# Uniformly Convergent Nonpolynomial Spline Method for Singularly Perturbed Robin-Type Boundary Value Problems with Discontinuous Source Term

Authors: Habtamu Garoma Debela; Gemechis File Duressa
Journal: Abstract and Applied Analysis
 (2021)

---

## Abstract

In this paper, a singularly perturbed second-order ordinary differential equation with discontinuous source term subject to mixed-type boundary conditions is considered. A fitted nonpolynomial spline method is suggested. The stability and parameter uniform convergence of the proposed method are proved. To validate the applicability of the scheme, two model problems are considered for numerical experimentation and solved for different values of the perturbation parameter,ε, and mesh size, h. The numerical results are tabulated in terms of maximum absolute errors and rate of convergence, and it is observed that the present method is more accurate and ε-uniformly convergent for h≥ε where the classical numerical methods fail to give good result and it also improves the results of the methods existing in the literature.


---

## Body

## 1. Introduction

Singular perturbation problems (SPPs) are differential equations with a small positive parameter multiplying the highest derivative term. These problems arise in several branches of applied mathematics, including fluid dynamics, quantum mechanics, elasticity, chemical reactor theory, and gas porous electrodes theory. Examples of SPPs include the Navier-Stokes equation of fluid flow at high Reynolds number, the equation governing flow in porous media, the drift-diffusion equation of semiconductor devices, physics and mathematical models of liquid crystal material, and the convection-diffusion and reaction-diffusion equations to mention but a few [1, 2].Such equations typically exhibit solutions with layers, which cause severe computational difficulties for standard numerical methods. Consequently, a variety of different numerical strategies have been devoted [1–3] and the references are therein to the construction and analysis of accurate numerical methods for SPPs. Recently, authors [4–8] have considered SPPs for second-/third-order ordinary differential equations (ODEs) with discontinuous source term and/or discontinuous convection coefficient. The novel aspect of the problem under consideration is that we take a source term in the differential equation which has a jump discontinuity at one or more points in the interior of the domain. This gives rise to an interior layer in the exact solution of the problem, in addition to the boundary layer at the outflow boundary point.A single discontinuity is assumed to occur at a pointd∈Ω=0,1. The solution of this problem has a boundary layer at x=0 and interior layers with different widths at x=d. It is convenient to introduce the notation Ω−=0,d and Ω+=d,1 and to denote the jump at d in any function with wd=wd+−wd−.Recently, Chandru and Shanthi [9] presented a fitted mesh method to solve singularly perturbed robin-type boundary value problems with discontinuous source term. But still, there is a room to increase the accuracy and develop its convergence for the developed scheme because the treatment of singular perturbation problem is not trivial distributions and the solution pended on perturbation parameter ε and mesh size h [10]. Due to this, numerical treatment of singularly perturbed boundary value problems needs improvement.Different authors in [11–23] have applied different families of spline method for second- and higher order singular and regular problems. As far as the researchers’ knowledge is concerned, numerical solution of fitted nonpolynomial cubic spline method for the problem under consideration is first being considered. Hence, in the present paper, motivated by the works of [9], we developed a fitted nonpolynomial cubic spline method for singularly perturbed robin-type boundary value problems with discontinuous source term. Therefore, it is important to develop more accurate and convergent numerical method for solving singularly perturbed robin boundary value problems. Thus, the purpose of this study is to develop stable, convergent, and more accurate numerical method for solving singularly perturbed robin boundary value problems with discontinuous source term.In this article, we consider the following class of problems:(1)Lyx≡εy″x+axy′x−bxyx=fx,x∈Ω−∪Ω+,subject to the boundary conditions,
(2)L1y0=α1y0−β1εy′0=ϕ0,L2y1=α2y1+β2y′1=φ,where α1,β1≥0,α1+β1>0,β2≥0,α2>0, and ε0<ε≪1 are a singular perturbation parameter; ax and bx are sufficiently smooth functions on Ω¯− and Ω¯+, respectively; and ax≥a>0 and bx≥b≥0. It is assumed that f is a sufficiently smooth function in Ω−∪Ω+∪0,1; the left and right limits of f and their derivatives are assumed to exist at x=d. The function f is assumed to have simple discontinuity at x=d. Hence, the solution y of equations (1) and (2) does not necessarily have a continuous second derivative at the point d; that is, y does not belong to the class of functions C2Ω. Hence, the class of functions, where y belongs to it, is taken as C0Ω¯∩C1Ω∩C2Ω−∪Ω+. Further, the SPPs in equations (1) and (2) have a unique solution y∈Y≡C0Ω¯∩C1Ω∩C2Ω−∪Ω+ [24].Boundary value problems of the type in equations (1) and (2) are model confinement of a plasma column by reaction pressure and geophysical fluid dynamics [25].


## 2. Properties of Continuous Solution

The following lemmas are necessary for the existence and uniqueness of the solution and for the problem to be well-posed.Lemma 1 (continues minimum principle).
Suppose that the functiony∈C1Ω¯∩C2Ω−∪Ω+ satisfies L1y0≥0 and L2y1≥0 and Lyx≤0,∀x∈Ω−∪Ω+, and y′d≤0,thenyx≥0,∀x∈Ω¯..Proof.
For the proof, refer to [26].Lemma 2 (stability result).
Consider the boundary value problem of equations (1) and (2) subject to the condition ax≥a>0,bx≥b≥0. If y∈C1Ω¯∩C2Ω−∪Ω+, then
(3)yΩ¯≤CmaxL1y0,L2y1,LyΩ−∪Ω+.Proof.
For the proof, refer to [26].Lemma 3.
For each integerk, satisfying 0≤k≤4, the solution y of equations (1) and (2) satisfy the bounds
(4)ykΩ¯\d≤Cε−k.Proof.
To establish the parameter robust properties of the numerical methods involved in this paper, the following decomposition ofY into smooth v and singular w components is required.The smooth componentv is defined as the solution of
(5)av0′−bv0=f,x∈Ω−∪Ω+,L1v00=q,av1′−bv1=−v0″,x∈Ω−∪Ω+,L1v11=0.Note thatv=v0=εv1+ε2v2, where v2∈C1Ω¯∩C2Ω−∪Ω+ and
(6)Lv2=−v2″,L1v20=v2d=L2v21=0.As in [27], we can obtain the following bounds on the derivatives of v for k=1,2,3:
(7)v≤CvkΩ−∪Ω+≤C1+ε2−k.Note also thatv′d≤C,v″d≤C. Define the singular components of the decomposition as follows. Find w∈C0Ω such that Lw=0,x∈Ω−∪Ω+,(8)L1w0=y0−v0,L2w1=0,w′d=−v′d,We can further decomposew as
(9)w=w1+w2,where w1 is the boundary layer function satisfying w1∈C2Ω,Lw2=0,x∈Ω,(10)L1w0=y0−v0,L2w1=0,and w2 is the interior layer function satisfying
(11)w2∈C0Ω,Lw2=0,x∈Ω−∪Ω+,L1w20=0,L2w20=0,w2′d=−v′d.(12)yεkx≤C1+ε−kexp−axε,forallx∈Ω−.Lemma 4.
Letyε be the solution of Pε. Then, for k=0,1,2,3,Proof.
For the proof, refer to [26].


## 3. Formulation of the Numerical Method

Consider the following singularly perturbed robin-type boundary value problems with discontinuous source terms of the form(13)εy″x+axy′x−bxyx=fx.Rearranging equation (13), we obtain
(14)εy″x+pxy′x+qxyx=fx,where px=ax and qx=−bx.Under boundary conditions,(15)y0=ϕ0,(16)y1=φ.Consider a uniform mesh with interval0,1 in which 0=x0<x1<x2<⋯<xN=1, whereh=1/Nandxi=ih,i=0,1,2,⋯,N.For each segmentxi,xi+1,i=1,2,⋯,N−1, the nonpolynomial cubic spline SΔx has the following form:
(17)SΔx=aisinwx−xi+bicoswx−xi+ciewx−xi−e−wx−xi+diewx−xi+e−wx−xi,where ai,bi,ci,anddi are unknown coefficients and w≠0 arbitrary parameter which will be used to increase the accuracy of the method of Bawa [28].To determine the unknown coefficients in equation (17) in terms of yi,yi+1,Mi,andMi+1, first, we define
(18)SΔx=yi,SΔxi+1=yi+1,SΔ″xi=Mi,SΔ″xi+1=Mi+1.By some algebra, the coefficients in equation ((17)) becomes
(19)ai=w2yi+1−Mi+1+Mi−yiw2cosθ2w2sinθ,bi=yiw2−Mi2w2,ci=2yi+1w2+Mi+1−yiw2+Mieθ+e−θ4w2eθ−e−θ,di=yiw2+Mi4w2,where θ=wh.Using the continuity condition of the first derivative atxi, SΔ−1′xi=SΔ′xi, we have
(20)ai−1cosθ−bi−1sinθ+ci−1eθ+e−θ+di−1eθ−e−θ=ai+2ci.Reducing indices of equation (19) by one and substituting into equation (20), we have
(21)−12sinθ−1eθ−e−θyi−1+cosθsinθ+eθ+e−θeθ−e−θyi+−12sinθ−1eθ−e−θyi+1=−h22θ2sinθ−h2θ2eθ−e−θMi−1+h2cosθθ2sinθ+h2eθ+e−θθ2eθ−e−θMi+−h22θ2sinθ−h2θ2eθ−e−θMi+1.Multiplying both side of equation (21) by −2sinθeθ−e−θ/eθ−e−θ+2sinθ, we obtain
(22)yi−1+kyi+yi+1=h2αMi−1+βMi+αMi+1,where
(23)k=−2eθcosθ+sinθ+e−θsinθ−cosθeθ−e−θ+2sinθ,α=eθ−e−θ−2sinθθ2eθ−e−θ+2sinθ,β=−2eθh2cosθ−sinθ−e−θsinθ+cosθθ2eθ−e−θ+2sinθ.Asθ⟶0,α,β,k⟶1/6,4/6,−2.SubstitutingSΔ″xi=yi″=Mi in equation (14), we get
(24)εMi=fi−piyi′−qiyi,εMi−1=fi−1−pi−1yi−1′−qi−1yi−1,εMi+1=fi+1−pi+1yi+1′−qi+1yi+1.Substituting equation (24) into equation (22), we obtain
(25)εh2yi−1+kyi+yi+1=αfi−1−pi−1yi−1′−qi−1yi−1+βfi−piyi′−qiyi+αfi+1−pi+1yi+1′−qi+1yi+1and the following approximations for the first derivatives of y by Bawa [28]:
(26)yi−1′≈−yi+1+4yi−3yi−12h,yi+1′≈3yi+1−4yi+yi−12h,yi′≈1−2h2ψqi+1−hψ3pi+1+pi−12hyi+1+2ψpi+1+pi−1yi,+1−2h2ψqi−1−hψpi+1+3pi−12hyi−1+hψri+1−ri−1,where ω is the parameter used to raise the accuracy of the method.Substituting equations (26) into equation (25) and rearranging, we get
(27)εh2yi−1+kyi+yi+1=3αpi−12h−αqi−1−βpi2h−1+2ψh2qi−1−ψhpi+1+3pi−1+αpi+12hyi−1+−2αpi−1h−2βpiψpi+1+pi−1−βpiqiαqi−1−2αpi+1hyi+3αpi+12h+αqi+1−βpi2h1−2ψh2qi+1−ψh3pi+1+pi−1+αpi−12hyi+1+α+βpiψhfi−1+βfi+α−βpiψhfi+1.From the theory of singular perturbation, it is known that the solution of equation (14) and equation (15) is of the form (O’Malley [29])
(28)yx=y0x+p0pxϕ0−y00exp−∫0xpxε−qxpxdx+Oε,where y0x is the solution of the reduced problem (obtained by setting ε=0 of equation (14) which is given by
(29)pxy0′x+qxy0x=fx,withy01=φ.Expandingpx and qx about point “0” up to the first term using Taylor’s series, the discretized form of equation (28) gives
(30)yxi=y0xi+ϕ0−y00exp−p0ε−q0p0xi+Oε.Therefore,(31)limh⟶∞yih=y00+ϕ0−y00exp−p20−εq0p0iρ+Oε, where ρ=h/ε.To handle the effect of the perturbation parameter in equation (27), artificial viscosity (σρ) is multiplied on the term containing the perturbation parameter as
(32)εσρh2yi−1+kyi+yi+1=3αpi−12h−αqi−1−βpi2h−1+2ψh2qi−1−ψhpi+1+3pi−1+αpi+12hyi−1+−2αpi−1h−2βpiψpi+1+pi−1−βpiqiαqi−1−2αpi+1hyi+3αpi+12h+αqi+1−βpi2h1−2ψh2qi+1−ψh3pi+1+pi−1+αpi−12hyi+1+α+βpiψhfi−1+βfi+α−βpiψhfi+1.Multiplying equation (32) by h and taking limit as h⟶0, we get the fitting factor
(33)σρ=ρp02ϕ0+φ2cothp20−εq0p0ρ2,which is a required fitting factor in the left boundary layer. From equation (32), we get
(34)εσ−h2ui+αh2qi−1yi−1−−ρεσ+h2vi−βh2qiyi+εσ−h2wi+αh2qi+1yi+1=h2α+βpiψhri−1+βri+α−βpiψhri+1,where
(35)ui=3αpi−1+βpiψhpi+1+3pi−1−βpiψh2qi−1−αpi+1+βpi,vi=−4αpi−1−4βhpiψpi+1+pi−1+4αpi+1,wi=αpi−1+βpiψh3pi+1+pi−1−3αpi+1−βpi+βpiψh2qi+1.Then, we get the three term recurrence relation of the form(36)Eiyi−1−Fiyi+Giyi+1=Hi,i=1,2,⋯,N−1,where
(37)Ei=εσ−h2ui+αh2qi−1,Fi=−ρεσ+h2vi−βh2qi,Gi=εσ−h2wi+αh2qi+1,Hi=h2α+βpiψhfi−1+βfi+α−βpiψhfi+1.To treat the boundary condition, we used forward finite difference formula fori=0 and backward difference formula for i=N, respectively, for the first derivative term.That is, fori=0, from equation (2), we have
(38)α1y0−β1εy′0=ϕ0⇒α1y0−β1εy0′=ϕ0,⇒α1+β1εhy0−β1εhy1=ϕ0.Similarly, fori=N, from equation (2), we have
(39)α2yN+β2y′N=φ⇒α2yN+β2yN′=φ,⇒−β2hyN−1+α2+β2hyN=φ.Therefore, the problem in equation (1) with given boundary condition in equation (2) can be solved using the scheme equations (36), (38), and (39) which form an N×N system of algebraic equations.


## 4. Truncation Error

From equation (26), we have
(40)ei+1′=Y′xi+1−yi+1′=h2yi33+h3yi412+h4yi530+Oh5,ei−1′=Y′xi−1−yi−1′=h2yi33−h3yi412+h4yi530+Oh5,ei′=Y′xi−yi′=−h216+2ψεσyi3−h4ψεσyi53+h4yi5120+Oh5.From equation (24) and equation (22),
(41)εσyi−1+εσkyi+εσyi+1=h2αfi−1−pi−1yi−1′−qi−1yi−1+βfi−piyi′−qiy+αfi+1−pi+1yi+1′−qi+1yi+1.Putting the exact solution in equation (41), we have
(42)εσYxi−1+εσkYxi+εσYxi+1=h2αfi−1−pi−1Y′xi−1−qi−1Yxi−1+βfi−piY′xi−qiYxi+h2αfi+1−pi+1Y′xi+1−qi+1Yxi+1+T0h,where
(43)T0h=−2+kεσyi+2α+β−1εσh2yi″+12α−1εσh4yi412+2α+β−1εσh6yi6.Subtract equation (41) from equation (42) and putting ej=Yxj−yj,j=i±1, we get
(44)εσ+h2αqi−1ei−1+εσk+βh2qiei+εσ+h2αqi+1ei+1=−h2αpi−1ei−1′+βpiei′+αpi+1ei+1′+T0h.Using equation (40), we have
(45)εσ+h2αqi−1ei−1+kσε+βh2qiei+εσ+h2αqi+1ei+1=−h4αpi−13+βpi16+2ψεσ+αpi+13yi3−h5−αpi−112+αpi+112yi4−h6130αpi−1+αpi+1y5+βpi−ψεσ3+1120yi5+T0h.Let(46)pi+1=pi+hpi′+h2pi″2+Oh3,pi−1=pi−hpi′+h2pi″2+Oh3.Substituting equation (46) into equation (45), we have
(47)εσ+h2αqi−1ei−1+kεσ+βh2qiei+εσ+h2αqi+1ei+1=Tih,where
(48)Tih=−2+kεσyi+2α+β−1εσh2yi″+−2α3+β16+2ψεσpih4yi3+12α−1εσh4yi412+Oh6.Tih is a local truncation error associated with the scheme developed in equation (36). Thus, for different values of k,α,β,andψ in the scheme of equation (36), the following different orders are indicated:
(1)
Fork=−2,2α+β=1, and for any value of ψ, the scheme of equation (36) gives the second-order method(2)
Fork=−2,α=1/12,β=10/12, and ψ=−1/20, from equation (36), the fourth-order method is derived


## 5. Stability and Convergence Analysis

Writing the tridiagonal system in equation (36) in matrix vector form, we have
(49)AY¯=C,where A=B0+BandB=B1+h2B2Q is a tridiagonal matrix of order N−1. Multiplying both sides of equation (34) by (-1), we get
(50)B0=−kσε−σε−σε−kσε−σε⋱⋱⋱−σε−kσε−σε−σε−kσε,B1=hv12hw12hu22hv22hw22⋱⋱⋱huN−22hvN−22hwN−22huN−12hvN−12,B2=−β−α−α−β−α⋱⋱⋱−α−β−α−α−β,Q=q1q2⋮qN−2qN−1.C=ci is the column vector, where
(51)ci=h2−α−βp1ψhf0−βf1−α−βp1ψhf2+εσϕ0−hϕ02u1,fori=1.ci=h2−α−βpiψhfi−1−βfi−α−βpiψhfi+1,for2≤i≤N−2.ci=h2−α−βpn−1ψhfN−2−βfN−1−α−βpN−1ψhfN+αh2qN+εσφ−hφ2wN,fori=N−1.Y¯=y¯1,y¯2,⋯,y¯N−1t is an approximation solution; Th=Oh6 for k=−2, α=1/12, β=10/12, and ψ=−1/20ε.Now, considering the above system with exact solutionY=yx1,yx2,⋯,yxN−1t, we have
(52)AY−Th=C,whereTh=T1h,T2h,⋯,TN−1ht is a local truncation error.From equation (49) and equation (52), we obtain
(53)AY−Y¯=Th.This gives the error equation(54)AE=Th,where
(55)E=Y−Y¯=e1,e2,⋯,eN−1t.Let(56)pi−1≤Z1,pi≤Z2,pi+1≤Z3,qi−1≤K1,qi≤K2,qi+1≤K3,where Z1,Z2,Z3,K1,K2,andK3 are positive constants.Letrij be the i,jth element of the matrix B; then, for i=1,2,⋯,N−1,
(57)ri,i+1=h2wi−αh2qi+1≤h2αZ1+βZ21+2ψh2K3+ψh3Z3+Z1+3αZ3+2αhK3,ri,i−1=h2ui−αh2qi−1≤h23αZ1+βZ21+2ωh2K1+ωhZ3+3Z1+αZ3+2αhK1..Thus, for sufficiently smallh, we have
(58)−εσ+ri,i+1≠0fori=1,2,⋯,N−2,−εσ+ri,i−1≠0fori=2,⋯,N−1.Hence, the matrixA is irreducible (Varga [30]).From equation (36), for sufficiently small handk=−2andε,0<ε≪1, and we get Ei+Gi<Fi.Hence,A is diagonally dominant. Under this condition, Thomas algorithm method is stable (Kadalbajoo and Reddy [31]).LetSi be the sum of the elements of the ith row of the matrix A; then, we have
(59)Si=−kσε−εσ+h2−3αpi−1−βpi+αpi+1+h22βψpipi+1+3pi−1−2βqi+αqi+1−h3βψpiqi+1,fori=1,Si=−kσε−2σε+h2−αqi−1−βqi−αqi+1+Oh3,for2≤i≤N−2,Si=−kσε−εσ+h2−αpi−1+βpi+3αpi+1+h22βψpipi+1+pi−1−2βqi+αqi−1+h3βψpiqi−1,fori=N−1.LetZ2∗=min1≤i≤N−1pi,Z2∗=max1≤i≤N−1pi,K2∗=min1≤i≤N−1qi,andK2∗=max1≤i≤N−1qi, then 0<Z2∗≤Z2≤Z2∗ and0<K2∗≤K2≤K2∗.For sufficiently smallh, A is monotone (Young [32]).Hence,A−1 exist and A−1≥0. From equation (63), we have
(60)E=A−1Th,E≤A−1Th.For sufficiently smallh, we have
(61)S1>h2K2∗D1,fori=1,Si>h2K2∗D2,for2≤i≤N−2,Sn−1>h2K2∗D1,fori=N−1,where D1=α+βandD2=2α+β.LetAi,k−1 be the i,kth element of A−1, and we define
(62)Ai,k−1=max1≤i≤N−1∑k=1n−1Ai,k−1,(63)Th=max1≤i≤N−1Tih.SinceAi,k−1≥0, then from the theory of matrices, we have
(64)∑k=1N−1Ai,k−1⋅Sk=1fori=1,2,⋯,N−1.Hence,(65)Ai,1−1≤1S1<1h2K2∗D1fork=1,(66)Ai,N−1−1≤1SN−1<1h2K2∗D1fork=N−1,(67)∑k=2N−2Ai,k−1≤1min2≤k≤N−2Sk<1h2K2∗D2fork=2,3,⋯,N−2.From equations (54)–(68) and (48), we get
(68)E≤1h2K2∗D1+1h2K2∗D1+1h2K2∗D2Oh6=D∗h4,where
(69)D∗=1K2∗D1+1K2∗D1+1K2∗D2,which is independent of mesh size h. This establishes that the present method is fourth-order convergent for k=−2,α=1/12,β=10/12,andψ=−1/20ε.


## 6. Numerical Examples and Results

To validate the applicability of the method, two model problems are considered for numerical experimentation for different values of the perturbation parameter and mesh points. The numerical results are tabulated in terms of maximum absolute errors, numerical rate of convergence, and uniform errors (see Tables1–4) and compared with the results of the previously developed numerical methods existing in the literature (Tables 2 and 4). Further, behavior of the numerical solution without and with fitted operator (Figures 1 and 2), point-wise absolute errors (Figure 3), and the ε-uniform convergence of the method are shown by the log-log plot of the ε-uniform error (Figure 4).Table 1
Maximum absolute errors for Example1 at different mesh points N.
εN=32N=64N=128N=256N=51210−41.0156e-025.0781e-032.5391e-031.2695e-036.2999e-0410−81.0156e-025.0781e-032.5391e-031.2695e-036.2999e-0410−121.0156e-025.0781e-032.5391e-031.2695e-036.2999e-0410−161.0156e-025.0781e-032.5391e-031.2695e-036.2999e-0410−201.0156e-025.0781e-032.5391e-031.2695e-036.2999e-04Table 2
Comparison of maximum absolute errors and order of convergence for Example1 at different of mesh points N.
N⟶641282565121024Present methodEN5.0781e-032.5382e-031.2467e-035.5910e-042.3203e-04RN1.00051.26671.15691.2688Method in [9]EN2.5658e−021.4128e−027.5359e−033.8225e−031.7536e−03RN0.860850.906710.979271.1242Table 3
Maximum absolute errors for Example2 at different mesh points.
ε↓N⟶N=32N=64N=128N=256N=51210−44.1730e-022.0931e-021.0482e-025.2450e-032.4856e-0310−84.1730e-022.0931e-021.0482e-025.2450e-032.4856e-0310−124.1730e-022.0931e-021.0482e-025.2450e-032.4856e-0310−164.1730e-022.0931e-021.0482e-025.2450e-032.4856e-0310−204.1730e-022.0931e-021.0482e-025.2450e-032.4856e-03Table 4
Comparison of maximum absolute errors and order of convergence for Example2 at number of mesh points N.
N⟶641282565121024Present methodEN2.0929e-021.0400e-024.7763e-031.9081e-037.6337e-04RN1.00891.12261.32381.3217Method in [9]EN9.6698e−015.8056e−013.2795e−011.7313e−018.1501e−02RN0.736040.823990.921611.08700Figure 1
Numerical behavior of Example1 without and with fitted operator using nonpolynomial cubic spline method at N=128 and ε=2−10.Figure 2
Numerical behavior of Example2 without and with fitted operator using nonpolynomial cubic spline method at N=128 and ε=2−10.Figure 3
Point-wise absolute error of Example1 and Example 2 with fitted operator using nonpolynomial cubic spline method, respectively, for N=128 and ε=2−10.Figure 4
ε-uniform convergence with log-log scale for Example 1 and Example 2, respectively.Havingyj≡yjN (the approx. solution obtained via the present method) for different values of N and ε since the exact solution is not available, the maximum errors (denoted by EεN) are evaluated using the formula given by the double-mesh principle:
(70)EεN≔max0≤j≤NyjN−y2j2N.Further, we will tabulate the errors(71)EN=max0<ε≤1EεN.The numerical rates of convergence are computed using the formula(72)rεN≔loglog2EεN/Eε2N,and the numerical rate of “ε-uniform convergence” is computed using
(73)RN=log2EN/E2N.Example 1.
(74)εy″x+y′x=fx,x∈Ω−∪Ω+,y0−εy′0=1,y1+y′1=−1,where
(75)fx=0.7,for0≤x≤0.5,−0.6,for0.5<x≤1.Example 2.
(76)εy″x+11+xy′x=fx,x∈Ω−∪Ω+,y0−εy′0=1,y1+y′1=1,where
(77)fx=1+x,for0≤x≤0.5,4.0,for0.5<x≤1.


## 7. Discussion and Conclusion

This study introduces uniformly convergent nonpolynomial cubic spline method based on exponential fitted operator for solving singularly perturbed second-order ODEs of Robin-type BVPs with discontinuous source term. Due to discontinuity in the source term, there is an interior layer occurring. To fit the interior and boundary layer, a suitable fitted operator finite difference method on uniform mesh is constructed. The behavior of the continuous solution of the problem is studied and shown that it satisfies the continuous stability estimate and the derivatives of the solution are also bounded. The numerical scheme is developed on uniform mesh. The Robin-type BVPs are treated using numerical finite difference techniques, and the results are compared accordingly. The stability of the developed scheme is established, and its uniform convergence is proved. Unlike other fitted operator finite difference methods constructed in standard ways, the method that we presented in this paper is fairly simple to construct.



---

*Source: 7569209-2021-10-22.xml*