p29 Test Problem

Gravity[{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};
TMax=17.06521656015796;
ySol=y/.NDSolve[ {
y'[t]==Gravity[y[t]],
y[0]==IC
},
y,
{t, 0, TMax}][[1]];
MmaPic=ParametricPlot[ySol[t][[{1,2}]],{t,0,TMax}]

Interpolation Order Conditions

Quadrature

Numerical integration has a lot of different versions.  The simplest version say I want to integrate a smooth function over a fixed interval using a bunch of evenly spaced points. The easiest way to get formulas like Simpson's rule etc. is to compute formulae that are exact on polynomials.  the Taylor's theorme tells you that you can approximate the function to a certain order with a polynomial.  You then know that you can integrate the polynomial exactly and so you know how accurate your "quadrature" scheme is on smooth functions.

Here is an example computation which produces some familiar formulas

Clear[p,f,x,a,b,h,n]
n=2;
{a,b}={0,h};
Data=Table[{x,f[x]}, {x, a, b, h/n}]
p[x_]=InterpolatingPolynomial[Data,x]
Integrate[p[x],{x,a,b}]

You can do the same thing with unequally spaced points.

What does this have to do with solving ODE?

Differentiation

Differentiating interpolating polynomials is easy.  Here is an extrapolated derivative from a bunch of equally spaced points

p[x_]=InterpolatingPolynomial[
Table[
{xn+i h, (Subscript[y, n+i])},
{i, -4,-1}],
x];
Simplify[(h^2)p''[xn+{0,1,2,3}h]]//Expand

Adams Type Schemes

Adams schemes come from interpolating back values of the fs with polynomials and using that polynomial to approximate the integral in
	y[(Subscript[t, n])]=y[(Subscript[t, n-1])]+(∫_(τ=(Subscript[t, n-1]))^(Subscript[t, n]))f[y[τ]]dτ
Note:  I am using t as the time variable and τ as the dummy time variable in the integral.  We are going to use equally spaced points with spacing h and a variety of back values.
As a result Adams schemes do NOT involve any back values of y.

Clear[f,h]
{StartCount,BackCount}={0,3}; (* StartCount=1 gives explicit =0 gives implicit *)
Data=Table[
{tn - i h, (Subscript[f, n-i])},
{i, StartCount, BackCount}];
TraditionalForm[Data]
fPoly[τ_]=InterpolatingPolynomial[Data,τ];
Integrate[fPoly[τ],{τ, tn-h, tn}]

This generates the implicit scheme
	(Subscript[y, n])=(Subscript[y, n-1])+ h ((3/8) (Subscript[f, n])+(19/24)(Subscript[f, n-1])-(5/24)(Subscript[f, n-2])+(1/24)(Subscript[f, n-3]))
A direct computation gives the order of this scheme

Series[
y[x]-(y[x-h]+h( (3/8)f[y[x]]+ (19/24)f[y[x-h]]-(5/24)f[y[x-2h]]+(1/24)f[y[x-3h]])),
{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 is defined by its coefficients by

AM3={(3/8), {1},{(19/24),-(5/24),(1/24)}};
AM2={(5/12), {1},{(2/3),-(1/12)}};
AM4={(251/720), {1},{(323/360),-(11/30),(53/360),(19/720)}};
AM5={(95/288), {1},{(1427/1440),-(133/240),(241/720),-(173/1440),(3/160)}};

To get explicit schemes you follow the same procedure but you do NOT include the value (Subscript[f, n]) in the interpolation.

Clear[f,h]
{StartCount,BackCount}={1,3}; (* StartCount=1 gives explicit =0 gives implicit *)
Data=Table[
{tn - i h, (Subscript[f, n-i])},
{i, StartCount, BackCount}];
TraditionalForm[Data]
fPoly[τ_]=InterpolatingPolynomial[Data,τ];
Integrate[fPoly[τ],{τ, tn-h, tn}]

This give the explicit scheme
	(Subscript[y, n])=(Subscript[y, n-1])+ h ((23/12)(Subscript[f, n-1])-(4/3)(Subscript[f, n-2])+(5/12)(Subscript[f, n-3]))
which has error term

Series[
y[x]-(y[x-h]+h( (23/12)f[y[x-h]]-(4/3)f[y[x-2h]]+(5/12)f[y[x-3h]])),
{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]]
}

If I define this in terms of its coefficients I am going to write

AB3={0, {1},{(23/12),-(4/3),(5/12)}};
AB2 ={0, {1},{(3/2),((-1)/2)}};
AB4={0, {1},{(55/24),-(59/24),(37/24),-((3 )/3)}};
AB5={0, {1},{(1901/720),(1387/360),(109/30),(637/360),(251/720)}};

BDF Type Schemes

BDF schemes do NOT involve any back values of f.  They can be derived by differentiating polynomials which interpolating back values of y.  We will look at them in two ways.  The first is just as a formula for an implicit linear multi-step method that does not involve any back values for f i.e.
	(Subscript[y, n])=h (Subscript[β, 0])(Subscript[f, n])+(∑_(i=1)^k)(Subscript[α, i])(Subscript[y, n-i])
and asking what values of (Subscript[α, 1]),… would make this have reasonable accuracy.  Assuming that all the old values (Subscript[y, n-i]) are exactly correct i.e.
	(Subscript[y, n-i])=y((Subscript[x, n])- i h)  for i=1,2,3,4,…
and that f((Subscript[y, n]))=(Subscript[f, n]) gives 
	(Subscript[y, n])-y((Subscript[x, n]))=h (Subscript[β, 0])f((Subscript[y, n]))+(∑_(i=1)^k)(Subscript[α, i])y((Subscript[x, n])-i h)-y((Subscript[x, n])) 
The expansion is 
	(Subscript[y, n])-y((Subscript[x, n]))=((∑_(i=1)^k)(Subscript[α, i])-1)y((Subscript[x, n]))+h ((Subscript[β, 0])f((Subscript[y, n]))-(∑_(i=1)^k)i (Subscript[α, i]) y'((Subscript[x, n])) )+(1/2)(h^(2 ))(∑_(i=1)^k)(i^2) (Subscript[α, i]) y''((Subscript[x, n])) -(1/(3!))(h^3)(∑_(i=1)^k)(i^3) (Subscript[α, i]) y'''((Subscript[x, n])) 
So I can generate schemes of this sort that match out to a specified order by setting these terms to zero.  I get as equations
	({{(∑_(i=1)^k)(Subscript[α, i]), =, 1}, {(∑_(i=1)^k)i (Subscript[α, i]), =, (Subscript[β, 0])}, {(∑_(i=1)^k)(i^2) (Subscript[α, i]) , =, 0}, {(∑_(i=1)^k)(i^3) (Subscript[α, i]) , =, 0}, {…,  ,  }})
I can also think about looking at the problem of advancing from (Subscript[y, n-1])to (Subscript[y, n])using only back values (Subscript[y, n-i]) and the forward value (Subscript[f, n]).  The simplest BDF scheme has k=1 and is simply implicit Euler ((Subscript[α, 1])=(Subscript[β, 0])=1)
	(Subscript[y, n])=(Subscript[y, n-1])+h (Subscript[f, n])
which has expansion (as above) 
	(Subscript[y, n])-y((Subscript[x, n]))=0y((Subscript[x, n])) + 0 h y'((Subscript[x, n]))+(1/2)(h^2)y''((Subscript[x, n]))
We could also get a similar expansion by combining the estimate (Subscript[f, n])  and the estimate ((Subscript[y, n])-(Subscript[y, n-1]))/h for the derivative y'((Subscript[x, n]))  to get
	(Subscript[y, n])=(Subscript[y, n-1])+h( θ (Subscript[f, n])+(1-θ) ((Subscript[y, n])-(Subscript[y, n-1]))/h)=(Subscript[y, n-1])+(1-θ)((Subscript[y, n])-(Subscript[y, n-1]))+θ h (Subscript[f, n])
where the most meaningfull range of values is 0<θ≤1 with θ=1 giving Implicit Euler.
We could fancy this up by using more back values of (Subscript[y, n-i]) to estimate the second derivative
	y''((Subscript[x, n]))≈(1/(h^2))((Subscript[y, n])-2 (Subscript[y, n-1])+(Subscript[y, n-2]))  
multiplying by 1/2 (h^2) and adding it in gives a k=2 BDF scheme 
	(Subscript[y, n])= (Subscript[y, n-1])+(1-θ)((Subscript[y, n])-(Subscript[y, n-1]))+θ h (Subscript[f, n])-(θ/2)((Subscript[y, n])-2(Subscript[y, n-1])+(Subscript[y, n-2]))
which has error estimate
	-(h^3)(1/3)θ y'''[(Subscript[x, n])]
which we can eliminate by including one more back value and using this to estimate the third derivative
	((-(Subscript[y, -3+n])+3 (Subscript[y, -2+n])-3 (Subscript[y, -1+n])+(Subscript[y, n]))/(h^3))
which gives a k=3 BDF scheme 
	(Subscript[y, n])= (Subscript[y, n-1])+(1-θ)((Subscript[y, n])-(Subscript[y, n-1]))+θ h (Subscript[f, n])-(θ/2)((Subscript[y, n])-2(Subscript[y, n-1])+(Subscript[y, n-2]))-(θ/3)((Subscript[y, n])-3(Subscript[y, n-1])+3 (Subscript[y, n-2])-(Subscript[y, n-3]))
with error estimate 
	-(h^4)(1/4)θ y''''[x]
and we can continue for ever if we want.   Since it is easy to estimate the error and change the order these codes are frequently run in a variable order versions.  We are going to look at what happens with these schemes when we run them on real problems.

Reaaranging these formula to get all the (Subscript[y, n]) on the same side the variants with θ=1 are 
	(Subscript[y, n])=(Subscript[y, n-1])+h (Subscript[f, n])
	(Subscript[y, n])=
	defined by

BDF1={1,{1},{}};
BDF2={2/3,{(4/3),-(1/3)},{}};
BDF3={6/11,{(18/11),-(9/11),(2/11)},{}};

Computation Details

Data= Table[
{xn -i h, (Subscript[y, n-i])},
{i,0,3}]

p[x_]=InterpolatingPolynomial[Data,x];
(h^3)Simplify[p'''[xn]]

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

Data= Table[
{xn -i h, (Subscript[y, n-i])},
{i,0,3}]

p[x_]=InterpolatingPolynomial[Data,x];
(h^3)Simplify[p'''[xn]]

θ=1;
Expand[
Solve[
(Subscript[y, n])== (Subscript[y, n-1])+(1-θ)((Subscript[y, n])-(Subscript[y, n-1]))+θ h (Subscript[f, n])-(θ/2)((Subscript[y, n])-2(Subscript[y, n-1])+(Subscript[y, n-2])) -(θ/3)((Subscript[y, n])-3(Subscript[y, n-1])+3 (Subscript[y, n-2])-(Subscript[y, n-3])),
(Subscript[y, n])
][[1]]
]

LMM Stepper: Predictor and Fixed Point

I am going to make a stepper for a LMM scheme defined by two lists α and β and a separate β0 entry.  It is going to take the "ADDRESS" of the list as the argument and insert the appropriate entry.  For an explicit scheme this is going to insert the values 
	(Subscript[y, n])=(∑_(i=0)^k)(Subscript[α, i])(Subscript[y, n-i])+h(∑_(i=0)^k)(Subscript[β, i])(Subscript[f, n-i])
	(Subscript[f, n])=f((Subscript[y, n])) 
 and (Subscript[f, n]) into the lists using the appropriate number of backvalues.  For an implicit scheme this is going to run a single step of the fixed point iteration
	(Subscript[y, n]^(m+1))=h (Subscript[β, 0])(Subscript[f, n]^m)+(∑_(i=0)^k)(Subscript[α, i])(Subscript[y, n-i]^m)+h(∑_(i=0)^k)(Subscript[β, i])(Subscript[f, n-i]^m)
	(Subscript[f, n]^(m+1))=f((Subscript[y, n]^(m+1)))
The historical process runs the explicit Predictor step and then a number of Corrector steps.  Note that each of these steps has a function Evaluation. The overall process is typically described as ((PE(CE))^m) where m indicates the number of corrector steps.  It is possible to skip the last evaluation in which case we have ((PE(CE))^(m-1))C.  It is also possible to continue the corrector steps to convergence ((PE(CE))^∞).

Clear[LMMStep]
LMMStep[{ys_,fs_},n_,h_, f_,{β0_,α_,β_}]:= Module[
{fOld=fs[[n]],yNew},
yNew=h*β0*fOld+Sum[α[[i]] ys[[n-i]] ,{i, Length[α]}]+h Sum[β[[i]] fs[[n-i]],{i, Length[β]}];
ys[[n]]=yNew;
fs[[n]]=f[yNew];
]
SetAttributes[LMMStep,HoldFirst];

Now we are going to test it on our gravitational problem.  First we need to get enough starting values and then work out what step size will get us round the circuit. Looks like it works.

h=0.1;
nMax=Floor[TMax/h]
BackVals=6;
IterCount=6;
f=Gravity;
ys=fs=ConstantArray[0,{nMax}];
Do[
yTemp=ySol[i*h];
ys[[i]]=yTemp;
fs[[i]]=f[yTemp],
{i,1,BackVals}]
Do[ 
LMMStep[{ys,fs},n,h,f,AB5];
Do[
LMMStep[{ys,fs},n,h,f,AM5],
{IterCount}],
{n, BackVals+1, nMax,1}
]
ListPlot[ys[[All,{1,2}]]]

LMM Stepper: Newton - Not Finished

The obvious question is what to do for implicit schemes!  Fixed point iteration is common using an explicit scheme as a predictor for non-stiff problems.  Newton Iterations are common for stiff problems.

LMMStepNewt[ {ys,fs}, n_, h_, {f_, Df_}, {β0_,α_,β_}]:= Module[
{RHS, i,dim, yVal, fVal,Matrix,IterMax=23},
(* Need to solve y - β0 h f[y] == Sum[(Subscript[α, i]) (Subscript[ys, n-i])] + h Sum[(Subscript[β, i])(Subscript[fs, n-1])] *)
(* First forming the Right Hand Side *)
RHS= Sum[ α[[i]] ys[[n-i]] ,{i, Length[α]}]+h Sum[ β[[i]] fs[[i]],{i, Length[β]}];
dim=Length[RHS];
(* Could use a better initial guess *)
yVal=ys[[n-1]]; fVal=f[yVal];
Do[
Matrix= IdentityMatrix[dim]-h β0  Df[yVal];
yVal=yVal-LinearSolve[Matrix,yVal-β0 fVal-RHS];
fVal=f[yVal],
{IterMax}
];
ys[[n]]=yVal;
fs[[n]]=fVal
]

Df[{y1_,y2_,y3_,y4_,y5_,y6_}]=D[Gravity[{y1,y2,y3,y4,y5,y6}],{{y1,y2,y3,y4,y5,y6}}];

Left overs

DGravity[{y1_,y2_,y3_,y4_,y5_,y6_}]:= Module[
{DTemp=SparseArray[{},{6,6}]},
DTemp[[1,4]]=DTemp[[2,5]]=DTemp[[3,6]]=DTemp[[4,1]]=DTemp[[5,2]]=1;
DTemp[[4,5]]=2; DTemp[[5,4]]=-2;
DTemp
]