WAVEFORM RELAXATION CODE

Team : Ahmad + Mohammad

Test Problems

Gravitational Problem

gravIC={0.994,0.0,0.0,0.0,-2.0015851063790825224,0.0};
gravTmin=0.0;
gravTmax=17.06521656015796;
grav[]={-∞,-∞};
grav[n_][{x1_,x2_,x3_,x4_,x5_,x6_}]:=Evaluate[
Module[{f,df,X1,X2,X3,X4,X5,X6,μ=1/81.45},
If[grav[]=={-∞,-∞},
f[X1_,X2_,X3_,X4_,X5_,X6_]={
X4,
X5,
0,
2X5+X1-((μ(X1+μ-1))/(((X2^2+(((X1+μ-1))^2)))^(3/2)))-(((1-μ)(X1+μ))/(((X2^2+(((X1+μ))^2)))^(3/2))),
-2X4+X2-((μ*X2)/(((X2^2+(((X1+μ-1))^2)))^(3/2)))-(((1-μ)X2)/(((X2^2+(((X1+μ))^2)))^(3/2))),
0
};
df[X1_,X2_,X3_,X4_,X5_,X6_]=D[f[X1,X2,X3,X4,X5,X6],{{X1,X2,X3,X4,X5,X6}}];
grav[]={0,0};
];
Switch[n,
0,grav[]=grav[]+{1,0};Evaluate[f[x1,x2,x3,x4,x5,x6]],
1,grav[]=grav[]+{0,1};Evaluate[df[x1,x2,x3,x4,x5,x6]]
]
]
]

Kepler Problem

e=3.0/4.0;
keplerIC={1-e,0.0,0.0,(((((1+e)/(1-e))))^(1/2))};
keplerTmin=0.0;
keplerTmax=2π;
kepler[]={-∞,-∞};
kepler[n_][{x1_,x2_,x3_,x4_}]:=Evaluate[
Module[
{f,df,X1,X2,X3,X4},
If[kepler[]=={-∞,-∞},
f[X1_,X2_,X3_,X4_]={X3,X4,-X1/((((X1^2)+(X2^2)))^(3/2)),-X2/((((X1^2)+(X2^2)))^(3/2))};
df[X1_,X2_,X3_,X4_]=D[f[X1,X2,X3,X4],{{X1,X2,X3,X4}}];
kepler[]={0,0};
];
Switch[n,
0,kepler[]=kepler[]+{1,0};Evaluate[f[x1,x2,x3,x4]],
1,kepler[]=kepler[]+{0,1};Evaluate[df[x1,x2,x3,x4]]
]
]
]

Dynamic Windowing WR Steppers and Driver

The Matrices A and B

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_]:=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;
];
];
Return[A]
]

Simpson stepper

DWRSIMP[n_,r_,y_,H0_,f_]:=Module[{A,B,Y0,h,k,F,g,Ynew},
A=matrixA[r];
B=4*matrixB[r];
Y0= ConstantArray[y,{r}];
Ynew=Y0;
h=H0/(r-1);
For[k=1,k≤n,k++,
F=Table[f[0][Ynew[[i]]],{i,1,r}];
g=Table[f[0][((Ynew[[i]]+Ynew[[i+1]])/2.0)+(1/8.0)h*(f[0][Ynew[[i]]]-f[0][Ynew[[i+1]]])],{i,1,r-1}];
Ynew=Y0+(h/6.0)*(A.F+B.g)
];
Return[Ynew]
]

Trapezoid Stepper

DWRTRAP[n_, r_, y_, H0_, f_] := Module[{A, Y0, Ynew, h, k, F}, 
   A = matrixA[r];
Y0 = ConstantArray[y, {r}]; 
Ynew = Y0;
h = H0/(r - 1); 
  For[k = 1, k <= n, k++, 
F = Table[f[0][Ynew[[i]]], {i, 1, r}];
Ynew = Y0 + (h/2.)*A . F
]; 
  Return[Ynew]
];

Midpoint Stepper

DWRMID[n_,r_,y_,H0_,f_]:=Module[{B,Y0,Ynew,h,k,g},
B=matrixB[r];
Y0= ConstantArray[y,{r}];
Ynew=Y0;
h=H0/(r-1);
For[k=1,k≤n,k++,
g=Table[f[0][((Ynew[[i]]+Ynew[[i+1]])/2.0)+(1/8.0)h*(f[0][Ynew[[i]]]-f[0][Ynew[[i+1]]])],{i,1,r-1}];
Ynew=Y0+h*B.g;
];
Return[Ynew]
]

Driver

DWRDriver[t0_,tmax_,y0_,r_,n_,Tol_,stepper_,f_]:=Module[{ya,λ,yData,t=t0,H,Ynew},
yData={};
ya=y0;
λ=Max[Abs[Eigenvalues[f[1][ya]]]];
While[t≤ tmax,
          H=0.5*(((((n+1)!*Tol))^(1/(n+1)))/λ); (* the 0.5 is a safety factor*)
          t=t+H;
         If[t<tmax,
           Ynew=stepper[n,r,ya,H,f];
           yData=Join[yData,Ynew];
           ya=Ynew//Last;
           λ=Max[Abs[Eigenvalues[f[1][ya]]]],
           H=tmax-(t-H);
          Ynew=stepper[n,r,ya,H,f];
          yData=Join[yData,Ynew];
          Goto[endwhile]
      ]
];
Label[endwhile];
yData
]

Fixed Window WR

Simpson stepper (fixed)

FWRSIMP[n_,r_,Tol_,y_,H0_,f_]:=Module[{A,B,Yold,Ynew,Y0,F,g,h,k},
A=matrixA[r];
B=4*matrixB[r];
Y0= ConstantArray[y,{r}];
Ynew=Y0;
h=H0/(r-1);
For[k=1,k≤n,k++,
F=Table[f[0][Ynew[[i]]],{i,1,r}];
g=Table[f[0][((Ynew[[i]]+Ynew[[i+1]])/2.0)+(1/8.0)h*(f[0][Ynew[[i]]]-f[0][Ynew[[i+1]]])],{i,1,r-1}];
Yold=Ynew;
Ynew=Y0+(h/6.0)*(A.F+B.g);
If[Norm[Yold-Ynew]≤Tol,Goto[endfor1]];
];
Label[endfor1];
Return[Ynew]
];

Trapezoid Stepper (fixed)

FWRTRAP[n_,r_,Tol_,y_,H0_,f_]:=Module[{A,Yold,Ynew,F,Y0,h,k},
A=matrixA[r];
Y0= ConstantArray[y,{r}];
Ynew=Y0;
h=H0/(r-1);
For[k=1,k≤n,k++,
F=Table[f[0][Ynew[[i]]],{i,1,r}];
Yold=Ynew;
Ynew=Y0+(h/2.0)*A.F;
If[Norm[Yold-Ynew]≤Tol,Goto[endfor2]];
];
Label[endfor2];
Return[Ynew]
];

Midpoint Stepper (fixed)

FWRMID[n_,r_,Tol_,y_,H0_,f_]:=Module[{B,Yold,Ynew,Y0,F,h,g,k},
B=matrixB[r];
Y0= ConstantArray[y,{r}];
Ynew=Y0;
h=H0/(r-1);
For[k=1,k≤n,k++,
F=Table[f[0][Ynew[[i]]],{i,1,r}];
g=Table[f[0][((Ynew[[i]]+Ynew[[i+1]])/2.0)+(1/8.0)h*(f[0][Ynew[[i]]]-f[0][Ynew[[i+1]]])],{i,1,r-1}];
Yold=Ynew;
Ynew=Y0+h*B.g;
If[Norm[Yold-Ynew]≤Tol,Goto[endfor3]];
];
Label[endfor3];
Return[Ynew]
];

Numerical Test

Test 1: Gravitational Problem

Mathematica Solution

MmaSolGrav=y/.NDSolve[ {
y'[t]==grav[0][y[t]],
y[gravTmin]==gravIC
},
y,
{t, gravTmin, gravTmax}][[1]];
MmaPicGrav=ParametricPlot[MmaSolGrav[t][[{1,2}]],{t,gravTmin,gravTmax}];

DWR-Simpson

grav[]={0,0};
Time1=AbsoluteTiming[
pic1=ListPlot[DWRDriver[gravTmin,gravTmax,gravIC,6,20,(10^(-12)),DWRSIMP,grav][[All,1;;2]],PlotStyle->{Black,PointSize[Medium]},Mesh->Full,AxesLabel-> {"(Subscript[y, 1])","(Subscript[y, 2])"},PlotLabel->"DWR-SIMP",PlotRange->All,AspectRatio->1];
][[1]];
count1=grav[][[1]];
Show[pic1,MmaPicGrav]

DWR-Trapezoid

grav[]={0,0};
Time2=AbsoluteTiming[
pic2=ListPlot[DWRDriver[gravTmin,gravTmax,gravIC,10,10,(10^(-12)),DWRTRAP,grav][[All,1;;2]],PlotStyle->{Black,PointSize[Medium]},Mesh->Full,AxesLabel-> {"(Subscript[y, 1])","(Subscript[y, 2])"},PlotLabel->"DWR-TRAP",PlotRange->All,AspectRatio->1];
][[1]];
count2=grav[][[1]];
Show[pic2,MmaPicGrav]

DWR-Midpoint

grav[]={0,0};
Time3=AbsoluteTiming[
pic3=ListPlot[DWRDriver[gravTmin,gravTmax,gravIC,10,10,(10^(-12)),DWRMID,grav][[All,1;;2]],PlotStyle->{Black,PointSize[Medium]},Mesh->Full,AxesLabel-> {"(Subscript[y, 1])","(Subscript[y, 2])"},PlotLabel->"WR-MIDPOINT",PlotRange->All,AspectRatio->1];
][[1]];
count3=grav[][[1]];
Show[pic3,MmaPicGrav]

FWR-Simpson

grav[]={0,0};
Time4=AbsoluteTiming[
pic4=ListPlot[FWRSIMP[200,200,(10^(-8)),MmaSolGrav[0.1],gravTmax-0.2,grav][[All,1;;2]],PlotStyle->{Black,PointSize[Medium]},Mesh->Full,AxesLabel-> {"(Subscript[y, 1])","(Subscript[y, 2])"},PlotLabel->"FWR-SIMP",PlotRange->All,AspectRatio->1];
][[1]];
count4=grav[][[1]];
Show[pic4,MmaPicGrav]

FWR-Trapezoid

grav[]={0,0};
Time5=AbsoluteTiming[
pic5=ListPlot[FWRTRAP[200,650,(10^(-8)),MmaSolGrav[0.1],gravTmax-0.2,grav][[All,1;;2]],PlotStyle->{Black,PointSize[Medium]},Mesh->Full,AxesLabel-> {"(Subscript[y, 1])","(Subscript[y, 2])"},PlotLabel->"FWR-TRAP",PlotRange->All,AspectRatio->1];
][[1]];
count5=grav[][[1]];
Show[pic5,MmaPicGrav]

FWR-Midpoint

grav[]={0,0};
Time6=AbsoluteTiming[
pic6=ListPlot[FWRMID[200,450,(10^(-8)),MmaSolGrav[0.1],gravTmax-0.2,grav][[All,1;;2]],PlotStyle->{Black,PointSize[Medium]},Mesh->Full,AxesLabel-> {"(Subscript[y, 1])","(Subscript[y, 2])"},PlotLabel->"FWR-MID",PlotRange->All,AspectRatio->1];
][[1]];
count6=grav[][[1]];
Show[pic6,MmaPicGrav]

Test 2: Kepler Problem

Mathematica Solution

MmaSolKepler=y/.NDSolve[ {
y'[t]==kepler[0][y[t]],
y[keplerTmin]==keplerIC
},
y,
{t, keplerTmin, keplerTmax}][[1]];
MmaPicKepler=ParametricPlot[MmaSolKepler[t][[{1,2}]],{t,keplerTmin,keplerTmax}];

DWR-Simpson

kepler[]={0,0};
Time7=AbsoluteTiming[
pic7=ListPlot[DWRDriver[keplerTmin,keplerTmax,keplerIC,5,20,(10^(-12)),DWRSIMP,kepler][[All,1;;2]],PlotStyle->{Black,PointSize[Medium]},Mesh->Full,AxesLabel-> {"(Subscript[y, 1])","(Subscript[y, 2])"},PlotRange->All,PlotLabel->"DWR-SIMP"];
][[1]];
count7=kepler[][[1]];
Show[pic7,MmaPicKepler]

DWR-Trapezoid

kepler[]={0,0};
Time8=AbsoluteTiming[
pic8=ListPlot[DWRDriver[keplerTmin,keplerTmax,keplerIC,20,20,(10^(-12)),DWRTRAP,kepler][[All,1;;2]],PlotStyle->{Black,PointSize[Medium]},Mesh->Full,AxesLabel-> {"(Subscript[y, 1])","(Subscript[y, 2])"},PlotLabel->"DWR-TRAP",PlotRange->All];
][[1]];
count8=kepler[][[1]];
Show[pic8,MmaPicKepler]

DWR-Midpoint

kepler[]={0,0};
Time9=AbsoluteTiming[
pic9=ListPlot[DWRDriver[keplerTmin,keplerTmax,keplerIC,20,20,(10^(-12)),DWRMID,kepler][[All,1;;2]],PlotStyle->{Black,PointSize[Medium]},Mesh->Full,AxesLabel-> {"(Subscript[y, 1])","(Subscript[y, 2])"},PlotLabel->"DWR-MIDPOINT",PlotRange->All];
][[1]];
count9=kepler[][[1]];
Show[pic9,MmaPicKepler]

FWR-Simpson

kepler[]={0,0};
Time10=AbsoluteTiming[
pic10=ListPlot[FWRSIMP[200,70,(10^(-8)),keplerIC,keplerTmax,kepler][[All,1;;2]],PlotStyle->{Black,PointSize[Medium]},Mesh->Full,AxesLabel-> {"(Subscript[y, 1])","(Subscript[y, 2])"},PlotLabel->"FWR-SIMP",PlotRange->All];
][[1]];
count10=kepler[][[1]];
Show[pic10,MmaPicKepler]

FWR-Trapezoid

kepler[]={0,0};
Time11=AbsoluteTiming[
pic11=ListPlot[FWRTRAP[200,2000,(10^(-8)),keplerIC,keplerTmax,kepler][[All,1;;2]],PlotStyle->{Black,PointSize[Medium]},Mesh->Full,AxesLabel-> {"(Subscript[y, 1])","(Subscript[y, 2])"},PlotLabel->"FWR-TRAP",PlotRange->All];
][[1]];
count11=kepler[][[1]];
Show[pic11,MmaPicKepler]

FWR-Midpoint

kepler[]={0,0};
Time12=AbsoluteTiming[
pic12=ListPlot[FWRMID[200,800,(10^(-8)),keplerIC,keplerTmax,kepler][[All,1;;2]],PlotStyle->{Black,PointSize[Medium]},Mesh->Full,AxesLabel-> {"(Subscript[y, 1])","(Subscript[y, 2])"},PlotLabel->"FWR-MID",PlotRange->All];
][[1]];
count12=kepler[][[1]];
Show[pic12,MmaPicKepler]

Conservation of energy (using Dynamic)

Simp=DWRDriver[keplerTmin,8π,keplerIC,20,10,(10^(-12)),DWRSIMP,kepler];
Trap=DWRDriver[keplerTmin,8π,keplerIC,20,10,(10^(-12)),DWRTRAP,kepler];
Mid=DWRDriver[keplerTmin,8π,keplerIC,20,10,(10^(-12)),DWRMID,kepler];
Energy[method_]:=Module[{r,y11,y22,y33,y44,energy,plot},
r=Dimensions[method][[1]];
y11=method[[All,1]];
y22=method[[All,2]];
y33=method[[All,3]];
y44=method[[All,4]];
energy=Table[(1/2)(((y33[[i]])^2)+((y44[[i]])^2))-(1/(Sqrt[((y11[[i]])^2)+((y22[[i]])^2)])),{i,1,r}]
];
e1=ListLinePlot[Energy[Simp],PlotRange->Full,AxesLabel-> {"Time","Energy"},PlotLabel->"Simpson",PlotRange->All];
e2=ListLinePlot[Energy[Trap],PlotRange->Full,AxesLabel-> {"Time","Energy"},PlotLabel->"Trapezoid",PlotRange->All];
e3=ListLinePlot[Energy[Mid],PlotRange->Full,AxesLabel-> {"Time","Energy"},PlotLabel->"MidPoint",PlotRange->All];
{e1,e2,e3}

Summary and conclusion

Results={Grid[{{"Gravitational"},{"Method","# Func Evals","Timings"},{"DWR-SIMP",count1,Time1},{"DWR-TRAP",count2,Time2},{"DWR-MID",count3,Time3},{"FWR-SIMP",count4,Time4},{"FWR-TRAP",count5,Time5},{"FWR-MID",count6,Time6}},Frame->All],
Grid[{{"Kepler"},{"Method","# Func Evals","Timings"},{"DWR-SIMP",count7,Time7},{"DWR-TRAP",count8,Time8},{"DWR-MID",count9,Time9},{"FWR-SIMP",count10,Time10},{"FWR-TRAP",count11,Time11},{"FWR-MID",count12,Time12}},Frame->All]}