Error estimates for splitting methods based on AMF-Runge-Kutta formulas for the time integration of advection diffusion reaction PDEs

0.32 MB·English

by  Severiano Gonzalez Pinto

Error estimates for splitting methods based
on AMF-Runge-Kutta formulas for the time
integration of advection diffusion reaction
PDEs.
5
1
0
2

n
S. Gonzalez-Pinto, D. Hernandez-Abreu and S. Perez-Rodriguez
a
J

2
Departamento de Ana´lisis Matem´atico. Universidad de La Laguna.
1
38071. La Laguna, Spain.

] email: [email protected], [email protected]
A
N
.
h
t
a
m Abstract
[

1 The convergence of a family of AMF-Runge-Kutta methods (in short AMF-RK)
v
forthetimeintegrationofevolutionaryPartialDifferentialEquations(PDEs)ofAd-
1
4 vectionDiffusionReactiontypesemi-discretizedinspaceisconsidered.Themethods
6 arebasedonveryfewinexactNewtonIterationsofAproximateMatrixFactorization
2
splitting-type (AMF) applied to the Implicit Runge-Kutta formulas, which allows
0
. very cheap and inexact implementations of the underlying Runge-Kutta formula.
1
0 Particular AMF-RK methods based on Radau IIA formulas are considered. These
5 methods have given very competitive results when compared with important for-
1
mulas in the literature for multidimensional systems of non-linear parabolic PDE
:
v problems. Uniform bounds for the global time-space errors on semi-linear PDEs
i
X when simultaneously the time step-size and the spatial grid resolution tend to zero
r are derived. Numerical illustrations supporting the theory are presented.
a
Key words: Evolutionary Advection-Diffusion-Reaction Partial Differential
equations, Approximate Matrix Factorization, Runge-Kutta Radau IIA methods,
Finite Differences, Stability and Convergence.
AMS subject classifications: 65M12, 65M15, 65M20.
Corresponding author: S. Gonzalez-Pinto ([email protected])
∗
1 This work has been supported by projects MTM2010-21630-C02-02 and
MTM2013-47318-C2-2
Preprint submitted to Elsevier 13 January 2015
1 Introduction
We consider numerical methods for the time integration of a family of Initial
Value Problems in ODEs
y′(t) = f (t,y (t)), y (0) = u∗ , 0 t t∗, y ,f Rm(h), h 0+,
h h h h 0,h ≤ ≤ h h ∈ →
(1.1)
coming from the spatial semi-discretization of an l dimensional Advection
−
Diffusion Reaction problem in time dependent Partial Differential Equations
(PDEs), with prescribed Boundary Conditions and an Initial Condition. Here
h denotes a small positive parameter associated with the spatial resolution
and usually l = 2,3,... .
The typical PDE problem with Dirichlet boundary conditions is given by (Ω is
a bounded openconnected region inRl,∂Ω its boundaryand isthe gradient
∇
operator)
¯
u (x,t) = (a(x,t)u(x,t))+ (d(x,t) u(x,t))+r(u,x,t),
t
−∇· ∇· ·∇
x Ω, t [0,t∗]; a(x,t) = (a (x,t))l Rl, d¯(x,t) = (d¯(x,t))l Rl,
∈ ∈ j j=1 ∈ j j=1 ∈
u(x,t) = g (x,t), (x,t) ∂Ω [0,t∗]; u(x,0) = g (x), x Ω,
1 2
∈ × ∈
(1.2)
¯
which is assumed to have some diffusion (d (x,t) d > 0, j = 1,...,l),
j 0
≥
namely that it is not of pure hyperbolic type, and it is also assumed that some
adequate spatial discretization based on Finite Differences or Finite Volume is
appliedtoobtainthesystem(1.1).Somestiffnesinthereactionpartr(u,x,t)is
also allowed. The treatment of Systems of PDEs do not involve more difficulty
forour analysis but forsimplicity ofpresentation we prefer to confineourselves
to the case of one PDE.
We denote by u (t) the solution of the PDE problem confined to the spatial
h
grid (or well to the h-space related). It will be tacitly assumed that the PDE
problem admits a smooth solution u(x,t) in the sense that continuous partial
derivatives in all variables up to some order p exist and are continuous and
uniform bounded on Ω [0,t∗] and that u(x,t) is continuous on Ω¯ [0,t∗]
¯ × ×
(Ω = Ω ∂Ω). It is also assumed that the spatial discretization errors
S σ (t) := u′(t) f (t,u (t)), (1.3)
h h − h h
satisfy in the norm considered,
σ (t) C hr, (C 0, r > 0), 0 t t∗, h 0. (1.4)
h
k k ≤ ≥ ≤ ≤ →
In general C, C′ or C∗ will refer to some constants that maybe different at
each occurrence but that all of them are independent of h 0 and from the
→
2
time-stepsize τ 0. The vector norm used is arbitrary as long as it is defined
→
for vectors of any dimension. For square matrices the norm used is the induced
operator norm, A = sup Av / v .
k k v6=0k k k k
In spite of most of our results apply in general, we will provide specific results
for weighted Euclidean norms of type
(v )N = N−1/2 (v )N .
k j j=1k k j j=1k2
It should be noted that in this case we have for any square matrix A that,
A = A , A RN,N, N = 1,2,3,....
2
k k k k ∀ ∈
We assume some natural splitting for f (directional or other),
h
d
f (t,y) = f (t,y), (1.5)
h j,h
j=1
X
which provides some natural splitting for the Jacobian matrix at the current
point (t ,y ),
n n
d ∂f (t ,y ) ∂f (t ,y )
h n n j,h n n
J = J , J := , J := . (1.6)
h j,h h j,h
∂y ∂y
j=1
X
This goal of the paper is to analyze the convergence order of the Method of
Lines (MoL) approach for time-dependent PDEs of Advection Reaction Dif-
fusion PDEs, with the main focuss on the time integration of the large ODE
systems resulting of the spatial PDE-semidiscretization, where some stiffness
is assumed (parabolic dominant problems with stiff reaction terms) and the
time integrators are based on very few iterations of splitting type (Approxi-
mate Matrix Factorizationand Newton-type schemes) applied to highly stable
Implicit Runge-Kutta methods. It should be remarked that the underlying
Implicit Runge-Kutta method is never solved up to convergence, hence the
convergence study does not follows from the results collected in classical ref-
erences about finite difference methods such as [14,4,10,17,13,9]. The kind of
approach to be considered here has interest since it is easily applicable to
general systems of PDEs as we will see later on and it is reasonably cheap
for non-linear problems in general (although we give convergence results for
semilinear problems only) when some splitting of the function f and its Ja-
h
cobian is available and the split terms can be handled efficiently. In particular
a method based on three AMF-iterations of the two-stage Radau IIA method
[1] has shown to be competitive [7] when compared with some standard PDE-
solvers such as VODPK [2,3] in some interesting non-linear diffusion reaction
problems widely considered in the literature. We also present two new meth-
ods based on the 2-stage Radau IIA, by performing just one or two iterations
3
of splitting type, respectively. The method based on two iterations is one of
the very few one-step methods of splitting type we have seen in the literature
that has order three in PDE-sense for the time integration.
The rest of the paper is organized as follows. In section 2 we introduce the
AMF -RK methods, and special attention is paid to some methods based on
q
Radau IIA formulas. In section 3, the convergence for semilinear PDEs is
studied in detail. The local and global errors are studied for the AMF -RK
q
splitting methods based on some general Runge-Kutta methods. Section 4
is devoted to some applications of the convergence results to 2D and 3D-
parabolic PDEs.
Henceforth, for simplicity in the notations, we omit in many cases the h-
dependence of some vectors such as f , f and of some matrices such as
h j,h
J and J (j = 1,...,d). It should be clear from the context which ones
h j,h
are h-dependent. Besides, we will refer to the identity matrix as I when its
dimension is clear from the context.
2 AMF-IRK methods
For the integration of the ODEs (1.1), we consider as a first step an implicit s-
stageRunge-Kuttamethodwithanonsingular coefficient matrixA = (a )s
ij i,j=1
and a weight vector b = (b )s . The method is given by the compact formula-
j j=1
tion(below denotestheKronecker product ofmatrices A B = (a B), A =
ij
⊗ ⊗
(a ), B = (b ))
ij ij
Y = e y +τ(A I )F(Y ),
n n m n
⊗ ⊗
y = ̟y +(ßT I )Y ,
n+1 n m n
⊗
c (c )s := Ae, e = (1,...,1)T Rs, ßT := bTA−1, ̟ = 1 ßTe,
≡ j j=1 ∈ −
Y = (Y )s Rms, F(Y ) = (f(t +τc ,Y ))s Rms.
n n,j j=1 ∈ n n j n,j j=1 ∈
(2.1)
It should be noted that we have replaced the usual formulation at the stepping
point y = y +τ(bT I )F(Y ) by the equivalent in (2.1), which has some
n+1 n m n
⊗
computational advantages for stiff problems when the algebraic system for the
stages is not exactly solved.
A typical Quasi-Newton iteration to solve the stage equations above is given
by (below, J = ∂f/∂y(t ,y ) is the exact Jacobian at the step-point (t ,y )),
n n n n
[I A τJ]∆ν = Dν−1, Yν = Yν−1 +∆ν, ν = 1,2,..., (2.2)
ms − ⊗ n n n
4
where
Dν−1 D(t ,τ,y ,Yν−1) := e y Yν−1 +τ(A I )F(Yν−1). (2.3)
n ≡ n n n ⊗ n − n ⊗ m n
A cheaper iteration of Newton-type when the matrix A has a multipoint spec-
trum has been considered in [6,12] (denoted as Single-Newton iteration)
[I T τJ]∆ν = Dν−1, Yν = Yν−1 +∆ν, ν = 1,2,...,q (2.4)
ms − ν ⊗ n n n
where
T = γS (I L )−1S−1, γ > 0,
ν ν − ν ν
S Rs,s are regular matrices and (2.5)
ν
∈
L Rs,s are strictly lower triangular matrices.
ν
∈
After some simple manipulations, by using standard properties of the Kro-
necker product, this iteration can be rewritten in the equivalent form,
[I (I γτJ)]Eν = ((I L )S−1 I )Dν−1 +(L I )Eν,
s ⊗ m − s − ν ν ⊗ m n ν ⊗ m (2.6)
Yν = Yν−1 +(S I )Eν, ν = 1,2,...,q.
n n ν ⊗ m
To reduce the algebra cost, we use the Approximate Matrix Factorization [8]
in short AMF, with J J and J J given in (1.6),
h j j,h
≡ ≡
d
Π := (I γτJ ) = (I γτJ)+ (τ2), (2.7)
d m j m
− − O
j=1
Y
and replace in (2.6) (I γτJ) by Π , which yields the AMF -RK method
m d q
−
based on the underlying Runge-Kutta method
[I Π ]Eν = ((I L )S−1 I )Dν−1 +(L I )Eν,
s ⊗ d s − ν ν ⊗ m n ν ⊗ m
Yν = Yν−1 +(S I )Eν, ν = 1,2,...,q
n n ν ⊗ m (2.8)
Y0 = e y (Predictor)
n ⊗ n
y = ̟y +(ßT I )Yq (Corrector).
n+1 n ⊗ m n
Our starting point for the convergence analysis in the next section takes into
account that the AMF -RK method can be rewritten in the equivalent form [5]
q
[I I T τP](Yν Yν−1) = D(t ,τ,y ,Yν−1), 1 ν q
⊗ − ν ⊗ n − n n n n ≤ ≤ (2.9)
Y0 = e y , y = ̟y +(ßT I )Yq,
n ⊗ n n+1 n ⊗ m n
5
where the matrix P plays a primary role
P := (γτ)−1(I Π )
d
−
= J +( γτ) J J +( γτ)2 J J J +...+( γτ)d−1J J J .
j k j k l 1 2 d
− − − ···
j<k j<k<l
X X
(2.10)
2.1 AMF -RK methods based on the 2 stage Radau IIA formula
q
We are going to deserve special attention to AMF -RK methods based on the
q
2 stage Radau IIA formula [1]. This formula has coefficient Butcher tableau
given by
1/3 5/12 1/12
−
c A
1 3/4 1/4
bT ≡
3/4 1/4
This is a collocation method (stage order is two) possessing good stability
properties, such asL-stability(i.e. A-stabilityplus R( ) = 0, withR(z) being
∞
thelinear stabilityfunctionofthemethod), andhasorder ofconvergence three
(in ODE sense), not only on non-stiff problems but also in many kinds of stiff
problems [4]. These properties for the underlying Runge-Kutta method are
convenient, since the family of ODEs (1.1) involves stiffness in most of cases,
due to the diffusion terms and possibly to the reaction part, and it is expected
that the methods to be built on inherit part of the good properties of the
original Runge-Kutta method.
The next three AMF -Rad methods have coefficient matrices (L , S and T )
q ν ν ν
and eigenvalue γ of the form
1 s 0 0
T = γS (I L )−1S−1, S = ν , L = , γ = det(A) = 1/√6.
ν ν 2− ν ν ν   ν  
0 1 l 0
ν q
   
    (2.11)
AMF -Radwas derived in [5] by looking for goodstability properties and order
1
two (ODE sense). In particular the method is A(π/2)-stable for a 2-splitting
(see in Definition 1 below, the concept of stability for a d-splitting), A(0)-
stable for any d-splitting and has stability wedges close to θ = π/(2(d 1))
d
−
for d = 3,4. The method is based on one iteration (q = 1) and was required
to fulfil (A T )c = 0 and it has coefficients given by
1
−
3+2√6 3
s = , l = ( 12+5√6). (2.12)
1 1
− 9 4 −
6
AMF -Radwas derived in [5] by looking for goodstability properties and order
2
three (ODE sense). The method is A(π/2)-stable for a 2-splitting, A(0)-stable
for any d-splitting and A(π/6)-stable for d = 3,4. The method is based on
two iterations (q = 2) and their matrices T and T were required to satisfy
1 2
(A T )c = 0andeTT−1(A T ) = 0T, eT = (0,1),respectively. Itscoefficients
− 1 2 2 − 2 2
are uniquely given by
3+2√6 3
s = , l = ( 12+5√6)
1 1
− 9 4 −
(2.13)
5 2√6 3√6
s = − , l = .
2 2
9 4
AMF -Rad was derived in [12,7] by looking for good stability properties and
3
order three (ODE sense). The method is A(π/2)-stable for a 2-splitting, A(0)-
stable for any d-splitting and close to A(θ )-stable for d = 3,4 with θ =
d d
π/(2(d 1)). The method is based on three iterations (q = 3) and their
−
matrices T = T = T = T were required to satisfy eTT−1(A T) = 0T. Its
1 2 3 2 −
coefficients are uniquely given by
5 2√6 3√6
s = s = s = − , l = l = l = . (2.14)
1 2 3 1 2 3
9 4
In [7], a variable-stepsize integrator based on the AMF -Rad method was suc-
3
cessfully tested on several interesting 2D and 3D advection diffusion reaction
PDEs by exhibiting good performances in comparison with state-of-the-art
codes like VODPK [2,3] and RKC [16,19] and its implicit-explicit counterpart,
IRKC [15,18]. The other two methods, AMF -Rad (q = 1,2), were introduced
q
later [5] after carefully analyzing the PDE errors on semilinear problems and
with the purpose of reducing the number of iterations w.r.t. AMF -Rad.
3
3 Convergence for semilinear problems
For our convergence analysis we consider AMF -RK methods applied to the
q
ODE problems coming from the spatial discretizations of semilinear PDE
problems of type (1.2) where the advection and diffusion vectors a(x,t) and
¯
d(x,t) are both constant and the reaction part has the form
r(u,x,t) = κu+g(x,t), κ being a constant, x Ω Rl. (3.1)
∈ ⊆
In this way, the ODE systems have the form
y′(t) = f (t,y ) := J y (t)+g (t), y (0) = u∗ , h 0+,
h h h h h h h 0,h → (3.2)
J = d J , t [0,t∗].
h j=1 j,h ∈
P
7
Here, theexactsolutionofthePDEconfinedtothespatial gridu (t) = u(x,t),
h
is assumed to satisfy (1.3) and (1.4). Thus, we focus on the global errors of
the MoL approach, where the spatial discretization is carried out first by
using finite differences (or finite volumes) and then the time discretization
is performed by using AMF -RK methods. It is important to remark that we
q
will not pursue the details of the spatial semidiscretizations but rather it is
assumed that the spatial semidiscretizations are stable and provides spatial
discretization errors satisfying (1.4). We shall provide uniform bounds for the
global errors of the MoL approach (y (t) henceforth denotes the numerical
h
solution of the MoL approach) in the sense
ǫ := u (t ) y (t ) = (τ)p1 + (hατp2), h 0+,τ 0+, (3.3)
n,h h n h n
− O O → →
which is meant that there exist constants C , C , p , p , α (all of them inde-
1 2 1 2
pendent on h and τ) so that in the norm considered,
ǫ C τp1 +C hατp2, h 0+,τ 0+ holds.
n,h 1 2
k k ≤ → →
Inourconvergence analysisweneedthatallthematricesJ pairwiseconmute
j,h
and that they can be brought to the following decomposition (it has some
resemblance with the Jordan’s decomposition, but it is a little more general)
J = Θ Λ Θ−1, Cond(Θ ) := Θ Θ−1 C, h 0+, 1 j d,
j,h h j,h h h k hk·k h k ≤ → ≤ ≤
Λ = BlockDiag(Λ(1),Λ(2),...,Λ(ϑh)), Λ(l) = λ(l)I +E(l), Re λ(l) 0,
j,h j,h j,h j,h j,h j,h h j,h ≤
dim(E(l)) N, E(l) C′, l = 1,2,...,ϑ (h 0+).
h ≤ k h k∞ ≤ h →
(l)
E are all of them strictly lower triangular matrices.
h
(3.4)
Another important approach for the convergence analysis of the MoL method
(mainly concerned with the time integration) is based on the pseudo-spectra
analysis of the matrix J [13] and the related matrices J . That analysis is
h j,h
of more general scope but it is much more difficult to make and as we will see
below, our analysis is enough for some interesting kind of semilinear problems
and it is expected that the results extend to most of the non-linear problems
of parabolic dominant type.
Next, we consider a standard 3D semilinear-PDEs problem where the assump-
tions in (3.4) are fulfilled.
3.1 An example
Consider the semilinear PDE-problem (1.2) with x Ω = (0,1)3, with con-
stant vectors, a(x,t) = (a )3 , d¯(x,t) = (d¯)3 , d∈¯ > 0 (j = 1,2,3) and
j j=1 j j=1 j
8
r(x,u,t) as in (3.1). Consider the spatial semidiscretization by using second
order central differences and spatial resolution h = 1/(N + 1). This yields a
semilinear ODE systems of dimension m = N3 of the form (3.2) for d = 3.
The matrices J are given by
j,h
J = I I , J = I I , J = I I
1,h N N 1 2,h N 2 N 3,h 3 N N
⊗ ⊗T ⊗T ⊗ T ⊗ ⊗
= Tridiag(α ,δ ,β) RN,N, l = 1,2,3,
l l l l
T ∈
α = h−2(d¯ 2−1ha ), β = h−2(d¯ +2−1ha ), δ = h−2( 2d¯ +h2κ),
l l l l l l l l
− −
(3.5)
and the vector g (t) includes the reaction part g(x,t) plus the boundary con-
h
ditions. It is straightforward to see that the J pairwise commute. Moreover,
l,h
by assuming Cell-P´eclet numbers [9, p. 67, formula (3.42) ]
¯
h a /d < 2, l = 1,2,3,
l l
| |
from [11, section 2] it follows that their spectral decomposition has the form
= Tridiag(α ,δ ,β ) = V Λ V−1, V = D U, l = 1,2,3,
Tl l l l l l l l l
kπ
Λ = Diag(λ )N , λ = δ +2 α β cos ,
l l,k k=1 l,k l l l N +1
q (3.6)
kjπ
U = ( 2 )1/2 sin is an orthogonal matrix and
N+1 N +1!k=1,N
j=1,N
N
D = (N+1)1/2Diag (α /β )k/2 .
l 2 l l k=1
(cid:16) (cid:17)
From here we conclude that all the matrices can be brought to the spectral
decomposition in (3.4) having negative eigenvalues and with matrix Θ =
h
V V V . Observe that
3 2 1
⊗ ⊗
Θ Θ−1 = 3 V V−1 = 3 D D−1
k hk2k h k2 l=1k lk2k l k2 l=1k lk2k l k2
Q3 2d¯ +h a N/2 Q 3 2d¯ +h a 1/(2h)
l l l l
= | | | |
¯ ¯
lY=1 2dl −h|al|! ≤ lY=1 2dl −h|al|!
3 a
l
exp | | as h 0.
¯
≃ l=1 2dl! →
X
3.2 Analysis of the Truncation Errors
TheAMF -RKmethodappliedonproblem(1.1) canbeexpressed inthesimple
q
one-step formaty = φ (t ,y ,τ), n 0.Thus, thetime-space globalerrors
n+1 f n n
≥
9
ǫ = u (t ) y satisfy
n h n n
−
ǫ := u (t ) φ (t ,y ,τ)
n+1 h n+1 f n n
−
= (u (t ) φ (t ,u (t ),τ))+(φ (t ,u (t ),τ) φ (t ,y ,τ))
h n+1 f n h n f n h n f n n
− −
= l(t ,τ,h)+[∂φ /∂y] (u (t ) y ),
n f n h n n
−
where
1 ∂φ
f
[∂φ /∂y] = (t ,u (t )+(θ 1)ǫ ,τ)dθ,
f n n h n n
∂y −
Z0
and the time-space local errors are defined by
l l(t ,τ,h) := u (t ) φ (t ,u (t ),τ). (3.7)
n n h n+1 f n h n
≡ −
Then, we have for the time-space global errors ǫ the recurrence
n
ǫ = [∂φ /∂y] ǫ +l , n = 0,1,2,...,t∗/τ 1. (3.8)
n+1 f n n n
· −
In order to get a better understanding of the latter recurrence, we next intro-
duce the following matrix operators (P is defined in (2.10))
Q = (I I T τP)−1, M = Q (A τJ T τP), ν 1, Q = I.
ν ν ν ν ν 0
⊗ − ⊗ ⊗ − ⊗ ≥
(3.9)
Lemma 1 The time-space global errors provided by the AMF -RK method
q
when applied to the problem (3.2) satisfy the recurrence
ǫ = R (τJ,τP) ǫ +l , n = 0,1,2,...,t∗/τ 1, (3.10)
n+1 q n n
· −
where l stands for the time-space local error defined in (3.7) and
n
1 j
R (τJ,τP) = ̟I +(ßT I) Q + ( M )Q (e I), (3.11)
q q i j−1
⊗   ⊗
j=q i=q
X Y
 
with Q ,M (ν 1) given by (3.9). Moreover, the function R (τJ,τP) fulfils
ν ν q
≥
1 j 1
R (τJ,τP) I = (ßT I) Q + ( M )Q M (c τJ). (3.12)
q q i j−1 i
− ⊗  −  ⊗
j=q i=q i=q
X Y Y
 
Remark 1 It must be observed that commutativity does not hold in general,
thus 1 M M M M . On the other hand, R ( ) can be seen as the
j=q j ≡ q q−1··· 1 q ·
linear stability function of the method. The identity (3.12) for the function
Q
R ( ) I will play a major role in a favourable propagation of the local errors
q
· −
in a similar way as indicated in Lemma 2.3 in [9, p.162].
Proof of Lemma 1. Our first step is to analyze the operator [∂φ /∂y] for
f n
the semilinear problem (3.2). Taking into account that the method is defined
10