# Directly Solving Special Second Order Delay Differential Equations Using Runge-Kutta-Nyström Method

Authors: M. Mechee; F. Ismail; N. Senu; Z. Siri
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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.


---

## Body

## 1. Introduction

A special second order differential equations (ODEs) of the form(1)y′′(t)=f(t,y(t)),y(t0)=α,y′(t0)=β,t≥t0
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 [1–6].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 [7–12]. 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:
(2)y′′(t)=f(t,y(t),y(t-τ1),y(t-τ2),…,y(t-τn)),y(t-τ2),…,y(t-τn)y((tτn)τ2),…,yt>t0
with initial conditions
(3)y(t)=φ(t),y′(t)=φ′(t),t≤t0.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:(4)y′′=f(t,y(t)),y(t0)=α,y′(t0)=β,t>t0
and are usually termed as Runge-Kutta-Nyström formula (RKN) since their introduction in 1925 by Nyström. An s-stage RKN method for the numerical integration of the IVP in (4) is given by
(5)yn+1=yn+hyn′+h2∑i=1s‍biki,yn+1′=yn′+h∑i=1s‍bi′ki,
where
(6)ki=f(tn+cih,yn+cihyn′+h2∑j=1s‍aijkj)
for i=1,2,…,s. The RKN parameters cj,aij,bj, and bj′ are assumed to be real for i,j=1,2,3,…,s and s is the number of stages of the method. Let us introduce the s-dimensional vectors c,b, and b′; moreover, the matrix A is s×s, where c=[c1,c2,c3,…,cs], b=[b1, b2, b3, …, bs],  b′=[b1′, b2′, b3′, …, bs′], and A=[aij], respectively. RKN method can be expressed in Butcher notation using the table of coefficients as follows:
(7)cAbTb′T
RKN methods can be divided into two classes explicit methods when aij=0 for i≤j and implicit methods otherwise.


## 3. Runge-Kutta-Nyström for DDEs

Consider second order DDE(8)y′′=f(y(t-τ3)t,y(t),y(t-τ1),y(t-τ2),y(t-τ3),…,y(t-τm)),t≥t0
with initial conditions
(9)y(t)=φ(t),y′(t)=φ′(t),t≤t0yn is the numerical solution of (8) for all tn>t0, n=1,2,3,….The RKN for second order ODE has been adapted to solve second order DDE and the formula can be written as follows:(10)yn+1=yn+hyn′+h2∑i=1s‍biki(tn+cih,Yi,y(tn+cih-τ1),y(tn+cih-τ2),…,y(tn+cih-τm)),yn+1′=yn′+h∑i=1s‍bi′ki(tn+cih,Yi,y(tn+cih-τ1),y(tn+cih-τ2),…,y(tn+cih-τm)),
where
(11)ki=f(tn+cih,yn+cihyn′+h2∑j=1s‍aijkj,y(xn+cih-τ1),y(tn+cih-τ2),…,y(tn+cih-τm)h2∑j=1s‍aijkj,)
for  i=1,2,…,s. It can be written as follows:
(12)yn+1=yn+hyn′+h2∑i=1s‍bif(tn+cih,Yi,zi),yn+1′=yn′+h∑i=1s‍bi′f(tn+cih,Yi,zi),
where
(13)Yi=yn+cihyi′+h2∑j=1s‍aijkj(tn+cih,Yj,zi),zi=(y(tn+cih-τ1),y(tn+cih-τ2),…,y(tn+cih-τm)).


### 3.1. Time Delay Interpolation

Let the interval of the definition of delay differential equation (8) be I=[a,b] and ti=a+ih for i=0,1,…,n and h=(b-a)/n, where n is the number of points in interval I. The numerical method approximates the solutions at point ti for i=1,2,…,n, to approximate the solution yi+1 at point ti+1 for i=0,1,…,n-1. Hence,
(14)yi+1=yi+1(ti,y(ti-τ1),y(ti-τ2),…,y(ti-τm)).
The method of interpolation has the following cases.Case 1.
If time delayτ1,τ2,…,τm are constants and suppose τi1<τi2<⋯<τim for i1,  i2,  …,  im=1,  2,  …,  m, therefore we can interpolate y(ti-τ1),y(ti-τ2),…,y(ti-τm) as follows:
(15)y(ti-τ1)=fτ1(ti-d,ti-d+1,…,ti,y(ti-d),y(ti-d+1),…,y(ti)),y(ti-τ2)=fτ2(ti-d+1,ti-d+2,…,ti,y(ti-d+1),y(ti-d+2),…,y(ti),y(ti-τ1)),y(ti-τ3)=fτ3(ti-d+2,ti-d+3,…,ti,y(ti-d+2),y(ti-d+3),…,y(ti),y(ti-τ1),y(ti-τ2)),⋮⋮y(ti-τd)=fτd(ti-1,ti,y(ti-1),y(ti),y(ti-τ1),y(ti-τ2),…,y(ti-τd-1)),y(ti-τj)=fτj(ti,y(ti),y(ti-τ1),y(ti-τ2),…,y(ti-τd)),m≥j≥d.The functionsfτj for j=1,2,…,m depend on the interpolation which is used in the numerical method which has degree d.Case 2.
Ifτ1,τ2,…,τm are the variables of time delays, that is, τj=τj(t) for j=1,2,…,m, however, we can consider
(16)y(ti-τj)=fτj(ti-d,ti-d+1,…,ti,y(ti-d),y(ti-d+1),…,y(ti))
for j=1,2,…,m, knowing that yk=y(tk) for k=i-d,i-d+1,…,i. The function fτj depends on the type and degree of the interpolation. In this paper, we consider the cubic interpolation to approximate DDEs with three time delays.


## 3.1. Time Delay Interpolation

Let the interval of the definition of delay differential equation (8) be I=[a,b] and ti=a+ih for i=0,1,…,n and h=(b-a)/n, where n is the number of points in interval I. The numerical method approximates the solutions at point ti for i=1,2,…,n, to approximate the solution yi+1 at point ti+1 for i=0,1,…,n-1. Hence,
(14)yi+1=yi+1(ti,y(ti-τ1),y(ti-τ2),…,y(ti-τm)).
The method of interpolation has the following cases.Case 1.
If time delayτ1,τ2,…,τm are constants and suppose τi1<τi2<⋯<τim for i1,  i2,  …,  im=1,  2,  …,  m, therefore we can interpolate y(ti-τ1),y(ti-τ2),…,y(ti-τm) as follows:
(15)y(ti-τ1)=fτ1(ti-d,ti-d+1,…,ti,y(ti-d),y(ti-d+1),…,y(ti)),y(ti-τ2)=fτ2(ti-d+1,ti-d+2,…,ti,y(ti-d+1),y(ti-d+2),…,y(ti),y(ti-τ1)),y(ti-τ3)=fτ3(ti-d+2,ti-d+3,…,ti,y(ti-d+2),y(ti-d+3),…,y(ti),y(ti-τ1),y(ti-τ2)),⋮⋮y(ti-τd)=fτd(ti-1,ti,y(ti-1),y(ti),y(ti-τ1),y(ti-τ2),…,y(ti-τd-1)),y(ti-τj)=fτj(ti,y(ti),y(ti-τ1),y(ti-τ2),…,y(ti-τd)),m≥j≥d.The functionsfτj for j=1,2,…,m depend on the interpolation which is used in the numerical method which has degree d.Case 2.
Ifτ1,τ2,…,τm are the variables of time delays, that is, τj=τj(t) for j=1,2,…,m, however, we can consider
(16)y(ti-τj)=fτj(ti-d,ti-d+1,…,ti,y(ti-d),y(ti-d+1),…,y(ti))
for j=1,2,…,m, knowing that yk=y(tk) for k=i-d,i-d+1,…,i. The function fτj 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 [13–17]. To study the stability of numerical method (10), consider the linear test equation
(17)y′′=λ2y(t)+μ2y(t-τ).When the method is applied to the linear test equation (17), we have
(18)yn+1=yn+hyn′+h2∑i=1s‍bif(tn+cih,Yi,y(tn+cih-τ))=yn+hyn′+h2∑i=1s‍bi(λ2Yi+μ2y(tn+cih-τ)),yn+1′=yn′+h∑i=1s‍bi′ki(tn+cih,Yi,y(tn+cih-τ))=yn′+h∑i=1s‍bi′ki(λ2Yi+μ2y(tn+cih-τ)),
where
(19)Yi=yn+cihyn′+h2∑j=1s‍aijf(tn+cih,Yj,y(tn+cih-τ))=yn+cihyn′+h2∑j=1s‍aij(λ2Yj+μ2y(tn+cih-τ)),Zn+1=TZn+λ2h2BY+μ2h2BZn(τ)
such that
(20)Zn=(ynhyn′),T=(1101),B=(b1b2⋯bsb1′b2′⋯bs′),Zn(τ)=(y(tn+c1h-τ)y(tn+c2h-τ)⋮y(tn+csh-τ)).
So
(21)(Y1Y2⋮Ys)=(ynyn⋮yn)+(c1c2⋮cs)hyn′+h2(a11a12⋯a1sa21a22⋯a2s⋮⋱⋮as1as2⋯ass)×(λ2Y1+μ2y(tn+c1h-τ)λ2Y2+μ2y(tn+c2h-τ)⋮λ2Ys+μ2y(tn+csh-τ))=(11⋮1)yn+(c1c2⋮cs)hyn′+h2λ2A(Y1Y2⋮Ys)+μ2h2A(μ2y(tn+c1h-τ)μ2y(tn+c2h-τ)⋮μ2y(xn+csh-τ)).
This implies that
(22)(I-HλA)(Y1Y2⋮Ys)=(1c11c2⋮⋮1cs)Zn+HμA(μ2y(tn+c1h-τ)μ2y(tn+c2h-τ)⋮μ2y(xn+csh-τ)),
where Hλ=(λh)2,  Hμ=(μh)2.So(23)(Y1Y2⋮Ys)=(I-HλA)-1×(CZn+HμA(μ2y(tn+c1h-τ)μ2y(tn+c2h-τ)⋮μ2y(tn+csh-τ))),
where
(24)A=(a11a12⋯a1sa21a22⋯a2s⋮⋱⋮as1as2⋯ass),C=(1c11c2⋮⋮1cs).
Hence, the stability polynomial of the method is
(25)Zn+1=TZn+HλB(I-HλA)-1×(CZn+HμAZn(τ))+HμBZn(τ),Zn+1=T1Zn+T2Zn(τ),Zn(τ)=(y(tn+c1h-τ)y(tn+c2h-τ)⋮y(xn+csh-τ)),T1=T+HλB(I-HλA)-1C,T2=HμB[Hλ(I-HλA)-1A+I].


## 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 Figures1, 2, 3, 4, 5, and 6:(i)
h: 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: Maxn|y(xn)-yn| Absolute value of the true solution minus the computed solution.Figure 1
The efficiency curves for all methods for Problem1 with tend=1 and h=1/4i,i=1,…,5.Figure 2
The efficiency curves for all methods for Problem2 with tend=1 and h=1/4i,i=1,…,5.Figure 3
The efficiency curves for all methods for Problem3 with tend=1 and h=1/4i,i=1,…,5.Figure 4
The efficiency curves for all methods for Problem4 with tend=1 and h=1/4i,i=1,…,5.Figure 5
The numerical solution for all methods for Linear Problem withtend=10 and h=1/10.Figure 6
The numerical solution for all methods for Nonlinear Problem withtend=10 and h=1/10.Problem 1 (nonlinear).
Consider the following:(26)y′′(t)=e-2τy2(t-τ)y(t),t≥t0,y(t)=e-t,y′(t)=-e-t,t≤t0,τ=h10.
Exact solution:y(t)=e-t.Problem 2 (nonlinear).
Consider the following:(27)y′′(t)=y2(t-τ)-14(1+t)3-(1+t)+τ,t≥t0,y(t)=1+t,y′(t)=121+t,t≤t0,τ=h10.
Exact solution:y(t)=1+t.Problem 3 (linear example).
Consider the following:(28)y′′(t)=y(t-τ)-1(1+t)2-ln(1+t-τ),t≥t0,y(t)=ln(1+t),y′(t)=11+t,t≤t0,τ=h10.
Exact solution:y(t)=ln(1+t).Problem 4 (nonlinear).
Consider the following:(29)y′′(t)=13(e-τ1y(t-τ1)+e-τ2y(t-τ2)+e-τ3y(t-τ3)),t≥t0,y(t)=e-t,y′(t)=-e-t,t≤t0,τ1=h10,τ2=h20,τ3=h30.
Exact solution:y(t)=e-t.


### 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:(30)y′′(t)+(δ+acost)y(t)+cy3(t)=by(t-T),
which is a nonlinear delay differential equation.Whereδ,a,b,c, and T are parameters.δ is the frequency squared of the simple harmonic oscillator, and a is the amplitude of the parametric resonance, and b is the amplitude of delay while c is the amplitude of the cubic nonlinearity and T 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δ=a=b=1 and c=0 we obtained the following linear Mathieu equation:
(31)y′′(t)=(1+cost)y(t)=y(t-T),t∈[0,10],y(t)=sin(t),y′(t)=cos(t),t≤0,
where T=τ=h/10 is the delay term, the exact solution does not exist.Whenδ=a=b=c=1, we obtained the following nonlinear Mathieu equation:
(32)y′′(t)=(1+cost)y(t)+y3(t)=y(t-T),t∈[0,10],y(t)=sin(t),y′(t)=cos(t),t≤0,
where T=τ=h/10 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 Figures5 and 6.


## 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:(30)y′′(t)+(δ+acost)y(t)+cy3(t)=by(t-T),
which is a nonlinear delay differential equation.Whereδ,a,b,c, and T are parameters.δ is the frequency squared of the simple harmonic oscillator, and a is the amplitude of the parametric resonance, and b is the amplitude of delay while c is the amplitude of the cubic nonlinearity and T 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δ=a=b=1 and c=0 we obtained the following linear Mathieu equation:
(31)y′′(t)=(1+cost)y(t)=y(t-T),t∈[0,10],y(t)=sin(t),y′(t)=cos(t),t≤0,
where T=τ=h/10 is the delay term, the exact solution does not exist.Whenδ=a=b=c=1, we obtained the following nonlinear Mathieu equation:
(32)y′′(t)=(1+cost)y(t)+y3(t)=y(t-T),t∈[0,10],y(t)=sin(t),y′(t)=cos(t),t≤0,
where T=τ=h/10 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 Figures5 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.



---

*Source: 830317-2013-05-27.xml*