Numerical Methods

Calculating Numerical Solutions Using the Euler Method

Find approximate values of the solution of the initial value problem
    y'[t] == t^2 - y[t]^2, y[0] == 1

at t == 1
    1) using the Euler method with h = 0.1
    2) using the Euler method with h = 0.05
Compare the results with NDSolve[] by superimposing the plots of the data.

Solution

EMStep[{t_,y_}] := N[{t+h,y+h*f[t,y]}];

f[t_,y_] = t^2 - y^2;
h = 0.1;

EMStep[{0,1}]

EMStep[%]

EMStep[%]

soln1=NestList[EMStep, {0,1}, 10]

p1 = ListPlot[soln1,PlotRange->{0,1},
     PlotJoined->True, PlotStyle->{RGBColor[0,1,0]}]

h = 0.05;

soln2 = NestList[EMStep,{0,1},20]

p2 = ListPlot[soln2,PlotRange->{0,1},
     PlotJoined->True, PlotStyle->{RGBColor[1,0,0]}]

soln = NDSolve[{y'[t] == t^2 - y[t]^2, y[0] == 1},
       y[t],{t,0,1}]

p = Plot[y[t]/.First[soln],{t,0,1},
PlotRange->{0,1}]

Show[p,p1,p2]

y[t_] = y[t] /. soln;

y[1]

Last[soln1]

Last[soln2]

Calculating Numerical Solutions Using the Euler Method

Find an approximate value of y[2] given that
    y'[t] == t^(1/2) + y[t]^(3/2), y[1] == 3
1) using the Euler method with h = 0.1
2) using the Euler method with h = 0.05
Compare the results with NDSolve[] by superimposing the plots of the data.

Solution

EMStep[{t_,y_}] := N[{t+h,y+h*f[t,y]}];

f[t_,y_] = t^(1/2) + y^(3/2);

h = 0.1;

EMStep[{1,3}]

EMStep[%]

EMStep[%]

soln1 = NestList[EMStep,{1,3},10]

p1 = ListPlot[soln1,PlotRange->{0,100},
     PlotJoined->True, PlotStyle->{RGBColor[0,1,0]}]

h = 0.05;

soln2 = NestList[EMStep,{1,3},20]

p2 = ListPlot[soln2,PlotRange->{0,100},
     PlotJoined->True, PlotStyle->{RGBColor[1,0,0]}]

soln = NDSolve[{y'[t] == t^(1/2)+y[t]^(3/2),y[1]==3},
       y[t],{t,1,2}]

y[t] /. soln /. t->2

p = Plot[y[t]/.First[soln],{t,1,2},
    PlotRange->{0,100}]

Show[p,p1,p2]

h = 0.025;
1/h

soln3 = NestList[EMStep,{1,3},40]

p3 = ListPlot[soln3,PlotRange->{0,200},
     PlotJoined->True, PlotStyle->{RGBColor[0,0,1]}]

Show[p,p2,p3,PlotRange->{0,400}]

Calculating Numerical Solutions Using the Runge-Kutta Method

Find an approximate value of y[2] given that
    y'[t] == t^(1/2) + y[t]^(3/2), y[1] == 3
1) using the Runge-Kutta method with h = 0.1
2) using the Runge-Kutta method with h = 0.05
Compre the results with NDSolve[] by superimposing the plots of the data.

Solution

The Runge-Kutta method is based on ideas similar to Simpson's rule for approximating integrals.
The procedure is basically the same as the Euler method except the formula for predicting the next
point is based on a special average of tangent slopes:
    t1 = t0 + h;
    y1 = y0 + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
where
    k1 = f[t0,y0];
    k2 = f[t0+h/2,y0+h*k1/2];
    k3 = f[t0+h/2,y0+h*k2/2];
    k2 = f[t0+h,y0+h*k3];

RKStep[{t_,y_}] :=
( k1 = N[f[t,y]];
  k2 = N[f[t+0.5*h,y+0.5*h*k1]];
  k3 = N[f[t+0.5*h,y+0.5*h*k2]];
  k4 = N[f[t+h,y+h*k3]];
  N[{t+h,y+(h/6)*(k1+2*k2+2*k3+k4)}] )

f[t_,y_] = t^(1/2) + y^(3/2);

h = 0.1;

RKStep[{1,3}]

RKStep[%]

RKStep[%]

soln1 = NestList[RKStep, {1,3}, 10]

p1 = ListPlot[soln1,PlotRange->{0,400},
     PlotJoined->True, PlotStyle->{RGBColor[0,1,0]}]

h = 0.05;

soln2 = NestList[RKStep, {1,3}, 20]

p2 = ListPlot[soln2,PlotRange->{0,400},
     PlotJoined->True, PlotStyle->{RGBColor[1,0,0]}]

Show[p1,p2,p]

h = 0.025;

soln3 = NestList[RKStep, {1,3}, 40]

p3 = ListPlot[soln3,PlotRange->{0,400},
     PlotJoined->True, PlotStyle->{RGBColor[0,0,1]}]

Show[p3,p]

soln = NDSolve[{y'[t] == t^(1/2)+y[t]^(3/2),y[1]==3},
       y[t],{t,1,3}]

Plot[y[t]/.%, {t,1,2.1}]

DSolve[{y'[t] == y[t]^(3/2),y[1]==3},y[t],t]

N[1+2/Sqrt[3]]