Numerical Methods for Differential Equations

f1[x_,y_]:=y;
f2[x_,y_]:=-y;
f3[x_,y_]:= -x (y^2);
f4[x_,y_]:= (x^2)-y;
f5[x_,y_]:= Cos[x] y;
f6[x_,y_]:=(y^2)(x-(x^3));
f7[x_,y_]:= -20 y+20 x+21;

combPair[x_,y_]:={x,y}
combList[x_List,y_List]:=MapThread[combPair,{x,y}]

exsol[a_,f_,y0_,t_]:=DSolve[{y'[t]==f[t,y[t]],y[a]==y0},y,t]
plotoptions={ PlotRangePadding->0.25,ImageSize->400};

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=combList[x,y]
]

eulerimp[{a_,b_},f_,y0_,n_]:=Module[{h,x, y,lsg,sol},
h=((b-a)/n);
x=Table[a+(i-1)h,{i,1,n+1}];
y=Table[0,{i,1,n+1}];
y[[1]]=y0;
Do [lsg=Solve[sol==y[[i-1]]+h f[ x[[i]] ,sol],sol];
y[[i]]=If[lsg==={},y[[i-1]],sol/.lsg[[1]]],
{i,2,n+1}];
lsg=combList[x,y]
]

rkgenStep[{y0_,x0_,h_,f_},{A_,b_,c_}]:=
Module[{s,k},
s=Length[b];
k=Table[0,{i,1,s}];
k[[1]]=f[x0,y0];
Do[k[[i]]=f[ x0+c[[i]]h, y0+h  Sum[A[[i-1,j]] k[[j]],{j,1,i-1}] ],{i,2,s}];
ynext= y0+h Sum[b[[i]] k[[i]],{i,1,s}]
]

rkgen[{a_,r_},f_,y0_,n_,{A_,b_,c_}]:=Module[{h,x, y,lsg,i},
h=((r-a)/n);
x=Table[a+(i-1)h,{i,1,n+1}];
y=Table[0,{i,1,n+1}];
y[[1]]=y0;
Do[y[[i]]=rkgenStep[ {y[[i-1]],x[[i-1]],h,f},{A,b,c}],{i,2,n+1}];
lsg=combList[x,y]
]

rkheun[{a_,r_},f_,y0_,n_]:=Module[{A,b,c},
A={{(1/2)}};
b={0,1};
c={0,(1/2)};
rkgen[{a,r},f,y0,n,{A,b,c}]
]

rk4[{a_,r_},f_,y0_,n_]:=Module[{A,b,c},
A={{(1/2)},{0,(1/2)},{0,0,1}};
b={(1/6),(1/3),(1/3),(1/6)};
c={0,(1/2),(1/2),1};
rkgen[{a,r},f,y0,n,{A,b,c}]
]

Manipulate[Show[{ Plot[Evaluate[y[t]/.exsol[a,f,c,t]],{t,a,b+0.1},PlotStyle->Thickness[0.005],Evaluate[plotoptions],AxesOrigin->{0,0},PlotRange->All,PlotLabel->Row[{"step size: ",((b-a)/nn)}]],
ListPlot [Tooltip[Check[method[{a,b},f,c,If[MemberQ[{f3,f5,f6},f],N[nn],nn]],{a}]],PlotStyle->{ColorData[1,2],PointSize[0.02]}],
(*ListPlot[Tooltip[{a,b}],PlotStyle->{ColorData[1,2],PointSize[0.01]}],*)
ListLinePlot[Check[method[{a,b},f,c,If[MemberQ[{f3,f5,f6},f],N[nn],nn]],{a,c}],
PlotStyle->{ColorData[1,2]},PlotRange->All]
},
ImageSize->{400,400},ImagePadding->{{10,10},{10,40}}],
{{f,f1,"y'="},{f1->"y",f2->"-y",f3->"-x (y^2)",f4->"(x^2)-y",f5->"cos(x) y",f6->"(y^2)(x-(x^3))",f7->"-20y+20x+21"},ControlType->PopupMenu},
{{a,0},-1/2,1/2,1/4,Appearance->"Labeled",ControlPlacement->Left,ImageSize->Tiny},
{{b,2},1,3,1/4,Appearance->"Labeled",ControlPlacement->Left,ImageSize->Tiny},
{{c,1,"y0"},1/4,3/2,1/4,Appearance->"Labeled",Exclusions->0,ControlPlacement->Left,ImageSize->Tiny},
{{nn,3,"steps"},1,50,1,Appearance->"Labeled",ControlPlacement->Left,ImageSize->Tiny},
{{method,eulerex,"method"},{eulerex->"explicit Euler",eulerimp->"implicit Euler",rkheun->"Heun",rk4->"classical Runge-Kutta"}},
TrackedSymbols->Manipulate,
SaveDefinitions->True]

This Demonstration shows the exact and the numerical solutions using a variety of simple numerical methods for ordinary differential equations. Use the sliders to vary the initial value or to change the number of steps or the method.

The Demonstration shows various methods for ODEs:

* Euler's method is the simplest method for the numerical solution of an ordinary differential equation y'=f(x,y(x)). Starting from an initial point ((Subscript[x, 0]), (Subscript[y, 0])) and dividing the interval [a, b] that is under consideration into n steps results in a step size h=(b-a)/n; the solution value at point (Subscript[x, i]) is recursively computed using (Subscript[y, i+1])=(Subscript[y, i])+h f((Subscript[x, i]),(Subscript[y, i])), i=1,…,n.

* Implicit Euler method

* Heun's method

* classical Runge-Kutta method of order 4

The last right-hand side given belongs to a stiff equation, such that the behavior of the method for this type of equation can be studied. See M. Heath, Scientific Computing: An Introductory Survey, New York: McGraw-Hill, 2002.

Note that Mathematica provides all of the methods outlined here and many others as part of the NDSolve framework. In contrast to the simple implementations used here, Mathematica uses more advanced methods which are e.g. equipped with error estimation and step size selection strategies as well as a stiffness switching; see Mathematica's advanced documentation for NDSolve.

Contributed by: Edda Eich-Soellner