Adaptive First-Order System Least-Squares Finite Element Methods for Second Order Elliptic Equations in Non-Divergence Form
This paper studies adaptive first-order least-squares finite element methods for second-order elliptic partial differential equations in non-divergence form. Unlike the classical finite element method which uses weak formulations of PDEs not applicable for the non-divergence equation, the first-order least-squares formulations naturally have stable weak forms without using integration by parts, allow simple finite element approximation spaces, and have build-in a posteriori error estimators for adaptive mesh refinements.
The non-divergence equation is first written as a system of first-order equations by introducing the gradient as a new variable. Then two versions of least-squares finite element methods using simple finite elements are developed in the paper, one is the -LSFEM which uses linear elements, the other is the weighted-LSFEM with a mesh-dependent weight to ensure the optimal convergence. Under a very mild assumption that the PDE has a unique solution, optimal a priori and a posteriori error estimates are proved. With an extra assumption on the operator regularity which is weaker than traditionally assumed, convergences in standard norms for the weighted-LSFEM are also discussed. -error estimates are derived for both formulations. We perform extensive numerical experiments for smooth, non-smooth, and even degenerate coefficients on smooth and singular solutions to test the accuracy and efficiency of the proposed methods.
Key words. non-divergence elliptic equation, first-order system least-squares, finite element method, a priori error estimate, a posteriori error estimate, adaptive method.
In this paper we consider finite element approximations of the following elliptic PDE in non-divergence form:
Here, the domain is an open and bounded polytope for or , the coefficient matrix is a positive definite matrix with eigenvalues bounded by below and above on , but not necessarily differentiable. The righthand side is assumed in .
Note that when , we have the following equation in divergence form, where is taken row-wise:
The elliptic PDE in non-divergence form arises in the linearization of fully nonlinear PDEs, for example, stochastic control problems, nonlinear elasticity, and mathematical finance. The matrix is not smooth nor even continuous in many such cases. For example, for fully nonlinear PDEs solving by finite element methods, the coefficient matrix of its linearization is possibly only element-wisely smooth if the coefficients containing derivatives of the numerical solution.
Since the matrix is not differentiable, the standard weak notion of elliptic equation is not applicable. The existence and uniqueness of equations of non-divergence form are often based on the classical or strong senses of the solutions, see discussions in [30, 29]. These PDE theories often assume that the domain is convex, the boundary is sufficiently smooth, or some other restrictive conditions on the smoothness of . For a discontinuous , there are possibilities that the solution is non-unique, see an example given in . It is worth to mention that these available theoretical PDE results are all sufficient theories. For example, since the Poisson equation is also an example of the non-divergence equation with , the existence and uniqueness condition of the equation (LABEL:eq_nondiv) dependence on the domain can be very weak.
There are several numerical methods are available for the problem in non-divergence form. Based on discrete Calderon-Zygmond estimates, Feng, Neilan, and co-authors developed finite element methods for problems with a continuous coefficient matrix in [18, 19, 28]. For equations with discontinuous coefficients satisfying the Cordes condition, a discontinuous Galerkin method , a mixed method , and a non-symmetric method  are developed. A weak Galerkin method is developed by Wang and Wang in . The analysis of these papers mostly assumes the full regularity of the operator and studies the -error estimates of the approximations. In some sense, these methods keep the non-divergence operator second order and borrow techniques from variational fourth-order problems. Nochetto and Zhang  studied a two-scale method, which is based on the integro-differential approach and focuses on error estimates.
Traditionally, the finite element method is based on the variational formulation of an elliptic equation, where the integration by parts plays an essential role. The integration by parts can shift a derivative from the trial variable to the test variable, thus reduces the differential order of the operator. For (LABEL:eq_nondiv), the integration by parts is not available. Luckily, there is another natural method to reduce the differential order of a PDE operator by introducing another auxiliary variable. We can reduce the second order equation into a system of a first order equation by using the new auxiliary variable. Normally, for the first-order system, we can two approaches. One is the mixed method which also involving the integration by parts and has difficulties to ensure the stability. The other method is the least-squares finite element method (LSFEM). The first-order system least-squares principle first re-write the PDE into a first-order system, then define an artificial, externally defined energy-type principle. The energy functional can be defined as summations of weighted residuals of the system. With the first-order least-squares functionals, corresponding LSFEMs can be defined. No integration by parts is needed to define the least-squares principle and thus the LSFEM, thus the first-order system LSFEM is ideal for the second order elliptic equation of non-divergence form.
Beside the obvious advantage of non-requirement of integration by parts, the LSFEM has other advantages. First, the least-squares weak formulation and its associated LSFEM using conforming finite element spaces are automatically coercive as long as the first-order system is well-posed. This is a significant advantage over other numerical methods since the well-posedness theory of the equation in non-divergence form is in general only sufficient. On the other hand, even without a rigorous mathematical proof, the elliptic equation in non-divergence form is often a result of some physical process that we are sure that a unique solution exists. Thus, in LSFEMs developed in this paper, we can reduce the condition of the PDE into a simple well-posedness without specifing the condition explicitly.
The other advantages of LSFEMs include conforming discretizations lead to stable and, ultimately, optimally accurate methods, the resulting algebraic problems are symmetric, positive definite, and can be solved by standard and robust iterative methods including multigrid methods.
The last important advantage of the LSFEMs is that it has a build-in a posteriori error estimator. The solution is probably singular due to the geometry of the domain or the coefficient matrix. Also, for problems like reaction diffusion equations, interior or boundary layers appear. To solve these problems efficiently, the a posteriori error estimator and adaptive mesh refinement algorithm are necessary.
In this paper, by introducing the gradient as an auxiliary variable, we first write the equation in non-divergence form into a system of first-order equations, then develop two least-squares minimization principles and two corresponding LSFEMs: one is based on an -norm square sum of the residuals and the other is based on a mesh-size weighted -norm square sum of the residuals. The two methods are called -LSFEM and weighted-LSFEM, respectively. For the -LSFEM, simpliest linear -finite elements are used to approximate both the solution and the gradient. For the weighted-LSFEM, the -finite element of degree , is used to approximate the solution, while the degree -finite element is used to approximate the gradient. Under the very weak assumption that the coefficient and domain is good enough to guarantee the existence and uniqueness of a solution, we show both continuous least-squares weak forms and their corresponding discrete problems are well-posed. A priori and a posteriori error estimates with respect to the least-squares norms are then discussed.
Numerical methods for non-divergence equation often use the following operator regularity assumption
to derive stability and error estimates. Unlike these papers, for the weighted-LSFEM, the error estimates of error in the -norm and the discrete broken -norm are investigated with a weaker assumption, see our discussion in section 4.2. Under stronger regularity assumptions, we show that the -norm of the error of the solution is one order higher than the least-squares norm of the error, providing the approximation degree for the solution is at least three. For the -LSFEM, we show the optimal and -error estimates with a solution regularity assumption. We perform extensive numerical experiments for smooth, non-smooth, and even degenerate coefficients on smooth and singular solutions to test the accuracy and efficiency of the proposed methods. With uniform refinements, we show the convergence orders match with the theory. With adaptive mesh refinements, optimal convergences results are obtained for singular solutions.
The LSFEM is well developed for the elliptic equations in divergence form, see for example, [9, 10, 5, 3, 8, 14]. A posteriori error estimates and adaptivity algorithms based on LSFEMs can be founded in [2, 13]. Compared the the LSFEMs for the elliptic equation in divergence form, the non-divergence equation has many differences in the stability analysis and choices of the finite element sub-spaces due to the non-divergence structure. We remark these differences in the various places of the paper as comparisons.
In a summary, the LSFEMs developed in this paper have several advantages compared to existing numerical methods: they are automatically stable under very mild assumption; they are easy to program due to that only simple Lagrange finite elements without jump terms are used; adaptive algorithms with the build-in a posteriori error estimators can handle problems with singular solutions or layers; under a condition on the operator regularity which is weaker than traditionally assumed, error estimates in standard norms are proved.
There are two least-squares finite element methods available for the non-divergence equation. None of them use a first-order reformulation. The paper  uses a second order least-squares formulation with -finite element approximations. The simple method developed by Ye and Mu in  uses -finite element spaces with orders higher than two and penalize the continuity of the solution and the normal component of the flux.
The remaining parts of this article are as follows: section 2 defines the first-order system least-squares weak problems and discusses their stabilities; section 3 presents the corresponding LSFEMs and their a priori and a posteriori error estimates in least-squares norms. Error estimates in other norms are discussed in sections 4 and 5 for the weighted and versions of methods, separately. Numerical experiments are presented in section 6.
Standard notation on function spaces applies throughout this article. Norms of functions in Lebesgue and Sobolev space () are denoted by . The subscript is omitted when . The inner product of real-valued matrices is denoted by . We use to denote the Hessian of .
2 First-Order System Least-Squares Weak Problems
2.1 Existence and uniqueness assumption
Define the solution space of (LABEL:eq_nondiv):
Notice that the space is weaker then , since we can expect that even though an individual , is not in , due to cancelation or good properties of , belongs to . For example, let and is the solution of the Poisson equation on an L-shaped domain , then but clearly .
We first state the assumption of the existence and uniqueness of the solution.
(Existence and uniqueness of the solution of the elliptic equation in non-divergence form) Assume that the coefficient matrix and the domain are nice enough, such that the equation (LABEL:eq_nondiv) has a unique solution for any .
There are various theories to ensure the existence and uniqueness of the equation, for example:
(Classical solution ) A classical solution exists if is Hölder continuous and if is sufficiently smooth.
(Strong solution ) If , a vanishing mean oscillation matrix with a uniform VMO-modulus of continuity, and if is of class , then there exists a unique solution to the problem.
( solution ) If the domain is convex and if satisfies the Cordes condition (LABEL:Cordes), then there exists a solution .
More detailed discussions can be found in the introduction of . We do find these theories are only sufficient theories. Besides many examples of Poisson equations on non-convex domains, the Test 3 problem from ( see also our numerical test 6.4) is an example that the matrix is not uniformly elliptic but a unique solution still exists.
2.2 Least-squares problems
Introduce a new gradient variable , we have the following first-order system:
It is clear that . For the gradient , the appropriate solution space is:
As a comparison, consider the equation in divergence form
It is well known that the flux for . The space is well studied [22, 4]. One of its property is that the normal component of its member vector function is continuous across the interfaces in the weak sense. Also, the negative divergence operator is the the dual operator of the gradient. These properties play important roles in the design of least-squares finite element methods for elliptic equations in divergence form. More importantly, -conforming finite element spaces such as Raviart-Thomas (RT) finite element space  with good approximation properties are also well known.
For the elliptic equation non-divergence form (LABEL:eq_nondiv), the property of the space is barely known. Similar to the space , for a vector function , we can expect that even though an individual , is not in , due to cancelation or good properties of , may belong to . But if we want to design a numerical method for a general non-divergence elliptic equation, we basically cannot assume or use any of these information, and we can not design an -intrinsic finite element subspace of as -conforming finite elements.
Let be a triangulation of using simplicial elements. The mesh is assumed to be shape-regular, but it does not to be quasi-uniform. Let be the diameter of the element .
We introduce two versions of least-squares functionals:
The functionals and are called the weighted version and the version, respectively. We use the notation to denote both and when two formulations can be presented in a unified framework and no confusion is caused.
The least-squares minimization problem is: seek , such that
The corresponding Euler-Lagrange formulations are: seek , such that
and find , such that
where for all and , the bilinear forms are defined as
The least-squares formulations can be easily extended to more general cases, for example, non-homogeneous Dirichlet boundary conditions and equations with convection and advection terms.
For example, for the general elliptic equation
with homogeneous boundary condition, let , the least-squares functional can be defined as:
To have a better robustness with respect to the coefficients, coefficient-weighted versions can also be used. For example, define the least-squares functional as:
where is a weight defined as a function of the coefficients , , and .
The -weighted least-squares functional can be defined similarly.
For the weighted functional , the -weight is on the term , similarly, we can also use,
as the weighted least-squares functional. For a uniform mesh, and are equivalent. But for an adaptively refined mesh, beahives more like a minimization problem with respect to the -norm of while is more like an optimization with respect to the -norm. We prefer the version in this paper since the minimum requirement of is not . Earlier discussion on the mesh-dependent least-squares methods can be found in .
The following are norms for :
We use to denote both versions when no confusion is caused.
To prove that defines a norm on , we only need to check conditions of a norm definition.
The linearity and the triangle inequality are obvious for .
If , due to the fact , we have , thus
in the sense. This means, , and
is true in the sense. By the existence and uniqueness of assumption of the solution Assumption LABEL:assump_EU, and . The norm is then well defined for both the weighted and versions of definition.
The condition is essential to the definition. This condition has the same role as the requirement of flux in for the equation in divergence form, which implicitly implies some weak continuity condition of its member functions.
It is also clear that
are semi-norms on an element .
The bilinear form or is continuous and coercive:
The lemma can be easily proved by a simple computation.
Assume that , the coefficient matrix and the domain are nice enough such that Assumption LABEL:assump_EU is true, then the least-squares problem (LABEL:LS) has a unique solution .
To prove the existence, for , by Assumption LABEL:assump_EU, there exists a unique solving the equation. Let , it is easy to that and , thus the least-squares functional (LABEL:LS) has a minimizer with minimum value zero. The minimizer is then the solution of the least-squares problem (LABEL:LS) and its corresponding Euler-Lagrange equation. The uniqueness is a simple consequence of the fact is a norm.
The above proof can be easily generalize to the case that the Drichelet boundary condition is not homogeneous.
The above argument to show the existence and uniqueness of the least-squares formulation is useful when the existences and uniqueness of the PDE is obtained from various non-variational techniques. A similar argument is used to prove the stability of least-squares formulations for the linear transport equation in .
Here, the assumption of the coefficient is quite weak. The matrix does not need to be in , it can even be degenerate as long as Assumption LABEL:assump_EU still holds.
A norm equivalence can be proved: there exists positive constant and , such that for :
Our least-squares functional essentially is a modification of . But due to the lack of the differentiability of , we cannot prove the following norm equivalence:
On the other hand, if is smooth enough, then (LABEL:eq_nondiv) can be written in the divergence form as (LABEL:eq_div), we do can prove (LABEL:equvalience) using the same technique for the equation divergence form with similar arguments in [8, 24].
However, for our least-squares functionals, we do have a one-sided bound, which can be easily proved for :
We do not use the minus- norm version in this paper due to its complicated discrete implementation, in stead, we choose a weighted mesh-dependent version to simplify the implementation and keep an optimal order of convergence.
3 Least-Squares Finite Element Methods
In this section, LSFEMs based on the least-squares minimization problems are developed. The a priori and a posteriori error estimates with respect to the least-squares norms are derived.
3.1 Least-squares finite element methods
For an element and an integer , let the space of polynomials with degrees less than or equal to . Define the finite element spaces and , , as follows:
We define the LSFEMs are follows.
(Weighted-LSFEM Problem) Seek , , such that
Or equivalently, find , , such that
(-LSFEM Problem) Seek , such that
Or equivalently, find , such that
The existence and uniqueness of the LSFEM problems are obvious from the facts that
For the approximation space of , the -conforming space is an obvious good choice. For the approximation space of , we use . The space is more restrictive than , but it has a simple conforming finite element space. The -conforming space is not the best choice if further information of is known. For example, in the case of the equation of divergence form, it is well known that and . Similarly for the non-divergence equation, for problems with low regularity and possible discontinuous coefficients, may not be in , then in some extreme case we will have even though is identical to the exact solution , cannot equal to at the same time, since may not be in . Such cases will pose problems of a posteriori error estimation and adaptive mesh refinements, see the discussions in  and  for the failure of the classical Zienkiewicz-Zhu error estimator, which recovers in to construct the a posteriori error estimator. Even though we do have this concern, due to the lack of information of and to keep the method suitable for a general coefficient matrix , the -conforming space to approximate is still a reasonable choice.
It is also worth to mention that, we also do not know whether or not for the non-divergence equation. In the the formulations suggested in  and , the normal jump of are used. Thus, these formulations have a similar possible inconsistency as our methods when applied to problems with less smoothness, where the exact jump of across element interfaces may not be zero for an exact , and have a possibility to introduce an extra error.
3.2 A priori error estimates
(Cea’s lemma type of result) Let be the solution of least-squares variational problem (LABEL:LS). Let , , be the solution of the weight-LSFEM problem (LABEL:LSFEM_h), the following best approximation result holds:
Let be the solution of the -LSFEM problem (LABEL:LSFEM_L2), the following best approximation result holds:
The proof of the best approximation result is standard.
(A priori error estimate for the weighted-LSFEM) Assume the solution , for some , and , be the solution of the weighted LSFEM problem (LABEL:LSFEM_h), then there exists a constant independent of the mesh size , such that
It is easy to see that
Then the a priori result is a direct consequence of Theorem LABEL:thm_bestapp and the approximation properties of functions in .
(A priori error estimate for the -LSFEM) Assume the solution and be the solution of the -LSFEM problem (LABEL:LSFEM_L2), then there exists a constant independent of the mesh size , such that
Then let be the interpolation of in and be the interpolation of in , by the approximation properties of and the fact that , we have
The theorem is proved.
From the a priori estimates, we can clearly see that with respect the least-squares norm, the weighted version is optimal when the regularity is high and a suitable high order finite element pair is used. For the -LSFEM, the optimal interpolation order for the -norm of is , which is one order high that the other two components, and thus sub-optimal, thus high order approximations of the -LSFEM is not suggested. But the -LSFEM can use the simplest linear conforming finite element space for , and has reasonable approximation orders, for example, we will find the -estimates of is the same order as the the weighted-LSFEM if assuming enough smoothness of the coefficient and the solution, see Theorems LABEL:L2_h and Theorems LABEL:L2_0 and our numerical tests.
In the a priori error estimates, in order to get a convergence order, we assume that the regularity of is at least or . To discuss the convergence without a high regularity, we first assume that the coefficient matrix is nice enough such that for any , there exists a , such that
Note that this condition is weaker than , since may not be in . If the assumption (LABEL:eps) is true, then a convergence result is easy to prove for only.
For the special case that , the restriction of the coefficient matrix on each , is a constant matrix, then contains a piecewise constant on each element , we have
This can be used to establish the convergence for a solution with a low regularity.
For the standard LSFEM for the divergence, such problems do not exist since Raviart-Thomas element is used to approximation , and the regularity requirement is on , which is weaker than the requirement on . Again, this is due to the special structure of the divergence equation, which is not available for the general non-divergence equation.
3.3 A posteriori error estimates
The least-squares functional can be used to define the following fully computable a posteriori local indicator and global error estimator.
Let be the solution of least-squares variational problem (LABEL:LS), and , be the solution of the weighted-LSFEM problem (LABEL:LSFEM_h), then define:
The a posteriori error estimator is exact with respect to the least-squares norm :
The following local efficiency bound is also true with a constant independent of the mesh size :
Using and , we obtain,
The proof of the local exactness is identical.
The locally efficiency (LABEL:eff_h) is a direct result of a local version of (LABEL:local_ineq_h).
For the -LSFEM, the a posteriori error estimator can be defined accordingly, and the corresponding results can be proved in a similar fashion.
Let be the solution of least-squares variational problem (LABEL:LS), and be the solution of the -LSFEM problem (LABEL:LSFEM_L2), define:
The a posteriori error estimator is exact with respect to the least-squares norm :
The following local efficiency bound is also true with independent of the mesh size :