Numerical Methods: Euler and Runge-Kutta

Abstract

Most real life phenomena change with time, hence dynamic. Differential equations are used in mathematical modeling of such scenarios. Linear differential equations can be solved analytically but most real life applications are nonlinear. Numerical solutions of nonlinear differential equations are approximate solutions. Euler and Runge-Kutta method of order four are derived, explained and illustrated as useful numerical methods for solving single and systems of linear and nonlinear differential equations. Accuracy of a numerical method depends on the step size used and degree of nonlinearity of the equations. Stiffness is another challenge with numerical solutions of nonlinear differential equations. Although better accuracy can be obtained with smaller step size, this takes more computational effort and time. Algorithms and codes can be written using available computer programming software to overcome this challenge and to avoid computational error. The Runge-Kutta method is more applicable and accurate for diverse classes of differential equations.

1. Introduction

In numerical analysis, mathematical methods are employed to produce numerical answers to mathematical expressions. It entails developing, examining, and putting into practice computer algorithms to solve continuous mathematics problems numerically [1]. These problems may originate from real-life applications in the natural sciences, social sciences, engineering, medical sciences and business where variables which vary continuously are involved.

Differential equations are the foundation of the majority of mathematical models used in the natural sciences and engineering. The numerical technique finds the solution function at a discrete set of points, approximating the derivatives or integrals in the relevant equation.

In an effort to obtain solutions that are more precise, more affordable, or more resilient to instabilities in the issue data, a number of techniques have been developed. Many of the differential equations that arise in practical problems are usually nonlinear. Nonlinear differential equations are nearly impossible to be solved precisely [2]. Most differential equations that arise in practical problems are typically nonlinear. The computational labour of solving a system of first order equations may be formidable even when the exact solution is possible [3]. It is for these reasons that techniques have been developed for computing to any degree of accuracy the numerical solution of almost any of such problems.

These techniques typically fall into “families” of increasing order of accuracy, with Euler’s method (or a close relative) frequently making up the lowest order. Adams-Bashforth techniques are “multistep” methods, whereas Runge-Kutta methods are “single-step” methods. These methods have been very successful and reliable.

In this chapter, the derivation of Euler and Runge Kutta methods are explained and illustrative examples given to explain how they are used. Algorithms/codes of the examples written with Maple mathematical programming software were also included for better comprehension and further practice.

2. The Euler method

The popular Euler method published in 1768 is attributed to Leonhard Euler (1707–1783). The method is one of the most straightforward and fundamental method of solving single and systems of ordinary differential eqs. [1]. The basic idea is as follows;

Consider initial value system:

$$\frac{\mathit{dy}}{\mathit{dt}}=f\left(t,x,y\right)y\left(0\right)=y_0$$E1

$$\frac{\mathit{dx}}{\mathit{dt}}=g\left(t,x,y\right)x\left(0\right)=x_0$$E2

By the definition of a derivative,

$$y^{\prime }\left(t\right)=\underset{h\to 0}{lim}\frac{f\left(t+h\right)-f\left(t\right)}{h}$$E3

$$x^{\prime }\left(t\right)=\underset{h\to 0}{lim}\frac{g\left(t+h\right)-f\left(t\right)}{h}$$E4

for small $h>0$, then, Eqs. (3) and (4) imply that a reasonable difference approximation for $y^{\prime }\left(t\right)$ and $x^{\prime }\left(t\right)$ [4] are.

$$y^{\prime }\left(t\right)=\frac{f\left(t+h\right)-f\left(t\right)}{h}$$E5

$$x^{\prime }\left(t\right)=\frac{g\left(t+h\right)-f\left(t\right)}{h}$$E6

Substituting Eq. (5) and Eq. (6) into Eq. (1) and Eq. (2) respectively yield the difference equations

$$\frac{f\left(t+h\right)-f\left(t\right)}{h}=f\left(t,x,y\right)$$E7

$$\frac{g\left(t+h\right)-f\left(t\right)}{h}=g\left(t,x,y\right)$$E8

which approximates the differential Eqs. (1) and (2).

According to [1], the Euler’s method approximates a solution (x, y) with step size h by:

$$y_{k+1}=y_k+\mathit{hf}\left(x_k,y_k\right)$$E9

$$x_{k+1}=x_k+\mathit{hg}\left(x_k,y_k\right)$$E10

Expressing Eq. (1) and Eq. (2) in vector notation:

$$\frac{\mathrm{dY}}{\mathrm{dt}}=F\left(Y\right),Y_0$$E11

Where $Y=\left(y,x\right)$, $Y_0=\left(y_0,x_0\right)$, $F\left(Y\right)=\left(f\left(x,y\right),g(x,y)\right)$ and $\frac{\mathrm{dY}}{\mathrm{dt}}=\left(\frac{\mathit{dy}}{\mathit{dt}},\frac{\mathit{dx}}{\mathit{dt}}\right)$.

The Euler’s method approximates a solution (x, y) by:

$$\left(x_{k+1},y_{k+1}\right)=\left(x_k,y_k\right)+\mathit{hF}\left(x_k,y_k\right)$$E12

Illustration 1

Given the initial value problem in [5]:

$$y^{\prime }=\frac{y-x}{y+x},y\left(0\right)=1.$$

Determine.

$y$ for $x=0.1$, using the Euler method with the step size $h=0.02$.

Solution

The Euler formula is.

$y_{n+1}=y_n+\mathit{hf}\left(x_n,y_n\right)$ where n is the number of iteration.

For the first iteration which is $n=0$, the Euler formula becomes

$$y_1=y_0+\mathit{hf}\left(x_0,y_0\right)$$

$x_0=0$ and $y_0=1$ as the given initial condition.

$$x_1=x_0+h$$

$$x_n=x_0+\mathit{nh}$$

$$n=\frac{x_n-x_0}{h}=\frac{0.1-0}{0.02}=5$$

$$y_1=1+0.02f\left(0,1\right)=1+0.02\left(\frac{1-0}{1+0}\right)=1+0.02=1.02$$

When $n=1$, $x_1=x_0+h=0+0.02=0.02$.

$y_2=y_1+0.02f\left(x_1,y_1\right)=1.02+0.02f\left(0.02,1.02\right)=1.02+$ $0.02\left(\frac{1.02-0.02}{1.02+0.02}\right)=1.039230769$ When $n=2$, $x_2=x_0+2h=0+2\left(0.02\right)=0.04$

$$y_3=y_2+0.02f\left(x_2,y_2\right)=1.039230769+0.02\left(\frac{1.039230769-0.04}{1.039230769+0.04}\right)=1.057748232$$

When $n=3$, $x_3=x_2+h=0.04+0.02=0.06$

$$y_4=y_3+0.02f\left(x_3,y_3\right)=1.057748232+0.02\left(\frac{1.057748232-0.06}{1.057748232+0.06}\right)=1.075601058$$

When $n=4$, $x_4=x_3+h=0.06+0.02=0.08$

$$y_5=y_4+0.02f\left(x_4,y_4\right)=1.075601058+0.02\left(\frac{1.075601058-0.08}{1.075601058+0.08}\right)=1.092831936$$

Table 1 gives the numerical solution of Illustration 1 using different step sizes. Obviously this takes more computational effort since more iteration is needed but gives more accurate result.

X	y(x) at h = 0.02, n = 5	y(x) at h = 0.01, n = 10	y(x) at h = 0.005, n = 20
0	1.0	1.0	1.0
0.02	$1.02000000$	1.019803922	1.019706799
0.04	$1.039230769$	1.038852886	1.038665609
0.06	$1.057748232$	1.057200704	1.056929219
0.08	$1.075601058$	1.074894334	1.074543762
0.10	$1.092831936$	1.091975065	1.091549847

Table 1.

Numerical solution of Example on Euler’s method.

The algorithm/ code given below can be adapted to various single differential equations for the Euler method (Figure 1).

Algorithm/ Codes using Maple programming software

  ## Illustration on Euler method#

 ode:=diff(y(x),x)=(y(x)-x)/(y(x)+x);

sol:dsolve({ode,y(0)=1},y(x));

  ##restart:y:=array(0..200):x:=array(0..200):n:=5.0:A:=0.0:h:=0.02:x[0]:=0.0:y[0]:=1.0:

  for m from 0 to n do x[m]:=A+m*h:y[0]:=1.0:od:

  for m from 0 to n do y[m+1]:=y[m]+h*((y[m]-x[m])/(y[m]+x[m])):od:

  for m from 0 to n do print (m, x[m], y[m]);od:

  0, 0., 1.0

  1, 0.02, 1.020000000

  2, 0.04, 1.039230769

  3, 0.06, 1.057748232

  4, 0.08, 1.075601058

  5, 0.10, 1.092831936

Figure 1.

Graph of numerical solution of Illustration 1 for six iterations.

For more on Maple programming codes and algorithms see [6].

Illustration 2

Consider the Lorenz dynamical system given as.

$$\frac{\mathit{dx}}{\mathit{dt}}=\sigma \left(y\left(t\right)-x\left(t\right)\right)\mathbf{,}x\left(0\right)=0.1$$

$$\frac{\mathit{dy}}{\mathit{dt}}=-x\left(t\right)z\left(t\right)+\mathit{rx}\left(t\right)-y\left(t\right)\mathbf{,}y\left(0\right)=0.1$$

$$\frac{\mathit{dz}}{\mathit{dt}}=x\left(t\right)y\left(t\right)-\mathit{bz}\left(t\right)\mathbf{,}z\left(0\right)=0.1$$

Where $\sigma =10$, $r=28$ and $b=\raisebox{1ex}{$8$}\!\left/ \!\raisebox{-1ex}{$3$}\right.$ with step size 0.01 determine the solution at $t=10$ using the Euler’s method_.

Solution

The Lorenz system is a classical nonlinear dynamical system. Due to the high degree of nonlinearity of the system and the number of iterations required, its numerical solution is only easily tractable using computer programming. The algorithm/Code in Maple mathematical software is given below after the table of numerical solution (Figures 2–4) and (Table 2).

Figure 2.

Graphical representation of numerical solution of Lorenz system X(t) using Euler method for 1000 iterations as in Table 2.

Figure 3.

Graphical representation of numerical solution of Lorenz system Y(t) using Euler method for 1000 iterations as in Table 2.

Figure 4.

Graphical representation of numerical solution of Lorenz system Z(t) using Euler method for 1000 iterations as in Table 2.

$t$	X(t)	Y(t)	Z(t)
1.0	$-$4.277716783	$-$1.835441598	26.77337837
2.0	3.307860445	5.707928748	12.82481891
3.0	$-$9.048257350	$-$3.394242770	33.60229026
4.0	2.906631532	4.606304372	17.27981641
5.0	$-$14.25707158	$-$18.85653324	29.67196074
6.0	0.1457099124	$-$1.419066178	21.71150940
7.0	8.807447039	11.68841199	23.40171689
8.0	4.038757258	0.5104950902	27.19520327
9.0	$-$6.260708216	$-$9.461245716	18.53277577
10.0	$-$16.65763702	$-$18.00805550	37.42332818

Table 2.

Numerical solution of the Lorenz system using Euler method at step size 0.01.

$$t_n=t_0+\mathit{nh}$$

$$n=\frac{t_n-t_0}{h}=\frac{10-0}{0.01}=1000$$

Algorithm/Maple programming Code of the numerical solution of Lorenz system using Euler’s method

## This is the numerical solution of Lorenz system using Euler method##

sys:=diff(x(t),t)=sigma*(y(t)-x(t)),diff(y(t),t)=-x(t)*z(t)+r*x(t)-y(t),

diff(z(t),t)=x(t)*y(t)-b*z(t);

##restart:X:=array(0..2000000):Y:=array(0..2000000):Z:=array

(0..2000000):T:=array(0..2000000):A:=0.0:h:=0.01:N:=1000:sigma:=

10.0:b:=(8/3):r:=28.0: for m from 0 to N do T[m]:=A+m*h:X[0]:=0.1:

Y[0]:=0.1:Z[0]:=0.1:od: for m from 0 to N do X[m+1]:=X[m]+

h*(sigma*(Y[m]-X[m])):Y[m+1]:=Y[m]+h*(-X[m]*Z[m]+r*X[m]-Y[m]):

Z[m+1]:=Z[m]+h*(X[m]*Y[m]-b*Z[m]):od:

for m from 0 by 100 to N do print (m, T[m], X[m],Y[m],Z[m]);od: with(plots):

func1:=listplot([seq(X[m],m=0..N)],style=line,color=black,thickness=1):

display(func1,labels=["t","X(t)"],axes=boxed,caption="Figure1:Euler solution for X(t)");

func2:=listplot([seq(Y[m],m=0..N)],style=line,color=blue,thickness=1):

display(func2,labels=["t","Y(t)"],axes=boxed,caption="Figure2:Euler solution for

Y(t)"); func3:=listplot([seq(Z[m],m=0..N)],style=line,color=red, thickness=1):

display(func3,labels=["t","Z(t)"],axes=boxed,caption="Figure

3:Euler solution for Z(t)");

For more on Lorenz system see [7].

3. The Runge-Kutta method of order four

One of the well-known numerical techniques for solving differential equations is the Runge-Kutta method. The Runge-Kutta method is a family of methods which depend on the order of derivatives involved. The Runge-Kutta method of order four is the classical form of the method in which four values of the derivatives are used for each iteration [1].

3.1 Derivation of Runge-Kutta method of order four

The Taylor series of functions of one variable is:

$$y\left(x\right)=y\left(x_0\right)+y^{\prime }\left(x_0\right)\left(x-x_0\right)+\frac{y^{″}\left(x_0\right)}{2!}{\left(x-x_0\right)}^2+\frac{y^{‴}\left(x_0\right)}{3!}{\left(x-x_0\right)}^3+\cdots$$E13

$$y\left(x\right)=\sum _{n=0}^{\infty }\frac{y^{\left(n\right)}}{n!}{\left(x-x_0\right)}^n$$E14

when $x=x_0+h$, $x-x_0=h$.

Then Eq. (13) becomes.

$$y\left(x\right)=y\left(x_0\right)+h{y}^{\prime }\left(x_0\right)+h^2\frac{y^{″}\left(x_0\right)}{2!}+h^3\frac{y^{‴}\left(x_0\right)}{3!}+\cdots$$E15

Taylor series of functions of two variables is:

$$f\left(x,y\right)=\sum _{i=0}^{\infty }\frac{1}{n!}{\left\{\left(x-x_0\right)\frac{\partial }{\partial x}+\left(y-y_0\right)\frac{\partial }{\partial y}\right\}}^{i}f\left(x,y\right)$$E16

when $x=x_0+\mathit{mh}$, $x-x_0=\mathit{mh}$and $y-y_0=\mathit{nh}$

$$f\left(x,y\right)=\sum _{i=0}^{\infty }\frac{1}{n!}{\left\{\mathit{mh}\frac{\partial }{\partial x}+\mathit{nh}\frac{\partial }{\partial y}\right\}}^{i}f\left(x,y\right)$$E17

$$f\left(x.y\right)=f\left(x_0,y_0\right)+f^{\text{'}}\left(x_0,y_0\right)\left(\mathit{mh}+\mathit{nh}\right)+\frac{f^{\text{'}\text{'}}\left(x_0,y_0\right){\left(\mathit{mh}+\mathit{nh}\right)}^2}{2!}$$

$$+\frac{f^{\text{'}\text{'}\text{'}}\left(x_0,y_0\right){\left(\mathit{mh}+\mathit{nh}\right)}^3}{3!}+\cdots$$E18

$$f\left(x.y\right)=f\left(x_0,y_0\right)+{\mathit{mhf}}^{\text{'}}\left(x_0,y_0\right)+{\mathit{nhf}}^{\text{'}}\left(x_0,y_0\right)+\frac{f^{\text{'}\text{'}}\left(x_0,y_0\right)\left({\left(\mathit{mh}\right)}^2+2\left(\mathit{mh}\right)\left(\mathit{nh}\right)+{\left(\mathit{nh}\right)}^2\right)}{2!}$$

$$+\frac{f^{\text{'}\text{'}\text{'}}\left(x_0,y_0\right)\left({\left(\mathit{mh}\right)}^3+3{\left(\mathit{mh}\right)}^2\left(\mathit{nh}\right)+3\left(\mathit{mh}\right){\left(\mathit{nh}\right)}^2+{\left(\mathit{nh}\right)}^3\right)}{3!}+\cdots$$E19

$$\begin{array}{c}f\left(x.y\right)=f\left(x_0,y_0\right)+h\left[{\mathit{mf}}_x+{\mathit{nf}}_y\right]+\frac{h^2}{2!}\left[m^2{f}_{\mathit{xx}}+2{\mathit{mnf}}_{\mathit{xy}}+n_2{f}_{\mathit{yy}}\right]\\ +\frac{h^3}{3!}\left[m^3{f}_{\mathit{xxx}}+3m^2{\mathit{nf}}_{\mathit{xxy}}+3{\mathit{mn}}^2{f}_{\mathit{xyy}}+n^3{f}_{\mathit{yyy}}\right]+\cdots \end{array}$$E20

where the partial derivatives are all evaluated at the point $\left(x_0,y_0\right)$.

Consider the differential equation

$$\frac{\mathit{dy}}{\mathit{dx}}=y^{\text{'}}=f\left(x,y\right)$$E21

$$\mathit{df}=\frac{\partial f}{\partial x}\mathit{dx}+\frac{\partial f}{\partial y}\mathit{dy}$$E22

$$\frac{d^{2}f}{{\mathit{dx}}^2}=\frac{\mathit{df}}{\mathit{dx}}=f^{\text{'}}=\frac{\partial f}{\mathit{dx}}\frac{\mathit{dx}}{\mathit{dx}}+\frac{\partial f}{\partial y}\frac{\mathit{dy}}{\mathit{dx}}$$E23

$$y^{\text{'}\text{'}}=f^{\text{'}}=\frac{\partial f}{\mathit{dx}}+\frac{\partial f}{\partial y}\frac{\mathit{dy}}{\mathit{dx}}$$E24

$$y^{\text{'}\text{'}}=f^{\text{'}}=\frac{\partial f}{\mathit{dx}}+\frac{\partial f}{\partial y}f$$E25

$$y^{\text{'}\text{'}}=f^{\text{'}}=f_x+{\mathit{ff}}_y$$E26

$$y^{‴}=f^{″}={\left[f_x+{\mathit{ff}}_y\right]}^{\prime }={\left[f_x\right]}^{\prime }+{\left[{\mathit{ff}}_y\right]}^{\prime }$$E27

$${\left(f_x\right)}^{\prime }=f_x+{\mathit{ff}}_{\mathit{xy}}$$E28

$${\left({\mathit{ff}}_y\right)}^{\prime }=f\left({f_y}^{\prime }\right)+f_y{f}^{\prime }$$E29

$$y^{‴}=f^{″}=f_{\mathit{xx}}+2{\mathit{ff}}_{\mathit{xy}}+f^2{f}_{\mathit{yy}}+f_x{f}_y+{{\mathit{ff}}_y}^2$$E30

From Eq. (15).

$$y\left(x\right)-y\left(x_0\right)=h{y}^{\prime }\left(x_0\right)+h^2\frac{y^{″}\left(x_0\right)}{2!}+h^3\frac{y^{‴}\left(x_0\right)}{3!}+\cdots$$E31

Hence

$$y\left(x\right)-y\left(x_0\right)=\mathit{hf}\left(x_0\right)+\frac{h^2}{2!}\left(f_x+{\mathit{ff}}_y\right)+\frac{h^3}{3!}\left(f_{\mathit{xx}}+2{\mathit{ff}}_{\mathit{xy}}+f^2{f}_{\mathit{yy}}+f_x{f}_y+{{\mathit{ff}}_y}^2\right)+\cdots$$E32

The main idea is to select several points so that the Taylor series of expansion Eq. (19), coincide with the terms on the right hand side of Eq. (32).

$$\text{Suppose}{m}_1=\mathit{hf}\left(x_0,y_0\right)$$E33

$$m_2=\mathit{hf}\left(x_0+\mathit{nh},y_0+{\mathit{nm}}_1\right)$$E34

$$m_3=\mathit{hf}\left(x_0+\mathit{ph},y_0+{\mathit{pm}}_2\right)$$E35

$$m_4=\mathit{hf}\left(x_0+\mathit{qh},y_0+{\mathit{qm}}_3\right)$$E36

Using Eq. (32), these values become;

$$m_1=\mathit{hf}$$E37

$$m_2=h\left[f+\mathit{nh}\left(f_x+{\mathit{ff}}_y\right)+\frac{{\left(\mathit{nh}\right)}^2}{2}\left(f_{\mathit{xx}}+2{\mathit{ff}}_{\mathit{xy}}+f^2{f}_{\mathit{yy}}\right)+\cdots \right]$$E38

$$m_3=h\left[f+\mathit{ph}\left(f_x+{\mathit{ff}}_y\right)+\frac{{\left(\mathit{ph}\right)}^2}{2}\left(f_{\mathit{xx}}+2{\mathit{ff}}_{\mathit{xy}}+f^2{f}_{\mathit{yy}}\right)+2\mathit{np}\left(f_x{f}_y+{{\mathit{ff}}_y}^2\right)+\cdots \right]$$E39

$$m_4=h\left[f+\mathit{qhc}+\frac{{\left(\mathit{qh}\right)}^2}{2}\left(f_{\mathit{xx}}+2{\mathit{ff}}_{\mathit{xy}}+f^2{f}_{\mathit{yy}}\right)+2\mathit{pq}\left(f_x{f}_y+{{\mathit{ff}}_y}^2\right)+\cdots \right]$$E40

where all functions are evaluated at the point (x₀, y₀).

Considering an expression of the form

$${\mathit{am}}_1+{\mathit{bm}}_2+{\mathit{cm}}_3+{\mathit{dm}}_4$$E41

Equating the expression to the right hand side of Eq. (37) through Eq. (40) to obtain four equations in seven unknowns.

$$\mathit{\text{Coefficient of}}\mathit{hf}:a+b+c+d=1$$E42

$$\mathit{\text{Coefficient of}}\mathit{hf}\left(f_x{f}_y+{{\mathit{ff}}_y}^2\right)\mathit{bn}+\mathit{cp}+\mathit{dq}=\frac{1}{2}$$E43

$$\text{Coefficient of}\mathit{hf}\left(f_{\mathit{xx}}+2{\mathit{ff}}_{\mathit{xy}}+f^2{f}_{\mathit{yy}}\right):{\mathit{bn}}^2+{\mathit{cp}}^2+{\mathit{dq}}^2=\frac{1}{3}$$E44

$$\mathit{\text{Coefficient of}}\mathit{hf}\left(f_x{f}_y+{{\mathit{ff}}_y}^2\right):\mathit{cnp}+\mathit{dpq}=\frac{1}{6}$$E45

Choosing three unknown quantities arbitrarily since there are four equations in seven unknowns.

Let $n=p=\frac{1}{2}$ and $q=1$ in Eq. (42) through Eq. (45). Then

$$a+b+c+d=1$$E46

$$b+c+2d=1$$E47

$$3b+3c+12d=4$$E48

$$3c+6d=2$$E49

Solving Eq. (46) through Eq. (47) simultaneously produced.

$$a=d=\frac{1}{6},b=c=\frac{1}{3}$$

Connecting Eq. (33) and the expression (41).

$$y\left(x\right)-y\left(x_0\right)={\mathit{am}}_1+{\mathit{bm}}_2+{\mathit{cm}}_3+{\mathit{dm}}_4$$E50

$$y\left(x\right)-y\left(x_0\right)=\frac{m_1}{6}+\frac{m_2}{3}+\frac{m_3}{3}+\frac{m_4}{6}$$E51

$$y\left(x\right)=y\left(x_0\right)+\frac{1}{6}\left(m_1+2m_2+2m_3+m_4\right)$$E52

In general

$$y_{n+1}=y_n+\frac{1}{6}\left(m_1+2m_2+2m_3+m_4\right)$$E53

Where

$$m_1=\mathit{hf}\left(x_n,y_n\right)$$

$$m_2=\mathit{hf}\left(x_n+\frac{1}{2}h,y_n+\frac{1}{2}{m}_1\right)$$

$$m_3=\mathit{hf}\left(x_n+\frac{1}{2}h,y_n+\frac{1}{2}{m}_2\right)$$

$$m_4=\mathit{hf}\left(x_n+h,y_n+m_3\right)$$

This is the Runge-Kutta formula of order four for a single first order differential equation.

The formula can be generalized for system of differential equations.

Consider the system of two initial-value differential equation:

$$\frac{\mathit{dx}}{\mathit{dt}}=f\left(t,x,y\right)x\left(0\right)=y_0$$E54

$$\frac{\mathit{dy}}{\mathit{dt}}=g\left(t,x,y\right)y\left(0\right)=x_0$$E55

The Runge-Kutta formula of order four for system of two initial-value differential equation of the form of Eq. (54) and Eq. (55) is

$$x_{n+1}=x_n+\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)$$E56

$$y_{n+1}=y_n+\frac{1}{6}\left(m_1+2m_2+2m_3+m_4\right)$$E57

Where

$$k_1=\mathit{hf}\left(t_n,x_n,y_n\right)$$

$$m_1=\mathit{hg}\left(t_n,x_n,y_n\right)$$

$$k_2=\mathit{hf}\left(t_n+\frac{1}{2}h,x_n+\frac{1}{2}{k}_1,y_n+\frac{1}{2}{m}_1\right)$$

$$m_2=\mathit{hg}\left(t_n+\frac{1}{2}h,x_n+\frac{1}{2}{k}_1,y_n+\frac{1}{2}{m}_1\right)$$

$$k_3=\mathit{hf}\left(t_n+\frac{1}{2}h,x_n+\frac{1}{2}{k}_2,y_n+\frac{1}{2}{m}_2\right)$$

$$m_3=\mathit{hg}\left(t_n+\frac{1}{2}h,x_n+\frac{1}{2}{k}_2,y_n+\frac{1}{2}{m}_2\right)$$

$$k_4=\mathit{hf}\left(t_n+h,x_n+k_3,y_n+m_3\right)$$

$$m_4=\mathit{hg}\left(t_n+h,x_n+k_3,y_n+m_3\right)$$

Generalizations of the formular to system of more than two equations follow a similar pattern.

Further insights and resources on the Runge- Kutta method can be found in [8, 9, 10].

Illustration 3

Given the initial value problem.

$$y^{\prime }=x+y,y\left(0\right)=1..$$

Determine.

$y$ for $x=0.2$, using the Runge-Kutta method of order four (RK4) with the step size $h=0.1$.

Solution

The RK4 formula is:

$$y_{n+1}=y_n+\frac{1}{6}\left(m_1+2m_2+2m_3+m_4\right)$$

Where

$$m_1=\mathit{hf}\left(x_n,y_n\right)$$

$$m_2=\mathit{hf}\left(x_n+\frac{1}{2}h,y_n+\frac{1}{2}{m}_1\right)$$

$$m_3=\mathit{hf}\left(x_n+\frac{1}{2}h,y_n+\frac{1}{2}{m}_2\right)$$

$$m_4=\mathit{hf}\left(x_n+h,y_n+m_3\right)$$

$$x_n=x_0+\mathit{nh}$$

$$n=\frac{x_n-x_0}{h}=\frac{0.2-0}{0.1}=2$$

Hence two iterations shall be needed.

When $n=0$

$$y_1=y_0+\frac{1}{6}\left(m_1+2m_2+2m_3+m_4\right)$$

$$m_1=\mathit{hf}\left(x_0,y_0\right)=0.1f\left(0,1\right)=0.1\left(0+1\right)=0.1$$

$$\begin{array}{l}{m}_2=\mathit{hf}\left(x_0+\frac{1}{2}h,y_0+\frac{1}{2}{m}_1\right)=0.1f\left(0+\frac{1}{2}\left(0.1\right),1+\frac{1}{2}\left(0.1\right)\right)\\ =0.1f\left(0.05,1.05\right)=0.1\left(0.05+1.05\right)=0.1\left(1.1\right)=0.11\end{array}$$

$$\begin{array}{l}{m}_3=\mathit{hf}\left(x_0+\frac{1}{2}h,y_0+\frac{1}{2}{m}_2\right)=0.1f\left(0+\frac{1}{2}\left(0.1\right),1+\frac{1}{2}\left(0.11\right)\right)\\ =0.1f\left(0.05,1.055\right)=0.1\left(0.05+1.055\right)=0.1\left(1.105\right)=0.1105\end{array}$$

$$\begin{array}{l}{m}_4=\mathit{hf}\left(x_0+h,y_0+m_3\right)=0.1f\left(0+0.1,1+0.1105\right)\\ =0.1f\left(0.1,1.1105\right)=0.1\left(0.1+1.1105\right)=0.1\left(1.2105\right)=0.12105\end{array}$$

$$\begin{array}{l}{y}_1=y_0+\frac{1}{6}\left(m_1+2m_2+2m_3+m_4\right)\\ =1+\frac{1}{6}\left(0.1+2\left(0.11\right)+2\left(0.1105\right)+0.12105\right)=1.110342\end{array}$$

For the second iteration, when n = 1:

$$y_2=y_1+\frac{1}{6}\left(m_1+2m_2+2m_3+m_4\right)$$

Where

$$x_1=x_0+h=0+0.1=0.1$$

$$\begin{array}{l}{m}_1=\mathit{hf}\left(x_1,y_1\right)=0.1f\left(0.1,1.110342\right)=0.1\left(0.1+1.110342\right)\\ =0.1\left(1.210342\right)=0.1210342\\ m_2=\mathit{hf}\left(x_1+\frac{1}{2}h,y_1+\frac{1}{2}{m}_1\right)=0.1f\left(0.1+\frac{1}{2}\left(0.1\right),1.110342+\frac{0.1210342}{2}\right)\\ =0.1f\left(0.15,\mathrm{1.1.708591}\right)=0.13208591\\ m_3=0.1f\left(0.15,1.176385\right)=0.1326385\\ m_4=0.1f\left(0.2,1.2429805\right)=0.144298049\end{array}$$

$$\begin{array}{l}{y}_2=1.110342+\frac{1}{6}\left(0.1210342+2\left(0.13208591\right)+2\left(0.1326385\right)+0.144298049\right)\\ =1.242805512\end{array}$$

Illustration 4

Consider the Lorenz dynamical system given as.

$$\frac{\mathit{dx}}{\mathit{dt}}=\sigma \left(y\left(t\right)-x\left(t\right)\right),x\left(0\right)=0.1.$$

$$\frac{\mathit{dy}}{\mathit{dt}}=-x\left(t\right)z\left(t\right)+\mathit{rx}\left(t\right)-y\left(t\right),y\left(0\right)=0.1.$$

$$\frac{\mathit{dz}}{\mathit{dt}}=x\left(t\right)y\left(t\right)-\mathit{bz}\left(t\right),z\left(0\right)=0.1.$$

Where$\sigma =10$, $r=28$ and $b=\raisebox{1ex}{$8$}\!\left/ \!\raisebox{-1ex}{$3$}\right.$with step size 0.1 determine the numerical solution at$t=10$ using the method of RK4 (Figures 5–7).

Figure 5.

Graphical representation of numerical solution of Lorenz system X(t) using Runge-Kutta method of order four for 100 iterations as in Table 3.

Figure 6.

Graphical representation of numerical solution of Lorenz system Y(t) using Runge-Kutta method of order four for 100 iterations as in Table 3.

Figure 7.

Graphical representation of numerical solution of Lorenz system X(t) using Runge-Kutta method of order four for 100 iterations as in Table 3.

Solution

The number of iteration needed is determine thus:

$$t_n=t_0+\mathit{nh}$$

$$n=\frac{t_n-t_0}{h}=\frac{10-0}{0.1}=100$$

The numerical solution is generated using Maple mathematical programming software. The codes and solution are presented below (Table 3).

Algorithm/Maple programming Code of the numerical solution of Lorenz system using RK4 method.

## This is the numerical solution of Lorenz system using RK 4method##    sys1:=diff(x(t),t)=sigma*(y(t)-x(t)),diff(y(t),t)=-x(t)*z(t)+r*x(t)-y(t),diff(z(t),t)=x(t)*y(t)-b*z(t);    ##    restart:X:=array(0..20000):Y:=array(0..20000):Z:=array(0..20000):T:=array(0..20000):A:=0.0:h:=0.1:N:=100:sigma:=10.0:b:=(8/3):r:=28.0:    for m from 0 to N do T[m]:=A+m*h:X[0]:=0.1:Y[0]:=0.1:Z[0]:=0.1:od:    for m from 0 to N-1 do K1:=h*(sigma*(Y[m]-X[m])): M1:=h*(-X[m]*Z    [m]+r*X[m]-Y[m]):N1:=h*(X[m]*Y[m]-b*Z[m]):K2:=h*(sigma*((Y[m]+    M1/2)-(X[m]+K1/2))): M2:=h*(-(X[m]+K1/2)*(Z[m]+N1/2)+r*(X[m]+    K1/2)-(Y[m]+M1/2)):N2:=h*((X[m]+K1/2)*(Y[m]+M1/2)-b*(Z[m]+N1/2))    :K3:=h*(sigma*((Y[m]+M2/2)-(X[m]+K2/2))): M3:=h*(-(X[m]+K2/2)*(Z    [m]+N2/2)+r*(X[m]+K2/2)-(Y[m]+M2/2)):N3:=h*((X[m]+K2/2)*(Y[m]+    M2/2)-b*(Z[m]+N2/2)):K4:=h*(sigma*((Y[m]+M3)-(X[m]+K3))):M4:=h*    (-(X[m]+K3)*(Z[m]+N3)+r*(X[m]+K3)-(Y[m]+M3)):N4:=h*((X[m]+K3)*(Y    [m]+M3)-b*(Z[m]+N3)):X[m+1]:=X[m]+(K1+2*K2+2*K3+K4)/6:Y[m+1]:=Y    [m]+(M1+2*M2+2*M3+M4)/6: Z[m+1]:=Z[m]+(N1+2*N2+2*N3+N4)/6:od:    for m from 0 by 10 to N do print (m, T[m], X[m],Y[m],Z[m]);od:

$t$	X(t)	Y(t)	Z(t)
1.0	$-$8.437853127	$-$9.711139533	25.75673332
2.0	$-$7.628715786	$-$6.966090160	26.86095558
3.0	$-$9.574336573	$-$9.690236844	28.21573058
4.0	$-$7.524494281	$-$8.107038156	24.99141661
5.0	$-$8.752300839	$-$7.562662622	28.79099030
6.0	$-$8.804061170	$-$9.973397587	25.79606864
7.0	$-$7.396786248	$-$6.824809818	26.47905713
8.0	$-$9.717608439	$-$9.622596776	28.65893980
9.0	$-$7.494744030	$-$8.276958009	24.64927991
10.0	$-$8.601421904	$-$7.222980039	28.83748475

Table 3.

Numerical solution of the Lorenz system using RK4 method at step size 0.1.

For further insight on programming with Maple see [6].

4. Conclusions

This chapter discusses the numerical solution of differential equation using Euler and Runge-Kutta methods. The formulas were derived and illustrations given to understand their applications. Algorithms written in Maple computational environment were provided for better understanding and further practice.