We consider a weakly coupled system of singularly perturbed convectiondiffusion equations with discontinuous source term. The diffusion term of each equation is associated with a small positive parameter of different magnitude. Presence of discontinuity and different parameters creates boundary and weak interior layers that overlap and interact. A numerical method is constructed for this problem which involves an appropriate piecewise uniform Shishkin mesh. The numerical approximations are proved to converge to the continuous solutions uniformly with respect to the singular perturbation parameters. Numerical results are presented which illustrates the theoretical results.
AMS Mathematics Subject Classification : 65L10, 65L11, 65L12.
1. Introduction
An extensive research had been done on numerical methods for a single singularly perturbed convectiondiffusion differential equation
[1]

[4]
, but for system of equations very few works had been done. The classical numerical methods fail to produce good approximations for singularly perturbed problems. Various nonclassical approaches produce better approximations and converge uniformly with respect to the small perturbation parameter. In the literature
[7]

[14]
methods were available to obtain numerical approximation for system of singularly perturbed convection diffusion differential equations the source term are smooth on the whole domain. Farrell et.al
[5]

[6]
considered scalar singularly perturbed convectiondiffusion equation with discontinuous source term. The interior layers in
[6]
were strong, in the sense that the solution was bounded but the magnitude of the first derivative grew unboundedly as
ε
→ 0, but in
[5]
they were weak, in the sense that the solution and the first derivative were bounded but the magnitude of the second derivative grows unboundedly as
ε
→ 0. In this work, we present a uniformly convergent numerical method for a weakly coupled system of singularly perturbed convectiondiffusion equations having discontinuous source term with different diffusion parameters. The solution to such equations has overlapping and interacting boundary and interior layers which makes the construction of numerical methods and analysis quite difficult. Tamilselvan and Ramanujam
[15]
considered the same problem but with equal diffusion parameters.
Consider a weakly coupled system of singularly perturbed convectiondiffusion equations with discontinuous source term on the unit interval Ω = (0, 1), having a single discontinuity in the source term at a point
d
∈ Ω. Let Ω
_{1}
= (0,
d
) and Ω
_{2}
= (
d
, 1). Let the jump in a function
ω
at a point
d
∈ Ω given as [
ω
](
d
) =
ω
(
d
+) −
ω
(
d
−). The corresponding boundary value problem is to find
u
_{1}
,
u
_{2}
∈
C
^{0}
(
) ∩
C
^{1}
(Ω) ∩
C
^{2}
(Ω
_{1}
∪ Ω
_{2}
), such that
where
E
= diag(
ε
_{1}
,
ε
_{2}
), the coupling matrix
A
= diag(
a
_{1}
,
a
_{2}
) and
B
= (
b_{ij}
)
_{2×2}
with 0 <
ε
_{1}
≤
ε
_{2}
≤ 1,
f
= (
f
_{1}
,
f
_{2}
)
^{T}
, and
u
= (
u
_{1}
,
u
_{2}
)
^{T}
. Assume for each
i
= 1, 2 and
x
∈
, the matrices
A
and
B
satisfy
Let
α
= min{
α
_{1}
,
α
_{2}
}. Further assume that the source terms
f
_{1}
,
f
_{2}
are sufficiently smooth on
and their derivatives have jump discontinuity at the same point.
Notations.
Throughout the paper, C denotes a generic positive constant and
C
= (
C
,
C
)
^{T}
denotes a generic positive constant vector, both are independent of perturbation parameters
ε
_{1}
,
ε
_{2}
and the discretization parameter N, but may not be same at each occurrence. Define
v
≤
w
if
v_{i}
≤
w_{i}
,
i
= 1, 2, and 
v
 = (
v
_{1}
, 
v
_{2}
)
^{T}
. We consider the maximum norm and denote it by║.║
_{S}
, where S is a closed and bounded subset in
. For a real valued function
v
∈
C
(
S
) and for a vector valued function
v
= (
v
_{1}
,
v
_{2}
)
^{T}
∈
C
(
S
)
^{2}
, we define
and ║
v
║
_{S}
= max{║
v
_{1}
║
_{S}
,║
v
_{2}
║
_{S}
}. Now let a mesh
be a set of points satisfying
x
_{0}
<
x
_{1}
<⋯<
x_{N}
= 1. A mesh function
is a realvalued function defined on Ω
^{N}
. Define the discrete maximum norm for such functions by
and for vector mesh functions
V
= (
V
_{1}
,
V
_{2}
)
^{T}
=
are used and define║
V
║
_{ΩN}
= max{║
V
_{1}
║
_{ΩN}
,║
V
_{2}
║
_{ΩN}
}.
 2. The continuous problem
Theorem 2.1
(Continuous maximum principle).
Suppose u
_{1}
,
u
_{2}
∈
C
^{0}
(
) ∩
C
^{1}
(Ω) ∩
C
^{2}
(Ω
_{1}
∪ Ω
_{2}
).
Further suppose that u
= (
u
_{1}
,
u
_{2}
)
^{T} satisfies u
(0) ≥
0, u
(1) ≥
0, Lu
(
x
) ≥
0 in
Ω
_{1}
∪ Ω
_{2}
and
[
u′
](
d
) ≤
0. Then u
(
x
) ≥
0, for all x
∈
.
Proof
. Let
u
:=
Now let
p
and
q
be the points at which
θ
_{1}
(
p
) :=
and
θ
_{2}
(
q
) :=
. Assume without loss of generality
θ
_{1}
(
p
) ≤
θ
_{2}
(
q
). If
θ
_{1}
(
p
) ≥ 0, then there is nothing to prove. Suppose that
θ
_{1}
(
p
) < 0, then proof is completed by showing that this leads to a contradiction. Note that
p
≠ {0, 1}. So either
p
∈ Ω
_{1}
∪ Ω
_{2}
or
p
=
d
.
In the first case for
x
∈ Ω
_{2}
,
In the second case, that is,
p
=
d
, we have [
u′
](
d
) =
, and at a negative minimum [Θ′](
d
) ≥ 0, which gives a contradiction. □
Lemma 2.2
(Stability Result).
Let u
= (
u
_{1}
,
u
_{2}
)
^{T}
be the solution of
(1) − (2),
then,
Proof
. Define the function Ψ
^{±}
(
x
) := max {║
u
(0)║,║
u
(1)║,║
Lu
║Ω
_{1}
∪Ω
_{2}
} (2 −
x
, 2 −
x
)
^{T}
±
u
(
x
). Then Ψ
^{±}
(0) ≥
0
, Ψ
^{±}
(1) ≥
0
,
L
Ψ
^{±}
(
x
) ≥
0
for each
x
∈ Ω
_{1}
∪ Ω
_{2}
, and [Ψ
^{±}
′ ](
d
) = ±[
u′
](
d
) =
0
, since
u
∈
C
^{1}
(Ω)
^{2}
. It follows from the maximum principle that Ψ
^{±}
(
x
) ≥
0
for all
x
∈ Ω, which leads to the required bound on
u
. Consequently, the problem (1) − (2) has a unique and stable solution. □
To derive sharper bounds on the derivatives of solution, the solution is decomposed into a sum, composed of a regular component
v
and a singular component
w
. That is,
u
=
v
+
w
. The regular component
v
, can be written in the form
v
=
, where
v
_{0}
= (
v
_{01}
,
v
_{02}
)
^{T}
,
v
_{1}
= (
v
_{11}
,
v
_{12}
)
^{T}
and
v
_{2}
= (
v
_{21}
,
v
_{22}
)
^{T}
are defined respectively to be the solutions of the problems
and
Thus the regular component
v
is the solution of
Further we decompose
w
as
w
=
w
_{1}
+
w
_{2}
where
w
_{1}
= (
w
_{11}
,
w
_{12}
)
^{T}
,
w
_{2}
= (
w
_{21}
,
w
_{22}
)
^{T}
. Thus
w
_{1}
=
w
_{11}
+
w
_{21}
and
w
_{2}
=
w
_{12}
+
w
_{22}
, where
w
_{1}
is the solution of
and
w
_{2}
is the solution of
Lemma 2.3.
For each integer k, satisfying
0 ≤
k
≤ 3
, the regular component v and its derivatives satisfy the bounds given by
Proof
. The proof follows from
[7]
and
[2]
. □
Lemma 2.4.
For each integer k, satisfying
0 ≤
k
≤ 3
, the singular component w
_{1}
and its derivatives satisfy the bounds given by
Proof
. The proof follows from
[7]
and
[2]
. □
Lemma 2.5.
For each integer k, satisfying
0 ≤
k
≤ 3
, the singular component w_{2} and its derivatives satisfy the bounds given by
Proof
. Consider the barrier function
ϕ
(1, 1)
^{T}
±
w
_{2}
, where
to bound
w
_{2}
. To bound derivatives of
w
_{2}
, use the technique used in
[7]
and bound on
w
_{2}
on the domain Ω
_{1}
and Ω
_{2}
. □
3. Discretization of the Problem
We use piecewise uniform Shishkin mesh which uses these transition parameters:
The interior points of the mesh are denoted by
Let
h_{i}
=
x_{i}
−
x
_{i}
_{−1}
be the
i^{th}
mesh step and
, clearly
= {
x_{i}
:
i
= 0, 1,...,
N
}. Let
N
= 2
^{l}
,
l
≥ 5 be any positive integer.
We divide
into three subintervals [0,
σ_{єl1}
], [
σ_{єl1}
,
σ_{єl2}
] and [
σ_{єl2}
,
d
] for some 0 <
σ_{єl1}
≤
σ_{єl2}
≤
. The subintervals [0,
σ_{єl1}
] and [
σ_{єl1}
,
σ_{єl2}
] are divided into
N
/8 equidistant elements and the subinterval [
σ_{єl2}
,
d
] is divided into
N
/4 equidistant elements. Similarly, in
the subintervals [
d
,
d
+
σ_{єr1}
] and [
d
+
σ_{єr1}
,
d
+
σ_{єr2}
] are divided into
N
/8 equidistant elements and the subinterval [
d
+
σ_{єr2}
, 1] is divided into
N
/4 equidistant elements, for some 0 <
σ_{єr1}
≤
σ_{єr2}
≤
.
Define the discrete finite difference operator
L^{N}
as follows
with boundary conditions
where
and at
x
_{N}
_{/2}
=
d
the scheme is given by
where
Lemma 3.1.
Suppose that a mesh function Z
(
x_{i}
)
satisfies Z
(
x
_{0}
) ≥
0,Z
(
x_{N}
) ≥
0, L^{N}Z
(
x_{i}
) ≥
0 for all x_{i}
∈ Ω
^{N}
and
≤
0, then Z
(
x_{i}
) ≥
0 for all x_{i}
∈
.
Lemma 3.2.
If Z
(
x_{i}
)
is any mesh function, then,
The discrete solution
U
can be decomposed into the sum
U
=
V
+
W
. The function
V
, is defined as the solution of the following problem:
The function
W
, is defined as the solution of the following problem:
where the jump in the discrete derivative of a mesh function
Z
at the point
x_{i}
=
d
is given by:
Further decompose
W
as
W
=
W
_{1}
+
W
_{2}
, where the function
W
_{1}
is defined as the solution of the following problem:
and the function
W
_{2}
is defined as the solution of the following problem:
4. Convergence analysis
By Taylor’s expansion and bounds on regular components defined in lemma 2.3 gives
Define the mesh function Ψ
^{±}
(
x_{i}
) as
Using discrete maximum principle, the error of the regular component satisfies the estimate
As in
[7]
, the error of the singular component satisfies the estimate
Lemma 4.1.
The following ε
_{1}
,
ε
_{2}
−
uniform bound
where W
_{2}
is the solution of
(14) − (15).
Proof
. At the point
x
=
d
we know that
First consider
From lemma 2.3 we have
Therefore, 
D
−
V
(
d
) ≤
Similarly, consider
Again from lemma 2.3 we have
Therefore, 
D
+
V
(
d
) ≤
On Ω
_{1}
, 
W
_{1}
(
x_{i}
) ≤
implies that 
D
^{−}
W
_{1}
(
d
) ≤
. On Ω
_{2}
,
D
^{+}
W
_{1}
(
d
) =
D
^{+}
(
W
_{1}
−
w
_{1}
)(
d
) +
D
^{+}
w
_{1}
(
d
). From lemma 2.4
and 
D
^{+}
(
W
_{1}
−
w
_{1}
)(
x_{i}
) ≤
.
Therefore,
Lemma 4.2.
The following ε
_{1}
,
ε
_{2}
−
uniform bound
is valid, where W
_{2}
is the solution of
(14) − (15).
Proof
. Consider the following function
= 1, 2 where
where Ψ = (
ψ
_{1}
,
ψ
_{2}
)
^{T}
is the solution of
Using the discrete maximum principle we get the required result.
Lemma 4.3.
The error of the singular component satisfies the estimate
Proof
.
Now
From lemma 4.1 we have
and
Hence
Likewise,
Using the bounds on derivatives of
w
_{2}
given in lemma 2.5, we have
Using the bounds on derivatives of
w
_{1}
given in lemma 2.4, we have
Also,
Collecting all the previous inequalities we get that
By the Taylor’s expansion and bounds on the derivatives of
w
_{2}
, given in lemma 2.5 we have
Case (i) For
x_{i}
∈ [
d
+
σ_{єr2}
, 1].
Similar arguments prove a similar result for the subinterval [
σ_{єl1}
,
d
).
Case (ii) For
x_{i}
∈ (0,
σ_{єl1}
).
Similar arguments prove a similar result for the subinterval (
d
,
d
+
σ_{єr1}
).
Case (iii) [
σ_{єl1}
,
σ_{єl2}
).
Using the inequality
Likewise, (
L^{N}
(
W
_{2}
−
w
_{2}
))
_{2}
(
x_{i}
) ≤
CN
^{−1}
ln
N
.
Similar arguments prove a similar result for the subinterval [
d
+
σ_{єr1}
, d+
σ_{єr2}
).
Combining all these gives,
Consider the following function
,
j
= 1, 2 where
where Ψ = (
ψ
_{1}
,
ψ
_{2}
)
^{T}
is the solution of
Using the discrete maximum principle we get the required result. □
Theorem 4.4.
Let u be the solution of given problem
(1) − (2)
and U is the solution of discrete problem on the Shishkin mesh defined in section 3, then
Proof
. Using the equation (16), (17) and lemma 4.3 we get the required result. □
5. Numerical Results
To illustrate the theoretical results the scheme in Section 3 is implemented on these test examples.
Example 5.1
Consider the following singularly perturbed convectiondiffusion problem with discontinuous source term:
where
For the construction of piecewiseuniform Shishkin mesh
, we take
α
= 0.8. The Exact solution of the examples are not known. Therefore we estimate the error for
U
by comparing it to the numerical solution
obtained on the mesh
that contains the mesh points of the original mesh and their midpoints, that is,
=
x_{j}
, j=0,...,N,
= (
x_{j}
+
x
_{j}
_{+1}
)/2, j=0,...,N1.
For different values of N and
ε
_{1}
,
ε
_{2}
, we compute
Maximum pointwise errorsandε1,ε2−uniform rate of convergencepNfor Example 5.1.
Maximum pointwise errors and ε_{1}, ε_{2}−uniform rate of convergence p^{N} for Example 5.1.
If
ε
_{1}
= 10
^{−j}
for some nonnegative integer j , set
Then the parameteruniform error is computed as
D^{N}
:=
and the order of convergence is calculated using the formula
p^{N}
:=
Finally, we want to show that similar results can be obtained for coupled system of
M
(> 2) singularly perturbed convection diffusion problem with discontinuous source term. Letting
N
= 2
^{M}
×
τ
, where
τ
is some positive power of 2, the mesh is defined using the following transition points
Then we divide the interval [0,
d
] into M+1 subintervals [0,
σ_{єl1}
], [
σ_{єl1}
,
σ_{єl2}
],..., [
σ_{єlM−1}
,
σ_{єlM}
], [
σ_{єlM}
,
d
]. On the subinterval [0,
σ_{єl1}
] a uniform mesh of
N
/2
^{M+1}
mesh intervals, on [
σ_{єlk}
,
σ_{єlk+1}
], 1 ≤
k
≤
M
− 1, a uniform mesh of
N
/2
^{M−k+2}
mesh intervals, and on [
σ_{єlM}
,
d
] a uniform mesh of N/4 mesh intervals are placed. Similarly, we divide the interval [
d
, 1] into subintervals [
d
,
d
+
σ_{єr1}
], [
d
+
σ_{єr1}
,
d
+
σ_{єr2}
],..., [
d
+
σ_{єrM−1}
,
d
+
σ_{єrM}
], [
d
+
σ_{єrM}
, 1]. On the subinterval [
d
,
d
+
σ_{єr1}
] a uniform mesh of
N
/2
^{M+1}
mesh intervals, on [
d
+
σ_{єrk}
,
d
+
σ_{єrk}
+1], 1 ≤ k ≤ M − 1, a uniform mesh of
N
/2
^{M−k+2}
mesh intervals, and on [
d
+
σ_{єrM}
, 1] a uniform mesh of
N
/4 mesh intervals are placed. Let
h
_{єl1}
and
h
_{єr1}
be the mesh lengths on [0,
σ_{єl1}
] and on [
d
,
d
+
σ_{єr1}
] respectively. Let
H
_{1}
and
H
_{2}
be the mesh lengths on [
σ_{єlM}
,
d
] and on [
d
+
σ_{єrM}
, 1] respectively;
h_{єlk}
and
h_{єrk}
be the mesh lengths on [
σ_{єlk}
,
σ_{єlk+1}
] and on [
d
+
σ_{єrk}
,
d
+
σ_{єrk+1}
],
k
= 2,...,
M
respectively.
In this case also, we obtain the scheme similar to (3.1), with
u
= (
u
_{1}
;
u
_{2}
,...,
u_{M}
)
^{T}
∈
∩
C
^{1}
(Ω)
^{M}
∩
C
^{2}
(Ω
_{1}
∪ Ω
_{2}
)
^{M}
and also expect the error bound
≤
C
(
N
^{−1}
ln
N
) to hold, although attempts of the authors have failed so far to provide a proof. To illustrate the order of uniform convergence of this method we consider the following test example.
Example 5.2
Consider the following singularly perturbed convectiondiffusion problem with discontinuous source term:
where
f
_{1}
(
x
) =
and
Maximum pointwise errors, withε2= 10−4,ε3= 10−1andε1,ε2,ε3−uniform rate of convergencepNfor Example 5.2.
Maximum pointwise errors , with ε_{2} = 10^{−4}, ε_{3} = 10^{−1} and ε_{1}, ε_{2}, ε_{3}−uniform rate of convergence p^{N}for Example 5.2.
BIO
Sheetal Chawla is doing her Ph.D. at Indian Institute of Technology, Delhi under the direction of Dr. S. Chandra Sekhara Rao. Her area of interests are Numerical Analysis and Differential Equations.
Department of Mathematics, Indian Institute of Technology Delhi, New Delhi110016.
email: chawlaasheetal@gmail.com
S. Chandra Sekhara Rao is working as a Professor, Department of Mathematics, Indian Institute of Technology, Delhi. His area of interests are Numerical Analysis and Parallel Computing.
Department of Mathematics, Indian Institute of Technology Delhi, New Delhi110016.
email: scsr@maths.iitd.ac.in
Miller J.J.H.
,
Riordan E.O.
,
Shishkin G.I.
1996
Fitted Numerical Methods for Singular Perturbation Problems: Error Estimates in the Maximum Norm for Linear Problems in One and Two Dimensions
World Scientific
Singapore
Farrell P.A.
,
Hegarty A.F.
,
Miller J.J.H.
,
Riordan E.O.
,
Shishkin G.I.
2000
Robust Computational Techniques for Boundary Layers
Chapman & Hall, CRC Press
Boca Raton, Florida
Roos H.G.
,
Stynes M.
,
Tobiska L.
1996
Numerical Methods for Singularly Perturbed Differential Equations
Springer
Berlin
Roos H.G.
,
Stynes M.
,
Tobiska L.
2008
Robust Numerical Methods for Singularly Perturbed Differential Equations
Springer
Berlin
Farrell P.A.
,
Hegarty A.F.
,
Miller J.J.H.
,
Riordan E.O
,
Shishkin G.I.
(2004)
Singularly perturbed convectiondiffusion problems with boundary and weak interior layers
J. Comput. Appl. Math.
166
133 
151
DOI : 10.1016/j.cam.2003.09.033
Farrell P.A.
,
Hegarty A.F.
,
Miller J.J.H.
,
Riordan E.O.
,
Shishkin G.I.
(2004)
Global Maximum Norm ParameterUniform Numerical Method for a Singularly Perturbed ConvectionDiffusion Problem with Discontinuous Convection Coefficient
Math. Comput. Modelling
40
1375 
1392
DOI : 10.1016/j.mcm.2005.01.025
Cen Z.
(2005)
Parameteruniform finite difference scheme for a system of coupled singularly perturbed convectiondiffusion equations
Intern. J. Comput. Math.
82
177 
192
DOI : 10.1080/0020716042000301798
Linss T.
(2009)
Analysis of a system of singularly perturbed convectiondiffusion equations with strong coupling
SIAM J. Numer. Anal.
47
1847 
1862
DOI : 10.1137/070683970
Clavero C.
,
Gracia J.L.
,
Lisbona F.
(1999)
High order methods on Shishkin meshes for singular perturbation problems of convectiondiffusion type
Numer. Algorithms
22
73 
97
DOI : 10.1023/A:1019150606200
Riordan E.O.
,
Stynes M.
(2009)
Numerical analysis of a strongly coupled system of two singularly preturbed convectiondiffusion problems
Adv. Comput. Math.
30
101 
121
DOI : 10.1007/s104440079058z
Bellew S.
,
Riordan E.O.
(2004)
A parameter robust numerical method for a system of two singularly perturbed convectiondiffusion equations
Appl. Numer. Math.
51
171 
186
DOI : 10.1016/j.apnum.2004.05.006
Munyakazi J.B.
(2014)
A uniformly convergent nonstandard finite difference scheme for a system of convectiondiffusion equations
Comp. Appl. Math.
DOI : 10.1007/s4031401401716
Priyadharshini R.M.
,
Ramanaujam N.
(2013)
Uniformlyconvergent numerical methods for a system of coupled singularly perturbed convectiondiffusion equations with mixed type boundary conditions
Math. Model. Anal.
18
577 
598
DOI : 10.3846/13926292.2013.851629
Tamilselva A.
,
Ramanujam N.
,
Priyadharshini R.M.
,
Valanarasu T.
(2010)
Parameter  uniform numerical method for a system of coupled singularly perturbed convectiondiffusion equations with mixed type boundary conditions
J. Appl. Math. Informatics
28
109 
130
Tamilselvan A.
,
Ramanujam N.
(2009)
A numerical method for singularly perturbed system of second order ordinary differential equations of convection diffusion type with a discontinuous source term
J. Appl. Math. Informatics
27
1279 
1292