Comparative Analysis of Second- and Fourth-Order Runge–Kutta Methods for Solving Chaotic Dynamical Systems

1. Introduction

The study focuses on investigating the efficiency and effectiveness of solving the Lorenz, Genesio Tesi and Rosler systems of equations by using the Midpoint, Improved Euler’s, Ralson methods and Runge Kutta method of order four. The Lorenz equations is a system of ordinary differential equations which were first studied by a mathematician and meteorologist Edward Lorenz. This system is used to model the unpredictable behavior of weather. They represent the convective motion of fluid cell that is warmed from below and cooled from above [1]. The Genesio-Tesi system is a system of ordinary differential equations which were implemented by Roberto Genesio and Alberto Tesi in 1992 as an electronic circuit. They are made of three-dimensional differential equations with one of them being quadratic. The Rossler system is a system of three nonlinear ordinary differential equations which were studied by Otto Rossler in the 1970s [2]. The ordinary differential equations in the system define a continuous-time dynamical system that exhibits chaotic dynamics associated with fractal properties of the system.

The chaotic system of equations will be solved by using the numerical methods, which are mathematical tools that are designed to solve numerical problems in the study of numerical analysis. In numerical analysis, these numerical methods attempt to find approximate solutions of problems than exact solutions [3]. In theory, a numerical method begins with an initial point and then moves forward in time to find the next solution point. The procedure is repeated in order to map out the solution. To determine the current value, single-step methods refer to only one previous point and its derivative. Runge-Kutta methods, for example, take some intermediate steps to obtain a higher order method, but then discard all previous information before proceeding to the next step, which in contrary to the Multi-step methods which improve efficiency by retaining and using information from previous steps rather than discarding it. In this project the selected numerical methods, the fourth order Runge-Kutta and second order Runge-Kutta methods are used to solve the mathematical problems that are chaotic in nature in order to investigate the accuracy and efficiency of the methods.

The Runge-Kutta methods were developed around the 1990s by two German mathematicians Carl Runge and Wilhelm Kutta [4]. They are a family of implicit and explicit iterative single step methods including the Euler’s method. Runge-Kutta methods are effective and widely used family of numerical methods for solving the initial-value problems of differential equations. The primary disadvantage of these methods is that they require significantly more computational time than multi-step methods whereas their main advantages are that they are easy to implement and stable [5]. We have single step methods such as the Euler, Heun’s, Midpoint and the Ralston methods, which are used to find approximate solutions of Initial Value Problems (IVPs). The main advantage of these other single step methods is that they are simple and direct. They mostly require less evaluations at a point than the Runge-Kutta methods that requires more than just one evaluation at a point to obtain a solution. These methods require less computational time than the Runge-Kutta methods. Although this is the case, the Runge-Kutta methods are more stable than these other single step methods. For instance, the Euler’s methods are conditionally stable. They are only stable when the step size is significantly small [7].

The theory of chaos dates back to 1963 when it was first reported by Lorenz for a system of ordinary differential equations modeling weather phenomena [11]. Since then, a large number of chaotic phenomena and chaotic behavior have been observed and reported in a wide range of systems, including electrical circuits, lasers, fluid dynamics, mechanical devices, population growth, and many other scientific applications [12]. The chaos theory is a study of random or unpredictable behavior in systems [13]. Such behaviour is mostly studied in natural systems such as the weather, biological and physical processes, and electromagnetic processes. Chaos theory states that a system is chaotic if ordered regular patterns can be seen to arise out of random, under certain conditions [14].

The selected chaotic systems are considered to be highly sensitive to their initial conditions, that is, a small change in their initial conditions results in a significant change in the solution. This motivates the investigation and applicability of the fourth order Runge-Kutta and three step Adams Moulton numerical methods. These study presents a novel comparisons of solutions obtained by the Midpoints, improved Eulers, Ralston and Runge Kutta method of order four when solving the Lorenz, Genessio Tesi and Rosler systems of chaotic equations.

The existing literature contains numerous studies on numerical approximations of chaotic systems using Runge–Kutta methods. However, limited studies provide a systematic comparison of multiple second-order Runge–Kutta variants against the fourth-order Runge–Kutta method across different chaotic systems while simultaneously assessing convergence properties, stability characteristics, and computational efficiency. This study addresses this gap through a unified comparative framework using the Lorenz, Genesio–Tesi, and Rössler systems.

Section 2 presents the literature review. Section 3 presents the methodology. Section 4 discusses numerical results, while Section 5 concludes the study.

2. Problem Formulation

The problems of interest in this research are the Lorenz, Genezio Tesi, and the Rossler chaotic systems. In this chapter, the stability analysis of these systems is done using the Lyapunov stability analysis in the first approximation [23]. The second and fourth orders Runge Kutta methods are derived and applied to the problems.

3. Theoretical Results

3.1. Theorem 1: Existence and Uniqueness of Solution

Theorem

Consider the initial value problem

$${\displaystyle \frac{dy}{dt}}=f(t,y),y\left(t_0\right)=y_0,$$

where $\left(f\right(t,y\left)\right)$ is continuous in a region

$$R=(t,y):|t-t_0|\le a,|y-y_0|\le b,$$

and satisfies a Lipschitz condition in (y), that is,

$$|f(t,y_1)-f(t,y_2)|\le L|y_1-y_2|,$$

for some positive constant $\left(L\right)$.

Then there exists a unique solution $\left(y\right(t\left)\right)$ of the initial value problem in some interval containing $\left(t_0\right)$.

Proof

The proof follows directly from the Picard–Lindelöf theorem. Since $\left(f\right)$ is continuous and satisfies the Lipschitz condition, successive approximations generated through Picard iteration converge uniformly to a unique solution.

3.2. Theorem 2: Consistency of the Runge–Kutta Methods

Theorem

The Midpoint, Improved Euler, Ralston, and RK4 methods are consistent numerical methods.

Proof

A numerical method is consistent if its local truncation error satisfies

$$\underset{h\to 0}{lim}{\displaystyle \frac{\tau \left(h\right)}{h}}=0.$$

For RK2 methods,

$$\tau \left(h\right)=O\left(h^3\right),$$

while for RK4,

$$\tau \left(h\right)=O\left(h^5\right).$$

Thus,

$$\underset{h\to 0}{lim}\tau \left(h\right)=0,$$

showing that the methods are consistent.

3.3. Theorem 3: Order of Accuracy of RK2 and RK4

Theorem

The Midpoint, Improved Euler, and Ralston methods possess second-order accuracy, while RK4 possesses fourth-order accuracy.

Proof

Expanding the exact solution using Taylor series gives

$$y(t+h)=y\left(t\right)+h{y}^{\prime }+{\displaystyle \frac{h^2}{2}}{y}^{″}+{\displaystyle \frac{h^3}{6}}{y}^{‴}+\dots$$

Comparing the Taylor expansion with the RK2 approximations shows agreement through terms of order $\left(h^2\right)$, producing local truncation error

$$O\left(h^3\right),$$

and global error

$$O\left(h^2\right).$$

Similarly, RK4 matches terms up to $\left(h^4\right)$, resulting in

$$O\left(h^5\right)$$

local truncation error and

$$O\left(h^4\right)$$

global error.

Hence the RK2 methods are second-order accurate and RK4 is fourth-order accurate.

3.4. Theorem 4: Convergence of the Runge–Kutta Methods

Theorem

If a numerical method is both consistent and stable, then it is convergent.

Proof

By the Lax Equivalence principle,

$$\mathrm{Consistency}+\mathrm{Stability}\Rightarrow \mathrm{Convergence}.$$

Since the considered Runge–Kutta methods are consistent and stable for sufficiently small step sizes, they are convergent.

3.5. Theorem 5: Global Error Bound

Theorem

Assume that (f(t,y)) satisfies the Lipschitz condition with constant (L). Then the global error satisfies

$$|y\left(t_n\right)-y_n|\le {\displaystyle \frac{e^{L(t_n-t_0)}-1}{L}}\tau \left(h\right),$$

where $\left(\tau \right(h\left)\right)$ denotes the local truncation error.

Proof

Using Gronwall’s inequality together with the local truncation error estimate gives the stated bound.

4. Chaotic Systems and Stability Analysis

4.1. Lorenz System

Edward Lorenz derived the system [24], given by the following system of ordinary differential equations

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\displaystyle \frac{du}{dt}}\hfill & =a_1(v-u),\hfill \\ {\displaystyle \frac{dv}{dt}}\hfill & =-us+a_{2}u-v,\hfill \\ {\displaystyle \frac{ds}{dt}}\hfill & =uv-a_{3}s,\hfill \end{array}\right.\end{array}$$

(1)

where $a_1$, $a_2$, and $a_3$ are greater than zero. These equations (1) were derived by Lorenz to model two-dimensional fluid cell between two parallel plates at different temperatures [25]. In this project we study the case where $a_1=10$, $a_2=28$, and $a_3={\displaystyle \frac{8}{3}}$ with initial conditions $u\left(0\right)=1$, $v\left(0\right)=5$, and $s\left(0\right)=10$, which is the case study that Lorenz conducted [26].

4.1.1. Stability Analysis of the Lorenz System

Given the Lorenz system (1), we set ${\displaystyle \frac{du}{dt}}=0$, ${\displaystyle \frac{dv}{dt}}=0$, and ${\displaystyle \frac{ds}{dt}}=0$ to determine the equilibrium points of the system [27]. In setting the Lorenz system equations to zero, we obtain the following equations

$$\begin{array}{cc}\hfill a_1(v-u)& =0,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill -us+a_{2}u-v& =0,\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill uv-a_{3}s& =0.\hfill \end{array}$$

(4)

Solving equation (84), (85), and (4) simultaneously, we learn that

$$v=u.$$

(5)

From equation (5) $v=u$, using this in equation (85) we obtain the following

$$\begin{array}{cc}\hfill s& =a_2-1.\hfill \end{array}$$

(6)

From equation (6), we learn that $u=0$. Using this in equation (5) gives the following

$$v=u=0.$$

(7)

By substituting $v=0$ and $u=0$ as determined in equation (7) into equation (4), we obtain the following

$$s=0.$$

(8)

From equation (7) and (8) we learn that the origin is one the the equilibrium points of the system, that is

$${\mathbf{x}}_1=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right].$$

(9)

By substituting equation (6) $s=a_2-1$, and equation (5) $v=u$ into equation (4) we obtain the following equation

$$\begin{array}{cc}\hfill u^2& =a_3\left(a_2-1\right),\hfill \\ \hfill u& =\pm \sqrt{a_3\left(a_2-1\right)}.\hfill \end{array}$$

(10)

It follows from equation (10) that when $s=a_2-1$, $u=\pm \sqrt{a_3\left(a_2-1\right)}$ and from equation (5), $v=u=\pm \sqrt{a_3\left(a_2-1\right)}$. From this we learn that the Lorenz system has two more equilibrium points which are

$${\mathbf{x}}_2=\left[\begin{array}{c}\sqrt{a_3\left(a_2-1\right)}\\ \sqrt{a_3\left(a_2-1\right)}\\ a_2-1\end{array}\right],$$

and

$${\mathbf{x}}_3=\left[\begin{array}{c}-\sqrt{a_3\left(a_2-1\right)}\\ -\sqrt{a_3\left(a_2-1\right)}\\ a_2-1\end{array}\right].$$

By substituting the parameters $a_1=10$, $a_2=28$, and $a_3={\displaystyle \frac{8}{3}}$, the equilibrium points of the system are

$${\mathbf{x}}_1=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right],$$

(11)

$${\mathbf{x}}_2=\left[\begin{array}{c}6\sqrt{2}\\ 6\sqrt{2}\\ 27\end{array}\right],$$

(12)

$${\mathbf{x}}_3=\left[\begin{array}{c}-6\sqrt{2}\\ -6\sqrt{2}\\ 27\end{array}\right].$$

(13)

To analyse the stability of the system, a Jacobian matrix [27] of the form $J\left(\mathbf{x}\right)={\displaystyle \frac{\partial f\left(x\right)}{\partial x}}$ is constructed. Where

$$\begin{array}{cc}\hfill f_1(t,u,v,s)& =a_1(v-u),\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill f_2(t,u,v,s)& =-us+a_{2}u-v,\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill f_3(t,u,v,s)& =uv-a_{3}s.\hfill \end{array}$$

(16)

such that

$$J\left(\mathbf{x}\right)=\left[\begin{array}{ccc}{\displaystyle \frac{\partial f_1}{\partial u}}& {\displaystyle \frac{\partial f_1}{\partial v}}& {\displaystyle \frac{\partial f_1}{\partial s}}\\ {\displaystyle \frac{\partial f_2}{\partial u}}& {\displaystyle \frac{\partial f_2}{\partial v}}& {\displaystyle \frac{\partial f_2}{\partial s}}\\ {\displaystyle \frac{\partial f_3}{\partial u}}& {\displaystyle \frac{\partial f_3}{\partial v}}& {\displaystyle \frac{\partial f_3}{\partial s}}\end{array}\right],$$

(17)

is the Jacobian matrix corresponding to the system.

By finding the first partial derivatives of equation (14), (15), and (16) with respect to u, v and s, and substituting them into the Jacobian matrix (17) we obtain the following Jacobian matrix

$$J\left(\mathbf{x}\right)=\left[\begin{array}{ccc}-a_1& a_1& 0\\ -s+a_2& -1& -u\\ v& u& -a_3\end{array}\right].$$

(18)

Using the given parameters of $a_1$, $a_2$ and $a_3$ for the Lorenz system (1) the Jacobian matrix becomes

$$\begin{array}{c}\hfill J\left(\mathbf{x}\right)=\left[\begin{array}{ccc}-10& 10& 0\\ -s+28& -1& -u\\ v& u& -{\displaystyle \frac{8}{3}}\end{array}\right].\end{array}$$

(19)

To make conclusions about the stability of the system, the eigenvalues of the Jacobian matrix (19) are determined with respect to each equilibrium point [27], that is

$$det\left(J\left(\mathbf{x}\right)-\lambda I\right)=0.$$

(20)

For the equilibrium point ${\mathbf{x}}_1\left(11\right)$

$$J\left({\mathbf{x}}_1\right)=\left[\begin{array}{ccc}-10& 10& 0\\ 28& -1& 0\\ 0& 0& -{\displaystyle \frac{8}{3}}\end{array}\right],$$

$$\begin{array}{cc}\hfill det\left(J\left({\mathbf{x}}_1\right)-\lambda I\right)& =0,\hfill \\ \hfill -\left(\lambda +{\displaystyle \frac{8}{3}}\right)\left({\lambda }^2+11\lambda -270\right)& =0.\hfill \end{array}$$

(21)

By solving equation (21), we learn that the eigenvalues are

$$\begin{array}{cc}\hfill {\lambda }_1& =-{\displaystyle \frac{8}{3}},\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\lambda }_2& =-{\displaystyle \frac{\sqrt{1201}-11}{2}},\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\lambda }_3& ={\displaystyle \frac{\sqrt{1201}-11}{2}}.\hfill \end{array}$$

(24)

${\lambda }_3$ is greater than zero, this implies that the equilibrium point ${\mathbf{x}}_1$ is unstable.

Subsequently, the same procedure is followed to determine the eigenvalues for the other equilibrium points ${\mathbf{x}}_2$ and ${\mathbf{x}}_3$.

For the equilibrium point ${\mathbf{x}}_2$ (12), we obtain the following Jacobian matrix and eigenvalues.

$$J\left({\mathbf{x}}_2\right)=\left[\begin{array}{ccc}-10& 10& 0\\ 1& -1& -6\sqrt{2}\\ 6\sqrt{2}& 6\sqrt{2}& -{\displaystyle \frac{8}{3}}\end{array}\right],$$

$$\begin{array}{cc}\hfill det\left(J\left({\mathbf{x}}_2\right)-\lambda I\right)& =0,\hfill \\ \hfill \left|\begin{array}{ccc}-10-\lambda & 10& 0,\\ 1& -1-\lambda & -6\sqrt{2},\\ 6\sqrt{2}& 6\sqrt{2}& -{\displaystyle \frac{8}{3}}-\lambda \end{array}\right|& =0,\hfill \\ \hfill -{\lambda }^3-{\displaystyle \frac{41}{3}}{\lambda }^2-{\displaystyle \frac{304}{3}}\lambda -1440& =0.\hfill \end{array}$$

(25)

Solving equation (25) using the quadratic formula we find the following eigenvalues

$$\begin{array}{cc}\hfill {\lambda }_1& =-13.8546,\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill {\lambda }_2& =0.09396+10.1945i,\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\lambda }_3& =0.09396-10.1945i.\hfill \end{array}$$

(28)

The real parts of ${\lambda }_2$ and ${\lambda }_3$ are all greater than zero, this suggests that ${\mathbf{x}}_2$ is an unstable equilibrium point.

For the equilibrium point ${\mathbf{x}}_3$ (13), we obtain

$$J\left({\mathbf{x}}_3\right)=\left[\begin{array}{ccc}-10& 10& 0\\ 1& -1& 6\sqrt{2}\\ -6\sqrt{2}& -6\sqrt{2}& -{\displaystyle \frac{8}{3}}\end{array}\right],$$

$$\begin{array}{cc}\hfill det\left(J\left({\mathbf{x}}_3\right)-\lambda I\right)& =0,\hfill \\ \hfill -{\lambda }^3-{\displaystyle \frac{41}{3}}{\lambda }^2-{\displaystyle \frac{304}{3}}\lambda -1440& =0,\hfill \end{array}$$

(29)

Solving equation (29) using the quadratic formula we find the following eigenvalues

$$\begin{array}{cc}\hfill {\lambda }_1& =-13.8546,\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\lambda }_2& =0.09396+10.1945i,\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\lambda }_3& =0.09396-10.1945i.\hfill \end{array}$$

(32)

The real parts of ${\lambda }_2$ and ${\lambda }_3$ are all greater than zero, this suggests that ${\mathbf{x}}_3$ is an unstable equilibrium point.

For the given values of the parameters of $a_1$, $a_2$ and $a_3$, none of the equilibrium points of the Lorenz system are stable. This implies that the system is unstable at all its equilibrium points.

4.2. Genesio-Tesi System

This system was proposed by Roberto Genesio and Alberto Tesi [28]. It is given by the following differential equations

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\displaystyle \frac{du}{dt}}\hfill & =v,\hfill \\ {\displaystyle \frac{dv}{dt}}\hfill & =s,\hfill \\ {\displaystyle \frac{ds}{dt}}\hfill & =a_{1}s+a_{2}v+a_{3}u+u^2,\hfill \end{array}\right.\end{array}$$

(33)

where the parameters $a_1=-1.2$, $a_2=-2.92$, and $a_3=-6$ with initial conditions $u\left(0\right)=0.2$, $v\left(0\right)=-0.3$ and $s\left(0\right)=0.1$. It is said that the three ordinary differential equations of this system depend on real negative parameters of $a_1$, $a_2$, and $a_3$ [29].

4.2.1. Stability Analysis of the Genesio Tesi System

Given the Genesio Tesi System (33), we set ${\displaystyle \frac{du}{dt}}=0$, ${\displaystyle \frac{dv}{dt}}=0$, and ${\displaystyle \frac{ds}{dt}}=0$ to determine the equilibrium points of the system. This leads to the following the equations

$$\begin{array}{cc}\hfill v& =0,\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill s& =0,\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill a_{1}s+a_{2}v+a_{3}u+u^2& =0.\hfill \end{array}$$

(36)

We learn from equation (34) and (35) that v=0 and s=0, using this in equation (36) we obtain

$$\begin{array}{cc}\hfill u& =-a_3.\hfill \end{array}$$

(37)

From equation () and (37), it follows that the system (33) has two equilibrium points which are

$${\mathbf{x}}_1=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right],$$

(38)

$${\mathbf{x}}_2=\left[\begin{array}{c}0\\ 0\\ -a_3\end{array}\right],$$

(39)

For the stated parameters of $a_1$, $a_2$, and $a_3$ of the system (33), we learn that the equilibrium point (39) is

$${\mathbf{x}}_2=\left[\begin{array}{c}0\\ 0\\ 6\end{array}\right],$$

(40)

By setting

$$\begin{array}{cc}\hfill f_1(t,u,v,s)& =v,\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill f_2(t,u,v,s)& =s,\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill f_3(t,u,v,s)& =a_{1}s+a_{2}v+a_{3}u+u^2,\hfill \end{array}$$

(43)

and finding the first partial derivatives with respect to u, v, s of equation (41), (42) and (43) to construct a Jacobian matrix (17) corresponding to the system , we obtain

$$J\left(\mathbf{x}\right)=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ a_3+2u& a_2& a_1\end{array}\right].$$

(44)

For the chosen parameters of $a_1=-1.2$, $a_2=-2.92$, and $a_3=-6$ for the system (33), the Jacobian matrix (44) becomes

$$J\left(\mathbf{x}\right)=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -6+2u& -2.92& -1.2\end{array}\right].$$

(45)

The eigenvalues of the Jacobian matrix (45) with respect to the equilibrium point (38) are given by

$$det(J\left({\mathbf{x}}_1\right)-\lambda I)=0,$$

which leads to the following equation.

$$\begin{array}{cc}\hfill \left|\begin{array}{ccc}-\lambda & 1& 0\\ 0& -\lambda & 1\\ -6& -2.92& -1.2-\lambda \end{array}\right|& =0,\hfill \\ \hfill -{\lambda }^3-1.2{\lambda }^2-2.92\lambda -6& =0.\hfill \end{array}$$

(46)

Solving equation (46) using the quadratic formula we find the following eigenvalues

$$\begin{array}{cc}\hfill {\lambda }_1& =0.22197+1.8975i,\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill {\lambda }_2& =0.22197-1.8975i,\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill {\lambda }_3& =-1.64393.\hfill \end{array}$$

(49)

The real parts of ${\lambda }_1$ and ${\lambda }_2$ are both greater than zero, this means that the equilibrium point ${\mathbf{x}}_1$ (38) is unstable.

The eigenvalues corresponding to the Jacobian matrix (45) with respect to the equilibrium point ${\mathbf{x}}_2$ (40) are given by

$$det(J\left({\mathbf{x}}_2\right)-\lambda I)=0.$$

(50)

Determining the eigenvalues from equation (50) we obtain

$$\begin{array}{cc}\hfill \left|\begin{array}{ccc}-\lambda & 1& 0\\ 0& -\lambda & 1\\ 6& -2.92& -1.2-\lambda \end{array}\right|& =0,\hfill \\ \hfill -{\lambda }^3-1.2{\lambda }^2-2.92\lambda +6& =0.\hfill \end{array}$$

(51)

Solving equation (51) using the quadratic formula we find the following eigenvalues

$$\begin{array}{cc}\hfill {\lambda }_1& =-1.15027+2.03193i,\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill {\lambda }_2& =-1.15027-2.03193i,\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill {\lambda }_3& =1.10054.\hfill \end{array}$$

(54)

We learn that ${\lambda }_3$ is greater than zero, this implies that the equilibrium point ${\mathbf{x}}_2$ of the system is also unstable.

The Genesio Tesi system (33) is unstable at both of its equilibrium points. This means that for the stated parameters of $a_1$, $a_2$ and $a_3$, the system is unstable.

4.3. Rossler System

The system was studied by Otto Rossler [2]. It is represented by the following set of ordinary differential equations

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\displaystyle \frac{du}{dt}}\hfill & =-(v+s),\hfill \\ {\displaystyle \frac{dv}{dt}}\hfill & =u+a_{1}v,\hfill \\ {\displaystyle \frac{ds}{dt}}\hfill & =a_2+us-a_{3}s,\hfill \end{array}\right.\end{array}$$

(55)

where the parameters $a_1$, $a_2$, and $a_3$ are greater than zero. The system (55) defines a continuous-time dynamical system that exhibits chaotic dynamics with the fractal properties of the attractor for selected values of $a_1$, $a_2$ and $a_3$ [17]. Rossler observed a chaotic behavior when conducting the study [30] with $a_1=0.2$, $a_2=0.2$, and $a_3=5.7$ subject to the initial conditions $u\left(0\right)=0$, $v\left(0\right)=0$, and $s\left(0\right)=0$.

4.3.1. Stability Analysis of the Rossler System

Given the Rossler system (55), by setting ${\displaystyle \frac{du}{dt}}=0$, ${\displaystyle \frac{dv}{dt}}=0$, and ${\displaystyle \frac{ds}{dt}}=0$ to find the equilibrium points of the system, we obtain the following equations

$$\begin{array}{cc}\hfill -(v+s)& =0,\hfill \end{array}$$

(56)

$$\begin{array}{cc}\hfill u+a_{1}v& =0,\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill a_2+us-a_{3}s& =0.\hfill \end{array}$$

(58)

Solving the equations (56)-(58) simultaneously, we find that

$$v=-s.$$

(59)

from equation (56).

By substituting $v=-s$ from equation (59) into equation (57) we obtain

$$u=a_{1}s.$$

(60)

From equation (60) $u=a_{1}s$, using this in equation (58), it follows that

$$\begin{array}{cc}\hfill a_1{s}^2-a_{3}s+a_2& =0,\hfill \\ \hfill s& ={\displaystyle \frac{a_3\pm \sqrt{a_3^2-4a_1{a}_2}}{2a_1}},\hfill \end{array}$$

(61)

By substituting $s={\displaystyle \frac{a_3\pm \sqrt{a_3^2-4a_1{a}_2}}{2a_1}}$ from equation (61) into equation (60) and (61), we obtain

$$\begin{array}{cc}\hfill u& =a_1\left({\displaystyle \frac{a_3\pm \sqrt{a_3^2-4a_1{a}_2}}{2a_1}}\right).\hfill \end{array}$$

(62)

$$\begin{array}{cc}\hfill v& =-\left({\displaystyle \frac{a_3\pm \sqrt{a_3^2-4a_1{a}_2}}{2a_1}}\right).\hfill \end{array}$$

(63)

Equation (61), (62), and (63) suggest that the system has two equilibrium points which are

$$\begin{array}{cc}\hfill {\mathbf{x}}_1& =\left[\begin{array}{c}{a}_1\left({\displaystyle \frac{a_3+\sqrt{a_3^2-4a_1{a}_2}}{2a_1}}\right)\\ -\left({\displaystyle \frac{a_3+\sqrt{a_3^2-4a_1{a}_2}}{2a_1}}\right)\\ {\displaystyle \frac{a_3+\sqrt{a_3^2-4a_1{a}_2}}{2a_1}}\end{array}\right].\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill {\mathbf{x}}_2& =\left[\begin{array}{c}{a}_1\left({\displaystyle \frac{a_3-\sqrt{a_3^2-4a_1{a}_2}}{2a_1}}\right)\\ -\left({\displaystyle \frac{a_3-\sqrt{a_3^2-4a_1{a}_2}}{2a_1}}\right)\\ {\displaystyle \frac{a_3-\sqrt{a_3^2-4a_1{a}_2}}{2a_1}}\end{array}\right].\hfill \end{array}$$

(65)

For the stated parameters of $a_1$, $a_2$, and $a_3$ of the system (55), the equilibrium points ${\mathbf{x}}_1$ (64) and ${\mathbf{x}}_2$ (65) of the system are

$$\begin{array}{cc}\hfill {\mathbf{x}}_1& =\left[\begin{array}{c}5.69297\\ -28.46487\\ 28.46487\end{array}\right].\hfill \end{array}$$

(66)

$$\begin{array}{cc}\hfill {\mathbf{x}}_2& =\left[\begin{array}{c}0.00703\\ -0.03513\\ 0.03513\end{array}\right].\hfill \end{array}$$

(67)

By letting

$$\begin{array}{cc}\hfill f_1(t,u,v,s)& =-(v+s),\hfill \end{array}$$

(68)

$$\begin{array}{cc}\hfill f_2(t,u,v,s)& =u+a_{1}v,\hfill \end{array}$$

(69)

$$\begin{array}{cc}\hfill f_3(t,u,v,s)& =a_2+us-a_{3}s,\hfill \end{array}$$

(70)

and finding the first partial derivatives of (68), (69) and (70) with respect to u, v, s to construct a Jacobian matrix (17) corresponding to the system we obtain

$$J\left(\mathbf{x}\right)=\left[\begin{array}{ccc}0& -1& -1\\ 1& 0.2& 0\\ s& 0& u-5.7\end{array}\right].$$

(71)

The eigenvalues of the Jacobian matrix with respect to the equilibrium point ${\mathbf{x}}_1$ (66) are determined by evaluating

$$det\left(J\left({\mathbf{x}}_1\right])-\lambda I\right)=0.$$

(72)

Evaluating equation (72) gives

$$\begin{array}{cc}\hfill \left|\begin{array}{ccc}-\lambda & -1& -1\\ 1& 0.2-\lambda & 0\\ 28.46487& 0& -0.00703-\lambda \end{array}\right|& =0,\hfill \\ \hfill -{\lambda }^3+0.19297{\lambda }^2-29.4635\lambda +5.6894& =0.\hfill \end{array}$$

(73)

Solving equation (73) using the quadratic formula we find the following eigenvalues

$$\begin{array}{cc}\hfill {\lambda }_1& =6.42954\times {10}^{-6}+5.42803i,\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill {\lambda }_2& =-6.42954\times {10}^{-6}-5.42803i,\hfill \end{array}$$

(75)

$$\begin{array}{cc}\hfill {\lambda }_3& =0.19298.\hfill \end{array}$$

(76)

The real parts of ${\lambda }_1$ and ${\lambda }_3$ are greater than zero, this indicates that the equilibrium point ${\mathbf{x}}_1$ is unstable.

The eigenvalues of the Jacobian matrix (71) with respect to the equilibrium point ${\mathbf{x}}_2$ (67). are determined by evaluating

$$det\left(J\left({\mathbf{x}}_2\right)-\lambda I\right)=0.$$

(77)

Evaluating (77) gives

$$\begin{array}{cc}\hfill \left|\begin{array}{ccc}-\lambda & -1& -1\\ 1& 0.2-\lambda & 0\\ 0.03513& 0& -5.69297-\lambda \end{array}\right|& =0,\hfill \\ \hfill -{\lambda }^3-5.49297{\lambda }^2+0.103464\lambda +5.68594& =0.\hfill \end{array}$$

(78)

Solving equation (78) using the quadratic formula we find the following eigenvalues

$$\begin{array}{cc}\hfill {\lambda }_1& =-5.68697,\hfill \end{array}$$

(79)

$$\begin{array}{cc}\hfill {\lambda }_2& =0.09700+0.99519i,\hfill \end{array}$$

(80)

$$\begin{array}{cc}\hfill {\lambda }_3& =0.09700-0.99519i.\hfill \end{array}$$

(81)

The real parts of ${\lambda }_2$ and ${\lambda }_3$ are greater than zero, this indicates that the equilibrium point ${\mathbf{x}}_2$ is also unstable.

The Rossler system (55) is unstable at all its equilibrium points for the stated parameters of $a_1$, $a_2$, and $a_3$.

4.4. Derivation of the Runge-Kutta Methods of Order Two and Four

The second order Runge-Kutta method forms part of the family of explicit iterative Runge-Kutta methods. It is used to solve ordinary differential equations of the form ${\displaystyle \frac{dy}{dx}}=f(x,y)$ given an initial condition $y\left(a\right)=\alpha$. It provides the approximate value of y for a given point x. It is considered to be more effective than the other explicit Runge-Kutta methods as its local truncation error is of order 5 (i.e., $O\left(h^2\right)$) [31].

The general s-stage Runge–Kutta recurrence formula is

$$y_{n+1}=y_n+h{\displaystyle \sum _{i=1}^s}{a}_i{k}_i,$$

(82)

where the stage values $k_i$ are defined by

$$k_i=f\left(t_n+p_{i}h,y_n+h{\displaystyle \sum _{j=1}^s}{q}_{ij}{k}_j\right),i=1,2,\dots ,s.$$

(83)

The coefficients $a_i$ ,$q_{ij}$, and $p_i$ define the specific Runge-Kutta method.

$$\begin{array}{c}\hfill y_{i+1}=y_i+h(a_1{k}_1+a_2{k}_2),\end{array}$$

(84)

with

$$\begin{array}{cc}\hfill k_1& =f(t_i,y_i),\hfill \\ \hfill k_2& =f(t_i+p_{1}h,y_i+q_{11}h{k}_1).\hfill \end{array}$$

(85)

The Taylor series expansion of $y(t_i+h)$ about $t_i$ is given as follows

$$\begin{array}{c}\hfill y(t_i+h)=y\left(t_i\right)+h{y}^{\prime }\left(t_i\right)+{\displaystyle \frac{h^2}{2}}{y}^{″}\left(t_i\right)+{\displaystyle \frac{h^3}{3!}}{y}^{‴}\left(t_i\right)+O\left(h^4\right).\end{array}$$

(86)

where

$$\begin{array}{cc}\hfill y^{\prime }& =f(t_i,y_i),\hfill \\ \hfill y^{″}& =f_t+f_{y}f.\hfill \end{array}$$

(87)

By substituting equation 87 into 86 we get the following

$$\begin{array}{cc}\hfill y(t_i+h)& =y\left(t_i\right)+hf+{\displaystyle \frac{h^2}{2}}\left(f_t+f_{y}f\right)+O\left(h^3\right).\hfill \end{array}$$

(88)

By substituting equation 85 into 84 with the use of Taylor expansion, we obtain the following equation

$$\begin{array}{cc}\hfill y_{i+1}& =y_i+h(a_1+a_2)f+{\displaystyle \frac{h^2}{2!}}\left(a_2(p_1{f}_t+q_{11}{f}_{y}f)\right).\hfill \end{array}$$

(89)

By comparing the coefficients on equation 88 and 89, the following system of equations is obtained

$$\begin{array}{c}\hfill a_1+a_2=1,\\ \hfill p_1{a}_2={\displaystyle \frac{1}{2}},\\ \hfill q_{11}{a}_2={\displaystyle \frac{1}{2}}.\end{array}$$

(90)

Since we have 3 equations and 4 unknowns, we can assume the value of one of the unknowns. The other three will then be determined from the three equations. Generally the value of $a_2$ is chosen to evaluate the other three constants. The three values generally used for $a_2$ are ${\displaystyle \frac{1}{2}}$, 1 and ${\displaystyle \frac{2}{3}}$, and are known as Heun’s Method, the midpoint method and Ralston’s method, respectively. The Improved Euler’s(Heun’s method) is given as follow:

$$\begin{array}{c}\hfill y_{i+1}=y_i+{\displaystyle \frac{h}{2}}\left(k_1+k_2\right),\end{array}$$

(91)

where

$$k_1=f(t_i,y_i),$$

(92)

$$k_2=f(t_i+h,y_i+h{k}_1).$$

(93)

The Midpoint is given as follow:

$$\begin{array}{c}\hfill y_{i+1}=y_i+h{k}_2,\end{array}$$

(94)

where

$$k_1=f(t_i,y_i),$$

(95)

$$k_2=f\left(t_i+{\displaystyle \frac{1}{2}}h,y_i+{\displaystyle \frac{1}{2}}h{k}_1\right).$$

(96)

The Ralston’s method is given as follow:

$$\begin{array}{c}\hfill y_{i+1}=y_i+h\left({\displaystyle \frac{1}{3}}{k}_1+{\displaystyle \frac{2}{3}}{k}_2\right),\end{array}$$

(97)

where

$$k_1=f(t_i,y_i),$$

(98)

$$k_2=f\left(t_i+{\displaystyle \frac{3}{4}}h,y_i+{\displaystyle \frac{3}{4}}h{k}_1\right).$$

(99)

4.5. Derivation of Runge-Kutta Method of Order Four

The fourth order Runge-Kutta method forms part of the family of explicit iterative Runge-Kutta methods. It is used to solve ordinary differential equations of the form ${\displaystyle \frac{dy}{dx}}=f(x,y)$ given an initial condition $y\left(a\right)=\alpha$. It provides the approximate value of y for a given point x. It is considered to be more effective than the other explicit Runge-Kutta methods as its local truncation error is of order 5 (i.e., $O\left(h^5\right)$) [31].

In general, the Runge-Kutta method of order is given by

$$y_{t+h}=y_t+h{\displaystyle \sum _{i=1}^s}{b}_i{k}_i,$$

(100)

where

$$k_i=f\left(t_n+c_{i}h,y_n+h{\displaystyle \sum _{j=1}^s}{a}_{ij}kj\right),\mathrm{for}\mathrm{i}=1,..,\mathrm{s}\&\mathrm{n}=0,...,\mathrm{N}$$

(101)

where the fixed scalars $a_{ij}$, $b_i$ and $c_i$ are the coefficients of the Runge Kutta formula.

In order to derive the 4th order Runge-Kutta method, we obtain the following recurrence equation for the 4th order Runge-Kutta method from (100)

$$\begin{array}{c}\hfill y_{i+1}=y_i+h(a_1{k}_1+a_2{k}_2+a_3{k}_3+a_4{k}_4),\end{array}$$

(102)

 

(103)

where,

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{k}_1\hfill & =f\left(t_i,y_i\right),\hfill \\ k_2\hfill & =f\left(t\left(i\right)+{\displaystyle \frac{h}{2}},y_i+{\displaystyle \frac{h}{2}}{k}_1\right),\hfill \\ k_3\hfill & =f\left(t\left(i\right)+{\displaystyle \frac{h}{2}},y_i+{\displaystyle \frac{h}{2}}{k}_2\right),\hfill \\ k_4\hfill & =f\left(t\left(i\right)+h,y_i+h{k}_3\right).\hfill \end{array}\right.\end{array}$$

(104)

Therefore, by using the Taylor method to expand the terms (104) of the Runge Kutta method (102) we obtain the following

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{k}_1\hfill & =f(t,y),\hfill \\ k_2\hfill & =f\left(t+{\displaystyle \frac{h}{2}},y+{\displaystyle \frac{h}{2}}\right),\hfill \\ \hfill & =f+{\displaystyle \frac{1}{2}}h{f}_t,\hfill \\ k_3\hfill & =f\left(t+{\displaystyle \frac{h}{2}},y+{\displaystyle \frac{h}{2}}\right),\hfill \\ \hfill & =f\left(t+{\displaystyle \frac{h}{2}},y+{\displaystyle \frac{h}{2}}f(y+{\displaystyle \frac{h}{2}}{k}_1)\right),\hfill \\ \hfill & =f+{\displaystyle \frac{h}{2}}\left(f_t+{\displaystyle \frac{h}{2}}{f}_{tt}\right)\hfill \\ k_4\hfill & =f(t+h,y+h{k}_3),\hfill \\ \hfill & =f+h{f}_t+{\displaystyle \frac{1}{2}}{h}^2{f}_{tt}+{\displaystyle \frac{1}{4}}{h}^3{f}_{ttt}.\hfill \end{array}\right.\end{array}$$

(105)

By substituting (105) into equation (102) gives the following equation

$$\begin{array}{cc}\hfill y_{i+1}=& y_i+h\left[\left(a_{1}f+a_2(f+{\displaystyle \frac{h}{2}}{f}_t)\right)+a_3\left(f+{\displaystyle \frac{h}{2}}(f_t+{\displaystyle \frac{h}{2}}{f}_{tt})\right)+\right.\hfill \\ \hfill & \left.a_4\left(f+h{f}_t+{\displaystyle \frac{h^2}{2}}{f}_{tt}+{\displaystyle \frac{h^3}{2}}{f}_{ttt})\right)\right],\hfill \\ \hfill =& y_i+hf\left(a_1+a_2+a_3+a_4\right)+h^2\left({\displaystyle \frac{1}{2}}{a}_3{f}_t+{\displaystyle \frac{1}{2}}{a}_3{f}_t+a{f}_t\right)+\hfill \\ \hfill & h^3\left({\displaystyle \frac{1}{4}}{a}_3{f}_{ttt}+{\displaystyle \frac{1}{2}}{a}_4{f}_{tt}\right)+h^4\left({\displaystyle \frac{1}{4}}{a}_4{f}_{tt}\right).\hfill \end{array}$$

(106)

Using the Taylor method to expand equation (106) on the right gives

$$\begin{array}{cc}\hfill y_{i+1}=& y_i+h{y}^{\prime }+{\displaystyle \frac{h^2}{2!}}{y}^{″}+{\displaystyle \frac{h^3}{3!}}{y}^{‴}+{\displaystyle \frac{h^4}{4!}}{y}^{″″}+O\left(h^5\right),\hfill \\ \hfill =& y_i+hf+{\displaystyle \frac{h^2}{2}}\left(f_t+f{f}_y\right)+{\displaystyle \frac{h^3}{6}}\left(f_{tt}+2f{f}_{ty}+f^2{f}_{yy}+f_y{f}_t+f{\left(f_y\right)}^2\right)+\hfill \\ \hfill & {\displaystyle \frac{h^4}{24}}\left(f_{ttt}+3f_{ty}{f}_t+2f_{ty}\right).\hfill \end{array}$$

(107)

Comparing the coefficients of equation (107) and equation (106) we obtain the following system of equations

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}1\hfill & =a_1+a_2+a_3+a_4,\hfill \\ {\displaystyle \frac{1}{2}}\hfill & ={\displaystyle \frac{1}{2}}{a}_2+{\displaystyle \frac{1}{2}}{a}_3+a_4,\hfill \\ {\displaystyle \frac{1}{6}}\hfill & ={\displaystyle \frac{1}{4}}{a}_3+{\displaystyle \frac{1}{4}}{a}_4,\hfill \\ {\displaystyle \frac{1}{24}}\hfill & ={\displaystyle \frac{1}{4}}{a}_4.\hfill \end{array}\right.\end{array}$$

(108)

Solving the above system of equation gives $a_1={\displaystyle \frac{1}{6}}$, $a_2={\displaystyle \frac{1}{3}}$, $a_3={\displaystyle \frac{1}{3}}$, and $a_4={\displaystyle \frac{1}{6}}$. It then follows that equation (102), which is the Runge Kutta method of order 4 is

$$\begin{array}{c}\hfill y_{i+1}=y_i+{\displaystyle \frac{h}{6}}(k_1+2k_2+2k_3+k_4).\end{array}$$

(109)

5. Algorithm Addition

Algorithm 1: Numerical Solution Procedure

1.

Define chaotic system parameters

2.

Specify initial conditions

3.

Select step size h

4.

Implement Midpoint method

5.

Implement Improved Euler method

6.

Implement Ralston method

7.

Implement RK4 method

8.

Compute numerical solutions

9.

Determine infinity norm errors

10.

Evaluate convergence rates

11.

Compute computational time

12.

Compare numerical results

13.

Plot trajectories and phase portraits

6. Results and Discussion

MATLAB implementation details of Runge–Kutta methods.

For the Lorenz system in Table 1, all second-order methods (Midpoint, Improved Euler, and Ralston) exhibit relatively large errors for coarse step sizes. As the step size decreases from $h=0.01000$ to $h=0.00063$, the errors reduce significantly, confirming convergence of the methods. Among the second-order schemes, the Improved Euler method initially produces slightly smaller errors than the Midpoint and Ralston methods, while the Ralston method shows competitive performance for smaller step sizes. However, the RK4 method clearly outperforms all second-order methods, producing errors several orders of magnitude smaller across all step sizes. For example, at h=0.00500, RK4 attains an error of $2.12\times {10}^{-3}$, whereas the second-order methods still produce errors of order ${10}^0$ to ${10}^1$ . This behaviour confirms the superior accuracy and stability of the fourth-order Runge–Kutta scheme for highly sensitive chaotic systems such as the Lorenz equations. Table 2 presents the convergence results for the Genesio system. The Midpoint, Improved Euler, and Ralston methods display nearly identical error values for all step sizes, indicating that these second-order methods possess comparable numerical behaviour for this problem. The errors decrease consistently as the step size is refined, which verifies the expected convergence property of the methods. The RK4 method again provides substantially smaller errors than the second-order methods. Interestingly, the RK4 error remains almost constant around $3.80\times {10}^{-5}$ even as the step size decreases. This suggests that the method has already reached a very high level of accuracy, where further refinement of the step size contributes only marginal improvement, possibly due to machine precision limitations or accumulated round-off effects.

Similarly, the results for the Rossler system in Table 3 indicate strong convergence properties for all methods. The Midpoint, Improved Euler, and Ralston methods again produce identical numerical errors, with the error decreasing rapidly when the step size is reduced from $0.01000$ to $0.00500$. However, beyond this point, the reduction becomes less pronounced, indicating that the solutions are approaching a limiting numerical accuracy. The RK4 method consistently achieves the smallest errors, approximately $3.05\times {10}^{-6}$, and remains highly stable for all step sizes considered. The near-constant RK4 error for smaller step sizes further demonstrates the robustness and efficiency of the fourth-order method for solving chaotic dynamical systems.

The numerical results presented in Table 1, Table 2 and Table 3 illustrate the convergence behaviour of the Midpoint, Improved Euler, Ralston, and RK4 methods when applied to the Lorenz, Genesio, and Rossler chaotic systems. In general, the results demonstrate that decreasing the step size h improves the numerical accuracy of all methods, although the rate of improvement differs among the schemes and dynamical systems considered. The numerical experiments confirm that all methods converge as the step size decreases, but the RK4 method provides the highest accuracy and best stability characteristics among the methods investigated. The second-order methods are computationally simpler and still produce acceptable approximations for sufficiently small step sizes, but their errors remain significantly larger than those of RK4, particularly for strongly nonlinear and chaotic systems such as the Lorenz model. These findings demonstrate the advantage of higher-order numerical integration techniques when solving sensitive dynamical systems where accuracy is critical.

The results presented in Table 4 illustrate that the computational results show that the performance of the numerical methods varies depending on the chaotic system being solved. RK4 achieved the best computational efficiency for the Lorenz system, while the Ralston method performed best for the Genesio–Tesi system, and the Midpoint method was most efficient for the Rössler system. Although RK4 is a higher-order method, it did not consistently produce the lowest computational cost across all systems. Overall, the findings indicate that computational efficiency is problem dependent, and lower-order methods can sometimes outperform higher-order methods for certain chaotic systems.

The infinity norm plots shown in Figure 1, Figure 2 and Figure 3 further validate the theoretical orders of convergence of the methods. The slopes corresponding to the second-order methods align closely with $O\left(h^2\right)$, while the RK4 method follows the expected $O\left(h^4\right)$ convergence trend. This agreement between theoretical and numerical convergence confirms the correctness of the implementation of the methods.

The time evolution plots and phase portraits presented in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19 illustrate the qualitative behaviour of the chaotic systems. The numerical solutions generated by all methods generally followed similar trajectories, especially for smaller step sizes. However, the RK4 solutions remained smoother and more stable over longer integration intervals, whereas slight deviations became visible in the second-order methods due to accumulated truncation errors. The phase portraits also successfully captured the chaotic attractors associated with the Lorenz, Genesio–Tesi, and Rossler systems, demonstrating that the numerical schemes were capable of reproducing the essential dynamics of the systems.

The results indicate that while the Midpoint, Improved Euler, and Ralston methods are computationally simpler and easier to implement, they are less accurate for highly sensitive chaotic systems. In contrast, the RK4 method provides significantly better numerical accuracy, stronger stability properties, and more reliable long-term approximations. These findings highlight the importance of using higher-order numerical integration techniques when solving nonlinear chaotic differential equations where precision and stability are critical.

The superior performance of RK4 can be attributed to its higher-order truncation accuracy and multiple intermediate evaluations within each time step. Chaotic systems are highly sensitive to small perturbations, and numerical errors propagate rapidly over time. The RK4 method effectively reduces the accumulation of these local errors, thereby providing more stable long-term approximations than second-order methods.

Future work may extend this analysis to adaptive Runge–Kutta methods, fractional chaotic systems, and hybrid machine-learning-assisted numerical solvers.

7. Conclusions

This study investigated the efficiency and effectiveness of the Midpoint, Improved Euler, Ralston, and fourth-order Runge Kutta numerical methods when solving the Lorenz, Genesio Tesi, and Rossler chaotic systems. The stability analysis showed that all selected systems are unstable at their equilibrium points, confirming their chaotic nature and sensitivity to initial conditions.

The numerical experiments demonstrated that all methods converge as the step size decreases. However, the RK4 method consistently produced the smallest errors and exhibited superior stability and convergence properties compared to the second-order methods. The Midpoint, Improved Euler, and Ralston methods produced acceptable approximations for sufficiently small step sizes, but their errors remained considerably larger than those of RK4, especially for the Lorenz system where chaotic effects are more pronounced.

The infinity norm analysis and graphical comparisons confirmed the theoretical convergence orders of the methods, with the second-order schemes following $O\left(h^2\right)$ convergence and RK4 following $O\left(h^4\right)$. Furthermore, the phase portraits and solution trajectories successfully reproduced the chaotic dynamics of the systems under consideration.

Based on the obtained results, it can be concluded that the fourth order Runge Kutta method is the most efficient and reliable method among the investigated schemes for solving chaotic systems of ordinary differential equations. Its high accuracy, improved stability, and robustness make it highly suitable for nonlinear and highly sensitive dynamical systems. Future studies may extend this work by investigating adaptive step-size Runge Kutta methods, implicit schemes, or machine learning-based numerical approaches for solving more complex chaotic and hyperchaotic systems.