# A New Way to Generate an Exponential Finite Difference Scheme for 2D Convection-Diffusion Equations

Authors: Caihua Wang
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

The idea of direction changing and order reducing is proposed to generate an exponential
difference scheme over a five-point stencil for solving two-dimensional (2D) convection-diffusion
equation with source term. During the derivation process, the higher order derivatives alongy-direction are removed to the derivatives along x-direction iteratively
using information given by the original differential equation (similarly from x-direction to
y-direction) and then instead of keeping finite terms in the Taylor series expansion, infinite
terms which constitute convergent series are kept on deriving the exponential coefficients
of the scheme. From the construction process one may gain more insight into the relations
among the stencil coefficients. The scheme is of positive type so it is unconditionally
stable and the convergence rate is proved to be of second-order. Fourth-order accuracy
can be obtained by applying Richardson extrapolation algorithm. Numerical results
show that the scheme is accurate, stable, and especially suitable for convection-dominated
problems with different kinds of boundary layers including elliptic and parabolic ones.
The idea of the method can be applied to a wide variety of differential equations.


---

## Body

## 1. Introduction

“Singular perturbation” means that a small perturbation may cause a large impact in mathematical or physical problems. This terminology was first used by Friedrichs and Wasow in their paper [1] in 1946. Despite this fairly long history, the subject of singular perturbations is not a settled one and there are still a lot of open problems to be investigated. Some surveys on the computational techniques for different kinds of singularly perturbed problems can be seen in [2–4].Among abundant mathematical models for singular perturbation problems, the convection-diffusion problems play a very important role as they arise in fluid flows, groundwater flows, reactive flows, traffic flows, and so forth. Stynes [5] presented a brief history of the development of numerical methods for steady-state convection-diffusion problems, and he also extensively discussed some common numerical methods in [6].In this paper, the 2D linear convection-diffusion model is considered:(1)
L
u
(
x
)
≡
-
u
x
x
-
u
y
y
+
p
(
x
,
y
)
u
x
+
q
(
x
,
y
)
u
y
=
f
(
x
,
y
)
,
for transport variable u defined in a bounded domain Ω
=
(
0,1
)
×
(
0,1
), with appropriate boundary condition:
(2)
u
(
x
,
y
)
|
Γ
=
γ
(
x
,
y
)
.
The second-order derivatives in (1) model diffusion while the first-order derivatives are associated with convective or transport process; p
(
x
,
y
)
,
q
(
x
,
y
), and f
(
x
,
y
) are sufficiently smooth functions in Ω. When the magnitudes of p
(
x
,
y
) and q
(
x
,
y
) are large (convection-dominated case), boundary and interior layers will normally emerge in the solutions of such kind of problems. The term “layers” means some thin regions in which the solution fluctuates rapidly.It is well known that global unphysical oscillations may occur if standard discretization schemes on general meshes are used [7, 8]. Hence, stabilized methods and/or a priori adapted meshes are widely developed in order to get discrete solutions with satisfactory accuracy. In this paper, we mainly focus on special finite difference (FD) methods on uniform meshes.The defect-correction method, as investigated in [9–11], made the early attempt to combine the accuracy of the central difference (CD) method and the stability of the upwind difference (UD) method. But according to Segal’s report [12], the method is not useful for convection-dominated problems because of its slow convergence.In recent years, high order compact FD methods have aroused renewed interest and a variety of techniques has been developed. A compact difference scheme is one that is restricted to cells surrounding any given node and does not extend further, so it is convenient for use since no special techniques are needed for points near the boundary,.A lot of compact FD methods may be classified as polynomial FD schemes for the influence coefficients are connected to polynomials functions of the coefficients of the differential operator; interested readers are referred to [13–20] and the references therein for more details. Numerical experiments have showed that some high order compact polynomial FD schemes can yield high accuracy results. But each scheme has its range of application, just as Stynes have stated in [5], some polynomial FD schemes are difficult to develop due to the need for extensive algebraic manipulation, and the local mesh size should be small enough to make sure that the coefficients matrix satisfies the maximum principle.Another kind of compact FD methods is the exponential ones that the coefficients of the schemes are connected to exponential functions of the coefficients of the differential operator. The exponential FD scheme was first introduced by Allen and Southwell [21] to solve the second-order partial differential equation governing the transport of vorticity. Later Il’in [22] derived in principle the same scheme. But both [21] and [22] were investigated only for 1D case. MacKinnon and Johnson [23] derived a fourth-order exponential FD scheme for one-dimensional convection-diffusion equation and extended the formulation to two dimensions. Chen et al. [24] developed a perturbational fourth-order exponential FD scheme with diagonally dominant coefficient matrix for the convection-diffusion equations based on a second-order exponential FD scheme. The authors in [25, 26] also developed a class of fourth-order compact exponential FD schemes for solving 1D and 2D convection-diffusion problems, respectively. Here we should point out that the coefficients of the high-order exponential schemes in [23–26] involve both exponential and polynomial functions so they are still conditionally stable. Most references lack theoretical analysis of the stability and convergence of the scheme.We notice that coefficients of most of those exponential type FD schemes [23–26] have close contact with a famous scheme Il’in [22] for 1D convection-diffusion problems. Roos [27] described ten different approaches to generate the Il’in and related schemes for the 1D problems, including the compact exponentially fitted method, collocation method, finite-volume method, and exponential Petrov-Galerkin finite-element method. But Roos also pointed out that “there are still a lot of open questions and technical difficulties” as to the generalization of most of those methods to the 2D case.As is known to all, the remainder reapproximation technique is an efficient procedure to increase the accuracy of approximations for many problems in numerical analysis. More specifically, a basic scheme (e.g., CDS) is given first and the truncation error is analyzed by using Taylor series expansions and then the high-order derivatives that appear in the leading term of the truncation error are reapproximated by using information given by the original differential equation. The details of the idea can be seen in [13, 14, 19, 23, 25, 26].Just inspired by the remainder reapproximation method, we aim to provide completely new way to generate an exponential FD scheme for 2D convection-diffusion problem directly without the aid of schemes for 1D case. The stencil coefficients of the FD scheme are to be determined on a five-point stencil. Different from the previous methods in the literature, the main idea is direction changing and order reducing. We keep as many high-order derivatives in the Taylor series expansions as we can at the very beginning of the construction of the scheme. During the derivation process the higher order derivatives alongy-direction are removed to the derivatives along x-direction (or from x-direction to y-direction) iteratively using information given by the original differential equation (direction changing) and result in terms with the orders of y-directional (or x-directional) derivatives less than second-order (order reducing). Then infinite terms which constitute convergent series are kept to form a system of linear equations to determine the final stencil coefficients of the scheme.Here I should emphasize that the meaning of the paper lies not merely on the scheme over a five-point stencil but also on the broad prospect of the idea in application: it can be used to construct the scheme of high-order convergence over a nine-point stencil, it can also be used to construct schemes on nonuniform meshes and to solve many other different kinds of problems.The outline of the rest of the paper is as follows. In Section2, a detailed derivation of a five-point exponential difference scheme is included for the convection-diffusion problem (1) together with a thorough discussion. Proofs of stability and convergence of the scheme are given in Section 3. Firstly, we show that our scheme is of positive type and then based on estimates of the remainder stability and convergence of the difference scheme are analyzed using barrier functions and the maximum principle for elliptic differential equation. The error is proved to be of second-order accuracy. In Section 4, different kinds of problems are included for numerical experiments, fourth-order accuracy can be obtained by Richardson extrapolation algorithm. Results show that the method can cope with problems with and without boundary layers and especially fit for convection-dominated ones. Some concluding remarks are summarized in Section 5.


## 2. A New Way to Generate Exponential FD Scheme and Some Discussions

### 2.1. Derivation of the FD Scheme

LetM and N be two positive integers. Partition the region Ω by the uniform rectangular grid of points with a spacing h in the x-direction and k in the y-direction, where h
=
1
/
M and k
=
1
/
N. Denote all the interior mesh points in Ω by Ω
h
,
k and the mesh points on the boundary by Γ
h
,
k.The value of a functionu
(
x
,
y
) at a reference mesh point (
x
,
y
) is denoted by u
0, and those values at its four immediate neighboring mesh points are denoted by u
i
,
i
=
1,2
,
3,4. The discretized values p
i
,
q
i
,
f
i, which will be used in the following part of the paper, have their obvious meanings too. Before we go on, we assume p
(
x
,
y
)
>
0 and q
(
x
,
y
)
>
0 firstly just for brevity, though we may see at the end of the section that this is not the necessary condition for our method.The compact five-grid point stencil is shown in Figure1.Figure 1
2D five-point stencil.The method will be based on a five-point approximation using Taylor series expansion at the reference mesh point(
x
,
y
) numbered 0 in Figure 1. For simplicity, we introduce denotation s
(
i
,
j
)
=
s
(
i
,
j
)
(
x
,
y
)
=
∂
(
i
+
j
)
s
(
x
,
y
)
/
∂
x
i
∂
y
j for any smooth function s
(
x
,
y
). Assume that the exact solution to differential equation (1) satisfies the following expression:
(3)
ω
0
u
0
-
(
ω
1
u
1
+
ω
2
u
2
+
ω
3
u
3
+
ω
4
u
4
)
=
f
0
+
γ
01
f
0
(
0,1
)
+
γ
10
f
0
(
1,0
)
+
R
0
,
where ω
i
(
i
=
0,1
,
2,3
,
4
) and γ
01
,
γ
10 are coefficients to be determined later and R
0 is the term connected with the local truncation error that also will be shown later. Sometimes we may also omit the subscript 0 and write u
0
=
u
,
p
0
=
p
,
q
0
=
q
,
f
0
=
f, and so forth, just for brevity.Set(4)
v
=
ω
0
u
0
-
(
ω
1
u
1
+
ω
2
u
2
+
ω
3
u
3
+
ω
4
u
4
)
.
u
1
,
u
2
,
u
3
,
u
4 can be defined at the reference mesh point 0, respectively, using Taylor’s theorem:
(5)
u
1
=
u
+
h
u
(
1,0
)
+
1
2
h
2
u
(
2,0
)
+
1
3
!
h
3
u
(
3,0
)
+
1
4
!
h
4
u
0
(
4,0
)
+
1
5
!
h
5
u
0
(
5,0
)
+
⋯
,
u
3
=
u
-
h
u
(
1,0
)
+
1
2
h
2
u
(
2,0
)
-
1
3
!
h
3
u
(
3,0
)
+
1
4
!
h
4
u
(
4,0
)
-
1
5
!
h
5
u
(
5,0
)
+
⋯
,
(6)
u
2
=
u
+
k
u
(
0,1
)
+
1
2
k
2
u
(
0,2
)
+
1
3
!
k
3
u
(
0,3
)
+
1
4
!
k
4
u
(
0,4
)
+
1
5
!
k
5
u
(
0,5
)
+
⋯
,
(7)
u
4
=
u
-
k
u
(
0,1
)
+
1
2
k
2
u
(
0,2
)
-
1
3
!
k
3
u
(
0,3
)
+
1
4
!
k
4
u
(
0,4
)
-
1
5
!
k
5
u
(
0,5
)
+
⋯
.
From the differential equation (1) one can have
(8)
u
y
y
=
-
u
x
x
+
p
(
x
,
y
)
u
x
+
q
(
x
,
y
)
u
y
-
f
(
x
,
y
)
.
In a next step we freeze the data and assume that the convection terms p
(
x
,
y
)
,
q
(
x
,
y
) are constant in the neighbouring meshes of grid 0; that is, p
(
x
,
y
)
=
p
0
,
q
(
x
,
y
)
=
q
0, and then using (8) as an auxiliary relation, we may remove all the higher order derivatives (equal to or more than second) of function u
(
x
,
y
) along y-direction to the derivatives of u along x-direction; for instance,
(9)
u
(
0,3
)
=
q
2
u
(
0,1
)
+
p
q
u
(
1,0
)
+
p
u
(
1,1
)
-
q
u
(
2,0
)
-
u
(
2,1
)
-
q
f
-
f
(
0,1
)
.
Equation (9) is really the key to learn the idea of the construction of the differencing scheme. We have removed the higher order (equal to 3) derivatives of function u
(
x
,
y
) along y-direction to the derivatives of u along x-direction (direction changing), and result in the right side of (9) with the orders of y-directional derivatives of function u
(
x
,
y
) less than 2 (order reducing). u
(
0,4
)
,
u
(
0,5
)
… can be treated in the same way. Generally,
(10)
u
(
0
,
n
)
=
∂
n
-
2
∂
y
n
-
2
(
-
u
x
x
+
p
u
x
+
q
u
y
-
f
)
,
for
n
≥
2
,
substituting (10) into (6) iteratively until all the derivatives of direction-y are less than second-order. The idea of direction changing and order reducing is just reflected in this procedure which can be carried out easily by means of Mathematica software, and then we obtain
(11)
u
2
=
u
+
(
k
+
1
2
q
k
2
+
1
3
!
q
2
k
3
+
1
4
!
q
3
k
4
+
⋯
)
u
(
0,1
)
+
(
1
2
k
2
+
1
3
!
q
k
3
+
1
4
!
q
2
k
4
+
1
5
!
q
3
k
5
+
⋯
)
p
u
(
1,0
)
+
(
1
3
!
k
3
+
2
4
!
q
k
4
+
3
5
!
q
2
k
5
+
4
6
!
q
3
k
6
+
⋯
)
p
u
(
1,1
)
+
⋯
+
[
-
(
1
2
k
2
+
1
3
!
q
k
3
+
1
4
!
q
2
k
4
+
1
5
!
q
3
k
5
+
⋯
)
]
f
+
[
-
(
1
3
!
k
3
+
1
4
!
q
k
4
+
1
5
!
q
2
k
5
+
1
6
!
q
3
k
6
+
⋯
)
]
f
(
0,1
)
+
⋯
.
Noticing that the infinite terms in the brackets of (11) constitute convergent series, we rewrite it as
(12)
u
2
=
u
+
α
01
u
(
0,1
)
+
α
10
u
(
1,0
)
+
α
11
u
(
1,1
)
+
⋯
+
ξ
00
f
+
ξ
01
f
(
0,1
)
+
⋯
,
with
(13)
α
01
=
∑
i
=
1
∞
‍
k
i
q
i
-
1
i
!
=
exp
⁡
⁡
(
q
k
)
-
1
q
,
α
10
=
p
∑
i
=
2
∞
‍
k
i
q
i
-
2
i
!
=
p
exp
⁡
⁡
(
q
k
)
-
q
k
-
1
q
2
,
α
11
=
p
∑
i
=
3
∞
‍
(
i
-
2
)
k
i
q
i
-
3
i
!
=
p
exp
⁡
(
q
k
)
(
q
k
-
2
)
+
q
k
+
2
q
3
,
ξ
00
=
-
∑
i
=
2
∞
‍
k
i
q
i
-
2
i
!
=
q
k
+
1
-
exp
⁡
(
q
k
)
q
2
,
ξ
01
=
-
∑
i
=
3
∞
‍
k
i
q
i
-
3
i
!
=
q
2
k
2
+
2
q
k
+
2
-
2
exp
⁡
(
q
k
)
2
q
3
.
Similarly, we get
(14)
u
4
=
u
+
β
01
u
(
0,1
)
+
β
10
u
(
1,0
)
+
β
11
u
(
1,1
)
+
⋯
+
η
00
f
+
η
01
f
(
0,1
)
+
⋯
,
with
(15)
β
01
=
∑
i
=
1
∞
‍
(
-
k
)
i
q
i
-
1
i
!
=
-
(
1
-
exp
⁡
(
-
q
k
)
)
q
,
β
10
=
p
∑
i
=
2
∞
‍
(
-
k
)
i
q
i
-
2
i
!
=
p
q
k
-
1
+
exp
⁡
(
-
q
k
)
q
2
,
β
11
=
p
∑
i
=
3
∞
‍
(
i
-
2
)
(
-
k
)
i
q
i
-
3
i
!
w
w
=
-
p
(
q
k
-
2
)
+
(
q
k
+
2
)
exp
⁡
(
-
q
k
)
q
3
,
η
00
=
-
∑
i
=
2
∞
‍
(
-
k
)
i
q
i
-
2
i
!
=
-
q
k
-
1
+
exp
⁡
(
-
q
k
)
q
2
,
η
01
=
-
∑
i
=
3
∞
‍
(
-
k
)
i
q
i
-
3
i
!
=
q
2
k
2
-
2
q
k
+
2
-
2
exp
⁡
(
-
q
k
)
2
q
3
.
Substituting (5), (12), and (14) into (4) yields
(16)
v
=
(
ω
0
-
ω
1
-
ω
2
-
ω
3
-
ω
4
)
u
+
(
-
α
01
ω
2
-
β
01
ω
4
)
u
(
0,1
)
+
(
-
h
ω
1
+
h
ω
3
-
α
10
ω
2
-
β
10
ω
4
)
u
(
1,0
)
+
(
-
α
11
ω
2
-
β
11
ω
4
)
u
(
1,1
)
+
⋯
+
(
-
ξ
00
ω
2
-
η
00
ω
4
)
f
+
(
-
ξ
01
ω
2
-
η
01
ω
4
)
f
(
0,1
)
+
⋯
.
Denoting v
=
σ
00
u
+
σ
01
u
(
0,1
)
+
σ
10
u
(
1,0
)
+
σ
11
u
(
1,1
)
+
⋯
+
μ
00
f
+
μ
01
f
(
0,1
)
+
⋯ and comparing it with (3), we know clearly that it is necessary to set σ
00
=
0
,
σ
01
=
0
,
σ
10
=
0, and μ
00
=
1
,
γ
01
=
μ
01 in order to get the numerical scheme of (1), and after some easy calculation, we get five equations:
(17)
ω
0
-
ω
1
-
ω
2
-
ω
3
-
ω
4
=
0
,
(18)
(
1
-
exp
⁡
(
-
q
k
)
)
ω
4
q
-
(
exp
⁡
(
q
k
)
-
1
)
ω
2
q
=
0
,
-
h
ω
1
+
h
ω
3
-
p
(
exp
⁡
(
q
k
)
-
q
k
-
1
)
ω
2
q
2
-
p
(
q
k
-
1
+
exp
⁡
(
-
q
k
)
)
ω
4
q
2
=
0
,
(19)
(
exp
⁡
(
q
k
)
-
q
k
-
1
)
ω
2
q
2
+
(
q
k
-
1
+
exp
⁡
(
-
q
k
)
)
ω
4
q
2
=
1
,
(20)
γ
01
=
(
2
exp
⁡
(
q
k
)
-
2
-
2
q
k
-
q
2
k
2
)
ω
2
2
q
3
-
(
q
2
k
2
-
2
q
k
+
2
-
2
exp
⁡
⁡
(
-
q
k
)
)
ω
4
2
q
3
.Instead of (8), if we change differential equation (3) into the form
(21)
u
x
x
=
-
u
y
y
+
p
(
x
,
y
)
u
x
+
q
(
x
,
y
)
u
y
-
f
(
x
,
y
)
,
similar procedure can be carried out as above, (17) will be obtained together with the following formulae:
(22)
(
1
-
exp
⁡
(
-
p
h
)
)
ω
3
p
-
(
exp
⁡
(
p
h
)
-
1
)
ω
1
p
=
0
,
(23)
-
k
ω
2
+
k
ω
4
-
q
(
exp
⁡
(
p
h
)
-
p
h
-
1
)
ω
1
p
2
-
q
(
p
h
-
1
+
exp
⁡
(
-
p
h
)
)
ω
3
p
2
=
0
,
(
exp
⁡
(
p
h
)
-
p
h
-
1
)
ω
1
p
2
+
(
p
h
-
1
+
exp
⁡
(
-
p
h
)
)
ω
3
p
2
=
1
,
(24)
γ
10
=
(
2
exp
⁡
(
p
h
)
-
2
-
2
p
h
-
p
2
h
2
)
ω
1
2
p
3
-
(
p
2
h
2
-
2
p
h
+
2
-
2
exp
⁡
(
-
p
h
)
)
ω
3
2
p
3
.These seven equations, (17)–(19) and (22)-(23) which contain five unknown values ω
i
,
i
=
0,1
,
2,3
,
4, are compatible; in fact (17)–(19) and (22) can be easily solved by means of Mathematica software to give
(25)
ω
1
=
exp
⁡
(
-
p
h
)
p
(
1
-
exp
⁡
(
-
p
h
)
)
h
,
(26)
ω
2
=
exp
⁡
(
-
q
k
)
q
(
1
-
exp
⁡
(
-
q
k
)
)
k
,
(27)
ω
3
=
p
(
1
-
exp
⁡
(
-
p
h
)
)
h
,
(28)
ω
4
=
q
(
1
-
exp
⁡
(
-
q
k
)
)
k
,
(29)
ω
0
=
p
(
1
+
exp
⁡
(
-
p
h
)
)
h
(
1
-
exp
⁡
(
-
p
h
)
)
+
q
(
1
+
exp
⁡
(
-
q
k
)
)
k
(
1
-
exp
⁡
(
-
q
k
)
)
,
and then (23) holds immediately. In fact, σ
11
=
0 holds too. Substituting (25)–(29) into (20) and (24), we get, respectively,
(30)
γ
01
=
-
(
2
+
q
k
)
exp
⁡
⁡
(
-
q
k
)
+
q
k
-
2
2
q
(
1
-
exp
⁡
⁡
(
-
q
k
)
)
,
(31)
γ
10
=
-
(
2
+
p
h
)
exp
⁡
⁡
(
-
p
h
)
+
p
h
-
2
2
p
(
1
-
exp
⁡
⁡
(
-
p
h
)
)
.
We denote the numerical solution at stencil points by U
i
(
i
=
0,1
,
2,3
,
4
); the FD approximation can then be obtained by dropping the truncation error R
0 from (3)
(32)
L
h
U
0
≡
ω
0
U
0
-
(
ω
1
U
1
+
ω
2
U
2
+
ω
3
U
3
+
ω
4
U
4
)
=
f
0
+
γ
01
f
0
(
0,1
)
+
γ
10
f
0
(
1,0
)
,
with ω
i
(
i
=
0,1
,
2,3
,
4
) defined in (25)–(29) and γ
01
,
γ
10 given in (30) and (31). And the boundary points values are directly given by the boundary condition (2) as
(33)
U
j
=
u
(
x
j
,
y
j
)
,
for
(
x
j
,
y
j
)
∈
Γ
h
,
k
.Numerical methods like these, whose coefficients involve exponential functions, are known collectively as compact exponentially difference schemes. For brevity we call the scheme EDS (exponential difference scheme) in the sequel. If we freeze the data off
(
x
,
y
) to be constant locally too, then in (32) we should let γ
01
=
0 and γ
10
=
0; we will call the scheme EDS0 as below.


### 2.2. Discussions of EDS and EDS0

(
1
) It can be easily seen from (25)–(29) that the coefficients matrix of difference operator L
h is five-diagonal and of positive type: the matrix is diagonally dominant with positive diagonal elements and negative off-diagonal elements. And since upwind effect is preserved in the scheme so hopefully it will serve well for convection-dominated problems.(
2
) The scheme EDS0 is equivalent to central differencing applied to a modified differential equation:
(34)
-
δ
(
p
h
2
)
u
x
x
-
δ
(
q
k
2
)
u
y
y
+
p
u
x
+
q
u
y
=
f
,
where δ
(
v
)
=
v
·
coth
⁡
⁡
v
. Since
(35)
δ
(
p
h
2
)
=
1
+
p
2
h
2
12
+
O
(
h
4
)
,
δ
(
q
k
2
)
=
1
+
q
2
k
2
12
+
O
(
k
4
)
,
the scheme EDS0 can be said to have artificial diffusion. Likewise, the scheme EDS is equivalent to central differencing applied to the following modified differential equation:
(36)
-
δ
(
p
h
2
)
u
x
x
-
δ
(
q
k
2
)
u
y
y
+
p
u
x
+
q
u
y
=
f
+
γ
01
f
(
0,1
)
+
γ
10
f
(
1,0
)
.(3) If we implement Taylor expansion forω
i
(
i
=
0,1
,
2,3
,
4
) with only two terms to be reserved, that is, ω
1
≈
(
1
/
h
2
)
-
(
p
/
2
h
)
,
ω
2
≈
(
1
/
k
2
)
-
(
q
/
2
k
)
,
ω
3
≈
(
1
/
h
2
)
+
(
p
/
2
h
)
,
ω
4
≈
(
1
/
k
2
)
+
(
q
/
2
k
), and ω
0
≈
(
2
/
h
2
)
+
(
2
/
k
2
), together with γ
01
=
0 and γ
10
=
0, then EDS (32) becomes CDS.(4) If consideringp
(
x
0
,
y
0
)
>
0 and q
(
x
0
,
y
0
)
=
0 on some reference point (
x
0
,
y
0
)
∈
Ω, we have no difficulty to follow the whole derivation process from the beginning of this section to arrive at the difference operator similar to L
h defined in (32) but with ω
2
=
1
/
k
2
,
ω
4
=
1
/
k
2
,
γ
01
=
0, and with ω
1
,
ω
3
,
γ
10 no change in (25), (27), and (31). Correspondingly, ω
0
=
(
p
(
1
+
exp
⁡
⁡
(
-
p
h
)
)
)
/
(
h
(
1
-
exp
⁡
⁡
(
-
p
h
)
)
)
+
(
2
/
k
2
). Easily we can see that EDS and EDS0 are still of positive type.Furthermore, ifp
(
x
,
y
)
>
0 and q
(
x
,
y
)
≡
0 are satisfied in Ω, this problem typically has an exponential boundary layer at x
=
1 and parabolic boundary layers at y
=
0 and y
=
1 (see [5, 28]). The parabolic layers raise interesting numerical issues. They cause numerical instabilities that are far less serious than those engendered by exponential layers; yet it is difficult for traditional methods to approximate them accurately. But EDS and EDS0 can cope with such kind of problems without difficulty, and we will provide numerical Example 10 in Section 4 to validate it.As forp
(
x
0
,
y
0
)
=
0 at some points (or p
(
x
,
y
)
≡
0) and q
(
x
0
,
y
0
)
>
0, similar formulae can be derived.(5) Ifp
(
x
,
y
)
≡
0 and q
(
x
,
y
)
≡
0 hold in Ω, EDS becomes CDS for Poisson’s equation.(6) We assumep
(
x
,
y
)
>
0 and q
(
x
,
y
)
>
0 at the beginning of the section just for brevity, as for p
(
x
,
y
)
<
0 or q
(
x
,
y
)
<
0, the difference scheme (32) still holds.


## 2.1. Derivation of the FD Scheme

LetM and N be two positive integers. Partition the region Ω by the uniform rectangular grid of points with a spacing h in the x-direction and k in the y-direction, where h
=
1
/
M and k
=
1
/
N. Denote all the interior mesh points in Ω by Ω
h
,
k and the mesh points on the boundary by Γ
h
,
k.The value of a functionu
(
x
,
y
) at a reference mesh point (
x
,
y
) is denoted by u
0, and those values at its four immediate neighboring mesh points are denoted by u
i
,
i
=
1,2
,
3,4. The discretized values p
i
,
q
i
,
f
i, which will be used in the following part of the paper, have their obvious meanings too. Before we go on, we assume p
(
x
,
y
)
>
0 and q
(
x
,
y
)
>
0 firstly just for brevity, though we may see at the end of the section that this is not the necessary condition for our method.The compact five-grid point stencil is shown in Figure1.Figure 1
2D five-point stencil.The method will be based on a five-point approximation using Taylor series expansion at the reference mesh point(
x
,
y
) numbered 0 in Figure 1. For simplicity, we introduce denotation s
(
i
,
j
)
=
s
(
i
,
j
)
(
x
,
y
)
=
∂
(
i
+
j
)
s
(
x
,
y
)
/
∂
x
i
∂
y
j for any smooth function s
(
x
,
y
). Assume that the exact solution to differential equation (1) satisfies the following expression:
(3)
ω
0
u
0
-
(
ω
1
u
1
+
ω
2
u
2
+
ω
3
u
3
+
ω
4
u
4
)
=
f
0
+
γ
01
f
0
(
0,1
)
+
γ
10
f
0
(
1,0
)
+
R
0
,
where ω
i
(
i
=
0,1
,
2,3
,
4
) and γ
01
,
γ
10 are coefficients to be determined later and R
0 is the term connected with the local truncation error that also will be shown later. Sometimes we may also omit the subscript 0 and write u
0
=
u
,
p
0
=
p
,
q
0
=
q
,
f
0
=
f, and so forth, just for brevity.Set(4)
v
=
ω
0
u
0
-
(
ω
1
u
1
+
ω
2
u
2
+
ω
3
u
3
+
ω
4
u
4
)
.
u
1
,
u
2
,
u
3
,
u
4 can be defined at the reference mesh point 0, respectively, using Taylor’s theorem:
(5)
u
1
=
u
+
h
u
(
1,0
)
+
1
2
h
2
u
(
2,0
)
+
1
3
!
h
3
u
(
3,0
)
+
1
4
!
h
4
u
0
(
4,0
)
+
1
5
!
h
5
u
0
(
5,0
)
+
⋯
,
u
3
=
u
-
h
u
(
1,0
)
+
1
2
h
2
u
(
2,0
)
-
1
3
!
h
3
u
(
3,0
)
+
1
4
!
h
4
u
(
4,0
)
-
1
5
!
h
5
u
(
5,0
)
+
⋯
,
(6)
u
2
=
u
+
k
u
(
0,1
)
+
1
2
k
2
u
(
0,2
)
+
1
3
!
k
3
u
(
0,3
)
+
1
4
!
k
4
u
(
0,4
)
+
1
5
!
k
5
u
(
0,5
)
+
⋯
,
(7)
u
4
=
u
-
k
u
(
0,1
)
+
1
2
k
2
u
(
0,2
)
-
1
3
!
k
3
u
(
0,3
)
+
1
4
!
k
4
u
(
0,4
)
-
1
5
!
k
5
u
(
0,5
)
+
⋯
.
From the differential equation (1) one can have
(8)
u
y
y
=
-
u
x
x
+
p
(
x
,
y
)
u
x
+
q
(
x
,
y
)
u
y
-
f
(
x
,
y
)
.
In a next step we freeze the data and assume that the convection terms p
(
x
,
y
)
,
q
(
x
,
y
) are constant in the neighbouring meshes of grid 0; that is, p
(
x
,
y
)
=
p
0
,
q
(
x
,
y
)
=
q
0, and then using (8) as an auxiliary relation, we may remove all the higher order derivatives (equal to or more than second) of function u
(
x
,
y
) along y-direction to the derivatives of u along x-direction; for instance,
(9)
u
(
0,3
)
=
q
2
u
(
0,1
)
+
p
q
u
(
1,0
)
+
p
u
(
1,1
)
-
q
u
(
2,0
)
-
u
(
2,1
)
-
q
f
-
f
(
0,1
)
.
Equation (9) is really the key to learn the idea of the construction of the differencing scheme. We have removed the higher order (equal to 3) derivatives of function u
(
x
,
y
) along y-direction to the derivatives of u along x-direction (direction changing), and result in the right side of (9) with the orders of y-directional derivatives of function u
(
x
,
y
) less than 2 (order reducing). u
(
0,4
)
,
u
(
0,5
)
… can be treated in the same way. Generally,
(10)
u
(
0
,
n
)
=
∂
n
-
2
∂
y
n
-
2
(
-
u
x
x
+
p
u
x
+
q
u
y
-
f
)
,
for
n
≥
2
,
substituting (10) into (6) iteratively until all the derivatives of direction-y are less than second-order. The idea of direction changing and order reducing is just reflected in this procedure which can be carried out easily by means of Mathematica software, and then we obtain
(11)
u
2
=
u
+
(
k
+
1
2
q
k
2
+
1
3
!
q
2
k
3
+
1
4
!
q
3
k
4
+
⋯
)
u
(
0,1
)
+
(
1
2
k
2
+
1
3
!
q
k
3
+
1
4
!
q
2
k
4
+
1
5
!
q
3
k
5
+
⋯
)
p
u
(
1,0
)
+
(
1
3
!
k
3
+
2
4
!
q
k
4
+
3
5
!
q
2
k
5
+
4
6
!
q
3
k
6
+
⋯
)
p
u
(
1,1
)
+
⋯
+
[
-
(
1
2
k
2
+
1
3
!
q
k
3
+
1
4
!
q
2
k
4
+
1
5
!
q
3
k
5
+
⋯
)
]
f
+
[
-
(
1
3
!
k
3
+
1
4
!
q
k
4
+
1
5
!
q
2
k
5
+
1
6
!
q
3
k
6
+
⋯
)
]
f
(
0,1
)
+
⋯
.
Noticing that the infinite terms in the brackets of (11) constitute convergent series, we rewrite it as
(12)
u
2
=
u
+
α
01
u
(
0,1
)
+
α
10
u
(
1,0
)
+
α
11
u
(
1,1
)
+
⋯
+
ξ
00
f
+
ξ
01
f
(
0,1
)
+
⋯
,
with
(13)
α
01
=
∑
i
=
1
∞
‍
k
i
q
i
-
1
i
!
=
exp
⁡
⁡
(
q
k
)
-
1
q
,
α
10
=
p
∑
i
=
2
∞
‍
k
i
q
i
-
2
i
!
=
p
exp
⁡
⁡
(
q
k
)
-
q
k
-
1
q
2
,
α
11
=
p
∑
i
=
3
∞
‍
(
i
-
2
)
k
i
q
i
-
3
i
!
=
p
exp
⁡
(
q
k
)
(
q
k
-
2
)
+
q
k
+
2
q
3
,
ξ
00
=
-
∑
i
=
2
∞
‍
k
i
q
i
-
2
i
!
=
q
k
+
1
-
exp
⁡
(
q
k
)
q
2
,
ξ
01
=
-
∑
i
=
3
∞
‍
k
i
q
i
-
3
i
!
=
q
2
k
2
+
2
q
k
+
2
-
2
exp
⁡
(
q
k
)
2
q
3
.
Similarly, we get
(14)
u
4
=
u
+
β
01
u
(
0,1
)
+
β
10
u
(
1,0
)
+
β
11
u
(
1,1
)
+
⋯
+
η
00
f
+
η
01
f
(
0,1
)
+
⋯
,
with
(15)
β
01
=
∑
i
=
1
∞
‍
(
-
k
)
i
q
i
-
1
i
!
=
-
(
1
-
exp
⁡
(
-
q
k
)
)
q
,
β
10
=
p
∑
i
=
2
∞
‍
(
-
k
)
i
q
i
-
2
i
!
=
p
q
k
-
1
+
exp
⁡
(
-
q
k
)
q
2
,
β
11
=
p
∑
i
=
3
∞
‍
(
i
-
2
)
(
-
k
)
i
q
i
-
3
i
!
w
w
=
-
p
(
q
k
-
2
)
+
(
q
k
+
2
)
exp
⁡
(
-
q
k
)
q
3
,
η
00
=
-
∑
i
=
2
∞
‍
(
-
k
)
i
q
i
-
2
i
!
=
-
q
k
-
1
+
exp
⁡
(
-
q
k
)
q
2
,
η
01
=
-
∑
i
=
3
∞
‍
(
-
k
)
i
q
i
-
3
i
!
=
q
2
k
2
-
2
q
k
+
2
-
2
exp
⁡
(
-
q
k
)
2
q
3
.
Substituting (5), (12), and (14) into (4) yields
(16)
v
=
(
ω
0
-
ω
1
-
ω
2
-
ω
3
-
ω
4
)
u
+
(
-
α
01
ω
2
-
β
01
ω
4
)
u
(
0,1
)
+
(
-
h
ω
1
+
h
ω
3
-
α
10
ω
2
-
β
10
ω
4
)
u
(
1,0
)
+
(
-
α
11
ω
2
-
β
11
ω
4
)
u
(
1,1
)
+
⋯
+
(
-
ξ
00
ω
2
-
η
00
ω
4
)
f
+
(
-
ξ
01
ω
2
-
η
01
ω
4
)
f
(
0,1
)
+
⋯
.
Denoting v
=
σ
00
u
+
σ
01
u
(
0,1
)
+
σ
10
u
(
1,0
)
+
σ
11
u
(
1,1
)
+
⋯
+
μ
00
f
+
μ
01
f
(
0,1
)
+
⋯ and comparing it with (3), we know clearly that it is necessary to set σ
00
=
0
,
σ
01
=
0
,
σ
10
=
0, and μ
00
=
1
,
γ
01
=
μ
01 in order to get the numerical scheme of (1), and after some easy calculation, we get five equations:
(17)
ω
0
-
ω
1
-
ω
2
-
ω
3
-
ω
4
=
0
,
(18)
(
1
-
exp
⁡
(
-
q
k
)
)
ω
4
q
-
(
exp
⁡
(
q
k
)
-
1
)
ω
2
q
=
0
,
-
h
ω
1
+
h
ω
3
-
p
(
exp
⁡
(
q
k
)
-
q
k
-
1
)
ω
2
q
2
-
p
(
q
k
-
1
+
exp
⁡
(
-
q
k
)
)
ω
4
q
2
=
0
,
(19)
(
exp
⁡
(
q
k
)
-
q
k
-
1
)
ω
2
q
2
+
(
q
k
-
1
+
exp
⁡
(
-
q
k
)
)
ω
4
q
2
=
1
,
(20)
γ
01
=
(
2
exp
⁡
(
q
k
)
-
2
-
2
q
k
-
q
2
k
2
)
ω
2
2
q
3
-
(
q
2
k
2
-
2
q
k
+
2
-
2
exp
⁡
⁡
(
-
q
k
)
)
ω
4
2
q
3
.Instead of (8), if we change differential equation (3) into the form
(21)
u
x
x
=
-
u
y
y
+
p
(
x
,
y
)
u
x
+
q
(
x
,
y
)
u
y
-
f
(
x
,
y
)
,
similar procedure can be carried out as above, (17) will be obtained together with the following formulae:
(22)
(
1
-
exp
⁡
(
-
p
h
)
)
ω
3
p
-
(
exp
⁡
(
p
h
)
-
1
)
ω
1
p
=
0
,
(23)
-
k
ω
2
+
k
ω
4
-
q
(
exp
⁡
(
p
h
)
-
p
h
-
1
)
ω
1
p
2
-
q
(
p
h
-
1
+
exp
⁡
(
-
p
h
)
)
ω
3
p
2
=
0
,
(
exp
⁡
(
p
h
)
-
p
h
-
1
)
ω
1
p
2
+
(
p
h
-
1
+
exp
⁡
(
-
p
h
)
)
ω
3
p
2
=
1
,
(24)
γ
10
=
(
2
exp
⁡
(
p
h
)
-
2
-
2
p
h
-
p
2
h
2
)
ω
1
2
p
3
-
(
p
2
h
2
-
2
p
h
+
2
-
2
exp
⁡
(
-
p
h
)
)
ω
3
2
p
3
.These seven equations, (17)–(19) and (22)-(23) which contain five unknown values ω
i
,
i
=
0,1
,
2,3
,
4, are compatible; in fact (17)–(19) and (22) can be easily solved by means of Mathematica software to give
(25)
ω
1
=
exp
⁡
(
-
p
h
)
p
(
1
-
exp
⁡
(
-
p
h
)
)
h
,
(26)
ω
2
=
exp
⁡
(
-
q
k
)
q
(
1
-
exp
⁡
(
-
q
k
)
)
k
,
(27)
ω
3
=
p
(
1
-
exp
⁡
(
-
p
h
)
)
h
,
(28)
ω
4
=
q
(
1
-
exp
⁡
(
-
q
k
)
)
k
,
(29)
ω
0
=
p
(
1
+
exp
⁡
(
-
p
h
)
)
h
(
1
-
exp
⁡
(
-
p
h
)
)
+
q
(
1
+
exp
⁡
(
-
q
k
)
)
k
(
1
-
exp
⁡
(
-
q
k
)
)
,
and then (23) holds immediately. In fact, σ
11
=
0 holds too. Substituting (25)–(29) into (20) and (24), we get, respectively,
(30)
γ
01
=
-
(
2
+
q
k
)
exp
⁡
⁡
(
-
q
k
)
+
q
k
-
2
2
q
(
1
-
exp
⁡
⁡
(
-
q
k
)
)
,
(31)
γ
10
=
-
(
2
+
p
h
)
exp
⁡
⁡
(
-
p
h
)
+
p
h
-
2
2
p
(
1
-
exp
⁡
⁡
(
-
p
h
)
)
.
We denote the numerical solution at stencil points by U
i
(
i
=
0,1
,
2,3
,
4
); the FD approximation can then be obtained by dropping the truncation error R
0 from (3)
(32)
L
h
U
0
≡
ω
0
U
0
-
(
ω
1
U
1
+
ω
2
U
2
+
ω
3
U
3
+
ω
4
U
4
)
=
f
0
+
γ
01
f
0
(
0,1
)
+
γ
10
f
0
(
1,0
)
,
with ω
i
(
i
=
0,1
,
2,3
,
4
) defined in (25)–(29) and γ
01
,
γ
10 given in (30) and (31). And the boundary points values are directly given by the boundary condition (2) as
(33)
U
j
=
u
(
x
j
,
y
j
)
,
for
(
x
j
,
y
j
)
∈
Γ
h
,
k
.Numerical methods like these, whose coefficients involve exponential functions, are known collectively as compact exponentially difference schemes. For brevity we call the scheme EDS (exponential difference scheme) in the sequel. If we freeze the data off
(
x
,
y
) to be constant locally too, then in (32) we should let γ
01
=
0 and γ
10
=
0; we will call the scheme EDS0 as below.


## 2.2. Discussions of EDS and EDS0

(
1
) It can be easily seen from (25)–(29) that the coefficients matrix of difference operator L
h is five-diagonal and of positive type: the matrix is diagonally dominant with positive diagonal elements and negative off-diagonal elements. And since upwind effect is preserved in the scheme so hopefully it will serve well for convection-dominated problems.(
2
) The scheme EDS0 is equivalent to central differencing applied to a modified differential equation:
(34)
-
δ
(
p
h
2
)
u
x
x
-
δ
(
q
k
2
)
u
y
y
+
p
u
x
+
q
u
y
=
f
,
where δ
(
v
)
=
v
·
coth
⁡
⁡
v
. Since
(35)
δ
(
p
h
2
)
=
1
+
p
2
h
2
12
+
O
(
h
4
)
,
δ
(
q
k
2
)
=
1
+
q
2
k
2
12
+
O
(
k
4
)
,
the scheme EDS0 can be said to have artificial diffusion. Likewise, the scheme EDS is equivalent to central differencing applied to the following modified differential equation:
(36)
-
δ
(
p
h
2
)
u
x
x
-
δ
(
q
k
2
)
u
y
y
+
p
u
x
+
q
u
y
=
f
+
γ
01
f
(
0,1
)
+
γ
10
f
(
1,0
)
.(3) If we implement Taylor expansion forω
i
(
i
=
0,1
,
2,3
,
4
) with only two terms to be reserved, that is, ω
1
≈
(
1
/
h
2
)
-
(
p
/
2
h
)
,
ω
2
≈
(
1
/
k
2
)
-
(
q
/
2
k
)
,
ω
3
≈
(
1
/
h
2
)
+
(
p
/
2
h
)
,
ω
4
≈
(
1
/
k
2
)
+
(
q
/
2
k
), and ω
0
≈
(
2
/
h
2
)
+
(
2
/
k
2
), together with γ
01
=
0 and γ
10
=
0, then EDS (32) becomes CDS.(4) If consideringp
(
x
0
,
y
0
)
>
0 and q
(
x
0
,
y
0
)
=
0 on some reference point (
x
0
,
y
0
)
∈
Ω, we have no difficulty to follow the whole derivation process from the beginning of this section to arrive at the difference operator similar to L
h defined in (32) but with ω
2
=
1
/
k
2
,
ω
4
=
1
/
k
2
,
γ
01
=
0, and with ω
1
,
ω
3
,
γ
10 no change in (25), (27), and (31). Correspondingly, ω
0
=
(
p
(
1
+
exp
⁡
⁡
(
-
p
h
)
)
)
/
(
h
(
1
-
exp
⁡
⁡
(
-
p
h
)
)
)
+
(
2
/
k
2
). Easily we can see that EDS and EDS0 are still of positive type.Furthermore, ifp
(
x
,
y
)
>
0 and q
(
x
,
y
)
≡
0 are satisfied in Ω, this problem typically has an exponential boundary layer at x
=
1 and parabolic boundary layers at y
=
0 and y
=
1 (see [5, 28]). The parabolic layers raise interesting numerical issues. They cause numerical instabilities that are far less serious than those engendered by exponential layers; yet it is difficult for traditional methods to approximate them accurately. But EDS and EDS0 can cope with such kind of problems without difficulty, and we will provide numerical Example 10 in Section 4 to validate it.As forp
(
x
0
,
y
0
)
=
0 at some points (or p
(
x
,
y
)
≡
0) and q
(
x
0
,
y
0
)
>
0, similar formulae can be derived.(5) Ifp
(
x
,
y
)
≡
0 and q
(
x
,
y
)
≡
0 hold in Ω, EDS becomes CDS for Poisson’s equation.(6) We assumep
(
x
,
y
)
>
0 and q
(
x
,
y
)
>
0 at the beginning of the section just for brevity, as for p
(
x
,
y
)
<
0 or q
(
x
,
y
)
<
0, the difference scheme (32) still holds.


## 3. Analysis of Stability and Convergence of the Difference Scheme

### 3.1.  Local Truncation Error

We begin with analyzing the truncation error of the scheme about EDS0. For brevity we only consider the case thatp
(
x
,
y
)
>
0 and q
(
x
,
y
)
>
0; similar analysis can be done for other cases.Estimating the truncation error of EDS0 from (3) yields
(37)
R
=
ω
0
u
0
-
(
ω
1
u
1
+
ω
2
u
2
+
ω
3
u
3
+
ω
4
u
4
)
-
f
.
Substituting Taylor expansions (5)–(7) together with values ω
i
,
(
i
=
0,1
,
2,3
,
4
) given in (25)–(29) into (37), and making some arrangement, we have
(38)
R
=
-
f
+
p
u
(
1,0
)
+
q
u
(
0,1
)
-
p
h
(
1
+
exp
⁡
⁡
(
-
p
h
)
)
2
(
1
-
exp
⁡
⁡
(
-
p
h
)
)
u
(
2,0
)
-
q
k
(
1
+
exp
⁡
⁡
(
-
q
k
)
)
2
(
1
-
exp
⁡
⁡
(
-
q
k
)
)
u
(
0,2
)
+
p
[
h
2
3
!
u
(
3,0
)
+
h
4
5
!
u
(
5,0
)
+
⋯
]
-
p
h
(
1
+
exp
⁡
⁡
(
-
p
h
)
)
(
1
-
exp
⁡
⁡
(
-
p
h
)
)
[
h
2
4
!
u
(
4,0
)
+
h
4
6
!
u
(
6,0
)
+
⋯
]
+
q
[
k
2
3
!
u
(
0,3
)
+
k
4
5
!
u
(
0,5
)
+
⋯
]
-
q
k
(
1
+
exp
⁡
⁡
(
-
q
k
)
)
(
1
-
exp
⁡
⁡
(
-
q
k
)
)
[
k
2
4
!
u
(
0,4
)
+
k
4
6
!
u
(
0,6
)
+
⋯
]
.
Using Taylor expansion again, we arrive at
(39)
p
h
(
1
+
exp
⁡
⁡
(
-
p
h
)
)
2
(
1
-
exp
⁡
⁡
(
-
p
h
)
)
=
1
+
p
2
h
2
12
+
O
(
h
4
)
,
q
k
(
1
+
exp
⁡
⁡
(
-
q
k
)
)
(
1
-
exp
⁡
⁡
(
-
q
k
)
)
=
1
+
q
2
k
2
12
+
O
(
k
4
)
.
Substituting the two formulae above into (38), we obtain
(40)
R
=
(
-
u
(
2,0
)
-
u
(
0,2
)
+
p
u
(
1,0
)
+
q
u
(
0,1
)
-
f
)
+
1
12
(
-
p
2
u
(
2,0
)
+
2
p
u
(
3,0
)
-
u
(
4,0
)
)
h
2
+
1
12
(
-
q
2
u
(
0,2
)
+
2
q
u
(
0,3
)
-
u
(
0,4
)
)
k
2
+
O
(
h
4
+
k
4
)
.
Assuming that the solution to differential equations (1) and (2) satisfies u
∈
C
4
(
Ω
-
), that is to say, u
(
i
,
j
)
(
x
,
y
)
,
(
i
,
j
=
0,1
,
2,3
,
4
,
i
+
j
≤
4
) are bounded by some constant, we use the norm
(41)
∥
u
∥
C
4
(
Ω
-
)
=
max
⁡
0
≤
i
,
j
≤
4
i
+
j
≤
4
sup
⁡
(
x
,
y
)
∈
Ω
|
u
(
i
,
j
)
(
x
,
y
)
|
,
and assuming p
(
x
,
y
)
≤
α
- and q
(
x
,
y
)
≤
β
-, then we can easily see from (40) that the truncation error of the difference approximation satisfies
(42)
R
≤
C
1
(
h
2
+
k
2
)
∥
u
∥
C
4
(
Ω
-
)
.
Here the constant C
1
=
max
⁡
⁡
{
(
α
-
2
+
2
α
-
+
1
)
/
12
,
(
β
-
2
+
2
β
-
+
1
)
/
12
}, which depends on the coefficients of differential operator L defined in (1) but not on mesh steps h
,
k, or U. We summarize the conclusion as below.Lemma 1.
Assuming that the solutionu
(
x
,
y
) of differential equations (1) and (2) satisfies u
∈
C
4
(
Ω
-
), and both p
(
x
,
y
)
≤
α
- and q
(
x
,
y
)
≤
β
- are satisfied in Ω, the truncation error of EDS0 is of second-order which is shown in (42).Remark 2.
In all cases, (42) will hold for EDS0. For instance, if p
(
x
,
y
)
≡
0 in Ω, from (40) we can see that (42) still holds with C
1
=
max
⁡
⁡
{
1
/
12
,
(
β
-
2
+
2
β
-
+
1
)
/
12
}. In particular, if p
(
x
,
y
)
≡
q
(
x
,
y
)
≡
0 in Ω, then (1) becomes Poisson’s equation, and (42) holds with C
1
=
max
⁡
⁡
{
1
/
12,1
/
12
}as well known to us and so do other cases.
Now we complete the local truncation error estimation for EDS0. As for EDS (32), if f
(
x
,
y
)
∈
C
1
(
Ω
-
), we can see that EDS is just a slight modification of EDS0 without changing the latter’s accuracy order since γ
01
=
-
(
q
k
2
/
12
)
+
(
q
3
k
4
/
720
)
-
(
q
5
k
6
/
30240
)
+
⋯
,
γ
10
=
-
(
p
h
2
/
12
)
+
(
p
3
h
4
/
720
)
-
(
p
5
h
6
/
30240
)
+
⋯.


### 3.2. Convergence Analysis

For continuing our analysis we show a discrete maximum principle next.Lemma 3 (maximum principle).
IfU is a mesh function such that L
h
U
j
≤
0
(
L
h
U
j
≥
0
) for any x
j
∈
Ω
h
,
k, then U attains its maximum (minimum) value for some x
j
∈
Γ
h
,
k.Proof.
From discussions in Section2, the difference operator L
h defined in (32) is of positive type, so it satisfies the discrete maximum principle. Detailed proof can be seen in [29].
The maximum principle leads to a stability estimate in the discrete maximum norm, as we will now demonstrate it. In this paper, we write discrete maximum norm for mesh functions:(43)
|
U
|
S
=
max
⁡
x
j
∈
S
|
U
j
|
.Lemma 4.
WithL
h defined in (32), if p
(
x
i
,
y
i
)
≥
α
_
>
0 (or q
(
x
i
,
y
i
)
≥
β
_
>
0) in Ω
h
,
k, we have, for any mesh function U,
(44)
|
U
|
Ω
-
h
,
k
≤
|
U
|
Γ
h
,
k
+
1
α
_
|
L
h
U
|
Ω
h
,
k
(
or
|
U
|
Ω
-
h
,
k
≤
|
U
|
Γ
h
,
k
+
1
β
_
|
L
h
U
|
Ω
h
,
k
)
.Proof.
We only consider the casep
(
x
,
y
)
≥
α
_
>
0. Set W
(
x
,
y
)
=
x and define the mesh function W
i
=
W
(
x
i
,
y
i
); then W
i
≥
0 and L
h
W
i
=
p
i
≥
α
_ in Ω
h
,
k. Furthermore, setting V
i
±
=
±
U
i
-
(
1
/
α
_
)
|
L
h
U
|
Ω
h
,
k
W
i, we get
(45)
L
h
V
i
±
=
±
L
h
U
i
-
1
α
_
|
L
h
U
|
Ω
h
,
k
L
h
W
i
=
±
L
h
U
i
-
p
i
α
_
|
L
h
U
|
Ω
h
,
k
≤
0
,
and since W
i
≥
0 for (
x
i
,
y
i
)
∈
Γ
h
,
k, it follows from Lemma 3 that V
i
±
≤
|
U
|
Γ
h
,
k. Since W
i
≤
1 in Ω
h
,
k, so
(46)
|
U
|
Ω
h
,
k
≤
|
U
|
Γ
h
,
k
+
1
α
_
|
L
h
U
|
Ω
h
,
k
|
W
|
Ω
h
,
k
≤
|
U
|
Γ
h
,
k
+
1
α
_
|
L
h
U
|
Ω
h
,
k
.
This proves Lemma 4.Lemma4 immediately shows the existence and uniqueness of the solution of the EDS0. We are now ready for an error estimate.Theorem 5.
LetU and u be the solutions of difference scheme (32)-(33) and differential equation (1)-(2), respectively, and assuming u
∈
C
4
(
Ω
-
) and conditions for Lemmas 1 and 4 are satisfied, then
(47)
|
U
-
u
|
Ω
-
h
,
k
≤
C
(
h
2
+
k
2
)
∥
u
∥
C
4
(
Ω
-
)
,
where C
=
C
1
/
α
_.Proof.
From (3) and (32), the error e
0
=
U
0
-
u
0 satisfies, at the interior reference mesh points,
(48)
L
h
e
0
=
-
R
0
.
Since R
0
≤
C
1
(
h
2
+
k
2
)
∥
u
∥
C
4
(
Ω
-
), from (42), the result therefore follows by application of Lemma 4 to e
j, since e
j
=
0 for (
x
j
,
y
j
)
∈
Γ
h
,
k.
Iff
(
x
,
y
)
∈
C
1
(
Ω
-
), that is, f
(
i
,
j
)
(
x
,
y
)
,
(
i
,
j
=
0,1
,
i
+
j
≤
1) are bounded by some constant, say by C
3. Since γ
01
,
γ
10 can be considered as just relatively minor modifications of EDS0 (equal to second-order), so the truncation error of EDS is still of second-order accuracy. The maximum principle and convergent theorem holds too for EDS.Remark 6.
So far all the derivation process and analysis of EDS are based on the fact thatp
(
x
,
y
)
≥
0 and q
(
x
,
y
)
≥
0, but, in fact, in case of p
(
x
,
y
)
≤
0
,
q
(
x
,
y
)
≤
0 (see that it will not change the sign of ω
i
(
i
=
0,1
,
2,3
,
4
) from (25)–(29)), the method is still effective and similar theoretical analysis as we have shown above could be done.


## 3.1.  Local Truncation Error

We begin with analyzing the truncation error of the scheme about EDS0. For brevity we only consider the case thatp
(
x
,
y
)
>
0 and q
(
x
,
y
)
>
0; similar analysis can be done for other cases.Estimating the truncation error of EDS0 from (3) yields
(37)
R
=
ω
0
u
0
-
(
ω
1
u
1
+
ω
2
u
2
+
ω
3
u
3
+
ω
4
u
4
)
-
f
.
Substituting Taylor expansions (5)–(7) together with values ω
i
,
(
i
=
0,1
,
2,3
,
4
) given in (25)–(29) into (37), and making some arrangement, we have
(38)
R
=
-
f
+
p
u
(
1,0
)
+
q
u
(
0,1
)
-
p
h
(
1
+
exp
⁡
⁡
(
-
p
h
)
)
2
(
1
-
exp
⁡
⁡
(
-
p
h
)
)
u
(
2,0
)
-
q
k
(
1
+
exp
⁡
⁡
(
-
q
k
)
)
2
(
1
-
exp
⁡
⁡
(
-
q
k
)
)
u
(
0,2
)
+
p
[
h
2
3
!
u
(
3,0
)
+
h
4
5
!
u
(
5,0
)
+
⋯
]
-
p
h
(
1
+
exp
⁡
⁡
(
-
p
h
)
)
(
1
-
exp
⁡
⁡
(
-
p
h
)
)
[
h
2
4
!
u
(
4,0
)
+
h
4
6
!
u
(
6,0
)
+
⋯
]
+
q
[
k
2
3
!
u
(
0,3
)
+
k
4
5
!
u
(
0,5
)
+
⋯
]
-
q
k
(
1
+
exp
⁡
⁡
(
-
q
k
)
)
(
1
-
exp
⁡
⁡
(
-
q
k
)
)
[
k
2
4
!
u
(
0,4
)
+
k
4
6
!
u
(
0,6
)
+
⋯
]
.
Using Taylor expansion again, we arrive at
(39)
p
h
(
1
+
exp
⁡
⁡
(
-
p
h
)
)
2
(
1
-
exp
⁡
⁡
(
-
p
h
)
)
=
1
+
p
2
h
2
12
+
O
(
h
4
)
,
q
k
(
1
+
exp
⁡
⁡
(
-
q
k
)
)
(
1
-
exp
⁡
⁡
(
-
q
k
)
)
=
1
+
q
2
k
2
12
+
O
(
k
4
)
.
Substituting the two formulae above into (38), we obtain
(40)
R
=
(
-
u
(
2,0
)
-
u
(
0,2
)
+
p
u
(
1,0
)
+
q
u
(
0,1
)
-
f
)
+
1
12
(
-
p
2
u
(
2,0
)
+
2
p
u
(
3,0
)
-
u
(
4,0
)
)
h
2
+
1
12
(
-
q
2
u
(
0,2
)
+
2
q
u
(
0,3
)
-
u
(
0,4
)
)
k
2
+
O
(
h
4
+
k
4
)
.
Assuming that the solution to differential equations (1) and (2) satisfies u
∈
C
4
(
Ω
-
), that is to say, u
(
i
,
j
)
(
x
,
y
)
,
(
i
,
j
=
0,1
,
2,3
,
4
,
i
+
j
≤
4
) are bounded by some constant, we use the norm
(41)
∥
u
∥
C
4
(
Ω
-
)
=
max
⁡
0
≤
i
,
j
≤
4
i
+
j
≤
4
sup
⁡
(
x
,
y
)
∈
Ω
|
u
(
i
,
j
)
(
x
,
y
)
|
,
and assuming p
(
x
,
y
)
≤
α
- and q
(
x
,
y
)
≤
β
-, then we can easily see from (40) that the truncation error of the difference approximation satisfies
(42)
R
≤
C
1
(
h
2
+
k
2
)
∥
u
∥
C
4
(
Ω
-
)
.
Here the constant C
1
=
max
⁡
⁡
{
(
α
-
2
+
2
α
-
+
1
)
/
12
,
(
β
-
2
+
2
β
-
+
1
)
/
12
}, which depends on the coefficients of differential operator L defined in (1) but not on mesh steps h
,
k, or U. We summarize the conclusion as below.Lemma 1.
Assuming that the solutionu
(
x
,
y
) of differential equations (1) and (2) satisfies u
∈
C
4
(
Ω
-
), and both p
(
x
,
y
)
≤
α
- and q
(
x
,
y
)
≤
β
- are satisfied in Ω, the truncation error of EDS0 is of second-order which is shown in (42).Remark 2.
In all cases, (42) will hold for EDS0. For instance, if p
(
x
,
y
)
≡
0 in Ω, from (40) we can see that (42) still holds with C
1
=
max
⁡
⁡
{
1
/
12
,
(
β
-
2
+
2
β
-
+
1
)
/
12
}. In particular, if p
(
x
,
y
)
≡
q
(
x
,
y
)
≡
0 in Ω, then (1) becomes Poisson’s equation, and (42) holds with C
1
=
max
⁡
⁡
{
1
/
12,1
/
12
}as well known to us and so do other cases.
Now we complete the local truncation error estimation for EDS0. As for EDS (32), if f
(
x
,
y
)
∈
C
1
(
Ω
-
), we can see that EDS is just a slight modification of EDS0 without changing the latter’s accuracy order since γ
01
=
-
(
q
k
2
/
12
)
+
(
q
3
k
4
/
720
)
-
(
q
5
k
6
/
30240
)
+
⋯
,
γ
10
=
-
(
p
h
2
/
12
)
+
(
p
3
h
4
/
720
)
-
(
p
5
h
6
/
30240
)
+
⋯.


## 3.2. Convergence Analysis

For continuing our analysis we show a discrete maximum principle next.Lemma 3 (maximum principle).
IfU is a mesh function such that L
h
U
j
≤
0
(
L
h
U
j
≥
0
) for any x
j
∈
Ω
h
,
k, then U attains its maximum (minimum) value for some x
j
∈
Γ
h
,
k.Proof.
From discussions in Section2, the difference operator L
h defined in (32) is of positive type, so it satisfies the discrete maximum principle. Detailed proof can be seen in [29].
The maximum principle leads to a stability estimate in the discrete maximum norm, as we will now demonstrate it. In this paper, we write discrete maximum norm for mesh functions:(43)
|
U
|
S
=
max
⁡
x
j
∈
S
|
U
j
|
.Lemma 4.
WithL
h defined in (32), if p
(
x
i
,
y
i
)
≥
α
_
>
0 (or q
(
x
i
,
y
i
)
≥
β
_
>
0) in Ω
h
,
k, we have, for any mesh function U,
(44)
|
U
|
Ω
-
h
,
k
≤
|
U
|
Γ
h
,
k
+
1
α
_
|
L
h
U
|
Ω
h
,
k
(
or
|
U
|
Ω
-
h
,
k
≤
|
U
|
Γ
h
,
k
+
1
β
_
|
L
h
U
|
Ω
h
,
k
)
.Proof.
We only consider the casep
(
x
,
y
)
≥
α
_
>
0. Set W
(
x
,
y
)
=
x and define the mesh function W
i
=
W
(
x
i
,
y
i
); then W
i
≥
0 and L
h
W
i
=
p
i
≥
α
_ in Ω
h
,
k. Furthermore, setting V
i
±
=
±
U
i
-
(
1
/
α
_
)
|
L
h
U
|
Ω
h
,
k
W
i, we get
(45)
L
h
V
i
±
=
±
L
h
U
i
-
1
α
_
|
L
h
U
|
Ω
h
,
k
L
h
W
i
=
±
L
h
U
i
-
p
i
α
_
|
L
h
U
|
Ω
h
,
k
≤
0
,
and since W
i
≥
0 for (
x
i
,
y
i
)
∈
Γ
h
,
k, it follows from Lemma 3 that V
i
±
≤
|
U
|
Γ
h
,
k. Since W
i
≤
1 in Ω
h
,
k, so
(46)
|
U
|
Ω
h
,
k
≤
|
U
|
Γ
h
,
k
+
1
α
_
|
L
h
U
|
Ω
h
,
k
|
W
|
Ω
h
,
k
≤
|
U
|
Γ
h
,
k
+
1
α
_
|
L
h
U
|
Ω
h
,
k
.
This proves Lemma 4.Lemma4 immediately shows the existence and uniqueness of the solution of the EDS0. We are now ready for an error estimate.Theorem 5.
LetU and u be the solutions of difference scheme (32)-(33) and differential equation (1)-(2), respectively, and assuming u
∈
C
4
(
Ω
-
) and conditions for Lemmas 1 and 4 are satisfied, then
(47)
|
U
-
u
|
Ω
-
h
,
k
≤
C
(
h
2
+
k
2
)
∥
u
∥
C
4
(
Ω
-
)
,
where C
=
C
1
/
α
_.Proof.
From (3) and (32), the error e
0
=
U
0
-
u
0 satisfies, at the interior reference mesh points,
(48)
L
h
e
0
=
-
R
0
.
Since R
0
≤
C
1
(
h
2
+
k
2
)
∥
u
∥
C
4
(
Ω
-
), from (42), the result therefore follows by application of Lemma 4 to e
j, since e
j
=
0 for (
x
j
,
y
j
)
∈
Γ
h
,
k.
Iff
(
x
,
y
)
∈
C
1
(
Ω
-
), that is, f
(
i
,
j
)
(
x
,
y
)
,
(
i
,
j
=
0,1
,
i
+
j
≤
1) are bounded by some constant, say by C
3. Since γ
01
,
γ
10 can be considered as just relatively minor modifications of EDS0 (equal to second-order), so the truncation error of EDS is still of second-order accuracy. The maximum principle and convergent theorem holds too for EDS.Remark 6.
So far all the derivation process and analysis of EDS are based on the fact thatp
(
x
,
y
)
≥
0 and q
(
x
,
y
)
≥
0, but, in fact, in case of p
(
x
,
y
)
≤
0
,
q
(
x
,
y
)
≤
0 (see that it will not change the sign of ω
i
(
i
=
0,1
,
2,3
,
4
) from (25)–(29)), the method is still effective and similar theoretical analysis as we have shown above could be done.


## 4. Numerical Experiments

In this section, we test our methods with four examples. Example7 is a convection-diffusion equation with constant coefficients and has exponential boundary layers. Example 8 is the one with constant coefficients and has no boundary layers, but the convection terms have singular lines on the boundary; that is, p
(
x
,
y
) and q
(
x
,
y
) vanish on some part of the boundary. Example 9 is convection-dominated equation with variable coefficients and has corner exponential boundary layer (Figure 5). At last, we provide Example 10 with both exponential and parabolic boundary layers.We discuss the accuracy and stability of the EDS in comparison with CDS and UDS methods and some results given by other authors.Still the computational domain for all the following problems is a unit squareΩ
=
(
0,1
)
×
(
0,1
). For the sake of simplicity, we assume that the step sizes h and k are equal. However, this is not a restriction at all for the methods proposed in this paper. The mesh points Ω
h
,
k and the boundary points Γ
h
,
k can be simplified to be denoted by Ω
h and Γ
h.The exact solutions for these problems are known, so the maximum pointwise error at all mesh points can be calculated using the formula(49)
E
h
:
=
max
⁡
i
∈
Ω
h
|
U
i
-
u
i
|
,
where u is the exact solution and U is the numerical solution. The rate of convergence is defined by
(50)
rate
=
log
⁡
⁡
(
E
2
h
/
E
h
)
log
⁡
⁡
2
.Example 7.
Consider convection-diffusion equation with constant coefficients:(51)
-
u
x
x
-
u
y
y
+
Re
(
u
x
+
u
y
)
=
0
,
(
x
,
y
)
∈
Ω
,
and with boundary condition:
(52)
u
(
x
,
0
)
=
0
,
u
(
x
,
1
)
=
exp
⁡
⁡
(
Re
x
)
-
1
exp
⁡
⁡
(
Re
)
-
1
,
u
(
0
,
y
)
=
0
,
u
(
1
,
y
)
=
exp
⁡
⁡
(
Re
y
)
-
1
exp
⁡
⁡
(
Re
)
-
1
.
In order to prevent overflow of the computer, we rewrite (52) as
(53)
u
(
x
,
0
)
=
0
,
u
(
x
,
1
)
=
exp
⁡
⁡
(
Re
(
x
-
1
)
)
-
exp
⁡
⁡
(
-
Re
)
1
-
exp
⁡
⁡
(
-
Re
)
u
(
0
,
y
)
=
0
,
u
(
1
,
y
)
=
exp
⁡
⁡
(
Re
(
y
-
1
)
)
-
exp
⁡
⁡
(
-
Re
)
1
-
exp
⁡
⁡
(
-
Re
)
.
Comparison of (51) with (1) shows that
(54)
p
(
x
,
y
)
=
Re
,
q
(
x
,
y
)
=
Re
,
f
(
x
,
y
)
=
0
.
The exact solution of this problem is
(55)
u
(
x
,
y
)
=
exp
⁡
⁡
(
-
Re
(
1
-
x
)
)
-
exp
⁡
⁡
(
-
Re
)
1
-
exp
⁡
⁡
(
-
Re
)
·
exp
⁡
⁡
(
-
Re
(
1
-
y
)
)
-
exp
⁡
⁡
(
-
Re
)
1
-
exp
⁡
⁡
(
-
Re
)
.
Figure2 shows that there is no boundary layer when Re
=
1, but when Re is large the boundary layers will emerge near x
=
1 and y
=
1. Notice that since the source term f
(
x
,
y
)
=
0, so EDS and EDS0 are virtually identical.
The Gauss-Seidel iteration is used to solve the resulting systems of equations. The convergence criterion for the iteration is chosen to be3
e
-
16, which is almost the minimum positive machine-precision number (the minimum positive machine-precision number is 2.22045
*
1
0
-
16 in double-precision computer system). The maximum absolute errors for different values of h and k are given in Table 1.
From Table1 we can see that EDS gives more accurate results compared with CDS and UDS. In fact, the result of the EDS is almost as accurate as the machine-precision just because EDS is the exact scheme for Example 7. An exact difference scheme [30] is one for which the computed solution agrees exactly with the true solution at the mesh points. It is easy to verify this fact because if we substitute the exact solution (55) and the values of the coefficients (25)–(29) into the truncation error (37) we will get
(56)
R
0
=
ω
0
u
0
-
(
ω
1
u
1
+
ω
2
u
2
+
ω
3
u
3
+
ω
4
u
4
)
=
ω
0
u
(
x
0
,
y
0
)
-
(
ω
1
u
(
x
0
+
h
,
y
0
)
+
ω
2
u
(
x
0
,
y
0
+
k
)
AAAAAAAAA
+
ω
3
u
(
x
0
-
h
,
y
0
)
+
ω
4
u
(
x
0
,
y
0
-
k
)
)
=
0
at any reference point.
Table2 shows clearly that with Re increases, the results become more and more accurate which validates the unconditional numerical stability in solving these linear system of equations of the EDS scheme. Figure 3 shows the absolute errors of EDS for Example 7, with mesh steps h
=
k
=
1
/
32.
We continue to implement our experiments; the convergence criterion for the Gauss Seidel iteration is chosen to be1
0
-
10 for the following problems.
We prefer to apply the Richardson extrapolation technique here to compute a fourth-order accurate solution since the technique will not destroy the unconditional stability of EDS0 and EDS. Now we denote byU
i
,
j
h
,
k the difference solution at any mesh point (
x
i
,
y
j
)
∈
Ω
h
,
k; the Richardson extrapolation can be written as
(57)
U
i
,
j
h
,
k
=
2
r
U
2
i
,
2
j
h
/
2
,
k
/
2
-
U
i
,
j
h
,
k
2
r
-
1
,
where r
=
2 is the order of accuracy before the extrapolation and the order of accuracy will be increased to r
+
2 after the extrapolation.Table 1
Maximum absolute errors with Re = 1 (left) and Re = 100 (right), Example7.
h
=
k
Re = 1
Re = 100
UDS
CDS
EDS
UDS
CDS
EDS
1/8
3.57
e
-
3
8.05
e
-
5
1.55
e
-
15
2.66
e
-
7
—
1.21
e
-
22
1/16
1.87
e
-
3
2.03
e
-
5
6.83
e
-
15
2.45
e
-
4
—
2.96
e
-
21
1/32
9.63
e
-
4
5.12
e
-
6
2.70
e
-
14
7.70
e
-
3
—
1.03
e
-
18
1/64
4.87
e
-
4
1.28
e
-
6
1.03
e
-
13
3.03
e
-
2
1.70
e
-
2
4.19
e
-
17
1/128
2.44
e
-
4
3.20
e
-
7
4.06
e
-
13
3.22
e
-
2
6.68
e
-
3
2.10
e
-
16
1/256
1.23
e
-
4
8.00
e
-
8
1.58
e
-
12
2.01
e
-
2
1.64
e
-
3
9.42
e
-
16
— means CDS diverges for Re = 100 when the mesh size is big.Table 2
Absolute errors for EDS solution of Example7 with different values of Re.
h
=
k
Re = 5
Re = 50
Re = 500
Re = 5000
1/8
8.60
e
-
16
5.36
e
-
19
1.86
e
-
82
0.00
e
+
0
1/16
4.30
e
-
15
7.74
e
-
18
9.78
e
-
42
0.00
e
+
0
1/32
1.82
e
-
14
2.71
e
-
16
4.49
e
-
28
1.33
e
-
204
1/64
7.43
e
-
14
1.57
e
-
15
1.01
e
-
20
2.06
e
-
84
1/128
2.87
e
-
13
6.73
e
-
15
3.20
e
-
18
6.39
e
-
52
1/256
2.72
e
-
14
2.72
e
-
14
9.58
e
-
17
1.79
e
-
26
1/512
1.08
e
-
13
1.08
e
-
13
8.79
e
-
17
6.81
e
-
22Exact solution of Example7 for different values of Re.
(a)
Re
=
1
(b)
Re
=
10
(c)
Re
=
100
(d)
Re
=
1000Max absolute errors of EDS, Example7.
(a)
Re
=
1
,
h
=
k
=
1
/
32
(b)
Re
=
100
,
h
=
k
=
1
/
32Example 8.
Consider(58)
-
u
x
x
-
u
y
y
+
P
x
u
x
+
P
y
u
y
=
f
(
x
,
y
)
,
(
x
,
y
)
∈
Ω
,
where boundary condition:
(59)
u
(
x
,
0
)
=
0
,
u
(
x
,
1
)
=
0
,
u
(
0
,
y
)
=
0
,
u
(
1
,
y
)
=
0
,
holds in Ω. The source term f
(
x
,
y
) is determined so as u
(
x
,
y
)
=
x
y
(
1
-
x
)
(
1
-
y
) being a solution of (58).
The figure of the exact solution is presented on the left in Figure4. The solution has no relation to P and has no boundary layer.
From Table3 we can see that EDS0 and EDS get better results than UDS for Example 8, and as step h decreases, the rate of EDS0 and EDS get closer and closer to second-order, whereas the rate of UDS gets closer to 1.
In order to get much more accurate results, we carry out Richardson’s extrapolation. In this paper, we denote the method of the EDS0 (EDS) with Richardson’s extrapolation by EDS0-RE (EDS-RE). The results are given in Table4. We can see that Richardson’s extrapolation really increases the accuracy a lot.
Table5 shows that EDS is a little bit more accurate than EDS0, but EDS-RE is more accurate than EDS0-RE.Table 3
Maximum absolute errors and rate of convergence, Example8 with P
=
100.
h
=
k
UDS
EDS0
EDS
Errors
Rate
Errors
Rate
Errors
Rate
1/8
3.07
e
-
2
2.37
e
-
2
1.00
e
-
2
1/16
1.75
e
-
2
0.81
8.64
e
-
3
1.46
5.07
e
-
3
0.98
1/32
9.31
e
-
3
0.91
2.51
e
-
3
1.78
1.68
e
-
3
1.60
1/64
4.80
e
-
3
0.96
6.59
e
-
4
1.93
4.48
e
-
4
1.90
1/128
2.44
e
-
3
0.98
1.67
e
-
4
1.98
1.14
e
-
4
1.97
1/256
1.23
e
-
3
0.99
4.19
e
-
5
1.99
2.88
e
-
5
1.98Table 4
Maximum absolute errors and rate of EDS0-RE and EDS-RE, Example8 with P
=
100.
EDS0-RE
h
=
k
1/8
1/16
1/32
1/64
1/128
error
3.20
e
-
3
4.71
e
-
4
4.03
e
-
5
2.83
e
-
6
2.62
e
-
7
Rate
2.76
3.55
3.83
3.43
EDS-RE
error
1.46
e
-
3
3.26
e
-
4
3.80
e
-
5
2.75
e
-
6
2.55
e
-
7
Rate
2.16
3.10
3.79
3.43Table 5
Maximum absolute errors in the computed solution, Example8 with P
=
1000.
h
=
k
UDS
EDS0
EDS
EDS0-RE
EDS-RE
1/8
3.51
e
-
2
3.45
e
-
2
2.20
e
-
2
1.16
e
-
2
1.46
e
-
3
1/16
1.89
e
-
2
1.78
e
-
2
1.56
e
-
2
5.17
e
-
3
3.26
e
-
4
1/32
9.96
e
-
3
8.39
e
-
3
7.64
e
-
3
1.93
e
-
3
3.80
e
-
5
1/64
5.08
e
-
3
3.56
e
-
3
3.14
e
-
3
5.29
e
-
4
2.75
e
-
6
1/128
2.57
e
-
3
1.29
e
-
3
1.14
e
-
3
8.82
e
-
5
2.55
e
-
7
1/256
1.29
e
-
3
3.88
e
-
4
3.54
e
-
4Exact and EDS numerical solution for Example8.
(a)
Exact solution
(b)
EDS numerical solution whenP
=
100
,
h
=
k
=
1
/
16Exact solution of Example9 for different values of P.
(a)
P
=
5
(b)
P
=
50Example 9.
Consider convention-dominated problem with variable convection coefficient:(60)
-
u
x
x
-
u
y
y
+
P
(
1
+
x
)
u
x
+
P
(
1
+
y
)
u
y
=
f
(
x
,
y
)
,
(
x
,
y
)
∈
Ω
,
with boundary condition
(61)
u
(
x
,
0
)
=
1
,
u
(
x
,
1
)
=
1
-
exp
⁡
⁡
(
-
P
(
1
-
x
)
)
1
-
exp
⁡
⁡
(
-
P
)
,
u
(
0
,
y
)
=
1
,
u
(
1
,
y
)
=
1
-
exp
⁡
⁡
(
-
P
(
1
-
y
)
)
1
-
exp
⁡
⁡
(
-
P
)
.
The source term f
(
x
,
y
) is determined so as u
(
x
,
y
)
=
(
1
-
exp
⁡
⁡
(
-
P
(
1
-
x
y
)
)
)
/
(
1
-
exp
⁡
⁡
(
-
P
)
) being a solution of (60).
Corner boundary layer occurs atx
=
1 and y
=
1. Table 6 lists maximum absolute errors in the computed solution of the problem. Table 6 again shows that EDS method gets much more accurate results for convection-dominated problems.Table 6
Maximum absolute errors in the computed solution of Example9.
h
=
k
P
=
10
P
=
100
P
=
1000
UDS
EDS
UDS
EDS
UDS
EDS
1/8
4.60
e
-
1
1.32
e
-
1
8.00
e
-
2
3.44
e
-
9
8.45
e
-
3
0.00
e
+
0
1/16
4.88
e
-
1
3.19
e
-
1
1.47
e
-
1
6.00
e
-
5
1.63
e
-
2
0.00
e
+
0
1/32
6.19
e
-
1
5.43
e
-
1
2.61
e
-
1
2.13
e
-
3
3.17
e
-
2
0.00
e
+
0
1/64
7.64
e
-
1
7.35
e
-
1
4.13
e
-
1
5.58
e
-
2
6.15
e
-
2
3.52
e
-
12
1/128
8.65
e
-
1
8.56
e
-
1
4.90
e
-
1
2.51
e
-
1
1.17
e
-
1
3.73
e
-
6Example 10.
Consider a constant coefficient convection-diffusion equation (1) with p
(
x
,
y
)
=
1
/
ϵ
,
q
(
x
,
y
)
=
0, and f
(
x
,
y
)
=
0; that is,
(62)
-
ϵ
(
u
x
x
+
u
y
y
)
+
u
x
=
0
,
(
x
,
y
)
∈
Ω
.
The boundary condition is prescribed as
(63)
u
(
x
,
0
)
=
0
,
u
(
x
,
1
)
=
0
,
u
(
0
,
y
)
=
sin
⁡
(
π
y
)
,
u
(
1
,
y
)
=
2
sin
⁡
(
π
y
)
.
The exact solution [31] is
(64)
u
(
x
,
y
)
=
(
exp
⁡
⁡
(
x
2
ϵ
)
sin
⁡
(
π
y
)
×
[
2
exp
⁡
⁡
(
-
1
2
ϵ
)
sinh
⁡
⁡
(
θ
x
)
+
sinh
⁡
⁡
(
θ
(
1
-
x
)
)
]
)
×
(
sinh
⁡
⁡
(
θ
)
)
-
1
,
where θ
2
=
π
2
+
1
/
(
4
ϵ
2
).
This problem represents a convection dominated flow and was used as one of the test problems by Gupta et al. [18] and many other authors [31, 32]. Figure 6 shows the exact solution with ϵ
=
0.1 and ϵ
=
0.001. For most part of the domain, the exact solution is smooth, but it has a steep boundary layer (exponential layer) along the downstream edge at x
=
1 and has parabolic layers along the edges at y
=
0 and y
=
1.
Results of Example10 from different methods are shown in Table 7. We see that EDS and Richardson’s extrapolation for EDS get more accurate results.
Figure7 draws errors of EDS for Example 10 with ϵ
=
0.01
,
h
=
k
=
1
/
32. The errors of different schemes along x
=
0.5 and y
=
0.5 are plotted in Figure 8. We can see that UDS cannot approximate the exact solution quite well near the exponential layer. With fixed coordinate y, the errors of all the methods increase with x which increases from x
=
0 to x
=
1 along the horizontal x axis, but EDS and EDS-RE get much better results (Table 8).Table 7
Maximum absolute errors and rate in the computed solution, Example10 with ϵ
=
0.1.
h
=
k
UDS
CDS
EDS
EDS-RE
Errors
Rate
Errors
Rate
Errors
Rate
Errors
Rate
1/8
2.41
e
-
1
7.50
e
-
2
7.93
e
-
3
3.51
e
-
5
1/16
1.39
e
-
1
0.79
1.65
e
-
2
2.19
2.02
e
-
3
1.96
2.16
e
-
6
3.93
1/32
7.76
e
-
2
0.84
4.19
e
-
3
1.98
5.09
e
-
4
1.99
1.32
e
-
7
3.98
1/64
4.11
e
-
2
0.92
1.04
e
-
3
2.01
1.28
e
-
4
2.00
8.16
e
-
9
4.00
1/128
2.12
e
-
2
0.96
2.59
e
-
4
2.00
3.19
e
-
5
2.00
Note. The convergence criterion of the Gauss Seidel iteration for EDS-RE scheme is chosen to be 3
×
1
0
-
16 to obtain highly accurate results.Table 8
Max absolute error in the computed solution, Example10 with ϵ
=
0.001.
h
=
k
UDS
EDS
EDS-RE
MDC
1/8
8.12
e
-
3
1.15
e
-
4
3.51
e
-
5
1/16
1.60
e
-
2
3.21
e
-
5
2.16
e
-
6
1.31
e
-
4
1/32
3.13
e
-
2
8.98
e
-
6
1.32
e
-
7
5.82
e
-
5
1/64
6.07
e
-
2
2.58
e
-
6
8.15
e
-
9
2.72
e
-
5
1/128
1.15
e
-
1
2.58
e
-
6
8.15
e
-
9
2.72
e
-
5
Note. The last column of data comes from [31] of the MDC finite-element solution.Exact solution of Example10.
(a)
ϵ
=
0.1
(b)
ϵ
=
0.001Figure 7
Errors of EDS for Example10 with ϵ
=
0.01
,
h
=
k
=
1
/
32.Errors of numerical results, Example10 with ϵ
=
0.01
,
h
=
k
=
1
/
32.
(a)
Errors of three methods forx
=
0.5
(b)
Errors of three methods fory
=
0.5
(c)
Errors of EDS forx
=
0.5
(d)
Errors of EDS fory
=
0.5


## 5. Conclusions

Idea of direction changing and order reducing has been used to construct a compact exponential FD scheme over a five-point stencil for 2D convection-diffusion problem in the paper. The scheme is of concise exponential type and unconditionally stable; the convergent accurate order is second. Four numerical examples show that EDS0 and EDS can get much highly accurate results for different kinds of problems, including those with boundary layers and without layers. Richardson’s extrapolation can improve the accurate order to four.We should mention that if one takes a scheme for its 1D analogue to form the “tensor product” of the scheme in 2D case, one can directly arrive at a scheme on five-point stencil but not on a nine-point stencil. But our method can be applied to construct much more accurate scheme over a nine-point stencil. Also, the idea can be used to form schemes on nonuniform grids. Theoretically, idea of direction changing and order reducing can also be used to solve varieties of problems, for example, parabolic equationu
t
+
b
(
x
)
u
(
x
)
=
ν
u
x
x, and so forth. Those will be investigated in future studies.



---

*Source: 457938-2014-06-29.xml*