Analysis of a finite element method for second order uniformly elliptic PDEs in non-divergence form

1. Problem Statement

The paper addresses the numerical analysis of two classes of second-order partial differential equations (PDEs) defined on Lipschitz polyhedral domains $\Omega \subset \mathbb{R}^d$ ($d=2, 3$):

1. Linear Elliptic PDE in Non-Divergence Form:
   $$A : D^2u + b \cdot \nabla u + cu = f \quad \text{in } \Omega, \quad u = 0 \text{ on } \partial\Omega$$
   where $A$ is a symmetric, uniformly positive definite matrix, and $b, c$ are bounded coefficients.

2. Uniformly Elliptic Hamilton-Jacobi-Bellman (HJB) Equation:
   $$\sup_{\alpha \in \Lambda} [A_\alpha : D^2u + b_\alpha \cdot \nabla u + c_\alpha u - f_\alpha] = 0 \quad \text{in } \Omega, \quad u = 0 \text{ on } \partial\Omega$$
   where $\Lambda$ is an index set, and the coefficients satisfy uniform ellipticity and boundedness conditions.

Key Challenges:

• Non-Divergence Form: Standard weak formulations (based on integration by parts) are not directly applicable because the second-order term $A:D^2u$ cannot be easily transferred to the test function without assuming high regularity of $A$ or the domain.
• Solution Space: The analysis focuses on strong solutions in the Sobolev space $W^{2,p}(\Omega)$ rather than viscosity solutions.
• Domain Geometry: The paper aims to handle Lipschitz polyhedral domains, including non-convex polygons in 2D, which are more general than the $C^{1,1}$ or convex domains often required in existing literature.
• Coefficient Regularity: The method must accommodate discontinuous coefficients (specifically for the HJB equation) while maintaining optimal convergence rates.

2. Methodology

The author proposes a unified Finite Element Method (FEM) for both the linear PDE and the HJB equation.

A. Discretization and Operators

• Mesh: A quasi-uniform, simplicial, conforming triangulation $\mathcal{T}_h$.
• Finite Element Space: $V_h = H^1_0(\Omega) \cap P_r(\mathcal{T}_h)$ with polynomial degree $r \geq 2$.
• Discrete Gradient ($\nabla_h$): The method utilizes a specific discrete gradient operator defined via integration by parts on the mesh elements and jump terms across interior faces. This operator allows the approximation of the Hessian $D^2u$ without requiring $C^1$ continuity of the basis functions.
  $$(\nabla_h v)_i = \partial_{h, x_i} v$$
• Weight Functions: To handle the non-divergence form, the equations are multiplied by a weight function $\gamma$ (or $\gamma_\alpha$ for HJB).
  • For the linear case: $\gamma(x) = \frac{\text{Tr} A(x)}{|A(x)|^2}$.
  • This transforms the equation into a form amenable to the discrete operator: $\tilde{A} : D^2 u + \dots = \tilde{f}$.

B. The Numerical Schemes

• Linear Case (Eq. 14): Find $u_h \in V_h$ such that for all $v_h \in V_h$:
  $$\int_\Omega (\tilde{A} : \nabla_h(\nabla u_h) + \tilde{b} \cdot \nabla u_h + \tilde{c} u_h) v_h \, dx = \int_\Omega \tilde{f} v_h \, dx$$
• HJB Case (Eq. 16): Find $u_h \in V_h$ such that for all $v_h \in V_h$:
  $$\left( \sup_{\alpha \in \Lambda} \gamma_\alpha [A_\alpha : \nabla_h(\nabla u_h) + b_\alpha \cdot \nabla u_h + c_\alpha u_h - f_\alpha], v_h \right)_\Omega = 0$$
  Note: If $\Lambda$ contains a single element, the HJB scheme reduces exactly to the linear scheme.

C. Analytical Framework

The analysis relies on discrete $W^{2,p}$-norms defined as:
$$\|v\|_{W^{2,p}_h(\Omega)} = \|D^2_h v\|_{L^p(\Omega)} + \left( \sum_{F \in \mathcal{F}_h^I} h_F^{1-p} \|\llbracket \nabla v \rrbracket\|_{L^p(F)}^p \right)^{1/p}$$
Key analytical tools include:

1. Global $W^{2,p}$ Estimates: Proving well-posedness of the strong solution for the continuous problem on Lipschitz polyhedra.
2. Discrete Compactness: Establishing that bounded sequences in the discrete $W^{2,p}_h$ norm have convergent subsequences in $W^{1,p}$.
3. Uniqueness Arguments: Using the uniqueness of the strong solution to prove stability via contradiction (avoiding specialized duality arguments used in previous works).
4. Measurability of Supremum: For the HJB equation, proving that the supremum over the index set $\Lambda$ is measurable and can be approximated by a countable subset $\tilde{\Lambda}$ (Assumption 6).

3. Key Contributions

A. Unified Framework

The paper demonstrates that a single FEM formulation can effectively solve both linear non-divergence PDEs and fully nonlinear HJB equations. The HJB scheme naturally degenerates to the linear scheme when the index set is a singleton.

B. Relaxation of Domain and Coefficient Assumptions

• Domains: The analysis extends optimal convergence results from $C^{1,1}$ or strictly convex domains to Lipschitz polyhedral domains.
  • For $d=3$, $\Omega$ must be convex.
  • For $d=2$, $\Omega$ can be non-convex (with restrictions on $p$).
• Coefficients: The method allows for discontinuous coefficients in the HJB equation, provided they satisfy a specific approximation condition (Assumption 6) involving a countable dense subset. This relaxes the continuity assumptions ($A \in C^0$) required by many existing DG methods.

C. Optimal Convergence in $W^{2,p}$

The paper proves optimal convergence rates in the discrete $W^{2,p}_h$-norm for a wide range of $p$:

• Convex Polyhedra ($d=2,3$): $p \in (1, 2]$.
• Non-convex Polygons ($d=2$): $p \in (4/3 - \epsilon_1, 4/3 + \epsilon_2)$.
  This is a significant improvement over methods that only guarantee convergence in $L^2$ or $H^2$ (where $p=2$).

D. Novel Stability Proof

Unlike previous works (e.g., Feng et al.) that rely on specialized duality arguments contingent on the well-posedness of the continuous problem, this paper establishes stability using discrete compactness and the uniqueness of the strong solution. This approach is more robust for non-convex domains and discontinuous coefficients.

4. Main Results

1. Well-Posedness (Theorems 3.3, 4.1):

   • Existence and uniqueness of strong solutions in $W^{2,p}(\Omega) \cap W^{1,p}_0(\Omega)$ are proven for both the linear PDE and the HJB equation under the specified conditions on $A, b, c$ and the domain.
   • For the HJB equation, the Cordes condition is generalized to allow discontinuous coefficients via the countable subset assumption.

2. Stability (Theorems 3.2, 4.3):

   • The discrete solution $u_h$ satisfies a stability estimate: $\|u_h\|_{W^{2,p}_h} \leq C \|f\|_{L^p}$.
   • This stability holds for the specific range of $p$ dependent on the domain geometry.

3. Optimal Convergence (Theorems 3.3, 4.3):

   • The error estimate is derived:
     $$\|u - u_h\|_{W^{1,p}} + \|u - u_h\|_{W^{2,p}_h} \leq C \left( \inf_{\chi_h \in V_h} (\|u - \chi_h\|_{W^{1,p}} + \|u - \chi_h\|_{W^{2,p}_h}) + \|P_h(D^2u) - D^2u\|_{L^p} \right)$$
   • This implies that if the exact solution is smooth enough, the method achieves optimal convergence rates in the discrete $W^{2,p}$ norm.

5. Significance

• Bridging Theory and Practice: The paper provides a rigorous theoretical foundation for using standard finite element spaces (with $C^0$ continuity) to solve non-divergence form PDEs, which are notoriously difficult to discretize.
• Handling Non-Convexity: By extending the analysis to non-convex 2D polygons, the method becomes applicable to a broader class of practical engineering and physical problems where domains are rarely perfectly convex.
• HJB Applications: The relaxation of coefficient continuity assumptions makes the method suitable for stochastic control problems where coefficients may be discontinuous (e.g., switching regimes), a scenario where previous $C^0$-conforming methods failed.
• Unified Approach: Demonstrating that one numerical scheme handles both linear and fully nonlinear cases simplifies implementation and theoretical analysis for complex systems.

In summary, Weifeng Qiu's work advances the state-of-the-art in numerical PDEs by establishing a robust, unified finite element framework for non-divergence form equations that achieves optimal convergence in strong norms on general polyhedral domains with relaxed regularity requirements.