A Finite Element Splitting Method for a Convection-Diffusion Problem

1 Introduction

In this paper, we shall consider a numerical method for the solution of the convection-diffusion problem in the square Ω = ( 0 , 1 ) × ( 0 , 1 ) ,

with periodic boundary conditions, where, with x = ( x 1 , x 2 ) , the initial function V = V ⁢ ( x ) , the positive definite 2 × 2 matrix a = a ⁢ ( x ) = ( a i ⁢ j ⁢ ( x ) ) , the vector b = b ⁢ ( x ) = ( b 1 ⁢ ( x ) , b 2 ⁢ ( x ) ) and the forcing term F = F ⁢ ( x , t ) are 1-periodic in x 1 and x 2 and smooth. Our method is an explicit-implicit time stepping method based on Lie splitting of a spatially discrete finite element version of (1.1).

With A ⁢ U = - ∇ ⋅ ( a ⁢ ∇ ⁡ U ) and B ⁢ U = b ⋅ ∇ ⁡ U , the exact solution of (1.1) may be formally expressed as

This representation in terms of the solution operator ℰ ⁢ ( t ) of the homogeneous case of (1.1) is the basis for our discretization method. Such a method was analyzed in [1] within the framework of finite differences, and here our purpose is to carry out the corresponding program with finite elements. We refer to [1] also for further references to work on splitting methods.

As a first step to define our finite element splitting method, we thus consider a spatially discrete finite element version of (1.1). Let 𝒯 h be a quasi-uniform family of triangulations of Ω, and let S h be the periodic, continuous piecewise linear functions on 𝒯 h . With ( ⋅ , ⋅ ) the inner product in L 2 ⁢ ( Ω ) , and

the standard Galerkin spatially semidiscrete version of (1.1) is then to find u ⁢ ( t ) ∈ S h for t ≥ 0 such that

This equation may also be expressed in matrix form. With { Φ j } j = 1 N the pyramid function basis for S h , let

With u h ⁢ ( t ) = ∑ j = 1 N 𝐮 j ⁢ ( t ) ⁢ Φ j and 𝐮 = ( 𝐮 1 , … , 𝐮 N ) T , we have

For our purposes, it will be more convenient to use instead a semidiscrete method based on the lumped mass method, employing the approximate inner product on S h defined by

where, for a triangle τ of the triangulation 𝒯 h , P τ , j , j = 1 , 2 , 3 , are its vertices and where | τ | = area ⁡ ( τ ) ; cf. [2, Chapter 15]. We shall then consider the semidiscrete problem, with u ⁢ ( t ) ∈ S h for t ≥ 0 ,

Introducing the operators A h , B h : S h → S h by

equation (1.3) may also be written as

where P ~ h : L 2 → S h is defined by ( P ~ h ⁢ F , χ ) h = ( F , χ ) for all χ ∈ S h . In matrix form, this means that the matrix ℳ in (1.2) is replaced by the diagonal matrix 𝒟 = ( d i ⁢ j ) , where d i ⁢ j = ( Φ i , Φ j ) h . The solution of (1.4) satisfies

We recall that E h ⁢ ( t ) , the solution operator of the homogeneous case of (1.4), is stable and that the solution of (1.4) satisfies an O ⁢ ( h 2 ) error estimate for suitable v.

To define our basic finite element splitting method, let k be a time step and t n = n ⁢ k for n ≥ 0 . On each time interval ( t n - 1 , t n ) , we then introduce the Lie splitting E h , k = e k ⁢ B h ⁢ e - k ⁢ A h of E h ⁢ ( k ) and define a time discrete solution of (1.4) by

We note that the two factors of E h , k are associated with the hyperbolic and parabolic parts of equation (1.1),

For a computable discretization in space and time, we then need to replace e - k ⁢ A h and e k ⁢ B h by rational functions of A h and B h , respectively. For the parabolic part, we shall use the implicit backward Euler operator

where I h is the identity operator on S h . For the hyperbolic part, we define

where the operator L h , defined by ( L h ⁢ ψ , χ ) h = ( ∇ ⁡ ψ , ∇ ⁡ χ ) for all ψ , χ ∈ S h , is added to secure stability, under a mesh-ratio condition for k / h . In matrix form, the application of H h , k takes the form

where the diagonal matrix on the left shows that H h , k is essentially explicit.

More generally, on each interval ( t n - 1 , t n ) , we shall allow the use m steps of the hyperbolic time stepping operator with step k m = k / m , i.e., instead of H h , k , we apply H h , k m m ≈ e k ⁢ B h . This will increase the accuracy and the bound for the mesh ratio but not significantly the cost of the computation since H h , k m is explicit. We consider thus the time discrete solution defined for n ≥ 1 by

Our main result is then that, for v appropriate, the error satisfies

We remark that the first term in this error bound comes from the spatial discretization, the second from the splitting and the backward Euler discretization of the parabolic part, and the third from the approximation of the hyperbolic part. In [1], numerical illustrations of these partial errors were presented in the finite difference case, in the present context essentially corresponding to uniform triangulations, and we note that the spatial discretization is then of order O ⁢ ( h 2 ) . In [1], the choice of m to balance the last two terms was also discussed.

2 Splitting of the Semidiscrete Problem

In this section, we will show that the result u n in (1.5) of Lie splitting of the spatially discrete problem (1.3) differs from the exact solution U ⁢ ( t n ) of (1.1) by O ⁢ ( h + k ) , under the appropriate regularity assumptions and choice of v.

For 1-periodic functions, we define ∥ V ∥ s = ∥ V ∥ H s ⁢ ( Ω ) for s ≥ 0 . Note that ∥ A ⁢ V ∥ s ≤ C ⁢ ∥ V ∥ s + 2 , ∥ B ⁢ v ∥ s ≤ C ⁢ ∥ V ∥ s + 1 and ∥ V ∥ s ≤ C ⁢ ( ∥ A s / 2 ⁢ V ∥ + ∥ V ∥ ) , where ∥ V ∥ = ∥ V ∥ L 2 ⁢ ( Ω ) . Further, for (1.1), we recall the stability estimates

and the smoothing property

We begin the analysis with the following stability result for the operator E h , k in (1.5). Here and below, we shall use the norm ∥ χ ∥ h = ( χ , χ ) h 1 / 2 on S h , which is equivalent to ∥ χ ∥ or, more precisely, satisfies

We next show an error estimate for one step of the Lie splitting of ℰ ⁢ ( k ) .

We shall need error bounds for the spatial discretizations of the parabolic and hyperbolic parts of E h , k . Here and below, we shall use the Ritz projection R h : H 1 → S h defined by A ^ ⁢ ( R h ⁢ V - V , χ ) = 0 for all χ ∈ S h , where A ^ ⁢ ( ψ , χ ) = A ⁢ ( ψ , χ ) + ( ψ , χ ) . Recall that

We also recall that, for ε h ⁢ ( ψ , χ ) = ( ψ , χ ) h - ( ψ , χ ) ,

where, in the last step, we have used the inverse inequality

valid since the family { 𝒯 h } is quasi-uniform.

This shows the following error estimate for the operator E h , k .

We now show a global error estimate for the homogeneous equation.

One possible choice for v with ∥ v - V ∥ ≤ C ⁢ h ⁢ ∥ V ∥ 1 which will be used below is v = P ~ h ⁢ V . In fact, since

we find ∥ P ~ h ⁢ V - P h ⁢ V ∥ h ≤ C ⁢ h ⁢ ∥ P h ⁢ V ∥ 1 ≤ C ⁢ h ⁢ ∥ V ∥ 1 , from which our claim follows.

We now show a complete error estimate for our basic splitting method.

3 Complete Discretization in Time and Space

We now turn to the analysis of the complete discretization using the operators Q h , k and H h , k m m in (1.6) and (1.7) for the parabolic and hyperbolic factors of E h , k = e k ⁢ B h ⁢ e - k ⁢ A h . We note that, since

we have ∥ B h ⁢ ψ ∥ h ≤ β 1 ⁢ ∥ ∇ ⁡ ψ ∥ for ψ ∈ S h . Further, by (2.9), ∥ L h ⁢ ψ ∥ h ≤ ν ⁢ h - 1 ⁢ ∥ ∇ ⁡ ψ ∥ .

We shall first consider the stability of E ~ h , k , m .