MONOTONE FINITE DIFFERENCE DOMAIN
DECOMPOSITION ALGORITHMS AND
APPLICATIONS TO NONLINEAR SINGULARLY
PERTURBED REACTION-DIFFUSION PROBLEMS
IGOR BOGLAEV AND MATTHEW HARDY
Received 16 September 2004; Revised 21 Decembe r 2004; Accepted 11 January 2005
This paper deals with monotone finite difference iterative algorithms for solving non-
linear singularly perturbed reaction-diffusion problems of elliptic and parabolic types.
Monotone domain decomposition algorithms based on a Schwarz alternating method
and on box-domain decomposition are constructed. These monotone algorithms solve
only linear discrete systems at each iterative step and converge monotonically to the ex-
act solution of the nonlinear discrete problems. The rate of convergence of the mono-
tone domain decomposition algorithms are estimated. Numerical experiments are pre-
sented.
Copyright © 2006 Hindawi Publishing Corporation. All rights reserved.
1. Introduction
We are interested in monotone discrete Schwarz alternating algorithms for solving non-
linear singularly perturbed reaction-diffusion problems.
The first problem considered corresponds to the singularly perturbed reaction-diffu-
sion problem of elliptic type
−μ
2

u
xx
+ u
yy

+ f (x, y,u) = 0, (x, y) ∈ ω,
u
= g on ∂ω, ω = ω

2

u
xx
+ u
yy

+ f (x, y,t,u)+u
t
= 0, (x, y) ∈ ω, t ∈ (0,T],
f
u
≥ 0, (x, y,t,u) ∈ ω × [0,T] × (−∞,∞),
(1.2)
where ω
={0 <x<1}×{0 <y<1} and μ is a small positive parameter. The initial-
boundary conditions are defined by
u(x, y,0)
= u
0
(x, y), (x, y) ∈ ω,
u(x, y, t)
= g(x, y,t), (x, y,t) ∈ ∂ω × (0,T].
(1.3)
The functions f , g,andu
0
are sufficiently smooth. Under suitable continuity and com-
patibility conditions on the data, a unique solution u of (1.2) exists (see [5] for details).
For μ
 1, problem (1.2) is singularly perturbed and characterized by the boundary lay-

and horizontal interfacial subdomain problems provide Dirichlet data on γ and ρ for the
problems on the nonoverlapping box-subdomains.
I. Boglaev and M. Hardy 3
In Section 2, we introduce the classical nonlinear finite difference schemes for the nu-
merical solution of (1.1)and(1.2). Iterative methods by which each of these schemes
may be solved are presented in [3, 4]. From an arbitrary initial mesh function, one may
construct a sequence of functions which converges monotonically to the exact solution
of the nonlinear difference scheme. Each function in the sequence is generated as the so-
lution of a linear difference problem. In Section 3, we consider the elliptic problem and
extend the monotone method to a box-decomposition of the computational domain.
We show that monotonic convergence is maintained under the proposed decomposition
and associated algorithm. Further, we develop estimates of the rate of convergence. The
box-decomposition of the spatial domain is applied to the parabolic nonlinear difference
scheme in Section 4. Numerical experiments are presented in Section 5. These confirm
the theoretical estimates of the earlier sections. Suggestions are made regarding future
parallel implementation.
2. Difference schemes for solving (1.1)and(1.2)
On
ω and [0,T] introduce nonuniform meshes ω
h
= ω
hx
× ω
hy
and ω
τ
:
ω
hx
=

y
= 1; h
yj
= y
j+1
− y
j

,
ω
τ
=

t
k
= kτ,0≤ k ≤ N
τ
, N
τ
τ = T

.
(2.1)
For approximation of the elliptic problem (1.1), we use the classical difference scheme
on nonuniform meshes
ᏸ
h
U + f (P,U) = 0, P ∈ ω
h
, U = g on ∂ω

=


xi

−1


U
i+1, j
− U
ij

h
xi

−1
−

U
ij
− U
i−1, j

h
x,i−1

−1

,

y,j−1

−1

,

xi
= 2
−1

h
x,i−1
+ h
xi

, 
yj
= 2
−1

h
y,j−1
+ h
yj

,
(2.4)
where P
= (x
i

= u
0
(P), P ∈ ω
h
, U(P,t) = g(P, t), (P,t) ∈ ∂ω
h
× ω
τ
,
(2.5)
where ᏸ
h
is defined in (2.3).
Consider the linear versions of problems (2.2)and(2.5)
ᏸW + c(P)W(P)
= F(P), P ∈ ω
h
,
W(P)
= W
0
(P), P ∈ ∂ω
h
, c(P) ≥ c
0
> 0, P ∈ ω
h
, c
0
= const,


∂ω
h
, F
ω
h
/

c
0
+ βτ
−1

,


W
0


∂ω
h
≡ max
P∈∂ω
h


W
0
(P)


x
m−1
,x
m

×

y
l−1
, y
l

, x
0
= 0, x
M
= 1, y
0
= 0, y
L
= 1.
(3.1)
Additionally, we introduce (M
− 1) interfacial subdomains θ
m
, m = 1, , M − 1(ver-
tical strips):
θ
m

, γ
e
m
=

x = x
e
m
,0≤ y ≤ 1

,
x
b
m
<x
m
<x
e
m
, γ
0
m
= ∂ω ∩ ∂θ
m
,
(3.2)
I. Boglaev and M. Hardy 5
θ
m−1
x

l
ϑ
l
θ
m
x
b
m
−1
x
e
m
−1
x
b
m
x
e
m
ω
m,l+l
Figure 3.1. Fragment of the domain decomposition.
and (L − 1) interfacial subdomains ϑ
l
, l = 1, ,L − 1 (hor izontal strips):
ϑ
l
= ω
x
× ϑ

0 ≤ x ≤ 1, y = y
e
l

,
y
b
l
<y
l
<y
e
l
, ρ
0
l
= ∂ω ∩ ∂ϑ
l
.
(3.3)
Figure 3.1 illustrates a fragment of the domain decomposition.
On
ω
ml
, m = 1, ,M, l = 1, ,L; θ
m
, m = 1, ,M − 1andϑ
l
, l = 1, ,L − 1, int ro-
duce meshes:

m

M−1
m
=1
∈ ω
hx
,

y
b
l
, y
l
, y
e
l

L−1
l
=1
∈ ω
hy
,
(3.4)
with
ω
hx
, ω
hy

h
satisfying the boundary conditions V
(0)
(P) = g(P)on∂ω
h
.
6 Monotone domain decomposition algorithms
Step 2. On subdomains
ω
h
ml
, m = 1, ,M, l = 1, ,L, compute mesh functions V
(n+1)
ml
(P)
(here the index n stands for a number of iterative steps) satisfy ing the following difference
problems

ᏸ
h
+ c
∗

Z
(n+1)
ml
=−G
(n)
(P), P ∈ ω
h

ml
.
(3.5)
Step 3. On the vertical interfacial subdomains
θ
h
m
, m = 1, , M − 1, compute the differ-
ence problems

ᏸ
h
+ c
∗

Z
(n+1)
m
=−G
(n)
(P), P ∈ θ
h
m
,
Z
(n+1)
m
(P) =
⎧
⎪

h
m+1,l
, l = 1, , L,
V
(n+1)
m
(P) = V
(n)
(P)+Z
(n+1)
m
(P), P ∈ θ
h
m
,
(3.6)
where we use the notation
γ
h0
m
= γ
0
m
∩ ∂ω
h
, γ
hb
m
= γ
b

=−G
(n)
(P), P ∈ ϑ
h
l
,

Z
(n+1)
l
(P) =
⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩


∩
ω
h
m,l+1
, m = 1, , M;
Z
(n+1)
m
(P), P ∈ ∂ϑ
h
l
∩ θ
h
m
, m = 1, , M − 1,

V
(n+1)
l
(P) = V
(n)
(P)+

Z
(n+1)
l
(P), P ∈ ϑ
h
l

, ρ
hb
l
= ρ
b
l
∩ ϑ
h
l
, ρ
he
l
= ρ
e
l
∩ ϑ
h
l
.
(3.9)
I. Boglaev and M. Hardy 7
Step 5. Compute the mesh function V
(n+1)
(P), P ∈ ω
h
by piecing together the solutions
on the subdomains
V
(n+1)
(P) =

m
\ ϑ
h
, m = 1, , M − 1;

V
(n+1)
l
(P), P ∈ ϑ
h
l
, l = 1, , L − 1.
(3.10)
Step 6. Stopping criterion: If a prescribed accuracy is reached, then stop; otherwise go to
Step 2.
Algorithm (3.5)–(3.10) can be carried out by parallel processing. Steps 2, 3,and4
must be performed sequentially, but on each step, the independent subproblems may be
assigned to different computational nodes.
Remark 3.1. We note that the original Schwarz alternating algorithm with overlapping
subdomains is a purely sequential algorithm. To obtain parallelism, one needs a subdo-
main colouring strategy, so that a set of independent subproblems can be introduced.
The modification of the Schwarz algorithm (3.5)–(3.10) can be considered as an additive
Schwarz algorithm.
3.2. Monotone convergence of algorithm (3.5)–(3.10). Additionally, we assume that f
from (1.1) satisfies the two-sided constraints
0 <c
∗
≤ f
u
≤ c

h
,
δV(P)
≥ 0, P ∈ ∂ω
h
,
(3.14)
where f
u
(P) ≡ f
u
[P,V(P)+Θ(P)δV(P)], 0 < Θ(P) < 1. In view of the maximum princi-
ple in Lemma 2.1,weconclude(3.13).
The following convergence property of algorithm (3.5)–(3.10)holdstrue.
8 Monotone domain decomposition algorithms
Theorem 3.2. Let
V
(0)
and V
(0)
beupperandlowersolutionsof(2.2), and let f (x, y,u) sat-
isfy (3.11). Then the upper sequence
{V
(n)
} generated by (3.5)–(3.10) converges monotoni-
cally from above to the unique solution U of (2.2), and the lower sequence
{V
(n)
} generated
by (3.5)–(3.10) converges monotonically from below to U:

(P), we obtain the difference
equation for V
(n+1)
ml
ᏸ
h
V
(n+1)
ml
+ f

P,V
(n+1)
ml

=−

c
∗
− f
(n)
u,ml
(P)

Z
(n+1)
ml
(P) ≥ 0, P ∈ ω
h
ml

(3.17)
where nonnegativeness of the right-hand side of the difference equation follows from
(3.11)and(3.16).
Taking into acco un t (3.16)and
V
(n)
is an upper solution, by the maximum pr inciple
in Lemma 2.1,from(3.6)and(3.8) it follows that
Z
(n+1)
m
(P) ≤ 0, P ∈ θ
h
m
, m = 1, , M − 1,

Z
(n+1)
l
(P) ≤ 0, P ∈ ϑ
h
l
, l = 1, , L − 1.
(3.18)
Similar to (3.17), we obtain the difference problems for V
(n+1)
m
ᏸ
h
V

⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎩
g(P), P ∈ γ
h0
m
;
V
(n+1)
ml
(P), P ∈ γ
hb
m
∩ ω
h
ml
, l = 1, , L;
V
(n+1)
m+1,l
(P), P ∈ γ
he
m
∩ ω

u,l
(P)


Z
(n+1)
l
(P) ≥ 0, P ∈ ϑ
h
l
,

V
(n+1)
l
(P) =
⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪

l
\ θ
h

∩
ω
h
m,l+1
, m = 1, , M;
V
(n+1)
m
(P), P ∈ ∂ϑ
h
l
∩ θ
h
m
, m = 1, , M − 1,
(3.20)
where nonnegativeness of the right-hand sides of the difference equations follows from
(3.11)and(3.18). Now we verify that the mesh function
V
(n+1)
defined by (3.10)isanup-
per solution. From the boundary conditions for V
(n+1)
ml
, V
(n+1)

γ
h
∪ ρ
h

,
γ
hb,e
ml
=

x
i
= x
b,e
m
, y
e
l
−1
<y
j
<y
b
l

, γ
hb,e
m
=

l
.
(3.21)
To p r ov e tha t
V
(n+1)
is an upper solution of problem (2.2), we have to verify only that
the last inequality holds true on the interfacial boundaries
γ
hb,e
ml
and ρ
hb,e
l
, m = 1, ,M − 1,
l
= 1, , L − 1.
We check this inequality in the case of the left interfacial boundary
γ
hb
ml
, since the case
with
γ
he
ml
is checked in a similar way. From (3.5), (3.6), and (3.18), we conclude that the
mesh function W
(n+1)
ml

⎧
⎨
⎩
0, P ∈ γ
hb
ml
= γ
hb
m
∩ ω
h
ml
;
≥ 0, P ∈ ∂θ
h
ml
\ γ
hb
ml
.
(3.22)
In view of the maximum principle in Lemma 2.1,
V
(n+1)
ml
(P) − V
(n+1)
m
(P) ≥ 0, P ∈ θ
h

,
−μ
2
Ᏸ
2
x
V
(n+1)
ml
(P) ≤−μ
2
Ᏸ
2
x
V
(n+1)
(P), P ∈ γ
hb
ml
.
(3.24)
10 Monotone domain decomposition algorithms
Thus, using (3.17), we conclude
G
(n+1)
(P) ≥ ᏸ
h
V
(n+1)
ml


V
(n+1)
l
satisfies the difference problem

ᏸ
h
+ c
∗


W
(n+1)
ml
= 0, P ∈ ϑ
h
ml
= ω
h
ml
∩ ϑ
h
l
,

W
(n+1)
ml
(P) =

.
(3.26)
By the maximum principle in Lemma 2.1,
V
(n+1)
ml
(P) −

V
(n+1)
l
(P) ≥ 0, P ∈ ϑ
h
ml
. (3.27)
By (3.8),

V
(n+1)
l
(P) = V
(n+1)
ml
(P), P ∈ ρ
hb
ml
∪{(x
e
m
−1

−μ
2
Ᏸ
2
y
V
(n+1)
ml
(P) ≤−μ
2
Ᏸ
2
y
V
(n+1)
(P), P ∈ ρ
hb
ml
.
(3.28)
Thus, using (3.17), we conclude
G
(n+1)
(P) ≥ ᏸ
h
V
(n+1)
ml
+ f



W
(n+1)
ml
= 0, P ∈ τ
h
ml
= θ
h
m
∩ ϑ
h
l
,

W
(n+1)
ml
(P) =
⎧
⎨
⎩
0, P ∈ ρ
hb,e
ml
=

x
b
m

(P) −

V
(n+1)
l
(P) ≥ 0, P ∈ τ
h
ml
. (3.31)
By (3.8),

V
(n+1)
l
(P)= V
(n+1)
m
(P), P ∈ ρ
hb
ml
∪{(x
e
m
, y
b
l
),(x
b
m
, y

m
(P) ≤−μ
2
Ᏸ
2
y
V
(n+1)
(P), P ∈ ρ
hb
ml
.
(3.32)
I. Boglaev and M. Hardy 11
Thus, using (3.19), we conclude
G
(n+1)
(P) ≥ ᏸ
h
V
(n+1)
m
+ f

P,V
(n+1)
m

≥
0, P ∈ ρ

b
l
), we have
V
(n+1)
ml

P
b
ml

=
V
(n+1)
m

P
b
ml

=

V
(n+1)
l

P
b
ml


V
(n+1)
ml

P
b
ml

h
b+
xm
−
V
(n+1)
ml

P
b
ml

−
V
(n+1)
ml

P
bx−
ml

h

ml

−
V
(n+1)
ml

P
b
ml

h
b+
yl
−
V
(n+1)
ml

P
b
ml

−
V
(n+1)
ml

P
by−

, y
b
l
± h
b±
yl

,

b
xm
= 2
−1

h
b−
xm
+ h
b+
xm

, 
b
yl
= 2
−1

h
b−
yl

h
V
(n+1)
ml
+ f

P,V
(n+1)
ml

≥
0, P = P
b
ml
. (3.36)
With a similar argument for mesh point P
e
l
= (x
e
m
, y
b
l
), we prove that V
(n+1)
is an upper
solution of problem (2.2) on the whole computational domain
ω
h

h
, (3.37)
is an exact solution to (2.2). The uniqueness of the solution to (2.2) follows from esti-
mate (2.8). Indeed, if by contradiction, we assume that there exist two solutions U
1
and
U
2
to (2.8), then by the mean-value theorem, the difference δU = U
1
− U
2
satisfies the
following difference problem
ᏸ
h
δU + f
u
δU = 0, P ∈ ω
h
, δU = 0, P ∈ ∂ω
h
. (3.38)
12 Monotone domain decomposition algorithms
By (2.8), δU
= 0 which leads to the uniqueness of the solution to (2.2). This proves the
theorem.

Remark 3.3. Consider the following approach for constructing initial upper and lower
solutions

,
Z
(0)
ν
(P) = 0, P ∈ ∂ω
h
, ν = 1, −1.
(3.39)
Then the functions
V
(0)
= R + Z
(0)
1
, V
(0)
= R + Z
(0)
−1
are upper and lower solutions, respec-
tively. The proof of this result can be found in [4].
Remark 3.4. Since the initial iteration in algorithm (3.5)–(3.10)iseitheranupperora
lower solution, which can be constructed directly from the difference equation without
any knowledge of the solution as we have suggested in the previous remark, this algorithm
eliminates the search for the initial iteration as is often needed in Newton’s method. This
gives a practical advantage in the computation of numerical solutions.
3.3. Convergence analysis of algorithm (3.5)–(3.10). We now establish convergence
properties of algorithm (3.5)–(3.10).
If we denote
Z

ml
\

θ
h
∪ ϑ
h

;
Z
(n+1)
m
(P), P ∈ θ
h
m
\ ϑ
h
, m = 1, , M − 1;

Z
(n+1)
l
(P), P ∈ ϑ
h
l
, l = 1, ,L − 1.
(3.41)
Introduce the following notation

b,e

m
and h
b±,e±
yl
are
the mesh step sizes on the top and bottom from points y
b,e
l
,and
κ
b
xm
≡
μ
2
c
∗

b
xm
h
b+
xm
, κ
e
xm
≡
μ
2
c

yl
h
b+
yl
, κ
e
yl
≡
μ
2
c
∗

e
yl
h
e−
yl
, q
II
= max
1≤l≤L−1

κ
b
yl
;κ
e
yl


= 1 − c
∗
/c
∗
.
Proof. Suppose that the sequence
{V
(n)
} is generated by algorithm (3.5)–(3.10). Using
(2.8), from (3.5) we get the following estimate on Z
(n+1)
ml


Z
(n+1)
ml


ω
h
ml
≤
1
c
∗


G
(n)

1≤l≤L



Z
(n+1)
ml


γ
hb
ml
;


Z
(n+1)
m+1,l


γ
he
ml


≤
1
c
∗




Z
(n+1)
l


ϑ
h
l
≤
1
c
∗


G
(n)


ω
h
. (3.47)
Thus, by the definition of Z
(n+1)
,wehave


Z
(n+1)

(n)
u
(P)

Z
(n)
(P), P ∈ ω
h
, ω
h
= ω
h
\

γ
h
∪ ρ
h

, (3.49)
where
γ
h
and ρ
h
are defined in (3.21). By (3.11),
1
c
∗


e
l
−1
<y
j
<y
b
l
},werepresentG
(n)
in
the form
G
(n)


P
b
m

=
ᏸ
h
V
(n)
ml
+ f


P


P
b+
m

,

P
b
m
=

x
b
m
, y
j

∈ 
γ
hb
ml
,

P
b+
m
=

x

m

≤
V
(n−1)


P
b+
m

−
V
(n)


P
b+
m

=−
Z
(n)


P
b+
m

. (3.52)


Z
(n)


ω
h
. (3.53)
Similarly, we can prove the estimate
1
c
∗


G
(n)


γ
he
ml
≤

q + κ
e
xm



Z



ω
h
. (3.55)
On
ρ
hb
ml
={x
e
m
−1
<x
i
<x
b
m
, y
j
= y
b
l
},werepresentG
(n)
in the form
G
(n)

P

(n)
l

P
b+
l

−
V
(n)
ml

P
b+
l

,
P
b
l
=

x
i
, y
b
l

∈ 
ρ


V
(n)
l

P
b+
l

≤
V
(n−1)

P
b+
l

−
V
(n)

P
b+
l

=−
Z
(n)

P

yl



Z
(n)


ω
h
. (3.58)
Similarly, we can prove the estimate
1
c
∗


G
(n)


ρ
he
ml
≤

q + κ
e
yl


(n)

P
b
l

=
ᏸ
h
V
(n)
m
+ f

P
b
l
,V
(n)
m

−
μ
2

b
yl
h
b+
yl


∈ 
ρ
hb
ml
, P
b+
l
=

x
i
, y
b
l
+ h
b+
yl

.
(3.60)
From ( 3.18) at the iterative step n and the definition of V
(n)
,wehave
V
(n)
m

P
b+

(n)

P
b+
l

. (3.61)
I. Boglaev and M. Hardy 15
From here and (3.19), and taking into account that Z
(n)
m
(P) = Z
(n)
(P), P ∈ ρ
hb
ml
,wegetthe
estimate
1
c
∗


G
(n)


ρ
hb
ml

q + κ
e
yl



Z
(n)


ω
h
. (3.63)
At P
b
ml
= (x
b
m
, y
b
l
), we represent G
(n)
in the form
G
(n)

P
b


P
bx+
ml

−
V
(n)
ml

P
bx+
ml

−
μ
2

b
yl
h
b+
yl


V
(n)
l

P

ml
=

x
b
m
, y
b
l
+ h
b+
yl

.
(3.64)
From ( 3.16) at the iterative step n and the definition of V
(n)
,wehave
V
(n)
ml

P
bx+
ml

−
V
(n)
m

V
(n)
ml

P
by+
ml

−

V
(n)
l

P
by+
ml

≤
V
(n−1)

P
by+
ml

−
V
(n)




G
(n)

P
b
ml



≤

q + κ
b
xm
+ κ
b
yl



Z
(n)


ω
h
. (3.66)
By the same reasonings, the following estimate holds true


ω
h
, P
e
m
−1,l
=

x
e
m
−1
, y
b
l

. (3.67)
On ρ
hb
l
,weconcludetheestimate
1
c
∗


G
(n)



G
(n)


ρ
h
≤

q + q
I
+ q
II



Z
(n)


ω
h
. (3.69)
16 Monotone domain decomposition algorithms
From here, (3.50)and(3.55), we conclude the estimate
1
c
∗



that is, g(P)
= 0. This assumption can always be obtained via a change of variables. Let
the initial function V
(0)
be chosen in the form of (3.39)withR(P) = 0, that is, V
(0)
is the
solution of the following difference problem

ᏸ
h
+ c
∗

V
(0)
= ν


f (P,0)


, P ∈ ω
h
,
V
(0)
(P) = 0, P ∈ ∂ω
h
, ν = 1,−1.


f (P,0)


ω
h
, c
0
=
3c
∗
+ c
∗
c
∗
c
∗
, (3.72)
where U is the solution to (2.2).
Proof. Using (3.44), we have


V
(n+k)
− V
(n)


ω
h



Z
(n)


ω
h
≤
(q)
n
1 − q


Z
(1)


ω
h
.
(3.73)
Taking into account that limV
(n+k)
= U as k →∞,whereU is the solution to (2.2), we
conclude the estimate


V
(n)

∗


G
(0)


ω
h
≤
1
c
∗


ᏸ
h
V
(0)


ω
h
+
1
c
∗


f

h

+
1
c
∗


f (P,0)


ω
h
+


V
(0)


ω
h
.
(3.75)
I. Boglaev and M. Hardy 17
From here and estimating V
(0)
from (3.71)with(2.8),



0


f (P,0)


ω
h
, (3.77)
where c
0
is defined in (3.72). Thus, from here and (3.74), we prove (3.72). 
Remark 3.8. In the next section, we present sufficientconditionstoguaranteethein-
equality
q<1requiredinTheorem 3.7.
3.4. Uniform convergence of the monotone domain decomposition algorithm (3.5)–
(3.10). Here we analyze a convergence rate of algorithm (3.5)–(3.10)appliedtothedif-
ference scheme (2.2) defined on the piecewise uniform mesh introduced in [7]. On this
mesh, the difference scheme (2.2)convergesμ-uniformly to the solution of (1.1).
The piecewise uniform mesh is formed in the following manner. We divide each of the
intervals
ω
x
= [0,1] and ω
y
= [0,1] into three parts each [0,σ
x
], [σ
x
,1− σ

/4 + 1 mesh points, respectively, and in the parts [σ
x
,1 − σ
x
], [σ
y
,1 − σ
y
]with
N
x
/2+1 andN
y
/2 + 1 mesh points, respectively. This defines the piecewise equidistant
mesh in the x-andy-directions condensed in the boundary layers at x
= 0,1 and y = 0,1:
x
i
=
⎧
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎩
ih
xμ

x
/4+1, ,N
x
,
y
j
=
⎧
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎩
jh
yμ
, j = 0,1, ,N
y
/4;
σ
y
+

j − N
y
/4

h

y
are the step sizes inside and outside the boundary layers, respec-
tively. The transition points σ
x
,(1− σ
x
)andσ
y
,(1− σ
y
)aredeterminedby
σ
x
= min

4
−1
,μc
−(1/2)
∗
lnN
x

, σ
y
= min

4
−1
,μc

lnN
x
, N
−1
x
<h
x
< 2N
−1
x
,
σ
y
= μc
−(1/2)
∗
lnN
y
, h
yμ
= 4μc
−(1/2)
∗
N
−1
y
lnN
y
, N
−1

,N
y

, (3.81)
where constant C is independent of μ and N. The proof of this result can be found in [7].
Theorem 3.9. Let the interfacial subdomains
θ
h
m
, m = 1, ,M − 1 and ϑ
h
l
, l = 1, ,L − 1
be located in the x-andy-directions, respectively, outside the boundary layers. Assume μ
≤
μ
0
 1, and the following condition
N ≤
α

c
∗
/2
μ
0
, N = max

N
x

Q)
n
(1 −

Q)


f (P,0)


ω
h
,

Q = 1 −

1 − α
2

c
∗
c
∗
< 1, c
0
=
3c
∗
+ c
∗

2
x
<

μ
0
N

2
c
∗
, q
II
=
μ
2
c
∗
h
2
y
<

μ
0
N

2
c
∗

∗
as follows.
Step 1. Initialization: on the mesh
ω
h
,chooseV
(0)
(P,t), P ∈ ω
h
satisfying the boundary
condition V
(0)
(P,t) = g(P, t)on∂ω
h
.
For n
= 0 to n
∗
− 1 do Steps 2–5
Step 2. On subdomains
ω
h
ml
, m = 1, , M, l = 1, ,L, compute mesh functions V
(n+1)
ml
(P,
t), m
= 1, ,M, l = 1, ,L satisfying the following difference problems


(n+1)
ml
(P,t) = 0, P ∈ ∂ω
h
ml
,
V
(n+1)
ml
(P,t) = V
(n)
(P,t)+Z
(n+1)
ml
(P,t), P ∈ ω
h
ml
.
(4.1)
Step 3. On the vertical interfacial subdomains
θ
h
m
, m = 1, , M − 1, compute the differ-
ence problems

ᏸ
hτ
+ c
∗

(n+1)
ml
(P,t), P ∈ γ
hb
m
∩ ω
h
ml
, l = 1, , L;
Z
(n+1)
m+1,l
(P,t), P ∈ γ
he
m
∩ ω
h
m+1,l
, l = 1, , L,
V
(n+1)
m
(P,t) = V
(n)
(P,t)+Z
(n+1)
m
(P,t), P ∈ θ
h
m

(P,t) =
⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩
0, P ∈ ρ
h0
l
;
Z
(n+1)
ml
(P,t), P ∈

ρ

h
l
∩ θ
h
m
, m = 1, , M − 1,

V
(n+1)
l
(P,t) = V
(n)
(P,t)+

Z
(n+1)
l
(P,t), P ∈ ϑ
h
l
,
(4.3)
where we use the notation from (3.8).
Step 5. Compute the mesh function V
(n+1)
(P,t), P ∈ ω
h
by piecing together the solutions
on the subdomains
V

(P,t), P ∈ θ
h
m
\ ϑ
h
, m = 1, , M − 1;

V
(n+1)
l
(P,t), P ∈ ϑ
h
l
, l = 1, , L − 1.
(4.4)
Step 6. Set up
V(P,t)
= V
(n
∗
)
(P,t), P ∈ ω
h
. (4.5)
4.2. Monotone convergence of algorithm (4.1)–(4.5). Additionally, we assume that f
from (1.2) satisfies the two-sided constraints
0
≤ f
u
≤ c

V
(P,t) ≤ V(P,t), P ∈ ω
h
. (4.8)
The proof of this result is similar to that of (3.13).
On each time level t
∈ ω
τ
, we have the following convergence property of algorithm
(4.1)–(4.5).
I. Boglaev and M. Hardy 21
Theorem 4.1. Let V(P,t
− τ) be given and V
(0)
(P,t), V
(0)
(P,t) beupperandlowersolu-
tions corresponding to V(P,t
− τ).Supposethat f satis fies (4.6). Then the upper sequence
{V
(n)
(P,t)} generated by (4.1)–(4.5) converges monotonically from above to the unique so-
lution ᐂ(P, t) of the problem
ᏸ
hτ
V(P,t)+ f (P,t,V) − τ
−1
V(P,t − τ) = 0, P ∈ ω
h
,

Proof. The proof of the theorem is similar to the proof of Theorem 3.2 and based on the
maximum pr inciple in Lemma 2.1 and the estimate (2.8)withβ
= 1forthedifference
operator ᏸ
hτ
. 
Remark 4.2. Inthecaseofalgorithm(4.1)–(4.5), Remarks 3.1–3.4 hold still tr ue at each
time step t
∈ ω
τ
. We only mention that the difference problem in (3.39)becomes
ᏸ
hτ
Z
(0)
ν
= ν


ᏸ
hτ
R(P,t)+ f (P,t,R) − τ
−1
V(P,t − τ)


, P ∈ ω
h
,
Z

≡
μ
2

c
∗
+ τ
−1


e
xm
h
e−
xm
,
υ
b
yl
≡
μ
2

c
∗
+ τ
−1


b

xm
;υ
e
xm

, r
II
= max
1≤l≤L−1

υ
b
yl
;υ
e
yl

,
(4.12)
where all the step sizes are defined in (3.42).
Similar to Theorem 3.5,oneachtimelevelt
∈ ω
τ
, we have the following convergence
property of algorithm (4.1)–(4.5).
Theorem 4.3. For algorithm (4.1)–(4.5), the following estimate holds true


Z
(n+1)

(n−1)
(P,t) and r = c
∗
/(c
∗
+ τ
−1
).
22 Monotone domain decomposition algorithms
Proof. The proof of the theorem is similar to the proof of Theorem 3.5 and based on the
maximum pr inciple in Lemma 2.1 and the estimate (2.8)withβ
= 1forthedifference
operator ᏸ
hτ
. 
Remark 4.4. In similar fashion to the proof of Theorem 3.5,theproofofTheorem 4.3
includes the result
r = r for the undecomposed algorithm.
Without loss of generality, we assume that for the parabolic problem (1.2), the bound-
ary condition g(P,t)
= 0. This assumption can always be obtained via a change of vari-
ables. Let on each t ime level the initial function V
(0)
(P,t) be chosen in the form of (4.11)
with R(P,t)
= 0, that is, V
(0)
(P,t) is the solution of the following difference problem
ᏸ
hτ

iterates n
∗
satisfies n
∗
≥ 2. Then the following estimate on convergence rate holds
max
t
k
∈ω
τ


V

t
k

−
U

t
k



≤
D

c
∗

(n)
(P,t)}
converges monotonically.
Proof. The difference problem for V(P,t
k
) = V
(n
∗
)
(P,t
k
)canberepresentedintheform
ᏸ
hτ
V

P,t
k

+ f

P,t
k
,V

−
τ
−1
V


From here, (2.5) and using the mean-value theorem, we get the difference pr oblem for
W(P,t
k
) = V(P,t
k
) − U(P,t
k
)

ᏸ
hτ
+ f
(n
∗
)
u

W

P,t
k

=
G
(n
∗
)

P,t
k

G
(n
∗
)

t
k



ω
h
≤

c
∗
+ η



Z
(n
∗
)

t
k





t
k



ω
h
+


W

t
k−1



ω
h
. (4.19)
By (2.8), from (4.1)–(4.5), we have


Z
(1)

t
k





ω
h
. (4.20)
By the mean-value theorem,
f

P,t
k
,V
(0)

=
f (P,t
k
,0)+ f
(0)
u
V
(0)

P,t
k

, (4.21)
and from (4.20), it follows that



u
V
(0)

t
k



ω
h
+ τ


f

P,t
k
,0

−
τ
−1
V

t
k−1




P,t
k
,0

−
τ
−1
V

t
k−1



ω
h
≤

2τ + c
∗
τ
2



f

P,t
k
,0

W

t
k



ω
h
≤

k

l=1
D
l

τ

c
∗
+ η

(r)
n
∗
−1
, k = 1, ,N
τ
. (4.24)

pendent of μ and τ.Fork
= 2, we have


Z
(1)

t
2



ω
h
≤

2τ + c
∗
τ
2



f

P,t
2
,0



(P,t
1
). As follows from Theorem 4.1, the monotone sequences
{V
(n)
(P,t
1
)} and {V
(n)
(P,t
1
)} are μ-uniformly bounded from above by V
(0)
(P,t
1
)and
from below by V
(0)
(P,t
1
). Applying (2.8)withβ = 1totheproblem(4.14)att = t
1
,we
have


V
(0)
(t
1

2
is independent of μ
and τ. Now by induction on k, we prove that all constants D
k
are independent of μ,and,
hence, constant D
= T max
1≤k≤N
τ
D
k
in (4.15)isindependentofμ and τ.Thus,weprove
the theorem.

Remark 4.6. In the next section, we present sufficientconditionstoguaranteethein-
equality
r<1requiredinTheorem 4.5.
4.4. Uniform convergence of the monotone domain decomposition algorithm (4.1)–
(4.5). Here we analyze a convergence rate of algorithm (4.1)–(4.5)appliedtothediffer-
ence scheme (2.5) defined on the piecewise uniform mesh (3.80).
The difference scheme (2.5) on the piecewise uniform mesh (3.80)convergesμ-uni-
formly to the solution of (1.2):
max
t
k
∈ω
τ


U

, (4.27)
where constant C is independent of μ, τ and N. The proof of this result can be found in
[7].
Similar to Theorem 3.9, we have the following uniform convergence property of algo-
rithm (4.1)–(4.5).
Theorem 4.7. Let the interfacial subdomains
θ
h
m
, m = 1, ,M − 1 and ϑ
h
l
, l = 1, ,L − 1
be located in the x-andy-directions, respectively, outside the boundary layers (unbalanced
decomposition). Suppose that μ
≤ μ
0
 1, and that the following conditions are satisfied
N ≤
1
μ
0
, N = max

N
x
,N
y

, τ<

, as assumed in the theorem, then r<1. From here, (4.15)and
(4.27), we conclude
max
t
k
∈ω
τ


V
(n)

t
k

−
u

t
k



ω
h
≤ C

N
−1
lnN

D are independent of μ, τ and N.Weprovethetheorem. 
Remark 4.8. The implicit two-level difference scheme (2.5) is of first order with respect
to τ.Since

Q = ᏻ(τ), one may choose n
∗
= 2 to keep the global er ror of algorithm (4.1)–
(4.5) consistent with the global error of the difference scheme (2.5).
5. Numerical experiments
Now the monotone domain decomposition algorithms (3.5)–(3.10)and(4.1)–(4.5)are
respectively applied to reaction-diffusion problems of elliptic and parabolic types. All ex-
periments are performed on a serial computer equipped with a 2.8 GHz Pentium 4 pro-
cessor. Some consequences for parallel implementation are also discussed. We consider
in turn the elliptic problem
−μ
2

u
xx
+ u
yy

+
u
− 4
5 − u
= 0, (x, y) ∈ ω ={0 <x<1}×{0 <y<1},
u(x, y)
= 1, (x, y) ∈ ∂ω,
(5.1)

= 4. For μ  1theproblem
is singularly perturbed and the solution increases sharply from u
= 1on∂ω to u = 4on
the interior. The solution to the parabolic problem approaches this steady state with time.
For the continuous problems (5.1)and(5.2), we solve the corresponding nonlinear
difference schemes (2.2)and(2.5) with the monotone domain decomposition algorithms
(3.5)–(3.10)and(4.1)–(4.5), respectively. We employ a piecewise uniform mesh (3.80)
and suppose that N
x
= N
y
= N. Because the mesh is only piecewise uniform, the linear
system arising from the difference problem on a given subdomain may be nonsymmetric.
Therefore, we solve all linear systems with the restarted GMRES algorithm from [10],
suitable for nonsymmetric systems.
5.1. The elliptic problem. Define upper and lower solutions
V
(0)
and V
(0)
by V
(0)
(ω
h
) =
4, V
(0)
(∂ω
h
) = 1andV