algorithms

Article
An Efficient Third-Order Scheme Based on Runge–Kutta and
Taylor Series Expansion for Solving Initial Value Problems
Noori Y. Abdul-Hassan 1 , Zainab J. Kadum 2 and Ali Hasan Ali 1,3, *

1 Department of Mathematics, College of Education for Pure Sciences, University of Basrah, Basrah 61001, Iraq
2 Basrah Education Directorate, Ministry of Education, Basrah 61001, Iraq
3 Institute of Mathematics, University of Debrecen, Pf. 400, H-4002 Debrecen, Hungary
* Correspondence: [email protected]

Abstract: In this paper, we propose a new numerical scheme based on a variation of the standard for-
mulation of the Runge–Kutta method using Taylor series expansion for solving initial value problems
(IVPs) in ordinary differential equations. Analytically, the accuracy, consistency, and absolute stability
of the new method are discussed. It is established that the new method is consistent and stable
and has third-order convergence. Numerically, we present two models involving applications from
physics and engineering to illustrate the efficiency and accuracy of our new method and compare it
with further pertinent techniques carried out in the same order.

Keywords: numerical methods; initial value problem; autonomous equation; local truncation error;
consistency; stability

MSC: 65L05; 65L07; 65L09; 65L12; 65L20

1. Introduction
Numerical analysis is the area of mathematics that deals with computational tech-
Citation: Abdul-Hassan, N.Y.; niques for studying and finding solutions to math problems. It is an offshoot of computer
Kadum, Z.J.; Ali, A.H. An Efficient science and mathematics that develops, analyzes, and generates algorithms for numerically
Third-Order Scheme Based on solving mathematical problems. Numerical methods are typically centered on the imple-
Runge–Kutta and Taylor Series mentation of numerical algorithms. The goal of these numerical methods is to provide
Expansion for Solving Initial Value systematic procedures for numerically solving mathematical problems. The ordinary differ-
Problems. Algorithms 2024, 17, 123. ential equation (ODE) is a mathematical equation that is used in natural physical processes,
https://doi.org/10.3390/a17030123 chemistry, engineering, and other sciences. Ordinary differential equations are notoriously
Academic Editors: Alicia Cordero and difficult to approximate analytically, so the numerical treatment is crucial because it offers
Juan Ramón Torregrosa Sánchez a powerful alternate solution method for resolving the differential equation.
We frequently use initial value problems (IVPs), such that
Received: 20 November 2023
Revised: 22 February 2024 y0 = f ( x, y( x )), y( x0 ) = y0 (1)
Accepted: 14 March 2024
Published: 16 March 2024 where x is the independent variable, which might also indicate time in physical problems,
and y( x ) is the solution. Furthermore, because y( x ) can be a vector-valued function with
N-dimensions, the domain and range of f and the solution y are given by
Copyright: © 2024 by the authors.
f : R×R N → R
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
y:R→RN
conditions of the Creative Commons Furthermore, Equation (1), in which f is a function both of x and y, is known as
Attribution (CC BY) license (https:// a “non-autonomous” equation. However, by simply adding an extra variable, which is
creativecommons.org/licenses/by/
4.0/).

Algorithms 2024, 17, 123. https://doi.org/10.3390/a17030123 https://www.mdpi.com/journal/algorithms

Algorithms 2024, 17, 123 2 of 12

always equal to x, it can be easily rewritten in an equivalent “autonomous” form below, in

which f is a function of y only:

y0 = f (y( x )), y( x0 ) = y0 (2)

There are numerous well-known numerical methods that can be used to solve the
IVP in Equation (1). One category of these is third-order methods, like Ralston’s third-
order Runge–Kutta method (Ralston’s scheme) [1], the third-order Runge–Kutta method
(RK3 scheme) [2], the third-order Runge–Kutta method based on the arithmetic mean
(ARK3 scheme) [3], the third-order Heun’s method (Heun’s scheme) [4], and several other
methods. These numerical methods have been constructed to solve Equations (1) and (2)
using different techniques such as the Taylor series expansion technique, homotopy analysis
technique, quadrature formulas technique, and Adomian decomposition technique; for
more details, see [5–7].
In this study, we construct a new numerical method based on a variation of the
standard formulation of the Runge–Kutta method using Taylor series expansion. The rest
of this paper is divided into the following sections. In Section 2, we recall some definitions
and auxiliary results that we will use in our work. The derivation of the new method is
described in Section 3. Section 4 provides details on the local truncation error. Section 5
discusses the stability analysis of the suggested technique. The consistency of the new
method and its convergence are detailed in Section 6. Several numerical models are shown
in Section 7 to show the efficiency of this method and compare it with other relevant
methods. Finally, Section 8 offers the discussion and the conclusions.

2. Preliminaries

Theorem 1 ([8]). (Existence and Uniqueness Theorem).

∂f
Let f ( x, y) and ∂y be continuous functions of x and y at all points ( x, y) in some
neighborhood of the initial point ( x0 , y0 ). Then, there is a unique function y( x ) defined on
some interval [ x0 − ε, x0 + ε] satisfying

y0 = f ( x, y( x )), y( x0 ) = y0 , x ∈ [ x0 − ε, x0 + ε], ε > 0 (3)

Definition 1 ([9]). A Taylor sum or Taylor series is a mathematical function representation in

the form of a series consisting of terms calculated using the values of the function’s derivation at a
specific point and given by
∞ r
hr

∂ ∂
ψT ( x n , y n ; h ) = ∑ (r + 1)! ∂x
+f
∂y
f ( x, y) (4)
r =0

Definition 2 ([10]). Let I ⊆ R be an open interval and f : I → R be a n-times differentiable

function at α ∈ I. The Taylor polynomial Tnα f of degree n at a point α of f is the polynomial

( x − α)2 00 ( x − α)3 000 ( x − α )n (n)

Tnα f (α) = f (α) + ( x − α) f 0(α) + f (α) + f (α) + · · · + f (α) (5)
2! 3! n!

Definition 3 ([11]). The difference between the numerical solution

yn+1 = yn + hψ( xn , yn ; h) (6)

and the exact solution y( xn+1 ) is called the local truncation error (L.T.E.) for a one-step method xn
with step size h is given by

L.T.E. = y( xn+1 ) − yn+1 (7)

Algorithms 2024, 17, 123 3 of 12

Definition 4 ([12]). The numerical method is said to be stable if there exists h > 0 for each
differential equation such that changing the initial values by a fixed amount produces a bounded
change in the numerical solution for all h

Definition 5 ([13]). A numerical method is said to be consistent when all step sizes h → 0 , as it
will converge to the differential equation. In other words, we say the method is consistent if and
only if

L.T.E. b−a
Lim = 0 f or h = (8)
h →0 h N

3. Construction of the New Scheme

In order to construct new single-step methods to solve the IVP Equation (1), we rely
on a variant of the standard form of the Runge–Kutta method given by

yn+1 = yn + hψ( xn , yn ; h)

where
2m
ψ( xn , yn ; h) = ∑ wi k i +1
i =1

and
 
i −1
k1 = f ( xn , yn ), k i = f  xn + ci h, yn + h∑ a k ,
j=1 ij j
i = 2, 3, . . . , 2m + 1
j 6 =i

and
i −1
∑ aij = ci , i = 2, 3, . . . , 2m + 1
j =1

By using the Taylor series expansion of any arbitrary function ψ( xn , yn ; h) for the
non-autonomous Equation (1), we have
∞  r
1 ∂ ∂
ψT ( x n , y n ; h ) = ∑ (r + 1)! ∂x
+f
∂y
f ( x, y)
r =0

and in the autonomous case Equation (2), ψT ( xn , yn ; h) becomes

∞  r
1 ∂
ψT ( y n ; h ) = ∑ (r + 1)! f
∂y
f (y)
r =0

We now start with the numerical method by using the family of explicit Runge–Kutta
methods listed below to solve the mentioned problem in (2).

y n +1 = y n + h ( w1 k 2 + w2 k 3 ) (9)

with

k1 = f ( xn , yn ) 

k2 = f ( xn + c2 h, yn + a1 hk1 )

(10)
k3 = f ( xn + c3 h, yn + h( a2 k1 + a3 k2 )) 


Algorithms 2024, 17, 123 4 of 12

where a1 = a21 , a2 = a31 , and a3 = a32 . We must first calculate the unknowns a1 , a2 ,
a3 , w2 and w3 . Using Taylor series expansion around f (y), we obtain

k1 = f (yn ) (11)

∞
d r
 
1
k2 = ∑ ha1 k1 f (yn ) (12)
r =0
r! dy

∞
d r
 
1
k3 = ∑ r! h ( a k
2 1 + a k
3 2 )
dy
f (yn ) (13)
r =0

By expanding k1 , k2 , and k3 in Equations (11)–(13), we have


k1 = f 
h2 h2

k2 = f + ha1 f f y + 2!  f a1 2 f yy +
2 a1 3 f 3 f yyy + · · ·




3!

(14)

h2
k3 = f + h( a2 + a3 ) f f y + 2! a2 + 2a2 a3 + a3 2 f 2 f yy + 2a1 a3 f f y 2
2


   


3 1 3 2 1 3 1 3 3 2 1 2 2 · · ·

+h 2 a 2 a 3 + a 2 a 3 + 6 a 3 + 6 a 2 f f yyy + a 1 a 2 a 3 f f y f yy + 2 a 1 a 3 f f y f yy + 

Substituting Equation (14) into Equation (9) yields

1   
yn+1 = yn + (w1 + w2 ) f h + (( a2 + a3 ) w2 + a1 w1 ) f f y h2 + ( a2 + a3 )2 w2 + a1 2 w1 f f yy + 2 w2 a1 a3 f y 2 f h3 + · · · (15)
2

The Taylor series expansion of an exact solution y( x n+1 ) is given by

 
1 1 2  1 3 
y ( x n +1 ) = y ( x n ) + f h + f f y h2 + f f yy + f f y 2 h3 + f f yy + 4 f 2 f y f yy + f f y 3 h4 + . . . (16)
2 6 24
The following system of equations is obtained by expressing k1 , k2 , and k3 in the
Taylor expansion, ignoring terms with powers of h higher than 3, and then substituting
them into formula (15) and comparing them to Equation (16):

h f : w1 + w2 = 1

1
h 2 f f y : ( a 2 + a 3 ) w2 + w1 a 1 =
2

2 1
h3 f f yy : ( a2 + a3 )2 w2 + a1 2 w1 =
3

1
h 3 f f y 2 : w2 a 1 a 3 =
6
This is a system with an infinite number of solutions comprising four equations and
five unknowns. Assuming that w1 = 73 , we obtain the optimal solution listed below:

3 4 1 7
w1 = , w2 = , a1 = , a2 = −1, a3 =
7 7 6 4

Additionally, from Equation (6), we obtain c2 = 16 and c3 = 34 .

Thus, substituting the above results in Equations (9) and (10), we present the new
method, and we call it the variation Runge–Kutta method of order three (VRK3), given
as follows:
h
yn+1 = yn + (3k2 + 4k3 ) (17)
7
with
Algorithms 2024, 17, 123 5 of 12


k1 = f ( xn , yn )  

k2 = f xn + 6h , yn + 6h k1



  (18)
k3 = f xn + 3h 7


4 , yn + h −k1 + 4 k2






4. Accuracy of the New Scheme

Here, the local truncation error of the proposed scheme is investigated as follows.
The set of Equation (18), when expanded using Taylor expansion, yields

k1 = f (19)

h h2 2 h3 3
k2 = f + f fy + f f yy + f f yyy + · · · (20)
6 72 1296
   
3h 9 2 7 9 3 77 2
k3 = f + f f y + h2 f f yy + f f y 2 + h3 f f yyy − f f y f yy + · · · (21)
4 32 24 128 288

Now, substituting Equations (19)–(21) into Equation (17) yields

   
1 1 2 1 2 35 3 11 2
y n +1 = yn + h f + h2 f f y + f f yy + f f y h +3
f f yyy − f f y f yy h4 + . . .
2 6 6 864 72
Hence, from Definition 3, we have

1
L.T.E. = ( f ( f 2 f yyy + 276 f f y f yy + 36 f y 3 ))h4 + O(h5 ) (22)
864
As per Equation (22), our proposed method (VRK3) is of third order, with an L.T.E. of
fourth order.

5. Stability Analysis of New Scheme

To test the absolute stability of the presented scheme (VRK3), we use the set of
Equations (18) to derive the following:

k1 = λyn (23)
 
hλ
k2 = λyn 1+ (24)
6
!
11hλ 7h2 λ2
k3 = λyn 1+ + (25)
4 24

By substituting Equations (23)–(25) into (17) and allowing z = hλ, we obtain

Z2 z3
 
y n +1 = y n 1 + z + + (26)
2 6

Then, from Equation (26), the stability polynomial is

Z2 z3
 
y n +1  
R(z) = = 1+z+ + + o z4 (27)
yn 2 6

Utilizing the MATLAB program, Figure 1 below graphically illustrates the absolute
stability region of the Formula (27):
Algorithms 2024, 17, 123
Algorithms 2024, 17, 123 6 of 12

Figure 1. The absolute stability region of the VRK3 method.

Figure 1. The absolute stability region of the VRK3 method.
6. Consistency of the New Scheme
To explain the consistency
6. Consistency of the Newproperty of the newly proposed scheme, we adopt Definition 5.
Scheme
Therefore, by substituting Equation (22) into Equation (8), we obtain
To explain the consistency property of the newly proposed scheme, w
nition 5.Lim L.T.E.
Therefore,
1 2
by(substituting
864 f ( f f yyy + 276
Equation f y 3 ))hinto
f f y f yy + 36(22) 4 + O ( h5 )
Equation (8),(28)
we obta
= Lim =0
h →0 h h →0 h
1
𝐿. 𝑇. 𝐸.
According to Lambert ( method
[13], a numerical 𝑓(𝑓 𝑓is consistent
+ 276𝑓𝑓 if the𝑓order
+is36𝑓
bigger)ℎ
than+ 𝑂(ℎ
one. Therefore, Lim
our new method= isLim 864
consistent since it is of order three.
→ ℎ a numerical
Lambert also defines
→
method as convergent if it is ℎ consistent and stable.
Following from Equations
According (27) and [13],
to Lambert (28), this method is consistent
a numerical methodand is stable. We con-
consistent if the o
clude that the new method (VRK3) is convergent because it satisfies the consistency and
than one.
stability Therefore, our new method is consistent since it is of order three.
properties.
Lambert also defines a numerical method as convergent if it is consiste
7. Numerical Examples
Following from Equations (27) and (28), this method is consistent and stable
In this section, we introduce two models of IVPs with varying step sizes h to compare
that the new
the efficiency andmethod (VRK3)
the accuracy of theisproposed
convergent because
new method it satisfies
(VRK3 the other
scheme) with consistenc
properties.
third-order methods, like Ralston’s scheme, RK3 scheme, ARK3 scheme, and Heun’s
scheme. Here, all calculations and figures are performed using MATLAB (R2022a) software.
The numerical results are presented in Tables 1–8, and the error analysis is illustrated in
7. Numerical Examples
Figures 2–7.
In this section, we introduce two models of IVPs with varying step si
7.1. Problem 1 [14]
pare the efficiency and the accuracy of the proposed new method (VRK3
Take into consideration the first order IVP y0 = y − x2 + 1, y(0) = 0.5, with the exact
other
solutionthird-order
y = x2 + 2x − methods,
0.5e x + 1, 0 ≤like
x ≤Ralston’s
1. Tables 1–3scheme, RK3 scheme,
show the results ARK3 schem
that were obtained.
scheme. Here, all calculations and figures are performed using MATLAB
The graphs of absolute errors are shown in Figures 2–4. The comparison of CPU time
between the VRK3 scheme and other relevant third-order schemes is shown in Table 7.
ware. The numerical results are presented in Tables 1–8, and the error an
trated in Figures 2–7.

7.1. Problem 1 [14]

Take into consideration the first order IVP 𝑦 = 𝑦 − 𝑥 + 1, 𝑦(0) = 0.5
act solution 𝑦 = 𝑥 + 2𝑥 − 0.5𝑒 + 1, 0 ≤ 𝑥 ≤ 1. Tables 1–3 show the results
tained. The graphs of absolute errors are shown in Figures 2–4. The compa
time between the VRK3 scheme and other relevant third-order schemes is s
Algorithms 2024, 17, 123 7 of 12

Table 1. Comparison of the absolute errors among third-order schemes in Problem 1, for h = 10−1 .

xi Ralston’s Scheme RK3 Scheme ARK3 Scheme Heun’s Scheme VRK3 Scheme
zero zero zero zero zero zero
0.1 6.2076 × 10−6 6.2076 × 10−6 8.9854 × 10−6 3.4299 × 10−6 6.5207 × 10−7
0.2 1.2845 × 10−5 1.2845 × 10−5 1.8692 × 10−5 6.9968 × 10−6 1.1492 × 10−6
0.3 1.9932 × 10−5 1.9932 × 10−5 2.9173 × 10−5 1.0692 × 10−5 1.4515 × 10−6
0.4 2.7492 × 10−5 2.7492 × 10−5 4.0482 × 10−5 1.4502 × 10−5 1.5125 × 10−6
0.5 3.5546 × 10−5 3.5546 × 10−5 5.2680 × 10−5 1.8412 × 10−5 1.2781 × 10−6
0.6 4.4113 × 10−5 4.4113 × 10−5 6.5826 × 10−5 2.2399 × 10−5 6.8563 × 10−7
0.7 5.3212 × 10−5 5.3212 × 10−5 7.9987 × 10−5 2.6437 × 10−5 3.3777 × 10−7
0.8 6.2861 × 10−5 6.2861 × 10−5 9.5229 × 10−5 3.0492 × 10−5 1.8762 × 10−6
0.9 7.3074 × 10−5 7.3074 × 10−5 1.1162 × 10−4 3.4524 × 10−5 4.0265 × 10−6
1 8.3864 × 10−5 8.3864 × 10−5 1.2925 × 10−4 3.8482 × 10−5 6.9006 × 10−6

Table 2. Comparison of the absolute errors among third-order schemes in Problem 1, for h = 5 × 10−2 .

xi Ralston’s Scheme RK3 Scheme ARK3 Scheme Heun’s Scheme VRK3 Scheme
zero zero zero zero zero zero
0.1 7.9184 × 10−7 7.9184 × 10−7 1.1480 × 10−6 4.3572 × 10−7 7.9594 × 10−8
0.2 1.6379 × 10−6 1.6379 × 10−6 2.3876 × 10−6 8.8818 × 10−7 1.3848 × 10−7
0.3 2.5408 × 10−6 2.5408 × 10−6 3.7254 × 10−6 1.3561 × 10−6 1.7141 × 10−7
0.4 3.5031 × 10−6 3.5031 × 10−6 5.1684 × 10−6 1.8377 × 10−6 1.7229 × 10−7
0.5 4.5273 × 10−6 4.5273 × 10−6 6.7240 × 10−6 2.3307 × 10−6 1.3400 × 10−7
0.6 5.6159 × 10−6 5.6159 × 10−6 8.3997 × 10−6 2.8321 × 10−6 4.8292 × 10−8
0.7 6.7710 × 10−6 6.7710 × 10−6 1.0204 × 10−5 3.3383 × 10−6 9.4375 × 10−8
0.8 7.9947 × 10−6 7.9947 × 10−6 1.2144 × 10−5 3.8448 × 10−6 3.0504 × 10−7
0.9 9.2884 × 10−6 9.2884 × 10−6 1.4231 × 10−5 4.3460 × 10−6 5.9642 × 10−7
1 1.0653 × 10−5 1.0653 × 10−5 1.6472 × 10−5 4.8351 × 10−6 9.8318 × 10−7

Table 3. Comparison of the absolute errors among third-order schemes in Problem 1, for h = 25 × 10−3 .

xi Ralston’s Scheme RK3 Scheme ARK3 Scheme Heun’s Scheme VRK3 Scheme
zero zero zero zero zero zero
0.1 9.9973 × 10−8 9.9973 × 10−8 1.4505 × 10−7 5.4894 × 10−8 9.8153 × 10−9
0.2 2.0675 × 10−7 2.0675 × 10−7 3.0165 × 10−7 1.1185 × 10−7 1.6954 × 10−8
0.3 3.2066 × 10−7 3.2066 × 10−7 4.7062 × 10−7 1.7070 × 10−7 2.0745 × 10−8
0.4 4.4202 × 10−7 4.4202 × 10−7 6.5283 × 10−7 2.3121 × 10−7 2.0405 × 10−8
0.5 5.7114 × 10−7 5.7114 × 10−7 8.4920 × 10−7 2.9308 × 10−7 1.5023 × 10−8
0.6 7.0830 × 10−7 7.0830 × 10−7 1.0607 × 10−6 3.5592 × 10−7 3.5422 × 10−9
0.7 8.5378 × 10−7 8.5378 × 10−7 1.2883 × 10−6 4.1926 × 10−7 1.5261 × 10−8
0.8 1.0078 × 10−6 1.0078 × 10−6 1.5331 × 10−6 4.8250 × 10−7 4.2799 × 10−8
0.9 1.1705 × 10−6 1.1705 × 10−6 1.7962 × 10−6 5.4492 × 10−7 8.0702 × 10−8
Algorithms
1 2024, 17, 123 1.3422 × 10−6 1.3422 × 10−6 2.0787 × 10−6 6.0565 × 10−7 8 of 13 1.3084 × 10−7

10-4
1.4

Ralston’s method
1.2 RK3 method
ARK3 method
Heun’s method
1 VRK3 method

0.8

0.6

0.4

0.2

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
xi

Figure
Figure2. The absolute
2. The errors for
absolute numerical
errors results in Table
for numerical 1.
results in Table 1.
 0.2

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
xi
Algorithms 2024, 17, 123 8 of 12
Figure 2. The absolute errors for numerical results in Table 1.

Absolute Error

Algorithms 2024, 17, 123 9 of 13

Figure 3.
Figure Theabsolute
3. The absoluteerrors
errorsfor
for numerical
numerical results
results in
in Table
Table 2.
2.
Absolute Error

Theabsolute
Figure4.4.The
Figure absoluteerrors
errorsfor
fornumerical
numericalresults
resultsininTable
Table3.3.

7.2.
7.2.Problem
Problem22(Mixture
(MixtureModel)
Model)
We
Weconsider here the
consider here theIVP
IVPproposed
proposedinin [15],
[15], which
which waswas a model
a model of aof a storage
storage tank tank
in an in
oil
an oil refinery
refinery that holds holdsgal2000
that 2000 gal of with
of gasoline gasoline
100 lbwith 100
of an lb of an
additive additive
mixed withinmixed
it. To within
prepare
it.forTo
winter weather,
prepare 40 gal/min
for winter weather, 40 gal/min
of gasoline that contains 2 lb of
of gasoline additive
that containsper 2gallon
lb ofisadditive
pumped
intogallon
per the storage tank. into
is pumped The the
well-mixed solution
storage tank. Theiswell-mixed
pumped out at a rateisof
solution 45 gal/min.
pumped out atLet
a
y beofthe45amount
rate gal/min. Let 𝑦 be(in
of additive thepounds)
amountinofthe tank at(in
additive time x. When
pounds) = tank
in xthe 0, weat know 𝑥.
timethat
y = 100.
When 𝑥 =The0,mixture
we know that 𝑦is =
process 100. The
modeled bymixture y0 = 80 −
the IVP, process 45
2000−5x) by the IVP, 𝑦 =
is (modeled y , y(0) = 100, and
80
the−analytic solution,
𝑦 , 𝑦(0)y(x= )=100 , and
(2000 − 5x){the 3900
2 −analytic solution, 𝑦(𝑥) = (2000 − 5𝑥){2 −
8
9 (2000 − 5x ) }, 0 ≤ x ≤ 1.
( ) (2000)
(2000 4–6
Tables − 5𝑥)
show}, 0the
≤ absolute
𝑥 ≤ 1. errors among third-order methods and the VRK3 scheme,
( )
withTables
different
4–6step sizesthe
show = 0.1, h errors
of habsolute = 0.05, among
and h =third-order
0.025. Figures 5–7 depict
methods andthe graphical
the VRK3
analysis used to support the numerical results in Tables 4–6. A comparison
scheme, with different step sizes of ℎ = 0.1, ℎ = 0.05, and ℎ = 0.025. Figures 5–7 depictof CPU time
between the new method and other third-order schemes is shown in Table 8.
the graphical analysis used to support the numerical results in Tables 4–6. A comparison
of CPU time between the new method and other third-order schemes is shown in Table 8.

Table 4. Comparison of the absolute errors among third-order schemes in Problem 2, for ℎ = 10 .

𝒙𝒊 Ralston’s Scheme RK3 Scheme ARK3 Scheme Heun’s Scheme VRK3 Scheme
zero zero zero zero zero zero
0.1 3.7246 × 10−9 4.1133 × 10−9 3.8694 × 10−9 3.3207 × 10−9 3.1607 × 10−9
Algorithms 2024, 17, 123 9 of 12

Table 4. Comparison of the absolute errors among third-order schemes in Problem 2, for h = 10−1 .

xi Ralston’s Scheme RK3 Scheme ARK3 Scheme Heun’s Scheme VRK3 Scheme
zero zero zero zero zero zero
0.1 3.7246 × 10−9 4.1133 × 10−9 3.8694 × 10−9 3.3207 × 10−9 3.1607 × 10−9
0.2 7.4366 × 10−9 8.2127 × 10−9 7.7257 × 10−9 6.6303 × 10−9 6.3109 × 10−9
0.3 1.1136 × 10−8 1.2298 × 10−8 1.1568 × 10−8 9.9282 × 10−9 9.4499 × 10−9
0.4 1.4821 × 10−8 1.6367 × 10−8 1.5397 × 10−8 1.3214 × 10−8 1.2577 × 10−8
0.5 1.8494 × 10−8 2.0424 × 10−8 1.9213 × 10−8 1.6488 × 10−8 1.5694 × 10−8
0.6 2.2154 × 10−8 2.4466 × 10−8 2.3015 × 10−8 1.9752 × 10−8 1.8800 × 10−8
0.7 2.5801 × 10−8 2.8493 × 10−8 2.6804 × 10−8 2.3003 × 10−8 2.1895 × 10−8
0.8 2.9435 × 10−8 3.2507 × 10−8 3.0579 × 10−8 2.6243 × 10−8 2.4979 × 10−8
0.9 3.3057 × 10−8 3.6506 × 10−8 3.4342 × 10−8 2.9472 × 10−8 2.8053 × 10−8
1 2024, 17, 123
Algorithms 3.6665 × 10−8 4.0491 × 10−8 3.8090 × 10−8 3.2689 × 10−8 10 of 13 3.1114 × 10−8

Comparison
Table5.5.Comparison
Table of absolute
of the the absolute
errorserrors
amongamong third-order
third-order schemes
schemes in for ℎ = 2, for h
Problemin2,Problem = 5 × 10−2 .
5 × 10 .
xi 𝒙 Ralston’s
Ralston’s Scheme
Scheme RK3 Scheme
RK3 Scheme ARK3 Scheme ARK3 Scheme
Heun’s Scheme Heun’s VRK3Scheme
Scheme VRK3 Scheme
𝒊
zerozero zero
zero zero zero zero zero zero zerozero zero
0.1 0.1 4.6501× 10
4.6501 10−10
× −10 × 10−10 × 10−104.8308 × 10−10
5.1354 5.1354 4.8308 × 10 −10
4.1457 × 10−10 4.1457 × 10−
3.9459 × 10
10−10 3.9459 × 10−10
0.2 0.2 × −10 × −9 × −10 × −10 −10 7.8832 × 10−10
9.2889 × 10 10
9.2889 −10 1.0258
1.0258 × 10 −9 10 9.6499 × 10 9.6499
−10 10
8.2821 × 10 −10 8.2821 10
7.8832 × 10
0.3 0.3 1.3910
1.3910 × −910−9
× 10 1.5361 1.5361
−9 10−9 × 10−9
× 10−9 × 10 1.4451 × 10−91.4451 × 1.2403 1.2403 × 10−
1.1805
9
× 10−9 1.1805 × 10−9
0.4 0.4 1.8506
1.8506 × −910−9
× 10 2.0437
2.0437 × 10−9 × 10 −9
1.9225 × 10−91.9225 × 10 −9
1.6499 × 10−9 1.6499 × 10 −9
1.5704 × 10−9 1.5704 × 10−9
0.5 2.3096 × −910−9 2.5506 × 10−9 2.3994 × 10−9 2.0592 × 10−9 −9 1.9600 × 10−9
0.5 2.3096 × 10 −9 2.5506 × 10−9 − 2.3994 × 10−9 2.0592
− × 10−9 1.9600 × 10
0.6 2.7670 × 10 3.0556 × 10 9 2.8745 × 10 9 2.4671 × 10−9 −9 2.3482 × 10−9
0.6 2.7670 × 10−9 −9 3.0556 × 10−9 − 2.8745 × 10−9 2.4671
− × 10−9 2.3482 × 10
0.7 3.2222 × 10 3.5584 × 10 9 3.3474 × 10 9 2.8729 × 10−9 2.7345 × 10−9
0.8 0.7 3.2222 × 10 3.5584 4.0598
× 10−9 × 10−93.3474 × 10−93.8191 × 2.8729
10−9 × 10 2.7345 × 910−9
−9 −9
3.6763 × 10−9 3.2778 × 10− 3.1199 × 10−9
0.9 0.8 3.6763 × 10 4.0598 4.5595
× 10−9 × 10−93.8191 × 10−94.2892 × 3.2778 × 10 3.6813 × 10 × 910
3.1199
−9 −9 −9
4.1288 × 10−9 10 − 9 − 3.5040 × 10−9
1 0.9 4.1288
4.5789 × −910−9
× 10 4.5595 5.0566
−9 −
× 10 × 10 4.2892 × 10 4.7568 × 3.6813
9 −9
10 − 9 × 10 −9 3.5040
4.0826 −
× 10 × 10−9
9 3.8859 × 10−9
1 4.5789 × 10−9 5.0566 × 10−9 4.7568 × 10−9 4.0826 × 10−9 3.8859 × 10−9

Table
Table6.6.Comparison
Comparisonof the absolute
of the errorserrors
absolute amongamong
third-order schemes schemes
third-order in Problem for ℎ = 2, for h
in2,Problem = 25 × 10−3 .
25 × 10 .

xi 𝒙𝒊 Ralston’s Scheme
Ralston’s Scheme RK3 Scheme
RK3 Scheme ARK3 Scheme Heun’s Scheme Heun’s
ARK3 Scheme VRK3 Scheme
Scheme VRK3 Scheme
zerozero zero
zero zero zero zero zero zero zerozero zero
0.1 0.1 5.7938
5.7938× 10 10−11
× −11 6.4006 6.4006
× 10−11 × 10−116.0197 × 10−11 5.1628
6.0197 × 10 −11 × 10−11 4.9127
5.1628 × 10− × 11
10−11 4.9127 × 10−11
0.2 0.2 1.1617
1.1617× 10 10−10
× −10 × 10−10 × 10−101.2069 × 10−10
1.2828 1.2828 −10 × 10−10
1.0360
1.2069 × 10 × 10−
9.8595
1.0360 × 10
10−11 9.8595 × 10−11
0.3 0.3 1.7408 × 10−10
1.7408× 10 −10 1.9220 1.9220
−10 − 10
× 10 × 10 1.8083 × 10 1.8083 × 10
−10 − 10
1.5525 × 10 −10 1.4776
1.5525 − 10
× 10 × 10−10 1.4776 × 10−10
0.4 0.4 2.3087× 10
2.3087 10−10
× −10 × 10−10 × 10−102.3985 × 10−10
2.5497 2.5497 2.3985 × 10 −10
2.0580 × 10−10 2.0580 × 10−
1.9583 × 10
10−10 1.9583 × 10−10
0.5 0.5 × 10−10
2.8851× 10
2.8851 −10 3.1858 3.1858
−10
− 10
× 10 × 10 2.9971 × 10 2.9971 × 10
−10
− 10
2.5724 × 10 −10 2.5724
2.4477
−
× 10 ×10 10−10 2.4477 × 10−10
0.6 0.6 3.4589
3.4589 × 10 10−10
× −10 3.8193
3.8193 × 10 −10 × 10 −10
3.5934 × 10 3.5934
−10 × 10 −10
3.0846 × 10 −10 3.0846 × 10 −10
2.9351 × 10−10 2.9351 × 10−10
0.7 0.7 4.0254 × 10−10
4.0254 × 10−10 4.4454 4.4454
−10
× 10−10 × 10−104.1820 × 10−10 4.1820 × 10 −10
3.5894 × 10 −10 3.5894 × 10−
3.4157
10
× 10−10 3.4157 × 10−10
0.8 10−10
4.5941 × −10 5.0736 × 10 4.7731 × 10 −10 4.0967 × 10 −10 3.8989 × 10−10
0.8 4.5941 × 10 −10 5.0736 × 10−10 4.7731 × 10−10 4.0967 × 10−10 3.8989 × 10−10
0.9 5.1622 × −10
10 5.7003 × 10−10 5.3629 × 10−10 4.6035 × 10−10 −10 4.3809 × 10−10
0.9 5.1622 × 10 −10 5.7003 × 10−10 − 5.3629 × 10−10 4.6035
− × 10−10 4.3809 × 10
1 5.7196 × 10 6.3162 × 10 10 5.9421 × 10 10 5.0994 × 10−10 4.8530 × 10−10
1 5.7196 × 10−10 6.3162 × 10−10 5.9421 × 10−10 5.0994 × 10−10 4.8530 × 10−10
Absolute Error

Figure
Figure5. The absolute
5. The errors for
absolute numerical
errors results in Table
for numerical 4.
results in Table 4.
Algorithms 2024, 17, 123 11 of 13

Algorithms 2024,
Algorithms 2024, 17, 123 11 of 1310 of 12

Absolute
Absolute Error
Error

Figure
Figure6.6.6.The
Figure
The
Theabsolute
absolute
absolute
errors
errors
errors
for
forfor
numerical
numerical
numerical
results
results
results
in Table
in Table
in Table
5.
5. 5.
Error
AbsoluteError
Absolute

Figure 7. The absolute errors for numerical results in Table 6.

Theabsolute
Figure7.7.The
Figure absolute errors
errors for
for numerical
numerical results
results in
in Table
Table 6.
6.
Table 7. Comparisons of CPU time in Problem 1, for different step sizes h.
Table
Table7.
7.Comparisons
Comparisons ofof CPU
CPU time in Problem
Problem 1,
1, for
for different
differentstep
stepsizes
sizesh.h.
CPU Time
Step Size
Ralston’s Scheme RK3 Scheme ARK3 Scheme
CPUCPU TimeHeun’s Scheme
Time VRK3 Scheme
Step Size
h = 0.1Step Size Ralston’s
0.003325 0.003508 0.004631 0.005653Scheme 0.001558
Scheme
Ralston’s Scheme RK3 Scheme
RK3 Scheme ARK3 Scheme
ARK3 Scheme Heun’sHeun’s Scheme VRK3
VRK3Scheme
Scheme
h =
h = 0.10.05 0.003407
0.003325 0.004827
0.003508 0.003017
0.004631 0.005066
0.005653 0.001027
0.001558
h = 0.1 0.003325 0.003508 0.004631 0.005653 0.001558
h =h 0.05
= 0.025 0.003655
0.003407 0.003144
0.004827 0.005021
0.003017 0.004364
0.005066 0.001273
0.001027
h = 0.05 0.003407 0.004827 0.003017 0.005066 0.001027
h = 0.025 0.003655 0.003144 0.005021 0.004364 0.001273
h = 0.025 Table 8. Comparisons
0.003655 of CPU time in Problem
0.003144 2, for different step0.004364
0.005021 sizes h. 0.001273
CPU
Table 8. Comparisons of CPU time Time 2, for different step sizes h.
in Problem
Step Size
Ralston’s Scheme RK3 Schemeof CPUARK3
Table 8. Comparisons time inScheme
Problem 2,HEUN’S SCHEME
step sizes h.VRK3 Scheme
hStep
= 0.1Size 0.004160 0.003116
CPU
0.004698
Time for different
0.004627 0.001358
h = 0.05 Ralston’s Scheme
0.003328 RK3 Scheme
0.004782 ARK3
0.004879Scheme
CPU Time HEUN’S
0.003469SCHEME 0.001003 VRK3 Scheme
h =h 0.1 Step Size 0.004160 0.003116 0.004698 0.004627 0.001358
= 0.025 0.003506 Scheme
Ralston’s 0.004392
RK3 Scheme 0.003230
ARK3 Scheme 0.004234
HEUN’S SCHEME 0.001296 VRK3 Scheme
h = 0.05 0.003328 0.004782 0.004879 0.003469 0.001003
h = 0.1 0.004160 0.003116 0.004698 0.004627 0.001358
h = 0.025 0.003506 0.004392 0.003230 0.004234 0.001296
h = 0.05 0.003328 0.004782 0.004879 0.003469 0.001003
h = 0.025 0.003506 0.004392 0.003230 0.004234 0.001296
Algorithms 2024, 17, 123 11 of 12

8. Discussion and Conclusions

In this study, we introduced an innovative third-order method designed for solving
initial value problems (IVPs). Our approach is rooted in a novel adaptation of the standard
formulation employed in Runge–Kutta methods, incorporating Taylor series expansion. To
validate the effectiveness of this new method, we employed two distinct numerical models,
effectively showcasing its fundamental capabilities. It is important to underscore that all our
numerical findings, including the accompanying tables and figures, were calculated using
MATLAB (R2022a) software on a dedicated computer system operating with Windows 11
Pro. The system uses an 11th Generation Intel(R) Core (TM) i7-11800H processor running
at 2.30 GHz, backed by 16.0 GB of RAM (15.7 GB usable).
A comprehensive numerical assessment was conducted using Tables 1–6, which
present an intricate comparison of absolute errors across various step sizes, specifically
h = 10−1 , h = 5 × 10−2 , and h = 25 × 10−3 . Through the graphical representations
found in Figures 2–7, we were able to discern that our novel method, referred to as VRK3,
consistently outperformed several benchmark techniques including Ralston’s scheme,
RK3 scheme, ARK3 scheme, and Heun’s scheme. This superiority primarily stems from
the reduced local truncation error of VRK3. Additionally, our investigation revealed
a significant insight regarding the impact of step size on accuracy. As we decreased
the step size, the error progressively approached zero, strongly indicating that precision
increased with smaller step sizes. This observation reinforces the importance of carefully
selecting step sizes to achieve higher levels of accuracy in numerical solutions. Turning
our attention to computational efficiency, Tables 7 and 8 provided valuable insights. The
VRK3 scheme consistently demonstrated reduced CPU time compared to its counterparts,
further validating its utility in practical applications. Furthermore, Figure 1 depicts the
stability region of our third-order VRK3 scheme, establishing its equivalence to similar
methodologies. Importantly, we substantiated the convergence of our VRK3 scheme, as it
satisfies both the consistency and stability criteria.
In conclusion, our newly proposed third-order method exhibits a commendable blend
of efficiency and reliability. The method’s stability and high accuracy render it partic-
ularly robust for a wide range of applications. This research contributes to the field of
numerical methods for IVPs by presenting an innovative approach that holds promise for
improving computational accuracy and efficiency. Future research directions might explore
the extension of this method to more complex problems or its integration into broader
computational frameworks.

Author Contributions: Conceptualization, N.Y.A.-H., Z.J.K. and A.H.A.; formal analysis, N.Y.A.-H.,
Z.J.K. and A.H.A.; investigation, N.Y.A.-H., Z.J.K. and A.H.A.; methodology, N.Y.A.-H., Z.J.K. and
A.H.A.; software, A.H.A.; supervision, N.Y.A.-H.; writing—original draft, N.Y.A.-H. and Z.J.K.;
writing—review and editing, A.H.A. All authors have read and agreed to the published version of
the manuscript.
Funding: The author received no direct funding for this work.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: Data are contained within the article.
Acknowledgments: The author would like to express gratitude to the anonymous referees for their
valuable comments and suggestions.
Conflicts of Interest: The authors declare that they have no competing interests.

References
1. Ralston, A. Runge-Kutta Methods with Minimum Error Bounds. Math. Comput. 1962, 16, 431–437. [CrossRef]
2. Ralston, A.; Rabinowitz, P. A First Course in Numerical Analysis, 2nd ed.; McGraw-Hill: New York, NY, USA, 1978.
Algorithms 2024, 17, 123 12 of 12

3. Wazwaz, A.M. A Comparison of Modified Runge-Kutta Formulas Based on a Variety of Means. Int. J. Comput. Math. 1994, 50,
105–112. [CrossRef]
4. Hatun, M.; Vatansever, F. Differential Equation Solver Simulator for Runge-Kutta Methods. Uludağ Univ. J. Fac. Eng. 2016, 21, 145–162.
5. Burden, R.L.; Faires, J.D. Numerical Analysis, 9th ed.; Brooks/Cole Publishing Company: Pacific Grove, CA, USA, 2011.
6. Ahmad, N.; Charan, S. A Comparative Study on Numerical Solution of Ordinary Differential Equation by Different Method with
Initial Value Problem. Int. J. Recent. Sci. Res. 2017, 8, 21134–21139.
7. Nhawu, G.; Mafuta, P.; Mushanyu, J. The Adomian Decomposition Method for Numerical Solution of First-Order Differential
Equations. J. Math. Comput. Sci. 2016, 6, 307–314.
8. Nagle, R.K.; Saff, E.B.; Snider, A.D. Fundamentals of Differential Equations, 9th ed.; Pearson: Boston, MA, USA, 2018.
9. Wusu, A.S.; Akanbi, M.A.; Okunuga, S.A. A Three-stage Multiderivative Explicit Runge-Kutta Method. Am. J. Comput. Math.
2013, 3, 121–126.
10. Ali, A.H.; Pales, Z.S. Taylor-type Expansions in Terms of Exponential Polynomials. Math. Inequalities Appl. 2022, 25, 1123–1141.
[CrossRef]
11. Kadum, Z.J.; Abdul-Hassan, N.Y. New Numerical Methods for Solving the Initial Value Problem Based on a Symmetrical
Quadrature Integration Formula Using Hybrid Functions. Symmetry 2023, 15, 631. [CrossRef]
12. Corless, R.M.; Kaya, C.Y.; Moir, R.H.C. Optimal Residuals and the Dahlquist Test Problem. Numer. Algorithms 2019, 81, 1253–1274.
[CrossRef]
13. Lambert, J.D. Computational Methods in Ordinary Differential Equations; John Wiley & Sons Inc.: New York, NY, USA, 1973.
14. Ram, T.; Solangi, M.A.; Sanghah, A.A. A Hybrid Numerical Method with Greater Efficiency for Solving Initial Value Problems.
Math. Theory Model. 2020, 10, 2224–5804.
15. Omar, Z.; Adeyeye, O. Numerical Solution of First Order Initial Value Problems Using a Self-Starting Implicit Two-Step
Obrechkoff-Type Block Method. J. Math. Stat. 2016, 12, 127–134. [CrossRef]

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual
author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to
people or property resulting from any ideas, methods, instructions or products referred to in the content.