Skip to content
BY 4.0 license Open Access Published by De Gruyter July 12, 2018

An Explicit-Implicit Splitting Method for a Convection-Diffusion Problem

Vidar Thomée EMAIL logo and A. S. Vasudeva Murthy


We analyze a second-order accurate finite difference method for a spatially periodic convection-diffusion problem. The method is a time stepping method based on the Strang splitting of the spatially semidiscrete solution, in which the diffusion part uses the Crank–Nicolson method and the convection part the explicit forward Euler approximation on a shorter time interval. When the diffusion coefficient is small, the forward Euler method may be used also for the diffusion term.

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)Ut=div(aU)+bUin Ωfor t>0with 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+Ufor t0with U(0)=V,

where 𝒜 and represent different physical processes, in our case 𝒜U=-div(aU), U=bU, 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(𝒜-)Vfor t0.

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),


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=e12ke-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(𝒜-)-e12ke-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}, 1jlM,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´ijiu)+j=1d12bj(j+¯j)uin Ωhfor t>0with 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+Bufor t>0with 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,


where d(xl)=a(xl+0.5h)+a(xl-0.5h) (recall a(xM+0.5)=a(x0.5)), and


The solution of (1.6), the spatially semidiscrete solution, is


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


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


Here and below, denotes the standard matrix norm induced by the 2 vector inner product. Stability therefore holds if k12h-2Ck1, or if k1Ch2. 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


Further, we set ju(x)=1h(u(x+hej)-u(x)) and αu=1α1dαdu for α=(α1,,αd), and define the discrete Sobolev norm, with |α|=α1++αd,

uh,s=(|α|sαuh2)12for s0.

We shall also use ¯ju(x)=1h(u(x)-u(x-hej)). We define


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 Uhh,sCUs.

Consider now the matrices A and B in our convection-diffusion problem (1.2). They satisfy, with C independent of h,


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,




expressing, in particular, that (1.5) is a second-order approximation of (1.1).

We note that, for |α|=s,


and we may conclude from (2.1) that


The matrix A is positive semidefinite, with


Further, it follows easily from the definition of A that

(2.7)Avhαh-2vh,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+B1with B0=12(B-BT),B1=12(B+BT).

Here B0 is skew-symmetric and


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 s0, with C=Cs independent of h,

E(t)vh,seCtvh,sfor t0.


Let s0 and let |α|=s. From (1.6) we find, for u(t)=E(t)v,


Here, by (2.4),


Further, since B=B0+B1 with B0 skew-symmetric,


Therefore, by (2.9),


Hence, using (2.6) and summing over |α|s we find, with c>0,


or, with C=Cs,


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)hCh2V4for nkT<.


Setting ω=u-Uh, we find

dωdt=-Aω+Bω+ρin Ωhfor t>0with ω(0)=0,

where ρ=((𝒜U)h-AUh)-((Uh)-BUh). Here, by (2.2) and (2.3),


and hence


where ρ(s)hCh2U(s)4. Hence, by Lemma 2.1, since U(s)4CV4,


Turning to the analysis of the time discretization, we first show the stability of etB.

Lemma 2.2.

For any s0, we have, with Cs independent of h,

(2.12)etBvh,seCstvh,sfor t0.

Here we may choose C0=β1, as in (2.8).


Let s0 and let |α|=s. Then for the solution of (1.6) with A=0,


where |Qα|Csuh,s2. Since (Bu,u)h=(B1u,u)hβ1uh2, 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,



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


Here, for sk, using (2.1) and Lemmas 2.1 and 2.2,


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)pe12μk,where μ=12(γβ0)2+β1.

Further, Ek=Bkr0(kA)Bk is stable, and

(2.14)EkneμTfor nkT.


Since B0 is skew-symmetric, we have




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 esMC for st. Then we have, for st,




If also (I+12sM)-1C for st, then



By Taylor expansion we have


and hence


Estimate (2.16) follows analogously. For (2.17), we use (2.16) together with


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γ,



Since pk2=12k, we may write


By Lemma 2.5, we have


By Lemma 2.2, e(p-j-1)k2Be12kβ1 for jp-1. Using also (2.15), we find


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γ,



In view of Lemma 2.3, it remains to show


We have


Here, by the above lemmas,


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)hCk2vh,6for nkT.


We write


Using Lemmas 2.4, 2.7 and 2.1, we obtain, for nkT,


Since Vhh,6CV6, 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)hCh2V6for nkT.

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,

Ut=εdiv(aU)+bUin Ωfor t>0with U(0)=V.

The corresponding semidiscrete system (1.6) may then be written

(3.1)dudt=-εAu+Bufor t>0with 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 huh,sCuh,s-1, and hence

(3.2)εAuh,sCuh,s+1for εγhif γ>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 s0, we have, with C=Cs independent of ε>0 and h>0,

E~(t)vh,seCtvh,sfor t0.


Following the steps in the proof of Lemma 2.1, we have, for u(t)=E~(t)v and |α|=s,


Here, as in the proof of Lemma 2.1, |(αBu,αu)h|Cuh,s2. Hence, by (3.3) and using (2.10),


and thus, by (2.6), with c>0,


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)hCh2V4for nkT.

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εAekBvhC(εk2+k3)vh,3for εγh.


With F(k)=e-k(εA-B)-e-kεAekB, we have F(0)=F(0)=0 and hence, by Taylor’s formula,




Using (3.2), (2.5) and the boundedness of the exponentials, we find

G1(s)v+G2(s)vhCεvh,3for sk.

Further, G3(0)=0 and hence

G3(s)vhssupσsG3(σ)vhCkvh,3for sk.

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


If also kγh, then E~k is stable, or, with μ as in Lemma 2.4,

E~kneμTfor nkT.


We note that, since A is positive semidefinite,

I-kεA1if kεA2,

and thus (3.4) holds by (2.7). Hence, by Lemma 2.4,


We now turn to the error analysis and show the following.

Lemma 3.4.

If max(ε,k)γh and kε2h2α, we have



In view of Lemma 3.2 it remains to show


We first note that by Lemma 2.5,


We write


Here by (3.5) and Lemma 2.5, and by (3.4) and Lemma 2.4,


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)hC(εk+k2)vh,3for nkT.


Using the analogue of (2.18), we find


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)hCh2V4for nkT.

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+bUxfor 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.12.3.

Table 1

Numerical results for constant coefficients. Here h=2πM, k=1N.


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


and find, after some effort,


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.12.3.

Table 2

Numerical results for variable coefficients. Here h=2πM, k=1N.


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 εk2h2α, 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.


Dedicated to Amiya K. Pani on the occasion of his 60th birthday and part of the special issue “Recent Advances in PDE: Theory, Computations and Applications” in his honour.


[1] M. Baldauf, Linear stability analysis of Runge–Kutta-based partial time-splitting schemes for the Euler equations, Monthly Weather Rev. 138 (2010), 4475–4496. 10.1175/2010MWR3355.1Search in Google Scholar

[2] D. Estep, V. Ginting, D. Ropp, J. N. Shadid and S. Tavener, An a posteriori-a priori analysis of multiscale operator splitting, SIAM J. Numer. Anal. 46 (2008), no. 3, 1116–1146. 10.1137/07068237XSearch in Google Scholar

[3] A. Gassmann and H.-J. Herzog, A consistent time-split numerical scheme applied to the nonhydrostatic compressible equations, Monthly Weather Rev. 135 (2007), 20–36. 10.1175/MWR3275.1Search in Google Scholar

[4] E. Hansen and A. Ostermann, Exponential splitting for unbounded operators, Math. Comp. 78 (2009), no. 267, 1485–1496. 10.1090/S0025-5718-09-02213-3Search in Google Scholar

[5] W. Hundsdorfer and J. Verwer, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, Springer Ser. Comput. Math. 33, Springer, Berlin, 2003. 10.1007/978-3-662-09017-6Search in Google Scholar

[6] T. Jahnke and C. Lubich, Error bounds for exponential operator splittings, BIT 40 (2000), no. 4, 735–744. 10.1023/A:1022396519656Search in Google Scholar

[7] S. MacNamara and G. Strang, Operator splitting, Splitting Methods in Communication, Imaging, Science, and Engineering, Sci. Comput., Springer, Cham (2016), 95–114. 10.1007/978-3-319-41589-5_3Search in Google Scholar

[8] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968), 506–517. 10.1137/0705041Search in Google Scholar

Received: 2017-11-17
Revised: 2018-04-16
Accepted: 2018-06-07
Published Online: 2018-07-12
Published in Print: 2019-04-01

© 2019 Walter de Gruyter GmbH, Berlin/Boston

This work is licensed under the Creative Commons Attribution 4.0 Public License.

Downloaded on 9.12.2022 from
Scroll Up Arrow