# Integration Processes of Delay Differential Equation Based on Modified Laguerre Functions

Authors: Yeguo Sun
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

We propose long-time convergent numerical integration processes for delay differential equations. We first construct an integration process based on modified Laguerre functions. Then we establish its global convergence in certain weighted Sobolev space. The proposed numerical integration processes can also be used for systems of delay differential equations. We also developed a technique for refinement of modified Laguerre-Radau interpolations. Lastly, numerical results demonstrate the spectral accuracy of the proposed method and coincide well with analysis.


---

## Body

## 1. Introduction

Time delay systems which are described by delay differential equations (DDEs) or more generally functional differential equations (FDEs) have been studied rather extensively in the past thirty years since time delays are often the sources of instability and encountered in various engineering systems such as chemical processes, economic markets, chemical reactions, and population dynamics [1, 2]. Nakagiri [1] studied the structural properties of linear autonomous functional differential equations in Banach spaces within the framework of linear operator theory. The system stability and the compactness of the operators describing the solution trajectories are investigated in [2]. Depending on whether the existence of time delays or not, stability criteria for time delay systems can be divided into two types: delay-dependent ones and delay-independent ones. De la Sen and Luo [3] deal with the global uniform exponential stability independent of delay of time-delay linear and time-invariant systems. By exploiting appropriate Lyapunov functional candidate, new delay-dependent robust stability criteria of uncertain time-delay systems are proposed in [4].However, most DDEs do not have analytic solution, so it is generally necessary to resort to numerical methods [5, 6]. It is well known that a numerical method which is convergent in a finite interval does not necessarily yield the same asymptotic behavior as the underlying differential equation. If the numerical solution defines a dynamical system, then we would study whether this dynamical system inherits the dynamics of the underlying system. Hence it is crucial to understand the behavior of numerical solution in order that we may interpret the data and facilitate the design of algorithms which provide correct qualitative information without being unduly expensive. The classical analysis of linear multistep methods [7–9] and Runge-Kutta methods [10–12] for delay differential equations involves assessment of stability and accuracy but has also been supplemented by considerable practical experience and experimentation.As a basic tool, the Runge-Kutta method plays an important role in numerical integrations of delay differential equations. We usually designed this kind of numerical schemes in two ways. The first way is based on Taylor's expansion coupled with other techniques. The next is to construct numerical schemes by using collocation approximation. For instance, Butcher [13, 14] provided some implicit Runge-Kutta processes based on the Radau Quadrature formulas; also see [15, 16] and the references therein.Recently, Guo et al. [17] proposed an integration process for ordinary differential equations based on modified Laguerre functions, which are very efficient for long-time numerical simulations of dynamical systems. But so far, to our knowledge, there is no work concerning the applications of Laguerre approximation to integration process for delay differential equations.In this paper, we construct a new integration processes for delay differential equations based on modified Laguerre functions and establish its global convergence in certain weighted Sobolev space. Numerical results demonstrate the spectral accuracy of the proposed method and coincide well with analysis.


## 2. Modified Laguerre-Radau Function for Delay Differential Equations

Letωβ(t)=e-βt, β>0, and define the weighted space Lωβ2(0,∞) as usual, with the following inner product and norm:
(2.1)(u,v)ωβ=∫0∞‍u(t)v(t)ωβ(t)dt,∥v∥ωβ=(v,v)ωβ1/2.
The modified Laguerre polynomial of degree l is defined by (cf. [18])
(2.2)ℒl(β)(t)=1l!eβtdldtl(tle-βt),l≥0.
For example,
(2.3)ℒ0(β)(t)=1,ℒ1(β)(t)=1-βt,ℒ2(β)(t)=-2βt+β22t2.
The modified Laguerre polynomials satisfy the recurrence relation
(2.4)ddtℒl(β)(t)=ddtℒl-1(β)(t)-βℒl-1(β)(t),l≥1.
The set of Laguerre polynomials is a complete Lωβ2(0,∞)-orthogonal system, namely,
(2.5)(ℒl(β),ℒm(β))ωβ=1βδł,m,
where δl,m is the Kronecker function. Thus, for any v∈Lωβ2(0,∞),
(2.6)v(t)=∑l=0∞‍v^lℒl(β)(t),v^l=β(v,ℒl(β))ωβ.
Let N be any positive integer and 𝒫N(0,∞) the set of all algebraic polynomials of degree at most N. We denote by tβ,jN the modified Laguerre-Radau interpolation points. Indeed, tβ,0N=0, and tβ,jN(1≤j≤N) are the distinct zeros of (d/dt)ℒN+1(β)(t). By using (2.1) and the formula (2.12) of [19], the corresponding Christoffel numbers are as follows:
(2.7)ωβ,0N=1β(N+1),ωβ,jN=1β(N+1)ℒN(β)(tβ,jN)ℒN+1(β)(tβ,jN),1≤j≤N.
For any ρ∈𝒫2N(0,∞),
(2.8)∑j=0N‍ρ(tβ,jN)ωβ,jN=∫0∞‍ρ(t)ωβ(t)dt.
Next, we define the following discrete inner product and norm:
(2.9)(u,v)ωβ,N=∑j=0N‍u(tβ,jN)v(tβ,jN)ωβ,jN,∥v∥ωβ,N=(v,v)ωβ,N1/2.
For any ϕ,ψ∈𝒫N(0,∞),
(2.10)(ϕ,ψ)ωβ=(ϕ,ψ)ωβ,N,∥ϕ∥ωβ=∥ϕ∥ωβ,N.
For all v∈Lωβ2(0,∞), the modified Laguerre-Radau interpolant Iβ,Nv∈𝒫N(0,∞) is determined by
(2.11)Iβ,Nv(tβ,jN)=v(tβ,jN),0≤j≤N.
By (2.10), for any ϕ∈𝒫N(0,∞),
(2.12)(Iβ,Nv,ϕ)ωβ=(Iβ,Nv,ϕ)ωβ,N=(v,ϕ)ωβ,N.
The interpolant Iβ,Nv can be expanded as
(2.13)Iβ,Nv(t)=∑l=0N‍v~β,lNℒl(β)(t).
By virtue of (2.5) and (2.10),
(2.14)v~β,lN=β(Iβ,Nv,ℒl(β))ωβ=β(v,ℒl(β))ωβ,N.
Define the modified Laguerre functions ℒ~l(β)(t)=e-(1/2)βtℒl(β)(t) as the base functions. According to (2.4), the functions ℒ~l(β)(t) satisfy the recurrence relation
(2.15)ddtℒ~l(β)(t)=ddtℒ~l-1(β)(t)-12βℒ~l(β)(t)-12βℒ~l-1(β)(t),l≥1.
Denote by (u,v) and ∥v∥ the inner product and the norm of the space L2(0,∞), respectively. The set of ℒ~l(β)(t) is a complete L2(0,∞)-orthogonal system, that is,
(2.16)〈ℒ~l(β),ℒ~m(β)〉=1βδl,m.
We now introduce the new Laguerre-Radau interpolation. Set
(2.17)QN(0,∞)=span{ℒ~0(β),ℒ~1(β),…,ℒ~N(β)}.
Let tβ,jN and ωβ,jN be the same as in (2.7), and take the nodes and weights of the new Laguerre-Radau interpolation as
(2.18)t~β,jN=tβ,jN,ω~β,jN=1ℒ~N(β)(tβ,jN)ℒ~N+1(β)(tβ,jN)=eβtβ,jNωβ,jN.
The discrete inner product and norm can be defined similarly as
(2.19)〈u,v〉β,N=∑j=0N‍u(tβ,jN)v(tβ,jN)ω~β,jN,∥v∥β,N=〈v,v〉β,N1/2.
For any ϕ1, ϕ2∈QN(0,∞), we have ϕ1=e-(1/2)βtψ1,  ϕ2=e-(1/2)βtψ2, and ψ1, ψ2∈𝒫N(0,∞). Thus by (2.10),
(2.20)〈ϕ1,ϕ2〉β,N=〈ψ1,ψ2〉ωβ,N=〈ψ1,ψ2〉ωβ=〈ϕ1,ϕ2〉.
The new Laguerre-Radau interpolant I~β,Nv∈QN(0,∞) is determined by
(2.21)I~β,Nv(tβ,jN)=v(tβ,jN),0≤j≤N.
Due to the equality (2.20), for any ϕ∈QN(0,∞),
(2.22)〈I~β,Nv,ϕ〉=〈I~β,Nv,ϕ〉β,N=〈v,ϕ〉β,N.
Let
(2.23)I~β,Nv(t)=∑l=0N‍v~β,lNℒ~l(β)(t).
Then, with the aid of (2.16) and (2.22), we derive that
(2.24)v~β,lN=β〈I~β,Nv,ℒ~l(β)〉=β〈I~β,Nv,ℒ~l(β)〉β,N=β〈v,ℒ~l(β)〉β,N.
There is a close relation between Iβ,N and I~β,N. From the previous two equalities, it follows that
(2.25)e(1/2)βtI~β,Nv(t)=∑l=0N‍v~β,lNℒl(β)(t)=β∑l=0N‍(v,ℒ~l(β))β,Nℒl(β)(t)=β∑l=0N‍(e(1/2)βtv,ℒl(β))ωβ,Nℒl(β)(t).
This with (2.14) implies
(2.26)I~β,Nv(t)=e-(1/2)βtIβ,N(e(1/2)βtv(t)).
Consider the following delay differential equation:
(2.27)ddtW(t)=f(W(t),W(t-τ)),t>0,W(t)=Φ(t),-τ≤t≤0.
For any fixed positive integer N, we define
(2.28)U(t)=W(tτtβ,NN),
and denote
(2.29)τ^=τtβ,NN.
Then system of (2.27) can be transformed into
(2.30)ddtU(t)=τ^f(U(t),U(t-tβ,NN)),t>0,U(t)=Φ(tτ^),-tβ,NN≤t≤0.
We suppose that U(t) is sufficiently continuously differentiable for t≥0. Let
(2.31)Gβ,1N(t)=ddtI~β,NU(t)-I~β,NddtU(t).
Then we obtain that
(2.32)ddtI~β,NU(tβ,kN)=τ^f(U(tβ,kN),U(tβ,kN-tβ,NN))+Gβ,1N(tβ,kN),1≤k≤N.
Now we derive an explicit expression for the left side of (2.32). Let U~β,lN be the coefficients of I~β,NU(t) in terms of ℒ~l(β)(t). Due to (2.15), we have
(2.33)ddtI~β,NU(t)=∑l=0N‍U~β,lNddtℒ~l(β)(t)=-12β∑l=1N‍U~β,lN(2∑m=0l-1‍ℒ~m(β)(t)+ℒ~l(β)(t))-12βU~β,0Nℒ~0(β)(t).
This equality and (2.24) imply that
(2.34)ddtI~β,NU(tβ,kN)=-12β2∑l=1N‍(∑j=0N‍U(tβ,jN)ℒ~l(β)(tβ,jN)ω~β,jN)(2∑m=0l-1‍ℒ~m(β)(tβ,kN)+ℒ~l(β)(tβ,kN))-12β2(∑j=0N‍U(tβ,jN)ℒ~0(β)(tβ,jN)ω~β,jN)ℒ~0(β)(tβ,kN).
Denote for 0≤j≤N and 1≤k≤N(2.35)aβ,k,jN=-12β2ω~β,jN(∑l=1N‍ℒ~l(β)(tβ,jN)(2∑m=0l-1‍ℒ~m(β)(tβ,kN)+ℒ~l(β)(tβ,kN))+ℒ~0(β)(tβ,jN)ℒ~0(β)(tβ,kN)).
Then
(2.36)ddtI~β,NU(tβ,kN)=∑j=0N‍aβ,k,jNU(tβ,jN).
Denote
(2.37)𝕌N=(U(0),U(tβ,1N),…,U(tβ,NN))T,𝔽βN(𝕌N)=(τ^f(U(tβ,0N),ϕ0),τ^f(U(tβ,1N),ϕ1),…,τ^f(U(tβ,NN),ϕN))T,ϕi=Φ(τ^(tβ,iN-tβ,NN)),i=0,1,…,N,𝔾β,1N=(Gβ,1N(tβ,1N),Gβ,1N(tβ,2N),…,Gβ,1N(tβ,NN))T,𝔸βN=(aβ,1,0Naβ,1,1N⋯aβ,1,NNaβ,2,0Naβ,2,1N⋯aβ,2,NN⋮⋮⋯⋮⋮⋮⋯⋮aβ,N,0Naβ,N,1N⋯aβ,N,NN).
Then, we obtain
(2.38)𝔸βN𝕌N=𝔽βN(𝕌ℕ)+𝔾β,1N,U(0)=Φ(0).
We now approximate U(t) by uN(t)∈QN(0,∞). Clearly, I~β,NuN(t)=uN(t). Furthermore, we set
(2.39)uN=(uN(0),uN(tβ,1N),…,uN(tβ,NN))T,𝔽βN(uN)=(τ^f(uN(tβ,0N),ϕ0),τ^f(uN(tβ,1N),ϕ1),…,τ^f(uN(tβ,NN),ϕN))T.
By replacing 𝕌N by uN and neglecting 𝔾β,1N in (2.38), we derive a new integration process by using the modified Laguerre functions. It is to seek uN such that
(2.40)𝔸βNuN=𝔽βN(uN),uN(0)=Φ(0).
The global numerical solution is
(2.41)uN(t)=∑l=0N‍u~β,lNℒ~l(β)(t),t≥0,
with
(2.42)u~β,lN=β(uN,ℒ~l(β))β,N=β∑j=0N‍uN(tβ,jN)ℒ~l(β)(tβ,jN)ω~β,jN.
Let
(2.43)u^N(t)={uN(t)∈QN(0,∞),t∈(0,∞),Φ(tτ^),t∈[-tβ,NN,0].
Then (2.40) is equivalent to the system
(2.44)ddtu^N(t)=τ^f(u^N(t),u^N(t-tβ,NN)),t>0,u^N(t)=Φ(tτ^),-tβ,NN≤t≤0.


## 3.  Convergence Analysis

In this section, we estimate the global error of numerical solution. For anyrth continuously differentiable function v(t), we set
(3.2)ℛN,r,β(1)(v)=β-1∥t(r-1)/2drvdtr∥ωβ+(1+β-(1/2))(lnN)1/2∥tr/2drvdtr∥ωβ,ℛN,r,β(2)(v)=β-1∥t(r+1)/2dr+2vdtr+2∥ωβ+N-1/2∥t(r+1)/2dr+2vdtr+2∥ωβ+(1+β-1/2)(lnN)1/2∥t(r+1)/2dr+2vdtr+2∥ωβ.
The following lemmas will play a key role in obtaining our main results.Lemma 3.1 (see [19]).
Ifv∈Lωβ2(0,∞), then for an integer r≥1,
(3.3)∥Iβ,Nv-v∥ωβ≤c(βN)1/2-r/2ℛN,r,β(1)(v).Lemma 3.2 (see [19]).
Ifv∈Lωβ2(0,∞), then for an integer r≥1,
(3.4)∥ddt(Iβ,Nv-v)∥ωβ≤c(βN)1/2-r/2ℛN,r,β(2)(v).Theorem 3.3.
Suppose that there exists a real numberγ0>0 such that
(3.5)(f(u1,v)-f(u2,v))(u1-u2)≤-γ0|u1-u2|2,∀u1,u2∈ℝ,
and ℛN,r,β(1)(U), ℛN,r,β(2)(U), and ℛN,r,β(1)(dU/dt) are finite. Then
(3.6)∥U-u^N∥≤cγ0τ^(βN)1/2-r/2×((γ0τ^+β)ℛN,r,β(1)(e(1/2)βtU)+ℛN,r,β(2)(e(1/2)βtU)+ℛN,r,β(1)(e(1/2)βtdUdt)).Proof.
LetEN(t)=u^N(t)-I~β,NU(t). Subtracting (2.32) from (2.44), we get
(3.7)ddtEN(tβ,kN)=Gβ,2N(tβ,kN)-Gβ,1N(tβ,kN),1≤k≤N,EN(0)=0,
where
(3.8)Gβ,2N(tβ,kN)=τ^[f(uN(tβ,kN),uN(tβ,kN-tβ,NN))-f(Iβ,NU(tβ,kN),U(tβ,kN-tβ,NN))].
We now multiply (3.7) by 2EN(tβ,kN)ω~β,kN and sum the result for 1≤k≤N. Due to EN(0)=0, we obtain that
(3.9)2〈EN,ddtEN〉β,N=Aβ,1N+Aβ,2N,
where
(3.10)Aβ,1N=-2〈Gβ,1N,EN〉β,N,Aβ,2N=2〈Gβ,2N,EN〉β,N.
Using (2.20) and the Cauchy inequality, we deduce that
(3.11)2〈EN,ddtEN〉β,N=2〈EN,ddtEN〉=|EN(∞)|2,|Aβ,1N|≤2∥Gβ,1N∥β,N∥EN∥β,N=2∥Gβ,1N∥∥EN∥.
Thus (3.9) reads
(3.12)|EN(∞)|2≤Aβ,2N+2∥Gβ,1N∥∥EN∥.
Since there exists γ0>0 such that
(3.13)(f(u1,v)-f(u2,v))(u1-u2)≤-γ0|u1-u2|2,∀u1,u2∈ℝ,Aβ,2N≤-2γ-0∥EN∥2,
where γ-0=γ0τ^. Therefore
(3.14)|EN(∞)|2+γ-0∥EN∥2≤1γ-0∥Gβ,1N∥2.
Hence it suffices to estimate ∥Gβ,1N∥2. With the aid of (2.26), Lemmas 3.1 and 3.2, we deduce that, for r≥1,
(3.15)∥I~β,Nv-v∥=∥Iβ,N(e(1/2)βtv)-e(1/2)βtv∥ωβ≤c(βN)1/2-r/2ℛN,r,β(1)(e(1/2)βtv).
On the other hand,
(3.16)ddt(I~β,Nv-v)=-12βe-(1/2)βt(Iβ,N(e(1/2)βtv)-e(1/2)βtv)+e-(1/2)βtddt(Iβ,N(e(1/2)βtv)-e(1/2)βtv).
It follows from the above result that for r≥1(3.17)∥ddt(I~β,Nv-v)∥≤c(βN)1/2-r/2(βℛN,r,β(1)(e(1/2)βtv)+ℛN,r,β(2)(e(1/2)βtv)).
Consequently,
(3.18)∥Gβ,1N∥≤∥ddt(I~β,NU-U)∥+∥ddtU-I~β,NddtU∥≤c(βN)1/2-r/2(βℛN,r,β(1)(e(1/2)βtU)+ℛN,r,β(2)(e(1/2)βtU)+ℛN,r,β(1)(e(1/2)βtdUdt)).
Thus, (3.14) implies that
(3.19)|EN(∞)|+γ-01/2∥EN∥≤cγ-0(1/2)(βN)(1/2)-(r/2)(βℛN,r,β(1)(e(1/2)βtU)+ℛN,r,β(2)(e(1/2)βtU)+ℛN,r,β(1)(e(1/2)βtdUdt)).
Furthermore, using (3.15) again, we obtain that, for any β>0,
(3.20)∥U-u^N∥≤cγ-0(βN)(1/2)-(r/2)×((γ-0+β)ℛN,r,β(1)(e(1/2)βtU)+ℛN,r,β(2)(e(1/2)βtU)+ℛN,r,β(1)(e(1/2)βtdUdt)),|U(∞)-u^N(∞)|≤|I~β,NU(∞)-U(∞)|+|EN(∞)|≤2∥I~β,NU-U∥1/2∥I~β,NU-U∥11/2+|EN(∞)|≤c(βN)1/2-r/2((βγ-0-1/2+β+1)ℛN,r,β(1)(e(1/2)βtU)+(1+γ-0-1/2)ℛN,r,β(2)(e(1/2)βtU)+γ-0-1/2ℛN,r,β(1)(e(1/2)βtdUdt)).
This completes the proof.Remark 3.4.
The norm∥U∥ is finite as long as that f(u,v) satisfies certain conditions. If f(u,v) satisfies conditions
(3.21)〈f(u1,v)-f(u2,v),u1-u2〉≤γ1∥u1-u2∥2,∀u1,u2∈ℝN,∥f(u,v1)-f(u,v2)∥≤γ2∥v1-v2∥,∀v1,v2∈ℝN,γ1<0,0<γ2<-γ1,
and f(0,0)=0, then |U(t)|=O(e-γ*t), γ*>0, see Tian [5]. Furthermore, if f(u,v) fulfills some additional conditions, then the norms appearing at the right sides of (3.20) are finite. Therefore, for certain positive constant c* depending only on β,
(3.22)∥U-u^N∥+|U(∞)-u^N(∞)|=c*(1+1γ0)(lnN)1/2N1/2-r/2.
Consequently, for r>1, the scheme (2.40) has the global convergence and the spectral accuracy in L2(0,∞). Moreover, at the infinity, the numerical solution has the same accuracy. This also indicates that the pointwise numerical error decays rapidly as the mode N increases, with the convergence rate as c*(lnN)1/2N1/2-r/2. On the other hand, for any fixed N, the norm ∥U-u^N∥ is bounded, and so the pointwise numerical error decays automatically as t→∞, at least less than cNt-1/2, cN being a small number. Hence, it is very efficient for long-time numerical simulations of dynamical systems, especially for stiff problems.Remark 3.5.
The method proposed is also available for solving systems of delay differential equations. Let(3.23)U→(t)=(U(1)(t),U(2)(t),…,U(m)(t))T,V→(t)=(V(1)(t),V(2)(t),…,V(m)(t))T,f→(U→,V→)=(f(1)(U→,V→),f(2)(U→,V→),…,f(m)(U→,V→))T.
We consider the system
(3.24)ddtU→(t)=f→(U→(t),V→(t)),t>0,U→(t)=Φ→(t),-τ≤t≤0.
We approximate U→(t) by u→N(t). We can derive a numerical algorithm which is similar to (2.38). Further, let |V→|E be the Euclidean norm of V→. Assume that
(3.25)(f→(Z→1,V→)-f→(Z→2,V→))·(Z→1-Z→2)≤-γ0|Z→1-Z→2|E2,γ0>0.
Then we can obtain an error estimate similar to (3.20).


## 4. Refinement and Numerical Results

In the preceding sections, we introduced an integration process for solving delay differential equations. Theoretically, their numerical errors with biggerN decrease faster. But in practical computation, it is not convenient to use them with very big N. On the other hand, the distance between the adjacent interpolation nodes tβ,jN and tβ,j-1N increases fast as N and j increase, especially for the nodes which are located far from the origin point t=0. This is one of advantages of the Laguerre interpolation approximation, since we can use moderate mode N to evaluate the unknown function at large t, but it is also its shortcoming. In fact, if the exact solution oscillates or changes very rapidly between two large adjacent interpolation nodes, then we may lose information about the structure of exact solution between those nodes. To remedy this deficiency, we may refine the numerical results as follows.LetN be a moderate positive integer, β>0, and the set of nodes {t0,β,jN}j=0N={tβ,jN}j=0N. We use (2.40) with the interpolation nodes {t0,β,jN}j=0N to obtain the original numerical solution
(4.1)u(0,N)(t)=uN(t),0≤t≤t0,β,NN.
Then we take t1,β,0N=t0,β,NN and consider the following delay differential equation:
(4.2)ddtU(1)(t)=τ^f(U(1)(t),U(1)(t-tβ,NN)),t>t1,β,0N,U(1)(t)=u(0,N)(t),t1,β,0N-τ≤t≤t1,β,0N.
We get the refined numerical solution u(1,N)(t) for t1,β,0N≤t≤t1,β,NN. Repeating the above procedure, we obtain the refined numerical solution u(m,N)(t) for tm,β,0N≤t≤tm,β,NN.


## 5. Numerical Results

We consider the following system of four homogeneous delay differential equations:(5.1)y˙1(t)=y3(t),y˙2(t)=y4(t),y˙3(t)=-2ny2(t)+(1+n2)(-1)ny1(t-π),y˙4(t)=-2ny1(t)+(1+n2)(-1)ny2(t-π).
The initial functions and solutions are given by
(5.2)y1(t)=sin(t)cos(nt),y2(t)=cos(t)sin(nt),y3(t)=y˙1(t),y4(t)=y˙2(t),t∈[-π,∞).
For n=1, Figures 1 and 2 show plots of the exact solution y1(t) and y3(t) on [0,2π], which controls the frequency of oscillation in the initial data and solution. Figures 3 and 4 are the exact and numerical solution on [0,π] with N=15, β=10. The refined numerical results are given in Figures 5 and 6 on [0,3π].Figure 1
Exact solution.Figure 2
Exact solution.Figure 3
Numerical solution and exact solution.Figure 4
Numerical solution and exact solution.Figure 5
Numerical solution of refined method and exact solution.Figure 6
Numerical solution of refined method and exact solution.The numerical experiments show that our numerical integration processes are efficient for numerically solving delay differential equations.


## 6. Conclusions

In this paper, we proposed new integration processes of delay differential equations, which have fascinating advantages. On the one hand, the suggested integration processes are based on the modified Laguerre functions on the half line; they provide the global numerical solution and the global convergence naturally and thus are available for long-time numerical simulations of dynamical systems. On the other hand, benefiting from the rapid convergence of the modified Laguerre functions, these processes possess the spectral accuracy. In particular, the numerical results fit the exact solutions well at the interpolation nodes. Furthermore, We also developed a technique for refinement of modified Laguerre-Radau interpolations. Lastly, numerical results demonstrate the spectral accuracy of the proposed method and coincide well with analysis.



---

*Source: 978729-2012-09-23.xml*