GK WR algorithm: version 1

So here is the simplest version of our algorithm.  All we do is implement a fixed point iteration using the matrix A.

Utilities

CreateWeightMatrix[absc_]:= Module[
{m=Length[absc],values,p,P,P0,t,i,j, ExtraPrecision=4},
Transpose[
Chop[
N[Table[
values=Table[ {absc[[i]],KroneckerDelta[i,j]},{i,m}];
p[t_]= Expand[InterpolatingPolynomial[ values, t]];
P[t_]=Integrate[p[t],t];
P0=P[0];
Table[P[absc[[i]]]-P0,{i,1,m}],
{j,1,m}
]
]
]
]
]
WriteMatrix[A_,n_,type_]:= Module[{name},
name=ToString[StringForm["WeightMat_``_``.dat",n,type]];
Export[name,A]
]

P29 Gravity Problem

Here is our standard test problem.

Clear[GravityP29]
GravityP29[{y1_,y2_,y3_,y4_,y5_,y6_}] :=Module[
{μ=100/8145, 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 = SetPrecision[{0.994,0,0,0,-2.001585106370825224,0},64];
IC=IC1;
TMax=17.06521656015796; 
ySol[t_]=y[t]/.NDSolve[{
y'[t]==GravityP29[y[t]],
y[0]==IC
},
y[t],
{t,0, TMax},
PrecisionGoal->16,
AccuracyGoal->16,
InterpolationOrder->16,
WorkingPrecision->64][[1]];
TStart=TMax/3;
H=TMax/3;
MmaPic=ParametricPlot[ySol[t][[{1,2}]],{t,0,TMax},
Epilog->{PointSize[0.02],Pink, 
Point[ySol[TStart][[{1,2}]]],
Point[ySol[TStart + H][[{1,2}]]]
}
]

GKWR Algorithm

Here is a test of a specific version of the agorithm.

First we build the matrix and interpolation points.

n=21;
precision=256;
absck=NIntegrate`GaussKronrodRuleData[n,precision][[1]];
(* absck=Partition[absck,2][[All,1]]; *)
ListPlot[absck]
A=CreateWeightMatrix[absck];
MatrixPlot[A]
ListPlot3D[A]

Now we test.

Y0=Table[ySol[TStart],{Length[absck]}];
Y=Y0;
IterCount=100;
Data=Table[
F=Map[GravityP29,Y];
Y=Y0+H A.F,
{IterCount}];
Show[MmaPic,
ListPlot[Y[[All,{1,2}]],
PlotStyle->Red]
]

An obvious question is does it really needs all the iterations.  Here is a picture of how it is converging.

Manipulate[
Show[
MmaPic, 
ListPlot[
Data[[i,All,{1,2}]],
Joined->True,
PlotStyle->Red,
Mesh->All
]
],
{i, 1,Length[Data],1}
]

Here is a pic of how far away it is from the “converged” value.  I apparently stall in my convergence at some point and wind up bouncing zeros.

ListContourPlot[
Table[
Log[10,Norm[Data[[i,j]]-Data[[-1,j]]]],
{i,Length[Data]-1},
{j, Length[absck]}],
InterpolationOrder->1,
ContourLabels->(Text[Framed[#3],{#1,#2},Background->White]&)
]

A second obvious question is how far away am I from the NDSolve answer?  In this case the answer is about (10^(-4)) which is not that impressive.

ListContourPlot[
Table[
Log[10,Norm[Data[[i,j]]-ySol[TStart+H*absck[[j]]]]],
{i,Length[Data]-1},
{j, Length[absck]}],
InterpolationOrder->1,
ContourLabels->(Text[Framed[#3],{#1,#2},Background->White]&)
]

There are convergence questions: Does it converege (and if so how fast) as a function of stiffness, n and H?    There are accuracy questions: How accurate is the converged (assuming it converges) result as a function of n and H?   There are efficiency questions:  If I have m concurrent processors is this more or less than efficient than other purely serial codes as a function of n, m, stiffness and H.  There are control questions:  How should I pick H, n, etc.

Some theory would seem to be in order.

Convergence Theory for linear iterative schemes.

For a linear function f[y]=λy we have the iteration
	(Y^(n+1)) = Y0 + H λ A.(Y^n)
which (if it converges) converges to
	(Id-H λA)(Y^∞)=Y0
which we might rewrite as 
	(Y^∞)=(((Id-Hλ A))^(-1))Y0 
where we would have to worry a little bit about what some of this means.  Subtracting the limiting value gives
	(Y^(n+1))-(Y^∞)=Y0+H λ A (Y^n)-(Y^∞)
which we can rewrite as
	(Y^(n+1))-(Y^∞)=((Y^∞)-H λ A (Y^∞))+H λ A (Y^n)-(Y^∞)=H λ A ((Y^n)-(Y^∞))
which is 
	(W^(n+1))=H λ A (W^n)
where 
	(W^n)=(Y^n)-(Y^∞).
So the convergence of this thing is governed by  the iteration matrix H λ A.  As in the linear stability arguments we put H λ=z and remember that it could well be complex.

The eigenvalues of A are relatively easy to plot.

{evals,evecs}=Eigensystem[A];
ListPlot[
evals/.{
z_Complex:>{Re[z],Im[z]},
z_Real:> {z,0}
},
AspectRatio->Automatic,
PlotRange->All,
PlotLabel->Max[Abs[evals]],
Epilog->{Red,Circle[{0,0},Max[Abs[evals]]]}
]
ListLogPlot[Abs[evals],
PlotRange->All]

We can see that we could multiply by about 43 (roughly 1/0.023) before we had an eigenvalue bigger than one.

evals[[1;;4]]
ListLogPlot[{Re[evecs[[1]]],-Re[evecs[[1]]],Im[evecs[[1]]],-Im[evecs[[1]]]},
PlotRange->All,
PlotStyle->{Red,Pink,Blue,Cyan}]
ListLogPlot[{Re[evecs[[3]]],-Re[evecs[[3]]],Im[evecs[[3]]],-Im[evecs[[3]]]},
PlotRange->All,
PlotStyle->{Pink,Red,Cyan,Blue}]
ListLogPlot[{Re[evecs[[5]]],-Re[evecs[[5]]],Im[evecs[[5]]],-Im[evecs[[5]]]},
PlotRange->All,
PlotStyle->{Pink,Red,Cyan,Blue}]

Here is the convergence for the simple linear case

fLinear[y_]=0.7 y;
Y0=Table[{1.,3.1,2.1},{Length[absck]}];
Y=Y0;
IterCount=32;
Data=Table[
F=Map[fLinear,Y];
Y=Y0+H A.F,
{IterCount}];
ListContourPlot[
ConvData=Table[
Log[10,Norm[Data[[i,j]]-Data[[-1,j]]]],
{i,Length[Data]-1},
{j, Length[absck]}],
InterpolationOrder->1,
ContourLabels->(Text[Framed[#3],{#1,#2},Background->White]&)
]