# Numerical Analysis for the Plateau Problem by the Method of Fundamental
Solutions

###### Abstract

Towards identifying the number of minimal surfaces sharing the same boundary
from the geometry of the boundary, we propose a numerical scheme with high
speed and high accuracy. Our numerical scheme is based on the method of
fundamental solutions. We establish the convergence analysis for Dirichlet
energy and $L^{\infty}$-error analysis for mean curvature. Each of the
approximate solutions in our scheme is a smooth surface, which is a
significant difference from previous studies that required mesh division.

##  1 Introduction

In 1762, Lagrange proposed the problem of finding a surface with the least
area spanned by a given closed Jordan curve [19]. Because the first variation
of the area of a surface gives the mean curvature vector, the mean curvature
of the surface with the least area is everywhere zero. Thus, a surface whose
mean curvature is everywhere zero is called a minimal surface. In a more
general sense, Lagrange’s problem can be rephrased as the problem of finding a
minimal surface spanned by a given closed Jordan curve. Later in 1873, Plateau
investigated the properties of minimal surfaces through an experiment on soap
films and pointed out that a single closed wire, regardless of its geometry,
bounded at least one soap film. Nowadays, the problem of the existence of
minimal surfaces bounded by a given closed Jordan curve is called Plateau’s
problem.

Although Douglas has established the existence of a solution to the Plateau
problem [9], it is difficult to determine whether the solution is unique or
whether there are finitely many solutions, even if the simple closed curve is
smooth and if the solutions are homeomorphic to the disk. The problem has not
been completely solved. For example, the solution is unique if a closed Jordan
curve has a one-to-one parallel projection onto a convex closed curve in the
plane [24] or if the total curvature is less than or equal to $4\pi$ [20]. For
the case where there is more than one solution, it is known that there are two
other minimal surfaces with the same boundary as the Enneper surface, which is
one of the exact solutions of the Plateau problem [22]. The finite solvability
is established for generic curves [5], stable solutions [18], and a curve
whose total curvature is less than $6\pi$ [21]. As seen from these results,
the problem of determining the number of solutions according to the geometry
of a closed Jordan curve is one of the most important issues in the
geometrical analysis of minimal surfaces.

To solve this problem, numerical analysis on minimal surfaces should help to
find a connection between the geometry of the curve and the number of minimal
surfaces sharing the curve as the boundary. Since minimal surfaces are
considered stationary points of functionals, such as area and Dirichlet
integrals, we can obtain minimal surfaces as the convergence limit of some
optimization problem. Therefore, one possible method to determine the number
of solutions corresponding to a given simple closed curve is to take a random
initial guess for the optimization problem and count the number of different
minimal surfaces obtained as the convergence limit of the optimization
problem.

In order to make the above approach practical, it is necessary to establish a
numerical scheme, which is

  1. 1.

so fast that many initial values can be handled in a reasonable time, and

  2. 2.

so accurate that the numerical solutions obtained can be distinguished.

Therefore, this paper aims to establish a numerical scheme with high speed and
high accuracy enough to achieve the above goals.

Numerical analysis of minimal surfaces has a long history, and since Douglas
[8] proposed a numerical scheme using the finite difference method, various
methods have been proposed up to the present day. Tsuchiya gave the first
theoretical convergence analysis [29, 30]. He considered the finite element
method and proved the existence of discrete minimal surfaces and convergence
for the $H^{1}$ norm. However, because he used non-direct arguments, he did
not derive the convergence order. The finite element method is used for the
analysis in the Dziuk–Huthinson paper [10, 11], as in Tsuchiya’s papers.
Although the target is somewhat limited to nondegenerate minimal surfaces, the
error for the $H^{1}$ norm is proved to be $\mathrm{O}(h)$. Their result is
the first analytical one that includes the order of convergence. Following
this work, Pozzi [23] gave the $L^{2}$ error estimate for finite element
solutions. In Dziuk–Hutchinson [12], a method for computing surfaces with a
specified mean curvature is given with error analysis. Numerical analysis of
minimal surfaces whose boundaries are given as polygons is considered by Hinze
[15], who proves convergence for the $H^{1}$ norm. Other than these, various
numerical methods have been proposed, for example, high-precision spatial
discretization using B-spline curves [14] and higher-order polynomials [28],
computation of mean curvature flow by finite volume method [27], improved
performance by mesh refinement [13], and methods based on algebraic topology
[26] and differential forms [31]. However, to our knowledge, no convergence
analysis has been given for any of them.

As seen in Definition 2.1, numerical computation of minimal surfaces requires
solving the Dirichlet boundary value problem for the Laplace equation on the
unit disk. Therefore, a fast and accurate solver for the Laplace equation can
be used to compute minimal surfaces more accurately than in previous studies.
Thus, in this paper, we employ the method of fundamental solutions (MFS for
short), a numerical method for the potential problem. The MFS is a mesh-free
numerical solver because it does not require meshing of the region, as with
the finite element and finite difference methods. It is so named because a
linear combination of the fundamental solution of the partial differential
operator of interest constructs the approximate solution. Under certain
conditions, the approximate solution by the MFS has the remarkable property of
exponentially converging to the exact solution with an increasing number of
approximation points (Theorem 2.6 holds in the setting of this paper).
Moreover, since a linear combination of fundamental solutions constructs the
approximate solution, the approximate solution is smooth, and differential
operations can be performed analytically, resulting in a smooth approximation
of minimal surfaces. When the Jordan curve is embedded into a two-dimensional
plane, the problem of finding a minimal surface is reduced to the problem of
finding a conformal map from the unit disk onto the simply-connected region
bounded by the Jordan curve. Amano and his colleagues have been known to
compute conformal mappings with high accuracy using the MFS [1, 2, 3, 4, 25].
However, to our knowledge, the MFS has no successful application to the
numerical computation of minimal surfaces.

By utilizing the convergence theorem for the Dirichlet boundary value problem
of the Laplace equation in the disk region (Theorem 2.6), we succeed in giving
the existence of approximate minimal surfaces (Theorem 4.1), and the
convergence of Dirichlet energy and evaluation of mean curvature of
approximate minimal surfaces (Theorem 4.2). These are the first results of the
MFS numerical analysis of minimal surfaces. As a result, the MFS can now
obtain smooth approximations of minimal surfaces, a significant difference
from previous studies that required mesh division.

This paper is organized as follows. In Section 2, we first briefly summarize
the geometry of surfaces, then formulate the Plateau problem and present known
results on the existence of minimal surfaces. Then, the MFS is explained in
detail, and the results of the convergence analysis used in this paper are
presented. In Section 3, we build the basis of an algorithm for solving the
Plateau problem and, in particular, present the existence theorem for
approximate minimal surfaces. Section 4 shows that the approximate minimal
surfaces obtained by our algorithm converge to the minimal surfaces. In
particular, we show the Dirichlet energy’s convergence and the mean
curvature’s $L^{\infty}$-error. In Section 5, after describing the algorithm
in actual numerical computations, we demonstrate the usefulness of the
algorithm proposed in this paper by giving various numerical examples. In
Section 6, we propose an algorithm for computing all minimal surfaces based on
the numerical scheme proposed in this paper and verify its usefulness through
numerical experiments. Finally, Section 7 summarizes the paper and provides
directions for future research.

##  2 Preliminaries

###  2.1 Geometry of surfaces

In this section, we briefly describe some basic notions from differential
geometry closely related to this work. Let $B_{\rho}$ be the disk with the
radius $\rho$ in the complex plane $z=x^{1}+\mathrm{i}x^{2}\in{\mathbb{C}}$
and $\partial B$ be its boundary:

$\displaystyle B_{\rho}:$ $\displaystyle=\\{z\in{\mathbb{C}}|\,|z|<\rho\\},$
(1) $\displaystyle\partial B_{\rho}:$
$\displaystyle=\\{z\in{\mathbb{C}}|\,|z|=\rho\\}.$ (2)

For notational simplicity, we write $B=B_{1}$. Let us denote the
differentiation by $\partial=\partial/\partial z$,
$\bar{\partial}=\partial/\partial z$, $\partial_{1}=\partial/\partial x^{1}$,
and $\partial_{2}=\partial/\partial x^{2}$, respectively. For each $f\in
C^{1}({\mathbb{C}})$, we now have

$\displaystyle\partial$
$\displaystyle=(\partial_{1}-\mathrm{i}\partial_{2})/2,$ (3)
$\displaystyle\bar{\partial}$
$\displaystyle=(\partial_{1}+\mathrm{i}\partial_{2})/2,$ (4)
$\displaystyle\triangle$
$\displaystyle=4\partial\bar{\partial}=\partial_{u}^{2}+\partial_{v}^{2},$ (5)
$\displaystyle 4|\partial f|^{2}$
$\displaystyle=(\partial_{1}f)^{2}+(\partial_{2}f)^{2},$ (6) $\displaystyle
4(\partial f)^{2}$
$\displaystyle=(\partial_{1}f)^{2}-(\partial_{2}f)^{2}+2\mathrm{i}(\partial_{1}f)(\partial_{2}f),$
(7)

Given a vector-valued function $X:z\in
B\to(X_{i}(z))_{i=1}^{3}\in{\mathbb{R}}^{3}$, we write
$\partial_{j}X=(\partial_{j}X_{i})_{i=1}^{3}\in{\mathbb{R}}^{3}$ and
$\triangle X=(\triangle X_{i})_{i=1}^{3}\in{\mathbb{R}}^{3}$.

A surface $(M,g)$ is a Riemann surface equipped with a Riemannian metric $g$,
which yields complex analysis, and Riemannian geometry can be utilized as the
surface. We now focus on some geometric structures defined on the surface.
Remenber that a complex function $f\colon{\mathbb{C}}\to{\mathbb{C}}$ is said
to be holomorphic if it satisfies the Cauchy-Riemann equation

$\displaystyle\bar{\partial}f=0.$ (8)

Moreover, if $f$ is bijective and the inverse is also holomorphic, then $f$ is
called a biholomorphism.

The Riemann surface $M$ is equipped with complex charts
${(U_{\alpha},\phi_{\alpha})}_{\alpha}$ such that for each
$(U_{\alpha},\phi_{\alpha})$ and $(U_{\beta},\phi_{\beta})$,
$\phi_{\alpha}\circ\phi_{\beta}^{-1}:\phi_{\beta}(U_{\alpha}\cap
U_{\beta})\to\phi_{\alpha}(U_{\alpha}\cap U_{\beta})\subset{\mathbb{C}}$ is a
biholomorphism, and $M=\bigcup_{\alpha}U_{\alpha}$ holds. The union of all
compatible complex charts defines the complex structure on the surface, and a
smooth $2$-manifold with a complex structure is called a Riemann surface.
Given two Riemann surfaces $M$ and $N$, a diffeomorphism $f:M\to N$ is a
biholomorphism if $f$ is a biholomorphism in each complex chart. In addition,
if a biholomorphism exists between two Riemann surfaces $M$ and $N$, then the
Riemann surface $M$ is said to be biholomorphic to $N$.

The Riemannian metric $g$ induces a conformal structure. A chart $(U,\phi)$ is
isothermal if the metric has the following local representation in the chart:

$\displaystyle g=\lambda^{2}|\operatorname{d}z|^{2},$ (9)

where $\lambda>0$ and $\lambda$ is called the conformal factor. A geometric
structure defined by the union of all compatible isothermal charts is called a
conformal structure, and a smooth $2$-manifold equipped with a conformal
structure is called a conformal surface. Given two conformal surfaces $M$ and
$N$, a diffeomorphism $f:M\to N$ is a conformal mapping if $f$ is conformal in
each isothermal chart. Moreover, if a conformal mapping exists between two
Riemann surfaces $M$ and $N$, then the Riemann surface $M$ is said to be
conformal to $N$. It is easy to see that every isothermal chart becomes
complex, and the converse holds using conventional identification
$z=x^{1}+\mathrm{i}x^{2}$. Hence, given surfaces $M$ and $N$, a diffeomorphism
$f:M\to N$ is a conformal mapping if and only if $f$ is a biholomorphism.

For instance, let a vector-valued function $X\colon B\to{\mathbb{R}}^{3}$ be a
$C^{1}$ immersion; that is, the the rank of the matrix
$(\partial_{i}X_{j})_{i,j}$ is $2$. Then, disk $B$ may have two different
Riemannian metrics: first, the standard 2D Euclidean metric
$g_{s}=(\operatorname{d}x^{1})^{2}+(\operatorname{d}x^{2})^{2}$, and second,
the induced metric
$g_{X}=\sum_{i,j}\langle\partial_{i}X,\partial_{j}X\rangle\operatorname{d}x^{i}\operatorname{d}x^{j}$,
where $\langle\cdotp,\cdotp\rangle$ is the 3D Euclidean metric. For the
identity map ${\rm id}:(B,g_{s})\to(B,g_{X})$ to be conformal, the immersion
$X$ must satisfy

$\displaystyle\langle\partial_{1}X,\partial_{1}X\rangle=\langle\partial_{2}X,\partial_{2}X\rangle,\quad\langle\partial_{1}X,\partial_{2}X\rangle=0.$
(10)

By using the complex coordinate, it is rewritten as

$\displaystyle\Phi_{X}(z):$ $\displaystyle=\sum_{i=1}^{3}(\partial
X_{i}(z))^{2}$ (11)
$\displaystyle=\frac{1}{4}\left(\langle\partial_{1}X,\partial_{1}X\rangle_{z}-\langle\partial_{2}X,\partial_{2}X\rangle_{z}+2\mathrm{i}\langle\partial_{1}X,\partial_{2}X\rangle_{z}\right)$
(12) $\displaystyle=0.$ (13)

The function $\Phi_{X}:B\to{\mathbb{C}}$ is called the complex dilatation. If
and only if the complex dilatation $\Phi_{X}$ vanishes, the induced metric
$g_{X}$ satisfies $g_{X}=\lambda^{2}g_{s}$, where
$\lambda^{2}=\langle\partial_{1}X,\partial_{1}X\rangle=\langle\partial_{2}X,\partial_{2}X\rangle$.
Moreover, if $\triangle X=0$, $\Phi_{X}$ is a holomorphic function since

$\displaystyle\bar{\partial}\Phi_{X}$ $\displaystyle=\sum_{i=1}^{3}2\partial
X_{i}(z)\bar{\partial}\partial X_{i}(z)$ (14)
$\displaystyle=\frac{1}{2}\sum_{i=1}^{3}\partial X_{i}(z)\triangle X_{i}(z)$
(15) $\displaystyle=0.$ (16)

###  2.2 Plateau problem

We review the Plateau problem in terms of a construction of solutions. The
following facts in this subsection is briefly introduced by the Dierkes,
Hildebrandt and Sauvigny [7] in Chapter 4.2. The Plateau problem is formulated
as finding an immersed surface spanned by a given Jordan curve that minimizes
the area. A variational argument yields that the minimizer of the area has
zero mean curvature. The mean curvature $H_{X}$ of a parametrized surface
$X:B\to{\mathbb{R}}^{3}$ satisfies

$\displaystyle H_{X}n=\triangle_{X}X,$ (17)

where $n=\partial_{1}X\times\partial_{2}X/|\partial_{1}X\times\partial_{2}X|$
is the unit normal vector and $\triangle_{X}$ is the Laplace-Beltrami operator
associated with the metric $g_{X}$ if $|\partial_{1}X\times\partial_{2}X|\neq
0$. In particular, if the complex dilatation $\Phi_{X}$ vanishes, we can
deduce

$\displaystyle\triangle_{X}=\lambda^{-2}\triangle,$ (18)

where $\triangle$ is the Laplacian
$\triangle=\partial_{1}^{2}+\partial_{2}^{2}$. Hence, the minimal surface is
formulated as a solution to the following problem.

######  Definition 2.1.

Given a closed Jordan curve $\Gamma\subset{\mathbb{R}}^{3}$,
$X:\bar{B}\to{\mathbb{R}}^{3}$ is called a minimal surface spanned by $\Gamma$
if the following four conditions are satisfied:

  1. 1.

$X\in C^{0}(\bar{B},{\mathbb{R}}^{3})\cap C^{2}(B,{\mathbb{R}}^{3});$

  2. 2.

$\triangle X=0;$

  3. 3.

the restriction $X|_{\partial B}:\partial B\to\Gamma$ is a homeomorphism;

  4. 4.

$\Phi_{X}=0.$

Courant [6] established the existence of a minimal surface spanned by a given
closed Jordan curve by minimizing the Dirichlet energy

$\displaystyle
D(X):=\frac{1}{2}\int_{B}(|\partial_{1}X|^{2}+|\partial_{2}X|^{2})\operatorname{d}x^{1}\operatorname{d}x^{2}$
(19)

in the Sobolev class $X\in H^{1}(B,{\mathbb{R}}^{3})$. More precisely, the
minimizing problem of the Dirichlet integral is performed in the following
space of admissible functions.

######  Definition 2.2.

Given a closed Jordan curve $\Gamma$ in ${\mathbb{R}}^{3}$, a mapping $X\in
H_{2}^{1}(B,{\mathbb{R}}^{3})$ is said to be of class $\mathcal{C}(\Gamma)$
for a fixed orientation if its Sobolev trace $X|_{\partial B}$ can be
represented by a weakly monotonic, continuous mapping $\varphi:\partial
B\to\Gamma$ onto $\Gamma$.

Unless otherwise stated, hereafter, $\mathcal{C}(\Gamma)$ is always defined
for a fixed orientation. Consequently, the minimizing problem of the Dirichlet
integral is formulated as follows:

$\displaystyle\mathcal{P}(\Gamma):\quad D(X)\to\min\quad\text{in the class
}\mathcal{C}(\Gamma).$ (20)

The existence of a solution to the problem $\mathcal{P}(\Gamma)$ is obtained
when $\mathcal{C}(\Gamma)$ is nonempty. In particular, it is satisfied if
$\Gamma$ is a closed Jordan curve of finite length.

######  Theorem 2.3 ([7, Chapter 4.3, Theorem 1]).

If $\mathcal{C}(\Gamma)$ is nonempty, then the minimizing problem
$\mathcal{P}(\Gamma)$ has at least one solution, continuous on $\bar{B}$ and
harmonic in $B$. In particular, $\mathcal{P}(\Gamma)$ has such a solution for
every rectifiable curve $\Gamma$.

It follows from Weyl’s lemma that all minimizers $X\in{\cal C}(\Gamma)$ are
actually, i.e., $\triangle X=0$. However, changing the coordinate, the
minimizer is not in general harmonic again in the new coordinate since the
Laplacian $\triangle$ and the Dirichlet energy $D(X)$ change as the
coordinates are changed. Courant [6] showed that a minimal surface is obtained
by taking the variation of the Dirichlet energy by changing variables.

######  Theorem 2.4 ([7, Chapter 4.5, Theorem 2]).

Every solution $X$ of the variational problem $\mathcal{P}(\Gamma)$ is a
minimal surface.

In particular, the following is helpful to show the existence of the limit for
an approximate solution for a minimal surface.

######  Theorem 2.5 ([7, Chapter 4.3, Theorem 3]).

Let $\\{\Gamma_{n}\\}_{n\geq 1}$ be a sequence of closed (oriented) Jordan
curves in ${\mathbb{R}}^{3}$, which converge in the sense of Fréchet to some
closed (oriented) Jordan curve $\Gamma$. Let $\\{X_{n}\\}_{n\geq 1}\subset
C^{0}(\bar{B},{\mathbb{R}}^{3})\cap C^{2}(B,{\mathbb{R}}^{3})$ be $\triangle
X_{n}=0$ and $X_{n}(\partial B)=\Gamma_{n}$. Then, there exists a subsequence
$\\{X_{n_{p}}\\}_{p\geq 1}$ and $X\subset C^{0}(\bar{B},{\mathbb{R}}^{3})\cap
C^{2}(B,{\mathbb{R}}^{3})$ with $\triangle X=0$ and $X(\partial B)=\Gamma$
such that $X_{n_{p}}\to X$ uniformly on $\bar{B}$ as $p\to\infty$.

The above facts is briefly introduced by the Dierkes, Hildebrandt and Sauvigny
[7] in Chapter 4.

###  2.3 Method of fundamental solutions

The method of fundamental solutions (MFS for short) is a mesh-free numerical
solver for linear partial differential equations such as the Laplace equation,
the Helmholtz equation, and the biharmonic equation. Its idea is quite simple,
and the algorithm is described for the Laplace equation, which is the subject
of this paper.

Let $\Omega$ be a bounded region in $\mathbb{C}$ with smooth boundary
$\partial\Omega$, and consider the Dirichlet boundary value problem for the
Laplace equation in $\Omega$ with a given boundary data
$f:\partial\Omega\to{\mathbb{R}}$:

$\displaystyle\begin{dcases*}\triangle u=0&in $\Omega$,\\\ u=f&on
$\partial\Omega$.\end{dcases*}$ (21)

The MFS constructs an approximate solution for this problem by the following
procedure.

  1. 1.

Take $N\in\mathbb{N}$ and fix it. Moreover, arrange $N$ points $\zeta_{k}$
($k=1,2,\ldots,N$) “suitably” in $\mathbb{C}\setminus\overline{\Omega}$, which
we call the singular points.

  2. 2.

Seek an approximate solution $u^{(N)}$ in the following form:

$\displaystyle u^{(N)}(z)=\sum_{k=1}^{N}Q_{k}G(z-\zeta_{k}),$ (22)

where $G(z)=(2\pi)^{-1}\log|z|$ is the fundamental solution of the Laplace
operator. Note that $u^{(N)}$ satisfies the Laplace equation exactly in
$\Omega$ since the singular points $\\{\zeta_{k}\\}_{k=1}^{N}$ are outside
$\Omega$.

  3. 3.

Determine coefficients $\\{Q_{k}\\}_{k=1}^{N}$ by the collocation method.
Namely, choose $N$ points $z_{j}$ ($j=1,2,\ldots,N$) “suitably” on
$\partial\Omega$, and impose the following conditions:

$\displaystyle u^{(N)}(z_{j})=f(z_{j}),\quad j=1,2,\ldots,N.$ (23)

Eq. (23) can be rewritten as a linear system called the collocation equations,

$\displaystyle\mathbf{G}\bm{Q}=\bm{f},$ (24)

where

$\displaystyle\mathbf{G}=(G(z_{j}-\zeta_{k}))_{j,k}\in\mathbb{R}^{N\times
N},\quad\bm{Q}=(Q_{k})_{k}\in\mathbb{R}^{N},\quad\bm{f}=(f(z_{j}))_{j}\in\mathbb{R}^{N}.$
(25)

As seen from the algorithm, the MFS does not require meshing the region, and
the approximate solution is constructed by choosing appropriate points on the
boundary and outside the region. However, what constitutes appropriate point
placement is still mathematically unsolved. In this paper, it is only
necessary to consider the case where the problem region is the unit disk. In
this case, it is natural to place the collocation points
$\\{z_{j}\\}_{j=1}^{N}$ and singular points $\\{\zeta_{k}\\}_{k=1}^{N}$
uniformly on concentric circles as follows:

$\displaystyle z_{j}$ $\displaystyle=\omega^{j},$ $\displaystyle
j=1,2,\ldots,N,$ (26) $\displaystyle\zeta_{k}$ $\displaystyle=R\omega^{k},$
$\displaystyle k=1,2,\ldots,N,$ (27)

where $R>1$ and $\omega=\exp(2\pi\mathrm{i}/N)$.

Then, we can solve the collocation equations (24) explicitly. Since the
coefficient matrix $\bm{G}$ is now circulant, its inverse
$\bm{G}^{-1}=(G^{-1}_{kj})_{kj}$ is presented by

$\displaystyle
G_{kj}^{-1}\coloneqq\frac{1}{N}\sum_{l=1}^{N}\frac{\omega^{(k-j)(l-1)}}{\varphi_{l-1}^{(N)}},$
(28)

where

$\displaystyle\varphi_{p}^{(N)}\coloneqq\sum_{k=1}^{N}\omega^{p(k-1)}G(1-\zeta_{k}),\quad
p\in\mathbb{Z}.$ (29)

As a result, the coefficients $Q_{k}$ are explicitly given by

$\displaystyle
Q_{k}=\sum_{j=1}^{N}G_{kj}^{-1}f(z_{j})=\frac{1}{N}\sum_{j=1}^{N}\sum_{l=1}^{N}\frac{\omega^{(k-j)(l-1)}}{\varphi_{l-1}^{(N)}}f(z_{j})$
(30)

for each $k=1,2,\ldots,N$. Hence, we find that an approximate solution exists
and that it can be concretely constructed.

Under the above setting, the following theorem holds.

######  Theorem 2.6 ([17, Theorem 2], [16, Theorem 2.3]).

  1. 1.

Suppose that the boundary data $f$ is real analytic. Then, there are constants
$C>0$ and $\tau\in(0,1)$, independent of $N$, such that

$\displaystyle\|u-u^{(N)}\|_{L^{\infty}(B)}\leq C\tau^{N}.$ (31)

  2. 2.

Let $\\{f_{n}\\}$ be the Fourier coefficients of $f$.

    1. (a)

If the Fourier series $\sum_{n}f_{n}\mathrm{e}^{\mathrm{i}n\theta}$ is
absolutely convergent, then the approximate solution uniformly converges to
the exact solution in $B$ as $N\to\infty$; that is,

$\displaystyle\|u-u^{(N)}\|_{L^{\infty}(B)}\longrightarrow 0\quad(N\to\infty)$
(32)

    2. (b)

If $f_{n}=\mathrm{O}(|n|^{-\alpha})$ for some $\alpha>1$ as $|n|\to\infty$,
then we have

$\displaystyle\|u-u^{(N)}\|_{L^{\infty}(B)}=\mathrm{O}(N^{-\alpha+1})\quad(N\to\infty).$
(33)

Since the approximate solution (22) by the MFS is analytic, its derivatives
can also be computed analytically. For instance,

$\displaystyle\partial u^{(N)}(z)=\sum_{k=1}^{N}Q_{k}\partial G(z-\zeta_{k}).$
(34)

As pointed out in [17, Section 4] and [16, Remark 4.1], under the same
situation in Theorem 2.6, we can also prove similar estimates for $\|\partial
u-\partial u^{(N)}\|_{L^{\infty}(B)}$; the $W^{1,\infty}$-error
$\|u-u^{(N)}\|_{W^{1,\infty}(B)}$ tends to $0$ as $N\to\infty$ under mild
assumptions on the boundary data $f$. The maximum principle for harmonic
functions implies that

$\displaystyle\|u-u^{(N)}\|_{L^{p}(B)}\leq\pi^{1/p}\|u-u^{(N)}\|_{L^{\infty}(B)},$
(35) $\displaystyle\|\partial u-\partial
u^{(N)}\|_{L^{p}(B)}\leq\pi^{1/p}\|\partial u-\partial
u^{(N)}\|_{L^{\infty}(B)}$ (36)

for $p\in[1,\infty)$. Repeating the same procedure, we find that, for
$m\in\mathbb{N}$ and $p\in[1,\infty]$, the $W^{m,p}$-error
$\|u-u^{(N)}\|_{W^{m,p}(B)}$ tends to $0$ as $N\to\infty$, depending on the
regularity of the solution.

##  3 Numerical scheme solving Plateau problem

In what follows, we construct a minimal surface $X\in
C^{0}(\bar{B},{\mathbb{R}}^{3})\cap C^{2}(B,{\mathbb{R}}^{3})$ spanned by a
given rectifiable closed Jordan curve $\Gamma\subset{\mathbb{R}}^{3}$. Let us
fix a homeomorphism $b:\partial B\to\Gamma$. Then, since $X|_{\partial
B}:\partial B\to\Gamma$ is also homeomorphism, we deduce that there exists a
homeomorphism $\phi:\partial B\to\partial B$ such that

  1. 1.

$\triangle X=0$ on $B$;

  2. 2.

$X=b\circ\phi$ on $\partial B$;

  3. 3.

$\Phi_{X}=0$ on $B$,

by taking $\phi=b^{-1}\circ X|_{\partial B}$. Note that for a given $\phi$,
$X$ is solved using Poisson kernel $P$ of the Dirichlet boundary value
problem; that is, $X$ is given by

$\displaystyle
X=X(z;\phi)=(P*(b_{1}\circ\phi),P*(b_{2}\circ\phi),P*(b_{3}\circ\phi)).$ (37)

From this point of view, we can say $\phi$ is chosen to attain
$\Phi_{X(\cdotp;\phi)}=0$. In other words, for a given boundary mapping
$b:\partial B\to\Gamma$, find $\phi:\partial B\to\partial B$ subject to Eq.
(37) and $||\Phi_{X(\cdotp;\phi)}||_{L^{\infty}(B)}=0$. Remember that
$\Phi_{X}$ is a holomorphic function if $\triangle X=0$. Owing to the maximum
principle, for minimizing $||\Phi_{X(\cdotp;\phi)}||_{L^{\infty}(B)}$, it
suffices to minimize $||\Phi_{X(\cdotp;\phi)}||_{L^{\infty}(\partial
B_{\rho})}$ for sufficiently close $\rho\in(0,1]$ to $1$. Hence, we can obtain
the desired $\phi$ by solving the minimization problem
$\min\\{||\Phi_{X(\cdotp;\phi)}||_{L^{\infty}(\partial
B_{\rho})}|\,\phi\in\mathrm{Homeo}(\partial B)\\}$, where
$\mathrm{Homeo}(\partial B)$ is the space of all homeomorphisms on $\partial
B$. In this paper, we solve Eq. (37) and
$\min\\{||\Phi_{X(\cdotp;\phi)}||_{L^{\infty}(\partial
B_{\rho})}|\,\phi\in\mathrm{Homeo}(\partial B)\\}$ by the MFS and Nesterov’s
accelerated gradient descent, respectively.

First, discretize Eq. (37) for a given $b:\partial B\to\Gamma$ and a given
$\phi:\partial B\to\partial B$ by the following procedure.

  1. 1.

Take $z_{j}$ and $\zeta_{k}$ given by Eq. (26) and Eq (27) for a given
$N\in{\mathbb{Z}}_{\geq 1}$.

  2. 2.

Take a transformation vector
$\bm{\phi}=(\phi(z_{j}))_{j=1}^{N}\in\mathbb{T}^{N}$, where $\mathbb{T}^{N}$
is the $N$-dimensional torus.

  3. 3.

Seek an approximate solution
$X_{i}^{(N)}(\cdotp;\bm{\phi}):\bar{B}\to{\mathbb{R}}$ for each $i=1,2,3$ in
the following form:

$\displaystyle
X_{i}^{(N)}(z;\bm{\phi})=\sum_{k=1}^{N}Q_{ik}(\bm{\phi})G(z-\zeta_{k})\quad\text{on
}B.$ (38)

  4. 4.

Determine the coefficients
$\bm{Q}_{i}(\bm{\phi})=(Q_{ik}(\bm{\phi}))_{k=1}^{N}\in{\mathbb{R}}^{N}$ by

$\displaystyle X_{i}^{(N)}(z_{j};\bm{\phi})$ $\displaystyle=b_{i}(\phi_{j}),$
(39)

which can be solved explicitly using Eq. (30).

  5. 5.

Make an approximate surface
$X^{(N)}(\cdotp;\bm{\phi}):\bar{B}\to{\mathbb{R}}^{3}$ by

$\displaystyle
X^{(N)}(z;\bm{\phi})=(X_{1}^{(N)},X_{2}^{(N)},X_{3}^{(N)})(z;\bm{\phi}).$ (40)

Second, minimize $||\Phi_{X^{(N)}(\cdotp;\bm{\phi})}||_{L^{\infty}(B)}$ in the
class $\phi\in\mathrm{Homeo}(\partial B)$. Since $X^{(N)}$ is determined by
$\bm{\phi}\in\mathbb{T}^{N}$, the admissible space for the minimization
problem can be reduced to the finite-dimensional space $\mathbb{T}^{N}$ from
the infinite-dimensional space $\mathrm{Homeo}(\partial B)$. We write
$\Phi_{X^{(N)}(\cdotp;\bm{\phi})}(z)$ by $\Phi^{(N)}(z;\bm{\phi})$ shortly. It
is worth noting that $\triangle X^{(N)}(z,\bm{\phi})=0$ for any
$\bm{\phi}\in\mathbb{T}^{N}$, which yields that $\Phi^{(N)}(z;\bm{\phi})$
becomes a holomorphic function. Hence, it is sufficient to minimize
$||\Phi^{(N)}(\cdotp;\bm{\phi})||_{L^{\infty}(\partial B_{\rho})}$ with a
given sufficiently close $\rho\in(0,1]$ to $1$. In particular, we discretize
$||\Phi^{(N)}(\cdotp;\bm{\phi})||_{L^{\infty}(\partial B_{\rho})}$ by

$\displaystyle E=E(\bm{\phi})=\sum_{j=1}^{N}|\Phi^{(N)}(\rho
z_{j};\bm{\phi})|^{2}.$ (41)

Hence, we deduce to minimize $E(\bm{\phi})$ in the class
$\bm{\phi}\in\mathbb{T}^{N}$ by the following procedure.

  1. 1.

Choose an initial vector $\bm{\phi}_{1}\in\mathbb{T}^{N}$ arbitrarily and a
small enough step size $\eta\in(0,1)$.

  2. 2.

Set $\bm{\varphi}_{1}=\bm{\phi}_{1}$

  3. 3.

Update

$\displaystyle\bm{\varphi}_{n+1}$ $\displaystyle=\bm{\phi}_{n}-\eta\nabla
E(\bm{\phi}_{n}),$ (42) $\displaystyle\bm{\phi}_{n+1}$
$\displaystyle=\bm{\varphi}_{n}+\frac{n-1}{n+2}\bm{\phi}_{n}.$ (43)

We can compute $\nabla E$ analytically as follows. The gradient of
$E(\bm{\phi})$ with respect to $\phi_{j}$ is given by

$\displaystyle\partial_{\phi_{j}}E=2\sum_{l=1}^{N}\mathrm{Re}\left(\partial_{\phi_{j}}\Phi^{(N)}(\rho
z_{l};\bm{\phi})\overline{\Phi^{(N)}(\rho z_{l};\bm{\phi})}\right),$ (44)

and the term $\partial_{\phi_{j}}\Phi^{(N)}(\rho z_{l};\bm{\phi})$ can be
computed using

$\displaystyle\partial_{\phi_{j}}\Phi^{(N)}(\rho
z_{l};\bm{\phi})=2\sum_{i=1}^{3}\partial X_{i}^{(N)}(\rho
z_{l};\bm{\phi})\partial_{\phi_{j}}\partial X_{i}^{(N)}(\rho
z_{l};\bm{\phi}).$ (45)

Moreover, we have

$\displaystyle\partial_{\phi_{j}}\partial X_{i}^{(N)}(\rho
z_{l};\bm{\phi})=\sum_{k=1}^{N}\partial_{\phi_{j}}Q_{ik}(\bm{\phi})\partial
G(\rho z_{l}-\zeta_{k}).$ (46)

Owing to Eq. (30), $\partial_{\phi_{j}}Q_{ik}(\bm{\phi})$ can be obtained
analytically as

$\displaystyle\partial_{\phi_{j}}Q_{ik}(\bm{\phi})=G_{kj}^{-1}b_{i}^{\prime}(\phi_{j}).$
(47)

Hence, we can obtain $\nabla E$ explicitly.

##  4 Convergence and error analysis

We first show the existence of an approximate solution for the minimization
problem with a given precision $\varepsilon>0$.

######  Theorem 4.1.

Let $\Gamma\subset{\mathbb{R}}^{3}$ be a rectifiable closed Jordan curve with
a fixed homeomorphism $b:\partial B\to\Gamma$. Suppose the Fourier series of
$b_{i}$ is absolutely convergent for each $i=1,2,3$. Let
$X^{(N)}(\cdotp;\bm{\phi}):\bar{B}\to{\mathbb{R}}^{3}$ be the approximate
surface for a given $N\in{\mathbb{Z}}_{\geq 1}$ and a given
$\bm{\phi}\in\mathbb{T}^{N}$. Let
$\Phi^{(N)}(\cdotp;\bm{\phi}):\bar{B}\to{\mathbb{C}}$ be the complex
dilatation of $X^{(N)}$.

  1. 1.

For any $\varepsilon>0$ and any sufficiently large $N\in{\mathbb{Z}}_{\geq
1}$, there exists $\bm{\varphi}=\bm{\varphi}(\varepsilon,N)\in\mathbb{T}^{N}$
such that

$\displaystyle\|\Phi^{(N)}(\cdotp;\bm{\varphi})\|_{L^{\infty}(B)}<\varepsilon.$
(48)

We call $\bm{\varphi}$ $\varepsilon$-conformal configulation.

  2. 2.

If the Fourier coefficient $b_{i,n}$ of $b_{i}$ satisfies
$b_{i,n}=\mathrm{O}(|n|^{-\alpha})$ for each $i=1,2,3$ and some $\alpha>2$,
then for any sufficiently large $N\in{\mathbb{Z}}_{\geq 1}$, there exists
$\bm{\varphi}=\bm{\varphi}(N)\in\mathbb{T}^{N}$ such that

$\displaystyle\|\Phi^{(N)}(\cdotp;\bm{\varphi})\|_{L^{\infty}(B)}=\mathrm{O}(N^{-\alpha+1})\quad(N\to\infty).$
(49)

  3. 3.

If $b_{i}$ is real analytic for each $i=1,2,3$, then there exists
$\bm{\varphi}=\bm{\varphi}(N)\in\mathbb{T}^{N}$ and constants $C>0$ and
$\tau\in(0,1)$, independent of $N$, such that

$\displaystyle\|\Phi^{(N)}(\cdotp;\bm{\varphi})\|_{L^{\infty}(B)}\leq
C\tau^{N}.$ (50)

###### Proof.

Take a minimal surface $X\in C^{0}(\bar{B},{\mathbb{R}}^{3})\cap
C^{2}(B,{\mathbb{R}}^{3})$ spanned by $\Gamma$ and define $\varphi:\partial
B\to\partial B$ by $\varphi=b^{-1}\circ X|_{\partial B}$. In particular, we
now have $\Phi_{X}=0$. Setting $\bm{\varphi}=(\varphi(z_{j}))_{j=1}^{N}$, we
show this is the desired object.

Since $X^{(N)}(\cdotp;\bm{\varphi})$ is constructed by the MFS, we deduce from
Theorem 2.6 that for any $i=1,2,3$,

$\displaystyle\|\partial X_{i}-\partial
X_{i}^{(N)}(\cdotp;\bm{\varphi})\|_{L^{\infty}(B)}<1.$ (51)

Hence, we see that for each $z\in B$,

$\displaystyle|\Phi^{(N)}(z;\bm{\varphi})|$
$\displaystyle=|\Phi_{X}(z)-\Phi^{(N)}(z;\bm{\varphi})|$ (52)
$\displaystyle\leq\sum_{i=1}^{3}|\partial X_{i}(z)^{2}-\partial
X_{i}^{(N)}(z;\bm{\varphi})^{2}|$ (53)
$\displaystyle\leq\sum_{i=1}^{3}|\partial X_{i}(z)-\partial
X_{i}^{(N)}(z;\bm{\varphi})|(2|\partial X_{i}(z)|+|\partial X_{i}(z)-\partial
X_{i}^{(N)}(z;\bm{\varphi})|)$ (54) $\displaystyle\leq\sum_{i=1}^{3}\|\partial
X_{i}-\partial X_{i}^{(N)}(\cdotp;\bm{\varphi})\|_{L^{\infty}(B)}(2\|\partial
X_{i}\|_{L^{\infty}(B)}+\|\partial X_{i}-\partial
X_{i}^{(N)}(\cdotp;\bm{\varphi})\|_{L^{\infty}(B)})$ (55) $\displaystyle\leq
C\|\partial X_{i}-\partial
X_{i}^{(N)}(\cdotp;\bm{\varphi})\|_{L^{\infty}(B)}.$ (56)

Hence, the decay of $\|\Phi^{(N)}(\cdotp;\bm{\varphi})\|_{L^{\infty}(B)}$
follows directly from the decay of $\|\partial X_{i}-\partial
X_{i}^{(N)}(\cdotp;\bm{\varphi})\|_{L^{\infty}(B)}$ given in Theorem 2.6. ∎

We next see that the approximate surface for an $\varepsilon$-conformal
configuration gives a minimal surface as $N\to\infty$.

######  Theorem 4.2.

Let $\Gamma\subset{\mathbb{R}}^{3}$ be a rectifiable closed Jordan curve with
a fixed homeomorphism $b:\partial B\to\Gamma$. Suppose the Fourier series of
$b_{i}$ is absolutely convergent for each $i=1,2,3$. Let
$\bm{\varphi}\in\mathbb{T}^{N}$ be an $\varepsilon$-conformal configuration
for a given $\varepsilon>0$ and a sufficiently large $N\in{\mathbb{Z}}_{\geq
1}$. Let $X^{(N)}(\cdotp;\bm{\varphi}):\bar{B}\to{\mathbb{R}}^{3}$ be the
approximate surface. Then, there exists a minimal surface $X\in
C^{0}(\bar{B},{\mathbb{R}}^{3})\cap C^{2}(B,{\mathbb{R}}^{3})$ spanned by
$\Gamma$ and a subsequence $\\{X^{(N_{p})}\\}_{p\geq 1}$ such that

$\displaystyle\lim_{p\to\infty}D(X^{(N_{p})})=D(X).$ (57)

Moreover, if $X$ is non-singular, i.e.,
$\det(\partial_{i}X\cdotp\partial_{j}X)\neq 0$, there exists a constant $C>0$,
independent of $p$ and $\varepsilon$, such that

$\displaystyle\|H_{X^{(N_{p})}}\|_{L^{\infty}(B)}<C\varepsilon.$ (58)

###### Proof.

Theorem 2.6 gives the uniform convergence of $\Gamma_{N}=X^{(N)}(\partial B)$
to $\Gamma$. Hence, applying Theorem 2.5 to $\\{X^{(N)}\\}$, we obtain the
existence of a subsequence $\\{X^{(N_{p})}\\}_{p\geq 1}$ and $X\subset
C^{0}(\bar{B},{\mathbb{R}}^{3})\cap C^{2}(B,{\mathbb{R}}^{3})$ with $\triangle
X=0$ and $X(\partial B)=\Gamma$ such that $X^{(N_{p})}\to X$ uniformly on
$\bar{B}$ as $p\to\infty$. Moreover, owing to the maximum principle, we see
that there exists a constant $C>0$, independent $p$, such that

$\displaystyle\|X-X^{(N_{p})}\|_{W^{1,\infty}(B)}\leq
C\|X-X^{(N_{p})}\|_{L^{\infty}(\partial B)}\to 0\quad\text{as }p\to\infty.$
(59)

We thus deduce that

$\displaystyle\|\Phi_{X}\|_{L^{\infty}(B)}$
$\displaystyle\leq\|\Phi_{X}-\Phi\|_{L^{\infty}(B)}+\|\Phi_{X^{(N)}}\|_{L^{\infty}(B)}$
(60) $\displaystyle\leq
C\|X-X^{(N_{p})}\|_{W^{1,\infty}(B)}+\|\Phi_{X^{(N)}}\|_{L^{\infty}(B)}$ (61)
$\displaystyle\leq 2\varepsilon.$ (62)

Since $\varepsilon>0$ is arbitrary, we obtain
$\|\Phi_{X}\|_{L^{\infty}(B)}=0$, which yields that $X$ is a minimal surface
spanned by $\Gamma$. We now have

$\displaystyle|D(X)-D(X^{(N_{p})}|$
$\displaystyle\leq\|X-X^{(N_{p})}\|_{W^{1,2}(B)}$ (63) $\displaystyle\leq
C\|X-X^{(N_{p})}\|_{W^{1,\infty}(B)}$ (64) $\displaystyle\leq
C\|X-X^{(N_{p})}\|_{L^{\infty}(\partial B)}\to 0\quad\text{as }p\to\infty.$
(65)

Lastly, we examine the convergence of the mean curvature. Define the first
fundamental forms $g_{ij}^{(N)}$ as

$\displaystyle
g_{ij}^{(N)}\coloneqq\langle\partial_{i}X^{(N)},\partial_{j}X^{(N)}\rangle,\quad
i,j=1,2.$ (66)

Since $\|\Phi_{X^{(N)}}\|_{L^{\infty}(B)}<\varepsilon$ holds, we have

$\displaystyle\|g_{11}^{(N)}-g_{22}^{(N)}\|_{L^{\infty}(B)}\leq
C\varepsilon,\quad\|g_{12}^{(N)}\|_{L^{\infty}(B)}\leq C\varepsilon$ (67)

for sufficiently large $N$. Here and hereafter, $C$ is a positive constant
that can change with each appearance and does not depend on $N$. Since
$g=(g_{ij})$ is non-singular, for sufficiently large $N$, we have

$\displaystyle\det g^{(N)}\geq C.$ (68)

A pointwise error of the mean curvature is given by

$\displaystyle|H_{X}^{(N)}|$ $\displaystyle=|H_{X}-H_{X^{(N)}}|$ (69)
$\displaystyle=\left|\frac{g_{11}h_{22}+g_{22}h_{11}-2g_{12}h_{12}}{2\det
g}-\frac{g_{11}^{(N)}h_{22}^{(N)}+g_{22}^{(N)}h_{11}^{(N)}-2g_{12}^{(N)}h_{12}^{(N)}}{2\det
g^{(N)}}\right|,$ (70)

where

$\displaystyle
h_{ij}^{(N)}=\langle\partial_{i}\partial_{j}X^{(N)},e^{(N)}\rangle,\quad
e^{(N)}=\frac{\partial_{1}X^{(N)}\times\partial_{2}X^{(N)}}{\|\partial_{1}X^{(N)}\times\partial_{2}X^{(N)}\|}.$
(71)

By repeatedly applying the very simple equation

$\displaystyle ac-bd=\frac{1}{2}\left[(a+b)(c-d)+(a-b)(c+d)\right],$ (72)

the error in the mean curvature can be estimated by evaluating
$\|g_{ij}-g_{ij}^{(N)}\|_{L^{\infty}(B)}$ and
$\|h_{ij}-h_{ij}^{(N)}\|_{L^{\infty}(B)}$. Since the estimate of
$\|g_{ij}-g_{ij}^{(N)}\|_{L^{\infty}(B)}$ have already been obtained, we will
consider $\|h_{ij}-h_{ij}^{(N)}\|$. By applying the product read as inner
product in equation (72) and the Cauchy–Schwarz inequality, we have

$\displaystyle|h_{ij}-h_{ij}^{(N)}|$
$\displaystyle\leq\frac{1}{2}\left[|\langle\partial_{i}\partial_{j}X-\partial_{i}\partial_{j}X^{(N)},e+e^{(N)}\rangle|+|\langle\partial_{i}\partial_{j}X+\partial_{i}\partial_{j}X^{(N)},e-e^{(N)}\rangle|\right]$
(73)
$\displaystyle\leq\frac{1}{2}\left[|\partial_{i}\partial_{j}X-\partial_{i}\partial_{j}X^{(N)}||e+e^{(N)}|+|\partial_{i}\partial_{j}X+\partial_{i}\partial_{j}X^{(N)}||e-e^{(N)}|\right]$
(74)
$\displaystyle\leq\frac{3}{2}\left[\|\partial_{i}\partial_{j}X-\partial_{i}\partial_{j}X^{(N)}\|_{L^{\infty}(B)}\|e+e^{(N)}\|_{L^{\infty}(B)}\right.$
(75) $\displaystyle\hskip
100.0pt\left.+\|\partial_{i}\partial_{j}X+\partial_{i}\partial_{j}X^{(N)}\|_{L^{\infty}(B)}\|e-e^{(N)}\|_{L^{\infty}(B)},\right]$
(76)

where $\|F\|_{L^{\infty}(B)}$ for an $\mathbb{R}^{3}$-valued function
$F=(F_{1},F_{2},F_{3})$ is defined as

$\displaystyle\|F\|_{L^{\infty}(B)}=\max\\{\|F_{1}\|_{L^{\infty}(B)},\|F_{2}\|_{L^{\infty}(B)},\|F_{3}\|_{L^{\infty}(B)}\\}.$
(77)

Since
$\|\partial_{i}\partial_{j}X-\partial_{i}\partial_{j}X^{(N)}\|_{L^{\infty}(B)}$
and $\|e-e^{(N)}\|_{L^{\infty}(B)}$ converge to $0$ as $N\to\infty$, we see
that

$\displaystyle\|h_{ij}-h_{ij}^{(N)}\|_{L^{\infty}(B)}\leq C\varepsilon$ (78)

for sufficiently large $N$. Combining these estimates, we conclude that

$\displaystyle\|H_{X^{(N)}}\|_{L^{\infty}(B)}<C\varepsilon$ (79)

holds for sufficiently large $N$. ∎

##  5 Numerical examples

In this section, we show several results of numerical experiments, which
exemplify the effectiveness of our method. We briefly explain how we obtain an
$\varepsilon$-conformal configuration $\phi^{(N)}$. In what follows, for
several rectifiable closed Jordan curves $b:\theta\in\partial
B\to\Gamma\subset{\mathbb{R}}^{3}$, we perform numerical experiments to
describe the behavior of the proposed optimization method. All were carried
out by using Julia 1.8.0 on a machine with 3.2 GHz Apple M1 Ultra 20 cores,
128 GB memory, OS X 12.5.1. In every computation, the Nesterov iteration is
performed $10^{5}$ times.

###  5.1 Jordan domain in the plane

As a first example, consider the case where the curve $\Gamma$ is embedded in
the plane. In this case, the problem of finding the minimal surface is nothing
more than finding an isometric map from the unit disk onto the Jordan domain
bounded by $\Gamma$.

We compute conformal mappings for the ellipse
$b(\theta)=(2\cos\theta,\sin\theta,0)$ and the Cassini oval
$b(\theta)=(r(\theta)\cos\theta,r(\theta)\sin\theta,0)$, where

$\displaystyle r(\theta)=\sqrt{\cos 2\theta+\sqrt{1.1^{4}-\sin^{2}2\theta}}.$
(80)

For both cases, we choose equidistant 150 points as initial data. It takes
124.28 seconds for the ellipse and 153.29 seconds for the oval to complete the
computation. As a result, $\varepsilon$-conformal configurations in Fig. 1 and
Fig. 2 are obtained. The contour plot of the dilatation $|\Phi_{X^{(N)}}(z)|$
in (b) of the figures confirms the accuracy of the computations. Since
dilatation is now a holomorphic function, it should obey the maximum
principle. For both figures, we see some focusing nodes around the circle
$\\{|z|=0.7\\}$. In the inner side of the focusing nodes, the dilatation is
smaller than $10^{-10}$.

Ideally, the number of the focusing nodes equals $N$ and they are placed on
circle $\\{|z|=\rho\\}$. However, since the computation is carried out with
$\rho=0.87$, the focusing nodes are formed inside $\\{|z|=\rho\\}$. In
addition, the numbers of nodes are 64 for the ellipse and 52 for the oval,
which is less than $N=150$. Hence, the structure of the focusing nodes yields
the discrete version of the maximum principle, but it is not the $N$-points on
circle $\\{|z|=\rho\\}$.

Figure 1: Ellipse.

Figure 2: Cassini oval.

###  5.2 Crown

Computations not included in the planar domain are as accurate and fast as
those included. A minimal surface spanned by a crown-shaped curve
$b(\theta)=(\cos\theta,\sin\theta,0.3\sin n\theta)$ is demonstrated in Fig. 3
for $n=5$. The focusing nodes are on circle $\\{|z|=0.80\\}$ with $\rho=0.9$,
and the dilatation inside the focusing nodes is smaller than $10^{-10}$. It
costs 143.71 seconds in our computation when we choose equidistant 150 points
as initial data. It is worth noting that the focusing nodes form two distinct
regular polygons: $8$-polygon inside and $10$-polygon outside. If we reduce
the number of times of the Nesterov iteration from $10^{5}$ to $10^{4}$, we
still have the same result, which implies the optimization is so fast and the
polygonal structure of the focusing nodes is an intrinsic nature of the crown-
shaped curve.

Figure 3: Crown for $n=5$.

###  5.3 Torus knot

Even if the resulting minimal surfaces can not be embedded but immersed, and
even if a given boundary curve is knotting, we can compute them with the same
accuracy and speed. The torus knot $b(\theta)=((2+\cos q\theta)\cos
p\theta,(2+\cos q\theta)\sin p\theta,-\sin q\theta)$ is a knot and gives such
a minimal surface. Fig. 4 shows the result for $(3,2)$-torus knot. Taking
equidistant 150 points as initial data, we obtain the result in Fig. 4 in
156.88 seconds for $(3,2)$-torus knot with $\rho=0.85$. The focusing nodes are
placed on $\\{|z|=0.70\\}$ with $10^{-10}$ accuracy.

Figure 4: Knot for $(p,q)=(3,2)$.

##  6 Searching methods for all solutions

We finally propose two methods to search for all solutions to the Plateau
problem for a given boundary curve. We demonstrate these methods by taking the
Enneper wire as an example, and the Enneper wire is a closed Jordan curve
given by

$\displaystyle b_{1}(\theta)$ $\displaystyle=r\cos\theta-\frac{r^{3}}{3}\cos
3\theta,$ (81) $\displaystyle b_{2}(\theta)$
$\displaystyle=-r\sin\theta-\frac{r^{3}}{3}\sin 3\theta,$ (82) $\displaystyle
b_{3}(\theta)$ $\displaystyle=r^{2}\cos 2\theta.$ (83)

for $r\in(0,\sqrt{3})$. The Enneper wire comes from the boundary of the
Enneper surface, which is an exact solution to the Plateau problem. The
homeomorphism $\phi:\partial B\to\partial B$ is then the identity map on
$\partial B$. Hence, the equidistant point set is the $\varepsilon$-conformal
configuration corresponding to this exact solution.

Fig. 5 shows the numerical solution for the equidistant points corresponding
to the Enneper surface. We see no focusing nodes in the contour plot of the
dilatation, which implies that focusing nodes cannot be formed without
optimization. It is known that the Enneper wire bounds two distinct minimal
surfaces other than the Enneper surface.

Our goals in this section are to find solutions and to confirm that there is
no other solution than what we have found. First, let us search for solutions
in a unified manner as follows. We take a one-parameter family of initial
configurations
$\bm{\phi}^{(N)}:s\in(S_{1},S_{2})\to\bm{\phi}^{(N)}(s)\in\mathbb{T}^{N}$ with

$\displaystyle\phi_{j}^{(N)}(s)=\frac{2\pi(j-1)}{N}+s\sin\frac{2\pi
m(j-1)}{N}$ (84)

for $s>0$ and $m\in{\mathbb{Z}}$, and $\phi^{(N)}$ is regarded as an
associated configuration of a perturbed identity map on $\partial B$ by a
Fourier mode $s\sin mx$. Fig. 6 indicates that each initial configurations
converges to either of the three minimizers for
$s\in\\{-2.95,-2.9,\ldots,2\\}$ and $m=2$.

Since $\bm{\phi}^{(N)}(s)$ varies continuously for $s$, the value of the
Dirichlet energy also varies continuously. Hence, we can distinguish the
minimizers from connected components of the value of the Dirichlet energy as
$s$ varies. Fig. 7 shows the value of Dirichlet energy for each $s$. The color
of the points is painted according to the accuracy of $E$. We see that the
distribution of the values is well-organized if $E$ is more accurate than
5-digit precision; otherwise, it is scattered. Therefore, it is suggested that
the numerical solutions must be more accurate than at least 5-digit precision
to distinguish minimizers.

We pick out three representatives of the connected components with 5-digit
precision, choosing $s=-1,0,1$. Fig. 8 shows their initial and
$\varepsilon$-conformal configurations for $s=-1,0,1$ in left, center, and
right, respectively. Indeed, it is worth noting that these representatives
give three distinct minimal surfaces in Fig 9. They attain 8-digit precision
by 160 seconds and form focusing nodes, as in Fig. 10. Hence, we can find
solutions by taking $\bm{\phi}^{(N)}(s)$.

Second, let us choose initial configurations at random, thereby check there is
no other solution than what we found above. Picking up points randomly and
interpolating them by B-spline interpolation, we obtain 5000 samples with
5-digit precision. Fig. 11 and Fig. 12 show the distribution and its
rearrangement of the Dirichlet energy according to the accuracy with 5-digit
(left), 6-digit (center), and 7-digit (right) precision.

We observe that there are at least two solutions by distinguishing these
samples by whether the value of the Dirichlet energy is greater than 12 or
not. When the Dirichlet energy of a sample is greater than 12, the sample
corresponds to the case $s=0$. Otherwise, they correspond to either the case
$s=1$ or $s=-1$. Note that the case for $s=1,-1$ cannot be distinguished by
the Dirichlet energy since the minimal surfaces for $s=1,-1$ coincide by
certain rigid rotation. However, we can see from the shape of the minimal
surface that there is no solution other than what we have found by taking
random initial configurations.

##  7 Concluding remarks

We propose a numerical scheme with such a high speed and high accuracy as to
find all minimal surfaces spanned by a closed Jordan curve. The numerical
scheme is based on the method of fundamental solutions (MFS) and the
Nesterov’s accelerated gradient decent. We can compute an approximate
solutions for the Dirichlet boundary value problem by the MFS, and the error
of the approximate solution to an exact solution decays exponentially for
sufficiently smooth boundary data. In the computation of a minimal surface,
the Dirichlet boundary value problem arises to compute the coordinate function
of the minimal surface. However, not all surfaces obtained as the solution for
the boundary value problem are minimal surfaces because the solution varies as
a change of variable of the boundary that gives a different boundary data.
Hence, it is necessary to compute minimal surfaces via the MFS to find a
suitable change of variable that yields minimal surfaces. We proposed a
minimization problem for discrete energy of the complex dilatation around the
boundary. The significant characteristic is that every MFS approximate
solution is smooth and harmonic on the whole domain, eliminating the need for
the integration on the whole space of the functional. Nesterov’s accelerated
gradient decent solves this minimization problem.

For the proposed numerical method, we proved the existence of a solution for
the minimization problem for the complex dilatation in Theorem 4.1. We
obtained the error estimate of the complex dilatation according to the
regularity of a given boundary. We next showed the existence of a subsequence
of a given sequence of approximate surfaces and the limit. In particular, the
limit became a minimal surface and the value of the Dirichlet energy
convergeed to the one for the minimal surface. We also obtained the
$L^{\infty}$-error estimate for the mean curvature.

We investigated the performance of numerical computation for the proposed
minimal surface. We saw from the computation for some boundaries that the
complex dilatation satisfied a discrete version of the maximum principle. The
errors of the complex dilatation were with 10-digit precision for at most 160
seconds for each of boundaries in our computation. Lastly, we proposed two
methods of finding all minimal surfaces spanned by a given closed Jordan
curve. We chose a one-parameter family of initial configurations in the first
method based on a perturbed identity map on the boundary by a Fourier mode.
Since the initial configurations changed continuously, we could find distinct
minimal surfaces as a minimizer of boundary mapping or connected components of
the value of the Dirichlet energy. In the second method, we selected the
initial configuration at random. We can see from the distribution of the
Dirichlet energy or the surface shape that there is no solution other than
what we found. It remains for future work to determine whether the above
methods cannot find some minimal surfaces or how many samples must be
calculated to ensure that all minimal surfaces have been obtained.

### Acknowledgements

This work was supported by JSPS KAKENHI Grant Numbers JP18K13455, JP22K03425,
JP21J00025.

Figure 5: Enneper surface for $R=1.1$.

Figure 6: $\bm{\phi}^{(N)}(s)$ (left) and $\varepsilon$-conformal
configuration (right). Figure 7: Dirichlet energy.

Figure 8: Initial and $\varepsilon$-conformal configurations.

Figure 9: Minimal surfaces.

Figure 10: Contour plots.

Figure 11: Distribution of the Dirichlet energy.

Figure 12: Rearranged Distribution of the Dirichlet energy.

## References

  * [1] K. Amano, A bidirectional method for numerical conformal mapping based on the charge simulation method, J. Inform. Process., 14 (1991), pp. 473–482.
  * [2] K. Amano, A charge simulation method for the numerical conformal mapping of interior, exterior and doubly-connected domains, J. Comput. Appl. Math., 53 (1994), pp. 353–370.
  * [3] K. Amano, A charge simulation method for numerical conformal mapping onto circular and radial slit domains, SIAM J. Sci. Comput., 19 (1998), pp. 1169–1187.
  * [4] K. Amano, D. Okano, H. Ogata, and M. Sugihara, Numerical conformal mappings onto the linear slit domain, Jpn. J. Ind. Appl. Math., 29 (2012), pp. 165–186.
  * [5] R. Böhme and A. J. Tromba, The index theorem for classical minimal surfaces, Ann. of Math. (2), 113 (1981), pp. 447–499.
  * [6] R. Courant, Plateau’s problem and dirichlet’s principle, Ann. of Math. (2), 38 (1937), pp. 679–724.
  * [7] U. Dierkes, S. Hildebrandt, and F. Sauvigny, Minimal surfaces, vol. 339 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], Springer, Heidelberg, second ed., 2010.  With assistance and contributions by A. Küster and R. Jakob.
  * [8] J. Douglas, A method of numerical solution of the problem of Plateau, Ann. of Math. (2), 29 (1927/28), pp. 180–188.
  * [9] , Solution of the problem of Plateau, Trans. Amer. Math. Soc., 33 (1931), pp. 263–321.
  * [10] G. Dziuk and J. E. Hutchinson, The discrete Plateau problem: algorithm and numerics, Math. Comp., 68 (1999), pp. 1–23.
  * [11] , The discrete Plateau problem: convergence results, Math. Comp., 68 (1999), pp. 519–546.
  * [12] , Finite element approximations to surfaces of prescribed variable mean curvature, Numer. Math., 102 (2006), pp. 611–648.
  * [13] A. Grodet and T. Tsuchiya, Finite element approximations of minimal surfaces: algorithms and mesh refinement, Jpn. J. Ind. Appl. Math., 35 (2018), pp. 707–725.
  * [14] Y.-X. Hao, C.-J. Li, and R.-H. Wang, An approximation method based on MRA for the quasi-Plateau problem, BIT, 53 (2013), pp. 411–442.
  * [15] M. Hinze, On the numerical approximation of unstable minimal surfaces with polygonal boundaries, Numer. Math., 73 (1996), pp. 95–118.
  * [16] M. Katsurada, A mathematical study of the charge simulation method ii, J. Fac. Sci. Univ. Tokyo Sect. IA Math., 36 (1989), pp. 135–162.
  * [17] M. Katsurada and H. Okamoto, A mathematical study of the charge simulation method i, J. Fac. Sci. Univ. Tokyo Sect. IA Math., 35 (1988), pp. 507–518.
  * [18] M. Koiso, On the finite solvability of Plateau’s problem for extreme curves, Osaka J. Math., 20 (1983), pp. 177–183.
  * [19] J. Lagrange, Essai dúne nouvelle méthode pour determiner les maxima et les minima des formules intégrales indéfinies., Misc. Philos.-Math. Soc. Priv. Taurinensis, 2 (1760–1762), pp. 173–195.
  * [20] J. C. C. Nitsche, A new uniqueness theorem for minimal surfaces, Arch. Rational Mech. Anal., 52 (1973), pp. 319–329.
  * [21] , Contours bounding at most finitely many solutions of Plateau’s problem, in Complex analysis and its applications (Russian), “Nauka”, Moscow, 1978, pp. 438–446, 670.
  * [22] , Lectures on minimal surfaces. Vol. 1, Cambridge University Press, Cambridge, 1989.  Introduction, fundamentals, geometry and basic boundary value problems, Translated from the German by Jerry M. Feinberg, With a German foreword.
  * [23] P. Pozzi, $L^{2}$-estimate for the discrete Plateau problem, Math. Comp., 73 (2004), pp. 1763–1777.
  * [24] T. Radó, Some remarks on the problem of plateau., Proc. Natl. Acad. Sci. USA, 16 (1930), pp. 242–248.
  * [25] K. Sakakibara, Bidirectional numerical conformal mapping based on the dipole simulation method, Eng. Anal. Bound. Elem., 114 (2020), pp. 45–57.
  * [26] H. Schumacher and M. Wardetzky, Variational convergence of discrete minimal surfaces, Numer. Math., 141 (2019), pp. 173–213.
  * [27] L. Tomek and K. Mikula, Discrete duality finite volume method with tangential redistribution of points for surfaces evolving by mean curvature, ESAIM Math. Model. Numer. Anal., 53 (2019), pp. 1797–1840.
  * [28] O. y. Trå sdahl and E. M. Rø nquist, High order numerical approximation of minimal surfaces, J. Comput. Phys., 230 (2011), pp. 4795–4810.
  * [29] T. Tsuchiya, Discrete solution of the Plateau problem and its convergence, Math. Comp., 49 (1987), pp. 157–165.
  * [30] , A note on discrete solutions of the Plateau problem, Math. Comp., 54 (1990), pp. 131–138.
  * [31] S. Wang and A. Chern, Computing minimal surfaces with differential forms, ACM Trans. Graph., 40 (2021).