1 Introduction↩︎

The Monge–Ampère equation is a fully non-linear partial differential equation, originating from the study of surfaces with prescribed curvature which goes back to Monge in the late 18th and Ampère in the early 19th century. Since then, variations of this equation have found a number of applications, for example in optimal transport [1], [2] and Riemannian geometry (the local isometric embedding problem [3]). The equation is also of great theoretical interest and is often used as a model example of a fully non-linear PDE. The complex version of the Monge–Ampère equation plays a central role in pluripotential theory, i.e.in the study of plurisubharmonic functions on domains in \(\mathbb{C}^n\) and more recently on Kähler manifolds [4].

In this paper, we introduce a novel numerical method for solving the Dirichlet problem for the real Monge–Ampère equation, based on Bellman’s principle: representing the non-linear operator as an infimum of elliptic linear partial differential operators. We construct a numeric iterative scheme for approximating the operator attaining this infimum. Under the assumption that our algorithm produces a convex approximant, we show monotonicity. More precisely, we show that the convex approximants always majorise the solution to the original Monge–Ampère equation.

Numerical experiments illustrate that the method works well in settings where the solution is smooth and strictly convex. In these cases, it outperforms recently published methods by a factor of 3–10 in execution time. For mildly singular equations, the speedup is even better: in some cases improving the running time by a factor of 100 or more. In the fully degenerate case, when the right hand side of the equation vanishes on a non-empty open set, our proposed method does not converge in general. We intend to return to the fully degenerate problem in a subsequent paper.

2 Background↩︎

The Dirichlet problem for the Monge–Ampère operator can be formulated as \[\label{eq:MA} \begin{cases} \det(D^2u(x)) = f(x), & x \in \Omega \\ u(x) = \phi(x) , & x \in \partial\Omega \end{cases}\tag{1}\] where \(\Omega \subset \mathbb{R}^n\) is a convex domain, \(f \ge 0\) and the boundary values of \(u\) are suitably interpreted. Here, \(D^2 u = (u_{jk})\) denotes the Hessian of \(u\). By restricting our attention to convex solutions of 1 , the equation is (degenerate) elliptic, and we will only consider such solutions. If we impose some extra conditions on \(\Omega\), \(f\) and \(\phi\), we can guarantee that the equation has a unique solution. Note that the left hand side, the Monge–Ampère operator, \(\mathop{\mathrm{MA}}(u) := \det(D^2u(x))\) makes perfect sense for \(C^2\) functions, and it can be interpreted pointwise also for functions in \(W_{\text{loc}}^{2,n}\). For functions of lower regularity, the Monge–Ampère operator must be viewed in a generalized sense, and since the operator is fully non-linear, there is no direct interpretation in the sense of distributions. However, convex functions are automatically locally Lipschitz and it turns out that we can define the Monge–Ampère measure (in the sense of Alexandrov) of any convex function \(u\), by \[\mu_u(E) = \mathop{\mathrm{Leb}}(\partial u(E))\] for any Borel subset \(E \subset \Omega\), where \(\mathop{\mathrm{Leb}}\) denotes the Lebesgue measure, \(\partial u(E) = \cup_{x\in E} \partial u(x)\) and \(\partial u(x)\) is the set of all slopes of the supporting hyperplanes to the graph of \(u\) at \(x\). If \(u \in C^2\) and \(D^2 u\) is positive semi-definite, then \(\mu_u(E) = \det(D^2 u)\). We should mention, that in contrast to the real setting, it is impossible to define the complex Monge–Ampère operator for every plurisubharmonic function. We refer to the monograph by Guedj and Zeriahi [4] for the theory of complex Monge–Ampère equations.

There is a vast number of analytic results in the literature concerning existence and regularity of solutions to 1 . A complete overview is outside the scope of this paper, but for example, it follows from results by Alexandrov [5] and Bakelman [6] from the 1950’s that if \(\Omega\) is assumed to be a bounded strictly convex domain, \(f\) is a finite Borel measure on \(\Omega\), and \(\phi \in C(\partial\Omega)\), then 1 has a unique convex solution \(u\) (in the sense of Alexandrov) with \(u \in C(\bar\Omega)\).

Furthermore, if \(\Omega\) is strictly convex (or uniformly convex in the unbounded setting) with \(C^{k+2, \alpha}\) boundary, \(0 < \phi \in C^{k+2, \alpha}(\partial\Omega)\) and \(f \in C^{k,\alpha}(\bar\Omega)\) for some integer \(k \ge 2\) and some \(0 < \alpha < 1\), then the Dirichlet problem 1 has a unique convex solution \(u \in C^{k+2, \alpha}(\bar\Omega)\), see [7]–[9].

If we merely assume that \(\Omega\) is convex, not necessarily strictly convex, the situation is a little less straight-forward, and we need to assume that \(\phi\) is convex on all line segments in \(\partial\Omega\) to guarantee existence of solutions to 1 . Similarly, if \(f\) is only assumed to be non-negative instead of strictly positive, the solutions to the Dirichlet problem will generally be less smooth, even in strictly convex domains. If \(f\) has zeros, the equation is called degenerate. See for example Błocki [10] for results in this direction, where \(C^{1,1}\)-regularity of the solution is shown under suitable technical assumptions. We will come back to this phenomenon in example 5.2.3.

For an accessible and thorough overview of the Alexandrov solutions as well as results on classical solutions including the ones mentioned, we refer to the recent monographs by Figalli [11] and Le [12].

There have been a number of suggested numerical methods for solving the Dirichlet problem for the Monge–Ampère operator. As far as we are aware, the first approach was due to Oliker and Prussner [13] who proposed a direct discretization of the Alexandrov interpretation of the Monge–Ampère operator. Since then, a number of other approaches to solving Monge–Ampère equations have been tried: finite difference discretization, often using wide stencils to obtain stable and consistent methods [14]–[18], variational methods [19], [20], fixpoint methods [14], [21] as well as finite elements [22], [23] only to mention a few.

3 Bellman’s principle↩︎

Our method relies on the following linear algebra result, sometimes called Bellman’s principle. In the setting of complex Monge–Ampère equations, this idea was introduced by Gaveau [24] and has later found several theoretical applications. We will formulate the result in the real setting and for convenience, we provide a proof of Bellman’s principle, similar to the approach taken by Wiklund [25].

Using Bellman’s principle, we can describe the Monge–Ampère operator (of a sufficiently smooth function) as an infimum over a class of elliptic partial differential operators. More precisely:

Note that, for every choice of \(\mathcal{B}\)-valued function \(B\), the operator \[u \mapsto \mathop{\mathrm{tr}}\big(B(x)D^2u(x)\big)\] is a linear elliptic partial differential operator, and if we somehow know the associated \(B(x)\) corresponding to the solution \(u\) of the Dirichlet problem 1 , then we can recover \(u\) by solving the linear elliptic Dirichlet problem \[\label{eq:linear95eq} \begin{cases} \mathop{\mathrm{tr}}(B(x)\,D^2u(x)) = n\,\big(f(x)\big)^{1/n}, & x \in \Omega \\ u(x) = \phi(x) , & x \in \partial\Omega, \end{cases}\tag{3}\] for which there are highly efficient numerical algorithms.

While we’re unable to prove that our method (as described in Section 4) converges, we remark that we do have a monotonicity result. More precisely,

Unfortunately, it seems like there are no known effective characterizations of when the solution to Equation  3 is convex, even in the case when \(B \equiv I\), i.e. for Poisson’s equation. Our numerical experiments suggest that our proposed algorithm produces a decreasing sequence of approximants converging to the solution of 1 , at least once one of the approximants is convex. So far, we have not been able to prove this.

We should point out that there is no proof of convergence for most methods that have been proposed for the Monge–Ampère equation. In particular, this is the case for the ones we compare our method to in Section 5. Notable exceptions are the direct discretization of the Alexandrov construction proposed in [13] which was shown to converge. Wide stencil finite difference discretizations [18] are also known to converge (to a viscocity solution).

In the remainder of the paper, we will restrict our attention to the two-dimensional case, \(n=2\), but we expect the method to be competitive also in higher dimensions.

4 Description of our algorithm↩︎

The above observation to Corollary 1 serves as the main underlying idea of the algorithm presented below, which we call the Bellman algorithm.

First, we discretize the domain \(\Omega\) using a rectangular \(N \times N\) mesh. For simplicity of the implementation, we only consider the case where \(\Omega\) is a bounded rectangle in \(\mathbb{R}^2\), despite the fact that these domains are not strictly convex, and the theoretical framework for solving Monge–Ampère equations on rectangles is less well-developed.

Next, we construct a sequence of functions \(u_k\) where \(u_k\) solves the discretized version of the following Dirichlet problem, using a narrow stencil second order finite difference solver: \[\label{eq:32trace32equation} \begin{cases} \mathop{\mathrm{tr}}(B_{k}(x)D^{2}u_{k}(x)) = 2\sqrt{f(x)}, & x \in \Omega \\ u_k(x) = \phi(x) , & x \in \partial\Omega. \end{cases}\tag{4}\] Here \(B_k\) is a sequence of \(\mathcal{B}\)-valued functions constructed as follows:

1. \(u_{0}\) is obtained by choosing \(B_{0}=I\), i.e.identity matrix. We observe that this renders equation 4 a Dirichlet problem for the Poisson equation.

2. for \(k\geq 1\) we loop over the mesh and construct the \(B_{k}\) matrices as follows:

   If the computed discrete Hessian \(D^2u_{k-1}(x)\) is positive definite, we apply the algorithm from [25]. Namely

   1. Diagonalise \(D^{2}u_{k-1}(x)=T_{k-1}(x)\,\Lambda_{k-1}(x)\,T^{T}_{k-1}(x)\),

   2. Construct \(A_{k-1}(x)\mathrel{\vcenter{:}}= \sqrt{\det(D^{2}u_{k-1})(x)}\,\Lambda^{-1}_{k-1}(x)\),

   3. Let \(B_{k}(x)\mathrel{\vcenter{:}}= T_{k-1}(x)\,A_{k-1}(x)\,T_{k-1}^{T}(x)\).

   If the computed \(D^2u_{k-1}\) is not positive definite, we mark the corresponding point \(x\), and temporarily take \(B_k(x)\) as the identity matrix.

3. We loop over the mesh a second time and replace \(B_k(x)\) at the marked points with the (normalized to determinant 1) convex combination of the computed Bellman matrices corresponding to the nearest grid points in horizontal and vertical direction for which the computed Hessian was positive definite.

We terminate the computation when \(\|u_k - u_{k-1}\|_{\infty} < 10^{-12}\) or when the number of iterations exceeds 10 000.

If the solution of the original Monge–Ampère equation is strictly convex, we expect the constructed functions \(u_k\) to be strictly convex as well, at least when \(k\) sufficiently large, but it can happen that convexity fails in the first few steps – see Example 5.1.2 – even when the final solution is strictly convex. The interpolation step [step:interpolation] is necessary to assure convergence and in practice it pushes the approximants \(u_k\) to become more and more convex as \(k\) increases. We will illustrate this phenomenon more clearly in Section 5.

If the solution of the original Monge–Ampère equation is not strictly convex, i.e. if \(f\) has zeros, we cannot expect our method to converge. In practice, our proposed algorithm still performs very well when the set of degeneracy is small, for example if \(f\) has an isolated zero, or even if the zero-set of \(f\) is one-dimensional (cf. Example 5.2.2).

In the fully degenerate case, where \(f\) vanishes on a non-empty open set or even vanishes identically, our proposed method fails to converge, or at least converges very slowly (cf. Example 5.3.1). We intend to return to the degenerate Monge–Ampère equations in a subsequent paper.

In our implementation of the Bellman method, we use the Python package FinDiff [26] to perform the necessary computations. Also, for efficiency reasons, the computations in step ([step:32B95matrices]) are performed in parallel, utilizing NumPy’s tensor operations, thus removing the need for an explicit loop over the mesh elements.

For comparison, we also implemented two methods (M1 and M2) proposed by Benamou et al. [14].

The M2 method is a fixed-point method which in spirit is somewhat similar to our method. The idea is to solve a sequence of Dirichlet problems for the Poisson equation, and update the right hand side in each iteration, whereas we instead change the partial differential operator at each step, keeping the right hand side constant. For detailed descriptions of these methods, see [14], but in short, for M2 we solve \[\label{eq:M2} \begin{cases} \Delta u_k = g_k(x), & x \in \Omega \\ u_k(x) = \phi(x), & x \in \partial\Omega. \end{cases}\tag{5}\] where \(g_0 = 2\sqrt{f}\) and \(g_{k+1} = T(g_k) = \sqrt{(\Delta u_k)^2 + 2(f-\det(D^2 u_k))}\).

It’s not too hard to show that the solution to the original Monge–Ampère equation 1 is obtained when \(g\) is a fixed point to the operator  \(T\). Note that our proposed algorithm can also be viewed as a fixed-point method: If \(u\) is a strictly convex solution of the Monge–Ampère equation 1 , and \(B(x)\) is the associated Bellman function, the solution of 4 reproduces \(u\).

Method M1 is different. Using a central difference discretization of the second order derivatives of \(u\), the discretized version of the Monge–Ampère equation \[(\mathcal{D}^2_{xx} u_{ij})(\mathcal{D}^2_{yy} u_{ij}) - (\mathcal{D}^2_{xy} u_{ij})^2 = f_{ij}\] is rewritten as a quadratic equation for \(u_{ij}\) as \[\label{eq:M1} 4(a_1 - u_{ij})(a_2 - u_{ij}) - \frac{1}{4}(a_3-a_4)^2 = h^4 f_{ij},\tag{6}\] where \[\begin{cases} a_1 = \frac{1}{2} ( u_{i+1,j} + u_{i-1,j} ) \\ a_2 = \frac{1}{2} ( u_{i,j+1} + u_{i,j-1} ) \\ a_3 = \frac{1}{2} ( u_{i+1,j+1} + u_{i-1,j-1} ) \\ a_4 = \frac{1}{2} ( u_{i-1,j+1} + u_{i+1,j-1} ). \end{cases}\] By solving 6 for \(u_{ij}\), selecting the smaller solution of the quadratic equation to impose local convexity, we have an iterative scheme for M1: \[\label{eq:M195iterative} u_{ij} = \frac{1}{2}\left( a_1 + a_2 + \sqrt{(a_1-a_2)^2 + \frac{1}{4}(a_3-a_4)^2 + h^4 f_{ij}} \right),\tag{7}\] which is then solved using Gauss–Seidel iteration, keeping the boundary values fixed and updating one \(u_{ij}\) at the time, using 7 . We should note that this is computationally expensive, especially in our NumPy implementation where explicit loops are slow, so M1, and to some extent our proposed method, are the ones that would benefit the most from an optimized implementation.

In our comparison we chose the same termination criterion for all methods, i.e.to terminate the computation when \(\|u_k - u_{k-1}\|_{\infty} < 10^{-12}\) or the number of iterations exceeds 10000 (for M1, we allow up to 300000 iterations, where one iteration means updating each \(u_{ij}\) once).

5 Numerical experiments↩︎

In this section, we show the results of testing the Bellman algorithm on several examples with various degrees of smoothness and regularity of the right-hand side of the Monge–Ampère equation, \(f\), on the square \(\Omega = [-1,1]^{2}\), unless otherwise stated. Moreover, we reproduce the algorithms M1 and M2 presented in [14].

Then we compare the performance of both algorithms in terms of the approximation of the exact solution (restricted to the mesh) \(\tilde{u}\) with the numerical solution \(u\), i.e.\(\|u-\tilde{u}\|_{\infty}\), \(\|u-\tilde{u}\|_{2}\), the number of iterations and CPU time needed for convergence. We also demonstrate the effectiveness of interpolating the Bellman functions, i.e. \(B_{k}(x)\) from Equation 4 , in our new algorithm on Example 5.1.2.

5.1 Smooth, strictly convex examples↩︎

5.1.1 Standard example↩︎

The first example that we will look at has the exact solution \[\label{eq:32standard95example} u(x,y)= e^{\frac{1}{2}(x^{2} +y^{2})}\tag{8}\] of the Monge–Ampère equation with a right-hand side given by \[\mathop{\mathrm{MA}}(u) = f(x,y)= (1+x^{2} +y^{2})e^{x^{2} +y^{2}}.\] This example has been used in several previous papers, see for example [14], [17], [21]. Clearly, \(u\) is \(C^\infty\)-smooth and strictly convex, so we expect any reasonable solver to give good results for this example.

In fact, the Bellman method produces a sequence of strictly convex approximants, and the interpolation step is never invoked. Comparing the Bellman method to the other methods, all approaches converge to the same solution regardless of mesh size, see Figure 1 (a). The Bellman method needs 6–7 iterations to converge, whereas M2 method requires 40–50 iterations, independent of mesh size. Even if each iteration in our method is computationally a little more demanding than in the M2 method, the overall running time is 5–6 times faster for the Bellman method, see Figure 1 (b). Both methods have similar running time complexity (\(\sim N^{2.8}\) and \(\sim N^{2.9}\), respectively).

For the M1 method, each iteration is very fast, but the number of iterations needed for convergence grows about linearly with \(N\), and 1500–90000 iterations are needed which results in slow total running times (\(\sim N^{4.0}\)). For a \(255 \times 255\) grid, the M1 method ran in just under 3 h, where the Bellman method took 24 s and the M2 method around 2 min.

Error in supremum (top) and \(\ell^{2}\) norm (bottom). The results for all three algorithms are almost identical and all graphs overlap. The results suggest that \(\|u-u_N\| \sim N^{-2}\) for the sup-norm as well as the \(\ell^2\)-norm

Running time (s). The results suggest that the running time is \(\sim N^{2.8}\) for the Bellman method, \(\sim N^{2.9}\) for M2 method and \(\sim N^{4.0}\) for M1

5.1.2 Degenerate example with regularisation↩︎

Error in supremum (top) and \(\ell^{2}\) norm (bottom). The results for the three algorithms are almost identical and suggest that \(\|u-u_N\| \sim N^{-2}\), for both norms. The Bellman algorithm without interpolation (Bellman2 in the graph) does not converge

Running time (s). The results suggest that the running time is \(\sim N^{2.7}\) for the Bellman method, \(\sim N^{2.8}\) for the M2 method and \(\sim N^{4.0}\) for M1

Iteration 1

Iteration 2

Iteration 3

Iteration 4

For this example, the exact solution to the Monge–Ampère equation is \[\label{eq:32degenerate95example95with95regularisation} u(x,y)=0.5(x-0.5)^{4} + \epsilon x^{2}+y^{2},\tag{9}\] where \(\epsilon \in \mathbb{R}_{+}\), and the right-hand side of the equation is \[\; \mathop{\mathrm{MA}}(u) = f(x,y)= 12(x-0.5)^{2}+4\epsilon.\] We note that the solution is smooth and strictly convex for each \(\epsilon > 0\), but the smaller the choice of \(\epsilon\), the less convex the exact solution is along the line \(x=0.5\). We solve the Monge–Ampère equation for \(\epsilon=0.1\) and observe that all three algorithms solve the Monge–Ampère equation successfully.

However, the Bellman algorithm requires the additional interpolation step (step [step:interpolation]) in order to converge, see Figure 2 (a). In fact, the solution to the corresponding Poisson equation, \(\Delta u = 2\sqrt{f}\) fails to be convex when \(0 \le \epsilon \lesssim 0.2\), thus the first approximant produced by the Bellman method (as well as M2) is not convex. After a few iterations in the Bellman method, convexity of \(u_k\) is achieved and the interpolation step is no longer needed, see Figure 3.

For this example, the difference in running time is more pronounced, and the Bellman method runs 10–13 times faster than the M2 method. The Bellman method converges in less than 10 iterations, whereas the M2 method require around 140 iterations. See Figure 2 for comparison. Again, the M1 method runs much slower and requires a large, increasing number of iterations as \(N\) grows. For a \(255 \times 255\) grid, the Bellman method took 31 s, the M2 needed 376 s, while M1 required around 3.5h.

5.2 Mildly singular examples↩︎

5.2.1 Trigonometric↩︎

On the domain \(\Omega=[0,1]^{2}\) we consider the example where the exact solution is given by \[u(x,y)=-\cos{(0.5\pi x)}-\cos{(0.5\pi y)}\] with the corresponding right hand side of the Monge–Ampère equation \[\mathop{\mathrm{MA}}(u) = f(x,y)=(0.5\pi)^{4}\cos{(0.5\pi x)}\cos{(0.5\pi y)}.\] We remark that \(f=0\) on the boundary where \(x=1\) or \(y=1\). This implies that for \(x\in \Omega\) near this part of the boundary, the Bellman function gives less and less strictly positive matrices \(B_{k}(x)\). Despite this complication, all three algorithms converge and produce the same result, see Figure 4 (a). For this example, the Bellman algorithm runs 3 times faster than the M2 method, see Figure 4 (b). The number of iterations needed for the convergence is 6–10 for the Bellman algorithm and 35 for the M2 algorithm, regardless of the mesh size. Yet again, M1 runs much slower: on a \(255 \times 255\) grid, the Bellman method needs 31 s, the M2 method just over 90 s and M1 needs a little over 3 h.

Error in supremum (top) and \(\ell^{2}\) norm (bottom). The results for all three algorithms are almost identical and suggest that \(\|u-u_N\| \sim N^{-2}\) for both norms. All three graphs overlap in the picture

Running time (s). The results suggest that the running time is \(\sim N^{2.7}\) for the Bellman method and \(\sim N^{2.8}\) for the M2, while the running time for M1 is \(\sim N^{4.0}\)

5.2.2 Degenerate example↩︎

Let us reconsider Example 5.1.2, i.e. Equation 9 , and remove the regularisation term by setting \(\epsilon=0\). In other words, we consider the exact solution \[u(x,y)=0.5(x-0.5)^{4} + y^{2},\] where \[\mathop{\mathrm{MA}}(u) = f(x,y) = 12(x-0.5)^{2}.\] Note that the right-hand side vanishes along the line \(x=0.5\), and while we can’t expect the Bellman method to produce approximants with positive definite Hessians along this line, nevertheless, all three algorithms successfully tackle the degeneracy and converge to the correct solution, again producing the same results.

As in the strictly convex examples, the Bellman method requires around 10 iterations to converge, regardless of mesh size, but for this example the M2 method struggles, and the number of iterations grows, at a roughly linear rate, as the mesh size \(N\) increases, see Table 1. This phenomenon heavily penalises the running time, which in the strictly convex settings was observed to be roughly proportional to \(N^\alpha\) with \(\alpha \approx 2.7\) for both the Bellman and M2 methods.

For this example, the running time of the Bellman method is still around \(N^{2.7}\), but for the M2 method, the running time increases to \(\sim N^{3.6}\) giving the Bellman method a huge advantage for large mesh sizes. For a \(511 \times 511\) mesh, the Bellman method ran in just over four minutes, and the method M2 needed almost eleven hours (on a reasonably fast personal computer). Here, the performance of the M1 method was closer to the M2 method, but still about three times slower.

Error in supremum (top) and \(\ell^{2}\) norm (bottom). The results for all three algorithms are once again almost identical and suggest that \(\|u-u_N\| \sim N^{-2}\) for both norms

Running time (s). The results suggest that the running time is \(\sim N^{2.7}\) for the Bellman method, \(\sim N^{3.5}\) for M2 method and \(\sim N^{4.0}\) for M1

5.2.3 Constant Monge–Ampère↩︎

Let us consider the Dirichlet problem \[\label{eq:constant95MA} \begin{cases} \mathop{\mathrm{MA}}(u) \equiv 1 & \text{on \Omega = [-1,1]^2} \\ u \equiv 1 & \text{on \partial\Omega} \end{cases}\tag{10}\] The exact solution of 10 is unknown, and even though the setup looks innocent enough at first glance, this example has shown to be a challenge to most if not all numerical algorithms that have been proposed for the Monge–Ampère equation. The problem is that the boundary conditions force \(u\) to be almost constant on lines parallel to and close to the boundary of the square, and the condition \(\det(D^2 u) = 1\) then forces the second order derivative in the normal direction to blow up. It follows from results by Błocki [10] that the solution is \(C^{1,1}(\Omega) \cap C(\bar\Omega)\), but by the argument above, we can’t expect higher order regularity of \(u\) up to \(\partial\Omega\).

Running the Bellman algorithm on Equation 10 produces a function that (barely) fails to be convex near the corners of \(\Omega\), but the same is true for Method M2, see Figure 6. This example was also studied by Dean and Glowinski [27] (using a variational formulation) and Oberman [18] (wide stencil finite difference discretization), all with similar results. In [14], they suggest comparing the different methods using the minimum (i.e.the value at the origin) of the solution. The wide stencils methods from Oberman [18] are known to converge monotonically to a viscocity supersolution of Equation 1 , and are thus too large. If the solution produced by the Bellman method had been convex, Proposition 2 would have implied that it too is too large, but since convexity fails near the corners of the square, we have no guarantee for this. Our method produces a solution which is slightly larger than the M2 and M1 methods, see Table 2.

The two methods used in [14] were observed to have running times of \(\sim N^{4.0}\) (M1) and \(\sim N^{4.6}\) (M2), respectively (results were reported for grid sizes \(21 \le N \le 141\)). When we ran those two methods on grid sizes up to \(255 \times 255\) (for M1), we observed similar running times.

The M2 method converges very slowly even for moderately large grid sizes, and was terminated after 10,000 iterations when \(N \ge 127\). The available results for smaller meshes suggest a running time of \(\sim N^{4.6}\). For this example, M1 outperforms M2 and runs at \(\sim N^{4.0}\). The Bellman method is an order of magnitude faster still, and clocks in at \(\sim N^{3.0}\), see Table 3. Already on a \(101 \times 101\) grid, the Bellman method runs more than 300 times faster than the M2 method and about 80 times faster than M1.

5.3 Singular examples↩︎

5.3.1 Circular degeneracy↩︎

Error in supremum and \(\ell^{2}\) norm. The errors for Bellman method are worse with \(\|u-u_N\|_{\infty} \sim N^{-0.8}\) and \(\|u-u_N\|_{\ell^2} \sim N^{-1.0}\). The other methods produce almost identical results with \(\|u-u_N\|_{\infty} \sim N^{-1.3}\) and \(\|u-u_N\|_{\ell^2} \sim N^{-1.5}\)

Running time (s). The results suggest that the running time is \(\sim N^{2.9}\) for the Bellman method, \(\sim N^{3.9}\) for the M2 method and \(~\sim N^{4.0}\) for M1

Difference \(u_{B}-u\) (light) and \(u_{M2}-u\) (dark), given by Bellman and Method M2, respectively

Comparing the running time to the achieved accuracy (\(\ell^2\) norm), the performance of the three methods are fairly similar

We consider \[\label{eq:32circular95degeneracy} u(x,y)= 0.5(\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}-0.2)^{+},\tag{11}\] which would be a solution to the Monge–Ampère equation with the following right-hand side \[\mathop{\mathrm{MA}}(u) = f(x,y)=\frac{(\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}-0.2)^{+}}{\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}}.\] In this example, \(u \in C^{1,1}\), but \(f = \mathop{\mathrm{MA}}(u)\) vanishes on a disc centred at \((0.5,0.5)\) with radius \(r=0.2\). Here, both the M1 and the M2 methods outperform the Bellman algorithm, see Figure 7 (a). As expected, the Bellman algorithm struggles on large areas where the exact solution is not strictly convex, that is on the aforementioned disc. This is the area where the solutions \(u_{\text{B}}\) and \(u_{M2}\) deviate the most from the exact solution \(u\), see Figure 7 (c). In fact, we see that the error in the supremum and \(\ell^{2}\) norms are proportional to \(N^{-0.8}\) and \(N^{-1.0}\), respectively, for the Bellman algorithm, whereas they are around \(N^{-1.3}\) and \(N^{-1.5}\) for the other methods.

While the accuracy of the Bellman method is notably worse than the other methods, it seems like all methods still converge to the exact solution as the mesh size \(N\) tends to \(\infty\), but the rate of convergence, for all methods, is slower than for the strictly convex examples, see Figure 7 (a).

Even though the Bellman method gives worse results, it still outperforms the other methods with respect to running time. The measured time complexity for the Bellman method is similar to the strictly convex examples (\(\sim N^{2.9}\)), whereas the running time of the M2 method is \(\sim N^{3.9}\) , resulting in a running time that is about 70 times slower than the Bellman method for a \(511 \times 511\) mesh. The number of iterations required is 6–14 for the Bellman method, where Method M2 needs an increasing number of iterations as the number of mesh points increases (1092 for the \(511 \times 511\) mesh). The M1 method runs even slower. If we look at accuracy as a function of computational time (see Figure 7 (d)), we see that all three methods perform fairly similarly.

5.3.2 Unbounded Monge–Ampère↩︎

We investigate the convergence for an unbounded right-hand side. For \(\Omega=[0,1]^{2}\) and \[u(x,y)=-\sqrt{2-x^{2}-y^{2}}.\] The right-hand side of the equation is \[\mathop{\mathrm{MA}}(u) = f(x,y)=\frac{2}{(2-x^{2}-y^{2})^{2}}\] and we note that \(f\) is unbounded near the point \((1,1)\). All three algorithms struggle, and only achieve an accuracy of order \(N^{-0.5}\) (measured in sup-norm), which is the same as observed by Benamou et al. [14] for M1 and M2. Most of the error is concentrated near the problematic point \((1,1)\), and if convergence is measured by average \(\ell^2\) norm instead, the accuracy is \(N^{-1.5}\). All algorithms run slowly, and the Bellman method and M2 runs at roughly the same speed, while M1 is about 20 times slower.

Curiously enough, the Bellman method seems to enjoy a relative speed-up as well as an improvement in accuracy for larger meshes, see Figure 8.

If we instead do the computations on a slightly smaller square, e.g. \([0, 0.99]^2\), the three methods run at full accuracy (\(\sim N^{-2}\)), and the Bellman method runs 2–3 times faster than M2.

Error in supremum (top) and \(\ell^{2}\) norm (bottom). The results for all three algorithms are fairly close, but not identical, and suggest that \(\|u-u_N\|_{\infty} \sim N^{-0.5}\) and \(\|u-u_N\|_{\ell^2} \sim N^{-1.5}\)

Running time (s). The results suggest that the running time is \(\sim N^{3.0}\) for the Bellman method, \(\sim N^{3.1}\) for M2 and \(\sim N^{4.0}\) for M1

Comparing the running time to achieved accuracy (\(\ell^2\) norm), the performance of the three methods are fairly similar, with the Bellman method winning for very fine meshes

6 Conclusions↩︎

The Bellman algorithm relies heavily on the assumption that the solution to the Monge–Ampère equation is strictly convex on the domain. When this assumption is satisfied, the performance of the algorithm is very good. Using an \(N\times N\) mesh, the discretization error in supremum and \(\ell^{2}\) norm is proportional to \(N^{-2}\), and is achieved with few iterations (typically fewer than \(10\)) resulting in a running time which is proportional to \(N^{2.7}\)–\(N^{2.8}\).

Similar performance was also observed in cases where the right hand side is mildly degenerate, such as the set where \(f=0\) is at most one dimensional. In this setting, our method requires a few more iterations to converge (typically 10–15), but this only has a minor effect on the running time. In comparison, the M2 method converges a lot slower for this class of examples, typically needing more and more iterations as the number of grid points increases, resulting in running times that vary between \(N^{2.8}\) and \(N^{3.6}\).

The net effect is that our method runs 3–10 times faster than the M2 method on smooth strictly convex examples, but as much as 20–100 (or more) times faster for the mildly degenerate cases. The M1 method is even slower, and has a running time proportional to around \(N^{4.0}\) for most examples. With the exception of the constant Monge–Ampère example, M1 runs substantially slower than M2.

In the most challenging examples, when \(f\) vanishes on an open set or when \(f\) is unbounded, the Bellman method begins to struggle. While the running time and the number of iterations is similar to the less challenging examples, the accuracy becomes worse resulting in supremum norm error between \(N^{-0.8}\) and \(N^{-0.5}\). In these examples, the errors are mostly concentrated near the zero set (or the unbounded locus) of \(f\), and measuring the error in \(\ell^2\) norm gives slightly better results, with an error between \(N^{-1.0}\) and \(N^{-1.5}\).

On the other hand, other published methods struggle with these more challenging examples as well. Method M2 outperforms the Bellman method in the circular degeneracy example (Example 5.3.1), at the cost of very slow running times. In all the other examples that we have considered, the accuracy of the Bellman method, M2 method and M1 are very similar, indicating that all methods converge to the optimal discretization of the underlying PDE on the given mesh, and in almost all examples the Bellman method runs considerably faster than the remaining methods.

References↩︎