Numerical ODE HW 1

Use a prebuilt ODE solver in any programming envornment matlab, Mathematica, the GNU Scientific Library, MathCad, etc. to recreate figure 120 in the text by numerically solving and ploting a numerical solution of an ODE.

f1[{y1_,y2_,y3_,y4_,y5_,y6_},μ_] :=Module[
{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; 
μVal =1/81.45;
ySol[t_]=y[t]/.NDSolve[{
y'[t]==f1[y[t],μVal],
y[0]==IC
},
y[t],
{t, 0, TMax}
][[1]];
ParametricPlot[ ySol[t],{t,0,TMax},
AxesLabel->{"y1","y2"}]

Repeat question one in a different programming environment.

Calculus texts (in the current MTU text it is section 11.3 ) describe Euler's method for numerically approximating solutions of initial value problems.   For y'=(y^2) with y(0)=1.

Plot an analytical solution on an appropriate domain.

Use Euler's method with stepsize h=0.1 to create a table of approximate solutions on the domain 0≤x≤1. Include the values and a plot comparing the approximate and analytical solutions.

yAnalytic[t_]:=1/(1-t)
f3[y_]:=(y^2)
Euler[y0_,f_,h_,n_]:=Module[
{y=y0,t=0},
Table[
y=y+h f[y];
t=t+h;
{t,y},
{n}]
]
Show[
Plot[yAnalytic[t],{t,0,1}],
ListPlot[Euler[1.0,f3,0.1,10]]
]

Use a number of different step sizes to examine the dependence of the Euler approximation to y(0.5) on the step size. Include a plot showing the dependence of the error on the stepsize.

ConvDataEuler=Table[
n=(2^p);
h=(0.5/n);
{t,yHalf}=Euler[1.0,f3,h,n][[-1]];
{h,Abs[yHalf-2]},
{p,4, 14}];
Show[
LogLogPlot[x,{x,(10^(-12)),1}],
ListLogLogPlot[ConvDataEuler],
AxesLabel->{"h","err"}
]

Note:  The offset in the convergence data is because it is C h and C > 1. The best estinmate for the power is 0.987 and the best estimate for the constant is 2.48.  See below.  I I use larger steps I will not get a power close to one:  the erorr results we have are as h->0.  If we use step sizes that are too small we may also get funny answers.

Clear[x]
g[x_]=E^Fit[Log[ConvDataEuler],{1,x},x];
g[Log[h]]