Int. J. Anal. Appl. (2023), 21:5 Computational Approach for a Singularly Perturbed Differential Equations With Mixed Shifts Using a Non-Polynomial Spline Kumar Ragula1, G. BSL Soujanya2,∗, D. Swarnakar3 1 Department Mathematics, Rajiv Gandhi University of Knowledge Technologies, Basar, India 2 Department Mathematics, University P.G College for Women,Kakatiya University, Warangal, India 3Department Mathematics, VNR Vignana Jyothi Institute of Engineering and Technology, Hyderabad, India ∗Corresponding author: gbslsoujanya@gmail.com Abstract. In this paper, a second order singularly perturbed differential difference equation with both the negative and positive shifts is considered. A fitted non-polynomial spline approach is applied to solve the problem. Taylor series expansion process is being used to produce an approximated form of the considered problem, and then a fitted non-polynomial spline approach is devised in the form of a three-term recurrence relation. The convergence of the method is examined, and a quadratic rate of convergence is achieved. The maximum absolute error with quadratic rate of convergence of the solution is recorded. Layer profile is examined using the graphs. 1. Introduction In applied science and engineering, there are many physical and biological problems with solutions that include boundary and interior layers for specific parameter values and these are known as singular perturbation problems(SPPs). A SPP is one in which no asymptotic expansion is always true over the interval as the perturbation parameter approaches zero. Due to the boundary or interior layer structure of the solutions,the numerical and analytical treatments of such SPPs are very difficult. There is a narrow transitional layer in which the solution changes most quickly, whereas the solution responds uniformly and slowly away from the layer. When the perturbation parameter approaches Received: Dec. 8, 2022. 2020 Mathematics Subject Classification. 65L11, 65L12. Key words and phrases. singularly perturbed differential-difference equation; boundary layer; non-polynomial spline; mixed shifts. https://doi.org/10.28924/2291-8639-21-2023-5 ISSN: 2291-8639 © 2023 the author(s). https://doi.org/10.28924/2291-8639-21-2023-5 2 Int. J. Anal. Appl. (2023), 21:5 zero, the regularity of the solution decreases. A singularly perturbed differential-difference equation (SPDDE) is formed by multiplying the highest order derivative term of a differential equation with a small positive perturbation parameter ε and involves a delay and advance terms. For the solutions such SPDDEs in the interval of boundary conditions is a challenging in the modeling of a variety of physical and biological problems like, the initial exit time problem in neuronal variability activation models [25], oscillations of the human pupil light reflex with delayed and mixed responses [15], red cell system evaluation [27] and bifurcated gap in a hybrid optical system [3] etc. For a detailed information on the SPDDEs is given in [2], [4], [5], [6], [17], [18]. A few numerical methods are proposed to obtain valid and realistic approaches for the solution for second order SPDDEs with mixed shifts (negative and positive) with boundary constraints. In the articles [11] - [14] the authors presented asymptotic analysis of SPDDEs. In [13], [14] Taylor series is used to analyze a term containing small shift. The same authors investigated the impact of small shifts on the oscillatory solution of the problem in [11]. In [1] the authors proposed an a nonstandard exponentially fitted finite difference method(FDM) to solve SPDDEs with boundary layer on each sides of the interval. To solve SPDDE with delay, an impactful Haar wavelet collocation methodology is developed in [20]. In the articles [21], [22] the authors suggested an exponentially fitted FDM to solve SPDDEs with delay and advanced terms and turning point problems. The authors advised a fourth order FDM with a fitting factor in [10], [16] to solve the considered problem. In [8], [9] the authors developed some numerical techniques to solve SPDDEs with mixed shifts. An asymptotic expansion estimation of the solution is designed to solve SPDDEs in [23]. The author suggested successive complementary expansion method (SCEM) for solving this model in [19]. A mixed FDM is proposed to solve SPDDEs with both the shifts delay and advanced in [24]. The approach given in this study gives a novel difference method to the SPDDE with mixed shifts. In Section 2, the problem description is given and the non-polynomial spline is discussed. In Section 3, a three term difference scheme is devised using the non-polynomial spline. Convergence analysis of this computational approach is discussed in Section 4. In Section 5, the computational solutions for the test problems to illustrate the technique, along with comparisons to alternative approach and conclusions are presented. 2. Definition of the problem and Non-Polynomial Spline Consider the following SPDDE with delay and advance terms εz′′ (θ)+p(θ)z′ (θ) + q (θ)z (θ−δ)+ r (θ)z (θ)+ s (θ)z (θ+η) = f (θ) , (2.1) for θ ∈ (0,1) with the boundary constraints z (θ)= ϕ (θ) , −δ ≤ θ ≤ 0; z (θ)= ψ(θ) ,1≤ θ ≤ 1+η. (2.2) Int. J. Anal. Appl. (2023), 21:5 3 where 0 < ε � 1 and p(θ), q (θ), r (θ), f (θ), ϕ(θ) and ψ(θ) are sufficiently smooth functions on (0,1) and 0 < δ, η = o(ε), δ, η are delay and advance shifts respectively. The solution of equation (2.1) with (2.2) represents the layer behavior at each end of the interval if p(θ)−δq (θ)+ηs (θ) > 0 and p(θ)−δq (θ)+ηs (θ) < 0 (left end and right end). Depending on the sign of q (θ)+ r (θ)+ s (θ), existence of boundary layers in two cases reported in [21]. The Taylor’s series expansions of z (θ−δ), z (θ+η) in the neighborhood of the point θ, we have z (θ−δ) ≈ z (θ)−δz′ (θ)+O ( δ2 ) and z (θ+η)≈ z (θ)−ηz′ (θ)+O ( η2 ) (2.3) Using (2.3) in (2.1), we get εz′′ (θ)+a(θ)z′ (θ)+b(θ)z (θ)= f (θ) , 0 < θ < 1 (2.4) with the boundary constraints z (0)= ϕ(0)= ϕ0, z (1)= ψ(1)= ψ1 (2.5) Consider the domain [0,1] and it is split into N equal length of sub-domains with constant step size h. Let 0 = θ0 < θ1 < ... < θN = 1 be the N grid points. Then we have, θi = ih and in each sub-domain [θi,θi+1] , i =1,2, . . . ,N −1 the non-polynomial spline is of the form Gi (θ)= aisinτ (θ−θi)+bi ( e−τ(θ−θi) +eτ(θ−θi) ) +ci (θ−θi)+di (2.6) where ai, bi, ci and di are constant coefficients, and τ 6= 0 arbitrary parameter. To derive the coefficients ai, bi, ci and di (2.6) in terms of zi, zi+1, Qi and Qi+1, we define Li (θi) = zi, Li (θi+1) = zi+1 (2.7) L′′i (θi) = Qi, L ′′ i (θi+1) = Qi+1 (2.8) By using the conditions (2.7) and (2.8), we calculate the coefficients in (2.6) as .   ai = Qih 2(eω+e−ω) 2ω2sinω − Qi+1h 2 2ω2sinω bi = Qih 2 2ω2 ci = zi+1− zi h + Qih 2(1 − eω−e−ω) hω2 + Qi+1h 2 hω2 di = yi − Qih 2 ω2 (2.9) where ω = τh. Using the first derivative of continuity, Li−1 (m) (θi)= Li (m) (θi) ,m =1, we obtain the relation (z i−1 −2zi +zi+1)= h 2(αQi−1 +βQi +γQi+1), i =1,2, . . . ,N −1 (2.10) where, α = ω ( e−ω +eω ) cosω +ω ( − e−ω +eω ) sinω +2 ( 1−eω −e−ω ) sinω 2ω2sinω , 4 Int. J. Anal. Appl. (2023), 21:5 β = −2ωcosω +2sinω−ω ( eω +e−ω ) −2 ( 1−eω −e−ω ) sinω 2ω2sinω , γ = ω − sinω ω2sinω If h → 0, then ω = hτ → 0, we have (α,β,γ)→ ( 1 6 , 4 6 , 1 6 ) , (2.10) reduces to cubic spline [7]. 3. Numerical Algorithm At the grid points θi, (2.4) can be written as εz′′i =−a(θi)z ′ i −b(θi)zi + f (θi) Using L′′i (θi) = Qi = z ′′ i in above equation, we get εQj =−aj (θi)z′i −bj (θi)zi + fj (θi) for j = i, i ±1 (3.1) Substitute (3.1) in (2.10) and then using z′j for j = i, i ±1 i.e z′i ≈ 1 2h (zi+1 −zi−1) , z′i+1 ≈ 1 2h (3zi+1 − 4zi +zi−1) and z′i−1 ≈ 1 2h (−zi+1 +4zi −3zi−1) ε h2 (z i+1 −2zi +zi−1) =−αai−1 (−zi+1 +4zi −3zi−1) 2h −βai (−zi−1 +zi+1) 2h −γai+1 (zi−1 −4zi +3zi+1) 2h −αbi−1zi−1 −βbizi −γbi+1zi+1+( αf i−1 +βfi +γf i+1 ) (3.2) We incorporate a fitting parameter σi (ρ) in the present approach to increase the solution’s accuracy and manage the layer behaviour. Then, we have εσi (ρ) h2 (z i+1 −2zi +zi−1) =−αai−1 (−zi+1 +4zi −3zi−1) 2h − βai (zi+1 −zi−1) 2h −γai+1 (zi−1 −4zi +3z i+1) 2h −αbi−1zi−1 −βbizi− γbi+1zi+1 + ( αf i−1 +βfi +γf i+1 ) (3.3) On simplification, we obtain the following tridiagonal system Eizi−1 +Fizi +Gizi+1 = Hi, i = 1,2, . . . , N −1 (3.4) where Ei = εσi − 3αhai−1 2 − βhai 2 + γhai+1 2 + h2αbi−1 , Fi = −2σiε + 2αhai−1 − 2γhai+1 + h2βbi , Int. J. Anal. Appl. (2023), 21:5 5 Gi = εσi − αhai−1 2 + βhai 2 + 3γhai+1 2 + h2γbi+1 , Hi = h 2 ( α f i−1 +βfi +γf i+1 ) The system (3.4) can be solved using Thomas algorithm with the constraints z (0)= ϕ0, z (1)= ψ1. To calculate the fitting parameter, we adopt the process given in [18]. The following is an approxi- mation for the solution of the homogeneous problem of (2.1) z (θ)= z0 (θ)+ a(0) a(θ) (α−z0 (0))e − ∫ θ 0 ( a(θ) ε −b(θ) a(θ) ) dθ +o (ε) (3.5) where z0 (θ) is the solution of a(θ)z0′ (θ)+b(θ)z0 (θ)= f (θ) , z0 (1)= ψ1 By using the expansion for a(θ) and b(θ) about the point zero, then (3.5) becomes z (θ)= z0 (θ)+(ϕ0 −z0 (0))e − ( a(θ) ε ) θ +o (ε) From (3.5), we have lim h→0 z (ih)= z0 (0)+(ϕ0 −z0 (0))e−a(θ)iρ, Using these limit values in (3.3), we obtain the following fitting factor σi (ρ)= ( α+ β 2 ) aiρ Coth (aiρ 2 ) ,where ρ = h ε . 4. Convergence Analysis The local error estimate for the numerical scheme of (3.4) is Ti (h)= h 2 [1− (α+β +γ)]εz′′i +h 3 (γ −α) [ b′izi +a ′ iz ′ i +biz ′ i − f ′ i ] +h4 (γ +α) 2 [ b′′i zi +a ′′ i z ′ i +2b ′ iz ′ i +2a ′ iz ′′ i +biz ′′ i − f ′′ i ] + (α−γ)O ( h5 ) +O ( h6 ) Hence, with (α,β,γ)= ( 1 6 , 4 6 , 1 6 ) , truncation error is of fourth order. With the help of (2.2), the matrix form of (3.4) is( R̃+ G̃ ) Z +M̃ +T (h)= O (4.1) where R̃ =   −2εσ εσ 0 0 ... 0 εσ −2εσ εσ .. ... 0 0 .. .. .. ... 0 .. .. .. .. .. .. .. .. .. .. .. .. 0 .. .. .. εσ −2εσ   , 6 Int. J. Anal. Appl. (2023), 21:5 G̃ = [xi,vi,wi] =   v1 w1 0 0 ........ 0 x2 v2 w2 0 ........ 0 0 x3 v3 w3 ........ 0 .. .. .. .. .. .. .. .. .. .. .. .. 0 .. .. 0 xN−1 vN−1   xi =− 3αhai−1 2 − βhai 2 + γhai+1 2 + h2αbi−1 , vi =2αhai−1 − 2γhai+1 + h2βbi , wi =− αhai−1 2 + βhai 2 + 3γhai+1 2 + γh2bi+1 , for all i =1,2, . . . ,N −1 M̃ = [m1 +(εσ +x1)φ(0) ,m2,m3 , . . . ,mN−2,mN−1 +(εσ +wN−1)γ] T where, mi = h 2 ( γf i+1 +βfi +αf i−1 ) , for i =1,2, . . . ,N −1, T (h)= o ( h4 ) , Z = [Z1,Z2, ...,ZN−1] T , T (h)= [t1,t2, ...,tN−1] T and O = [0,0, ...,0] T are corresponding vectors of (4.1). Let z = [z1,z2, ...,zN−1] T ∼= z which satisfies the equation( R̃+ G̃ ) z +M̃ = 0 (4.2) Let ei = zi−Zi, i =1,2 . . . ,N−1 denote the discretized error so that Ẽ = [e1,e2, ...,eN−1] T = z−Z. From (4.1) and (4.2), we obtain the error equation( R̃+ G̃ ) Ẽ = T(h) (4.3) Let |ai | ≤ K1, | bi| ≤ K2 so that, if Q̃i,j is the (i, j) th element of matrix G̃, then ∣∣Q̃i, i+1∣∣ = |wi| ≤ ε+h(α+β +3γ)K1 +h2αK2, i =1,2, . . . , N −2 (4.4a) ∣∣Q̃i, i−1∣∣ = |ui| ≤ ε+h(3α+β +γ)K1 +h2αK2, i =2,3, . . . ,N −1 (4.4b) As a result, for relatively small h (h → 0), we see that∣∣Q̃i,i+1∣∣ < ε, ∀ i = 1,2, . . . , N −1, ∣∣Q̃i,i−1∣∣ < ε,∀ i = 2,3, . . . ,N −1 (4.4c) Hence, ( R̃+ G̃ ) is irreducible [26]. Let the sum of elements of ith row of ( R̃+ G̃ ) be Si, then we have Si =−εσ + 3αhai−1 2 + βhai 2 − γhai+1 2 +h2 (γbi+1 +βbi) , for i =1 Si = h2 ( αbi−1 +βbi +γbi+1 ) , for i = 2,3, . . . , N −2 Si =−εσ + αhai−1 2 − βhai 2 − 3γhai+1 2 +h2 (αbi−1 +βbi) , for i = N −1 Int. J. Anal. Appl. (2023), 21:5 7 Let K1∗ = min 1 ≤i≤ N−1 | ai|, K1∗ = max 1 ≤i ≤N |ai| ,K2∗ = min 1 ≤i ≤N−1 | bi|, K2 ∗ = max 1 ≤i ≤N | bi| then 0≤ K1∗ ≤ K1 ≤ K1∗, 0≤ K2∗ ≤ K2 ≤ K2∗ For relatively small h, ( R̃+ G̃ ) is monotone [26]. Hence ( R̃+ G̃ )−1 exists and ( R̃+ G̃ )−1 ≥ 0. Thus from (4.3), we have ||Ẽ|| ≤ ||R̃+ G̃||−1||T || (4.5) For relatively small h, we have Let ( R̃+ G̃ )−1 i,k be the (i,k)th element of ( R̃+ G̃ )−1 and define ||R̃+ G̃||−1 = max 1≤ i≤ N−1 N−1∑ k=1 ( R̃+ G̃ )−1 i,k , and ||T(h)||= max 1≤ i≤ N−1 |Ti|. Since (R̃+ G̃)−1 i,k ≥ 0 and ∑N−1 k=1 (R̃+ G̃) −1 i,k .Sk =1, for i =1(1)N −1( R̃+ G̃ )−1 i, k ≤ 1 Si < 1 h2K2 , i =1. (4.6a) ( R̃+ G̃ )−1 i, k ≤ 1 Si < 1 h2K2 , i = N −1 (4.6b) Further more, N−1∑ k=1 ( R̃+ G̃ )−1 i,k ≤ 1 min 2≤i≤N−2 Si < 1 h2K2 , for i =2,3, . . . , N −2 (4.6c) With the help of (4.6a), (4.6b), (4.6c) and using (4.5), we have ||Ẽ|| ≤ O (h2) This illustrates the second-order convergence for the scheme (3.4) with (α,β,γ)= ( 1 6 , 4 6 , 1 6 ) 5. Numerical Examples To show the validity and robustness of the suggested technique, we reported the computational results of the four example problems in terms of maximum absolute errors (MAEs) with calculated rates of convergence (ROC) in the tables. Since the exact solutions are not known for considered examples, the MAEs are estimated using the double mesh approach using the formula EN = max 0≤i≤N ∣∣ziN −z2i2N∣∣ where ziN and z2i2N are the computational solutions of the example problem for N and 2N grid points respectively. Further, the rate of convergence is determine by the formula RN = log2 ∣∣∣∣ ENE2N ∣∣∣∣ . 8 Int. J. Anal. Appl. (2023), 21:5 Example 5.1. εz′′ (θ)+z′ (θ)−z (θ−δ)+z (θ)−z (θ+η)=−1, subject to boundary constraints z (θ)=   1, − δ ≤ θ ≤ 0 1, 1≤ θ ≤ 1+η. Example 5.2. εz′′ (θ)+2.5z′ (θ)−2exp(θ)z (θ−δ)−z (θ)−θz (θ+η)=1, subject to boundary constraints z (θ)=   1, − δ ≤ θ ≤ 0 1, 1≤ θ ≤ 1+η. Example 5.3. εz′′ (θ)− ( 1+exp ( −θ2 )) z ′ (θ)−θz (θ−δ)−θ2z (θ)− (1.5−exp(−θ))z (θ+η)=1 , with boundary constraints z (θ)=   1, − δ ≤ θ ≤ 0 1, 1≤ θ ≤ 1+η. Example 5.4. εz′′ (θ)− ( 1+exp ( θ2 )) z ′ (θ)−θz (θ−δ)+θ2z (θ)− (1−exp(−θ))z (θ+η)=1, with boundary constraints z (θ)=   1, − δ ≤ θ ≤ 0−1, 1≤ θ ≤ 1+η. 6. Conclusion A novel finite difference algorithm is suggested for solving SPDDE of second order with mixed shifts using a non-polynomial cubic spline with fitting factor. To represent the validity and efficiency of the method, we solved four test problems for different values N and with δ = 0.5ε = η and recorded computational results in the form of MAEs and ROCs. Using MATLAB, the MAEs in the solutions of the Examples 5.1, 5.2 and 5.3 are listed in comparison to the method given in [21] in Tables 1, 2, 3 and 4. Tables 5 and 6 compare the MAEs in Example 5.4 solution to the method described in [16]. The mixed shifts have no significant impact on the layer behaviour of the problems with boundary layers at the left-side and right-side of the points in the given interval shown in Figures. 1, 2, 3, 4, 5, 6, 7 and 8. Based on the results, we observe that the thickness of the layer increases as the size of the delay parameter increases and decreases as the size of the advance parameter increases. The proposed method is simple and can be easily implemented on a computer. Int. J. Anal. Appl. (2023), 21:5 9 Table 1. MAEs of Example 5.1 with various values of ε ε ↓ N → 23 24 25 26 27 28 Present method η = δ =0.5ε 10−1 3.769e-03 9.401e-04 2.209e-04 5.439e-05 1.359e-05 3.396e-06 2.0034 2.0889 2.0223 2.0002 1.7786 10−2 8.914e-03 4.578e-03 2.006e-03 7.512e-04 1.893e-04 4.425e-05 0.9611 1.1905 1.4168 1.9885 2.0968 10−3 8.968e-03 5.034e-03 2.683e-03 1.386e-03 6.778e-04 2.835e-04 0.8330 0.9076 0.9528 1.0324 1.2573 10−4 8.970e-03 5.036e-03 2.684e-03 1.388e-03 7.060e-04 3.561e-04 0.8334 0.9075 0.9517 0.9752 0.9874 10−5 8.970e-03 5.036e-03 2.684e-03 1.388e-03 7.060e-04 3.561e-04 0.8334 0.9075 0.9516 0.9752 0.9874 10−6 8.970e-03 5.036e-03 2.684e-03 1.388e-03 7.061e-04 3.561e-04 0.8334 0.9075 0.9516 0.9752 0.9874 Results in [21] η = δ =0.5ε 10−1 3.658e-03 9.595e-04 2.409e-04 6.759e-05 1.776e-05 1.232e-05 10−2 1.695e-02 7.297e-03 2.486e-03 6.964e-04 1.776e-04 2.616e-05 10−3 2.020e-02 1.047e-02 5.210e-03 2.461e-03 1.057e-03 3.771e-04 10−4 2.052e-02 1.079e-02 5.520e-03 2.769e-03 1.363e-03 6.539e-04 10−5 2.061e-02 1.088e-02 5.608e-03 2.858e-03 1.453e-03 7.417e-04 10−6 1.951e-02 9.783e-03 4.513e-03 1.762e-03 3.577e-04 3.729e-04 10 Int. J. Anal. Appl. (2023), 21:5 Table 2. MAEs in Example 5.2 ε ↓ N → 101 102 103 104 Present method η =0.5ε δ =0.7ε 10−1 1.243e-02 1.581e-04 1.585e-06 1.585e-08 10−2 2.466e-02 1.819e-03 2.077e-05 2.080e-07 10−3 2.488e-02 3.337e-03 1.914e-04 2.157e-06 10−4 2.490e-02 3.341e-03 3.458e-04 1.924e-05 Results in [21] η =0.5ε δ =0.7ε 10−1 1.533e-02 1.917e-04 1.921e-06 1.917e-08 10−2 2.817e-02 1.865e-03 2.024e-05 2.026e-07 10−3 2.853e-02 3.389e-03 1.919e-04 2.162e-06 10−4 2.857e-02 3.395e-03 3.463e-04 1.925e-05 Int. J. Anal. Appl. (2023), 21:5 11 Table 3. MAEs and ROCs in Example 5.2 ε ↓ N → 25 26 27 28 29 210 Present method η = δ =0.5ε 2−3 1.150e-03 2.919e-04 7.324e-05 1.832e-05 4.583e-06 1.145e-06 1.9789 1.9948 1.9987 1.9996 1.9999 2−4 2.582e-03 6.718e-04 1.697e-04 4.253e-05 1.064e-05 2.660e-06 1.9423 1.9851 1.9961 1.9990 1.9997 2−5 5.108e-03 1.437e-03 3.712e-04 9.358e-05 2.344e-05 5.864e-06 1.8290 1.9533 1.9880 1.9969 1.9992 2−6 8.048e-03 2.759e-03 7.653e-04 1.968e-04 4.955e-05 1.241e-05 1.5444 1.8501 1.9592 1.9896 1.9973 2−7 9.448e-03 4.293e-03 1.438e-03 3.959e-04 1.016e-04 2.556e-05 1.1381 1.5751 1.8615 1.9623 1.9904 2−8 9.619e-03 5.019e-03 2.221e-03 7.353e-04 2.015e-04 5.165e-05 0.9384 1.1762 1.5946 1.8674 1.9640 2−9 9.641e-03 5.099e-03 2.590e-03 1.130e-03 3.718e-04 1.017e-04 0.9191 0.9770 1.1965 1.6038 1.8703 Results in [21] η = δ =0.5ε 2−3 1.378e-03 3.486e-04 8.742e-05 2.187e-05 5.469e-06 1.367e-06 2−4 2.880e-03 7.458e-04 1.881e-04 4.714e-05 1.179e-05 2.948e-06 2−5 5.477e-03 1.526e-03 3.930e-04 9.902e-05 2.480e-05 6.204e-06 2−6 8.487e-03 2.862e-03 7.898e-04 2.028e-04 5.105e-05 1.278e-05 2−7 9.922e-03 4.413e-03 1.466e-03 4.024e-04 1.031e-04 2.596e-05 2−8 1.009e-02 5.148e-03 2.252e-03 7.424e-04 2.032e-04 5.206e-05 2−9 1.011e-02 5.228e-03 2.624e-03 1.138e-03 3.736e-04 1.021e-04 12 Int. J. Anal. Appl. (2023), 21:5 Table 4. MAEs and ROCs in Example 5.3 ε ↓ N → 25 26 27 28 29 210 Present method η = δ =0.5ε 2−3 5.102e-04 1.261e-04 3.161e-05 7.897e-06 1.973e-06 4.934e-07 2.0165 1.9960 2.0010 2.0002 2.0000 2−4 1.174e-03 2.837e-04 7.153e-05 1.784e-05 4.457e-06 1.114e-06 2.0489 1.9878 2.0034 2.0008 2.0002 2−5 2.676e-03 6.223e-04 1.550e-04 3.864e-05 9.656e-06 2.413e-06 2.1045 2.0049 2.0042 2.0007 2.0002 2−6 4.201e-03 1.413e-03 3.222e-04 8.159e-05 2.024e-05 5.067e-06 1.5715 2.1330 1.9815 2.0111 1.9978 2−7 4.818e-03 2.205e-03 7.283e-04 1.642e-04 4.195e-05 1.041e-05 1.1270 1.5986 2.1483 1.9692 2.0107 2−8 4.980e-03 2.520e-03 1.132e-03 3.700e-04 8.300e-05 2.129e-05 0.9823 1.1550 1.6131 2.1563 1.9629 2−9 4.997e-03 2.604e-03 1.290e-03 5.736e-04 1.865e-04 4.172e-05 0.9401 1.0127 1.1700 1.6206 2.1604 Results in [21] η = δ =0.5ε 2−3 8.434e-04 2.112e-04 5.284e-05 1.321e-05 3.303e-06 8.260e-07 2−4 4.172e-03 1.047e-03 2.640e-04 6.602e-05 1.650e-05 4.127e-06 2−5 1.858e-02 4.743e-03 1.190e-03 2.980e-04 7.452e-05 1.864e-05 2−6 6.074e-02 1.988e-02 5.080e-03 1.275e-03 3.192e-04 7.981e-05 2−7 1.111e-01 6.451e-02 2.061e-02 5.270e-03 1.323e-03 3.311e-04 2−8 1.297e-01 1.176e-01 6.658e-02 2.101e-02 5.372e-03 1.349e-03 2−9 1.310e-01 1.372e-01 1.212e-01 6.766e-02 2.122e-02 5.425e-03 Int. J. Anal. Appl. (2023), 21:5 13 Table 5. MAEs and ROCs in Example 5.4 ε ↓ N → 25 26 27 28 29 Present method η = δ =0.5ε 2−3 4.896e-03 1.129e-03 2.774e-04 6.934e-05 1.731e-05 2.1154 2.0257 2.0005 2.0018 2−4 1.114e-02 2.544e-03 5.859e-04 1.434e-04 3.587e-05 2.1312 2.1187 2.0302 1.9992 2−5 1.898e-02 5.677e-03 1.295e-03 2.979e-04 7.289e-05 1.7412 2.1315 2.1208 2.0309 2−6 2.371e-02 9.579e-03 2.864e-03 6.537e-04 1.501e-04 1.3076 1.7414 2.1317 2.1221 2−7 2.448e-02 1.188e-02 4.811e-03 1.438e-03 3.282e-04 1.0425 1.3047 1.7415 2.1319 2−8 2.449e-02 1.224e-02 5.950e-03 2.411e-03 7.210e-04 1.0000 1.0414 1.3032 1.7416 Results in [16] η = δ =0.5ε 2−3 8.354e-03 2.013e-03 4.986e-04 1.249e-04 3.121e-05 2−4 1.719e-02 4.378e-03 1.041e-03 2.571e-04 6.429e-05 2−5 2.517e-02 8.889e-03 2.238e-03 5.290e-04 1.303e-04 2−6 3.154e-02 1.294e-02 4.516e-03 1.131e-03 2.664e-04 2−7 4.478e-02 1.622e-02 6.559e-03 2.276e-03 5.686e-04 2−8 7.878e-02 2.317e-02 8.224e-03 3.301e-03 1.142e-03 14 Int. J. Anal. Appl. (2023), 21:5 Table 6. MAE in Example 5.4 with ε =0.1 N → 101 102 103 104 Present method δ ↓ η =0.5ε 0.00 8.123e-02 6.690e-03 3.573e-04 2.225e-05 0.05 8.066e-02 6.582e-03 3.529e-04 2.194e-05 0.09 8.019e-02 6.495e-03 3.494e-04 2.171e-05 η ↓ δ =0.5ε 0.00 8.051e-02 6.527e-03 3.508e-04 2.180e-05 0.05 8.066e-02 6.582e-03 3.529e-04 2.194e-05 0.09 8.077e-02 6.626e-03 3.546e-04 2.206e-05 Results in [16] δ ↓ η =0.5ε 0.00 9.109e-02 1.112e-02 6.382e-04 4.004e-05 0.05 9.047e-02 1.095e-02 6.306e-04 3.950e-05 0.09 8.996e-02 1.082e-02 6.244e-04 3.906e-05 η ↓ δ =0.5ε 0.00 9.604e-02 1.116e-02 6.458e-04 3.924e-05 0.05 9.621e-02 1.124e-02 6.494e-04 3.950e-05 0.09 9.634e-02 1.131e-02 6.522e-04 3.970e-05 Int. J. Anal. Appl. (2023), 21:5 15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 N um er ic al s ol ut io n of z ( ) =0.00 =0.05 =0.09 Figure 1. Layer profile in Example 5.1 with δ , N =25, ε =10−1 and η =0.5ε 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 N um er ic al s ol ut io n of z ( ) =0.00 =0.05 =0.09 Figure 2. Layer profile of Example 5.1 with η , N =25, ε =10−1 and δ =0.5ε 16 Int. J. Anal. Appl. (2023), 21:5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 N um er ic al s ol ut io n of z ( ) =0.00 =0.05 =0.09 Figure 3. Layer profile of Example 5.2 with δ , N =25, ε =10−1 and η =0.5ε 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 N um er ic al s ol ut io n of z ( ) =0.00 =0.05 =0.09 Figure 4. Layer profile of Example 5.2 with η , N =25, ε =10−1 and δ =0.5ε Int. J. Anal. Appl. (2023), 21:5 17 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 N um er ic al s ol ut io n of z ( ) =0.00 =0.05 =0.09 Figure 5. Layer profile of Example 5.3 with δ , N =25, ε =10−1 and η =0.5ε 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 N um er ic al s ol ut io n of z ( ) =0.00 =0.05 =0.09 Figure 6. Layer profile of Example 5.3 with η , N =25, ε =10−1 and δ =0.5ε 18 Int. J. Anal. Appl. (2023), 21:5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 N um er ic al s ol ut io n of z ( ) =0.00 =0.05 =0.09 Figure 7. Layer profile of Example 5.4 with δ , N =25, ε =10−1 and η =0.5ε 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 N um er ic al s ol ut io n of z ( ) =0.00 =0.05 =0.09 Figure 8. Layer profile of Example 5.4 with η , N =25, ε =10−1 and δ =0.5ε Int. J. Anal. Appl. (2023), 21:5 19 Conflicts of Interest: The authors declare that there are no conflicts of interest regarding the publi- cation of this paper. References [1] M. Adilaxmi, D. Bhargavi, Y. Reddy, An Initial Value Technique Using Exponentially Fitted Nonstandard Finite Difference Method for Singularly Perturbed Differential-Difference Equations, Appl. Math. 14 (2019), 245–269. [2] R.E. Bellman, K.L. Cooke, Differential-Difference Equations, Academic Press, New York, 1963. [3] M.W. Derstine, H.M. Gibbs, F.A. Hopf, D.L. Kaplan, Bifurcation Gap in a Hybrid Optically Bistable System, Phys. Rev. A. 26 (1982), 3720–3722. https://doi.org/10.1103/physreva.26.3720. [4] E.P. Doolan, J.J.H Miller, W.H. Schilders, Uniform Numerical Methods for Problems with Initial and Boundary Layers, Boole Press, Dublin, 1980. [5] R.D. Driver, Ordinary and Delay Differential Equations, Springer, New York, 1977. [6] M. Holmes, Introduction to Perturbation Methods, Springer, Berlin, 1995. [7] M.K. Kadalbajoo, R.K. Bawa, Variable-Mesh Difference Scheme for Singularly-Perturbed Boundary-Value Problems Using Splines, J. Optim. Theory Appl. 90 (1996), 405–416. https://doi.org/10.1007/bf02190005. [8] M.K. Kadalbajoo, K.K. Sharma, Numerical Analysis of Boundary-Value Problems for Singularly-Perturbed Differential-Difference Equations with Small Shifts of Mixed Type, J. Optim. Theory Appl. 115 (2002), 145–163. https://doi.org/10.1023/a:1019681130824. [9] M.K. Kadalbajoo, K.K. Sharma, Numerical Treatment of a Mathematical Model Arising From a Model of Neuronal Variability, J. Math. Anal. Appl. 307 (2005), 606–627. https://doi.org/10.1016/j.jmaa.2005.02.014. [10] D. Kumara Swamy, K. Phaneendra, Y.N. Reddy, Accurate Numerical Method for Singularly Perturbed Differential- Difference Equations with Mixed Shifts, Khayyam J. Math. 4 (2018), 110–122. https://doi.org/10.22034/kjm. 2018.57949. [11] C.G. Lange, R.M. Miura, Singular Perturbation Analysis of Boundary Value Problems for Differential-Difference Equations, SIAM J. Appl. Math. 42 (1982), 502–531. https://doi.org/10.1137/0142036. [12] C.G. Lange, R.M. Miura, Singular Perturbation Analysis of Boundary Value Problems for Differential-Difference Equations III. Turning Point Problems, SIAM J. Appl. Math. 45 (1985), 708–734. https://doi.org/10.1137/ 0145042. [13] C.G. Lange, R.M. Miura, Singular Perturbation Analysis of Boundary Value Problems for Differential-Difference Equations. V. Small Shifts with Layer Behavior, SIAM J. Appl. Math. 54 (1994), 249–272. https://doi.org/10. 1137/s0036139992228120. [14] C.G. Lange, R.M. Miura, Singular Perturbation Analysis of Boundary-Value Problems for Differential-Difference Equations. VI. Small Shifts with Rapid Oscillations, SIAM J. Appl. Math. 54 (1994), 273–283. https://doi.org/ 10.1137/s0036139992228119. [15] A. Longtin, J.G. Milton, Complex Oscillations in the Human Pupil Light Reflex With "Mixed" and Delayed Feedback, Math. Biosci. 90 (1988), 183–199. https://doi.org/10.1016/0025-5564(88)90064-8. [16] M. Ayalew, G.G. Kiltu, G.F. Duressa, Fitted Numerical Scheme for Second-Order Singularly Perturbed Differential- Difference Equations with Mixed Shifts, Abstr. Appl. Anal. 2021 (2021), 4573847. https://doi.org/10.1155/ 2021/4573847. [17] A.H. Nayfeh, Perturbation Methods, Wiley, New York, 1979. [18] R.E. O’Malley Jr, Introduction to Singular Perturbations, Academic Press, New York, 1974. [19] S. Priyadarshana, S.R. Sahu, J. Mohapatra, Asymptotic and Numerical Methods for Solving Singularly Perturbed Differential Difference Equations With Mixed Shifts, Iran. J. Numer. Anal. Optim. 12 (2022), 55–72. https: //doi.org/10.22067/ijnao.2021.70731.1038. https://doi.org/10.1103/physreva.26.3720 https://doi.org/10.1007/bf02190005 https://doi.org/10.1023/a:1019681130824 https://doi.org/10.1016/j.jmaa.2005.02.014 https://doi.org/10.22034/kjm.2018.57949 https://doi.org/10.22034/kjm.2018.57949 https://doi.org/10.1137/0142036 https://doi.org/10.1137/0145042 https://doi.org/10.1137/0145042 https://doi.org/10.1137/s0036139992228120 https://doi.org/10.1137/s0036139992228120 https://doi.org/10.1137/s0036139992228119 https://doi.org/10.1137/s0036139992228119 https://doi.org/10.1016/0025-5564(88)90064-8 https://doi.org/10.1155/2021/4573847 https://doi.org/10.1155/2021/4573847 https://doi.org/10.22067/ijnao.2021.70731.1038 https://doi.org/10.22067/ijnao.2021.70731.1038 20 Int. J. Anal. Appl. (2023), 21:5 [20] A. Raza, A. Khan, Treatment of Singularly Perturbed Differential Equations with Delay and Shift Using Haar Wavelet Collocation Method, Tamkang J. Math. 53 (2021), 303–322. https://doi.org/10.5556/j.tkjm.53.2022.3250. [21] R. Ranjan, H.S. Prasad, A Novel Approach for the Numerical Approximation to the Solution of Singularly Perturbed Differential-Difference Equations With Small Shifts, J. Appl. Math. Comput. 65 (2020), 403–427. https://doi. org/10.1007/s12190-020-01397-6. [22] R. Ranjan, H.S. Prasad, A Novel Exponentially Fitted Finite Difference Method for a Class of 2nd Order Singularly Perturbed Boundary Value Problems With a Simple Turning Point Exhibiting Twin Boundary Layers, J. Ambient Intell. Human Comput. 13 (2022), 4207–4221. https://doi.org/10.1007/s12652-022-03902-0. [23] L.S. Senthilkumar, R. Mahendran, V. Subburayan, A Second Order Convergent Initial Value Method for Singularly Perturbed System of Differential-Difference Equations of Convection Diffusion Type, J. Math. Computer Sci. 25 (2021), 73–83. https://doi.org/10.22436/jmcs.025.01.06. [24] L. Sirisha, K. Phaneendra, Y.N. Reddy, Mixed Finite Difference Method for Singularly Perturbed Differential Difference Equations With Mixed Shifts via Domain Decomposition, Ain Shams Eng. J. 9 (2018), 647–654. https://doi.org/10.1016/j.asej.2016.03.009. [25] R.B. Stein, Some Models of Neuronal Variability, Biophys. J. 7 (1967), 37–68. https://doi.org/10.1016/ s0006-3495(67)86574-3. [26] R.S. Varga, Matrix Iterative Analysis, Prentice-Hall, Englewood Cliffs, 1962. [27] M. Wazewska, A. Lasota, Mathematical Models of the Red Cell System, Mat. Stosowana, 6 (1976), 25–40. https://doi.org/10.5556/j.tkjm.53.2022.3250 https://doi.org/10.1007/s12190-020-01397-6 https://doi.org/10.1007/s12190-020-01397-6 https://doi.org/10.1007/s12652-022-03902-0 https://doi.org/10.22436/jmcs.025.01.06 https://doi.org/10.1016/j.asej.2016.03.009 https://doi.org/10.1016/s0006-3495(67)86574-3 https://doi.org/10.1016/s0006-3495(67)86574-3 1. Introduction 2. Definition of the problem and Non-Polynomial Spline 3. Numerical Algorithm 4. Convergence Analysis 5. Numerical Examples 6. Conclusion References