Computer Graphics of Solutions of the Generalized Monge-Ampère Equation

Alfred Baldes and Ortwin Wohlrab

Abstract

Using the method of Minkowski or Alexandrov one finds simple discretizations of elliptic Monge-Ampère equations, including the equation of graphs with prescribed positive Gaussian curvature. It is shown how these discrete problems can be solved numerically, and computer graphics of the piecewise linear, convex solutions are presented.

1 Introduction

Provided a sufficiently regular surface S ⊂ ^{3} can be represented as a graph

S = S_{z} = {(x, y, z(x, y)) : (x, y) ∈ Ω}

of a function z = z(x, y) defined in some domain Ω ⊂ ^{2}, its most important intrinsic geometric invariant, the Gauss curvature, can be computed as

κ = \frac{z_{xx}z_{yy} − z_{xy}^{2}}{(1 + z_{x}^{2} + z_{y}^{2})^{2}} .

Thus, the inverse problem of trying to find a graph over Ω of prescribed Gauss curvature amounts to the solution of the Monge-Ampère equation

z_{xx}z_{yy} − z_{xy}^{2} = φ(x, y)(1 + z_{x}^{2} + z_{y}^{2})^{2} . (1.1)

Only slightly more general we here consider the equation

z_{xx}z_{yy} − z_{xy}^{2} = φ(x, y)R(z_{x}, z_{y}) , (1.2)

with R = R(p, q) defined in ^{2}, but immediately restrict ourselves to the elliptic case φ = φ(x, y) > 0 on Ω and R = R(p, q) > 0 in ^{2}. Without further restriction z = z(x, y) can therefore assumed to be convex. Because of its widespread occurrence in differential geometry the resulting Monge-Ampère equation serves as an important model in the study of the Dirichlet problem for a fully nonlinear partial differential equation of second order

(1.3) \{ \begin{matrix}z_{xx}z_{yy} − z_{xy}^{2} = φ(x, y) R(z_{x}, z_{y}) in Ω \\ z(x, y) = h(x, y) on ∂Ω . \end{matrix}

In this context it is natural to assume also the domain Ω to be convex, in order not to exclude even trivial boundary data h = h(x, y) right from the beginning.

2 Numerical Procedure

We want to approximate classical solutions of (1.3) by piecewise linear convex functions solving corresponding discrete problems first introduced by Minkowski [4] and later revived by Alexandrov [1].

Theorem 1 . Let Ω ⊂ ^{2} be a convex polygon with corners \underline{B} = [B_{s}, 1 ≤ s ≤ k] and distinguished interior points \underline{A} = [A_{j}, 1 ≤ j ≤ N]. Then for all boundary values \underline{h} = [h_{s}, 1 ≤ s ≤ k] and positive weights \underline{μ} = [μ^{j}, 1 ≤ j ≤ N] with

∑_{j=1}^{N} μ^{j} < A(R): = ∫_{2} \frac{dpdq}{R(p, q)} (2.1)

there exists a unique convex, piecewise linear function z = z(x, y) on Ω , such that

ω(R, z, A_{j}): = ∫_{χ_{z}(A_{j})} \frac{dpdq}{R(p, q)} = μ^{j} , 1 ≤ j ≤ N,

(2.2)

z(B_{s}) = h_{s} , 1 ≤ s ≤ k ,

and such that all interior vertices of the graph S_{z} project onto \underline{A}, all boundary vertices onto \underline{B}.

Here χ_{z} denotes the generalized gradient map with

χ_{z}(x_{∘}, y_{∘}) : = \{ \begin{matrix}(p, q) ∈ ^{2}: z(x, y) ≥ z(x_{∘}, y_{∘}) + p(x − x_{∘}) + q(y − y_{∘}) \\ for all (x, y) ∈ Ω \end{matrix} \}

Thus χ_{z}(x_{∘}, y_{∘}) is the set of all supporting directions in (x_{∘}, y_{∘}) ∈ Ω, and it is implicitly assumed that R = R(p, q) is at least locally integrable.

Theorem 2 . Let Ω be an arbitrary uniformly convex, bounded domain in 2, φ = φ(x, y) ∈ C^{α}(Ω), R = R(p, q) ∈ C^{α}(^{2}), both positive, α ∈ (0, 1). Assume

∫_{Ω} φ(x, y)dxdy < A(R) = ∫_{2} \frac{dpdq}{R(p, q)} (2.3)

and φ(x, y) ≤ κd^{β} , R(p, q) ≤ C(1+p^{2}+q^{2})^{γ} with positive constants κ, c, γ, β, 2γ ≤ 3 + β , d : = dist_{∂Ω}(x, y). Moreover, suppose we choose a sequence of convex polygons Ω_{m} , m ≥ 1 , Ω_{1} ⊂ Ω_{2} ⊂ …, with ∪_{m=1}^{∞} Ω_{m} = Ω, of boundary points \underline{B}^{m} = [B_{s}^{m}, 1 ≤ s ≤ k(m)] , B_{s}^{m} ∈ ∂Ω_{m} , of interior points \underline{A}^{m} = [A_{j}^{m}, 1 ≤ j ≤ N(m)] , A_{j}^{m} ∈ Ω_{m} , lim_{m→∞} k(m) = lim_{m→∞} N(m) = ∞, and of positive weights \underline{μ}_{m} = [μ_{m}^{j}, 1 ≤ j ≤ N(m)] , ∑_{j=1}^{N(m)} μ_{m}^{j} < A(R), such that for every disk D = D_{ρ}(x_{∘}, y_{∘}) with radius ρ > 0 and center (x_{∘}, y_{∘}) ∈ ^{2}

lim_{m→∞} ∑_{A_{m}^{j}∈D} μ_{m}^{j} = ∫_{D∩Ω} φ(x, y)dxdy.

Then, given any sequence of boundary values \underline{h}^{m} = [h_{s}^{m}, 1 ≤ s ≤ k(m)] with piecewise linear extensions on ∂Ω_{m} converging uniformly to h = h(x, y) on ∂Ω, the piecewise linear, convex solutions z_{m} = z_{m}(x, y) , (x, y) ∈ Ω_{m}, of the corresponding discrete problems (2.2) converge uniformly on Ω to a unique solution z ∈ C^{2,α}(Ω) ∩ C^{∘}(\overline{Ω}) of the classical problem (1.3).

For arbitrary Ω uniform convexity means that through every boundary point there exists a circle enclosing Ω with uniformly bounded radii. Each function z_{m} being defined only on Ω_{m}, uniform convergence on Ω of course means that for all ε > 0 these exists m_{∘} with

|z_{m}(x, y) − z(x, y)| < ε for m ≥ m_{∘} and (x, y) ∈ Ω_{m} .

Convergence of the boundary values may be defined via homeomorphisms between ∂Ω_{m}, m ≥ 1, and ∂Ω, provided by projection along rays from a fixed origin. For a more detailed discussion of the results above we refer to [2] and the references therein.

Given any bounded uniformly convex Ω ⊂ ^{2}, we might assume its center of mass being the origin, and think of ∂Ω given in polar coordinates. Then, Ω_{m}, m ≥ 1, can be taken to be the polygon determined by its boundary vertices \underline{B}^{m} = [B_{s}^{m}, 1 ≤ s ≤ k] on ∂Ω, where

B_{s}^{m} = ( r(θ_{s}) cos θ_{s} , r(θ_{s}) sin θ_{s}) , θ_{s} = \frac{2π}{k} (s − 1) ,

and we always take k = 6m. Then, for any boundary function h ∈ C^{∘}(∂Ω) we can fix \underline{h}^{m} = [h_{s}^{m}, 1 ≤ s ≤ k] by setting h_{s}^{m} = h(B_{s}^{m}). We chose the number of interior points in \underline{A}^{m} = [A_{j}^{m}, 1 ≤ j ≤ N] as N = 1+3m(m−1) = 1+∑_{i=1}^{m−1} 6i with A_{∘}^{m} the origin and the remaining 3m(m−1) points distributed on (m−1) equidistant layers as

A_{j}^{m} = ( \frac{i}{m}r(θ_{s}) cos θ_{s} , \frac{i}{m}r(θ_{s}) sin θ_{s}) , 1 ≤ i ≤ m − 1 ,

1 ≤ s(i) ≤ 6i , θ_{s} = π \frac{s − 1}{3i} , j = 1 + 3i(i − 1) + s .

Starting from a central hexagon it is easy to determine a regular triangulation T_{m} of Ω_{m} with interior vertices \underline{A}^{m} and boundary vertices \underline{B}^{m}. Given a function φ = φ(x, y) ∈ C^{α}(Ω), the weights \underline{μ}_{m} = [μ_{m}^{j}, 1 ≤ j ≤ N] are taken to be

μ_{m}^{j} = \frac{1}{3} ∑_{i=1}^{6} ∫_{Δ_{i}(A_{j})} φ(x, y)dxdy , where Δ_{i}(A_{j}) , 1 ≤ i ≤ 6 ,

are the triangles of our triangulation T_{m} having A_{j} as vertex. Then, if with some R = R(p, q) ∈ C^{α}(^{2}) inequality (2.3) and the growth conditions on φ and R in Theorem 2 are satisfied, the corresponding sequence of discrete solutions z_{m} = z_{m}(x, y) will indeed converge to the desired solution of (1.3). As explained in [3], taking the boundary conditions into account, the function z_{m} = z_{m}(x, y) is uniquely determined by its values in the vertices z_{m,j} = z_{m}(A_{j}^{m}) , 1 ≤ j ≤ N. Thus a problem as (2.2) with Ω, z, \underline{A}, \underline{μ}, \underline{B}, \underline{h} replaced by Ω_{m}, z_{m}, \underline{A}^{m}, \underline{μ}_{m}, \underline{B}^{m}, \underline{h}^{m} respectively constitutes a nonlinear system in N(m) variables, which we solved numerically by a Newton method with variable step length.

Difficulties arise from the fact, that, as necessary in each step, in order to compute the values ω(R, z_{m}, A_{j}) , 1 ≤ j ≤ N, and the derivatives with respect to z_{m,j} with all vertices (A_{j}^{m}, z_{m}(A_{j})) given, we have to find the unique possible connectivity satisfying the convexity condition. We have to determine the whole graph of z_{m} with its facets and edges as part of the convex hull of both all vertices in the interior and on the boundary. The most important use of a variable step length is to ensure that all points (A_{j}^{m}, z_{m}(A_{j})) are in fact true vertices of the corresponding graph, and all sets χ_{z_{m}}(A_{j}) in gradient space have nonvoid interior. For more details we have to refer to [2] again.

3 Computer Graphics

The first numerical computations of functions z_{m} = z_{m}(x, y) approximating solutions of (1.3) where done on an IBM mainframe at the University of California at Santa Cruz. There we also could get reasonable wireframe and raster computer graphics, but there was no fast and immediate connection between the computing engine and the graphics device. We continued our work at the SFB 256 in Bonn, where the final version of our program has been developed on the Silicon-Graphics IRIS-3D workstation. For debugging and performance tests on our algorithms the immediate availability of even simple wireframe graphics at any given moment during program execution turned out to be extremely helpful. Because often one is not so much interested in the solution of a single particular Dirichlet problem but rather would like to study the behavior of the solution as boundary values, parameters in the equation, or even the domain are changing, our program allows the onruntime specification of several parameters which pick the desired functions φ, R, h and the domain Ω out of one-parameter families encoded in corresponding subroutines. One more parameter controls the graphics display; you might want to see a model of the final result only, or alternatively get a glimpse of the graph S_{z} after each single Newton step.

As starting point of this numerical imitation of the continuation method mentioned above it is useful to have some explicitly known solutions. For Ω an ellipse,

Ω = \{(x, y) ∈ ^{2}: \frac{x^{2}}{a^{2}} + y^{2} < 1\},

and constant right hand side c > 0, we have

z(x, y) = b (\frac{x^{2}}{a^{2}} + y^{2}) − b , b = \frac{a\sqrt{c}}{2},

as solution of z_{xx}z_{yy} − z_{xy}^{2} = c with vanishing boundary data. A spherical cap yields a graph with constant Gauss curvature; thus for the unit disk D = D_{1}(0) and constant c , 0 < c ≤ 1, we get

z(x, y) = \sqrt{R^{2} − r^{2}} − \sqrt{R^{2} − 1} , R^{2} = \frac{1}{c} , r^{2} = x^{2} + y^{2};

as solution of

z_{xx}z_{yy} − z_{xy}^{2} = c (1 + z_{x}^{2} + z_{y}^{2})^{2}

with trivial boundary data again.

For example to study the one-parameter family of functions solving

z_{xx}z_{yy} − z_{xy}^{2} = 4 in the unit circle D,
z(x, y) = λ cos θ for 0 ≤ θ < 2π , x = cos θ, y = sin θ,

we start with the paraboloid z(x, y) = x^{2} + y^{2} − 1 solving the problem for λ = 0; then we change the parameter step by step, always using the solution found in the step before to construct reasonable initial values for our Newton solver. The shaded images Ia, Ib, Ic show the solutions for λ = 0.2, λ = 0.6, λ = 1.0 respectively. All graphs shown in this section were computed with the choice m = 20 , k = 120 , N = 1141 , and consist of 6m^{2} = 2400 triangular facets. II shows the solution of (1.3) for

φ(x, y) = 1 + x + 2y^{2} , R(p, q) = (1 + p^{2} + q^{2})^{\frac{3}{2}} ,

and constant boundary data. The family corresponding to φ(x, y) = 1 + μ(x + 2y^{2}) , 0 ≤ μ ≤ 1 had been considered.

Returning to the case of prescribed Gauss curvature R(p, q) = (1+p^{2}+q^{2})^{2} we must take β ≥ 1 in Theorem 2. Thus, φ = φ(x, y) has to vanish Lipschitz continuously on ∂Ω, for us to be sure that a solution of (1.3) exists for all boundary values. Nevertheless even if this last conditions is not met, all discrete solutions z_{m} = z_{m}(x, y) will still exist provided inequality (2.3) hence condition (2.1) is satisfied. The functions solving

z_{xx}z_{yy} − z_{xy}^{2} = \frac{3}{4}(1 + z_{x}^{2} + z_{y}^{2})^{2} in D

z(x, y) = λ cos 3θ for 0 ≤ θ < 2π , x = cos θ, y = sin θ

for λ = 0.2 and λ = 0.4 are shown by III. All graphics in this paper has been produced at the University of Massachusetts at Amherst using the VPL graphics programming environment, as described in [3]. Output device had been the Apple LaserWriter and for color display the RasterTechnologies Model I/380.

After fixing a coloring of the lower hemisphere, in RGB - terms for example

R = −z , G = \frac{1 + x}{2} , B = \frac{1 + y}{2} ,

x^{2} + y^{2} + z^{2} = 1 , z ≤ 0,

we could color each facet differently according to its unit outward normal; we used simple flat shading at this point, for Gouraud shading later resulted in less interesting images. Some specular parallel light and point light sources were added to produce a more realistic impression. IVa, IVb, IVc show solutions on the ellipse

Ω = \{(x, y) ∈ ^{2} : \frac{x^{2}}{a^{2}} + y^{2} < 1\}

with a = \frac{3}{2}, corresponding to

φ(x, y) = \frac{1}{2} + μ(2 − x^{2}) , R(p, q) = (1 + p^{2} + q^{2})^{2}

z(x, y) = h(θ) = λ cos 2θ on ∂Ω

for (λ, μ) = (0, 0) , (λ, μ) = (0, \frac{1}{8}) , (λ, μ) = (\frac{3}{10}, \frac{1}{8}) respectively. Here the deformation of a spherical cap into the graph IVa had to be computed first.

REFERENCES

l. A. D. Alexandrov, Dirichlet’s Problem for the Equation Det ||u_{ij}|| = Φ(z_{1} … z_{n}, z, x_{1} … x_{n}) , Vestnik Leningrad Univ. Math. Mekh. Astronom 13 (1958), 5-24.

2. A. Baldes and O. Wohlrab, Convex Discrete Solutions of Generalized Monge-Ampère Equations, Preprint SFB 256, Bonn (to appear).

3. M. J. Callahan, D. Hoffman and J. T. Hoffman, Computer Graphics Tools for the Study of Minimal Surfaces, Communications of the ACM 31 (1988), 648-661.

4. H. Minkowski, Volumen and Oberfläche, Math. Ann. 57 (1903), 447-495.

A.B., Mathematisches Institut, Beringstr. 4, 53 Bonn, West Germany O.W., Sonderforschungsbereich 256, Wegelerstr. 8, 53 Bonn, West Germany.

Research supported in part by DOE Grant DE-FG02-86ER25015.