Integration Processes of Delay Differential Equation Based on Modified Laguerre Functions

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.

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 𝜔_𝛽(𝑡)=𝑒^−𝛽𝑡, 𝛽>0, and define the weighted space 𝐿²_𝜔𝛽(0,∞) as usual, with the following inner product and norm: (𝑢,𝑣)_𝜔𝛽=^∞₀𝑢(𝑡)𝑣(𝑡)𝜔_𝛽(𝑡)𝑑𝑡,‖𝑣‖_𝜔𝛽=(𝑣,𝑣)_𝜔^1/2_𝛽. (2.1) The modified Laguerre polynomial of degree 𝑙 is defined by (cf. [18]) ℒ_𝑙^(𝛽)(𝑡)=1𝑒𝑙!^𝛽𝑡𝑑^𝑙𝑑𝑡^𝑙𝑡^𝑙𝑒^−𝛽𝑡,𝑙≥0. (2.2) For example, ℒ₀^(𝛽)(𝑡)ℒ=1,_1(𝛽)(𝑡)ℒ=1−𝛽𝑡,_2(𝛽)(𝑡)𝛽=−2𝛽𝑡+²2𝑡². (2.3) The modified Laguerre polynomials satisfy the recurrence relation 𝑑ℒ𝑑𝑡^𝑙(𝛽)(𝑡)=𝑑ℒ𝑑𝑡^{(𝛽)𝑙−1}(𝑡)−𝛽ℒ^(𝛽)_𝑙−1(𝑡),𝑙≥1. (2.4) The set of Laguerre polynomials is a complete 𝐿²_𝜔𝛽(0,∞)-orthogonal system, namely, ℒ_𝑙^(𝛽),ℒ_𝑚^(𝛽)_𝜔𝛽=1𝛽𝛿_,𝑚, (2.5) where 𝛿_𝑙,𝑚 is the Kronecker function. Thus, for any 𝑣∈𝐿²_𝜔𝛽(0,∞), 𝑣(𝑡)=^∞_𝑙=0𝑣_𝑙ℒ_𝑙^(𝛽)(𝑡),𝑣_𝑙=𝛽𝑣,ℒ_𝑙^(𝛽)_𝜔𝛽. (2.6) Let 𝑁 be any positive integer and 𝒫_𝑁(0,∞) the set of all algebraic polynomials of degree at most 𝑁. We denote by 𝑡^𝑁_𝛽,𝑗 the modified Laguerre-Radau interpolation points. Indeed, 𝑡^𝑁_𝛽,0=0, and 𝑡^𝑁_𝛽,𝑗(1≤𝑗≤𝑁) are the distinct zeros of (𝑑/𝑑𝑡)ℒ^(𝛽)_𝑁+1(𝑡). By using (2.1) and the formula (2.12) of [19], the corresponding Christoffel numbers are as follows: 𝜔^𝑁_𝛽,0=1𝛽(𝑁+1),𝜔^𝑁_𝛽,𝑗=1𝛽(𝑁+1)ℒ_𝑁^(𝛽)𝑡^𝑁_𝛽,𝑗ℒ^(𝛽)_𝑁+1𝑡^𝑁_𝛽,𝑗,1≤𝑗≤𝑁. (2.7) For any 𝜌∈𝒫_2𝑁(0,∞), ^𝑁_𝑗=0𝜌𝑡^𝑁_𝛽,𝑗𝜔^𝑁_𝛽,𝑗=^∞₀𝜌(𝑡)𝜔_𝛽(𝑡)𝑑𝑡. (2.8) Next, we define the following discrete inner product and norm: (𝑢,𝑣)_{𝜔𝛽,𝑁}=^𝑁_𝑗=0𝑢𝑡^𝑁_𝛽,𝑗𝑣𝑡^𝑁_𝛽,𝑗𝜔^𝑁_𝛽,𝑗,‖𝑣‖_{𝜔𝛽,𝑁}=(𝑣,𝑣)_𝜔^1/2_𝛽,𝑁. (2.9) For any 𝜙,𝜓∈𝒫_𝑁(0,∞), 𝜙,𝜓_𝜔𝛽=𝜙,𝜓_{𝜔𝛽,𝑁},𝜙_𝜔𝛽=𝜙_{𝜔𝛽,𝑁}. (2.10) For all 𝑣∈𝐿²_𝜔𝛽(0,∞), the modified Laguerre-Radau interpolant 𝐼_𝛽,𝑁𝑣∈𝒫_𝑁(0,∞) is determined by 𝐼_𝛽,𝑁𝑣𝑡^𝑁_𝛽,𝑗=𝑣𝑡^𝑁_𝛽,𝑗,0≤𝑗≤𝑁. (2.11) By (2.10), for any 𝜙∈𝒫_𝑁(0,∞), 𝐼_𝛽,𝑁𝑣,𝜙_𝜔𝛽=𝐼_𝛽,𝑁𝑣,𝜙_{𝜔𝛽,𝑁}=𝑣,𝜙_{𝜔𝛽,𝑁}. (2.12) The interpolant 𝐼_𝛽,𝑁𝑣 can be expanded as 𝐼_𝛽,𝑁𝑣(𝑡)=^𝑁_𝑙=0𝑣^𝑁_𝛽,𝑙ℒ_𝑙^(𝛽)(𝑡). (2.13) By virtue of (2.5) and (2.10), 𝑣^𝑁_𝛽,𝑙=𝛽𝐼_𝛽,𝑁𝑣,ℒ_𝑙^(𝛽)_𝜔𝛽=𝛽𝑣,ℒ_𝑙^(𝛽)_{𝜔𝛽,𝑁}. (2.14) Define the modified Laguerre functions ℒ_𝑙^(𝛽)(𝑡)=𝑒^{−(1/2)𝛽𝑡}ℒ_𝑙^(𝛽)(𝑡) as the base functions. According to (2.4), the functions ℒ_𝑙^(𝛽)(𝑡) satisfy the recurrence relation 𝑑ℒ𝑑𝑡^𝑙(𝛽)(𝑡)=𝑑ℒ𝑑𝑡^{(𝛽)𝑙−1}(𝑡)−12𝛽ℒ_𝑙^(𝛽)(𝑡)−12𝛽ℒ^(𝛽)_𝑙−1(𝑡),𝑙≥1. (2.15) Denote by (𝑢,𝑣) and ‖𝑣‖ the inner product and the norm of the space 𝐿²(0,∞), respectively. The set of ℒ_𝑙^(𝛽)(𝑡) is a complete 𝐿²(0,∞)-orthogonal system, that is, ℒ_𝑙^(𝛽),ℒ_𝑚^(𝛽)=1𝛽𝛿_𝑙,𝑚. (2.16) We now introduce the new Laguerre-Radau interpolation. Set 𝑄_𝑁(0,∞)=spanℒ₀^(𝛽),ℒ₁^(𝛽),…,ℒ_𝑁^(𝛽). (2.17) Let 𝑡^𝑁_𝛽,𝑗 and 𝜔^𝑁_𝛽,𝑗 be the same as in (2.7), and take the nodes and weights of the new Laguerre-Radau interpolation as 𝑡̃^𝑁_𝛽,𝑗=𝑡^𝑁_𝛽,𝑗,𝜔^𝑁_𝛽,𝑗=1ℒ_𝑁^(𝛽)𝑡^𝑁_𝛽,𝑗ℒ^(𝛽)_𝑁+1𝑡^𝑁_𝛽,𝑗=𝑒^{𝛽𝑡𝑁𝛽,𝑗}𝜔^𝑁_𝛽,𝑗. (2.18) The discrete inner product and norm can be defined similarly as ⟨𝑢,𝑣⟩_𝛽,𝑁=^𝑁_𝑗=0𝑢𝑡^𝑁_𝛽,𝑗𝑣𝑡^𝑁_𝛽,𝑗𝜔^𝑁_𝛽,𝑗,‖𝑣‖_𝛽,𝑁=⟨𝑣,𝑣⟩^1/2_𝛽,𝑁. (2.19) For any 𝜙₁, 𝜙₂∈𝑄_𝑁(0,∞), we have 𝜙₁=𝑒^{−(1/2)𝛽𝑡}𝜓₁,  𝜙₂=𝑒^{−(1/2)𝛽𝑡}𝜓₂, and 𝜓₁, 𝜓₂∈𝒫_𝑁(0,∞). Thus by (2.10), 𝜙₁,𝜙_2𝛽,𝑁=𝜓₁,𝜓_{2𝜔𝛽,𝑁}=𝜓₁,𝜓_2𝜔𝛽=𝜙₁,𝜙₂. (2.20) The new Laguerre-Radau interpolant 𝐼̃_𝛽,𝑁𝑣∈𝑄_𝑁(0,∞) is determined by 𝐼̃_𝛽,𝑁𝑣𝑡^𝑁_𝛽,𝑗=𝑣𝑡^𝑁_𝛽,𝑗,0≤𝑗≤𝑁. (2.21) Due to the equality (2.20), for any 𝜙∈𝑄_𝑁(0,∞), 𝐼̃_𝛽,𝑁𝑣,𝜙=𝐼̃_𝛽,𝑁𝑣,𝜙_𝛽,𝑁=𝑣,𝜙_𝛽,𝑁. (2.22) Let 𝐼̃_𝛽,𝑁𝑣(𝑡)=^𝑁_𝑙=0𝑣^𝑁_𝛽,𝑙ℒ_𝑙^(𝛽)(𝑡). (2.23) Then, with the aid of (2.16) and (2.22), we derive that 𝑣^𝑁_𝛽,𝑙=𝛽𝐼̃_𝛽,𝑁𝑣,ℒ_𝑙^(𝛽)=𝛽𝐼̃_𝛽,𝑁𝑣,ℒ_𝑙^(𝛽)_𝛽,𝑁=𝛽𝑣,ℒ_𝑙^(𝛽)_𝛽,𝑁. (2.24) There is a close relation between 𝐼_𝛽,𝑁 and 𝐼̃_𝛽,𝑁. From the previous two equalities, it follows that 𝑒^{(1/2)𝛽𝑡}𝐼̃_𝛽,𝑁𝑣(𝑡)=^𝑁_𝑙=0𝑣^𝑁_𝛽,𝑙ℒ_𝑙^(𝛽)(𝑡)=𝛽^𝑁_𝑙=0𝑣,ℒ_𝑙^(𝛽)_𝛽,𝑁ℒ_𝑙^(𝛽)(𝑡)=𝛽^𝑁_𝑙=0𝑒^{(1/2)𝛽𝑡}𝑣,ℒ_𝑙^(𝛽)_{𝜔𝛽,𝑁}ℒ_𝑙^(𝛽)(𝑡). (2.25) This with (2.14) implies 𝐼̃_𝛽,𝑁𝑣(𝑡)=𝑒^{−(1/2)𝛽𝑡}𝐼_𝛽,𝑁𝑒^{(1/2)𝛽𝑡}𝑣(𝑡). (2.26) Consider the following delay differential equation: 𝑑𝑊(𝑡)(𝑊(𝑡)(𝑊(𝑡)(𝑡)𝑑𝑡=𝑓,𝑊𝑡−𝜏)),𝑡>0,=Φ,−𝜏≤𝑡≤0. (2.27) For any fixed positive integer 𝑁, we define 𝑈(𝑡)=𝑊𝑡𝜏𝑡^𝑁_𝛽,𝑁, (2.28) and denote 𝜏𝜏=𝑡^𝑁_𝛽,𝑁. (2.29) Then system of (2.27) can be transformed into 𝑑𝑈(𝑡)=𝑈(𝑡)𝑑𝑡𝜏𝑓,𝑈𝑡−𝑡^𝑁_𝛽,𝑁𝑈(𝑡)(𝑡𝜏),𝑡>0,=Φ,−𝑡^𝑁_𝛽,𝑁≤𝑡≤0. (2.30) We suppose that 𝑈(𝑡) is sufficiently continuously differentiable for 𝑡≥0. Let 𝐺^𝑁_𝛽,1(𝑡)=𝑑̃𝐼𝑑𝑡^𝛽,𝑁𝑈(𝑡)−𝐼̃_𝛽,𝑁𝑑𝑈(𝑡).𝑑𝑡 (2.31) Then we obtain that 𝑑̃𝐼𝑑𝑡^𝛽,𝑁𝑈𝑡^𝑁_𝛽,𝑘=𝜏𝑓𝑈𝑡^𝑁_𝛽,𝑘,𝑈𝑡^𝑁_𝛽,𝑘−𝑡^𝑁_𝛽,𝑁+𝐺^𝑁_𝛽,1𝑡^𝑁_𝛽,𝑘,1≤𝑘≤𝑁. (2.32) Now we derive an explicit expression for the left side of (2.32). Let 𝑈^𝑁_𝛽,𝑙 be the coefficients of 𝐼̃_𝛽,𝑁𝑈(𝑡) in terms of ℒ_𝑙^(𝛽)(𝑡). Due to (2.15), we have 𝑑̃𝐼𝑑𝑡^𝛽,𝑁𝑈(𝑡)=^𝑁_𝑙=0𝑈^𝑁_𝛽,𝑙𝑑ℒ𝑑𝑡^𝑙(𝛽)(𝑡)1=−2𝛽^𝑁_𝑙=1𝑈^𝑁_𝛽,𝑙2^𝑙−1_𝑚=0ℒ_𝑚^(𝛽)(𝑡)+ℒ_𝑙^(𝛽)(𝑡)−12𝛽𝑈^𝑁_𝛽,0ℒ₀^(𝛽)(𝑡). (2.33) This equality and (2.24) imply that 𝑑̃𝐼𝑑𝑡^𝛽,𝑁𝑈𝑡^𝑁_𝛽,𝑘1=−2𝛽^2𝑁_𝑙=1^𝑁_𝑗=0𝑈𝑡^𝑁_𝛽,𝑗ℒ_𝑙^(𝛽)𝑡^𝑁_𝛽,𝑗𝜔^𝑁_𝛽,𝑗2^𝑙−1_𝑚=0ℒ_𝑚^(𝛽)𝑡^𝑁_𝛽,𝑘+ℒ_𝑙^(𝛽)𝑡^𝑁_𝛽,𝑘−12𝛽^2𝑁_𝑗=0𝑈𝑡^𝑁_𝛽,𝑗ℒ₀^(𝛽)𝑡^𝑁_𝛽,𝑗𝜔^𝑁_𝛽,𝑗ℒ₀^(𝛽)𝑡^𝑁_𝛽,𝑘. (2.34) Denote for 0≤𝑗≤𝑁 and 1≤𝑘≤𝑁𝑎^𝑁_{𝛽,𝑘,𝑗}1=−2𝛽²𝜔^𝑁_𝛽,𝑗^𝑁_𝑙=1ℒ_𝑙^(𝛽)𝑡^𝑁_𝛽,𝑗2^𝑙−1_𝑚=0ℒ_𝑚^(𝛽)𝑡^𝑁_𝛽,𝑘+ℒ_𝑙^(𝛽)𝑡^𝑁_𝛽,𝑘+ℒ₀^(𝛽)𝑡^𝑁_𝛽,𝑗ℒ₀^(𝛽)𝑡^𝑁_𝛽,𝑘. (2.35) Then 𝑑̃𝐼𝑑𝑡^𝛽,𝑁𝑈𝑡^𝑁_𝛽,𝑘=^𝑁_𝑗=0𝑎^𝑁_{𝛽,𝑘,𝑗}𝑈𝑡^𝑁_𝛽,𝑗. (2.36) Denote 𝕌^𝑁=𝑈(0),𝑈𝑡^𝑁_𝛽,1,…,𝑈𝑡^𝑁_𝛽,𝑁^𝑇,𝔽^𝑁_𝛽𝕌^𝑁=𝜏𝑓𝑈𝑡^𝑁_𝛽,0,𝜙₀,𝜏𝑓𝑈𝑡^𝑁_𝛽,1,𝜙₁,…,𝜏𝑓𝑈𝑡^𝑁_𝛽,𝑁,𝜙_𝑁^𝑇,𝜙_𝑖=Φ𝜏𝑡^𝑁_𝛽,𝑖−𝑡^𝑁_𝛽,𝑁𝔾,𝑖=0,1,…,𝑁,_𝑁𝛽,1=𝐺^𝑁_𝛽,1𝑡^𝑁_𝛽,1,𝐺^𝑁_𝛽,1𝑡^𝑁_𝛽,2,…,𝐺^𝑁_𝛽,1𝑡^𝑁_𝛽,𝑁^𝑇,𝔸^𝑁_𝛽=𝑎^𝑁_𝛽,1,0𝑎^𝑁_𝛽,1,1⋯𝑎^𝑁_𝛽,1,𝑁𝑎^𝑁_𝛽,2,0𝑎^𝑁_𝛽,2,1⋯𝑎^𝑁_𝛽,2,𝑁𝑎⋮⋮⋯⋮⋮⋮⋯⋮_{𝑁𝛽,𝑁,0}𝑎^𝑁_𝛽,𝑁,1⋯𝑎^𝑁_{𝛽,𝑁,𝑁}. (2.37) Then, we obtain 𝔸^𝑁_𝛽𝕌^𝑁=𝔽^𝑁_𝛽𝕌^ℕ+𝔾^𝑁_𝛽,1,𝑈(0)=Φ(0). (2.38) We now approximate 𝑈(𝑡) by 𝑢^𝑁(𝑡)∈𝑄_𝑁(0,∞). Clearly, 𝐼̃_𝛽,𝑁𝑢^𝑁(𝑡)=𝑢^𝑁(𝑡). Furthermore, we set 𝐮^𝑁=𝑢^𝑁(0),𝑢^𝑁𝑡^𝑁_𝛽,1,…,𝑢^𝑁𝑡^𝑁_𝛽,𝑁^𝑇,𝔽^𝑁_𝛽𝐮^𝑁=𝜏𝑓𝑢^𝑁𝑡^𝑁_𝛽,0,𝜙₀,𝜏𝑓𝑢^𝑁𝑡^𝑁_𝛽,1,𝜙₁,…,𝜏𝑓𝑢^𝑁𝑡^𝑁_𝛽,𝑁,𝜙_𝑁^𝑇. (2.39) By replacing 𝕌^𝑁 by 𝐮^𝑁 and neglecting 𝔾^𝑁_𝛽,1 in (2.38), we derive a new integration process by using the modified Laguerre functions. It is to seek 𝐮^𝑁 such that 𝔸^𝑁_𝛽𝐮^𝑁=𝔽^𝑁_𝛽𝐮^𝑁,𝑢^𝑁(0)=Φ(0). (2.40) The global numerical solution is 𝑢^𝑁(𝑡)=^𝑁_𝑙=0𝑢^𝑁_𝛽,𝑙ℒ_𝑙^(𝛽)(𝑡),𝑡≥0, (2.41) with 𝑢^𝑁_𝛽,𝑙=𝛽𝑢^𝑁,ℒ_𝑙^(𝛽)_𝛽,𝑁=𝛽^𝑁_𝑗=0𝑢^𝑁𝑡^𝑁_𝛽,𝑗ℒ_𝑙^(𝛽)𝑡^𝑁_𝛽,𝑗𝜔^𝑁_𝛽,𝑗. (2.42) Let 𝑢^𝑁(𝑡)=𝑢^𝑁(𝑡)∈𝑄_𝑁()(),Φ(𝑡𝜏)0,∞,𝑡∈0,∞,𝑡∈−𝑡^𝑁_𝛽,𝑁,0. (2.43) Then (2.40) is equivalent to the system 𝑑𝑢𝑑𝑡^𝑁(𝑡)=𝜏𝑓𝑢^𝑁(𝑡),𝑢^𝑁𝑡−𝑡^𝑁_𝛽,𝑁𝑢,𝑡>0,_𝑁(𝑡)=Φ(𝑡𝜏),−𝑡^𝑁_𝛽,𝑁≤𝑡≤0. (2.44)

3. Convergence Analysis

In this section, we estimate the global error of numerical solution. For any 𝑟th continuously differentiable function 𝑣(𝑡), we set ℛ⁽¹⁾_{𝑁,𝑟,𝛽}(𝑣)=𝛽⁻¹𝑡^(𝑟−1)/2𝑑^𝑟𝑣𝑑𝑡^𝑟_𝜔𝛽+1+𝛽^−(1/2)(ln𝑁)^1/2𝑡^𝑟/2𝑑^𝑟𝑣𝑑𝑡^𝑟_𝜔𝛽,ℛ⁽²⁾_{𝑁,𝑟,𝛽}(𝑣)=𝛽⁻¹𝑡^(𝑟+1)/2𝑑^𝑟+2𝑣𝑑𝑡^𝑟+2_𝜔𝛽+𝑁^−1/2𝑡^(𝑟+1)/2𝑑^𝑟+2𝑣𝑑𝑡^𝑟+2_𝜔𝛽+1+𝛽^−1/2(ln𝑁)^1/2𝑡^(𝑟+1)/2𝑑^𝑟+2𝑣𝑑𝑡^𝑟+2_𝜔𝛽. (3.2) The following lemmas will play a key role in obtaining our main results.

Lemma 3.1 (see [19]). If 𝑣∈𝐿²_𝜔𝛽(0,∞), then for an integer 𝑟≥1, 𝐼_𝛽,𝑁𝑣−𝑣_𝜔𝛽≤𝑐𝛽𝑁^1/2−𝑟/2ℛ⁽¹⁾_{𝑁,𝑟,𝛽}(𝑣). (3.3)

Lemma 3.2 (see [19]). If 𝑣∈𝐿²_𝜔𝛽(0,∞), then for an integer 𝑟≥1, 𝑑𝐼𝑑𝑡^𝛽,𝑁𝑣−𝑣_𝜔𝛽≤𝑐𝛽𝑁^1/2−𝑟/2ℛ⁽²⁾_{𝑁,𝑟,𝛽}(𝑣). (3.4)

Theorem 3.3. Suppose that there exists a real number 𝛾₀>0 such that 𝑓𝑢₁,𝑣−𝑓𝑢₂,𝑣𝑢₁−𝑢₂≤−𝛾₀𝑢₁−𝑢₂²,∀𝑢₁,𝑢₂∈ℝ, (3.5) and ℛ⁽¹⁾_{𝑁,𝑟,𝛽}(𝑈), ℛ⁽²⁾_{𝑁,𝑟,𝛽}(𝑈), and ℛ⁽¹⁾_{𝑁,𝑟,𝛽}(𝑑𝑈/𝑑𝑡) are finite. Then 𝑈−𝑢^𝑁≤𝑐𝛾₀𝜏𝛽𝑁^1/2−𝑟/2×𝛾₀𝜏+𝛽ℛ⁽¹⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑈+ℛ⁽²⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑈+ℛ⁽¹⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑑𝑈.𝑑𝑡 (3.6)

Proof. Let 𝐸^𝑁(𝑡)=𝑢^𝑁(𝑡)−𝐼̃_𝛽,𝑁𝑈(𝑡). Subtracting (2.32) from (2.44), we get 𝑑𝐸𝑑𝑡^𝑁𝑡^𝑁_𝛽,𝑘=𝐺^𝑁_𝛽,2𝑡^𝑁_𝛽,𝑘−𝐺^𝑁_𝛽,1𝑡^𝑁_𝛽,𝑘𝐸,1≤𝑘≤𝑁,_𝑁(0)=0, (3.7) where 𝐺^𝑁_𝛽,2𝑡^𝑁_𝛽,𝑘=𝜏𝑓𝑢^𝑁𝑡^𝑁_𝛽,𝑘,𝑢^𝑁𝑡^𝑁_𝛽,𝑘−𝑡^𝑁_𝛽,𝑁−𝑓𝐼_𝛽,𝑁𝑈𝑡^𝑁_𝛽,𝑘,𝑈𝑡^𝑁_𝛽,𝑘−𝑡^𝑁_𝛽,𝑁. (3.8) We now multiply (3.7) by 2𝐸^𝑁(𝑡^𝑁_𝛽,𝑘)𝜔^𝑁_𝛽,𝑘 and sum the result for 1≤𝑘≤𝑁. Due to 𝐸^𝑁(0)=0, we obtain that 2𝐸^𝑁,𝑑𝐸𝑑𝑡^𝑁_𝛽,𝑁=𝐴^𝑁_𝛽,1+𝐴^𝑁_𝛽,2, (3.9) where 𝐴^𝑁_𝛽,1=−2𝐺^𝑁_𝛽,1,𝐸^𝑁_𝛽,𝑁,𝐴^𝑁_𝛽,2=2𝐺^𝑁_𝛽,2,𝐸^𝑁_𝛽,𝑁. (3.10) Using (2.20) and the Cauchy inequality, we deduce that 2𝐸^𝑁,𝑑𝐸𝑑𝑡^𝑁_𝛽,𝑁=2𝐸^𝑁,𝑑𝐸𝑑𝑡^𝑁=𝐸^𝑁(∞)²,𝐴^𝑁_𝛽,1≤2𝐺^𝑁_{𝛽,1𝛽,𝑁}𝐸^𝑁_𝛽,𝑁=2𝐺^𝑁_𝛽,1𝐸^𝑁. (3.11) Thus (3.9) reads 𝐸^𝑁(∞)²≤𝐴^𝑁_𝛽,2+2𝐺^𝑁_𝛽,1𝐸^𝑁. (3.12) Since there exists 𝛾₀>0 such that 𝑓𝑢₁,𝑣−𝑓𝑢₂,𝑣𝑢₁−𝑢₂≤−𝛾₀𝑢₁−𝑢₂²,∀𝑢₁,𝑢₂𝐴∈ℝ,_𝑁𝛽,2≤−2𝛾₀𝐸^𝑁2, (3.13) where 𝛾₀=𝛾₀𝜏. Therefore 𝐸^𝑁(∞)²+𝛾₀𝐸^𝑁2≤1𝛾₀𝐺^𝑁_𝛽,1². (3.14) Hence it suffices to estimate ‖𝐺^𝑁_𝛽,1‖². With the aid of (2.26), Lemmas 3.1 and 3.2, we deduce that, for 𝑟≥1, 𝐼̃_𝛽,𝑁𝑣−𝑣=𝐼_𝛽,𝑁𝑒^{(1/2)𝛽𝑡}𝑣−𝑒^{(1/2)𝛽𝑡}𝑣_𝜔𝛽≤𝑐𝛽𝑁^1/2−𝑟/2ℛ⁽¹⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑣. (3.15) On the other hand, 𝑑̃𝐼𝑑𝑡^𝛽,𝑁1𝑣−𝑣=−2𝛽𝑒^{−(1/2)𝛽𝑡}𝐼_𝛽,𝑁𝑒^{(1/2)𝛽𝑡}𝑣−𝑒^{(1/2)𝛽𝑡}𝑣+𝑒^{−(1/2)𝛽𝑡}𝑑𝐼𝑑𝑡^𝛽,𝑁𝑒^{(1/2)𝛽𝑡}𝑣−𝑒^{(1/2)𝛽𝑡}𝑣. (3.16) It follows from the above result that for 𝑟≥1𝑑̃𝐼𝑑𝑡^𝛽,𝑁𝑣−𝑣≤𝑐𝛽𝑁^1/2−𝑟/2𝛽ℛ⁽¹⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑣+ℛ⁽²⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑣. (3.17) Consequently, 𝐺^𝑁_𝛽,1≤𝑑̃𝐼𝑑𝑡^𝛽,𝑁+𝑑𝑈−𝑈𝐼̃𝑑𝑡𝑈−_𝛽,𝑁𝑑𝑈𝑑𝑡≤𝑐𝛽𝑁^1/2−𝑟/2𝛽ℛ⁽¹⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑈+ℛ⁽²⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑈+ℛ⁽¹⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑑𝑈.𝑑𝑡 (3.18) Thus, (3.14) implies that 𝐸^𝑁(∞)+𝛾₀^1/2𝐸^𝑁≤𝑐𝛾₀^(1/2)𝛽𝑁^{(1/2)−(𝑟/2)}𝛽ℛ⁽¹⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑈+ℛ⁽²⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑈+ℛ⁽¹⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑑𝑈.𝑑𝑡 (3.19) Furthermore, using (3.15) again, we obtain that, for any 𝛽>0, 𝑈−𝑢^𝑁≤𝑐𝛾₀𝛽𝑁^{(1/2)−(𝑟/2)}×𝛾₀+𝛽ℛ⁽¹⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑈+ℛ⁽²⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑈+ℛ⁽¹⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑑𝑈,𝑈(∞)−𝑢𝑑𝑡_𝑁(∞)≤𝐼̃_𝛽,𝑁𝑈(∞)−𝑈(∞)+𝐸^𝑁(∞𝐼̃)≤2_𝛽,𝑁𝑈−𝑈^1/2𝐼̃_𝛽,𝑁𝑈−𝑈₁^1/2+𝐸^𝑁(∞)≤𝑐𝛽𝑁^1/2−𝑟/2𝛽𝛾₀^−1/2+𝛽+1ℛ⁽¹⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑈+1+𝛾₀^−1/2ℛ⁽²⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑈+𝛾₀^−1/2ℛ⁽¹⁾_{𝑁,𝑟,𝛽}𝑒^{(1/2)𝛽𝑡}𝑑𝑈.𝑑𝑡 (3.20) This completes the proof.

Remark 3.4. The norm ‖𝑈‖ is finite as long as that 𝑓(𝑢,𝑣) satisfies certain conditions. If 𝑓(𝑢,𝑣) satisfies conditions 𝑓𝑢₁,𝑣−𝑓𝑢₂,𝑣,𝑢₁−𝑢₂≤𝛾₁𝑢₁−𝑢₂²,∀𝑢₁,𝑢₂∈ℝ^𝑁,𝑓𝑢,𝑣₁−𝑓𝑢,𝑣₂≤𝛾₂𝑣₁−𝑣₂,∀𝑣₁,𝑣₂∈ℝ^𝑁,𝛾₁<0,0<𝛾₂<−𝛾₁, (3.21) and 𝑓(0,0)=0, then |𝑈(𝑡)|=𝑂(𝑒^−𝛾∗𝑡), 𝛾_∗>0, see Tian [5]. Furthermore, if 𝑓(𝑢,𝑣) fulfills some additional conditions, then the norms appearing at the right sides of (3.20) are finite. Therefore, for certain positive constant 𝑐_∗ depending only on 𝛽, 𝑈−𝑢^𝑁+𝑈(∞)−𝑢^𝑁(∞)=𝑐_∗11+𝛾₀(ln𝑁)^1/2𝑁^1/2−𝑟/2. (3.22) Consequently, for 𝑟>1, the scheme (2.40) has the global convergence and the spectral accuracy in 𝐿²(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 𝑁 increases, with the convergence rate as 𝑐_∗(ln𝑁)^1/2𝑁^1/2−𝑟/2. On the other hand, for any fixed 𝑁, the norm ‖𝑈−𝑢^𝑁‖ is bounded, and so the pointwise numerical error decays automatically as 𝑡→∞, at least less than 𝑐_𝑁𝑡^−1/2, 𝑐_𝑁 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 𝑈^→(𝑡)=𝑈⁽¹⁾(𝑡),𝑈⁽²⁾(𝑡),…,𝑈^(𝑚)(𝑡)^𝑇,_→𝑉(𝑡)=𝑉⁽¹⁾(𝑡),𝑉⁽²⁾(𝑡),…,𝑉^(𝑚)(𝑡)^𝑇,_→𝑓𝑈^→,𝑉^→=𝑓𝑈^(1)→,𝑉^→,𝑓𝑈^(2)→,𝑉^→,…,𝑓𝑈^(𝑚)→,𝑉^→𝑇. (3.23) We consider the system 𝑑𝑑𝑡^→𝑈(𝑡)=^→𝑓𝑈^→(𝑡),𝑉^→(𝑡),𝑡>0,_→𝑈(𝑡)=Φ^→(𝑡),−𝜏≤𝑡≤0. (3.24) We approximate 𝑈^→(𝑡) by ^→𝑢^𝑁(𝑡). We can derive a numerical algorithm which is similar to (2.38). Further, let |𝑉^→|_𝐸 be the Euclidean norm of 𝑉^→. Assume that 𝑓𝑍^→→₁,𝑉^→−^→𝑓𝑍^→₂,𝑉^→⋅𝑍^→₁−𝑍^→₂≤−𝛾𝑍₀^→₁−𝑍^→₂²_𝐸,𝛾₀>0. (3.25) 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 bigger 𝑁 decrease faster. But in practical computation, it is not convenient to use them with very big 𝑁. On the other hand, the distance between the adjacent interpolation nodes 𝑡^𝑁_𝛽,𝑗 and 𝑡^𝑁_{𝛽,𝑗−1} increases fast as 𝑁 and 𝑗 increase, especially for the nodes which are located far from the origin point 𝑡=0. This is one of advantages of the Laguerre interpolation approximation, since we can use moderate mode 𝑁 to evaluate the unknown function at large 𝑡, 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.

Let 𝑁 be a moderate positive integer, 𝛽>0, and the set of nodes {𝑡^𝑁_0,𝛽,𝑗}^𝑁_𝑗=0={𝑡^𝑁_𝛽,𝑗}^𝑁_𝑗=0. We use (2.40) with the interpolation nodes {𝑡^𝑁_0,𝛽,𝑗}^𝑁_𝑗=0 to obtain the original numerical solution 𝑢^(0,𝑁)(𝑡)=𝑢^𝑁(𝑡),0≤𝑡≤𝑡^𝑁_0,𝛽,𝑁. (4.1) Then we take 𝑡^𝑁_1,𝛽,0=𝑡^𝑁_0,𝛽,𝑁 and consider the following delay differential equation: 𝑑𝑈𝑑𝑡⁽¹⁾(𝑡)=𝜏𝑓𝑈⁽¹⁾(𝑡),𝑈⁽¹⁾𝑡−𝑡^𝑁_𝛽,𝑁,𝑡>𝑡^𝑁_1,𝛽,0,𝑈⁽¹⁾(𝑡)=𝑢^(0,𝑁)(𝑡),𝑡^𝑁_1,𝛽,0−𝜏≤𝑡≤𝑡^𝑁_1,𝛽,0. (4.2) We get the refined numerical solution 𝑢^(1,𝑁)(𝑡) for 𝑡^𝑁_1,𝛽,0≤𝑡≤𝑡^𝑁_1,𝛽,𝑁. Repeating the above procedure, we obtain the refined numerical solution 𝑢^(𝑚,𝑁)(𝑡) for 𝑡^𝑁_𝑚,𝛽,0≤𝑡≤𝑡^𝑁_{𝑚,𝛽,𝑁}.

5. Numerical Results

We consider the following system of four homogeneous delay differential equations: 𝑦̇₁(𝑡)=𝑦₃(𝑡),𝑦̇₂(𝑡)=𝑦₄(𝑡),𝑦̇₃(𝑡)=−2𝑛𝑦₂(𝑡)+1+𝑛²(−1)^𝑛𝑦₁(),𝑦̇𝑡−𝜋₄(𝑡)=−2𝑛𝑦₁(𝑡)+1+𝑛²(−1)^𝑛𝑦₂(𝑡−𝜋). (5.1) The initial functions and solutions are given by 𝑦₁(𝑡)(𝑡)(),𝑦=sincos𝑛𝑡₂(𝑡)(𝑡)(),𝑦=cossin𝑛𝑡₃(𝑡)=𝑦̇₁(𝑡),𝑦₄(𝑡)=𝑦̇₂(𝑡),𝑡∈[−𝜋,∞). (5.2) For 𝑛=1, Figures 1 and 2 show plots of the exact solution 𝑦₁(𝑡) and 𝑦₃(𝑡) 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 𝑁=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.

Acknowledgments

This work is supported by the Program of Science and Technology of Huainan no. 2011A08005. The authors wish to thank the anonymous referees for carefully correcting the preliminary version of the paper. Special and warm thanks are addressed to Professor Ram N. Mohapatra for his valuable comments and suggestions for improving this paper.