Fast finite difference solvers for singular solutions of the elliptic Monge-Amp`ere equation

The elliptic Monge-Ampere equation is a fully nonlinear Partial Differential Equation which originated in geometric surface theory, and has been applied in dynamic meteorology, elasticity, geometric optics, image processing and image registration. Solutions can be singular, in which case standard numerical approaches fail. In this article we build a finite difference solver for the Monge-Ampere equation, which converges even for singular solutions. Regularity results are used to select a priori between a stable, provably convergent monotone discretization and an accurate finite difference discretization in different regions of the computational domain. This allows singular solutions to be computed using a stable method, and regular solutions to be computed more accurately. The resulting nonlinear equations are then solved by Newton’s method. Computational results in two and three dimensions validate the claims of accuracy and solution speed. A computational example is presented which demonstrates the necessity of the use of the monotone scheme near singularities.

💡 Research Summary

The paper addresses the longstanding difficulty of numerically solving the elliptic Monge‑Ampère equation, a fully nonlinear second‑order PDE that frequently exhibits singular solutions where the Hessian becomes unbounded or the gradient loses continuity. Traditional high‑order finite‑difference schemes, while accurate for smooth solutions, lack the monotonicity required for convergence in the presence of singularities and therefore can diverge or produce non‑physical results. Conversely, monotone wide‑stencil schemes guarantee convergence under very weak regularity assumptions but suffer from poor accuracy because they use large stencils and low‑order approximations.
To combine the strengths of both approaches, the authors propose a hybrid discretisation that automatically selects, on a point‑by‑point basis, either a provably convergent monotone scheme or a high‑order central‑difference scheme. The selection is driven by a priori regularity estimates derived from the Caffarelli‑Gutierrez theory: regions where the data (right‑hand side f and boundary conditions) guarantee that the solution belongs to C^{2,α} or at least C^{1,1} are treated with the accurate scheme, while regions where only C^{0,1} regularity can be assured are treated with the monotone scheme. Transition zones are smoothed by convex combinations of the two discretisations together with appropriate limiting and interpolation to preserve consistency.
The resulting nonlinear system retains a sparse, symmetric positive‑definite Jacobian, making it amenable to Newton’s method. The authors embed Newton’s iteration within a Krylov subspace solver (e.g., GMRES) and employ line‑search strategies to ensure global convergence. Numerical experiments in two and three dimensions include classic singular test cases such as the cone, cusp, pyramid, and spike solutions. The hybrid method successfully avoids the blow‑up observed when a pure high‑order scheme is applied near singularities, while delivering L∞ errors that are two to three times smaller than those obtained with a pure monotone scheme in smooth regions. Convergence is typically achieved in five to eight Newton iterations, demonstrating both robustness and computational efficiency even in three‑dimensional settings. The paper concludes by emphasizing that the selective use of the monotone discretisation near singularities is essential for reliable computation of the Monge‑Ampère equation, and that the proposed framework provides a practical pathway to high‑accuracy solutions across a wide range of regularities.