Numerical Method, Scalar Equations
by R Clark Robinson

ClearSystemCache;

This notebook applis the three numerical methods, Euler, Heun, and Runge-Kutta, to a scalar differential equation. The function that is the right-hand side of the differential equation is

f := Function[{t,x}, x*(1-x^2)]

The step size is

h := 0.1

The initial conditions are

T[0] =0 ;
x0 = 2.0;

Euler Method

The slope at the initial point {t = 0, x=2} is

T[0] = 0;xe[0] = 2 ;
f[T[0],xe[0]]

After one step the values of  t, x, and the slope are approximately

T[1] = T[0] + h
xe[1] = xe[0] + h*f[T[0],xe[0]]
f[T[1],xe[1]]

Ater two steps the values of t,x, and the slope are approximately

T[2] = T[1] + h
xe[2] = xe[1] + h*f[T[1],xe[1]]
f[T[2],xe[2]]

Repeating the process for steps 3, 4, and 5, we get

T[3] = T[2] + h
xe[3] = xe[2] + h*f[T[2],xe[2]]
f[T[3],xe[3]]

T[4] = T[3] + h
xe[4] = xe[3] + h*f[T[3],xe[3]]
f[T[4],xe[4]]

T[5] = T[4] + h
xe[5] = xe[4] + h*f[T[4],xe[4]]
f[T[5],xe[5]]

Combining the steps into a table and calculating 10 terms up to t = 10

h= 0.1; Num := 10; xe[0] = 2 ;T[0] = 0;
f := Function[{t,x}, x*(1-x^2)]
k[0] = f[T[0],xe[0]];
Print["       Euler Method"]
Print[" t_i    x_i         k_i"]
Print[
PaddedForm[T[0],{2,1}],
"    ",
PaddedForm[xe[0],{6,5},NumberPadding -> {"","0"}],
"    ",
PaddedForm[k[0],{6,5},NumberPadding -> {"","0"}]]
Do[  
T[i]   =  T[i-1] +0.1 ;
 xe[i] = xe[i-1]+ h*f[T[i-1],xe[i-1]];
k[i] =  f[T[i],xe[i]];
Print[
PaddedForm[T[i],{2,1}],
"    ",
PaddedForm[xe[i],{6,5},NumberPadding -> {"","0"}],
"    ",
PaddedForm[k[i],{6,5},NumberPadding -> {"","0"}]], 
{i,1,Num}]

To understand the approximation, we plot the Euler solution with the direction graph and also with the “true” solution (actually also numerical)

OdeScalar := {x'[t] == x[t]*(1-x[t]^2)}
xmin := 1.0
xmax := 2.0
Tmin := 0
Tmax := 1
dirgraph=VectorPlot[{1,f[t,x]},{t,Tmin,Tmax},{x,xmin,xmax},VectorScale  -> {0.05,0.8}];
SolnOde := NDSolve[ {x'[t] == x[t]*(1-x[t]^2), x[0]== 2.0},x,{t,Tmin,Tmax}]
phaseplot = Plot[Evaluate[x[t]/.SolnOde],{t,Tmin,Tmax},
PlotRange -> {xmin,xmax}, PlotStyle -> Blue ];
eulerplot = ListPlot[Table[{T[i],xe[i]},{i,0,Num}], Joined -> True,PlotStyle -> Red];
Show[dirgraph,eulerplot,phaseplot]

The red plot is the Euler approximation,  the blue curve is the true solution.  
Observe how the computed solution closely follows the direction field, but the Euler solution is not close.

Heun or Improved Euler

First we calculate the table of the first 10 steps, with the two slopes displayed

h= 0.1; Num := 10; xh[0] = 2 ;T[0] = 0;
f := Function[{t,x}, x*(1-x^2)]
Print["       Heun Method"]
Print[" t_i    x_i         k1_i         k2_i"]
Do[  
T[i+1] = T[i] + h;
kh1[i] =  f[T[i],xh[i]];
zh[i] = xh[i] + h*kh1[i];
kh2[i] = f[T[i+1],zh[i]];
xh[i+1] = xh[i] + 0.5*h*(kh1[i]+kh2[i]);
Print[
PaddedForm[T[i],{2,1}],
"    ",
PaddedForm[xh[i],{6,5},NumberPadding -> {"","0"}],
"    ",
PaddedForm[kh1[i],{6,5},NumberPadding -> {"","0"}], 
"    ",
PaddedForm[kh2[i],{6,5},NumberPadding -> {"","0"}]],
{i,0,Num}]

We plot the Heun numerical solution in magenta and the true solution in blue

heunplot = ListPlot[Table[{T[i],xh[i]},{i,0,Num}], Joined -> True,PlotStyle -> Magenta];
Show[dirgraph,heunplot,phaseplot]

Runge-Kutta

First we calculate the table of the first 10 steps

h= 0.1; Num := 10; xrk[0] = 2 ;T[0] = 0;
f := Function[{t,x}, x*(1-x^2)]
h2 = 0.5*h;
Print["  Runge-Kutta Method"]
Print[" t_i    x_i"]
Do[  
T[i+1] = T[i] + h;
krk1[i] =  f[T[i],xrk[i]];
zrk1[i] = xrk[i] + h2*krk1[i];
krk2[i] = f[T[i]+h2,zrk1[i]];
zrk2[i] = xrk[i] + h2*krk2[i];
krk3[i] = f[T[i]+h2,zrk2[i]];
zrk3[i] = xrk[i] + h*krk3[i];
krk4[i] = f[T[i+1],zrk3[i]];
xrk[i+1] = xrk[i] + h*(krk1[i]+2*krk2[i]+2*krk3[i]+krk4[i])/6.0;
Print[
PaddedForm[T[i],{2,1}],
"    ",
PaddedForm[xrk[i],{6,5},NumberPadding -> {"","0"}]],
{i,0,Num}]

Show the plots of the Runge-Kutta approximation in green and true solution in blue

RungeKuttaplot = ListPlot[Table[{T[i],xrk[i]},{i,0,Num}], Joined -> True,PlotStyle -> Green];
Show[dirgraph,RungeKuttaplot,phaseplot]

Compare the three numerical solutions and the “true” solution:  Euler is in red, Heun is in Magenta, Runge-Kutta is in green, and the “true” solution is in blue

Show the plots of the three numerical methods: Euler in red, Heun in magenta, Runge-Kutta in green, and “true” solution in blue.

Show[{eulerplot,heunplot,RungeKuttaplot,phaseplot},PlotRange -> {{Tmin,Tmax},{xmin,xmax}}]

Compare the numerical values of the three methods

The exact solution is  x(t) = sqrt( 4/(4-3e^{-2t}) )

N[e

Print["\  t      Euler       Heun   Runge-Kutta   Exact"]
Do[
Print[
PaddedForm[T[i],{2,1}],
"    ",
PaddedForm[xe[i],{6,5},NumberPadding -> {"","0"}],
"    ",
PaddedForm[xh[i],{6,5},NumberPadding -> {"","0"}],
"    ",
PaddedForm[xrk[i],{6,5},NumberPadding -> {"","0"}],
"    ",
PaddedForm[Sqrt[4/(4-3 E^(-2 T[i]))],{6,5},NumberPadding -> {"","0"}]],
{i,0,Num}]