Optimisation approach for the Monge-Ampere equation

0.55 MB·

by  Fethi Ben Belgacem

OPTIMIZATION APPROACH FOR THE MONGE-AMPÈRE
EQUATION
7
1
0 fethi ben belgacem 1
2

n
a Abstract. This paper studies the numerical approximation of solution of
J theDirichletproblemforthefullynonlinearMonge-Ampèreequation. Inthis
approach,wetaketheadvantageofreformulationtheMonge-Ampèreproblem
9
as an optimization problem, to which we associate a well defined functional
1
whoseminimumprovidesuswiththesolutiontotheMonge-Ampèreproblem

after resolvingaPoissonproblem bythe finiteelementGalerkinmethod. We
]
P presentsomenumericalexamples,forwhichagoodapproximationisobtained
A in68iterations.
Keywords. ellipticMonge-Ampèreequation,gradientconjugatemethod,
h. finiteelementGalerkinmethod.
t AMS subject classifications. 35J60,65K10,65N30.
a
m
[ 1. Introduction

1 In this paper, we give a numerical solution for the following Monge-Ampère
v
problem
8
8 det[D2u]=f(x) x Ω,
3 (1.1) ∈
u =0, uconvex,
5 (cid:26) |Γ
0 where Ω is a smooth convex and bounded domain in R2, D2u is the Hessian
.
1 of u and f C∞(Ω), f >0.
0 ∈ (cid:2) (cid:3)
Equation (1.1) belongs to the class of fully nonlinear elliptic equation. The
7
mathematical analysis of real Monge-Ampère and related equations has been a
1
: source of intense investigations in the last decades; let us mention the following
v
references ( among many others and in addition to [7], [9], [15]): [10], [8], [17,
i
X chapter 4], [2], [28], [11]-[14]. Applications to Mechanics and Physics can be found
r in [27], [4], [5], [18], [24], [26],[31], (see also the references therein).
a
The numericalapproximationsofthe Monge-Ampèreequationaswellasrelated
equations have recently been reported in the literature. Let us mention the ref-
erences [4], [29], [39], [26 ], [11], [32], [25],[28], [33]; the method discussed in [11],
[32],[25] is very geometrical in nature. In contrast with the method introduced by
Dean and Glowinski in [19 ] [20] [21], which is of the variational type.
On the existence of smooth solution for (1.1), we recall that if f C∞(Ω)
∈
equation (1.1) has a unique strictly convex solution u C∞(Ω) (see [14]).
∈
1
Department of Mathematics, HigherInstitute of Computer Sciences and Mathematics of
Monastir (ISIMM), Avenue de la Korniche - BP 223 - Monastir - 5000, TUNISIA.
([email protected]).
1
OPTIMIZATION APPROACH FOR THE MONGE-AMPÈRE EQUATION 2
To obtain a numerical solution for (1.1), we propose a least-square formulation
of (1.1). In this approach, we take the advantage of reformulation of the Monge-
Ampère problem as a well defined optimization problem, to which we associate a
wellfunctionalwhoseminimumprovidesuswiththesolutiontotheMonge-Ampère
problem after resolving a Poisson problem by the finite element Galerkin method.
The minimum is computed by the conjugate gradient method.
The remainder of this article is organized as follows. In section 2, We introduce
the optimization problem. In section 3, we discuss a conjugate gradient algorithm
for the resolution of the optimization problem. The finite element implementation
oftheabovealgorithmisdiscussedinsection4. Finally,insection5,weshowsome
numerical results.
2. Formulation of the Dirichlet problem for the elliptic
Monge-Ampère equation
Let u be the solution of (1.1). Let λ and λ be the eigenvalues of the matrix
I 1 2
[D2u ]. We have
I
λ +λ = ∆u ,
1 2 I
λ λ = det[D2u ]=f.
1 2 I
(cid:26)
Then λ and λ are the solutions of the equation
1 2
X2 ∆u X+f =0.
I
−
So
(∆u )2 4f 0.
I
− ≥
Then
∆u 2 f 0.
I
− ≥
Let us set
p
∆u 2 f =g C∞(Ω).
I
− ∈
We conclude that u is solution of the following Dirichlet Poisson problem
I p
e
∆u=2√f +g,
Peg u|Γ =0.
(cid:26)
To compute g, we consider the least-squares functeional J defined on
E = ϕ C∞(Ω), ϕ 0 ,
∈ ≥
e
as follows: (cid:8) (cid:9)
J(g)= 1 det D2ug f 2dx,
2ˆ −
Ω
where ug is the solution of the Dirich(cid:0)let P(cid:0)oisson(cid:1)prob(cid:1)lem
∆u=2√f +g,
g
P u|Γ =0.
(cid:26)
The minimization problem
g E,
(2.1) ∈
J(g) J(g) g E,
(cid:26) ≤ ∀ ∈
is thus a least-squares formulationeof (1.1).
e
Theorem 1. u is the strictly convex solution of (1.1) if and only if there exist a
I
unique solution g of (2.1) such that u =ug˜.
I
e
OPTIMIZATION APPROACH FOR THE MONGE-AMPÈRE EQUATION 3
Proof. Since u is solutionof( g˜), wehaveu =ug˜. So J(g)=0andg˜is aunique
I I
P
solution of (2.1).
Conversely, let g be a solution of (2.1). Since (1.1) has a solution u , we can
e I
deduce immediatly that J(g)=0 and so, J(g)=0. It follows that
det[D2ug¯]=f,
e ug¯ =0.
(cid:26) |Γ
We have ∆ug¯ = 2f +g > 0 and det[D2ug¯] > 0, we can deduce that ug¯ is strictly
convex and from the uniqueness of solution for (1.1) we get ug¯ =u . (cid:3)
I
3. Iterative solution for the minimisation problem
3.1. Descriptionofthealgorithm. Thealgorithmweconsidertosolvetheprob-
lem (2.1) which is basedon the PRP (Polak-Ribière-Polyak[36,37])conjugate gra-
dient method reads:
Given g0 E;
∈
then, for k 0, gk being known in E, solve
≥
gk ∆u=2√f +gk,
P u|Γ =0.
(cid:26)
Compute, J(gk),
∇
J(gk)T( J(gk) J(gk−1)
If k 1, βk = ∇ ∇ −∇ ;
≥ J(gk−1) 2
k∇ k2
J(g0) ifk=0
dk = −∇
J(gk)+βkdk−1 ifk 1;
(cid:26) −∇ ≥
and update gk by
gk+1 =gk+αkdk.
Where α is computed with the Armijo-type line search.
k
3.2. Solution of sub-problem ( g).. We consider first the variational formula-
P
tion of ( g)
P
Find ug H1(Ω), such that,
(3.1) ∈ 0
a(ug,v)=L(v), v H1(Ω),
(cid:26) ∀ ∈ 0
where
(3.2) a(u,v)= u vdx
ˆ ∇ ∇
Ω
and
(3.3) L(v)= 2 f +g vdx,
−ˆ
Ω
(cid:16) p (cid:17)
a in (3.2) is coercive on H1(Ω). For f L2(Ω)we have √f L2(Ω). Since Ω is
0 ∈ ∈
bounded and for g L2(Ω), L in (3.3) is continuous, then by the Lax-Milgram
∈
theorem ( g) has a unique solution ug.
V
P
OPTIMIZATION APPROACH FOR THE MONGE-AMPÈRE EQUATION 4
4. Finite element approximation of the minimization problem
For simplicity, we assume that Ω is a bounded polygonal domain of R2. Let
h
T
a finite triangulation of Ω (like those discussed in e.g, [16]).
We introduce a
with P the space of the two-variable polynomials of degree 1. A function ϕ
1
being given in H2(Ω) we denote ∂2ϕ by ≤
∂xixj
∂2ϕ ∂ϕ ∂v
(4.1) vdx= dx, v H1(Ω), i=1,2,
ˆ ∂x2 −ˆ ∂x ∂x ∀ ∈ 0 ∀
Ω i Ω i i
∂2ϕ 1 ∂ϕ ∂v ∂ϕ ∂v
(4.2) vdx= + dx, v H1(Ω).
ˆ ∂x x −2ˆ ∂x ∂x ∂x ∂x ∀ ∈ 0
Ω 1 2 Ω(cid:20) 1 2 2 1(cid:21)
Let ϕ V ; taking advantage of relations (4.1) and (4.2) we define the discrete
h
∈
analogues of the differential operators D2 by
ij
i=1,2, D2 (ϕ) V ,
∀ hii ∈ 0h
(4.3) ∂ϕ ∂v
 ´ΩDh2ii(ϕ)vdx=−ˆΩ ∂xi∂xidx, ∀v ∈V0h,
 D2 (ϕ) V ,
h12 ∈ 0h
(4.4) 1 ∂ϕ ∂v ∂ϕ ∂v
 ´ΩDh212(ϕ)vdx=−2ˆΩ(cid:20)∂x1∂x2 + ∂x2∂x1(cid:21)dx, ∀v ∈V0h.
To compute the above discrete second order partial derivatives we will use the

trapezoidal rule to evaluate the inegrals in the left hand sides of (4.3) and (4.4).
We consider the set Ph of the vertices of Th and P0h={P|P ∈Ph,P ∈/ Γ}. We
define the integers Nh and N0h by Nh = Card(Ph) and N0h = Card(P0h). So
dimV =N and dimV =N .
h h 0h 0h
For Pk ∈Ph we associate the function wk uniquely defined by
w V , w (P )=1, w (P )=0, ifl=1,...N l=k.
k h k k k l h,
∈ 6
It is wellknown(e.g., [16])that the sets B = w Nh andB = w N0h are
h { k}k=1 0h { k}k=1
vector bases of V and V , respectively.
h 0h
We denote by A the area of the polygonalwhich is the union of those triangles
k
of which have P as a common vertex. By applying the trapezoidal rule to the
h k
T
integrals in the left hand side of relations (4.3) and (4.4) we obtain:
i=1,2, D2 (ϕ) V ,
∀ hii ∈ 0h
(4.5)  Dh2ii(ϕ)(Pk)=−A3k ˆΩ ∂∂xϕi∂∂wxki dx, ∀k =1,2,...,N0h,
 D2 (ϕ) =D2 (ϕ) V ,
h12 h21 ∈ 0h
(4.6)  Dh212(ϕ)(Pk)=−2A3(cid:0)k ˆΩ(cid:20)∂∂xϕ1∂∂(cid:1)wx2k + ∂∂xϕ2∂∂wx1k(cid:21)dx, ∀k =1,2,...,N0h.
Computingthe integralsinthe righthandsidesof(4.5)and(4.6)isquitesimple

since the first order derivatives of ϕ and w are piecewise constant.
k
Taking the above relations into account. We approximate the space E by
E = ϕ V , ϕ 0 ,
h h
{ ∈ ≥ }
OPTIMIZATION APPROACH FOR THE MONGE-AMPÈRE EQUATION 5
and then the minimization problem (2.1) by
g E ,
h h
∈
J (g ) J (g ) g E ,
h h h h h h
(cid:26) ≤ ∀ ∈
e
Where
e
J (g )= 1Nh0A D2 (ugh)(P )D2 (ugh)(P ) D2 (ugh(P ) 2 f (P ) 2, g E ,
h h 6 k h11 h k h22 h k − h12 h k − h k ∀ h ∈ h
Xk=1 (cid:12)(cid:12) (cid:0) (cid:1) (cid:12)(cid:12)
andf ,g ,g arer(cid:12)espectivelyacontinuousapproximationsoffunctionsf,g,(cid:12)g and
h h h
ugh is the solution of the discret variant of the Dirichlet Poisson problem ( g).
h P
e e
4.1. Discrete variant of the algorithm. We will discuss now the solution of
(2.1) by a discrete variant of algorithm 3.1.
Given g0 E ;
h ∈ h
then, for k ≥0, ghk being known in Eh, solve Pghk,
Compute, J(gk), γk = J(gk) 2;
∇ h h ∇ h 2
If k 1, βk =γk/γ(cid:13)k−1; (cid:13)
≥ h h (cid:13)h (cid:13)
J(g0) ifk =0
dk = −∇
h J(gk)+βkdk−1 ifk 1;
(cid:26) −∇ h h h ≥
and update gk by
h
gk+1 =gk+αkdk.
h h h h
Remark 2. There are many approches for finding an avaible step size αk. Among
h
themtheexactlinesearchisanidealone,butiscost-consumingorevenimpossible
to use to find the step size. Some inexact line searches are sometimes useful and
effective in practical computation, such as Armijo line search [1], Goldstein line
search and Wolfe line search [24,38].
The Armijo line search is commonly used and easy to implement in practical
computation.
Armijo line search
Let s > 0 be a constant, ρ (0,1) and µ (0,1). Choose α to be the largest
k
∈ ∈
α in s, sρ, sρ2,..., such that
(cid:8) J(cid:9) (gk) J (x +αdk) αµ J (gk)Tdk.
h h − h k h ≥− ∇ h h h
However, this line search cannot guarantee the global convergence of the PRP
method and even cannot guarantee d to be descent direction of J at gk.
k
fi4.n1i.t1e.dSimoleuntisoionnoaflvthaeriasutibonparolblilnemearPpgrhko.blAemnywshuibc-hprroeabdlesmas(Pfogllhko)w,sis: Feqinudivualgehnt tVo a
h ∈ 0h
such that
(4.7) a(ugh,v )=L(v ), v V .
h h h ∀ h ∈ 0h
By the Lax-milgram theorem we can easily show that (4.7) has a unique solution
ugh V .
h ∈ 0h
OPTIMIZATION APPROACH FOR THE MONGE-AMPÈRE EQUATION 6
5. Numerical experiments
Inthissectionwearegoingtoapplythemethoddiscussedintheprevioussection
to the solution of some test problems. For all these test problems we shall assume
that Ω is the unit disk. We first approximate Ω by a polygonal domain Ω . We
h
consider a finite triangulation of Ω .
h h
T
The first test problem is expressed as follows
det[D2u]=4 1+2 x2 e2(|x|2−1) x Ω,
(5.1) | | ∈
( (cid:16)u|Γ =(cid:16)0, (cid:17)(cid:17) uconvex.
with x2 =x2+x2. The exact solution u C∞(Ω¯) to problem (5.1) is given by
| | 1 2 ∈
u(x)=e(|x|2−1) 1.
−
Remark 3. When computing the approximate solutions of these problems, we
stopped the iterations of the algorithm as soon as J (g ) 10−6.
h h
| |≤
We have discretized the optimization problem associated to the problem (5.1).
We solved the Poisson problem encountred at each iteration of the algorithm by a
fast Poisson solvers.
We have used as initial guess three different constant values for g0. The results
h
obtainedafter 68iterationsaresummarizedinTable 1(where uc denotes the com-
h
puted approximate solution and . = . ).
kk0,Ω kkL2(Ω)
Thegraphofuchanditscontourplotobtained,forh=1/128hasbeenrespectively
visualized on Figure 2 and Figure 3.
Table 1. First test problem : Convergence of the approximate solution.
h g0 u uc
h k − hk0,Ω
1/32 0.1 0.8861 10−4
×
1/32 0.2 0.5497 10−4
×
1/32 0.3 0.3720 10−4
×
1/64 0.1 0.3416 10−4
×
1/64 0.2 0.9121 10−5
×
1/64 0.3 0.7554 10−5
×
1/128 0.1 0.6305 10−5
×
1/128 0.2 0.4981 10−5
×
1/128 0.3 0.7203 10−6
×
We conclude from the results in Table 1 that the value g0 = 0.3 is optimal and
quite accurate approximations of the exact solutions are obtained.
Remark 4. Wedidnottrytofindtheoptimalvalueofg0(itseemsthatisadifficult
problem).
In the second test problem we take
2
4 π π
f(x)= π2 cos2 1 x2 + x2 sin π 1 x2 .
5 2 −| | 2 | | −| |
(cid:18) (cid:19) h (cid:16) (cid:16) (cid:17)(cid:17) (cid:16) (cid:17) (cid:16) (cid:16) (cid:17)(cid:17)i
OPTIMIZATION APPROACH FOR THE MONGE-AMPÈRE EQUATION 7
Figure 5.1. First test problem : Graph of uc.
h
Figure 5.2. First test problem : Contour plot of uc.
h
The solution to the corresponding Monge-Ampère problem is the function u
C∞(Ω¯) defined by ∈
4 π
u(x)= sin 1 x2 .
−5 2 −| |
The method provides after 64 iterations(cid:16)the(cid:16)results su(cid:17)m(cid:17)murized in Table 2.
The value g0 =0.3 is again optimal.
The third test problem is defined as follows
det[D2u]=1 x Ω,
(5.2) ∈
u =0, uconvex.
|Γ
(cid:26)
The function u given by
1
u(x)= x2 1
2 | | −
is the solution of (5.2) and u C∞(Ω¯). (cid:16) (cid:17)
∈
We deduce from Table 3 that g =0.2 is an optimal value.
0
UnfortunalyIdidnotfindanyotherinitialvaluethatgivesmoreaccurateresults.
Even for g =0.4 the results are not satisfied.
0
OPTIMIZATION APPROACH FOR THE MONGE-AMPÈRE EQUATION 8
Table 2. second test problem : Convergence of the approximate solution.
h g0 u uc
h k − hk0,Ω
1/32 0.1 0.6466 10−4
×
1/32 0.2 0.4510 10−4
×
1/32 0.3 0.2983 10−4
×
1/64 0.1 0.1749 10−4
×
1/64 0.2 0.8507 10−5
×
1/64 0.3 0.6221 10−5
×
1/128 0.1 0.3743 10−5
×
1/128 0.2 0.1180 10−5
×
1/128 0.3 0.5591 10−6
×
Figure 5.3. Second test problem : Graph of uc.
h
Figure 5.4. Second test problem : Contour plot of uc.
h
References
[1] L.Armijo,MinimizationoffunctionshavingLipschitscontinuouspartialderivatives,Pacific
J.Math. 16(1966) 1-3.
OPTIMIZATION APPROACH FOR THE MONGE-AMPÈRE EQUATION 9
Table 3. Third test problem : Convergence of the approximate solution.
h g0 u uc
h k − hk0,Ω
1/32 0.1 0.3830 10−3
×
1/32 0.2 0.2564 10−3
×
1/32 0.3 0.2971 10−3
×
1/64 0.1 0.7448 10−4
×
1/64 0.2 0.8529 10−6
×
1/64 0.3 0.9193 10−5
×
1/128 0.1 0.6215 10−5
×
1/128 0.2 0.5837 10−6
×
1/128 0.3 0.3806 10−5
×
Figure 5.5. Third test problem : Graph of uc.
h
Figure 5.6. Third test problem : Contour plot of uc.
h
[2] T. Aubin, Nonlinear analysis on manifolds, Monge-Ampère equations, Springer-Verlag,
Berlin,1982.
[3] J. D. BenamouandY. Brenier,The Monge-Kantorovich mass transfer and itscomputa-
tional fluid mechanics formulation, Int.J.Numer.Meth.Fluids,2002.
OPTIMIZATION APPROACH FOR THE MONGE-AMPÈRE EQUATION 10
[4] J-D. Benamou, B. D. Frose, A. M. Oberman, Two numerical methods for the elliptic
Monge-Ampere equation, ESSAIM:M2AN44(2010)737-758.
[5] O. Bokanowski, B. Grébert, Deformations of density functions in molecular quantum
chemistry,J.Math.Phys.37(1996), no.4,1553-1573.
[6] K.Böhmer,onfiniteelementmethodsforfullynonlinearellipticequationsofsecondorder,
SIAMJ.Numer.Anal.,46(3):1212-1249, 2008.
[7] Y. Brenier, Some geometric PDEs related to hydrodynamics and electrodynamics, inPro-
ceedings of the International Congress of Mathematicians, Beijing 2002, August 20-28, Vol.
III,HigherEductionPress,Beijing,2002,pp.761-771.
[8] X.Cabré, Topicsinregularity andqualitativespropertiesof solutionsofnonlinearelliptic
equations, Current developments inpartial differential equations (Temucooo,1999). Discrete
Contin.Dyn.Syst.8(2002). pp.289-302.
[9] L. Caffarelli, Non linear elliptic theory and the Monge-Ampère equation, inProceedings
oftheInternational Congress ofMathematicians, Beijing2002, August20-28, Vol.I,Higher
EductionPress,Beijing,2002,pp.179-187.
[10] L.CaffarelliandX.Cabré,FullynonlinearEllipticEquations,AmericanMathematical
society,Providence, RI, 1995.
[11] L.Caffarelli,S.A.Kochengin,V.I.Oliker,OntheNumericalsolutionoftheproblem
of reflector design with given far-field scattering data, Contemporary Mathematics, 226,
(1999), 13-32.(Eds.L. A. CaffarelliandM. Milman)
[12] L. Caffarelli, and Y. Y. Li, A liouville theorem for solutions of the Monge-Ampère
equation with periodic data, Ann. Inst. H. Poincarè Anal. Non Linéaire, 21 (2004), pp. 97-
120.
[13] L. Caffarelli, L. Nirenberg, J. Spruck, The Dirichlet problem for nonlinear second-
order elliptic equation I. Monge-Ampère equation, Comm. on Pure and Appl. Math., 17,
(1984), 396-402.
[14] L.Caffarelli,LNirenberg,J.Spruck,TheDirichletproblemforthedegenerateMonge-
Ampère equation,Rev.Mat.Ibero.,2,(1985), 19-27.
[15] S-Y. A. Chang and P. C. Yang, Nonlinear partial differential equations in conformal
geometry, in Proceedings of the International Congress of Mathematicians, Beijing 2002,
August20-28,Vol.I,HigherEductionPress,Beijing,2002, pp.189-207.
[16] P. G. Ciarlet, The FiniteElement Method for EllipticProblems,North-Holland, Amster-
dam(1978).
[17] R.Courant,D.Hilbert,Methodsofmathematicalphysics,Vol.II,NewYorkInterscience
Publishers,Wiley(1962).
[18] M.J.PCullenandR.J.Douglas,applicationsoftheMonge-AmpèreequationandMonge
transportproblemtometorologyandoceanog-raphy,contemporaryMathematics,226,(1999),
33-53.
[19] 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.Acad.Sci.Paris,Ser.1336(2003) 779-784.
[20] 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.Acad.
Sci.Paris,Ser.1339(2004)887-892.
[21] E.J.Dean,R. Glowinski,Numericalmethodsforfullynonlinearellipticequationsof the
Monge-Ampère type, Comput. Methods Appl. Mech. Engrg., 195(13-16):1344-1386, 2006.
[22] E.J. Dean,R. Glowinski, Onthenumericalsolutionof theellipticMonge-Ampère equa-
tionindimensiontwo: aleast-squaresapproach. InPartialdifferentialequations,volume16
of Comput. Methods Appl. sci., pages 43-63, Springer, Dordrecht, 2008.
[23] X.FengandM.Neilan,Vanishingmomentmethodandmomentsolutionsforsecondorder
fullynonlinearpartialdifferentialequations,J.Scient.Comp.,DOI10.1007/s10915-008-9221-
9,2008.
[24] A. A. Goldstein,Onsteepest descent,SIAMJ.Control 3(1965) 147-151.
[25] G. J. Haltiner,Numerical weather prediction, NewYork,Wiley(1971).
[26] S.A. Kochengin,V.L. Oliker,Determinationofreflectorsurfacesfromnear-filedscatter-
ingdata, InverseProblems13(1997), no.2,363-373.
[27] F. X. Le Dimet and M. Ouberdous, Retrieval of balanced fields: an optimal control
method,Tellus,45A(1993), pp.449-461.