Numerical Method for Systems of Differential Equations
by R Clark Robinson

ClearSystemCache;

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

fx := Function[{x,y}, y]
fy := Function[{x,y}, -4 x]

The step size is

h := 0.1

The initial conditions are

T[0] =0 ;
x0 = 0.0;
y0= 1.0;

Euler Method

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

T[0] = 0;xe[0] = x0 ; ye[0] = y0;
fx[xe[0],ye[0]]
fy[xe[0],ye[0]]

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

T[1] = T[0] + h
xe[1] = xe[0] + h*fx[xe[0],ye[0]]
ye[1] = ye[0] + h*fy[xe[0],ye[0]]

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

T[2] = T[1] + h
xe[2] = xe[1] + h*fx[xe[1],ye[1]]
ye[2] = ye[1] + h*fy[xe[1],ye[1]]

T[3] = T[2] + h
xe[3] = xe[2] + h*fx[xe[2],ye[2]]
ye[3] = ye[2] + h*fy[xe[2],ye[2]]

T[4] = T[3] + h
xe[4] = xe[3] + h*fx[xe[3],ye[3]]
ye[4] = ye[3] + h*fy[xe[3],ye[3]]

T[5] = T[4] + h
xe[5] = xe[4] + h*fx[xe[4],ye[4]]
ye[5] = ye[4] + h*fy[xe[4],ye[4]]

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

h= 0.1; Num := 10;T[0] = 0;
 
xe[0] = x0;
ye[0] = y0;
Print["       Euler Method"]
Print[" t_i    x_i    y_i        kx_i       ky_i"]
Do[  
kx[i] =  fx[xe[i],ye[i]];
ky[i] =  fy[xe[i],ye[i]];
T[i+1]   =  T[i] + h ;
xe[i+1] = xe[i]+ h*kx[i];
ye[i+1] = ye[i] + h*ky[i];
Print[
PaddedForm[T[i],{2,1}],
"    ",
PaddedForm[xe[i],{6,5},NumberPadding -> {"","0"}],
"    ",
PaddedForm[ye[i],{6,5},NumberPadding -> {"","0"}],
"    ",
PaddedForm[kx[i],{6,5},NumberPadding -> {"","0"}],
"    ",
PaddedForm[ky[i],{6,5},NumberPadding -> {"","0"}]],
{i,0,Num}]

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

OdeSystem := {x'[t] == y[t], y'[t] == -4 x[t]}
xmin := 0.0;
xmax := 1.0;
ymin := -1.0;
ymax := 1.0;
Tmin := 0
Tmax := 1
OdeVectorField := {y,-4 x}
dirgraph=VectorPlot[{fx[x,y],fy[x,y]},{x,xmin,xmax},{y,ymin,ymax},VectorScale  -> {0.07,0.5}];
SolnOde := NDSolve[ {OdeSystem[[1]], OdeSystem[[2]], x[0]== x0, y[0]==y0},{x,y},{t,Tmin,Tmax}]
phaseplot = ParametricPlot[Evaluate[{x[t],y[t]}/.SolnOde],{t,Tmin,Tmax},
PlotRange -> {{xmin,xmax},{ymin,ymax}}, PlotStyle -> Blue];
eulerplot = ListPlot[Table[{xe[i],ye[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 vectors displayed

h= 0.1; Num := 10; xh[0] = x0;yh[0]=y0;T[0] = 0;
Print["       Heun Method"]
Print[" t_i    x_i         y_i"]
Do[  
k1hx[i] =  fx[xh[i],yh[i]];
k1hy[i] =  fy[xh[i],yh[i]];
zhx[i] = xh[i] + h*k1hx[i];
zhy[i] = yh[i] + h*k1hy[i];
k2hx[i] = fx[zhx[i],zhy[i]];
k2hy[i] = fy[zhx[i],zhy[i]];
T[i+1] = T[i] + h;
xh[i+1] = xh[i] + 0.5*h*(k1hx[i]+k2hx[i]);
yh[i+1] = yh[i] + 0.5*h*(k1hy[i]+k2hy[i]);
Print[
PaddedForm[T[i],{2,1}],
"    ",
PaddedForm[xh[i],{6,5},NumberPadding -> {"","0"}],
"    ",
PaddedForm[yh[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[{xh[i],yh[i]},{i,0,Num}], Joined -> True,PlotStyle -> Magenta];
Show[dirgraph,heunplot]

Runge-Kutta

First we calculate the table of the first 10 steps

h= 0.1; Num := 10;T[0] = 0; xrk[0] = x0;yrk[0] = y0;
h2 = 0.5*h;
Print["  Runge-Kutta Method"]
Print[" t_i    x_i   y_i"]
Do[  
k1rkx[i] =  fx[xrk[i],yrk[i]];
k1rky[i] =  fy[xrk[i],yrk[i]];
z1rkx[i] = xrk[i] + h2*k1rkx[i];
z1rky[i] = yrk[i] + h2*k1rky[i];
k2rkx[i] = fx[z1rkx[i],z1rky[i]];
k2rky[i] = fy[z1rkx[i],z1rky[i]];
z2rkx[i] = xrk[i] + h2*k2rkx[i];
z2rky[i] = yrk[i] + h2*k2rky[i];
k3rkx[i] = fx[z2rkx[i],z2rky[i]];
k3rky[i] = fy[z2rkx[i],z2rky[i]];
z3rkx[i] = xrk[i] + h*k3rkx[i];
z3rky[i] = yrk[i] + h*k3rky[i];
k4rkx[i] = fx[z3rkx[i],z3rky[i]];
k4rky[i] = fy[z3rkx[i],z3rky[i]];
T[i+1] = T[i] + h;
xrk[i+1] = xrk[i] + h*(k1rkx[i]+2*k2rkx[i]+2*k3rkx[i]+k4rkx[i])/6.0;
yrk[i+1] = yrk[i] + h*(k1rky[i]+2*k2rky[i]+2*k3rky[i]+k4rky[i])/6.0;
Print[
PaddedForm[T[i],{2,1}],
"    ",
PaddedForm[xrk[i],{6,5},NumberPadding -> {"","0"}],
"    ",
PaddedForm[yrk[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[{xrk[i],yrk[i]},{i,0,Num}], Joined -> True,PlotStyle -> Green];
Show[dirgraph,RungeKuttaplot]

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[{eulerplot,heunplot,RungeKuttaplot,phaseplot},PlotRange -> {{Tmin,Tmax},{xmin,xmax}}]

Compare the numerical values

The solution in closed form is x(t) = sin(2t), y(t) = cos(2t).  We print the y values of the methods

Print["\  t    Euler    Heun   Runge-Kutta   Exact"]
Do[
Print[
PaddedForm[T[i],{2,1}],
"    ",
PaddedForm[ye[i],{6,5},NumberPadding -> {"","0"}],
"    ",
PaddedForm[yh[i],{6,5},NumberPadding -> {"","0"}],
"    ",
PaddedForm[yrk[i],{6,5},NumberPadding -> {"","0"}],
"    ",
PaddedForm[Cos[2 T[i]],{6,5},NumberPadding -> {"","0"}]],
{i,0,Num}]