Global and Local Errors in Runge-Kutta Methods

Manipulate[
Quiet[
x0=0;
poligonal[x0,xf,f,nn,y0,y];
Show[
Plot[Evaluate[y[t]/.teosol[x0,f,y0,t]],{t,x0,xf*1.01},PlotStyle->{Thickness[0.0025],RGBColor[0.1,0.,0.1]},Evaluate[plotoptions],AxesOrigin->{0,0},PlotRange->All,PlotLabel->Row[{"step size: ",(((xf-x0))/nn)//N}]],
ListPlot[Table[{numsol[[i]]},{i,1,Length[numsol]}],PlotStyle->Directive[PointSize[Medium],Blue]],

ListLinePlot[Tooltip[Table[{numsol[[i]],numsol[[i+1]]},{i,1,Length[numsol]-1}],"Numerical solution"],PlotStyle->RGBColor[0,0,1]],
(* local error *)
If[el,
{ListLinePlot[Tooltip[Table[{teo[[i]],solocal[[i+1]]},{i,1,Length[numsol]-1}],"Local numerical solution"],PlotStyle->RGBColor[1,0,0]],
ListLinePlot[Tooltip[Table[{{teo[[i,1]]-0.01,teo[[i,2]]},{solocal[[i,1]]-0.01,solocal[[i,2]]}},{i,1,Length[numsol]}],"local error"],PlotStyle-> Directive[Thickness[0.0025],RGBColor[1,0.1,0.1]]]},Graphics[{PointSize[Large],Point[{x0,y0}]}]]
,
(* global error *)
If[eg,ListLinePlot[Tooltip[Table[{{teo[[i,1]]+0.01,teo[[i,2]]},{numsol[[i,1]]+0.01,numsol[[i,2]]}},{i,1,Length[numsol]}],"global error"],
PlotStyle-> Directive[Thickness[0.0025],RGBColor[0,0.5,0]]],
Graphics[{PointSize[Large],Point[{x0,y0}]}]],
ImageSize->{400,400},ImagePadding->{{10,10},{10,40}}]],
{{f,f1,"y'"},{f1->"-y + cos(x)",f2->"-y",f3->"x + 1",f4->"cos(x) y",f5->"(y^2)(1 - x)",f6->"-5 y-2 x + 1",f7->"y (1 - y/10)",f8->"-y + U(x - 1)"},ControlType->PopupMenu,ControlPlacement->Left},
{{nn,4,"steps"},1,50,1,Appearance->"Labeled",ControlPlacement->Left,ImageSize->Tiny},
{{y0,1,"(Subscript[y, 0])"},1/4,3/2,1/4,Appearance->"Labeled",Exclusions->0,ControlPlacement->Left,ImageSize->Tiny},
{{xf,1,"(Subscript[x, e])"},1,10,1,Appearance->"Labeled",Exclusions->0,ControlPlacement->Left,ImageSize->Tiny},
{{method,eulerex,"method"},{eulerex->"Expl. Euler (order 1)",eulerex1->"Heun (order 2)",eulerex2->"RK2 (order 2)"},ControlType->PopupMenu,ControlPlacement->Left,ImageSize->Small},
Delimiter,
{{eg,False,"global error"},{True,False},ControlPlacement->Left},
{{el,False,"local error"},{True,False},ControlPlacement->Left},
{{numsol,1},ControlType->None},
{{teo,1},ControlType->None},
{{solocal,1},ControlType->None},
Delimiter,
Dynamic[If[eg,ListPlot[Table[{numsol[[i,1]],teo[[i,2]]-numsol[[i,2]]},{i,1,Length[numsol]}],Filling->Axis,PlotLabel->"global error",ImageSize->{200,100},ImagePadding->{{55,10},{15,10}}],Pane["no global error", ImageSize->{200,100}]],SynchronousUpdating->False],
Delimiter,
Dynamic[If[el,ListPlot[Table[{numsol[[i,1]],teo[[i,2]]-solocal[[i,2]]},{i,1,Length[numsol]}],Filling->Axis,PlotLabel->"local error",ImageSize->{200,100},ImagePadding->{{55,10},{15,10}}],Pane["no local error",ImageSize->{200,100}]],
SynchronousUpdating->False],ControlPlacement->Left,
AutorunSequencing->{1,4},
TrackedSymbols:>{f,nn,y0,xf,method,eg,el},
Initialization:>(
(* Differential Equations *)
f1[x_,y_]:=-y+Cos[x];
f2[x_,y_]:=-y;
f3[x_,y_]:= x+1;
f4[x_,y_]:= Cos[x] y;
f5[x_,y_]:= y^2(1-x);
f6[x_,y_]:= -5 y-2 x+1;
f7[x_,y_]:=   y(1-y/10);
f8[x_,y_]:=   -y+UnitStep[x-1];
(* Explicit Euler's method, order 1*)
eulerex[{a_,b_},f_,y0_,n_]:=Module[{h,x, y,lsg},
h=((b-a)/n);
x=Table[a+(i-1)h,{i,1,n+1}];
y=Table[0,{n+1}];
y[[1]]=y0;
Do [y[[i]]=y[[i-1]]+h f[x[[i-1]],y[[i-1]]],{i,2,n+1}];
lsg=Transpose[{x,y}]
];
(* Hein's method, order 2*)
eulerex1[{a_,b_},f_,y0_,n_]:=Module[{h,x, y,k,lsg},
h=((b-a)/n);
x=Table[a+(i-1)h,{i,1,n+1}];
y=Table[0,{n+1}];
y[[1]]=y0;
Do [k=f[x[[i-1]],y[[i-1]]];y[[i]]=y[[i-1]]+h f[x[[i-1]]+h/2,y[[i-1]]+h/2*k],{i,2,n+1}];
lsg=Transpose[{x,y}]
];
(* two-order, two-stage RK method *)
eulerex2[{a_,b_},f_,y0_,n_]:=Module[{h,x, y,k1,k2,k3,lsg},
h=((b-a)/n);
x=Table[a+(i-1)h,{i,1,n+1}];
y=Table[0,{n+1}];
y[[1]]=y0;
Do [k1=f[x[[i-1]],y[[i-1]]];y[[i]]=y[[i-1]]+h/2(k1+ f[x[[i-1]]+h,y[[i-1]]+h*k1]),{i,2,n+1}];
lsg=Transpose[{x,y}]
];
teosol[x0_,f_,y0_,t_]:=DSolve[{y'[t]==f[t,y[t]],y[x0]==y0},y,t];
plotoptions={ PlotRangePadding->0.10,ImageSize->475};poligonal[x0_,xf_,f_,nn_,y0_,y_]:=Module[{nn1,h,i,t,pp,aux},
nn1=N[nn];
h=(xf-x0)/nn1;
(* Compute the numerical solution at the mesh points *) 
numsol=method[{x0,xf},f,y0,nn1]//N;
teo=Table[{(i-1)*h,0},{i,1,Length[numsol]}]//N;
teo[[1]]={x0,y0};
(* Evaluate the exact solution at the mesh points *)
pp=Evaluate[y[t]/.teosol[x0,f,y0,t]];
For[i=1,i≤Length[numsol]-1,i++,
aux=pp/.t->numsol[[i+1,1]];
teo[[i+1,2]]=aux[[1]];
];
solocal=Table[{i*h,0},{i,0,Length[numsol]-1}];
solocal[[1,2]]=y0;
(* Compute the local solution *)
For[i=1,i≤Length[numsol]-1,i++,
aux=method[{teo[[i,1]],teo[[i+1,1]]},f,teo[[i,2]],1]//N;
solocal[[i+1,2]]=aux[[2,2]];
];
];)
]

CAPTION

This Demonstration shows the global and local errors generated by a one-step Runge–Kutta method in the numerical solution of initial value problems. The local errors are from the exact solution. For several initial value problems, you can select the number of steps, the initial value (Subscript[y, 0]), the end point of the interval of integration (Subscript[x, e]), and an explicit Runge–Kutta method with order less than or equal to 2. You can also plot the local and/or global errors.

DETAILS

For problems 2 and 6, some stability problems can appear when you vary the number of steps and the integration interval.

Problem 8 is a discontinuous initial value problem; U is the unit-step function, U(x)=0 for x<0 and U(x)=1 for x≥0. Therefore all the considered integration methods present some difficulties.

THIS NOTEBOOK IS THE SOURCE CODE FROM

"Global and Local Errors in Runge-Kutta Methods" from the Wolfram Demonstrations Project
 http://demonstrations.wolfram.com/GlobalAndLocalErrorsInRungeKuttaMethods/

Contributed by: Luis Rández

Based on a program by: Edda Eich-Soellner

A full-function Wolfram Mathematica system (Version 6 or higher) is required to edit this notebook.
GET WOLFRAM MATHEMATICA »

© Wolfram Demonstrations Project & Contributors   |   Terms of Use   |   Make a new version of this Demonstration »