Stability and guaranteed error control of approximations to the Monge–Ampère equation

Abstract

This paper analyzes a regularization scheme of the Monge–Ampère equation by uniformly elliptic Hamilton–Jacobi–Bellman equations. The main tools are stability estimates in the $L^\infty$ norm from the theory of viscosity solutions which are independent of the regularization parameter $\varepsilon$. They allow for the uniform convergence of the solution $u_\varepsilon$ to the regularized problem towards the Alexandrov solution u to the Monge–Ampère equation for any nonnegative $L^n$ right-hand side and continuous Dirichlet data. The main application are guaranteed a posteriori error bounds in the $L^\infty$ norm for continuously differentiable finite element approximations of u or $u_\varepsilon$.

1 Introduction

1.1 Overview

Let $\varOmega \subset {\mathbb {R}}^n$, $n\ge 2$, be a bounded and convex domain. Given a nonnegative function $0 \le f \in L^n(\varOmega )$ and continuous Dirichlet data $g \in C(\partial \varOmega )$, the Monge–Ampère equation seeks the unique (convex) Alexandrov solution $u \in C(\overline{\varOmega })$ to

$$\begin{aligned} \det \textrm{D}^2 u = (f/n)^n \text { in } \varOmega \quad \text {and}\quad u = g \text { on } \partial \varOmega . \end{aligned}$$

(1)

If the Dirichlet data $g \ne 0$ is non-homogenous, then we additionally assume that $\varOmega$ is strictly convex. The re-scaling ${\tilde{f}} :=(f/n)^n$ of the right-hand side is not essential, but turns out convenient for purposes of notation. By the Alexandrov solution u to (1) we mean a convex function $u \in C(\varOmega )$ with $u=g$ on $\partial \varOmega$ and

$$\begin{aligned} {\mathcal {L}}^n(\partial v(\omega )) = \int _\omega {\tilde{f}}\,\textrm{d}x \quad \text {for any Borel subset }\omega \subset \varOmega . \end{aligned}$$

The left-hand side denotes the Monge–Ampère measure of $\omega$, i.e., the n-dimensional Lebesgue measure of all vectors in the subdifferential $\partial v(\omega ) :=\cup _{x \in \omega } \partial v(x)$ where $\partial v(x)$ is the usual subdifferential of v in a point x. We remark that this solution concept admits more general right-hand sides (Borel measures), which are, however, disregarded in this work. For further details, we refer to the monographs [14, 17]. It is known [1] that the Alexandrov solution to (1) exists and is unique. In addition, it was shown [4] that if $f \in C^{0,\alpha }(\varOmega )$, $0 < \lambda \le f \le \varLambda$, and $g \in C^{1,\beta }(\partial \varOmega )$ with positive constants $0< \alpha ,\beta < 1$ and $0 < \lambda \le \varLambda$, then $u \in C(\overline{\varOmega }) \cap C^{2,\alpha }_{\text {loc}}(\varOmega )$.

It is known [12, 19] that (1) can be equivalently formulated as a Hamilton–Jacobi–Bellman (HJB) equation, a property that turned out useful for the numerical solution of (1) [12, 15]; one of the reasons being that the latter is elliptic on the whole space of symmetric matrices ${\mathbb {S}} \subset {\mathbb {R}}^{n \times n}$ and, therefore, the convexity constraint is automatically enforced by the HJB formulation. For nonnegative continuous right-hand sides $0 \le f \in C(\varOmega )$, the Monge–Ampère equation (1) is equivalent to

$$\begin{aligned} F_0(f;x,\textrm{D}^2 u) = 0 \text { in } \varOmega \quad \text {and}\quad u = g \text { on } \partial \varOmega \end{aligned}$$

with $F_0(f;x,M) :=\sup _{A \in {\mathbb {S}}(0)}(-A:M + f\root n \of {\det A})$ for any $x \in \varOmega$ and $M \in {\mathbb {R}}^{n\times n}$. Here, ${\mathbb {S}}(0) :=\{A \in {\mathbb {S}}: A \ge 0 \text { and } \textrm{tr}\,A = 1\}$ denotes the set of positive semidefinite symmetric matrices A with unit trace $\textrm{tr}\,A = 1$. Since $F_0$ is only degenerate elliptic, the regularization scheme proposed in [15] replaces ${\mathbb {S}}(0)$ by a compact subset ${\mathbb {S}}(\varepsilon ) :=\{A \in {\mathbb {S}}(0): A \ge \varepsilon \textrm{I}_{n \times n}\} \subset {\mathbb {S}}(0)$ of matrices with eigenvalues bounded from below by the regularization parameter $0 < \varepsilon \le 1/n$. (We note that $\varepsilon > 1/n$ would lead to an empty set ${\mathbb {S}}(\varepsilon )$ because of the unit trace condition.) The solution $u_\varepsilon$ to the regularized PDE solves

$$\begin{aligned} F_\varepsilon (f;x,\textrm{D}^2 u_\varepsilon ) = 0 \text { in } \varOmega \quad \text {and}\quad u_\varepsilon = g \text { on } \partial \varOmega \end{aligned}$$

(2)

where, for any $x \in \varOmega$ and $M \in {\mathbb {R}}^{n \times n}$, the function $F_\varepsilon$ is defined as

$$\begin{aligned} F_\varepsilon (f;x,M) :=\sup \nolimits _{A \in {\mathbb {S}}(\varepsilon )} (-A:M + f\root n \of {\det A}). \end{aligned}$$

(3)

In two space dimensions $n = 2$, uniformly elliptic HJB equations satisfy the Cordes condition [20] and this allows for a variational setting for (2) with a unique strong solution $u_\varepsilon \in H^2(\varOmega )$ in the sense that $F_\varepsilon (f;x,\textrm{D}^2 u_\varepsilon ) = 0$ holds a.e. in $\varOmega$ [25, 26]. The paper [15] establishes uniform convergence of $u_\varepsilon$ towards the generalized solution u to the Monge–Ampère equation (1) as $\varepsilon \searrow 0$ under the assumption $g \in H^2(\varOmega ) \cap C^{1,\alpha }(\overline{\varOmega })$ and that $0 \le f \in L^2(\varOmega )$ can be approximated from below by a pointwise monotone sequence of positive continuous functions.

1.2 Contributions of this paper

The variational setting of (2) in two space dimensions leads to $H^2$ stability estimates that deteriorate with $\varepsilon ^{-1} \rightarrow \infty$ as the regularization parameter $\varepsilon \rightarrow 0$ vanishes. This can be explained by the regularity of Alexandrov solutions to the Monge–Ampère equation (1) as they are, in general, not in $H^2(\varOmega )$ without additional assumptions on the domain $\varOmega$ and the data f, g. Consequently, error estimates in the $H^2$ norm may not be of interest, and the focus is on error estimates in the $L^\infty$ norm.

The analysis departs from the following $L^\infty$ stability estimate that arises from the Alexandrov maximum principle. If $v_1,v_2 \in C(\overline{\varOmega })$ are viscosity solutions to $F_\varepsilon (f_j; x, \textrm{D}^2 v_j) = 0$ in $\varOmega$ with $0 \le \varepsilon \le 1/n$ and $f_1,f_2 \in C(\overline{\varOmega })$, then

$$\begin{aligned} \Vert v_1-v_2\Vert _{L^\infty (\varOmega )} \le \Vert v_1-v_2\Vert _{L^\infty (\partial \varOmega )} + C(n,\textrm{diam}(\varOmega ))\Vert f_1-f_2\Vert _{L^n(\varOmega )}. \end{aligned}$$

(4)

(For $\varepsilon = 0$, we additionally assume $f_1,f_2 \ge 0$.) The constant $C(n,\textrm{diam}(\varOmega ))$ exclusively depends on the dimension n and the diameter $\textrm{diam}(\varOmega )$ of $\varOmega$, but not on the ellipticity constant of (2) or on the regularization parameter $\varepsilon$. Consequently, this allows for control of the $L^\infty$ error even as $\varepsilon \rightarrow 0$. By density of $C(\overline{\varOmega })$ in $L^n(\varOmega )$, the $L^\infty$ stability estimate (4) leads to a generalized viscosity solution concept with $L^n$ right-hand sides, which coincides with $L^n$ viscosity solutions from [3, 18]. Furthermore, (4) also holds for these solutions as well with two applications highlighted below. First, this paper establishes, in extension to [15], uniform convergence of (generalized) viscosity solutions $u_\varepsilon$ of the regularized PDE (2) to the Alexandrov solution $u \in C(\overline{\varOmega })$ of the Monge–Ampère equation (2) under the (essentially) minimal assumptions $0 \le f \in L^n(\varOmega )$ and $g \in C(\partial \varOmega )$ on the data. Second, (4) provides guaranteed error control in the $L^\infty$ norm (even for inexact solve) for $H^2$ conforming FEM. The observation that continuous dependence on the data can lead to a posteriori error bounds was used in other contexts as well, see, e.g., [8, 13].

1.3 Outline

The principal tool we use for establishing our results is the celebrated Alexandrov maximum principle. It provides an upper bound for the $L^\infty$ norm of any convex function in dependence of its Monge–Ampère measure.

Lemma 1

(Alexandrov maximum principle) There exists a constant $c_n$ solely depending on the dimension n such that any convex function $v \in C(\overline{\varOmega })$ with homogenous boundary data $v|_{\partial \varOmega } = 0$ over an open bounded convex domain $\varOmega$ satisfies

$$\begin{aligned} |v(x)|^n \le c_n^n\textrm{dist}(x,\partial \varOmega ) \textrm{diam}(\varOmega )^{n-1}{\mathcal {L}}^n(\partial v(\varOmega )) \quad \text {for any } x \in \varOmega . \end{aligned}$$

(5)

Proof

This is [14, Theorem 2.8] and the constant $c_n :=(2(2\pi )^{n/2-1}/((n-1)!!n)$ arises therein from the n-dimensional volume formula for a cone ${\mathcal {C}} \subset \partial v(\varOmega )$. If $n = 2$, then $c_2 = 1$. We note that, under additional assumptions, the scaling with respect to $\textrm{dist}(x,\partial \varOmega )$ in (5) can be improved at the cost of control over the involved constants, cf. [5] for $n \ge 3$ and [16] for $n = 2$. $\square$

The remaining parts of this paper are organized as follows. Section 2 establishes $L^\infty$ stability estimates for viscosity solutions to the HJB equation (2) for all parameters $0 \le \varepsilon \le 1/n$ in any space dimension. Section 3 provides a proof of convergence of the regularization scheme. A posteriori error estimates for the discretization error in the $L^\infty$ norm for $H^2$-conforming FEM are presented in Sect. 4. The three numerical experiments in Sect. 5 conclude this paper.

Standard notation for function spaces applies throughout this paper. Let $C^k(\varOmega )$ for $k \in {\mathbb {N}}$ denote the space of scalar-valued k-times continuously differentiable functions. Given a positive parameter $0 < \alpha \le 1$, the Hölder space $C^{k,\alpha }(\varOmega )$ is the subspace of $C^k(\overline{\varOmega })$ such that all partial derivates of order k are Hölder continuous with exponent $\alpha$. For any set $\omega \subset {\mathbb {R}}^n$, $\chi _\omega$ denotes the indicator function associated with $\omega$. For $A,B \in {\mathbb {R}}^{n \times n}$, the Euclidean scalar product $A:B :=\sum _{j,k = 1}^{n} A_{jk} B_{jk}$ induces the Frobenius norm $|A| :=\sqrt{A:A}$ in ${\mathbb {R}}^{n \times n}$. The notation $|\cdot |$ also denotes the absolute value of a scalar or the length of a vector. The relation $A \le B$ of symmetric matrices $A,B \in {\mathbb {S}}$ holds whenever $B - A$ is positive semidefinite. The notation $a \lesssim b$ abbreviates $a \le Cb$ for a generic constant C independent of the mesh-size as well as the regularization parameter $\varepsilon$ and $a \approx b$ abbreviates $a \lesssim b \lesssim a$.

2 Stability estimate

We first recall the concept of viscosity solutions to the HJB equation (2).

Definition 1

(Viscosity solution) Given $f \in C(\varOmega )$ and $0 \le \varepsilon \le 1/n$, a function $v \in C(\overline{\varOmega })$ is a viscosity subsolution (resp. supersolution) to $F_\varepsilon (f;x,\textrm{D}^2 v) = 0$ if, for all $x_0 \in \varOmega$ and $\varphi \in C^2(\varOmega )$ such that $v - \varphi$ has a local maximum (resp. minimum) at $x_0$, $F_\varepsilon (f;x_0,\textrm{D}^2 \varphi (x_0)) \le 0$ (resp. $F_\varepsilon (f;x_0,\textrm{D}^2 \varphi (x_0)) \ge 0$). If v is viscosity sub- and supersolution, then v is called viscosity solution to $F_\varepsilon (f;x,\textrm{D}^2 v) = 0$.

The following result provides the first tool in the analysis of this section.

Lemma 2

(Classical comparison principle) Given $0 \le \varepsilon \le 1/n$ and a continuous right-hand side $f \in C(\varOmega )$, where we assume $f\ge 0$ if $\varepsilon =0$, let $v^* \in C(\overline{\varOmega })$ resp. $v_* \in C(\overline{\varOmega })$ be a super- resp. subsolution to the PDE

$$\begin{aligned} F_\varepsilon (f;x,\textrm{D}^2 v) = 0 \text { in } \varOmega . \end{aligned}$$

(6)

If $v_* \le v^*$ on $\partial \varOmega$, then $v_* \le v^*$ in $\overline{\varOmega }$.

Proof

The proof applies the arguments from [9, Section 3] to the PDE (6) and can follow [12, Lemma 3.6] with straightforward modifications; further details are therefore omitted. $\square$

An extended version of Lemma 2 is the following.

Lemma 3

(Comparison principle) Given any $0 \le \varepsilon _* \le \varepsilon ^* \le 1/n$ and $f_*,f^* \in C(\varOmega )$ with $f_* \le f^*$ in $\varOmega$, where we assume $f_*\ge 0$ if $\varepsilon _*=0$, let $v_*, v^* \in C(\overline{\varOmega })$ be viscosity solutions to

$$\begin{aligned} F_{\varepsilon ^*}(f_*; x, \textrm{D}^2 v^*) = 0 \text { in } \varOmega \quad \text {and}\quad F_{\varepsilon _*}(f^*; x, \textrm{D}^2 v_*) = 0 \text { in } \varOmega . \end{aligned}$$

If $v_* \le v^*$ on $\partial \varOmega$, then $v_* \le v^*$ in $\overline{\varOmega }$.

Proof

Given any test function $\varphi \in C^2(\varOmega )$ and $x \in \varOmega$ such that $v^* - \varphi$ has a local minimum at x, then $F_{\varepsilon ^*}(f_*; x, \textrm{D}^2 v^*) = 0$ in the sense of viscosity solutions implies $0 \le F_{\varepsilon ^*}(f_*; x, \textrm{D}^2 \varphi (x))$. This, $f_* \le f^*$ in $\varOmega$, and ${\mathbb {S}}(\varepsilon ^*) \subset {\mathbb {S}}(\varepsilon _*)$ show

$$\begin{aligned} 0 \le F_{\varepsilon ^*}(f_*; x, \textrm{D}^2 \varphi (x)) \le F_{\varepsilon _*}(f^*; x, \textrm{D}^2 \varphi (x)), \end{aligned}$$

whence $v^*$ is viscosity supersolution to the PDE $F_{\varepsilon _*}(f^*; x, \textrm{D}^2 v_*) = 0$. Therefore, the comparison principle from Lemma 2 with $v_* \le v^*$ on $\partial \varOmega$ concludes $v_* \le v^*$ in $\overline{\varOmega }$. $\square$

The comparison principle from Lemma 2 allows for the existence and uniqueness of viscosity solutions (2) by Perron’s method. In the next proposition we summarize some results from the literature [6, 12, 17, 19, 24] applied to our context.

Proposition 1

(Properties of HJB equation) Given any $0 \le \varepsilon \le 1/n$, $f \in C(\varOmega ) \cap L^n(\varOmega )$, where we assume $f\ge 0$ if $\varepsilon =0$, and $g \in C(\partial \varOmega )$, there exists a unique viscosity solution $u \in C(\overline{\varOmega })$ to the HJB equation (2). It satisfies (a)–(b):

• (a) (viscosity = Alexandrov) If $\varepsilon = 0$ and $f \ge 0$ is nonnegative, then the viscosity solution to the HJB equation (2) and the Alexandrov solution to the Monge–Ampère equation (1) coincide.

• (b) (interior regularity for HJB) If $\varepsilon > 0$ and $f \in C^{0,\alpha }(\varOmega )$ with $0< \alpha < 1$, then $u \in C(\overline{\varOmega }) \cap C^{2,\kappa }_{\text {loc}}(\varOmega )$ with a constant $0< \kappa < 1$ that solely depends on $\alpha$ and $\varepsilon$.

• (c) (interior regularity for Monge–Ampère) If $\varepsilon = 0$, $f \in C^{0,\alpha }(\varOmega )$ with $0< \alpha < 1$, $f > 0$ in $\overline{\varOmega }$, and $g \in C^{1,\beta }(\partial \varOmega )$ with $\beta > 1 - 2/n$, then $u \in C(\overline{\varOmega }) \cap C^{2,\alpha }_{\text {loc}}(\varOmega )$.

Proof

On the one hand, an elementary reasoning as in the proof of Lemma 3 proves that the viscosity solution $v^*$ to the Poisson equation $F_{\varepsilon ^*}(f_*; x, \textrm{D}^2 v^*) = 0$ with $\varepsilon ^* :=1/n$, $f_* :=f$, and Dirichlet data $v^* = g$ on $\partial \varOmega$ is a viscosity supersolution to (2). On the other hand, the Alexandrov solution $v_*$ to the Monge–Ampère equation (1) with the right-hand side |f| [14, Theorem 2.14] is the viscosity solution to the HJB equation $F_{\varepsilon _*}(f^*; x, \textrm{D}^2 v_*) = 0$ with $\varepsilon _* :=0$, $f^* :=|f|$, and Dirichlet data $v_* = g$ on $\partial \varOmega$ [17, Proposition 1.3.4]. Hence, the function $v_*$ is viscosity subsolution to (2). Therefore, Perron’s method [9, Theorem 4.1] and the comparison principle from Lemma 2 conclude the existence and uniqueness of viscosity solutions to (2). The combination of [12, Theorem 3.3 and Theorem 3.5] with [17, Proposition 1.3.4] implies the assertion in (a). The interior regularity in (b) is a classical result from [6, 24]. For the Monge–Ampère equation, the interior regularity in (c) holds under the assumption that the Alexandrov solution u is strictly convex [14, Corollary 4.43]. Sufficient conditions for this are that $f > 0$ is bounded away from zero and $g \in C^{1,\beta }(\partial \varOmega )$ is sufficiently smooth [14, Corollary 4.11]. $\square$

Some comments are in order, before we state a precise version of the $L^\infty$ stability estimate (4) from the introduction. In general, these estimates arise from the Alexandrov–Bakelman–Pucci maximum principle for the uniform elliptic Pucci operator, cf. [3] and the references therein for further details. However, the constant therein may depend on the ellipticity constant of $F_\varepsilon$ and therefore, on $\varepsilon$. In the case of the HJB equation (2) that approximates the Monge–Ampère equation (1) as $\varepsilon \rightarrow 0$, the Alexandrov maximum principle is the key argument to avoid a dependency on $\varepsilon$. Recall the constant $c_n$ from Lemma 1.

Theorem 1

($L^\infty$ stability) Given a nonnegative parameter $0 \le \varepsilon \le 1/n$ and right-hand sides $f_1, f_2 \in C(\overline{\varOmega })$, where we assume $f_1,f_2\ge 0$ if $\varepsilon =0$, let $v_1, v_2 \in C(\overline{\varOmega })$ be viscosity solutions to the HJB equation $F_\varepsilon (f_j; x, \textrm{D}^2 v_j) = 0 \text { in } \varOmega$ for $j \in \{1,2\}$. Then, for any subset $\omega \subset \varOmega$,

$$\begin{aligned} \Vert v_1 - v_2\Vert _{L^\infty (\omega )} \le \Vert v_1 - v_2\Vert _{L^\infty (\partial \varOmega )} + \frac{C}{n} \max _{x \in \overline{\omega }} \textrm{dist}(x,\partial \varOmega )^{1/n}\Vert f_1-f_2\Vert _{L^n(\varOmega )} \end{aligned}$$

(7)

with the constant $C :=c_n \textrm{diam}(\varOmega )^{(n-1)/n}$. In particular,

$$\begin{aligned} \Vert v_1 - v_2\Vert _{L^\infty (\varOmega )} \le \Vert v_1 - v_2\Vert _{L^\infty (\partial \varOmega )} + \frac{C}{n} (\textrm{diam}(\varOmega )/2)^{1/n}\Vert f_1-f_2\Vert _{L^n(\varOmega )}. \end{aligned}$$

(8)

Proof

The proof is divided into two steps.

Step 1: The first step establishes (7) under the assumptions $f_2 \le f_1$ in $\overline{\varOmega }$ and $v_1 \le v_2$ on $\partial \varOmega$. For $f_\varDelta :=f_1 - f_2 \ge 0$, let the sequence $(f_{\varDelta ,k})_{k \in {\mathbb {N}}}$ of smooth functions $f_{\varDelta ,k} \in C^\infty (\overline{\varOmega })$ approximate $f_\varDelta \in C(\overline{\varOmega })$ from above such that $f_\varDelta \le f_{\varDelta ,k}$ and $0 < f_{\varDelta ,k}$ in $\overline{\varOmega }$ for all $k \in {\mathbb {N}}$ and $\lim _{k \rightarrow \infty } \Vert f_\varDelta - f_{\varDelta ,k}\Vert _{L^\infty (\varOmega )} = 0$. Let $w_k \in C(\overline{\varOmega })$ be viscosity solutions to the PDE, for all $k \in {\mathbb {N}}$,

$$\begin{aligned} F_\varepsilon (f_{\varDelta ,k}; x, \textrm{D}^2 w_k) = 0 \text { in } \varOmega \quad&\text {and}\quad w_k = 0 \text { on } \partial \varOmega . \end{aligned}$$

(9)

Since $v_1 \le v_2$ on $\partial \varOmega$ and $f_2 \le f_1$ by assumption of Step 1, Lemma 3 proves

$$\begin{aligned} v_1 \le v_2 \text { in } \overline{\varOmega }. \end{aligned}$$

(10)

Proposition 1(b)–(c) provides the interior regularity $w_k \in C^{2,\alpha }_\textrm{loc}(\varOmega )$ for some positive parameter $\alpha$ that (possibly) depends on $\varepsilon$. In particular, $w_k \in C^2(\varOmega )$ is a classical solution to the PDE (9). We define the continuous function $v_* :=v_2 - \Vert v_1 - v_2\Vert _{L^\infty (\partial \varOmega )} + w_k \in C(\overline{\varOmega })$. Given any $x \in \varOmega$ and $\varphi \in C^2(\varOmega )$ such that $v_* - \varphi = v_2 - (\Vert v_1 - v_2\Vert _{L^\infty (\partial \varOmega )} - w_k + \varphi )$ has a local maximum at x, the function $\psi :=\Vert v_1 - v_2\Vert _{L^\infty (\partial \varOmega )} - w_k + \varphi \in C^2(\varOmega )$ is smooth and, therefore, an admissible test function in the definition of viscosity solutions. Since $v_2$ is viscosity solution to $F_\varepsilon (f_2; x, \textrm{D}^2 v_2) = 0$, $F_\varepsilon (f_2; x, \textrm{D}^2 \psi (x)) \le 0$ follows. This, $\textrm{D}^2 \psi = \textrm{D}^2(\varphi - w_k)$, the sub-additivity $\sup (X + Y) \le \sup X + \sup Y$ of the supremum, $f_\varDelta \le f_{\varDelta ,k}$, and $F_\varepsilon (f_{\varDelta ,k}; x, \textrm{D}^2 w_k(x)) = 0$ from (9) lead to

$$\begin{aligned} F_\varepsilon (f_1; x, \textrm{D}^2 \varphi (x))&\le F_\varepsilon (f_2; x, \textrm{D}^2 \psi (x)) + F_\varepsilon (f_\varDelta ; x, \textrm{D}^2 w_k(x))\\&\le F_\varepsilon (f_2; x, \textrm{D}^2 \psi (x)) + F_\varepsilon (f_{\varDelta ,k}; x, \textrm{D}^2 w_k(x)) \le 0, \end{aligned}$$

whence $v_*$ is viscosity subsolution to the PDE $F_\varepsilon (f_1; x, \textrm{D}^2 v) = 0$ in $\varOmega$. Therefore, $v_* \le v_1$ on $\partial \varOmega$ by design and the comparison principle from Lemma 2 provide

$$\begin{aligned} v_* \le v_1 \text { in } \overline{\varOmega }. \end{aligned}$$

(11)

On the one hand, the zero function with $F_\varepsilon (f_{\varDelta ,k}; x, 0) \ge 0$ is a viscosity supersolution to $F_\varepsilon (f_{\varDelta ,k}; x, \textrm{D}^2 w_k) = 0$. Hence, the comparison principle from Lemma 2 shows $w_k \le 0$ in $\overline{\varOmega }$. On the other hand, Proposition 1(a) proves that the Alexandrov solution $z_k \in C(\overline{\varOmega })$ to $\det \textrm{D}^2 z_k = (f_{\varDelta ,k}/n)^n$ with homogenous boundary is viscosity solution to $F_0(f_{\varDelta ,k}; x, \textrm{D}^2 z_k) = 0$ and Lemma 3 reveals $z_k \le w_k$, whence $z_k \le w_k \le 0$ in $\overline{\varOmega }$. Consequently, the Alexandrov maximum principle from Lemma 1 and ${\mathcal {L}}^n(\partial z_k(\varOmega ))^{1/n} = \Vert (f_{\varDelta ,k}/n)^n\Vert _{L^1(\varOmega )}^{1/n} = \Vert f_{\varDelta ,k}\Vert _{L^n(\varOmega )}/n$ imply

$$\begin{aligned} 0 \le - w_k \le - z_k \le \frac{C}{n}\max _{x \in \overline{\omega }} \textrm{dist}(x,\partial \varOmega )^{1/n}\Vert f_{\varDelta ,k}\Vert _{L^n(\varOmega )} \quad \text {in } \overline{\omega } \end{aligned}$$

(12)

for any subset $\omega \subset \varOmega$. The combination of (10)–(12) with $v_* = v_2 - \Vert v_1 - v_2\Vert _{L^\infty (\partial \varOmega )} + w_k$ results in

$$\begin{aligned} \Vert v_1 - v_2\Vert _{L^\infty (\omega )}&\le \Vert v_2 - v_*\Vert _{L^\infty (\omega )} = \Vert v_1 - v_2\Vert _{L^\infty (\partial \varOmega )} + \Vert w_k\Vert _{L^\infty (\omega )}\\&\le \Vert v_1 - v_2\Vert _{L^\infty (\partial \varOmega )} + \frac{C}{n} \max _{x \in \overline{\omega }} \textrm{dist}(x,\partial \varOmega )^{1/n}\Vert f_{\varDelta ,k}\Vert _{L^n(\varOmega )}. \end{aligned}$$

A passage of the right-hand side to the limit as $k \rightarrow \infty$ and $\lim _{k \rightarrow \infty } \Vert f_{\varDelta ,k}\Vert _{L^n(\varOmega )} = \Vert f_\varDelta \Vert _{L^n(\varOmega )}$ conclude (7).

Step 2: The second step establishes (7) without the additional assumptions from Step 1. For the functions $f_* :=\min \{f_1,f_2\}$, $f^* :=\max \{f_1,f_2\}$, and $f_\varDelta :=f^* - f_* = |f_1 - f_2| \ge 0$, let $v^*,v_* \in C(\overline{\varOmega })$ be viscosity solutions to the PDE

$$\begin{aligned}&F_\varepsilon (f_*; x, \textrm{D}^2 v^*) = 0 \text { in } \varOmega \quad \text {and}\quad v^* = \max \{v_1,v_2\} \text { on } \partial \varOmega , \end{aligned}$$

(13)

$$\begin{aligned}&F_\varepsilon (f^*; x, \textrm{D}^2 v_*) = 0 \text { in } \varOmega \quad \text {and}\quad v_* = \min \{v_1,v_2\} \text { on } \partial \varOmega , \end{aligned}$$

(14)

Since $f_* \le f_j \le f^*$ and $v_* \le v_j \le v^*$ on $\partial \varOmega$ for $j \in \{1,2\}$, Lemma 3 verifies $v_* \le \{v_1,v_2\} \le v^*$ in $\overline{\varOmega }$, whence

$$\begin{aligned} \Vert v_1 - v_2\Vert _{L^\infty (\omega )} \le \Vert v^* - v_*\Vert _{L^\infty (\omega )}\quad \text {for any open subset } \omega \subset \varOmega . \end{aligned}$$

(15)

The application of Step 1 to the viscosity solutions $v^*,v_*$ of (13)–(14) with $f_* \le f^*$ and $v_* \le v^*$ on $\partial \varOmega$, and the identity $\max \{a,b\} - \min \{a,b\} = |a - b|$ reveal

$$\begin{aligned} \Vert v^* - v_*\Vert _{L^\infty (\omega )} \le \Vert v_1 - v_2\Vert _{L^\infty (\partial \varOmega )} + \frac{C}{n} \max _{x \in \overline{\omega }}\textrm{dist}(x,\partial \varOmega )^{1/n}\Vert f_1 - f_2\Vert _{L^n(\varOmega )}. \end{aligned}$$

The combination of this with (15) concludes (7). $\square$

The stability estimate from Theorem 1 motivates a solution concept for the HJB equation (2) with $L^n$ (rather than continuous) right-hand sides.

Lemma 4

(Generalized viscosity solution) Given $f \in L^n(\varOmega )$, $g \in C(\partial \varOmega )$ and $0 \le \varepsilon \le 1/n$, where we assume $f\ge 0$ if $\varepsilon =0$, there exists a unique function $u \in C(\overline{\varOmega })$ such that u is the uniform limit of any sequence $(u_j)_{j \in {\mathbb {N}}}$ of viscosity solutions $u_j \in C(\overline{\varOmega })$ to

$$\begin{aligned} F_\varepsilon (f_j; x, \textrm{D}^2 u_j) = 0 \text { in } \varOmega \quad \text {and}\quad u_j = g_j \text { on } \partial \varOmega \end{aligned}$$

(16)

for right-hand sides $f_j \in C(\overline{\varOmega })$ and Dirichlet data $g_j \in C(\overline{\varOmega })$ with $\lim _{j \rightarrow \infty } \Vert f - f_j\Vert _{L^n(\varOmega )} = 0$ and $\lim _{j \rightarrow \infty } \Vert g - g_j\Vert _{L^\infty (\partial \varOmega )} = 0$. The function u is called generalized viscosity solution to (2). If $\varepsilon = 0$ and $f \ge 0$, then the generalized viscosity solution to (2) and the Alexandrov solution to (1) coincide.

Proof

Let $(f_j)_{j \in {\mathbb {N}}} \subset C(\overline{\varOmega })$ (resp. $(g_j)_{j \in {\mathbb {N}}} \subset C(\overline{\varOmega })$) approximate f in $L^n(\varOmega )$ (resp. g in $C(\partial \varOmega )$). For any index $j,k \in {\mathbb {N}}$, the stability estimate (8) from Theorem 1 provides

$$\begin{aligned} \Vert u_j - u_k\Vert _{L^\infty (\varOmega )} \le \Vert g_j - g_k\Vert _{L^\infty (\partial \varOmega )} + \frac{C}{n} (\textrm{diam}(\varOmega )/2)^{1/n}\Vert f_j - f_k\Vert _{L^n(\varOmega )}. \end{aligned}$$

Since $(f_j)_{j \in {\mathbb {N}}}$ (resp. $(g_j)_{j \in {\mathbb {N}}}$) is a Cauchy sequence in $L^n(\varOmega )$ (resp. $C(\partial \varOmega )$), this implies that $(u_j)_{j \in {\mathbb {N}}}$ is a Cauchy sequence in the Banach space $C(\overline{\varOmega })$ endowed with the $L^\infty$ norm. Therefore, there exists $u \in C(\overline{\varOmega })$ with $\lim _{j \rightarrow \infty } \Vert u - u_j\Vert _{L^\infty (\varOmega )} = 0$. It remains to prove that u is independent of the choice of the approximation sequences for f and g. To this end, let $(\widetilde{f}_j)_{j \in {\mathbb {N}}}$ be another sequence of continuous functions $\widetilde{f}_j \in C(\overline{\varOmega })$ with $\lim _{j \rightarrow \infty } \Vert f - \widetilde{f}_j\Vert _{L^n(\varOmega )} = 0$. Then the sequence $(\widetilde{u}_j)_{j \in {\mathbb {N}}}$ of viscosity solutions $\widetilde{u}_j \in C(\overline{\varOmega })$ to (16) with $f_j$ replaced by $\widetilde{f}_j$ converges uniformly to some $\widetilde{u} \in C(\overline{\varOmega })$. The stability estimate (8) from Theorem 1 shows

$$\begin{aligned} \Vert u_j - \widetilde{u}_j\Vert _{L^\infty (\varOmega )} \le \frac{C}{n}(\textrm{diam}(\varOmega )/2)^{1/n}\Vert f_j - \widetilde{f}_j\Vert _{L^n(\varOmega )} \end{aligned}$$

for any $j \in {\mathbb {N}}$. The right-hand side of this vanishes in the limit and the left-hand side converges to $\Vert u - \widetilde{u}\Vert _{L^\infty (\varOmega )}$ as $j \rightarrow \infty$, whence $u = \widetilde{u}$ in $\overline{\varOmega }$. If $f \ge 0$, then there exists a sequence $(f_j)_{j \in {\mathbb {N}}}$ of nonnegative continuous functions $0 \le f_j \in C(\overline{\varOmega })$ with $\lim _{j \rightarrow \infty } \Vert f - f_j\Vert _{L^\infty (\varOmega )}$ (e.g., from convolution with a nonnegative mollifier). Proposition 1(a) provides, for all $j \in {\mathbb {N}}$, that the viscosity solution $u_j$ to (16) with $\varepsilon = 0$ is the Alexandrov solution to $\det \textrm{D}^2 u_j = f_j$ in $\varOmega$. Since $u_j$ converges uniformly to the generalized viscosity solution u to (2), the stability of Alexandrov solutions [14, Corollary 2.12 and Proposition 2.16] concludes that u is the Alexandrov solution to (1). $\square$

We note that, for the HJB equation (2) with $\varepsilon > 0$, the generalized viscosity solutions coincide with the $L^n$ viscosity solutions from the literature [3, 18]. While our approach is limited to a smaller class of operators, it requires less technical effort. By approximation of the right-hand sides, the stability estimates from Theorem 1 also applies to generalized viscosity solutions to the HJB equation (2).

Corollary 1

(Extended $L^\infty$ stability) Given any $0 \le \varepsilon \le 1/n$, $f_j \in L^n(\varOmega )$, where we assume $f_j\ge 0$ if $\varepsilon =0$, the generalized viscosity solutions $v_j \in C(\overline{\varOmega })$ to $F_\varepsilon (f_j;x,\textrm{D}^2 v_j) = 0$ in $\varOmega$ for $j \in \{1,2\}$ satisfy (7)–(8).

Proof

For any index $j \in \{1,2\}$, there exists a sequence $(f_{j,k})_{j \in {\mathbb {N}}}$ of smooth functions $f_{j,k} \in C^\infty (\overline{\varOmega })$ that approximates $f_j$ in $L^n(\varOmega )$, i.e., $\lim _{k \rightarrow \infty } \Vert f_j - f_{j,k}\Vert _{L^n(\varOmega )} = 0$. Given any $j \in \{1,2\}$ and $k \in {\mathbb {N}}$, let $v_{j,k} \in C(\overline{\varOmega })$ denote the viscosity solution to the HJB equation $F_\varepsilon (f_{j,k}; x, \textrm{D}^2 v_{j,k}) = 0$ in $\varOmega$ and $v_{j,k} = v_j$ on $\partial \varOmega$. The $L^\infty$ stability estimate (7) from Theorem 1 shows, for any $k \in {\mathbb {N}}$, that

$$\begin{aligned} \Vert v_{1,k} - v_{2,k}\Vert _{L^\infty (\omega )}&\le \Vert v_1 - v_2\Vert _{L^\infty (\partial \varOmega )} + \frac{C}{n} \max _{x \in \overline{\omega }} \textrm{dist}(x,\partial \varOmega )^{1/n}\Vert f_{1,k} - f_{2,k}\Vert _{L^n(\varOmega )}. \end{aligned}$$

The left-hand side of this converges to $\Vert v_1 - v_2\Vert _{L^\infty (\varOmega )}$ by the definition of generalized viscosity solutions in Lemma 4. Hence, $\lim _{k \rightarrow \infty } \Vert f_{1,k} - f_{2,k}\Vert _{L^n(\varOmega )} = \Vert f_1 - f_2\Vert _{L^n(\varOmega )}$ concludes the proof. $\square$

Remark 1

($L^\infty$ stability for Alexandrov solutions) If the right-hand sides $0 \le f_1,f_2 \in L^n(\varOmega )$ are nonnegative, then the generalized solutions $v_1,v_2$ from Corollary 1 are Alexandrov solutions to $\det \textrm{D}^2 v_j = (f_j/n)^n$, cf. Lemma 4. Therefore, Corollary 1 provides $L^\infty$ stability estimates for Alexandrov solutions.

The convexity and uniform ellipticity of the differential operator $F_\varepsilon$ in ${\mathbb {S}}$ lead to existence (and uniqueness) of strong solutions $u_\varepsilon \in C(\overline{\varOmega }) \cap W^{2,n}_\textrm{loc}(\varOmega )$ to (2) for any $\varepsilon > 0$, $f \in L^n(\varOmega )$, and $g \in C(\partial \varOmega )$ [3]. It turns out that strong solutions are generalized viscosity solutions. For the purpose of this paper, we only provide a weaker result.

Theorem 2

(Strong solution implies generalized viscosity solution) Let $0 < \varepsilon \le 1/n$, $f \in L^n(\varOmega )$, and $g \in C(\partial \varOmega )$ be given. Suppose that $u_\varepsilon \in W^{2,n}(\varOmega )$ is a strong solution to (2) in the sense that (2) is satisfied a.e. in $\varOmega$. Then this strong solution $u_\varepsilon$ is the unique generalized viscosity solution to (2).

The proof of Theorem 2 utilizes the following elementary result.

Lemma 5

(Computation and stability of right-hand side) Let $\varepsilon > 0$ be given. For any $M \in {\mathbb {S}}$, there exists a unique $\xi (M) \in {\mathbb {R}}$ such that $\max _{A \in {\mathbb {S}}(\varepsilon )}(- A: M + \xi (M)\root n \of {\det A}) = 0$. Furthermore, any $M, N \in {\mathbb {S}}$ satisfy the stability $|\xi (M) - \xi (N)| \le C(\varepsilon )|M - N|$ with a constant $C(\varepsilon )$ depending on the regularization parameter $\varepsilon$.

Proof

Given a symmetric matrix $M \in {\mathbb {S}}$, define the continuous real-valued function

$$\begin{aligned} \varPsi _{M}(\xi ) :=\max _{A \in {\mathbb {S}}(\varepsilon )} (-A:M + \xi \root n \of {\det A}). \end{aligned}$$

(17)

Since $\varPsi _M$ is strictly monotonically increasing with the limits $\lim _{\xi \rightarrow -\infty } \varPsi _M = -\infty$ and $\lim _{\xi \rightarrow \infty } \varPsi _M = +\infty$, there exists a unique root $\xi (M)$ such that $\varPsi _M(\xi (M)) = 0$. For any $M,N \in {\mathbb {S}}$, the inequality $\max X - \max Y \le \max (X - Y)$ shows

$$\begin{aligned} 0 = \varPsi _{M}(\xi (M)) - \varPsi _{N}(\xi (N)) \le \varPsi _{M - N}(\xi (M) - \xi (N)). \end{aligned}$$

(18)

Let $A \in {\mathbb {S}}(\varepsilon )$ be chosen such that $\varPsi _{M - N}(\xi (M) - \xi (N)) = -A:(M-N) + (\xi (M) - \xi (N))\root n \of {\det A}$. Then it follows from (18) that

$$\begin{aligned} \xi (N) - \xi (M) \le A:(N-M)/\root n \of {\det A}\le |A||M-N|/\root n \of {\det A}. \end{aligned}$$

(19)

Exchanging the roles of M and N in (19) leads to $\xi (M) - \xi (N) \le |B||M-N|/\root n \of {\det B}$ for some $B \in {\mathbb {S}}(\varepsilon )$. Since $|A|/\root n \of {\det A} \le 1/(\sqrt{\varepsilon ^{n-1}(1-(n-1)\varepsilon )})$ holds for any $A \in {\mathbb {S}}(\varepsilon )$, the combination of this with (19) concludes $|\xi (N) - \xi (M)| \le |M-N|/\root n \of {\varepsilon ^{n-1}(1-(n-1)\varepsilon )}$. $\square$

Proof of Theorem 2

Let $v_j \in C^2(\overline{\varOmega })$ be a sequence of smooth functions that approximate $u_\varepsilon$ with $\lim _{j \rightarrow \infty } \Vert u_\varepsilon - v_j\Vert _{W^{2,n}(\varOmega )} = 0$. Lemma 5 proves that there exists a (unique) function $f_j :=\xi (\textrm{D}^2 v_j)$ with $F_\varepsilon (f_j;x,\textrm{D}^2 v_j) = 0$ in $\varOmega$. We apply the stability from Lemma 5 twice. First, $|f_j(x) - f_j(y)| \le C(\varepsilon )|\textrm{D}^2 v_j(x) - \textrm{D}^2 v_j(y)|$ for any $x,y \in \varOmega$ implies continuity $f_j \in C(\overline{\varOmega })$ of $f_j$ and second, $|f(x) - f_j(x)| \le C(\varepsilon )|\textrm{D}^2 u_\varepsilon (x) - \textrm{D}^2 v_j(x)|$ for a.e. $x \in \varOmega$ implies the convergence $\lim _{j \rightarrow \infty } \Vert f - f_j\Vert _{L^n(\varOmega )} = 0$. Notice from the Sobolev embedding that $v_j$ converges uniformly to $u_\varepsilon$ in $\overline{\varOmega }$ as $j \rightarrow \infty$. In conclusion, $u_\varepsilon$ is the uniform limit of classical (and in particular, viscosity) solutions $v_j$ such that the corresponding right-hand sides and Dirichlet data converge in the correct norm, i.e., $\lim _{j \rightarrow \infty } \Vert f - f_j\Vert _{L^n(\varOmega )} = 0$ and $\lim _{j \rightarrow \infty }\Vert g - v_j\Vert _{L^\infty (\partial \varOmega )} = 0$. Lemma 4 proves that $u_\varepsilon$ is the unique (generalized) viscosity solution. $\square$

3 Convergence of the regularization

This section establishes the uniform convergence of the generalized viscosity solution $u_\varepsilon$ of the regularized HJB equation (2) to the Alexandrov solution u of the Monge–Ampère equation (1) for any nonnegative right-hand side $0 \le f \in L^n(\varOmega )$. The proof is carried out in any space dimension n and does not rely on the concept of strong solutions in two space dimensions from [25, 26]. It departs from a main result of [15].

Theorem 3

(Convergence of regularization for smooth data) Let $f \in C^{0,\alpha }(\varOmega )$, $0 < \lambda \le f \le \varLambda$, and $g \in C^{1,\beta }(\partial \varOmega )$ with positive constants $0< \alpha < 1$, $1 - 2/n< \beta < 1$, and $0 < \lambda \le \varLambda$ be given. Let $u \in C(\overline{\varOmega }) \cap C^{2,\alpha }_{\textrm{loc}}(\varOmega )$ be the unique classical solution to (1) from Proposition 1(c).

• (a) For any sequence $0<(\varepsilon _j)_{j \in {\mathbb {N}}}\le 1/n$ with $\lim _{j \rightarrow \infty } \varepsilon _j = 0$, the sequence $(u_{\varepsilon _j})_{j \in {\mathbb {N}}}$ of classical solutions $u_{\varepsilon _j} \in C(\overline{\varOmega }) \cap C^2(\varOmega )$ to (2) with $\varepsilon :=\varepsilon _j$ from Proposition 1(b) converges uniformly to u in $\varOmega$ as $j \rightarrow \infty$.

• (b) If $g\equiv 0$, $f \in C^{2,\alpha }(\varOmega )$, and $f > 0$ in $\overline{\varOmega }$, then, for some constant C and all $0 < \varepsilon \le 1/n$, the generalized viscosity solution $u_\varepsilon$ to (2) satisfies

  $$\begin{aligned} \Vert u - u_\varepsilon \Vert _{L^\infty (\varOmega )} \lesssim \varepsilon ^{1/(n^2(2n+3))}. \end{aligned}$$

Proof

The proof of Theorem 3 can follow the lines of the proof of [15, Theorem 4.1], where Lemma 6 below replaces its counterpart [15, Lemma 4.2] in two space dimensions. We note that the assumption $g \in H^2(\varOmega )$ in [15, Theorem 4.1] is only required for the existence of strong solutions $u_\varepsilon \in H^2(\varOmega )$ and can be dropped. Further details of the proof are omitted. $\square$

Lemma 6

(Effect of regularization) Given $0 < \varepsilon \le 1/n$, $M \in {\mathbb {S}}$, and $\xi > 0$, suppose that $|M|_n^n \le \xi ^n(1/\varepsilon - (n-1))/n^n$ and $\max _{A \in {\mathbb {S}}(0)} (-A:M + \xi \sqrt{\det A}) = 0$, then $\max _{A \in {\mathbb {S}}(\varepsilon )} (-A:M + \xi \sqrt{\det A}) = 0$.

Proof

The assumption $\max _{A \in {\mathbb {S}}(0)} (-A:M + \xi \sqrt{\det A}) = 0$ implies that $M > 0$ is positive definite and $\det M = (\xi /n)^n$ [19, p. 51]. Let $\varrho _1,\dots ,\varrho _n$ denote the positive eigenvalues of M and $t_j :=\varrho _j^{-1}/(\sum _{k=1}^n \varrho _k^{-1})$ for $j = 1,\dots ,n$. By design of $t_j$,

$$\begin{aligned} \varrho _j^{-1} = t_j\left( \frac{\varrho _1^{-1} \dots \varrho _n^{-1}}{t_1 \dots t_n} \right) ^{1/n}, \end{aligned}$$

whence $\varrho _j = \xi (t_1 \dots t_n)^{1/n}/(n t_j)$. Without loss of generality, suppose that $t_1 \le t_2 \le \dots \le t_n$. The elementary bound $t_1 \dots t_n \ge t_1^{n-1}(1-(n-1)t_1)$ proves

$$\begin{aligned} \xi ^n(1-(n-1)t_1)/t_1 \le \xi ^n(t_1 \dots t_n)/(nt_1)^n = n^n \varrho _1^n \le n^n|M|_n^n. \end{aligned}$$

Hence, $1/t_1 \le n^n|M|_n^n/\xi ^n + (n-1) \le 1/\varepsilon$ by assumption and so, $t_1 \ge \varepsilon$. In particular, $\varepsilon \le t_1 \le \dots \le t_n$ and $t_1 + \dots + t_n = 1$. Notice that $t :=(t_1,\dots ,t_n) \in {\mathbb {R}}^n$ maximizes the scalar-valued function $g : {\mathbb {R}}^{n} \rightarrow {\mathbb {R}}$ with

$$\begin{aligned} \psi (s) :=-s_1 \varrho _1 - \dots - s_n \varrho _n + \xi \root n \of {s_1 \dots s_n} \end{aligned}$$

among $s \in S(0)$ with $S(\varepsilon ) :=\{s = (s_1,\dots ,s_n) : s \ge \varepsilon \text { and } s_1 + \dots + s_n = 1\}$. Since $\psi (t) = \max _{s \in S(0)} \psi (s) = \max _{A \in {\mathbb {S}}(0)} (-A:M + \xi \sqrt{\det A})$ [19, pp. 51–52] and $t \in S(\varepsilon )$, this implies that $0 = \psi (t) = \max _{A \in {\mathbb {S}}(\varepsilon )} (-A:M + \xi \sqrt{\det A})$. $\square$

The approximation of nonsmooth data leads to the following convergence result under (almost) minimal assumptions (general Borel measures as right-hand sides are excluded).

Theorem 4

(Convergence of regularization) Let a sequence $(\varepsilon _j)_{j \in {\mathbb {N}}} \subset (0,1/n]$ with $\lim _{j \rightarrow \infty } \varepsilon _j = 0$, a nonnegative right-hand side $0 \le f \in L^n(\varOmega )$, and Dirichlet data $g \in C(\partial \varOmega )$ be given. Then the sequence $(u_j)_{j \in {\mathbb {N}}}$ of generalized viscosity solutions $u_j \in C(\overline{\varOmega })$ to

$$\begin{aligned} F_{\varepsilon _j}(f; x, \textrm{D}^2 u_j) = 0 \text { in } \varOmega \quad \text {and}\quad u_j = g \text { on } \partial \varOmega \end{aligned}$$

converges uniformly $\lim _{j \rightarrow \infty } \Vert u - u_j\Vert _{L^\infty (\varOmega )} = 0$ to the Alexandrov solution u to the Monge–Ampère equation (1).

Proof

Recall the constant $c_n$ from Lemma 1 and $C :=c_n \textrm{diam}(\varOmega )^{(n-1)/n}$. Given $\delta > 0$, there exist smooth functions $f_\delta , g_\delta \in C^\infty (\overline{\varOmega })$ such that

• (i) $f_\delta > 0$ in $\overline{\varOmega }$ and $\Vert f-f_\delta \Vert _{L^n(\varOmega )} \le n\delta /(8C(\textrm{diam}(\varOmega )/2)^{1/n})$ (the approximation $f_\delta$ can be constructed by the convolution of f with a nonnegative mollifier plus an additional small constant),

• (ii) $\Vert g - g_\delta \Vert _{L^\infty (\partial \varOmega )} \le \delta /4$.

Notice that the bound $f_\delta > 0$ in $\overline{\varOmega }$ and the smoothness of the Dirichlet data $g_\delta \in C^{\infty }(\partial \varOmega )$ allow for strict convexity of the Alexandrov solution $u_\delta$ to the Monge–Ampère equation $\det \textrm{D}^2 u_\delta = (f_\delta /n)^n$ with Dirichlet data $u_\delta = g_\delta$ on $\partial \varOmega$ [14, Corollary 4.11]. This is a crucial assumption in Theorem 3, which leads to the uniform convergence of the sequence $(u_{\delta ,j})_{j \in {\mathbb {N}}}$ of viscosity solutions $u_{\delta ,j} \in C(\overline{\varOmega })$ to the HJB equation

$$\begin{aligned} F_{\varepsilon _j}(f_\delta ; x, \textrm{D}^2 u_{\delta ,j}) = 0 \text { a.e.~in } \varOmega \quad \text {and}\quad u_{\delta ,j} = g_\delta \text { on } \partial \varOmega \end{aligned}$$

towards $u_\delta$ as $j \rightarrow \infty$. Therefore, there exists a $j_0 \in {\mathbb {N}}$ such that $\Vert u_\delta - u_{\delta ,j}\Vert _{L^\infty (\varOmega )} \le \delta /4$ for all $j \ge j_0$. The stability estimate (8) from Corollary 1 and (i)–(ii) provide

$$\begin{aligned}&\Vert u - u_\delta \Vert _{L^\infty (\varOmega )} + \Vert u_j - u_{\delta ,j}\Vert _{L^\infty (\varOmega )}\\&\quad \le 2\Vert g-g_\delta \Vert _{L^\infty (\partial \varOmega )} + \frac{2C}{n} (\textrm{diam}(\varOmega )/2)^{1/n}\Vert f - f_\delta \Vert _{L^n(\varOmega )} \le 3\delta /4. \end{aligned}$$

This, the triangle inequality, and $\Vert u_\delta - u_{\delta ,j}\Vert _{L^\infty (\varOmega )} \le \delta /4$ verify, for all $j \ge j_0$, that $\Vert u - u_j\Vert _{L^\infty (\varOmega )} \le \delta$, whence $u_j$ converges uniformly to u as $j \rightarrow \infty$. $\square$

4 A posteriori error estimate

In this section we prove an a posteriori error bound for a given approximation $v_h$ to the viscosity solution $u_\varepsilon$ of the regularized PDE (2) as well as the Alexandrov solution u of the Monge–Ampère equation. In what follows we assume a given finite partition ${\mathcal {T}}$ of $\overline{\varOmega }$ of closed polytopes such that the interiors of any distinct $T,K\in {\mathcal {T}}$ are disjoint and the union over ${\mathcal {T}}$ equals $\overline{\varOmega }$. Let $V_h\subset C^{1,1}(\overline{\varOmega }) \cap W^{2,\infty }(\varOmega )$ be a subspace of functions in $C^{2}(T)$ when restricted to any set $T\in {\mathcal {T}}$ of the partition. (Here, $C^2$ up to the boundary of T means that there exists a sufficiently smooth extension of the function $v_h|_{\textrm{int}(T)}$ to T for $v_h \in V_h$.) In practical examples, we think of $V_h$ as a space of $C^1$-regular finite element functions. Given any $v \in C(\varOmega )$, its convex envelope is defined as

$$\begin{aligned} \varGamma _{v} (x) :=\sup _{\begin{array}{c} w : {\mathbb {R}}^n \rightarrow {\mathbb {R}} \text { affine}\\ w\le v \end{array}} w (x) \quad \text {for any }x\in \varOmega . \end{aligned}$$

(20)

Let ${\mathcal {C}}_v :=\{x\in \varOmega : v (x) = \varGamma _{v}(x)\}$ denote the contact set of v.

Theorem 5

(Guaranteed error control for Monge–Ampère) Given a nonnegative right-hand side $f \in L^n(\varOmega )$ and $g \in C(\partial \varOmega )$, let $u \in C(\overline{\varOmega })$ be the Alexandrov solution to (1). Let $v_h\in V_h$ with its convex envelope $\varGamma _{v_h}$ be given and define $f_h :=\chi _{{\mathcal {C}}_{v_h}} n (\det \textrm{D}^2 v_h)^{1/n}$. For any convex subset $\varOmega ' \subset \varOmega$, we have, for the constant $c_n$ from Lemma 1, that

$$\begin{aligned} \Vert u - \varGamma _{v_h}\Vert _{L^\infty (\varOmega )}&\le \limsup _{x \rightarrow \partial \varOmega } |(g - \varGamma _{v_h})(x)| + \frac{c_n}{2^{1/n}n}\textrm{diam}(\varOmega ')\Vert f - f_h\Vert _{L^n(\varOmega ')}\nonumber \\&\quad + \frac{c_n}{n}\textrm{diam}(\varOmega )^{(n-1)/n} \max _{x\in \overline{\varOmega \setminus \varOmega '}} {\text {dist}}(x,\partial \varOmega )^{1/n} \Vert f - f_h\Vert _{L^n(\varOmega )}=:\textrm{RHS}_0. \end{aligned}$$

(21)

The practical evaluation of $f_h$ is described in Sect. 5.1. The proof of Theorem 5 requires the following result on the Monge–Ampère measure of the convex envelope $\varGamma _{v_h}$.

Lemma 7

(MA measure of the convex envelope) The convex envelope $\varGamma _{v_h}$ of any $v_h\in V_h$ satisfies $\det \textrm{D}^2 \varGamma _{v_h} = \widetilde{f}_h \,\textrm{d}x$ in the sense of Monge–Ampère measure with the nonnegative function $\widetilde{f}_h :=\chi _{{\mathcal {C}}_{v_h}}\det \textrm{D}^2 v_h \in L^\infty (\varOmega )$.

Proof

We first claim that $\partial \varGamma _{v_h}(x) = \partial v_h(x) = \{\nabla v_h(x)\}$ holds for all $x \in \varOmega \cap {\mathcal {C}}_{v_h}$. In fact, if $p \in \partial \varGamma _{v_h}(x)$, then $\ell _{x,p}(z) :=\varGamma _{v_h}(x) + p \cdot (z - x)$ is a supporting hyperplane touching $\varGamma _{v_h}$ from below at x. By design of the convex envelope $\varGamma _{v_h}$, $\ell _{x,p} \le v_h$. Since $\ell _{x,p}(x) = v_h(x)$ because $x \in \varOmega \cap {\mathcal {C}}_{v_h}$, $\ell _{x,p}$ touches $v_h$ at x from below. We deduce $p = \nabla v_h(x)$ from the differentiability of $v_h$. The claim then follows from the fact that the subdifferential $\partial \varGamma _{v_h}$ is nonempty in $\varOmega$ [23, Theorem 23.4]. The set $\partial \varGamma _{v_h}(\varOmega {\setminus } {\mathcal {C}}_{v_h})$ has Lebesgue measure zero [10, p. 995] and $\partial \varGamma _{v_h}(x) = \partial v_h(x) = \{\nabla v_h(x)\}$ holds for all $x \in \varOmega \cap {\mathcal {C}}_{v_h}$. Therefore, the area formula [14, Theorem A.31] implies, for any Borel set $\omega \subset \varOmega$, that

$$\begin{aligned} \mu _{\varGamma _{v_h}}(\omega ) = {\mathcal {L}}^n(\partial \varGamma _{v_h}(\omega )) = {\mathcal {L}}^n(\nabla v_h(\omega \cap {\mathcal {C}}_{v_h})) = \int _{\omega \cap {\mathcal {C}}_{v_h}} \det \textrm{D}^2 v_h \,\textrm{d}x. \end{aligned}$$

This formula implies that $\chi _{{\mathcal {C}}_{v_h}} \det \textrm{D}^2 v_h \ge 0$ is a nonnegative function a.e. in $\varOmega$. Consequently, $\mu _{\varGamma _{v_h}} = \widetilde{f}_h \,\textrm{d}x$ with ${\tilde{f}}_h :=\chi _{{\mathcal {C}}_{v_h}} \det \textrm{D}^2 v_h \ge 0$. $\square$

Proof of Theorem 5

Lemma 7 proves that the Monge–Ampère measure $\mu _{\varGamma _{v_h}} = (f_h/n)^n {\d x}$ of $\varGamma _{v_h}$ can be expressed by the $L^1$ density function $(f_h/n)^n$. In particular, $\varGamma _{v_h}$ is the generalized viscosity solution to $F_0(f_h;x,\textrm{D}^2 \varGamma _{v_h}) = 0$ in $\varOmega$. The application of the stability estimate (8) from Corollary 1 on the convex subset $\varOmega ' \subset \varOmega$ instead of $\varOmega$ leads to

$$\begin{aligned} \Vert u - \varGamma _{v_h}\Vert _{L^\infty (\varOmega ')} \le \Vert u - \varGamma _{v_h}\Vert _{L^\infty (\partial \varOmega ')} + \frac{c_n}{2^{1/n} n}\textrm{diam}(\varOmega ') \Vert f - f_h\Vert _{L^n(\varOmega ')}. \end{aligned}$$

The unknown error $\Vert u - \varGamma _{v_h}\Vert _{L^\infty (\partial \varOmega ')} \le \Vert u - \varGamma _{v_h}\Vert _{L^\infty (\varOmega {\setminus } \varOmega ')}$ can be bounded by the local estimate (7) from Corollary 1 with $\omega :=\varOmega {\setminus } \varOmega '$. If $\varGamma _{v_h} \in C(\overline{\varOmega })$ is continuous up to the boundary $\partial \varOmega$ of $\varOmega$, this reads

$$\begin{aligned}&\Vert u - \varGamma _{v_h}\Vert _{L^\infty (\varOmega \setminus \varOmega ')} \le \Vert g - \varGamma _{v_h}\Vert _{L^\infty (\partial \varOmega )}\\&\quad + \frac{c_n}{n} \textrm{diam}(\varOmega )^{(n-1)/n} \max _{x\in \overline{\varOmega \setminus \varOmega '}} {\text {dist}}(x,\partial \varOmega )^{1/n} \Vert f - f_h\Vert _{L^n(\varOmega )}. \end{aligned}$$

Since $\varGamma _{v_h}$ may only be continuous in the domain $\varOmega$, $\Vert g - \varGamma _{v_h}\Vert _{L^\infty (\partial \varOmega )}$ is replaced by $\limsup _{x \rightarrow \partial \varOmega } |(g - \varGamma _{v_h})(x)|$ in general. The combination of the two previously displayed formula concludes the proof. $\square$

We note that, for certain examples, the convex envelope $\varGamma _{v_h}$ of an approximation $v_h$ is continuous up to the boundary.

Proposition 2

(Continuity at boundary) Let $v \in C^{0,1}(\overline{\varOmega })$ be Lipschitz continuous such that $v|_{\partial \varOmega }$ can be extended to a Lipschitz-continuous convex function $g \in C^{0,1}(\overline{\varOmega })$. Then $\varGamma _{v} \in C(\overline{\varOmega })$ and $\varGamma _{v} = v$ on $\partial \varOmega$.

Proof

We first prove the assertion for homogenous boundary condition $v|_{\partial \varOmega } = 0$. Given any point $x \in \varOmega$, let $x' \in \partial \varOmega$ denote a best approximation of x onto the boundary $\partial \varOmega$ so that $|x - x'| = \textrm{dist}(x,\partial \varOmega )$. Define the affine function $a_x(z) :=L(z - x') \cdot (x' - x)/|x - x'|$ for $z \in \varOmega$, where L denotes the Lipschitz constant of the function $v \in C^{0,1}(\overline{\varOmega })$. It is straight-forward to verify that $a_x \le v$ in $\overline{\varOmega }$ [17, p. 12]. Therefore, $-L\textrm{dist}(x,\partial \varOmega ) = a_x(x) \le \varGamma _{v}(x) \le 0$ by definition of the convex envelope. This shows $\varGamma _{v} \in C(\overline{\varOmega })$ with $\varGamma _{v} \equiv 0$ on $\partial \varOmega$. In the general case, we observe that $v - g \in C^{0,1}(\overline{\varOmega })$ is Lipschitz continuous. The first case proves $\varGamma _{v - g} \in C(\overline{\varOmega })$ with $\varGamma _{v - g} = v - g$ on $\partial \varOmega$. We deduce that $w :=g + \varGamma _{v - g} \in C(\overline{\varOmega })$ is a convex function with $w \le v$ in $\varOmega$ and $w = v$ on $\partial \varOmega$. Let $(x_j)_j \subset \varOmega$ be a sequence of points converging to some point $x \in \partial \varOmega$ on the boundary. For a given $\gamma > 0$, there exists, from the uniform continuity of $v - w$ in the compact set $\overline{\varOmega }$, a $\delta > 0$ such that $|(v - w)(x_j) - (v - w)(x)| \le \gamma$ whenever $|x - x_j| \le \delta$. Since $w \le \varGamma _{v} \le v$ in $\varOmega$, this implies $|(v - \varGamma _{v})(x_j)| \le \gamma$ for sufficiently large j. In combination with the triangle inequality and the Lipschitz continuity of v, we conclude $|v(x) - \varGamma _{v}(x_j)| \le \gamma + |v(x) - v(x_j)| \le \gamma + L|x - x_j|$. Therefore, $\lim _{j \rightarrow \infty } \varGamma _{v}(x_j) = v(x)$. $\square$

The theory of this paper also allows for an a posteriori error control for the regularized HJB equation (2). We state this for the sake of completeness as, in general, it is difficult to quantify the regularization error $\Vert u - u_\varepsilon \Vert _{L^\infty (\varOmega )}$.

Theorem 6

(Guaranteed $L^\infty$ error control for uniform elliptic HJB) Given a positive parameter $0 < \varepsilon \le 1/n$ and a $C^1$ conforming finite element function $v_h \in V_h$, there exists a unique $f_h \in L^\infty (\varOmega )$ such that

$$\begin{aligned} F_\varepsilon (f_h; x, \textrm{D}^2 v_h) = 0 \text { a.e.~in } \varOmega . \end{aligned}$$

(22)

The viscosity solution $u_\varepsilon$ to (2) with right-hand side $f \in L^n(\varOmega )$ and Dirichlet data $g \in C(\partial \varOmega )$ satisfies, for any convex subset $\varOmega ' \Subset \varOmega$, that

$$\begin{aligned} \Vert u_\varepsilon - v_h\Vert _{L^\infty (\varOmega )}&\le \Vert g - v_h\Vert _{L^\infty (\partial \varOmega )} + \frac{c_n}{2^{1/n}n}\textrm{diam}(\varOmega ')\Vert f - f_h\Vert _{L^n(\varOmega ')} \nonumber \\&\quad + \frac{c_n}{n}\textrm{diam}(\varOmega )^{(n-1)/n} \max _{x\in \overline{\varOmega \setminus \varOmega '}} {\text {dist}}(x,\partial \varOmega )^{1/n} \Vert f - f_h\Vert _{L^n(\varOmega )} =:\textrm{RHS}_{\varepsilon }. \end{aligned}$$

(23)

Proof

As in the proof of Theorem 2, Lemma 5 provides a (unique) piecewise continuous and essentially bounded function $f_h :=\xi (\textrm{D}^2_{\textrm{pw}} v_h) \in L^\infty (\varOmega )$ with (22). Theorem 2 shows that $v_h$ is the generalized viscosity solution to (22). Therefore, the stability estimates from Corollary 1 can be applied to $u_\varepsilon$ and $v_h$. First, the application of (8) to the subdomain $\varOmega '$ instead $\varOmega$ leads to

$$\begin{aligned} \Vert u_\varepsilon - v_h\Vert _{L^\infty (\varOmega ')} \le \Vert u_\varepsilon - v_h\Vert _{L^\infty (\partial \varOmega ')} + \frac{c_n}{2^{1/n}n}\textrm{diam}(\varOmega ')\Vert f - f_h\Vert _{L^n(\varOmega ')}. \end{aligned}$$

Second, the local estimate (7) with $\omega :=\varOmega \setminus \varOmega '$ implies

$$\begin{aligned}&\Vert u_\varepsilon - v_h\Vert _{L^\infty (\varOmega \setminus \varOmega ')} \le \Vert g - v_h\Vert _{L^\infty (\partial \varOmega )}\\&\quad + \frac{c_n}{n}\textrm{diam}(\varOmega )^{(n-1)/n} \max _{x\in \overline{\varOmega \setminus \varOmega '}} {\text {dist}}(x,\partial \varOmega )^{1/n} \Vert f - f_h\Vert _{L^n(\varOmega )}. \end{aligned}$$

Since $\Vert u_\varepsilon - v_h\Vert _{L^\infty (\partial \varOmega ')} \le \Vert u_\varepsilon - v_h\Vert _{L^\infty (\varOmega {\setminus } \varOmega ')}$, the combination of the two previously displayed formulas concludes the proof. $\square$

In Theorems 5 and 6, it is possible to apply the stability estimate (7) to further subsets of $\varOmega$ to localize the error estimator. We proceed with a discussion on the efficiency of the proposed error estimators. The point of departure is the following efficiency estimate for the uniformly elliptic HJB equation (2), which bounds the error estimator from above by the error in the $W^{2,n}$ norm. It is proven in [3] that, for any $\varepsilon > 0$, the (generalized) viscosity solution $u_\varepsilon$ to (2) satisfies the regularity $u_\varepsilon \in C(\overline{\varOmega }) \cap W^{2,n}_\textrm{loc}(\varOmega )$, i.e., $u_\varepsilon \in W^{2,n}(\omega )$ for any open subset $\omega \Subset \varOmega$ of $\varOmega$.

Theorem 7

(Local efficiency for HJB with respect to $W^{2,n}$) Under the assumptions of Theorem 6, the difference $f-f_h$ satisfies

$$\begin{aligned} \Vert f - f_h\Vert _{L^n(\omega )} \lesssim \varepsilon ^{(1-n)/n}\Vert \textrm{D}^2 (u_\varepsilon - v_h)\Vert _{L^n(\omega )} \end{aligned}$$

(24)

for any open subset $\omega \Subset \varOmega$ of $\varOmega$. If $u_\varepsilon \in W^{2,n}(\varOmega )$, then (24) holds with $\omega = \varOmega$.

Proof

For a.e. $x \in \varOmega$, the solution properties (2), (22), and the elementary bound $\sup X - \sup Y \le \sup (X - Y)$ reveal

$$\begin{aligned} \begin{aligned} 0&= F_\varepsilon (f;x,\textrm{D}^2 u_\varepsilon (x)) - F_\varepsilon (f_h;x,\textrm{D}^2 v_h(x))\\&\le \sup _{A \in {\mathbb {S}}(\varepsilon )} \big (-A:\textrm{D}^2(u_\varepsilon (x) - v_h(x)) + (f(x)-f_h(x))\root n \of {\det A}\big ). \end{aligned} \end{aligned}$$

(25)

The set ${\mathbb {S}}(\varepsilon )$ is compact, whence the supremum on the right-hand side is attained a.e. in $\varOmega$. Let $A=A(x) \in {\mathbb {S}}(\varepsilon )$ denote the maximizer. Since the eigenvalues of A(x) are nonnegative and sum up to 1, the Frobenius norm is bounded as $|A(x)|\le 1$. As in the proof of Lemma 6, the product of the eigenvalues satisfies the bound $\det A(x)\ge \varepsilon ^{n-1}(1-(n-1)\varepsilon )$, so that

$$\begin{aligned} |A(x)|/\root n \of {\det A(x)} \le 1/\root n \of { \varepsilon ^{n-1}(1-(n-1)\varepsilon )}. \end{aligned}$$

Therefore, (25) and the Cauchy inequality lead to

$$\begin{aligned} f_h(x) - f(x) \le \frac{A(x):\textrm{D}^2 (v_h(x) - u_\varepsilon (x))}{\root n \of {\det A(x)} } \lesssim \varepsilon ^{(1-n)/n}|\textrm{D}^2 (u_\varepsilon (x) - v_h(x))|. \end{aligned}$$

Exchanging the roles of $u_\varepsilon$ and $v_h$ leads to the same upper bound for $|f(x) - f_h(x)|$. The assertion then follows from integration over $\omega$. $\square$

We note the discrepancy of norms in the reliability ($L^\infty$) of Theorem 6 and the efficiency estimate ($W^{2,n}$) of Theorem 7. This is to be expected as pointed out by the following remark.

Remark 2

(Equivalence in 2d) Suppose that $n = 2$ and $g \equiv 0$, then it is known from [26] that $u_\varepsilon \in H^2(\varOmega )$ for any fixed $\varepsilon > 0$ and so, $\Vert f - f_h\Vert _{L^2(\varOmega )} \lesssim \varepsilon ^{-1/2}\Vert \textrm{D}^2(u_\varepsilon - v_h)\Vert _{L^2(\varOmega )}$. On the other hand, the stability estimate from [15, Theorem 3.3(a)] provides $\Vert u_\varepsilon - v_h\Vert _{H^2(\varOmega )} \lesssim \Vert f - f_h\Vert _{L^2(\varOmega )}/\varepsilon$ under the natural assumption $v_h = 0$ on $\partial \varOmega$. In conclusion, $\varepsilon ^{1/2}\Vert f - f_h\Vert _{L^2(\varOmega )} \lesssim \Vert u - v_h\Vert _{H^2(\varOmega )} \lesssim \Vert f - f_h\Vert _{L^2(\varOmega )}/\varepsilon$.

Proving the efficiency of the error estimator of Theorem 5 turns out much more challenging: in constrast to what is used in the proof of Theorem 7, in the case of the Monge–Ampère equation it is not known whether the maximizer A(x) in ${\mathbb {S}}(0)$ belongs to a prescribed subset ${\mathbb {S}}(\varepsilon )$. We are therefore formulate a weaker and partly conditional efficiency estimate.

Theorem 8

(Local efficiency for Monge–Ampère with respect to $W^{2,n}$) Let the assumptions of Theorem 5 hold. Assume additionally that $f \in C^{0,\alpha }(\varOmega )$, $0 < \lambda \le f \le \varLambda$, and $g \in C^{1,\beta }(\partial \varOmega )$ with positive parameters $0 < \lambda \le \varLambda$, $0< \alpha < 1$, and $1 - 2/n< \beta < 1$. For any open subset $\omega \Subset \varOmega$ such that $\omega \subset {\mathcal {C}}_{v_h}$ and $|\textrm{det}\,\textrm{D}^2 v_h| \ge \mu$ a.e. in $\omega$ for some $\mu > 0$, there exists $\varepsilon > 0$ depending on $\Vert \textrm{D}^2 u\Vert _{L^\infty (\omega )}$, $\Vert \textrm{D}^2 v_h\Vert _{L^\infty (\omega )}$, $\lambda$, $\mu$ and n such that the difference $f - f_h$ satisfies

$$\begin{aligned} \Vert f - f_h\Vert _{L^n(\omega )} \lesssim \varepsilon ^{(1-n)/n}\Vert \textrm{D}_{\textrm{pw}}^2(u - v_h)\Vert _{L^n(\omega )}. \end{aligned}$$

Proof

The assumptions on f and g allow for interior regularity results [4] and leads to $u\in C^{2,\alpha }(\omega )$. In particular, $\textrm{D}^2 u$ and $\textrm{D}^2 v_h$ are uniformly bounded in $\omega$. From Lemma 6, $F_0(f;x,\textrm{D}^2 u) = 0$, and $F_0(f_h;x,\textrm{D}^2 v_h) = 0$ a.e. in $\omega$, we deduce that there exists a $\varepsilon > 0$, which solely depends on $\Vert \textrm{D}^2 u\Vert _{L^\infty (\omega )}$, $\Vert \textrm{D}^2 v_h\Vert _{L^\infty (\omega )}$, $\lambda$, $\mu$, and n, such that $F_\varepsilon (f;x,\textrm{D}^2 u) = 0$ and $F_\varepsilon (f_h;x,\textrm{D}^2 v_h) = 0$ a.e. in $\omega$. The remaining parts of the proof follow the same lines as for Theorem 7. $\square$

Although the term $\Vert \textrm{D}^2 v_h\Vert _{L^\infty (\omega )}$ can be computed a posteriori for a given discrete approximation $v_h$ of u, we remark that the efficiency in Theorem 8 is mainly of theoretical interest due to the lack of bounds on the $L^\infty$ norm of $\textrm{D}^2 u$.

5 Numerical examples

In this section, we apply the theory from Sect. 4 to numerical benchmarks on the (two-dimensional) unit square domain $\varOmega :=(0,1)^2$.

5.1 Implementation

Some remarks on the practical realization precede the numerical benchmarks of this section.

5.1.1 Setup

Given ${\mathcal {T}}$ as a rectangular partition of the domain $\varOmega$ with the set ${\mathcal {E}}$ of edges, we choose $V_h$ to be the Bogner–Fox–Schmit finite element space [7]. It is the space of global $C^{1,1}(\overline{\varOmega })$ functions that are bicubic when restricted to any element $T\in {\mathcal {T}}$. We compute the discrete approximation in $V_h$ by approximating the regularized problem (3) with a Galerkin method. In the two-dimensional setting, this yields a strongly monotone problem with a unique discrete solution $u_{h,\varepsilon }$ [15]. Since $v_h :=u_{h,\varepsilon }$ is a $C^{1,1}(\overline{\varOmega })$ function, we can apply Theorem 5 to obtain error bounds for $\Vert u - \varGamma _{v_h}\Vert _{L^\infty (\varOmega )}$, which motivates an adaptive scheme as outlined below.

5.1.2 Evaluation of the upper bound of Theorem 5

We proceed as follows for the computation of the right-hand side $\textrm{RHS}_0$ of (21).

Integration of $f - f_h$ for $f_h :=2 \chi _{{\mathcal {C}}_{v_h}} (\det \textrm{D}^2 v_h)^{1/2}$. The integral $\Vert f - f_h\Vert _{L^2(\omega )}$ for any subset $\omega \subset \varOmega$ is computed via numerical integration. Given a set of Gauss points ${\mathcal {N}}_\ell$ associated to the degree of exact integration $\ell$, this reads

$$\begin{aligned} \sum _{T \in {\mathcal {T}}} \sum _{x \in {\mathcal {N}}_\ell \cap T \cap \omega } {\text {meas}}(T) w_{\ell ,T}(x)(f(x) - 2\chi _{{\mathcal {C}}_{v_h}}(x) (\det \textrm{D}^2 v_h(x))^{1/2})^2 \end{aligned}$$

(26)

with some positive weight function $w_{\ell ,T} \in L^\infty (T)$. In the evaluation of (26), the term $\chi _{{\mathcal {C}}_{v_h}}(x)$ determines whether the quadrature point x belongs to the contact set. A point $x \in {\mathcal {N}}_\ell$ is in the contact set ${\mathcal {C}}_{v_h}$ of $v_h$ if (and only if)

$$\begin{aligned} 0 \le v_h(z) - v_h(x) - \nabla v_h(x) \cdot (z - x) \quad \text {for all } z \in \varOmega \end{aligned}$$

(27)

(because $\partial \varGamma _{v_h} (x) = \{\nabla v_h(x)\}$ for any $x \in \varOmega \cap {\mathcal {C}}_{v_h}$ from the proof of Theorem 5). While this condition can be checked explicitly, it is a global condition and, therefore, leads to a global problem for each Gauss point, which may become rather expensive. Instead, (27) is verified at only a finite number of points, e.g., $z \in {\mathcal {V}}_\ell :={\mathcal {N}}_\ell \cup {\mathcal {N}}_\ell ^b$, where ${\mathcal {N}}_\ell ^b \subset \partial \varOmega$ is a discrete subset of $\partial \varOmega$. The set of points ${\mathcal {V}}_\ell$ create a quasi-uniform refinement ${\mathcal {T}}_\ell$ of the partition ${\mathcal {T}}$ into triangles and we assume that the mesh-size of ${\mathcal {T}}_\ell$ tends to zero as $\ell \rightarrow \infty$. Let $\textrm{I}_\ell v_h$ denote the nodal interpolation of $v_h$ w.r.t. the mesh ${\mathcal {T}}_\ell$. We replace the function $\chi _{{\mathcal {C}}_{v_h}}$ in (26) by the indicator function $\chi _{{\mathcal {C}}^\ell _{v_h}}$ of the set

$$\begin{aligned} {\mathcal {C}}^\ell _{v_h} :={\mathcal {C}}_{\textrm{I}_\ell v_h} \cap \{x \in \varOmega {\setminus } \cup {\mathcal {E}}: \textrm{D}^2 v_h(x) \ge 0 \text { is positive semi-definite}\}. \end{aligned}$$

Similar discrete contact sets were considered in [22] with application to the Monge–Ampère equation in [21]. In practice, the numerical integration formula for $\Vert f - f_h\Vert _{L^2(\omega )}$ reads

$$\begin{aligned} \sum _{T \in {\mathcal {T}}} \sum _{x \in {\mathcal {N}}_\ell \cap T \cap \omega } {\text {meas}}(T) w_{\ell ,T}(x)(f(x) - 2\chi _{{\mathcal {C}}^\ell _{v_h}}(x) (\det \textrm{D}^2 v_h(x))^{1/2})^2. \end{aligned}$$

(28)

The convex envelope $\varGamma _{\textrm{I}_\ell v_h}$ of $\textrm{I}_\ell v_h$ can be computed, for instance, by the Quickhull algorithm [2]. Therefore, it is straight-forward to compute (28). We note that if $x \in {\mathcal {C}}_{v_h} \cap {\mathcal {N}}_\ell$, then (27) holds for any $z \in {\mathcal {V}}_\ell$. Since the convex envelope of the continuous piecewise affine function $\textrm{I}_\ell v_h$ only depends on the nodal values of $v_h$, this implies $x \in {\mathcal {C}}^\ell _{v_h} \cap {\mathcal {N}}_\ell$. However, the reverse is not true. Hence, (28) and (26) may not coincide. From the uniform convergence of $\textrm{I}_\ell v_h$ to $v_h$ as $\ell \rightarrow \infty$, we deduce

$$\begin{aligned} \limsup _{\ell \rightarrow \infty } {\mathcal {C}}^\ell _{v_h} :=\cap _{\ell \in {\mathbb {N}}} \cup _{k \ge \ell } {\mathcal {C}}^\ell _{v_h} \subset {\mathcal {C}}_{v_h}, \end{aligned}$$

cf. [3, Lemma A.1]. Given any $\delta > 0$, this implies ${\mathcal {C}}^\ell _{v_h} {\setminus } {\mathcal {C}}_{v_h} \subset \{x \in \varOmega {\setminus } {\mathcal {C}}_{v_h} : \textrm{dist}(x,{\mathcal {C}}_{v_h}) \le \delta \}$ for sufficiently large $\ell$. Therefore, the set of all points $x \in {\mathcal {N}}_\ell$ with $\chi _{{\mathcal {C}}_{v_h}} \ne \chi _{{\mathcal {C}}^\ell _{v_h}}(x)$ is a subset of ${\mathcal {C}}^\ell _{v_h} {\setminus } {\mathcal {C}}_{v_h}$, whose Lebesgue measure vanishes in the limit as $\ell \rightarrow \infty$. In conclusion, the limits of (26) and (28) coincide.

Computation of $\mu :=\limsup _{x \rightarrow \partial \varOmega } |(g - \varGamma _{v_h})(x)|$. The boundary residual $\mu$ is approximated by $\Vert g - \varGamma _{\textrm{I}_\ell v_h}\Vert _{L^\infty (\partial \varOmega )}$. Since $\varGamma _{v_h} \le \textrm{I}_\ell v_h$ and $\textrm{I}_\ell v_h$ is piecewise affine, $\varGamma _{v_h} \le \varGamma _{\textrm{I}_\ell v_h}$ holds in $\varOmega$. On the other hand, we have $\lim _{\ell \rightarrow \infty } \Vert v_h - \textrm{I}_\ell v_h\Vert _{L^\infty (\varOmega )} = 0$. Hence, any supporting hyperplane $a_{x}$ of $\varGamma _{\textrm{I}_\ell v_h}$ at $x \in \varOmega$ satisfies $a_x - \delta _\ell \le v_h$ in $\varOmega$ with $\delta _\ell :=\Vert v_h - \textrm{I}_\ell v_h\Vert _{L^\infty (\varOmega )}$. Since $a_x - \delta _\ell$ is an affine function, $\varGamma _{\textrm{I}_\ell v_h}(x) - \delta _\ell = a_x(x) - \delta _\ell \le \varGamma _{v_h}(x)$. We conclude $\varGamma _{\textrm{I}_\ell v_h} - \delta _\ell \le \varGamma _{v_h} \le \varGamma _{\textrm{I}_\ell v_h}$ in $\varOmega$. In particular, $\lim _{\ell \rightarrow \infty } \Vert g - \varGamma _{\textrm{I}_\ell v_h}\Vert _{L^\infty (\partial \varOmega )} = \mu$.

Choice of $\varOmega '$. Let $\delta :=\min _{E \in {\mathcal {E}}} h_E$ denote the minimal edge length of the mesh ${\mathcal {T}}$. For all integers $0 \le j < 1/(2\delta )$, define $\varOmega _{j\delta } :=\{x \in \varOmega : \textrm{dist}(x,\partial \varOmega ) \ge j\delta \}$. It seems canonical to choose $\varOmega ' :=\varOmega _{j\delta }$, where j is the index that minimizes $\textrm{RHS}_0$. However, this choice may lead to significant computational effort. From the interior regularity of Alexandrov solutions [4], we can expect that the error is concentrated on the boundary and so, the best j will be close to one. Accordingly, the smallest $j \ge 0$ is chosen so that $\textrm{RHS}_0$ with $\varOmega ' :=\varOmega _{(j+1)\delta }$ is larger than $\textrm{RHS}_0$ with $\varOmega ' :=\varOmega _{j\delta }$.

5.1.3 Adaptive marking strategy

We define the refinement indicator

$$\begin{aligned} \eta (T) :=j\delta \sqrt{2}\Vert f-f_h\Vert _{L^2(T)}^2 + (1-2j\delta )^2\Vert f-f_h\Vert _{L^2(T \cap \varOmega _{j\delta })}^2 \end{aligned}$$

for any $T \in {\mathcal {T}}$, where the scaling in $\delta$ arises from (21) with $n = 2$. Let $\sigma :=\textrm{RHS}_0 - \mu$ denote the remaining contributions of $\textrm{RHS}_0$, where $\mu = \limsup _{x \rightarrow \partial \varOmega } |(g - \varGamma _{u_{h,\varepsilon }})(x)|$ from above. If $\sigma /10 < \Vert g - u_{h,\varepsilon }\Vert _{L^\infty (\partial \varOmega )}$, then we mark one fifth of all boundary edges $E \in {\mathcal {E}}$ with the largest contributions $\Vert g - u_{h,\varepsilon }\Vert _{L^\infty (E)}$. Otherwise, we mark a set ${\mathcal {M}}$ of rectangles with minimal cardinality so that

$$\begin{aligned} \frac{1}{2}\sum _{T \in {\mathcal {T}}} \eta (T) \le \sum _{T \in {\mathcal {M}}} \eta (T). \end{aligned}$$

5.1.4 Displayed quantities

The convergence history plots display the errors $\Vert u - u_{h,\varepsilon }\Vert _{L^\infty (\varOmega )}$, $\textrm{LHS} :=\Vert u - \varGamma _{u_{h,\varepsilon }}\Vert _{L^\infty (\varOmega )}$ as well as the error estimator $\textrm{RHS}_0$ against the number of degrees of freedom $\textrm{ndof}$ in a log–log plot. (Note that $\textrm{ndof}$ scales like $h_\mathrm {\max }^{-2}$ on uniformly refined meshes.) We utilize $\Vert u - \varGamma _{\textrm{I}_\ell u_{h,\varepsilon }}\Vert _{L^\infty (\varOmega )}$ as an approximation of $\Vert u - \varGamma _{u_{h,\varepsilon }}\Vert _{L^\infty (\varOmega )}$, where $\varGamma _{\textrm{I}_\ell u_{h,\varepsilon }}$ is computed by the Quickhull algorithm [2]. Whenever the solution u is sufficiently smooth, the errors $\Vert u - u_{h,\varepsilon }\Vert _{H^1(\varOmega )}$ and $\Vert u - u_{h,\varepsilon }\Vert _{H^2(\varOmega )}$ are also displayed. Solid lines in the convergence history plots indicate adaptive mesh-refinements, while dashed lines are associated with uniform mesh-refinements. The experiments are carried out for the regularization parameters $\varepsilon = 10^{-3}$ in the first two experiments and $\varepsilon = 10^{-4}$ for the third experiment. For a numerical comparison of various $\varepsilon$, we refer to [15].

5.2 Regular solution

In this example from [11], the exact solution u is given by

$$\begin{aligned} u(x) = \frac{(2|x|)^{3/2}}{3} \end{aligned}$$

with $f(x) = 1/|x|$. The solution belongs to $H^{5/2-\nu }(\varOmega )$ for any $\nu >0$, but not to $C^2(\overline{\varOmega })$. It is proven in [15] that u is the viscosity solution to $F_\varepsilon (f;x,\textrm{D}^2 u) = 0$ in $\varOmega$ for any regularization parameter $0 < \varepsilon \le 1/3$. Accordingly, we observed no visual differences in the convergence history plots for different $0 < \varepsilon \le 1/3$. Figure 1 displays the convergence rates 0.8 for $\Vert u - u_{h,\varepsilon }\Vert _{L^\infty (\varOmega )}$ and $\textrm{RHS}$, 3/4 for $\Vert u - u_{h,\varepsilon }\Vert _{H^1(\varOmega )}$, and 1/4 for $\Vert u-u_{h,\varepsilon }\Vert _{H^2(\varOmega )}$ on uniform meshes. The adaptive algorithm refines towards the singularity of u at 0 and leads to improved convergence rates for all displayed quantities. We observe the rate 1.75 for $\Vert u - u_{h,\varepsilon }\Vert _{L^\infty (\varOmega )}$, 1 for $\textrm{LHS}$, $\textrm{RHS}_0$, and $\Vert u - u_{h,\varepsilon }\Vert _{H^2(\varOmega )}$, and 1.5 for $\Vert u-u_{h,\varepsilon }\Vert _{H^1(\varOmega )}$. It is also worth noting that $\textrm{RHS}_0$ seems to be efficient on adaptive meshes.

Fig. 1

Convergence history for the first experiment with $\varepsilon = 10^{-3}$

5.3 Convex envelope of boundary data

In the second example, we approximate the exact solution

$$\begin{aligned} u(x,y) :=|x - 1/2| \end{aligned}$$

to $\det \textrm{D}^2 u = 0$ in $\varOmega$, which is the largest convex function with prescribed boundary data. The solution belongs to $H^{3/2-\delta }(\varOmega )$ for any $\delta > 0$, but not to $H^{3/2}(\varOmega )$. It was observed in [15] that the regularization error of $u-u_\varepsilon$ dominates the discretization error $u - u_{h,\varepsilon }$ on finer meshes. Therefore, the errors $\Vert u - u_{h,\varepsilon }\Vert _{L^\infty (\varOmega )}$ and $\Vert u - u_{h,\varepsilon }\Vert _{H^1(\varOmega )}$ stagnate at a certain value (depending on $\varepsilon$) as displayed in Fig. 2. However, $\textrm{LHS}$ converges with convergence rate 1/2 on uniform meshes even for fixed $\varepsilon$. At first glance on the discrete solution shown in Fig. 3, we can expect that the maximum of $|u - u_{h,\varepsilon }|$ is attained along the line $\textrm{conv}\{(1/2,0),(1/2,1)\}$. This error depends on the regularization parameter and only vanishes in the limit as $\varepsilon \rightarrow 0$, but the convex envelope of $u_{h,\varepsilon }$ provides an accurate approximation of u along this line. In fact, Fig. 4 shows that the adaptive algorithm refines towards the points (1/2, 0) and (1/2, 1), but the whole line $\textrm{conv}\{(1/2,0),(1/2,1)\}$ is only of minor interest. We observe the improved convergence rate 2.5 for $\textrm{LHS}$ on adaptive meshes. The guaranteed upper bound $\textrm{RHS}_0$ can provide an accurate estimate of $\textrm{LHS}$, but seems to oscillate due to the nature of the problem. The goal of the adaptive algorithm is the reduction of $\textrm{RHS}_0$, which consists of the error $\Vert f - f_h\Vert _{L^2(\varOmega )}$ in the Monge–Ampère measures and of some boundary data approximation error. Thanks to the additional regularization provided by the convex envelope, $\Vert f - f_h\Vert _{L^2(\varOmega )}$ is concentrated at the points (1/2, 0) and (1/2, 1), but becomes very small after some mesh-refining steps. We even observed in Fig. 2 that $\textrm{LHS} = \textrm{RHS}_0$ on two meshes, i.e., $\Vert f - f_h\Vert _{L^2(\varOmega )} = 0$. Then $\textrm{RHS}_0$ is dominated by the data boundary approximation error and leads to mesh refinements on the boundary. This may result in significant changes in the Monge–Ampère measure of $\varGamma _{u_{h,\varepsilon }}$, because the convex envelope of the discrete function $u_{h,\varepsilon }$ depends heavily on its values on the boundary in this class of problems.

Fig. 2

Convergence history for the second experiment with $\varepsilon = 10^{-4}$

Fig. 3

Discrete solution on a uniform mesh with 4225 nodes

Fig. 4

Adaptive mesh with 1907 nodes for the second experiment

5.4 Nonsmooth exact solution

In this example, the function

$$\begin{aligned} u(x,y) :=-\big (\sin (\pi x)^{-1} + \sin (\pi y)^{-1}\big )^{-1} \end{aligned}$$

is the solution to the Monge–Ampère equation (1) with homogenous boundary data and right-hand side

$$\begin{aligned} f(x,y) = \frac{4\pi ^2\sin (\pi x)^2\sin (\pi y)^2(2-\sin (\pi x)\sin (\pi y))}{(\sin (\pi x) + \sin (\pi y))^4}. \end{aligned}$$

The function u belongs to $C^2(\varOmega ) \cap H^{2-\delta }(\varOmega )$ for all $\delta > 0$, but neither to $H^{2}(\varOmega )$ nor $C^2(\overline{\varOmega })$. The convergence history is displayed in Fig. 5. Notice from Proposition 2 that $\textrm{RHS}_0$ consists solely of the error in the Monge–Ampère measures. In this example, f exhibits strong oscillations at the four corners of the domain $\varOmega$ and the adaptive algorithm seems to solely refine towards these corners as displayed in Fig. 6. While $\textrm{RHS}_0$ converges on uniform meshes (although with a slow rate), there is only a marginal reduction of $\textrm{RHS}_0$ for adaptive computation. We can conclude that the discrete approximation cannot resolve the infinitesimal oscillation of the Monge–Ampère measure of u properly. This results in the stagnation of $\Vert u - u_{h,\varepsilon }\Vert _{L^\infty (\varOmega )}$ and $\textrm{LHS}$ at an early level in comparison to uniform mesh refinements. However, we also observed that the stagnation point depends on the maximal mesh-size. In fact, if we start from an initial uniform mesh with a small mesh-size $h_0$, significant improvements of $\textrm{RHS}_0$ are obtained on the first levels of adaptive mesh refinements as displayed in Fig. 7. Undisplayed experiments show the same behaviour for $\Vert u - u_{h,\varepsilon }\Vert _{L^\infty (\varOmega )}$. This leads us to believe that, in this example, a combination of uniform and adaptive mesh-refining strategy provides the best results.

Fig. 5

Convergence history for the third experiment with $\varepsilon = 10^{-4}$

Fig. 6

Adaptive mesh with 1351 nodes for the third experiment

Fig. 7

Convergence history of $\textrm{LHS}$ for the third experiment with $\varepsilon = 10^{-4}$ and different initial meshes