210 6 Linear Algebra
ofstream ofile;
// Main program only, no other functions
intmain(intargc,charargv[])
{
charoutfilename;
inti, j, n;
doubleh, btemp;
// Read in output file, abort if there are too few command-line arguments
if( argc <= 1 ){
cout <<"Bad Usage: "<< argv[0] <<
" read also output file on same line"<< endl;
exit(1);
}
else{
outfilename=argv[1];
}
ofile.open(outfilename);
cout <<"Read in number of mesh points"<< endl;
cin >> n;
h = 1.0/( (double) n+1);
// Use Blitz to allocate arrays
// Use range to change default arrays from 0:n-1 to 1:n
Range r(1,n);
Array<double,1> a(r), b(r), c(r), y(r), f(r), temp(r);
// set up the matrix defined by three arrays, diagonal, upperand lower diagonal band
b = 2.0; a = -1.0 ; c = -1.0;
// Then define the value of the right hand side f (multiplied by hh)
for(i=1; i <= n; i++){
// Explicit expression for f, could code as separate function
f(i) = hh(ih3.0+(ih)(ih))exp(ih);
}
// solve the tridiagonal system, first forward substitution
btemp = b(1);
for(i = 2; i <= n; i++){
temp(i) = c(i-1) / btemp;
btemp = b(i) - a(i)temp(i);
y(i) = (f(i) - a(i)y(i-1)) / btemp;
}
// then backward substitution, the solution is in y()
for(i = n-1; i >= 1; i--){
y(i) -= temp(i+1)y(i+1);
}
// write results to the output file
for(i = 1; i <= n; i++){
ofile << setiosflags(ios::showpoint | ios::uppercase);
ofile << setw(15) << setprecision(8) << ih;
ofile << setw(15) << setprecision(8) << y(i);
ofile << setw(15) << setprecision(8) << ih(1.0-ih)exp(i*h) <<endl;
}
ofile.close();
}
The program writes also the exact solution to file. In Fig. 6.4we show the results obtained
withn= 10. Even with so few points, the numerical solution is very close to the analytic
answer. Withn= 100 it is almost impossible to distinguish the numerical solution from the
analytical one, as shown in Fig. 6.5. It is therefore instructive to study the relative error,
which we display in Table 6.2 as function of the step lengthh= 1 /(n+ 1 ).
The mathematical truncation we made when computing the second derivative goes like
O(h^2 ). Our results fornfromn= 10 to somewhere betweenn= 104 andn= 105 result in a
slope which is almost exactly equal 2 ,in good agreement with the mathematical truncation