1 Introduction
In this paper we shall
consider the numerical solution of the following convection-diffusion problem
in the cube
Ω=(0,2π)d:
(1.1)∂U∂t=div(a∇U)+b⋅∇U in Ω for t>0 with U(0)=V,
under periodic boundary conditions,
where the positive definite d×d matrix a(x)=(aij(x))
and the vector b=b(x)=(b1,…,bd) are periodic and smooth.
Equation (1.1)
is a special case of the initial-value problem for the operator equation
(1.2)dUdt=-𝒜U+ℬU for t≥0 with U(0)=V,
where 𝒜 and ℬ represent different
physical processes, in our case
𝒜U=-div(a∇U), ℬU=b⋅∇U,
with 𝒜 representing a slow and
ℬ a fast physical process.
The solution of (1.2)
may be formally expressed as
U(t)=ℰ(t)V=e-t(𝒜-ℬ)V for t≥0.
To discretize such an equation in time,
a common approach is to split
𝒜-ℬ into 𝒜 and -ℬ.
With k a time step one introduces tn=nk and
one may then use
the
second-order
symmetric Strang splitting [8, 7]
on each time interval (tn-1,tn),
(1.3)ℰ(k)=e-k(𝒜-ℬ)≈e12kℬe-k𝒜e12kℬ,
which thus locally involves solutions of
Ut=-𝒜U and Ut=ℬU, see e.g. Hundsdorfer and Verwer
[5] and references therein.
The three exponentials on the right in (1.3) are then approximated
by rational functions of 𝒜 and ℬ, respectively, such
as
the Crank–Nicolson method
r0(k𝒜), with r0(λ)=(1-12λ)/(1+12λ),
for the middle factor. The choice of approximation involving ℬ
is not so obvious.
It has been suggested in the context of numerical weather prediction,
e.g. in Baldauf [1], that time steps
of different length
could be used for the two processes, with shorter time intervals for
the fast process and longer for the slow one.
Also Gassmann and
Herzog [3] discuss the difficulties associated with splitting in
such situations.
In the case of reaction-diffusion equation, see
Estep, Ginting, Ropp, Shadid and Tavener [2] and references therein. Our aim in this work is to discuss
this problem in a somewhat rigorous fashion for our simple model problem.
We note that if 𝒜 and ℬ commute, which holds for (1.1)
when a and b
are independent of x, then
e-k(𝒜-ℬ)=e-k𝒜ekℬ=e12kℬe-k𝒜e12kℬ
so that the error in
(1.3) is zero. When 𝒜 and ℬ do not commute, then
formally, by Taylor expansion,
e-k(𝒜-ℬ)-e-k𝒜ekℬ=O(k2) and
e-k(𝒜-ℬ)-e12kℬe-k𝒜e12kℬ=O(k3).
Error estimates for the time splitting, depending on the regularity of the initial values may be
found in Jahnke and Lubich [6],
Hansen and Ostermann [4] and references therein.
In the method that we study in this paper,
we begin by discretizing (1.1) in the spatial variables.
We let h=2πM, where M is a positive integer, and define
a corresponding uniform mesh
(1.4)Ωh={x=xj=jh:j=(j1,…,jd}, 1≤jl≤M,l=1,…,d}.
For M-periodic vectors u with elements uj, corresponding
to the mesh-point xj,
and with uj+Mel=uj,
we consider the following simple second-order finite difference
approximation of (1.1):
(1.5)dudt=∑i,j=1d∂¯j(a´ij∂iu)+∑j=1d12bj(∂j+∂¯j)u in Ωh for t>0 with u(0)=v.
Here ∂j, ∂¯j are forward and
backward finite difference quotients in the direction
of xj,
a´ij(xl)=aij(xl+12hei), and v the restriction
of V to Ωh.
Problem (1.5)
may be written as a system of ODEs in time,
(1.6)dudt=-Au+Bu for t>0 with u(0)=v,
where the Md×Md matrices A and B
correspond to the differential
operators 𝒜 and ℬ. It is then to (1.6) that we will
apply the splitting approach.
In the one-dimensional case we may take, with h=2πM the mesh-width
and xl=lh,
A=1h2[d(x1)-a(x1.5)0……-a(x0.5)-a(x1.5)d(x2)-a(x2.5)……00-a(x2.5)d(x3)……0⋮⋮⋮⋱⋱-a(xM-0.5)-a(x0.5)00…-a(xM-0.5)d(xM)],
where d(xl)=a(xl+0.5h)+a(xl-0.5h)
(recall a(xM+0.5)=a(x0.5)),
and
B=12h[0b(x1)……-b(x1)-b(x2)0b(x2)…0⋮⋮⋱⋮⋮b(xM-1)b(xM)0…-b(xM)0].
The solution of (1.6), the spatially semidiscrete solution, is
u(t)=E(t)v=e-t(A-B)v,
and we shall see that the error in this approximation is O(h2), under
the appropriate regularity assumptions.
For the time discretization we shall work with basic time
intervals of length k=1N,
where N=2p is an even positive integer,
then
apply the Strang splitting (1.3) on each of the
interval, so that
(1.7)E(k)=e-k(A-B)≈e12kBe-kAe12kB,
and finally
approximate the three exponential factors.
We would like the time discretization
error to match that of the discretization in space. For a
second-order time discretization method, this will require
k=O(h).
For k≤γh2, with γ appropriate, the problem
may be solved by explicit approximations but since we prefer
N to be relatively small, we will consider methods with h and k of the same order.
For the approximation of
e-kA in (1.7) we shall use the Crank–Nicolson method.
Then, in order to approximate e12kB, we would like to use the
forward Euler method on a time
interval of length k1<k.
Assuming for the moment that b is constant and thus
B skew-symmetric, we note that
∥I+k1B∥=(1+k12∥B∥2)12≤1+Ck12h-2.
Here and below, ∥⋅∥ denotes the standard matrix norm
induced by the ℓ2 vector inner product.
Stability therefore holds if k12h-2≤Ck1, or if k1≤Ch2.
Since k should be of the same order as h, this makes it natural to choose k1=k2.
We thus
subdivide the time intervals of length k into N
subintervals of length
k2=kN and apply an explicit forward Euler approximation on
each of these.
As we shall see, this approximation
matches the second-order of the Crank–Nicolson method.
Thus the diffusion part of the equation is approximated
on intervals of length k and the convection part on intervals
of length k2.
We consider thus the time discrete solution at time tn=nk,
(1.8)un=Eknv,where Ek=ℬkr0(kA)ℬk with ℬk=(I+k2B)p.
In fact, in the successive time stepping,
only the matrix ℬ~k=ℬk2=(I+k2B)N is used,
except in the first and last steps.
The method proposed thus replaces the solution at each time
step of
a nonsymmetric problem,
by the solution of a symmetric problem,
plus applications of
an explicit method, but
successively repeated N2=p times before and after the diffusion
approximation, thus covering the interval of length Nk2=k.
We re-emphasize that the splitting is applied to the spatially
semidiscrete problem and not to the continuous problem.
The analysis sketched above is carried out in Section 2.
The analysis
will use discrete Sobolev norms.
After this, in Section 2,
we discuss the case when equation (1.1)
contains a small diffusion coefficient ε.
In this case we are able to show
that if ε≤γh with γ sufficiently small,
then the approximation of
e-kεA can be done by the forward Euler method,
and, with the convection part as before, we have a purely explicit
second-order approximation method.
We close the paper by presenting some numerical illustrations
in Section 4.
2 Basic Error Analysis
For the periodic problem in Ω=(0,2π)d and with
Ωh as in (1.4),
we introduce the
discrete inner product and norm
(v,w)h=hd∑xj∈Ωhvjwj and ∥v∥h=(v,v)h12.
Further, we set
∂ju(x)=1h(u(x+hej)-u(x)) and
∂αu=∂1α1…∂dαdu
for α=(α1,…,αd),
and define the discrete Sobolev norm,
with |α|=α1+…+αd,
∥u∥h,s=(∑|α|≤s∥∂αu∥h2)12 for s≥0.
We shall also use ∂¯ju(x)=1h(u(x)-u(x-hej)).
We define
∥w∥ℂs(R)=max|α|≤ssupx∈R|Dαw(x)|,Dα=(∂∂x1)α1⋯(∂∂xd)αd.
We shall write ℂs for ℂs(Ω), and ℂ for ℂ0.
We note that defining Uh to be the restriction to the mesh
Ωh of a smooth
function U, i.e. by (Uh)j=U(xj), we have
∥Uh∥h,s≤C∥U∥ℂs.
Consider now the matrices A and B in our convection-diffusion problem
(1.2).
They satisfy, with C independent of h,
(2.1)∥Au∥h≤C∥u∥h,2 and ∥Bu∥h≤C∥u∥h,1.
Setting Q(x)={y:|ys-xs|≤h,s=1,…,d},
we have, since the terms in AUh are symmetric difference
quotients of U at the mesh-points xj,
and the terms in (𝒜U)h are the corresponding
derivatives of U at xj,
(2.2)|(AUh)(xj)-(𝒜U)(xj)|≤Ch2∥U∥ℂ4(Q(xj))
and
(2.3)|(BUh)(xj)-(ℬU)(xj)|≤Ch2∥U∥ℂ3(Q(xj)),
expressing, in particular, that
(1.5) is a second-order approximation of
(1.1).
We note that, for |α|=s,
(2.4)∥∂αAu-A∂αu∥h≤C∥u∥h,s+1 and ∥∂αBu-B∂αu∥h≤C∥u∥h,s,
and we may conclude from (2.1) that
(2.5)∥Au∥h,s≤C∥u∥h,s+2 and ∥Bu∥h,s≤C∥u∥h,s+1.
The matrix A is positive semidefinite, with
(2.6)∥v∥h,12≤C((Av,v)h+∥v∥h2).
Further, it follows easily from the definition of A that
(2.7)∥Av∥h≤αh-2∥v∥h,where α=4dλmax(a).
For B=(bil) we have
bil={±bj(xi)if l=i±ej,j=1,…,d,0for other l,
and that then bi+ej,i=-bj(xi+ej). It follows from this
that
B=B0+B1 with B0=12(B-BT),B1=12(B+BT).
Here
B0 is skew-symmetric
and
(2.8)∥B0∥≤h-1β0,β0=∑j=1d∥bj∥ℂ,∥B1∥≤β1=∑j=1d∥∂bj∂xj∥ℂ.
Note also that (B0v,v)=0 for all v.
We begin with the straightforward standard analysis
of the spatially semidiscrete problem
(1.6),
which we include for completeness. We first show
the stability of
the solution operator of (1.6)
in discrete Sobolev norms.
Lemma 2.1.
Let E(t)=e-t(A-B). Then, for any s≥0,
with C=Cs independent of h,
∥E(t)v∥h,s≤eCt∥v∥h,s for t≥0.
Proof.
Let s≥0 and let |α|=s. From (1.6) we find,
for u(t)=E(t)v,
(2.9)(∂αut,∂αu)h+(∂αAu,∂αu)h=(∂αBu,∂αu)h.
Here, by (2.4),
(2.10)|(∂αAu,∂αu)h-(A∂αu,∂αu)h|≤C∥u∥h,s+1∥u∥h,s.
Further, since B=B0+B1 with B0 skew-symmetric,
|(∂αBu,∂αu)h|≤|(B∂αu,∂αu)h|+C∥u∥h,s≤C∥u∥h,s2.
Therefore, by (2.9),
12ddt∥∂αu∥h2+(A∂αu,∂αu)h+∥∂αu∥h2≤C∥u∥h,s+1∥u∥h,s.
Hence, using (2.6) and summing over |α|≤s we find,
with c>0,
ddt∥u∥h,s2+c∥u∥h,s+12≤C∥u∥h,s+1∥u∥h,s≤c∥u∥h,s+12+C∥u∥h,s2,
or, with C=Cs,
(2.11)ddt∥u∥h,s2≤C∥u∥h,s2,
from which the lemma follows.
∎
Note that the special case of e-tA is included for B=0.
As a consequence, we have the following second-order
error estimate.
Theorem 2.1.
We have, for the solutions of (1.6) and (1.1),
with v=Vh,
and C=CT,
∥u(t)-Uh(t)∥h≤Ch2∥V∥ℂ4 for nk≤T<∞.
Proof.
Setting ω=u-Uh, we find
dωdt=-Aω+Bω+ρ in Ωh for t>0 with ω(0)=0,
where ρ=((𝒜U)h-AUh)-((ℬUh)-BUh).
Here, by (2.2) and (2.3),
∥ρ(t)∥h≤∥((𝒜U)h-AUh)(t)∥h+∥((ℬU)h-BUh)(t)∥h≤Ch2∥U(t)∥ℂ4,
and hence
ω(t)=∫0tE(t-s)ρ(s)𝑑s,
where ∥ρ(s)∥h≤Ch2∥U(s)∥ℂ4.
Hence, by Lemma 2.1,
since ∥U(s)∥ℂ4≤C∥V∥ℂ4,
∥ω(t)∥h≤C∫0t∥ρ(s)∥h𝑑s≤Ch2∫0t∥U(s)∥ℂ4𝑑s≤Ch2∥V∥ℂ4.∎
Turning to the analysis of the time discretization, we first
show the stability of etB.
Lemma 2.2.
For any s≥0, we have, with Cs independent of h,
(2.12)∥etBv∥h,s≤eCst∥v∥h,s for t≥0.
Here we may choose C0=β1, as in (2.8).
Proof.
Let s≥0 and let |α|=s.
Then
for the solution of (1.6) with A=0,
(∂αut,∂αu)h=(∂αBu,∂αu)h=(B∂αu,∂αu)h+Qα,
where |Qα|≤Cs∥u∥h,s2.
Since
(Bu,u)h=(B1u,u)h≤β1∥u∥h2,
we conclude that (2.11) holds,
which shows (2.12).
For s=0 we have Qα=0 and hence (2.11) holds with C=2β1
which implies (2.12), with C0=β1.
∎
We now show the following error estimate for the Strang splitting.
Lemma 2.3.
We have, with C independent of h and k,
∥e-k(A-B)v-e12kBe-kAe12kBv∥h≤Ck3∥v∥h,6.
Proof.
Setting
F(k)=e-k(A-B)-e12kBe-kAe12kB
and noting that
F(0)=F′(0)=F′′(0)=0, we may use
Taylor’s formula to obtain
∥F(k)∥=∥(F(k)-F(0)-kF′(0)-12k2F′′(0))v∥h≤16k3sups≤k∥F′′′(s)v∥h.
Here,
for s≤k,
using (2.1) and Lemmas 2.1
and 2.2,
∥F′′′(s)v∥h≤∥(A-B)3e-k(A-B)v∥h+C∑i1+i2+i3=3∥Bi1e12sBAi2e-sABi3e12sBv∥h
≤C(∥v∥h,6+∑i1+i2+i3=3∥v∥h,i1+2i2+i3)≤C∥v∥h,6,
which shows the lemma.
∎
We now turn to the time stepping operator Ek defined in
(1.8) and begin with the following
stability result.
Lemma 2.4.
Let k≤γh. Then, with β0,β1 as in (2.8),
(2.13)∥ℬk∥=∥(I+k2B)p∥≤e12μk,where μ=12(γβ0)2+β1.
Further, Ek=Bkr0(kA)Bk is stable, and
(2.14)∥Ekn∥≤eμT for nk≤T.
Proof.
Since B0 is skew-symmetric,
we have
∥I+k2B0∥2=1+k4∥B0∥2≤1+k4h-2β02≤1+(γβ0)2k2≤e(γβ0)2k2.
Hence
(2.15)∥I+k2B∥≤∥I+k2B0∥+k2∥B1∥≤e12(γβ0)2k2+β1k2≤eμk2,
which shows (2.13)
since 2pk=1. Hence
(2.14) follows
by ∥r0(kA)∥≤1.
∎
We start the analysis of the time discretization error with the following.
Lemma 2.5.
Let M be a square matrix, and
assume ∥esM∥≤C for s≤t. Then
we have, for s≤t,
∥(etM-(I+tM))v∥h≤Ct2∥M2v∥h
and
(2.16)∥(etM-(I+tM+12t2M2))v∥h≤Ct3∥M3v∥h.
If also ∥(I+12sM)-1∥≤C for s≤t, then
(2.17)∥(e-tM-r0(tM))v∥h≤Ct3(∥M3v∥h+∥v∥h).
Proof.
By Taylor expansion we have
etM=I+tM+∫0t(t-s)eMsM2𝑑s,
and hence
∥(etM-(I+tM))v∥h≤C∫0t(t-s)∥M2v∥h𝑑s=12Ct2∥M2v∥h.
Estimate (2.16) follows analogously. For (2.17), we use (2.16)
together with
r0(tM)=I-tM+12t2M2+12∫0t(t-s)2r0′′′(sM)𝑑s,
where
r0′′′(λ)=-32(1+12λ)-4,
to complete the proof.
∎
We now show the following error estimate for ℬk.
Lemma 2.6.
Let k≤γh. Then we have, with C=Cγ,
∥(e12kB-ℬk)v∥h≤Ck3∥v∥h,2.
Proof.
Since pk2=12k, we may write
e12kBv-(1+k2B)pv=∑j=0p-1e(p-j-1)k2B(ek2B-(I+k2B))(I+k2B)jv.
By Lemma 2.5, we have
∥(ek2B-(1+k2B))v∥h≤Ck4∥B2v∥h.
By Lemma 2.2,
∥e(p-j-1)k2B∥≤e12kβ1
for j≤p-1. Using also
(2.15),
we find
∥e12kBv-(1+k2B)pv∥h≤Ck4∑j=0p-1∥B2(I+k2B)jv∥h
≤Ck4∑j=0p-1eμjk2∥B2v∥h≤Ck4pe12μk∥B2v∥h≤Ck3∥v∥h,2
which completes the proof.
∎
We now show the following error estimate for Ek.
Lemma 2.7.
Let k≤γh. Then we have, with C=Cγ,
∥E(k)v-Ekv∥h≤Ck3∥v∥h,6.
Proof.
In view of Lemma 2.3,
it remains to show
∥e12kBe-kAe12kBv-Ekv∥h≤Ck3∥v∥h,6.
We have
e12kBe-kAe12kB-ℬkr0(kA)ℬk=(e12kB-ℬk)e-kAe12kB+ℬk(e-kA-r0(kA))e12kB+ℬkr0(kA)(e12kB-ℬk)
=J1+J2+J3.
Here, by the above lemmas,
∥J1v∥h≤Ck3∥e-kAe12kBv∥h,2≤Ck3∥v∥h,2,
∥J2v∥h≤Ck3(∥A3e12kBv∥h+∥e12kBv∥h)≤Ck3∥e12kBv∥h,6≤Ck3∥v∥h,6,
∥J3v∥h≤Ck3∥v∥h,2,
which completes the proof.
∎
We can now prove the following error estimate:
Theorem 2.2.
Let k≤γh. Then we have
for the solutions of (1.8) and (1.6),
with C=Cγ,T,
∥un-u(nk)∥h≤Ck2∥v∥h,6 for nk≤T.
Proof.
We write
(2.18)un-u(nk)=Eknv-E(nk)v=∑j=0n-1Ekn-1-j(Ek-E(k))E(jk).
Using Lemmas 2.4,
2.7 and 2.1, we obtain, for nk≤T,
∥Eknv-E(nk)v∥h≤C∑j=0n-1∥(Ek-E(k))E(jk)v∥h
≤Ck3∑j=0n-1∥E(jk)v∥h,6≤Cnk3∥v∥h,6≤Ck2∥v∥h,6.∎
Since ∥Vh∥h,6≤C∥V∥ℂ6, we immediately obtain from
Theorems 2.1 and 2.2
the following
total error estimate.
Theorem 2.3.
Let v=Vh and
k≤γh. Then we have
for the solutions of (1.8) and (1.1),
with C=Cγ,T,
∥un-Uh(nk)∥h≤Ch2∥V∥ℂ6 for nk≤T.
3 The Case of a Small Diffusion Coefficient
In this section,
we consider the variant of problem (1.1) with
a small diffusion coefficient ε>0,
∂U∂t=εdiv(a∇U)+b⋅∇U in Ω for t>0 with U(0)=V.
The corresponding semidiscrete system (1.6) may then be written
(3.1)dudt=-εAu+Bu for t>0 with u(0)=v,
where
A and B are as before.
We shall see that (3.1) is
stable, and satisfies an O(h2) error estimate, independently of ε.
Further,
for ε and k small, or more precisely, if
max(ε,k)≤γh and kε≤2h2α,
we will be able to show
an O(k2) estimate for the time discretization
error, even when we use the less
accurate forward Euler method for the A part of the
time stepping operator, and with weaker
regularity requirements than earlier. Also, we do not need to use the
symmetric Strang splitting, and consider now, with
r1(λ)=1-λ,
Un=E~knv,where E~k=r1(kεA)ℬ~k with ℬ~k=ℬk2=(I+k2B)2p.
We note the inverse inequality
h∥u∥h,s≤C∥u∥h,s-1,
and hence
(3.2)ε∥Au∥h,s≤C∥u∥h,s+1 for ε≤γh if γ>0.
As in Section 2, we first attend to the spatially semidiscrete problem.
Lemma 3.1.
Let E~(t)=e-t(εA-B). Then,
for any s≥0, we have,
with C=Cs independent of ε>0 and h>0,
∥E~(t)v∥h,s≤eCt∥v∥h,s for t≥0.
Proof.
Following the steps in the proof of Lemma 2.1, we have,
for u(t)=E~(t)v and
|α|=s,
(3.3)(∂αut,∂αu)h+ε(∂αAu,∂αu)h=(∂αBu,∂αu)h.
Here, as in the proof of Lemma 2.1,
|(∂αBu,∂αu)h|≤C∥u∥h,s2.
Hence, by (3.3) and using (2.10),
12ddt∥∂αu∥h2+ε(A∂αu,∂αu)h+ε∥∂αu∥h2≤Cε∥u∥h,s+1∥u∥h,s+C∥u∥h,s2,
and thus, by (2.6), with c>0,
ddt∥u∥h,s2+εc∥u∥h,s+12≤εc∥u∥h,s+12+C∥u∥h,s2.
This implies (2.11), with C independent of ε
and h, and thus
completes the proof.
∎
In the same way as in Section 2, the stability shows the following
error estimate.
Theorem 3.1.
We have, for the solutions of (3.1) and (1.1),
with C=CT independent of ε,
∥u(t)-Uh(t)∥h≤Ch2∥V∥ℂ4 for nk≤T.
In the analysis of the time discretization we begin with the analogue of Lemma 2.3.
Lemma 3.2.
We have, with C=Cγ independent of h and k,
∥E(k)v-e-kεAekBv∥h≤C(εk2+k3)∥v∥h,3 for ε≤γh.
Proof.
With
F(k)=e-k(εA-B)-e-kεAekB, we
have
F(0)=F′(0)=0 and hence, by
Taylor’s formula,
∥F(k)v∥h=∥(F(k)-F(0)-kF′(0))v∥h≤12k2sups≤k∥F′′(s)v∥h.
Here
F′′(s)=e-s(εA-B)(εA-B)2-e-sεA(ε2A2-2εAB+B2)esB
=e-s(εA-B)(ε2A2-εAB-εBA)-e-sεA(ε2A2-2εAB)esB+(e-s(εA-B)-e-sεAesB)B2
=G1(s)+G2(s)+G3(s).
Using (3.2), (2.5) and the boundedness of the exponentials, we find
∥G1(s)v+G2(s)v∥h≤Cε∥v∥h,3 for s≤k.
Further, G3(0)=0 and hence
∥G3(s)v∥h≤ssupσ≤s∥G3′(σ)v∥h≤Ck∥v∥h,3 for s≤k.
Together these estimates complete the proof of the lemma.
∎
We now turn to the time stepping operator E~k and begin
with its stability.
Lemma 3.3.
If kε≤2h2α, then
(3.4)∥r1(kεA)∥=∥I-kεA∥≤1.
If also k≤γh, then
E~k is stable, or, with μ as in Lemma 2.4,
∥E~kn∥≤eμT for nk≤T.
Proof.
We note that, since A is positive semidefinite,
∥I-kεA∥≤1 if kε∥A∥≤2,
and thus (3.4) holds
by (2.7).
Hence,
by Lemma 2.4,
∥E~k∥≤∥r1(kA)∥∥ℬk∥2≤eμk.∎
We now turn to the error analysis and show the following.
Lemma 3.4.
If max(ε,k)≤γh and kε≤2h2α, we have
∥E(k)v-E~kv∥h≤C(εk2+k3)∥v∥h,3.
Proof.
In view of Lemma 3.2 it remains to show
∥e-kεAekBv-E~kv∥h≤C(εk2+k3)∥v∥h,3.
We first note that by Lemma 2.5,
(3.5)∥(e-kεA-r1(kεA))v∥h≤Cε2k2∥A2v∥h≤Cεk2∥v∥h,3.
We write
e-kεAekB-E~k=(e-kεA-r1(kεA))ekB+r1(kεA)(ekB-ℬ~k)=J1+J2.
Here by (3.5) and Lemma 2.5, and by (3.4)
and Lemma 2.4,
∥J1v∥h≤Cεk2∥ekBv∥h,3≤Cεk2∥v∥h,3 and ∥J2v∥h≤Ck3∥v∥h,2,
which completes the proof
∎
The following is the resulting error estimate.
Theorem 3.2.
If max(ε,k)≤γh and kε≤2h2/α we have,
with C=Cγ,T independent of h,k and ε,
∥un-u(nk)∥h≤C(εk+k2)∥v∥h,3 for nk≤T.
Proof.
Using the analogue of (2.18), we find
∥E~knv-E(nk)v∥h≤C∑j=0n-1∥(E~k-E(k))E(jk)v∥h≤C(εk2+k3)∑j=0n-1∥E(jk)v∥h,3≤CT(εk+k2)∥v∥h,3.∎
As in Section 2, our error estimates in
Theorems 3.1 and 3.2 together show
a total error estimate.
Theorem 3.3.
With v=Vh and for
max(ε,k)≤γh and kε≤2h2α we have,
with C=Cγ,T independent of h,k and ε,
∥un-Uh(nk)∥h≤Ch2∥V∥ℂ4 for nk≤T.
4 Numerical Illustrations
In this section we present some numerical computations
to illustrate our error estimates. We restrict ourselves to
the one-dimensional version of
(1.1),
(4.1)Ut=(aUx)x+bUx for x∈(0,2π),t>0,with U(x,0)=sinx.
As before we shall
choose h=2πM, k=1N,
where M and N are positive integers, and study the
effect of doubling these integers.
We begin with the simple case a=b=1, in which
case the exact solution is
U(x,t)=e-tsin(x+t).
In Table 1 we compile the errors in the numerical solution
at t=1, first
in the spatial discretization, then
when our time stepping method is applied to the semidiscrete solution,
and finally the total error. We use N=4,8,16,32,64, and M=5N,
so that kh=52π=0.7958.
The successive ratios of the total errors are given in the
last column and confirm the second-order convergence
estimates resulting from Theorems 2.1–2.3.
Table 1
Numerical results for constant coefficients. Here h=2πM, k=1N.
| M | N | ∥(u-Uh)(⋅,1)∥h | ∥uN-u(⋅,1)∥h | ∥uN-Uh(⋅,1)∥h | Ratio |
| 20 | 4 | 0.01670 | 0.01199 | 0.02494 | |
| 40 | 8 | 0.00423 | 0.00300 | 0.00621 | 4.01 |
| 80 | 16 | 0.00106 | 0.00075 | 0.00155 | 4.01 |
| 160 | 32 | 0.00027 | 0.00019 | 0.00039 | 3.97 |
| 320 | 64 | 0.00007 | 0.00005 | 0.00010 | 3.90 |
We recall that in the case that a and b are constant,
the matrices A and B involved in our method
commute, and consequently the splitting error
given in Lemma 2.3 vanishes. In order to also consider
a situation when this does not happen, we let
a(x)=1+12cosx and b(x)=1+12sinx.
To indicate that the matrices A and B do not commute in this case,
we consider the corresponding continuous operators
𝒜U=-((1+12cosx)Ux)x,ℬU=(1+12sinx)Ux,
and find, after some effort,
(𝒜ℬ-ℬ𝒜)U=-(a(bUx)x)x+b(aUx)xx
=(1-142-sinx2+cosx)Ux-(12+12sinx-14sin2x+cosx)Uxx.
Thus 𝒜 and ℬ do not commute, and therefore neither could A and B.
The exact solution U is taken to be the
semidiscrete solution with M=2560,N=512.
The errors are presented in Table 2. Again we see that
the errors are of second order, which agrees with the error bounds in
Theorems 2.1–2.3.
Table 2
Numerical results for variable coefficients. Here h=2πM, k=1N.
| M | N | ∥(u-Uh)(⋅,1)∥h | ∥uN-u(⋅,1)∥h | ∥uN-Uh(⋅,1)∥h | Ratio |
| 20 | 4 | 0.02641 | 0.01419 | 0.03323 | |
| 40 | 8 | 0.00651 | 0.00356 | 0.00817 | 4.07 |
| 80 | 16 | 0.00163 | 0.00089 | 0.00203 | 4.02 |
| 160 | 32 | 0.00041 | 0.00022 | 0.00051 | 3.98 |
| 320 | 64 | 0.00010 | 0.00005 | 0.00013 | 3.92 |
We finally consider a numerical example for Section 3,
for which we use (4.1) with a=ε=0.01, b=1.
Here
U(x,t)=e-εtsin(x+t),
and un=E~knv with E~k=r1(kA)ℬ~k.
Note that the condition εk≤2h2α, with α=4, now
reduces to ε≤π5h<0.8h, or ε≤π2800=0.0123
for N=64, which is satisfied for our choice of ε. The results are given in Table 3
and agree with the error bounds of Section 3.
Table 3
Numerical results for constant coefficients and with small diffusion. Here a=0.01, b=1, h=2πM, k=1N.
| M | N | ∥(u-Uh)(⋅,1)∥h | ∥uN-u(⋅,1)∥h | ∥uN-Uh(⋅,1)∥h | Ratio |
| 20 | 4 | 0.02872 | 0.05381 | 0.06237 | |
| 40 | 8 | 0.00721 | 0.01365 | 0.01555 | 4.01 |
| 80 | 16 | 0.00180 | 0.00342 | 0.00388 | 4.00 |
| 160 | 32 | 0.00045 | 0.00085 | 0.00097 | 4.00 |
| 320 | 64 | 0.00011 | 0.00021 | 0.00024 | 4.04 |
