Novel approach for solving higher-order differential equations with applications to the Van der Pol and Van der Pol–Duffing equations

Background Numerical methods are used to solve differential equations, but few are effective for nonlinear ordinary differential equations (ODEs) of order higher than one. This paper proposes a new method for such ODEs, based on Taylor series expansion. The new method is a second-order method for second-order ODEs, and it is equivalent to the central difference method, a well-known method for solving differential equations. The new method is also simple to implement for higher-order differential equations. The proposed technique was applied to solve the Van der Pol and Van der Pol–Duffing equations. It is stable over a wide range of nonlinearity and produces accurate and reliable results. For the self-excitation Van der Pol equation, the proposed technique was applied with different values of nonlinear damping. Results The results were compared with those obtained using the ODE15s solver in MATLAB. The two sets of results showed excellent agreement. For the forced Van der Pol–Duffing equation, the proposed technique was applied with different values of exciting force amplitude and frequency. It was found that for certain conditions, the solution obtained using the proposed technique differed from that obtained using ODE15s. Conclusions The solution obtained using the proposed technique showed good agreement with the solutions obtained using ODE45 and Runge–Kutta fourth order. The results show that the proposed approach is very simple to apply and produces acceptable error. It is a powerful and versatile tool for solving of high-order nonlinear differential equations accurately.


Background
Ordinary differential equations (ODEs) are used to model a broad range of problems, such as the motion of objects, the spread of diseases, and the behavior of financial markets.Many ODEs are nonlinear.Nonlinear ODEs can be difficult to solve analytically, so numerical methods are often used to approximate their solutions.However, there are few effective numerical methods for solving nonlinear ODEs, especially those of higher order.The greatest popular methods are Runge-Kutta [1, 2] and predictor-corrector methods based on Adams-Bashforth and Adams-Moulton methods [3,4], or their variations [5][6][7][8].
Runge-Kutta methods are a family of numerical methods that are commonly utilized to solve ODEs.They are known for their accuracy and efficiency, and they can be used to solve both linear and nonlinear ODEs.Predictor-corrector methods are another family of numerical methods that can be used to solve ODEs.They are based on using a predictor method to estimate the solution at the next time step, and then using a corrector method to improve the estimate.Predictor-corrector methods are often more accurate than Runge-Kutta methods, but they can also be more computationally expensive.
Adams-Bashforth and Adams-Moulton methods are two specific types of predictor-corrector methods.They are often used together, with Adams-Bashforth methods used for prediction and Adams-Moulton methods used for correction.
"Adams-Bashforth and Adams-Moulton" methods can be used to solve both linear and nonlinear ODEs.However, they are particularly well suited for solving nonlinear ODEs of order higher than one.
This paper proposes a new approach to solve higherorder differential equations emergent from initial value problems (IVPs).The new method is efficient and accurate and can be used to solve any nth-order nonlinear ordinary differential equation emergent from IVPs.This is a momentous advance in the study of numerical methods and could be used in a variety of applications, such as physics, engineering, and finance due to its ease of application.
The new approach has been developed to solve the differential equations of the self-excited and nonlinear oscillators, which are two types of nonlinear differential equations that arise in a variety of physical and engineering applications.These equations are known to be difficult to solve, especially for higher orders.The new approach is able to solve these equations with greater accuracy and efficiency in a simple manner than previous methods.
The solution of the free response of self-excited equation considering different damping effects and with different initial conditions as well as the forced response with different values of exciting force amplitude and frequency is obtained.
Results are compared with that obtained by the ODE15s algorithm for self-excitation Van der Pol oscillator and they are in excellent agreement.The absolute error is of order 10 −3 or less.Results for forced nonlinear oscillator are compared with ODE15s, ODE45 in MATLAB, and Runge-Kutta fourth order.The comparison shows that the present technique and the Runge-Kutta fourth-order (RK4) method produce very similar results with minimum relative error compared to the other techniques.

Methodology of the proposed technique
The initial value problem can be written in a general form as: Subject to the initial conditions In Eq. (1), derivatives of the function f can be denoted by either superscripts or dots.The function f depends on the time variable t, the unknown func- tion y, and its derivatives up to order n.The values of the first n − 1 derivatives of the function y at the initial value t 0 (which is assumed to be 0 unless stated other- wise) are given by the constants a i (where i = 1, 2, …, n − 1).
The new method for solving ODEs requires that the ODE must be rewritten so that the highest-order derivative term is alone on one side of the equation, and the lower-order derivative terms and all other terms are on the other side of the equation.
Since the n th derivative at the initial condition a n is unknown, it can be found by substituting the initial conditions (2) into Eq.(3).Therefore, we will now describe a solution procedure to solve equations of the form (1), subject to initial conditions of the form (2).
The given text is a description of the new approach for solving ordinary differential equations (ODEs) using Taylor series expansion.The method begins by discretizing the domain of the variable t , i.e., divid- ing it into a number of subintervals of the length of �t = (T − t 0 )/N .Then, at each grid point t j , the func- tion y and its derivatives up to order n − 1 are approxi- mated using Taylor series expansion.Finally, the nth derivative is obtained using Eq.(3).
This process is repeated for j = 1 to N , to obtain the approximate values of y and its derivatives at all grid points.The procedure of the numerical solution of ODEs using Taylor series expansion is as follows: 1. Find the nth derivative at the initial condition using Eq.(3). 2. Discretize the domain of the independent variable t. 3.At each grid point t j , approximate y and its deriva- tives up to order n − 1 using Taylor series expansion using Eq.(4).
(1) f y (n) (t), y (n−1) (t), . . ., ẏ(t), y(t), t = 0 (2) 4. Obtain the nth derivative using Eq.(3). 5. Repeat steps 3 and 4 for j = 1 to N. This method is a powerful tool for solving ODEs numerically, but it is essential to remember that the accuracy of the solution is based on the number of subintervals N and the order of the Taylor series expansion n.
The new method expresses the derivatives of a function in a simplified manner using the Taylor expansion directly.This means that instead of using complex mathematical formulas, the derivatives can be calculated using a few simple steps.
For example, when solving a second-order differential equation, the first derivative of the function can be calculated using the central difference method.The central difference method can be used to approximate the first derivative of a function at a point, which can be useful for solving second-order differential equations.This method is based on the Taylor expansion of the function at two points, one slightly ahead of the point of interest and one slightly behind.
In the new method, the Taylor expansion for the first derivative is applied directly in simplified form.This means that the derivatives can be calculated more quickly and easily, and with less chance of error.
The Taylor expansion of the function y at time t + ∆t is: The first derivative of y at time t + ∆t using the central difference method is: Substituting the Taylor expansion of y(t + 2∆t) into the central difference equation, we get: (4) Simplifying, we get: This is the expression for the first derivative of y at time t + ∆t using the new method.The two forms are equivalent, but the second form is simpler and easier to understand.It is also easier to implement in computer code.
The simplified form of the central difference method is particularly useful for students and early-career professionals who are new to numerical methods.

The problem of Van der Pol oscillator
Van der Pol was a Dutch engineer who studied how things move.In 1927, he came up with a new equation to describe how electrical circuits oscillate.This equation is now used to study many different kinds of systems, including human hearts and the Earth's climate.The Van der Pol oscillator is a special type of system that can oscillate on its own.This means that it does not need any external input to keep moving.Van der Pol oscillators are found in many different natural systems, and they are also used in many different technologies.It is a nonlinear differential equation, which means that there is no easy solution, but it can be solved numerically.The selfexcited oscillator equation is a resourceful model that can be utilized to represent a broad range of oscillatory phenomena in many different fields [9][10][11][12].
The self-excited oscillator equation generally takes the form Equation ( 8) is a nonlinear inhomogeneous ordinary differential equation (ODE) with order 2 that describes the dynamics of the self-excited oscillator.It requires two initial conditions to solve.The parameter µ represents the strength of the oscillator damping.The right-hand side of Eq. ( 8) stands for the forced manner of the oscillator and is given as f(t) = Fcos(ωt), where F is the amplitude of the exciting force and ω is the frequency of the (7) oscillating force.If f(t) = 0, then Eq. ( 8) reduces to a nonlinear homogeneous second-order independent ODE, which is Eq. ( 9).
The unforced self-excited oscillator equation is a nonlinear differential equation that is given by Eq. ( 9).Truncated expansions can be utilized to find approximate solutions to this equation, although there are no known exact solutions [13].
Since then many people have tried to solve Eq. (9) using both analytical and numerical methods to see if it has a limit cycle [16].Equation ( 9) has a well-developed theory of existence, uniqueness, and stability, as extensively discussed in [17].Many authors have found the solution to Eq. ( 9), including in [18] using the collocation method and in [19] using MATLAB Ode15s and Ode45 built-in functions, and in [20] using the damped Fourier series method.Also, in [11] using (9) ÿ − µ 1 − y 2 ẏ + y = 0 Subjected to y(t 0 ) = y 0 and ẏ(t 0 ) = ẏ0 predictor-corrector Adam-Bashforth-Moulton method, in [12] using restarted A domain decomposition method, and in [21] utilizing the segmenting recursion method.
The literature explores the Van der Pol oscillator through diverse methods, ranging from classical techniques like harmonic balance and Krylov-Bogoliubov-Mitropolsky to modern approaches like neural networks.Several studies tackle specific aspects: references [22] and [23] focus on different stiffness conditions, [24] introduces a new variant method, [25] utilizes a two-point block method, [26] leverages the KBMM method, [27] employs hybrid functions, and [28,29] applies a meshfree neural network algorithm to a more complex variant of the oscillator.This showcases the breadth and ongoing development of techniques used to analyze and understand this important nonlinear system.
Advances in computational mathematics have made it possible to test the accuracy and performance of new or modified numerical methods for solving nonlinear differential equations by applying them to well-developed physical models, such as the unforced self-excited equation.This model is a classic test problem that is utilized to assess the efficiency and reliability of new methods [11,20].The system is characterized by its linear and nonlinear damping coefficients (ε and μ), its linear and nonlinear stiffness coefficients (α and β), and the excitation amplitude (F) and frequency (ω).The initial displacement and velocity of the system are y 0 and ẏ0 , respectively.

Application of the proposed technique
In this technique, the differential Eq. ( 10) is rearranged as follows: By direct substitution of the initial conditions given in (10), the acceleration at starting point can be written as: The approximate displacement function at a time t + t is obtained using the first of (4): The approximate velocity function at a time t + t is obtained using (4): The approximate acceleration function at a time t + t is obtained using Eq.(11) as: So, the first iteration is obtained from Eqs. ( 13) through (15) as: (10) ÿ + ε + µy 2 ẏ + αy + βy 3 = F cos (ωt) Subjected to y(t 0 ) = y 0 and ẏ(t 0 ) = ẏ0 The recurrence formula can be written as: The recurrence formula is simple, direct, and straightforward.

Results
Presented below are the results of this research.

Free vibration
The recursive relations in Sect. 5 are developed to solve Van der Pol's oscillator Eq. ( 9) with the non-exciting force for different values of the nonlinearity parameter µ .The initial conditions are y(0) = 2 and ẏ(0) = 0.
The solution obtained using the present technique is compared with the solutions obtained using Euler's method and the ODE15s solver in MATLAB.Figure 1 shows the comparison for µ = 20.
Tables 1, 2, 3, 4, 5 and 6 show the absolute error in the present solution compared to the ODE15s solution at a different time to t = 50 s for µ = 0.1, 1, 10, 20, 50, and 100, respectively.The error is of order 10 −3 or less.
This shows that the present technique is a powerful and easy way to solve "Van der Pol's" oscillator problem.
Example 2 is with the following data ε = −µ = 0.2, F = 0.53 and ω = 1 with the same initial conditions as example 1.The results are presented in Fig. 9a.It was found that the results have a big difference from that obtained using ODE15s in the period from t = 30 to t = 35.This prompted me to compare the results of the present solution with the results of other solutions.The results of the present solution were compared with the results of ODE15s, RK4, and ODE45.It has been observed that all techniques are identical in the period from t = 0 to t = 8, as shown in Fig. 9b.The results in the periods from t = 10: 25, t = 25:35, and t = 35: 50 are shown in Fig. 9c-e, respectively.The big difference is in the period t = 25:35.Table 7 presents a comparison of Van der Pol-Duffing oscillator response obtained by different techniques and their relative errors compared to RK4.It is clear that the results of the present technique are in good agreement with RK4 over the whole range.

Conclusions
This paper presents a new technique to solve ordinary differential equations (ODEs) which is based on Taylor expansion.The technique is developed to solve nonlinear Van der Pol and Van der Pol-Duffing oscillators with damping effects under different initial conditions.The results are compared to three other well-known ODE solvers: ODE15s, ODE45, and fourthorder Runge-Kutta (RK4).The comparison shows that the new technique is just as accurate as the other three solvers, but it is simpler to understand and implement in computer code.The new technique is proved to be equivalent to the central difference method, but the authors argue that their simplified form of the central difference method is easier to use, especially for      Van der Pol Oscilator = 100

Table 3
Comparison of the Van der Pol oscillator response for the present solution and the ode15s solution with μ = 10

Table 4
Comparison of the Van der Pol oscillator response for the present solution and the ode15s solution with μ = 20

Table 2
Comparison of the Van der Pol oscillator response for the present solution and the ode15s solution with μ = 1

Table 1
Comparison of the Van der Pol oscillator response for the present solution and the ode15s solution with μ = 0.1

Table 5
Comparison of the Van der Pol oscillator response for the present solution and the ode15s solution with μ = 50

Table 6
Comparison of the Van der Pol oscillator response for the present solution and the ode15s solution with μ = 100