Comparison of the Selected Methods Used for Solving the Ordinary Differential Equations and Their Systems

: Ordinary differential equations (ODEs), and the systems of such equations, are used for describing many essential physical phenomena. Therefore, the ability to efﬁciently solve such tasks is important and desired. The goal of this paper is to compare three methods devoted to solving ODEs and their systems, with respect to the quality of obtained solutions, as well as the speed and reliability of working. These approaches are the classical and often applied Runge–Kutta method of order 4 (RK4), the method developed on the ground of the Taylor series, the differential transformation method (DTM), and the routine available in the Mathematica software (Mat).


Introduction
Many problems in mathematics, physics or technics are modeled with the use of ordinary differential equations (ODEs) and the systems of such equations. However, these models are sometimes so complex that it is impossible to apply the analytical methods for solving them. In such cases, the only possibility lies in using some approximate methods. In this paper, we intend to compare three selected approaches dedicated for solving the ODEs. The first of them belong to the group of Runge-Kutta methods, very popular because of the ease of implementation and good speed of convergence of the obtained approximate solution. The most often used method, belonging to this family, is the classical Runge-Kutta method of order 4 (RK4), which is applied in the presented elaboration. This method is very popular and described in many sources; therefore, we omit here its detailed description. Instead, we recommend to the interested readers, for example, Ref. [1]. Let us only recall that in order to solve the ODE task with the aid of the RK4 method, we need to present the problem in the form of equation y = f (x, y), x ∈ [a, b], with the initial condition y(a) = y 0 . Function y(x) is the sought function, whereas f is the given function of two variables, defined in some region Ω ⊂ R 2 . In this method, we obtain the discrete approximate solution from the relation y i+1 = y i + h 6 (k 1 + 2k 2 + 2k 3 + k 4 ), k = 0, 1, . . . , n − 1, where y 0 is known from the initial condition, h = b−a n , whereas where x i = a + ih, i = 0, 1, . . . , n.

Differential Transformation Method
Since the differential transformation method (DTM) is strictly connected with the expansion of a function into the Taylor series, we assume that we take into consideration only such functions, which can be expanded into the Maclaurin series (we can discuss the Maclaurin series, because the simple transformation x = x − α reduces the problem of expanding function f into the Taylor series to the expansion of function f into the Maclaurin series for variable x ). Such functions are called "the originals". So, if function f is an original, then the following equality holds which results directly from the theorem about the unique expansion of a function into the Maclaurin power series. This theorem states that any expansion of a function into the power series at point zero is equivalent to the expansion of this function into the Maclaurin power series. Every original f corresponds to function F with nonnegative, integer arguments k, k = 0, 1, 2, . . ., according to the formula Function F is called the T-function of function f , whereas the discussed transformation is called the Taylor transformation.
By having the T-function F, one can find the corresponding original in the form of its expansion into the Maclaurin series, according to the following formula The Taylor transformation, introduced in [5], possesses many useful properties. Some of them result directly from the form of the Maclaurin series of the given originals, while the others must be derived in more complicated way. The discussed properties are of great practical importance and make this tool quite simple to apply. In Section 2.1, we present the properties of DTM, which will be used in the current elaboration. Proofs of these properties can be found, among others, in papers [6,7]. The authors of the current paper have also some achievements in this matter, which will be presented in detail in future papers. The DTM method quickly gained recognition and was used for many different problems, such as for solving ODE systems [8], partial differential equations [9], differentialalgebraic equations [7], selected types of equations, such as the Schrödinger equation [10], Riccati equation [11] or Bratu problem [12], integral equations and integro-differential equations [13,14], fuzzy differential equations [15], differential-difference equations [16], fractional differential equations [17]. The authors successfully used this method also in the calculus of variations, the difference equation, the difference-differential equation and various types of systems of the equations mentioned above (see, for example, the application of a similar method in [18] for the systems of equations).

Properties of DTM
In all the investigated examples, we assume that the considered functions are the originals and that k = 0, 1, 2, . . .

Initial Assumptions and Remarks
We investigate the differential equations of order n (or the systems of such equations) possible to present in the form where y(x) is the sought function of variable x, x ∈ X ⊂ R, F is the function of n + 2 variables defined in some region Ω ⊂ R n+2 . Many methods, including the RK method, require Equation (12) to be presented in the normal form, which means the form of Equation (12) that is possible to solve with respect to variable y (n) : This is quite a strong requirement since many nonlinear differential equations cannot be presented in the normal form; therefore, the applicability of the RK methods is significantly limited.
The next requirement for applicability of the RK methods is the form of conditions ensuring the uniqueness of solutions. The most frequent form is the form of the Cauchy problem in which the conditions are presented as follows where x 0 ∈ X and y 0i ∈ R for i = 0, 1, . . . , n − 1. It happens, however, that we look for the solution of Equation (12) by having the conditions written in the form where x i ∈ X, y i ∈ R, i = 0, 1, . . . , n − 1, or having the combination of conditions (14) and (15). The next problem, which can be met while using the RK methods, concerns the differential equations of the higher orders. In a case when such a kind of equation cannot be presented in the normal form, it is possible sometimes to transform it into the system of n differential equations of the first order (the number of equations in the system is equal to the order of input equation).

Examples
In this section, we present a few selected examples, on the base of which we compare the investigated approaches.

Example 1
The first discussed equation is the Riccati equation of the following form One can easy check that the solution of Equation (16) is given by function y(x) = sinh x. Mathematica software finds the analytical form of the general solution of this equation, but it contains the integral that is impossible to determine, and therefore, such a solution is useless in practice. However, there is no trouble to obtain the numerical solution in the Mathematica software. Similarly easy is determining the approximate solution with the use of the RK4 method. By applying the DTM method for solving Equation (16), we use the Taylor transformation. For this purpose, we apply, among others, Formulas (3), (6)-(9), and we obtain the relation, given below, thanks to which we can determine recursively the values of the successive coefficients Y k for k ≥ 10: From the initial condition, we have Y(0) = 0 and by substituting the successive values k ≥ 0 to the relation (17), we obtain in turn Thanks to this, we can construct the following approximate solution based on eight components or based on ten components as follows In Figure 1, we present the comparison of solutions obtained with the aid of discussed methods together with the absolute errors ∆ of the received approximations. In this figure, the solution obtained from the RK4 method with 21 discretization nodes is denoted by the blue dots and from the RK4 method with 41 discretization nodes denoted by brown dots. Solution (18) is denoted by yellow line and solution (19) by the red line, whereas the solution calculated in the Mathematica software is denoted by the green dashed line. It is worth emphasizing that by using the DTM method, we obtain the exact values of the successive elements of the expansion of function sinh x into the Maclaurin series. Noting this fact, we can easily find the exact solution.

Example 2
In the second example, we discuss the system of nonlinear differential equations which can be presented in the normal form x ∈ [0, π], with conditions One can easily check that the exact solution of this system is represented by functions It is worth noticing that this exact solution cannot be found by using the Mathematica software; however, there is no problem to obtain the approximate solution of the above defined system in this way.
Since the system (20) is presented in the normal form, we can apply the RK4 method, but before we do that, we have to modify this method properly in order to adapt it to the systems of differential equations. In this case (in the case of systems with more equations, the modification is similar), we have (y 0 and z 0 are known from the initial condition (21) and h = b−a n ): where x i = a + ih, i = 0, 1, . . . , n. Thus, we obtain for i = 0, 1, . . . , n − 1.
Similar to the previous example, the greatest attention should be paid to the DTM method. This time, by using, among others, the properties (3), (4), (6)-(10), we transform the system (20) with conditions (21) into the following system for k = 0, 1, 2, . . .: where Substituting the successive values k = 0, 1, 2, . . . in the system (22), we obtain Generating in this way the successive values of T-functions Y i and Z i , i ≤ 13, we can construct the approximate solution based on eleven components .  ). In the plots, the blue dots denote the solutions obtained by using the RK4 method with 51 discretization nodes, brown dots-by using RK4 method with 101 discretization nodes, yellow line-solutions y 15 and z 15 obtained by using the DTM method, red linesolutions y 19 and z 19 received by using the DTM method, and finally the green dashed line represents the solution calculated by applying the Mathematica software. It should be emphasized that by applying the DTM method, we get the exact values of the successive elements of the expansions of functions y and z into the Maclaurin series. This time, however, it is not easy to notice what the analytical forms of these functions are. Therefore, we use the routine FindGeneratingFunction, implemented in the Mathematica software, thanks to which we can obtain the exact solutions (20) and (21) of the investigated problem.

Example 3
In this example, we intend to find the solution of the initial problem described by the system of equations for x ∈ [0, π 2 ], with conditions In order to apply the RK4 method for solving problems (23) and (24), we need to transform the discussed problem into the system of equations of the first order. Having the n-order differential equation of the form y (n) (x) = f x, y(x), y (x), . . . , y (n−1) (x) , for x ∈ [a, b], with conditions y(a) = y 0 , y (a) = y 1 , y (a) = y 2 , . . . , y (n−1) (a) = y n−1 , we can introduce new variables z 1 (x) = y(x), z 2 (x) = y (x), z 3 (x) = y (x), . . . , z n (x) = y (n−1) (x), thanks to which the given equation can be written in the form of the system of the first-order differential equations, as follows . . , z n (x)), with conditions z 1 (a) = y 0 , z 2 (a) = y 1 , z 3 (a) = y 2 , . . . , z n (a) = y (n−1) .
By acting in this way, we transform problems (23) and (24) into the form of the system of four differential equations of the first order, and next, similar to example 4.2, we adapt the formulas describing the RK4 method to the dimension of the task.
Again, similar to the previous examples, Mathematica software is not able to give the solution of the investigated problem in the analytical form, but it can easily give the approximate solution.
Solving this problem with the aid of the DTM method, we assume that the functions Y and Z represent the images of functions y and z. Basing on the properties (4)- (7) and (9), we obtain the following system of equations from the system of Equations (23) and conditions (24), for k = 0, 1, 2, . . .: where Substituting in the system (25) the successive values k = 0, 1, 2, . . ., we obtain , . . .
Generating in this way the successive values of T-functions Y i and Z i , for i ≤ 16, and applying next, once again, the routine FindGeneratingFunction, built in the Mathematica software, we receive the exact solution of the examined problem y 1 = − sin 2 x, y 2 = 1 + sin 2x.
Correctness of this result can be easily verified. In Figure 4, we show the comparison of solutions for functions y (left figure) and z (right figure) obtained by using all the examined approaches, whereas Figure 5 presents the absolute errors ∆ of these approximate solutions for function y (left figure) and z (right figure). In the plots, the blue dots denote the solutions obtained by using the RK4 method with 51 discretization nodes, brown dots-by using the RK4 method with 101 discretization nodes, yellow line-solutions y 15 and z 15 received by using the DTM method, red linesolutions y 19 and z 19 obtained by using the DTM method, and finally the green dashed line represents the approximate solution given by the Mathematica software.

Example 4
As the next example, in this section, we consider the system of differential equations with the boundary conditions, different from the previous example when the system is completed by the initial conditions. So, we have the system of differential equations with boundary conditions of the first kind Because of the form of conditions (27), this problem cannot be solved by means of the RK4 method. It is possible to solve this problem by applying the Mathematica software, no matter the form of conditions (27). In this case, like in the previous examples, the useful analytical solution cannot be found, but there is no problem with obtaining the approximate solution.
In order to use the DTM method, we need to notice that we do not know the values of the first-order derivatives of functions y and z for x = 0. Therefore, we assume temporarily that y (0) = α ∈ R and z (0) = β ∈ R. Hence, and from the half of conditions (27), we obtain If we assume that functions Y k and Z k are the images of functions y(x) and z(x), respectively, then from Equation (27), on the basis of properties (3), (4), (6)-(9), we receive the following relations for k = 0, 1, 2, . . . : Assuming that we determine five successive components from the above relations, that is, by taking 0 ≤ k ≤ 4, we obtain -for k = 0: -for k = 3: -for k = 4: Thus we have the approximate solution depending on the unknown parameters α and β. To determine their values we need to use the other half of the defined boundary conditions and solve the system of equations y 6 (π) = e −π − 1, z 6 (π) = π − 1.
It appears that three solutions of this systems exist: In order to select the best one from among these solutions, we construct three groups of the approximate solutions, y 6,i , z 6,i , and i = 1, 2, 3, and we examine the errors of these solutions, ∆y 6,i = |y 6,i − P (y 6,i )|, ∆z 6,i = |z 6,i − P (z 6,i )|, for i = 1, 2, 3, where P means the right-hand sides of the system (26) of equations.
We can observe that the best solution is the one provided by the third pair of parameters. Figure 6 shows the exact solution and the approximate solution y 6 = y 6,3 , z 6 = y 6,3 received for the selected pair of parameters together with the errors of these approximations. If we increase the number of the determined components and create in this way the approximate solution y 8 , z 8 , then we obtain the result displayed in Figure 7 (in this case, we have five groups of solutions and the best one is for α = −0.0028424 and β = 0.000929434chosen in the same way as described above). Figures 6 and 7 could be made on the basis of the known exact solutions, which are the functions whereas the parameters α and β then take the values α = β = 0.

Example 5
As the last example in this section, let us consider the following differential equation with condition y(0) = 0.
Equation (28) cannot be solved by using the RK method (this equation cannot be presented in the normal form). Mathematica software is also helpless in this case (even the approximate solution cannot be found). The last hope is the DTM method. In order to apply the DTM method, let us assume that function Y k is the image of function y(x). Then, from Equation (28), basing, among others, on properties (3) and (6)-(10), we obtain for k = 0, 1, 2, . . .
-for k = 1: -for k ≥ 2: Solving the equations, described by means of relation (30), for the successive values k ≥ 2, we obtain the sought values Y k for k ≥ 3: On this ground, we can suppose that Hence we are able to predict the solution of problem (28): The above function is in fact the exact solution of task (28), which can be easily verified.

Conclusions
In the paper, we have presented the comparison of three approaches used for solving ordinary differential equations of any order, as well as the systems of such equations. The main goal of this elaboration was to indicate and justify the advantages of the differential transformation method (DTM), which is less known and not often applied. For this purpose, we investigated this method on the background of two other popular methods: the Runge-Kutta method (we decided to use the Runge-Kutta method of order 4 to emphasize better the advantages of the DTM method) and the methods built in the Mathematica software (the discussed examples and the figures illustrating the obtained solutions were realized in version 12.2 of the Mathematica program). The choice of these methods was dictated by their popularity (considering especially the RK4 method) or their simplicity-Mathematica software offers the special routines dedicated for solving the differential equations and their systems, based on the implementation of selected methods and numerical treatment of differential equations. The offered routines are the obvious and default tools for Mathematica users. It is worth mentioning that one does not need to purchase the license for using Mathematica because the tools of this software are available freely on the web page https://www.wolframalpha.com (last access date: 18 December 2021).
The exemplary problems, presented in this paper, were not accidentally chosen. The examples describe the problems that are possible to solve by using all the three discussed approaches (Examples 1-3), impossible to solve by using the RK4 method (Example 4because of the form of conditions; and Example 5-because of the implicit form of the equation) and impossible to solve with the aid of the Mathematica program (Example 5). Despite these problems, all the presented examples could be successfully solved by applying DTM method which shows universality of this method-this method can be used for various forms of initial conditions and for various forms of the solved equation, including the implicit forms. Moreover, the errors of approximate solutions, obtained in the DTM method, are comparable, or very often even lower, than the errors generated by the other methods. What is more and worth emphasizing, the DTM method often leads to the exact solution of the problem, if it exists.
The usefulness of the DTM method, shown in this paper, encourages undertaking a further study on some other applications of this method. At present, the authors are investigating the possibility of applying the discussed method for solving the differentialintegral-algebraic equations and the integral equations with retarded argument. Many interesting and important problems are also described by means of the fractional differential equations (see for example [19,20]), possible to solve by applying the DTM method as well.