Comparing total error against truncation error at each grid point resulting from the numerical solution of u''(x)=exp(x) using finite difference discretization

Nasser M. Abbasi, 10/25/2010

Introduction

The differential equation u''[x]==(e^x) with Dirichlet boundary conditions u[0]=0,u[1]=0 is solved numerically by approximating u''[x] by (((Subscript[u, i-1])-2(Subscript[u, i])+(Subscript[u, i+1]))/(h^2))where h is the spacing between grid points. The truncation error τ at each (Subscript[x, i]) is given by ((h^2)/12)(u^((4)))((Subscript[x, i]))while the total error e at the same point is given by (Subscript[U, i])-(Subscript[u, i]) where (Subscript[U, i]) is the numerical solution found at that point and (Subscript[u, i]) is the exact solution u(x) evaluated at the same point.
A small program is written to show the distribution of e and τ along the length of the domain [0,1] as h is made smaller and smaller.
This was run for few iterations, where h was divided by half each time, the system was solved for U and the result is shown in a series of plots.

Results

Get["G:\\nabbasi\\data\\nabbasi_web_Page\\my_notes\\my_dsp_mathematica_package\\dsp.m"];
Get["G:\\nabbasi\\data\\mma_package\\plot_custom_ticks\\LevelScheme-3.45\\CustomTicks.m"];

matrixLaplace1DRobin[n_,h_,α_]:=Module[{A,i,j},
A=Table[If[i==j,-2,If[i==j||i==j-1,1,0]],{i,n},{j,n}];
A[[1,1]]=-1-α*h;
A
];
matrixLaplace1DNeumannScheme1[n_]:=Module[{A,i,j},
A=Table[If[i==j,-2,If[i==j+1||i==j-1,1,0]],{i,n},{j,n}];
A[[1,2]]=2;
A[[-1,-2]]=2;
A
];
matrixLaplace1DNeumannScheme2[n_]:=Module[{A,i,j},
A=Table[If[i==j,-2,If[i==j+1||i==j-1,1,0]],{i,n},{j,n}];
A[[1,1]]=3/2;
A[[1,2]]=-2;
A[[1,3]]=1/2;
A[[-1,-1]]=3/2;
A[[-1,-2]]=-2;
A[[-1,-3]]=1/2;

A
];
matrixLaplace1DDirichletn[n_]:=Module[{A,i,j},
A=Table[If[i==j,-2,If[i==j+1||i==j-1,1,0]],{i,n},{j,n}]
];

errorVectorProblem1Parta[points_,h_]:=Module[{U,u,e,xcoordinates,α,β,len,A,f,force,exact,numberInternalPoints,x,τ},
force[x_]:=Exp[x];
exact[x_]:=-1+Exp[x]+2*x-Exp[1]*x;
truncationError[x_,spacing_]:=((Exp[x])/12)*spacing^2;  (* this is - u''''(x)/12*)
α=0;
β=1;
len=1;

numberInternalPoints=points-2;
xcoordinates=Range[h,len-h,h];(*xcoordinates at internal points*)

A=matrixLaplace1DDirichletn[numberInternalPoints];
A=(A/(h^2));

(*build the f vector,make sure to adjust the first and last entries*)
f=force[xcoordinates];
f[[1]]=f[[1]]-(α/(h^2));
f[[-1]]=f[[-1]]-(β/(h^2));
τ=truncationError[xcoordinates,h];
U=LinearSolve[A,f];
u=exact[xcoordinates];(*% find exact solution at the grid points*)e=U-u;
N[{e,U,u,xcoordinates,τ}]
];

prob1Parta[]:=Module[{nIter,points,h,e,U,u,x,i,imSize=250,τ,POINTS=1,SPACE=2,EMAX=3,EMAXRATIO=4,EONE=5,EONERATIO=6,PLOTS=7,TAOINF=8,EINF=9,data,p1,p2,p3,p4,tbl},

nIter=7;
points=Table[(2^i)+1,{i,nIter}];
h=N@Table[(1/(2^i)),{i,nIter}];
data=Table[0,{i,nIter},{j,9}];


Do[{e,U,u,x,τ}=errorVectorProblem1Parta[points[[i]],h[[i]]];
data[[i,POINTS]]=points[[i]];
data[[i,SPACE]]=h[[i]];
data[[i,EMAX]]=Norm[e,Infinity];
data[[i,EONE]]=h[[i]] Norm[e,1];
If[i>1,
{
data[[i,EMAXRATIO]]=((data[[i-1,EMAX]])/(data[[i,EMAX]]));
data[[i,EONERATIO]]=((data[[i-1,EONE]])/(data[[i,EONE]]));
}];
data[[i,TAOINF]]=(Sqrt[h[[i]]])Norm[τ,2];
data[[i,EINF]]=(Sqrt[h[[i]]])Norm[e,2];


p1=ListPlot[Inner[List,x,e,List],Joined->False,AxesOrigin->{0,0},PlotRange->All,ImageSize->imSize,Filling->Axis,Frame->True,BaseStyle->{FontFamily->"Times"},ImagePadding->{{60,10},{40,40}},FrameLabel->{{None,None},{"x","total error e"}}(*FrameTicks-> {CustomTicks`LinTicks,CustomTicks`LinTicks,None,None}*)];
p2=ListPlot[Inner[List,x,τ,List],Joined->False,AxesOrigin->{0,0},PlotRange->All,ImageSize->imSize,Filling->Axis,Frame->True,BaseStyle->{FontFamily->"Times"},ImagePadding->{{60,10},{40,40}},FrameLabel->{{None,None},{"x","truncation error τ"}}(*FrameTicks-> {CustomTicks`LinTicks,CustomTicks`LinTicks,None,None}*)];

data[[i,PLOTS]]=Labeled[Framed@Grid[{{p1,p2}},Spacings->2],Style[Row[{"number of grid points = ",points[[i]],", spacing h=",h[[i]]}],Bold],Top];
,{i,nIter}
];


tbl=data[[All,1;;6]];
For[i=1,i≤nIter,i++,{
tbl[[i,2]]=dsp`numItse[tbl[[i,2]],10,8,12];
tbl[[i,3]]=dsp`numItse[tbl[[i,3]],10,8,12];
tbl[[i,4]]=dsp`numItse[tbl[[i,4]],10,8,12];
tbl[[i,5]]=dsp`numItse[tbl[[i,5]],10,8,12];
tbl[[i,6]]=dsp`numItse[tbl[[i,6]],10,8,12];
}];

p1=ListPlot[Inner[List,x,u,List],Joined->True];
p2=ListPlot[Inner[List,x,U,List],PlotStyle->Red];

Do[ Print@data[[i,7]],{i,3,7}];

p3=ListPlot[Inner[List,h,data[[All,8]],List],Joined-> True,Frame->True,PlotStyle->Red,FrameLabel->{{None,None},{"h",None}}];
p4=ListPlot[Inner[List,h,data[[All,9]],List],Joined-> True,Frame->True];
Print@Labeled[Show[p3,p4],Style["||e(Subscript[||, 2]) compared to ||τ(Subscript[||, 2]) linear scale",Bold],Top];

p3=ListPlot[Inner[List,Log[10,h],Log[10,data[[All,8]]],List],Joined-> True,Frame->True,PlotStyle->Red,FrameLabel->{{None,None},{"h",None}},FrameTicks-> {CustomTicks`LogTicks,CustomTicks`LogTicks,None,None}];
p4=ListPlot[Inner[List,Log[10,h],Log[10,data[[All,9]]],List],Joined-> True,Frame->True,FrameTicks-> {CustomTicks`LogTicks,CustomTicks`LogTicks,None,None}];
Print@Labeled[Show[p3,p4],Style["||e(Subscript[||, 2]) compared to ||τ(Subscript[||, 2]) log scale",Bold],Top];


Grid[{

{
TableForm[tbl,TableHeadings->{None,{"N","h","e-max","ratio","e-1","ratio"}},TableAlignments->Center]
},

{
ListPlot[Inner[List,Log[10,data[[All,2]]],Log[10,data[[All,5]]],List],Joined->True,
FrameTicks-> {CustomTicks`LogTicks,CustomTicks`LogTicks,None,None},
Frame->True,
BaseStyle->{FontFamily->"Times"},
GridLines->Automatic,
FrameLabel->{{Style["log(error)",10],None},{Style["log(h)",10],Style["error accuracy by slope method",12]}},ImageSize->300
]
}

}]
]

prob1Parta[]

Observations on the above plots

The total error e was largest in the middle of the grid, and smallest at the edges, while the truncation error was smallest at the left edge and largest at the right edge. 
Comparing e to τ on each grid point is useful to understand the distribution of errors along the domain, but for verification that the numerical solution converges to the exact solution as h becomes smaller, the grid norm of e is the one examined to verify that it is less than or equal to the grid norm of τ. 
Therefore, the following plots are generated which shows how ||e(Subscript[||, 2]) compares to ||τ(Subscript[||, 2])  as h becomes smaller. These norms are grid norms and not the standard norms.

Observations on the above plots

The total error has smaller grid norm. This shows convergence, since this implies that ||e||≤||(A^(-1))|| ||τ||, therefore showing that ||(A^(-1))||=O(1) which is the condition for convergence. The above results confirms this for this problem.