# A Numerical Approach for Singularly Perturbed Nonlinear Delay Differential Equations Using a Trigonometric Spline

Authors: M. Lalu; K. Phaneendra
Journal: Computational and Mathematical Methods
 (2022)

---

## Abstract

In this paper, a computational procedure for solving singularly perturbed nonlinear delay differentiation equations (SPNDDEs) is proposed. Initially, the SPNDDE is reduced into a series of singularly perturbed linear delay differential equations (SPLDDEs) using the quasilinearization technique. A trigonometric spline approach is suggested to solve the sequence of SPLDDEs. Convergence of the method is addressed. The efficiency and applicability of the proposed method are demonstrated by the numerical examples.


---

## Body

## 1. Introduction

Consider a nonlinear singularly perturbed delay differential equation in the form(1)εθ′′=Fs,θ,θ′s−δon0,1,under the interval and boundary conditions
(2)θs=μson−δ≤s≤0,θ1=γ,where 0<ε≪1 is a perturbation parameter and δ is a delay parameter of oε. The solution θs of Equations (1) and (2) must be continuous on [0,1] and differentiable continuously on (0,1). Suppose F is smooth function and follows the conditions
(3)∂∂wFs,θ,w≥0,∂∂wFs,θ,w≤0.(1)
∂/∂θ−∂/∂wFs,θ,w≥φ>0,φ is a positive constant(2)
The growth conditionFs,θ,w=Οw2asw⟶∞foralls∈0,1 and all real θ and w.Under the above conditions, Equations (1) and (2) have a unique solution [1] for δ=0.The proposed equation usually plays an important role in illustrating different applications, such as theory of nonpremixed combustion [2], geodynamics [3], oceanic and atmosphere circulation [4], and chemical reactions [5]. More attention has been given in the past to the computational analysis of SPDDEs [6–9]. However, motivation for the research and solution of the SPNDDE has been increasing in the last few years. These problems may have steep exponential boundary layers as a solution. Classical methods for solving such types of problems are ineffective since a boundary layer structure is present when the perturbation parameter goes to zero. For these equations, effective numerical methods should be established, the accuracy of which does not depend on ε. Hence, in this work, we proposed a higher order numerical scheme using a trigonometric spline which gives more accuracy with a smaller number of mesh points. The existence and originality of the solutions of a SPNDDE with shift were studied by Lange and Miura [10]. The authors in [11] presented a fixed-point strategy to solve a second order SPDDE. The authors in [12] assemble two methodical spectral Legendre’s derivative methods to solve numerically the Lane-Emden, Bratu’s, and singularly perturbed type equations. For generating numerical spectrum solutions to linear and nonlinear second-order boundary value problems, a new operational matrix approach based on shifted Legendre polynomials is introduced and studied in [13].In [14], the authors proposed schemes with finite differences for solving the system of SPNDDE. In [15], a B-spline collocation method is constructed to solve Equations (1) and (2). In [16], the authors used shifted Legendre polynomials for studying the spectral collocation approach to solve neutral functional-differential equations with proportional delays. In [17], the Legendre spectral collocation approach is suggested by the authors for handling multipantograph delay boundary value problems. In [18], a new numerical method is proposed for solving a class of delay time-fractional partial differential equations. The fractional partial differential equations are reduced into an associated system of algebraic equations that may be solved by some robust iterative solvers using the localization method, which is based on space-time collocation in some appropriate points. In [19], the authors developed a numerical technique for nonlinear singly perturbed two-point boundary value problems based on a noniterative integration method with a modest deviation argument.The following is a concise summary of the contents of the paper. In Section2, the approach of quasilinearization and the analysis of convergence are discussed. The continuous problem is discussed in Section 3. In Section 4, the procedure using a trigonometric spline for the solution of the problem is derived. Error estimates of the proposed scheme are discussed in Section 5. Numerical examples and computational results are shown in Section 6. Finally, the Section 7 ends with the conclusion.


## 2. The Method of Quasilinearization

Using the method of quasilinearization [20], the given nonlinear differential Equations (1) and (2) are reduced into a sequence of SPLDDEs. We take the initial approximation θ0s which serves as a starting point for the function θs in F and expand Fs,θs,θ′s−δ, around the function θ0s; we get
(4)Fs,θ1,θ′1s−δ=Fs,θ0,θ′0s−δ+θ1−θ0∂F∂θs,θ0,θ′0+θ′1−θ′0∂F∂θ′s,θ0,θ′0+⋯.In general, we can write(5)Fs,θν+1,θ′ν+1s−δ=Fs,θ0,θ′0s−δ+θν+1−θν∂F∂θs,θ0,θ′0+θ′ν+1−θ′ν∂F∂θ′s,θ0,θ′0+⋯,for ν=0,1,2,⋯Using the quasilinearization technique, Equations (1) and (2) become
(6)εθ′′ν+1s+−∂Fν∂θ′θ′ν+1s−δ+−∂Fν∂θθν+1s=Fν−−∂Fν∂θ−∂Fν∂θ′θ′νs−δ,ν=0,1,⋯,with
(7)θν+1s=μs,−δ≤s≤0,θν+11=γ,and Fν=Fs,θν,θ′νs−δ. Thus, Equation (6) with Equation (7) is linear in θν+1s. Now, we solve the problems given by Equations (6) and (7) using the nonpolynomial spline method.Theoretically, the solution to the nonlinear problem satisfies(8)limn⟶∞θνs=θ∗s,0≤s≤1,where θ∗s is the solution of the nonlinear problem. Computationally, we require
(9)θν+1s−θνs<Tol.,0≤s≤1.Here, Tol. is a prescribed small tolerance. Once the tolerance test is achieved, the iteration is terminated.


## 3. Convergence Analysis

The convergence of the sequence of solutionsθν is obtained as follows. For convenience purpose, we refer Fs,θ,θ′s−δ as Fθ in the entire convergence part. Consider the problem
(10)εθ′′=Fθ,with
(11)θs=φson−δ≤s≤0,θ1=γ.After quasilinearization, we have a sequenceθν of linear equations defined by the following recurrence relation:
(12)εθ′′ν+1s=Fθν+θν+1s−θνsF′θν,(13)θν+1s=φs,−δ≤s≤0,θν+11=γ,where F′θ=∂Fθ/∂θ.Letθ0s be an initial approximation; then using Equation (12), we have
(14)εθ′′νs=Fθν−1+θνs−θν−1sF′θν−1.Using Equations (12) and (14), we have
(15)εθν+1s−θνs′′=Fθν−Fθν−1−θνs−θν−1sF′θν−1+θν+1s−θνsF′θν.Equation (15) is a differential equation of second order in θν+1s−θνs. Thus, by using Green’s function, the integral form of Equation (15) is
(16)εθν+1s−θνs=∫01Gs,tFθν−Fθν−1−θνs−θν−1sF′θν−1+θν+1s−θνsF′θνds,where the Gs,t is the Green’s function and determined by [21]
(17)Gs,t=s−1t,0≤t≤s≤1,st−1,0≤s≤t≤1,where maxs,tGs,t=1/4. By using the mean value theorem, we have
(18)Fθν−Fθν−1=θνs−θν−1sF′θν−1+θνs−θν−1s22F′′t,where θν−1s≤t≤θνs. Substituting Equation (18) into Equation (16), we get
(19)εθν+1s−θνs=∫01Gs,tFθν−Fθν−12F′′θ2+θν+1s−θνsF′θνds.Letmaxw≤1F′θ=a1,maxw≤1F′′θ=a2. On both sides of Equation (18), taking the maximum of the moduli over the region of interest, we get
(20)max0≤s≤1θν+1s−θνs≤14ε∫01max0≤s≤1θνs−θν−1s2max0≤s≤1F′′θds+14ε∫01max0≤s≤1θνs−θν−1smax0≤s≤1F′θνds.A simplification yields(21)max0≤s≤1θν+1s−θνs≤K1max0≤s≤1θνs−θν−1s2,where K1=a2/8ε1−a1/4ε<1. This shows that, given K1<1, the sequence θνs of linear equations converges quadratically. As a result, to get the approximate solution of Equation (1) with Equation (2), it is required to estimate the solution of the sequence of SPLDDEs of the form
(22)εθ′′ν+1s+aνsθ′ν+1s−δ+bνsθν+1s=cνs,ν=0,1,⋯,with
(23)θν+1s=μs,−δ≤s≤0,θν+11=γ,where
(24)aνs=−∂Fν∂θ′,bνs=−∂Fν∂θ,cνs=Fν−θν∂Fν∂θ−∂Fν∂θ′θ′νs−δ.


## 4. Continuous Problem

When the delay argumentδ is oε, sequential expanding for the term θ′ν+1s−δ in Equation (22) yields
(25)Lπ=εθ′′ν+1s+pνsθ′ν+1s+qνsθν+1s=rνs,ν=0,1,⋯,where
(26)pνs=aνs1−aνsδ/ε,qns=bνs1−bνsδ/ε,rνs=Fνs−θν∂Fν/∂θ−θ′νs−δθ′′νs∂Fν/∂θ′1−aνsδ/ε,with
(27)θν+10=μ0=μ0,θν+11=γ.The boundary layer appears on the left or right side of the interval depending on the sign of the coefficientpνs, i.e., as pνs>0 or pνs<0, respectively.The operatorLπ=εd2/ds2+pνsd/ds+qνsI in Equation (25) satisfies Lemma 1.Lemma 1.
AssumeGs is a smooth function satisfying Gν+10≥0, Gν+11≥0. Then, LπGν+1s≤0,∀x∈0,1 implies Gν+1s≥0,∀s∈0,1,k=0,1,⋯Proof.
Letk∈0,1 be such that Gν+1k<0 and
(28)Gν+1k=mins∈0,1Gν+1s
Clearly,k∉0,1; therefore, G′ν+1k=0 and G′′ν+1k≥0.☐Now, consider(29)LπGν+1ms=k=εG′′ν+1k+pνkG′ν+1k+qνkGν+1k>0,which is opposite to our assumption.Hence,Gν+1s≥0,∀s∈0,1.Lemma 2.
Letθν+1s be the solution of Equation (25) with Equation (27); then, we have
(30)θν+1≤β∗−1rν+maxμ0,γ,∀ν=0,1,⋯,where . is the L∞ norm given by θν+1=max0≤s≤1θν+1s.Proof.
Letψ±ν+1s be two barrier functions specified by
(31)ψ±ν+1s=β∗−1rν+maxμ0,γ±θν+1s.☐Then, this implies(32)ψ±ν+1s=β∗−1rν+maxμ0,γ±θν+10=β∗−1rν+maxμ0,γ±μ0≥0,ψ±ν+1s=β∗−1rν+max∣μ0∣,∣γ∣±θν+11=β∗−1rν+max∣μ0∣,∣γ∣±γ≥0,Lπψν+1s=εψ±ν+1s′′s+pνsψ±ν+1s′+qνsψ±ν+1s=qνsβ∗−1rν+max∣μ0∣,∣γ∣±Lπθν+1s=qνsβ∗−1rν+max∣μ0∣,∣γ∣±rνs.Asqνs≤−β∗<0 implies qνs≤β∗−1≤−1 and since rν≥rνs, we have Lπϕν+1s≤0,∀s∈0,1. Using Lemma 1, the required estimate is obtained.Lemma 3.
Ifθν+1s satisfies Equations (25)–(27), then
(33)θν+1s=uν+1s+vν+1s+zν+1s,ν=0,1,2,⋯,where
(34)uν+1s=−εθ′ν+10pν0exp−pν0sε,vν+1s=−εθ′ν+11pν1exp−pν11−sε,zν+1s≤M1+ε−s+1exp−c11−sε,s=0,1,2,3,4,5,6,where c1 is constant and M is positive constant independent of h and ε.


## 5. Trigonometric Spline

The integration domain [0,1] with mesh size h=1/N is decomposed into N equal subintervals, so that si=ih,i=0,1,2,⋯,N are the nodes with 0=s0,1=sN. Let θs be the exact solution and θi be an approximation to θsi by the trigonometric spline Sis passing through the points si,θi and si+1,θi+1. Here, Sis meets the interpolatory conditions at si and si+1; also the continuity of first derivatives at the common nodes si,θi is fulfilled. For each ith subinterval, the trigonometric spline function Sis has the form
(35)Sis=ai+bis−si+cisinτs−si+dicosτs−si,i=0,1,⋯,N−1,where ai,bi,ci, and di are constants and τ is a free parameter.To get the values of the coefficients of Equation (35) in terms of θi,θi+1,Mi,andMi+1, the following are defined:
(36)Ssi=θi,Ssi+1=θi+1,s′′si=Mi,s′′si+1=Mi+1.Using simple calculations, the following expressions are obtained for the coefficients:(37)ai=θi+Miτ2,bi=θi−θi+1h+Mi+1−Miτθ,ci=Micosθ−Mi+1τ2sinθ,di=−Miτ2,where θ=τh, for i=0,1,2,⋯,N−1. Using the first derivative continuity at si,θi, that is Si+1′si=Si′si, we get the following relation:
(38)αMi+1+2βMi+αMi−1=θi−1−2θi+θi+1h2,fori=1,2,3,⋯,N−1.Here,α=−1/θ2+1/θsinθ,β=1/θ2−cosθ/θsinθ, Mj=θ′′sj,j=i−1,i.i+1 and θ=τh.


## 6. Method of Solution

At the grid pointssi, Equation (25) may be discretised by
(39)εMjν+1=rνsj−pνsjθ′ν+1sj−qνsjθν+1sj,∀j=i,i±1.Using Equation (39) in Equation (38) and utilising the first derivatives of θ using the following estimations:
(40)θi+1′ν+1=θi−1−4θi+3θi+12h,θi−1′ν+1=−3θi+1+4θi−θi−12h,θi′ν+1=1+2ωh2qi+1+ωh3pi+1+pi−12hθi+1−2ωpi+1+pi−1θi−1+2ωh2qi−1−ωhpi+1+3pi−12hθi−1+ωhri+1−ri−1,we get
(41)ε+3αpi−1h2−βpiwh2pi+1+3pi−1+2βpiωh3qi−1−αhpi+12−h2αqi−1+hβpiθi−1+−2ε−2hαpi−1+4βpi2ωpi+1+pi−1+2hαpi+1−2h2βqiθi+ε+αhpi−12−βpih2ω3pi+1+pi−1−2βpih3ωqi+1−3hαpi+12−h2αqi+1−hβpiθi+1=h2α−2βhωri−1+2βri+α+2βhωri+1.Using Equation (41), we have the following tridiagonal system:
(42)Ei−1θi−1+Fiθi+Gi+1θi+1=Hi,fori=1,2,3,⋯,N−1.Here,(43)Ei−1=ε+3αpi−1h2−βpiwh2pi+1+3pi−1+2βpiωh3qi−1−αhpi+12−h2αqi−1+hβpi,Fi=−2ε−2hαpi−1+4βpi2ωpi+1+pi−1+2hαpi+1−2h2βqiGi+1=ε+αhpi−12−βpih2ω3pi+1+pi−1−2βpih3ωqi+1−3hαpi+12−h2αqi+1−hβpiHi=h2α−2βhωri−1+2βri+α+2βhωri+1.


## 7. Error Estimate

The truncation error in the proposed numerical scheme is given by(44)μih=1−2α+βh2εθi″+4ωε+13β−2α3Piθi‴+−1−12αεθi412h4+oh6.Thus, for different values ofω, αandβ in the approach (Equation (42)), the following different orders are indicated:
(1)
The scheme of Equation (44) yields the second order approach for any arbitrary choice of α and β with α+β=1/2 and for any value of ω(2)
Forα=1/12,β=5/12,ω=−1/20ε, the scheme Equation (44) gives the sixth orderTheorem 4.
LetWjν+1,j=0,1,2,⋯.N be an approximate solution to θν+1s of Equation (1) obtained by the proposed method. Then, the required estimate is
(45)maxjθjν+1−Wjν+1≤Mh41+e−k1sj/ε+e−k21−sj/ε,for α=1/12,β=5/12,ω=−1/20ε. Here, k1,k2,M are positive constants, independent of h and ε.Proof.
Using Lemma3, we have
(46)Uj6ν+1=−pν0ε5exp−pν0sjε.
Therefore(47)μjuν+1≤Mh6ε5exp−pν0sjε.
Similarly,(48)μjϑν+1≤Mh6ε5exp−pν11−sjε,zj6ν+1≤M1+ε−5exp−c11−sjε,μizν+1≤Mh61+ε−5exp−c11−sjε.
Now,(49)μiθν+1=μiuk+1+μiϑν+1+μizν+1,μiθν+1≤Mh6ε51+exp−pν0sjε+exp−pν11−sjε.☐The matrix form of the system Equation (42) is
(50)Aθν+1=B,l=0,1,⋯,where A is the matrix of the system Equation (42), θν+1 and B are the corresponding vectors, and μiθν+1 is the local truncation error. Thus,
(51)maxjθjν+1−Wjν+1≤A−1maxjμiθν+1,ν=0,1,⋯By using a result in [22], we have
(52)A−1≤max1≤i≤N−1Eiν+Fiν+Giν,ν=0,1,⋯,giving A−1≤1/ch2.Therefore,(53)maxjθjν+1−Wjν+1≤Mh41+e−pν0sj/ε+e−pν11−sj/ε,for α=1/12,β=5/12,ω=−1/20ε. Here, k1,k2,M are positive constants, independent of h and ε.


## 8. Numerical Examples

To show the relevance and validity of the approach, it was implemented for the following problems. The maximum pointwise errors (MAEs) (EN,εK) are determined by using the double mesh principle [3]:
(54)EN,εK=max0≤s≤1θliN−θl2i2N.Example 1.
εθ′′s+θsθ′s−δ−θs=0 with the interval and boundary constraints
(55)θs=1,−δ≤s≤0,θ1=1.Example 2.
εθ′′s+2θ′s−δ−eθs=0 with the interval and boundary constraints
(56)θs=0,−δ≤s≤0,θ1=0.


## 9. Conclusion

To solve a singularly perturbed nonlinear delay differentiation equation, a computational technique is proposed using a trigonometric spline. The SPNDDE is reduced into a series of linear SPDDEs using quasilinearization. A trigonometric spline approach is suggested to solve the sequence of linear SPDDEs. The scheme was implemented on two problems. The values of the maximum absolute errors produced by the suggested scheme are compared to the results in [15, 23] presented in Tables 1–4. Comparisons reveal that the suggested scheme outperforms the methods given in [15, 23] in terms of maximum error. Results of simulation have shown that as we increase the value of the parameter N, the accuracy of the computed approximate solutions is significantly improved.Table 1
MAEs in Example1 with δ=0.8ε.
ε/N6412825651210242048Proposed method10−11.1651e−033.2884e−048.7186e−052.2435e−055.6893e−061.4325e−0610−23.3515e−031.0682e−032.9959e−047.9176e−052.0341e−055.1542e−0610−33.8819e−031.2730e−033.6165e−049.6155e−052.4775e−056.2866e−0610−43.9415e−031.2967e−033.6893e−049.8159e−052.5299e−056.4208e−0610−53.9476e−031.2991e−033.6967e−049.8363e−052.5352e−056.4344e−0610−63.9482e−031.2994e−033.6974e−049.8383e−052.5358e−056.4358e−0610−73.9482e−031.2994e−033.6975e−049.8385e−052.5358e−056.4359e−0610−83.9482e−031.2994e−033.6975e−049.8385e−052.5358e−056.4359e−06Results in [15]10−16.036e−021.973E−025.134e−031.276e−033.185e−047.963e−0510−25.894e−021.902e−025.945e−031.812e−035.279e−041.444e−0410−35.891e−021.895e−025.906e−031.791e−035.164e−041.382e−0410−45.893e−021.902e−025.921e−031.789e−035.153e−041.376e−0410−55.892e−021.901e−025.943e−031.806e−035.172e−041.376e−0410−65.899e−021.901e−025.938e−031.806e−035.250e−041.409e−0410−75.973e−021.902e−025.938e−031.805e−035.237e−041.420e−0410−86.713e−021.915e−025.940e−031.805e−035.234e−041.416e−04Table 2
MAEs in Example1 withδ=0.5ε.
ε/N3264128256512Proposed method2−13.6533e−049.5648e−052.4452e−056.1804e−061.5535e−062−25.6399e−041.4855e−043.8081e−059.6382e−062.4243e−062−49.0467e−042.4051e−046.1920e−051.5704e−053.9539e−062−61.0520e−032.8073e−047.2400e−051.8377e−054.6288e−062−81.0954e−032.9264e−047.5510e−051.9171e−054.8293e−062−121.1097e−032.9655e−047.6529e−051.9431e−054.8950e−062−161.1106e−032.9679e−047.6594e−051.9447e−054.8992e−062−321.1106e−032.9681e−047.6598e−051.9449e−054.8995e−06Result in [23]2−14.84e−41.31e−43.39e−58.66e−62.18e−62−21.54e−34.28e−41.13e−42.90e−57.37e−62−48.86e−32.82e−38.17e−42.21e−45.78e−52−63.87e−21.23e−32.91e−36.93e−41.68e−42−83.68e−28.95e−32.12e−35.22e−41.30e−42−124.05e−28.33e−31.87e−34.46e−41.11e−42−164.84e−28.33e−31.87e−34.45e−41.09e−42−324.84e−28.33e−31.87e−34.45e−41.09e−4Table 3
MAEs in Example2 with δ=0.4ε.
ε/N6412825651210242048Proposed method10−21.9998e−045.2245e−051.3345e−3.3720e−068.4746e−072.1242e−0710−32.0404e−045.3327e−051.3624e−3.4429e−068.6532e−072.1690e−0710−42.0446e−045.3438e−051.3653e−3.4501e−068.6714e−072.1736e−0710−52.0450e−045.3449e−051.3656e−3.4508e−068.6732e−072.1741e−0710−62.0450e−045.3450e−051.3656e−3.4509e−068.6734e−072.1741e−0710−72.0450e−045.3450e−051.3656e−3.4509e−068.6734e−072.1741e−0710−82.0450e−045.3450e−051.3656e−3.4509e−068.6734e−072.1741e−07Results in [15]10−25.623e−021.717e−025.387e−031.630e−034.750e−041.300e−0410−35.641e−021.718e−025.367e−031.620e−034.692e−041.269e−0410−45.640e−021.727e−025.418e−031.623e−034.687e−041.266e−0410−55.638e−021.725e−025.426e−031.650e−034.776e−041.267e−0410−65.638e−021.725e−025.421e−031.646e−034.836e−041.338e−0410−75.640e−021.725e−025.420e−031.645e−034.820e−041.335e−0410−85.655e−021.725e−025.420e−031.645e−034.817e−041.331e−04Table 4
MAEs in Example2 with δ=0.5ε.
ε/N3264128256512Proposed method2−12.1589e−069.4190e−072.3548e−075.8869e−081.4717e−082−22.2793e−054.9244e−061.2344e−063.0893e−075.0693e−082−42.1081e−051.2064e−063.3601e−062.9526e−075.3741e−082−82.7909e−057.0930e−061.7870e−063.0396e−077.5799e−082−122.7908e−057.0926e−061.7869e−063.0461e−077.5949e−082−162.7908e−057.0926e−061.7869e−063.0461e−077.5949e−082−182.7908e−057.0926e−061.7869e−063.0461e−077.5949e−08Result in [23]2−17.32e−61.71e−64.07e−79.85e−82.42e−82−22.88e−56.10e−61.37e−63.21e−77.75e−82−44.98e−59.33e−61.92e−64.26e−79.96e−82−82.34e−44.54e−56.89e−69.22e−71.20e−72−122.41e−45.85e−51.43e−53.54e−68.48e−72−162.41e−45.85e−51.43e−53.54e−68.80e−72−182.41e−45.85e−51.43e−53.54e−68.80e−7In addition, while the error values generally increase as the perturbation parameterε decreases, they are usually within reasonable limits even for small values of it. It is also worth noting that the approach works well even when h≥ε is used. Figures 1–4 depict the layer behaviour at various δ values. It has been noticed that when the delay value increases, the thickness of the boundary layer increases as well. The simulation results show that the computational method proposed in this study is capable of giving accurate results for SPNDDE.Figure 1
Layer profile in Example1 withε=10−1.Figure 2
Layer profile in Example1 with ε=10−2.Figure 3
Layer profile in Example2 withε=10−1.Figure 4
Layer profile in Example2 with ε=10−2.



---

*Source: 8338661-2022-02-13.xml*