# An Iterative Least-Squares Method for the Hyperbolic Monge-Ampère Equation
with Transport Boundary Condition

M.W.M.C. Bertens  CASA, Department of Mathematics and Computer Science,
Eindhoven University of Technology, PO Box 513, 5600 MB Eindhoven, The
Netherlands  Corresponding author<EMAIL_ADDRESS>M.J.H. Anthonissen
CASA, Department of Mathematics and Computer Science, Eindhoven University of
Technology, PO Box 513, 5600 MB Eindhoven, The Netherlands
J.H.M. ten Thije Boonkkamp  CASA, Department of Mathematics and Computer
Science, Eindhoven University of Technology, PO Box 513, 5600 MB Eindhoven,
The Netherlands  W.L. IJzerman  CASA, Department of Mathematics and Computer
Science, Eindhoven University of Technology, PO Box 513, 5600 MB Eindhoven,
The Netherlands  Signify Research, High Tech Campus 7, 5656 AE Eindhoven, The
Netherlands

###### Abstract

A least-squares method for solving the hyperbolic Monge-Ampère equation with
transport boundary condition is introduced. The method relies on an iterative
procedure for the gradient of the solution, the so-called mapping. By
formulating error functionals for the interior domain, the boundary, both
separately and as linear combination, three minimization problems are solved
iteratively to compute the mapping. After convergence, a fourth minimization
problem, to compute the solution of the Monge-Ampère equation, is solved. The
approach is based on a least-squares method for the elliptic Monge-Ampère
equation[1], and is improved upon by the addition of analytical solutions for
the minimization on the interior domain and by the introduction of two new
boundary methods. Lastly, the iterative method is tested on a variety of
examples. It is shown that, when the iterative method converges, second-order
global convergence as function of the spatial discretization is obtained.

Keywords— Hyperbolic Monge-Ampère equation, Transport boundary condition,
Iterative least-squares method

##  1 Introduction

In this paper we introduce a least-squares method for the hyperbolic Monge-
Ampère equation with transport boundary condition. We are motivated by
applications to optical design. In [2] it was found that designing lenses and
reflectors for some single-optical-surface systems is equivalent to solving
the Monge-Ampère (MA) equation with transport boundary condition. The optical
surface satisfies either the elliptic (+) or hyperbolic (-) MA equation, i.e.,

$\displaystyle\det(\mathrm{\mathbf{D}}^{2}u(\mathbf{x}))=\pm\frac{E(\mathbf{x})}{I(\nabla
u(\mathbf{x}))},\qquad\mathbf{x}\in\mathcal{X},$ (1a) where
$\mathrm{\mathbf{D}}^{2}u$ is the Hessian matrix of the optical surface
$z=u(\mathbf{x})$, $E\geq 0$ the emittance of the source with domain
$\mathcal{X}\subset\mathbb{R}^{2}$ and $I>0$ the illuminance on the target
with domain $\mathcal{Y}=\nabla u(\mathcal{X})\subset\mathbb{R}^{2}$. The
accompanying transport boundary condition is given by $\displaystyle\nabla
u(\partial\mathcal{X})=\partial\mathcal{Y}.$ (1b)

In the elliptic case this boundary condition follows from convexity or
concavity arguments of the optical surface [2]. For the hyperbolic MA
equation, the transport boundary condition cannot be derived in the same
manner as $z=u(\mathbf{x})$ is a saddle surface, and one instead relies on
optical arguments, e.g., the edge-ray principle [3]. Further disparities
between the elliptic and hyperbolic variant are common. This for one is due to
the connection of the elliptic variant to the rich field of optimal transport
(OT) [4, p. 282], which is absent for the hyperbolic MA equation. OT was
originally established by Monge, who was concerned with rearranging mass from
one distribution to another[5, p. xiv]. Brenier proved that, assuming
regularity conditions [6], the corresponding optimal (point-wise) transport
map $\mathbf{m}=\nabla u$ satisfies the elliptic Monge-Ampère equation (1a),
where the plus sign is assumed and $E$ and $I$ should be interpreted as
densities. See [4, p. 323-332], for example, for regularity, uniqueness and
existence of solutions to the elliptic MA equation. Assumptions on energy
conservation allow design of optical systems to be cast in the framework of
OT, but only for the elliptic variant and not for the hyperbolic MA equation.
As a consequence, results regarding regularity, uniqueness and existence for
the hyperbolic MA equation are scarce. The most relevant results for the
hyperbolic MA equation

$\displaystyle\det(\mathrm{\mathbf{D}}^{2}u(\mathbf{x}))=-f^{2}(\mathbf{x}),\qquad\mathbf{x}\in\mathcal{X}\subset\mathbb{R}^{2},$
(2)

with $f:\mathbb{R}^{2}\mapsto(0,\infty)$, follow from the method of
characteristics [7, 8] and only hold for Cauchy boundary conditions. Equation
(2) has received little interest in numerical and computational journals. To
the best of our knowledge, the equation has been solved twice, once on a
triangular computational domain [9, p. 614] and secondly on a rectangular
computational domain by Bertens et al. [10]. The former method, based on a
finite difference scheme, assumes Cauchy boundary conditions on an initial
curve, and does not treat boundary conditions on the rest of the domain. The
latter method, derived using the method of characteristics, shows that the
requirements on the boundary conditions are strenuous. The method of
characteristics shows, by parametrizing the characteristics with the
$x$-coordinate and assuming Cauchy boundary conditions on an initial curve
$x=$ const, that the remaining boundary conditions depend on the location of
the characteristics. Consequently, the required boundary conditions are
fundamentally different from the transport boundary condition.

In this work, we therefore do not consider the method of characteristics, but
instead resort to a least-squares method which has been proven to work for
various elliptic problems, among which the Monge-Ampère equation [1], the
generalized Monge-Ampère equation and the generated Jacobian equation [11].
The least-squares algorithm is an iterative method which does not directly
solve for the unknown $u$, but instead first constructs the mapping
$\mathbf{m}=\nabla u$ and afterwards approximates $u$. The general outline of
the least-squares method is as follows: first, we approximate the Jacobi
matrix of $\mathbf{m}$ in the interior of the domain by minimizing an error
functional. Secondly, $\mathbf{m}$ restricted to the boundary of the domain is
approximated. By minimizing another error functional involving the newly found
Jacobi matrix and the boundary approximation, we obtain a new approximation
for the mapping. We proceed by repeating these three steps iteratively until
$\mathbf{m}$ no longer changes and subsequently calculate $u$ by minimizing a
fourth functional. One of the benefits of this method is that each of the
three stages can be adapted for the problem at hand. For example, the
minimization for $\mathbf{m}$ in [1] relies on a finite difference scheme
while in [12] it uses a finite volume scheme. This three-stage approach allows
us to introduce two new boundary methods, viz. a segmented projection method
and a segmented arc length method, which both lead to better results and
higher computational efficiency than the original projection method[1]. Even
more importantly, the iterative method using the segmented projection method
converges in some cases when the original projection method does not. And, as
we will show, the segmented arc length method converges for all examples.
Furthermore, we improve upon the first minimization procedures, viz. the
procedure for approximating the Jacobi matrix in the interior of the domain.
Numerical experiments have shown that grid lines in target space can
intersect, preventing proper numerical convergence of our algorithm.
Therefore, we introduce a method to prevent these so-called grid shocks.

The content of this paper is as follows. We discuss the theory of the least-
squares method for the Monge-Ampère equation in Section 2. In Section 2.1 the
least-squares method is introduced. Afterwards, we adapt part of the method,
viz. the optimization in the interior domain, in Section 2.2. Next, in Section
2.3, we introduce various boundary methods to replace the existing projection
method and in Section 2.4 we introduce a grid shock correction method. In
Section 3 we compare the boundary methods, show their weaknesses and strengths
and elaborate on the convergence of the algorithm for various test cases.
Lastly, we end with a discussion of the results followed by conclusions in
Section 4.

##  2 The least-squares formulation

We are interested in the two-dimensional hyperbolic Monge-Ampère equation with
transport boundary condition, given by

$\displaystyle\det{\left(\mathrm{\mathbf{D}}^{2}u(\mathbf{x})\right)}+f^{2}(\mathbf{x},\nabla
u(\mathbf{x}))$ $\displaystyle=0,\quad\mathbf{x}\in\mathcal{X},$ (3a)
$\displaystyle\nabla u(\partial\mathcal{X})$
$\displaystyle=\partial\mathcal{Y},$ (3b)

where $u=u(\mathbf{x})$ is the unknown, $\mathrm{\mathbf{D}}^{2}u$ the Hessian
matrix of $u$, $f^{2}>0$ and $\mathcal{X},\mathcal{Y}\subset\mathbb{R}^{2}$
connected domains. We require the boundaries $\partial\mathcal{X}$ and
$\partial\mathcal{Y}$ to be orientable. The transport boundary condition (3b)
can be interpreted as

$\displaystyle\begin{cases}\forall\mathbf{x}\in\partial\mathcal{X}:\nabla
u(\mathbf{x})\in\partial\mathcal{Y},\\\
\forall\mathbf{y}\in\partial\mathcal{Y}\,\,\,\exists\mathbf{x}\in\partial\mathcal{X}:\nabla
u(\mathbf{x})=\mathbf{y},\end{cases}$ (4a)

where the latter condition is recognized as surjectivity of $\nabla u$.
Bijectivity is generally not implied, not even when restricted to the
boundary, as will become apparent by the example discussed in Section 3.4.
Hyperbolicity of (3a) follows from the discriminant of the characteristic
condition, which can be obtained by rewriting (3a) as

$\displaystyle F(\mathbf{x},u,p,q,r,s,t)=rt-s^{2}+f^{2}=0,$ (5)

where $p=u_{x_{1}}$, $q=u_{x_{2}}$, $r=u_{x_{1}x_{1}}$, $s=u_{x_{1}x_{2}}$ and
$t=u_{x_{2}x_{2}}$. The characteristic condition is given by [10, p. 10]

$\displaystyle F_{r}\mu^{2}-F_{s}\mu+F_{t}=0,$ (6)

for the unknown function $\mu$, representing the slope of the characteristics.
For the MA equation to be hyperbolic, two real characteristics need to exist
for every point in the domain, hence the slopes of the two characteristics,
and thus the roots of (6), need to be real and distinct. Henceforth, the
discriminant of (6) should be strictly positive. It follows that the
discriminant $\Delta$ of (6) is given by

$\displaystyle\Delta=F_{s}^{2}-4F_{r}F_{t}=4s^{2}-4tr=4f^{2},$ (7)

which is, by assumption, strictly positive. Hence, equation (3a) is
hyperbolic.

###  2.1 Least-squares approach

In [1] a least-squares method was introduced to solve the elliptic Monge-
Ampère equation given by
$\det{\left(\mathrm{\mathbf{D}}\mathbf{m}\right)}=f^{2}(\mathbf{x},\nabla
u(\mathbf{x}))$ for $\mathbf{x}\in\mathcal{X}$ and
$\mathrm{\mathbf{D}}\mathbf{m}$ the Jacobi matrix of $\mathbf{m}$. The main
idea of the least-squares method is to reformulate the Monge-Ampère equation
in terms of the mapping $\mathbf{m}:\mathcal{X}\rightarrow\mathcal{Y}$,
representing $\nabla u$, and solve for $\mathbf{m}$. Subsequently, $u$ is
reconstructed from $\mathbf{m}$ in a least-squares sense. To solve the
hyperbolic problem we replace the right-hand side of the elliptic Monge-Ampère
equation by $-f^{2}(\mathbf{x},\nabla u(\mathbf{x}))$ and substitute
$\mathrm{\mathbf{D}}\mathbf{m}=\mathrm{\mathbf{D}}^{2}u$, thus obtaining

$\displaystyle\det{\left(\mathrm{\mathbf{D}}\mathbf{m}(\mathbf{x})\right)}+f^{2}(\mathbf{x},\mathbf{m}(\mathbf{x}))$
$\displaystyle=0,\quad\mathbf{x}\in\mathcal{X},$ (8a)
$\displaystyle\mathbf{m}(\partial\mathcal{X})$
$\displaystyle=\partial\mathcal{Y}.$ (8b)

We formulate a minimization problem for $\mathbf{m}$ which we solve
numerically. For this, we introduce the auxiliary functions
$\mathbf{P}:\mathcal{X}\rightarrow\mathbb{R}^{2\times 2}$ and
$\mathbf{b}:\partial\mathcal{X}\rightarrow\partial\mathcal{Y}$ which are used
to approximate $\mathrm{\mathbf{D}}\mathbf{m}$ on the whole domain and
$\mathbf{m}$ on the boundary, respectively. This is achieved by the least-
squares method, i.e., subsequently minimizing three separate functionals given
by

$\displaystyle J_{\textrm{I}}(\mathbf{m},\mathbf{P})$
$\displaystyle=\frac{1}{2}\iint_{\mathcal{X}}\|\mathrm{\mathbf{D}}\mathbf{m}-\mathbf{P}\|^{2}\mathop{}\\!\mathrm{d}\mathbf{x},$
(9a) $\displaystyle J_{\textrm{B}}(\mathbf{m},\mathbf{b})$
$\displaystyle=\frac{1}{2}\oint_{\partial\mathcal{X}}|\mathbf{m}-\mathbf{b}|^{2}\mathop{}\\!\mathrm{d}s,$
(9b) $\displaystyle J(\mathbf{m},\mathbf{P},\mathbf{b})$ $\displaystyle=\alpha
J_{\textrm{I}}(\mathbf{m},\mathbf{P})+(1-\alpha)J_{\textrm{B}}(\mathbf{m},\mathbf{b}),$
(9c)

where $|\cdot|$ is the standard 2-norm, $\|\cdot\|$ is the Frobenius norm
defined by $\|\mathbf{A}\|^{2}=\mathrm{Tr}(\mathbf{A}\mathbf{A}^{\text{T}})$
for a matrix $\mathbf{A}$ and $0<\alpha<1$ is a control parameter to either
place weights on the boundary and the interior. Starting with an initial guess
$\mathbf{m}^{0}$, the iterative optimization procedure for $n=0,1,2,\dots$
reads

$\displaystyle\mathbf{P}^{n+1}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{P}\in\mathcal{P}(\mathbf{m}^{n})}J_{\textrm{I}}(\mathbf{m}^{n},\mathbf{P}),$
(10a) $\displaystyle\mathbf{b}^{n+1}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{b}\in\mathcal{B}}J_{\textrm{B}}(\mathbf{m}^{n},\mathbf{b}),$
(10b) $\displaystyle\mathbf{m}^{n+1}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{m}\in\mathcal{V}}J(\mathbf{m},\mathbf{P}^{n+1},\mathbf{b}^{n+1}).$
(10c)

The spaces $\mathcal{P}(\mathbf{m}^{n})$, $\mathcal{B}$ and $\mathcal{V}$
follow from three key observations. First, because $\mathbf{m}=\nabla u$, the
Jacobi matrix $\mathrm{\mathbf{D}}\mathbf{m}=\mathrm{\mathbf{D}}^{2}u$ is
symmetric and
$\det(\mathrm{\mathbf{D}}\mathbf{m})=-f^{2}(\mathbf{x},\mathbf{m}(\mathbf{x}))$.
Secondly, by the transport boundary condition, for all
$\mathbf{x}\in\partial\mathcal{X}:\mathbf{m}(\mathbf{x})\in\partial\mathcal{Y}$.
As we require $\mathbf{m}$ to be twice continuously differentiable later on,
we impose this requirement. The three sets are then given by

$\displaystyle\mathcal{P}(\mathbf{m})$
$\displaystyle=\left\\{\mathbf{P}\in[C^{1}(\mathcal{X})]^{2\times
2}\mid\det(\mathbf{P}(\mathbf{x}))=-f^{2}(\mathbf{x},\mathbf{m}(\mathbf{x})),\mathbf{P}=\mathbf{P}^{\text{T}}\right\\},$
(11a) $\displaystyle\mathcal{B}$
$\displaystyle=\left\\{\mathbf{b}\in[C(\mathcal{\partial
X})]^{2}\mid\mathbf{b}(\mathbf{x})\in\partial\mathcal{Y}\right\\},$ (11b)
$\displaystyle\mathcal{V}$ $\displaystyle=[C^{2}(\mathcal{X})]^{2}.$ (11c)

We first outline the minimization of $J$, as it remains unchanged w.r.t. [1],
and in the next sections we elaborate on the minimization of $J_{\mathrm{I}}$
and $J_{\mathrm{B}}$. Taking the variational derivative of (9c) and applying
the fundamental lemma of calculus of variations [13, p. 185] yields that for
the optimal $\mathbf{m}$, each of its components should satisfy a Poisson
equation with Robin boundary condition given by

$\displaystyle\Delta m_{1}$
$\displaystyle=\nabla\boldsymbol{\cdot}\mathbf{p}_{1},$
$\displaystyle\mathbf{x}\in\mathcal{X},$ (12a)
$\displaystyle(1-\alpha)m_{1}+\alpha\nabla
m_{1}\boldsymbol{\cdot}\hat{\mathbf{n}}$
$\displaystyle=(1-\alpha)b_{1}+\alpha\mathbf{p}_{1}\boldsymbol{\cdot}\hat{\mathbf{n}},$
$\displaystyle\mathbf{x}\in\partial\mathcal{X},$ (12b)

for the first component $m_{1}$ and

$\displaystyle\Delta m_{2}$
$\displaystyle=\nabla\boldsymbol{\cdot}\mathbf{p}_{2},$
$\displaystyle\mathbf{x}\in\mathcal{X},$ (13a)
$\displaystyle(1-\alpha)m_{2}+\alpha\nabla
m_{2}\boldsymbol{\cdot}\hat{\mathbf{n}}$
$\displaystyle=(1-\alpha)b_{2}+\alpha\mathbf{p}_{2}\boldsymbol{\cdot}\hat{\mathbf{n}},$
$\displaystyle\mathbf{x}\in\partial\mathcal{X},$ (13b)

for the second component $m_{2}$. The functions $\mathbf{p}_{i}$ ($i=1,2$)
denote the $i^{\text{th}}$ column of the matrix $\mathbf{P}$ and
$\hat{\mathbf{n}}$ is the unit outward normal vector to $\partial\mathcal{X}$.

Upon convergence of (10) we reconstruct $u$ from $\mathbf{m}$ by minimizing
another least-squares functional, viz.

$\displaystyle u=\operatorname*{argmin}_{\psi\in
C^{2}(\mathcal{X})}\frac{1}{2}\iint_{\mathcal{X}}|\nabla\psi-\mathbf{m}|^{2}\mathop{}\\!\mathrm{d}\mathbf{x}.$
(14)

Using calculus of variations once more, we obtain the Poisson equation with
Neumann boundary conditions for $u$, which reads

$\displaystyle\Delta u$ $\displaystyle=\nabla\boldsymbol{\cdot}\mathbf{m},$
$\displaystyle\qquad\mathbf{x}\in\mathcal{X},$ (15a) $\displaystyle\nabla
u\boldsymbol{\cdot}\hat{\mathbf{n}}$
$\displaystyle=\mathbf{m}\boldsymbol{\cdot}\hat{\mathbf{n}},$
$\displaystyle\qquad\mathbf{x}\in\partial\mathcal{X}.$ (15b)

For (15) to admit a solution, the compatibility condition [14, p. 184]

$\displaystyle\iint_{\mathcal{X}}\nabla\boldsymbol{\cdot}\mathbf{m}\mathop{}\\!\mathrm{d}\mathbf{x}-\oint_{\partial\mathcal{X}}\mathbf{m}\boldsymbol{\cdot}\hat{\mathbf{n}}\mathop{}\\!\mathrm{d}s=0,$
(16)

is automatically satisfied due to the divergence theorem.

We solve the three Poisson equations using finite differences (FD), more
specifically, standard second-order central differences for both the first and
second order derivatives. For grid points on the boundary we introduce ghost
points, which we eliminate using the normal derivatives in the Robin boundary
condition. The system we obtain from discretizing (12) and (13) needs to be
solved in each iteration. In order to increase computational efficiency, we
compute the LU-decomposition in the initialization of the algorithm. Note that
the solution for $u$ is not unique due to the (transport) boundary condition
[15, p. A1438], which is also reflected by (15b), so we enforce uniqueness by
fixing one function value of $u$, i.e., let $\mathbf{x}\in\mathcal{X}$ be
arbitrary, we then impose the condition $u(\mathbf{x})=0$. In practice we
assume
$\mathcal{X}=[x_{\mathrm{m}},x_{\mathrm{M}}]\times[y_{\mathrm{m}},y_{\mathrm{M}}]$
and we impose $u(x_{\mathrm{m}},y_{\mathrm{m}})=0$. Alternatively, one could
prescribe the average value of $u$ on the domain [16, p. 177].

###  2.2 $\mathbf{P}$-optimization

The matrix $\mathrm{\mathbf{D}}\mathbf{m}$ cannot be determined exactly during
the iterative process. Because the integrand of $J_{\textrm{I}}$, i.e.,
$\|\mathrm{\mathbf{D}}\mathbf{m}-\mathbf{P}\|^{2}$, does not depend on
derivatives of $\mathbf{P}$ we employ a piece-wise minimization. To this end
we approximate $\mathrm{\mathbf{D}}\mathbf{m}$ using standard finite
difference. Let $\mathbf{x}_{ij}=((x_{1})_{i},(x_{2})_{j})\in\mathcal{X}$ be
the grid points of a Cartesian grid with $i=1,\dots,N_{x_{1}}$ and $j=1,\dots
N_{x_{2}}$ denoting the first and second coordinate, respectively. We write
$\mathbf{m}_{ij}\approx\mathbf{m}(\mathbf{x}_{ij})$ and similar for the other
variables. We approximate $(\mathrm{\mathbf{D}}\mathbf{m})_{ij}$ by
$\mathbf{D}_{ij}$ using central and one-sided second-order finite differences
in the interior and at the boundary, respectively. This implies that
$\mathbf{D}$ is in general not symmetric, while
$\mathrm{\mathbf{D}}\mathbf{m}$ and $\mathbf{P}$ are. By virtue of the point-
wise minimization we proceed to drop the subscripts, e.g., we write
$\mathbf{m}$ instead of $\mathbf{m}_{ij}$, for brevity.

Let $F(p_{11},p_{22},p_{12})=\|\mathbf{D}-\mathbf{P}\|^{2}$; expanding it
yields

$\displaystyle
F(p_{11},p_{22},p_{12})=\frac{1}{2}\Big{(}(p_{11}-d_{11})^{2}+(p_{12}-d_{12})^{2}+(p_{12}-d_{21})^{2}+(p_{22}-d_{22})^{2}\Big{)}.$
(17)

We replace $\mathbf{D}$ by its symmetric part
$\mathbf{D}_{\mathrm{s}}=\tfrac{1}{2}(\mathbf{D}+\mathbf{D}^{\text{T}})$, or
written in its components, we introduce
$d_{\textrm{s}}=\tfrac{1}{2}(d_{12}+d_{21})$ and

$\displaystyle\mathbf{D}_{\textrm{s}}=\begin{pmatrix}d_{11}&d_{\textrm{s}}\\\
d_{\textrm{s}}&d_{22}\end{pmatrix}.$ (18)

Furthermore, we replace $F$ by
$F_{\textrm{s}}=\tfrac{1}{2}\|\mathbf{P}-\mathbf{D}_{\textrm{s}}\|^{2}$, i.e.,

$\displaystyle
F_{\textrm{s}}(p_{11},p_{22},p_{12})=\frac{1}{2}\Big{(}(p_{11}-d_{11})^{2}+2(p_{12}-d_{\textrm{s}})^{2}+(p_{22}-d_{22})^{2}\Big{)}.$
(19)

To justify the replacement, note that

$\displaystyle\begin{split}F(p_{11},p_{22},p_{12})=F_{\textrm{s}}(p_{11},p_{22},p_{12})+\frac{1}{4}(d_{12}-d_{21})^{2},\end{split}$
(20)

hence, $(p_{11},p_{22},p_{12})$ minimizes $F$ if and only if it minimizes
$F_{\textrm{s}}$. To obtain the minimizers, we minimize $F_{\textrm{s}}$ under
the condition $\mathbf{P}\in\mathcal{P}(\mathbf{m})$ using Lagrange
multipliers. The Lagrangian is thus given by

$\displaystyle\Lambda(p_{11},p_{22},p_{12},\lambda)=F_{\textrm{s}}(p_{11},p_{22},p_{12})+\lambda\left(p_{11}p_{22}-p_{12}^{2}+f^{2}\right).$
(21)

By setting the partial derivatives of $\Lambda$ with respect to $p_{11}$,
$p_{22}$, $p_{12}$ and $\lambda$ to zero, we find that the critical points of
$\Lambda$ have to satisfy

$\displaystyle p_{11}+\lambda p_{22}$ $\displaystyle=d_{11},$ (22a)
$\displaystyle\lambda p_{11}+p_{22}$ $\displaystyle=d_{22},$ (22b)
$\displaystyle(1-\lambda)p_{12}$ $\displaystyle=d_{\textrm{s}},$ (22c)
$\displaystyle p_{11}p_{22}-p_{12}^{2}$ $\displaystyle=-f^{2}.$ (22d)

This system can be solved analytically and the results are given by Prins et
al. [1, p. B942-B947] for the elliptic Monge-Ampère equation, with $-f^{2}$
replaced by $f^{2}$ in (22d). Unfortunately, the list of solutions is not
complete as for the case $d_{11}=-d_{22}$, two roots of (22) are missing. We
propose a different solution strategy here. First, two remarks are in place.
While minimizing $F_{\textrm{s}}$, the matrix $\mathbf{D}_{\mathrm{s}}$ and
the function value of $f$ are given and both $\mathbf{P}$ and $\lambda$ have
to be computed. Hence, we provide a classification in terms of
$\mathbf{D}_{\mathrm{s}}$ and the corresponding solutions of (22).
Furthermore, because the matrix $\mathbf{D}_{\mathrm{s}}$ is an approximation,
$\det(\mathbf{D}_{\mathrm{s}})\neq-f^{2}$ and in general
$\det(\mathbf{D}_{\mathrm{s}})\geq 0$ could possibly occur. We first write the
linear equations of (22) as

$\displaystyle\boldsymbol{\Lambda}\mathbf{p}=\mathbf{d},\quad\boldsymbol{\Lambda}=\begin{pmatrix}1&\lambda&0\\\
\lambda&1&0\\\
0&0&1-\lambda\end{pmatrix},\quad\mathbf{p}=\begin{pmatrix}p_{11}\\\ p_{22}\\\
p_{12}\end{pmatrix},\quad\mathbf{d}=\begin{pmatrix}d_{11}\\\ d_{22}\\\
d_{12}\end{pmatrix}.$ (23)

The vector $\mathbf{p}$ is uniquely determined when $\boldsymbol{\Lambda}$ is
regular, i.e., when
$0\neq\det(\boldsymbol{\Lambda})=(1-\lambda)^{2}(1+\lambda)$. We should
therefore distinguish between the cases $\lambda=1$, $\lambda=-1$ and
$\lambda\neq\pm 1$.

Although we should consider the cases $\lambda=1$, $\lambda=-1$ and
$\lambda\neq\pm 1$ separately, $\mathbf{d}$ and $f$ are given and $\lambda$
and $\mathbf{p}$ are to be calculated. Therefore we consider three cases based
on $\mathbf{D}_{\mathrm{s}}$, viz., Case 1: $d_{11}=d_{22}$ and
$d_{\mathrm{s}}=0$, Case 2: $d_{11}=-d_{22}$ and Case 3: all other
$\mathbf{D}_{\mathrm{s}}$. We consider $\mathbf{D}_{\mathrm{s}}=\mathbf{0}$ as
a special case of $d_{11}=d_{22}$ and $d_{\mathrm{s}}=0$.

We start with some general results, to be used in the subsequent derivations.
First, let $\mathrm{Tr}(\mathbf{A})$ denote the trace of a matrix
$\mathbf{A}$. Using (22) we find

$\displaystyle\delta_{\mathrm{s}}:=\det(\mathbf{D}_{\mathrm{s}})$
$\displaystyle=\lambda\mathrm{Tr}(\mathbf{P})^{2}-(\lambda-1)^{2}f^{2},$ (24a)
$\displaystyle\mathrm{Tr}(\mathbf{D}_{\mathrm{s}})$
$\displaystyle=(\lambda+1)\mathrm{Tr}(\mathbf{P}).$ (24b)

Solving the second equation for $\mathrm{Tr}(\mathbf{P})$ and subsequently
substituting it in the first equation yields

$\displaystyle
f^{2}(\lambda^{2}-1)^{2}+\delta_{\mathrm{s}}(\lambda+1)^{2}-\mathrm{Tr}(\mathbf{D}_{\mathrm{s}})^{2}\lambda=0.$
(25)

Next, we consider the roots of (25) and the corresponding solutions
$\mathbf{P}$.

Case 1: $d_{11}=d_{22}$ and $d_{\mathrm{s}}=0$, which we write as
$\mathbf{D}_{\mathrm{s}}=d\mathbf{I}$ with $d\in\mathbb{R}$. We will show that
this condition is equivalent with $\lambda=1$. So, let
$\mathbf{D}_{\mathrm{s}}=d\mathbf{I}$. We show that $\lambda=1$ by forcing a
contradiction, so, assume $\lambda\neq 1$. Then subtracting (22a) from (22b)
gives $p_{11}=p_{22}$ and by (22c) we have $p_{12}=0$. Substitution of
$p_{11}=p_{22}$ and $p_{12}=0$ in (22d) yields $p_{11}^{2}=-f^{2}<0$, being a
contradiction. Therefore $\lambda=1$. Conversely, substitution of $\lambda=1$
in $\boldsymbol{\Lambda}$ gives

$\displaystyle\boldsymbol{\Lambda}=\begin{pmatrix}1&1&0\\\ 1&1&0\\\
0&0&0\end{pmatrix}.$ (26)

In this case the null space of $\boldsymbol{\Lambda}$ is given by
$\mathcal{N}(\boldsymbol{\Lambda})=\langle\mathbf{v}_{1},\mathbf{v}_{2}\rangle$
with $\mathbf{v}_{1}=(1,-1,0)^{\text{T}}$ and
$\mathbf{v}_{2}=(0,0,1)^{\text{T}}$. So
$\boldsymbol{\Lambda}\mathbf{p}=\mathbf{d}$ only has a solution if
$\mathbf{d}$ lies in the column space of $\boldsymbol{\Lambda}$, i.e., if
$d_{11}=d_{22}$ and $d_{\mathrm{s}}=0$ or
$\mathbf{D}_{\mathrm{s}}=d\mathbf{I}$ with $d\in\mathbb{R}$. Henceforth we
have that $\lambda=1$ is equivalent with $\mathbf{D}_{\mathrm{s}}=d\mathbf{I}$
and thus $\lambda=1$ only occurs in Case 1. The general solution to
$\boldsymbol{\Lambda}\mathbf{p}=\mathbf{d}$ is now given by

$\displaystyle\mathbf{p}=(p,\,\,d-p,\,\,0)^{\text{T}}+\mu_{1}\mathbf{v}_{1}+\mu_{2}\mathbf{v}_{2},\qquad
p,\mu_{1},\mu_{2}\in\mathbb{R}.$ (27)

We aim to minimize $F_{s}$. Substitution of (27) in $F_{s}$ gives

$\displaystyle
F_{\textrm{s}}(p_{11},p_{22},p_{12})=\tfrac{1}{2}\Big{(}(p+\mu_{1}-d)^{2}+2\mu_{2}^{2}+(p+\mu_{1})^{2}\Big{)},$
(28)

thus showing $\mu_{2}=0$. Furthermore, the minimum is independent of the
choice for $\mu_{1}$ as can be seen by writing
$\mathbf{p}=(\tilde{p},\,\,d-\tilde{p},\,\,0)^{\text{T}}$ with
$\tilde{p}=p+\mu_{1}$. For simplicity we choose $\mu_{1}=0$. Subsequent
substitution of $\mathbf{p}$ into (22d) gives $p(d-p)=-f^{2}$. This second
order polynomial in $p$ has two real roots, viz.

$\displaystyle p=\frac{1}{2}\left(d\pm\sqrt{d^{2}+4f^{2}}\right).$ (29)

So in total we find the two solutions

$\displaystyle p_{11}=\frac{1}{2}\left(d\pm\sqrt{d^{2}+4f^{2}}\right),\qquad
p_{22}=d-p_{11},\qquad p_{12}=0.$ (30)

In case $d=0$, i.e., in case $\mathbf{D}_{\mathrm{s}}=\mathbf{0}$, the above
derivation still holds so we consider $\mathbf{D}_{\mathrm{s}}=\mathbf{0}$ an
instance of Case 1.

Case 2: $d_{11}=-d_{22}$, which we write as
$\mathbf{d}=(d,\,\,-d,\,\,d_{\mathrm{s}})^{\text{T}}$ with
$d,d_{\mathrm{s}}\in\mathbb{R}$. We have that
$\mathrm{Tr}(\mathbf{D}_{\mathrm{s}})=0$ and
$\delta_{\mathrm{s}}=-(d^{2}+d_{\mathrm{s}}^{2})$. For this case the fourth
order polynomial (25) can be written as

$\displaystyle(\lambda+1)^{2}\left((\lambda-1)^{2}+\frac{\delta_{\mathrm{s}}}{f^{2}}\right)=0.$
(31)

It follows that we have the three unique roots, $\lambda=-1$ (with
multiplicity 2) and $\lambda=1\pm\sqrt{|\delta_{\mathrm{s}}|}/f$.

  * •

In case $\lambda=-1$ we have

$\displaystyle\boldsymbol{\Lambda}=\begin{pmatrix}1&-1&0\\\ -1&1&0\\\
0&0&2\end{pmatrix},$ (32)

and the corresponding null space
$\mathcal{N}(\boldsymbol{\Lambda})=\langle\mathbf{v}_{3}\rangle$ with
$\mathbf{v}_{3}=(1,1,0)^{\text{T}}$. For $\mathbf{p}$ to be a solution to
$\boldsymbol{\Lambda}\mathbf{p}=\mathbf{d}$ we require $\mathbf{d}$ to be in
the column space of $\boldsymbol{\Lambda}$. It follows that
$\mathbf{d}=(d,\,\,-d,\,\,d_{\mathrm{s}})^{\text{T}}$,
$d,d_{\mathrm{s}}\in\mathbb{R}$. Henceforth $\lambda=-1$ only occurs for Case
2.

The general solution to $\boldsymbol{\Lambda}\mathbf{p}=\mathbf{d}$ is
therefore given by

$\displaystyle\mathbf{p}=(p,\,\,p-d,\,\,\tfrac{1}{2}d_{\mathrm{s}})^{\text{T}}+\mu_{3}\mathbf{v}_{3},\qquad
p,\mu_{3}\in\mathbb{R}.$ (33)

Writing
$\mathbf{p}=(\tilde{p},\,\,\tilde{p}-d,\,\,\tfrac{1}{2}d_{\mathrm{s}})^{\text{T}}$
with $\tilde{p}=p+\mu_{3}$ shows that the actual solution $\mathbf{p}$ does
not change by choosing $\mu_{3}$, so we simply choose $\mu_{3}=0$. By (22d) it
follows that

$\displaystyle p(p-d)-\frac{1}{4}d_{\textrm{s}}^{2}+f^{2}=0.$ (34)

Consequently, solving for $p$ we find that for
$|\delta_{\mathrm{s}}|-4f^{2}\geq 0$ we have

$\displaystyle
p_{11}=\tfrac{1}{2}\left(d\pm\sqrt{|\delta_{\mathrm{s}}|-4f^{2}}\right),\qquad
p_{22}=p_{11}-d,\qquad p_{12}=\tfrac{1}{2}d_{\textrm{s}}.$ (35)

When $|\delta_{\mathrm{s}}|-4f^{2}<0$ the solution $\mathbf{p}$ is complex.
Because we are only interested in real-valued solutions, we do not consider
$\lambda=-1$ when $|\delta_{\mathrm{s}}|-4f^{2}<0$.

  * •

In the case $\lambda=1\pm\sqrt{|\delta_{\mathrm{s}}|}/f$, the matrix
$\boldsymbol{\Lambda}^{-1}$ is uniquely defined (see (37) for an explicit
expression) and by $\mathbf{p}=\boldsymbol{\Lambda}^{-1}\mathbf{d}$ we obtain

$\displaystyle p_{11}=\mp\frac{df}{\sqrt{|\delta_{\mathrm{s}}|}},\qquad
p_{22}=-p_{11},\qquad
p_{12}=\mp\frac{d_{\textrm{s}}f}{\sqrt{|\delta_{\mathrm{s}}|}}.$ (36)

The solutions (36) are new with respect to those found by Prins et al. [1] and
are not specific to the hyperbolic Monge-Ampère equation.

Case 3: All other $\mathbf{D}_{\mathrm{s}}$, i.e., both
$\mathbf{D}_{\mathrm{s}}\neq d\mathbf{I}$ and
$\mathbf{d}_{\mathrm{s}}\neq(d,\,\,-d,\,\,d_{\mathrm{s}})^{\text{T}}$ for all
$d,d_{\mathrm{s}}\in\mathbb{R}$. By Case 1 we have $\lambda\neq 1$ and by Case
2 we have $\lambda\neq-1$. Therefore
$\det\boldsymbol{\Lambda}=(1-\lambda)^{2}(1+\lambda)\neq 0$. Consequently
$\boldsymbol{\Lambda}$ is invertible and its inverse is given by

$\displaystyle\boldsymbol{\Lambda}^{-1}=\frac{1}{1-\lambda^{2}}\begin{pmatrix}1&-\lambda&0\\\
-\lambda&1&0\\\ 0&0&1+\lambda\end{pmatrix}.$ (37)

The values for $\lambda$ are obtained by solving (25). The roots of this
fourth order polynomial can be determined analytically using Ferrari’s method
[17, p. 22] and are given in [1, p. B945]. For $\mathbf{p}$ we subsequently
find $\mathbf{p}=\boldsymbol{\Lambda}^{-1}\mathbf{d}$, or more explicitly

$\displaystyle p_{11}=\frac{\lambda d_{22}-d_{11}}{\lambda^{2}-1},\qquad
p_{22}=\frac{\lambda d_{11}-d_{22}}{\lambda^{2}-1},\qquad
p_{12}=\frac{d_{\textrm{s}}}{1-\lambda}.$ (38)

###  2.3 $J_{\mathrm{B}}$-Optimization

In [2, p. 131-133] a projection method (PM) has been proposed for the
minimization of (9b). As our numerical results will show, this method proves
insufficient for some examples. Therefore we developed two improved methods,
viz., the segmented projection method (SPM) and the segmented arc length
method (SALM). Before we introduce the boundary methods, we first introduce
some notation.

Let $\mathbf{x}_{ij}=((x_{1})_{i},(x_{2})_{j})\in\mathcal{X}$ be the grid
points of a Cartesian grid with $i=1,\dots,N_{x_{1}}$ and $j=1,\dots
N_{x_{2}}$ denoting the first and second coordinate, respectively. Let
$\mathbf{x}_{l}$ be the grid points restricted to $\partial\mathcal{X}$ for
$l=1,\dots N$. We index $\mathbf{x}_{l}=1,\dots N$ in the clockwise direction
such that $\mathbf{x}_{1}=\mathbf{x}_{1,1}$, i.e. the first point on
$\partial\mathcal{X}$ equals the point $\mathbf{x}_{ij}$ with $i=j=1$. We
approximate $\mathbf{m}_{ij}\approx\mathbf{m}(\mathbf{x}_{ij})$,
$\mathbf{m}_{l}\approx\mathbf{m}(\mathbf{x}_{l})$ and similarly for the other
variables.

The main idea behind SPM and SALM is to partition the boundaries of the source
and target domains in segments. We then uniquely enforce one source segment to
be mapped to one target segment. We then distribute $\mathbf{b}_{l}$,
corresponding to $\mathbf{m}_{l}$ by either a projection (SPM) or by a ratio
of arc lengths (SALM).

Let the boundary segments of $\mathcal{X}$ be the curves
$\Gamma_{k}^{\mathcal{X}}\subset\partial\mathcal{X}$ such that for
$N_{\Gamma}$ boundary segments we have
$\cup_{k=1}^{N_{\Gamma}}\Gamma_{k}^{\mathcal{X}}=\partial\mathcal{X}$. We
denote $\Gamma_{N_{\Gamma}+1}^{\mathcal{X}}=\Gamma_{1}^{\mathcal{X}}$ and
assume the intersections
$\Gamma_{k_{1}}^{\mathcal{X}}\cap\Gamma_{k_{2}}^{\mathcal{X}}$ contain
precisely one element if $k_{2}=k_{1}+1$ and no elements otherwise.
Furthermore, we require each $\Gamma_{k}^{\mathcal{X}}$ to be parametrizable.
We assume similar properties for $\Gamma^{\mathcal{Y}}_{k}$. We aim to map
each boundary segment of $\partial\mathcal{X}$ to a boundary segment of
$\partial\mathcal{Y}$, hence we enforce
$\mathbf{m}(\Gamma_{k}^{\mathcal{X}})=\Gamma_{k}^{\mathcal{Y}}$ for
$k=1,\dots,N_{\Gamma}$, from which it follows that

$\displaystyle\mathbf{m}(\partial\mathcal{X})=\mathbf{m}(\cup_{k=1}^{N_{\Gamma}}\Gamma^{\mathcal{X}}_{k})=\cup_{k=1}^{N_{\Gamma}}\mathbf{m}(\Gamma^{\mathcal{X}}_{k})=\cup_{k=1}^{N_{\Gamma}}\Gamma^{\mathcal{Y}}_{k}=\partial\mathcal{Y},$
(39)

which is the required transport boundary condition.

Figure 1 shows a part of $\partial\mathcal{Y}$ and (parts of) three boundary
segments. In the following we fix $k$ and for brevity drop the subscript in
$\Gamma_{k}^{\mathcal{X}}$ and $\Gamma_{k}^{\mathcal{Y}}$.

Let $\mathbf{y}_{i}$ for $i=1,\dots,N_{\mathrm{b}}$ be a discretization of the
boundary segment $\Gamma^{\mathcal{Y}}$ such that for a given counter
clockwise parametrization
$\mathbf{y}(s):[0,1]\rightarrow\Gamma^{\mathcal{Y}}$, we have
$\mathbf{y}_{1}=\mathbf{y}(0)$ and
$\mathbf{y}_{N_{\mathrm{b}}}=\mathbf{y}(1)$. We choose to parametrize
$\partial\mathcal{X}$ and $\partial\mathcal{Y}$ in opposite directions because
$\mathbf{m}$ reverts the direction if it is a solution to the hyperbolic
Monge-Ampère.

$\Gamma^{\mathcal{Y}}_{k}$$\mathbf{y}_{1}$$\mathbf{y}_{i}$$\mathbf{y}_{i+1}$$\mathbf{y}_{N_{\textrm{b}}}$$\Gamma^{\mathcal{Y}}_{k-1}$$\Gamma^{\mathcal{Y}}_{k+1}$$\tau_{i}$
Figure 1: Schematic overview of the discretization of $\partial\mathcal{Y}$.

Projection method. We briefly explain PM as introduced by [2, p. 131-133]. We
perform the following for each approximation $\mathbf{m}_{l}$ individually.
Let $N_{\Gamma}=1$, i.e., we consider the whole boundary as one boundary
segment. Furthermore, let $\mathbf{y}_{N_{\textrm{b}}+1}=\mathbf{y}_{1}$, we
connect adjacent points $\mathbf{y}_{i}$ and $\mathbf{y}_{i+1}$ by straight
line segments. The projection of $\mathbf{m}_{l}$ onto the line connecting
$\mathbf{y}_{i}$ and $\mathbf{y}_{i+1}$ is given by

$\displaystyle\mathbf{m}^{\mathrm{P}}_{i}(t_{i})$
$\displaystyle=\mathbf{y}_{i}+t_{i}(\mathbf{y}_{i+1}-\mathbf{y}_{i}),$ (40a)
$\displaystyle t_{i}$
$\displaystyle=\frac{(\mathbf{m}_{l}-\mathbf{y}_{i})\boldsymbol{\cdot}(\mathbf{y}_{i+1}-\mathbf{y}_{i})}{|\mathbf{y}_{i+1}-\mathbf{y}_{i}|^{2}}.$
(40b)

As only $0\leq t_{i}\leq 1$ corresponds to a point on the line segment between
$\mathbf{y}_{i}$ and $\mathbf{y}_{i+1}$, we limit $t_{i}$ according to
$\hat{t}_{i}=\min(1,\max(0,t_{i}))$. Among all possible line segments, we
choose $\mathbf{b}_{l}$ corresponding to $\mathbf{m}_{l}$, such that the
distance is smallest, i.e.,

$\displaystyle i_{\mathrm{min}}$
$\displaystyle=\operatorname*{argmin}_{i}\\{|\mathbf{m}^{\mathrm{P}}_{i}(\hat{t}_{i})-\mathbf{m}_{l}|\\},$
(41a) $\displaystyle\mathbf{b}_{l}$
$\displaystyle=\mathbf{m}^{\mathrm{P}}_{i_{\mathrm{min}}}(\hat{t}_{\mathrm{min}}).$
(41b)

Segmented projection method. Let
$\Gamma_{k}^{\mathcal{X}}\subseteq\partial\mathcal{X}$ and
$\Gamma_{k}^{\mathcal{Y}}\subseteq\partial\mathcal{Y}$ be boundary segments of
the source and target, respectively. For this method, we apply PM to the
individual boundary segments instead of the whole boundary at once.
Furthermore, we set

$\displaystyle\mathbf{b}(\Gamma_{k}^{\mathcal{X}}\cap\Gamma_{k+1}^{\mathcal{X}})=\Gamma_{k}^{\mathcal{Y}}\cap\Gamma_{k+1}^{\mathcal{Y}},\qquad
1\leq k\leq N_{\Gamma}$ (42)

meaning, we map the end points of the source segments to the end points of the
corresponding target segments. In practice, these end points are the corners
of the source and target domains.

Segmented arc length method. The core idea of this method is as follows: if
$\mathbf{m}(\Gamma^{\mathcal{X}})=\Gamma^{\mathcal{Y}}$, then the arc length
of the curve $\mathbf{m}(\Gamma^{\mathcal{X}})$ should be equal to the arc
length of the curve $\Gamma^{\mathcal{Y}}$. Numerically we approximate this
condition by approximating the arc length of both $\Gamma^{\mathcal{Y}}$ and
the distance between the points
$\\{\mathbf{m}(\mathbf{x}_{l})\mid\mathbf{x}_{l}\in\Gamma^{\mathcal{X}}\\}$.

We start with the arc length of the curve $\Gamma^{\mathcal{Y}}$. We
approximate the arc length between $\mathbf{y}_{i}$ and $\mathbf{y}_{i+1}$
along $\Gamma^{\mathcal{Y}}$ by the length of the line segment connecting
$\mathbf{y}_{i}$ and $\mathbf{y}_{i+1}$. We denote the approximation by

$\displaystyle\tau_{i}=|\mathbf{y}_{i+1}-\mathbf{y}_{i}|,\quad i=1,\dots
N_{\mathrm{b}}-1.$ (43)

The approximate cumulative arc length between $\mathbf{y}_{1}$ and
$\mathbf{y}_{i}$ in the direction of increasing $s$ is then given by

$\displaystyle t_{i}=\sum_{j=1}^{i-1}\tau_{j},\quad i=1,\dots,N_{\mathrm{b}}.$
(44)

The total arc length from $\mathbf{y}_{1}$ to $\mathbf{y}_{N_{\mathrm{b}}}$ is
then approximated by $L=t_{N_{\textrm{b}}}$. We use the cumulative arc lengths
to introduce a piece-wise linear interpolation $\mathbf{b}_{\mathrm{int}}$
approximating $\mathbf{y}(s)$, viz.

$\displaystyle\mathbf{b}_{\mathrm{int}}(t)=\mathbf{y}_{i}+\frac{t-t_{i}}{t_{i+1}-t_{i}}(\mathbf{y}_{i+1}-\mathbf{y}_{i}),\quad
t_{i}\leq t\leq t_{i+1},$ (45)

where the scalar factor is a scaled coordinate between $\mathbf{y}_{i}$ and
$\mathbf{y}_{i+1}$. Note that by construction $\mathbf{b}_{\mathrm{int}}$
satisfies

$\displaystyle\mathbf{b}_{\mathrm{int}}(t_{i})=\mathbf{y}_{i},\quad
i=1,\dots,N_{\mathrm{b}},$ (46)

and is an approximation of $\partial\mathcal{Y}$.

Next we consider the points
$\\{\mathbf{m}(\mathbf{x}_{l})\mid\mathbf{x}_{l}\in\Gamma^{\mathcal{X}}\\}$.
Let $N_{\mathrm{m}}$ be the number of grid points on $\Gamma^{\mathcal{X}}$
such that $\mathbf{x}_{l}\in\Gamma^{\mathcal{X}}$ for
$l=1,\dots,N_{\mathrm{m}}$. Furthermore, let

$\displaystyle\sigma_{l}=|\mathbf{m}_{l+1}-\mathbf{m}_{l}|,\quad
l,\dots,N_{\textrm{m}}-1,$ (47)

be an approximation of the arc length from $\mathbf{m}(\mathbf{x}_{l})$ to
$\mathbf{m}(\mathbf{x}_{l+1})$ along $\partial\mathcal{Y}$. This again
introduces a cumulative arc length and a total arc length, respectively, given
by

$\displaystyle s_{l}=\sum_{j=1}^{l-1}\sigma_{j},\quad
l=1,\dots,N_{\mathrm{m}},\qquad\tilde{L}=s_{N_{\mathrm{m}}}.$ (48)

Because $\mathbf{m}_{l}$ is an approximation and $\Gamma^{\mathcal{Y}}$ is
approximated by straight line segments, $\tilde{L}\neq L$ in general. Hence,
$s_{l}\neq L$ may occur such that the end points of $\Gamma^{\mathcal{X}}$ may
not be mapped to the end points of $\Gamma^{\mathcal{Y}}$. We fix this by
letting

$\displaystyle\tilde{s}_{l}=\frac{L}{\tilde{L}}s_{l}.$ (49)

It follows that $\mathbf{b}_{\mathrm{int}}(\tilde{s}_{l})$ forms a proper
approximation for $\mathbf{m}_{l}$ restricted to $\Gamma^{\mathcal{Y}}$, viz.

$\displaystyle\mathbf{b}_{l}=\mathbf{b}_{\mathrm{int}}(\tilde{s}_{l}),\quad
l=1,\dots,N_{\mathrm{m}}.$ (50)

###  2.4 Grid shock correction

Using the methods outlined above, it is possible that the approximation
$\mathbf{m}^{n}$ of $\mathbf{m}$ contains crossing grid lines, also known as
grid shocks [18]. This phenomenon is shown in Figure 2 for an example we
discuss in Section 3.2, with grid parameters $N_{x_{1}}=N_{x_{2}}=321$ after
$n=15,000$ iterations. Though the solution on the left may look visually
correct, the grid shock, as seen on the right, prevents proper numerical
convergence of our algorithm.

Figure 2: Example of grid shock. The global numerical approximation is shown
on the left, and a zoomed-in version on the right.

$\partial\mathcal{X}$$\mathbf{x}_{i-1,1}$$\mathbf{x}_{i,1}$$\mathbf{x}_{i+1,1}$$\mathbf{x}_{i-1,2}$$\mathbf{x}_{i,2}$$\mathbf{x}_{i+1,2}$$\mathbf{x}_{i,3}$

$\partial\mathcal{Y}$$\mathbf{y}_{k}$$\mathbf{y}_{k+1}$$\mathbf{y}_{k+2}$$\mathbf{m}_{i,3}$$\mathbf{m}_{i+1,2}$$\mathbf{m}_{i,2}$$\mathbf{m}_{i-1,2}$$\mathbf{m}_{i,1}$$\mathbf{b}_{i,1}$$\mathbf{m}_{i+1,1}$$\mathbf{m}_{i-1,1}$

Figure 3: Schematic overview of stencil used for detecting grid shocks. On the
right $\mathbf{m}_{i-1,2}$ is closer to $\mathbf{b}_{i,1}$ than
$\mathbf{m}_{i,1}$ is, so a grid shock occurs.

To resolve grid shocks, consider a point
$\mathbf{x}_{ij}\in\partial\mathcal{X}$ as shown for $j=1$ in Figure 3 on the
left, and the corresponding image $\mathbf{m}_{ij}$ shown on the right. If
both $\mathbf{m}_{ij}$ and $\mathbf{b}_{ij}$ are exact, then
$|\mathbf{m}_{ij}-\mathbf{b}_{ij}|=0$ and
$|\mathbf{m}_{kl}-\mathbf{b}_{ij}|>0$ for all $(k,l)\neq(i,j)$. Because both
$\mathbf{m}$ and $\mathbf{b}$ are approximated,
$|\mathbf{m}_{ij}-\mathbf{b}_{ij}|\neq 0$ in general. To detect grid shocks,
we compute the distance $|\mathbf{m}_{ij}-\mathbf{b}_{ij}|$ for all $(k,l)$
such that $|(k,l)^{\text{T}}-(i,j)^{\text{T}}|\leq 2$. If the minimum distance
if found for $(k,l)\neq(i,j)$, then we assume a grid shock orrcurs and we
recompute $\mathbf{m}$. We do so by making $\alpha$ in (9c) dependent on the
coordinate, i.e., $\alpha=\alpha(\mathbf{x})$ and
$\alpha_{ij}=\alpha(\mathbf{x}_{ij})$ and subsequently reduce $\alpha_{ij}$ on
the boundary which puts more emphasis on the minimization of
$|\mathbf{m}_{ij}-\mathbf{b}_{ij}|$.

Introducing the $\mathbf{x}$-dependency, the coefficients $\alpha$ and
$1-\alpha$ in (9c) formally have to be moved inside the integrals of
$J_{\textrm{I}}$ and $J_{\textrm{B}}$. After doing so, we compute the first
variation of
$J(\mathbf{m},\mathbf{P},\mathbf{b},\alpha)=\frac{1}{2}\iint_{\mathcal{X}}\alpha\|\mathrm{\mathbf{D}}\mathbf{m}-\mathbf{P}\|^{2}\mathop{}\\!\mathrm{d}\mathbf{x}+\frac{1}{2}\oint_{\partial\mathcal{X}}(1-\alpha)|\mathbf{m}-\mathbf{b}|^{2}\mathop{}\\!\mathrm{d}s$
w.r.t. $\mathbf{m}$ and apply the fundamental lemma of calculus of variations.
Consequently we obtain

$\displaystyle\nabla\alpha\boldsymbol{\cdot}\nabla m_{1}+\alpha\Delta m_{1}$
$\displaystyle=\mathbf{p}_{1}\boldsymbol{\cdot}\nabla\alpha+\alpha\nabla\boldsymbol{\cdot}\mathbf{p}_{1},$
$\displaystyle\mathbf{x}\in\mathcal{X},$ (51a)
$\displaystyle(1-\alpha)m_{1}+\alpha\nabla m_{1}$
$\displaystyle=(1-\alpha)b_{1}+\alpha\mathbf{p}_{1}\cdot\hat{\mathbf{n}},$
$\displaystyle\mathbf{x}\in\partial\mathcal{X},$ (51b)

for the first component of $\mathbf{m}$. If $\alpha$ is constant in the
interior, we have $\nabla\alpha=\mathbf{0}$ in the interior and equations (51)
reduce to (12). By analogy, we have (13) for the second component of
$\mathbf{m}$. Let $\alpha_{1}\in(0,1)$. We set
$\alpha(\mathbf{x}_{ij})=\alpha_{1}$ for $\mathbf{x}_{ij}$ in the interior of
$\mathcal{X}$. For a boundary point $\mathbf{x}_{ij}\in\partial\mathcal{X}$ we
instead set

$\displaystyle\alpha_{ij}=\begin{cases}\alpha_{2},&\quad\text{if}\quad\min\limits_{(k,l)\in\mathcal{I}(i,j)}|\mathbf{m}_{kl}-\mathbf{b}_{ij}|<|\mathbf{m}_{ij}-\mathbf{b}_{ij}|,\\\
\alpha_{1}&\quad\text{otherwise},\end{cases}$ (52)

where
$\mathcal{I}(i,j)=\\{(k,l)\mid\mathbf{x}_{kl}\in\mathrm{int}(\mathcal{X}),\,|(k,l)^{\text{T}}-(i,j)^{\text{T}}|\leq
2\\}$ is the space over which we minimize, $\mathrm{int}(\mathcal{X})$ the
interior of $\mathcal{X}$ and $\alpha_{2}\in(0,\alpha_{1})$ a constant. We
choose a distance of 2 and the values $\alpha_{2}=0.005$ and $\alpha_{1}=0.2$
since these have proven to work well in practice.

If in the $n^{\text{th}}$ iteration we obtain
$\alpha(\mathbf{x})\not\equiv\alpha_{1}$, we solve (12) and (13) for a second
time with the updated $\alpha$ to perform a correction.

Recall, we use a finite difference method for inverting the Poisson equations
(12) and (13), yielding a system of equations. This system of equations
depends on $\alpha$. Without grid shock correction, a LU-factorization can be
calculated once and used for each subsequent iteration making the inversion of
the system efficient. In case $\alpha_{ij}\neq\alpha_{1}$ for any $(i,j)$, the
same LU-factorization can no longer be used due to the component $\alpha\nabla
m_{k}\boldsymbol{\cdot}\hat{\mathbf{n}}$ in the Robin boundary conditions and
a new LU-factorization has to be calculated for the iteration.

##  3 Numerical results

In this section we present numerical results for five examples. For each
example we know the exact solution and compare the numerical methods. We
choose $\mathcal{X}=[x_{1}^{m},x_{1}^{M}]\times[x_{2}^{m},x_{2}^{M}]$, a
rectangle which may vary per case. For each example we choose $\mathcal{Y}$
such that it has a unique feature to it. We measure the residual

$\displaystyle\epsilon_{r}$
$\displaystyle=\Big{|}D_{x_{1}}[m_{1}]_{ij}D_{x_{2}}[m_{2}]_{ij}-D_{x_{1}}[m_{2}]_{ij}D_{x_{2}}[m_{1}]_{ij}+f^{2}(\mathbf{x}_{ij},\mathbf{m}_{ij})\Big{|}_{\infty},$
(53)

with $D_{x_{1}}$ and $D_{x_{2}}$ standard second-order (central in interior
and one-sided on boundary) finite difference operators for the first-order
derivatives with respect to $x_{1}$ and $x_{2}$, respectively. Furthermore, we
measure the global discretization errors $\epsilon_{u}$, $\epsilon_{m_{1}}$
and $\epsilon_{m_{2}}$ defined by

$\displaystyle\begin{split}\epsilon_{u}&=\big{|}\big{(}u_{ij}-u_{11}\big{)}-\big{(}u(\mathbf{x}_{ij})-u(\mathbf{x}_{11})\big{)}\big{|}_{\infty},\\\
\epsilon_{m_{1}}&=\big{|}(m_{1})_{ij}-m_{1}(\mathbf{x}_{ij})\big{|}_{\infty},\\\
\epsilon_{m_{2}}&=\big{|}(m_{2})_{ij}-m_{2}(\mathbf{x}_{ij})\big{|}_{\infty},\\\
\end{split}$ (54)

where the terms $u_{11}$ and $u(\mathbf{x}_{11})$ are introduced due to the
nonuniqueness of $u$ given $\mathbf{m}$; see the discussion following
equations (15). Any fixed grid point could be used, here $\mathbf{x}_{11}$ is
chosen. The choice for the $\infty$-norm is arbitrary in the sense that any
standard norm would give similar results. However, the $\infty$-norm is more
sensitive to differences in the local errors than, for example, the standard
2-norm. Starting the least-squares algorithm requires an initial guess
$\mathbf{m}^{0}$, so we introduce
$\widetilde{\mathcal{Y}}=[y_{1}^{m},y_{1}^{M}]\times[y_{2}^{m},y_{2}^{M}]$,
the smallest bounding box of $\mathcal{Y}$. We then choose
$\partial\mathbf{m}^{0}(\mathcal{X})=\partial\widetilde{\mathcal{Y}}$, such
that $\mathbf{m}^{0}_{ij}$ is equidistantly distributed, i.e.,
$\mathbf{m}^{0}_{ij}$ is the result a bilinear uniform interpolation of the
bounding box of $\mathcal{Y}$ with
$\det(\mathrm{\mathbf{D}}\mathbf{m}^{0})<0$. The initial guess then reads

$\displaystyle(m^{0}_{1})_{ij}$
$\displaystyle=\frac{(x_{1})_{ij}-x_{1}^{m}}{x_{1}^{M}-x_{1}^{m}}y_{1}^{M}+\frac{x_{1}^{M}-(x_{1})_{ij}}{x_{1}^{M}-x_{1}^{m}}y_{1}^{m},$
(55a) $\displaystyle(m^{0}_{2})_{ij}$
$\displaystyle=\frac{x_{2}^{M}-(x_{2})_{ij}}{x_{2}^{M}-x_{2}^{m}}y_{2}^{M}+\frac{(x_{2})_{ij}-x_{2}^{m}}{x_{2}^{M}-x_{2}^{m}}y_{2}^{m}.$
(55b)

The initial guess is a (discretized) solution of the hyperbolic Monge-Ampère
equation with
$f^{2}=\mathrm{area}(\widetilde{\mathcal{Y}})/\mathrm{area}(\mathcal{X})$,
$\mathbf{m}=\nabla u$ and $u=\tfrac{1}{2}(x^{2}-y^{2})f$. Three more such
initial guesses exist, viz., $u=\tfrac{1}{2}(y^{2}-x^{2})f$ and $u=\pm xyf$.

We segmentate the source boundary in segments, clockwise, according to

$\displaystyle\begin{split}\Gamma_{1}^{\mathcal{X}}&=\\{x_{1}^{m}\\}\times[x_{2}^{m},x_{2}^{M}],\qquad\Gamma_{2}^{\mathcal{X}}=[x_{1}^{m},x_{1}^{M}]\times\\{x_{2}^{M}\\},\\\
\Gamma_{3}^{\mathcal{X}}&=\\{x_{1}^{M}\\}\times[x_{2}^{m},x_{2}^{M}],\qquad\Gamma_{4}^{\mathcal{X}}=[x_{1}^{m},x_{1}^{M}]\times\\{x_{2}^{m}\\},\end{split}$
(56)

and we write

$\displaystyle\Gamma^{\mathcal{Y}}_{k}=\left\\{\Gamma^{\mathcal{Y}}_{k}(s)\,\Big{|}\,s\in[0,1]\right\\},$
(57)

for $k=1,\dots,4$. Furthermore, we apply grid shock correction only for
iteration step $n\geq 100$, as the distance between the boundary of the
initial guess and the boundary of the target may be large for small $n$.

Lastly, we stop the iteration (10) based on the update of $\mathbf{m}^{n}$,
i.e., based on

$\displaystyle\Delta m^{n}=|\mathbf{m}^{n}-\mathbf{m}^{n-1}|,$ (58)

instead of, the already introduced measures, $J_{\textrm{I}}$ and
$J_{\textrm{B}}$. This is because the values for $J_{\textrm{I}}$ and
$J_{\textrm{B}}$ may stagnate over the iterations while $\Delta m^{n}$ is
still changing. Conversely, if $\Delta m^{n}$ has stagnated, then so have the
functionals $J_{\textrm{I}}$ and $J_{\textrm{B}}$. We stop the iterative
process when $\Delta m^{n}$ reaches floating-point precision.

###  3.1 Annulus segment

For the first example we consider $\mathcal{X}=[0,1]\times[-1/2,1/2]$,
$\partial\mathcal{Y}=\cup_{k=1}^{4}\Gamma_{k}^{\mathcal{Y}}$ with

$\displaystyle\Gamma_{1}^{\mathcal{Y}}(s)$
$\displaystyle=(\cos(\tfrac{1}{2}-s),\,\sin(\tfrac{1}{2}-s)),$ (59a)
$\displaystyle\Gamma_{2}^{\mathcal{Y}}(s)$
$\displaystyle=(e^{s}\cos(\tfrac{1}{2}),\,-e^{s}\sin(\tfrac{1}{2})),$ (59b)
$\displaystyle\Gamma_{3}^{\mathcal{Y}}(s)$
$\displaystyle=(e\cos(\tfrac{1}{2}-s),\,e\sin(s-\tfrac{1}{2})),$ (59c)
$\displaystyle\Gamma_{4}^{\mathcal{Y}}(s)$
$\displaystyle=(e^{1-s}\cos(\tfrac{1}{2}),\,e^{1-s}\sin(\tfrac{1}{2})),$ (59d)

Figure 4:  The exact mapping $\mathbf{m}$ (left) and solution $u$ (right) on a
$21\times 21$ grid.

as shown together with the exact mapping on a $21\times 21$ grid in Figure 4
on the left. We choose $\Gamma^{\mathcal{Y}}_{k}=\nabla
u(\Gamma^{\mathcal{X}}_{k})$ for all examples. Furthermore, this choice of
$\Gamma_{k}^{\mathcal{Y}}$ implies that for SALM and SPM the corners of
$\partial\mathcal{X}$ are mapped to the corners of $\partial\mathcal{Y}$. Let
$f^{2}(x_{1},x_{2})=e^{2x_{1}}$. The solution is then given by

$\displaystyle u(x_{1},x_{2})=e^{x_{1}}\cos(x_{2}),$ (60)

which is symmetric in $x_{2}=0$ as can be seen in Figure 4. Unless specified
otherwise, we take $N_{\mathrm{b}}=10^{4}$ and for each target segment
$\Gamma^{\mathcal{Y}}_{k}$, with $k=1,\dots,4$, we construct
$\mathbf{y}_{i}=\Gamma^{\mathcal{Y}}_{k}(s_{i})$ with
$s_{i}=(i-1)/(N_{\mathrm{b}}-1)$ and $i=1,\dots,N_{\mathrm{b}}$. The results
for PM, SPM and SALM are shown in Figure 5 for varying grid configurations
with $N_{x_{1}}=N_{x_{2}}$. The three figures clearly show second-order
convergence of the relevant errors and residual, which is in accordance with
the discretization error of the finite difference approximations used. In
terms of $\epsilon_{u}$, PM and SPM ($3\cdot 10^{-6}$) slightly outperform
SALM ($6\cdot 10^{-6}$) for $N_{x_{1}}=N_{x_{2}}=473$, though the difference
is small.

Figure 5:  Global error and the residual for PM (left), SPM (middle) and SALM
(right).

Figure 6 shows the $J$-errors over the iterations on a grid of
$N_{x_{1}}=N_{x_{2}}=473$. For PM and SPM we obtained $J_{\mathrm{I}}\approx
2\cdot 10^{-12}$, $J_{\textrm{B}}\approx 3\cdot 10^{-13}$ in approximately
60,000 iterations. SALM gave $J_{\mathrm{I}}\approx 4\cdot 10^{-12}$,
$J_{\textrm{B}}\approx 7\cdot 10^{-13}$ in 40,000, iterations. SALM
consistently requires less iterations as seen on the left in Figure 7, where
the number of required iterations ($n_{\mathrm{max}}$) for various
$N_{x_{1}}=N_{x_{2}}$ is shown.

Convergence of $u$ with respect to $N_{\mathrm{b}}$ is shown on the right of
Figure 7 for $N_{x_{1}}=N_{x_{2}}=473$. Two observations are in place. First,
for increasing $N_{\mathrm{b}}$, the error $\epsilon_{u}$ reaches an
asymptotic value (dashed black line). This phenomenon is to be expected and
occurs when the discretization errors in $\mathbf{m}_{ij}$, $\mathbf{P}_{ij}$
and $u_{ij}$, and the finite differences $\mathbf{D}_{ij}$ become dominant,
i.e., when the discretization errors due to the choice of $N_{x_{1}}$ and
$N_{x_{2}}$ dominate the errors due to discretizing the boundary. Secondly, in
the regime prior to the asymptote, the discretization error in $u$ due to the
boundary discretization is second-order accurate for all three boundary
methods.

Figure 6:  $J$-errors and $\Delta m$ over the iterations for PM (left), SPM
(middle) and SALM (right).

Finally, the computational cost per iteration for SALM is lowest, second comes
SPM and third PM. The projection methods calculate
$C\max(N_{x_{1}},N_{x_{2}})N_{\mathrm{b}}$ projections and performs
$C\max(N_{x_{1}},N_{x_{2}})$ searches over $N_{\mathrm{b}}$ points each, with
$C\in\mathbb{N}_{+}$. Similarly, SALM performs a linear interpolation of
$C\max(N_{x_{1}},N_{x_{2}})$ points over $N_{\mathrm{b}}-1$ segments.
Therefore, one would expect the average time per iteration to scale linearly
in $N_{x},N_{y}$ and $N_{\mathrm{b}}$ for PM and SPM when
$N_{x_{1}}=N_{x_{2}}$, and $N_{\mathrm{b}}$ is fixed. For SALM, a linear
relation is also expected, with a possible asymptote when either the
computational load due to $\max(N_{x_{1}},N_{x_{2}})$ or $N_{\mathrm{b}}$
dominates. This is also shown in Figure 8, where, from left to right,
$N_{x_{1}}=N_{x_{2}}=501$ is fixed while $N_{\textrm{b}}$ varies,
$N_{\textrm{b}}=37$ is fixed and $N_{x_{1}}=N_{x_{2}}$ varies, and lastly,
$N_{\textrm{b}}=10007$ is fixed and $N_{x_{1}}=N_{x_{2}}$ varies.
Additionally, it is observed that SALM, on average, significantly outperforms
the projection methods. Lastly, SPM is approximately four times faster than PM
because SPM projects one source segment on a target segment (four times)
instead of the whole source boundary on the whole target boundary.

Figure 7:  The total number of iterations needed for convergence (left) and
the influence of the boundary discretization on $\epsilon_{u}$ (right).

Figure 8:  Timing results in seconds for the boundary procedures. On the left,
$N_{x_{1}}=N_{x_{2}}=501$ is fixed and in the middle and on the right
$N_{\textrm{b}}=37$ and $N_{\textrm{b}}=10007$, respectively.

###  3.2 Deformed square

As a second example we consider a mapping and surface with no symmetries, viz.
$\mathcal{X}=[0,1]\times[-1/2,1/2]$,
$\partial\mathcal{Y}=\cup_{k=1}^{4}\Gamma_{k}^{\mathcal{Y}}$ with

$\displaystyle\Gamma_{1}^{\mathcal{Y}}(s)$
$\displaystyle=(s-\tfrac{1}{2},\,-s+\tfrac{1}{2}),$ (61a)
$\displaystyle\Gamma_{2}^{\mathcal{Y}}(s)$
$\displaystyle=(\tfrac{1}{3}s^{3}+s^{2}+\tfrac{1}{2},\,s-\tfrac{1}{2}),$ (61b)
$\displaystyle\Gamma_{3}^{\mathcal{Y}}(s)$
$\displaystyle=(-s+\tfrac{11}{6},\,s+\tfrac{1}{2}),$ (61c)
$\displaystyle\Gamma_{4}^{\mathcal{Y}}(s)$
$\displaystyle=(-\tfrac{1}{3}s^{3}+2s^{2}-3s+\tfrac{5}{6},\,\tfrac{3}{2}-s),$
(61d)

as shown on the left of Figure 9. Let $f^{2}(x_{1},x_{2})=(x_{1}+1)^{2}$, then
the exact solution is given by

$\displaystyle u(x_{1},x_{2})$
$\displaystyle=\frac{1}{12}x_{1}^{4}+\frac{1}{3}x_{1}^{3}+x_{1}x_{2}-\frac{1}{2}x_{2}^{2},$
(62)

which is shown on the right of Figure 9. In Figure 10 errors and residuals are
shown for SPM (left) and SALM (right) for varying grid configurations with
$N_{x_{1}}=N_{x_{2}}$. Both figures show second-order convergence which is in
accordance with the discretization error of the finite difference
approximations used.

Figure 9:  The exact mapping $\mathbf{m}$ (left) and solution $u$ (right) on a
$21\times 21$ grid.

Figure 10:  Convergence of the global error and the residual for SPM (left)
and SALM (right).

Figure 11 shows the behaviour of $J_{\textrm{I}}$, $J_{\textrm{B}}$ and
$\Delta m$ for SPM and SALM. For SPM, on the left $\Delta m$ exhibit
oscillations starting at 100 iterations. This is due to the grid shock
correction, enabled in the 100${}^{\text{th}}$ iteration. As it turns out,
this is one example where grid shocks occur using SPM. Without the grid shock
correction, $\Delta m$ would still go to computer precision but the grid shock
(as visualized on the right of Figure 2) would remain and subsequently
$J_{\textrm{I}}$, $J_{\textrm{B}}$ and the errors $\epsilon_{u}$,
$\epsilon_{m_{1}}$ and $\epsilon_{m_{2}}$ and the residual $\epsilon_{r}$
would be three orders of magnitude higher.

Figure 11:  Comparison of $J_{\textrm{I}}$, $J_{\textrm{B}}$ and $\Delta m$
for SPM (left) and SALM (right) for a grid of $295\times 295$.

Lastly we discuss the results for PM. In Figure 12 the errors $J_{\textrm{I}}$
and $J_{\textrm{B}}$ are shown, both with grid shock correction (left) and
without (right). Clearly the example with grid shock correction does not
converge, the method actually oscillates between intermediate solutions. One
may be tempted to think that without grid shock correction the method does
work, as $\Delta m$ goes to machine precision, but this is not the case as
shown in Figure 13. The two leftmost figures show the mapping for a $29\times
29$ grid after 50,000 iterations. Clearly, neither of the methods work as
intended as there are gaps between the mesh spanned by $\mathbf{m}_{l}$ and
$\partial\mathcal{Y}$, i.e., the transport boundary condition has not been
satisfied. The reason why the algorithm with PM does not converge is that the
projection of $\mathbf{m}$ onto $\partial\mathcal{Y}$ does not distribute
$\mathbf{b}$ well. In particular, no points $\mathbf{b}_{l}$ near
$(0.5,-0.5)\in\partial\mathcal{Y}$ are obtained, as can be seen in Figure 13
on the right, where the blue circles represent $\mathbf{m}_{l}$, the red
squares $\mathbf{b}_{l}$ and the thin black lines connect $\mathbf{m}_{l}$ to
$\mathbf{b}_{l}$ for $l=1,\dots,N$.

Figure 12:  The errors $J_{\textrm{I}}$, $J_{\textrm{B}}$ and the update
$\Delta m$ for PM with shock correction (left) and without (right).

Figure 13:  The mapping after 50,000 iterations for PM with shock correction
(left) and without (middle) and the accompanying projection of
$\mathbf{m}_{l}$ onto $\partial\mathcal{Y}$ for construction of
$\mathbf{b}_{l}$ without shock correction on the right.

###  3.3 Inward fold

For this example we consider the target as illustrated in Figure 14, for the
exact solution on a $61\times 61$ grid. On the right a zoomed-in version of
the target is shown. In the figure we have marked two points, one by a solid
circle, and one by an asterisk. The former is a point for which the boundary
of the target is not differentiable, while for the latter it is. We will come
back to this.

Figure 14: The target domain and mapping, with a zoomed-in version of the
Inward fold example.

The example shown in Figure 14 corresponds to
$\mathcal{X}=[0,1]\times[-1/2,1/2]$ and
$\partial\mathcal{Y}=\cup_{k=1}^{4}\Gamma_{k}^{\mathcal{Y}}$ with

$\displaystyle\Gamma_{1}^{\mathcal{Y}}(s)$
$\displaystyle=(0,-\tfrac{1}{4}s^{4}+\tfrac{1}{2}s^{3}-\tfrac{3}{8}s^{2}-\tfrac{7}{8}s+\tfrac{31}{64}),$
(63a) $\displaystyle\Gamma_{2}^{\mathcal{Y}}(s)$
$\displaystyle=(\tfrac{9}{8}s-\tfrac{1}{2}s^{3},-\tfrac{1}{4}s^{4}+\tfrac{3}{8}s^{2}-\tfrac{33}{64}),$
(63b) $\displaystyle\Gamma_{3}^{\mathcal{Y}}(s)$
$\displaystyle=(-s^{3}+\tfrac{3}{2}s^{2}+\tfrac{1}{4}s+\tfrac{5}{8},-\tfrac{1}{4}s^{4}+\tfrac{1}{2}s^{3}+\tfrac{9}{8}s^{2}-\tfrac{3}{8}s-\tfrac{25}{64}),$
(63c) $\displaystyle\Gamma_{4}^{\mathcal{Y}}(s)$
$\displaystyle=(-\tfrac{1}{2}s^{3}+\tfrac{3}{2}s^{2}-\tfrac{19}{8}s+\tfrac{11}{8},-\tfrac{1}{4}s^{4}+s^{3}-\tfrac{9}{8}s^{2}+\tfrac{1}{4}s+\tfrac{39}{64}).$
(63d)

Furthermore we have

$\displaystyle f^{2}(x_{1},x_{2})$
$\displaystyle=x_{1}^{6}+3x_{1}^{4}x_{2}^{2}+3x_{1}^{2}x_{2}(x_{2}^{3}-2)+(1+x_{2}^{3})^{2},$
(64a) $\displaystyle u(x_{1},x_{2})$
$\displaystyle=\frac{x_{1}^{2}}{2}-\frac{x_{1}^{4}x_{2}}{4}-\frac{x_{2}^{2}}{2}+\frac{x_{1}^{2}x_{2}^{3}}{2}-\frac{x_{2}^{5}}{20}.$
(64b)

Taking derivatives of $u$ yields the mapping, i.e.,

$\displaystyle m_{1}(x_{1},x_{2})$
$\displaystyle=u_{x_{1}}(x_{1},x_{2})=x_{1}-x_{1}^{3}x_{2}+x_{1}x_{2}^{3},$
(65a) $\displaystyle m_{2}(x_{1},x_{2})$
$\displaystyle=u_{x_{2}}(x_{1},x_{2})=-\frac{x_{1}^{4}}{4}-x_{2}+\frac{3x_{1}^{2}x_{2}^{2}}{2}-\frac{x_{2}^{4}}{4}.$
(65b)

A straightforward calculation shows that $\mathbf{m}(1,1/2)=(5/8,-25/64)$
which corresponds to the solid circle in Figure 14. Henceforth, $\mathbf{m}$
is not differentiable in the point $(1,1/2)$ as it is the image of a
nondifferentiable (corner) point in $\mathcal{X}$ under a continuously
differentiable map. The point depicted by the asterisk originates from the
source boundary segment
$\Gamma^{\mathcal{X}}_{2}=[0,1]\times\\{\tfrac{1}{2}\\}$. Let $m_{1}$ and
$m_{2}$ along the boundary be parametrized by $s$. Then in the point indicated
by the asterisk, both
$\tfrac{\mathop{}\\!\mathrm{d}m_{1}(s)}{\mathop{}\\!\mathrm{d}s}$ and
$\tfrac{\mathop{}\\!\mathrm{d}m_{2}(s)}{\mathop{}\\!\mathrm{d}s}$ change sign.
Henceforth, the location of the asterisk can be obtained by solving
$\tfrac{\mathop{}\\!\mathrm{d}m_{1}(s)}{\mathop{}\\!\mathrm{d}s}=\tfrac{\mathop{}\\!\mathrm{d}m_{2}(s)}{\mathop{}\\!\mathrm{d}s}=0$,
which is equivalent to

$\displaystyle\left.\frac{\partial m_{1}}{\partial
x_{1}}\right|_{x_{2}=1/2}=0,\qquad\left.\frac{\partial m_{2}}{\partial
x_{1}}\right|_{x_{2}=\tfrac{1}{2}}=0,\qquad 0\leq x_{1}\leq 1.$ (66)

Indeed, doing so one uniquely obtains $x_{1}=\sqrt{3}/2$ such that
$\mathbf{m}(\sqrt{3}/2,1/2)=(3\sqrt{3}/8,-3/8)$, which corresponds to the
point indicated by an asterisk in Figure 14. Furthermore, smoothness of the
boundary in said point is implied.

PM and SPM do not yield converging numerical approximations. Figure 15 shows a
zoomed-in version of two numerical solutions for $161\times 161$ grids. The
sharp inward fold seems to be the culprit for the boundary method, as is also
seen in Figure 16, which shows the projection of $\mathbf{m}_{l}$ onto
$\partial\mathcal{Y}$. The figure clearly shows that the method does not pick
$\mathbf{b}_{l}$ deep within the fold and, consequently, the optimization for
$\mathbf{m}$ does not produce a mapping with such a sharp fold.

Because SALM does force points $\mathbf{b}_{l}$ to be located along the whole
boundary, naturally points will end up in the fold. This can be seen in Figure
17, where on the left the first iteration of applying SALM to the result of
SPM is shown. The $50^{\text{th}}$ iteration of SALM is shown in the middle,
showing $\mathbf{m}(\partial\mathcal{X})$ being positioned in the fold.
Continuation using SALM yields similar results to using SALM starting from the
default initial guess. SALM shows approximately second-order convergence, as
graphed on the right of Figure 17, when using (55) as initial guess. SALM does
not show any visual distortions, in contrast to PM and SPM.

For the remaining results we will not discuss PM, as it performs, at best, as
good as SPM while being more computationally expensive.

Figure 15: Zoomed-in results for PM (left) and SPM (right) on a $161\times
161$ grid.

Figure 16: Projection step after convergence for PM (left) and for SPM (right)
on a $161\times 161$ grid.

Figure 17: The first (left) and 50th (middle) iteration of continuation by
SALM after SPM has converged on a $161\times 161$ grid and the convergence of
the errors and residual for SALM starting from a uniform initial guess
(right).

###  3.4 Annulus

Figure 18: The target and the exact mapping on the left, and the solution
surface on the right, both shown on a $51\times 51$ grid.

For this example we consider the target given in Figure 18, where the central
part near $\mathbf{m}=(0,0)$ is not part of $\mathcal{Y}$. Let
$\mathcal{X}=[0,2\pi]\times[-1/2,1/2]$,
$\partial\mathcal{Y}=\cup_{k=1}^{4}\Gamma_{k}^{\mathcal{Y}}$ with

$\displaystyle\Gamma_{1}^{\mathcal{Y}}(s)$
$\displaystyle=(0,\,e^{s-\frac{1}{2}}),$ (67a)
$\displaystyle\Gamma_{2}^{\mathcal{Y}}(s)$ $\displaystyle=\sqrt{e}(-\sin(2\pi
s),\,\cos(2\pi s)),$ (67b) $\displaystyle\Gamma_{3}^{\mathcal{Y}}(s)$
$\displaystyle=(0,\,e^{\frac{1}{2}-s}),$ (67c)
$\displaystyle\Gamma_{4}^{\mathcal{Y}}(s)$
$\displaystyle=\frac{1}{\sqrt{e}}(\sin(2\pi s),\,\cos(2\pi s)),$ (67d)

and $f^{2}(x_{1},x_{2})=e^{2x_{2}}$, such that the exact solution is given by

$\displaystyle u(x_{1},x_{2})$ $\displaystyle=e^{x_{2}}\cos(x_{1}),$ (68)

as shown on the right of Figure 18. Observe that
$\mathbf{m}|_{\partial\mathcal{X}}$ is not bijective, as
$\Gamma^{\mathcal{Y}}_{1}=\Gamma^{\mathcal{Y}}_{3}$. Nevertheless, we
introduce both $\Gamma^{\mathcal{Y}}_{1}$ and $\Gamma^{\mathcal{Y}}_{3}$ as
the orientation, i.e., the parametrization of the segments, matters for SALM.

Figure 19 shows results for SPM. On the left the mapping after the algorithm
has converged for a grid with $N_{x_{1}}=115$ and $N_{x_{2}}=19$. The grid
parameters are chosen such that $h_{x_{1}}\approx h_{x_{2}}$ as
$N_{x_{1}}/N_{x_{2}}\approx(x_{1}^{M}-x_{1}^{m})/(x_{2}^{M}-x_{2}^{m})=2\pi$.
Although the figure on the right clearly shows $J_{\textrm{I}}$ and
$J_{\textrm{B}}$ have converged, and that $\Delta m$ reached computer
precision, the algorithm does not yield a correct solution, as it does not
satisfy the transport boundary condition because there are points
$\mathbf{m}_{ij}$ which lie outside $\mathcal{Y}$, nor does it solve the
hyperbolic Monge-Ampère equation as is shown by the residual $\epsilon_{r}$ in
Figure 20 on the left.

Figure 19:  The numerical mapping $\mathbf{m}$ after convergence (left) and
the history of $J_{\textrm{I}}$, $J_{\textrm{B}}$ and $\Delta m$ for SPM.

Figure 20:  Residual for SPM (left) and SALM (right). No convergence for SPM,
and second-order convergence for SALM is observed.

The reason why SPM does not produce accurate solutions is easiest demonstrated
by visualizing a few iterations. To this end, consider the boundary routine
for the first, third and tenth iteration as shown in Figure 21. We focus on
one segment of the mapping of the boundary, i.e, $\mathbf{m}_{ij}$ with
$\mathbf{x}_{ij}\in\Gamma^{\mathcal{X}}_{4}=[0,2\pi]\times\\{-\tfrac{1}{2}\\}$,
which corresponds to $m_{2}(0,-\tfrac{1}{2})=\exp(-1/2)\approx 0.6$ for the
initial guess in Figure 21. For the exact solution, $\Gamma^{\mathcal{X}}_{4}$
needs to be mapped to the entire inner circle of the target, i.e.,
$\Gamma^{\mathcal{Y}}_{4}$. As shown for the first iteration,
$\Gamma^{\mathcal{X}}_{4}$ is mapped to only part of
$\Gamma^{\mathcal{Y}}_{4}$, viz., the accompanying $\mathbf{b}_{l}$ lies on
the northern part of $\Gamma^{\mathcal{Y}}_{4}$ (the inner circle). In
subsequent iterations, shown in the middle and on the right in Figure 21,
$\Gamma^{\mathcal{X}}_{4}$ will again not be mapped to the whole of
$\Gamma^{\mathcal{Y}}_{4}$, as the distance to the northern part of the inner
circle remains minimal. This process continues indefinitely.

For SALM such accumulation of $\mathbf{b}_{l}$ does not occur, as by
construction, $\mathbf{b}_{l}$ is distributed over the boundary segments. The
results of the first, second and third iteration of the
$\mathbf{b}$-minimization are shown in Figure 22. Clearly, SALM does not
suffer from the same flaws as SPM. As such, the convergence is expected to
behave as for the other examples, which is confirmed by the results shown in
Figure 20 on the right.

Figure 21: From left to right, the first, third and tenth iteration for SPM.

Figure 22: From left to right, the first, second and third iteration for SALM.

###  3.5 Gradient dependent problem

Lastly we consider an example with $f$ dependent on the gradient of the
solution, i.e., $f=f(x_{1},x_{2},\nabla u)$, viz.

$\displaystyle
f^{2}(\mathbf{x},\mathbf{m})=3x_{2}^{2}-m_{1}\sin(x_{1})-\frac{1}{4}m_{2}^{2}.$
(69)

We consider the domain $\mathcal{X}=[-1,1]\times[1,3/2]$ and
$\partial\mathcal{Y}=\cup_{k=1}^{4}\Gamma_{k}^{\mathcal{Y}}$ with

$\displaystyle\Gamma_{1}^{\mathcal{Y}}(s)$
$\displaystyle=(\tfrac{\sin(1)}{4}s^{2}+\sin(1)s+\sin(1),\,\cos(1)s+2\cos(1)),$
(70a) $\displaystyle\Gamma_{2}^{\mathcal{Y}}(s)$
$\displaystyle=(-\tfrac{9}{4}\sin(2s-1),\,3\cos(2s-1)),$ (70b)
$\displaystyle\Gamma_{3}^{\mathcal{Y}}(s)$
$\displaystyle=(-\tfrac{\sin(1)}{4}s^{2}+\tfrac{3}{2}s-\tfrac{9\sin(1)}{4},\,-\cos(1)s+3\cos(1)),$
(70c) $\displaystyle\Gamma_{4}^{\mathcal{Y}}(s)$
$\displaystyle=(\sin(2s-1),\,2\cos(2s-1)).$ (70d)

The exact solution is given by

$\displaystyle u(x_{1},x_{2})$ $\displaystyle=x_{2}^{2}\cos(x_{1}),$ (71)

and is, together with the mapping and target domain, shown in Figure 23.

By construction of the algorithm, little effort is required for $f$ to be
dependent on the mapping $\mathbf{m}$. The difference being that during the
$n^{\text{th}}$ iteration, $f^{2}(\mathbf{x}_{ij},\mathbf{m}^{n}_{ij})$ has to
be evaluated instead of $f^{2}(\mathbf{x}_{ij})$ in the optimization of
$\mathbf{P}$. The results for SPM and SALM with $N_{x_{1}}=N_{x_{2}}$ are
given in Figure 24, showing second-order convergence for both methods. In this
case grid shock correction is needed for SPM to ensure proper convergence.

Figure 23:  The target and the exact mapping on the left, and the solution
surface on the right, both shown on a $51\times 51$ grid.

Figure 24:  Convergence of SPM (left) and SALM (right).

##  4 Conclusion

We have introduced a least-squares solver for the hyperbolic Monge-Ampère
equation with transport boundary condition. The algorithm, originally
introduced by Prins et al. [1] for the elliptic Monge-Ampère equation, has
been improved to encompass a more complete description of the roots for the
$P$-optimization. Furthermore, we introduced two new boundary methods. All
three boundary methods, if convergent, show second-order convergence of the
residual and the global discretization errors as function of the mesh size,
and also second-order convergence as function of the number of boundary
points. Of the three boundary methods, the segmented arc length method is both
the only method to converge for all examples and is most computationally
efficient, both in terms of computation time per iteration, as in total number
of iterations required.

As far as the authors are aware, the least-squares method paired with the
segmented arc length method for the boundary, is the first method to solve the
hyperbolic Monge-Ampère equation with transport boundary conditions.

## References

  * [1] C. R. Prins, R. Beltman, J. H. M. ten Thije Boonkkamp, W. L. IJzerman, and T. W. Tukker.  A least-squares method for Optimal Transport using the Monge–Ampère equation.  SIAM Journal on Scientific Computing, 37(6):B937–B961, jan 2015\.
  * [2] C. R. Prins.  Inverse Methods for Illumination Optics.  PhD thesis, Eindhoven University of Technology, 2014.
  * [3] H. Ries and A. Rabl.  Edge-ray principle of nonimaging optics.  Journal of the Optical Society of America A, 11(10):2627, October 1994.
  * [4] C. Villani.  Optimal Transport: Old and New.  Number 338 in Grundlehren der Mathematischen Wissenschaften. Springer, Berlin, 2009.
  * [5] F. Santambrogio.  Optimal Transport for Applied Mathematicians.  Springer International Publishing, 2015.
  * [6] Y. Brenier.  Décomposition polaire et réarrangement monotone des champs de vecteurs. (French).  R. Acad. Sci. Paris Sér. I Math., 305(19):805–808, 1987.
  * [7] D. V. Tunitskii.  On the global solubility of the Monge–Ampère hyperbolic equations.  Izvestiya: Mathematics, 61(5):1069–1111, oct 1997.
  * [8] J. Hong.  The global smooth solutions of Cauchy problems for hyperbolic equation of Monge-Ampère type.  Nonlinear Analysis: Theory, Methods & Applications, 24(12):1649–1663, June 1995.
  * [9] B. S. Westcott and F. Brickell.  Computation of reflector surfaces for two-variable beam shaping in the hyperbolic case.  Journal of Physics A: Mathematical and General, 9(4):611–625, April 1976.
  * [10] M. W. M. C. Bertens, E. M. T. Vugts, M. J. H. Anthonissen, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman.  Numerical methods for the hyperbolic Monge–Ampère equation based on the method of characteristics.  Partial Differential Equations and Applications, 3(4):52, Jul 2022\.
  * [11] L. B. Romijn, J. H. M. ten Thije Boonkkamp, M. J. H. Anthonissen, and W. L. IJzerman.  An iterative least-squares method for generated Jacobian equations in freeform optical design.  SIAM Journal on Scientific Computing, 43(2):B298–B322, January 2021.
  * [12] N. K. Yadav, J. H. M. ten Thije Boonkkamp, and W. L IJzerman.  A Monge–Ampère problem with non-quadratic cost function to compute freeform lens surfaces.  Journal of Scientific Computing, 80(1):475–499, 2019.
  * [13] R. Courant and D. Hilbert.  Methods of Mathematical Physics, volume 1.  John Wiley & Sons, Ltd, 1989.
  * [14] S Abdallah.  Numerical solutions for the pressure Poisson equation with Neumann boundary conditions using a non-staggered grid, I.  J. Comput. Phys., 70(1):182–192, May 1987.
  * [15] B. D. Froese.  A numerical method for the elliptic Monge–Ampère equation with transport boundary conditions.  SIAM Journal on Scientific Computing, 34(3):A1432–A1459, January 2012.
  * [16] L. B. Romijn.  Generated Jacobian Equations in Freeform Optical Design: Mathematical Theory and Numerics.  PhD thesis, Eindhoven University of Technology, 2021.
  * [17] Jean-Pierre Tignol.  Galois’ theory of algebraic equations.  World Scientific Publishing, Singapore, Singapore, April 2001.
  * [18] J. Cordova and T. Barth.  Grid generation for general 2-D regions using hyperbolic equations.  In 26th Aerospace Sciences Meeting. American Institute of Aeronautics and Astronautics, January 1988.