Euler Improved with “Runge One” routines

MAT360
Andy Long
Spring, 2017

Clear["Global `*"]

euler[f_,y0_,{a_,b_,n_}]:=
Module[
(* create some local variables: *)
{h=N[(b-a)/n],t=a,y=y0,soln},
(* initialize the solution list *)
soln={{t,y}};
(* iterate, updating y each time: *)
Do[y=y+h N[f[t,y]];
(* step time with a constant step-size: *)
t=t+h;
(* update soln: *)
AppendTo[soln,{t,y}
],
(* do this n times: *)
{n}
];
(* return the solution list of times and y values: *)
soln
]
(* Now we'll recreate Figure 8.11, p. 323 *)
f[t_,y_]:=1+ y^2
n=10;
(* Note: the text mentions that the points are connected by piecewise linear interpolation (bottom of page 323). You can change that using the InterpolationOrder option to ListPlot: *)
p0=ListPlot[euler[f,0,{0,N[Pi]/4,n}],Joined->True,InterpolationOrder->1];
p1=ListPlot[euler[f,0,{0,N[Pi]/4,n}],PlotMarkers->{Automatic,Medium}];
p2=Plot[Tan[x],{x,0,N[Pi]/4}];
Show[p0,p1,p2]

Note: we can use this Euler routine with a list of functions, however, so it handles systems! Lynx and Hares....

f[t_,y_]:=
{
2y[[1]]-0.01y[[1]]y[[2]],
-y[[2]]+0.01y[[1]]y[[2]]
}
n=100; (* Results are crazy for 100! We need a smaller h....*)
y0={300,150};
soln=euler[f,y0,{0,20,n}];
s=Transpose[soln[[Range[1,n],2]]];
ListPlot[{s[[1]],s[[2]]},AspectRatio->1.6, Joined-> True]
ListPlot[soln[[Range[1,n],2]],AspectRatio->1.6]

Note: we can turn second order DEs into systems of first-order DEs, so we can handle those, too, with this Euler routine. Harmonic oscillators:

f[t_,y_]:=
{
y[[2]],
-y[[1]]
}
n=100; (* Results get crazy a too small h.... we oscillate out of control!*)
y0={1,0};
soln=euler[f,y0,{0,2,n}];
s=Transpose[soln[[Range[1,n],2]]];
ListPlot[{s[[1]],s[[2]]},AspectRatio->1.6, Joined-> True]
ListPlot[soln[[Range[1,n],2]],AspectRatio->1.6]

Improving Euler: use the midpoint approximation, sort of!

Clear[t,halfway]
eulerRungeMid[f_,y0_,{a_,b_,n_}]:=
Module[
(* create some local variables: *)
{h=N[(b-a)/n],t=a,y=y0,soln,halfway},
(* initialize the solution list *)
soln={{t,y}};
(* iterate, updating y each time: *)
Do[halfway=y+h/2* N[f[t,y]];
(* step time with a constant step-size. But we follow Euler to what he thinks is the midpoint, compute the slope there, and then back up to follow that slope from our current position (rather than the slope from our starting point): *)
y=y+h N[f[t+h/2,halfway]];
t=t+h;
(* update soln: *)
AppendTo[soln,{t,y}
],
(* do this n times: *)
{n}
];
(* return the solution list of times and y values: *)
soln
]

(* Now we'll recreate Figure 8.11, p. 323 *)
f[t_,y_]:=1+ y^2
n=10;
(* Note: the text mentions that the points are connected by piecewise linear interpolation (bottom of page 323). You can change that using the InterpolationOrder option to ListPlot: *)
p0=ListPlot[eulerRungeMid[f,0,{0,N[Pi]/4,n}],Joined->True,InterpolationOrder->1];
p1=ListPlot[eulerRungeMid[f,0,{0,N[Pi]/4,n}],PlotMarkers->{Automatic,Medium}];
p2=Plot[Tan[x],{x,0,N[Pi]/4}];
Show[p0,p1,p2]

erunge=eulerRungeMid[f,0,{0,N[Pi]/4,n}];
erunge[[n+1]]
Tan[N[Pi]/4]

Improving Euler: use the trapezoidal approximation, sort of!

eulerRungeTrap[f_,y0_,{a_,b_,n_}]:=
Module[
(* create some local variables: *)
{h=N[(b-a)/n],t=a,y=y0,soln,eulerslope},
(* initialize the solution list *)
soln={{t,y}};
(* iterate, updating y each time: *)
Do[eulerslope= N[f[t,y]];
(* step time with a constant step-size. But we follow Euler to what he thinks is the midpoint, compute the slope there, and then back up to follow that slope from our current position (rather than the slope from our starting point): *)
y=y+h/2(eulerslope+ N[f[t+h,y+h*eulerslope]]);
t=t+h;
(* update soln: *)
AppendTo[soln,{t,y}
],
(* do this n times: *)
{n}
];
(* return the solution list of times and y values: *)
soln
]

(* Now we'll recreate Figure 8.11, p. 323, using this improved scheme: *)
p0=ListPlot[eulerRungeTrap[f,0,{0,N[Pi]/4,n}],Joined->True,InterpolationOrder->1];
p1=ListPlot[eulerRungeTrap[f,0,{0,N[Pi]/4,n}],PlotMarkers->{Automatic,Medium}];
p2=Plot[Tan[x],{x,0,N[Pi]/4}];
Show[p0,p1,p2]

erunge=eulerRungeTrap[f,0,{0,N[Pi]/4,n}];
erunge[[n+1]]
Tan[N[Pi]/4]

Improving Euler: Let’s compare to Taylor2, which is another way of improving Euler (which we can think of as Taylor1):

Taylor2[f_,y0_,{a_,b_,n_}]:=
Module[
(* create some local variables: *)
{h=N[(b-a)/n],t=a,y=y0,soln,eulerslope},
(* initialize the solution list *)
soln={{t,y}};
(* iterate, updating y each time: *)
Do[y=y+h(1+y^2 +h/2(2 y (1+y^2)));
(* step time with a constant step-size *)
t=t+h;
(* update soln: *)
AppendTo[soln,{t,y}
],
(* do this n times: *)
{n}
];
(* return the solution list of times and y values: *)
soln
]

(* Now we'll recreate Figure 8.11, p. 323, using this improved scheme: *)
p0=ListPlot[Taylor2[f,0,{0,N[Pi]/4,n}],Joined->True,InterpolationOrder->1];
p1=ListPlot[Taylor2[f,0,{0,N[Pi]/4,n}],PlotMarkers->{Automatic,Medium}];
p2=Plot[Tan[x],{x,0,N[Pi]/4}];
Show[p0,p1,p2]

erunge=Taylor2[f,0,{0,N[Pi]/4,n}];
erunge[[n+1]]
Tan[N[Pi]/4]