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.
In this paper we shall consider the numerical solution of the following convection-diffusion problem in the cube :
under periodic boundary conditions, where the positive definite matrix and the vector 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 , with representing a slow and a fast physical process.
The solution of (1.2) may be formally expressed as
To discretize such an equation in time, a common approach is to split into and . With k a time step one introduces and one may then use the second-order symmetric Strang splitting [8, 7] on each time interval ,
which thus locally involves solutions of and , see e.g. Hundsdorfer and Verwer  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 , with , 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 , 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  discuss the difficulties associated with splitting in such situations. In the case of reaction-diffusion equation, see Estep, Ginting, Ropp, Shadid and Tavener  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 so that the error in (1.3) is zero. When and do not commute, then formally, by Taylor expansion, and . Error estimates for the time splitting, depending on the regularity of the initial values may be found in Jahnke and Lubich , Hansen and Ostermann  and references therein.
In the method that we study in this paper, we begin by discretizing (1.1) in the spatial variables. We let , where M is a positive integer, and define a corresponding uniform mesh
For M-periodic vectors u with elements , corresponding to the mesh-point , and with , we consider the following simple second-order finite difference approximation of (1.1):
Here , are forward and backward finite difference quotients in the direction of , , and v the restriction of V to . Problem (1.5) may be written as a system of ODEs in time,
where the 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 the mesh-width and ,
where (recall ), and
The solution of (1.6), the spatially semidiscrete solution, is
and we shall see that the error in this approximation is , under the appropriate regularity assumptions. For the time discretization we shall work with basic time intervals of length , where 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 . For , 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 in (1.7) we shall use the Crank–Nicolson method. Then, in order to approximate , we would like to use the forward Euler method on a time interval of length . 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 vector inner product. Stability therefore holds if , or if . Since k should be of the same order as h, this makes it natural to choose . We thus subdivide the time intervals of length k into N subintervals of length 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 . We consider thus the time discrete solution at time ,
In fact, in the successive time stepping, only the matrix 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 times before and after the diffusion approximation, thus covering the interval of length . 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 with γ sufficiently small, then the approximation of 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 and with as in (1.4), we introduce the discrete inner product and norm
Further, we set and for and define the discrete Sobolev norm, with
We shall also use . We define
We shall write for , and for . We note that defining to be the restriction to the mesh of a smooth function U, i.e. by , we have .
Consider now the matrices A and B in our convection-diffusion problem (1.2). They satisfy, with C independent of h,
Setting we have, since the terms in are symmetric difference quotients of U at the mesh-points , and the terms in are the corresponding derivatives of U at ,
We note that, for ,
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
For we have
and that then . It follows from this that
Here is skew-symmetric and
Note also that 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.
Let . Then, for any with independent of h,
Let and let . From (1.6) we find, for ,
Here, by (2.4),
Further, since with skew-symmetric,
Therefore, by (2.9),
Hence, using (2.6) and summing over we find, with ,
or, with ,
from which the lemma follows. ∎
Note that the special case of is included for .
As a consequence, we have the following second-order error estimate.
Setting , we find
where . Hence, by Lemma 2.1, since ,
Turning to the analysis of the time discretization, we first show the stability of .
For any , we have, with independent of h,
Here we may choose , as in (2.8).
Let and let . Then for the solution of (1.6) with ,
We now show the following error estimate for the Strang splitting.
We have, with C independent of h and k,
Setting and noting that , we may use Taylor’s formula to obtain
which shows the lemma. ∎
We now turn to the time stepping operator defined in (1.8) and begin with the following stability result.
Let . Then, with as in (2.8),
Further, is stable, and
Since is skew-symmetric, we have
We start the analysis of the time discretization error with the following.
Let M be a square matrix, and assume for . Then we have, for ,
If also for , then
By Taylor expansion we have
where , to complete the proof. ∎
We now show the following error estimate for .
Let . Then we have, with ,
Since , we may write
By Lemma 2.5, we have
which completes the proof. ∎
We now show the following error estimate for .
Let . Then we have, with ,
In view of Lemma 2.3, it remains to show
Here, by the above lemmas,
which completes the proof. ∎
We can now prove the following error estimate:
3 The Case of a Small Diffusion Coefficient
In this section, we consider the variant of problem (1.1) with a small diffusion coefficient ,
The corresponding semidiscrete system (1.6) may then be written
where A and B are as before. We shall see that (3.1) is stable, and satisfies an error estimate, independently of ε. Further, for ε and k small, or more precisely, if and , we will be able to show an 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 ,
We note the inverse inequality , and hence
As in Section 2, we first attend to the spatially semidiscrete problem.
Let . Then, for any , we have, with independent of and ,
Following the steps in the proof of Lemma 2.1, we have, for and ,
and thus, by (2.6), with ,
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.
In the analysis of the time discretization we begin with the analogue of Lemma 2.3.
We have, with independent of h and k,
With , we have and hence, by Taylor’s formula,
Further, and hence
Together these estimates complete the proof of the lemma. ∎
We now turn to the time stepping operator and begin with its stability.
If , then
If also then is stable, or, with μ as in Lemma 2.4,
We note that, since A is positive semidefinite,
We now turn to the error analysis and show the following.
If and , we have
In view of Lemma 3.2 it remains to show
We first note that by Lemma 2.5,
which completes the proof ∎
The following is the resulting error estimate.
If and we have, with independent of and ε,
Using the analogue of (2.18), we find
With and for and we have, with independent of and ε,
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),
As before we shall choose , , where M and N are positive integers, and study the effect of doubling these integers.
We begin with the simple case , in which case the exact solution is . In Table 1 we compile the errors in the numerical solution at , first in the spatial discretization, then when our time stepping method is applied to the semidiscrete solution, and finally the total error. We use , and , so that . 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.
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 and 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 . 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.
We finally consider a numerical example for Section 3, for which we use (4.1) with , . Here and with . Note that the condition , with , now reduces to , or for , which is satisfied for our choice of ε. The results are given in Table 3 and agree with the error bounds of Section 3.
 M. Baldauf, Linear stability analysis of Runge–Kutta-based partial time-splitting schemes for the Euler equations, Monthly Weather Rev. 138 (2010), 4475–4496. Search in Google Scholar
 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. Search in Google Scholar
 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. Search in Google Scholar
 E. Hansen and A. Ostermann, Exponential splitting for unbounded operators, Math. Comp. 78 (2009), no. 267, 1485–1496. Search in Google Scholar
 W. Hundsdorfer and J. Verwer, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, Springer Ser. Comput. Math. 33, Springer, Berlin, 2003. Search in Google Scholar
 T. Jahnke and C. Lubich, Error bounds for exponential operator splittings, BIT 40 (2000), no. 4, 735–744. Search in Google Scholar
 S. MacNamara and G. Strang, Operator splitting, Splitting Methods in Communication, Imaging, Science, and Engineering, Sci. Comput., Springer, Cham (2016), 95–114. Search in Google Scholar
 G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968), 506–517. Search in Google Scholar
© 2019 Walter de Gruyter GmbH, Berlin/Boston
This work is licensed under the Creative Commons Attribution 4.0 Public License.