Higher-order Adaptive Finite Difference Methods for Fully Nonlinear Elliptic Equations

TL;DR

This paper presents innovative higher-order adaptive finite difference methods for solving fully nonlinear elliptic PDEs, emphasizing flexibility, efficiency, and convergence on complex geometries.

Contribution

It introduces generalized finite difference schemes on adaptive, piecewise Cartesian meshes with boundary points, enabling higher-order accuracy and efficient computation.

Findings

Methods demonstrate high accuracy in computational examples.

Adaptive meshes improve solution flexibility for complex geometries.

Algorithms efficiently compute non-traditional finite difference stencils.

Abstract

We introduce generalised finite difference methods for solving fully nonlinear elliptic partial differential equations. Methods are based on piecewise Cartesian meshes augmented by additional points along the boundary. This allows for adaptive meshes and complicated geometries, while still ensuring consistency, monotonicity, and convergence. We describe an algorithm for efficiently computing the non-traditional finite difference stencils. We also present a strategy for computing formally higher-order convergent methods. Computational examples demonstrate the efficiency, accuracy, and flexibility of the methods.

Figures

Tables

Table 1. Table 1. Convergence results for the C 2 superscript 𝐶 2 C^{2} solution of the Monge-Ampère equation.

	Errors and order, Example 5.1 (circle)
$N$	$h$	Monotone		Filtered
101	$2.500\text{×}{10}^{-01}$	$1.739\text{×}{10}^{-02}$	-	$7.008\text{×}{10}^{-03}$	-
349	$1.250\text{×}{10}^{-01}$	$5.690\text{×}{10}^{-03}$	1.61	$2.219\text{×}{10}^{-03}$	1.66
1137	$6.250\text{×}{10}^{-02}$	$2.780\text{×}{10}^{-03}$	1.03	$6.184\text{×}{10}^{-04}$	1.84
4289	$3.125\text{×}{10}^{-02}$	$1.876\text{×}{10}^{-03}$	0.57	$1.708\text{×}{10}^{-04}$	1.86
15685	$1.562\text{×}{10}^{-02}$	$1.630\text{×}{10}^{-03}$	0.20	$4.320\text{×}{10}^{-05}$	1.98
58449	$7.812\text{×}{10}^{-03}$	$1.567\text{×}{10}^{-03}$	0.06	$1.082\text{×}{10}^{-05}$	2.00
	Errors and order, Example 5.1 (ellipse)
$N$	$h$	Monotone		Filtered
79	$2.500\text{×}{10}^{-01}$	$9.467\text{×}{10}^{-03}$	-	$8.884\text{×}{10}^{-03}$	-
257	$1.250\text{×}{10}^{-01}$	$3.026\text{×}{10}^{-03}$	1.65	$1.233\text{×}{10}^{-03}$	2.85
893	$6.250\text{×}{10}^{-02}$	$1.262\text{×}{10}^{-03}$	1.26	$3.628\text{×}{10}^{-04}$	1.76
3143	$3.125\text{×}{10}^{-02}$	$8.277\text{×}{10}^{-04}$	0.61	$1.012\text{×}{10}^{-04}$	1.84
11257	$1.562\text{×}{10}^{-02}$	$7.443\text{×}{10}^{-04}$	0.15	$2.602\text{×}{10}^{-05}$	1.96
42863	$7.812\text{×}{10}^{-03}$	$7.238\text{×}{10}^{-04}$	0.04	$6.512\text{×}{10}^{-06}$	2.00
	Errors and order, Example 5.1 (diamond)
$N$	$h$	Monotone		Filtered
57	$2.500\text{×}{10}^{-01}$	$6.162\text{×}{10}^{-03}$	-	$7.244\text{×}{10}^{-03}$	-
177	$1.250\text{×}{10}^{-01}$	$2.332\text{×}{10}^{-03}$	1.40	$1.478\text{×}{10}^{-03}$	2.29
673	$6.250\text{×}{10}^{-02}$	$9.716\text{×}{10}^{-04}$	1.26	$3.561\text{×}{10}^{-04}$	2.05
2369	$3.125\text{×}{10}^{-02}$	$3.413\text{×}{10}^{-04}$	1.51	$8.881\text{×}{10}^{-05}$	2.00
9089	$1.562\text{×}{10}^{-02}$	$2.771\text{×}{10}^{-04}$	0.30	$2.218\text{×}{10}^{-05}$	2.00
35073	$7.812\text{×}{10}^{-03}$	$1.878\text{×}{10}^{-04}$	0.56	$5.544\text{×}{10}^{-06}$	2.00
	Errors and order, Example 5.1 (diamond stretched)
$N$	$h$	Monotone		Filtered
61	$2.500\text{×}{10}^{-01}$	$2.535\text{×}{10}^{-03}$	-	$1.888\text{×}{10}^{-03}$	-
153	$1.250\text{×}{10}^{-01}$	$7.891\text{×}{10}^{-04}$	1.68	$5.212\text{×}{10}^{-04}$	1.86
497	$6.250\text{×}{10}^{-02}$	$2.802\text{×}{10}^{-04}$	1.49	$1.356\text{×}{10}^{-04}$	1.94
1633	$3.125\text{×}{10}^{-02}$	$1.516\text{×}{10}^{-04}$	0.89	$3.428\text{×}{10}^{-05}$	1.98
5569	$1.562\text{×}{10}^{-02}$	$7.684\text{×}{10}^{-05}$	0.98	$8.597\text{×}{10}^{-06}$	2.00
20353	$7.812\text{×}{10}^{-03}$	$2.357\text{×}{10}^{-05}$	1.70	$2.151\text{×}{10}^{-06}$	2.00

Table 2. Table 2. Convergence results for the C 1 superscript 𝐶 1 C^{1} solution of the Monge-Ampère equation.

	Errors and order, Example 5.2 (circle)
$N$	$h$	Monotone		Filtered
101	$2.500\text{×}{10}^{-01}$	$1.606\text{×}{10}^{-02}$	-	$4.142\text{×}{10}^{-03}$	-
349	$1.250\text{×}{10}^{-01}$	$9.106\text{×}{10}^{-03}$	0.82	$2.384\text{×}{10}^{-03}$	0.80
1137	$6.250\text{×}{10}^{-02}$	$5.706\text{×}{10}^{-03}$	0.67	$1.763\text{×}{10}^{-03}$	0.44
4289	$3.125\text{×}{10}^{-02}$	$4.827\text{×}{10}^{-03}$	0.24	$7.871\text{×}{10}^{-04}$	1.16
15685	$1.562\text{×}{10}^{-02}$	$4.330\text{×}{10}^{-03}$	0.16	$3.674\text{×}{10}^{-04}$	1.10
58449	$7.812\text{×}{10}^{-03}$	$4.300\text{×}{10}^{-03}$	0.01	$1.839\text{×}{10}^{-04}$	1.00
	Errors and order, Example 5.2 (ellipse)
$N$	$h$	Monotone		Filtered
79	$2.500\text{×}{10}^{-01}$	$1.193\text{×}{10}^{-02}$	-	$4.276\text{×}{10}^{-03}$	-
257	$1.250\text{×}{10}^{-01}$	$7.126\text{×}{10}^{-03}$	0.74	$3.201\text{×}{10}^{-03}$	0.42
893	$6.250\text{×}{10}^{-02}$	$4.569\text{×}{10}^{-03}$	0.64	$1.579\text{×}{10}^{-03}$	1.02
3143	$3.125\text{×}{10}^{-02}$	$3.924\text{×}{10}^{-03}$	0.22	$8.118\text{×}{10}^{-04}$	0.96
11257	$1.562\text{×}{10}^{-02}$	$3.519\text{×}{10}^{-03}$	0.16	$3.764\text{×}{10}^{-04}$	1.11
42863	$7.812\text{×}{10}^{-03}$	$3.497\text{×}{10}^{-03}$	0.01	$2.224\text{×}{10}^{-04}$	0.76
	Errors and order, Example 5.2 (diamond)
$N$	$h$	Monotone		Filtered
57	$2.500\text{×}{10}^{-01}$	$8.503\text{×}{10}^{-03}$	-	$6.576\text{×}{10}^{-03}$	-
177	$1.250\text{×}{10}^{-01}$	$6.670\text{×}{10}^{-03}$	0.35	$2.155\text{×}{10}^{-03}$	1.61
673	$6.250\text{×}{10}^{-02}$	$4.265\text{×}{10}^{-03}$	0.65	$1.555\text{×}{10}^{-03}$	0.47
2369	$3.125\text{×}{10}^{-02}$	$3.442\text{×}{10}^{-03}$	0.31	$7.661\text{×}{10}^{-04}$	1.02
9089	$1.562\text{×}{10}^{-02}$	$2.303\text{×}{10}^{-03}$	0.58	$3.597\text{×}{10}^{-04}$	1.09
35073	$7.812\text{×}{10}^{-03}$	$1.109\text{×}{10}^{-03}$	1.05	$2.142\text{×}{10}^{-04}$	0.75
	Errors and order, Example 5.2 (diamond stretched)
$N$	$h$	Monotone		Filtered
61	$2.500\text{×}{10}^{-01}$	$2.653\text{×}{10}^{-03}$	-	$4.060\text{×}{10}^{-03}$	-
153	$1.250\text{×}{10}^{-01}$	$3.780\text{×}{10}^{-03}$	-0.51	$3.026\text{×}{10}^{-03}$	0.42
497	$6.250\text{×}{10}^{-02}$	$2.426\text{×}{10}^{-03}$	0.64	$1.213\text{×}{10}^{-03}$	1.32
1633	$3.125\text{×}{10}^{-02}$	$1.820\text{×}{10}^{-03}$	0.41	$1.228\text{×}{10}^{-03}$	-0.02
5569	$1.562\text{×}{10}^{-02}$	$1.203\text{×}{10}^{-03}$	0.60	$3.726\text{×}{10}^{-04}$	1.72
20353	$7.812\text{×}{10}^{-03}$	$6.089\text{×}{10}^{-04}$	0.98	$2.015\text{×}{10}^{-04}$	0.89

Equations

$$F(x,u(x),D^{2}u(x))\equiv \underset{\theta \in [0,2\pi )}{\max }{F}_{\theta }(x,u(x),u_{\theta \theta }(x))=0,$$

$$\underset{\theta \in \mathcal{A}\subset [0,2\pi )}{\max }{F}_{\theta }\left(x_i,u(x_i),\sum _{j\in \mathcal{N}(i,\theta )}{a}_{i,j,\theta }(u(x_i)-u(x_j))\right)=0,x_i\in \mathcal{G}.$$

$$F(x,u(x),D^{2}u(x))=0,x\in \bar{\Omega }\subset {\mathbb{R}}^2.$$

$$F(x,u,X)\le F(x,v,Y)$$

$$F(x,u(x),D^{2}u(x))=u(x)-g(x),x\in \partial \Omega .$$

$$u^{\ast }(x)=\underset{y\to x}{\mathrm{limsup}}u(y),$$

$$u_{\ast }(x)=\underset{y\to x}{\mathrm{liminf}}u(y).$$

$$F_{\ast }^{(\ast )}(x,u(x),D^2\varphi (x))\le (\ge )0.$$

$$F^{ϵ}(x,u(x),u(x)-u(\cdot ))=0$$

$$\underset{ϵ\to 0^{+},y\to x,\xi \to 0}{\mathrm{limsup}}{F}^{ϵ}(y,\varphi (y)+\xi ,\varphi (y)-\varphi (\cdot ))\le F^{\ast }(x,\varphi (x),\nabla \varphi (x),D^2\varphi (x)),$$

$$\underset{ϵ\to 0^{+},y\to x,\xi \to 0}{\mathrm{liminf}}{F}^{ϵ}(y,\varphi (y)+\xi ,\varphi (y)-\varphi (\cdot ))\ge F_{\ast }(x,\varphi (x),\nabla \varphi (x),D^2\varphi (x)).$$

$$x_j\in \mathrm{argmin}\left\{{\sin }^2\varphi \mid (r,\varphi )\in {\mathcal{G}}^h\cap B(x_0,ϵ)\text{ is in the jth quadrant}\right\}$$

$${\mathcal{D}}_{\theta \theta }^{h}u(x_0)=\sum _{j=1}^4{a}_j(u(x_j)-u(x_0))$$

$$S_j=r_j\sin {\varphi }_j,C_j=r_j\cos {\varphi }_j$$

$$\begin{array}{rl}{a}_1& =\frac{2S_4(C_3{S}_2-C_2{S}_3)}{(C_3{S}_2-C_2{S}_3)(C_1^2{S}_4-C_4^2{S}_1)-(C_1{S}_4-C_4{S}_1)(C_3^2{S}_2-C_2^2{S}_3)}\\ a_2& =\frac{2S_3(C_1{S}_4-C_4{S}_1)}{(C_3{S}_2-C_2{S}_3)(C_1^2{S}_4-C_4^2{S}_1)-(C_1{S}_4-C_4{S}_1)(C_3^2{S}_2-C_2^2{S}_3)}\\ a_3& =\frac{-2S_2(C_1{S}_4-C_4{S}_1)}{(C_3{S}_2-C_2{S}_3)(C_1^2{S}_4-C_4^2{S}_1)-(C_1{S}_4-C_4{S}_1)(C_3^2{S}_2-C_2^2{S}_3)}\\ a_4& =\frac{-2S_1(C_3{S}_2-C_2{S}_3)}{(C_3{S}_2-C_2{S}_3)(C_1^2{S}_4-C_4^2{S}_1)-(C_1{S}_4-C_4{S}_1)(C_3^2{S}_2-C_2^2{S}_3)}.\end{array}$$

$$F_i[u]\equiv \underset{\theta \in \mathcal{A}}{\max }{F}_{\theta }\left(x_i,u(x_i),\sum _{j\in \mathcal{N}(i,\theta )}{a}_{i,j,\theta }(u(x_i)-u(x_j))\right)=0,x_i\in \mathcal{G}.$$

$$h_{B,S}^n=\underset{x\in \partial \Omega \cap S}{\sup }\underset{y\in {\mathcal{G}}^n\cap \partial \Omega \cap S}{\min }|x-y|,{\delta }_S^n=\underset{x\in \Omega \cap {\mathcal{G}}^n}{\min }\underset{y\in \partial \Omega \cap S}{\inf }|x-y|.$$

$$F_F[u]=F_M[u]+h^{\alpha }S\left(\frac{F_A[u]-F_M[u]}{h^{\alpha }}\right),$$

$$S(x)=\{\begin{array}{ll}x,& |x|\le 1,\\ 0,& |x|\ge 2,\\ -x+2,& 1\le x\le 2,\\ -x-2,& -2\le x\le -1.\end{array}$$

$$u_{xx}$$

$$u_{yy}$$

$$u_{xy}$$

$$u_{xx}$$

$$u_{yy}$$

$$u_{xy}$$

$$u_{xx}$$

$$u_{yy}$$

$$u_{xy}$$

$$u(x_i)-u(x_0)=\sum _{0<|\alpha |\le 3}\frac{(x_i-x_0{)}^{\alpha }}{\alpha !}({\partial }^{\alpha }u)(x_0)+\mathcal{O}(h^4),$$

$$\sum _{i=1}^N{a}_i(u(x_i)-u(x_0))={\partial }^{\beta }u(x_0)+\mathcal{O}(h^2)$$

$$\sum _{i=1}^N\frac{(x_i-x_0{)}^{\alpha }}{\alpha !}{a}_i=1_{\{\alpha =\beta \}}$$

$$\{\begin{array}{ll}-\det (D^{2}u)+f=0,& x\in \Omega \\ u=g,& x\in \partial \Omega \\ \text{u is convex.}& \end{array}$$

$$-\underset{\theta \in [0,\pi /2)}{\min }\left\{\max \left\{\frac{{\partial }^{2}u}{\partial e_{\theta }^2},0\right\}\max \left\{\frac{{\partial }^{2}u}{\partial e_{\theta +\pi /2}^2},0\right\}+\min \left\{\frac{{\partial }^{2}u}{\partial e_{\theta }^2},0\right\}+\max \left\{\frac{{\partial }^{2}u}{\partial e_{\theta +\pi /2}^2},0\right\}\right\}+f=0.$$

Full text

The first author was partially supported by NSF DMS-1619807. The second author was partially supported by NSERC Discovery grant RGPIN-2016-03922 and by Fundação para a Ciência e Tecnologia (FCT) doctoral grant (SFRH/BD/84041/2012).

1. Introduction

In this article we design generalised finite difference methods for a large class of fully nonlinear degenerate elliptic partial differential equations (PDEs). The approximation schemes are almost-monotone, which allows us to exploit the Barles-Souganidis convergence framework. A key feature of these methods is the use of piecewise Cartesian grids, augmented by additional discretisation points along the boundary. Because of the underlying structure of the grids, the methods we design overcome several hurdles that exist for current numerical methods. 1) The non-standard finite difference stencils can be constructed efficiently. 2) The monotone approximation schemes preserve consistency near the boundary. 3) Higher-order schemes are easily designed. 4) Complicated geometries are easily handled.

1.1. Background

Fully nonlinear elliptic partial differential equations (PDEs) arise in numerous applications including reflector/refractor design [GO03], meteorology [CNP91], differential geometry [CM99], astrophysics [FMMS02], seismology [EF14], mesh generation [BW09], computer graphics [OP03], and mathematical finance [FS06]. Realistic applications often involve complicated domains and highly non-smooth data. Thus the development of robust numerical methods that are compatible with adaptive mesh refinement is a priority.

In recent years, the numerical solution of these equations has received a great deal of attention, and several new methods have been developed including finite difference methods [BFO10, FDC08, LR05, SAK15, SWR11], finite element methods [Awa15, Böh08, BGNS11, SS14], least squares methods [DG06], and methods involving fourth-order regularisation terms [FN09]. However, these methods are not designed to compute weak solutions. When the ellipticity of the equation is degenerate or no smooth solution exists, methods become very slow, are unstable, or converge to an incorrect solution.

Using a framework developed by Barles and Souganidis [BS91], provably convergent (monotone) methods have recently been constructed for several fully nonlinear equations using wide finite difference stencils [Obe08a]. Recently, the idea of wide stencil schemes has been adapted to produce monotone approximations for a large class of fully nonlinear elliptic operators posed on very general meshes or points clouds [Fro17b]. However, constructing the necessary finite difference stencils is a very expensive process, and the resulting approximations have only $\mathcal{O}(\sqrt{h})$ accuracy.

1.2. Contributions of this work

The goal of this article is to design higher-order, adaptive, convergent finite difference methods for the solution of the degenerate elliptic PDE

[TABLE]

where $u_{\theta\theta}$ denotes the second directional derivative of $u$ in the direction $e_{\theta}=(\cos\theta,\sin\theta)\in\mathbb{R}^{2}$. This includes a wide range of PDE operators including monotone functions of the eigenvalues of the Hessian matrix $D^{2}u$ and general convex functions of the Hessian matrix [Eva82, Proposition 5.3]. We note that the maxima can also be taken over a subset of directions by simply setting $F_{\theta}=-1$ for directions that are not active in the PDE operator. The methods described in this article can also be trivially adapted to include minima and other monotone functions of the operators $F_{\theta}$.

We remark that the PDE operator can also include Lipschitz continuous dependence on the gradient $\nabla u$. Monotone approximation of first-order operators is fairly well-established and does not require the wide-stencil structure that is necessary to correctly discretise second-order operators. In this article, we omit terms involving the gradient for the sake of compactness and clarity.

We take as a starting point the meshfree finite difference approximations developed in [Fro17b]. Given a set of discretisation points $\mathcal{G}$, that work introduced a generalised finite difference approximation of (1) of the form

[TABLE]

Above, the set of neighbours $\mathcal{N}(i,\theta)$ gives indices of several points that are near $x_{i}$ and align closely with the $e_{\theta}$ direction. In the original paper, these sets of neighbours are computed through brute force, by finding and inspecting all nodes lying within a distance $\sqrt{h}$ of each other.

In this article, we describe the construction of piecewise Cartesian meshes using a quadtree structure, which is augmented by additional discretisation points along the boundary in order to preserve consistency. We demonstrate that the resulting set of discretisation points satisfies the conditions required by [Fro17b, Theorem 13]. This ensures that it is possible to build monotone (convergent) approximations. Moreover, the structure of these quadtrees allows for easy mesh adaptation, with refinement criteria that can either be specified a priori or determined automatically from the quality of the solution of (2).

By exploiting the underyling structure of the quadtree meshes, we design an efficient method for constructing the set of neighbours $\mathcal{N}(i,\theta)$ required by our generalised finite difference stencils. This leads to a much faster numerical method for approximating solutions of (1) that can easily handle singular solutions, complicated geometries, and highly non-uniform meshes.

Finally, we describe a strategy for producing higher-order almost-monotone methods using the framework of [FO13]. This requires combining the monotone scheme with a formally higher-order approximation. By utilising the structure and symmetry within the quadtree mesh, we obtain a simple strategy for constructing higher-order schemes in the interior of the domain. Near the boundary, where these simple schemes no longer exist, we propose a least-squares approach to constructing higher-order schemes.

1.3. Contents

In section 2, we review the theory of generalised finite difference approximations for fully nonlinear elliptic equations. In section 3, we describe our strategy for constructing meshes and finite difference stencils. In section 4, we describe a higher-order implementation of these methods. In section 5, we present several computational examples. Finally, in section 6, we provide conclusions and perspectives.

2. Generalised Finite Difference Schemes

In this section we review existing results on the construction and convergence of numerical methods for fully nonlinear elliptic equations.

2.1. Weak solutions

One of the challenges associated with the approximation of fully nonlinear PDEs is the fact that classical (smooth) solutions may not exist. It thus becomes necessary to interpret PDEs using some notion of weak solution, and the numerical methods that are used need to respect this notion of weak solution. The most common concept of weak solution for this class of PDEs is the viscosity solution, which involves transferring derivatives onto smooth test functions via a maximum principle argument [CIL92].

The PDEs we consider in this work belong to the class of degenerate elliptic equations,

[TABLE]

Definition 2.1** (Degenerate elliptic).**

The operator $F:\bar{\Omega}\times\mathbb{R}\times\mathcal{S}^{2}\to\mathbb{R}$ is degenerate elliptic if

[TABLE]

whenever $u\leq v$ and $X\geq Y$.

We note that the operator is also defined on the boundary $\partial\Omega$ of the domain, which allows both the PDE and the boundary conditions to be contained within equation (3). For example, if Dirichlet data $u(x)=g(x)$ is given on the boundary, the elliptic operator at the boundary will be defined by

[TABLE]

The PDE operators (1) that we consider in this work are degenerate elliptic if they are non-decreasing functions of their second argument ($u$) and non-increasing functions of all subsequent arguments (which involve second directional derivatives).

Since degenerate elliptic equations need not have classical solutions, solutions need to be interpreted in a weak sense. The numerical methods developed in this article are guided by the very powerful concept of the viscosity solution [CIL92]. Checking the definition of the viscosity solution requires checking the value of the PDE operator for smooth test functions lying above or below the semi-continuous envelopes of the candidate solution.

Definition 2.2** (Upper and Lower Semi-Continuous Envelopes).**

The upper and lower semi-continuous envelopes of a function $u(x)$ are defined, respectively, by

[TABLE]

[TABLE]

Definition 2.3** (Viscosity subsolution (supersolution)).**

An upper (lower) semi-continuous function $u$ is a viscosity subsolution (supersolution) of (1) if for every $\phi\in C^{2}(\bar{\Omega})$, whenever $u-\phi$ has a local maximum (minimum) at $x\in\bar{\Omega}$, then

[TABLE]

Definition 2.4** (Viscosity solution).**

A function $u$ is a viscosity solution of (1) if $u^{*}$ is a subsolution and $u_{*}$ a supersolution.

An important property of many elliptic equations is the comparison principle, which immediately implies uniqueness of the solution.

Definition 2.5** (Comparison principle).**

A PDE has a comparison principle if whenever $u$ is an upper semi-continuous subsolution and $v$ a lower semi-continuous supersolution of the equation, then $u\leq v$ on $\bar{\Omega}$.

2.2. Convergence of elliptic schemes

In order to construct convergent approximations of elliptic operators, we will rely on the framework provided by Barles and Souganidis [BS91] and further developed by Oberman [Obe06].

We consider finite difference schemes that have the form

[TABLE]

where $\epsilon$ is a small parameter.

The convergence framework requires notions of consistency and monotonicity, which we define below.

Definition 2.6** (Consistency).**

The scheme (4) is consistent with the equation (1) if for any smooth function $\phi$ and $x\in\bar{\Omega}$,

[TABLE]

[TABLE]

Definition 2.7** (Monotonicity).**

The scheme (4) is monotone if $F^{\epsilon}$ is a non-decreasing function of its final two arguments.

Schemes that satisfy these two properties respect the notion of the viscosity solution at the discrete level. In particular, these schemes preserve the maximum principle and are guaranteed to converge to the solution of the underlying PDE.

Theorem 2.8** (Convergence [Obe06]).**

Let $u$ be the unique viscosity solution of the PDE (1), where $F$ is a degenerate elliptic operator with a comparison principle. Let the finite difference approximation $F^{\epsilon}$ be consistent and monotone and let $u^{\epsilon}$ be any solution of the scheme (4), with bounds independent of $\epsilon$. Then $u^{\epsilon}$ converges uniformly to $u$ as $\epsilon\to 0$.

The above theorem assumes existence of a bounded solution to the approximation scheme. This is typically straightforward to show for a consistent, monotone approximation of a well-posed PDE, though the precise details can vary slightly and rely on available well-posedness theory for the PDE in question.

Theorem 2.9** (Existence and Stability [Fro17a, Lemmas 35-36]).**

Let $F^{\epsilon}$ be a consistent, monotone scheme that is Lipschitz in its last two arguments. Suppose also that there exist strict classical sub- and super-solutions to the PDE (3). Then for small enough $\epsilon>0$, the scheme (4) has a solution $u^{\epsilon}$. Moreover, there exists a constant $M>0$ such that $\|u^{\epsilon}\|_{\infty}\leq M$ for sufficiently small $\epsilon>0$.

In many cases, simple quadratic functions will serve as the sub- and super-solutions required by Theorem 2.9. For more complicated PDE operators, particularly those with a non-trivial dependence on the gradient $\nabla u$, the theory of classical solutions of the equation can often be used to show the existence of these sub- and super-solutions.

2.3. Meshfree finite difference approximations

In [Fro17b], a new generalised finite difference method was introduced for approximating fully nonlinear second order elliptic operators on point clouds. We review the key results of that work, which will be foundational to the higher-order adaptive methods that will be developed in the remainder of this article.

Definition 2.10** (Notation).**

1. (N1)

$\Omega\subset\mathbb{R}^{2}$ is a bounded domain with Lipschitz boundary $\partial\Omega$. 2. (N2)

$\mathcal{G}\subset\bar{\Omega}$ is a point cloud consisting of the points $x_{i}$, $i=1,\ldots,N$. 3. (N3)

$h=\sup\limits_{x\in{\Omega}}\min\limits_{y\in\mathcal{G}}|x-y|$ is the spatial resolution of the point cloud. In particular, every ball of radius $h$ contained in $\bar{\Omega}$ contains at least one discretisation point $x_{i}$. 4. (N4)

$h_{B}=\sup\limits_{x\in{\partial\Omega}}\min\limits_{y\in\mathcal{G}\cap\partial\Omega}|x-y|$ is the resolution of the point cloud on the boundary. In particular, every ball of radius $h_{B}$ centred at a boundary point $x\in\partial\Omega$ contains at least one discretisation point $x_{i}\in\mathcal{G}\cap\partial\Omega$ on the boundary. 5. (N5)

$\delta=\min\limits_{x\in\Omega\cap\mathcal{G}}\inf\limits_{y\in\partial\Omega}|x-y|$ is the distance between the set of interior discretisation points and the boundary. In particular, if $x_{i}\in\mathcal{G}\cap\Omega$ and $x_{j}\in\partial\Omega$, then the distance between $x_{i}$ and $x_{j}$ is at least $\delta$. 6. (N6)

$d\phi$ is the angular resolution used to approximate the second directional derivatives $u_{\theta\theta}$. 7. (N7)

$d\theta$ is the angular resolution used to approximate the nonlinear operator. 8. (N8)

$\epsilon$ is the search radius associated with the point cloud.

Discretising the PDE requires approximating second directional derivatives $u_{\theta\theta}$ at each interior discretisation point $x_{i}\in\mathcal{G}$. To accomplish this, we consider all points $x_{j}\in\mathcal{G}\cap B(x_{i},\epsilon)$ within a search neighbourhood of radius $\epsilon$ centred at $x_{i}$. Discretisation points within this neighbourhood can be written in polar coordinates $(r,\phi)$ with respect to the axes defined by the lines $x_{0}+t(\cos\theta,\sin\theta)$, $x_{0}+t(-\sin\theta,\cos\theta)$. We seek one neighbouring discretisation point in each quadrant described by these axes, with each neighbour aligning as closely as possible with the line $x_{0}+t\nu$, where $\nu=(\cos\theta,\sin\theta)$. That is, we select the neighbours

[TABLE]

for $j=1,\ldots,4$. See Figure 1. We say that a stencil with angular resolution $d\phi$ exists for the point cloud $\mathcal{G}$ if for all interior discretisation points, the four discretisation points $x_{j}\in\mathcal{G}$ defined by (5) exist and satisfy $d\phi=\max\{\phi_{j}\}$.

Because of the “wide-stencil” nature of these approximations (since the search radius $\epsilon\gg h$), care must be taken near the boundary. In order to preserve consistency up to the boundary, it is necessary that the boundary be more highly resolved than the interior ($h_{B}\ll h$). In particular, this means that a simple Cartesian mesh (or piecewise Cartesian mesh) is not sufficient for producing consistent schemes up to the boundary.

Then a consistent, monotone approximation of $u_{\theta\theta}$ is

[TABLE]

where we use the polar coordinate characterisation of the neighbours to define

[TABLE]

and the coefficients are given by

[TABLE]

In general, the PDE requires evaluating second directional derivatives in all possible directions. Instead, we consider a finite subset $\mathcal{A}=\left\{jd\theta\mid j=0,\ldots,\lfloor\frac{2\pi}{d\theta}\rfloor\right\}\subset[0,2\pi)$ with a resolution $d\theta$.

Then we can substitute these coefficients into (2) to obtain the scheme:

[TABLE]

We recall the convergence result from [Fro17b, Theorem 18].

Theorem 2.11** (Convergence).**

Let $F$ be a degenerate elliptic operator with a comparison principle that is Lipschitz continuous in $u_{\theta\theta}$ for each $\theta\in[0,2\pi)$ and let $u$ be the unique viscosity solution of the PDE (1). Suppose also that (1) has a strict classical sub- and super-solution. Consider a sequence of point clouds $\mathcal{G}^{n}$, with parameters defined as in Definition 2.10, which satisfy the following conditions.

• •

The spatial resolution $h^{n}\to 0$ as $n\to\infty$.

• •

The boundary resolution satisfies $h_{B}^{n}/\delta^{n}\to 0$ as $n\to\infty$.

• •

The search radius satisfies both $\epsilon^{n}\to 0$ and $h^{n}/\epsilon^{n}\to 0$ as $n\to\infty$.

• •

The angular resolution $d\theta^{n}\to 0$ as $h^{n}\to 0$.

Then for sufficiently large $n$, the approximation scheme (6) admits a solution $u^{n}$ and $u^{n}$ converges uniformly to $u$ as $n\to\infty$.

We note that the angular resolution that emerges from the scheme satisfies $d\phi=\mathcal{O}(\max\{h/\epsilon,h_{B}/\delta\})$ (Figure 2). For a uniform grid, a natural choice of parameters is $\epsilon=\mathcal{O}(\sqrt{h})$, $h_{B}=\mathcal{O}(h^{3/2})$, $\delta=\mathcal{O}(h)$, $d\theta=\mathcal{O}(\sqrt{h})$. This leads to a formally optimal discretisation error of $\mathcal{O}(\sqrt{h})$.

We remark also that these parameters can be defined locally instead of globally in order to accommodate highly non-uniform meshes.

3. Construction of Meshes and Stencils

In this section, we explain how we use augmented quadtrees [dBCvKO08, Chapter 14] to build piecewise Cartesian meshes with additional discretisation points on the boundary. We organise this section as follows. In subsection 3.1, we recall the basic structure of a quadtree. In subsection 3.2 we explain how we augment the quadtree to deal with the boundary. In subsection 3.4, we explain how the quadtrees are used to efficiently find the stencils. Finally, in subsection 3.3 we discuss mesh adaptation.

3.1. Quadtrees

Quadtrees are based on a simple idea: a square can be divided into four smaller squares which correspond to the four quadrants of the square. A quadtree is then a rooted tree in which every internal node has four children and every node in the tree corresponds to a square. A square with no children is called a leaf square. See Figure 3.

Quadtrees can then easily be used to build uniform and non-uniform meshes: the squares’ vertices are the mesh points. This structure is appealing because it is general enough to allow for local mesh adaptation, while still maintaining enough structure to efficiently build the finite difference stencils. Indeed, as we will see in subsection 3.4, the quadtree structure allow us to significantly reduce the number of mesh points inspected when constructing our stencils. However, the quadtree in and of itself is not ideal for handling complicated geometries as the mesh points are restricted to be vertices of the squares. We observe that the global spatial resolution $h$ of a quadtree corresponds to the length-scale of the largest leaf square. However, the quadtree can be highly non-uniform and the local spatial resolution near a particular point may be much less than $h$.

3.2. Meshing the boundary

Quadtrees alone are not enough for the schemes proposed here: the boundary requires additional treatment. As discussed in section 2, the boundary must be more highly resolved than the interior to maintain consistency of the numerical method. As a result, we cannot restrict the mesh points to be the vertices of the squares in the quadtree.

To overcome this, we build augmented quadtrees: each leaf square that intersects the boundary is marked as such and additional mesh points that lie on the boundary are added and associated with this boundary leaf square. The immediate advantage of this approach is that mesh points may lie exactly on the boundary, which allows us to handle complicated geometries with ease. In addition, by keeping track of which boundary leaf square the mesh points belong to, we preserve one of the key properties of the quadtree: knowledge of the relative position of the mesh points. This allows for efficient construction of the finite difference stencils.

We make the following general assumption: each edge of a leaf square intersects the boundary at most once. This is a reasonable assumption that simply entails that our quadtree must be sufficiently refined near the boundary. See Figure 4, where each edge of the grey squares intersects the boundary at most once.

The only question left to address is exactly how many additional boundary mesh points one must add to guarantee the existence of a consistent stencil. This is addressed in Theorem 2.11, which requires that the boundary resolution go to zero more quickly than the resolution of the “standard” quadtree, $h_{B}=o(h)$, and more quickly than the gap between the boundary and the interior nodes, $h_{B}=o(\delta)$. Note that these conditions need only be satisfied locally rather than globally.

We define a simple algorithm that enlarges a given point cloud so that the condition $h_{B}\leq 2\delta\tan(d\theta/2)$ is satisfied locally (see Algorithm 1). This ensures that the angular resolution of the finite difference approximations is commensurate with the angular resolution used to approximate the nonlinear operator: $d\phi\lesssim 2\tan^{-1}\left(\dfrac{h_{B}}{2\delta}\right)\leq d\theta$ (Figure 2).

In Figure 4, the meshes obtained by applying Algorithm 1 to the point cloud obtained from a quadtree of depth $4$ are displayed. The fan-shaped domain illustrates the advantages of the local criteria: the boundary is only highly resolved when there are interior mesh points nearby.

These augmented quadtrees enable us to construct convergent (consistent and monotone) finite difference approximations.

Lemma 3.1** (Approximation with augmented quadtrees).**

Consider a sequence of augmented quadtrees $\mathcal{G}^{n}$ constructed via Algorithm 1 with spatial resolution $h_{n}\to 0$. Consider also a sequence of search radii $\epsilon^{n}=\mathcal{O}(\sqrt{h^{n}})$ and angular resolutions $d\theta^{n}=\mathcal{O}(\sqrt{h^{n}})$. Then $\mathcal{G}^{n}$ satisfies the hypotheses of Theorem 2.11.

Proof.

We need only verify that $h_{B}^{n}/\delta^{n}\to 0$; the remaining conditions of Theorem 2.11 are trivially satisfied.

We recall that both $h_{B}^{n}$ and $\delta^{n}$ can be defined locally. Indeed, for each boundary leaf square $S$ we can let

[TABLE]

By construction, Algorithm 1 ensures that $h_{B,S}^{n}/\delta_{S}^{n}=\mathcal{O}(d\theta_{S}^{n})\to 0$.

Moreover, it is sufficient to verify these conditions at boundary leaf squares; other interior squares will produce larger values of $\delta^{n}_{S}$ and smaller values of $h_{B,S}^{n}/\delta^{n}_{S}$. ∎

3.3. Refinement, adaptivity and balance

The use of quadtrees also provides a natural means of doing mesh adaptation. A refinement criteria can either be specified a priori or determined automatically from the quality of the solution.

In Figure 5, we provide an example of a priori refinement: the mesh is refined near the corners of the domain.

Simply refining the quadtree can lead to a very unbalanced quadtree when large squares adjoin several smaller squares. This is an undesirable property for our meshes as it makes the construction of high-order schemes significantly more difficult. Therefore, we always maintain a balanced quadtree: any two neighbouring squares differ by at most a factor of two in length scale (see Figure 5). Balancing a quadtree can be done efficiently; we refer to [dBCvKO08, Theorem 14.4] for details.

3.4. Generating the stencil

We explain how quadtrees are used to efficiently find the neighbours for each interior mesh point. The main idea is the following: given the quadtree structure we know the relative position of the mesh points and can significantly restrict the number of nodes we examine.

Consider a direction $\nu=(\cos\theta,\sin\theta)$ and the line $x_{0}+t\nu$. Without loss of generality, assume the line has positive slope as in Figure 6. We describe the procedure for finding the required neighbours of $x_{0}$ lying in the first and fourth quadrants.

From the list of potential neighbours in each quadrant, the precise neighbours used in the stencil are determined via (5). See Figure 6, for a close-up of the neighbours search in a uniform (left) and non-uniform (right) mesh.

Referring to Figure 6 (a), we provide a rough estimate on the improvement this algorithm yields for a uniform $N\times N$ grid with grid spacing $h=\mathcal{O}(1/N)$. Recall that the search region is a disc of radius $\epsilon$. The brute force algorithm used in [Fro17b] examines $\mathcal{O}((\epsilon/h)^{2})$ neighbours, while the algorithm proposed above using quadtrees examines only $\mathcal{O}(\epsilon/h)$ neighbours. Given the typical choice $\epsilon=\mathcal{O}(\sqrt{h})$, the cost of constructing the stencil at each point is reduced from $\mathcal{O}(N)$ to $\mathcal{O}(\sqrt{N})$. A similar speed-up is seen for the piecewise Cartesian meshes produced by the quadtree.

4. Higher-Order Methods

We also introduce a technique for building formally higher-order approximations on highly non-uniform/unstructured grids. We focus on second-order schemes, which is typically sufficient for applications, but the same ideas can be easily be extended to higher-order schemes.

4.1. Filtered schemes

The meshfree finite difference approximation discussed in section 2 is low accuracy; formally it is at best $\mathcal{O}(\sqrt{h})$. However it can be used as the foundation for higher-order convergent filtered schemes as in [FO13]. The main idea is to blend a monotone convergent scheme with a non-monotone accurate scheme and retain the advantages of both: stability and convergence of the former, and higher accuracy of the latter.

To accomplish this, we let $F_{A}[u]$ be any higher-order scheme and $F_{M}[u]$ be a monotone approximation scheme, both defined on the same mesh. The filtered scheme is then defined as

[TABLE]

where the filter S is given by

[TABLE]

As long as $\alpha>0$, this approximation converges to the viscosity solution of the PDE under the same conditions as the monotone scheme converges. The underlying reason is that this scheme is a small perturbation of a monotone scheme and the proof of the Barles-Souganidis theorem is easily modified to accommodate this. Moreover, if $h^{\alpha}$ is larger than the discretisation error of the monotone scheme, the formal accuracy of the filtered scheme is the same as the formal accuracy of the non-monotone scheme.

4.2. Higher-order schemes in interior

We discuss how to build high-order schemes for the non-uniform meshes proposed in section 3. In this section, we focus on interior mesh points away from the boundary.

Defining higher-order schemes for (1) reduces to defining higher-order approximations to $u_{xx}$, $u_{yy}$ and $u_{xy}$. We will focus on building second order approximations, although the ideas are easily generalised. For a uniform Cartesian grid, such as in Figure 7, these are widely known and are given by

[TABLE]

In a uniform cartesian grid, all the nodes are regular nodes; i.e., each node is the vertex of four different squares like the one depicted in Figure 7. However, the meshes proposed here are non-uniform and in general not all nodes will be regular nodes. We need also consider dangling nodes, which occur midway along the shared edge of two equally-sized squares, one of which is subdivided. Thus additional work is required to define the higher-order schemes.

As explained in section 3, the meshes are generated using quadtrees that are kept balanced (the lengths of neighbouring squares differ by at most a factor of two). Thus each interior mesh point can be associated to one of five different configurations. These are depicted in Figures 7, 9, and 9. The generic element chosen to represent each configuration is one where each square is a leaf square. In general, one or more of the smaller squares may have children in the quadtree; i.e., they may be subdivided into smaller squares. We consider these to be redundant when constructing the high-order schemes. Considering all possible different configurations would only increase the complexity of the schemes with no additional benefits as the schemes would remain asymptotically second order; only the asymptotic error constant could be improved.

For the configurations in Figure 9, we use the following approximations

[TABLE]

As for the configurations in Figure 9, we have

[TABLE]

The standard Taylor expansion argument shows that the above expressions are second order accurate.

Finally, we explain how one can efficiently determine the configuration of each interior mesh point. As mentioned in subsection 3.4, for each interior mesh point a record is kept of the four leaf squares that the interior mesh point as a vertex. Thus the configuration is easily determined by determining the depth of the neighbouring squares and respective parent squares in the quadtree.

4.3. Least-squares constructions near boundary

In this section, we discuss how to construct second order schemes at interior points near the boundary. When near the boundary, the construction of the schemes cannot reduce to the cases considered in the previous section: in general, not all the neighbouring mesh points will be the vertices of squares and some will lie on the boundary, which we allow to have a complicated geometry. Thus additional care is needed. Here we describe a general strategy for building high-order schemes.

Let $\{x_{i}\}_{i=1}^{N}$ denote neighbouring mesh points to the interior mesh point $x_{0}$ with $\left\|x_{i}-x_{0}\right\|_{\infty}=\mathcal{O}(h)$. Using Taylor expansion we obtain, for each $i=1,\ldots,N$,

[TABLE]

where we are using the multi-index notation. Hence

[TABLE]

if the $\{a_{i}\}_{i=1}^{N}$ solve the linear system

[TABLE]

for $0<|\alpha|\leq 3$.

To approximate second derivatives to second order, we expect to require $N=9$ neighbours.

Designing second order schemes is now reduced to determining the neighbouring mesh points and solving the respective linear system. However, since we are particularly interested in the case where some of the neighbouring mesh points lie on the boundary, which may have a complicated geometry, it is hard to make any a priori claim regarding the invertability and conditioning of the linear system. It is important to point out that we are interested in obtaining any particular solution. As we saw in the previous section, depending on the derivative being approximated and the location of the neighbouring mesh points, the number of neighbouring mesh points required changes.

We are now ready to describe the strategy implemented to construct the higher-order schemes. First, we determine which configuration we are in. If all the vertices of squares neighbouring $x_{0}$ are mesh points, we use the approximations described in the previous section. Otherwise, we build the linear system above using all mesh points $x_{i}$ that lie within squares adjoining $x_{0}$ (some will lie on the boundary of the domain and will not be vertices of the squares). In general, we will have $N\geq 9$, but the linear system may still have no solution or be ill-conditioned. If that is the case, we consider additional neighbouring mesh points by adding the mesh points that belong to neighbouring squares. In general, we end up with an under-determined system and we select the least squares solution. In practice, adding additional neighbouring mesh points was not always required, and it was never required more than once. Thus the high-order scheme has stencil width $\mathcal{O}(h)$ and preserves the formal discretisation error of $\mathcal{O}(h^{2})$.

5. Computational Examples

5.1. Monge-Ampère equation

We consider the Monge-Ampère equation

[TABLE]

The PDE is only elliptic in the space of convex functions. However, as in [Fro12], we can use the globally elliptic extension

[TABLE]

We will consider four different domains given by $\Omega=\{(x,y)\in\mathbb{R}^{2}\mid\phi(x,y)<0\}$ where $\phi$ is given by

1. (a)

(circle) $\phi(x,y)=x^{2}+y^{2}-1$, 2. (b)

(ellipse) $\phi(x,y)=x^{2}+2y^{2}-1$, 3. (c)

(diamond) $\phi(x,y)=|x|+|y|-1$, 4. (d)

(diamond stretched) $\phi(x,y)=|x|+|2y|-1$.

Example 5.1**.**

We consider first the following $C^{2}$ solution of the Monge-Ampère equation

[TABLE]

Results are displayed in Table 1. On each domain, the filtered implementation recovers the desired second-order accuracy even though the boundary nodes do not belong to the structured piecewise Cartesian mesh.

Example 5.2**.**

We consider also a $C^{1}$ solution of the Monge-Ampère equation, for which the ellipticity is degenerate in an open set.

[TABLE]

This solution is singular, and there is thus no realistic hope of attaining the formal second-order discretisation error obtained from Taylor’s Theorem. Nevertheless, the method converges and we observe roughly first-order accuracy (Table 2).

5.2. Computation time

One of the reasons to use quadtrees to build non-uniform meshes was to efficiently find the neighbours for the finite difference schemes. Here we compare the implementation with quadtrees, discussed in depth in section 3, with the simple but inefficient brute force approach of [Fro17b]. Reduced CPU time is the ultimate goal, but not necessarily a robust measure of efficiency as it depends on both hardware and software implementations. Here both implementations were vectorised whenever possible in order to optimise the code for MATLAB. Although it is still possible to further optimise the code, the CPU time should be a fair indication of the improvement gained by using quadtrees.

In Figure 10, we compare the number of mesh points versus the CPU time required to generate the mesh and find the stencil. We let the domain $\Omega$ be the ellipse given by $\Omega=\{(x,y)\in\mathbb{R}^{2}\mid x^{2}+2y^{2}<1\}$. The results indicate that our new approach is optimal, with the computation time required to construct the stencils being roughly proportional to the number of mesh points. This represents an improvement of roughly one order over the original brute force approach. In practice, we can generate a mesh and find the stencil based on a uniform $256\times 256$ grid in roughly $50$ seconds with the use of quadtrees, instead of over $6$ minutes with the previous brute force approach.

5.3. Adaptivity

In our next example, we demonstrate the improvements possible with adaptivity using our generalised finite difference approximations.

Example 5.3**.**

We consider the fully nonlinear convex envelope equation

[TABLE]

where

[TABLE]

The equation is posed on an ellipse with semi-major axis equal to one and semi-minor axis equal to one-half, which is rotated through an angle of $\phi=\pi/6$. The obstacle $g$ consists of two cones,

[TABLE]

and the exact solution is

[TABLE]

We note that this solution is only Lipschitz continuous, and the equation must be understood in a weak sense.

We defined the following refinement strategy. Given a solution a solution $u^{h}$,

1. (1)

Compute numerically $\left\|D^{2}u\right\|=\sqrt{u_{xx}^{2}+2u_{xy}^{2}+u_{yy}^{2}}$. 2. (2)

Refine twice every square with a vertex such that $h\left\|D^{2}u\right\|>0.5$.

This refinement strategy is inspired by the hybrid scheme proposed in [Obe08b].

A convergence plot is displayed in Figure 11. We observe an improvement in accuracy from roughly $\mathcal{O}(N^{-0.38})$ to $\mathcal{O}(N^{-0.64})$. Moreover, the use of adaptivity leads to a clear qualitative improvement in the computed solutions. This is evident in Figure 12, which shows that the solution obtained with the uniform mesh is visibly non-convex along the singularity, which does not align with the grid. By resolving this singularity, the adaptive method produces a solution that has a dramatically better quality.

6. Conclusions

In this article, we described generalised finite difference methods for solving a large class of fully nonlinear elliptic partial differential equations. These methods were inspired by the meshfree methods described in [Fro17b], which are flexible and convergent, but very expensive to implement.

Our meshes used a modified quadtree structure that used piecewise Cartesian grids in the interior of the domain, augmented by a set of points lying exactly on the boundary. The inclusion of the additional boundary points was needed to ensure convergence of the numerical methods. This type of mesh also allows us to deal easily with complicated boundary geometries.

By relying on the underlying quadtree structure, we developed an algorithm for efficiently constructing the mesh and finite difference stencils. This led to a dramatic improvement in computational efficiency as compared to the original brute force approach.

We also described a strategy for constructing higher-order approximations, which still fit within the convergence framework. In the interior of the domain, the quadtree structure allows us to explicitly write out the higher-order finite difference schemes. This strategy fails near the boundary, which can have a very complicated geometry. At these points, we employed a least-squares approach to construct higher-order schemes.

We also used these methods to perform automatic mesh adaptation, which refines the mesh near singularities in the computed solutions. This led to a dramatic qualitative and quantitative improvement in computed solutions.

Acknowledgments

We thank Adam Oberman for helpful discussions and support of this project.