A Comparison Study of Two High Accuracy Numerical Methods for a Parabolic System in Air Pollution Modelling

0.62 MB·

by  Ivan Dimov

A Comparison Study of Two High Accuracy
Numerical Methods for a Parabolic System in
Air Pollution Modelling
7
1
0 I. Dimova, J. Kandilarovb, V. Todorova, L. Vulkovb
2

n aInstitute of Information and Communication Technologies, BAS, Bulgaria
a
J bCenter of Applied Mathematics and Informatics, University of Ruse, Ruse 7017,

1 Bulgaria
1

]
A
N Abstract
.
h
We present two approaches for enhancing the accuracy of second order finite dif-
t
a ference approximations of two-dimensional semilinear parabolic systems. These are
m
the fourth order compact difference scheme and the fourth order scheme based on
[
Richardson extrapolation. Our interest is concentrated on a system of ten parabolic

1 partial differential equations in air pollution modeling. We analyze numerical ex-
v
periments to compare the two approaches with respect to accuracy, computational
9
4 complexity, non-negativity preserving and etc. Sixth-order approximation based on
0 thefourth-ordercompactdifferenceschemecombinedwithRichardsonextrapolation
3
is also discussed numerically.
0
.
1
0 Key words: air pollution model, semilinear parabolic systems, compact finite
7
difference schemes, Richardson extrapolation.
1
:
v
i
X
r
a
1 Introduction
In many fields of sciences and engineering parabolic equations are always used
to describe many phenomena, so that the finite-difference method that solves
the parabolic equation is always a focus of concern, see e.g. [2,11,14,15]. In
the context of the finite difference discretization, the standard second-order
discretization schemes may need fine griddings to yield approximate solutions
of acceptable accuracy. The resulting large size systems have to be solved,
Email addresses: [email protected], [email protected],
[email protected], [email protected] (L. Vulkov).
Preprint submitted to Elsevier January 12, 2017
which may consume a lot of memory space and CPU cycles even on present
generation supercomputers.
One approach to reduce computational cost in very large-scale modelings and
simulations is to used higher-order discretization methods. Other important
factor affecting the computational efficacy of a discretized method is to solve
the resulting linear and nonlinear systems of algebraic equations. The higher-
order methods usually generate algebraic systems of much smaller size, com-
pared to the lower-order methods.
Because of this and other advantages of high-order methods, there has been
growinginterest ofdeveloping andusing highlyaccuratenumerical schemes for
solving partial differential equations, leading to renewed interest in high-order
compact difference schemes [7,10,13,17,18,19].
Compact schemes, proposed by Kreiss and Oliger [8] use similar stencil, but
requires a scalar tridiagonal or pentadiagonal matrix inversion. In this paper
we use another idea to obtain high-order compact schemes, namely, to operate
on the differential equations as auxiliary relations in order to express hight or-
der derivatives in the truncation error [16,20]. More details and discussions on
construction of compact difference schemes for convection-diffusion problems
can be found in [10,17,20].
In [9] the air pollution problem, which is the base of the mathematical model
of the present paper is stated. A preconditioned iterative solution method for
nonlinear parabolic transport system is done. The ingredients of the method
are implicit Euler discretization in time and FEM discretization in space, then
an outer-inner iteration and preconditioning via an ℓ-tuple of independent
elliptic operators.
Another known approach for increasing the order of accuracy of the finite
difference schemes is the use of Richardson extrapolation [11]. Fourth order
compact difference scheme for a system of two semilinear toy 1D parabolic
equations is derived in [4].
This article is arranged as follows. In Section 2 we present the two dimensional
model problem. In Section 3 the second order central difference scheme (CDS)
is presented and the application of the Richardson extrapolation for higher-
order approximations is analyzed. In Section 4 the fourth-order compact finite
difference schemes (CFDS) for general weakly coupled parabolic system of
two equations is introduce. In Section 5 numerical results and comparisons
are presented and analysed. Concluding remarks are included in Section 6.
2
2 The Two Dimensional Model Problem of Air Pollution
The simulation of various processes in chemistry, physics and engineering uses
models of systems of coupled parabolic problems. In this work we construct
compact high-order finite difference schemes for semilinear parabolic systems
and propose fast algorithms for solution of the nonlinear algebraic equations.
Problems of air pollution transport with coupling in the nonlinear reactions
terms are of our main consideration, namely,
∂u
l −K△u +b ∇u = R (x,y,u ,...,u ), (x,y,t) ∈ Ω×(0,T], (1)
l l l l 1 L
∂t
u = 0, (x,y,t) ∈ ∂Ω×(0,T], (2)
u = u (x,y), (x,y) ∈ Ω, (3)
0
where u = (u ,u ,...,u ), u = u (x,y,t), l = 1,...,L are the concentrations of
1 2 L l l
L chemical species (pollutants) and K > 0 is the diffusion coefficient and Ω ∈
R2 isaboundeddomain.Theassumption regardingconstant K := K = K is
x y
not a restriction for developing our numerical approach. This just corresponds
to the physical model described in [5,6,9].
The main goal of the paper is the application and numerical illustration of
above-mentioned difference approximations to the following real-life parabolic
transport system described in [6]. Following [6,9,21] the advection part in (1)
may be presented in the following form:
∂u ∂u
bl.∇ul = µ(y −yc) l +µ(xc −x) l,
∂x ∂y
where x ∈ (0,X), y ∈ (0,Y), x = X/2, y = Y/2. The nonlinear chemical
c c
part of the model is (see [9]):
R (u ,...,u )=k u −(k u +k u +k u )u ,
1 1 10 5 2 6 5 4 7 3 8 1
R (u ,...,u )=(k u +k u +k u )u −(k +k u )u ,
2 1 10 6 5 4 7 3 8 1 5 9 9 2
R (u ,...,u )=−k u u ,
3 1 10 1 3 9
R (u ,...,u )=2k u u +k u u −k u ,
4 1 10 1 3 9 3 1 8 2 4
R (u ,...,u )=k u (4)
5 1 10 2 5
R (u ,...,u )=k u u ,
6 1 10 9 2 9
R (u ,...,u )=2k u +k u u +k u −k u u ,
7 1 10 2 4 3 1 8 10 9 4 1 7
R (u ,...,u )=4k u u −k u u ,
8 1 10 1 3 9 3 1 8
R (u ,...,u )=k u u +2k u −(k u −k u +k )u ,
9 1 10 4 1 7 8 10 1 3 9 2 10 9
R (u ,...,u )=k u −k u .
10 1 10 7 5 8 10
ThechemicalpartofthemodelisgiveninTable1forthesakeofcompleteness.
The rate coefficients can be found in Table 2. Some of the coefficients belong
3
Table 1
The chemical reactions of the model
1 HC +OH → 4RO +2ALD 6 NO+O → NO +O
2 3 2 2
2 ALD+hν → 2HO +CO 7 O +hν → O +O(1D)
2 3 2
3 RO +NO → NO +ALD+HO 8 O(1D)+H O → 2OH
2 2 2 2
4 NO+HO → NO +OH 9 NO +OH → HNO
2 2 2 3
5 NO +hν → NO+O 3 10 CO+OH → CO +HO
2 3 2 2
Table 2
The coefficients of the chemical reactions
k 6.0e−12 k 1.6e−14
1 6
k 7.8e−05.exp(−0.87/cosθ) k 1.6e−04.exp(−1.9/cosθ)
2 7
k 8.0e−12 k 2.3e−10
3 8
k 8.0e−12 k 1.0e−11
4 9
k 1.0e−02.exp(−0.39/cosθ) k 2.9e−13
5 10
to photochemical reactions (the ones with term hν), which means that this
reactionsdependonthelight,morepreciselyonthepositionoftheSunrelative
tothehorizon:ink ,k andk theangleθ denotesthesolarzenithangle,which
2 5 7
is the angle of the Sun measured from vertical. The chemical species involved
in the simplified reactions are written in Table 3.
Table 3
The chemical species in the model
u u u u u u u u u u
1 2 3 4 5 6 7 8 9 10
NO NO HC ALD O HNO HO RO OH O(1D)
2 3 3 2 2
Fromboththepracticalandmathematicalpointofview,oneisnaturallyinter-
ested in the existence and qualitative of the solutions to the problem (1)-(4).
The well-posedness of initial boundary value problems for a system more gen-
eralthan(1)isobtainedin[12].Throughout oftherest ofthepaper weassume
existence and uniqueness of classical solution of problem (1)-(4) which means
a function that belongs to C([0,T]×Ω) C1((0,T);C(Ω)) (C(0,T);C2(Ω))
and satisfies the equations (1)-(3) pointwise. Moreover, at the finite difference
T T
approximations in Sections 3, 4 we assume fourth in time and sixth in space
derivatives.
Since we are interested in systems describing chemical concentrations, the
nonnegativity of the solutions has to be preserved. It is proved in [1], that if:
1. u (x,y) ≥ 0;
0
2. R (x,y,u), l = 1,...,L is Lipshitz continuous with respect to the concentra-
l
4
tions u ,u ,...,u and it satisfies the inequality R (x,y,u) ≥ 0, whenever
1 2 L l
u = 0, and u ∈ RL ≡ {u ≥ 0, k = 1,...,L},
l + k
than u ≥ 0 for all (x,y) ∈ Ω and t ∈ [0,T].
It is easily to check that the chemical reactions R (u ,u ,...,u ), l = 1,...,10
l 1 2 10
given by (4) satisfy the point 2. and the solution of problem (1)-(3) with (4)
is nonnegative in time t > 0 if the initial data u (x,y) ≥ 0.
0
3 Central Difference Schemes and Richardson Extrapolation
Inthissection,forclarityexpositionwedescribe theconstructionofthesecond
order CDS for the weakly coupled system of two equations
∂u ∂2u ∂2u ∂u ∂u
−a(x,y) −b(x,y) +c(x,y) +d(x,y) = r(x,y,t,u,v),
∂t ∂x2 ∂y2 ∂x ∂y
(5a)
∂v ∂2v ∂2v ∂v ∂v
−e(x,y) −f(x,y) +g(x,y) +h(x,y) = s(x,y,t,u,v),
∂t ∂x2 ∂y2 ∂x ∂y
(5b)
defined on the cylindric domain Q = Ω×(0,T], where Ω ⊂ R2 is a bounded
T
domain with Lipshitz boundary. The nonlinear functions r and s are suffi-
ciently smooth of their arguments. The coefficients a(x,y), b(x,y), e(x,y) and
f(x,y) are positive in Ω. We consider Dirichlet boundary conditions
¯ ¯¯
u(x,y,t) = φ(x,y,t), v(x,y,t) = φ(x,y,t), (x,y,t) ∈ ∂Ω×(0,T] (6)
and initial conditions
¯ ¯¯
u(x,y,0) = ψ(x,y), v(x,y,0) = ψ(x,y), (x,y) ∈ Ω, (7)
¯ ¯¯ ¯ ¯¯
where φ, φ, ψ and ψ are given and smooth data and compatibility of the
boundary and initial data is ensured.
Let for simplicity the domain Ω is a rectangle Ω = [0,X] × [0,Y]. We in-
troduce uniform meshes in the following way: ω = {x = ih , i =
h,x i x
0,1,...,M , h = X/M }, ω = {y = jh , j = 0,1,...,M , h =
x x x h,y j y y y
Y/M } and then Ω = ω × ω , Ω = Ω ∪ ∂Ω , where Ω consist of all
y h h,x h,y h h h h
interior mesh points and ∂Ω - of all boundary mesh points.
h
We will used the index pair (i,j) to represent the mesh point (x ,y ) and
i j
define
u = u(x ,y ,t), v = v(x ,y ,t), r = r(x ,y ,t,u ,v ), ect.
i,j i j i,j i j i,j i j i,j i,j
5
For w = u,v we introduce the central difference operators
δ w = (w −w )/(2h ), δ2w = (w −2w +w )/h2,
x i,j i+1,j i−1,j x x i,j i+1,j i,j i−1,j x
(8)
δ w = (w −w )/(2h ), δ2w = (w −2w +w )/h2.
y i,j i+1,j i−1,j y y i,j i,j+1 i,j i,j−1 y
3.1 Second-order space semidiscretization
Application of the difference operators (8) into the system (5) for every point
(i,j) ∈ Ω leads to
h
∂u
−a δ2u −b δ2u +c δ u +d δ u +χ =r ,
∂t(cid:12) i,j x i,j i,j y i,j i,j x i,j i,j y i,j i,j,1 i,j
(cid:12)(xi,yj)
(cid:12)
∂v(cid:12)
(cid:12) −e δ2v −f δ2v +g δ v +h δ v +χ =s ,
∂t(cid:12) i,j x i,j i,j y i,j i,j x i,j i,j y i,j i,j,2 i,j
(cid:12)(xi,yj)
(cid:12)
(cid:12)
(cid:12)
where the truncation errors χ and χ are
i,j,1 i,j,2
χ = h2x 2c∂3u −a∂4u + h2y 2d∂3u −b∂4u +O(h4 +h4),
i,j,1 12 ∂x3 ∂x4 i,j 12 ∂y3 ∂y4 i,j x y (9)
χ = h2x (cid:16)2g∂3v −e∂4v(cid:17) + h2y (cid:16)2h∂3v −f∂4v(cid:17) +O(h4 +h4).
i,j,2 12 ∂x3 ∂x4 i,j 12 ∂y3 ∂y4 i,j x y
(cid:16) (cid:17) (cid:16) (cid:17)
After dropping the truncation error terms a semi-discrete second-order central
difference approximation of (5) is obtained:
∂uh −a δ2uh −b δ2uh +c δ uh +d δ uh = rh ,
∂t (xi,yj) i,j x i,j i,j y i,j i,j x i,j i,j y i,j i,j (10)
(cid:12)
∂vh(cid:12) −e δ2vh −f δ2vh +g δ vh +h δ vh = sh ,
∂t (cid:12)(xi,yj) i,j x i,j i,j y i,j i,j x i,j i,j y i,j i,j
(cid:12)
(cid:12)
(cid:12)
where for (i,j) ∈ Ω
h
uh ≈ u(x ,y ,t), vh ≈ v(x ,y ,t),
i,j i j i,j i j
rh ≈ r(x ,y ,t,uh ,vh ), sh ≈ s(x ,y ,t,uh ,vh ).
i,j i j i,j i,j i,j i j i,j i,j
Now we introduce the matrix representation for the system (10). We order
the mesh points lexicographically from left to right in x direction and from
the bottom to the top in y direction. Excluding the boundary mesh points
(i,j) ∈ ∂Ω ,forj = 1,2,...,M −1wedefinethefollowing(M −1)dimensional
h y x
vectors:
6
Uh = uh ,uh ,...,uh , Vh = vh ,vh ,...,vh ,
j 1,j 2,j Mx−1,j j 1,j 2,j Mx−1,j
R (Uh,Vh) =(cid:16)(R ,R ,...,R (cid:17) ), S (Uh(cid:16),Vh) = (S , S ,..(cid:17).,S )
j j j 1,j 2,j Mx−1,j j j j 1,j 2,j Mx−1,j
and then
T T
U = Uh,Uh,...,Uh , V = Vh,Vh,...,Vh ,
1 2 My−1 1 2 My−1
(cid:16) (cid:17)T (cid:16) T(cid:17)
R = R ,R ,...,R , S = S ,S ,...,S .
1 2 My−1 1 2 My−1
(cid:16) (cid:17) (cid:16) (cid:17)
We then rewrite the system (10) as a system of ordinary differential equations
d
U +P¯U =R+Φ¯, t ∈ (0,T], (11)
dt
d
V +P¯¯V =S +Φ¯¯, t ∈ (0,T] (12)
dt
¯ ¯¯
with initial conditions U(0) and V(0) obtaining from ψ and ψ for (i,j) ∈ Ω
h
after the reordering. In (11) the matrix P¯ is (M − 1) × (M − 1) block-
y y
tridiagonal matrix P¯ = tridiag(P¯ ,P¯ ,P¯ ) and P¯ , l = k −1,k,k +
k,k−1 k,k k,k+1 k,l
1 are tridiagonal matrixes for l = k and diagonal for l = k ± 1 of order
(M −1)×(M −1). Let for two natural numbers m and M, m < M denote
x x
m : M = m,m+1,...,M and assume that p is a vector with entrances
k,m:M
p = (p ,p ,...,p ) . Then from (10) and (8) the entrances of P¯
k,m:M k,m k,m+1 k,M k,l
are
P¯ = tridiag(p(−1,ε) ,p(0,ε) ,p(1,ε) ) l = k +ε, ε = 0,±1, (13)
k,l k,2:Mx−1 k,2:Mx k,1:Mx−2
where
c(i,j) a(i,j)
(±1,0)
p =± − ,
i,j 2h h2
x x
d(i,j) b(i,j)
(0,±1)
p =± − , (14)
i,j 2h h2
y y
a(i,j) b(i,j)
(0,0)
p =2 +2 .
i,j h2 h2
x y
Replacing a ↔ e, b ↔ f, c ↔ g and d ↔ h in a similar way we obtain the
entrances of the matrix P¯¯.
The vectors Φ¯ and Φ¯¯ in (11)-(12) are associated with the boundary functions
and also depend on time t.
7
3.2 Full discretization
For discretization in time the so called θ-weight method is used. Let ω =
τ
{t = nτ, n = 0,1,...,N, τ = T/N} be uniform mesh in time with time
n
step τ. Then the weight θ-discretization of (11), (12) may be written in the
following way:,
Un+1 −Un
+P¯Un,θ=Rn,θ +Φ¯n,θ, t ∈ (0,T), (15)
τ
Vn+1 −Vn
+P¯¯Vn,θ=Sn,θ +Φ¯¯n,θ, t ∈ (0,T),
τ
where Zn,θ = θZn+1 + (1 − θ)Zn for Z = U,V,R,S,Φ¯,Φ¯¯, Zn ≈ Z(t ) and
n
0 ≤ θ ≤ 1, n = 0,1,...N − 1. For θ = 1 one obtain the fully implicit finite
difference scheme, for θ = 0 - explicit and for θ = 1/2 - the Crank-Nicolson
scheme. The last case has an advantage that the scheme is of second order
in time and as we want to derive schemes of higher order, in the numerical
experiments we use mainly θ = 1/2.
For θ > 0 the finite difference schemes requires solving of nonlinear alge-
braic systems. We briefly discuss the application of the Newton method on
the problem (15). To apply the classical Newton method the system (15) is
rewritten in the form Υ(W) = 0, where W = [UT,VT]T is a vector of length
0
2(M −1)(M −1).We set Wn+1 as initialguess on the new time layer t = t
x y n+1
tobethenumerical solutionontheprevious time layer t = t .Then tofindthe
n
solution on t = t the iterative process with appropriate stopping criteria is
n+1
used:
k k k
Υ′(Wn+1) ∆= −Υ(Wn+1),
 (16)
Wk+n1+1=Wkn+1 + ∆k .
Here ∆k is a vector of the increments and the Jacobian matrix Υ′(Wkn+1) for
θ = 1/2 is
k ∂Υ 1I + 1P¯ − 1∂R 1∂R
Υ′(Wn+1) = = τ 2 2∂U 2∂V , (17)
∂W  1∂S 1I + 1P¯¯ − 1∂S (cid:12)(cid:12)
 2∂U τ 2 2∂V (cid:12)(cid:12)(U,V)=(Uk,Vk)
 (cid:12)
(cid:12)
(cid:12)
where I is the identity matrix and P, P - as defined by (13), (14). In the
numerical experiments to solve the first line in (16) which is a linear system of
2(M −1)(M −1) equations we use the so called inexact Newton method [3],
x y
i.e. we solve this system approximately using the MatLab function bicgstab(l)
(biconjugate gradients stabilized (l) method) that gives better results for our
examples in sense of convergence of the inner iterations and the CPU time.
8
3.3 Richardson extrapolation
Richardson extrapolation is a powerful computational tool which can success-
fully be used in the efforts to improve the accuracy of the of the approximate
solutions of the systems of partial differential equations (PDEs) obtained by
finite difference methods.
Therefor, another way for obtaining the difference schemes of higher order is
to use the Richardson extrapolation method. The main idea [11] is to solve
the difference scheme on two or more consecutive meshes and then to combine
the obtained numerical solutions with appropriate weights. Let assume that
h = h = handforthenumericalsolutiononthen-thtimelayerthefollowing
x y
expression is true:
Uτ = Un = u(x ,y ,tn)+C hσ +χ(h,τ), (x ,y ,t ) ∈ Ω , (18)
h (i,j) i j 1 i j n h,τ
where function χ(h,τ) is a remainder term and C does not depend on h , h
1 x y
and τ. If we want to eliminate the term C hσ, we do the following steps:
1
• solve the difference scheme on two consecutive meshes: coarse one Ω and
h,τ
fine one Ω and let the corresponding numerical solutions be Uτ and
h/2,τ h
Uτ ;
h/2
• find the weights γ and γ from the system
1 2
γ +γ =1 (19)
1 2
γ
2
γ + =0
1
2σ
• obtain a new numerical solution on the coarse mesh
U = γ Uτ +γ Uτ (x ,y ,t ) ∈ Ω .
extr 1 h 2 h/2 i j n h,τ
From (19) we have for the case of central Crank-Nicolson Scheme (σ = 2) that
the coefficients for the Richardson extrapolation are
γ = −1/3 γ = 4/3. (20)
1 2
In the case of CFDS and Richardson Extrapolation (σ = 4) the corresponding
weight coefficients are
γ = −1/15 γ = 16/15. (21)
1 2
If in (18) the more detailed analysis of the LTE is done, then the prolongation
of the idea of space-time Richardson extrapolation [13] can be applied.
9
4 Compact Difference Schemes
In this section, just for clarity we describe the construction of the CFDS again
for the system of two equations (2).
4.1 Space discretization
InordertoeliminatethetermsofO(h2+h2)in(9)wedifferentiatetheequation
x y
(5a) twice with respect to x obtaining expressions for ∂3u, ∂4u, and twice with
∂x3 ∂x4
respect to y for ∂3u, ∂4u.
∂y3 ∂y4
Let
˜
a˜ = (c +2δ a )/a , b = (d +2δ b )/b , (i,j) ∈ Ω .
i,j i,j x i,j i,j i,j i,j y i,j i,j h
Let also
h2 h2
α =a + x δ2a −a˜ (δ a −c )−2δ c + y δ2a −˜b δ a ,
i,j i,j 12 x i,j i,j x i,j i,j x i,j 12 y i,j i,j y i,j
(cid:16) (cid:17) (cid:16) (cid:17)
h2 h2
β =b + x δ2b −a˜ δ b + y δ2b −˜b (δ b −d )−2δ d ,
i,j i,j 12 x i,j i,j x i,j 12 y i,j i,j y i,j i,j y i,j
(cid:16) (cid:17) (cid:16) (cid:17)
h2 h2
α˜ =c + x δ2c −a˜ δ c + y δ2c −˜b δ c ,
i,j i,j 12 x i,j i,j x i,j 12 y i,j i,j y i,j
(cid:16) (cid:17) (cid:16) (cid:17)
h2 h2
β˜ =d + x δ2d −a˜ δ d + y δ2d −˜b δ d ,
i,j i,j 12 x i,j i,j x i,j 12 y i,j i,j y i,j
(cid:16) (cid:17) (cid:16) (cid:17)
and
h2 h2 h2 h2
θ = yc − x(2δ b −a˜ b ), θ˜ = xd − y(2δ a −˜b a ),
i,j i,j x i,j i,j i,j i,j i,j y i,j i,j i,j
12 12 12 12
h2 h2 h2 h2
γ = xb + ya , γ˜ = x(2δ −a˜ d )+ y(2δ c −˜b c ).
i,j i,j i,j i,j x i,j i,j y i,j i,j i,j
12 12 12 12
Define the following difference operators
lh =−α δ2 −β δ2 +α˜ δ +β˜ δ −γ δ2δ2 +θ δ δ2 +θ˜ δ2δ +γ˜ δ δ
i,j i,j x i,j y i,j x i,j y i,j x y i,j x y i,j x y i,j x y
h2 h2
νh =1+ x(δ2 −a˜ δ )+ y(δ2 −˜b δ ).
i,j 12 x i,j x 12 y i,j y
Applying these operators to (5a) we have
lh u = νh (r −u )+O(h4 +h2h2 +h4). (22)
i,j i,j i,j i,j t,i,j x x y y
10