General Linear Methods

A general linear method combines elements of a multi-value scheme and an RK scheme into one procedure.  To be an Adams type scheme the multi-value passed from one time step to the next are derivatives of the fs.  To be a BDF type scheme the multi-values passed from one time to the next are derivatives of y.  There are lots of other possibilities:  for instance there are schemes that pass prior RK stage derivatives from one time to the next.

As always we are solving 
	y'=f[y]
with y[t]∈(ℝ^m).  We assume that we are going to pass r values at each step and have s internal stage valus and derivatives.

If r=1 we have an RK scheme.  If s=0 we have a LMM scheme.  A scheme consists of a description of the mapping
	(({{(Subscript[y, 1]^old)}, {(Subscript[y, 2]^old)}, {⋮}, {(Subscript[y, r]^old)}}))⟶(({{(({{(Subscript[y, 1]^int)}, {(Subscript[y, 2]^int)}, {(Subscript[y, 3]^int)}, {⋮}, {(Subscript[y, s]^int)}})),, (({{(Subscript[f, 1]^int)}, {(Subscript[f, 2]^int)}, {(Subscript[f, 3]^int)}, {⋮}, {(Subscript[f, s]^int)}}))}}))⟶(({{(Subscript[y, 1]^new)}, {(Subscript[y, 2]^new)}, {⋮}, {(Subscript[y, r]^new)}}))
As always it is desirable to have a consistent description of all such schemes to permit us to build consistent tools to implement error control, variable step sizes, orders, etc.

The overall plan is to build an extended Butcher tableau (({{A, U}, {B, V}})) of matrices A∈(ℝ^(s×s)), U∈(ℝ^sxr), B∈(ℝ^rxs), V∈(ℝ^rxr) which describes the scheme.  Consistent with the old Butcher table the top row describes the construction of the internal or stage values from the old input values while the bottom row describes the construction of the new output values from the internal stage values.   The only real changes are the matrix U that describes how the r  old y values {(Subscript[y, 1]^old), (Subscript[y, 2]^old), …,(Subscript[y, r]^old)} enter the internal computation and the matrix V which describes how the r old y values enter the computation of the new output values.

As always there are possibilities for including a lot of additional bells and whistles but the basic internal formula is
(Subscript[(GLM), in])	(Subscript[y, i]^int)=h(Underoverscript[ ∑, (j=1), s])(Subscript[A, i,j])(Subscript[f, j]^int)+(Underoverscript[ ∑, (j=1), r])(Subscript[U, i,j])(Subscript[y, j]^old)
for 1≤i≤s and the basic output formula is
(Subscript[(GLM), out])	(Subscript[y, i]^out)=h(Underoverscript[ ∑, (j=1), s])(Subscript[B, i,j])(Subscript[f, j]^int)+(Underoverscript[ ∑, (j=1), r])(Subscript[V, i,j])(Subscript[y, j]^old)
We can immediately write all of our LMM and RK schemes in this form. The schemes are explicit if A is strictly lower triangular and the schemes inherit much of their terminology from RK schemes. To code an implicit scheme we need to describe how we are going to solve (or approximately solve) the implicit equation (or equations) described by the stage computation in (Subscript[(GLM), in]).

There is considerable freedom as to what the "extra" values represent.  They could be errors, estimated derivatives, values at prior time steps, etc.  The scheme descibed in (GLM) does not specify what they are.

ExplicitGLMStepper[ y_,f_,h_,GLMTableau_]:= Module[
{A,U,B,V,r,s,yInt,fInt},
{{A,U},{B,V}}=GLMTableau[[1]];
s=Length[A];
r=Length[B];
yInt=fInt=ConstantArray[0,s];
Do[
yInt[[i]]=h Sum[A[[i,j]] fInt[[j]],{j,1,i-1}] + Sum[U[[i,j]] y[[j]],{j,1,r}];
fInt[[i]]=f[yInt[[i]]],
{i,1, s}
];
Table[
h Sum[ B[[i,j]] fInt[[j]],{j,1,s}] +Sum[V[[i,j]] y[[j]],{j,1,r}],
{i,1,r}
]
]

There are two schemes described on p126 which are not LMM or RK schemes. Their tableaus are

GLM261A={(({{(({{0, 0, 0}, {(1/2), 0, 0}, {-(1/3), (4/3), 0}})), (({{1, 0, 0, 0}, {1, 0, 0, 0}, {1, 0, 0, 0}}))}, {(({{(11/12), (1/3), (1/4)}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}})), (({{1, (1/12), -(1/3), -(1/4)}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}))}})),{"p126","{(y^n),h ((Subscript[f, 1])^n),h ((Subscript[f, 2])^n),h ((Subscript[f, 3])^n)}"}};
GLM262={(({{(({{0, 0, 0}, {(189/92), 0, 0}, {(3375/5152), (25/168), 0}})), (({{-(529/3375), (3904/3375), (4232/3375), (1472/3375)}, {(152/25), -(127/25), -(419/100), -(1118/575)}, {1, 0, (19/96), -(1/552)}}))}, {(({{(3375/5152), (25/168), 0}, {0, 0, 0}, {0, 0, 1}, {0, 0, 0}})), (({{1, 0, (19/96), -(1/552)}, {1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}))}})),{"p126","{(y^n),(y^(n - 1)),h y((')^n),h y((')^((n - 1)))}"}};

There are also a collection of RK schemes on p126.

GLM232A={(({{(({{0, 0}, {1, 0}})), (({{1}, {1}}))}, {(({{(1/2), (1/2)}})), {{1}}}})),{"p125","{(y^old)}"}};
GLM233f={(({{(({{0, 0, 0}, {(1/2), 0, 0}, {-1, 2, 0}})), (({{1}, {1}, {1}}))}, {(({{(1/6), (2/3), (1/6)}})), {{1}}}})),{"p125","{(y^old)}"}};
GLM235I={(({{(({{0, 0, 0, 0}, {(1/2), 0, 0, 0}, {0, (1/2), 0, 0}, {0, 0, 1, 0}})), (({{1}, {1}, {1}, {1}}))}, {(({{(1/6), (1/3), (1/3), (1/6)}})), {{1}}}})),{"p125","{(y^old)}"}};

Finally on p126 there is a AB/ABM PECE schemes in GLM form.

GLMPECE126={(({{(({{0, 0}, {(1/2), 0}})), (({{1, (3/2), -(1/2)}, {1, (1/2), 0}}))}, {(({{(1/2), 0}, {0, 1}, {0, 0}})), (({{1, (1/2), 0}, {0, 0, 0}, {0, 1, 0}}))}})),{"p126","{(y^n),h y((')^n), h y((')^((n - 1)))}"}};

None of these are currently equipped with error estimates but we can run them with a fixed stepsize on our p29 orbit problem and compare how large a step size we can take and still get visually plausible results.  We can also check the order of the error.   Here is a basic fixed stepsize driver.

ExplicitGLMDriverFixedStep[IC_, f_,{h_,n_}, GLMTableau_]:= Module[
{},
NestList[ ExplicitGLMStepper[ #,f,h,GLMTableau]&,IC,n]
]

The problem is that we need to have some mechanism for getting any genuinely multi-value scheme started.  We need to create the multi-value stuff.  I am going to do this manually for each scheme at the moment!

p29 Test Problem

Here is the orbit problem from p29

Gravityp129[{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)
}
]
ICp129 = {0.994,0,0,0,-2.001585106370825224,0};
TMax=17.06521656015796;
ySol=y/.NDSolve[ {
y'[t]==Gravityp129[y[t]],
y[0]==ICp129
},
y,
{t, 0, TMax},
PrecisionGoal->12,
AccuracyGoal->12][[1]];
MmaPic=ParametricPlot[ySol[t][[{1,2}]],{t,0,TMax}]

Testing Correctness and Order

First the GLM RK schemes.  For these the multi-value is trivial:  it is just a list of the single initial value

Clear[h]
IC={ySol[0]}
{a,pMin, pMax}={1.3,12,43};
Data=Table[
Table[
h=(a^(-p));
{h, Norm[ExplicitGLMStepper[IC,Gravity,h,Stepper][[1]]-ySol[h]]},
{p, pMin,pMax}
],
{Stepper, {GLM232A,GLM233f,GLM235I}}];
CompData=Table[ {{(a^(-pMin)),((((a^(-pMin))))^n)},{(a^(-pMax)),((((a^(-pMax))))^n)}},{n, 1, 7}];
Show[
ListLogLogPlot[Data],
ListLogLogPlot[CompData,
Joined->True]
]

Next the GLM LMM scheme.  For this the multi-value is non trivial

Clear[h,IC]
IC[h_]:={ySol[h], h ySol'[h], h ySol'[0]}
{a,pMin, pMax}={1.2,12,88};
Data = Table[
h=(a^(-p));
{h, Norm[ExplicitGLMStepper[IC[h],Gravity,h,GLMPECE126][[1]]-ySol[2h]]},
{p, pMin,pMax}];
CompData=Table[ {{(a^(-pMin)),((((a^(-pMin))))^n)},{(a^(-pMax)),((((a^(-pMax))))^n)}},{n, 1, 7}];
Show[
ListLogLogPlot[Data],
ListLogLogPlot[CompData,
Joined->True]
]

Next the GLM262

Clear[h,IC]
IC[h_]:={ySol[h], ySol[0],h ySol'[h], h ySol'[0]};
{a,pMin, pMax}={1.2,12,74};
Data = Table[
h=(a^(-p));
{h, Norm[ExplicitGLMStepper[IC[h],Gravity,h,GLM262][[1]]-ySol[2h]]},
{p, pMin,pMax}];
CompData=Table[ {{(a^(-pMin)),((((a^(-pMin))))^n)},{(a^(-pMax)),((((a^(-pMax))))^n)}},{n, 1, 7}];
Show[
ListLogLogPlot[Data],
ListLogLogPlot[CompData,
Joined->True]
]

Note that when the error is below (10^(-12)) our Mathematica numerical scheme is no longer sufficiently accurate.

Improved Order Computation and Picture: Not finished

It would be good to fit the central portion of the data and see how good the fit is.  This is not obvious with the current lines. The problem is that all the schemes need different starting methods.  So I can not really automate the set up of the error function I can however automate the extracting of the data, some fitting, and some plotting.

Clear[h]
IC={ySol[0]}
{a,pMin, pMax}={1.3,12,43};
Data=Table[
Table[
h=(a^(-p));
{h, Norm[ExplicitGLMStepper[IC,Gravity,h,Stepper][[1]]-ySol[h]]},
{p, pMin,pMax}
],
{Stepper, {GLM232A,GLM233f,GLM235I}}];
CompData=Table[ {{(a^(-pMin)),((((a^(-pMin))))^n)},{(a^(-pMax)),((((a^(-pMax))))^n)}},{n, 1, 7}];
Show[
ListLogLogPlot[Data],
ListLogLogPlot[CompData,
Joined->True]
]

RandomReal[]

Clear[OrderPlot]
OrderPlot[ErrFun_,{hMin_,hMax_},n_:100]:= Module[
{a=N[(((hMax/hMin))^(1/n))],Data, dataFit,h,hh,color=Hue[RandomReal[]]},
Data=Table[
{hMin (a^i), Norm[ErrFun[hMin (a^i)]]},
{i, 0,n}
];
dataFit[h_]=Simplify[(E^(Fit[Log[Data],{1,Log[h]},Log[h]]))];
Show[
ListLogLogPlot[Data,
GridLines->Automatic,
FrameLabel->{"h","err"},
Frame->True,
PlotStyle->color,
PlotLabel->dataFit["h"],
FrameTicks->{Select[Table[ (10.0^(-i)),{i,-100,100,1}] ,hMin<#<hMax&],Automatic,None, None}
],
LogLogPlot[dataFit[h],{h,hMin,hMax},PlotStyle->color]
]
]

Clear[hMin,hMax,n, ErrFun,h]
{hMin,hMax}={(10^(-6.5)),(10^(-2.5))};
n=30;
ErrFun1[h_]:=ExplicitGLMStepper[IC,Gravity,h,Stepper][[1]]-ySol[h];
OrderPlot[ErrFun1,{hMin,hMax},n]

Stability Computations (Single Step Format)

Clear[h,λ]
GLMTableau=GLM262;
{{A,U},{B,V}}=GLMTableau[[1]];
MValueLength=Length[V];
IC=Array[v,{MValueLength}];
fTest[y_]:=λ y

SingleStep=Simplify[ ExplicitGLMDriverFixedStep[IC,fTest,{1,1}, GLMTableau][[-1]] ]
{b,A}=Normal[CoefficientArrays[SingleStep, IC]];
G262Stab[λ_]=Max[Map[Abs,Eigenvalues[A]]];

Clear[h,λ]
GLMTableau=GLM233f;
{{A,U},{B,V}}=GLMTableau[[1]];
MValueLength=Length[V];
IC=Array[v,{MValueLength}];
fTest[y_]:=λ y

SingleStep=Simplify[ ExplicitGLMDriverFixedStep[IC,fTest,{1,1}, GLMTableau][[-1]] ]
{b,A}=Normal[CoefficientArrays[SingleStep, IC]];
G233fStab[λ_]=Max[Map[Abs,Eigenvalues[A]]]

Clear[h,λ]
GLMTableau=GLMPECE126;
{{A,U},{B,V}}=GLMTableau[[1]];
MValueLength=Length[V];
IC=Array[v,{MValueLength}];
fTest[y_]:=λ y

SingleStep=Simplify[ ExplicitGLMDriverFixedStep[IC,fTest,{1,1}, GLMTableau][[-1]] ]
{b,A}=Normal[CoefficientArrays[SingleStep, IC]];
GLMPECE126Stab[λ_]=Max[Map[Abs,Eigenvalues[A]]];

ContourPlot[{
G262Stab[x+I y]==1,
GLMPECE126Stab[x+I y]==1,
G233fStab[x+I y]==1
},
{x, -3,0.2},
{y,-3,3}
]