An efficient approach for the numerical solution of the Monge–Ampère equation

Applied Numerical Mathematics 61 (2011) 298–307

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

An efficient approach for the numerical solution of the Monge–Ampère equation Mohamed M. Sulman ∗ , J.F. Williams, Robert D. Russell Department of Mathematics, Simon Fraser University, Burnaby, British Columbia, V5A 1S6 Canada

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 1 May 2009 Received in revised form 6 March 2010 Accepted 13 October 2010 Available online 23 October 2010

In this paper, we present a new method to compute the numerical solution of the elliptic Monge–Ampère equation. This method is based on solving a parabolic Monge– Ampère equation for the steady state solution. We study the problem of global existence, uniqueness, and convergence of the solution of the fully nonlinear parabolic PDE to the unique solution of the elliptic Monge–Ampère equation. Some numerical experiments are presented to show the convergence and the regularity of the numerical solution. © 2010 IMACS. Published by Elsevier B.V. All rights reserved.

Keywords: Monge–Ampère equation Global existence Convergence

1. Introduction In this work, we study the numerical solution of the Monge–Ampère equation

det ∇ 2 Ψ = g (ξ , ∇Ψ ),

ξ ∈ Ωc

(1)

where ∇ 2 Ψ is the Hessian matrix of Ψ , Ωc is a bounded domain of Rd , and g is a given positive function. As seen below, for grid generation we are interested in the numerical solution of a particular form of (1) written as

det ∇ 2 Ψ (ξ ) =

f 1 (ξ ) f 2 (∇Ψ )

,

ξ ∈ Ωc , ∇Ψ (ξ ) ∈ Ω ⊂ Rd ,

(2)

where Ψ is convex, ∇Ψ : Ωc → Ω with Ω is a bounded domain in Rd , and f 1 , f 2 are given smooth positive functions satisfying the compatibility condition





f 1 (ξ ) dξ = Ωc

f 2 (x) dx.

(3)

Ω

The numerical solution of the elliptic Monge–Ampère equation (2) has been of considerable interest in the contexts of a wide range of applications in different areas of science and engineering, such as image processing, meteorology, fluid dynamics and adaptive grid generation. The existence, uniqueness and regularity of the solution of (2) have been studied extensively, for example, see [10–12,7,5,4,6,3,19,2]. Aleksandrov [1] showed that there exists a unique generalized solution of (2) with Dirichlet boundary conditions in the class of convex functions. Brenier [7] showed the existence and uniqueness of the solution of (2) under some suitable conditions. In [9] Caffarelli constructed interior W 2, p and C 2,α estimates for the solution of (2). Dean and Glowinski [14,15] proposed an augmented Lagrangian and a least-squares approach for the

*

Corresponding author. E-mail address: [email protected] (M.M. Sulman).

0168-9274/$30.00 © 2010 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2010.10.006

M.M. Sulman et al. / Applied Numerical Mathematics 61 (2011) 298–307

299

numerical solution of (2) with Dirichlet boundary condition. In [21], a 2D numerical solution of (2) with periodic boundary conditions is computed by Newton’s algorithm. In this paper, we present a new approach of computing a numerical solution of the elliptic Monge–Ampère equation (2) based on solving an appropriate parabolic Monge–Ampère equation (PMAE). We then develop a theory of global existence, uniqueness and convergence of the numerical solution of the PMAE to the solution of elliptic equation (2) defined on square or cubic domains. The paper is organized as follows. In Section 2, we review the derivation of the elliptic Monge–Ampère equation from the theory of the L 2 Monge–Kantorovich problem (MKP). In Section 3 we describe the PMA adaptive mesh method based on solving a parabolic Monge–Ampère equation. In Section 4, we present some theoretical results for global existence, uniqueness and convergence of the solution of the parabolic Monge–Ampère equation to the solution of the elliptic Monge– Ampère equation. In Section 5, we present some numerical experiments to illustrate the global existence and convergence results. Finally, we conclude with some remarks and discussion in Section 6. 2. The L 2 Monge–Kantorovich problem The L 2 Monge–Kantorovich problem is stated as follows: Given ρ0 and ρ1 defined in Rd , d  1, which are two positive density functions with equal masses of a given material, find the mapping x = φ(ξ ) : Rd −→ Rd that transfers the density ρ0 to ρ1 and minimizes the transport cost



C (φ) =

  φ(ξ ) − ξ 2 ρ0 (ξ ) dξ .

(4)

Rd

Mathematically, the mapping φ realizes the transfer of





ρ1 (x) dx, ∀Ω ⊂ Rd ,

ρ0 (ξ ) dξ = φ −1 (Ω)

ρ0 to ρ1 if (5)

Ω

where | · | is the standard Euclidean norm defined on the space Rd . If φ is a smooth one-to-one map, then by the change of variables x = φ(ξ ), (5) leads to





J φ(ξ )





ρ1 φ(ξ ) = ρ0 (ξ ),

(6)

where J(φ(ξ )) is the determinant of the Jacobian matrix of φ . The existence and the characterization of the solution for the L 2 Monge–Kantorovich problem is addressed in the following theorem (see [20,7]). Theorem 2.1 (Brenier). There is a unique optimal solution φ = φˆ of the L 2 Monge–Kantorovich problem. Moreover, φˆ is characterized as the gradient of a convex potential Ψ , i.e.,

φˆ = ∇Ψ.

(7)

Substituting (7) into (6) shows that Ψ is a solution of the elliptic Monge–Ampère equation (MAE)





ρ1 ∇Ψ (ξ ) det ∇ 2 Ψ (ξ ) = ρ0 (ξ ).

(8)

3. The PMA adaptive mesh method For a standard adaptive mesh method, the adaptive mesh is obtained as the image of a mesh in a computational domain

Ωc ⊂ Rd under a coordinate transformation φ mapping Ωc into a physical domain Ω ⊂ Rd , i.e.,

φ : Ωc → Ω,

x = φ(ξ ),

ξ ∈ Ωc , x ∈ Ω.

For the PMA adaptive mesh method the transformation φ is determined as the solution of the L problem for which the density functions are defined as

ρ0 (ξ ) := 1,

ρ1 (x) := ρ (x), ξ ∈ Ωc , x ∈ Ω,

(9) 2

Monge–Kantorovich

(10)

where ρ (x) is a user defined monitor function. In this case, the transformation φ for the mesh adaptation is determined as the gradient of the solution of the Monge– Ampère equation (8) (one equation in one scalar potential) with some suitable boundary conditions. However, the Monge– Ampère equation (8) is fully nonlinear, so developing a fast and robust numerical method is quite challenging. Avoiding solving (8) directly, Brenier and Benamou [4] introduce a fluid dynamic framework to compute a numerical solution of the

300

M.M. Sulman et al. / Applied Numerical Mathematics 61 (2011) 298–307

L 2 Monge–Kantorovich problem. Our alternative approach for computing a numerical solution for (8) is based on finding a steady-state solution of the logarithmic parabolic Monge–Ampère equation

  ∂Ψ ρ1 (∇Ψ (ξ )) det ∇ 2 Ψ , = log ∂t ρ0 (ξ )

(11)

together with some suitable temporal initial and spatial boundary conditions. The logarithmic term in (11) is introduced in order to preserve the convexity of the solution. For an approach which uses a different form of PMA, see [8]. If Ψ∞ is the solution of (11) as t → ∞, then Ψ∞ solves the Monge–Ampère equation (8). Thus, the solution φˆ of the L 2 Monge–Kantorovich problem is determined by taking the spatial gradient of Ψ∞ with respect to ξ , i.e.,

ˆ ) = ∇Ψ∞ (ξ ), φ(ξ

∀ξ ∈ Ωc ,

(12)

ˆ ), ξ ∈ Ωc , x ∈ Ω. and the adaptive mesh in Ω is given from x = φ(ξ We consider the case where both the computational and the physical domains Ωc and Ω are squares (in 2D) or cubes (in 3D). Eq. (11) is solved using the Neumann boundary condition

∇Ψ · n = ξ · n,

for ξ ∈ ∂Ωc ,

(13)

and the initial condition

1

Ψ (ξ , 0) = Ψ0 (ξ ) = ξ · ξ T , 2

(14)

where n is the outward normal to ∂Ωc , the boundary of Ωc . The condition (13) enforces the boundary points of Ωc to be mapped to the boundary points of Ω . 4. Existence, uniqueness and convergence results In this section, we study global existence, uniqueness and convergence of the solution of the parabolic Monge–Ampère equation (11) to the solution of the Monge–Ampère equation (8). Eq. (8) is a particular case of the general form of the Monge–Ampère equation

det ∇ 2 Ψ = g (ξ , ∇Ψ ).

(15)

Existence of globally smooth convex solutions of the parabolic Monge–Ampère equation of the form (11) and its convergence to the solution of (15) have been investigated in [24,23,22] under suitable regularity and structure conditions on Ωc , Ω , and g. In what follows we describe these results. Let Ωc and Ω be smooth strictly convex domains in Rd , d  2, and consider the initial boundary value problem

⎧   ∂Ψ ⎪ ⎪ = log det ∇ 2 Ψ − log g (ξ , ∇Ψ ) in Ωc × [0, T ), ⎨ ∂t ∇Ψ (Ωc ) = Ω, ⎪ ⎪ ⎩ Ψ (·, 0) = Ψ0 in Ωc ,

(16)

on a time interval [0, T ), T > 0. The boundary condition ∇Ψ (Ωc ) = Ω means that ∇Ψ maps Ωc onto Ω . Assume that Ψ0 : Ωc → R, is a smooth strictly convex function, ∇Ψ0 (Ωc ) = Ω ( = 0th compatibility condition). The boundary condition ∇Ψ (Ωc ) = Ω is equivalent to (for more details see [24,27])

h(∇Ψ ) = 0,

on ∂Ωc

(17)

for a smooth strictly convex function Ψ , where h : Rd → R is a smooth strictly concave function such that h|∂Ωc = 0 and |∇ h| = 1 on ∂Ωc . The compatibility conditions are defined as (cf. [24])



d dt

m

 

h(∇Ψ )

t =0

= 0,

m ∈ N.

(18)

Moreover, assume that

g (ξ , ∇Ψ ) =

g 1 (ξ ) g 2 (∇Ψ )

(19)

for smooth positive functions g 1 and g 2 satisfying





g 1 (ξ ) dξ = Ωc

g 2 (x) dx, Ω

(20)

M.M. Sulman et al. / Applied Numerical Mathematics 61 (2011) 298–307

301

Fig. 1. Modifying a square domain Ω to a smooth strictly convex domain Ω n .

for all smooth convex domains Ωc and Ω . In the context of adaptive grid generation, we notice that the function g is of the form (19). In particular, for the PMA adaptive mesh method, we have

g (ξ , ∇Ψ ) =

ρ0 (ξ ) , ρ1 (∇Ψ )

(21)

moreover, ρ0 and ρ1 , defined as in (10), satisfy the condition (20), which is simply equivalent to the equidistribution principle (6) (for example, see [8]). The existence, uniqueness, and convergence of the solution of (16) to the solution of (15) follow from the following theorem. Theorem 4.1. Assume that Ωc , Ω , g and Ψ0 are defined as above. Furthermore, assume that



0  det ∇ 2 Ψ0 − log g (ξ , ∇Ψ0 ) in Ωc , 1st compatibility condition

(22)

on ∂Ωc

is satisfied. Then there exists a smooth strictly convex function, Ψ : Ω × (0, ∞) → R, that solves (16) for all times, i.e., T = ∞, and as t → ∞, Ψ converges smoothly to the solution Ψ∞ of



det ∇ 2 Ψ∞ (ξ ) = g (ξ , ∇Ψ∞ ),

in Ωc

(23)

∇Ψ∞ (Ωc ) = Ω. Furthermore, if (18) is satisfied for all m, then Ψ is smooth in [0, T ) and it converges exponentially to Ψ∞ .

Theorem 4.1 follows directly from Theorem A.1 of [22] and [24]. We remark that for the PMA adaptive mesh method, it is straightforward to show that the conditions (22) and (18) of Theorem 4.1 follow directly from the initial and boundary conditions (14) and (13). In this case, the function h is defined as

h(∇Ψ ) = ∇Ψ · n − ξ · n,

on ∂Ωc .

(24)

In what follows we extend the results of Theorem 4.1 to include the case where the domains Ωc and Ω are squares or cubes (for simplicity, we take Ωc = Ω = [0, 1]d ⊂ Rd , d = 2, 3, . . .). 

Definition 4.1. Let {n } be a sequence of small positive real numbers such that n → 0 as n → ∞. Define Ωc n and Ω n as smooth strictly convex domains obtained by modifying the boundaries of Ωc and Ω (see Fig. 1), where

n = sup |x1 − x2 |: x1 ∈ ∂Ω, x2 ∈ ∂Ω n



x1 ,x2

and | · | is the standard Euclidean norm defined in Rd . 

For each n > 0, define Ψ n : Ωc n → R to be a solution to the initial-boundary value problem (11), (13), (14) and assume that ρ0 and ρ1 satisfy the condition (20) and g has the form (21). Then from Theorem 4.1, as t → ∞, Ψ n converges smoothly to the unique solution of

    ρ1n ∇Ψ n (ξ ) det ∇ 2 Ψ n = ρ0n (ξ ) in Ωcn ,   ∇Ψ n Ωcn = Ω n .

(25)



Remark 4.1. We can extend the definition of Ψ n on Ωc n to Rd by using the Legendre transform convention

 Ψ n (ξ ) = sup ξ · x − Φ n (x) , x∈Ω n

ξ ∈ Rd ,

302

M.M. Sulman et al. / Applied Numerical Mathematics 61 (2011) 298–307

where Φ n is defined as



Φ n (x) = sup ξ · x − Ψ n (ξ ) , ξ ∈Ωcn

x ∈ Rd .

Theorem 4.2. Let 0 ∈ Ω, 0 ∈ Ω n ∀n and Ψ n , Ψ be defined as in (25), (8) such that Ψ n (0) = 0 = Ψ (0). Suppose that as n → 0, ρ1n → ρ1 , ρ0n → ρ0 weakly in L 1 (Rd ) and the support of ρ1n , ρ0 be a subset of B R , a disk of radius R. Then as n → 0, 1, p 

∇Ψ n → ∇Ψ



uniquely in W loc Rd ,

(26)

for 1  p < ∞. The proof of Theorem 4.2 is a direct consequence of the following facts. Lemma 4.1. If Ψ n : Rd → R is convex for all n > 0 and for each compact set K ,



   Ψ n   C K

K

where C K is a constant depending only on K , then there exist Ψ : Rd → R convex such that |∇Ψ |  M K on each K and a subsequence 1, p {Ψ nk } such that Ψ nk → Ψ in W loc (Rd ). Proof. To prove this result, it suffices to bound



  ∇Ψ n  dx and

K



  det ∇ 2 Ψ n  dx.

K

The first bound is given by Theorem 1 of Evans and Gariepy [18]. In fact, that theorem gives the stronger result

     Ψ n , ∇Ψ n   C K .

(27)

By (27) and the Ascoli–Arzela theorem (e.g., see [17]), up to a subsequence Ψ n converges uniformly to a convex function Ψ and ∇Ψ n converges weakly to ∇Ψ . Next, let an ,1 , . . . , an ,d be the eigenvalues of ∇ 2 Ψ n . Then

  det ∇ 2 Ψ n 2 = a2

d

n ,1 + · · · + an ,d

and

Ψ n = an ,1 + · · · + an ,d , where ∇ 2 is the Hessian matrix and is the Laplacian operator. Because Ψ n is convex, the eigenvalues are nonnegative. Thus,

  √ det ∇ 2 Ψ n   d Ψ n .

(28)

But if B is a ball then by the divergence theorem





Ψ n = B



∇Ψ n · n dσ  ∂B

  ∇Ψ n  dσ  C K σ (∂ B ),

(29)

∂B

where n is the unit outward normal to the boundary ∂ B. To obtain the last inequality, we have used (27). By (28) and (29), we have



  √ det ∇ 2 Ψ n   dC K σ (∂ B ).

(30)

K

Using (27) and (30), we have proven that ∇Ψ n is bounded in B V ( K ) (functions of bounded variation) for each compact set. By Theorem 4 of [18], up to a subsequence, we have that ∇Ψ n converges in L 1 ( K ) to ∇Ψ. Since by (27) ∇Ψ n is bounded in L ∞ ( K ), we conclude that up to that subsequence, ∇Ψ n converges in L p ( K ) to ∇Ψ . 2 Remark 4.2. The uniqueness of the solution of the Monge–Ampère equation (8) implies that the whole sequence {Ψ n } converges to Ψ . To see this, assume that is not true, i.e., there exists some positive number  ∗ such that

   Ψ n − Ψ  1, p  ∗, W (Q ) loc

Q ⊂ Rd ,

(31)

M.M. Sulman et al. / Applied Numerical Mathematics 61 (2011) 298–307

303

for infinitely many values of n. Let {Ψ nq } be the subsequence of {Ψ n } consisting of all elements of the sequence that satisfy (31), i.e.,    Ψ nq − Ψ  1, p   ∗ , ∀q. (32) W loc ( Q )

nq

Since {Ψ } satisfies (27) then by the Ascoli–Arzela theorem, there exists a subsequence {Ψ Ψ is the unique solution of the Monge–Ampère equation (8)), which contradicts (32).

nql

} that converges to Ψ (since

5. Numerical experiments In this section we discuss the details of the numerical solution of (11) and give some examples which illustrate the global existence and convergence results for the solution. We use a standard centered second order finite difference scheme for the spatial discretization of (11) and a forward Euler scheme for the temporal discretization. The simple algorithm to compute a numerical solution of (11) is the following: 1. 2. 3. 4.

Given Ψ n . Compute the approximations for ∇ 2 Ψ n and ∇Ψ n . Compute ρ1 (∇Ψ n ). Set



F n = log



ρ1 (∇Ψ n ) det ∇ 2 Ψ n . ρ0

5. Update Ψ n : Ψ n+1 = Ψ n + t F n . 6. If ∇ F n 2 < Tol, where Tol is the user defined tolerance, STOP; else set n = n + 1 and go to 2. Example 1. We begin with the simple problem of computing a mapping from one circle to another concentric one using the PMA method. For radially symmetric densities, the solution of the Monge–Ampère equation (8) can be found exactly by direct integration. Specifically, defining the polar coordinates

ξ = r cos(θ),

η = r sin(θ),

ξ = (ξ, η) ∈ Ωc ,

the Monge–Ampère equation (8) simplifies to

1 r

Ψr Ψrr =

ρ0 (r ) ρ1 (Ψr (r ))

(33)

where Ψr and Ψrr are the first and second derivatives with respect to the radial variable r. In this case, the density functions ρ0 and ρ1 are defined as



ρ0 (r ) =



ξ 2 + η2  r0 ,

c1 ,

if

c2 ,

otherwise,



ρ1 (Ψr ) =



x2 + y 2  r1 ,

d1 ,

if

d2 ,

otherwise

where x = (x, y ) ∈ Ω , r0 , r1 are the radii of the circles defined in Ωc and Ω respectively, and c 1 , c 2 , d1 , d1 are some positive real numbers. Integrating (33) once gives

Ψr =

√ 

a r2 + c

where a is the ratio ρ0 /ρ1 and c is the integration constant. Note that our primary goal is to compute the optimal mapping that transfers ρ0 to ρ1 and not the potential Ψ (r ) in its entirety. Therefore, using (7) we obtain the optimal mapping x = (x, y ) in terms of Ψr as

x=

ξ r

Ψr ,

y=

η r

Ψr .

The exact solution for the optimal transformation is shown in Fig. 2(b). We implement the PMA method to compute the numerical solution on a grid of size 33 × 33. In Fig. 2(a) we show the numerical solution for the optimal transformation (i.e., ∇Ψ∞ ). The maximal densities being transported are concentric circles; therefore, the grid points can be expected to be pushed radially towards the center of the circle. Example 2. To illustrate the global convergence of the solutions of the parabolic Monge–Ampère equation (11), we study the numerical solution of the initial boundary value problem

⎧   ρ0 (ξ ) ∂Ψ ⎪ ⎪ = log det ∇ 2 Ψ − log ⎪ ⎪ ρ1 (∇Ψ ) ⎨ ∂t ∇Ψ · n = ξ · n on ∂Ωc , ⎪ ⎪ ⎪   ⎪ ⎩ Ψ (ξ, η, 0) = 1 ξ 2 + η2 in Ωc , 2

in Ωc × [0, ∞), (34)

304

M.M. Sulman et al. / Applied Numerical Mathematics 61 (2011) 298–307

Fig. 2. Example 1. Shown here are (a) numerical solution for the optimal mapping computed using the PMA method, and (b) exact solution for the optimal mapping.

Fig. 3. Example 2. Gray-scale plots of the time evolution | ∂Ψ ∂ t | on a unit square domain at different times.

with Ωc = Ω = [0, 1]2 ⊂ R2 , where we define

ρ0 ≡ 1 and ρ1 (∇Ψ ) = ρ1 (x, y ) as   2  ρ1 (x, y ) = 1 + 5 exp −50 (x − 0.5) + ( y − 0.5)2 − 0.09 . 

(35)

In order to guarantee the convergence, we normalize the function ρ1 so that the compatibility condition (20) is satisfied. For the numerical solutions of (34), we first discretize the flow equation in space on a grid of size 41 × 41, and then integrate in time using an explicit Euler scheme. Fig. 3 presents gray-scale plots of the evolution of | ∂Ψ ∂ t | at different times. −5 throughout the unit square domain, and this reflects the global From Fig. 3(c) we observe that the | ∂Ψ ∂ t | is of order 10

convergence of Ψ to the stationary solution Ψ∞ . In Fig. 4(a), we plot δ(t ) = ∂Ψ ∂ t L 2 (Ωc ) , t > 0, and in Fig. 4(b) we plot −1

log δ(t ) which shows the expected exponential convergence. In grid application we use the steady state solution of (34) to generate the adaptive mesh for the given monitor function (35). We notice that the function ρ1 has its maximum values on a circle of radius 0.3, so using it as a monitor function with the PMA method, the adaptive mesh is expected to cluster around the circle. In Figs. 5(a) and 5(b) we plot respectively the generated adaptive mesh and the surface of the equidistribution measure E (x, y ) which is defined as t

E (x, y ) = ρ1 (x, y )J|Ωc |,

(x, y ) ∈ Ω

(36)

where |Ωc | is the area of Ωc and J is the determinant of the Jacobian matrix of the coordinate transformation x = ∇Ψ∞ . Note that E (x, y ), which measures how close the mesh comes to equidistributing the monitor function ρ (x, y ), is very close to 1 throughout Ω . Example 3. In this example, we consider the evolution problem

∂Ψ = log det ∇ 2 Ψ − log ˜f , ∂t

(37)

M.M. Sulman et al. / Applied Numerical Mathematics 61 (2011) 298–307

305

−1 Fig. 4. Example 2. Plots of (a) the decay δ(t ) = ∂Ψ ∂ t L 2 (Ωc ) and (b) the convergence rate t log δ(t ). The expected exponential convergence can be seen from plot (b).

Fig. 5. Example 2. Plots of (a) the adaptive mesh and (b) the equidistribution measure function E (x, y ) using the monitor function (35). The surface of E (x, y ) shows a very accurate equidistribution of the mesh.

where

˜f (x, y ) = 

f (x, y ) Ω f (x, y ) dx d y

.

To investigate the convergence of the solution in space and to compare the convergence results with the method of [16], we choose

  2 2   f x(ξ, η), y (ξ, η) = 1 + ξ 2 + η2 e ξ +η ,

(ξ, η) ∈ Ωc = [0, 1]2 .

The steady state solution of (37) solves the Monge–Ampère equation

det ∇ 2 Ψ = ˜f . The boundary conditions are chosen so that

Ψe = e (ξ

2

+η2 )/2

is the exact solution of (38).

(38)

306

M.M. Sulman et al. / Applied Numerical Mathematics 61 (2011) 298–307

Fig. 6. Example 3. Shown here are (a) exact solution of elliptic Monge–Ampère equation (b) numerical solution computed using PMA method, and (c) pointwise error surface.

Fig. 7. Example 3. Convergence study of the PMA method: semilog plot of L 2 -norm error in numerical solution of (37) computed using PMA method for N = 16, 32, 64, 128.

To study the convergence of the solution in space, we integrate the discretization of (37) on grids of sizes N × N, N = 16, 32, 64, and 128. In Fig. 6, we plot (a) the exact solution, (b) the computed numerical solution and (c) the pointwise error |Ψe − Ψc | measured as the absolute value of the difference between the numerical and the exact solutions on a grid of size 32 × 32. In Fig. 7 we show a semilog plot of the L 2 -norm error which is defined as

 L 2 -error =

|Ψe − Ψ |2 dξ dη

 12 .

Ωc

For the integral approximation, we use

 F (ξ, η) dξ dη =

N −1 N −1  

F˜ (ξi , η j ) ξ · η,

i =1 j =1

Ωc

where

F˜ (ξi , η j ) =

ξ = η =

1

N −1

1 4



F (ξi , η j ) + F (ξi +1, j , ηi +1, j ) + F (ξi , j +1 , ηi , j +1 ) + F (ξi +1, j +1 , ηi +1, j +1 ) ,

for an N × N grid defined on a unit square.

In Table 1, the L 2 -norm error and rate of convergence measured by the ratio log2 (error N /error2N ) are shown for N = 16, 32, 64, 128. The results shown in Table 1 are comparable to the results of the same example given in [16].

M.M. Sulman et al. / Applied Numerical Mathematics 61 (2011) 298–307

307

Table 1 Example 3. Convergence study of the PMA method. Shown are the L 2 error and convergence rate for N = 16, 32, 64, 128. N

L 2 -error

Convergence rate

16 32 64 128

8.5129 × 10−5

2.0866 2.0534 2.0620

2.0042 × 10−5 4.8285 × 10−6 1.1563 × 10−6

6. Discussion The global existence and convergence results given in this paper support numerical evidence of the robustness of the PMA method for adaptive mesh generation. Indeed, the numerical results show that the PMA method described here maintains a global equidistribution property and compares favorably with some other adaptive mesh methods, including the methods in [13,8] (for further details see [25]). The PMA method is fairly simple and straightforward to implement, and hence makes it attractive for applying the PMA method for other similar applications such as image processing [26], with some confidence. Acknowledgements The first author would like to thank the Institute of Pure and Applied Mathematics (IPAM) at UCLA for the opportunity to attend the Optimal Mass Transport Workshop and to Dr. Wilfrid Gangbo for many fruitful discussions at the workshop. The authors are grateful to a referee who made many valuable comments on an earlier version of this paper. References [1] A.D. Aleksandrov, Certain estimates for the Dirichlet problem, Soviet Math. Dokl. 1 (1961) 1151–1154. [2] L. Ambrosio, L.A. Caffarelli, Y. Brenier, G. Buttazzo, C. Villani, Optimal Transportation and Applications, Lecture Notes in Mathematics, vol. 1813, Springer-Verlag, Berlin, 2003, lectures from the C.I.M.E. Summer School held in Martina Franca, September 2–8, 2001, edited by Caffarelli and S. Salsa. [3] J.-D. Benamou, Y. Brenier, A domain decomposition method for the polar factorization of vector fields, in: Domain Decomposition Methods in Science and Engineering, Como, 1992, in: Contemp. Math., vol. 157, Amer. Math. Soc., Providence, RI, 1994, pp. 231–236. [4] J.-D. Benamou, Y. Brenier, A computational fluid mechanics solution to the Monge–Kantorovich mass transfer problem, Numer. Math. 84 (3) (2000) 375–393. [5] J.-D. Benamou, Y. Brenier, K. Guittet, The Monge–Kantorovich mass transfer and its computational fluid mechanics formulation, Internat. J. Numer. Methods Fluids 40 (1–2) (2002) 21–30, iCFD Conference on Numerical Methods for Fluid Dynamics (Oxford, 2001). [6] J.-D. Benamou, Y. Brenier, K. Guittet, The Monge–Kantorovich mass transfer and its computational fluid mechanics formulation, in: iCFD Conference on Numerical Methods for Fluid Dynamics, Oxford, 2001, Internat. J. Numer. Methods Fluids 40 (1–2) (2002) 21–30. [7] Y. Brenier, Polar factorization and monotone rearrangement of vector-valued functions, Comm. Pure Appl. Math. 44 (4) (1991) 375–417. [8] C.J. Budd, J.F. Williams, Parabolic Monge–Ampère methods for blow-up problems in several spatial dimensions, J. Phys. A 39 (19) (2006) 5425–5444. [9] L.A. Caffarelli, Interior W 2, p estimates for solutions of the Monge–Ampère equation, Ann. of Math. (2) 131 (1) (1990) 135–150. [10] L. Caffarelli, The regularity of mappings with convex potentials, in: Partial Differential Equations and Related Subjects, Trento, 1990, in: Pitman Res. Notes Math. Ser., vol. 269, Longman Sci. Tech., Harlow, 1992, pp. 53–58. [11] L.A. Caffarelli, Boundary regularity of maps with convex potentials, Comm. Pure Appl. Math. 45 (9) (1992) 1141–1151. [12] L.A. Caffarelli, The regularity of mappings with a convex potential, J. Amer. Math. Soc. 5 (1) (1992) 99–104. [13] W. Cao, W. Huang, R.D. Russell, A moving mesh method based on the geometric conservation law, SIAM J. Sci. Comput. 24 (1) (2002) 118–142 (electronic). [14] E.J. Dean, R. Glowinski, Numerical solution of the two-dimensional elliptic Monge–Ampère equation with Dirichlet boundary conditions: an augmented Lagrangian approach, C. R. Math. Acad. Sci. Paris 336 (9) (2003) 779–784. [15] E.J. Dean, R. Glowinski, Numerical solution of the two-dimensional elliptic Monge–Ampère equation with Dirichlet boundary conditions: a least-squares approach, C. R. Math. Acad. Sci. Paris 339 (12) (2004) 887–892. [16] E.J. Dean, R. Glowinski, Numerical methods for fully nonlinear elliptic equations of the Monge–Ampère type, Comput. Methods Appl. Mech. Engrg. 195 (13–16) (2006) 1344–1386. [17] J.H. Dshalalow, Real Analysis, Studies in Advanced Mathematics, Chapman & Hall/CRC, Boca Raton, FL, 2001, an introduction to the theory of real functions and integration. [18] L.C. Evans, R.F. Gariepy, Measure Theory and Fine Properties of Functions, Studies in Advanced Mathematics, CRC Press, Boca Raton, FL, 1992. [19] K. Guittet, On the time-continuous mass transport problem and its approximation by augmented Lagrangian techniques, SIAM J. Numer. Anal. 41 (1) (2003) 382–399, (electronic). [20] M. Knott, C.S. Smith, On the optimal mapping of distributions, J. Optim. Theory Appl. 43 (1) (1984) 39–49. [21] G. Loeper, F. Rapetti, Numerical solution of the Monge–Ampère equation by a Newton’s algorithm, C. R. Math. Acad. Sci. Paris 340 (4) (2005) 319–324. [22] O.C. Schnürer, Translating solutions to the second boundary value problem for curvature flows, Manuscripta Math. 108 (3) (2002) 319–347. [23] O.C. Schnürer, H.R. Schwetlick, Translating solutions for Gauss curvature flows with Neumann boundary conditions, Pacific J. Math. 213 (1) (2004) 89–109. [24] O.C. Schnürer, K. Smoczyk, Neumann and second boundary value problems for Hessian and Gauss curvature flows, Ann. Inst. H. Poincaré Anal. Non Linéaire 20 (6) (2003) 1043–1073. [25] M. Sulman, J. Williams, R.D. Russell, Optimal mass transport for higher dimensional adaptive grid generation, submitted for publication. [26] M.H. Sulman, J. Williams, R.D. Russell, F.M. Beg, Volumetric image registration methods based on solving the Monge–Ampère equation, manuscript. [27] J. Urbas, On the second boundary value problem for equations of Monge–Ampère type, J. Reine Angew. Math. 487 (1997) 115–124.