ODE example for MAT360

Andy Long, Spring, 2024

Tea Time Numerical Analysis, Taylor Methods, starting on p. 209
Leave it to Leo to go backwards in time in his first example...:)

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 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 6.2.2, p. 210  *)
f[t_,y_]:=-y/t+t^2
y0=20;
n=8;
a=4;
b=2;
g[x_]:=x^3/4+16/x (* exact solution *)

soln=euler[f,y0,{a,b,n}];
MatrixForm[soln]
p0=ListPlot[soln,Joined->True,InterpolationOrder->1];
p2=Plot[g[x],{x,a,b}];
p1=ListPlot[soln];
Show[p0,p1,p2]

taylor2[f_,fp_,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( f[t,y]+h/2*fp[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
]

f[t_,y_]:=-y/t+t^2
(* This calculation is done on page 209 *)
fp[t_,y_]:=2y/t^2+t
y0=20;
n=8;
a=4;
b=2;
g[x_]:=x^3/4+16/x (* exact solution *)

soln=taylor2[f,fp,y0,{a,b,n}];
MatrixForm[soln]
p0=ListPlot[soln,Joined->True,InterpolationOrder->1];
p2=Plot[g[x],{x,a,b}];
p1=ListPlot[soln];
Show[p0,p1,p2]

taylor3[f_,fp_,fpp_,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( f[t,y]+h/2*(fp[t,y]+h/3fpp[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 Table 6.2, p. 212 *)
f[t_,y_]:=-y/t+t^2
fp[t_,y_]:=2y/t^2+t
(*  I did the following calculation, and hope that it's right...:) *)
fpp[t_,y_]:=2(t^3-3y)/t^3+1
y0=20;
n=8;
a=4;
b=2;
g[x_]:=x^3/4+16/x (* exact solution *)

soln=taylor3[f,fp,fpp,y0,{a,b,n}];
MatrixForm[soln]
p0=ListPlot[soln,Joined->True,InterpolationOrder->1];
p2=Plot[g[x],{x,a,b}];
p1=ListPlot[soln];
Show[p0,p1,p2]

#1-3, ac

f[t_,y_]:=3t-2y
fp[t_,y_]:=3-2(3t-2y)
(*  I did the following calculation, and hope that it's right...:)
	Notice that I cheat and don't simplify y', leaving it as f[t,y]. *)
fpp[t_,y_]:=-6+4f[t,y]
y0=1;
n=2;
a=1;
b=2;
soln=euler[f,y0,{a,b,n}];
MatrixForm[soln]
soln=taylor2[f,fp,y0,{a,b,n}];
MatrixForm[soln]
soln=taylor3[f,fp,fpp,y0,{a,b,n}];
MatrixForm[soln]

It’s easy to operate on systems:

Here’s the pendulum problem:

Clear[x,y,z,f,c,m,g,l]
(* We're thinking of f as a vector-valued function, obtained by reducing the second order equation for the pendulum to a first order system. The arguments a time t and a location (u,v): *)   
f[t_,{u_,v_}]:={
-c/m u-g/l Sin[v],
u
}
g=9.81; (* m/s2 *)
c=1.1;
l=.31; (* m *)
m=70; (* kg *)
theta0=Pi/3;
thetaprime0=0;

(* We just give an initial value that is also a list, of the proper size: *)
soln=euler[f,{thetaprime0,theta0},{0,10,40}];
curve=Table[soln[[k]][[2]],{k,1,Length[soln]}];
Show[Graphics[{Line[curve]}],Axes->Automatic]

Let’s try that again, with Taylor 3:

soln=taylor3[f,fp,fpp,{thetaprime0,theta0},{0,10,40}];
curve=Table[soln[[k]][[2]],{k,1,Length[soln]}];
Show[Graphics[{Line[curve]}],Axes->Automatic]

Clear[g]

Here is the lovely Lorenz System, exhibiting chaos:

Clear[x,y,z,f]
f[t_,{x_,y_,z_}]:={
-3.(x-y),
-x z+26.5 x-y,
x y-z
}
soln=euler[f,{0.0,1.0,0.0},{0,50,40000}];
curve=Table[soln[[k]][[2]],{k,1,Length[soln]}];
Show[Graphics3D[{Line[curve]}],Axes->Automatic]

soln=taylor3[f,fp,fpp,{0.0,1.0,0.0},{0,50,40000}];
curve=Table[soln[[k]][[2]],{k,1,Length[soln]}];
Show[Graphics3D[{Line[curve]}],Axes->Automatic]