Enhanced Knowledge in Sciences and Technology Vol. 2 No.

1 (2022) 517-526
© Universiti Tun Hussein Onn Malaysia Publisher’s Office

EKST
Homepage: http://publisher.uthm.edu.my/periodicals/index.php/ekst
e-ISSN : 2773-6385

Fourth-Order Runge-Kutta Method for Solving

Applications of System of First-Order Ordinary
Differential Equations
Mimi Aida Razali1, Syahirbanun Isa1*
1
Department of Mathematics and Statistics, Faculty of Applied Sciences and
Technology,
Universiti Tun Hussein Onn Malaysia. 84600 Pagoh, Johor, MALAYSIA

*Corresponding Author Designation

DOI: https://doi.org/10.30880/ekst.2022.02.01.056
Received 02 January 2022; Accepted 02 March 2022; Available online 1 August 2022

Abstract: This paper study the applicability and the efficiency of the fourth-order
Runge-Kutta method (RK4), resulting in solving the application of system of first-
order ordinary differential equations (ODEs). Two applications of system of first-
order ODEs have been considered with initial condition which are arm race model
and drug diffusion into human body model. The real-life application of these model
has been discussed and applied in real-world. These applications have been solved by
RK4 method by using different step size with the help of MATLAB. The absolute
error analysed are also carried out explicitly in the framework. At each point of
interval, the value of the system has been calculated and compared with its exact value
at that point. Moreover, for arm race model, a comparison of the computation is
compared to existing method is worked out to illustrate the general advantage of
proposed method.

Keywords: Fourth-order Runge-Kutta method, system of Ordinary Differential

Equations, Linear Application

1. Introduction
Ordinary differential equations (ODEs) are essential to solve problems or to gain insight into the
phenomenon that arises outside mathematics. Although the basics mathematical model is discussed, the
model has solved a variety and numerous complex issues. The primary method that derives is constantly
improved and used by mathematicians to provide equations applicable in the modern world. Systems
of ODEs play vital roles not just in mathematics but also apply in another professional field. In biology,
system of ODEs is used to identify and predict the spread of disease while the diffusion of the drug into
human body uses the system of ODEs to solve the pharmacology problem.

*Corresponding author: syahir@uthm.edu.my

2022 UTHM Publisher. All rights reserved.
publisher.uthm.edu.my/periodicals/index.php/ekst
Razali et al., Enhanced Knowledge in Sciences and Technology Vol. 2 No. 1 (2022) p. 517-526

The fourth-order Runge-Kutta (RK4) is the numerical method widely used to compute an
approximation at each step of the sequences. Ten ODEs of the first-order equation with boundary
conditions are solved with three methods. Comparison results show that the RK4 method is better in all
cases among all in [1]. RK methods became important in studying explicit and implicit methods for
solving ODEs through time discretization by [3]. RK4 technique is more accurate and the numerical
solution derived by this approach converges faster. The numerical solution of the two first-order ODEs
was also calculated more efficiently using the RK4 technique than by the forward Euler technique [3].
Meanwhile, [4] showed RK5 and RK8 are less efficient than RK4 since RK4 involves less
computational time to compute truncation global error in the numerical solution.
1.1 Applications of system of first-order ODEs
The main focus of this paper is to study and understand the application that used a linear system of
first-order ODEs and solved it numerically. Therefore, in this sub-section, we will introduce two
relevant applications appropriate in the real world nowadays. In this case, we study the drug diffusion
in the gastrointestinal track (GI tract) and blood. A two-compartment model for drug absorption and
circulation through the GI tract and blood been formulated. Pharmacokinetics used to introduce
biomathematical modelling to the students. The goals for the model are to achieve the desired effect for
a medication and how long it takes to reach the desired effect in [5].
Next application that considers system of linear of first-order ODEs in the equation is economic
arm race model. Arm race models can be defined as can be thought of as enduring rivalries between
pairs of hostile powers which prompt competitive acquisition of military capability [6]. The phrase
"arms race" was coined in the nineteenth century and was widely cited as one of the causes of World
War I. He quotes Lord Grey, Britain's Foreign Secretary at the time of the conflict, as saying after the
war, "Great weapons lead inexorably to war”. If there are armaments on one side, there must be
armaments on the other sides [7]. In [8] MacKay combines Richardson's arms race equations with
Lanchester's attritional dynamics.
The objectives of this study are to study application of system of first-order ODEs by solving two
mathematical models which are arm race model and drug diffusion into human body model. Next
objective is solving the applications by using exact method and RK4 method. Then, the result of
applications of system of first order ODEs with RK4 method is compares with its exact solution and
existing method which is ZDTM method for arm race model.
2. Materials and Methods
The RK4 method is used to solve two problems regarding applications of system of first-order ODE
which is drug diffusion into human body model and arm race model.
2.1 Mathematical Model
2.1.1 Drug diffusion into human body model
Theoretically regarded as a continuous function, the flow diffuses from the GI tract, denoted as
c1 (t ) into the bloodstream, c2 (t ) before eliminated from the body. Assume the initial concentration of
drug as 𝑐𝑜 . The general formula of the equation is given as

dc1
 -k1c1
dt
𝐸𝑞. 1
dc2
 k1c1  k 2 c2
dt

with initial conditions c1 (0)  c0 and c2 (0)  0 .

518
 Razali et al., Enhanced Knowledge in Sciences and Technology Vol. 2 No. 1 (2022) p. 517-526

2.1.2 Arm race model

In this application of system of first order ODEs will consider two countries denoted as x and y
for country 1 and country 2, respectively. Let assume that the defense spending would be decreasing at
a rate proportional to the amount being spent. The two-time functions constitute the solution of the
following system of ODEs as below

dx
  kx   y  g
dt
𝐸𝑞. 2
dy
 lx   y  h
dt

where x(t) and y(t) represent the amount of weaponry at time t for country 1 and country 2 respectively,
k and l represent the efficiency increasing the armaments of x(t) and y(t) respectively. While 𝛼 and 𝛽
represent the restraining factor; g and h are grievance constants.

2.2 Numerical method

2.2.1 Fourth-order Runge-Kutta Method

The RK4 method is based on [9] and described as the following

1
y j  y j 1  h ( k1  2 k 2  2 k 3  k 4 )
6
k1  f ( t j  1 , y j  1 )

1 1
k 2  f (t j 1  h , y j 1  hk1 ) 𝐸𝑞. 3
2 2
1 1
k 3  f ( t j 1  h , y j 1  hk 2 )
2 2
k 4  f (t j 1  h, y j 1  hk 3 )

RK4 is a four-stage explicit Runge-Kutta method, which means that four function evaluations are
needed to advance the numerical solution in a one-time step.

2.3 Exact solution

Homogeneous system of ODEs can be written as follow

x 1 (t )  a11 x1  a22 x2
'

𝐸𝑞. 4
x 2 (t )  a21 x1  a22 x2
'

where xn , n  1, 2 is a function of t .

To solve homogenous system of ODEs, it can be written in a matrix form x '  Ax . Where

519
Razali et al., Enhanced Knowledge in Sciences and Technology Vol. 2 No. 1 (2022) p. 517-526

 x '1   a11 a12  x1 

 '   a   𝐸𝑞. 5
 x 2   21 a22  x2 

and A is known as the coefficient of matrix. Then, the eigenvalues are determined using

A    0 𝐸𝑞. 6

with I is an identity matrix. The roots of the characteristic’s polynomial, 1 and 2 are the
eigenvalues of A. To find the corresponding eigenvectors, we put each eigenvalue above into

( A  )v  0 𝐸𝑞. 7

 v11   v21 
Let   and   be the eigenvectors. Therefore, the general solutions become
 v12   v22 

 v11   t  v21   t
x (t )  c1   e  c2   e
1 2
𝐸𝑞. 8
 v12   v22 

Next, consider the system of non-homogenous as follow

x 1 (t )  a11 x1  a12 x2  f1
'

𝐸𝑞. 9
x 2 (t )  a21 x1  a22 x2  f 2
'

it can be written in matrix form x '  Ax  f (t ) , where

 x '1   a11 a12  x1   f1 

 '       𝐸𝑞. 10
 x 2   a21 a22  x2   f 2 

Therefore, the complementary solution is given as

 v11   t  v21 
x(t )  c1   e  c2   e  t
1 2
𝐸𝑞. 11
 v12   v22 

Then we set

x p  M (t )u (t ) 𝐸𝑞. 12

 u1 (t ) 
 u (t ) 
where u ( t )   2  is a fundamental matrix of a system x '  Ax, and let
 
 
 u n (t ) 

M '(t )  AM (t ) 𝐸𝑞. 13

520
 Razali et al., Enhanced Knowledge in Sciences and Technology Vol. 2 No. 1 (2022) p. 517-526

From Eq 12, by product rules we get

x ' p  M (t )u '(t )  u (t ) M '(t ) 𝐸𝑞. 14

1
substitute Eq. 14 into Eq. 10 and consider Eq. 12 and Eq. 13. Then multiplying both sides by M (t )

and integrate with respect to t, we conclude that a particular solution is


1
x p  M (t ) M (t ) f (t )dt 𝐸𝑞. 15

Therefore, the general solution to the system is

 v11   t  v21   t
x (t )  c1   e  c 2   e  x p (t )
1 2
𝐸𝑞. 16
 v12   v22 

3. Results
The applications with initial conditions are discussed to solved the applications problem. Drug
diffusion into human body model is denoted as test problem (a) which is the homogenous equation
while arm race model is denoted as test problem (b) which is the non-homogenous system. The exact
solution of the applications will be obtain using the method in section 2.3. The application then will be
solved numerically using RK4 method. The result for RK4 method is tested using different step size, h
which is h=0.01 and h=0.001 to determine the right step size to obtain the accuracy of RK4 method.
3.1 Test problem (a)
Consider the homogenous system

dc1
 -k1c1
dt
𝐸𝑞. 17
dc2
 k1c1  k 2 c2
dt

with initial conditions c1 (0)  500 and c2 (0)  0 .

3.1.1 Exact solution

First write the system in matrix form as follow

 c1 '    k1 0   c1 
 c '   k   𝐸𝑞. 18
 2   1 ke   c2 

Obtain the eigenvalue by 𝐸𝑞. 6 we get (  k1 )(  ke )  0 . Where, 1   k1 and 2  ke . For
case 1   k1 the system is

0 0   v1   0 
k      .
 k e  K 1   v2   0 
𝐸𝑞. 19
 1

521
Razali et al., Enhanced Knowledge in Sciences and Technology Vol. 2 No. 1 (2022) p. 517-526

k e  k1
Solving the bottom equation k1v1    ke  k1  v2  0 , and let v2  1 we get v1  .
k1

For case 2   ke , the system is

  k1  ke 0   v1  0
 k      . 𝐸𝑞. 20
 1 0   v2   0 

Solving the bottom equation ( k1  ke )v1  0v2  0 , and let v2  1 we get v1  0 . The general solution is

 ke  k1 
 c1 (t )     k1t  0   ket
   ca  k1  e  cb   e 𝐸𝑞. 21
 c2 (t )   1  1
 

ke  k1  k1
We find that ca  500 and cb  500 hence the particular solution is
k1 ke  k1

 500e
(k t )1
k1 
 1  
c ( t ) 
    k1 500(e(  k t )  e(  k t ) )  𝐸𝑞. 22
 c2 (t )   
e 1

 k1  ke 

3.1.2 Numerical results

The RK4 in 𝐸𝑞. 3 is used to find the numerical results of test problem (a). The result are calculated
using Matlab software. The performance of the fourth-order Runge Kutta method is presented in Table
1 for step size h  0.01 and Table 2 for step size h  0.001 and graphically depicted in Figure 1 and
Figure 2.
Table 1: Numerical results of RK4 method for solving test problem (a) h  0.01 .

Fourth-order Runge Kutta Method

ti Exact Solution
Numerical Solution Absolute Error
/hour
c1 c2 c1 ' c2 ' c1  c1 ' c2  c2 '
0 500 0 500 0 0 0
1 188.1064636 274.8495153 188.1064636 274.8495152 1.41113E-08 3.0914E-08
2 70.76808331 323.6872861 70.76808333 323.6872861 1.06177E-08 2.86398E-08
3 26.62386778 298.3288033 26.62386779 298.3288033 5.99182E-09 2.14817E-08
4 10.01624323 253.7385503 10.01624324 253.7385503 3.00556E-09 1.55594E-08
5 3.768240186 208.8713287 3.768240188 208.8713286 1.41342E-09 1.14278E-08
6 1.417660671 169.476797 1.417660672 169.476797 6.38094E-10 8.60857E-09
7 0.533342271 136.6109062 0.533342271 136.6109062 2.80071E-10 6.62996E-09
8 0.200650257 109.7835679 0.200650257 109.7835679 1.20418E-10 5.18806E-09
9 0.075487221 88.09921423 0.075487221 88.09921422 5.09659E-11 4.10039E-09
10 0.028399268 70.65091881 0.028399268 70.65091881 2.13044E-11 3.25981E-09
11 0.010684172 56.64064108 0.010684172 56.64064108 8.8165E-12 2.60025E-09

522
 Razali et al., Enhanced Knowledge in Sciences and Technology Vol. 2 No. 1 (2022) p. 517-526

12 0.004019524 45.40199894 0.004019524 45.40199893 3.61843E-12 2.07786E-09

13 0.001512197 36.39083117 0.001512197 36.39083117 1.47474E-12 1.66193E-09
14 0.000568908 29.16721638 0.000568908 29.16721638 5.97494E-13 1.32997E-09
15 0.000214031 23.37714252 0.000214031 23.37714252 2.40841E-13 1.06463E-09
16 8.05211E-05 18.73634104 8.05E-05 18.73634104 9.66479E-14 8.52332E-10

Table 2: Numerical results of RK4 method for solving test problem (a) with h  0.001

Fourth-order Runge Kutta Method

ti Exact Solution
Numerical Solution Absolute Error
/hour
c1 c2 c1 ' c2 ' c1  c1 ' c2  c2 '
0 500 0 500 0 0 0
1 188.1064636 274.8495153 188.1064636 274.8495153 1.27898E-12 1.26989E-08
2 70.76808331 323.6872861 70.76808331 323.6872861 1.09424E-12 1.49537E-08
3 26.62386778 298.3288033 26.62386778 298.3288033 6.18172E-13 1.37827E-08
4 10.01624323 253.7385503 10.01624323 253.7385503 2.59348E-13 1.17213E-08
5 3.768240186 208.8713287 3.768240186 208.8713286 1.37668E-13 9.64877E-09
6 1.417660671 169.476797 1.417660671 169.476797 6.35048E-14 7.82856E-09
7 0.533342271 136.6109062 0.533342271 136.6109062 2.88658E-14 6.30996E-09
8 0.200650257 109.7835679 0.200650257 109.7835679 1.24623E-14 5.07106E-09
9 0.075487221 88.09921423 0.075487221 88.09921422 5.31519E-15 4.06968E-09
10 0.028399268 70.65091881 0.028399268 70.65091881 2.21004E-15 3.26372E-09
11 0.010684172 56.64064108 0.010684172 56.64064108 8.98587E-16 2.61664E-09
12 0.004019524 45.40199894 0.004019524 45.40199893 3.84241E-16 2.09746E-09
13 0.001512197 36.39083117 0.001512197 36.39083117 1.56125E-16 1.68124E-09
14 0.000568908 29.16721638 0.000568908 29.16721638 6.15827E-17 1.34747E-09
15 0.000214031 23.37714252 0.000214031 23.37714252 2.41777E-17 1.07993E-09
16 8.05211E-05 18.73634104 8.05E-05 18.73634104 9.97466E-18 8.6553E-10

C1
C1 C1

C2 C2

Figure 1: Graph of Drug Diffusion into human body Figure 2: Graph of Drug Diffusion into human
model for exact solution and RK4 method with step body model for exact solution and RK4 method
size h=0.01 with step size h=0.001

The drug initial doses are given as 500 and the time for the drug to completely eliminated from GI
tract is calculated 16 hours. The graph in Figure 1 and Figure 2 shows that the drug absorbs into GI
tract that denoted as c1 , we can see the graph decreasing over time. Meanwhile the drug eliminated

523
Razali et al., Enhanced Knowledge in Sciences and Technology Vol. 2 No. 1 (2022) p. 517-526

from GI tract and entered bloodstream which is denoted as c 2 we can see the graph increase and at
some point, it decreases and approaching zero. At last the amount of drug is totally eliminated from the
body when the graph of c1 and c 2 is equal to zero.

Meanwhile in the Table 1 and Table 2 showed the result of exact solution and RK4 method. In
comparing the two solutions, we can see clearly that RK4 method is more preferable by using smaller
step size.
3.2 Test problem (b)
Consider the non-homogenous system

dx
 4 y  3 x  2,
dt
𝐸𝑞. 23
dy
 2x  y  2
dt

with initial conditions of x(0)  4 and y(0)  1 .

3.2.1 Exact solution

First obtain the eigenvalue from 𝐸𝑞. 6. The eigenvalues are 1  1 and 2  5 . For case 1  1 ,
 4  v  0
4
     
1
𝐸𝑞. 24
2 2   v   0  2

Solving the top equation v1  v2  0 and let v1  1 we get v1  v2 . For case 2  5 ,

 2 4   v1   0 
 2 4v   0 𝐸𝑞. 25
  2   

Solving the top equation v1  2v2  0 , we get v1  2v2 . Let v 2  1 then, v1  2 . Therefore, the
complementary solution is

1 e  2  e
   
5 t
f (t )  C C 𝐸𝑞. 26
t

1 1
1 2

 et 2e 5t 
M (t )    𝐸𝑞. 27
e
t
e 5t 

 2
Find the inverse and x p as in 𝐸𝑞. 15. we get x p    . Therefore, the general solution to the system is
 2

 x(t )   4e  2e t ( 5 t )
 2
   𝐸𝑞. 28
 y (t )   4e  e 
( 5 t )
2
t

3.2.2 Numerical results

The performance of the RK4 method is presented in Table 3 for step size h  0.01 and Table 4 for
step size h  0.001 and graphically represented in Figure 3 and Figure 4. The RK4 method will be

524
 Razali et al., Enhanced Knowledge in Sciences and Technology Vol. 2 No. 1 (2022) p. 517-526

compare to exact solution and ZDTM method that obtain from [10]. In [10], the arm race model has
been solved using ZDTM method to determine the efficiency of ZDTM method. Therefore, we compare
the result of RK4 with the ZDTM method.
Table 3: Numerical result of RK4 method for solving test problem (b) with h=0.01.

Fourth-order Runge Kutta Method ZDTM Method

Exact Solution
ti Numerical Solution Absolute Error Absolute Error
x y x' y' x  x' y  y' x  x' y  y'
0 4 1 4 1 0 0 0 0
0.1 3.633744992 1.814153013 3.633745025 1.814152996 3.28985E-08 1.6504E-08 5.28E-03 3.87E-02
0.2 3.621369915 2.517731591 3.621369955 2.517731571 3.98714E-08 2.00568E-08 4.84E-02 1.47E-01
0.3 3.845695551 3.17630507 3.845695587 3.176305052 3.62145E-08 1.8308E-08 2.01E-01 2.98E-01
0.4 4.237969357 3.831963507 4.237969386 3.831963492 2.91979E-08 1.48948E-08 6.00E-01 4.33E-01
0.5 4.75905508 4.512800084 4.759055102 4.512800073 2.20138E-08 1.14157E-08 1.46E+00 4.53E-01

Table 4: Numerical result of RK4 method for solving application test problem (b) with h=0.001.

Fourth-order Runge Kutta Method ZDTM Method

Exact Solution
ti Numerical Solution Absolute Error Absolute Error
x y x' y' x  x' y - y' x  x' y  y'
0 4 1 4 1 0 0 0 0
0.1 3.633744992 1.814153013 3.633744992 1.814153013 3.17257E-12 1.58762E-12 5.28E-03 3.87E-02
0.2 3.621369915 2.517731591 3.621369915 2.517731591 3.83604E-12 1.93756E-12 4.84E-02 1.47E-01
0.3 3.845695551 3.17630507 3.845695551 3.17630507 3.48743E-12 1.76303E-12 2.01E-01 2.98E-01
0.4 4.237969357 3.831963507 4.237969357 3.831963507 2.80309E-12 1.43885E-12 6.00E-01 4.33E-01
0.5 4.75905508 4.512800084 4.75905508 4.512800084 2.12008E-12 1.09335E-12 1.46E+00 4.53E-01

x x
y y

Figure 3: Graph of Arm race model Figure 4: Graph of Arm Race model
for exact solution, RK4 method and for exact solution, RK4 method and
ZDTM method for x and y. ZDTM method for x and y.

The arm race model solution in Figure 3 and Figure 4 shows that the graph the country 1
denoted as x increasing their armaments. In a beginning, country 2 denoted as y will decrease their
defences spending. When x increase the spending for defences, y will see the needs to increase the
spending. Meanwhile for the numerical results in Table 3 and Table 4, comparing the two methods, we
can see clearly that RK4 method is more preferable than ZDTM method and the result of RK4 method

525
Razali et al., Enhanced Knowledge in Sciences and Technology Vol. 2 No. 1 (2022) p. 517-526

with smaller step size gives better approximations to the exact solutions. If only ZDTM increasing the
order of approximate more accurate solution but its time consuming and tedious to apply.
4. Conclusion
The RK4 method has very effective approach for many applications in various field of science and
engineering. By comparing exact solution and RK4 method for both applications of system of ODEs,
RK4 method is approximately accurate since its approaching to exact solution. The result with smaller
step sizes gives better result and less absolute error compare to larger step size. Meanwhile, in
application of system of ODEs of arm race model the result of RK4 is more accurate approaching to
exact value than ZDTM method. To obtain the approximate result, ZDTM method need to increasing
the order which is tedious and time consuming for computational work. In this research, RK4 solved
the arm race model and drug diffusion into human body model which reduces the computational work
than other traditional method.
Acknowledgement
The authors would also like to thank the Faculty of Applied Sciences and Technology,
Universiti Tun Hussein Onn Malaysia for its support.
References
[1] N. Ahmad, S. Charan, and V. P. Singh, “Study of Numerical Accuracy of Runge-Kutta Second
, Third and Fourth-order Method,” vol. 4, no. 6, pp. 111–118, 2015.
[2] K. A Taher (2020). Comparison of numerical methods for solving a system of ordinary
differential equations: accuracy, stability and efficiency.
[3] V. Chauhan and P. K. Srivastava, “Computational techniques based on runge-kutta method of
various order and type for solving differential equations,” Int. J. Math. Eng. Manag. Sci., vol. 4,
no. 2, pp. 375–386, 2019, doi: 10.33889/ijmems.2019.4.2-030.
[4] M. A. Islam, “Accurate Solutions of Initial Value Problems for Ordinary Differential Equations
with the Fourth-order Runge Kutta Method,” J. Math. Res., vol. 7, no. 3, pp. 41–45, 2015, doi:
10.5539/jmr.v7n3p41.
[5] G. A. Koch-Noble, “Drugs in the classroom: Using pharmacokinetics to introduce
biomathematical modeling,” Math. Model. Nat. Phenom., vol. 6, no. 6, pp. 227–244, 2011, doi:
10.1051/mmnp/20116612.
[6] Smith, R. P. (2019). The Influence of the Richardson Arms Race Model. Pioneers in Arts,
Humanities, Science, Engineering, Practice, 25–34.
[7] Maiolo, Joseph (2016) Introduction. Ch. 1 in: Thomas Mahnken, Joseph Maiolo & David
Stevenson (eds) Arms Races in International Politics: From the Nineteenth to the Twenty-first
Century. Oxford: Oxford University Press, 1–10
[8] MacKay, Niall (2020) When Lanchester met Richardson: The interaction of warfare with
psychology. Ch. 9 in this volume.
[9] U. A. M. Roslan, Z. Salleh, and A. Kiliçman, “Solving Zhou chaotic system using fourth-order
Runge-Kutta method,” World Appl. Sci. J., vol. 21, no. 6, pp. 939–944, 2013
[10] Onkar Warade, N., & Rastogi, P. (2018). Zhou’s Differential Transformation Method for Study
of Arm Race Richardson Model. International Journal of Mathematics Trends and Technology,
57(6), 382–387. https://doi.org/10.14445/22315373/ijmtt-v57p552

526