A UNIFORMLY CONVERGENT NUMERICAL METHOD FOR A WEAKLY COUPLED SYSTEM OF SINGULARLY PERTURBED CONVECTION-DIFFUSION PROBLEMS WITH BOUNDARY AND WEAK INTERIOR LAYERS

Abstract

We consider a weakly coupled system of singularly perturbed convection-diffusion 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.

1. Introduction

An extensive research had been done on numerical methods for a single singularly perturbed convection-diffusion 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 non-classical 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 convection-diffusion 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 convection-diffusion 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 convection-diffusion 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 u1, u2 ∈ C0() ∩ C1(Ω) ∩ C2(Ω1 ∪ Ω2), such that

where E = diag(ε1, ε2), the coupling matrix A = diag(a1, a2) and B = (bij)2×2 with 0 < ε1 ≤ ε2 ≤ 1, f = (f1, f2)T , and u = (u1, u2)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 f1, f2 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 vi ≤ wi, i = 1, 2, and |v| = (|v1|, |v2|)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 = (v1, v2)T ∈ C(S)2, we define and ║v║S= max{║v1║S,║v2║S}. Now let a mesh be a set of points satisfying x0 < x1 <⋯< xN = 1. A mesh function is a real-valued function defined on ΩN. Define the discrete maximum norm for such functions by and for vector mesh functions V = (V1, V2)T = are used and define║V║ΩN = max{║V1║ΩN,║V2║ΩN}.

2. The continuous problem

Theorem 2.1 (Continuous maximum principle). Suppose u1, u2 ∈ C0() ∩ C1(Ω) ∩ C2(Ω1 ∪ Ω2). Further suppose that u = (u1, u2)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 = (u1, u2)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 ∈ C1(Ω)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 v0 = (v01, v02)T, v1 = (v11, v12)T and v2 = (v21, v22)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 = w1 + w2 where w1 = (w11,w12)T , w2 = (w21,w22)T . Thus w1 = w11+w21 and w2 = w12+w22, where w1 is the solution of

and w2 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 w1 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 w2 and its derivatives satisfy the bounds given by

Proof. Consider the barrier function ϕ(1, 1)T ± w2, where

to bound w2. To bound derivatives of w2, use the technique used in [7] and bound on w2 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 hi = xi − xi−1 be the ith mesh step and , clearly = {xi : i = 0, 1,...,N}. Let N = 2l , l ≥ 5 be any positive integer.

We divide into three sub-intervals [0, σєl1 ], [σєl1 , σєl2 ] and [σєl2 , d] for some 0 < σєl1 ≤ σєl2 ≤ . The sub-intervals [0, σєl1 ] and [σєl1 , σєl2 ] are divided into N/8 equidistant elements and the sub-interval [σєl2 , d] is divided into N/4 equidistant elements. Similarly, in the sub-intervals [d, d + σєr1 ] and [d + σєr1 , d + σєr2 ] are divided into N/8 equidistant elements and the sub-interval [d+σєr2 , 1] is divided into N/4 equidistant elements, for some 0 < σєr1 ≤ σєr2 ≤ .

Define the discrete finite difference operator LN as follows

with boundary conditions

where

and at xN/2 = d the scheme is given by

where

Lemma 3.1. Suppose that a mesh function Z(xi) satisfies Z(x0) ≥ 0,Z(xN) ≥ 0, LNZ(xi) ≥ 0 for all xi ∈ ΩN and ≤ 0, then Z(xi) ≥ 0 for all xi ∈ .

Lemma 3.2. If Z(xi) 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 xi = d is given by:

Further decompose W as W = W1 +W2, where the function W1 is defined as the solution of the following problem:

and the function W2 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 Ψ±(xi) 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 W2 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, |W1(xi)| ≤ implies that |D−W1(d)| ≤ . On Ω2, D+W1(d) = D+(W1 − w1)(d) + D+w1(d). From lemma 2.4 and |D+(W1 − w1)(xi)| ≤ .

Therefore,

Lemma 4.2. The following ε1, ε2− uniform bound

is valid, where W2 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 w2 given in lemma 2.5, we have

Using the bounds on derivatives of w1 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 w2, given in lemma 2.5 we have

Case (i) For xi ∈ [d + σєr2 , 1].

Similar arguments prove a similar result for the sub-interval [σєl1 , d).

Case (ii) For xi ∈ (0, σєl1 ).

Similar arguments prove a similar result for the sub-interval (d, d + σєr1 ).

Case (iii) [σєl1 , σєl2 ).

Using the inequality

Likewise, |(LN(W2 − w2))2(xi)| ≤ CN−1 ln N.

Similar arguments prove a similar result for the sub-interval [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 convection-diffusion problem with discontinuous source term:

where

For the construction of piecewise-uniform 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, = xj , j=0,...,N, = (xj + xj+1)/2, j=0,...,N-1.

For different values of N and ε1, ε2, we compute

Table 1.Maximum point-wise errors and ε1, ε2−uniform rate of convergence pN for Example 5.1.

If ε1 = 10−j for some non-negative integer j , set

Then the parameter-uniform error is computed as DN := and the order of convergence is calculated using the formula pN :=

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 = 2M ×τ , 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/2M+1 mesh intervals, on [σєlk, σєlk+1], 1 ≤ k ≤ M − 1, a uniform mesh of N/2M−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/2M+1 mesh intervals, on [d + σєrk, d + σєrk+1], 1 ≤ k ≤ M − 1, a uniform mesh of N/2M−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 H1 and H2 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 = (u1; u2,...,uM)T ∈ ∩ C1(Ω)M ∩ C2(Ω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 convection-diffusion problem with discontinuous source term:

where f1(x) =

and

Table 2.Maximum point-wise errors , with ε2 = 10−4, ε3 = 10−1 and ε1, ε2, ε3−uniform rate of convergence pNfor Example 5.2.