Low-Storage Runge-Kutta Schemes
Chad- Elaheh- Jane

Introduction

Storage requirements are typically one of the main considerations, especially if a high order ODE solver is required. Many problem of mathematical physics result in the numerical approximation of systems of coupled ordinary equations. The technique for integrating the equations depends on, among other things, the desired accuracy, efficiency, robustness and simplicity. For the ODE's that result from the direct simulation of fluid dynamics, the overriding consideration is often the availability of fast memory. A numerical integration technique that minimizes memory storage is essential and can be formulated with a Runge-Kutta (RK) methodology.  
The purpose of this project is to describe a set of low storage Runge-Kutta schemes namely 2N-storage. 
A family of three-stage third-order Runge-Kutta schemes is implemented, these schemes requires only two storage locations.
These schemes are competitive with high-storage Runge-Kutta schemes and are considerably more efficient and accurate.

Williamson (2N) methods

Williamson first presented the relationship between the general RK scheme and the 2N-storage scheme. Specifically the following equation provide the relationship between the 2N-storage variables Aj , Bj and the conventional RK variables (Subscript[a, i,j]) , (Subscript[b, j]) and (Subscript[c, i])



And according to these equalities, The original Butcher array is expanded in terms of the 2N-storage variables Aj and Bj to obtain:

Condition for casting a scheme into 2N-storage scheme

The low-storage algorithm can be used if we have:

(1- c2- b3)  a3,2 =   (c3 -c2)  b2 
 
Or based on c2 and c3: 

(c3^2)(1-c2)+c3((c2^2)+(1/2)c2-1)+((1/3)-(1/2)c2)=0

Butcher principal error for three-stage third-order Runge Kutta schemes

<<NumericalDifferentialEquationAnalysis`
RungeKuttaOrderConditions[3,3,
ButcherRowSum->True,
RungeKuttaMethod->Explicit
]

Eqs={(Subscript[b, 1])+(Subscript[b, 2])+(Subscript[b, 3])==1,(Subscript[b, 2]) (Subscript[c, 2])+(Subscript[b, 3]) (Subscript[c, 3])==(1/2),(Subscript[a, 3,2]) (Subscript[b, 3]) (Subscript[c, 2])==(1/6),(Subscript[b, 2]) (Subscript[c, 2]^2)+(Subscript[b, 3]) (Subscript[c, 3]^2)==(1/3),(Subscript[a, 2,1])==(Subscript[c, 2]),(Subscript[a, 3,1])+(Subscript[a, 3,2])==(Subscript[c, 3])}

sols= Solve[Eqs]

pErr=ButcherPrincipalError[3,3,
ButcherRowSum->True,
RungeKuttaMethod->Explicit]

pErr=ButcherPrincipalError[3,3,
ButcherRowSum->True,
RungeKuttaMethod->Explicit]/.sols[[1]]
Goodness[a32_,b2_,b3_,c2_,c3_]=Simplify[pErr.pErr/.{
(Subscript[a, 3,2])->a32,
(Subscript[b, 2])->b2,
(Subscript[b, 3])->b3,
(Subscript[c, 2])->c2,
(Subscript[c, 3])->c3}]

pErr=ButcherPrincipalError[3,3,
ButcherRowSum->True,
RungeKuttaMethod->Explicit]/.sols[[2]]
Goodness[a32_,b2_,b3_,c2_,c3_]=Simplify[pErr.pErr/.{
(Subscript[a, 3,2])->a32,
(Subscript[b, 2])->b2,
(Subscript[b, 3])->b3,
(Subscript[c, 2])->c2,
(Subscript[c, 3])->c3}]

pErr=ButcherPrincipalError[3,3,
ButcherRowSum->True,
RungeKuttaMethod->Explicit]/.sols[[3]]
Goodness[a32_,b2_,b3_,c2_,c3_]=Simplify[pErr.pErr/.{
(Subscript[a, 3,2])->a32,
(Subscript[b, 2])->b2,
(Subscript[b, 3])->b3,
(Subscript[c, 2])->c2,
(Subscript[c, 3])->c3}]

ContourPlot[(1/(7776 b3))(80-60 (Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))])+3 b3 (54-90 c2+83 (c2^2)-24 (c2^3)+9 (c2^4))+9 (c2^3) (-3+(Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))]))-6 (c2^2) (-9+4 (Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))]))+3 c2 (-48+19 (Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))]))),
{b3,0,1},
{c2,0,1},
Contours->40,
ColorFunction->Hue,
FrameLabel->{"b3","c2"}
]

FindMinimum[{
(1/(7776 b3))(80-60 (Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))])+3 b3 (54-90 c2+83 (c2^2)-24 (c2^3)+9 (c2^4))+9 (c2^3) (-3+(Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))]))-6 (c2^2) (-9+4 (Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))]))+3 c2 (-48+19 (Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))])))},
{b3,0.4},{c2,0.4}]

pErr=ButcherPrincipalError[3,3,
ButcherRowSum->True,
RungeKuttaMethod->Explicit]/.sols[[4]]
Goodness[a32_,b2_,b3_,c2_,c3_]=Simplify[pErr.pErr/.{
(Subscript[a, 3,2])->a32,
(Subscript[b, 2])->b2,
(Subscript[b, 3])->b3,
(Subscript[c, 2])->c2,
(Subscript[c, 3])->c3}]

ContourPlot[(1/(7776 b3))(80+60 (Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))])+3 b3 (54-90 c2+83 (c2^2)-24 (c2^3)+9 (c2^4))-9 (c2^3) (3+(Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))]))+6 (c2^2) (9+4 (Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))]))-3 c2 (48+19 (Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))]))),
{b3,0,1},
{c2,0,1},
Contours->40,
ColorFunction->Hue,
FrameLabel->{"b3","c2"}
]

FindMinimum[{
(1/(7776 b3))(80+60 (Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))])+3 b3 (54-90 c2+83 (c2^2)-24 (c2^3)+9 (c2^4))-9 (c2^3) (3+(Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))]))+6 (c2^2) (9+4 (Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))]))-3 c2 (48+19 (Sqrt[3]) (Sqrt[b3 (4-6 c2+3 b3 (c2^2))])))},
{b3,0.4},{c2,0.4}]

Runge Kutta Order Conditions

For third-order accuracy, seven constrains between the coefficients ai,j , bj and ci must be satisfied:

<<NumericalDifferentialEquationAnalysis`
RungeKuttaOrderConditions[3,3,
ButcherRowSum->True,
RungeKuttaMethod->Explicit
]

Eqs={{(Subscript[b, 1])+(Subscript[b, 2])+(Subscript[b, 3])==1},{(Subscript[b, 2]) (Subscript[c, 2])+(Subscript[b, 3]) (Subscript[c, 3])==(1/2)},{(Subscript[a, 3,2]) (Subscript[b, 3]) (Subscript[c, 2])==(1/6),(Subscript[b, 2]) (Subscript[c, 2]^2)+(Subscript[b, 3]) (Subscript[c, 3]^2)==(1/3)},{(Subscript[a, 2,1])==(Subscript[c, 2]),(Subscript[a, 3,1])+(Subscript[a, 3,2])==(Subscript[c, 3])}};
sol=Solve[Flatten[{Eqs}],{(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[a, 3,1]),(Subscript[a, 3,2]),(Subscript[a, 2,1])}]

Clear[c2,c3]
Solve[(c3^2)(1-c2)+c3((c2^2)+(1/2)c2-1)+((1/3)-(1/2)c2)==0,{c2,c3}]

LTE=ButcherPrincipalError[3,3,
ButcherRowSum->True,
RungeKuttaMethod->Explicit
]/.sol

BPE=ButcherPrincipalError[3,3,ButcherRowSum->True,RungeKuttaMethod->Explicit]/.sol;
BPE2=Flatten[BPE];
BPE3=Sqrt[BPE2.BPE2]/.{(Subscript[c, 2])->c2,(Subscript[c, 3])->c3}
FindMinimum[{BPE3,(c3^2)(1-c2)+c3((c2^2)+(1/2)c2-1)+((1/3)-(1/2)c2)==0},{{c2,0.5},{c3,0.75}}]
Show[
ContourPlot[BPE3,{c2,0,1},{c3,0,1},Contours->50,ColorFunction->Hue],
ContourPlot[(c3^2)(1-c2)+c3((c2^2)+(1/2)c2-1)+((1/3)-(1/2)c2)==0,{c2,0,1},{c3,0,1},ContourStyle->Black]
]

The curve denotes values of the coefficients c2 and c3 for third-order low-storage schemes are possible.

Comparison of stability

RK33[y_,f_,h_]:= Module[
{k1,k2,k3,k4,k5,k6,yNew,zNew},
k1=f[y];
k2=f[y+(1/3) h k1];
k3=f[y-(3/16) h k1+(15/16) h k2];
yNew=y + h ((1/6) k1 + (3/10)k2 + (8/15) k3);
zNew=y+h k1;
yNew
];

Williamson[y_,f_,h_]:= Module[
{k1,k2,k3,k4,k5,k6,yNew,zNew},
k1=f[y];
k2=f[y+0.45736 h k1];
k3=f[y-0.13268h k1+0.92530 h k2];
yNew=y + h (0.19546 k1 + 0.41071k2 + 0.39383 k3);
zNew=y+h k1;
yNew
];

RKClassic[y_,f_,h_]:= Module[
{k1,k2,k3,k4,k5,k6,yNew,zNew},
k1=f[y];
k2=f[y+(1/2) h k1];
k3=f[y-1 h k1+2 h k2];
yNew=y + h ((1/6) k1 + (2/3)k2 + (1/6) k3);
zNew=y+h k1;
yNew
];

CEJ[y_,f_,h_]:= Module[
{k1,k2,k3,k4,k5,k6,yNew,zNew},
k1=f[y];
k2=f[y+0.539 h k1];
k3=f[y-0.115 h k1+0.890 h k2];
yNew=y + h (0.225 k1 + 0.428k2 + 0.347 k3);
zNew=y+h k1;
yNew
];

RK43[y_,f_,h_]:= Module[
{k1,k2,k3,k4,k5,k6,yNew},
k1=f[y];
k2=f[y+(1/3) h k1];
k3=f[y-(5/12) h k1+(3/4) h k2];
k4=f[y+(1/4) h k1+(1/12) h k2+(2/3)h k3];
yNew=y + h ((1/3) k2 + (5/12)k3 + (1/4) k4);
yNew
];

Williamson2N[y_,f_,h_]:= Module[
{yNew,s1,s2,A,B,i=1},
A={0,-0.63768,-1.30663};
B={0.45736,0.92530,0.39383};
s1=y;
While[
i≤3,
s2=A[[i]]s2+h f[s1];
s1=s1+B[[i]]s2;
i++
];
yNew=s1;
yNew
];

CEJ2N[y_,f_,h_]:= Module[
{yNew,s1,s2,A,B,i=1},
A={0,-0.734,-1.331};
B={0.539,0.890,0.347};
s1=y;
While[
i≤3,
s2=A[[i]]s2+h f[s1];
s1=s1+B[[i]]s2;
i++
];
yNew=s1;
yNew
];

RK332N[y_,f_,h_]:= Module[
{yNew,s1,s2,A,B,i=1},
A={0,((-5)/9),((-153)/128)};
B={(1/3),(15/16),(8/15)};
s1=y;
While[
i≤3,
s2=A[[i]]s2+h f[s1];
s1=s1+B[[i]]s2;
i++
];
yNew=s1;
yNew
];

RK43[y_,f_,h_]:= Module[
{k1,k2,k3,k4,k5,k6,yNew,zNew},
k1=f[y];
k2=f[y+(1/3) h k1];
k3=f[y-(5/12) h k1+(3/4) h k2];
k4=f[y+(1/4) h k1+(1/12) h k2+(2/3)h k3];
yNew=y + h ((1/3) k2 + (5/12)k3 + (1/4) k4);
zNew=y+h k1;
yNew
]  ;

y0=1;
M=3;
g[z_]:=λ z;
Stab1[z_]=Simplify[RK33[1,g,(z/λ)]]

y0=1;
M=3;
g[z_]:=λ z;
Stab2[z_]=Simplify[CEJ[1,g,(z/λ)]]

y0=1;
M=3;
g[z_]:=λ z;
Stab3[z_]=Simplify[Williamson[1,g,(z/λ)]]

y0=1;
M=3;
g[z_]:=λ z;
Stab4[z_]=Simplify[RKClassic[1,g,(z/λ)]]

y0=1;
M=3;
g[z_]:=λ z;
Stab5[z_]=Simplify[CEJ2N[1,g,(z/λ)]]

y0=1;
M=3;
g[z_]:=λ z;
Stab6[z_]=Simplify[Williamson2N[1,g,(z/λ)]]

y0=1;
M=3;
g[z_]:=λ z;
Stab7[z_]=Simplify[RK332N[1,g,(z/λ)]]

y0=1;
M=3;
g[z_]:=λ z;
Stab8[z_]=Simplify[RK43[1,g,(z/λ)]]

Show[
ContourPlot[
Abs[Stab1[x+I y]]==1,
{x,-M,M},{y,-M,M},ContourStyle->Dashed]
,
ContourPlot[
Abs[Stab2[x+I y]]==1,
{x,-M,M},{y,-M,M},ContourStyle->Dashed]
,
ContourPlot[
Abs[Stab3[x+I y]]==1,
{x,-M,M},{y,-M,M},ContourStyle->Dashed]
,
ContourPlot[
Abs[Stab4[x+I y]]==1,
{x,-M,M},{y,-M,M},ContourStyle->Dashed]
,
ContourPlot[
Abs[Stab5[x+I y]]==1,
{x,-M,M},{y,-M,M},ContourStyle->Dashed]
,
ContourPlot[
Abs[Stab6[x+I y]]==1,
{x,-M,M},{y,-M,M},ContourStyle->Dashed]
,
ContourPlot[
Abs[Stab7[x+I y]]==1,
{x,-M,M},{y,-M,M},ContourStyle->Dashed]
,
ContourPlot[
Abs[Stab8[x+I y]]==1,
{x,-M,M},{y,-M,M}]
]

As we can see, four-stage third-order runge kutta scheme is more stable than three-stage third-order Runge Kutta Scheme.

Comparison of all the schemes in terms of Timing and Local Trunction Error

RK33[y_,f_,h_]:= Module[
{k1,k2,k3,k4,k5,k6,yNew},
k1=f[y];
k2=f[y+(1/3) h k1];
k3=f[y-(3/16) h k1+(15/16) h k2];
yNew=y + h ((1/6) k1 + (3/10)k2 + (8/15) k3);
{yNew}
];

Williamson[y_,f_,h_]:= Module[
{k1,k2,k3,k4,k5,k6,yNew},
k1=f[y];
k2=f[y+0.45736 h k1];
k3=f[y-0.13268h k1+0.92530 h k2];
yNew=y + h (0.19546 k1 + 0.41071k2 + 0.39383 k3);
{yNew}
];

Williamson2N[y_,f_,h_]:= Module[
{yNew,s1,s2,A,B,i=1},
A={0,-0.63768,-1.30663};
B={0.45736,0.92530,0.39383};
s1=y;
While[
i≤3,
s2=A[[i]]s2+h f[s1];
s1=s1+B[[i]]s2;
i++
];
yNew=s1;
{yNew}
];

CEJ2N[y_,f_,h_]:= Module[
{yNew,s1,s2,A,B,i=1},
A={0,-0.73421,-1.33302};
B={0.53918,0.89055,0.34710};
s1=y;
While[
i≤3,
s2=A[[i]]s2+h f[s1];
s1=s1+B[[i]]s2;
i++
];
yNew=s1;
{yNew}
];

RKClassic[y_,f_,h_]:= Module[
{k1,k2,k3,k4,k5,k6,yNew},
k1=f[y];
k2=f[y+(1/2) h k1];
k3=f[y-1 h k1+2 h k2];
yNew=y + h ((1/6) k1 + (2/3)k2 + (1/6) k3);
{yNew}
];

CEJ[y_,f_,h_]:= Module[
{k1,k2,k3,k4,k5,k6,yNew},
k1=f[y];
k2=f[y+0.53918 h k1];
k3=f[y-0.11468 h k1+0.89055 h k2];
yNew=y + h (0.22504 k1 + 0.42786k2 + 0.34710 k3);
{yNew}
];

RK43[y_,f_,h_]:= Module[
{k1,k2,k3,k4,k5,k6,yNew},
k1=f[y];
k2=f[y+(1/3) h k1];
k3=f[y-(5/12) h k1+(3/4) h k2];
k4=f[y+(1/4) h k1+(1/12) h k2+(2/3)h k3];
yNew=y + h ((1/3) k2 + (5/12)k3 + (1/4) k4);
{yNew}
];

RK332N[y_,f_,h_]:= Module[
{yNew,s1,s2,A,B,i=1},
A={0,((-5)/9),((-153)/128)};
B={(1/3),(15/16),(8/15)};
s1=y;
While[
i≤3,
s2=A[[i]]s2+h f[s1];
s1=s1+B[[i]]s2;
i++
];
yNew=s1;
{yNew}
];

RK432N[y_,f_,h_]:= Module[
{yNew,s1,s2,A,B,i=1},
A={0,-1,-1,-1};
B={(1/3),(3/4),(2/3),(1/4)};
s1=y;
While[
i≤4,
s2=A[[i]]s2+h f[s1];
s1=s1+B[[i]]s2;
i++
];
yNew=s1;
{yNew}
];

RKF45[y_,f_,h_]:= Module[
{k1,k2,k3,k4,k5,k6,yNew},
k1=f[y];
k2=f[y+ 1/4 h k1];
k3=f[y+3/32 h k1+9/32 h k2];
k4=f[y+1932/2197 h k1−7200/2197 h k2+7296/2197 h k3];
k5=f[y+439/216 h k1-8 h k2+3680/513 h k3−845/4104 h k4];
k6=f[y-8/27 h k1+2h k2−3544/2565 h k3+1859/4104 h k4−11/40 h k5];
yNew=y + h (25/216 k1 + 1408/2565 k3 + 2197/4104 k4−1/5 k5);
{yNew}
];

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)
}
];
IC1 = {0.994,0,0,0,-2.0015851063790825224,0};
IC2 = {0.87978,0,0,0,-0.3797,0};
TMax1=17.06521656015796;
TMax2=19.14045706162071;

FixedStep[y0_,h0_,f_,TMax_,Stepper_,NumFn_]:=Module[
{yy=y0,ConvData},
ConvData=Table[
yy=Stepper[yy,f,h0][[1]];
{yy[[1]],yy[[2]]},
{TMax/h0}];
Print["Total Steps: ",TMax/h0];
Print["Function Evaluations: ",TMax/h0 NumFn];
ConvData
];

h0=0.0002;
{Time1,Data1}=Timing[FixedStep[IC1,h0,f1,TMax1,RK33,3]];
Print["Time of computation is ",Time1," seconds."]
ListPlot[Data1,
PlotRange->All,
AxesLabel->{"y1","y2"}]

h0=0.0002;
{Time1,Data1}=Timing[FixedStep[IC1,h0,f1,TMax1,RKClassic,3]];
Print["Time of computation is ",Time1," seconds."]
ListPlot[Data1,
PlotRange->All,
AxesLabel->{"y1","y2"}]

h0=0.0002;
{Time1,Data1}=Timing[FixedStep[IC1,h0,f1,TMax1,CEJ,3]];
Print["Time of computation is ",Time1," seconds."]
ListPlot[Data1,
PlotRange->All,
AxesLabel->{"y1","y2"}]

h0=0.0016;
{Time1,Data1}=Timing[FixedStep[IC1,h0,f1,TMax1,Williamson2N,3]];
Print["Time of computation is ",Time1," seconds."]
ListPlot[Data1,
PlotRange->All,
AxesLabel->{"y1","y2"}]

h0=0.0002;
{Time1,Data1}=Timing[FixedStep[IC1,h0,f1,TMax1,RK332N,3]];
Print["Time of computation is ",Time1," seconds."]
ListPlot[Data1,
PlotRange->All,
AxesLabel->{"y1","y2"}]

h0=0.0002;
{Time1,Data1}=Timing[FixedStep[IC1,h0,f1,TMax1,CEJ2N,3]];
Print["Time of computation is ",Time1," seconds."]
ListPlot[Data1,
PlotRange->All,
AxesLabel->{"y1","y2"}]

h0=0.001;
{Time1,Data1}=Timing[FixedStep[IC1,h0,f1,TMax1,RK43,3]];
Print["Time of computation is ",Time1," seconds."]
ListPlot[Data1,
PlotRange->All,
AxesLabel->{"y1","y2"}]

h0=0.0002;
{Time1,Data1}=Timing[FixedStep[IC1,h0,f1,TMax1,Williamson2N,3]];
Print["Time of computation is ",Time1," seconds."]
ListPlot[Data1,
PlotRange->All,
AxesLabel->{"y1","y2"}]

h0=0.001;
{Time1,Data1}=Timing[FixedStep[IC1,h0,f1,TMax1,RK432N,3]];
Print["Time of computation is ",Time1," seconds."]
ListPlot[Data1,
PlotRange->All,
AxesLabel->{"y1","y2"}]

ySol=y/.NDSolve[ {
y'[t]==f1[y[t]],
y[0]==IC1
},
y,
{t, 0, TMax1},Method->"ExplicitRungeKutta"][[1]];
MmaPic=ParametricPlot[ySol[t][[{1,2}]],{t,0,TMax1}]

LTE[IC_,f_,h0_,TMax_,Stepper_]:=Module[
{y,t=0,err,i=1,j},
y=IC;
While[
And[t≤TMax,i≤(TMax/h0)], 
t=t+h0;
err[i]=Norm[Stepper[y,f,h0][[1]][[{1,2}]]-ySol[t][[{1,2}]]];
y=Stepper[y,f,h0][[1]];
i++
];
Table[err[j],{j,1,i-1}]
];

h0=0.0016
Show[
ListLogPlot[LTE[IC1,f1,h0,TMax1,RK33](*,PlotRange->{{5000,7000},{0.0000001,0.001}}*)],
ListLogPlot[LTE[IC1,f1,h0,TMax1,RKClassic],PlotStyle->Green(*,PlotRange->{{5000,7000},{0.0000001,0.001}}*)],
ListLogPlot[LTE[IC1,f1,h0,TMax1,Williamson],PlotStyle->Pink(*,PlotRange->{{5000,7000},{0.0000001,0.001}}*)],
ListLogPlot[LTE[IC1,f1,h0,TMax1,CEJ],PlotStyle->Red(*,PlotRange->{{5000,7000},{0.0000001,0.001}}*)],
AxesLabel->{"Steps","Error"}
]

h0=0.0004;
Show[
ListLogPlot[LTE[IC1,f1,h0,TMax1,RK33](*,PlotRange->{{5000,7000},{0.0000001,0.001}}*)],
ListLogPlot[LTE[IC1,f1,h0,TMax1,RKClassic],PlotStyle->Green(*,PlotRange->{{5000,7000},{0.0000001,0.001}}*)],
ListLogPlot[LTE[IC1,f1,h0,TMax1,Williamson],PlotStyle->Pink(*,PlotRange->{{5000,7000},{0.0000001,0.001}}*)],
ListLogPlot[LTE[IC1,f1,h0,TMax1,CEJ],PlotStyle->Red(*,PlotRange->{{5000,7000},{0.0000001,0.001}}*)],
AxesLabel->{"Steps","Error"}
]

h0=0.0004;
Show[
ListLogPlot[LTE[IC1,f1,h0,TMax1,RK332N]],
ListLogPlot[LTE[IC1,f1,h0,TMax1,Williamson2N],PlotStyle->Pink],
ListLogPlot[LTE[IC1,f1,h0,TMax1,CEJ2N],PlotStyle->Red],
AxesLabel->{"Steps","Error"}
]

h0=0.0016;
Show[
ListLogPlot[LTE[IC1,f1,h0,TMax1,CEJ]],
ListLogPlot[LTE[IC1,f1,h0,TMax1,CEJ2N],PlotStyle->Red],
AxesLabel->{"Steps","Error"}
]

h0=0.0016;
Show[
ListLogPlot[LTE[IC1,f1,h0,TMax1,Williamson]],
ListLogPlot[LTE[IC1,f1,h0,TMax1,Williamson2N],PlotStyle->Red],
AxesLabel->{"Steps","Error"}
]

h0=0.0016;
Show[
ListLogPlot[LTE[IC1,f1,h0,TMax1,RK33]],
ListLogPlot[LTE[IC1,f1,h0,TMax1,RK332N],PlotStyle->Red],
AxesLabel->{"Steps","Error"}
]

Timing  and Errors comparison tables

Timing

TableForm[{{2.402,2.433,2.324,2.324},{4.805,4.836,4.633,4.634},{9.61,9.719,9.276,9.282},{19.157,19.422,18.439,18.549}},TableHeadings->{{"0.0016","0.0008","0.0004","0.0002"},{"RKClassic","RK33","Williamson","CEJ"}}]

Error

TableForm[{{0.04483296,0.0611457,0.0344986,0.0176908},{0.00985629,0.00969983,0.00620298,0.00226482},{0.00143659,0.00131171,0.000875916,0.000222103},{0.00019233,0.000169525,0.000115496,0.00030449}},TableHeadings->{{"0.0016","0.0008","0.0004","0.0002"},{"RKClassic","RK33","Williamson","CEJ"}}]

Timing

TableForm[{{2.433,2.356},{4.836,4.68},{9.719,9.407},{19.422,18.72}},TableHeadings->{{"0.0016","0.0008","0.0004","0.0002"},{"RK33","RK332N"}}]

Timing

TableForm[{{2.324,2.262},{4.633,4.509},{9.267,8.986},{18.439,17.959}},TableHeadings->{{"0.0016","0.0008","0.0004","0.0002"},{"Williamson","Williamson2N"}}]

Timing

TableForm[{{2.324,2.262},{4.634,4.54},{9.282,9.142},{18.549,19.018}},TableHeadings->{{"0.0016","0.0008","0.0004","0.0002"},{"CEJ33","CEJ2N2N"}}]

Error

TableForm[{{0.0611457,0.0611457},{0.00969983,0.00969983},{0.00131171,0.00131171},{0.000169525,0.000169525}},TableHeadings->{{"0.0016","0.0008","0.0004","0.0002"},{"RK33","RK332N"}}]

Error

TableForm[{{0.0344986,0.0345098},{0.00620293,0.00620664},{0.000875916,0.000878345},{0.000115496,0.000117896}},TableHeadings->{{"0.0016","0.0008","0.0004","0.0002"},{"Williamson","Williamson2N"}}]

Error

TableForm[{{0.0176908,0.0192416},{0.00220482,0.00368141},{0.000222103,0.000896391},{0.00030449,0.000645435}},TableHeadings->{{"0.0016","0.0008","0.0004","0.0002"},{"CEJ33","CEJ2N2N"}}]

Conclusion

Low-storage schemes are possible up to the fourth order, that nearly all schemes can be immunized against round-off error, and that their truncation error bounds are comparable to those of the more familiar Runge-Kutta methods.
2N-Storage schemes are better than conventional ones in terms of timing.
Our schemes CEJ and CEJ(2N) are the best ones among all (3,3) RK and  (3,3) RK 2N-Storage schemes in terms of LTE, but a little worse than Williamson and Williamson(2N) in terms of timing.
Stability region of our scheme CEJ is the same as the other (3,3) schemes.

References

[1] Low-Storage Runge-Kutta Schems, J.H. Williamson, JCP 35, 48-56 (1980)
[2] Fourth-Order 2N-Storage Runge-Kutta Schemes, Mark H. Carpenter, Christopher A. Kennedy, NASA Technical     Memorandum 109112
[3] Third-Order 2N-Storage Runge-Kutta Schemes with Error Control, Mark H. Carpenter, Christopher A. Kennedy, NASA Technical Memorandum 109111
[4] Runge-Kutta methods with minimum-storage implementations, David I. Ketcheson