WIDE STENCIL FINITE DIFFERENCE SCHEMES FOR THE ELLIPTIC MONGE-AMPÈRE EQUATION AND FUNCTIONS OF THE EIGENVALUES OF THE HESSIAN

Certain fully nonlinear elliptic Partial Differential Equations can be written as functions of the eigenvalues of the Hessian. These include: the Monge-Ampère equation, Pucci’s Maximal and Minimal equations, and the equation for the convex envelope. In this article we build convergent monotone finite difference schemes for the aforementioned equations. Numerical results are presented.


Introduction
Certain fully nonlinear elliptic Partial Differential Equations can be written as functions of the eigenvalues of the Hessian. These include: the Monge-Ampère equation, the Pucci maximal and minimal equations, and the equations for the convex and concave envelope of the boundary data. In this article we present convergent approximation schemes for such equations, in the form of wide stencil finite difference schemes. We prove that numerical solutions converge uniformly to the unique viscosity solution of the equation without any regularity requirement on the solution.
To keep the ideas clear, we work in two dimensions and limit ourselves to the examples listed above. Our methods generalize to other equations which are functions of the eigenvalues, to equations with lower order terms, and to higher dimensions.
1.1. Presentation of the Equations. We begin by writing down the general form of the equations, and then show how the specific equations mentioned above fit into the framework. Write D 2 u(x) for the Hessian of u : R 2 → R and write λ − [u](x) ≤ λ + [u](x) for the smaller and larger eigenvalues of the Hessian of u, respectively. Let Ω ⊂ R 2 be a bounded domain, and consider the Dirichlet problem for the fully nonlinear, elliptic Partial Differential Equation Since F is elliptic, G is a non-decreasing function of the eigenvalues of the Hessian. Next we show how different choices of G which lead to the different equations. Key words and phrases. finite difference schemes, partial differential equations, viscosity solutions, Monge-Ampère equation, Pucci maximal equation, convex functions, convex hull.
(i) Convex functions satisfy λ − ≥ 0, and concave functions satisfy λ + ≤ 0 in the viscosity sense [17]. Setting f = 0 in (PDE) and setting (CE) G CE (M ) = λ − we get the equation for the convex envelope of the boundary data. Likewise setting G(M ) = λ + we get the equation for the concave hull of the boundary data.
(ii) The Monge-Ampère equation [19,12] is obtained when we set G(M ) ≡ det(M ) = λ − λ + . It has been the subject of much recent attention, in part because of its appearance in optimal mass transportation [11,3]. In order for the equation to be elliptic, the additional requirement that solution be convex needs to be imposed. This condition can be imposed using λ − ≥ 0 in the viscosity sense [17]. The result is the single equation and G m = Γλ − + γλ + where 0 < γ ≤ Γ. Solutions of these equations are the best underestimator and overestimators of solutions of elliptic equations with ellipticity constants related to γ, Γ, see [2, sec. 2.2] and [12, p. 442]. As such, they can be used to give pointwise upper and lower bounds on solutions of elliptic equations with unknown coefficients.
For each of the equations, the challenge is to discretize the nonlinear second order differential operator; adding zeroth or first order terms to the equation can be accomplished without too much difficultly, we refer to [16] for details.
1.2. Related work. We know of little early numerical work on these equations. Oliker and Prussner presented a convergent scheme for the Monge-Ampère equation in two dimensions [18]. It uses a natural discretization of the geometric Alexandrov-Bakelman interpretation of solutions.
Recently, Dean and Glowinksi have investigated numerical methods for the Monge-Ampère equation [8,7,6,10] and Pucci's equations [9]. As the authors are careful to point out, in order to converge, their method requires more regularity than can be expected for solutions.

1.3.
A semi-discrete scheme. To illustrate ideas, we begin with a semi-discrete scheme. The theoretical properties we will prove hold for the semi-discrete scheme, so it can be used to sketch out the main ideas in a simpler setting.
Consider, for the sake of illustration, the equation Using the classical Rayleigh-Ritz characterization of the eigenvalues, we can write where v θ = (cos θ, sin θ) is the unit vector in the direction of the angle θ. We substitute the standard centered difference approximation to the second derivative in the direction of v θ , to arrive at the semi-discrete scheme 1.4. Higher Dimensions. The discretization for the eigenvalues of the Hessian presented in Section 1.3 using the classical Rayleigh-Ritz formula for the eigenvalues of a symmetric matrix. This formula simplifies in two dimensions, since the maximum or a minimum over all directions to gives the largest and smallest eigenvalues, respectively. In order to compute the range of approximate eigenvalues in dimensions higher than three, the full formula must be used. This means finding the largest eigenvalue, using the max formula, then finding the next largest by taking a maximum over (approximately) orthogonal directions, and so on.
1.5. Wide stencil finite difference schemes. Monotonicity, along with consistency and stability is sufficient for convergence of finite difference schemes to the unique viscosity solution of a nonlinear, possibly degenerate elliptic PDE [1]. For second order equations, even if they are linear, it is may not be easy to build monotone approximation schemes without making further assumptions on the structure of the equation. In fact, Motzkin and Wasow showed that for a given fixed stencil, it is possible to find a linear elliptic equation for which no second order accurate, monotone scheme can be built [13]. Their example was non-constructive, but later work by Crandall and Lions [5] makes the problem clear. Suppose the equation is d 2 u/dv 2 , where v is a direction vector which doesn't line up with the grid. Then it is impossible to find a monotone second order accurate scheme. Since nonlinear PDEs often reduce to two derivatives in a given direction, the same difficulties apply in this case. Crandall and Lions resorted to using grid with impractically large stencils. An alternative is to use a wider, but still reasonable sized grid, and resulting in lower accuracy. The idea is to replace derivatives of the form d 2 u/dv 2 , for general directions vectors v, with d 2 u/dw 2 , where w is the closest direction vector to v which lines up on the grid. This incurs an error of which depends on the directional resolution dθ, which is a measure of the distance from v to w. If this can be done in a monotone way, then a convergent scheme can be built.
Wide stencil schemes which have previously been used by the author to approximate the equation for Motion of Level Sets by Mean Curvature [14] and the Infinity Laplace equation [15]. These schemes required using specific properties of the equations for which a discrete analogue could be found. The result (a median scheme for Mean Curvature, and a discretely locally minimizing Lipschitz extension for Infinity Laplace) led to an expression which could be shown to be consistent and monotone. For equations which are functions of the eigenvalues, we have a natural discrete characterization of the eigenvalues, which allows schemes for a wider class of equations to be built.
The schemes converge as the spatial resolution, h, and the directional resolution, dθ, go to zero. In practice, we work with reasonably small stencils (up to four times the grid spacing). The computational time for the wide stencil schemes scales no worse than linearly with the stencil size, depending on the particular choice of stencil. Wider stencil schemes have the advantage of better capturing the geometry of the solution, and enjoy approximate versions of the properties (e.g. rotational invariance) of the exact solutions, as well as being necessary to obtain theoretical convergence. The main disadvantage is the lower accuracy, which also depends on the directional resolution, for example, O(h 2 + dθ 2 ), as in (5.1).
It is natural to try to extend these schemes to obtain higher order accuracy. One approach to building higher order schemes would be to pursue an analogy with WENO schemes for first order equations. Using a foundation of a monotone scheme, it might be possible to selectively use a higher order scheme when the solution appears regular enough, and fall back on the lower order monotone scheme when necessary.
1.6. Properties of the semi-discrete scheme. While the scheme above is only semi-discrete, we set aside that consideration aside for now and study its properties. Many of these properties will also hold for the fully discrete schemes which follow.
This scheme is second order accurate, by which we mean, , for every twice differentiable function φ(x). Begin by substituting (1.2) into (1.1) and solving for the reference value u(x), thus rewriting the scheme in the explicit form, The scheme is defined locally, but by extending this definition to all points in the domain (with possible adjustments near the boundary) we can arrive at a scheme for the equation in a domain. While we make local arguments, the properties below apply globally. More precisely, the properties below apply to the solution operator of the scheme, when it is well-defined. The scheme is monotone: increasing the value at any of the neighboring values u(x + hv θ ) cannot decrease the value at u(x).
The scheme is degenerate: changing the value at some of the neighbors might have no effect on the value at u(x).
The scheme is stable, in fact it is non-expansive in the uniform norm. Consider a perturbation δ(x), writeũ(x + hv θ ) = u(x + hv θ ) + δ(x + hv θ ), and call the resulting valueũ(x). It follows that The scheme is nonlinear, but the nonlinearity is of a particular kind: it consists of a minimum of linear terms.
The consistency error doesn't give a global solution error. As we observed above, the consistency error for the scheme is second order. Since the equations support solutions which are not twice differentiable, (for example, we can take u(x, y) = |x| as solution for the convex envelope, or the Monge-Ampère equation). the interpolation error on an exact solution will not even be second order accurate. Thus the lack of a higher order convergence rate is not a deficiency of the scheme. Rather, the scheme is robust enough to handle non-smooth solutions.

Convergence Theory
The advantage of using viscosity solutions [4] is that the theory applies to whole classes of equations; the machinery for proving existence, uniqueness and stability results is the same for each equation. We review the corresponding theory for convergence of numerical schemes in the next sections.
2.1. Convergence of approximation schemes. One advantage of the viscosity solutions formulation is that we need only check consistency on twice continuously differentiable functions, even though we allow for non-differentiable solutions.
The fundamental result for the convergence of numerical schemes for fully nonlinear, degenerate equations is due to Barles and Souganidis.
Theorem 1 (Convergence of Approximation Schemes [1]). Consider a degenerate elliptic equation, (PDE), for which there exist unique viscosity solutions. A consistent, stable approximation scheme converges uniformly on compact subsets to the viscosity solution, provided it is monotone.
Note the theorem only gives a uniform convergence rate. As mentioned above, the lack of regularity of solutions in the general setting prevents asserting convergence rates like O(h 2 ), despite second order consistency for schemes.
Monotonicity is an essential property which we discuss in depth below. The definition of consistency is in terms of twice differentiable test functions: it amounts to checking Taylor series computations, see Definition 2. The theorem requires only a mild form of stability: our schemes always satisfy a stronger notion of stability: they are nonexpansive in the uniform norm.

2.2.
The comparison principle and monotonicity. The comparison principle is a nonlinear generalization of the maximum principle. When applied to (PDE) it can be written in Ω and u 1 ≤ u 2 on ∂Ω then u 1 ≤ u 2 in Ω Some authors prefer to avoid the minus sign in front of F [u] in (PDE), in which case the first inequality would be switched. A consequence of the comparison principle is the monotonicity of the solution operator. If u, v are solutions of the same equations, (PDE), but with different boundary conditions then We prefer to write this property more abstractly as follows. For a given equation F , as in (PDE), let S F : C(∂Ω) → C(Ω) be the solution operator, which takes boundary data g, to the solution u. Thus The the monotonicity property (2.2) becomes In this form, monotonicity of the solution operator takes the same form for the solution operator of the numerical scheme. Embracing a duplication of notation, we write S F for the solution operator of the corresponding numerical scheme, and take (2.3) as the definition of monotonicity.

Degenerate elliptic schemes.
Taking into account the Barles-Souganidis result, the problem of building convergent schemes for fully nonlinear equations is reduced to the problem of building consistent, stable, monotone schemes. However, building monotone schemes remains a challenge for many interesting equations. While monotonicity is a global property of the solution operator, we can take advantage of local conditions for the scheme which ensure monotonicity. Degenerate elliptic schemes are defined by a local structure condition which is the discrete analogue of the degenerate ellipticity condition for PDEs. The definition follows.
A finite difference scheme consists of a grid, G, and a finite difference equation, F which is defined on grid functions u. The grid is a graph consisting of a set of points, x i ∈ Ω, i = 1, . . . N , each endowed with a list of neighbors, N (i). The direction vectors x(j) − x(i) associated with the neighbors j ∈ N (i), defines the stencil at x i . A grid function is simply a real-valued function defined on the grid, with values u i = u(x i ).
Write the finite difference equation in the form where u j is shorthand for the list of neighbors u j | j=N (i) . Using this form emphasizes that each term u i − u j (which usually occurs with the factor 1/|x i − x j |) represents the first order approximation to the derivatives in the direction of the neighbor, Second derivatives in the direction v ij are approximated by averaging du/dv ij and A boundary point is a grid point with no neighbors. Dirichlet boundary conditions u = g may be imposed at boundary points by setting F i [u] = u i − g(x i ). A solution is a grid function which satisfies F [u] = 0.
We write u = S F (g) when we wish to specify the equation, for the solution operator, provided it is well-defined. We regard the solution operator as a mapping from the Dirichlet data on the boundary points to grid functions.
We now define degenerate elliptic schemes.
Remark. We emphasize that each component of the scheme F i is a nondecreasing function of u i and the differences u i − u j for j = 1, . . . , N (i).

Theorem 2 ([16]
). Under mild analytic conditions, degenerate elliptic schemes are monotone, and stable (in fact they are non-expansive in the uniform norm). They can be solved by an iterative method with time step dt = K −1 , where K is the Lipschitz constant of the scheme, regarded as a function from R N → R N .
In addition, we can combine degenerate elliptic schemes for simple equations into schemes for more complicated ones. In particular, once we have built degenerate elliptic schemes for λ − , λ + , we build schemes for functions of the eigenvalues by simple substitution, see §4.2 Level Neighbors of (i,j) in First Quadrant Table 1. Neighbors of reference point (i,j) in the first quadrant. The neighbors in the other quadrants are given by rotation. . In order to fully discretize the scheme, we need to take second derivatives only in directions θ which are available on the computational grid. The limited directional resolution, of a finite grid requires the introduction of a second parameter dθ, in addition to the spatial resolution, h. The spatial resolution is improved by using more grid points, the directional resolution is improved by increasing the size of the stencil.
The indices of neighboring grid points for successively wider stencil schemes are indicated in Table 1. The first two levels are displayed in Figure 1(a).
3.1. Boundary grid points. Near the boundary of the domain, the stencil width is limited by available grid points. Second derivatives in a given direction are still computed, but to lower accuracy. This is achieved in two steps. First the boundary value at the intermediate point is determined, either by using the given value, if it is available, or by quadratic interpolation, using the three nearest neighbors. Next, the intermediate value is used to compute a first order accurate second derivative. To give a concrete example, in a typical calculation, we might compute Near the boundary, if the i − 2 value is not available, we compute instead Here the u(i − 1, j − 1/2) value is computed using quadratic interpolation at the indices (i − 1, k), k = j − 1, j, j + 1 via We note that this gives different values for the intermediate value depending on the direction of the computed derivate, but this is acceptable.

3.2.
Local spatial and directional resolution. On a regular cartesian grid, let the stencil at the reference point x 0 consist of the neighbors x 1 , . . . , x N (illustrated by open and closed circles respectively in Figure 1). Write the direction vectors in polar coordinates We assume the stencil is symmetric, if v i is a direction vector, then so is −v i and define the local spatial resolution For the purpose of calculating the time step for iterative schemes, we also define

The Schemes
In this section we present the schemes. From the basic building block of the discrete second derivative, we build monotone schemes for λ − , λ + on a wide stencil uniform grid. These operators are combined to give schemes for general functions of the eigenvalues.
We show the scheme is degenerate elliptic, which allows us to establish consistency, stability, and monotonicity, and thereby prove convergence using Theorem 1.

4.1.
Definition of the schemes for λ − , λ + . Now we can discretize the eigenvalues by the following simple formula.
Proof. Each discrete second derivative in the direction v i is the average of the terms u(x 0 ± v i ) − u(x 0 ) which have the form u j − u i , so the are certainly non-decreasing in u j − u i . Taking a minimum (or maximum) of several non-decreasing functions yields a non-decreasing function.

4.2.
Definition of the schemes for functions of the eigenvalues. Given the schemes (4.1), (4.2) for λ − , and λ + , respectively, we build schemes for any nondecreasing function G(λ − , λ + ) by simple substitution . These schemes are clearly consistent, with local accuracy determined by substitution. In addition, they are degenerate elliptic, as we show below.

Convergence
Next we prove convergence of the schemes. By Theorem 1, we need only to verify that the schemes are consistent and monotone. The first property involves a Taylor series calculation. The effect of the parameter dθ which reflects the directional resolution of the grid means requires an additional argument beyond the standard Taylor series arguments used for consistency in h.

Monotonicity.
Lemma 2 (Monotonicity). The scheme G h,dθ (λ − , λ + ), defined by (4.1),(4.2), (4.3) is monotone, provided G : R 2 → R is a nondecreasing function in each variable. In particular, the schemes for G M A , G CE , G M , G m are monotone. In addition, the schemes for G − f (x) are monotone. More generally, if G 2 is a degenerate elliptic scheme, then so is G + G 2 .
The last part of the theorem allows the schemes to be used as building blocks for more complicated equations. For example, to obtain schemes for the PDE G − f (Du, x) = 0 we simply combine the scheme for G with a degenerate elliptic scheme for f .
Proof. By virtue of Lemma 1 the schemes (4.1)(4.2) are degenerate elliptic. As was already shown in [16], it follows from the properties of nondecreasing functions, that if If F 1 and F 2 are degenerate elliptic finite difference schemes, then so is F = G(F 1 , F 2 ), provided G : R 2 → R is a nondecreasing function. It follows from Theorem 2 that degenerate elliptic schemes are monotone.
Furthermore, combining the scheme for G with the trivial scheme for f i = f (x i ) does affect the ordering properties, so the result is still degenerate elliptic. The last remark follows from the fact that a sum of non-decreasing functions is still non-decreasing.

Consistency.
Definition 2 (Consistency). We say the scheme F h,dθ is consistent with the equation (PDE) at x 0 if for every twice continuously differentiable function φ(x) defined in a neighborhood of x 0 , The global scheme defined on Ω is consistent if the limit above holds uniformly for all x ∈ Ω. (The domain is assumed to be closed and bounded).
Note that in practice we test consistency only at grid points. We implicitly assume the existence of an interpolation function that gives the value of the scheme at points in the domain which are not grid points. Typically, we can use linear interpolation, over simplices which partition the domain. Other types of monotone interpolation are also permitted.
Next we define the local spatial and directional resolution, and compute the consistency error of the scheme. The consistency error for the scheme gives a formal accuracy estimate for the scheme, which as we remarked above, is not a rigorous convergence rate.
Lemma 3 (Consistency). Let x 0 be a reference point with neighbors x 1 , . . . , x N , and direction vectors v i = x i − x 0 , i = 1, . . . , N , arranged symmetrically so that (3.1) holds. Use the schemes for λ − , λ + given by (4.1), (4.2), respectively. If φ(x) is a twice continuously differentiable function defined in a neighborhood of the grid, then and likewise for λ + . Consistency also holds for the scheme (4.3).
Proof. The standard Taylor series computation gives Let M be a given symmetric 2×2 matrix. Without loss of generality we may choose coordinates so that M is diagonal. For the unit vector This last equation recovers the Rayleigh-Ritz characterization of the eigenvalues: minimizing over θ yields λ − , and maximizing yields λ + . Next minimize over a discrete set θ ∈ {θ 1 . . . θ N }, with resolution dθ, defined by (3.3). Calculate which follows from the fact that cos(π − θ) = −1 + O(θ 2 ). A similar argument works for λ + , with minimum replaced by maximum.
Combing the h error with the dθ error gives the desired result. Inserting (5.1) into the scheme (4.3) shows that G h,dθ (λ − , λ + ) → G(λ − , λ + ) as h, dθ → 0. For G lin , the accuracy is the same. For G M A expanding the product shows that the accuracy is multiplied by a factor of (λ + + λ − ).
Remark. The factor of λ + − λ − in the error term shows that the dθ error arises differently from the h error. The spatial consistency error arises as the difference between the function and a quadratic function. The directional consistency error arises as the difference between the function an one whose Hessian is a multiple of the identity.

Convergence.
Theorem 3 (Convergence). Suppose G is a non-decreasing function of the two variables, and that unique viscosity solutions exist for the equation G(λ − , λ + ) = f (x). Then the finite difference scheme given by (4.3) converges uniformly on compacts subsets of Ω to the unique viscosity solution of the equation.
Proof. By virtue of Theorem 1, we need only show that the scheme is consistent and monotone. Consistency follows from Lemma 3 and monotonicity follows from Lemma 2.

5.4.
Explicit optimal time steps. The schemes yield a fully nonlinear equation defined on grid functions. In order to solve F h,dθ [u] = 0 we perform the iteration By Theorem 2, this iteration will converge to a fixed point which is a solution of the equation, provided the nonlinear CFL condition is respected, where K is the Lipschitz constant of the scheme, regarding F map from grid functions to grid functions. Here dt is a parameter, so-named because it corresponds to the time step in the explicit Euler discretization of u t = F [u]. Next we compute the constant K for each of the schemes.
In the above, we can compute K(i) which is the Lipschitz constant of F i at the grid point x i . We can set K = max i K(i), and set dt accordingly. Alternatively, we can update at each grid point, using a local dt, and choosing the ordering of the grid points for the update.
We can do this by noting that the Lipschitz constant for λ h,dθ − , λ h,dθ + , defined by (4.1) (4.2), is given by 2h −2 . We write this as, where h is given by (3.4). This follows since the Lipschitz constant of a maximum or minimum is bounded by the maximum of the Lipschitz constants. Finally, we have Taking the worst possible value for |v i |, and using the definition of h (3.4) gives the result.
When the equation is linear in the eigenvalues of the Hessian, as is the case for the convex envelope equation, G CE , and the maximal and minimal Pucci equations, G M , G m , we calculate the time step as follows. Write . Therefore the optimal time step is  Table 2. Errors for the exact solution of Monge-Ampère , u(x) = exp(|x| 2 /2), f (x) = (1 + |x| 2 ) exp(|x| 2 ), using the 9 and 17 point schemes.
The resulting allowable time step in this case can be significantly less restrictive than the one given by the the dimensional scaling O(h 4 ).

Computations
In the section we perform several convergence tests for the Monge-Ampère equation and the Pucci equations, using examples as a basis of comparison from the works of Dean and Glowinski. In addition, we perform computations to test the accuracy dependence on the dθ parameter. Numerical solutions were computed for this example, see Table 2.
A second example involved solved the Monge-Ampère equation G M A [u] = 1 with constant boundary data 1 on the unit square. Since the function is constant on the boundary, which means one of u xx or u yy is zero on the boundary, classical solutions don't exist in this case. In particular, along the boundary, one of the eigenvalues is zero, so the larger eigenvalue must approach infinity as the boundary is approached. However the method used converges independently of the regularity of the solution.
Solutions were computed solutions using the 9 point and 17 point stencil. The difference in accuracy between the smaller and the wider stencil was not significant. See Table 3. Since the exact solution is unknown, we compared the minimum value of the solution, as well as the difference between solutions in the maximum norm. The solution is plotted in Figure 2. The figure was coarsened for the plot, since too many lines appear using the full grid. The contour plot shows that the level sets go from circles, which develop high curvature along the diagonals until they approach a square. The computed maximum and minimum eigenvalues, ranged from on the order of .1 and 10 to .001 and 100 as n increased.   Table 4. Errors for the exact solution of the Pucci equation, using the radial solution, as a function of of α and n.
6.2. Study in the directional resolution error. The consistency error of the scheme is given by (5.1). This section studies the numerical error of the scheme, using an example designed to measure the error in dθ. For this purpose, we take a quadratic solution of the form u α,θ (x, y) = (cos 2 θ + α sin 2 θ) x 2 2 + (1 − α) cos θ sin θxy + (cos 2 θα + sin 2 θ) y 2 2 with 1 < α, so that λ − [u](x, y) = 1 and λ + [u](x, y) = α. Then u α,θ is a solution of Monge-Ampère with f (x, y) = α. Solutions were computed for θ = i 20 π 4 , i = 1, . . . , 20, with α = 5 and n = 80. The computed solution was numerically exact whenever the direction θ lined up with the grid. However, for the 9 point stencil, the error was large, and reached a maximum of .13 with θ = π/8, which was midway between the grid directions. For the 17 and 33 point stencils, the error was acceptable, with a maximum of .004 for the 17 point scheme, and .002 for the 33 point scheme, again reaching a maximum when θ was between the grid directions. See Figure 3.
We can conclude that when we have a combination of a large ratio between λ − and λ + along directions which are not favorable to the grid, the dθ error can dominate. In more tame situations, the dθ error does not.
We solved on [−1, 1] 2 , with n = 256, 128, ..., for α = 2, 2.5, 3, using the 17 point scheme with interpolation on the boundary. Each time the spatial resolution was increased, there was a large error between the previous solution near the discontinuity at the boundary, which was better resolved. The solutions are increasing as a function of α, by the comparison principle. See Figure 4. By Theorem 2, we expect that numerical solutions of the difference equations solved on the same grid respect the comparison principle as well. We verified that this held numerically, for each of the equations, and also using coefficients which were functions of x. For example the variable coefficient Pucci equations (A(x)λ − + λ + ), where 1 ≤ A(x) ≤ 2. See Figure 5.
6.5. Boundary continuity of solutions. We remark that solutions need not be continuous up to the boundary. For example, If the boundary data is not convex, then there is no way to build a convex function which agrees at the boundary. Since both solutions of (MA) and (CE) are required to be convex, we can't expect boundary continuity. However this is a feature of the equation, not the scheme, so it shows the robustness of the scheme.
6.6. Accelerating Iterations. The iteration (5.2) is a simple, explicit, convergent method to find the solution of the difference equation. While it may thousands of iterations to converge, on the largest grids used, solution time was a few minutes.  Table 5. Number of iterations required for convergence comparing initial data given by a rough guess with the interpolated solution from a coarser grid.
A simple method to accelerate iterations was to interpolate the solution on a coarser grid. Given a grid size n × n, we solved the equation on a coarser grid of n/2 × n/2, and interpolated the solution onto the the finer grid. This proved to be faster than solving with a less accurate initial guess.
We compared the number of iterations to get a solution using the two methods. See Table 5, which records data from the solution of the convex envelope. The results suggest that the number of iterations using the coarse interpolant appears to O(n), compared to superlinear growth using less accurate initial data. Each iteration (5.2) is an explicit application of the finite difference scheme at each grid point.