Euler Compensated Summation Comparison

Implementing the code fragments from Butcher p82.

In native Mathematica

f[y_]:=y
ySol[t_]:=(E^t)
{pMax, base}={6,4};
EData=Table[
n=(base^p);
{t,y,h}={0.0,1.0,1.0/n};
Do[
t=t+h;
term = h*f[y];
y=y+term,
{n}];
{h, y},
{p, 1,pMax}]

CSEData=Table[
n=(base^p);
{t,y,h}={0.0,1.0,1.0/n};
z=0;
Do[
t=t+h;
term =h*f[y]+z;
newy=y+term;
z=term-(newy-y);
y=newy,
{n}];
{h, y},
{p, 1,pMax}]

The ouput from the two variants is slightly different.

FullForm[{CSEData[[-1,2]],EData[[-1,2]]}-E]

The question is which is better.

With native machine arithmetic.

Here they are with a compiled functions

Clear[f,y,ySol,t]
f= Compile[{{y, _Real}}, y]
ySol[t_]:=(E^t)
{pMax, base}={6,4};
EData=Table[
n=(base^p);
{t,y,h}={0.0,1.0,1.0/n};
Do[
t=t+h;
term = h*f[y];
y=y+term,
{n}];
{h, y},
{p, 1,pMax}]

CSEData=Table[
n=(base^p);
{t,y,h}={0.0,1.0,1.0/n};
z=0;
Do[
t=t+h;
term =h*f[y]+z;
newy=y+term;
z=term-(newy-y);
y=newy,
{n}];
{h, y},
{p, 1,pMax}]

Again the answers are slightly different

FullForm[{CSEData[[-1,2]],EData[[-1,2]]}-E]

Implementing low precision arithmetic

y=12.345678;
Digits=4;
{mantissa,exponent}=RealDigits[y]
a=N[FromDigits[{Take[mantissa,Digits],exponent}]

LowPrec[y_]:= Module[ {Digits=5,mantissa,exponent},
{mantissa,exponent}=RealDigits[y];
N[FromDigits[{Take[mantissa,Digits],exponent}]]
]
LowPrec[0.0]=0.0

Clear[f,y,ySol,t]
f[y_]:=LowPrec[y];
ySol[t_]:=(E^t)
{pMax, base}={14,2};
EData=Table[
n=(base^p);
{t,y,h}={0.0,1.0,1.0/n};
h = LowPrec[h];
Do[
t=t+h;
term =LowPrec[ h*f[y]];
y=LowPrec[LowPrec[y]+LowPrec[term]],
{n}];
{h, y},
{p, 1,pMax}]

CSEData=Table[
n=(base^p);
{t,y,h}={0.0,1.0,1.0/n};
h=LowPrec[h];
z=0.0;
Do[
t=t+h;
term =LowPrec[LowPrec[h*f[y]]+LowPrec[z]];
newy=LowPrec[LowPrec[y]+LowPrec[term]];
z=LowPrec[LowPrec[term]-LowPrec[LowPrec[newy]-LowPrec[y]]];
y=LowPrec[newy],
{n}];
{h, y},
{p, 1,pMax}]
ListLogLogPlot[{CSEData,EData},
PlotStyle->{Red, Blue},
PlotRange->All]

Stage Order

Each stage level in an RK scheme has its own order. For instance consider our simple Kutta 3rd order scheme contains two internal stage schemes.

ODEK3[y_,f_,h_]:= Module[
{k1,k2,k3,k4},
k1=f[y];                      
k2=f[y+ 1/2 h k1]; (* Euler Scheme (Subscript[Y, (1/2)])= (Subscript[y, 0])+ 1/2 k f[y] - 1st order scheme *)
k3=f[y+h (2 k2 -k1)]; (* (Subscript[Y, 1]) = (Subscript[y, 0]) + h (2 k2-k1) 2 stage explicit RK scheme *) 
y+(h/6)(k1+4k2 + k3)
]

The accuracy of the overall scheme relies on delicate cancellation (involving high derivatives of the solutions) between errors in the internal schemes.

Clear[f,y,h,x]
LTE=Series[y[x+h]-ODEK3[y[x],f, h],{h, 0,4}]/.{
y'[x]->f[y[x]],
y''[x]->f'[y[x]]f[y[x]],
y'''[x]->f''[y[x]]f[y[x]]f[y[x]]+f'[y[x]] f'[y[x]] f[y[x]]
}

This cancellation is not a problem for non-stiff problems and we will see the 3rd order behavior because the high derivatives are the same order of magnitude as the solution.  Howevere for stiff problems high derivatives are much larger than the solution and it is important that all stages of the computation are equally accurate.

The first stage of the above is basically an Euler half step with LTE

Clear[f,y,h,x]
Series[y[x+(1/2)h]-(y[x]+ 1/2 h f[y[x]]),{h,0, 3}]/.{
y'[x]->f[y[x]]
}

So the stage order of the first stage is only 1:  Not a surprise it is Euler after all.

The second stage LTE is

Series[
y[x+h]-(y[x]+h (2 f[y[x]+(h/2)f[y[x]]] -f[y[x]])),
{h, 0, 3}]/.{
y'[x]->f[y[x]],
y''[x]->f'[y[x]]f[y[x]]
}

which also only has order 1.

The stage order of an RK scheme is the lowest order of any of the internal stages.  In practice for a stiff problem you do not get the overall order but only the stage order on stiff problems.   So for stiff problems this scheme (which has Euler buried inside) is not very good.

Stability (reminder)

While we are at it lets look at the linear stability of this scheme.  We look at our linear test problem (reduced to a scalar using our eigenvalue trick)  when the real solution should decay and ask when the approximate solutions decays.

f[y_]:=λ y
Simplify[ODEK3[y0,f,h]] 
Stab[ξ_]=Simplify[ODEK3[1,f,ξ/λ]]  (* scale so that ξ=hλ and y0=1*)
M=3;
ContourPlot[
Abs[Stab[x+ I y]]==1,
{x, -M,0},
{y, 0,M}
]

The only difference with our fancier schemes is that we have more than one amplification factor.

Definition:  The scheme is A stable if the stability region includes the entire left half plane.

Definition:  The scheme is L stable if it is A stable and the amplification factor goes to zero as |ξ|->∞ in the left hand plane.

Definition:  The scheme is A[α] stable if the stability region includes the wedge with opening angle α centered on the negative real axis.