First order least squares method for convection diffusion problems

First order least squares method with weakly imposed boundary condition for convection dominated diffusion problems

Huangxin Chen School of Mathematical Sciences, Xiamen University, Xiamen, 361005, People’s Republic of China chx@xmu.edu.cn Guosheng Fu School of Mathematics, University of Minnesota, Minneapolis, MN 55455, USA fuxxx165@math.umn.edu Jingzhi Li Department of Financial Mathematics and Financial Engineering, South University of Science and Technology of China, Shenzhen, 518055, People’s Republic of China li.jz@sustc.edu.cn  and  Weifeng Qiu Department of Mathematics, City University of Hong Kong, 83 Tat Chee Avenue, Kowloon, Hong Kong, China weifeqiu@cityu.edu.hk
Abstract.

We present and analyze a first order least squares method for convection dominated diffusion problems, which provides robust a priori error estimate for the scalar variable even if the given data . The novel theoretical approach is to rewrite the method in the framework of discontinuous Petrov–-Galerkin (DPG) method, and then show numerical stability by using a key equation discovered by J. Gopalakrishnan and W. Qiu [Math. Comp. 83(2014), pp. 537-552]. This new approach gives an alternative way to do numerical analysis for least squares methods for a large class of differential equations. We also show that the condition number of the global matrix is independent of the diffusion coefficient. A key feature of the method is that there is no stabilization parameter chosen empirically. In addition, Dirichlet boundary condition is weakly imposed. Numerical experiments verify our theoretical results and, in particular, show our way of weakly imposing Dirichlet boundary condition is essential to the design of least squares methods - numerical solutions on subdomains away from interior layers or boundary layers have remarkable accuracy even on coarse meshes, which are unstructured quasi-uniform.

Key words and phrases:
error estimate; least squares method; DPG method; convection diffusion problems.
2000 Mathematics Subject Classification:
65N30, 65L12, 35L15
Corresponding author: Weifeng Qiu

1. Introduction

In this paper, we present a robust a priori analysis of first order least squares method with weakly imposed boundary condition for the following convection dominated diffusion equation

(1.1a)
(1.1b)

where () is a polyhedral domain, , a function in , a function in and a function in . Here, the variable flux satisfies the following assumption:

(1.2a)
(1.2b)

According to [1], the assumption (1.2a) is satisfied if

Least squares methods have been frequently used to simulate solutions of partial differential equations arising from fluid dynamics and continuum mechanics. We refer to [7, 43] for comprehensive summary. It is well known that least squares methods have the following desirable features: it leads to a minimization problem; its numerical stability is not sensitive to the choice of finite element space or meshes; the resulting global stiffness matrices are symmetric and positive definite; a practical a posteriori error estimator can be given without any additional cost, and so on (see [5, 6, 14, 15, 17, 16, 18, 21, 22, 23, 26, 37, 48]).

Unfortunately, primitive least squares methods for convection dominated diffusion problems (1.1) have the following drawbacks. Firstly, if the term is not uniformly bounded from below by a positive constant, a priori error estimate of primitive least squares methods will deteriorate as the diffusion coefficient goes to zero, even when the exact solution has no interior layers or boundary layers (see error estimates in [19, 41, 45]). Secondly, primitive least squares methods show a very poor performance for convection diffusion problem (1.1) with a sufficiently small diffusion coefficient, because large spurious oscillations are observed (see numerical experiments in [41]). We notice that in [41], residual-free bubble strategy is used to address the second drawback. But, the least squares method in [41] needs to compute basis functions element-wisely, which is relatively not easy to implement.

It is well known that streamline diffusion method [10], residual free bubble methods [8, 9, 13], and DG methods [1, 40, 42, 32, 33] do not suffer from the above two drawbacks of primitive least squares methods. We refer to [50, 52] as comprehensive summaries of numerical methods suitable for convection dominated diffusion problems. We would like to emphasize that none of these numerical methods (streamline diffusion method, residual free bubble methods, DG methods) results in symmetric global stiffness matrices. Hence from the point of view of solver design, the least squares method is more attractive than the other methods mentioned before and many works have been contributed to this subject (e.g. [16, 46]). Moreover, we derive that the condition number of linear system from our first order least squares method is at most , where is the mesh size. In particular, the condition number is independent of the diffusion coefficient. This property is important for designing efficient solver, e.g., multilevel method, for the first order least squares approximation of convection dominated diffusion equation.

In this paper, we propose and analyze our first order least squares method to address these two drawbacks for primitive least squares methods. In fact, it is difficult to provide robust error estimate by the traditional approach of numerical analysis for least squares methods in [7]. So, it is necessary to look for an alternative approach. We notice that a weighted test function was used in [44] to obtain the stability of the original DG method [49] for the transportation reaction equation, and this idea was generalized to convection-diffusion-reaction equation in [1] using the IP–DG methods. In this paper, we rewrite our method in the framework of discontinuous Petrov–-Galerkin (DPG) method, then show numerical stability by using a key equation discovered in [39]. The advantage of this new approach is that the weight function in [1] is shown to stay in some “equivalent” test function space (see (3.7) in section 3) such that numerical stability can be obtained without using any projection as in [1]. This approach is novel and useful to numerical analysis of least squares methods for a large class of differential equations. This new approach of numerical analysis is also different from traditional ones used for DPG method in [11, 27, 28, 29, 30]. We show that, roughly speaking, using polynomials of degree ,

(1.3)

if for any ,

Here, the constant is independent of . Thus, we can conclude that a priori error estimate in (1.3) is robust with respect to the diffusion coefficient , which addresses the first draw back. We also want to emphasize that the convergence result (1.3) shows our method has convergence rate even if , which means our method does not have excessive smoothness requirements than other methods. In order to overcome the second drawback, we impose Dirichlet boundary condition in an weak way, such that the error along the boundary layers will not propagate into the whole domain. We show the advantage of imposing boundary condition weakly by numerical experiments. We notice that our way of imposing boundary condition is similar to the weak imposition of Dirichlet boundary condition in [12, 2, 51], which belongs to Nitsche’s method in [47]. However, we do not have to choose any penalty terms empirically while [2] needs. We would like to emphasize that weakly imposing boundary condition is essential to least squares methods while it is incrementally helpful to streamline diffusion method and DG methods (see numerical experiments in [2]). If boundary condition is imposed strongly, the numerical solutions produced by streamline diffusion method and DG methods may have artificial oscillation along boundary layers, while the accuracy in subdomains away from boundary layers is still remarkable. However, according to our numerical experiments, if we impose boundary condition strongly, then numerical solution of least squares methods will be polluted on almost the whole domain by boundary layers. We have tried to add several stabilization terms, which have been utilized by streamline diffusion method or DG methods, into least squares methods to prevent propagation of error from boundary layers. None of them works except weakly imposing boundary condition.

This paper is the first one which addresses the two drawbacks of least squares methods for convection dominated diffusion problems. Now, we would like to compare with DPG methods, which can be considered as a special class of least squares methods. We notice that all DPG methods [3, 20, 31] need to compute optimal test function space. on each element such that the implementation is more complicated than ours (Methods in [24, 25] are similar to DPG methods.) Next, we would like to compare our results with IP–DG method [1]. [1] is the first paper which gives robust a priori error estimate for variable flux . However, IP–DG method [1] has to choose stabilization parameter empirically while our method does not. In addition, IP–DG method [1] can have robust convergence only when the mesh size where is a positive constant depending on (see Theorem in [1]). On the contrast, the convergence result of our method does not have this restriction.

The remainder of this paper is organized as follows. In section , we introduce our first order least squares method and the main theoretical results. In section , we show a novel approach to do numerical analysis for least squares methods (not restricted to our first order least squares method for convection dominated diffusion problems). In section , we prove a priori error estimates by using the approach introduced in section . In section , we prove the estimate of condition number of global stiffness matrix provided by our method. In section , we give numerical experiments which verify our theoretical results. In section , we extend our first order least squares method for transportation reaction problems.

2. First order least squares method and main theoretical result

In this section, we present setting of meshes, first order least squares method and the main theoretical results.

2.1. Setting of meshes

Let be a conforming triangulation of the domain made of shape-regular simplexes . For each element , we set and for each of its faces , , where denotes the Lebesgue measure in or dimensions. We associate to the set of faces as well as those of interior faces and boundary faces . We say that if there are two elements and in such that , and we say that if there is an element in such that . It is obvious that .

2.2. First order least squares method

We define . We can rewrite (1.1) as the following first-order equations:

(2.1a)
(2.1b)
(2.1c)

We define the finite element space , where

(2.2a)
(2.2b)

Here, is the space of polynomials on the domain of total degree at most , a non-negative integer. Obviously, there is a positive constant ,

(2.3a)
(2.3b)

The first order least squares method is to find satisfying

(2.4)

In (2.4), the inner products are defined in the following natural manner:

with the obvious modifications for vector-valued functions.

Remark 2.1.

Notice that the boundary condition imposed on outflow boundary () is . This is the reason we say Dirichlet boundary condition is weakly imposed. If we ignore the weight function used in [19], the first order least squares method in [19] is to find satisfying

(2.5)

where and . The main difference between (2.4) and (2.5) is that in (2.5), Dirichlet boundary condition is imposed strongly while our method (2.4) uses weakly imposed boundary condition. Our convergence analysis in section 4 is also valid for least squares method (2.5). But, if the solution of (1.1) has boundary layers or interior layers, then our way of weakly imposing boundary condition in (2.4) can improve the accuracy of numerical solution in subdomains away from layers dramatically. We refer to section for detailed description.

Remark 2.2.

In (2.4), we weakly imposed the Dirichlet boundary condition by

In fact, we can also weakly imposed the Dirichlet boundary condition as

Then, the first order least squares method is to find satisfying

(2.6)

We can still obtain the stability estimate similar to Lemma 4.1. However, the convergence rate of (2.6) is the same as that of (2.4), since the approximation to in -norm gets involved with the error analysis (see the proof of Theorem 2.3).

2.3. Main theoretical results

We state our main theoretical results on a priori error analysis and condition number estimation of the first order least squares method as the following Theorems.

We denote by a positive constant, which is independent of and . We assume the assumptions (1.2) on and hold.

Theorem 2.3.
(2.7)

In addition, if we further assume and ,

(2.8)
Theorem 2.4.

If for any ,

(2.9)

In addition, if we further assume and ,

(2.10)
Theorem 2.5.

We denote by the condition number of global stiffness matrix of first order least squares method (2.4). If we assume the meshes are quasi-uniform, then

(2.11)
Remark 2.6.

In Theorem 2.4, we obtain optimal convergence of if for any . The restriction on mesh size is due to the energy norm of in Lemma 4.1 is

Thus, we need to utilize inverse inequality and the above restriction on mesh size to obtain upper bound on .

3. The approach to analysis

In this section, an alternative approach is introduced for numerical analysis of least squares methods in general (not restricted to the first order least squares method (2.4)). First of all, we show least squares methods and DPG method in abstract settings, respectively. Then, we rewrite least squares methods in the framework of DPG method. A key equation in [39] is then used to show how to achieve numerical stability of least squares methods.

Suppose we want to approximate solution satisfying

Here, is a Banach space with norm , and is a Hilbert space with norm . We assume that the linear operator is in , which is the space of bounded linear operators from to . Then, least squares methods are to find satisfying

(3.1)

Here, is a finite dimensional trial subspace of (where denotes a parameter determining the finite dimension).

On the other hand, DPG method can be described in the following general context. Suppose we want to approximate satisfying

(3.2)

Here is a Banach space with norm and is a Hilbert space under an inner product with corresponding norm . We assume that the bi-linear form is continuous and the linear form is also continuous. Define by

(3.3)

Then, the DPG approximation to , lies in a finite dimensional trial subspace . It satisfies

(3.4)

where . Since in general, this is a Petrov-Galerkin approximation. The method (3.4) is the DPG method. The excellent stability and approximation properties of this method are well known [28, 29].

Now, we can rewrite the least squares methods (3.1) in the framework of DPG method in the following way. The corresponding DPG method is to find such that

(3.5)

where the bi-linear form

Here, for any ,

We say the corresponding DPG method (3.5) is numerical stable if there is a constant independent of mesh size , such that for any ,

(3.6)

In general, it is not easy to know what are elements in the finite dimensional space . So, it is usually not easy to estimate the right hand side of the inequality (3.6). But, according to the following Theorem 3.1, we have that for any ,

(3.7)

(3.7) is the key equation discovered in work for DPG methods. We call the “equivalent” test function space. By (3.7), in order to obtain (3.6), it is sufficient to achieve

(3.8)

The inequality (3.8) is usually easier to be obtained than (3.6). Notice that the least squares methods (3.1) is actually the same as the corresponding DPG method (3.5). This explains why it is relatively easier to obtain numerical stability of least squares methods. In fact, first order least squares methods based on first order Friedrichs’ systems in [34, 35, 36] can be shown to have numerical stability by using the key equation (3.7).

The proof of following Theorem 3.1 is included in that of Theorem  in [39]. In [39], the solution space is restricted to Hilbert spaces. We put the proof here in order to make this paper more self-contained.

Theorem 3.1.

Let be a continuous bi-linear mapping. Here, is a normed linear space, and is a Hilbert space with associate norm .

We define a linear operator by

(3.9)

Then, for any , we have that

(3.10)
Proof.

Since the bi-linear mapping is continuous and is a Hilbert space, the operator is well-defined.

We take arbitrarily. It is straightforward to see that

since .

If , then (3.10) is obviously true due to (3.9). If , then by (3.9) and the fact that is a Hilbert space, we have

We can conclude that for any , (3.10) is true. ∎

4. Convergence analysis

In this section, we shall prove Theorem 2.3 and Theorem 2.4. According to the approach of analysis in section 3, we rewrite first order least squares method (2.4) in the framework of DPG method. Then, we show the numerical stability of corresponding DPG method. Finally, we prove Theorem 2.3 and Theorem 2.4.

Throughout this section, we define and . The inner product of is defined by

(4.1)

Obviously, is a Hilbert space with respect to the inner product in (4.1). For any , we define

(4.2)

The norm of is

(4.3)

The norm of is defined by

(4.4)

4.1. Rewrite least squares method in DPG framework

The first order least squares method (2.4) is the same as the corresponding DPG method (4.5) in the following.

Equivalently, the first order least squares method (2.4) is to find satisfying

(4.5)

Here, the bi-linear form is defined by

(4.6)

for any , . Obviously, the bi-linear form (4.6) is continuous on with respect to the norm (4.4) and norm (4.3). The operator is defined by

(4.7)

It is easy to see that for any ,

(4.8)

4.2. Numerical stability

In the following Lemma 4.1, we show numerical stability of the corresponding DPG method (4.5). Notice that DPG method (4.5) is the same as first order least squares method (2.4).

Lemma 4.1.

If the assumptions (1.2) hold, then there is a positive constant , which is independent of , such that for any ,

(4.9)

In addition, if for any , then we have

(4.10)
Proof.

According to Theorem 3.1 and (4.7), we have

(4.11)

Given , we define

(4.12)

Recall that is introduced in (1.2a). Let be a non-negative number, which will be determined below. We take

(4.13)

where for any . It is easy to see that , and

Applying integration by parts on the term , we have

Applying integration by parts on the term and assumption (1.2a) and assumption (1.2b), we have

If we take big enough (which depends only on and ),

(4.14)

According to (4.2) and the definition of in (4.13), we have

By (2.3a) and shape regularity assumption, we have

We recall . So, we have

(4.15)

According to (4.14, 4.15), we have

It is easy to see . So, we have

By the definition (4.4) of norm of , we have

Then, by (4.11), we have (4.9) immediately.

By (2.3b), if for any , we have

Combing the above inequality with (4.9), we can conclude that (4.10) holds if for any . ∎

Remark 4.2.

We want to emphasize that the proof of Lemma 4.1 holds for any subspace of , which satisfies (2.3a) and (2.3b). It implies that Lemma 4.1 is still true if we impose Dirichlet boundary condition strongly, like in [19].

Lemma 4.3.

If the assumptions (1.2) hold,

(4.16)

Here, the constants and are independent of and .

Proof.

According to Lemma 4.1 and explicit formulation of operator in (4.8),