Gauss-Kronrod and other Integration

Notes: 
A) To get full double precision I need to form the matrices in higher precision.  This would be true for any of the higher order integration schemes. 
B) Somehow I need to use a version other than the default load of 7.0.  My current 8.0.1 load from Feb 23 2011 is fine to high order. 
C) I need to investigate the Gauss-Seidel implementation as well as the current Jacobi iteration.  On a GPU it is probably a block intermediate that is actually running. 
D) Note: lack of smoothness affects the accuracy badly as does an integrable singularity.  This is not much of a surprise!
E) There must be equivalent integration formulas based on Fourier series integration.  I need to find out/remember the names. 
F) Note:  The full interpolation formula for large n does take a little while with an analytic integrattion but it is a onesy computation and could easily be acclerated.  the easy solution is to use NIntegrate with high precision.Completely excessive but n=32 is taking a while!  n=7 takes 6 seconds and n=8 takes 9 seconds.

$Version

Making matrices for a bunch of dimensions and all the available integration formats is.

I have access to a few one dimensional integration rules.  See the list below.

CreateWeightMatrix[absc_]:= Module[
{m=Length[absc],values,p,P,P0,t,i,j, ExtraPrecision=4},
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}
],
$MachinePrecision+ExtraPrecision
],
(10^(-$MachinePrecision))
]
]
WriteMatrix[A_,n_,type_]:= Module[{name},
name=ToString[StringForm["WeightMat_``_``.dat",n,type]];
Export[name,A]
]

SetDirectory["/home/faculty/struther/Desktop/Link to Courses/5628_11/Data"];
precision=64;  (* Note:  Some extra precision is neccessary but this is prob excessive *)
{nMin, nMax}={3,12}; (* Note:  Actual array sizes are different e.g GK and GBE are 2n+1 square while CC and NC and LK and LP are 2n-1 square*)
TimingData=Table[
Timing[
Do[
absc=ToExpression[ToString[StringForm["NIntegrate`.```RuleData[``,``][[1]]",ToExpression[type],n,precision]]];
(* Print[StringForm["For n = ``: Length of `` is ``",n, type, Length[absc]]]; *)
A=CreateWeightMatrix[absc];
WriteMatrix[A,n,type],
{type, {"ClenshawCurtis","GaussBerntsenEspelid","GaussKronrod","LobattoKronrod","LobattoPeano","NewtonCotes"}}
]
][[1]],
{n,nMin, nMax}
];
ListLogPlot[TimingData]

F[n_]=Exp[Fit[Transpose[{Range[nMin,nMax],Log[TimingData]}],{1,n},n]]
F[15]

Data={{1.9916959999999335`,Null},{7.183908000000031`,Null},{25.754085000000032`,Null},{105.02903299999997`,Null}};
ListLogPlot[Transpose[{{3,4,5,6},Data[[All,1]]}]]

The prediction is it takes about an hour to make the n=31 data sets.  This is perfectly feasible.

On the interval (∫_0^1)f[τ]dτ

n=7;
precision=$MachinePrecision;
precision=64; 
{abscK,weightsK,errweightsK}=NIntegrate`GaussKronrodRuleData[n,precision];

f[t_]:= ((t^3)/(1+(t^2)));
FK=f[abscK];
err=Abs[NIntegrate[f[t],{t,0,1}]-weightsK.FK]

On the interval (∫_0^H)f[τ]dτ

All we have to do is scale the evaluation points and the integral

n=7;
precision=$MachinePrecision;
{abscK,weightsK,errweightsK}=NIntegrate`GaussKronrodRuleData[n,precision];
f[t_]:= ((Sin[(t^2)]Abs[t+2](((Abs[t-2]))^(11/3)))/(((Abs[t+1]))^0.3));
H=3.4;
FK=f[H*abscK];
Show[
Plot[f[t],{t,0,H},
Filling->Axis,
PlotRange->All],
ListPlot[Transpose[{H*abscK,FK}],
AxesLabel->{"t","f"}],
PlotRange->All
]
{
{GKInt=H*FK.weightsK,GKErr=H*FK.errweightsK}, 
{MmaInt= NIntegrate[f[t],{t,0,H}],(1-GKInt/MmaInt),(MmaInt-GKInt)}
}

Interpolation formulae.

These integration formulas are based on polynomial interpolation.  I need the equivalent formula for the integration up to each of the absicae. Fortunately, this is easy.  Note the extended precision in the computation to get full double precision in the output matrix A..  It is striking that the matrix is almost triangular with a low amplitude ± checkerboard pattern on the “small” side.  Note these were timed and reported at the top.  Do not run for n=32.  I am using quite excessive extra precision.

The dimensions that make sense here are simple: stuff that is less than or eaul to 16≥8+7, 32≥15+16, and 64≥31+32. 
I need to see what makes sense for the other schemes.

n=31;
ExtraPrecision=128;
precision=$MachinePrecision+ExtraPrecision; 
{abscK,weightsK,errweightsK}=NIntegrate`GaussKronrodRuleData[n,precision];
Clear[f,F]
poly[t_]=InterpolatingPolynomial[Table[ {abscK[[i]],(Subscript[f, i])},{i, Length[abscK]}],t];

Poly[t_,j_]:= Expand[InterpolatingPolynomial[Table[{abscK[[i]],KroneckerDelta[i,j]},{i, Length[abscK]}],t]]
{t1,prec}=Timing[
Precision[a=Integrate[Poly[t,2 n+1],{t, 0,abscK[[2n+1 ]]}]]]

A31=Table[ 
Integrate[ Poly[t,j],{t,0, abscK[[i]]}],
{i, Length[abscK]},
{j,Length[abscK]}
];

A=Transpose[A31];
ListContourPlot[Map[Precision,A,{2}],
PlotLabel->"Precision of A31",
FrameLabel->{"i","j"}]
{ListPlot3D[A,
PlotRange->All,
AxesLabel->{"i","j"},
PlotLabel->"A31"],
MatrixPlot[A,
FrameLabel->{"i","j"},
PlotLabel->"A31"]
}

SetDirectory["/home/faculty/struther/Desktop/Link to Courses/5628_11"]
A=N[Import["AGK31.dat"]];
Dimensions[A];
n=Length[A];
Id = IdentityMatrix[n];
ListLogPlot[{Abs[Eigenvalues[A]],SingularValueList[A]}]
ListContourPlot[A.Transpose[A]]
Min[A.Transpose[A]]

F=Table[(Subscript[f, i]),{i,1,Length[abscK]}];
A=Table[
 Coefficient[Simplify[Integrate[poly[t],{t, 0, abscK[[i]]}]],F],
{i,1, Length[abscK]}];
ListContourPlot[Map[Precision,A,{2}],
PlotLabel->"Precision of A",
FrameLabel->{"i","j"}]
{ListPlot3D[A,
PlotRange->All,
AxesLabel->{"i","j"},
PlotLabel->"A"],
MatrixPlot[A,
FrameLabel->{"i","j"},
PlotLabel->"A"]
}

Testing the integration formulae is pretty easy. All I need to do is do the integration for each absicca and plug in

Clear[f,H,t,FK,errs]
H=3.0;
f[t_]:= ((Sin[t+ Cos[ t+2]])/(1+(t^2)))
FK=f[H*abscK];
errs=H A.FK-Table[ 
NIntegrate[ f[t],{t,0, H*abscK[[i]]},PrecisionGoal->16,AccuracyGoal->16],
{i, Length[abscK]}];
ListLogPlot[{errs,-errs},
PlotLabel->StringForm["`` points",2 n+1],
AxesLabel->{"point","Int Err"},
PlotRange->All
]

The errors reported look remarkably flat.  for large H they are typically ± patterned!

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.

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 fixed iteration pass through the fixed point iteration.  I think this is pretty impressive.  It is doing a large section of the loop with a few points

I ned to improved this till it satisfies a tolerance

Possible “rules” are

NIntegrate::bdmtd: The Method option should be one of Automatic, CartesianRule, ClenshawCurtisRule, GaussBerntsenEspelidRule, GaussKronrodRule, LobattoKronrodRule, LobattoPeanoRule, MultiPanelRule, NewtonCotesRule, TrapezoidalRule, AdaptiveMonteCarlo, AdaptiveQuasiMonteCarlo, DoubleExponential, DuffyCoordinates, ExtrapolatingOscillatory, MonteCarlo, QuasiMonteCarlo, Trapezoidal, EvenOddSubdivision, SymbolicPiecewiseSubdivision, OscillatorySelection, UnitCubeRescaling, or a user-implemented integration strategy or rule.

Rule

({{"BooleRule", one-dimensional, weighted sum}, {"ClenshawCurtisRule", one-dimensional, weighted sum}, {"GaussBerntsenEspelidRule", one-dimensional, weighted sum}, {"GaussKronrodRule", one-dimensional, weighted sum}, {"LobattoKronrodRule", one-dimensional, weighted sum}, {"LobattoPeanoRule", one-dimensional, weighted sum}, {"NewtonCotesRule", one-dimensional, weighted sum}, {"TrapezoidalRule", one-dimensional, weighted sum}, {"ClenshawCurtisOscillatoryRule", one-dimensional, specialized oscillatory rule}, {"LevinRule", one- or multidimensional, general oscillatory rule}, {"MultipanelRule", one-dimensional, weighted sum, combination of 1D rules}, {"CartesianRule", multidimensional, weighted sum, product of 1D rules}, {"MultidimensionalRule", multidimensional, weighted sum}})

RuleList={"ClenshawCurtisRule","GaussBerntsenEspelidRule","GaussKronrodRule","LobattoKronrodRule","LobattoPeanoRule","NewtonCotesRule"}

n=5;
ExtraPrecision=128;
precision=$MachinePrecision+ExtraPrecision;
{ 
{absc,weights,errweights}=NIntegrate`ClenshawCurtisRuleData[n,precision];
{"CC",Length[absc]},
{absc,weights,errweights}=NIntegrate`GaussBerntsenEspelidRuleData[n-1,precision];
{"GBE",Length[absc]},
{absc,weights,errweights}=NIntegrate`GaussKronrodRuleData[n-1,precision];
{"GK",Length[absc]},
{absc,weights,errweights}=NIntegrate`LobattoKronrodRuleData[n,precision];
{"LKR",Length[absc]},
{absc,weights,errweights}=NIntegrate`NewtonCotesRuleData[n,precision];
{"NC",Length[absc]}
}

n=42;
NC=NIntegrate`NewtonCotesRuleData[n,precision][[1]];
ListPlot[
{
NC-NIntegrate`ClenshawCurtisRuleData[n,precision][[1]],
NC-NIntegrate`GaussBerntsenEspelidRuleData[n-1,precision][[1]],
NC-NIntegrate`GaussKronrodRuleData[n-1,precision][[1]],
NC-NIntegrate`LobattoKronrodRuleData[n,precision][[1]]
}
]

f[t_]= t;
NIntegrate[f[t],{t,0,1},
Method->dump]

Y0=Table[ySol[TStart],{Length[abscK]}];
Y=Y0;
IterCount=300;
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]&)
]

g[x_]=Fit[ Log[Take[ConvData,{10,50}]],{1,x},x]

Simplify[(E^(g[x]))]

ListPlot[ConvData]

Show[
MmaPic,
ListPlot[{Y[[All,{1,2}]], Table[ySol[(TMax/3)+H*abscK[[i]]][[{1,2}]],{i,1,Length[abscK]}]}]
]