Numerical Solution of Initial Value Problem

Our goal is to solve the initial value problem 
                            x' = f(t,x),    x(a) = (Subscript[x, 0]).
Suppose we want to solve the equation on the interval [a, b].  The solution is then a function x = x(t) on [a, b]. Instead of finding the function x(t) analytically, we find numerical values of x(t) on a collection of points (Subscript[t, i]) = a + i h, i = 0, 1, ..., n,  and compute iteratively the values (Subscript[x, i]), i = 0, 1,.... n, which approximate x((Subscript[t, i])),  i = 0, 1,.... n.

Euler's Method

For solving the initial value problem on [a, b], we choose  h = ((b-a)/n). The formula of Euler's method reads       
              (Subscript[t, k]) = a + k h,     (Subscript[x, k+1])= (Subscript[x, k])+f((Subscript[t, k]), (Subscript[x, k])),    k = 0, 1, ..., n,   
 where (Subscript[t, 0]) =a and (Subscript[x, 0]) is the given initial value.

Algorithm

Euler[f0_,a_,b_,x0_,n_] :=
Module[{f,h,k,T, X }, Set@@{f[t_,x_],f0};
  h = (b - a)/ n //N ;
  T = Table[0,{n+1}];
  X = Table[0,{n+1}];
  T[[1]] =  a //N;  X[[1]] = x0 //N;
  For[k=1, k<=n, k++, 
    X[[k+1]] = X[[k]] + h f[T[[k]],X[[k]]] //N; T[[k+1]] = a + h k //N];
 Return[ Transpose[{T,X}] ] ];

The function of the algorithm is of the form:  Euler[function, a, b, (Subscript[x, 0]), n], where n is the number of steps, so that h = ((b-a)/n).

Example

Example 1. Find the solution of x' = cos t - sin x + (t^2) and x(-1) =3 on [-1,1].

f[t_,x_] = Cos[t] - Sin[x] + t^2;
Euler[f[t,x],-1,1,3,10]

pts = Euler[f[t,x],-1,1,3,10];
  dots= ListPlot[pts,PlotStyle-> {PointSize[0.015]}, DisplayFunction->  Identity];
  lines= ListLinePlot[pts, AspectRatio -> 0.5 , DisplayFunction->  Identity ]; 
  Show[ lines, dots, DisplayFunction-> $DisplayFunction]

Here is the computation with 200 points:

pts = Euler[f[t,x],-1,1,3,200];
  dots= ListPlot[pts,PlotStyle-> {PointSize[0.005]}, DisplayFunction->  Identity];
  lines= ListLinePlot[pts,  AspectRatio -> 0.5 , DisplayFunction->  Identity ]; 
  Show[ lines, dots, DisplayFunction-> $DisplayFunction]

Compare with Mathematica solution:

ss=NDSolve[{x'[t]==f[t,x[t]], x[-1]==3},x[t],{t,-1,1}];
Plot[Evaluate[x[t]/.ss],{t,-1,1},AspectRatio -> 0.5];

According to Mathematica description, this solution is good up to MachinePrecision -10 = 6 digits.  Let us look at the error on the first 10 evaluation points, while  h = 2/200 = 0.01:

xx[t_]= x[t] /. First[ss]; 
Error = Table[pts[[i,2]]- xx[-1+ (i-1) /100], {i,1,10}]

This solution is a very flat curve and h is small, the accumulation of errors in the first 10 step is not very prominent.

Example 2. Find the solution of x' =(t x -  (x^2))/(t^2) and x(1) =2 on [1, 3]. The analytic solution of the this initial value problem is given by x(t) = ((2 t)/(1+2 log t)).

xx[t_]:=(2 t)/(1 + 2 Log[t]);

Clear[f,t,x];
f[t_,x_] = (t*x - x^2) / t^2;
Euler[f[t,x],1,3,2,20];

pts = Euler[f[t,x],1,3,2,20]; 
solu= Plot[xx[t], {t,1,3},DisplayFunction->Identity];
dots= ListPlot[pts,PlotStyle-> {PointSize[0.01]}, DisplayFunction->  Identity];
Show[  dots, solu, DisplayFunction-> $DisplayFunction]

Error = Table[pts[[i,2]]- xx[1+(i-1)/10], {i,1,21}]

pts = Euler[f[t,x],1,3,2,200];
dots= ListPlot[pts,PlotStyle-> {PointSize[0.01]}, DisplayFunction->  Identity];
Show[realsolu,  dots, DisplayFunction-> $DisplayFunction]

Let us look at the error on the first 11 points being evaluated, where h = 2/100 = 0.05.

xx[t_]=(2 t)/(1 + 2 Log[t]);
Error1 = Table[pts[[i,2]]- xx[1+ (i-1) /100], {i,1,10}]//N

We know the local errors accumulating from step to step, so that the error should get worse as k get large. Let us look at the error on the last 11 points being evaluated.

Error2=Table[pts[[i,2]]- xx[1+ (i-1) /100], {i,90,100}]//N

The accumulation of the errors is clearly visible in this case.

Heun's Method

For solving the initial value problem on [a, b], we choose  h = ((b-a)/n). The formula of Heun's method reads       
                    (Subscript[t, k])= a + k h,     (Subscript[x, k+1])= (Subscript[x, k])+(h/2)  [f((Subscript[t, k]), (Subscript[x, k]))+ f((Subscript[t, k+1]), (Subscript[t, k])+h (Subscript[x, k]))],   k = 0, 1, ..., n,   
 where (Subscript[t, 0]) =a and (Subscript[x, 0]) is the given initial value.

Algorithm

Heun[f0_,a_,b_,x0_,n_] := Module[{f,h,j},Set@@{f[t_,x_],f0};
  h = (b - a)/n //N;
  T = Table[0,{n+1}];
  X = Table[0,{n+1}];
  T[[1]] = a //N;  X[[1]] = x0 //N;
  For[j=1, j<=n,j++,
    F1 = f[T[[j]],X[[j]]] //N; T[[j+1]] = a + h j //N;
    F2 = f[T[[j+1]],X[[j]] + h F1 ] //N;
    X[[j+1]] = X[[j]] + h (F1 + F2)/2 //N];
  Return[Transpose[{T,X}]] ];

The function of the algorithm is of the form:  Heun[function, a, b, (Subscript[x, 0]), n], where n is the number of steps, so that h = ((b-a)/n).

Example

Example 2. Find the solution of x' =(t x -  (x^2))/(t^2) and x(1) =2 on [1, 3]. The analytic solution of the this initial value problem is given by x(t) = ((2 t)/(1+2 log t)).

Clear[f,t,x];
f[t_,x_] = (t *x - x^2) / t^2;
Heun[f[t,x],1,3,2,200];

Clear[pts,dots]; pts = Heun[f[t,x],1,3,2,20];
dots= ListPlot[pts,PlotStyle-> {PointSize[0.01]}, DisplayFunction->  Identity];
Show[  dots,solu, DisplayFunction-> $DisplayFunction]

Clear[pts,dots]; pts = Heun[f[t,x],1,3,2,200];
dots= ListPlot[pts,PlotStyle-> {PointSize[0.01]}, DisplayFunction->  Identity];
Show[  dots,solu, DisplayFunction-> $DisplayFunction]

Let us look at the error on the first 11 points being evaluated, where h = 2/100 = 0.05.

xx[t_]=(2 t)/(1 + 2 Log[t]);
Error1 = Table[pts[[i,2]]- xx[1+ (i-1) /100], {i,1,10}]//N

We know the local errors accumulating from step to step, so that the error should get worse as k get large. Let us look at the error on the last 11 points being evaluated.

Error2=Table[pts[[i,2]]- xx[1+ (i-1) /100], {i,90,100}]//N

Runge-Kutta's Method of Order 4

For solving the initial value problem on [a, b], we choose  h = ((b-a)/n). The formula of Heun's method reads       
       (Subscript[t, k]) = a + k h,     (Subscript[x, k+1])= (Subscript[x, k])+(h/6)  [(Subscript[F, 1]) + 2 (Subscript[F, 2 ]) + 2 (Subscript[F, 3 ])+ (Subscript[F, 4 ])],     k = 0, 1, ..., n,    where 
            (Subscript[F, 1]) = f((Subscript[t, k]), (Subscript[x, k])),  (Subscript[F, 2]) = f((Subscript[t, k])+(h/2), (Subscript[t, k])+h ((Subscript[F, 1])/2)),    (Subscript[F, 3]) = f((Subscript[t, k])+(h/2), (Subscript[t, k])+h ((Subscript[F, 2])/2)),    (Subscript[F, 4]) = f((Subscript[t, k])+h, (Subscript[t, k])+h (Subscript[F, 3])),
 with (Subscript[t, 0]) =a and (Subscript[x, 0]) being the given initial value.

Algorithm

RungeKutta[f0_,a_,b_,x0_,n_] :=
Module[{f,h,j,F1,F2,F3,F4,T,X},
  Set@@{f[t_,x_],f0};
  h = (b - a)/n //N;
  T = Table[0,{n+1}];
  X = Table[0,{n+1}];
  T[[1]] = a //N;
  X[[1]] = x0 //N;
  For[j=1, j<=n,j++,
    F1 = h f[T[[j]],X[[j]]] //N;
    F2=h f[T[[j]]+h/2,X[[j]]+F1/2]//N;
    F3=h f[T[[j]]+h/2,X[[j]]+F2/2]//N;
    F4=h f[T[[j]]+ h, X[[j]]+F3] //N;
    X[[j+1]] = X[[j]] + (F1 + 2 F2 + 2 F3 + F4)/6 //N;
    T[[j+1]] = a + h j //N];
 Return[Transpose[{T,X}]] ];

The function of the algorithm is of the form:  RungeKutta[function, a, b, (Subscript[x, 0]), n], where n is the number of steps, so that h = ((b-a)/n).

Example

Example 2. Find the solution of x' =(t x -  (x^2))/(t^2) and x(1) =2 on [1, 3]. The analytic solution of the this initial value problem is given by x(t) = ((2 t)/(1+2 log t)).

Clear[f,t,x];
f[t_,x_] = (t *x - x^2) / t^2;

Clear[pts,dots]; pts = RungeKutta[f[t,x],1,3,2,20];
dots= ListPlot[pts,PlotStyle-> {PointSize[0.01]}, DisplayFunction->  Identity];
Show[solu,  dots, DisplayFunction-> $DisplayFunction]

Clear[pts,dots]; pts = RungeKutta[f[t,x],1,3,2,200];
dots= ListPlot[pts,PlotStyle-> {PointSize[0.01]}, DisplayFunction->  Identity];
Show[solu,  dots, DisplayFunction-> $DisplayFunction]

Let us look at the error on the first 11 points being evaluated, where h = 2/100 = 0.05.

xx[t_]=(2 t)/(1 + 2 Log[t]);
Error1 = Table[pts[[i,2]]- xx[1+ (i-1) /100], {i,1,10}]//N

We know the local errors accumulating from step to step, so that the error should get worse as k get large. Let us look at the error on the last 11 points being evaluated.