Numerical Investigation on System of Ordinary Differential Equations Absolute Time Inference with Mathematica®

Abstract—The purpose of this research is to perform a comparative numerical analysis of an efficient numerical methods for second-order ordinary differential equations, by reducing the second-order ODE to a system of first-order differential equations. Then we obtain approximate solutions to the system of ODE. To validate the accuracy of the algorithm, a comparison between Euler’s method and the Runge-Kutta method or order four was carried out and an exact solution was found to verify the efficiency, accuracy of the methods. Graphical representations of the parametric plots were also presented. Time inference analysis is taken to check the time taken to executes the algorithm in Mathematica®12.2.0. The obtained approximate solution using the algorithm shows that the Runge-Kutta method of order four is more efficient for solving system of linear ordinary differential equations.


I. INTRODUCTION
In recent times, computer algebraic system (CAS) have been employed to investigate the use of built-in functions and construction of algorithm to obtain solutions to initial value problems. Scientific problems arises in stem fields as Biology, Physics, Chemistry ,and Engineering [1], Ecology with respect to inverse problems [2]. They also arise in nonstem fields as humanities, virtual art, designs and films [3], [4]. Differential equations recently have been understood, that its application model life realities and numerical methods are used to investigate ODEs. Some problems have been modeled to explore dynamics of music, literature, poetry, prey-predator models, decay radiation, and numerical methods [5], have been used to understand its dynamics in terms of approximation. Some mathematical real life situations modeled by system of ODEs do not have exact solutions, hence, approximation and/or numerical techniques must be used. In this paper, we consider second-order ODEs, convert them to a system of ODEs and compare the exact and numerical solution. This exploration is investigated through evaluation using Mathematica® [6], [7], [8]. Euler's and Runge-Kutta [9], [5], [10] algorithm are implemented with its built-in function. Other researchers have implemented several methods to solve initial value problems and the analysis of the accuracy in areas such as Epidemiology [11], stability and efficiency to system of ODEs [12] but not with the algorithm implemented in this paper and the use of time inference obtained with Mathematica®. The Euler's method has large errors which is illustrated using Taylor's series expansion. Taylor's series converges within a small range and as a result, if the step size is not relatively small, it diverges. From Taylor's expansion, the first two terms represent the Euler's method. Through this, its imperative enough to know if h is not small enough, the method is not accurate and will not be a good use for practical implementation. The Runge-Kutta method of order four gives a better approximation than the Euler's method but its complicated in computing. However, it converges faster to the exact solution. This procedure explore the comparative analysis between the exact and numerical solution of the ODEs [13], [14]. The algorithm works robustly for obtaining solution to system of first order ODEs [15] and comparing the solution by computation and investigating the absolute timing it takes for each algorithm to compute, hence the algorithm can be applied to biological models. This algorithm will be applied to investigate nonlinear system of ODEs for further studies.
In this paper, we are more concerned about the accuracy of the algorithm and its reliability after comparing the exact solution, Euler's method and Runge-Kutta method on the system of first-order linear ODEs.

II. METHODOLOGY
In this section, we will consider a second-order ODEs by recasting it to system of first-order ODEs. In considering the Euler's method and Runge-Kutta method of order four, this methods investigated by coding an algorithm using the CAS, Mathematica and its built-in function. The approximate solution of both methods would be obtained and compared by the result of algorithm and the built-in function of Math-ematica®. For Equation (1) and Equation (2), the explicit formula for the solution can't be obtained for an initial value problem of such form but can approximate the solution using the numerical methods which is based on the tangent line approximation. In this paper, the Euler's method would be used to approximate the system of ODEs with initial conditions. This technique will investigate numerically with the aid of the Mathematica® software.

III. EULER'S METHOD
Consider the initial value problem in with initial condition x(0) = 0, and y(0) = 0 valid for a ≤ t ≤ b Let Such that where ∆t = t n − t 0 n is known as the step size. The rule of the thumb to the above algorithm is the smaller the step size (over the increased number of length interval), the more accurate the approximate solution. However, even when extremely small step sizes are used, over a large number of steps the error starts to accumulate and the estimate diverges from the actual functional value which denotes its limitation and the requires more time to compute the approximate solution using the CAS.

IV. RUNGE-KUTTA METHOD OF ORDER 4(RK4)
The Runge-Kutta method are an important family of methods for approximate solutions of ODEs which were developed by mathematics duo C Runge (1856-1927) and M.W Kutta (1867-1944). In this paper, the Runge-Kutta method (RK4) was considered for a initial value problem (IVP), for a system of the first order ODEs, such that with initial condition x(0) = x 0 . To implement the RK4 method of solution of x(t) and y(t) of the IVP over the time interval t ∈ [a, b]. The time interval was subdivided into n equal intervals and we selected the step size such that where h is called the step size. Consider a problem for the implementation of RK4 To obtain the solution through investigation of built-in algorithm of Mathematica ® 12.2.0, it imperative to note that k s obtains the solution for x(t) and l s for solutions to y(t).

V. NUMERICAL EXAMPLES
A. Using Euler's Algorithm Consider the system of ODE The exact solution was obtained as x(t) = 4e −t sin(t), y(t) = 4e −t cos(t). Results generated using the Euler's algorithm were represented on the Table I Let x1 = 4e −t sin(t),y1 = 4e −t cos(t) be the exact solution, and X n ⇒ x = −z + y, Y n ⇒ y = −x − y be the numerical solution. The graphs in Figure 1 represents the exact solution, numerical solution of x(t) & y(t) and the error plots. Although, from a graphical perspective, Figure 1 (c) shows better agreement (between the exact and numerical solutions) as compared to Figure 1 The graph in Figure 2 represents the exact solution and numerical solution of x(t) & y(t)

B. Using Runge-Kutta Algorithm
Now, let us consider the system of ODEs from equations in Equation (11)-Equation (12). The exact solution was obtained as x(t) = 4e −t sin(t), y(t) = 4e −t cos(t). The results generated by the Runge-Kutta algorithm are represented in the Table II. Let L = 4e −t sin(t), P = 4e −t cos(t) be the exact solution, and X n ⇒ x = −x + y, Y n ⇒ y = −x − y be the numerical solution The graph in Figure 3 represents the exact solution and numerical solution of x(t) & y(t) and error plots.
The same phenomenon we observed under the Euler's algorithm section is also observed in this section. Figure 3 (b) and (d) are error plots having order of magnitudes 10 −1 and 10 −1 respectively. However, Figure 3 (a) and (c) shows that the agreement between the exact and numerical solutions for y(t) is better compared to the agreements between the exact and numerical solutions for x(t).
The graph in Figure 4 represents the exact solution and numerical solution of x(t)

C. Example 2
Consider the IVP of the form In what follows, we reduce Equation (13) to a system of firstorder ODEs The exact solution of the system was obtained as The result generated using Euler's algorithm by comparative analysis is presented in Table III. We define the following and www.ijacsa.thesai.org    The graph in Figure 5 represents the exact solution, numerical solution of x(t) & y(t) and error plots.
For exmaple 2, we can observe a stark contrast to example 1. Figure 5 (a) and (c) shows that the agreement between the exact and numerical solutions for x(t) is better compared to the agreements between the exact and numerical solutions for y(t). The error plots for the graph in Figure 5 (b) shows that the approximation is quite good, as the error has order of magnitude 10 −2 . The error plots for the solutions of y(t) is not quite good (it has dimensions of order 10 0 ). The graph in Figure 6 represents the exact solution and numerical solution of x(t)

D. Using Runge-Kutta Algorithm
Consider the system of ODEs in Equation (14)-Equation (15), result generated by RK4. The numerical solution is presented in Equation (20). The comparative Runge-Kutta algorithm is presented in the Table IV   TABLE IV. RESULTS OBTAINED USING RUNGE-KUTTA ALGORITHM    The graph in Figure 7 represents the exact solution, numerical solution of x(t) & y (t), and error plots.
For this section, the error plots for x(t) has order of magnitude 10 −3 (see Figure 7 (b)). The magnitude of order for the error plots of y(t) still maintains the same value as it was in the section for Euler's algorithm (see Figure 7 (d)).
The graph in Figure 8 represents the exact solution and numerical solution of x(t)

E. Example 3
In what follows, we will consider the well-known model of a vibrating particle attached to a stationary wall through a spring mx (t) + dx (t) + kx(t) = 0, where m is the mass of the particle, d is the damping factor, and k is the stiffness constant of the spring. Equation (21) is rewritten as The results generated from the Euler's algorithm is presented in     The graph in Figure 9 represents the exact solution, numerical solution of x(t) & y(t) and error plots. The graph in Figure 11 represents the exact solution and numerical solution of x(t)

F. Using Runge-Kutta Algorithm
Now, let us consider the system of ODEs in Equation (23)-Equation (24), for the results generated using RK4 and its numerical solution. The comparative analysis of RK4 algorithm is presented in the Table VI. The graph in Figure 10 represents the exact solution, numerical solution of x(t) & y(t), and error plots.
The graph in Figure 12 represents the exact solution and numerical solution of x(t)

VII. CONCLUSION
In this paper, we considered an algorithm which was simulated using Mathematica® 12. absolute timing was used to generate a better approximation in the three examples considered in preceding sections. The accuracy of this algorithm depends on the system of ODEs under consideration. The algorithm was built to simulate system of ODEs by converting second-order ODE to a system of first-order IVP, with graphical representation of exact and numerical solution. The error plot with the same step size h = 0.1 was also graphically illustrated. Considering the parametric plot of example 1, the plot shows that the solution presented by the both algorithm, that RK4 algorithm gives a better approximation than the Euler's algorithm. This was also the case in example 2, the Runge-Kutta algorithm showed that the solution obtained from the exact and numerical solution by comparative analysis obtained accurate result and the same goes for example 3. It has been observed by the results shown in Figure 13, through the parametric plot that the algorithm works well. Table VII shows how the Euler's algorithm computes faster but does not give accurate approximations, like that of RK4 which computes slower but gives better approximation.
In future work, we intend to harness the algorithm to solve higher orders of ODE that can be recasted to system of ODEs.
We also aim to adopt the algorithms to solve other differential equations such as Lotka-Volterra model, and stiff equations, and investigate the algorithm time taken in execution, using Mathematica®. By investigation, it is shown in Figure 9 that the parametric plot of the system of ODEs using the (Euler's and RK4) algorithms concludes that RK4 algorithm works well and gives better approximation. In Equation (25) and Equation (26), we observed the appearance of complex conjugate values in the exact solution. This was due to the low value of damping (δ = 0.01) used in the equation of motion, describing the particle attached to a stationary wall. This explanation also covers for the complex values observed in Table V and  Table VI. The numerical values in the aforementioned tables showed close approximations, despite the complex part. The values of the complex parts are also negligible as they are really small (having dimensions of 10 −34 , 10 −19 and 10 −18 ). Comparing Figure 9 (b) and (d) against Figure 10 (b) and (d), we can observe that the Runge-Kutta algorithm also gave a better approximation to the Euler algorithm. The errors for the Runge-Kutta algorithm was in order 5 × 10 −2 while Euler algorithm was in order 1 × 10 −1 .