# An Analysis of the Theta-Method for Pantograph-Type Delay Differential Equations

Authors: Fathalla A. Rihan; Ahmed F. Rihan
Journal: Complexity
 (2022)

---

## Abstract

The pantograph equation arises in electrodynamics as a delay differential equation (DDE). In this article, we provide theϑ-method for numerical solutions of pantograph equations. We investigate the stability conditions for the numerical schemes. The theoretical results are verified by numerical simulations. The theoretical results and numerical simulations show that implicit or partially implicit ϑ-methods, with ϑ>1/2, are effective in resolving stiff pantograph problems.


---

## Body

## 1. Introduction

Pantograph equations are a special kind of functional differential equations with proportional delays [1]. The name pantograph was derived in 1971 from the work of Ockendon and Tayler [2]. Over the past few years, pantograph equations have gained increasing importance in the investigation of various scientific models. These arise in industrial applications and in a variety of fields of pure and applied mathematics, such as electrodynamics, control systems, number theory, probability, and quantum mechanics. Since most of these types of differential equations cannot be solved analytically, researchers have developed numerical methods to solve them. There are a lot of research articles in the solution methods of pantograph equations, such as the exponential approximation [3], Taylor method [4], collocation method using Hermite polynomials [5], and improved Morgan-Voyce collocation method [6]. Legendre wavelet solutions are discussed in [7], and homotopy perturbation method is examined in Ref. [8]. The variational iteration method is discussed in Ref. [9], and the homotopy analysis method is discussed in Ref. [10]. The authors in Ref. [11] discuss the modified Chebyshev collocation method for pantograph equations. A pseudospectral method based on the Legendre principle is examined in Ref. [12], and a modified procedure based on the residual power series method is explored in [13]. Other spectral methods and recent research about numerical treatments of pantograph equations exist in Refs. [14–17].Despite the fact that there are several polynomial approximation methods, such as spectral and pseudospectral methods, a comprehensive study of step methods for solving pantograph delay differential equations is lacking. In this article, we present some numerical schemes using theϑ-method to solve pantograph equations numerically and stability conditions for such schemes. Although the maximum order of ϑ-methods is 2, the numerical schemes are very suitable and effective for pantograph equations in a long-run time interval [18–20].The general delay differential equations (DDEs) take the following form:(1)y′t=ℱt,yt,yγt,t≥t0,yt=ψt,t≤t0,where ℱ:t0,+∞×ℂn×ℂn⟶ℂn and yt:t0,+∞×ℝ+⟶ℂn. Time delay can be a constant, time-dependent, or state-dependent. A classic case that has been discussed in many papers is when γt=t−τ, and τ is positive. The pantograph equations are functional differential equations whose delays are proportional, where 0≤γt≤t,γ0=0. Differential equations with variable or proportional delays are more difficult to solve numerically than equations with constant delays.The purpose of this article is to consider a nonlinear pantograph type of DDE:(2)y′t=ℱt,yt,yλt,t≥t0,λ∈0,1,yt=ψt,t∈λt0,t0.The lag functionγt=λt, and ψ∈λt0,t0 is a continuous initial function. We remark that this kind of equation plays an interesting role in modeling many phenomena including biological and nonlinear dynamical systems. Considering the numerical aspect, it is important to identify how the numerical methods preserve the qualitative behavior of the analytic solutions. Herein, we investigate the stability properties of ϑ-methods when applied to pantograph delay differential equations.This study is organized as follows. Analytical stability conditions are discussed in Section2. In Section 3, we provide continuous ϑ-methods for the pantograph equation. Numerical stability analysis and conditions of the asymptotic stability are discussed in Section 4. Some numerical simulations showing the effectiveness of the theoretical results are provided in Section 5, and concluding remarks are given in Section 6.


## 2. Nonlinear Stability

One of the important characteristics of differential equation (2) is the sensitivity of a particular solution to small changes in the “parameters,” which occur in the equation, or in the “initial conditions,” that gives rise to a definition of stability. If we consider another system, defined by the same function ℱt,ut,uλt of (2) but with another initial condition, then(3)u′t=ℱt,ut,uλt,t≥t0,ut=ϕt,t∈λt0,t0.y0≠u0. Assume that u,v is the inner product of vectors u,v∈ℂN, such that u=u,u1/2.Definition 1.
The solution of (1) is asymptotically stable, with respect to perturbing the initial function, if yt−ut⟶0 for t⟶∞. The solution ξ-exponentially stable (or algebraically decay [19]) if yt−ut≤Ke−ξt−t0, where ξ>0 (see [1], Chapter 3).Theorem 1 (Halanay inequality [21]).
Suppose thatα>β>0 and pt be a continuous nonnegative function on λt0,t0 satisfying the inequality:(4)p′t≤−αpt+βsupλt≤s≤tps,t≥t0.
Then,(5)pt≤Ge−ξt−t0,t≥t0,where G=supλt0≤t≤t0pt (depends on the initial condition), and ξ is the unique positive solution of(6)ℋξ≔ξ−α+βeξ1−λ.Theorem 2.
The solution of (1) is ξ-exponentially stable (algebraically decay) if(7)ρ1t+qρ2t≤0,0≤q<1,where(8)ρ1t≥supu,y,y˜∈ℂny≠y˜Reℱt,y,u−ℱt,y˜,u,y−y˜y−y˜2,ρ2t≥supy,u,u˜∈ℂnu≠u˜ℱt,y,u−ℱt,y,u˜u−u˜.Proof.
Assume for the inner product.,. on ℂn such that (6) holds. For every t≥t0, we have (see [22])(9)12ddtyt−ut2=Rey′t−u′t,yt−ut=Reℱt,yt,yλt−ℱt,ut,uλt,yt−ut=Reℱt,yt,yλt−ℱt,yt,uλt,yt−ut+Reℱt,yt,uλt−ℱt,ut,uλt,yt−ut.
It follows from the definitions ofρ1t and ρ2t, and from the Schwartz inequality that(10)12ddtyt−ut2≤ρ1tyt−ut2+ℱt,yt,yλt−ℱt,yt,uλtyt−ut≤ρ1tyt−ut2+ρ2tyλt−uλtyt−ut.
Define(11)Yt:=yt−ut.
Note thatYt>0 for every t>t0 because we assume that the function f is such that (2) has unique solution y(t) for every initial condition yt0=y0. Then,(12)12ddtyt−ut2=yt−utddtyt−ut=YtY′t,so we have(13)YtY′t≤ρ1tY2t+ρ2tYλtYt,and hence,(14)Y′t≤ρ1tYt+ρ2tsupλt≤s≤tYs.
Applying Theorem1 to the above inequality yields the following:(15)Yt≤Ge−ξt−t0,t≥t0.
Here,G depends only on the initial conditions ψt−ϕt.Corollary 1.
For the scalar linear case, letℱt,yt,yλt=αyt+βyλt. Then, the solution of y′t=αyt+βyλt, with ρ1t=Reα and ρ2t=β, is algebraically asymptotically stable (see [23]).


## 3. Continuousϑ-Methods for Pantograph Equations

Here, we study a discretization of pantograph equation (2). Let 0<λ<1 be a positive real and T≡t0,t1,…,tn,… be the assigned mesh points, with the property that tn⟶∞ as n⟶∞. Let hn=tn+1−tn, where n=0,1,…. be the nonconstant step-size. We divide the whole interval if the solution 0,T into bounded subintervals(16)Ik=λ−kr,λ−k−1r,k≥0,ris a positive integer.We then partition any of them into a fixed numberm of subintervals of the same size using the grid points, such that(17)0=t0<tl<t2<⋯<tp=r<tp+1<⋯<tp+m=λ−1r,tp+k+1m=λ−k−1r,k=0,1,…Lety˜n denote an approximate to the exact solution yt of (2) at tn. Then, applying the linear ϑ-method to (2) yields(18)y˜n+1=y˜n+hnϑℱtn+1,y˜n+1,y^λtn+1+1−ϑℱtn,y˜n,y^λtn,where y^t extends the values of the numerical solution to nonmesh points. Since, for λ∈0,1, it is usually the case that λtn∉T, we shall then require a densely defined approximation via linear interpolation:(19)y^t=tn+1−thny˜n+t−tnhny˜n+1,tn≤t<tn+1,n=0,1.….To find the solution at nonmesh points, assuming a uniform constant step-sizehn=h, we have dense outputs for yλt:(20)y^λtn=λn+1−λnyλn+λn−λnλnyλn+1n=0,1.…,where . is denoting the integer part, so that tλn<λtn<tλn+1. Furthermore, we define the global mesh by partitioning every primary interval into a fixed number m of subintervals of the same size, so that the mesh points are defined by the following recursion formula:(21)hn+1=1−λmλn/m+1,n=0,1,…,with gridtn=λ−1tn−m,n>m.Therefore, whenϑ=1, it corresponds to the so-called backward (or implicit) Euler method. However, when ϑ=0, it corresponds to the forward method. Both of them are of 1− stage RK methods. Normally, the order of ϑ methods is 1, but for ϑ=1/2, which corresponds to the trapezium rule method, the order is 2.Consider the following scalar linear pantograph equation:(22)y′t=αyt+βyλt,λ∈0,1,t>t0,yt=ψt,λt0≤t≤t0.Applying theϑ− method for (22), with a constant step-size h=tn+1−tn, yields(23)y˜n+1−y˜n=αhϑy˜n+1+1−ϑy˜n+βhϑy˜n+1−m+1−ϑy˜n−m.This generates the sequencey˜nn≥1, when given y˜−ℓ=ψλℓh−t0,forℓ∈0,1,…,N.We also should mention that, in both analytical and numerical terms, delay differential equations with infinite lags and those with finite lags exhibit remarkable differences [24]. Let us compare (22) with the DDE:(24)y′t=αyt+βyt−τ,τ∈0,∞.Solution to (22) is an analytic function on 0,∞, while solution to (24) is initially nonsmooth and becomes smoother as t increases. The solution of (22) decays algebraically, whereas the solution of (24) decays exponentially (see [18, 25, 26]).


## 4. Numerical Stability

In this section, we can use the discrete Halanay inequality [21] to examine the stability conditions of the discrete schemes for linear pantograph (22).Theorem 3.
Givenα+β<0 and ξh is the solution of the characteristic equations (see [21]):(25)ξhm+1−1−αhξhm−βh=0,where m∈N is a positive integer number, and 0<αh<1. Assume that p˜n−m∞ is a sequence, satisfying numbers satisfies the recurrence relation:(26)z˜n+1−z˜nh≤−αhz˜+βhmaxℓ∈Tz˜n−ℓ,forn∈ℕ,where 0≤βh<αh,0<αhh<1. Then, z˜n≤Ghe−ξhtn−t0, where ξh>0 are the values that occur in (25) and Ghmaxℓ∈Tp˜ℓ.


### 4.1. Stability ofϑ−Method

From the auxiliary polynomial, it is possible to compute (asymptotic) stability regions of recurrence relations (23):(27)ξm+1−ϱ0ξm−ϱ1ξ−ϱ2=0,where(28)ϱ0=1+1−ϑαh1−ϑαh,ϱ1=ϑβh1−ϑαh,ϱ2=1−ϑβh1−ϑαh,1−ϑαh≠0,such that the roots of (27) must satisfy the condition ξi<1.A small perturbation in the initial conditionsy˜−ℓℓ∈T in (23) leads to a consequence variation δy˜ℓℓ≥1, satisfying(29)δy˜n+1=ϱ0δy˜n+ϱ1δy˜n+1−m+ϱ2δy˜n−m,n≥0.One can get [22].(30)δy˜n+12−δy˜n2=ϱ02−1δy˜n2+2ϱ0ϱ1δy˜n2δy˜n+1−m2+2ϱ0ϱ2δy˜n2δy˜n−m2+2ϱ1ϱ2δy˜n+1−m2δy˜n−m2+ϱ12δy˜n+1−m2+ϱ22δy˜n−m2.IfUV≠0, then SUV≤1/2U2+S2V2, with equally if S=U/Vv. Thus, UV=infS∈0,∞1/2SU2+S2V2≤U2/r+rV2 for all r∈0,∞. Therefore,(31)ϱiϱjδy˜lδy˜m≤ϱiϱj2δy˜l2rlm+rlmδy˜m2,for all arbitrary rlm∈0,∞.From (30) and (31), one can get(32)δy˜n+12−δy˜n2≤ϱ02−1δy˜n2+ϱ0ϱ1δy˜n2r1+r1δy˜n+1−m2+ϱ0ϱ2δy˜n2r2+r2δy˜n−m2+ϱ1ϱ2δy˜n−m+12r3+r3δy˜n−m2+ϱ12δy˜n+1−m2+ϱ22δy˜n−m2,for the arbitrary numbers r1,r2,r3∈0,∞. Therefore,(33)δy˜n+12−δy˜n2h≤−Ah♮δy˜n2+Bh♮maxℓ∈Tδy˜n−ℓ2,where(34)Ah♮≡Ah♮r1,r2=−1hϱ02−1+ϱ0ϱ1r1+ϱ0ϱ2r2.Forr3=1,(35)Bh♮≡Bh♮r1,r2=1hϱ0ϱ1r1+ϱ0ϱ2r2+ϱ1+ϱ22.We have following theorem and remarks.Theorem 4.
From inequality (33), using equation (26) of Theorem 3, if 0≤Bh♮<Ah♮ and hAh♮∈0,1 with any positive values of r1,r2, the recurrence relation (23) is ξ− exponentially stable.Remark 1.
Theϑ-method, as applied to equation (2) is asymptotically stable if and only ϑ>1/2.
It is important to mention that explicit methods are easy to implement and are recommended for solving nonstiff problems. For stiff problems, however, the standard explicit methods with poor stability properties suffer a lot from step-size reduction and turn out to be inefficient in terms of overall computational costs [27–29].Remark 2.
In order to solve stiff problems numerically, we should use implicit or partially implicitϑ-method’s schemes.


## 4.1. Stability ofϑ−Method

From the auxiliary polynomial, it is possible to compute (asymptotic) stability regions of recurrence relations (23):(27)ξm+1−ϱ0ξm−ϱ1ξ−ϱ2=0,where(28)ϱ0=1+1−ϑαh1−ϑαh,ϱ1=ϑβh1−ϑαh,ϱ2=1−ϑβh1−ϑαh,1−ϑαh≠0,such that the roots of (27) must satisfy the condition ξi<1.A small perturbation in the initial conditionsy˜−ℓℓ∈T in (23) leads to a consequence variation δy˜ℓℓ≥1, satisfying(29)δy˜n+1=ϱ0δy˜n+ϱ1δy˜n+1−m+ϱ2δy˜n−m,n≥0.One can get [22].(30)δy˜n+12−δy˜n2=ϱ02−1δy˜n2+2ϱ0ϱ1δy˜n2δy˜n+1−m2+2ϱ0ϱ2δy˜n2δy˜n−m2+2ϱ1ϱ2δy˜n+1−m2δy˜n−m2+ϱ12δy˜n+1−m2+ϱ22δy˜n−m2.IfUV≠0, then SUV≤1/2U2+S2V2, with equally if S=U/Vv. Thus, UV=infS∈0,∞1/2SU2+S2V2≤U2/r+rV2 for all r∈0,∞. Therefore,(31)ϱiϱjδy˜lδy˜m≤ϱiϱj2δy˜l2rlm+rlmδy˜m2,for all arbitrary rlm∈0,∞.From (30) and (31), one can get(32)δy˜n+12−δy˜n2≤ϱ02−1δy˜n2+ϱ0ϱ1δy˜n2r1+r1δy˜n+1−m2+ϱ0ϱ2δy˜n2r2+r2δy˜n−m2+ϱ1ϱ2δy˜n−m+12r3+r3δy˜n−m2+ϱ12δy˜n+1−m2+ϱ22δy˜n−m2,for the arbitrary numbers r1,r2,r3∈0,∞. Therefore,(33)δy˜n+12−δy˜n2h≤−Ah♮δy˜n2+Bh♮maxℓ∈Tδy˜n−ℓ2,where(34)Ah♮≡Ah♮r1,r2=−1hϱ02−1+ϱ0ϱ1r1+ϱ0ϱ2r2.Forr3=1,(35)Bh♮≡Bh♮r1,r2=1hϱ0ϱ1r1+ϱ0ϱ2r2+ϱ1+ϱ22.We have following theorem and remarks.Theorem 4.
From inequality (33), using equation (26) of Theorem 3, if 0≤Bh♮<Ah♮ and hAh♮∈0,1 with any positive values of r1,r2, the recurrence relation (23) is ξ− exponentially stable.Remark 1.
Theϑ-method, as applied to equation (2) is asymptotically stable if and only ϑ>1/2.
It is important to mention that explicit methods are easy to implement and are recommended for solving nonstiff problems. For stiff problems, however, the standard explicit methods with poor stability properties suffer a lot from step-size reduction and turn out to be inefficient in terms of overall computational costs [27–29].Remark 2.
In order to solve stiff problems numerically, we should use implicit or partially implicitϑ-method’s schemes.


## 5. Numerical Examples and Simulations

We present an example with some numerical simulations to illustrate the previous results.Example 1.
Consider the initial-value problem (IVP):(36)y′t=−40yt+20yλt,yt=1fort∈λt0,t0.
This IVP is considered a stiff problem. This problem can also be considered a singularly perturbed and stiff problemεy′t=αyt+βyλt, with ε=0.01,α=−0.4,β=0.2. Since α+β<0, where α=−40,β=20, the analytical solution is asymptotically stable. We apply the ϑ-method scheme to (36) with different values of ϑ=0.0,0.2,0.6. The mesh points are given by(37)0=t0<tl<t2<⋯<tp=r<tp+1<⋯<tp+m=λ−1r,tp+k+1m=λ−k−1r,k=0,1,…
Figures1 and 2 show the numerical simulations of model (36) with a full explicit scheme ϑ=0 and partially explicit with ϑ=0.2, for λ=0.75,0,8, 0.9. Although α+β<0, the numerical solutions when ϑ=0 are unstable. Usually, the asymptotic and stability of explicit numerical schemes are restricted by step-size, which is given in terms of λ. As the step-size decreases, the numerical approximation solutions become more asymptotic. Figure 3 shows the numerical simulations with fixed ϑ and different values of delays. Time delays significantly increase the stability of differential equations.
Figure4 displays numerical simulations when ϑ<0.5=0.2 with λ=0.75 and λ=0.8. The asymptotic and stability behavior are improved with bigger λ, and numerical solution is getting smoother with long-run time. However, Figures 5 and 6 show the numerical simulation with partially implicit scheme, with ϑ=0.6,0.9>0.5, with different values of λ=0.75,0.8,0.9. The solution is asymptotically stable regardless of the size of the step-size. The numerical schemes are unconditionally asymptotically stable with ϑ>1/2.Figure 1
Numerical simulations of model (36) when ϑ=0, and λ=0.75 (a), λ=0.8 (b).
(a)(b)Figure 2
Numerical simulations of model (36) with ϑ=0.0 (a) and ϑ=0.2 (b), when λ=0.9.
(a)(b)Figure 3
Numerical simulations of model (36) with ϑ=0.4 and different values of λ=0.7 (a) and λ=0.9 (b).
(a)(b)Figure 4
Numerical simulations of model (36) with ϑ=0.2 and λ=0.75 (a), λ=0.8 (b).
(a)(b)Figure 5
Numerical simulations of model (36) with ϑ=0.6, with λ=0.75 (a), λ=0.8 (b).
(a)(b)Figure 6
Numerical simulations of model (36) with ϑ=0.6 (a) and ϑ=0.9 (b), when λ=0.9.
(a)(b)


## 6. Conclusion

Time delays significantly increase the complexity and stability of differential equations. It is possible for a delay to stabilize an otherwise unstable system or unstabilize a stable one. We developed and discussed an efficient numerical scheme for solving pantograph-type delay differential equations using theϑ-method. Some stability conditions of analytical and numerical solutions of the problem have been investigated. From the theoretical results, numerical simulations can conclude that the ϑ-method, applied to pantograph equations, is a reliable method and asymptotically stable if and only if ϑ>1/2, regardless of the step-size of the scheme, especially with respect to the long-time behavior of the solutions. For ϑ<1/2, the stability of the numerical solution is almost unstable. Most of pantograph equations are stiff problems, and it is recommended to use an implicit or partially implicit sachems.We can extend the stability analysis, in the coming research, to other classes of higher-order Runge-Kutta methods with pantograph equations. The schemes can also be extended to stochastic pantograph delay differential equations [30, 31].



---

*Source: 8961352-2022-08-13.xml*