Adaptive FirstOrder System LeastSquares Finite Element Methods for Second Order Elliptic Equations in NonDivergence Form
Abstract
This paper studies adaptive firstorder leastsquares finite element methods for secondorder elliptic partial differential equations in nondivergence form. Unlike the classical finite element method which uses weak formulations of PDEs not applicable for the nondivergence equation, the firstorder leastsquares formulations naturally have stable weak forms without using integration by parts, allow simple finite element approximation spaces, and have buildin a posteriori error estimators for adaptive mesh refinements.
The nondivergence equation is first written as a system of firstorder equations by introducing the gradient as a new variable. Then two versions of leastsquares finite element methods using simple finite elements are developed in the paper, one is the LSFEM which uses linear elements, the other is the weightedLSFEM with a meshdependent 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 weightedLSFEM are also discussed. error estimates are derived for both formulations. We perform extensive numerical experiments for smooth, nonsmooth, and even degenerate coefficients on smooth and singular solutions to test the accuracy and efficiency of the proposed methods.
Key words. nondivergence elliptic equation, firstorder system leastsquares, finite element method, a priori error estimate, a posteriori error estimate, adaptive method.
1 Introduction
In this paper we consider finite element approximations of the following elliptic PDE in nondivergence form:
(1.1)  
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 rowwise:
(1.2)  
The elliptic PDE in nondivergence 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 elementwisely 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 nondivergence 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 nonunique, see an example given in [30]. 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 nondivergence 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 nondivergence form. Based on discrete CalderonZygmond estimates, Feng, Neilan, and coauthors 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 [31], a mixed method [20], and a nonsymmetric method [29] are developed. A weak Galerkin method is developed by Wang and Wang in [32]. 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 nondivergence operator second order and borrow techniques from variational fourthorder problems. Nochetto and Zhang [30] studied a twoscale method, which is based on the integrodifferential 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 firstorder 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 leastsquares finite element method (LSFEM). The firstorder system leastsquares principle first rewrite the PDE into a firstorder system, then define an artificial, externally defined energytype principle. The energy functional can be defined as summations of weighted residuals of the system. With the firstorder leastsquares functionals, corresponding LSFEMs can be defined. No integration by parts is needed to define the leastsquares principle and thus the LSFEM, thus the firstorder system LSFEM is ideal for the second order elliptic equation of nondivergence form.
Beside the obvious advantage of nonrequirement of integration by parts, the LSFEM has other advantages. First, the leastsquares weak formulation and its associated LSFEM using conforming finite element spaces are automatically coercive as long as the firstorder system is wellposed. This is a significant advantage over other numerical methods since the wellposedness theory of the equation in nondivergence form is in general only sufficient. On the other hand, even without a rigorous mathematical proof, the elliptic equation in nondivergence 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 wellposedness 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 buildin 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 nondivergence form into a system of firstorder equations, then develop two leastsquares minimization principles and two corresponding LSFEMs: one is based on an norm square sum of the residuals and the other is based on a meshsize weighted norm square sum of the residuals. The two methods are called LSFEM and weightedLSFEM, respectively. For the LSFEM, simpliest linear finite elements are used to approximate both the solution and the gradient. For the weightedLSFEM, 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 leastsquares weak forms and their corresponding discrete problems are wellposed. A priori and a posteriori error estimates with respect to the leastsquares norms are then discussed.
Numerical methods for nondivergence equation often use the following operator regularity assumption
to derive stability and error estimates. Unlike these papers, for the weightedLSFEM, 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 leastsquares 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, nonsmooth, 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 nondivergence equation has many differences in the stability analysis and choices of the finite element subspaces due to the nondivergence 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 buildin 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 leastsquares finite element methods available for the nondivergence equation. None of them use a firstorder reformulation. The paper [20] uses a second order leastsquares formulation with finite element approximations. The simple method developed by Ye and Mu in [27] 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 firstorder system leastsquares weak problems and discusses their stabilities; section 3 presents the corresponding LSFEMs and their a priori and a posteriori error estimates in leastsquares 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 realvalued matrices is denoted by . We use to denote the Hessian of .
2 FirstOrder System LeastSquares Weak Problems
2.1 Existence and uniqueness assumption
Define the solution space of (LABEL:eq_nondiv):
(2.1) 
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 Lshaped domain , then but clearly .
We first state the assumption of the existence and uniqueness of the solution.
Assumption 2.1.
(Existence and uniqueness of the solution of the elliptic equation in nondivergence 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 .
Remark 2.2.
There are various theories to ensure the existence and uniqueness of the equation, for example:

(Classical solution [23]) A classical solution exists if is Hölder continuous and if is sufficiently smooth.

(Strong solution [16]) If , a vanishing mean oscillation matrix with a uniform VMOmodulus of continuity, and if is of class , then there exists a unique solution to the problem.

( solution [26]) 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 [30]. We do find these theories are only sufficient theories. Besides many examples of Poisson equations on nonconvex domains, the Test 3 problem from [18]( 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 Leastsquares problems
Introduce a new gradient variable , we have the following firstorder system:
(2.2) 
It is clear that . For the gradient , the appropriate solution space is:
Remark 2.3.
As a comparison, consider the equation in divergence form
(2.3)  
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 leastsquares finite element methods for elliptic equations in divergence form. More importantly, conforming finite element spaces such as RaviartThomas (RT) finite element space [4] with good approximation properties are also well known.
For the elliptic equation nondivergence 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 nondivergence 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 shaperegular, but it does not to be quasiuniform. Let be the diameter of the element .
We introduce two versions of leastsquares functionals:
(2.4)  
(2.5) 
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 leastsquares minimization problem is: seek , such that
(2.6) 
The corresponding EulerLagrange formulations are: seek , such that
(2.7) 
and find , such that
(2.8) 
where for all and , the bilinear forms are defined as
Remark 2.4.
The leastsquares formulations can be easily extended to more general cases, for example, nonhomogeneous Dirichlet boundary conditions and equations with convection and advection terms.
For example, for the general elliptic equation
(2.9) 
with homogeneous boundary condition, let , the leastsquares functional can be defined as:
To have a better robustness with respect to the coefficients, coefficientweighted versions can also be used. For example, define the leastsquares functional as:
where is a weight defined as a function of the coefficients , , and .
The weighted leastsquares functional can be defined similarly.
Remark 2.5.
For the weighted functional , the weight is on the term , similarly, we can also use,
as the weighted leastsquares 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 meshdependent leastsquares methods can be found in [1].
Lemma 2.6.
The following are norms for :
We use to denote both versions when no confusion is caused.
Proof.
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.
Remark 2.7.
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.
Remark 2.8.
It is also clear that
are seminorms on an element .
Lemma 2.9.
The bilinear form or is continuous and coercive:
(2.10)  
(2.11) 
The lemma can be easily proved by a simple computation.
Theorem 2.10.
Assume that , the coefficient matrix and the domain are nice enough such that Assumption LABEL:assump_EU is true, then the leastsquares problem (LABEL:LS) has a unique solution .
Proof.
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 leastsquares functional (LABEL:LS) has a minimizer with minimum value zero. The minimizer is then the solution of the leastsquares problem (LABEL:LS) and its corresponding EulerLagrange equation. The uniqueness is a simple consequence of the fact is a norm.
Remark 2.11.
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 leastsquares formulation is useful when the existences and uniqueness of the PDE is obtained from various nonvariational techniques. A similar argument is used to prove the stability of leastsquares formulations for the linear transport equation in [25].
Remark 2.12.
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.
Remark 2.13.
For the elliptic equation in divergence form (LABEL:eqdiv), traditionally there are two forms on leastsquares functionals [9, 5]:
A norm equivalence can be proved: there exists positive constant and , such that for :
and
Our leastsquares functional essentially is a modification of . But due to the lack of the differentiability of , we cannot prove the following norm equivalence:
(2.12) 
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 leastsquares functionals, we do have a onesided bound, which can be easily proved for :
(2.13)  
(2.14) 
We do not use the minus norm version in this paper due to its complicated discrete implementation, in stead, we choose a weighted meshdependent version to simplify the implementation and keep an optimal order of convergence.
3 LeastSquares Finite Element Methods
In this section, LSFEMs based on the leastsquares minimization problems are developed. The a priori and a posteriori error estimates with respect to the leastsquares norms are derived.
3.1 Leastsquares 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.
(WeightedLSFEM Problem) Seek , , such that
(3.1) 
Or equivalently, find , , such that
(3.2) 
(LSFEM Problem) Seek , such that
(3.3) 
Or equivalently, find , such that
(3.4) 
The existence and uniqueness of the LSFEM problems are obvious from the facts that
and .
Remark 3.1.
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 nondivergence 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 [12] and [7] for the failure of the classical ZienkiewiczZhu 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 nondivergence equation. In the the formulations suggested in [18] and [27], 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
Theorem 3.2.
(Cea’s lemma type of result) Let be the solution of leastsquares variational problem (LABEL:LS). Let , , be the solution of the weightLSFEM problem (LABEL:LSFEM_h), the following best approximation result holds:
(3.5) 
Let be the solution of the LSFEM problem (LABEL:LSFEM_L2), the following best approximation result holds:
(3.6) 
Proof.
The proof of the best approximation result is standard.
Theorem 3.3.
(A priori error estimate for the weightedLSFEM) 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
(3.7) 
Proof.
It is easy to see that
(3.8) 
Then the a priori result is a direct consequence of Theorem LABEL:thm_bestapp and the approximation properties of functions in .
Theorem 3.4.
(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
(3.9) 
Proof.
We have
(3.10) 
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.
Remark 3.5.
From the a priori estimates, we can clearly see that with respect the leastsquares 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 suboptimal, 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 weightedLSFEM if assuming enough smoothness of the coefficient and the solution, see Theorems LABEL:L2_h and Theorems LABEL:L2_0 and our numerical tests.
Remark 3.6.
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
(3.11) 
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 RaviartThomas 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 nondivergence equation.
3.3 A posteriori error estimates
The leastsquares functional can be used to define the following fully computable a posteriori local indicator and global error estimator.
3.3.1 WeightedLSFEM
Let be the solution of leastsquares variational problem (LABEL:LS), and , be the solution of the weightedLSFEM problem (LABEL:LSFEM_h), then define:
Theorem 3.7.
The a posteriori error estimator is exact with respect to the leastsquares norm :
The following local efficiency bound is also true with a constant independent of the mesh size :
(3.12) 
Proof.
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).
3.3.2 Lsfem
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 leastsquares variational problem (LABEL:LS), and be the solution of the LSFEM problem (LABEL:LSFEM_L2), define:
and
Theorem 3.8.
The a posteriori error estimator is exact with respect to the leastsquares norm :
The following local efficiency bound is also true with independent of the mesh size :