Comparing Leapfrog Methods with Other Numerical Methods for Differential Equations

f1[x_,y_]:=y;
f2[x_,y_]:=y-2x;
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[

(* leapfrog method *)
(* second order explicit method *)
(* one f-evaluation per time step *)
lf[{a_,b_},f_,y0_,n_]:=Module[{h,h2,x,y},
h=((b-a)/n); h2=2 h;
x=Table[0,{n+1}];y=Table[0,{n+1}]; 
x[[1]]=a; y[[1]]=y0;
x[[2]]=a+h; y[[2]]=y0+f[a,y0] h;
Do [
x[[i]]=x[[i-2]]+h2; 
y[[i]]=y[[i-2]]+h2 f[x[[i-1]],y[[i-1]]],
{i,3,n+1}
];
combList[x,y]
];

(* densified leapfrog method *)
(* second order explicit method *)
(* two f-evaluations per time step *)
dlf[{a_,b_},f_,y0_,n_]:=Module[{h,hh,xh,yh,xp,yp,x,y},
h=((b-a)/n);hh=h/2;
x=Table[0,{n+1}];y=Table[0,{n+1}]; 
x[[1]]=a; y[[1]]=y0;
xh=x[[1]]+hh; yh=y[[1]]+hh f[ x[[1]],y[[1]] ] ;
x[[2]]=x[[1]]+h; y[[2]]=y[[1]]+h f[xh,yh];
Do [
xp=xh+h; yp=yh+h f[x[[i-1]],y[[i-1]]];
x[[i]]=xp+hh; y[[i]]=y[[i-1]]+h f[xp,yp]; 
xh=xp; yh=yp,
{i,3,n+1}
];
combList[x,y]
];

(* asynchronous leapfrog method *)
(* second order explicit method *)
(* one f-evaluation per time step *)
(* additional data (velocity v) and one additional f-evaluation for *)
(* the first time step *)
alf[{a_,b_},f_,y0_,n_]:=Module[{τ,x,y,v,xi,yi,vi},
τ=((b-a)/(2*n));
x=Table[0,{n+1}];y=Table[0,{n+1}]; v=Table[0,{n+1}];
x[[1]]=a; y[[1]]=y0; v[[1]]=f[a,y0];(* additional data v and f-evaluation *)
Do [
xi=x[[i-1]]; yi=y[[i-1]];vi=v[[i-1]];
xi+=τ;yi+=τ*vi;
vi+=2*( f[xi,yi]-vi);
xi+=τ;yi+=τ*vi;
x[[i]]=xi;y[[i]]=yi;v[[i]]=vi,
(* inelegant data re-shuffling *)
{i,2,n+1}];
 combList[x,y]
];

(* densified asynchronous leapfrog method *)
(* second order explicit method *)
(* two f-evaluations per time step *)
(* additional data (velocity v) and one additional f-evaluation for *)
(* the first time step *)
dalf[{a_,b_},f_,y0_,n_]:=Module[{τ,τ2,x,y,v,xi,yi,vi},
τ=((b-a)/(4*n)); τ2=2 τ;
x=Table[0,{n+1}];y=Table[0,{n+1}]; v=Table[0,{n+1}];
x[[1]]=a; y[[1]]=y0; v[[1]]=f[a,y0];
Do [
xi=x[[i-1]]; yi=y[[i-1]];vi=v[[i-1]];
xi+=τ;yi+=τ*vi;
vi+=2*( f[xi,yi]-vi);
xi+=τ2;yi+=τ2*vi;(* two τ-steps contracted *)
vi+=2*( f[xi,yi]-vi);
xi+=τ;yi+=τ*vi;
x[[i]]=xi;y[[i]]=yi;v[[i]]=vi,
{i,2,n+1}];
  combList[x,y]
];
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]}],
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,f4,"y' ="},{f1->"y",f2->"y-2x",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,"(Subscript[y, 0])"},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,dalf,"method"},{
eulerex->"explicit Euler",
eulerimp->"implicit Euler",
rkheun->"Heun",
rk4->"classical Runge-Kutta",
lf->"leapfrog",
alf->"aynchronous leapfrog",
dlf->"densified leapfrog",
dalf->"densified asynchronous leapfrog"
}},
TrackedSymbols->True,
AutorunSequencing->{6},
SaveDefinitions->True]

CAPTION

This Demonstration shows the exact and the numerical solution of some ODEs using a variety of numerical methods. Four of these methods are well-known simple standard methods. Additional methods are the well-known leapfrog method and the less-known asynchronous leapfrog method. Further, a straightforward derivative of each of these two methods is considered. Use the dropdown menus to vary the differential equation and the method. Use the sliders to vary the initial value or to change the number and size of the integration steps.

DETAILS

In comparing the performance of numerical methods for differential equations, one has to take into account the computational complexity of the underlying algorithms. A simple measure (and for many practical applications the most relevant one) counts how often the right-hand side of the differential equation has to be evaluated in one integration step. The leapfrog and asynchronous leapfrog methods excel in needing only one such evaluation and thus compare directly only with the explicit Euler method, which they outperform in a spectacular way. Since the leapfrog methods are second-order methods, they are most naturally compared with the second-order Runge–Kutta method, which is also known as Heun's method. Since Heun's method does two right-hand side evaluations, it is fair to compare one Heun step with a composition of two leapfrog steps of half the step size. To allow such a fair comparison without having to change the step size, such combined half-steps are implemented. They are referred to as densified methods.

Snapshots 1 and 2: behavior for a stiff differential equation; here, the Runge–Kutta methods are superior to the leapfrog methods

Snapshots 3 and 4: behavior for non-stiff differential equations; the densified asynchronous leapfrog method seems to be in most cases more accurate than Heun's method

The leapfrog methods were added by the present author to a program by Edda Eich-Soellner (see Related Links). Apart from additions, virtually no changes have been made to the program of this instructive Demonstration.

References

[1] U. Mutze. "An Asynchronous Leap-Frog Method." (2008) http://www.ma.utexas.edu/mp_arc/c/08/08-197.pdf.

[2] Alternate link for [1]: http://www.ulrichmutze.de/articles/leapfrog4.pdf.

THIS NOTEBOOK IS THE SOURCE CODE FROM

"Comparing Leapfrog Methods with Other Numerical Methods for Differential Equations" from the Wolfram Demonstrations Project
 http://demonstrations.wolfram.com/ComparingLeapfrogMethodsWithOtherNumericalMethodsForDifferen/

Contributed by: Ulrich Mutze

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 »