Abstract

Runge-Kutta-Nyström (RKN) method is adapted for solving the special second order delay differential equations (DDEs). The stability polynomial is obtained when this method is used for solving linear second order delay differential equation. A standard set of test problems is solved using the method together with a cubic interpolation for evaluating the delay terms. The same set of problems is reduced to a system of first order delay differential equations and then solved using the existing Runge-Kutta (RK) method. Numerical results show that the RKN method is more efficient in terms of accuracy and computational time when compared to RK method. The methods are applied to a well-known problem involving delay differential equations, that is, the Mathieu problem. The numerical comparison shows that both methods are in a good agreement.

1. Introduction

A special second order differential equations (ODEs) of the form which is not explicitly dependent on the first derivative of the solution are frequently found in many physical problems such as electromagnetic waves, thin film flow, and gravity driven flow. Most researchers, scientists, and engineers used to solve (1) by converting the second order differential equations to a system of first order equations twice the dimension. However there are also studies on numerical methods which directly solve (1) using one-step methods or multistep methods. Such work can be seen in [16].

Most of the methods for solving special second order ODEs can be adapted for solving special second order delay differential equations (DDEs). In recent years there has been a growing interest in numerical solutions of DDEs. This is due to the appearance of such equations in various areas such as neural network theory, epidemiology, and time lag control processes. DDEs also provide us with realistic model of many phenomena arising in real world problems. For example, DDEs can be used in modelling of population dynamics and spread of infectious diseases and two body problems of electrodynamics [712]. Most of the work concerning DDEs in the literature involved first order delay differential equations. Hence, in this research, we are going to focus on numerical methods for solving special second order DDEs. Special second order delay differential equations with multiple delays can be written in the following form: with initial conditions

Runge-Kutta-Nyström method of order four will be adapted for solving second order DDEs (2). Stability polynomial of the method is also presented when applied to linear second order DDE. Numerical results on a set of test problems are given and compared with the numerical results when the problems are reduced to a system of first order DDEs and solved using Runge-Kutta methods. We are also going to solve a well-known problem in engineering which involves second order delay differential equations, that is, the Mathieu problem.

2. Numerical Methods for Second Order ODEs

Runge-Kutta methods are designed for special second order differential equations: and are usually termed as Runge-Kutta-Nyström formula (RKN) since their introduction in 1925 by Nyström. An -stage RKN method for the numerical integration of the IVP in (4) is given by where for . The RKN parameters , and are assumed to be real for and is the number of stages of the method. Let us introduce the -dimensional vectors , and ; moreover, the matrix is , where , ,?,?,?,?,? ,?,?,?,?, and , respectively. RKN method can be expressed in Butcher notation using the table of coefficients as follows: RKN methods can be divided into two classes explicit methods when for and implicit methods otherwise.

3. Runge-Kutta-Nyström for DDEs

Consider second order DDE with initial conditions is the numerical solution of (8) for all , .

The RKN for second order ODE has been adapted to solve second order DDE and the formula can be written as follows: where for . It can be written as follows: where

3.1. Time Delay Interpolation

Let the interval of the definition of delay differential equation (8) be and for and , where is the number of points in interval . The numerical method approximates the solutions at point for , to approximate the solution at point for . Hence, The method of interpolation has the following cases.

Case 1. If time delay are constants and suppose for ,??,??,??,??,??,??, therefore we can interpolate as follows:

The functions for depend on the interpolation which is used in the numerical method which has degree .

Case 2. If are the variables of time delays, that is, for , however, we can consider for , knowing that for . The function depends on the type and degree of the interpolation. In this paper, we consider the cubic interpolation to approximate DDEs with three time delays.

4. Stability of the Method

Stability aspect of numerical methods for delay differential equations has been introduced in [1317]. To study the stability of numerical method (10), consider the linear test equation

When the method is applied to the linear test equation (17), we have where such that So This implies that where ,? .

So where Hence, the stability polynomial of the method is

5. Numerical Results

In this section, some of the problems involving second order DDEs are solved using RKN methods. Then the same set of problems is reduced to a first order DDEs system and solved using RK methods of the same order. Then numerical results are given.

The following notations are used in Figures 1, 2, 3, 4, 5, and 6:(i): stepsize used. (ii)RKN4-N: the existing Rung-Kutta-Nyrström method of order four derived by Senu et al. [18].(iii)RKN4-D: the existing Rung-Kutta-Nyrström method of order four as in [19].(iv)RK4: existing Runge-Kutta method order four as given in [20].(v)DOPRI: existing Runge-Kutta method order fifth as in Dormand [19].(vi)Total time: the total time in second to solve the problems. (vii)MAX ERROR: Absolute value of the true solution minus the computed solution.

Problem 1 (nonlinear). Consider the following:
Exact solution: .

Problem 2 (nonlinear). Consider the following:
Exact solution: .

Problem 3 (linear example). Consider the following:
Exact solution: .

Problem 4 (nonlinear). Consider the following:
Exact solution: .

5.1. An Application to Mathieu Equation

In this section we will apply the RKN and RK methods to solve a well-known equation in engineering, the Matheiu's equation, which is defined as follows: which is a nonlinear delay differential equation.

Where , and are parameters.

is the frequency squared of the simple harmonic oscillator, and is the amplitude of the parametric resonance, and is the amplitude of delay while is the amplitude of the cubic nonlinearity and is the time delay.

Equation (30) is a model for high speed milling, a kind of parametrically interrupted cutting as opposed to the self-interrupted cutting arising in an unstable turning process. More information on the problem can be found in [21]. Various special cases of (30) have been studied, depending on which parameters are zero.

When and we obtained the following linear Mathieu equation: where is the delay term, the exact solution does not exist.

When , we obtained the following nonlinear Mathieu equation: where is the delay term, the exact solution does not exist.

Both the linear and nonlinear Mathieu equations are solved using RKN and RK methods and the results are plotted in Figures 5 and 6.

6. Discussion and Conclusion

In this paper we adapted the Runge-Kutta-Nyström method for solving special second order delay differential equations in which cubic interpolation is used to evaluate the delay term. We also presented the stability of the method when applied to linear second order DDEs. We solved a set of DDEs using RKN4-N by Senu et al. [6] and RKN4-D by Dormand [19]. For comparison purposes the same set of problems a reduced to a system of first order DDEs and solved using classical fourth order Runge-Kutta method and fifth order Runge-Kutta method by Dormand [19]. The log of maximum error is plotted against time. From the numerical results, we observed that both RKN methods are just as efficient. However they are more efficient compared to the fourth order classical Runge-Kutta method followed by the fifth order Runge-Kutta method by Dormand [19]. This is because RKN method directly solved the equations whereas in RK method the equations are reduced to a system of first order DDEs. The fifth order RK method has more stages compared to the fourth order RK method which required more function evaluations at each step. Both RKN and RK methods are also used to solve the linear and nonlinear Mathieu equations. Since they do not have the exact solution, we just plot the numerical results in which both methods show a good agreement.

Acknowledgment

The authors would like to thank the reviewers for their constructive comments and their careful reading of the paper.