Numerical Solutions of Boundary Value Problems for K-Surfaces in R^{3} Frank E. Baginski^{1} The George Washington University Nathaniel Whitaker^{2} The University of Massachusetts at Amherst Abstract

A K-surface is a surface whose Gauss curvature K is equal to a positive constant. In this paper, we will consider K-surfaces that are defined by a nonlinear boundary value problem. In this setting, existence follows from some recent results on nonlinear second-order elliptic partial differential equations. The analytical techniques used to establish these results motivate effective numerical methods for computing K-surfaces. In theory, the solvability of the boundary value problem reduces to the existence of a subsolution. In an analogous way, if an approximate numerical subsolution can be determined, then the corresponding K-surface can be computed. We will consider two boundary value problems. In the first problem, the K-surface is a graph over a plane. In the second problem, the K-surface is a radial graph over a sphere. From certain geometrical considerations, it follows that there is a maximum allowable Gauss curvature K_{max} for these problems. The principal results in this paper are numerical estimates of K_{max} for a variety of geometries and boundary data. Using a continuation method, we determine numerically the unique one-parameter family of K-surfaces that exist for K ∈ (0, K_{max}). We can compare our numerical estimates for K_{max} to the true value when the Ksurface is a subset of a hyperbolic spherical surface of revolution. In this case, we find that our numerical estimates for K_{max} are in close agreement with the expected values.

Typeset by AMS-TEX

1Research supported in part by NASA under grant NAG5-697. This work was initiated while the author was on sabbatical leave at the Center for Geometry Analysis Numerics & Graphics.

2This work was suppported in part by ....

1. Introduction

A K-surface is a surface whose Gauss curvature K is equal to a positive constant. In this paper, we will consider K-surfaces in R^{3} with boundary. The study of K-surfaces of this type can be reduced to the study of certain Dirichlet problems. The existence of a K-surface with a given boundary can be reduced to the solvability of a nonlinear second-order elliptic partial differential equation of the Monge-Amperé type. These problems have received considerable attention in recent years (see, e.g., [3], [4], [7], [8], [9], [11], [14], [15], [16]). In this paper, we are concerned with the computation of numerical solutions of two nonlinear boundary value problems involving K-surfaces. In the first problem, the boundary of the K-surface is a curve in R^{3} that projects to a rectangle in the plane. In the second problem, the boundary of the K-surface consists of two disjoint closed curves over the sphere.

In Section 2, we formulate the problem of determining a K-surface in R^{3} that can be expressed as a graph over a rectangle in the plane. The problem is equivalent to finding a solution u that solves the boundary value problem: Given a positive constant K and a smooth function Φ, find u satisfying,

(P1)

det(u_{αβ}) = K(1 + |∇u|^{2})^{2}, on Ω ⊂ R^{2},
u = Φ, on Γ = ∂Ω.

The existence of a solution of (P 1) is a consequence of a more general result by Caffarelli, Nirenberg, and Spruck [3] and Krylov [14]. For a smooth and strictly convex boundary, the existence of a solution of (P 1) can be reduced to the existence of a convex subsolution \underline{u} ∈ C^{∞}(Ω̄). Problem (P 1) is replaced by its discretization (see Eqn. (2.5)) which is solved numerically in a way that parallels the analytic techniques used to prove the solvability of (P 1). In particular, a continuity method is used to determine the solutions of (2.5) for a fixed K. A continuation method is used to determine the range of K for which a solution of (2.5) exists. Equation (2.5) is derived in Appendix A.

One question of interest is to determine the maximal interval for which a K-surface with fixed boundary conditions can exist. We illustrate this with a simple example. Let Ω be a disk of radius R in the plane and let Φ = 0 on Γ. Suppose we wish to find the largest interval (0, K_{max}) for which (P 1) has a solution. Since a solution of (P 1) is a graph, it follows from simple geometrical considerations that (P 1) has a solution for every K ∈ (0, 1/R^{2}). A solution is a spherical cap of radius r, where K = 1/r^{2} and r > R. Besides a few simple geometries (as the previous example), there are no theoretical estimates on K_{max} for arbitrary domains. However, we can compute numerical estimates for K_{max} for each boundary value problem considered here. Numerical results for solutions of (P 1) are presented in Section 3.

In Section 4, we formulate the problem of determining a K-surface that is expressible as a radial graph over the sphere S^{2}. This problem is equivalent to the following boundary value problem: Given a positive constant K and a smooth function Φ, find u satisfying,

(P2)

a^{−1} det(u_{αβ}) = K(1 + |∇u|^{2})^{2}, on Ω ⊂ S^{2},
u = Φ, on Γ = ∂Ω.

Guan and Spruck [9] proved a result on the existence of K-surfaces that are expressible as radial graphs. The main theorem in [9] states: if Ω does not contain a hemisphere and Γ = radial-graph_{∂Ω}(Φ) bounds a locally strictly convex radial graph M over Ω, then for any K < K(M), Γ bounds a K-surface that is a radial graph. The discretization of (P 2) is discussed in Section 4 and Appendix A.

In Section 5, we present numerical results for K-surfaces expressible as radial graphs. For the examples presented in Section 5, the boundary Γ consists of two disconnected components Γ = Γ_{1} ∪ Γ_{2}. We consider a number of different cases, including Case 5.3 where Γ_{1} is a circle and Γ_{2} is a “saw-tooth” curve that lies in the sphere.

It should be noted that a rectangular domain, does not satisfy the smoothness properties that are assumed in the analytical work of [3]. However, if one considered (P 1) on a domain Ω′ ⊂ Ω, where Ω′ is obtained by smoothing the corners of Ω, one would expect to find results similar to those presented in Section 3. Although boundaries for the examples in Case 5.1 - 5.2 are smooth, the sawtooth boundary in Case 5.3 is not smooth. A result for K > 1 presented in Case 5.3 should be applied to a K-surface that is the radial graph for some Ω′ ⊂ Ω, where Ω′ is obtained by smoothing the corners of the sawtooth boundary. See Section 6 for a further discussion of these surfaces.

For each boundary value problem considered, there corresponds a one-parameter family of K-surfaces,

(K, u(K; Φ)), for K ∈ (0, K_{max}).

To determine the interval (0, K_{max}), we first determine some solution (K_{1}, u(K_{1}; Φ)) where K_{1} > 0. Theoretically, once (K_{1}, u(K_{1}; Φ)) is known, K-surfaces also must exist for every K ∈ (0, K_{1}), because u(K_{1}; Φ) can be used as a subsolution for any K ∈ (0, K_{1}). These K-surfaces are found by using the solution (K_{1}, u(K_{1}; Φ)) as an initial approximation to the solution u(K_{0}; Φ) where K_{1} − K_{0} = ΔK > 0 is sufficiently small. Typically, ΔK was on the order of 10^{−2}. Using a continuation method, we are able to march from K = K_{1} to K = 0, and determine the curve (K, u(K; Φ)) for 0 < K ≤ K_{1}. While a K-surface by definition is one with positive Gauss curvature, in Case 5.1 we are able to compute solutions of (P 2) for K < 0 (see Case 5.1). This was the only case where we found surfaces with K < 0.

We found that we were able to continue solutions for increasing values of K, enabling us to estimate K_{max}. In particular, the known solution (K_{1}, u(K_{1}; Φ)) could be used to compute the solution (K_{2}, u(K_{2}; Φ)) where 0 < K_{1} < K_{2}. This worked for most values of K in Problems (P 1) and (P 2). However, near the maximal value K_{max}, we found it useful to introduce a small positive perturbation ϵ and use u(K_{1}; Φ) + ϵ as an initial guess for the solution (K_{2}, u(K_{2}; Φ)) for solutions of (P 1). By adding this slight perturbation, we were able to approach K_{max} more closely. For solutions of (P 2), the positive perturbation is added to ρ (see Eqn. (4.1)), resulting in a decrease in the corresponding u, since u = 1/ρ.

For a number of cases, we are able to compare our numerical solutions to a known solution, a subset of a sphere of radius 1. For Cases 3.1, 5.1, and 5.3, with K = 1.00, we find that the corresponding residuals are on the order of 10^{−10} in absolute value.

When the boundary of the radial graph consists of two parallel circles on a sphere, the K-surface is a hyperbolic spherical surface of revolution (see [13, p. 153]). Its defining equation is an ordinary differential equation that can be solved exactly (see Eqn. (B.4)) and the true value of K_{max} can be computed. In order to determine the reliability of the computed value of K_{max}, we compared our computed values with the true values for Case 5.1. In Case 5.1a, we find that the computed value of K_{max} is 2.63, while the true value is 2.608. In Case 5.1b, the computed value of K_{max} is 1.15, while the true value is 1.147. These results give us confidence that our computed values of K_{max} are accurate to 0.03 and that our procedure for the computation of (K, u(K; Φ)) with increasing values of K is valid.

2. K-Surfaces that are graphs

In the following section, partial derivatives of a function will be indicated by a subscript, i.e., u_{1} = u_{x}, u_{2} = u_{y}, u_{11} = u_{xx}, u_{12} = u_{xy}, u_{22} = u_{yy}. A general partial derivative of a function will be denoted by a Greek subscript. For a surface parametrized as a graph,

S = {(x, y, u(x, y)) | (x, y) ∈ Ω ⊂ R^{2}},

the Gauss curvature is given by

(2.1)  K = \frac{det(u_{αβ})}{(1 + |∇u|^{2})^{2}} ,

where det(u_{αβ}) = u_{11}u_{22} − u_{12}u_{21} and |∇u|^{2} = u_{1}^{2} + u_{2}^{2}.

We are led to the following nonlinear boundary value problem for a K-surface. Given a positive constant K and a smooth function Φ, find the function u(x, y) satisfying

(P1)  \begin{aligned}det(u_{αβ}) & = K(1 + |∇u|^{2})^{2}, (x, y) ∈ Ω, \\ u & = Φ, (x, y) ∈ Γ = ∂Ω. \end{aligned}

The existence of a solution of (P 1) is equivalent to the existence of a strictly convex, strict subsolution (see [3]). We say that \underline{u} is a strictly convex, strict subsolution of (P 1) if

(2.2)  \begin{aligned}det(\underline{u}_{αβ}) & ≥ K(1 + |∇\underline{u}|^{2})^{2} + δ_{0}, (x, y) ∈ Ω, \\ u & = Φ, (x, y) ∈ Γ = ∂Ω, \end{aligned}

for some δ_{0} > 0. If δ_{0} = 0 in (2.2), then ū is called a strictly convex subsolution.

If \underline{u} satisfies (2.2) and the boundary is smooth, it follows from the implicit function theorem that there exists a function u^{0} with the properties,

(2.3)

det(\underline{u}_{αβ}) ≥ det(u_{αβ}^{0}) + ϵ_{0}, (x, y) ∈ Ω,
u^{0} = Φ, (x, y) ∈ Γ = ∂Ω,

where 0 < ϵ_{0} < δ_{0} (see [11]).

An existence proof for (P 1) requires establishing an apriori bound and the application of the continuity method. In the continuity method, one considers the following problem: find u^{t} such that

(2.4)

det(u_{αβ}^{t}) = tK(1 + |∇u^{t}|^{2})^{2} + (1 − t) det(u_{αβ}^{0}), (x, y) ∈ Ω,
u^{t} = Φ, (x, y) ∈ Γ = ∂Ω,

where 0 ≤ t ≤ 1. By applying the implicit function theorem and classical Schauder theory, one can show that the set of t for which (2.4) is solvable is open. If one can establish an apriori estimate ||u^{t}||_{2+α} ≤ C, then the set of such t is also closed and therefore the whole interval. The solution of (P 1) is obtained by setting t = 1. The reader is referred to [7] for a more thorough discussion of apriori bounds and the continuity method.

The continuity method suggests a numerical approach to solving (P 1). The numerical method by which we approximate a solution of (P 1) parallels the steps we have just outlined. First, we discretize (P 1), replacing partial derivatives by finite difference approximations, Ω by Ω^{h}, and Φ by Φ^{h} (see Appendix A). We are lead to the following discrete analogue of (P 1),

(2.5)  \begin{aligned}det(u_{αβ}^{h}) & = K(1 + |∇u^{h}|^{2})^{2}, (x_{i}, y_{j}) ∈ Ω^{h}, \\ u^{h} & = Φ^{h}, (x_{i}, y_{j}) ∈ ∂Ω^{h}, \end{aligned}

where i = 2, …, n_{x} − 1 and j = 2, …, n_{y} − 1 (see Eqn. (A.2)). A solution of Eqn. (2.5) is a set of u_{i,j}^{h}’s, where u_{i,j}^{h} denotes the value of u^{h} at the grid point (x_{i}, y_{j}). Ω^{h} is the set of interior grid points,

Ω^{h} = {(x_{i}, y_{j}) | i = 2, …, n_{x} − 1, j = 2, …, n_{y} − 1},
∂Ω^{h} = {(x_{i}, y_{j}) | i ∈ {1, n_{x}} and 1 ≤ j ≤ n_{y} or 1 ≤ i ≤ n_{x} and j ∈ {1, n_{y}}}.

We follow the convention that u^{h} denotes a numerical solution whose accuracy depends on the grid size h (see Appendix A). See Appendix A for a further discussion of the discretization.

The second step in determining a numerical solution is to find \underline{u}^{h}, an approximate subsolution satisfying

\begin{aligned}det(\underline{u}_{αβ}^{h}) & ≥ K(1 + |∇\underline{u}^{h}|^{2})^{2} + δ_{0}, (x_{i}, y_{j}) ∈ Ω^{h}, \\ \underline{u}^{h} & = Φ^{h}, (x_{i}, y_{j}) ∈ ∂Ω^{h}, \end{aligned}

i = 2, …, n_{x} − 1 and j = 2, …, n_{y} − 1. Then, we determine u^{0,h} such that

(2.6)

det(\underline{u}_{αβ}^{h}) ≥ det(u_{αβ}^{0,h}) + ϵ_{0}, (x_{i}, y_{j}) ∈ Ω^{h},
u^{0,h} = Φ^{h}, (x_{i}, y_{j}) ∈ ∂Ω^{h}.

In practice, u^{0,h} = \underline{u}^{h} for our computations. While in theory δ_{0} is strictly positive, numerically δ_{0} was allowed to be negative, but on the order of 10^{−10} in magnitude. When K was not near the maximum value, we found that we were able to compute an accurate K-surface even if we started with a nonconvex approximate subsolution. See Section 6 for a further discussion.

After u^{0,h} is determined, we solve the system

(2.7)

det(u_{αβ}^{h,t}) = tK(1 + |∇u^{h,t}|^{2})^{2} + (1 − t) det(u_{αβ}^{0,h}), (x_{i}, y_{j}) ∈ Ω^{h},
u^{h,t} = Φ^{h}, (x_{i}, y_{j}) ∈ ∂Ω^{h},

for a sequence t = t_{n}, where 0 = t_{0} < t_{1} < t_{2} < ⋯ < t_{n} = 1. The desired solution corresponds to t_{n} = 1. Eqn. (2.7) reduces to Eqn. (2.5) when t = 1. It was found that n = 5 was sufficient for calculating K-surfaces that were graphs.

We require that the numerical solution u^{h,t} satisfies Eqns. (2.7) at each of the internal grid points in Ω^{h}. Eqn. (2.7) is a system of polynomial equations of degree 4 in the unknowns u_{i,j}^{h,t}. We used n_{x} = 25 and n_{y} = 25 for the computations presented in Section 3. At each stage, we solve a system of 529 polynomial equations.

The actual method used to solve (2.7) was a Levenberg–Marquardt algorithm, a variation of Newton’s method (see IMSL subroutine DNEQNJ [12, p. 776]). Subroutine DNEQNJ requires the analytic Jacobian of (2.7). The solution process requires an initial guess for a solution of (2.7) and u^{0,h} was used for this purpose. It took approximately one minute of CPU time on a Silicon Graphics Indigo workstation with an R4000 processor to compute one solution of Eqn. (2.7) for some t = t_{k}, when (n_{x}, n_{y}) = (25, 25).

In practice, it was not very difficult to determine a subsolution. This was done by extending Φ to the interior of Ω. The boundary data Φ was chosen so that this could always be easily done. If \underline{u}^{0} is a subsolution for K = K_{0} > 0, then it is also a subsolution for all K ∈ (0, K_{0}). It follows that we could compute a solution of (2.5) for K = K_{0} (call it u_{K_{0}}^{h}) and then use u_{K_{0}}^{h} as the initial subsolution for K = K_{0} − ΔK. Continuing in this fashion, we generate a solution curve (K, u_{K}^{h}) for K ∈ (0, K_{0}).

We found that it was possible to compute solutions for increasing values of K, if the size of ΔK was not too large. That is, if u_{K_{0}} was a solution of (2.5) for K = K_{0}, then u_{K_{0}}^{h} could be used as the initial guess to solve (2.5) for K = K_{0} + ΔK. Continuing in this fashion, we are able to compute an estimate for K_{max}. When the maximal value of K was approached, we found that convergence was aided by adding a small positive constant to the solution computed from the previous step. Typically, ΔK was on the order of 10^{−2}. If the solution process converged for some K′ and diverged for K′ + ΔK, we defined K_{max} = K′.

The initial subsolution u^{h,0} was chosen to be the natural extention of the boundary data Φ into the interior of Ω. In order to evaluate the effectiveness of our computational scheme, we computed (K, u_{K}^{h}) for K = 1 and steadily decreased K. Reversing the process, we then computed a solution for K = 0.01 and steadily increased K. We found that the solutions determined by starting at K = 0.01 and increasing K were essentially identical to the corresponding solutions that were computed by starting at K = 1 and decreasing K. For the results that are presented in Section 3, we began the continuation at K = 0.01 and then increased K until we obtained an estimate for the maximum Gauss curvature.

While K-surfaces by definition have positive Gauss curvature, we found that we were able to compute a solution of Eqn. (2.5) for K = 0 for the boundary data considered in Section 3. This is not surprising, since it is not difficult to construct a ruled surface that matches the boundary data. We were unable numerically to continue (K, u_{K}) to negative values of K. In Section 5, we present an example of boundary data where it was possible to use our solution process to compute solutions with K < 0 (see Case 5.1).

3. Numerical Results for Graphs

In this section, we will consider (P 1) for a number of different boundary conditions and rectangular domains. The rectangular domain will be denoted by

(3.1)  R_{a,b} = {(x, y) | − a ≤ x ≤ a, −b ≤ y ≤ b}.

Figure 3.1a: Rectangular domain - R_{a,b}

For Cases 3.1a, 3.2a, 3.3a, a = 0.57, b = 0.57, n_{x} = 25, n_{y} = 25, and (h_{x}, h_{y}) = (.0471, 0.0471). For Cases 3.1b, 3.2b, 3.3b, a = 0.72, b = 0.36, n_{x} = 25, n_{y} = 25, (h_{x}, h_{y}) = (0.0596, 0.0298). Sup-norms of the numerical solutions will be computed as follows. For a solution u^{h} on Ω^{h},

||u^{h}||_{∞} = max_{(x_{i},y_{j})∈Ω^{h}} |u^{h}(x_{i}, y_{j})|.

Case 3.1a.

Let a = 0.57 and b = 0.57, Ω = R_{0.57,0.57} and
(3.2)  Φ = \sqrt{1 − x^{2} − y^{2}}, (x, y) ∈ Γ = ∂R_{0.57,0.57}.

In this case, we found that a K-surface exists for every K ∈ (0, K_{max}) where K_{max} = 2.10. To denote explicit dependence on K and Φ, we write u^{h}(x, y; K, Φ). However, to simplify notation, we will not always write out this dependence. In Figure 3.2, we present a graph of the curve (K, û(K)), where û(K) denotes the sup-norm of a numerical solution for a particular value of K and Φ, i.e.,

(3.3)  û(K) = ||u^{h}(K, Φ)||_{∞}, K ∈ (0, K_{max}).
The C^{2}-norm of u^{h}(x, y) on Ω^{h} is given by
(3.4)  ||u^{h}||_{2} = ||u^{h}||_{∞} + ||u_{1}^{h}||_{∞} + ||u_{2}^{h}||_{∞} + ||u_{11}^{h}||_{∞} + ||u_{12}^{h}||_{∞} + ||u_{22}^{h}||_{∞}.

The terms u_{1}^{h}, u_{2}^{h}, u_{11}^{h}, u_{12}^{h}, u_{22}^{h} are defined in Appendix A. In Figure 3.3, we present a graph of (K, η̂(K)) where

(3.5)  η̂(K) = ||u^{h}(K, Φ)||_{2}, K ∈ (0, K_{max}).

In Figures 3.4a - 3.4d, we present the solutions u^{h}(x, y; K, Φ) for K = 0.01, 0.70, 1.30, 2.10.

R_{0.57,0.57}, Φ = \sqrt{1 − x^{2} − y^{2}}

11

Case 3.1b.

Let a = 0.72 and b = 0.36, Ω = R_{0.72,0.36}, and
Φ = \sqrt{1 − x^{2} − y^{2}}, (x, y) ∈ Γ = ∂R_{0.72,0.36}.

Figures 3.5, 3.6, 3.7 show the results of Case 3.1b. We find that a K-surface exists for all K ∈ (0, K_{max}) where K_{max} = 2.61. In Figures 3.7a - 3.7d, we present a graph of the solution u(x, y; K, Φ) for K = 0.01, 1.00, 2.00, 2.61.

R_{0.72,0.36}, Φ = \sqrt{1 − x^{2} − y^{2}}

13

Case 3.2a.

Let a = 0.57 and b = 0.57, Ω = R_{0.57,0.57}, and
Φ = 1 − x^{2} − y^{2}, (x, y) ∈ Γ = ∂R_{0.57,0.57}.

In this case, we found that a K-surfaces exists for every K ∈ (0, K_{max}) where K_{max} = 2.24. Figure 3.8 contains a sketch of (K, û(K)) for K ∈ (0, K_{max}) and Figure 3.9 contains a sketch of (K, η̂(K)) for K ∈ (0, K_{max}). In Figures 3.10a - 3.10d, we present a graph of the solution u(x, y; K, Φ) for K = 0.01, 0.80, 1.60, 2.24.

R_{0.57,0.57}, Φ = 1 − x^{2} − y^{2}

15

Case 3.2b.

Let a = 0.72, b = 0.36, Ω = R_{0.72,0.36}, and
Φ = 1 − x^{2} − y^{2}, (x, y) ∈ Γ = ∂R_{0.72,0.36}.

Figure 3.11 contains a sketch of (K, û(K)) for K ∈ (0, K_{max}) and Figure 3.12 contains a sketch of (K, η̂(K)) for K ∈ (0, K_{max}). We find that K_{max} = 2.73. In Figures 3.13a - 3.13d, we present a graph of the solution u(x, y; K, Φ) for K = 0.01, 1.00, 1.90, 2.73.

R_{0.72,0.36}, Φ = 1 − x^{2} − y^{2}

17

Case 3.3a

Let a = 0.57, b = 0.57, Ω = R_{0.57,0.57}, and
Φ = 1 − (x − 0.075)^{2} − (y − 0.15)^{2}, (x, y) ∈ Γ = R_{0.57,0.57}.

In this case, we found that a K-surface exists for every K ∈ (0, K_{max}) where K_{max} = 1.85. Figure 3.14 contains a sketch of (K, û(K)) for K ∈ (0, K_{max}) and Figure 3.15 contains a sketch of (K, η̂(K)) for K ∈ (0, K_{max}). In Figures 3.16a - 3.16d, we present a graph of the solution u(x, y; K, Φ) for K = 0.01, 0.58, 1.20, 1.85.

R_{0.57,0.57} Φ = 1 − (x − 0.075)^{2} − (y − 0.15)^{2}

19

Figure 3.16a - 3.16d: K-surfaces, R_{0.57,0.57}, Φ = 1 − (x − 0.075)^{2} − (y − 0.15)^{2} Case 3.3b

Let a = 0.72, b = 0.36, Ω = R_{0.72,0.36}, and

Φ = 1 − (x − 0.075)^{2} − (y − 0.15)^{2}, (x, y) ∈ Γ = ∂R_{0.72,0.36}.

In this case, we found that a K-surface exists for every K ∈ (0, K_{max}) where K_{max} = 2.27. Figure 3.17 contains a sketch of (K, û(K)) for K ∈ (0, K_{max}) and Figure 3.18 contains a sketch of (K, η̂(K)) for K ∈ (0, K_{max}). In Figures 3.19a - 3.19d, we present a graph of the solution u(x, y; K, Φ) for K = 0.01, 0.75, 1.50, 2.27.

R_{0.72,0.36}, Φ = 1 − (x − 0.075)^{2} − (y − 0.15)^{2}

21

Figure 3.19a - 3.19d: K-surfaces - R_{0.72,0.36}, Φ = 1 − (x − 0.075)^{2} − (y − 0.15)^{2}

In Table 3.1, we summarize the results of our numerical computations for K-surfaces that are expressible as graphs. A solution is accepted if the relative error between two successive approximations is less than 5.0 × 10^{−8}. The residual of a solution of (2.5) is defined in Appendix A (see Eqn. (A.3)-(A.4)). The residuals for solutions corresponding to K_{max} are presented in Table 3.1.

Case  Domain  Boundary data Φ  (0.0, K_{max})  Residual R^{h}
3.1a  R_{0.57,0.57}  Φ = \sqrt{1 − x^{2} − y^{2}}  (0.00, 2.10)  0.6719 × 10^{−10}
3.1b  R_{0.72,0.36}  Φ = \sqrt{1 − x^{2} − y^{2}}  (0.00,2.61)  0.9656 × 10^{−11}
3.2a  R_{0.57,0.57}  Φ = 1 − x^{2} − y^{2}  (0.00, 2.24)  0.7883 × 10^{−10}
3.2b  R_{0.72,0.36}  Φ = 1 − x^{2} − y^{2}  (0.00, 2.73)  0.53433 × 10^{−8}
3.3a  R_{0.57,0.57}  Φ = 1 − (x − 0.075)^{2} − (y − 0.15)^{2}  (0.00, 1.85)  0.7883 × 10^{−10}
3.3b  R_{0.72,0.36}  Φ = 1 − (x − 0.075)^{2} − (y − 0.15)^{2}  (0.00, 2.27)  0.2311 × 10^{−9}

Table 3.1: Maximal Intervals for which a K-surface exists 4. Radial Graphs

A surface S is said to be a radial graph over a sphere S^{2} if there is a mapping X : Ω → S ⊂ R^{3} where

(4.1)  X(x) = ρ(x)x, x ∈ Ω ⊂ S^{2}.

Assuming ρ(x) > 0, we can let u(x) = 1/ρ(x) and it follows that the Gauss curvature of S is (see [7])

(4.2)  K = \frac{a^{−1} det(u_{αβ})}{(1 + |∇u|^{2})^{2}},

where

(4.3)  u_{αβ} = ∇_{αβ}u + a_{αβ}u,

det(u_{αβ}) = u_{11}u_{22} − u_{12}u_{21}, and ∇_{α}, ∇_{αβ} are covariant derivatives on S^{2}. The metric on S^{2} is denoted by a_{αβ} and the inverse of a_{αβ} is denoted by a^{αβ}; a = det(a_{αβ}), the surface gradient ∇u is calculated with respect to the metric a_{αβ}, and |∇u|^{2} = a^{αβ}∇_{α}u∇_{β}u (repeated indices are summed).

Parametrizing the unit sphere in terms of the usual spherical coordinates, we find that a point x ∈ S^{2} can be represented in the form

(4.4)  x = (cos ϕ sin θ, sin ϕ sin θ, cos θ), (θ, ϕ) ∈ Ω,

where

Ω = {(θ, ϕ) | 0 ≤ θ ≤ π, 0 ≤ ϕ ≤ 2π}.

The metric tensor on S^{2} relative to (4.4) is given by a_{11} = 1, a_{12} = a_{21} = 0, a_{22} = sin^{2} θ. Covariant derivatives are given by ∇_{1}u = u_{θ}, ∇_{2}u = u_{ϕ},

∇_{11}u = u_{θθ},
∇_{12}u = ∇_{21}u = u_{θϕ} − cot θu_{ϕ},
∇_{22}u = u_{ϕϕ} + sin θ cos θu_{θ},

and |∇u|^{2} = |u_{θ}|^{2} + (sin θ)^{−2}|u_{ϕ}|^{2}.

For K-surfaces given by radial graphs, we are led to the following boundary value problem. Given a positive constant K and a smooth function Φ, find the function u(θ, ϕ) satisfying

(P2)  \begin{matrix}a^{−1} det(u_{αβ}) = & K(1 + |∇u|^{2})^{2}, (θ, ϕ) ∈ Ω, \\ u = & Φ, (θ, ϕ) ∈ Γ = ∂Ω. \end{matrix}

The discrete version of (P 2) has the same form as Eqn. (2.5). In particular, we have

(4.6)  \begin{matrix}det(u_{αβ}^{h}) = & K(1 + |∇u^{h}|^{2})^{2}, (θ_{i}, ϕ_{j}) ∈ Ω^{h}, \\ u^{h} = & Φ^{h}, (θ_{i}, ϕ_{j}) ∈ Γ^{h} = ∂Ω^{h}. \end{matrix}

See Appendix A for the definitions of Ω^{h}, u^{h}, Φ^{h}, (θ_{i}, ϕ_{j}) and other terms introduced in (4.6). As in Section 3, a solution of (4.6) is a set of u_{i,j}^{h}’s, where u_{i,j}^{h} represents the value of the approximate solution at the interior grid point (θ_{i}, ϕ_{j}).

The steps in the numerical solution process of (P 2) are similar to those used to solve (P 1). We first determine an approximate subsolution and denote it by \underline{u}^{h}. Following the steps that are outlined in Section 2, we determine a function u^{0,h} that satisfies

(4.7)

det(\underline{u}_{αβ}^{h}) ≥ det(u_{αβ}^{0,h}) + ϵ_{0}, (θ_{i}, ϕ_{j}) ∈ Ω^{h},
u^{0,h} = Φ^{h}, (θ_{i}, ϕ_{j}) ∈ Γ^{h} = ∂Ω^{h}.

Finally, after introducing the continuation parameter t, we are led to

(4.8)

det(u_{αβ}^{h,t}) = tK(1 + |∇u^{h,t}|^{2})^{2} + (1 − t) det(u_{αβ}^{0,h}), (θ_{i}, ϕ_{j}) ∈ Ω^{h},
u^{h,t} = Φ^{h}, (θ_{i}, ϕ_{j}) = Γ^{h} ∈ ∂Ω^{h}.

Eqn. (4.8) is solved for a sequence t = t_{n}, where 0 = t_{0} < t_{1} < t_{2} < ⋯ < t_{n} = 1, where t_{n} = 1 corresponds to the desired solution. Eqn. (4.8) reduces to (4.7) when t = 1. A summary of our numerical results, which will be discussed below, is given in Table 6.1.

In the following, we will consider K-surfaces over domains in the sphere. In Cases 5.1 - 5.2, we will consider domains in the form,

(4.9)  Ω_{ϑ_{1},ϑ_{2}} = {(θ, ϕ) | ϑ_{1} ≤ θ ≤ ϑ_{2}, 0 ≤ ϕ ≤ 2π}.

See Figure 4.1a for a sketch of Ω_{ϑ_{1},ϑ_{2}} in (θ, ϕ)-space. In order to define the sawtooth boundary, first we define the following sets,

(4.10)  B_{k} = \{(θ, ϕ) | ϑ_{2} − δ ≤ θ ≤ ϑ_{2}, \frac{1}{2}πk − \frac{1}{8}π ≤ ϕ ≤ \frac{1}{2}πk + \frac{1}{8}π\}. k = 0, 1, 2, 3.

In Case 5.3, we consider radial graphs over sawtooth domains. The sawtooth domain Ω_{ϑ_{1},ϑ_{2},δ} is formed by removing the B_{k}’s from Ω_{ϑ_{1},ϑ_{2}}, i.e.,

(4.11)  Ω_{ϑ_{1},ϑ_{2},δ} = Ω_{ϑ_{1},ϑ_{2}} \ ⋃_{k=0}^{3} B_{k}.

For the sawtooth boundary, ϑ_{2} = \frac{1}{2}π and δ = \frac{1}{4}(ϑ_{2} − ϑ_{1}). A sketch of Ω_{ϑ_{1},ϑ_{2},δ} is presented in Figure 4.2.

Figure 4.1a: Ω_{ϑ_{1},ϑ_{2}} Figure 4.1b Ω_{ϑ_{1},ϑ_{2},δ} Sawtooth domain

Remark. We found that for radial graphs, it was possible solve Eqn. (4.6) directly without using the continuation parameter t. This is equivalent to solving Eqn. (4.8) with t = 1.

To solve Eqn. (4.6) (or equivalently, Eqn. (4.8) with t = 1), we used a variation of Newton’s method. We define the residual R_{i,j}^{h} of the approximate solution u^{h} at (θ_{i}, ϕ_{j}) by subtracting the right hand side of (4.6) from the left hand side of (4.6), .i.e.,

(4.12)  R_{i,j}^{h} = det(u_{αβ}^{h}) − K(1 + |∇u^{h}|^{2})^{2}, (θ_{i}, ϕ_{j}) ∈ Ω^{h}.

The goal is to find a u^{h} such that

(4.13)  R_{i,j}^{h} = 0, (θ_{i}, ϕ_{j}) ∈ Ω^{h}.

Eqn. (4.13) is a system of polynomial equations of degree four. The system is linearized about an approximate subsolution, u^{h,0} and the Jacobian matrix is approximated by finite differences. At each Newton step, the solution of this linear system is approximated using a conjugate gradient method for the least squares problem (see, [10, p. 295]). The iterative conjugate gradient method is very fast since our matrix is sparse. We use typically 200 iterations of the conjugate gradient method to approximate the solution to the linear system. Newton’s method typically converges in less than 15 iterations. Our convergence criteria is that the absolute value of the residual be less than 10^{−10}.

5. Numerical Results

In this section, Ω is a subset of the sphere with boundary Γ = ∂Ω = Γ_{1} ∪ Γ_{2}, where Γ_{1} ∩ Γ_{2} = ∅. Γ_{1} is a circle parallel to the xy-plane. Γ_{2} is either a circle parallel to the xy-plane (see Case 5.1-2) or a “sawtooth” curve (see Case 5.3).

Case 5.1

In Case 5.1, we consider domains in the form Ω_{ϑ_{1},ϑ_{2}}, with (ϑ_{1}, ϑ_{2}) = (\frac{1}{3}π, \frac{1}{2}π) and (ϑ_{1}, ϑ_{2}) = (\frac{1}{8}π, \frac{1}{2}π) where

\begin{aligned}Φ & = 1, (θ, ϕ) ∈ Γ = Γ_{1} ∪ Γ_{2}, \\ Γ_{i} & = \{(θ, ϕ) | θ = ϑ_{i}, 0 ≤ ϕ ≤ 2π\}, i = 1, 2. \end{aligned}

In Figure 5.1a, we present a graph of the curve (K, ρ̂(K)) where ρ̂(K) is the sup-norm of ρ^{h} = 1/u^{h}, u^{h} is a solution of (4.6), and ρ is defined in (4.1). In particular,

(5.1)  ρ̂(K) = max_{(θ_{i},ϕ_{j})∈Ω_{ϑ_{1},ϑ_{2}}^{h}} |ρ^{h}(θ_{i}, ϕ_{j}; K, Φ)| = max_{(θ_{i},ϕ_{j})∈Ω_{ϑ_{1},ϑ_{2}}^{h}} \frac{1}{|u^{h}(θ_{i}, ϕ_{j}; K, Φ)|}.

In Figure 5.1b, we present a graph of (K, η̂(K)) where

(5.2)  η̂(K) = ||u^{h}(K, Φ)||_{2}, K ∈ (0, K_{max}),

and u^{h} is a solution of Eqn. (4.6). Figures 5.1 - 5.3 correspond to Case 5.1a where (ϑ_{1}, ϑ_{2}) = (\frac{1}{3}π, \frac{1}{2}π). Figures 5.4 - 5.6 correspond to Case 5.1b where (ϑ_{1}, ϑ_{2}) = (\frac{1}{8}π, \frac{1}{2}π).

Figure 5.1: (K, ρ̂(K)) - Ω\frac{1}{3}π, \frac{1}{2}π Figure 5.2: (K, η̂(K)) - Ω\frac{1}{3}π, \frac{1}{2}π

Figure 5.3a: Radial graphs - Ω\frac{1}{3}π, \frac{1}{2}π, K = −2.99

28

Figure 5.3b: Radial graphs - Ω\frac{1}{3}π, \frac{1}{2}π, K_{max} = 2.63 Figure 5.4: (K, ρ̂(K)) - Ω\frac{1}{8}π, \frac{1}{2}π Figure 5.5: (K, η̂(K)) - Ω\frac{1}{8}π, \frac{1}{2}π

Figure 5.6a: Radial graph - Ω\frac{1}{8}π, \frac{1}{2}π, K = −1.27

30

Case 5.2

In Case 5.2, the boundary of the K-surface consists of two disjoint components formed by taking two small circles and mapping them into curves in R^{3}. In order to define the boundary in this case, let

Φ(θ, ϕ) = r(ϕ)x(θ, ϕ), (θ, ϕ) ∈ Γ = Γ_{1} ∪ Γ_{2},

where

Γ_{1} = {(θ, ϕ) | θ = ϑ_{1}, 0 ≤ ϕ ≤ 2π},
Γ_{2} = {(θ, ϕ) | ϑ_{2} = π − ϑ_{1}, 0 ≤ ϕ ≤ 2π},
r(ϕ) = \sqrt{a^{2} cos^{2} ϕ + b^{2} sin^{2} ϕ},

and a = 1.0, and b = 1.3. Note, the image of Γ under Φ lies inside V_{ϑ_{1}} ∪ V_{π−ϑ_{1}}, where V_{ϑ} is a cone with vertex angle ϑ, i.e., in spherical coordinates V_{ϑ} = {(r, θ, ϕ) | 0 < r < ∞, θ = ϑ, 0 ≤ ϕ ≤ 2π}.

Case 5.2a corresponds to (ϑ_{1}, ϑ_{2}) = (\frac{5}{12}π, \frac{7}{12}π) and the results are presented in Figures 5.7 - 5.9. In Figures 5.10 - 5.12, we present the results of Case 5.2b, (ϑ_{1}, ϑ_{2}) = (\frac{1}{3}π, \frac{2}{3}π). In Figures 5.13 - 5.15, we present the results of Case 5.2c, (ϑ_{1}, ϑ_{2}) = (\frac{1}{6}π, \frac{5}{6}π).

Figure 5.7: (K, ρ̂(K)) - Ω\frac{5}{12}π, \frac{7}{12}π Figure 5.8: (K, η̂(K)) - Ω\frac{5}{12}π, \frac{7}{12}π

32

Figure 5.9a Radial graph - Ω\frac{5}{12}π, \frac{7}{12}π, K = 0.00

33

Figure 5.9b Radial graph - Ω\frac{5}{12}π, \frac{7}{12}π, K_{max} = 1.48 Figure 5.10: (K, ρ̂(K)) - Ω\frac{1}{3}π, \frac{2}{3}π Figure 5.11: (K, η̂(K)) - Ω\frac{1}{3}π, \frac{2}{3}π

35

Figure 5.13: (K, ρ̂(K)) - Ω\frac{1}{6}π, \frac{5}{6}π Figure 5.14: (K, η̂(K)) - Ω\frac{1}{6}π, \frac{5}{6}π

37

Case 5.3

The boundary Γ for Case 5.3 consists of two disconnected components, Γ_{1} and Γ_{2}. Γ_{1} is a small circle in the unit sphere that is parallel to the xy-plane, Γ_{2} is the “sawtooth” component shown in Figure 4.2, and

Φ = 1, (θ, ϕ) ∈ Γ = ∂Ω_{ϑ_{1},ϑ_{2},δ}.

In Case 5.3a, (ϑ_{1}, ϑ_{2}) = (\frac{1}{8}π, \frac{1}{2}π), δ = \frac{1}{4}(\frac{1}{2}π − \frac{1}{8}π) = \frac{3}{32}π, and Ω = Ω_{ϑ_{1},ϑ_{2}, \frac{3}{32} π}. The numerical results are presented in Figures 5.16 - 5.18. The appropriate modifications need to made in the definitions of ρ̂ and û to take into account the sawtooth domain. In Case 5.3b, (ϑ_{1}, ϑ_{2}) = (\frac{1}{3}π, \frac{1}{2}π), δ = \frac{1}{4}(\frac{1}{2}π − \frac{1}{3}π) = \frac{1}{24}π, and Ω = Ω_{ϑ_{1},ϑ_{2}, \frac{1}{24} π}. The numerical results for Case 5.3b are presented in Figures 5.19 - 5.21.

Figure 5.16: (K, ρ̂(K)) - Ω_{ϑ_{1},ϑ_{2}, \frac{3}{32} π} Figure 5.17: (K, η̂(K)) - Ω_{ϑ_{1},ϑ_{2}, \frac{3}{32} π}

39 40

Figure 5.19: (K, ρ̂(K)) - Ω_{ϑ_{1},ϑ_{2}, \frac{1}{24} π} Figure 5.20: (K, η̂(K)) - Ω_{ϑ_{1},ϑ_{2}, \frac{1}{24} π}

42

In Case 5.1, we choose n_{θ} = 30 and n_{ϕ} = 30. In Case 5.2, we choose n_{θ} = 24 and n_{ϕ} = 48. The grid size (h_{θ}, h_{ϕ}) for each case is presented in Table 5.1. All residuals R^{h} are less than 10^{−10} in absolute value. A summary of the results of the radial graph K-surface computations is presented in Table 5.1.

Case  Boundary Data  (h_{θ}, h_{ϕ})  (0.0, K_{max})
5.1a  Ω\frac{1}{3}π, \frac{1}{2}π  (0.017, 0.210)  (0.00, 2.63)
5.1b  Ω\frac{1}{8}π, \frac{1}{2}π  (0.039, 0.210)  (0.00, 1.15)
5.2a  Ω\frac{5}{12}π, \frac{7}{12}π  (0.017, 0.210)  (0.00, 1.48)
5.2b  Ω\frac{1}{3}π, \frac{2}{3}π  (0.035, 0.210)  (0.00, 0.62)
5.2c  Ω\frac{1}{6}π, \frac{5}{6}π  (0.070, 0.210)  (0.00, 0.14)
5.3a  Ω\frac{1}{8}π, \frac{1}{2}π, \frac{3}{32}π  (0.049, 0.131)  (0.00, 1.23)^{∗}
5.3b  Ω\frac{1}{3}π, \frac{1}{2}π, \frac{1}{24}π  (0.022, 0.131)  (0.00, 2.90)^{∗}

Table 5.1: Maximal Intervals for which a radial graph K-surface exists

* The values of K_{max} for these cases applicable to sawtooth-like domains Ω′ ⊂ Ω where Ω′ is obtained by smoothing the corners of Ω. We found that by using a finer grid (h_{θ}, h_{ϕ}) = (0.011, 0.0655), the maximum Gauss curvature in Case 5.3b was reduced from 2.9 to 1.9. We conjecture that in the limit as (h_{θ}, h_{ϕ}) → (0, 0), K_{max} → 1. See the discussion in section 6.

Remark. Although the analytical theory is concerned only with surfaces of positive Gauss curvature, we were able to compute solutions of negative Gauss curvature in Case 5.1a, where we found a minimum value of the Gauss curvature to be −2.99. In Case 5.1b, we found a minimum value of the Gauss curvature to be −1.27. Our algorithm seemed to work well for these cases. For Cases 5.2 and 5.3, we tried to compute surfaces of negative curvature, but the algorithm failed to converge in these cases.

6. Concluding Remarks

The continuity method is a fundamental tool that is used to prove the solvability of a certain class of Dirichlet problems (see, e.g., [3] and [9]). The continuity method can be adapted to a scheme that can be used to compute K-surfaces. The solvability of Problems (P 1) and (P 2) essentially reduce to the existence of a subsolution. Computationally, if we can find an approximate subsolution, we can compute the corresponding K-surface. This suggests the possibility of establishing a rigorous convergence proof for numerical boundary value problems for K-surfaces.

When computing K-surfaces for values of K that were not close to K_{max}, our algorithm was robust. The initial subsolution did not need to be “close” to the K-surface in order for the numerical scheme to converge. However, the algorithm was sensitive to the initial guess when K was near K_{max}. This was not a serious problem, because we were able to compute K-surfaces for increasing values of K without any difficulty as long as the step size for ΔK was not too large. When K = 1 and the K-surface is a subset of the unit sphere (see Case 3.1, Case 5.1, Case 5.3), we find that the residuals are on the order of 10^{−10} in absolute value.

Given a solution (K_{1}, u(K_{1})), the natural way to compute a family of K-surfaces is to use u(K_{1}) as the subsolution for the continuity method for every 0 < K < K_{1}. However, we also found that we were able to compute a sequence of solutions for increasing values of K. In particular, we were able to compute (K_{k+1}, u^{h}(K_{k+1})) for K_{k+1} > K_{k} by using u^{h}(K_{k}) as the initial guess for the Newton iterations for K_{k+1} = K_{k} + ΔK. This could be done if the step size ΔK = K_{k+1} − K_{k} was not too large. For our applications, we chose ΔK = 0.01.

We are able to determine the accuracy of our computed values of K_{max} in Case 5.1, because the corresponding boundary value problem can be solved exactly and the true value of K_{max} can be computed. In Case 5.1, the problem of determining a K-surface reduces to solving an ordinary differential equation (see, Eqn. (B.4)). K_{max} can be computed by solving a system of scalar equations (see Eqn. (B.8) − (B.11)). The derivation of K_{max} for Case 5.1 is outlined in Appendix B and the results are summarized in Table 6.1, indicating that the computed values of K_{max} are within 0.03 of the true solution.

Case Boundary data  Computed K_{max} True K_{max}
5.1a  Ω\frac{1}{3}π, \frac{1}{2}π  2.63  2.608
5.1b  Ω\frac{1}{8}π, \frac{1}{2}π  1.15  1.147
Table 6.1: Computed K_{max} vs. true K_{max}

As was pointed out earlier, convergence of our algorithm hinged on the existence of a numerical subsolution. We carried out a number of numerical experiments to investigate what happens when u^{h,0} is not convex. Under certain conditions in Problem (P 1), we were able to compute a K-surface starting with a nonconvex u^{h,0}. For these cases K was not near K_{max}. We formed a nonconvex u^{h,0}(x_{i}, y_{j}) by adding a perturbation to the approximate subsolution \underline{u}^{h}, i.e.,

u^{h,0}(x_{i}, y_{j}) = \underline{u}_{i,j}^{h} + ∑_{i,j∈I} c_{i,j},

where I is a collection of indices for a number of internal grid points and |c_{i,j}| ≤ 0.025. Under these circumstances, the numerical continuity method would smooth out the perturbations, so that the resulting solution u^{h} was convex. We emphasize that the conditions K << K_{max} and |c_{i,j}| << 1 were important. If these conditions were violated, spurious solution could appear.

When K is large, we observed that the algorithm at times could converge to a spurious solution. That is, the residual for the solution is small, yet the surface is clearly not convex. In Figure 6.1, we present a K-surface that was computed using Case 5.1 boundary conditions and K = 1.03. The residual for this solution is 10^{−10}. However, the close-up of the K-surface presented in Figure 6.2 shows a dimple near the boundary Γ_{1}, illustrating that the solution is not convex. This happened when the initial guess for the Newton iteration was too far away from the solution. Initially, the residual increased, but did not diverge. Then we observed the residuals become smaller than our prescribed residual. Thus, even though the residuals were small, dimples such as the one shown in Figure 6.2 can develop in the solution.

Remark. The dimpling effect we have just described need not be localized near the boundary. We have observed also dimpling of numerical solutions over the interior of Ω^{h}. We should point out that the appearance of a dimple in a solution was atypical. If a solution exists for a certain value of K, then the algorithm would converge rapidly to that solution and the K-surface was convex. Only with a poor choice of an approximate subsolution were we able to produce these spurious solutions.

Figure 6.1: A spurious solution Figure 6.2: A close-up of the dimple in Figure 6.1

In an effort to accelerate the method by which we numerically tracked solution curves, we considered employing a pseudo-arc-length continuation method (PCM) (see [1] and [2]). In this method, the nonlinear equations are supplemented by a pseudo-arc-length condition, and the solution curve is tracked by first taking an Euler step in the direction of the tangent, and then solving the supplemented nonlinear system. For a variety of problems, it offers the advantages of being able to track solution curves around turning points and allows one to take larger step sizes. For the applications that we considered here, the capability of tracking curves around turning points is not necessary, because if a K-surface exists for K = K′, then it is the unique K-surface for that value of K′. This means there are no turning points. However, when the PCM is applied to an equation in the form (2.7) and an Euler step is taken, there is no assurance that the resulting guess will yield a good subsolution, unless the step size is small. If the step size is small, then there is no advantage to introducing the pseudo-arc-length parameter. If one compares the solution curve (K, η̂(K)) such as in Figure 3.2 with one generated by using a pseudo-arc-length continuation method (with a large step size), one notices that the solution curve generated by pseudo-arc-length continuation gradually drifts away from the curve (K, η̂(K)) and in fact turns around! These are not true solutions even though they have small residuals. We should point out that if the step size is sufficiently small, the turn-around exhibited by the PCM does not occur.

As K → K_{max} for a family of solutions K-surfaces (K, u_{K}), the corresponding K-surfaces are characterized by the behavior of their tangent planes near the boundary Γ. For a K-surface parametrized as a graph over a rectangle, the tangent planes become nearly vertical (see, e.g., Figures 3.4d, 3.7d, 3.10d, 3.13d, 3.16d, 3.19d). This behavior cannot be attributed to the sharp corners in the rectangle, since this is observed along the smooth components of the boundary. A very similar behavior is observed for radial graphs. In Figures 5.3b and 5.6b, one can see that tangent planes near the boundary Γ_{1} become nearly parallel to the xy-plane. In fact, this conditions is used in Appendix B to compute the true value of K_{max} (see Eqn. (B.11)). For the boundaries used in Case 5.2 (see Figures 5.9b and 5.12b), we see that there are four areas near the boundary Γ where the tangent plane becomes nearly parallel to the xy-plane (two of the regions on Γ_{2} cannot be seen). Because the mesh size in Figure 5.15b is larger than those in Figures 5.9b and 5.12b, one would expect to see small regions with horizontal tangent planes on a finer grid. This suggests the possiblity of establishing an estimate of K_{max} based on the boundary data (Γ, Φ).

Certain boundary conditions in Sections 3 and 5 did not satsify the smoothness properties that are assumed in the analytical work. However, this is not a serious problem. In particular, if one considered (P 1) on a domain Ω′ ⊂ Ω, where Ω′ is obtained by smoothing the corners of Ω, one would expect to find results similar to those presented in Section 3. We feel that the results of Section 3 are reasonable estimates of K_{max} for K-surfaces over such domains. The numerical results for Case 5.1 and Case 5.2 can be applied directly to problem (P 2) because the boundaries are smooth. The case of the sawtooth boundary is a more delicate one, because of the corner points. We should interpret the results for K_{max} as applicable to problem (P 2) for domains Ω′ ⊂ Ω that are obtained by smoothing the corners of Ω.

We graphed all the K-surfaces in Sections 3 and 5 over Ω^{h} ∪ ∂Ω^{h}. In doing so, the K-surfaces appear to have a singular behavior near the boundary. In Figure 6.3a, we present a K-surface with boundary conditions similar to those used in Case 5.3b (a finer grid was used also). The solution u^{h} on Ω^{h} with boundary data Φ = 1 on Γ = ∂Ω^{h} are plotted. One notices that the K-surface appears to be nonconvex in a band adjacent to the boundary Γ_{2} (See Figure 6.3a). This becomes more apparent when we focus on the area near a corner of the sawtooth boundary (see Figure 6.3b). On closer examination one observes that if we restrict our attention to the surface defined by u^{h} on Ω^{h}, then the corresponding K-surface appears convex (See Figure 6.4). This is not entirely unexpected given that the analytical problem defines K on Ω, and the approximate problem (4.6) defines K on Ω^{h}. On the other hand, the width of the band can be made smaller by refining the grid size.

In the limit, one would expect to obtain a more accurate estimate of K_{max}. However, based on the computations in Case 5.3, we cannot exclude the possibility that K_{max} = 1 for the actual nonsmooth sawtooth boundary.

However, refining the mesh size reduces the corresponding computed value of K_{max}. In particular, we found that be reducing (h_{θ}, h_{ϕ}) from () to (0.011, 0.0655), the maximum Gauss curvature changed from 2.9 to 1.9. We conjecture that as (h_{θ}, h_{ϕ}) approaches (0, 0), K_{max} approaches 1.

Figures 6.3a: Radial graph over Ω^{h} ∪ ∂Ω^{h} Figures 6.3b: Magnification of corner region of Figure 6.3a.

50

Figures 6.4: Radial graph over Ω^{h} Appendix A. Discrete Approximations

In this section, we present the discrete versions of (P 1) and (P 2), respectively. The discrete problems are obtained by replacing partial derivatives by their finite difference approximations and discretizing the domains.

First, we consider the approximation for a K-surface that is a graph over a rectangle Ω = R_{a,b}. Let x_{i} = −a+(i−1)h_{x} for i = 1, …, n_{x} and let y_{j} = −b+(j −1)h_{y} for j = 1, …, n_{y}. The grid size is characterized by h = (h_{x}, h_{y}) where the size of a typical rectangle in the grid is h_{x} × h_{y} and (h_{x}, h_{y}) = (2a/n_{x}, 2b/n_{y}). The domain Ω^{h} is given by

Ω^{h} = {(x_{i}, y_{j}) | i = 2, …, n_{x} − 1, j = 2, …, n_{y} − 1} .

The boundary ∂Ω^{h} is the set

∂Ω^{h} = {(x_{i}, y_{j}) | i ∈ {1, n_{x}} and 1 ≤ j ≤ n_{y} or 1 ≤ i ≤ n_{x} and j ∈ {1, n_{y}}} .

If u(x_{i}, y_{j}) denotes the value of the true solution at (x_{i}, y_{j}), the value of its approximation is given by u_{i,j}^{h} or u^{h}(x_{i}, y_{j}), i.e.,

u(x_{i}, y_{i}) ˙= u_{i,j}^{h}(x_{i}, y_{j}) = u_{i,j}^{h}

The finite difference approximations are as follows:

(A.1a)  \begin{aligned}u_{x}^{h}(x_{i}, y_{j}) & = \frac{u_{i+1,j}^{h} − u_{i−1,j}^{h}}{2h_{x}} , \\ u_{y}^{h}(x_{i}, y_{j}) & = \frac{u_{i,j+1}^{h} − u_{i,j−1}^{h}}{2h_{y}} , \\ u_{xx}^{h}(x_{i}, y_{j}) & = \frac{u_{i+1,j}^{h} − 2u_{i,j}^{h} + u_{i−1,j}^{h}}{h_{x}^{2}} , \\ u_{yy}^{h}(x_{i}, y_{j}) & = \frac{u_{i,j+1}^{h} − 2u_{i,j}^{h} + u_{i,j−1}^{h}}{h_{y}^{2}} , \\ u_{xy}^{h}(x_{i}, y_{j}) = u_{yx}^{h}(x_{i}, y_{j}) & = \frac{u_{i+1,j}^{h} − u_{i−1,j+1}^{h} − u_{i+1,j−1}^{h} + u_{i−1,j−1}^{h}}{4h_{x}h_{y}} . \end{aligned}

It follows that

(A.1b)  det(u_{αβ}^{h}) = u_{xx}^{h}(x_{i}, y_{j})u_{yy}^{h}(x_{i}, y_{j}) − (u_{xy}^{h}(x_{i}, y_{j}))^{2},

and

(A.1c)  |∇u_{i,j}^{h}(x_{i}, y_{j})|^{2} = |u_{x}^{h}(x_{i}, y_{j})|^{2} + |u_{y}^{h}(x_{i}, y_{j})|^{2}.

The finite difference approximation to (P 1) is given by

(A.2)  u_{xx}^{h}(x_{i}, y_{j})u_{yy}^{h}(x_{i}, y_{j}) − (u_{xy}^{h}(x_{i}, y_{j}))^{2} = K(1 + |u_{x}^{h}(x_{i}, y_{j})|^{2} + |u_{y}^{h}(x_{i}, y_{j})|^{2})^{2},

for i = 2, …, n_{x} − 1, j = 2, …, n_{y} − 1. Eqn. (A.2) is equivalent to Eqn. (2.5). Equations (A.2) are a system of (n_{x} − 2) × (n_{y} − 2) fourth-order polynomial equations in the unknowns u_{i,j}^{h}. For the cases discussed in Sections 2-3, we chose n_{x} = 25 and n_{y} = 25. The total number of equations is 529.

In order to describe how well the numerical solution satisfies the discrete equations, we define the residual R_{i,j}^{h} of a solution u^{h} at (x_{i}, y_{j}) as follows,

(A.3)  R_{i,j}^{h} = u_{xx}^{h}(x_{i}, y_{j})u_{yy}^{h}(x_{i}, y_{j}) − (u_{xy}^{h}(x_{i}, y_{j}))^{2} − K(1 + |u_{x}^{h}(x_{i}, y_{j})|^{2} + |u_{y}^{h}(x_{i}, y_{j})|^{2})^{2},

where 2 ≤ i ≤ n_{x} − 1, 2 ≤ j ≤ n_{y} − 1. Note, R_{i,j}^{h} depends on the value of the numerical solution at the 8 grid points that surround (x_{i}, y_{j}) through the relations (A.1). The residual of the solution is denoted by

(A.4)  R^{h} = max_{(i,j)∈J} |R_{i,j}|.

where J = {(i, j) | 2 ≤ i ≤ n_{x} − 1, 2 ≤ j ≤ n_{y} − 1}.

For surfaces that are expressible as radial graphs, we can derive a similar finite difference approximation. We begin by describing the regions considered in Cases 5.1 - 5.2. Let h = (h_{θ}, h_{ϕ}). Let θ_{i} = ϑ_{1} + (i − 1)h_{θ} for i = 1, 2, …, n_{θ} where h_{θ} = (ϑ_{2} − ϑ_{1})/(n_{θ} − 1). Let ϕ_{j} = (j − 1)h_{ϕ} for j = 1, 2, …, n_{ϕ} − 1 where h_{ϕ} = 2π/n_{ϕ}. For a domain that is a subset of the sphere, Ω^{h} consists of all (θ_{i}, ϕ_{j}),

Ω^{h} = {(θ_{i}, ϕ_{j}) | i = 2, …, n_{θ} − 1, j = 1, …, n_{ϕ} − 1}.

The boundary ∂Ω^{h} is the set

∂Ω^{h} = {(θ_{i}, ϕ_{j}) | i ∈ {1, n_{θ}} and 1 ≤ j ≤ n_{ϕ} − 1} .

The approximate solution at the points (θ_{i}, ϕ_{j}) is denoted by u_{i,j}^{h}, i.e.,

u(θ_{i}, ϕ_{i}) ˙= u_{i,j}^{h}(θ_{i}, ϕ_{j}) = u_{i,j}^{h}.

The finite difference approximations of the derivatives are given by

(A.5a)

u_{θ}^{h}(θ_{i}, ϕ_{j}) = \frac{u_{i+1,j}^{h} − u_{i−1,j}^{h}}{2h_{θ}} ,
u_{ϕ}^{h}(θ_{i}, ϕ_{j}) = \frac{u_{i,j+1}^{h} − u_{i,j−1}^{h}}{2h_{ϕ}} ,
∇_{11}u^{h}(θ_{i}, ϕ_{j}) = \frac{u_{i+1,j}^{h} − 2u_{i,j}^{h} + u_{i−1,j}^{h}}{h_{θ}^{2}} ,
∇_{12}u^{h}(θ_{i}, ϕ_{j}) = ∇_{21}u^{h}(θ_{i}, ϕ_{j}) = \frac{u_{i+1,j}^{h} − u_{i−1,j+1}^{h} − u_{i+1,j−1}^{h} + u_{i−1,j−1}^{h}}{4h_{θ}h_{ϕ}} − cot θ_{i}u_{ϕ}^{h}(θ_{i}, ϕ_{j}),
∇_{22}u^{h}(θ_{i}, ϕ_{j}) = \frac{u_{i,j+1}^{h} − 2u_{i,j}^{h} + u_{i,j−1}^{h}}{h_{ϕ}^{2}} + sin θ_{i} cos θ_{i}u_{θ}^{h}(θ_{i}, ϕ_{j}),

respectively. The approximations to the quantities u_{11}, u_{12}, u_{21}, u_{22} are given by

(A.5b)  \begin{aligned}u_{11}^{h}(θ_{i}, ϕ_{j}) & = ∇_{11}u^{h}(θ_{i}, ϕ_{j}) + u_{i,j}^{h}, \\ u_{12}^{h}(θ_{i}, ϕ_{j}) = u_{21}^{h}(θ_{i}, ϕ_{j}) & = ∇_{12}u^{h}(θ_{i}, ϕ_{j}), \\ u_{22}^{h}(θ_{i}, ϕ_{j}) & = ∇_{22}u^{h}(θ_{i}, ϕ_{j}) + sin^{2}θ_{i}u_{i,j}^{h}, \end{aligned}

respectively. The above definitions hold for all i = 2, …, n_{θ} − 1 and j = 1, …, n_{ϕ} − 1 with the obvious modifications to ensure periodicity in ϕ.

The finite difference approximation to (P 2) for Cases 5.1 and 5.2 is given by

(A.6) (sin θ_{i})^{−2}(u_{11}^{h}(θ_{i}, ϕ_{j})u_{22}^{h}(θ_{i}, ϕ_{j}) − (u_{12}^{h}(θ_{i}, ϕ_{j}))^{2}) = K(1 + |u_{θ}^{h}(θ_{i}, ϕ_{j})|^{2} + (sin θ_{i})^{−2}|u_{ϕ}^{h}(θ_{i}, ϕ_{j})|^{2})^{2},

where i = 2, …, n_{θ} − 1 and j = 1, …, n_{ϕ} − 1. Equations (A.6) and (4.6) are equivalent.

The sawtooth domain and its boundary (see Case 5.3) is cumbersome to write out explicitly. To form Ω_{ϑ_{1},ϑ_{2},δ}^{h}, the domain points (θ_{i}, ϕ_{j}) that fall inside any of the B_{k}’s are removed according to the definitions in Eqn. (4.11) (see Figure 4.1b).

Appendix B. Hyperbolic spherical surfaces of revolution

Surfaces of revolution with constant Gauss curvature can be determined by solving an ordinary differential equation. Depending on the values of certain parameters, these surfaces are called hyperbolic spherical, elliptic spherical, or pseudo-spherical surfaces of revolution (see [13]). The examples in Case 5.1 turn out to be hyperbolic spherical surfaces of revolution and so we can compute the true value of K_{max} in these cases.

For a surface of revolution,

(B.1)  X(u, v) = (ϕ(v) cos u, ϕ(v) sin u, ψ(v)), 0 < u < 2π, v_{1} < v < v_{2},

it follows that the Gauss curvature is given by (see [6])

(B.2)  K = −\frac{ϕ′′(v)}{ϕ(v)},

where v is arc length, i.e.,

(B.3)  |ϕ′(v)|^{2} + |ψ′(v)|^{2} = 1.

When K is a constant, we see that Eqn. (B.2) yields the ordinary differential equation,

(B.4)  ϕ′′ + Kϕ = 0,

whose general solution is

(B.5)  ϕ(v) = A cos \sqrt{K}v + B sin \sqrt{K}v.

From Eqn. (B.3), it follows that

(B.6)  ψ(v) − ψ(v_{1}) = ∫_{v_{1}}^{v} \sqrt{1 − |ϕ′(v)|^{2}}dv.

For the examples considered in Case 5.1, the radius of the circle Γ_{2} is 1. Let r^{∗} denote the radius of the circle Γ_{1} and let z^{∗} denote the distance from Γ_{1} to the xy-plane. It follows that

(B.7)  \begin{matrix}r^{∗} = & cos(\frac{1}{2}π − ϑ_{1}), \\ z^{∗} = & sin(\frac{1}{2}π − ϑ_{1}). \end{matrix}

The parameters A and B introduced in Eqn. (B.5) are determined by the conditions (B.7). In particular, v_{1} = 0 and ϕ(0) = 1 implies that B = 1. At v = v_{2}, we have

(B.8)  r^{∗} = A sin \sqrt{K}v_{2} + cos \sqrt{K}v_{2}.

Solving for A, we find

(B.9)  A = \frac{r^{∗} − cos \sqrt{K}v_{2}}{sin \sqrt{K}v_{2}} .
Finally, ψ(v_{1}) = 0 and ψ(v_{2}) = z^{∗}, so that
(B.10)  z^{∗} = ∫_{0}^{v_{2}} \sqrt{1 − K(A cos \sqrt{K}v − sin \sqrt{K}v)^{2}}dv.

Eqn. (B.10) defines a one-parameter family of K-surfaces. For every K ∈ (0, K_{max}), one can solve Eqn. (B.9) -(B.10) and determine the arc length of the corresponding generating curve, i.e., v_{2} = v_{2}(K). To determine the maximum allowable Gauss curvature K_{max}, we observe that if ψ′(v_{2}) < 0, the corresponding surface of revolution would have a band of negative curvature. Thus, K_{max} occurs when ψ′(v_{2}) = 0. It follows from (B.5) − (B.6) and (B.8) − (B.9) that we need to find v_{2} such that

(B.11)  1 − K(A cos \sqrt{K}v_{2} − sin \sqrt{K}v_{2})^{2} = 0.

Solving Equations (B.9)−(B.11), uniquely determines K = K_{max} and the corresponding arc length v_{2} = v_{max}. For Case 5.1a, we find that the computed value of K_{max} is 2.63, compared to the true value of 2.608. For Case 5.1b, we find that the computed value of K_{max} is 1.15, compared to the true value of 1.147. These results are summarized in Table 6.1. The case of a pseudo-spherical surface of revolution can be handled in a similar manner, with the condition ψ′(v_{2}) replaced by ψ′(0) = 0.

Acknowledgement. The authors wish to thank Bo Guan and Joel Spruck for a number of helpful discussions.

References

[1] F. Baginski, The computation of one-parameter families of bifurcating elastic surfaces, SIAM J. Appl. Math (to appear).

[2] J. H. Bolstad & H. B. Keller, A multigrid continuation method for elliptic problems with folds, SIAM J. Sci. Stat. Comput. 7 (4) (1986), 1081–1104.

[3] L. Caffarelli, L. Nirenberg, & J. Spruck, The Dirichlet problem for nonlinear second order elliptic equations I: Monge-Ampére equations, Comm. Pure Appl. Math 37 (1984), 369–402.

[4] L. Caffarelli, J.J. Kohn, L. Nirenberg, & J. Spruck, The Dirichlet problem for nonlinear second order elliptic equations II: Complex Monge-Ampére and uniformally elliptic equations, Comm. Pure Appl. Math 38 (1985), 209–252.

[5] L. Caffarelli, L. Nirenberg, & J. Spruck, The Dirichlet problem for nonlinear second order elliptic equations III: Functions of the eigenvalues of the Hessian, Acta Mathematica 155 (1985), 261–301.

[6] Manfredo P. Do Carmo, Differential Geometry of Curves and Surfaces, Prentice-Hall, Englewood Cliffs, New Jersey, 1976.

[7] Bo Guan, Boundary Value Problems for Surfaces of Prescribed Gauss Curvature and Monge-Ampere Equations, Thesis, University of Massachusetts, Amherst, Massachusetts, September 1992.

[8] Bo Guan, The Dirichlet problem for a class of fully nonlinear elliptic equations, (preprint).

[9] B. Guan & J. Spruck, Boundary value problems on S^{n} for surfaces of constant Gauss curvature, Annals of Math. (to appear).

[10] M. Hestenes, Conjugate Direction Methods in Optimization, Springer-Verlag, 1980.

[11] D. Hoffman, H. Rosenberg, & J. Spruck, Boundary value problems for surfaces of constant Gauss curvature, Comm. Pure Appl. Math. 45 (1992), 1051-1062.

[12], IMSL Stat/Library User’s Manual, Version 1.1, International Mathematical and Statistical Libraries, Houston, 1989.

[14] N.V. Krylov, Boundedly nonhomogeneous elliptic and parabolic equations in a domain, Izvestia. Math. Ser. 47 (1983), 75-108.

[15] N.V. Krylov, Nonlinear Elliptic and Parabolic Equations of the Second Order, Moscow Nauk Glavnaya Red. Fisico-Math. Lit., 1983.

[13] E. Kreyszig, An Introduction to Differential Geometry and Riemannian Geometry, University of Toronto Press, 1968.

[16] H. Rosenberg, Hypersurfaces of Constant Curvature in Space Forms (preprint).