Taylor Expansions about t=0

yhTaylor[n_] = y[0]+Sum[  (D[f[t],{t,k-1}] /.t->0)((h^k)/(k!))  , {k,n}] (* y(h) in terms of f(t) *)

yhTaylor[5]

ftTaylor[t_,n_]:=Sum[  (D[f[τ],{τ,k}] /.τ->0)((t^k)/(k!))  , {k,0,n}]  (* f(t) expansion *)

ftTaylor[h,5]

Forward Euler O((h^2))

yhFwdEuler=y[0]+h f[0]

yhFwdEuler-yhTaylor[2]

Backward Euler O((h^2))

yhBwdEuler = Expand[y[0]+h ftTaylor[h,1]]

yhBwdEuler-yhTaylor[2]

Trapezoidal Method O((h^3))

ftTaylor[h,2]

yhTrapezoid=Expand[y[0]+(h/2)(f[0]+ftTaylor[h,2])]

Series[  yhTrapezoid - yhTaylor[5], {h,0,3}]

Midpoint Method O((h^3))

yhMidpoint = Expand[y[0] +h ftTaylor[(h/2),5]]

Series[ yhMidpoint-yhTaylor[5], {h,0,3}]

Runge-Kutta Methods

Use degree-n Taylor series expansion of F[y] about (Subscript[y, 0]):

Fy0Taylor[Y_,n_]:=F[(Subscript[y, 0])]+Sum[  (D[F[X],{X,k}] /.X->(Subscript[y, 0]))((((Y-(Subscript[y, 0])))^k)/(k!))  , {k,n}]

Fy0Taylor[(Subscript[y, 0]),4]

Fy0Taylor[(Subscript[y, 0])+h v,4]

Expanding (y^((n)))(t) derivatives in terms of F and its gradients

dyF[1]=F[v]  (* (y^((1))) *)

dyF[2] = D[dyF[1],v] F[v]  (* (y^((2))) *)

dyF[3] = D[dyF[2],v] F[v]   (* (y^((3))) *)

dyF[4] = D[dyF[3],v] F[v]   (* (y^((4))) *)

dyF[5] = D[dyF[4],v] F[v]   (* (y^((5))) *)

dyF[6] = D[dyF[5],v] F[v]   (* (y^((6))) *)

Taylor series of y(h) about zero using F & its gradients for (y^((n)))(t=0) derivatives:

yhFExact[n_] := ((Subscript[y, 0]) + Sum[dyF[i] ((h^i)/(i!)), {i,1,n}] )

yhF[n_]:=yhFExact[n]/.v->(Subscript[y, 0])  (* Taylor series expansion of y(h) about t=0, y=(Subscript[y, 0]) *)

yhF[3]  (* test *)

Note: For each Runge-Kutta integrator, we will now compare its Taylor expansion (using F & gradients) to yhF[n].

Fehlberg3 / Heun’s method O((h^3))

yhF3[n_]:=(Subscript[y, 0])+(h/2)(F[(Subscript[y, 0])]+Fy0Taylor[(Subscript[y, 0])+h F[(Subscript[y, 0])],n-1])

Expand[yhF3[4]]

Series[yhF3[4]-yhF[4], {h,0,3}]  (* Verified: O((h^3)) accurate *)

Fehlberg4 O((h^4))

yhF4[n_]:=(Subscript[y, 0])+(h/6)(F[(Subscript[y, 0])]+Fy0Taylor[(Subscript[y, 0])+h F[(Subscript[y, 0])],n-1]+4 Fy0Taylor[(Subscript[y, 0])+(h/4)(F[(Subscript[y, 0])]+Fy0Taylor[(Subscript[y, 0])+h F[(Subscript[y, 0])],n-2]) ,n-1])

Series[yhF4[4]-yhF[4], {h,0,4}]  (* Verified: O((h^4)) accurate *)

RK4 O((h^5))

k1[n_]:=Fy0Taylor[(Subscript[y, 0]),n]

k2[n_]:=Fy0Taylor[(Subscript[y, 0])+(h/2)k1[n-1],n]  (* O((h^n)) accurate expansion -- NB: Only O((h^(n-1))) need for inner gradient since multiplied by h *)

k3[n_]:=Fy0Taylor[(Subscript[y, 0])+(h/2)k2[n-1],n]

k4[n_]:=Fy0Taylor[(Subscript[y, 0])+h k3[n-1],n]

k1[3]

k2[3]

k3[3]

k4[3]

yhRK4 = (Subscript[y, 0])+(h/6)(k1[4]+2 k2[4]+2 k3[4] + k4[4])  (* O((h^5)) accurate RK4 expansion *)

Simplify[Series[yhRK4-yhF[5],{h,0,5}]]  (* Verified: O((h^5)) accurate!! *)

Stability of Trapezoid Method on Model Equation (x=ha)

Rtrap[x_]:=Abs[((1+x/2)/(1-x/2))]

Plot[{Rtrap[x],1},{x,-10, 2}]

Stability of Heun’s Method on Model Equation (x=ha)

RHeun[x_]:=Abs[1+x+(1/2)(x^2)]

Plot[{RHeun[x],1},{x,-3,3}]