Code snippets for ODE Driver/Stepper

EulerStep[y_,f_,h_]:= Module[
{k1,k2,k3,k4,yNew,zNew},
k1=f[y];
yNew=y + h k1;
k2=f[yNew];
{yNew, Norm[h k1-h k2]}
]  

ODE23Step[y_,f_,h_]:= Module[
{k1,k2,k3,k4,yNew,zNew},
k1=f[y];
k2=f[y+ 1/2 h k1];
k3=f[y+3/4 h k2];
yNew=y + h ( 2/9 k1 + 1/3 k2 + 4/9 k3);
k4=f[yNew];
zNew=y +h*(7/24 k1 + 1/4 k2 + 1/3 k3 + 1/8 k4);
{yNew, Norm[zNew-yNew]}
]  

DriverVariableStep1[y0_,h0_,f_,n_, ErrorOrder_,LocalTol_,Stepper_]:=Module[
{y=y0,t=0.0,h=h0,err,Data, SafetyFactor=0.9, yNew},
  Data = Table[
{yNew, err} = Stepper[y, f, h];
If[err<LocalTol, {t,y}={t+h,yNew}, Print["Failed Step"]];
h=h*SafetyFactor*((((LocalTol/err)))^(1/ErrorOrder));
{t, y},
{n}];
(* fixes it so that it comes out in the form {ts, ys} *)
Transpose[Data]
]

Testing on our gravitational problem

Testing this code on the gravitaional problem

f1[{y1_,y2_,y3_,y4_,y5_,y6_}] :=Module[
{μ=1/81.45,d1,d2 },
d1 =((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2)); 
d2=((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2));
{y4,y5,y6,
2y5 + y1 - ((μ(y1+μ-1))/d1) -(((1-μ)(y1+μ))/d2),
-2y4 + y2 - ((μ y2)/d1) -(((1-μ)y2)/d2),
- ((μ y3)/d1) -(((1-μ)y3)/d2)
}
]
IC = {0.994,0,0,0,-2.001585106370825224,0};

n=1000;
h0=0.001;
LocalTol=(10^(-3));
ErrorOrder =1;
{tVs, yVs}=DriverVariableStep1[IC,h0,f1,n,ErrorOrder,LocalTol,EulerStep];
ListPlot[ yVs[[All,{1,2}]],
PlotStyle->{Red, Green},
PlotRange->All,
AxesLabel->{"y1","y2"}]
ListLogPlot[Drop[RotateLeft[tVs]-tVs,-1],
PlotRange->All]

B2=Simplify[D[f1[{y1,y2,y3,y4,y5,y6}], {{y1,y2,y3,y4,y5,y6}}].f1[{y1,y2,y3,y4,y5,y6}]];
B3=Simplify[D[B2,{{y1,y2,y3,y4,y5,y6}}].f1[{y1,y2,y3,y4,y5,y6}] ];
B4=D[B3,{{y1,y2,y3,y4,y5,y6}}].f1[{y1,y2,y3,y4,y5,y6}] ;
B5=D[B4,{{y1,y2,y3,y4,y5,y6}}].f1[{y1,y2,y3,y4,y5,y6}] ;

?f1

ButchD[f,n_]:=ButchD[{y1_,y2_,y3_,y4_,y5_,y6_},n]

Testing on our gravitational problem

Testing this code on the gravitaional problem

f2[{t_,y1_,y2_}] :={1,
-16 y1 + 12 y2 + 16 Cos[t] - 13 Sin[t],
12 y1 - 9 y2 -11 Cos[t] + 9 Sin[t]
}
IC = {0,1,0};

n=100;
h0=1.0;
LocalTol=(10^(-5));
ErrorOrder = 2;
{tVs, yVs}=DriverVariableStep1[IC,h0,f2,n,ErrorOrder,LocalTol,EulerStep];
ListPlot[ yVs[[All,{1,2}]],
PlotStyle->{Red, Green},
PlotRange->All,
AxesLabel->{"y1","y2"}]
ListPlot[Drop[RotateLeft[tVs]-tVs,-1],
PlotRange->All]

Butcher Trees

<<NumericalDifferentialEquationAnalysis`
RungeKuttaOrderConditions[3,3]

ButcherPlot[
ButcherTrees[5]
]

ButcherPhi[ ButcherTrees[3],s]