A First-Order Explicit-Implicit Splitting Method for a Convection-Diffusion Problem

1 Introduction

In this paper, we shall consider the numerical solution by finite differences and splitting of the convection-diffusion equation

in the cube Ω = ( 0 , 2 ⁢ π ) d , with d = 1 , 2 , 3 , under periodic boundary conditions. With x = ( x 1 , … , x d ) , we assume that the positive definite d × d matrix a = a ⁢ ( x ) = ( a i ⁢ j ⁢ ( x ) ) and the vector b = b ⁢ ( x ) = ( b 1 ⁢ ( x ) , … , b d ⁢ ( x ) ) as well as the forcing term F = F ⁢ ( x , t ) and the initial data V = V ⁢ ( x ) are periodic and smooth.

Equation (1.1) is a special case of the initial-value problem for the operator equation

where 𝒜 and ℬ represent different physical processes, in our case diffusion and convection, respectively, or

The solution of (1.2) may be formally expressed as

To discretize (1.2) in time, let k be a time step, and set t n = n ⁢ k . We then have

We may then find an approximate solution U ˘ n ≈ U ⁢ ( t n ) , n ≥ 1 , by approximating the integrand by its value at t n , or U ˘ n = ℰ ⁢ ( k ) ⁢ U ˘ n - 1 + k ⁢ F n , where F n = F ⁢ ( t n ) for n ≥ 1 , with U ˘ 0 = V .

In the next step, we shall approximate ℰ ⁢ ( k ) = e - k ⁢ ( 𝒜 - ℬ ) by splitting 𝒜 - ℬ into 𝒜 and - ℬ , and replace ℰ ⁢ ( k ) by using the Lie splitting

so that the time discrete solution U n , n ≥ 1 , is defined by

Applying ℰ k thus involves solving the parabolic equation U t = - 𝒜 ⁢ U and the hyperbolic equation U t = ℬ ⁢ U on ( t n - 1 , t n ) ; see e.g. Hundsdorfer and Verwer [5], Descombes [1] and references therein. The two exponential factors of ℰ k may then be further approximated by rational functions of 𝒜 and ℬ , respectively, such as the backward or forward Euler approximations.

We note that if 𝒜 and ℬ commute, which holds for (1.1) when a and b are independent of x, then e - k ⁢ ( 𝒜 - ℬ ) = e k ⁢ ℬ ⁢ e - k ⁢ 𝒜 so that the error in (1.4) is zero. When 𝒜 and ℬ do not commute, then formally, by Taylor expansion, e k ⁢ ℬ ⁢ e - k ⁢ 𝒜 - e - k ⁢ ( 𝒜 - ℬ ) = O ⁢ ( k 2 ) . Under the appropriate assumptions, this will lead to an O ⁢ ( k ) error bound in (1.5) for t n bounded.

Another possible approximation of ℰ ⁢ ( k ) is the Strang splitting ℰ ⁢ ( k ) ≈ ℰ k = e 1 2 ⁢ k ⁢ ℬ ⁢ e - k ⁢ 𝒜 ⁢ e 1 2 ⁢ k ⁢ ℬ for which the symmetry gives local accuracy of order O ⁢ ( k 3 ) , resulting in an error estimate in (1.5) of order O ⁢ ( k 2 ) for t n bounded, cf. [8]. In [9], a finite difference method based on the Strang splitting was analyzed for the homogeneous case of (1.1). The present paper using the less accurate Lie splitting may be thought of as a complement to [9]. It is intended to provide background for future work on more general problems, considering for instance maximum-norm estimates, reaction-diffusion equations, problems with nonsmooth data and also the combination of splitting with the finite element method. As a first step, the case of nonhomogeneous equations is included here.

Error estimates for time splittings may be found in [6, 3, 4, 2] and references therein. For an extensive list of related work, see MacNamara and Strang [7].

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

We denote by 𝒫 h the M-periodic vectors u with elements u ω , corresponding to the mesh points x ω , and with u ω + M ⁢ e l = u ω , l = 1 , … , d , where e l is the unit vector in the x l direction. Note that we use capital letters for functions in Ω and lowercase letters for vectors in 𝒫 h . We consider now the second-order spatially discrete, continuous in time, finite difference ODE system for u ⁢ ( t ) ∈ 𝒫 h ,

where the M d × M d matrices A and B corresponding to the differential operators 𝒜 and ℬ in (1.3) are defined by

Here ∂ j and ∂ ¯ j are forward and backward finite difference quotients in the direction of x j ,

and v and f the restrictions of V and F to Ω h . It is then to (1.7) that we will apply the Lie splitting.

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

and we shall see that the error in this approximation is of order O ⁢ ( h 2 ) , under the appropriate regularity assumptions.

For the time discretization of (1.7), an obvious first choice would be the backward Euler method, defining the discrete solution u ˘ n at t n by

The basis of the method we study here, however, will be the splitting method analogous to (1.4), (1.5), defined by

Similarly to the continuous problem, when A and B commute, E k = E ⁢ ( k ) .

For a practical numerical method, we then need to approximate the two exponential factors in E k by rational functions. For the approximation of e - k ⁢ A , we shall use the unconditionally stable backward Euler operator Q k = ( I + k ⁢ A ) - 1 ≈ e - k ⁢ A for k > 0 . To approximate e k ⁢ B , we would like to use the forward Euler method, or e k ⁢ B ≈ I + k ⁢ B . For stability reasons, we shall need to add an artificial viscosity term, so we set L ⁢ u = - ∑ j = 1 d ∂ ¯ j ⁢ ∂ j ⁡ u and define H k = I + k ⁢ B - γ ⁢ k 2 ⁢ L for k h ≤ ρ 0 , where the positive constants γ and ρ 0 will be defined later.

In order to increase the limit ρ 0 for the mesh ratio, and to increase the accuracy of the approximation in the hyperbolic part, one may replace the operator H k by H k m m , with k m = k m for some positive integer m, and the stability requirement will then be reduced to k h ≤ m ⁢ ρ 0 . This means that the approximate solution of the hyperbolic part of E k is obtained in m steps of length k m = k m .

We consider thus the time discrete solution at time t n = n ⁢ k , defined by

This method, which we will refer to as our fully discrete method, thus replaces at each time step the backward Euler solution of our nonsymmetric problem (1.7) by a backward Euler approximation of a symmetric parabolic problem, followed by an explicit finite difference solution of a hyperbolic problem.

After the introduction of notation and some preliminary observations in Section 2, the spatially semidiscrete problem (1.7), (1.8), the backward Euler method (1.9) and the corresponding basic time splitting method (1.10) will be analyzed in Section 3, where we will show O ⁢ ( h 2 ) and O ⁢ ( h 2 + k ) error bounds for these methods, respectively. The analysis of the fully discrete method (1.11) will be carried out in Section 4, using discrete Sobolev norms. The time discretization error can be divided into an O ⁢ ( k ) part associated with the parabolic part of the equation and an O ⁢ ( k m ) term for the hyperbolic part. In Section 5, we illustrate our theoretical results with computations for examples with one and two spatial dimensions.

2 Notation and Preliminaries

With Ω h as in (1.6), we introduce the discrete inner product and norm

Further, for x ∈ Ω h , we set

and define the discrete Sobolev inner product and norm, with | α | = α 1 + … + α d , by

We shall also use

For a domain 𝒢 ⊂ R d , we define the norm on the Sobolev space H s ⁢ ( 𝒢 ) by

For U periodic, we write ∥ U ∥ s for its norm on H s = H s ⁢ ( Ω ) . We note that, defining U h to be the restriction of U to the mesh Ω h , i.e., by ( U h ) ω = U ⁢ ( x ω ) , we have

In fact, let

Then, for | α | = s , ℓ ⁢ ( U ) = ( ∂ α ⁡ U h ) ω is a bounded linear functional on ℂ ⁢ ( Q s ⁢ ( x ω ) ) , and therefore, by the Sobolev embedding theorem, if s > d 2 , also on H s ⁢ ( Q s ⁢ ( x ω ) ) , which vanishes for polynomials of degree ≤ s - 1 . A simple argument using the Bramble–Hilbert lemma therefore shows

Hence | U h | h , s ≤ C ⁢ ∥ U ∥ s , where | v | h , s 2 = ∑ | α | = s ∥ ∂ α ⁡ v ∥ h 2 . Since ∥ U h ∥ h , s ≤ C ⁢ ( ∥ U h ∥ h + | U h | h , s ) , this shows (2.1). In particular, ∥ U h ∥ h , 2 ≤ C ⁢ ∥ U ∥ 2 for d = 1 , 2 , 3 ; it is our need for this inequality which is the reason for our restriction on d in (1.1).

Consider now the matrices A and B in (1.8), and note that, for | α | ≤ s ,

In fact, ∂ α ⁡ A and A ⁢ ∂ α are linear combinations of difference quotients of orders ≤ s + 2 , and in the difference, the highest-order terms cancel so that the orders are ≤ s + 1 , which shows the first inequality. The second inequality is shown analogously. Further, with C independent of h,

The matrix A is positive semidefinite, with

Since the terms in A ⁢ U h are symmetric difference quotients of U at the mesh points x ω , and the terms in ( 𝒜 ⁢ U ) h are the corresponding derivatives of U at x ω , we find, similarly to (2.2),

and similarly for B,

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

For B = ( b ω , χ ) , we have

and that then b ω + e j , ω = - 1 2 ⁢ h ⁢ b j ⁢ ( x ω + e j ) . We may write

Here B 0 is skew-symmetric and B 1 symmetric, with

In fact, for d = 1 , with l = ω 1 , we have b l , l ± 1 = ± 1 2 ⁢ b ⁢ ( x l ) h , and hence

so that ∥ B 1 ⁢ v ∥ h ≤ 1 2 ⁢ ∥ b ′ ∥ ℂ ⁢ ∥ v ∥ h . The case d > 1 is treated analogously.

We note that ( B 0 ⁢ v , v ) = 0 for all v so that

We shall also use

In fact, since ( B ⁢ v ) ω = 1 2 ⁢ ( ∑ j = 1 d b j ⁢ ( ∂ j + ∂ ¯ j ) ⁢ v ) ω , we find

Thus ∥ B ⁢ v ∥ h 2 ≤ β ~ ⁢ ∑ j = 1 d ∥ ∂ j ⁡ v ∥ h 2 = β ~ ⁢ ( L ⁢ v , v ) h .

3 The Semidiscrete Problem, the Backward Euler Method and the Basic Splitting

We begin with the straightforward standard error analysis of the spatially semidiscrete problem (1.7), which we include for completeness. We first show stability and a smoothing property of the solution operator of (1.7) for f = 0 , in discrete Sobolev norms.

Note that the special case of E ⁢ ( t ) = e - t ⁢ A is included for B = 0 .

As a consequence of Lemma 3.1, we have the following second-order error estimate for the spatially semidiscrete solution u ⁢ ( t ) of (1.7).

We now turn to the backward Euler method (1.9), the analysis of which will be the basis for the analysis of our splitting method (1.11). We begin with stability.

We can now prove the following error estimate for the time discretization of the homogeneous equation.