Numerical Solution of the Simple Monge-Amp\`ere Equation with Non-convex Dirichlet Data on Non-convex Domains

TL;DR

This paper provides numerical evidence that the semi-Lagrangian method for the Monge-Ampère equation converges even on non-convex domains and with non-convex boundary data, extending known convergence results.

Contribution

It demonstrates convergence of the semi-Lagrangian method without convexity assumptions, supported by numerical experiments and analysis.

Findings

Numerical solutions converge on non-convex domains.

Multi-valued functions can approximate solutions satisfying boundary conditions.

Convergence observed without strict convexity assumptions.

Abstract

The existence of a unique numerical solution of the semi-Lagrangian method for the simple Monge-Amp\`ere equation is known independently of the convexity of the domain or Dirichlet boundary data -- when the Monge-Amp\`ere equation is posed as Bellman problem. However, the convergence to the viscosity solution has only been proved on strictly convex domains. In this paper we provide numerical evidence that convergence of numerical solutions is observed more generally without convexity assumptions. We illustrate how in the limit multi-valued functions may be approximated to satisfy the Dirichlet conditions on the boundary as well as local convexity in the interior of the domain.

Equations

$$\det (D^{2}u)$$

$$u(x)$$

$$H(A,f)$$

$$\displaystyle H\bigl{(}D^{2}u(x),f(x)\bigr{)}$$

$$u(x)$$

$$\{v\in C(\overline{\Omega }):\text{v is convex and a viscosity solution of eqrefMA}\}$$

$$=$$

$$B:D^2\varphi (x)$$

$$=\sum _i{\lambda }_i\frac{\varphi (x+k{e}_i)-2\varphi (x)+\varphi (x-k{e}_i)}{k^2}+\mathcal{O}(k^2).$$

$$\det ({\underline{u}}_{ij})$$

$$\underline{u}$$

$$\det (u_{ij})$$

$$u$$

$$\det (u_{ij})=0\text{in}\Omega ,u=\varphi \text{on}\partial \Omega$$

$$\Omega=\bigl{[}(0,1)\times(-1,1)\bigr{]}\cup\bigl{[}(-1,1)\times(0,1)\bigr{]}.$$

$$\Omega=B_{1}(0)\cup\bigr{[}(0,1)\times(0,1)\bigr{]}=\{x\in\mathbf{R}^{2}:\|x\|<1\}\cup\bigr{[}(0,1)\times(0,1)\bigr{]}$$

$${\nu }_v:\partial \Omega \to \mathbf{R},x\mapsto \inf \{h\in \mathbf{R}:(x,h)\in \overline{Co}({\Gamma }_v)\}$$

$${\nu }_v(x)\le g(x)\forall x\in \partial \Omega .$$

$$\Omega =$$

$$\displaystyle\,\cup\bigl{[}B_{1}((0,0))\cap\{x\in\mathbf{R}^{2}:x_{1}\leq 0\}\bigr{]},$$

Full text

\headlinetitle

Monge–Ampère equation with non-convex data and domains \lastnameoneJensen \firstnameoneMax \nameshortoneM. Jensen \addressoneDepartment of Mathematics, University of Sussex, Brighton BN1 9QF \countryoneEngland \emailonem.jensen@sussex.ac.uk

\researchsupported

1 Introduction

The paper is concerned with the numerical computation of solutions of the so-called simple Monge–Ampère equation

[TABLE]

where $\Omega\subset\mathbf{R}^{2}$, $g\in C(\mathbf{R}^{2})$ and $f\in C(\Omega)$ is non-negative: $f\geq 0$.

This raises immediately the question how the notion of solution for (1) should be defined. Beyond classical solutions, the generalisations in the Aleksandrov and viscosity sense provide settings of less smooth solutions. Any such definition in the literature imposes implicitly or explicitly convexity properties on $u$.

If $\Omega$ is non-convex but it is known that there exists a locally convex subsolution which which attains the boundary data pointwise, then the works [7, 8] give criteria for the existence of a unique solution. If $\Omega$ is convex (possibly not strictly) and $g$ non-convex then [1] examines in the context of the closely related Gauss curvature problem the possibility of solutions which are in a certain sense multi-valued on $\partial\Omega$. We shall return to these results in sections 3 and 4.

To our knowledge there has been no systematic analysis of the well-posedness of (1) for the combination of a non-convex domain and general non-convex boundary data. It is therefore interesting if numerical methods can provide an insight into (1) in this case. The basis for our study is [6] where a semi-Lagrangian numerical scheme is proposed, which uses Krylov’s Bellman reformulation of (1). While the convergence proof to viscosity solutions in [6] requires strict convexity of $\Omega$ in order to make a comparison principle available, we highlight two results in that work which do not impose any form of convexity on $\Omega$ or $g$:

1. (a)

For any finite element mesh there exists a unique numerical solution $u_{i}$, which is the limit of a globally converging semi-smooth Newton method. 2. (b)

There exists a constant $C$ (depending on $\Omega$, $f$ and $g$ only) such that $\|u_{i}\|_{L^{\infty}}\leq C$ for all numerical solutions $u_{i}$ independently of the mesh.

Thus at this point we know that as the finite element mesh size $h$ approaches [math] some subsequences of numerical solutions converge weakly∗ to functions in $L^{\infty}$, which may be candidate solutions of (1) of some form.

The aims of this paper are

1. (a)

to provide numerical evidence that indeed not just subsequences but the whole sequence of numerical solutions converges, provided that the stencil size is scaled appropriately with the mesh size; 2. (b)

to highlight that the Bakelman interpretation of solutions might extend to the non-convex setting, consistently with the findings in [7, 8]; 3. (c)

to study the performance of the numerical scheme in the non-convex setting, e.g. the robustness of the semi-smooth Newton solver; 4. (d)

to present computations on general domains, e.g. which are not Lipschitz.

2 The Bellman formulation and its semi-Lagrangian approximation

A difficulty when solving (1) is that the Monge–Ampère operator $v\mapsto\det(D^{2}u)$ is only elliptic on the set of convex functions. It is therefore convenient to reformulate the problem in such a way that the set of solutions remains unchanged but ellipticity is established on the whole function space.

For that purpose we define the Bellman operator

[TABLE]

where $\mathbf{S}$ is the set of symmetric $2\times 2$ matrices, $\mathbf{S}_{+}:=\{A\in\mathbf{S};\,A\geq 0\}$ and $\mathbf{S}_{1}:=\{B\in\mathbf{S}_{+};\,\operatorname{tr}B=1\}$, and consider the boundary value problem

[TABLE]

It was shown [6] that

[TABLE]

This statement does not require boundedness or convexity of $\Omega$ and does in this form not refer to the boundary conditions.

To discretize (3a) we write $B=\sum_{i}\lambda_{i}\,e_{i}\,e_{i}^{\sf T}\in\mathbf{S}_{1}$ with the eigenvalues $\lambda_{i}$ and normalised eigenvectors $e_{i}$, so that for smooth functions $\phi$

[TABLE]

Evaluating these central differences on a P1 finite element space at interior nodes, combined with nodal interpolation of $g$ on $\partial\Omega$, gives the numerical scheme of [6].

Remark 2.1**.**

A comparison principle on the set of semi-continuous functions in respect to $H$ holds when classical boundary conditions are considered [6], where $\Omega$ may be non-convex. For the proof of convergence this comparison principle is applied to the upper and lower semi-continuous envelopes of the sequence of numerical solutions—in the setting of strictly convex domains where it is ensured that the envelopes attain the boundary conditions in the classical sense.

The numerical experiments in this paper highlight that on non-convex domains the envelopes may only satisfy the boundary conditions in a generalised form. It would be most convenient if the convergence proof could be translated to viscosity boundary conditions as they are defined either in [2] or alternatively in the form of [4]. However, in [12] counterexamples were given to show that comparison principles with these kinds of viscosity boundary conditions do in general not hold on the spaces of semi-continuous functions (even if $\Omega$ is convex), screening out trivial extensions of the proof.

3 The Monge–Ampère equation on non-convex domains

In the light of Remark 2.1, we wish to identify settings in which the exact solution of (1) satisfies the boundary conditions classically. If was shown in [7] that this is guaranteed provided one knows of the existence of a strict subsolution and sufficient smoothness of data and domain. Indeed [7] covers a more general equation with dependence on first-order derivatives.

Theorem 3.1** (Guan and Spruck ‘93).**

Let $\Omega\subset\mathbf{R}^{n}$ be a smooth bounded domain with boundary components $\partial\Omega=(\gamma_{1},\ldots,\gamma_{m})$. Assume that there is a smooth, strictly locally convex function $\underline{u}$ in $\overline{\Omega}$ satisfying

[TABLE]

where $\phi$ and $\psi$ are smooth, $\psi>0$ and $\psi^{1/n}(x,u,p)$ is convex in $p$. Then there is a smooth locally convex solution $u\in C^{\infty}(\overline{\Omega})$ to the Monge–Ampère boundary-value problem:

[TABLE]

If $\psi_{u}\geq 0$, then the solution is unique.

A similar result can be found in [8], with vanishing right-hand side but a more precise specification of the required regularity.

Theorem 3.2** (Guan ‘98).**

Assume that $\partial\Omega$ is in $C^{3,1}$ and $\phi\in C^{3,1}(\partial\Omega)$. Suppose there exists a locally strictly convex function $\underline{u}\in C^{2}(\overline{\Omega})$ with $\underline{u}=\phi$ on $\partial\Omega$. There there is a unique locally convex weak solution of

[TABLE]

in $C^{1,1}(\overline{\Omega})$.

Our first computational experiment concerns the convergence of the numerical approximations to an exact solution, for which we know from the above that the boundary conditions are admitted classically.

Experiment 3.1** (Quartic problem on L-shape).**

We approximate the exact solution $u(x)=(x_{1}^{2}+x_{2}^{2})^{2}$ on the L-shape:

[TABLE]

The quasi-uniform grid has at the coarsest level $58$ nodes and at the finest level after $6$ uniform refinements $168,961$ nodes. The stencil diameter $k$ is, away from the boundary, represented through $k=m\cdot h$ by a fixed positive factor $m$ and the (average) mesh size $h$. Near the boundary, so where $m\cdot h$ is larger than the distance to $\partial\Omega$, the stencil is reduced in size to remain within $\Omega$, cf. [6].

Figure 1 shows the decay of the relative $L^{\infty}(\Omega)$ approximation error for different choices of the multiplier $m$. Importantly, this gives numerical evidence of the convergence of the scheme on a non-convex domain. The table shows the multiplier $m$ which achieves the smallest relative $L^{\infty}$ error for a given number of degrees of freedom, and the number of Newton iterations of the respective computations to obtain a Newton step size of less than 5e-8 in the $\infty$-norm:

[TABLE]

Overall the largest number of Newton iterations in this computational experiment is $9$, which occurs on the finest mesh with $m=32$.

4 The simple Monge–Ampère equation with non-convex boundary data

In computational experiments where the domain $\Omega$ is not strictly convex and there is no subsolution $\underline{u}$ as in the above Theorems 3.1 and 3.2 we observe numerical solutions which appear to approximate multi-valued boundary data.

Experiment 4.1** (Non-convex boundary data on convex domain).**

Now

[TABLE]

is a convex but not strictly convex domain. On the vertical straight boundary segment \Gamma:=\partial\Omega\cap\bigl{[}\{1\}\times(0,1)\bigr{]} we impose the non-convex boundary data $g(x)=x_{2}(1-x_{2})$ while on the remainder $\partial\Omega\setminus\Gamma$ we set $g(x)=x_{1}^{4}-1$. We select $f=0$ on $\Omega$, implying that the graph of the solution $u$ is a surface of vanishing Gauss curvature.

A numerical solution is shown in Figure 2. In the vicinity of $\Gamma$ on the far side of the plot the numerical solution is nearly vertical, interpolating on $\Gamma$ the data $g(x)=x_{1}(1-x_{1})$ while attaining at the interior nodes neighbouring $\Gamma$ approximately the value $x_{1}^{4}-1\approx 0$.

Recalling from (4) that the Bellman formulation enforces convexity on the domain $\Omega$ of the differential operator, it appears that the boundary data is only extended into the interior of the domain in as far as convexity permits. Indeed, according Figure 3 the numerical solutions converge on the subdomain $\omega:=\{x\in\Omega:x_{1}<0.95\}$ under mesh refinement to $u(x)=x_{1}^{4}-1$; we see a relative $L^{\infty}(\omega)$ error of $1.35\cdot 10^{-4}$ on a mesh with $304129$ DoFs and $m=64$.

The table gives the numerical values of the error for the best choices of $m$ on a given mesh, and the associated number of Newton iterations:

[TABLE]

Overall the largest number of Newton iterations in this computational experiment is $19$, which occurs on the finest mesh with $m=64$.

To give an interpretation of Experiment 4.1 we review a result due to Bakelman [1]. Let $v$ be a bounded convex function on $\Omega$ and $\overline{Co}(\Gamma_{v})$ be the closed convex hull of the graph $\Gamma_{v}\subset\mathbf{R}^{2+1}$ of $v$. Then the function

[TABLE]

is called the border of the function $u$. We denote by $\mathcal{L}$ the set of functions $v$ which satisfy the Monge-Ampère differential equation in the Aleksandrov sense and

[TABLE]

If the domain $\Omega$ is strictly convex then Aleksandrov and viscosity solutions coincide for the simple Monge-Ampère equation: See [9] for the proof when the continuous $f$ is positive. This argument can be extended to the case of non-negative $f$ with the tools of [5]—alternatively one can use the uniqueness of solutions and [6].

In the context of the Gauss curvature problem, Bakelman shows under a data condition that on bounded strictly convex domains $\Omega$ the set $\mathcal{L}$ is non-empty. Moreover, there is a unique $u\in\mathcal{L}$ such that $\nu_{u}\geq\nu_{v}$ for all $v\in\mathcal{L}$. This $u$ is considered to be the generalised solution of the boundary value problem.

The numerical solutions of Experiment 4.1 evidently converge to $u(x)=x_{1}^{4}-1$ in $L^{\infty}_{\rm loc}(\Omega)$, so on relatively compact subsets of $\Omega$. The border of $u$ is its trace. As the numerical solution attains the continuum of values between $g(x)$ and $\nu_{u}(x)=u(x)$ in the vicinity of any $x\in\partial\Omega$, the visual impression is that of a convergence to a multi-valued limit which attains the interval $[\nu_{u}(x),g(x)]$ at $x$. It is therefore appealing to adopt this multi-valued interpretation for the purposes of this text.

Starting from [1], Froese recently studied a wide-stencil method for the Gauss curvature problem with non-classical boundary conditions on convex domains [10]. There, besides the definition with the border function $\nu_{u}$, also viscosity boundary conditions as in [4] are discussed.

5 The simple Monge–Ampère equation on non-convex domains without a classical subsolution

As far as we are aware there is no systematic analysis in the literature which extends [1, 7, 8] to boundary value problems where non-convex boundary data is imposed on a non-convex domain.

Yet, also in this setting the numerical scheme of [6] is guaranteed to have a unique solution for any mesh. From our point of view this makes numerical experiments interesting.

Indeed, equations of Monge–Ampère type have been proposed for physical and biological models where the application does not justify to impose convexity on $\Omega$ or $g$. An example is the Monge–Ampère Keller–Segel system of chemotaxis [3, 11]. There the solution of a Monge–Ampère equation describes the density of a chemical substance, whose physical domain must not necessarily be convex.

Experiment 5.1** (Non-convex boundary data on the L-shape).**

Combining in spirit the domain of Experiment 3.1 with the boundary data of Experiment 4.1, we set $g(x)=x_{2}(1-x_{2})$ on left boundary segment \Gamma:=\partial\Omega\cap\bigl{[}\{-1\}\times(0,1)\bigr{]} of the L-shape \Omega=\bigl{[}(0,1)\times(-1,1)\bigr{]}\cup\bigl{[}(-1,1)\times(0,1)\bigr{]}. On the remainder $\partial\Omega\setminus\Gamma$ we set $g(x)=x_{1}^{4}-1$ and on $\Omega$ we select $f=0$.

Similarly to the previous experiment, numerical solutions approximate a multi-valued solution in the vicinity of $\Gamma$, while on the subdomain $\omega:=\{x\in\Omega:x_{1}>-0.95\}$ the numerical solutions converge under mesh refinement to $u(x)=x_{1}^{4}-1$. We see a relative $L^{\infty}(\omega)$ error of $4.8\cdot 10^{-5}$ on a mesh with $168961$ DoFs and $m=64$, cf. Figure 4. Overall the largest number of Newton iterations in this computational experiment is $25$, which occurs on the finest mesh with $m=32$.

The multi-valued behaviour may also be exhibited by problems with convex or even vanishing boundary data, when no smooth classical subsolution exists.

Experiment 5.2** (Hölder domain).**

We now consider the heart-shaped domain

[TABLE]

where again $B_{\delta}(x)$ denotes the ball with radius $\delta$ and centre $x$. This domain does not satisfy a Lipschitz condition at the origin, as can also be seen in Figure 5 (b). In this example the boundary data $g=0$ vanishes, while $f=1$.

The plots of Figure 5 give the impression that at the origin the numerical solution approximates a multi-valued function but then transitions to a classical single-valued boundary condition. The Newton method stopped the computation with 166,176 DoFs after 16 iterations, when a Newton step size of $5\cdot 10^{-8}$ was reached.

A note about the stencil diameter: As the cartoon stencil in Figure 5 (b) illustrates, a too simple implementation a wide second-order central difference with the centre on the right of the re-entrant boundary could have a node in $\Omega$ which is on the left of the re-entrant boundary. On non-convex domains stencils should be scaled so that not just the nodes of the stencil belong to $\Omega$, but also any points between them.

Experiment 5.3** (Bent square).**

In this last experiment we solve the degenerate equation on a sequence of curved polyhedra, starting from a uniformly convex geometry and leading to polyhedra with concave faces:

To be precise: The first three domains are formed from the intersection of four circles with the midpoints $(c,0)$, $(0,c)$, $(-c,0)$, $(0,-c)$ where $c$ is $3$, $10$, $100$ respectively and where the radii are chosen so that the circles intersect at the points $(1,1)$, $(-1,1)$, $(1,-1)$ and $(-1,-1)$. The fourth domain is the square centred at the origin with side length $2$. The remaining three domains are obtained by subtracting from this square the points which lie in circles with the same centres and intersection points as above, but smaller radii.

With $f=0$ the graphs of the exact solutions are (when restricted to the interior) surfaces of vanishing Gauss curvature. The function $g$ is equal to the convex function $x_{2}^{2}-1$ where $x_{1}>x_{2}$ and equal to the concave function $1-x_{1}^{2}$ where $x_{1}\leq x_{2}$. We note that $g$ is continuous on $\partial\Omega$. Two of the resulting solutions are plotted in Figure 6. In the interior the numerical solutions do not significantly vary from the function $x_{2}^{2}-1$, which already defined the convex part of the boundary data. This is notable because on strictly convex domains the boundary function $g$ is attained on all of $\partial\Omega$ in the classical single-valued sense, while on the other domains the multi-valued cut-off mechanisms appears to screen out concave sections of $g$. Indeed on $\omega=(-3/4,3/4)^{2}$ we find the following relative $L^{\infty}$ errors, taking $x_{2}^{2}-1$ as the reference solution:

[TABLE]

The table also shows an increase in the required Newton iterations to achieve a Newton step size of $5\cdot 10^{-8}$ when the faces of $\Omega$ are concave. The number of DoFs vary in this experiment between 108,124 and 140,324, depending on the domain.