# Solving Systems of Singularly Perturbed Convection Diffusion Problems via Initial Value Method

Authors: Wondwosen Gebeyaw Melesse; Awoke Andargie Tiruneh; Getachew Adamu Derese
Journal: Journal of Applied Mathematics
 (2020)

---

## Abstract

In this paper, an initial value method for solving a weakly coupled system of two second-order singularly perturbed Convection–diffusion problems exhibiting a boundary layer at one end is proposed. In this approach, the approximate solution for the given problem is obtained by solving, a coupled system of initial value problem (namely, the reduced system), and two decoupled initial value problems (namely, the layer correction problems), which are easily deduced from the given system of equations. Both the reduced system and the layer correction problems are independent of perturbation parameter,ε. These problems are then solved analytically and/or numerically, and those solutions are combined to give an approximate solution to the problem. Further, error estimates are derived and examples are provided to illustrate the method.


---

## Body

## 1. Introduction

Singular perturbation problems (SPPs) arise most frequently in diversified fields of applied mathematics. For instance, in fluid mechanics, elasticity, aerodynamics, quantum mechanics, chemical-reactor theory, oceanography, meteorology, modeling of semiconductor devices, and many others in the area. A well known fact is that, the solutions of such problems have a multiscale character, i.e., there are thin transition layer(s) where the solution varies very rapidly, while away from the layer(s) the solution behaves regularly and varies slowly. This leads to boundary and/or interior layer(s) in the solution of the problems. For a detailed discussion on the analytical and numerical treatment of SPPs one can refer the books of Miller et al. [1], O’Malley [2] and Roos et al. [3].Due to the presence of the layer regions, it has been shown that the classical numerical methods fails to produce good approximations for singular perturbation problems (SPPs). In fact, some numerical techniques employed for solving second-order singularly perturbed boundary value problem (SPBVPs) are based on the idea of replacing this problems by suitable initial value problems (IVPs). The reason for that is, the numerical treatment of a boundary value problem is much more demanding than the treatment of the corresponding IVPs. There are different initial value methods in the literature of SPPs developed for solving SPBVPs, for the detail discussions of such methods one can refer the papers [4–8] and the references there in.In the past few decades, a considerable amount of works have been reported in the literature of SPPs. However, most of the works connected with the computational aspects are confined to second-order equation. Only few results are reported for higher order and systems of equations. The systems of SPPs have applications in electro analytical chemistry, predator prey population dynamics, modeling of optimal control situations and resistance-capacitor electrical circuits [9]. In recent years, few scholars developed non-classical methods for different classes of systems of singularly perturbed differential equations. A class of systems of singularly perturbed reaction-diffusion equations have been examined by the authors in [10, 11–14]. In the papers [15, 16], a class of strongly coupled systems of singularly perturbed convection–diffusion equations are examined. The scholars in [17–21], considered weakly coupled systems of singularly perturbed convection–diffusion equations with equal or different diffusion parameters. A brief survey of article on the current progress about the numerical treatment of systems of singularly perturbed differential equations is also discussed in [22]. However, most of the methods developed for systems of singularly perturbed problems focus on fitted mesh method, so it is natural to look for an alternative approach for such problems.In this paper, an initial value method for solving a weakly coupled system of two second-order singularly perturbed convection–diffusion equations exhibiting a boundary layer at one end is proposed. The technique, used in this work, is the careful factorization of original problem into a system of IVPs and two explicit IVPs which are independent of perturbation parameter. First, a system of IVPs is obtained by putting the perturbation parameter to zero, namely the reduced system, which corresponds to the outer solution. Next, using reduction of order together with stretching of variable gives two decoupled IVPs, namely the boundary layer correction problems, which corresponds to the inner solution. And then, the reduced system is solved numerically using fourth-order Runge–Kutta method, whereas, the boundary layer correction problems, are solved analytically. Finally, combining the above two solutions we obtain an approximation for the original problem. In addition, error estimates are derived and examples are provided to illustrate the method.


## 2. Statement of the Problem

Consider the problem of findingy1,y2∈C0Ω¯∩C2Ω such that(1a)L1y¯≡−εy1″x−a1xy1′x+b11xy1x+b12xy2x=f1x,∀x∈Ω,(1b)L2y¯≡−εy2″x−a2xy2′x+b21xy1x+b22xy2x=f2x,∀x∈Ω,with the boundary conditions(2a)y10=p1,y11=q1,(2b)y20=p2,y21=q2,whereΩ=0,1, Ω¯=0,1 and y¯=y1,y2T, p1, p2, q1 and q2 are given constants, and 0<ε≪1 is the singular perturbation parameter. The coefficient functions are taken to be sufficiently smooth on Ω¯ and satisfying the following conditions:(3a)a1x≥α1>0,a2x≥α2>0,(3b)b11x+b12x≥β1>0,b22x+b21x≥β2>0,∀x∈Ω¯,(3c)b12x≤0,b21x≤0.Under these assumptions the system (1a)–(2b) has a unique solution y¯x which exhibits a boundary layer of width Oε on the left side x=0 of the underlying interval. The case aix≤αi<0 for i=1,2, can be put in to (1a) and (1b)by the change of independent variable from x to 1−x.The above coupled system of Equations (1a)–(2b) can also be written in vector form as(3d)Ly¯x=−εd2dx200−εd2dx2y¯x−Axy′¯x+Bxy¯x=f¯x,∀x∈Ω,with the boundary conditions(3e)y¯0=p1p2,y¯1=q1q2,where(3f)y¯x=y1xy2x,Ax=a1x00a2x,Bx=b11xb12xb21xb22xf¯x=f1xf2x.Remark 1.
In this paper, we consider only the case where there is one boundary layer at the left end of the interval. The case when the layer occurs at right end, can be analyzed similarly.Notations. Let y:D→R, the appropriate norm for studying the convergence of the approximate solution to the exact solution is the maximum norm yD=maxx∈Dyx. In case of vectors y¯=y1,y2T, we define (3g)y¯=y1,y2Ty¯D=maxx∈Dy1D,y2D.Throughout this paper,C (sometimes sub-scripted) denotes generic positive constants independent of the singular perturbation parameter ε and in the case of discrete problems, also independent of the mesh parameter N, these constants may assume different values but remains to be constant.


## 3. Analytic Results

In this section, a maximum principle, a stability result, and estimates of the derivatives of the system of Equations (1a)–(2b) are presented. First, we consider the following property of the operators L1 and L2.Lemma 2 (Maximum principle).
Assume thatπ¯x is any sufficiently smooth function such that π¯0≥0¯, π¯1≥0¯ and L1π¯x≥0¯, L2π¯x≥0¯ in Ω, then π¯x≥0¯ in Ω¯.Proof.
Letx∗ and y∗ be arbitrary points in 0,1 such that π1x∗=minx∈Ω¯π1x and π2y∗=minx∈Ω¯π2x. Without loss of generality, we assume that π1x∗≤π2y∗ and suppose π1x∗<0. Clearly x∗∉0,1 and π1′x∗=0, π1″x∗≥0. Moreover,(4)L1π¯x∗=−επ1″x∗−a1x∗π1′x∗+b11x∗π1x∗+b12x∗π2x∗=−επ1″x∗+b11x∗+b12x∗π1x∗+π2x∗−π1x∗b12x∗<0,
which is a contradiction. It follows that our assumptionπ1x∗<0 is wrong, so that π1x∗≥0. Similarly, L2π¯x can be dealt. Hence, π¯x≥0¯ in Ω¯.An immediate consequence of the maximum principle is the following stability result.Lemma 3 (Stability result).
ify1,y2∈C0Ω¯∩C2Ω, then for i=1,2,(5)yix≤maxy¯0,y¯1+1βf¯,∀x∈Ω¯,
whereβ=minβ1,β2.Proof.
Define two barrier functionsφ¯±=φ1±,φ2±T as(6)
φ
¯
±
x
=
C
e
¯
±
y
¯
x
,wheree¯=1,1T is the unit column vector and C=max‖y¯0‖,y¯1+1/βf¯.Clearly,φ¯±0≥0¯ and φ¯±1≥1¯.Also(7)L1φ¯±x=−εφ1±x″−a1xφ1±x′+b11xφ1±x+b12xφ2±x=Cb11x+b12x±L1y¯x≥β1C±f1x≥βmaxy¯0,y¯1+f1±f1x≥0.Similarly, it can be proved thatL1φ¯±x≥0. Therefore, by the maximum principle, we obtain φ¯±x≥0¯ for all x∈Ω¯, which gives the required estimate. □Now we give bounds on the derivatives of the exact solutiony¯x for system (1a)–(2b).Lemma 4.
Lety¯x be the solution of (1a)–(2b), then ∀x∈Ω and for k=1,2,(8)
y
1
k
x
≤
C
1
+
ε
−
k
exp
−
α
1
x
ε
,
y
2
k
x
≤
C
1
+
ε
−
k
exp
−
α
2
x
ε
.Proof.
By following the approach used in the proof of Theorem 2 of [19] and the technique from [23], we can easily derive this Lemma. □The solutiony¯x of the problem (1a)–(2b) can be decomposed into smooth and singular components v¯ and w¯ respectively, as(9)y¯x=v¯x+w¯x,wherev¯=v1,v2T and w¯=w1,w2T. Further, the regular component can be written in the form of v¯=u¯0+εu¯1+ε2u¯2, where u¯0=u01,u02T is the solution of the following system(10)−a1xu01′x+b11xu01x+b12xu02x=f1x,−a2xu02x+b21xu01x+b22xu02x=f2x,u011=q1andu021=q2.andu¯1=u11,u12T is the solution of the following system(11)−a1xu11′x+b11xu11x+b12xu12x=u01″x,−a2xu12′x+b21xu11x+b22xu12x=u02″x,u111=0andu121=q2,andu¯2=u21,u22T is the solution of the following system(12)−εu21″x−a1xu21′x+b11xu21x+b12xu22x=u11″x,−εu22″x−a2xu22′x+b21xu21x+b22xu22x=u12″x,u210=r1,u220=r2,u211=0andu221=0,,where r1 and r2 are constants to be chosen such that ri≤C, for i=1,2. Thus the regular component v¯ is the solution of (13)Lv¯x=f¯x,v¯1=q¯,v¯0suitablychosen.and the singular componentw¯ is the solution of (14)Lw¯x=0¯,w¯0=p¯−v¯0,w¯1=0¯.The following lemma provides the bound on the derivatives of the regular and singular components of the solutiony¯.Lemma 5.
For0≤k≤2 the smooth component v¯ and the singular component w¯ and their derivatives satisfy the bounds:(15)νikx≤C1+ε2−kexp−αixε,∀x∈Ω¯,fori=1,2,(16)wikx≤Cε−kexp−αixε,∀x∈Ω¯fori=1,2,
respectively.Proof.
Using appropriate barrier functions, applying Lemma2 and adopting the method of proof used in [[5] page 44], the present Lemma can be proved.


## 4. Description of the Method

In this section, we will obtain the solution of the (1a)–(2b) as a combination of two solutions: outer solution and inner solution.


### 4.1. Outer Solution

Letu¯0x=u01,u02T, u0i∈C1Ω∩C0Ω¯ for i=1,2, be the solution of the reduced problem of (1a)–(2b) given by(17)−a1xu′01x+b11xu01x+b12xu02x=f1x,−a2xu′02x+b21xu01x+b22xu02x=f2x,u011=q1u021=q2.Since, the degenerate equation does not satisfy the condition atx=0, therefore, its contribution to the solution of (1a)–(2b) is for those values of x which are away from x=0. The problem (17) is therefore termed as outer problem.For the exact solution of the reduced problem the following theorem gives an error bound.Theorem 6.
Lety¯x be the solution of (1a)–(2b) and u¯0x be its reduced problem solution defined by (17). Then(18)yix−u0ix≤Cε+exp−αxε,∀x∈Ω¯fori=1,2,
whereα=minα1,α2.Proof.
Consider the following barrier functionφ¯±=φ1±,φ2±T, where(19)
φ
i
±
x
=
C
1
ε
1
−
x
2
+
C
1
exp
−
α
x
ε
±
y
i
x
−
u
0
i
x
,
x
∈
Ω
¯
,
i
=
1
,
2
.Its easy to see thatφi±∈C1Ω∩C0Ω¯ further,(20)φi±0=C1ε+C1±yi0−u0i0≥C1±yi0−u0i0≥0,for an appropriate choice ofC1, and(21)φi±1=C1ε12+C1exp−αε±yi1−u0i1=C112ε+exp−αε±qi−qi≥0,next, for the operatorL1 we have(22)L1φ¯±x=−εφ1±x′′−a1xφ1±x′+b11xφ1±x+b12xφ2±x=C1ε12a1+b11+b121−x2+C1α2ε+a1xαε+b11+b12exp−αxε±L1y¯x−u¯0x≥C1ε12α1+12β1+C1αεa1x−α+β1exp−αxε±−εu01″x≥C1ε12α1+β1+C1αεα1−α+β1exp−αxε∓εu01″x≥C2ε∓εu0″x=εC2∓u01″x≥0,sinceu01″x≤C, for an appropriate choice of C2. Similarly, we can prove that L2φ¯±x≥0, for all x∈Ω¯. Therefore, from the maximum principle of Lemma 2, we obtain φi±x≥0, ∀x∈Ω¯ and for i=1,2. Hence the proof of the theorem.Remark 7.
From the above theorem, it is clear that the solutiony¯ of problem (1a)–(2b) exhibits a strong boundary layer at x=0 and further away from the boundary layer region and in particular on kε,1, where k≥−lnε/α, for sufficiently small values of ε, we have.(23)
y
i
x
−
u
0
i
x
≤
C
ε
,
∀
x
∈
Ω
¯
and
for
i
=
1
,
2
.For the numerical solution of the reduced problem (17) we employ fourth-order Runge–Kutta method for system. Suppose U¯0=U01,U02 be the numerical solution of the reduced problem obtained from fourth-order Runge–Kutta method, then the maximum error becomes(24)u0ixj−U0ixj≤Ch4,forj=0,1,…,Nandfori=1,2,whereh=1/N is the equal mesh spacing of the domain of the problem.


### 4.2. Inner Solution

To obtain the inner solution for (1a)–(2b) we will use reduction of order together with streching of variable as follows:First we rewrite the given problem (1a) for i=1,2 equivalently as(25)−εd2yidx2−ddxaixyix=Fix,y1,y2,x∈Ω,where(26)Fix,y1,y2=fix−ai′xyix−∑j=12bijxyjx.From Theorem6 that the solution u¯0x satisfies (1a)–(2b) on most part of the interval 0,1 and away from x=0. Thus by replacing the solution y¯x by u¯0x on the right part of (25), we obtain an asymptotically equivalent approximation as:(27)−εd2yidx2−ddxaixyix=Fix,u01,u02+Oε,x∈Ω,where(28)Fix,u01,u02=fix−ai′xu0ix−∑j=12bijxu0jx.Integrating both sides of (27) with respect to x, gives(29)−εdyidx−aixux=∫Fix,u01,u02dx+Oε,x∈0,1,where(30)∫Fix,u01,u02dx=∫fix−ai′xu0ix−∑j=12bijxu0jxdx.Using the reduced problem (16) in the above integral, yields(31)∫Fix,u01,u02dx=∫fix−ai′xu0ix−∑j=12bijxu0jxdx=∫−ai′xu0ix−aixu0i′xdx=−∫ddxaixu0ixdx=−aixu0ix+ki.Then, substituting this in to (28), gives us(32)−εdyidx−aixyix=−aixu0ix+ki+Oε,wherek1 and k2 are integration constants. In order to determine ki’s, we introduce the condition that the reduced equations of (32) should satisfy the boundary condition at x=1. Thus, we get k1=k2=0.Hence, by substitutingk1=k2=0 in (32), a first-order initial value problem which is asymptotically equivalent to the second-order system of boundary value problems (1a)–(2b) is obtained, and written as:(33)εdwidx+aixwix=aixu0ix,fori=1,2,withaninitialconditions,w10=p1,w20=p2.Next, to compute the solution for the layer part, a new inner variable is introduced by stretching the spatial coordinatex, as(34)t=xε⇒x=εtanddtdx=1ε.Using this stretching transformation in to (33), we obtain(35)dwidt+aiεtwi=aiεtu0iεt.In spite of the simplification, these equations remains first-order differential equation and also regularly perturbed. Thus, forε=0, (24) becomes(36)dwidt+ai0wi=ai0u0i0.These are differential equation for the solution of the layer regions. Moreover, the solutions of (36) are supposed to counter act for the fact that the solutions of the reduced problem do not satisfy the boundary condition at x=0.Further, using the substitutions,Wit=wit−u0i0 in to (36), we obtain the following boundary layer correction problems(37)dWidt+ai0Wi=0,withWi0=pi−u0i0,fori=1,2,Since these equations are separable linear initial value problems with constant coefficients which can easily be solved analytically, thus and gives(38)Wixε=pi−u0i0exp−ai0xε,fori=1,2.Finally, from standard singular perturbation theory it follows that the solution of the (33) admits the representation in terms of the solutions of the reduced problem (16) and boundary layer correction problem (38), which approximates the solution of the system (1a)–(2b); that is,(39)yix=u0ix+wix+Oε=u0ix+Wixε+Oε,yix=u0ix+pi−u0i0exp−ai0xε+Oε,∀x∈Ω¯andfori=1,2.The numerical error of the present method has two sources: one from the asymptotic approximation of the modified problem (33) and the other from the numerical approximation of the reduced system (16). We can summarize the results of this section in the form of the following theorem.Theorem 8.
Lety¯x be the solution of the problem (1a)–(2b) given by (27) and Y¯x=Y1x,Y2xT be its approximate then
(40)
y
¯
x
−
Y
¯
x
=
O
ε
+
h
4
,
∀
x
∈
Ω
¯
,
whereh is the equal mesh spacing for the domain of the problem.Proof.
Assume thatU¯0 be the numerical solution of the outer problem determined by making use of the Runge–Kutta method and W¯ be the solution of the inner problem whose left boundary conditions are affected by U¯0, such that Y¯=U¯0+W¯ is the approximation for the exact solution of (1a)–(2b) given by y¯=u¯0+w¯+C1ε. Now using (17) and (27), we obtain,(41)
y
¯
x
−
Y
¯
x
=
u
¯
0
+
w
¯
+
C
1
ε
−
U
¯
0
+
W
¯
=
u
¯
0
−
U
¯
0
+
p
¯
−
u
¯
0
0
−
p
¯
+
U
¯
0
0
exp
−
a
¯
0
x
/
ε
+
C
1
ε
=
u
¯
0
−
U
¯
0
−
u
¯
0
0
−
U
¯
0
0
exp
−
a
¯
0
x
/
ε
+
C
1
ε
≤
C
2
h
4
1
−
exp
−
a
¯
0
x
/
ε
+
C
1
ε
≤
C
ε
+
h
4
,for an appropriate choice ofC.Therefore, we can conclude thaty¯x−Y¯x=Oε+h4, ∀x∈Ω¯. □


## 4.1. Outer Solution

Letu¯0x=u01,u02T, u0i∈C1Ω∩C0Ω¯ for i=1,2, be the solution of the reduced problem of (1a)–(2b) given by(17)−a1xu′01x+b11xu01x+b12xu02x=f1x,−a2xu′02x+b21xu01x+b22xu02x=f2x,u011=q1u021=q2.Since, the degenerate equation does not satisfy the condition atx=0, therefore, its contribution to the solution of (1a)–(2b) is for those values of x which are away from x=0. The problem (17) is therefore termed as outer problem.For the exact solution of the reduced problem the following theorem gives an error bound.Theorem 6.
Lety¯x be the solution of (1a)–(2b) and u¯0x be its reduced problem solution defined by (17). Then(18)yix−u0ix≤Cε+exp−αxε,∀x∈Ω¯fori=1,2,
whereα=minα1,α2.Proof.
Consider the following barrier functionφ¯±=φ1±,φ2±T, where(19)
φ
i
±
x
=
C
1
ε
1
−
x
2
+
C
1
exp
−
α
x
ε
±
y
i
x
−
u
0
i
x
,
x
∈
Ω
¯
,
i
=
1
,
2
.Its easy to see thatφi±∈C1Ω∩C0Ω¯ further,(20)φi±0=C1ε+C1±yi0−u0i0≥C1±yi0−u0i0≥0,for an appropriate choice ofC1, and(21)φi±1=C1ε12+C1exp−αε±yi1−u0i1=C112ε+exp−αε±qi−qi≥0,next, for the operatorL1 we have(22)L1φ¯±x=−εφ1±x′′−a1xφ1±x′+b11xφ1±x+b12xφ2±x=C1ε12a1+b11+b121−x2+C1α2ε+a1xαε+b11+b12exp−αxε±L1y¯x−u¯0x≥C1ε12α1+12β1+C1αεa1x−α+β1exp−αxε±−εu01″x≥C1ε12α1+β1+C1αεα1−α+β1exp−αxε∓εu01″x≥C2ε∓εu0″x=εC2∓u01″x≥0,sinceu01″x≤C, for an appropriate choice of C2. Similarly, we can prove that L2φ¯±x≥0, for all x∈Ω¯. Therefore, from the maximum principle of Lemma 2, we obtain φi±x≥0, ∀x∈Ω¯ and for i=1,2. Hence the proof of the theorem.Remark 7.
From the above theorem, it is clear that the solutiony¯ of problem (1a)–(2b) exhibits a strong boundary layer at x=0 and further away from the boundary layer region and in particular on kε,1, where k≥−lnε/α, for sufficiently small values of ε, we have.(23)
y
i
x
−
u
0
i
x
≤
C
ε
,
∀
x
∈
Ω
¯
and
for
i
=
1
,
2
.For the numerical solution of the reduced problem (17) we employ fourth-order Runge–Kutta method for system. Suppose U¯0=U01,U02 be the numerical solution of the reduced problem obtained from fourth-order Runge–Kutta method, then the maximum error becomes(24)u0ixj−U0ixj≤Ch4,forj=0,1,…,Nandfori=1,2,whereh=1/N is the equal mesh spacing of the domain of the problem.


## 4.2. Inner Solution

To obtain the inner solution for (1a)–(2b) we will use reduction of order together with streching of variable as follows:First we rewrite the given problem (1a) for i=1,2 equivalently as(25)−εd2yidx2−ddxaixyix=Fix,y1,y2,x∈Ω,where(26)Fix,y1,y2=fix−ai′xyix−∑j=12bijxyjx.From Theorem6 that the solution u¯0x satisfies (1a)–(2b) on most part of the interval 0,1 and away from x=0. Thus by replacing the solution y¯x by u¯0x on the right part of (25), we obtain an asymptotically equivalent approximation as:(27)−εd2yidx2−ddxaixyix=Fix,u01,u02+Oε,x∈Ω,where(28)Fix,u01,u02=fix−ai′xu0ix−∑j=12bijxu0jx.Integrating both sides of (27) with respect to x, gives(29)−εdyidx−aixux=∫Fix,u01,u02dx+Oε,x∈0,1,where(30)∫Fix,u01,u02dx=∫fix−ai′xu0ix−∑j=12bijxu0jxdx.Using the reduced problem (16) in the above integral, yields(31)∫Fix,u01,u02dx=∫fix−ai′xu0ix−∑j=12bijxu0jxdx=∫−ai′xu0ix−aixu0i′xdx=−∫ddxaixu0ixdx=−aixu0ix+ki.Then, substituting this in to (28), gives us(32)−εdyidx−aixyix=−aixu0ix+ki+Oε,wherek1 and k2 are integration constants. In order to determine ki’s, we introduce the condition that the reduced equations of (32) should satisfy the boundary condition at x=1. Thus, we get k1=k2=0.Hence, by substitutingk1=k2=0 in (32), a first-order initial value problem which is asymptotically equivalent to the second-order system of boundary value problems (1a)–(2b) is obtained, and written as:(33)εdwidx+aixwix=aixu0ix,fori=1,2,withaninitialconditions,w10=p1,w20=p2.Next, to compute the solution for the layer part, a new inner variable is introduced by stretching the spatial coordinatex, as(34)t=xε⇒x=εtanddtdx=1ε.Using this stretching transformation in to (33), we obtain(35)dwidt+aiεtwi=aiεtu0iεt.In spite of the simplification, these equations remains first-order differential equation and also regularly perturbed. Thus, forε=0, (24) becomes(36)dwidt+ai0wi=ai0u0i0.These are differential equation for the solution of the layer regions. Moreover, the solutions of (36) are supposed to counter act for the fact that the solutions of the reduced problem do not satisfy the boundary condition at x=0.Further, using the substitutions,Wit=wit−u0i0 in to (36), we obtain the following boundary layer correction problems(37)dWidt+ai0Wi=0,withWi0=pi−u0i0,fori=1,2,Since these equations are separable linear initial value problems with constant coefficients which can easily be solved analytically, thus and gives(38)Wixε=pi−u0i0exp−ai0xε,fori=1,2.Finally, from standard singular perturbation theory it follows that the solution of the (33) admits the representation in terms of the solutions of the reduced problem (16) and boundary layer correction problem (38), which approximates the solution of the system (1a)–(2b); that is,(39)yix=u0ix+wix+Oε=u0ix+Wixε+Oε,yix=u0ix+pi−u0i0exp−ai0xε+Oε,∀x∈Ω¯andfori=1,2.The numerical error of the present method has two sources: one from the asymptotic approximation of the modified problem (33) and the other from the numerical approximation of the reduced system (16). We can summarize the results of this section in the form of the following theorem.Theorem 8.
Lety¯x be the solution of the problem (1a)–(2b) given by (27) and Y¯x=Y1x,Y2xT be its approximate then
(40)
y
¯
x
−
Y
¯
x
=
O
ε
+
h
4
,
∀
x
∈
Ω
¯
,
whereh is the equal mesh spacing for the domain of the problem.Proof.
Assume thatU¯0 be the numerical solution of the outer problem determined by making use of the Runge–Kutta method and W¯ be the solution of the inner problem whose left boundary conditions are affected by U¯0, such that Y¯=U¯0+W¯ is the approximation for the exact solution of (1a)–(2b) given by y¯=u¯0+w¯+C1ε. Now using (17) and (27), we obtain,(41)
y
¯
x
−
Y
¯
x
=
u
¯
0
+
w
¯
+
C
1
ε
−
U
¯
0
+
W
¯
=
u
¯
0
−
U
¯
0
+
p
¯
−
u
¯
0
0
−
p
¯
+
U
¯
0
0
exp
−
a
¯
0
x
/
ε
+
C
1
ε
=
u
¯
0
−
U
¯
0
−
u
¯
0
0
−
U
¯
0
0
exp
−
a
¯
0
x
/
ε
+
C
1
ε
≤
C
2
h
4
1
−
exp
−
a
¯
0
x
/
ε
+
C
1
ε
≤
C
ε
+
h
4
,for an appropriate choice ofC.Therefore, we can conclude thaty¯x−Y¯x=Oε+h4, ∀x∈Ω¯. □


## 5. Test Problems and Numerical Results

To demonstrate the efficiency and applicability of the proposed method we considered the following two test problems:Example 5.1.
Consider the following boundary value problems for the systems of convention–diffusion equations on0,1:(42)
−
ε
y
1
″
x
−
y
1
′
x
+
4
y
1
x
−
y
2
x
=
4
,
−
ε
y
2
″
x
−
y
2
′
x
−
y
1
x
+
4
y
2
x
=
2
,
y
1
0
=
1
,
y
1
1
=
0
,
y
2
0
=
1
,
y
2
1
=
0
.The exact solution of this problem is(43)y1x=65−1expm1−expm2expm1x+1expm1−expm2expm2x−151−expm4expm3−expm4expm3x−15expm3−1expm3−expm4expm4x(44)y2x=45−1expm1−expm2expm1x+1expm1−expm2expm2x+151−expm4expm3−expm4expm3x+15expm3−1expm3−expm4expm4x,where,m
1
=
−
1
+
1
+
12
ε
/
2
ε, m2=−1−1+12ε/2ε, m3=−1+1+20ε/2ε and m4=−1−1+20ε/2ε.Since the exact solution is known, we calculate maximum point wise error by(45)Eε,iN=max0≤j≤Nyixj−Yi,jN,fori=1,2,whereyixj is the exact solution and Yi,jN is the numerical solution obtained by using N mesh intervals. Tables 1 and 2 display, respectively, the maximum point-wise errors for y1 and y2 for different values of N and ε. The plots of the exact and the approximate solution components for N=1024 and ε=0.01 are shown in Figure 1.Table 1
Maximum point wise errorEε,1N, for the solution y1 of Example 9.
ε
↓
N
=
64
N
=
128
N
=
256
N
=
512
N
=
1024
1
0
−
1
1.001
E
−
01
1.001
E
−
01
1.001
E
−
01
1.001
E
−
01
1.001
E
−
01
1
0
−
2
1.366
E
−
02
1.366
E
−
02
1.366
E
−
02
1.366
E
−
02
1.366
E
−
02
1
0
−
3
1.426
E
−
03
1.426
E
−
03
1.426
E
−
03
1.426
E
−
03
1.426
E
−
03
1
0
−
4
1.432
E
−
04
1.433
E
−
04
1.433
E
−
04
1.433
E
−
04
1.433
E
−
04
1
0
−
5
1.429
E
−
05
1.433
E
−
05
1.434
E
−
05
1.434
E
−
05
1.434
E
−
05
1
0
−
6
1.395
E
−
06
1.431
E
−
06
1.433
E
−
06
1.434
E
−
06
1.434
E
−
06
1
0
−
7
1.358
E
−
07
1.410
E
−
07
1.432
E
−
07
1.433
E
−
07
1.410
E
−
07
1
0
−
8
2.393
E
−
08
1.280
E
−
08
1.495
E
−
08
1.508
E
−
08
1.509
E
−
08Table 2
Maximum point wise errorEε,2N, for the solution y2 of Example 9.
ε
↓
N
=
64
N
=
128
N
=
256
N
=
512
N
=
1024
1
0
−
1
6.133
E
−
02
6.133
E
−
02
6.133
E
−
02
6.133
E
−
02
6.133
E
−
02
1
0
−
2
7.848
E
−
03
7.848
E
−
03
7.848
E
−
03
7.848
E
−
03
7.848
E
−
03
1
0
−
3
8.119
E
−
04
8.119
E
−
04
8.119
E
−
04
8.119
E
−
04
8.119
E
−
04
1
0
−
4
8.148
E
−
05
8.147
E
−
05
8.147
E
−
05
8.147
E
−
05
8.147
E
−
05
1
0
−
5
8.152
E
−
06
8.150
E
−
06
8.150
E
−
06
8.150
E
−
06
8.150
E
−
06
1
0
−
6
8.171
E
−
07
8.152
E
−
07
8.151
E
−
07
8.151
E
−
07
8.151
E
−
07
1
0
−
7
8.462
E
−
08
8.159
E
−
08
8.149
E
−
08
8.148
E
−
08
8.148
E
−
08
1
0
−
8
1.746
E
−
08
8.996
E
−
09
8.842
E
−
09
8.835
E
−
09
8.835
E
−
09Figure 1
Plot of the exact and approximate solutions of Example9 for ε=0.01 and N=1024.Example 10.
Consider the following boundry value problems for the systems of convention–diffusion equations on(0,1):0,1:(46)
−
ε
y
1
″
x
−
1
+
x
2
y
1
′
x
+
4
+
s
i
n
x
y
1
x
−
2
y
2
x
=
e
x
p
x
,
−
ε
y
2
″
x
−
2
+
x
y
2
′
x
−
y
1
x
+
2
+
c
o
s
x
y
2
x
=
x
2
,
w
i
t
h
y
1
0
=
3
,
y
1
1
=
1
,
y
2
0
=
3
,
y
2
1
=
1
.Since the exact solution is not known, we calculate maximum point wise error by using the double mesh principle defined by(47)Dε,iN=max0≤j≤NYi,jN−Yi,j2N,fori=1,2,whereYi,jN and Yi,j2N denotes the ith and 2ith components of the numerical solutions obtained by using N and 2N meshes points, respectively. Tables 3 and 4 display, respectively, the maximum point-wise errors for Y1 and Y2 for different values of N and ε. Figure 2 represents the numerical solutions of Example 10, for N=1024 and ε=0.05.Table 3
Maximum point wise errorDε,1N, for the solution Y1 of Example 9.
ε
↓
N
=
8
N
=
16
N
=
32
N
=
64
N
=
128
2
−
1
8.339
E
−
06
4.267
E
−
07
2.409
E
−
08
1.433
E
−
9
8.734
E
−
11
2
−
2
7.914
E
−
06
4.018
E
−
07
2.260
E
−
08
1.343
E
−
9
8.181
E
−
11
2
−
4
7.743
E
−
06
3.919
E
−
07
2.207
E
−
08
1.310
E
−
9
7.975
E
−
11
2
−
6
7.743
E
−
06
3.918
E
−
07
2.207
E
−
08
1.310
E
−
9
7.975
E
−
11
2
−
8
7.743
E
−
06
3.918
E
−
07
2.207
E
−
08
1.310
E
−
9
7.975
E
−
11
2
−
12
7.743
E
−
06
3.918
E
−
07
2.207
E
−
08
1.310
E
−
9
7.975
E
−
11
2
−
16
7.743
E
−
06
3.918
E
−
07
2.207
E
−
08
1.310
E
−
9
7.975
E
−
11
2
−
20
7.743
E
−
06
3.918
E
−
07
2.207
E
−
08
1.310
E
−
9
7.975
E
−
11
2
−
24
7.743
E
−
06
3.918
E
−
07
2.207
E
−
08
1.310
E
−
9
7.975
E
−
11Table 4
Maximum point wise errorDε,2N, for the solution Y2 of Example 9.
ε
↓
N
=
8
N
=
16
N
=
32
N
=
64
N
=
128
2
−
1
9.740
E
−
07
4.915
E
−
08
2.771
E
−
09
1.643
E
−
10
9.995
E
−
12
2
−
2
9.829
E
−
07
4.973
E
−
08
2.824
E
−
09
1.678
E
−
10
1.022
E
−
11
2
−
4
9.833
E
−
07
4.977
E
−
08
2.828
E
−
09
1.680
E
−
10
1.024
E
−
11
2
−
6
9.833
E
−
07
4.977
E
−
08
2.828
E
−
09
1.680
E
−
10
1.024
E
−
11
2
−
8
9.833
E
−
07
4.977
E
−
08
2.828
E
−
09
1.680
E
−
10
1.024
E
−
11
2
−
12
9.833
E
−
07
4.977
E
−
08
2.828
E
−
09
1.680
E
−
10
1.024
E
−
11
2
−
16
9.833
E
−
07
4.977
E
−
08
2.828
E
−
09
1.680
E
−
10
1.024
E
−
11
2
−
20
9.833
E
−
07
4.977
E
−
08
2.828
E
−
09
1.680
E
−
10
1.024
E
−
11
2
−
24
9.833
E
−
07
4.977
E
−
08
2.828
E
−
09
1.680
E
−
10
1.024
E
−
11Figure 2
Plot of the solutions of Example10 for ε=0.05 and N=1024.


## 6. Discussion

In this paper, an initial value method for solving a weakly coupled system of two linear second-order singularly perturbed convection–diffusion equations exhibiting a boundary layer at one end is proposed. The method is some how similar to the asymptotic expansion methods, but differs in detail. The approximate solution of the given problem is obtained by solving a coupled system of initial value problem (namely, the reduced system) and two decoupled initial value problems with constant coefficients (namely, the layer correction problems), which are easily deduced from the the original problem. Both the reduced system and the layer correction problems are independent of perturbation parameter,ε and therefore, we get easily the numerical solution by solving the reduced system using fourth-order Runge–Kutta method and solving the layer correction problems analytically. The method is simple to apply, very easy to implement on a computer and offers a relatively simple tool for obtaining approximate solution.We have implemented the method on two test problems to illustrate the theoretical results, and presented the computational results for different values ofε and N in Tables 1–4 and Figures 1 and 2. From the results it is observed that, for very small ε the present method approximates the exact solution of the problems very well.



---

*Source: 1062025-2020-02-07.xml*