ODE examples for MAT360
Andy Long
Spring, 2015

Figure 8.11, p. 323:

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=4;
(* 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]

Computer problem #1, p. 330 -- redo figure 8.11, with more steps:

n=8;
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}];
Show[p0,p1,p2]

n=16;
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}];
Show[p0,p1,p2]

n=32;
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}];
Show[p0,p1,p2]

Computer problem #2, p. 330

f[t_,y_]:=Sin[t] y
y0=1
{ListPlot[euler[f,y0,{0,20,20}],Joined->True,PlotRange->{{0,20},{0,7}}],
ListPlot[euler[f,y0,{0,20,200}],Joined->True],
ListPlot[euler[f,y0,{0,20,2000}],Joined->True]}
p0=Plot[Exp[1]Exp[-Cos[x]],{x,0,20},PlotStyle->Red];
p1=ListPlot[euler[f,y0,{0,20,20}],Joined->True,PlotRange->{{0,20},{0,8}}];
p2=ListPlot[euler[f,y0,{0,20,200}],Joined->True];
p3=ListPlot[euler[f,y0,{0,20,2000}],Joined->True];
Show[p0,p1,p2,p3,PlotRange->{{0,20},{0,8}}]

Systems of ODE: Rabbits and foxes. First using “NDSolve” (a black box, at this point), and then using Euler’s method with NestList.

Clear[r,f,t];
sol=NDSolve[
{
r'[t]==2r[t]-0.01r[t]f[t],
f'[t]==-f[t]+0.01r[t]f[t],
r[0]==300,
f[0]==150},
{r,f},
{t,20}]
ParametricPlot[Evaluate[{r[t],f[t]}/.sol],{t,0,20}]
Plot[{Evaluate[r[t]/.sol],Evaluate[f[t]/.sol]},{t,0,20}]

Computer problem #1, p. 337:

(* This is Euler's method on a system of ODEs -- a piece of cake! *)
step[{t_,r_,f_}]:={t+h,r+h r (2.-0.01 f),f+h f (0.01 r-1.)} 
h=0.01;
n=Floor[20/h];
soln=NestList[step,{0.,300.,150.},n];
ListPlot[soln[[Range[1,n],{2,3}]],AspectRatio->1.6]

(* We could also use our old Euler plot routine with a list of functions, however: : *)
f[t_,y_]:=
{
2y[[1]]-0.01y[[1]]y[[2]],
-y[[2]]+0.01y[[1]]y[[2]]
}
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]

Computer problem #2, p. 337

(*part a*)
solveode[f_,yinit_,ypinit_,tfinal_,h_]:=
Module[{n=Ceiling[tfinal/h],y=yinit,yp=ypinit,ypp},ypp=f[yinit];
{{0,y,yp}}~Join~Table[yp=yp+0.5 h ypp;
y=y+h yp;
ypp=f[y];
yp=yp+0.5 h ypp;
{(i-1)h,y,yp},{i,2,n+1}]]
f[y_]:=-y
solveode[f,1.,0.,6. N[Pi],0.1];
ListPlot[%[[Range[Length[%]],{2,3}]]]
(*part b*)
f[y_]=.
f[r_]:=34.7*(3.4/r)^7*(2*(3.4/r)^6-1)
soln=solveode[f,6.46,0,12.64,0.05];
Map[#[[{1,2}]]&,soln];
ListPlot[%,PlotLabel->"waveform"]
(*part c*)
Map[#[[{2,3}]]&,soln];
ListPlot[%,PlotLabel->"trajectory"]

The Lorenz equations provide us with a system exhibiting chaos:

Clear[s,x,y,z,t]
s=NDSolve[
{
x'[t]==-3 (x[t]-y[t]),
y'[t]==-x[t] z[t]+26.5 x[t]-y[t],
z'[t]==x[t] y[t]-z[t],
x[0]==z[0]==0,
y[0]==1
},{x,y,z},{t,0,200},MaxSteps->∞]
ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.s],{t,0,200},PlotPoints->1000,ColorFunction->(Hue[#4]&)]