Holomorphic potentials, symplectic integrators and CMC surfaces

David Lerner and Ivan Sterling ^{∗}

July 30, 1995

1 Introduction

The Gauss map of a surface Σ in R^{3} sends a point on the surface to the corresponding unit normal vector n ∈ S^{2}. In terms of a conformal coordinate z on the surface, if the surface has constant mean curvature (CMC) then its Gauss map is harmonic [RV]:

∂_{zz̄}^{2}n = \frac{2n̄}{1 + nn̄}∂_{z}n∂z̄n.  (1.1)

Moreover, given a harmonic map n from the simply connected domain D to S^{2}, one can construct a (possibly branched) CMC immersion with n as its Gauss map [K]. For a fixed value of the mean curvature, the construction is unique up to rigid motions of R^{3} and conformal transformations of the surface. So, to construct all CMC surfaces with domain D, it suffices to construct all harmonic maps from D to S^{2}. The recently discovered method of holomorphic potentials due to Dorfmeister, Pedit and Wu [DPW] gives a method for accomplishing this which is, in principle, constructive. The procedure, described below, involves first solving a differential equation with coefficients in a twisted loop group, followed by a “Gram-Schmidt” factorization of the resulting loop. Implementing these two fundamentally infinitedimensional processes on the computer requires further theoretical and numerical development, which forms the subject of this paper. Prior work on this subject has either been too CPU intensive [LS] or of such questionable accuracy that the computer renditions have been suspect. Our goal here is to demonstrate the existence of a class of numerical methods that may, in the future, allow us to accurately construct and visualize these surfaces in a routine manner.

∗Partially supported by NSF grant DMS-9204801, and KITCS at the University of Kansas

Following the introductory material in Sections 2 and 3, we establish the notion of a “perfect” integrator in Section 4, and show that these numerical schemes permit subsequent construction of the surface in an amount of time which is polynomial (rather than exponential) in the number of grid points. In Section 5, we give some focus to the computational possibilities by viewing the system of differential equations as a (complex) Hamiltonian system. Any perfect integrator is necessarily symplectic [S] in the usual sense, but the converse is false. We then establish the existence of these integrators in a number of important cases, including second and higher order methods when the holomorphic potential takes a particular, common form.

The existence of perfect integrators of sufficiently high accuracy makes it possible, in principle, to use the method of holomorphic potentials to construct explicit examples of smoothly immersed CMC surfaces of higher genus. This is now under investigation.

We would like to thank U. Pinkall and C. Gunn; the results in Section 5 were motivated by their (unpublished) theoretical and computer algorithms.

2 CMC Surfaces

Let Σ be a smoothly immersed surface in R^{3}; let D be a simply connected domain in R^{2}, and let f : D → Σ be a conformal chart on the surface. If z = u + iv is a complex coordinate in D, then the induced metric in Σ|_{f(D)} has the form ds^{2} = e^{ω}(du^{2} + dv^{2}) = e^{ω}dzd̄z. We define quantities Q and H by

Q = f_{zz} · n, and \frac{1}{2}He^{ω} = f_{zz̄} · n,  (2.1)

where n is the surface normal, and subscripts denote partial derivatives. If ∇ is the covariant derivative in Σ, and X is tangent to the surface, the Weingarten map W : X → ∇_{X}n is given in the (u, v) coordinates by

W = ( \begin{matrix}H + 2e^{−ω}Re(Q) & −2e^{−ω}Im(Q) \\ −2e^{−ω}Im(Q) & H − 2e^{−ω}Re(Q) \end{matrix} ) .  (2.2)

The mean curvature H is then half the trace of W. The Gauss-Codazzi equations of the surface in this local coordinate system are

Qz̄ = \frac{1}{2}H_{z}e^{ω}  (2.3)
ω_{zz̄} + \frac{H^{2}}{2} e^{ω} − 2|Q|^{2}e^{−ω} = 0.  (2.4)

If H is constant, then equation (2.3) implies that Q is holomorphic, and the above equations are equivalent to the harmonic map equation (1.1) for n. In this case, Qdz^{2} is the local representative of a holomorphic quadratic differential on Σ, known as the Hopf differential. At the zeros of this differential, known as umbilics, W becomes a scalar multiple of I. It is well-known that in a neighborhood of any nonumbilic point, a new complex coordinate exists with respect to which (2.4) reduces to the elliptic sinh-Gordon equation (see, for example Eisenhart [E]). However, on a compact, smoothly immersed CMC surface of genus g, the Riemann-Roch theorem gives a net count of 4g −4 zeros for the Hopf differential, so any CMC surface of g ≥ 2 necessarily has umbilics. Since one of our ultimate goals is the construction of CMC surfaces with arbitrary genus, we shall restrict ourselves to the equations above.

It will be convenient to identify R^{3} with the Lie algebra of su(2). If σ = (σ_{1}, σ_{2}, σ_{3}) are the Pauli matrices, let J = −iσ, and identify x ∈ R^{3} with

X = x · J = ( \begin{matrix}−ix_{3} & −x_{2} − ix_{1} \\ x_{2} − ix_{1} & ix_{3} \end{matrix} ) .  (2.5)

This is an isometry under the ad-invariant metric (X, Y) = \frac{1}{2}T r(XY). The unit sphere is now identified with the orbit of J_{3} under the adjoint action of SU(2); hence we regard it as the symmetric space SU(2)/K, where

K = \{( \begin{matrix}e^{iθ} & 0 \\ 0 & e^{−iθ} \end{matrix} ) : θ ∈ (0, 2π]\} .  (2.6)

The corresponding Cartan decomposition of su(2) as K ⊕ P is given by

K = \{( \begin{matrix}ir & 0 \\ 0 & −ir \end{matrix} ) : r ∈ R\}  (2.7)
P = \{( \begin{matrix}0 & b \\ −b̄ & 0 \end{matrix} ) : b ∈ C\} .  (2.8)

For f : D → R^{3} as above, let F_{u}, N_{z},… denote the corresponding matrices in su(2). Since D is simply connected, the Gauss map N can be lifted to an Ñ : D → SU(2), which satisfies

(adÑ)(z, z̄) : J_{3} → N(z, z̄).  (2.9)

If N is harmonic, then Ñ will be called a harmonic lift. Notice that [J_{3}, K] = 0 implies that (adÑ)(z, z̄) · J_{3} is unchanged under the replacement of Ñ by Ñ · k, where k : D → K is any smooth map. In particular, Ñ · k is another harmonic lift of N. Each such lift maps the ordered pair (J_{1}, J_{2}) to an orthonormal frame at each point of

the immersion. Moreover, if k(z, z̄) = ( \begin{matrix}e^{iθ(z,z̄)} & 0 \\ 0 & e^{−iθ(z,z̄)} \end{matrix} ), it is easily checked that multiplying Ñ by k rotates this frame through the angle 2θ(z, z̄) at the point labeled by (z, z̄). So we have

Proposition 2.1 Under the above circumstances, the Gauss map of the surface has a smooth lift Φ to SU(2) satisfying

(adΦ)(z, z̄) : J → \{e^{− \frac{ω}{2}} F_{u}, e^{− \frac{ω}{2}} F_{v}, N\}.

Φ is unique up to multiplication by −1.

We shall refer to Φ as a normalized lift of the Gauss map. Clearly it contains all the information to be found in the moving frame; in particular, we recover the Gauss-Codazzi equations (2.3),(2.4) as the integrability conditions for the overdetermined linear system

Φ^{−1}∂_{z}Φ = U = ( \begin{matrix}−ω_{z}/4 & Qe^{− \frac{ω}{2}} \\ −\frac{H}{2} e^{\frac{ω}{2}} & ω_{z}/4 \end{matrix} )  (2.10)

(2.11)

Φ^{−1}∂z̄Φ = V = −Ū^{t}.  (2.12)

If H is a constant, as we shall suppose from now on, rescaling the metric in R^{3} by a (constant) factor r^{−2} rescales H by the factor r. So without loss of generality, we can take H ≡ 2.

When Φ is replaced by Φ·k, the quantities U and V undergo the corresponding partial gauge transformations; in particular,

U → Û := k^{−1}Uk + k^{−1}∂_{z}k.  (2.13)

If we write k = ( \begin{matrix}e^{i \frac{τ}{4}} & 0 \\ 0 & e^{−i \frac{τ}{4}} \end{matrix} ), and define ω̂ = ω + iτ, then we find that:

Û = ( \begin{matrix}−¯ω̂_{z}/4 & Qe^{− \frac{ω̂}{2}} \\ −e^{\frac{ω̂}{2}} & ¯ω̂_{z}/4 \end{matrix} ) .  (2.14)

We therefore have

Proposition 2.2 Φ is a normalized harmonic lift if and only if U_{21} is real and negative. The quantity Q = −U_{12}U_{21} is gauge-invariant and can be computed from any lift of N.

3 The Scattering System

It is a classical result that any CMC surface is naturally embedded in a oneparameter family of such surfaces. From the soliton theoretic viewpoint, this parameter is (essentially) the spectral parameter in a scattering system: For any map Φ : D → SU(2), the su(2)-valued 1-form Φ^{−1}dΦ decomposes as α_{K} + α_{P} in K ⊕ P. If we further decompose α_{P} into its (1,0) and (0,1) parts, we may form the loop (for |λ| = 1)

α_{λ} = λ^{−1}U_{−1}dz + α_{K} + λU_{1}dz̄.  (3.1)

The α_{λ} may be regarded as a one-parameter family of connections in the trivial twoplane bundle over the surface Σ. Clearly α_{1} is flat, since by construction it has the form Φ^{−1}dΦ. Conversely, we have the following fundamental result, due originally to Pohlmeyer [Po] in a somewhat different context:

Proposition 3.1 If α_{λ} has the form given in (3.1), then dα_{λ} + \frac{1}{2}[α_{λ}, α_{λ}] = 0 (α_{λ} is flat) for all λ ∈ S^{1} if and only if there is a one-parameter family of harmonic lifts ϕ_{λ} such that ϕ_{λ}^{−1}dϕ_{λ} = α_{λ}.

Proof: (We shall use the notation ϕ_{λ} whenever we are thinking of λ as a fixed parameter; for example, dϕ_{λ} = ∂_{z}ϕdz + ∂z̄ϕdz̄ contains no derivatives with respect to λ.) Write α_{K} = U_{0}dz + V_{0}dz̄, and let functions a, b, and c be defined by

U_{−1} = ( \begin{matrix}0 & b \\ c & 0 \end{matrix} ) , U_{0} = ( \begin{matrix}a & 0 \\ 0 & −a \end{matrix} ) .

Since D is simply connected, the loop of connections is flat if and only if there exists a unitary loop ϕ_{λ}(z, z̄) such that

ϕ_{λ}^{−1}∂ϕ_{λ} = λ^{−1}U_{−1} + U_{0}, and ϕ_{λ}^{−1} ¯∂ϕ_{λ} = λV_{1} + V_{0},

where V_{1} = −Ū_{−1}^{t} and V_{0} = −Ū_{0}^{t}. We must show that the Gauss-Codazzi equations are satisfied.

The integrability condition [∂, ¯∂] = 0 now reads, in terms of the components,

∂z̄b − 2āb = 0  (3.2)
∂z̄c + 2āc = 0  (3.3)
∂z̄a + ∂_{z}ā + |b|^{2} − |c|^{2} = 0.  (3.4)

From the first two, we see that c∂z̄b + b∂z̄c = 0, so that the product bc is holomorphic in z. Introducing the functions μ and Q via

b = Qe^{− \frac{μ}{2}} , c = −e^{\frac{μ}{2}} ,

we find that (3.2) implies a = −\frac{μ̄_{z}}{4} ; and if we set u = Re(μ), then (3.4) reads

u_{zz̄} + 2(e^{u} − |Q|^{2}e^{−u}) = 0.

Now gauging with e^{−iarg(μ)/2} gives μ̂ real, as required, so that ϕ_{1} is a harmonic lift. To

see that this holds for all λ, notice that multiplying ϕ_{λ} on the right by ( \begin{matrix}\sqrt{λ} & 0 \\ 0 & 1/\sqrt{λ} \end{matrix} )

is equivalent to replacing Q by λ^{2}Q in equation (2.10). Since the Gauss-Codazzi equations only involve |Q|^{2}, we see that by integrating the linear system we actually obtain a one-parameter family of CMC surfaces, all locally isometric, and with the same mean and Gaussian curvature. This is the classical result of Lie and Enneper mentioned above.

Remark 3.2 Notice that gauging by elements of K only changes the phase of c, so we can always take

ω = 2 ln |c|.  (3.5)

4 Loop Groups and Holomorphic Potentials

For any z ∈ D, the loop α_{λ}(z, z̄) defined in (3.1) is invariant under the involution

(σα)(λ) = ( \begin{matrix}1 & 0 \\ 0 & −1 \end{matrix} ) α(−λ) ( \begin{matrix}1 & 0 \\ 0 & −1 \end{matrix} ) ,  (4.1)

so we can regard α as a smooth map from D to the twisted loop algebra

Λ^{σ}su(2) = {f : S^{1} ⟶ su(2); σf = f}.  (4.2)

A loop f is twisted in the sense of this paper if, when expanded as a Fourier series in λ, the coefficients of odd powers of λ belong to P, while the even ones belong to K. If α is a loop of flat connections, then there’s a smooth map ϕ from D to the corresponding group of unitary loops

Λ^{σ}SU(2) = {ϕ : S^{1} ⟶ SU(2); σϕ = ϕ}  (4.3)

satisfying ϕ_{λ}^{−1}dϕ_{λ} = α_{λ}. The important part of the preceding proposition can therefore be restated as

Proposition 4.1 Let ϕ : D ⟶ Λ^{σ}SU(2) be differentiable and suppose ϕ_{λ}^{−1}dϕ_{λ} is a Laurent polynomial of degree 1. Then ϕ_{λ} is a harmonic lift, for all λ in S^{1}.

We shall call such a map ϕ_{λ} an extended harmonic lift. The construction in [DPW] gives a method for constructing these lifts from essentially free, holomorphic data:

Let ξ(z, λ) = ∑_{n=−1}^{∞} ξ_{n}(z)λ^{n} ∈ Λ^{σ}sl(2, C) be given arbitrarily with coefficients which are holomorphic in z (the set of such loops is denoted by Λ_{−1,∞}^{σ}). In this context, ξ is called a holomorphic potential. Treating λ as a parameter, solve the initial value problem

\frac{dg}{dz} = gξ, g(0, λ) = g_{0}(λ)  (4.4)

to obtain a map g : D ⟶ Λ^{σ}SL(2, C) which is holomorphic in z. For each fixed z ∈ D, the loop g(z, ·) can be factored as a product

g = ϕh,  (4.5)

where ϕ : D ⟶ Λ^{σ}SU(2), and h : D ⟶ Λ_{+}^{σ}SL(2, C). The subscript “+” means that the Fourier expansion

h(z, z̄, λ) = ∑_{n≥0} h_{n}(z, z̄)λ^{n}  (4.6)

contains no negative powers of λ. In particular, as a function of λ, h is non-singular and holomorphic for |λ| ≤ 1.

We shall construct this factorization for algebraic loops in the next section. Observe that (4.5) implies

ϕ_{λ}^{−1}∂ϕ_{λ} = hg^{−1}∂(gh^{−1}) = hξh^{−1} + h∂h^{−1}  (4.7)

is necessarily of the form ∑_{n≥−1} u_{n}λ^{n}. On the other hand,

ϕ_{λ}^{−1} ¯∂ϕ_{λ} = h_{λ} ¯∂h_{λ}^{−1} = ∑_{n≥0} v_{n}λ^{n}.  (4.8)

But ϕ is unitary, so \overline{ϕ_{λ}^{−1}∂ϕ_{λ}}^{t} = −ϕ_{λ}^{−1} ¯∂ϕ_{λ}. Since λ̄ = λ^{−1}, comparing the above expressions leads to the conclusion that ϕ_{λ}^{−1}dϕ_{λ} is a Laurent polynomial of degree one. Thus it satisfies the hypotheses of Proposition 4.1. So we obtain an extended harmonic lift from the holomorphic potential, and hence a 1-parameter family of CMC surfaces. The proof that each such lift can be obtained in this way may be found in [DPW].

Now suppose the holomorphic potential ξ has the leading term

ξ_{−1} = ( \begin{matrix}0 & F (z) \\ G(z) & 0 \end{matrix} ) .  (4.9)

Then letting h_{0} = ( \begin{matrix}a(z, z̄) & 0 \\ 0 & a^{−1}(z, z̄) \end{matrix} ) in the factorization (4.7) implies

U_{−1}(z, z̄) = h_{0}ξh_{0}^{−1} = ( \begin{matrix}0 & a^{2}F (z) \\ a^{−2}G(z) & 0 \end{matrix} ) .  (4.10)

Note that h|_{λ=0} ∈ SL(2, C) ⟹ h_{0} ∈ SL(2, C), ∀(z, z̄) ∈ D. Comparing this with the proof of Proposition 3.1 above, we see that the Hopf differential Q = −F G can be read off directly from ξ and that the solution to the second Gauss-Codazzi equation is given by

ω = 2 ln |\frac{G}{a^{2}} |.  (4.11)

The factorization (4.5) is of course not unique since we can replace ϕ_{λ} by ϕ_{λ}k, (and h_{λ} by k^{−1}h_{λ}) where k : D ⟶ K is smooth. This gives us just enough freedom in the splitting ϕh so the k : D ⟶ K can be chosen to make the (2, 1) entry of U_{−1} real and negative. This is one way of achieving uniqueness in the factorization. Summarizing, given the holomorphic potential ξ as above, the solution to g′ = gξ satisfying g(z_{0}, λ) = g_{0}(λ) has a unique (up to ϕ → −ϕ) factorization g = ϕh, where ϕ is a normalized, extended harmonic lift. Our goal in the remainder of the paper is to address the computational issues involved in (approximately) solving the differential equation in such a manner that the subsequent factorization can be performed accurately and efficiently.

5 Numerical methods I: Rapid Factorizations and Perfect Integrators

We begin with some facts about twisted algebraic loops; in practice, we’re interested in polynomials, of course.

Proposition 5.1 Let g = ∑_{n=−L}^{∞} g_{n}λ^{n} ∈ Λ^{σ}SL(2, C). Then there exists ϕ_{L} ∈ Λ^{σ}SU(2), a Laurent polynomial of degree 1, such that

ϕ_{L}^{∗}g = ∑_{n=−L+1}^{∞} ĝ_{n}λ^{n}.

ϕ_{L} is unique up to right multiplication by an element of K.

Proof: There are four essentially identical cases, of which we’ll do one. Suppose, then, that L is odd, so g_{−L} is off-diagonal. Thus

g(λ) = ( \begin{matrix}aλ^{−L+1} + O(λ^{−L+3}) & bλ^{−L} + O(λ^{−L+2}) \\ cλ^{−L} + O(λ^{−L+2}) & dλ^{−L+1} + O(λ^{−L+3}) \end{matrix} ) .

Then det g(λ) ≡ 1 ⟹ bc = 0, so either b or c (but not both) must vanish. We suppose c = 0, so that the (2, 1) entry of g has the form eλ^{−L+2} + O(λ^{−L+4}).

In order that ϕ_{L} have degree 1, ϕ_{L}^{∗} must either be of the form

ϕ_{L}^{∗} = ( \begin{matrix}α & γλ^{−1} \\ −γ̄λ & ᾱ \end{matrix} )  (5.1)

or of the same form with λ replaced by λ^{−1}. Computing with the matrix in (5.1), one finds that ϕ_{L}^{∗}g has negative degree −L + 1 ⟺

α = −\frac{γd}{b} and
|γ|^{2} = \frac{|b|^{2}}{|b|^{2} + |d|^{2}}

Choosing γ = \sqrt{|γ|^{2}} real and positive gives existence, and multiplying the resulting ϕ_{L}^{∗} on the left by an element of K gives the rest. The other three cases (c ≠ 0 and the corresponding two with L even) are entirely similar.

Proceeding inductively, we obtain:

Theorem 5.2. An algebraic element g ∈ Λ^{σ}SL(2, C) can be factored in the form

g = ϕh,  (5.2)

where ϕ = ∏_{i=1}^{L} ϕ_{i}(λ) is the product of Laurent polynomials ϕ_{i} ∈ Λ^{σ}SU(2) of degree 1 and h ∈ Λ_{+}^{σ}SL(2, C). The factorization is unique up to ϕh ⟶ (ϕk) · (k^{−1}h) = ϕ̃h̃, where k is a λ-independent element of K.

Recall that our goal is to solve g′ = gξ; since the equation is linear, we seek numerical methods of the general form

g_{n+1} = g_{n}A_{n}.  (5.3)

Suppose, in addition, that for some such scheme, A_{n} is a polynomial of degree d in λ^{−1} belonging to Λ^{σ}SL(2, C). Assume that g_{n} has been factored already as in the theorem:

g_{n+1} = ϕ_{n}h_{n}A_{n}.

Then h_{n}A_{n} has degree at most d in λ^{−1} so at most d operations of the type described in Proposition 5.1 are required to yield

h_{n}A_{n} = ϕ̂h_{n+1}

The “operations” of course take an amount of time that depends on n, since h_{n} has, in general, d more terms than h_{n−1}. Nevertheless, one easily checks that the time required to produce the entire list \{ϕ_{0}, ϕ_{1}, ⋯, ϕ_{n}\} is \underline{polynomial} in n. From a practical, computational viewpoint, this is a critical observation.

Definition 5.3. A consistent numerical method of the form g_{n+1} = g_{n}A_{n}, with A_{n} ∈ Λ^{σ}SL(2, C) a polynomial element of degree d (independent of n) in λ^{−1} will be called “perfect of degree d”. A consistent method with A_{n} a polynomial in Λ^{σ}GL(2, C) will be called “admissible”.

Remark 5.4 If the method is admissible, then by using a sufficiently small stepsize we can certainly arrange for det(A_{n}) = 1 ± ϵ, where ϵ is the machine precision. But this is by no means sufficient to force the conclusion that the matrix g_{−L} in (5.1) has precisely one nonzero entry. Since this is the key to rapid factorization, it is not true that admissible methods are effectively perfect. Indeed, our first computations for CMC surfaces done in the spirit of the present work [LS] used a factorization algorithm that required exponential time because the numerical method employed was not perfect in the above sense. The first perfect method to be used in this context was discovered by U. Pinkall. Once one adopts the “Hamiltonian” viewpoint of the next section, it will become clear that many such methods exist.

6 Numerical methods II: The Hamiltonian Metaphor

Let A(z, λ) = ξ^{t}, where ξ ∈ Λ_{−1,∞}^{σ} is a holomorphic potential and consider the equation

X′ = AX  (6.1)

for a two component vector X = ( \begin{matrix}p \\ q \end{matrix} ). Let us write A = ( \begin{matrix}a & c \\ b & −a \end{matrix} ) where the

entries are holomorphic in D, a contains only even powers of λ and b, c only odd powers. We set

H = ( \begin{matrix}b & −a \\ −a & −c \end{matrix} )  (6.2)
H = \frac{1}{2}X^{t}HX,  (6.3)

and then (6.1) assumes the Hamiltonian form

p′ = −\frac{∂H}{∂q}
q′ = \frac{∂H}{∂p} .  (6.4)

Of course, (6.4) is only formally Hamiltonian, since ′ = \frac{d}{dz} and not \frac{d}{dt} for some real variable t, not to mention the spectral parameter. Even so, we can make substantial use of both the formalism of classical mechanics [Lan] and the recent numerical work on symplectic integrators [S],[SC].

Recall that a one-step numerical scheme for integrating (6.4) will have the general form

p_{n+1} = α(p_{n}, q_{n}, z_{n})
q_{n+1} = β(p_{n}, q_{n}, z_{n}).  (6.5)

The scheme is symplectic if (p_{n}, q_{n}) ⟶ (p_{n+1}, q_{n+1}) (computed at time z_{n}) is a symplectic map for all n: dp_{n+1} ∧ dq_{n+1} = dp_{n} ∧ dq_{n}.

Symplectic integrators have received a lot of attention during the past few years, primarily because the numerical solutions often exhibit the correct qualitative behavior for relatively long times and relatively large time steps [SC] But this is not always the case as Herbst and Ablowitz have shown [HA].

For long time integrations, the most widely used methods have tended to be implicit (i.e., at each step, one must solve algebraically for (p_{n+1}, q_{n+1}) in terms of the other variables). But these methods are generally unsuited to our purposes because they do not preserve any of the desired λ behavior in the solutions.

Since our system is linear, we seek explicit methods of the form X_{n+1} = A_{n}X_{n}, and these are symplectic precisely when det(A_{n}) = 1. In particular, perfect integrators are necessarily symplectic. The converse is false, unfortunately, as there is no a priori reason that A_{n} should be a polynomial in λ or even that it should belong to the twisted loop group. However, of the many known symplectic schemes, some are perfect; and it’s not too difficult to invent new ones. The key idea for our purposes is simply to retain the condition det(A_{n}) = 1.

When the Hamiltonian can be written as H(p, q, z) = T (p, z) + V (q, z) it is said to be separable. In this case, Hamilton’s equations (6.4) take the simpler form

p′ = r(q, z)
q′ = s(p, z)  (6.6)

and there are many known explicit methods (most of which, however, only apply to the autonomous case ∂H/∂z = 0). Perhaps the simplest is the symplectic Euler method:

p_{n+1} = p_{n} + δr(q_{n}, z_{n})
q_{n+1} = q_{n} + δs(p_{n+1}, z_{n}).  (6.7)

Here δ is the step size. For our system (6.1), separable means a = 0, and the equation (6.7) reads

X_{n+1} = ( \begin{matrix}1 & δb_{n} \\ δc_{n} & 1 + δ^{2}b_{n}c_{n} \end{matrix} ) X_{n}  (6.8)

where b_{n} = b(z_{n}), etc.

In writing down an expression such as this, it is tacitly assumed that some grid has been chosen in the domain D, that z_{n+1} and z_{n} are adjacent grid points, and that δ = z_{n+1} − z_{n}. Note in particular that δ will depend on n, and that different “paths” through the grid from z_{ℓ} to z_{m} will generally yield different results. Since the functions b and c are at worst linear in λ^{−1}, and a contains no λ^{−1} terms, the symplectic Euler method is perfect of degree 2. The symplectic Euler method was first used in this context by U. Pinkall. While it has the virtue of being relatively fast, it is of rather limited use due to its first-order accuracy.

There are consistent, perfect, first order methods in the non-separable case as well. Briefly, the idea is as follows:

1. Write

( \begin{matrix}a & c \\ b & −a \end{matrix} ) = ( \begin{matrix}a & 0 \\ 0 & −a \end{matrix} ) + ( \begin{matrix}0 & c \\ b & 0 \end{matrix} ) =: A + B.

2. B is “separated” and we can choose any of the explicit methods alluded to above. Pick one such and write

X^{∗} = F (δ, B)X,

where F is symplectic (det = 1) and δ represents the step size Δz.

3. Let Y^{∗} = P^{−1}X^{∗}, where

P^{−1} = \frac{1}{\sqrt{2}} ( \begin{matrix}1 & −1 \\ 1 & 1 \end{matrix} ) .

Note that P^{−1}AP corresponds to a separated system of the form Ẏ = ÂY, and that P is orthogonal, hence symplectic.

4. Put Ŷ := F (δ, Â)Y and then X̂ = P Ŷ.

So all in all, we have the scheme

X_{n+1} = P F (δ, Â)P^{−1}F (δ, B)X_{n}.  (6.9)

Notice that it is not necessary to use the same “F ” in steps 2 and 4 above; we have just written down one of the simpler possibilities. Since

F (δ, Â) = I + δÂ + O(δ^{2}),

we have

P F (δ, Â)P^{−1} = I + δA + O(δ^{2}),

and therefore

P F (δ, Â)P^{−1}F (δ, B) = I + δ(A + B) + O(δ^{2}),

so the method is consistent. But there is a problem with the spectral parameter, as we discover if we try to implement this for the symplectic Euler method. Following the algorithm outlined above, we obtain

P F (δ, Â)P^{−1} = ( \begin{matrix}1 + δa + \frac{δ^{2}}{2} a^{2} & \frac{δ^{2}}{2} a^{2} \\ \frac{δ^{2}}{2} a^{2} & 1 − δa + \frac{δ^{2}}{2} a^{2} \end{matrix} ) .  (6.10)

The off-diagonal terms are of even degree in λ, so this matrix does not belong to the twisted loop group. (There is no problem with the matrix F (δ, B) of course, as we’ve seen above.) To cure this defect (in an admittedly ad hoc fashion), we simply replace the off-diagonal entries in this matrix by \frac{δ^{2}}{2} a^{2}λ and \frac{δ^{2}}{2} a^{2}λ^{−1} respectively. The resulting numerical method is still consistent.

For the sake of completeness, we can write out the whole expression:

X̂ = ( \begin{matrix}1 + δa + \frac{δ^{2}}{2} a^{2} + \frac{δ^{3}}{2} a^{2}cλ & δb + δ^{2}(\frac{a^{2}}{2} λ + ab) + \frac{δ^{3}}{2} a^{2}b + \frac{δ^{4}}{2} a^{2}bcλ \\ δc + \frac{δ^{2}}{2} (a^{2}λ^{−1} − 2ac) + \frac{δ^{3}}{2} a^{2}c & 1 − δa + \frac{δ^{2}}{2} (a^{2} + 2bc) \\ & +\frac{δ^{3}}{2} (a^{2}bλ^{−1} − 2abc) + \frac{δ^{4}}{2} a^{2}bc \end{matrix} ) X

(6.11)

This method is first order accurate, and perfect of degree 2, but it suffers from the same problems as Euler’s method. For reasonable accuracy and speed, we require higher order schemes.

In the very special case where H is both separable and autonomous, it is relatively easy to construct explicit higher order methods. Unfortunately, this is much too restrictive for our purposes. A holomorphic potential which is a polynomial in λ and independent of z yields a CMC surface of finite type [PiSt], provided that it satisfies the reality condition ξ̄_{i}^{t} = −ξ_{L−1−i} for −1 ≤ i ≤ \frac{L−1}{2} [BP]. The characteristic polynomial of the potential is the spectral curve of the resulting finite type surface. From this it is easy to check that most of the interesting known surfaces (e.g. Delaunay surfaces, Wente Tori, Bubbletons, Dobriner surfaces, etc.) have non-separable Hamiltonians.

Away from a discrete subset of D, every CMC surface can be obtained from a single separable potential of the form

ξ = λ^{−1} ( \begin{matrix}0 & F (z) \\ G(z) & 0 \end{matrix} ) .

For this and other reasons the separable case has received more attention. The reduction to this form depends on solving a Riemann-Hilbert problem in D [D], and is not explicit. But a local reduction to separable is straightforward:

Proposition 6.1 For ξ ∈ Λ_{−1,∞}^{σ} there is a guage transformation

ξ̂ = P^{−1}ξP + P^{−1}dP/dz

with P : D ⟶ Λ_{+}^{σ}SL(2, C) holomorphic in z, such that ξ̂ is off-diagonal.

Proof: Write ξ = ξ_{even} +ξ_{odd}, where “even”(odd) refers to all the terms of even (odd) degree in λ. Then

ξ̂ = P^{−1}ξ_{even}P + P^{−1}ξ_{odd}P + P^{−1}dP/dz.

If P is taken to have only even powers of λ, then P^{−1}ξ_{odd}P has only odd powers of λ. Assuming this, we want P such that

P^{−1}ξ_{even}P + P^{−1}dP/dz = 0

or

dP/dz + ξ_{even}P = 0

Since ξ_{even}, P and dP/dz are all diagonal, this is an uncoupled pair of scalar equations with the solution

P = Ae^{− ∫ ξ_{even}}.
♢

Note that this algorithm introduces z-dependence into the potential if it was not there already and typically gauges polynomial potential to non-polynomial ones.

In general, the system (6.1) is non-autonomous, and higher order explicit methods are more difficult to come by, even when H is separable. In the following section, we prove the existence of a perfect, second order method in a particular (very useful) case.

7 Numerical methods III: Non-autonomous Systems

Assuming that ξ is a Laurent polynomial, we recall that except for the finite type potentials, which are well understood, our Hamiltonian is formally nonautonomous, since a, b and c depend on z. Using a standard trick, we can remedy this by increasing the number of dependent variables by two: The new pair is p_{2} = s, q_{2} = z, and the new Hamiltonian is H̃ = s + H. Hamilton’s equations for the augmented system now read

˙p⃗ = \{ \begin{matrix}ṗ = ap + cq \\ ṡ = −\frac{∂H̃}{∂z} = −\frac{p^{2}}{2} \frac{db}{dz} + pq \frac{da}{dz} + \frac{q^{2}}{2} \frac{dc}{dz}, \end{matrix}
˙q⃗ = \{ \begin{matrix}q̇ = bp − aq \\ ż = 1 \end{matrix} .  (7.1)

The system is now autonomous with, formally, · = \frac{d}{dt}, but it’s non-linear and generally not separable. Nevertheless, it is possible, under some circumstances, to use an explicit method on (7.1) and have the projection on to the original 2-dimensional space be linear and even perfect. In fact, this rather far-fetched idea turns out to work in one of the most important cases which we now describe:

For an arbitrary holomorphic function f, let c = −λ^{−1}f(z), b = λ^{−1}, a = 0. (Recall that the zeros of the product bc determine the umbilic points on the surface, so by choosing f correctly, one can assign umbilics of any order at arbitrary locations.) Now

H̃ = (s + \frac{p^{2}λ^{−1}}{2} ) + (f(z)q^{2} \frac{λ^{−1}}{2} ).  (7.2)

is separable, and Hamilton’s equations are

˙p⃗ = G⃗(q⃗) = \{ \begin{matrix}ṗ = λ^{−1}f(z)q \\ ṡ = −\frac{λ^{−1}}{2} \frac{df}{dz}q^{2} \end{matrix}
˙q⃗ = H⃗(p⃗) = \{ \begin{matrix}q̇ = λ^{−1}p \\ ż = 1 \end{matrix} .  (7.3)

We now have

Theorem 7.1. In the above setup, the numerical scheme

q^{∗} = q⃗_{n} + \frac{δ}{2}H⃗(p⃗_{n})
p⃗_{n+1} = p⃗_{n} + δG⃗(q^{∗})
q⃗_{n+1} = q^{∗} + \frac{δ}{2}H⃗(p⃗_{n+1})

(known as the “leapfrog” scheme [S]) is perfect and second order for the original p and q.

Proof: A routine computation gives

p⃗_{n+1} = ( \begin{matrix}p_{n} + δλ^{−1}f(z + \frac{δ}{2})(q_{n} + \frac{δλ^{−1}}{2} p_{n}) \\ s_{n} − \frac{δλ^{−1}}{2} f′(z + \frac{δ}{2})(q_{n} + \frac{δλ^{−1}}{2} p_{n})^{2} \end{matrix} )
q⃗_{n+1} = ( \begin{matrix}q_{n} + \frac{δλ^{−1}}{2} (2p_{n} + δλ^{−1}f(z + \frac{δ}{2})(q_{n} + \frac{δλ^{−1}}{2} p_{n})) \\ z_{n} + δ \end{matrix} )

and in particular X_{n+1} = A_{n}X_{n}, with

A_{n} = ( \begin{matrix}1 + \frac{δ^{2}λ^{−2}}{2} f(z_{n} + \frac{δ}{2}) & λ^{−1}δf(z_{n} + \frac{δ}{2}) \\ δλ^{−1} + \frac{δ^{3}λ^{−3}}{4} f(z_{n} + \frac{δ}{2}) & 1 + \frac{δ^{2}λ^{−2}}{2} f(z_{n} + \frac{δ}{2}) \end{matrix} ) .

Evidently, A_{n} ∈ Λ^{σ}SL(2, C), is polynomial and, writing f(z_{n} + \frac{δ}{2}) = f(z_{n})+ \frac{δ}{2}f′(z_{n})+ O(δ^{2}), we find that

A_{n} = I + δ ( \begin{matrix}0 & λ^{−1}f(z_{n}) \\ λ^{−1} & 0 \end{matrix} ) +

\frac{δ^{2}}{2} ( \begin{matrix}λ^{−2}f(z_{n}) & λ^{−1}f′(z_{n}) \\ 0 & λ^{−2}f(z_{n}) \end{matrix} ) + O(δ^{3}).

So X_{n+1} = (I + δA + \frac{δ^{2}}{2!} (A^{2} + A′) + O(δ^{3}))X_{n} and the method is second order, as claimed. ♢

Higher order schemes belonging to the general class of partitioned Runge-Kutta methods can be found for this case as well. As far as we know, this is the first perfect method with better accuracy than the symplectic Euler method. It seems quite likely that by combining the techniques outlined in the preceding section with these, that 4th and higher order perfect methods can be obtained. This is currently under investigation.

References

[BP] F. E. Burstall and F. Pedit Harmonic maps via Adler-Kostant-Symes theory, in Harmonic Maps and Integrable Systems, Vieweg (1994), 221–272.

[D] J. Dorfmeister (private communication).

[DPW] J. Dorfmeister, F. Pedit and H. Wu, Weierstrass type representations of harmonic maps into symmetric spaces, U. Mass. Amherst GANG Preprint 1994.

[E] L. P. Eisenhart, A treatise on the differential geometry of curves and surfaces, Dover Publications, Inc., New York 1909.

[HA] B. M. Herbst and M. J. Ablowitz, Numerical homoclinic instabilities in the sine-Gordon equation, Quaest. Math. 15 (1992), 345–363.

[K] K. Kenmotsu, Weierstrass Formula for Surfaces of Prescribed Mean Curvature, Math. Ann. 245 (1979), 89–99.

[Lan] C. Lanczos, Variational principles of classical mechanics, Dover, New York (1986).

[LS] D. Lerner, and I. Sterling, Construction of constant mean curvature surfaces using the DPW representation of harmonic maps, (to appear in the Proceedings of the Symposium on Differential Geometry, Hamiltonian Systems, and Operator Theory, Univ. of the West Indies, Mona, Jamaica, Feb. 7–11 (1994)).

[PiSt] U. Pinkall and I. Sterling, Classification of constant mean curvature tori, Ann. Mat. 130 (1989), 407-451.

[Po] K. Pohlmeyer, Integrable Hamiltonian systems and interaction through quadratic constraints, Comm. Math. Phys. 46 (1976), 207–221.

[PrSe] A. Pressley and G. Segal, Loop Groups, Oxford (1988).

[RV] E. Ruh, and J. Vilms, The tension field of the Gauss map, Trans. Amer. Math. Soc. 149 (1970), 569–573.

[S] J. M. Sanz-Serna, Symplectic integrators for Hamiltonian problems: an overview, Acta Numerica (1991), 243–286.

[SC] J. M Sanz-Serna and M. Calvo, Numerical Hamiltonian Problems, Chapman and Hall, London (1994).

David Lerner Ivan Sterling Dept. of Mathematics Dept. of Mathematics University of Kansas University of Toledo Lawrence, Kansas 66045-2142 Toledo, Ohio 43606-3390