FWR-SIMP- Kepler problem

Kepler Problem

Clear[e,IC,f,t0,u0,nIt,p,MyPic,MmaPic,tmin,tmaxold]

e=3.0/4.0;(* value of e*)
IC={1-e,0.0,0.0,(Sqrt[((1+e)/(1-e))])};(*initail condition for y'=f*)
tmin=0.0;
tmaxold=2π;
f[{y1_,y2_,y3_,y4_}] :=Module[
{d},(*to find the function for y'=f for the orbit problem*)
d =((((y1^2)+(y2^2)))^(3/2)); 
{y3,y4,-y1/d,-y2/d}
]

The coefficient matrices for the integral

matrixA[r_]:=Module[{A,i,j},
A=ConstantArray[0,{r,r}];
For[i=2,i≤r-1,i++,
A[[i,1]]=1.0;
For[j=2,j≤i,j++,
A[[i+1,j]]=2.0;
];
A[[i,i]]=1.0;
];
A[[r,1]]=A[[r,r]]=1.0;
Return[A]
]

matrixB[r_,p_]:=Module[{A,i,j},
A=ConstantArray[0,{r,r-1}];
For[i=2,i≤r,i++,
For[j=2,j≤i,j++,
A[[i,j-1]]=1.0*p;
];
];
Return[A]
]

WR using Simpson' s rule

WRSIMP[r_,y_,e_]:=Module[{A,B,Yold,Ynew,F,Y0,h,a,g,k},
A=matrixA[r];
B=matrixB[r,4];
Y0= ConstantArray[y,{r}];
Ynew=Y0;
h=(1.0*e)/(2*r-2);
For[k=1,k≤150,k++,
F=Table[f[Ynew[[i]]],{i,1,r}];
g=Table[f[((Ynew[[i]]+Ynew[[i+1]])/2.0)+((h*(f[Ynew[[i]]]-f[Ynew[[i+1]]]))/4.0)],{i,1,r-1}];
Yold=Ynew;
Ynew=Y0+(h/3.0)*(A.F+B.g);
If[Norm[Yold-Ynew]≤(10^(-6)),Goto[endfor]];
];
Label[endfor];
Return[Ynew]
]

Mathematica solution using NDSolve

MmaSol=y/.NDSolve[ {
y'[t]==f[y[t]],
y[0]==IC
},
y,
{t, 0, tmaxold}][[1]];
MmaPic=ParametricPlot[MmaSol[t][[{1,2}]],{t,0,tmaxold},AxesLabel->{"(Subscript[y, 1])","(Subscript[y, 2])"}];

Run on the test problem:

Example2=Manipulate[
yin=MmaSol[tmin];
eval=Max[Abs[Eigenvalues[D[f[{y1,y2,y3,y4}],{{y1,y2,y3,y4}}]/.y1->yin[[1]]/.y2->yin[[2]]]]];
myPlot=ListLinePlot[WRSIMP[nodes,yin,H][[All,1;;2]],PlotStyle->{Red,PointSize[Large]},Mesh->All];
Grid[{{Show[MmaPic,myPlot]," "},{"","Max Eigenvalue:",eval}}],
{nodes,10,100,2},
{tmin,0,tmaxold},
{H,0.1,tmaxold}
];

Coservation of energy property

EnergeyExample=Manipulate[
Sol=WRSIMP[nodes,IC,H];
y11=Sol[[All,1]];
y22=Sol[[All,2]];
y33=Sol[[All,3]];
y44=Sol[[All,4]];
energy=Table[(1/2)(((y33[[i]])^2)+((y44[[i]])^2))-(1/(Sqrt[((y11[[i]])^2)+((y22[[i]])^2)])),{i,2,nodes-3}];
ListLinePlot[energy,PlotRange->Full,AxesLabel-> {"Time","Energy"},PlotLabel->"WR: Energy Conservation"],
{nodes,10,400,2},
{H,π,10π,π}
];

Example2

EnergeyExample