First order least squares method with weakly imposed boundary condition for convection dominated diffusion problems
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
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
where () is a polyhedral domain, , a function in , a function in and a function in . Here, the variable flux satisfies the following assumption:
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 ). We notice that in , residual-free bubble strategy is used to address the second drawback. But, the least squares method in  needs to compute basis functions element-wisely, which is relatively not easy to implement.
It is well known that streamline diffusion method , 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 . So, it is necessary to look for an alternative approach. We notice that a weighted test function was used in  to obtain the stability of the original DG method  for the transportation reaction equation, and this idea was generalized to convection-diffusion-reaction equation in  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 . The advantage of this new approach is that the weight function in  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 . 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 ,
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 . However, we do not have to choose any penalty terms empirically while  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 ). 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 .  is the first paper which gives robust a priori error estimate for variable flux . However, IP–DG method  has to choose stabilization parameter empirically while our method does not. In addition, IP–DG method  can have robust convergence only when the mesh size where is a positive constant depending on (see Theorem in ). 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:
We define the finite element space , where
Here, is the space of polynomials on the domain of total degree at most , a non-negative integer. Obviously, there is a positive constant ,
The first order least squares method is to find satisfying
In (2.4), the inner products are defined in the following natural manner:
with the obvious modifications for vector-valued functions.
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 , the first order least squares method in  is to find satisfying
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.
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
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.
In addition, if we further assume and ,
If for any ,
In addition, if we further assume and ,
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
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  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
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
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
Then, the DPG approximation to , lies in a finite dimensional trial subspace . It satisfies
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
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 ,
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 ,
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 . In , the solution space is restricted to Hilbert spaces. We put the proof here in order to make this paper more self-contained.
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
Then, for any , we have that
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
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
Obviously, is a Hilbert space with respect to the inner product in (4.1). For any , we define
The norm of is
The norm of is defined by
4.1. Rewrite least squares method in DPG framework
Equivalently, the first order least squares method (2.4) is to find satisfying
Here, the bi-linear form is defined by
It is easy to see that for any ,
4.2. Numerical stability
If the assumptions (1.2) hold, then there is a positive constant , which is independent of , such that for any ,
In addition, if for any , then we have
Given , we define
Recall that is introduced in (1.2a). Let be a non-negative number, which will be determined below. We take
where for any . It is easy to see that , and
Applying integration by parts on the term , we have
If we take big enough (which depends only on and ),
By (2.3a) and shape regularity assumption, we have
We recall . So, we have
It is easy to see . So, we have
By the definition (4.4) of norm of , we have
If the assumptions (1.2) hold,
Here, the constants and are independent of and .