General Explicit RK Scheme Code

As always the best code does as many things as I can reasonably do.  Here code for a general RK stepper for an autonomous f with an error estimate.  The scheme is described using the Butcher tableau A and B where B includes the output prediction weights in a first vector, the error estimate weights in a second vector and an integer m giving the order of the LTE i.e.  EstimateLTE ~ C h^m

ExplicitRKStep[y_,h_,f_,ButcherTableau_]:= Module[
{A,B,m,name,s,K,k,yNew,EstimateLTE},
{A,B,m,name}=ButcherTableau;
s=Length[A];
K=ConstantArray[0,s];
If[A≠LowerTriangularize[A,-1], Return[Print["A does not describe an Explicit RK scheme!"]]];
K[[1]] = f[y];
Do[
K[[i]]=f[y+h*Sum[A[[i,j]] K[[j]],{j,1,i-1}]],
{i,2,s}
];
yNew=y+h*Sum[B[[1,i]] K[[i]],{i,1,s}];
EstimateLTE =Norm[h*Sum[B[[2,i]] K[[i]],{i,1,s}]];
{yNew, EstimateLTE}
]

Here are (I think) all the exlicit schemes with error estimates we have examined.

BS23={(({{0, 0, 0, 0}, {(1/2), 0, 0, 0}, {0, (3/4), 0, 0}, {(2/9), (1/3), (4/9), 0}})),(({{(2/9), (1/3), (4/9), 0}, {-(5/72), (1/12), (1/9), -(1/8)}})),3,"BS23"};
Euler={(({{0, 0}, {1, 0}})),(({{1, 0}, {1, -1}})),2,"Euler"};
Heun={(({{0, 0}, {1, 0}})),(({{(1/2), (1/2)}, {1, -1}})),2, "Heun"};
CK45={(({{0, 0, 0, 0, 0, 0}, {(1/5), 0, 0, 0, 0, 0}, {(3/40), (9/40), 0, 0, 0, 0}, {(3/10), -(9/10), (6/5), 0, 0, 0}, {-(11/54), (5/2), -(70/27), (35/27), 0, 0}, {(1631/55296), (175/512), (575/13824), (44275/110592), (253/4096), 0}})),(({{(37/378), 0, (250/621), (125/594), 0, (512/1771)}, {(277/64512), 0, -(6925/370944), (6925/202752), (277/14336), -(277/7084)}})),5, "CashKarp45"};
RKF45={(({{0, 0, 0, 0, 0, 0}, {(1/4), 0, 0, 0, 0, 0}, {(3/32), (9/32), 0, 0, 0, 0}, {(1932/2197), -(7200/2197), (7296/2197), 0, 0, 0}, {(439/216), -8, (3680/513), -(845/4104), 0, 0}, {-(8/27), 2, -(3544/2565), (1859/4104), -(11/40), 0}})),(({{(25/216), 0, (1408/2565), (2197/4104), -(1/5), 0}, {(1/360), 0, -(128/4275), -(2197/75240), (1/50), (2/55)}})),5, "RKF45"};
DP45={(({{0, 0, 0, 0, 0, 0, 0}, {(1/5), 0, 0, 0, 0, 0, 0}, {(3/40), (9/40), 0, 0, 0, 0, 0}, {(44/45), -(56/15), (32/9), 0, 0, 0, 0}, {(19372/6561), -(25360/2187), (64448/6561), -(212/729), 0, 0, 0}, {(9017/3168), -(355/33), (46732/5247), (49/176), -(5103/18656), 0, 0}, {(35/384), 0, (500/1113), (125/192), -(2187/6784), (11/84), 0}})),(({{(35/384), 0, (500/1113), (125/192), -(2187/6784), (11/84), 0}, {(71/57600), 0, -(71/16695), (71/1920), -(17253/339200), (22/525), -(1/40)}})),5,"DOPRI45"};

Here is a test that the error orders are as predicted.

AsymptoticRKErrorTest[y_,f_, Scheme_,{hMin_,hMax_}] := Module[
{n,nMax=12,Data,h, yOut, Err,LTEPower=Scheme[[3]],g,logh},
Data=Table[
h=hMin*((((hMax/hMin)))^(n/nMax));
{yOut, Err}=ExplicitRKStep[y,h,f,Scheme];
{h,Err},
{n,0,nMax}
];
g[logh_]=Fit[Log[Data],{1, logh},logh];
ListLogLogPlot[Data,
PlotLabel->StringForm["``: Claimed error (Ch^(``));   Err = ``",Scheme[[4]],LTEPower,(E^(g[Log["h"]]))],
FrameLabel->{"h","err"},
Frame->True
]
]

Clear[f1,y1,y2]
f1[{y1_,y2_,y3_,y4_,y5_,y6_}] :=Module[
{μ=1/81.75},
{y4,y5,y6,
2y5 + y1 - ((μ(y1+μ-1))/((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2))) -(((1-μ)(y1+μ))/((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2))),
-2y4 + y2 - ((μ y2)/((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2))) -(((1-μ)y2)/((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2))),
- ((μ y3)/((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2))) -(((1-μ)y3)/((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2)))
}
]

yIn=RandomReal[{-1,1},6];
AsymptoticRKErrorTest[yIn,f1,RKF45,{(10^(-3)),1.9}]
AsymptoticRKErrorTest[yIn,f1,DP45,{(10^(-3)),1.9}]
AsymptoticRKErrorTest[yIn,f1,CK45,{(10^(-3)),1.9}]

General Implicit RK Scheme Code

I want to set up a general Implicit RK scheme.  It should use a butcher tableau just like the explicit solver and I would like to implement as tests at least the three schemes 237a,b,c on p99.  They have Tableaus (with no error estimates so I put  ∞ in the locations )

RK237A={(({{(1/3), 0}, {1, 0}})),(({{(3/4), (1/4)}, {∞, ∞}})),∞,"Imp_RK237A"};
RK237B={(({{((5-3(Sqrt[2]))/4), ((7-5(Sqrt[2]))/4)}, {((1+(Sqrt[2]))/4), ((3-(Sqrt[2]))/4)}})),(({{((1+(Sqrt[2]))/4), ((3-(Sqrt[2]))/4)}, {∞, ∞}})),∞,"ImplicitRK237B"};
RK237C={(({{(1/4), (1/4)-((Sqrt[3])/6)}, {(1/4)+((Sqrt[3])/6), (1/4)}})),(({{(1/2), (1/2)}, {∞, ∞}})),∞,"ImplicitRK237C"};

The first thing to do is work out how to set up the implicit system that needs to be solved.  
For "A" I need to solve
	({{(Subscript[k, 1]), =, f[y+(h/3)(Subscript[k, 1])]}, {(Subscript[k, 2]), =, f[y+h (Subscript[k, 1])]}})
For "B" I need to solve
	({{(Subscript[k, 1]), =, f[y+h ( ((5-3(Sqrt[2]))/4)(Subscript[k, 1])+((7-5(Sqrt[2]))/4) (Subscript[k, 2]))]}, {(Subscript[k, 2]), =, f[y+h ( ((1+(Sqrt[2]))/4)(Subscript[k, 1])+((3-(Sqrt[2]))/4) (Subscript[k, 2]))]}})
For "C" I need to solve
	({{(Subscript[k, 1]), =, f[y+h ( (1/4)(Subscript[k, 1])+((1/4)-((Sqrt[3])/6)) (Subscript[k, 2]))]}, {(Subscript[k, 2]), =, f[y+h (( (1/4)+((Sqrt[3])/6))(Subscript[k, 1])+(1/4) (Subscript[k, 2]))]}})		
"A" is simpler because it only has one implicit equation. "B" and  "C" have two implicit equations.   My iterative fixed point implicit solver needs the A matrix, the step size h, the point y, the function f,  and an initial guess for  ALL the k values.  The solver is going to return the converged values of ALL the k values.   Here is a preliminary version.

ImplicitFixedPoint1[y_,h_,K0_,f_,A_]:= Module[
{s=Length[K0],KNew=K0,i=0, MaxIter=12,yCopies},
yCopies=Table[y,{s}];
Table[
KNew=Map[f,yCopies+h A.KNew], 
{MaxIter}
]
]

Clear[f,k1,k2,h,y,A]
n=4;
f[{y1_,y2_,y3_,y4_}]:= {y1 + (y2^2) Sin[y2 y1],  y4 + Cos[y1 y2],y3 +y2 y1, y1 - y4};
y0=RandomReal[{-1,1},n];
K0=RandomReal[ {-1,1},{2,n}];
A= RK237C[[1]];
h=0.1;
(* For 237A K0=Map[f,{y0+ (1/3)h f[y0],y0+h f[y0]}]; *)
K0=RandomReal[ {-1,1},{2,n}]
Data=ImplicitFixedPoint1[y0,h,K0,f,A];
MatrixForm[Data]

This is seriously inefficient.  What should I do to improve the fixed point iteration?

ImplicitFixedPoint2[y_,h_,K0_,f_,A_]:= Module[
{s=Length[K0],yCopies,KOld=K0,KNew,Tol=1.0 (10^(-7)),i=0, MaxIter=200},
yCopies=Table[y,{s}];
While[
And[
++i<MaxIter,
KNew=Map[f,yCopies+h A.KOld];Norm[KNew-KOld]>Tol],
KOld=KNew
];
Switch[i,
MaxIter,Print["Did not converge with tolerance ", Tol, "!"],
_,Print["Converged after ", i, " iterations!"]
];
KNew
]

Clear[f,k1,k2,h,y,A]
n=4;
f[{y1_,y2_,y3_,y4_}]:= {y1 + (y2^2) Sin[y2 y1],  y4 + Cos[y1 y2],y3 +y2 y1, y1 - y4};
y0=RandomReal[{-1,1},n];
K0=RandomReal[ {-1,1},{2,n}];
A= RK237C[[1]];
h=0.1;
(* For 237A K0=Map[f,{y0+ (1/3)h f[y0],y0+h f[y0]}]; *)
K0=RandomReal[ {-1,1},{2,n}]
Data=ImplicitFixedPoint2[y0,h,K0,f,A];

What would a Newton version look like? The first thing to do is work out what everything is. The equations are 
	({{(Subscript[k, 1]), =, f[y+h (Subscript[a, 1,1])(Subscript[k, 1])+h (Subscript[a, 1,2]) (Subscript[k, 2])]}, {(Subscript[k, 2]), =, f[y+h (Subscript[a, 2,1])(Subscript[k, 1])+h (Subscript[a, 2,2]) (Subscript[k, 2])]}})
and I need to remember that (Subscript[k, i]) are vectors in our vector space (ℝ^n).  Rewriting the equations gives 
	({{(Subscript[g, 1])[(Subscript[k, 1]),(Subscript[k, 2])], :=, (Subscript[k, 1]), -, f[y+h (Subscript[a, 1,1])(Subscript[k, 1])+h (Subscript[a, 1,2]) (Subscript[k, 2])], =, 0}, {(Subscript[g, 2])[(Subscript[k, 1]),(Subscript[k, 2])], :=, (Subscript[k, 2]), -, f[y+h (Subscript[a, 2,1])(Subscript[k, 1])+h (Subscript[a, 2,2]) (Subscript[k, 2])], =, 0}})
and taking the gradient w.r.t (Subscript[k, 1])and (Subscript[k, 2]) gives
	(Subscript[({{(Subscript[∇, 1])(Subscript[g, 1]), =, Id, -, h (Subscript[a, 1,1])∇(Subscript[f, 1])}, {(Subscript[∇, 2])(Subscript[g, 1]), =,  , -, h (Subscript[a, 1,2])∇(Subscript[f, 1])}, {(Subscript[∇, 1])(Subscript[g, 2]), =,  , -, h (Subscript[a, 2,1])∇(Subscript[f, 2])}, {(Subscript[∇, 2])(Subscript[g, 2]), =, Id, -, h (Subscript[a, 2,2])∇f}}), 2])
The structure in the expressions above is worth examining.
*  If (Subscript[a, 11])=(Subscript[a, 22]) then the complicted blocks are similar! Schemes with this characteristic are called Singly Diagonal Implicit
*  As h gets small the matrices become easy to "invert" because the diagonal terms look like the identity and the off-diagonal terms vanish.
* The diagonal matrices are invertible provided h<||(Subscript[a, 1,1])∇(Subscript[f, 1])(Subscript[||, 2]) and h<||(Subscript[a, 2,2])∇(Subscript[f, 2])(Subscript[||, 2]) 
In the fixed point iteration we organized the two vectors into a matrix K={(Subscript[k, 1]),(Subscript[k, 2])}∈(ℝ^(2×n)) .  In the Newton iteration we organize the two vectors into a single long vector 𝓀=(Subscript[k, 1]):(Subscript[k, 2])∈(ℝ^(2n)) and the two vector functions ℊ={(Subscript[g, 1]):(Subscript[g, 2])} into a single long vector because we need to solve linear systems. with this setup our overall gradient becomes
	∇ℊ=(({{Id-h (Subscript[a, 11])∇(Subscript[f, 1]), -h (Subscript[a, 12])∇(Subscript[f, 1])}, {-h (Subscript[a, 21])∇(Subscript[f, 2]), Id-h (Subscript[a, 22])∇(Subscript[f, 2])}}))
and our Newton iteration is
	(Subscript[𝓀, n+1])=(Subscript[𝓀, n])-LinearSolve[∇(Subscript[ℊ, n]),(Subscript[ℊ, n])]
followed by evaluations of ℊ and ∇ℊ.  The complications here are BOOKKEEPING!

ImplicitNewton1[y_,h_,K0_,{f_,df_},A_]:= Module[
{s=Length[K0],n=Length[y],Id,yCopies,KNew=Flatten[K0],MaxIter=12,
g,fder,gDer},
yCopies=Table[y,{s}];
Id=IdentityMatrix[s];
g[k_]:=k- Flatten[Map[f,yCopies+h A.Partition[k,n]]];
fder[k_]:= Map[df, yCopies+h A.Partition[k,n]];
gDer[k_]:=Module[
{dfVals=fder[k]},
IdentityMatrix[s*n]-h ArrayFlatten[
Table[A[[i,j]] dfVals[[j]],{i,s},{j,s}]]
];
Table[
KNew=KNew-LinearSolve[gDer[KNew],g[KNew]],
{MaxIter}
]
]

Clear[f,k1,k2,h,y,A]
n=4;
f[{y1_,y2_,y3_,y4_}]:= {y1 + (y2^2) Sin[y2 y1],  y4 + Cos[y1 y2],y3 +y2 y1, y1 - y4};
df[{y1_,y2_,y3_,y4_}]=D[f[{y1,y2,y3,y4}],{{y1,y2,y3,y4}}];
y0=RandomReal[{-1,1},n];
K0=RandomReal[ {-1,1},{2,n}];
A= RK237C[[1]];
h=0.8;
(* For 237A K0=Map[f,{y0+ (1/3)h f[y0],y0+h f[y0]}]; *)
Fred=ImplicitNewton1[y0,h,K0,{f,df},A]

RK Order Equations

Mathematica has the capability of working out the Butcher conditions and solving for low stae order for us.  Here is an explicit computation and all the 3 order 3 stage RK schemes.

Needs["NumericalDifferentialEquationAnalysis`"]
{order,s}={3,3};
Skeletons={LowerTriangularize[Array[a,{s,s}],-1],Array[b,{s}]}
Eqs=Flatten[
RungeKuttaOrderConditions[order,s,
RungeKuttaMethod->Explicit,
ButcherRowSum->True]]
Sols=Simplify[Skeletons/.Solve[Eqs]];
TableForm[Sols]

The error terms are combinations of derivatives multiplied by the following expressions.

ButcherPrincipalError[order,s,
RungeKuttaMethod->Explicit]

The fancy looking schemes minimize these coefficients and almost always need to round to fractions.

The things these coefficeints multiply are computable.