Monge–Ampère Equations with Applications in Optic Design
Von der Fakultät für Mathematik, Informatik und Naturwissenschaften der RWTH Aachen University zur Erlangung des akademischen Grades
eines Doktors der Naturwissenschaften genehmigte Dissertation
vorgelegt von
Elisa Kremer, Master of Science aus Arnstadt, Deutschland
Berichter: apl. Prof. Dr. Siegfried Müller Univ.-Prof. Dr. Wolfgang Dahmen
Tag der mündlichen Prüfung: 26. Mai 2020
Diese Dissertation ist auf den Internetseiten der Universitätsbibliothek online verfügbar.
Declaration of Authorship
I, Elisa Kremer, declare that this thesis and the work presented in it are my own and has been generated by me as the result of my own original research. I do solemnely swear that:
1. This work was done wholly or mainly while in candidature for the doctoral degree at this faculty and university;
2. Where any part of this thesis has previously been submitted for a degree or any other qualification at this university or any other institution, this has been clearly stated;
3. Where I have consulted the published work of others or myself, this is always clearly attributed;
4. Where I have quoted from the work of others or myself, the source is always given. This thesis is entirely my own work, with the exception of such quotations;
5. I have acknowledged all major sources of assistance;
6. Where the thesis is based on work done by myself jointly with others, I have made clear exactly what was done by others and what I have contributed myself;
7. No part of this work has been published before.
Aachen, May 17
iii
Danksagungen
Ich bedanke mich herzlich bei meinen Betreuern Prof. Müller und Prof.
Dahmen für Ihre Unterstützung bei der Erstellung dieser Arbeit. Sie hat- ten immer ein offenes Ohr und lieferten wertvolle Anregungen, ließen mir auch die Freiräume und Zeit die Arbeit nach meinen eigenen Vorstellun- gen umzusetzen. Ich bedanke mich bei Prof. Dahmen dafür, dass ich diese Arbeit durchführen konnte, obwohl er nur für eine begrenzte Zeit während der Erstellung in Aachen lehrte und bei Prof. Müller, dass er sich meiner fachfremden Dissertation annahm und die notwendigen Impulse für die Fer- tigstellung gab.
Ich bedanke mich bei der DFG für die unkomplizierte Förderung (Pro- jektnummer 259180742). Ohne die vorhergehenden Arbeiten von Andreas Platen und Kolja Brix wäre die Realisierung nicht möglich gewesen. Ins- besondere Koljas unermüdlicher Einsatz und immer währendes Interesse haben das Projekt ungemein bereichert.
Ich bin für die Zusammenarbeit bei allen Projektpartnern vom Fraunhofer ILT und Lehrstuhl für Technologie optischer Systeme dankbar, besonders für Ihren Einsatz das Projekt meinetwegen zu verlängern und die Diskussio- nen mit Rolf Wester, Annika Völl und Michael Berens.
Für weitreichenden Austausch in Workshops und Symposien über Fragen zu Optimal Transport und Monge–Ampère-Gleichungen in der Optik möchte ich mich bei Prof. ten thije Boonkkamp und Prof. Benamou bedanken.
Dann bedanke ich mich für die Unterstützung durch Freunde und Kolle- gen. Sie sind nicht ein wegzudenkender Bestandteil meiner Promotionszeit, versüßten Kaffee- und Mittagspausen und bildeten übergreifende Leidens- gemeinschaften. Dabei war vor allem Valentina König immer für mich da, glücklicherweise durfte ich sowohl Oberstufe, Studium, als auch eine Pro- motion am selben Lehrstuhl mit ihr teilen. Genauso standen mir meine guten Freunde Jan Rosendahl, Katharina Plückers und meine Mädels Sab- rina Kielmann, Maren Köhne, Ragna Soennecken und Iris Trosien, immer zur Seite. Für ein erfülltes Büroleben und eine tiefere Reflektion über die Welt danke ich Dima Moser.
Zu guter letzter möchte ich bei meiner Familie, inklusive meiner engagierten Schwiegerfamilie bedanken. Mein Mann war es nie leid, sich für meine Ar- beit zu interessieren und Code zu debuggen und meine beiden Kinder haben tapfer viele Stunden Mamas langweiliger Arbeit ertragen.
v
Abstract
Since antiquity humans are interested in the design of mirrors and lenses to create light distributions. Mathematically the problem is strongly con- nected to the optimal transport problem, where for an efficient transport between two different distributions is asked for. In both, illumination op- tics and optimal transport, the problem can be transferred to solving a fully nonlinear partial differential equation of Monge–Ampère type equipped with non-standard transport boundary conditions.
In this thesis we develop a versatile, robust numerical method to solve par- tial differential equations of such a type. Versatile means that the method works for most illumination problems of practical interest, especially the difficult near field case, and robust means that the algorithm contains as few parameters as possible and does not have to be tweaked in a long trial and error process for every new configuration. The method is based on aC1 finite element projection method designed to solve fully nonlinear partial differential equations with Dirichlet conditions. We discuss the key ele- ments implementing such complicated elements and how to incorporate the nonlinear transport boundary conditions in the finite element method. In order to solve practically relevant problems, we identify the critical parts of the algorithm and introduce stabilisation techniques and inevitable regular- isation within a nested iteration process. We confirm with a wide range of benchmark problems that the resulting algorithm yields numerical accuracy, geometric flexibility and high order of approximation for smooth solutions in both, classical optimal transport problems and optic design problems.
In a second part we discuss potential exploitation in the context of extended light sources. Although mathematically the hereto considered illumination system already provides major challenges, during the design of compact op- tics the size of the source is not negligible and, hence, these systems cannot be modelled by a partial differential equation of Monge–Ampère type. Cur- rent design techniques like the phase space method are based on least-squares methods and often yield suboptimal results. Typically these methods are ini- tialised with surfaces designed assuming a point source approximation. We develop a merging process that takes several different point source approxi- mations into account and could thereby improve the overall performance of the phase space method. The guiding idea for the merging procedure is a convex combination of optical surfaces based on a unified parametrisation.
1 Introduction 1 1.1 Brief history of illumination problems and optimal transport 1
1.2 Goal and structure of the thesis . . . 3
2 Optimal transport and PDEs of MA type 5 2.1 The optimal transport problem . . . 5
2.2 PDE of optimal transport . . . 9
2.3 Numerical challenges . . . 12
2.4 Methods for optimal transport . . . 14
2.5 Galerkin methods for PDEs of MA type . . . 15
3 A FE method for PDEs of MA type 19 3.1 Outline of the derivation . . . 19
3.2 Infinite-dimensional formulation . . . 20
3.3 Discretisation . . . 26
3.4 Regularisation and stabilisation . . . 31
3.5 Implementation remarks . . . 38
4 Numerical results 41 4.1 Interpretation of results . . . 41
4.2 Simple test problem . . . 45
4.3 Mapping an ellipse to an ellipse . . . 49
4.4 Quadrature of the circle . . . 57
4.5 Image densities . . . 63
4.6 Comparison and conclusion . . . 82
5 Extension to optics 87 5.1 The inverse reflector/refractor problem . . . 88
5.2 Cases with zero etendue . . . 91
5.3 State of the art and related work . . . 91
5.4 PDEs describing optical surfaces . . . 95
5.5 Applying the conforming finite element method . . . 104
CONTENTS vii
6 Numerical results 107
6.1 Reflective surface with parallel light . . . 107
6.2 Refractive surface with parallel light . . . 114
6.3 Reflective surfaces with a point source . . . 116
6.4 Refractive surfaces with a point source . . . 121
6.5 Comparing with other approaches . . . 124
6.6 Conclusion . . . 125
7 Extended Light Sources 127 7.1 Differences in the case of an extended source . . . 127
7.2 Problem statement . . . 130
7.3 Methods for extended light sources . . . 130
7.4 Initial guess for extended light systems . . . 133
7.5 Conclusion . . . 144
8 Conclusion and Outlook 147 8.1 Conclusion . . . 147
8.2 Further research . . . 149 A Strong boundary conditions for OT 151
Bibliography 155
Chapter 1
Introduction
1.1 Brief history of illumination problems and optimal transport
The shaping of light through lenses and mirrors is a question scientists are interested in for centuries. Greek mathematicians are said to have con- structed parabolic mirrors and lenses focussing light beams in the second century BC. Detailed recording of the design of such parabolic mirrors are given by the Arabian physicist Ibn Sahl in the tenth century [Ras90]. In modern history, since the 19th century light houses have used mirrors in- tensifying light beams in certain directions. In the beginning these were fabricated by piecing together small silver coated facets, but later on also smooth mirrors were manufactured [Soc]. In 1823 Fresnel and Arago in- vented prismatic lenses which replaced the mirrors needing less maintenance [Wat99]. The first success to derive equations governing the reflector shape has been made in the 1930s [Bol32] and shortly after in 1941 the connec- tion between certain reflectors and the Monge–Ampère (MA) equation has been discovered [Kom41]. In the following decades there have been several extensions to include different kind of geometries which are summarised in the book [Elm80], at this point all design techniques were limited to point sources and rotational or translational symmetric systems.
Mathematically illumination problems are strongly connected to the opti- mal transport problem. The original optimal transport (OT) problem asked how to best move given piles of sand to fill up given holes? In economy this translated directly to the question how to efficiently transport goods be- tween different places, for example production plants and costumer, where efficiency is evaluated according to transport cost function. This problem on its own is very interesting and has attracted a lot of attention in the last decades. On the one hand, although it is such an easily explainable problem,
even the existence of solutions poses difficult mathematical challenges due to the nonlinearity in cost functions. Furthermore the problem shows many connections with other branches of mathematics, for example geometric and functional inequalities and partial differential equations (PDEs). But also the problem formulation itself is very close to reality. It can be formulated either discrete or continuous and directly models problems occurring in spa- tial economics, traffic, networks, image processing, geophysics and many more. Nevertheless, after the first optimal transport problem statement from Monge in 1781 [Mon81] the OT problem has been mainly investigated from the theoretical point of view. Only recently, due to new results and contemporary computational possibilities, the research for numerical meth- ods became very active. A particular break through by Brenier [Bre91]
connected the solution of an OT problem to the solution of an equation of MA type and thus, gave rise to solve OT problems with PDE solvers.
How is OT connected to the propagation of light in an illumination sys- tem? The continuous version of an OT problem can be simply transferred to the question, how can light of a light source be transported, i.e. redi- rected, by an optical component, such that it creates a specific pattern on a target? With this basis, the design of a reflector can be attributed to solving PDEs of MA type [Kar10]. While this yields direct methods, the nature of the design problem is that of an inverse problem and it is referred to as the inverse reflector/refractor problem. Its solution in the point source setting is known to be non-unique and is unstable under small deviations of the desired light distributions. The connected forward problem simply is the simulation of light propagation in a lighting system, also known as ray-tracing. With ray tracing techniques at hand, designer design mirrors and refractors often in a trial and error process. Due to the many degrees of freedom in free form surfaces, this optimisation process cannot only rely on computational power, but has to be tweaked manually. Solving related PDEs instead is a huge speed up and can be done without much insight in optics. However, “despite its growing list of applications, and in contrast to its extensive and mature PDE theory, the construction and analysis of com- putational methods [for MA equations] is still a relatively new and emerging field in numerical analysis”[NSZ19] and best suited algorithms have yet to be discovered.
Due to the rise of computational power in the beginning of this century, design of general optical surfaces for extended sources has started. Unfortu- nately, it is not known whether systems with extended sources are governed by certain PDEs. Current industrial design strategies produce optical sur- faces far from optimal which are not only unnecessarily large in term of materials, but contribute to the increasing light pollution. In particular, with the advent of LEDs, which are often clustered and, thus, behaving like
1.2. GOAL AND STRUCTURE OF THE THESIS 3 an extended source, lightning is added to more and more situations and should be tackled efficiently.
Additionally to the theoretical design of optical components the achieve- ments of the manufacturing have been significant. While previously only certain radial symmetric mirrors and lenses could be produced, today, sur- faces represented by NURBS-surfaces can be precisely manufactured. Fur- thermore, LEDs produce not as much waste heat as conventional light bulbs, so today optical components may be made of plastic instead of glass or metal.
This has simplified the manufacturing process, such that today the limit- ing factor for luminaires is not the manufacturing, but the previous design process.
1.2 Goal and structure of the thesis
Although studying illumination systems has a long history, even today many design solutions are tailored to specific situations. A common simplification is to consider only settings of rotational or translational symmetry to ob- tain a two-dimensional problem. Still, the design of optical surfaces often is a slow process, even if surfaces for similar configurations are known, and depends heavily on the adjustments of an expert. Thus, our goal is to find a robust, versatile numerical method which can be applied to a gen- eral illumination problem. Robust means the parameters in the algorithm do not have to be tweaked in a long trial and error process for every new problem instance, but a working configuration is found quickly and mostly automatically. Versatile means the method should work for most illumina- tion systems of practical interest, and especially should work for the most difficult near field case. Due to possible modelling of illumination systems with point sources by PDEs of MA type, we first review the simpler, yet strongly related OT problem. The goal is to develop a stable and general method for this application, always keeping in mind to apply it later to the inverse reflector/refractor problem. This method directly gives rise to a robust point source algorithm, and we want to discuss its potential in the context of extended light sources. The general design practice for extended light sources is a least-squares optimisation. We hope to improve this pro- cess by combining several solutions considering point sources.
Due to these foci the thesis is divided into three major parts: The OT problem and the connected PDE of MA type, the inverse refractor/reflector problem for point sources and illumination systems with extended sources.
The thesis starts in Chapter 2 with the mathematical introduction to OT problems and their connection to PDEs of MA type. We summarise rel- evant results to get familiar with the problem and associated challenges.
Especially, we focus on problems relevant for the design of numerical algo-
rithms and then review existing approaches. In Chapter 3 we develop a new conforming finite element method to solve the PDEs arising in the context of OT problems. The main idea is to transfer a method of solving nonlin- ear second order PDEs to these particular PDE problems. This cannot be done trivially, as the boundary conditions from the OT problem need spe- cial treatment, and the PDEs of MA type are only elliptic and well-posed if restricted to the space of convex functions. Within the derivation we al- ways keep in mind that the solver should also be applicable to a wider range of PDE problems of MA type, in particular those arising in the context of optics. Hence, the method is stabilised to ensure good performance on ex- amples of practical interest. We conclude the part on classical OT problems in Chapter 4 with numerical results of the newly introduced method.
Afterwards, Chapter 5 introduces the inverse reflector/refractor problem model, analyses the similarities to the OT problem and revises solution strategies for the restriction to zero etendue sources. Due to the connec- tion between the problems we are able to apply the algorithm from Chapter 3 to the just introduced illumination problem. We present empirical evi- dence that this method is suitable for illumination problems in Chapter 6.
Although the restriction to zero etendue already provides a lot of mathemat- ical challenges, we discuss the difficulties when introducing extended light sources to illumination systems in Chapter 7. There, we introduce an idea how to improve the initial guess for the standard optimisation method for extended sources, the phase space method. While currently often a point source is used to start the optimisation method with, we merge the outcome of several runs of the point source algorithm, each considering a differently located point source. The driving idea for the merging procedure hereby is a convex combination of optical surfaces based on a unified parametrisation.
The practicability of this approach is verified with a numerical example. We conclude with a short resume of accomplishments of the thesis and give an outlook on further research opportunities in Chapter 8.
Chapter 2
Optimal transport and
partial differential equations of Monge–Ampère type
In this chapter we introduce and characterise the OT problem. We start with its definition, derive the related PDE of MA type and briefly repeat main results on existence, uniqueness and regularity of solutions. In the second part we identify challenges for numerical PDE solvers and review already developed solution procedures with an emphasis on state of the art Galerkin methods.
2.1 The optimal transport problem
2.1.1 The Monge formulation of the optimal transport problem
Originally Monge stated in [Mon81] the following problem. How can a fixed set of wholes be filled with rubble with minimum effort? For the mathe- matical formalisation we represent the rubble and the area of wholes by two densities f :X → R>0, g :Y → R>0 given on convex sets X, Y ⊂R3. An example forf and gwith X=Y = [0,1] are given in Figure 2.1.
Naturally, these densities, describing either the rubble or the moulds, should have the same amount of mass, i.e.
Z
X
f dx=Z
Y
gdy. (2.1)
The solution of the problem is now a distribution of rubble to wholes, i.e.
0 0.2 0.4 0.6 0.8 1 0
0.2 0.4 0.6 0.8 1
(a) densityf, a pile of sand
0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2 0
(b) density−g, a hole Figure 2.1: Example densities for a Monge problem.
an mapτ :X→Y with Z
τ−1(A)f dx=Z
Agdy for any Borel subset A⊂Y, (2.2a) which minimises the transportation costs modelled by
Z
Xf(x)kx−τ(x)k2 dx. (2.2b) Since τ models a transport, this problem is called the optimal transport (OT) problem.
Nowadays the following more general formulation is common [San15, Intro- duction].
2.1.2 The modern formulation of the OT problem
Let X and Y be separable, metric spaces, and P(X) and P(Y) the set of Borel measures onXandY, respectively. The optimal transport problem in the most general form reads as follows: Given two Borel measuresµ∈ P(X) and ν ∈ P(Y) find a transport mapτ :X→Y, such that
ν(A) =µ(τ−1(A)) for allA⊆Y, measurable (2.3) while it minimises the cost
J(τ) =Z c(x, τ(x))dµ(x).
induced by a cost functionc:X×Y →R+.
2.1. THE OPTIMAL TRANSPORT PROBLEM 7 In our application we assume that we have the spaces X and Y given by Ω ⊂ R2 and Σ ⊂ R2, as well as two densities f : Ω → R>0, g : Σ → R>0
withµ(dx) =f(x)dxandν(dy) =g(y)dy. Furthermore we consider the cost functionc to be the squared Euclidean distance function
c(x, y) :=kx−yk22 forx, y∈R2.
Hence, whenever we talk of the solution to an OT problem, the solution has to satisfy
Z
τ−1(A)f dx=Z
A
gdy for any Borel subsetA⊂Σ (2.4a) and
τ = arg min
υ:Ω→Σ
Z
Ωkx−υ(x)k22f(x) dx. (2.4b) This means an OT problem is completely defined by the starting domain Ω, the target area Σ, as well as a density on the start domainf : Ω→Rand a density on the targetg: Σ→R. Note that, sincef and g are densities, the conservation of mass in (2.1) automatically holds.
In the literature, especially in the context of economic problems, one also finds the discrete version. Here we have k different locations containing n different goods represented byxj ∈Rn, j= 1. . . kand we want to transport them according to a distribution yj ∈ Rn, j = 1. . . k. The transport of goods from xi toyj, 1 ≤i, j ≤k costsc(xi,yj). Here, the question for the cheapest delivery scheme, the discrete counter part of τ, arises.
The continuous formulation can be additionally weakened if one asks for a transport plan instead of a map, going back to an idea of Kantorovich. He describes the relaxation in [Kan42] in Russian, an translation to English is given in [Kan06b]. In essence, a plan does not transport the whole mass f(x) from a pointx∈Ω to some spotτ(x)∈Σ, but may split the transport to different places. Vice versa, the amount at y ∈ Σ can originate from different x ∈ Ω. A transportation plan t(x, y) then describes how much from x ∈ Ω is transferred to y ∈ Σ. For a solution to the OT problem it has to be ensured that all mass is moved to the target with minimal effort, i.e. a solution t with the marginals
f(x) =Z
Σt(x, y) dy, g(y) =Z
Ωt(x, y) dx
minimises Z
Ω×Σt(x, y)c(x, y) d(x, y).
2.1.3 Boundary conditions of the optimal transport problem
Note that (2.3) implies the surjectivity of the transport map τ. This may be reduced to the boundary equation [TW09]
τ(∂Ω) =∂Σ
which is known as the second boundary value problem.
We characterise the target boundary ∂Σ implicitly by a signed distance function H. To that end we define for a given Σ ⊂ R2 the function H : R2 →Rby
H(y) =
+dist(y, ∂Σ) = min
z∈∂Σkz−yk2, y /∈Σ
−dist(y, ∂Σ), y ∈Σ. (2.5)
After the introduction of H, we formulate for the OT mapτ the boundary condition
H(τ(x)) = 0 ∀x∈∂Ω. (2.6)
SinceH is an implicit characterisation of the manifold ∂Σ, the derivative of H aty is the normal aty0. In case of a convex boundary ∂Σ as presented in Figure 2.2, we use the representation ofH described in [BFO14; BFO12].
It holds that
H(y) = sup
knk2=1{y·n−H∗(n)} (2.7) withH∗ being the Legendre–Fenchel transformation which simplifies to
H∗(n) = sup
z∈∂Σ{z·n}. (2.8)
due to the convexity of Σ [BFO12, Section 2.1].
y6∈∂Σ
y0= arg min
z∈∂Σ kz−yk2
∂Σ
Figure 2.2: Visualisation of the signed distance functionH.
2.2. PDE OF OPTIMAL TRANSPORT 9 These expressions in (2.8) and (2.7) assume maxima for the boundary point y0 = arg minz∈∂Σkz−yk2andnybeing the normal of∂Σ aty0, respectively.
Hence, we obtain for the derivative
∇yH(y) = arg max
knk2=1 {y·n−H∗(n)}=ny. (2.9)
2.2 Partial differential equation of optimal transport
Condition (2.3) asserts thatτ ensures a mass conservation for any setA⊂Ω with respect to the densities f and g. One has to keep in mind that this does not imply that the density f and the concatenation of the transport map τ and g are equal pointwise, i.e. in general we have for x∈Ω
f(x)6=g(τ(x)).
Instead, for any sufficiently smooth transport map τ we can apply the change of variables on (2.3) and we infinitesimally get the Jacobian equation, cf. [DF14, Section 3] or [Caf03]
f(x) =g(τ(x))|det(∇τ(x))|. (2.10) Brenier showed that the solutionτ of the OT problem with Euclidean cost is always the gradient of a convex functionu[Bre91]. Two nice motivations of the idea of the proof can be found in [Caf03; Eva97], a more formal, recent proof can be found in [DF14, Theorem 3.1]. With (2.10) and (2.5) we have the following problem corresponding to the OT problem: We search for a convex u: Ω→R with∇u=τ : Ω→Σ satisfying
F(u) := detD2u− f
g(∇u) = 0 in Ω, (2.11a)
G(u) :=H(∇u(x)) = 0 for all x∈∂Ω, (2.11b) hui:= 1
|Ω|
Z
Ωudx=um0 . (2.11c)
Here, D2u denotes the Hessian of u and um0 ∈ R is constant. The partial differential equation (2.11a) is of Monge–Ampère (MA) type, characterised by its linearity in the determinant of the Hessian.
A solutionτ =∇uto the OT problem (2.4) is already fixed by the equations (2.11a)-(2.11b). However, they determine the MA solution u only up to an additive constant. To fix the additive constant, we additionally specify the mean value of u by a given value um0 and, hence, introduce the auxiliary constraint (2.11c). In the latter algorithm we find it convenient to set um0
to the mean hu0i of an initial guessu0.
We note that we sketched the connection between the OT problem (2.4) and the MA problem(2.11) only for smooth τ. For a review of the connection between both problems in cases of weak solutions we refer to the survey in [DF14]. Regarding the equivalence of the problems, we have that the gradient of any convex solution uof (2.11) yields an OT map solving (2.4), whereas the converse is not true. There exist f, g,Ω,Σ such that problem (2.4) is solvable, while (2.11) cannot be solved.
2.2.1 Existence, uniqueness and regularity
We will focus in this thesis on how to solve the MA problem (2.11) numer- ically. Nevertheless, to understand arising problems, we briefly summarise some theoretical aspects. Considering the underlying OT problem (2.4) a detailed overview is given in [Vil09], or specially written for applied math- ematicians in [San15]. While the problem formulation goes back to Monge, for a long time only certain facts about a solution were known. The exis- tence and regularity remained open. A first breakthrough was the relaxation to transport plans [Kan42](translated to English in [Kan06b]) and the addi- tional introduction of the dual problem as known from linear programming [Kan48] (translated to English in [Kan06a]). In the first instance, this set- tled the question of uniqueness and existence for the relaxed problem. For- tunately, the answer extends in the case of quadratic cost, cf. (2.4b): Brenier showed that the optimal transportation plan does not split any mass and, thus, is a valid solution to (2.4) [Bre91]. The only demand for the existence of such a unique solution map τ is that f is absolutely continuous with respect to the Lebesgue measure [San15, Theorem 1.17].
Let us now consider the substitute MA problem (2.11). PDEs of MA type, which are the classical representative of fully nonlinear equations, have been an active field of research in the past decades [GT83; Gut01]. The combina- tion with the second boundary value problem, i.e. the boundary condition (2.11b), has drawn less attention, but is studied in [Urb97; TW08; San15].
Compared to the OT problem (2.4), the formulation of the MA problem (2.11) has higher demands on the existence of a solution due to the appear- ance of second derivates. The most complete existence theory for classical solutions is restricted to the elliptic case: The operator F in (2.11a) is re- ferred to be elliptic if its linearisation is elliptic. This holds if the Hessian is positive definite, coinciding with the demand of convexity for solutions in (2.11). We restrict ourselves to the case that Σ is convex. Then the following result has been formally proven by Caffarelli [Caf92; Caf96] and is nicely presented in [DF14, Theorem 3.3]
• If f, g > α for an α >0 andf ∈Clock,β(Ω), as well as g ∈Clock,β(Σ) for
2.2. PDE OF OPTIMAL TRANSPORT 11 someβ ∈(0,1), it holds thatτ =∇u∈Ck+1,β(Ω).
• Furthermore, if f ∈ Ck,β(Ω), g ∈ Ck,β(Σ) and Ω and Σ are smooth and uniformly convex, then τ =∇u∈Ck+1,β(Ω).
This means the transport mapτ =∇uis locally ’one derivative better’ than our initial data f and g. In particular, for continuous, boundedf and g we can expect a solution in a classical sense.
Apart from the classical solution there are several weak solution notions for PDEs of MA type. We sum up only the most popular notions. A solution u is said to satisfy (2.11a)
• in a Brenier sense, if∇usatisfies the corresponding OT problem (2.4).
• in an Alexandrov sense, ifuis convex and for any Borel setA⊂Ω the measure
Mu(A) :=
[
x∈A
∂u(x)
is absolutely continuous and its density equals the right-hand side of (2.11a). Here
∂u(x) :=np∈R2 |u(z)≥u(x) +p·(z−x)∀z∈Ωo
denotes the subdifferential ofu atx∈Ω and|·|the Lebesgue measure of a set.
• in a viscosity sense, if u is a viscosity subsolution and a viscosity supersolution of (2.11a), i.e. if for every x ∈ Ω and every convex functionϕ∈C2(Ω) withu≤(≥)ϕin Ω, there holdsF(ϕ)(x)≤(≥)0.
While the first notion applies only to these PDEs of MA type which arise from the OT context, Alexandrov solutions arise in the context of differen- tial geometry and viscocity solutions can be defined for any fully nonlinear second order PDE [FGN13, Definition 1.1]. For convex Σ the notion of Bre- nier and Alexandrov solutions coincide, for non-convex Σ there exist exam- ples showing that Brenier solutions do not solve the PDE in an Alexandrov sense. As the regularity of these problems is quite delicate and complex, we refer the interested reader to [Gut01; FGN13; DF14]. For the context of the thesis we emphasis that none of these weak solution notions relies on a variational formulation based on partial integration. The nonlinearity of the MA operator prevents that the higher derivatives can be thrown on test functions.
Note that all these results apply to the OT problem with quadratic cost stated in (2.4). For different cost functionsc, the solutionτ can be expressed by means of a gradient of another potential function, and a PDE of MA type similar to (2.11a) derived. The understanding of these PDEs has been greatly advanced in [MTW05].
Before discussing numerical methods, let us have a look at the severe chal- lenges that have to be faced during a solution procedure.
2.3 Numerical challenges
In this section we first enumerate all challenges we need to pay to attention to during algorithm design. In its second part we explicate the impact of the hardest aspect, namely the nonlinearity in detail.
First, the PDE (2.11a) is equipped with nonlinear and non-natural bound- ary conditions (2.11b). The great majority of solvers for equations of MA type are developed for Dirichlet boundary conditions. Second, we know the regularity of the solution depends on the boundedness of f, g from below.
So, in particular we have to monitor the behaviour of the algorithm for small values of f and g. The third challenge is especially important in our latter application on similar MA problems arising in optics. There we of- ten solve the problem on circular domains Σ, such that we are particularly interested in algorithms working on general domains. Next to the aforemen- tioned problems, our main challenge lies in the full nonlinearity of the PDE (2.11a) itself. In general, solution procedures resolve the nonlinearity of an equation into a sequence of linear problems, e.g. by Newton’s method. The question is, how to intertwine this with the discretisation. We take a closer look at this general dilemma.
2.3.1 Concatenation of discretisation and linearisation Either one linearises the PDE first and discretises it afterwards, or one starts by a discretisation and embeds the result in the linearisation process. We will refer to the first strategy as the linearisation-discretisation strategy and, accordingly, to the second as discretisation-linearisation. Both approaches are sketched in Figure 2.3. For a more exemplary look, we consider a finite element (FE) discretisation and the linearisation as performed by a plain (idealised) Newton’s method. Furthermore we assume a smooth solution in H2(Ω) exists.
At first, we follow the linearisation-discretisation path. We need to choose a linearisation point ¯u, which directly yields the question of the space that ¯u originates from. Since we restrict ourselves to smooth solutions, we postulate
2.3. NUMERICAL CHALLENGES 13
F(u) = 0
Fh(uh) = 0 discretise
DF[u]w=−F(u) linearise at ¯u
ALD(wh) =lLD discretise
DFh[uh](wh) =−Fh(uh) linearise atuh
¯
u u¯+wh
uh←uh+wh
Figure 2.3: Schematic view on the interplay between linearisation and dis- cretisation.
¯
u∈H2(Ω). The Newton updatew∈H2(Ω) is determined
DF[u]w=−F(¯u) (2.12)
with DF[u]w being the Gateaux derivative of the PDE operator F atu in the direction w. In our case the linearisation of F reads
∇ ·cofD2u¯∇w+ f
g(∇¯u)2∇g(∇¯u)· ∇w, (2.13) where the matrix cof (M) is called the cofactor matrix ofM defined by
cof
m00 m01 m10 m11
!!
:= m11 −m10
−m01 m00
!
. (2.14)
Equation (2.13) is well-known, it is a convection-diffusion equation. Its discretisation has been extensively investigated and FE methods based on C0 orDG elements are applicable [ABC+02; ESW14]. Thus, performing a Galerkin approach we obtain the finite-dimensional problem
ALD(wh) =lLD
searching for wh ∈ VhLD, whereALD :VhLD → VhLD0 and lLD ∈VhLD0 with the Galerkin discretisation space VhLD and its dual space VhLD0. The crux now is the performance of another Newton step. We would usually update (2.12) with ¯u being the intermediate solution, here ¯u ← u¯+wh. But the formulation of (2.13) which is the left-hand side of (2.13) requires the guess ¯u to be twice differentiable. This holds true for the original linearisation point
¯
u in H2(Ω), but not necessarily for intermediate solutions in the Galerkin space VhLD.
On the other hand, in the discretisation-linearisation order one starts by defining a discrete space VhDL and a discrete form Fh :VhDL →Rsuch that the numerical solutionu∈VhDL is determined by
Fh(uh) = 0.
First, it is not obvious how to design Fh. Even if we imagine that we have found a candidate forFh, some aspects indicate severe challenges validating the method. First, we miss the chance of exploiting the (advanced) under- standing of the linear case (2.12). Thus, we have to make the analysis all from scratch, i.e. prove the existence of roots of Fh and that these roots really approximate the solutionu∈H2(Ω). Moreover, it has to be ensured the roots of Fh can be calculated numerically. Again we image that our Fh satisfies all this and moreover, we can find its roots by a plain Newton’s method. In this case, we can update an initial guessuh ∈VhDLby an update wh∈VhDL satisfying
DFh[uh](wh) =−Fh(uh). (2.15) We emphasise that in general (2.15) is a different discrete problem than (2.12) obtained by he linearisation-discretisation path, even if the discreti- sation on both paths is done in the same fashion.
So, we are facing various challenges solving nonlinear PDEs. Recall that we introduced the PDE problem (2.11) as a substitute for the OT problem (2.4). The original problem (2.4) can also be solved directly.
2.4 Methods for optimal transport
Until recently, the analysis of the OT problem (2.4) got more attention than its simulation. The most familiar solution procedure goes back to Kantorovich and formulates the problem as a linear programming problem [Vil09]. This problem can be solved directly by standard linear program- ming techniques such as the simplex algorithm. However, the size of the linear program, i.e. the problem formulation, grows quadratically with the number of unknowns. Thus, solving the problem becomes infeasible for large problems [UR00]. Further methods are summarised in [San15, Chapter 6], most famous being the Benamou-Brenier method [BB00] and the procedure proposed in [HZT+04]. Other approaches reformulate the OT problem as a convex optimisation problem [Mér11; LR17].
Recently, numerical methods solving the substitute MA problem (2.11) with its unusual boundary conditions became popular. Several finite difference (FD) schemes are developed in [LR05; BFO14; BCM16]. They fit in the framework of Barles and Sougadinis which provides a general proof of their
2.5. GALERKIN METHODS FOR PDES OF MA TYPE 15 convergence [BS91]. For some of those FD schemes, the convergence to singular (weak) solutions can additionally be proven [BFO12; BF14]. Yet, the convergence order is only linear, and the construction is not possible for narrow stencils [MW53]. The latter problem can be overcome by so-called wide-stencil finite difference schemes [Obe08]. However, to obtain highly accurate solutions, in every grid refinement the size of the stencil needs to be increased, making them of limited use for high resolution problems [FGN13, Section 2.1.5.]. In [PBt+15] a three step procedure for the PDE problem (2.11) based on the work in [CGS13] is proposed. The idea is to use three kind of degrees of freedoms to detach the interior gradient
∇u, the Hessian D2u and ∇u |∂Σ, i.e. the interior transport mapping, the Jacobian of the transport mapping and boundary information, respectively.
If one kind of degrees of freedom is considered variable, while the others are temporarily fixed, the MA problem (2.11) leads to a simple minimisation problem. The variable degrees are updated accordingly and the process is reiterated with the next kind of degrees of freedom. The procedure is finished, if no significant change is registered any more.
Despite the variety of numerical methods, so far there is no robust, higher- order solver on general domains. Due to the connection to MA problems, recent developments in numerical methods for fully nonlinear PDEs can also be applied to the substitute PDE problem (2.11). For a summary of numerical methods we refer the reader to [FGN13], as we only review Galerkin methods for equations of MA type in the next section.
2.5 Galerkin methods for partial differential equations of Monge–Ampère type
For an overview about relevant methods, we now do not consider (2.11a), but instead the simpler problem of the standard MA equation equipped with Dirichlet boundary conditions, i.e.
detD2u=φin Ω, (2.16a)
u= 0 at ∂Ω (2.16b)
withφ∈L2(Ω).
The canonical variational formulation candidate is the following L2 projec- tion: Find u∈H02(Ω) such that for allv∈L2(Ω) it holds
Z
Ω
detD2u−φvdx= 0. (2.17) This problem indeed is solvable if and only if u is a classical solution to the MA equation (2.16a) [Böh08]. While for linear PDEs the second derivatives
of u could be shifted to the test functions v by integration of parts, this is not longer possible due to the nonlinearity in D2u. Hence, (2.17) is not well-defined for standard finite elements fromC0(Th(Ω)) for a triangulation Th(Ω) of Ω, since they are non-differentiable between neighbouring simplices.
To understand the main ideas behind existing finite element methods, we recall the diagram in Figure 2.3 which is also applicable to the standard MA equation (2.16a). As we already stated in Section 2.3.1, on the linearisation- discretisation path a conflict between possible linearisation points and the discrete candidates arises. On the discretisation-linearisation path however, the linear problem DFh[uh](wh) =−Fh(uh) is, in general, not consistent with the infinite-dimensional, linearised problemDF[u]w=−F(u).
The first finite element method described in [Böh08] suggests a discretisa- tion of the L2 projection (2.17) with an ansatz space contained in C1(Ω).
As such DFh[uh] provides a valid discretisation of the infinite, linearised operator DF[uh] and thus, the diagram actually commutes. We also refer to finite elements with this property as conforming elements, but stress that their implementations is not trivial. We revisit this method in more detail in Section 2.5.1.
The vanishing moment method proposed in [FN09; FN11] tackles arising problems already before the discretisation and linearisation stage. Instead of handling the PDE (2.16a) directly, the method adds a regularisation term to generate a sequence of easy-to-handle higher order PDEs. The discreti- sation for those higher order PDEs may be of Galerkin type, but this is not exclusive.
Recently, further Galerkin methods with non-conforming elements, i.e. finite element spaces not contained in C1(Ω), have been developed in [BGN+11;
LP13; Nei14; NNZ18a]. These non-conforming methods have one idea in common. Although they start by proposing a discrete Fh, i.e. derive the method along the discretisation-linearisation path, their goal is to keep con- sistency to the linearisation-discretisation path. To do so, they have to take care thatDFh[uh]is a valid discretisation of the infinite, linearised operator DF[u]. In [BGN+11] the authors achieve this by introducing additional stabilisation parameters to the nonlinear, discrete Fh. It is known that this method fails to converge except for smooth examples with solution inH2(Ω) [Nei14, p.366]. In [LP13; Nei14] a kind of mixed method is introduced: For an intermediate discrete solutionuk,h, which is not necessarily contained in H2(Ω), the notion of a discrete Hessian as a natural extension of second derivatives is established. Thus, the linearisation step may also be formed based on discrete numerical solutions. Although not proven, numerical re- sults show that this approach also converges in cases where only a weak solution exists [Fri14, Chapter 5]. In [KLP18] this procedure is adapted to fit the OT boundary conditions (2.11b), and thus they solve the system
2.5. GALERKIN METHODS FOR PDES OF MA TYPE 17 (2.11a) and (2.11b). The work in [NNZ18a; NNZ18b] embeds the ideas of FD discretisations from [FO11] in a finite element setting. Whenever the values of the Hessian is required, the FD operator developed in [FO11] is evaluated using the FE solution. The stencil width of the embedded FD operator is chosen differently from the grid width of the FE discretisation, thus, yielding a method on two different scales. Currently, it is the only finite element method proven to converge to weak viscosity solutions [NNZ18b].
While the non-conforming methods are promising, they have some disadvan- tages. The stabilisation parameters from [BGN+11] have to be calculated according to the original PDE operator and, thus, are not easily adapted to the MA problem (2.11). The mixed method introduces additional degrees of freedom for the Hessian, yielding nearly a fivefold increase of the degrees of freedom. Since the method includes a simple assembly process compared to conforming methods, one could argue there is a trade-off between the complexity of the assembly process and the size of the linear system. The method described in [NNZ18a; NNZ18b] is based on linear elements, thus does not exceed first order convergence in case of smooth test problems.
Although the conforming finite element method is inherently connected to classical solutions and not suitable for weak solution concepts, i.e. solutions not in H2(Ω), it is in advantage with its mathematical simplicity and the possibility to extend it to complex PDEs of MA type. In the following section we describe the conforming finite element method in more detail.
2.5.1 Conforming C1 method
The approach of [Böh08] may not only be applied to the standard MA equa- tion but works for a general second order PDE with homogeneous boundary conditions. However, for simplicity we consider problem (2.16) on a polyhe- dral domain Ω. By testing the strong (nonlinear) formulation of the PDE (2.16a) in the L2-sense, we obtain the following variational formulation:
Find uh∈Vh,0 such that Z
Ω
detD2uh−φvh dx= 0 for allvh ∈Vh,0, (2.18) where Vh,0 is a trial space embedding the homogeneous boundary condi- tions. We emphasis that this variational formulation (2.18) does not define any weak solution concept, but requires the same regularity as the original PDE. To avoid wrong associations with the linear case we refer to (2.18) as a strong variational formulation. Hence, (2.18) is only applicable to conform- ing elements, where the trial spaceVh,0 is contained inC1(Ωh). A candidate with the necessary approximation quality isVh,0 being a proper subspace of
Sh1,d :=nf ∈C1(Ωh)|f |T∈Πd forT ∈ Th
o,
where Πd is the space of all polynomials of total degree≤ d. Furthermore Vh,0 must incorporate the Dirichlet boundary conditions, i.e.
Vh,0⊂Sh1,d∩ {f |f |∂Ω= 0}.
To solve (2.18) numerically Böhmer proposes the Newton method. Hence, the convergence of the method depends on the quality of the initial guess.
While (2.18) can be formulated for every nonlinear PDE, optimal error esti- mates are only proven if the original PDE provides existence, local unique- ness of a solution u and that the Gateaux derivative of F is bounded and invertible at u. These optimal convergence rates are not achieved for non- polyhedral domains. To this end, consider the domain Ω to be such that
∂Ω can be described by a function in C2. Here, if we choose elements of Sh1,d which have support on triangular meshes Th instead of Ω, the conver- gence rate is bounded by two. However, the work in [Böh08] focuses on the theoretical point of view in case of classical solutions and lacks of any numerical experiments. An implementation described in [DS13] fails to con- verge for a problem with a solution not in H2(Ω), but verifies the higher order convergence for smooth examples with solutions in H2(Ω).
Chapter 3
A conforming finite element method for partial
differential equations of Monge–Ampère type
After having discussed the general problem and known solution procedures, in this chapter we develop a new method for MA problems (2.11) which is loosely based on the idea of the conformingC1 method described in Section 2.5.1.
3.1 Outline of the derivation
The algorithm consists of a delicate intertwining of procedures each tackling one of the challenges described in Section 2.3. Before going into detail we shortly anticipate the derivation on an abstract level.
First, the core PDE (2.11a) is processed along the linearisation-discretisation path as described in Section 2.3.1. Thus, in the beginning we linearise prob- lem (2.11) with the aid of an idealised Newton’s method. Then we formulate the variational framework for the linear problem. To gain control of the accompanying conditions resulting of (2.11b) and (2.11c) we introduce La- grangian multipliers. Afterwards, we formulate a computable, correspond- ing problem in suitable finite discretisation spaces. In the subsequent part, we tailor several adjustments for real world applications. In particular, we introduce a nested iteration, see [Bra13, §4], also sometimes called multires- olution technique or multilevel approach. This nested approach ensures that the initial guess needs to be of good quality only on a coarse resolution. Ad- ditionally, we regularise data in several places to speed up the convergence
process, as detailed in Section 3.4.1. Although this latter subtle changes may not be in the spirit of an irregular weak solution, the algorithm pro- duces satisfactory results in empirical examples of interest. Interestingly, in the related optics context this regularisation is even beneficial for the man- ufacturing process. At the end of this chapter we conclude with remarks on the actual implementation.
3.2 Infinite-dimensional formulation
We tackle the nonlinearity in F and G in problem (2.11) by Newton’s method.
3.2.1 A Newton step
Let V :=H2(Ω). Given that a Newton update wfor an initial guessu∈V exists, this update whas to formally satisfy
DF[u]w=−F(u) in Ω, (3.1a)
DG[u]w=−G(u) on ∂Ω, (3.1b)
hu+wi=um0 . (3.1c)
Here, the left-hand sides are determined by the Gateaux-derivatives of F andGevaluated atuin the direction ofw. ForF a formal derivation yields the Gateaux-derivative
DF[u]w=∇ ·cofD2u∇w+ f
g(∇u)2∇g(∇u)· ∇w. (3.2) Thus, we have in (3.1a) a convection diffusion equation and we denote the convection term by
b(u) := f
g(∇u)2∇g(∇u). (3.3)
Concerning (3.1b), we have
DG[u]w=∇yH(∇u)· ∇w.
For the gradient of H, we shortly recall the results of Section 2.1.3. If we evaluate∇yH at∇u(x) forx∈Ω, we know that its value corresponds to the normal projecting the transported point∇u(x) to the boundary∂Σ, see also Figure 2.2. More formally, ∇yH(∇u(x)) is given by the normalny(∇u)(x) which is the normal of ∂Σ aty0 ∈∂Σ. Here,y0 is the closest toy =∇u(x) on the desired boundary∂Σ. Hence, we obtain
DG[u]w=ny(∇u)· ∇w.
3.2. INFINITE-DIMENSIONAL FORMULATION 21 Hence, we observe that (3.1b) states non-natural Neumann boundary condi- tions. These conditions do not fix a solution uniquely such that the auxiliary condition (3.1) is necessary.
3.2.2 Weak formulation
We construct a weak solution to (3.1a) in analogy to the classical approach for the Dirichlet problem for second order elliptic equations. For that end, we introduce the subspace (depending onu)
W0,uBC :={z∈H1(Ω)|ny(∇u)· ∇z= 0 on ∂Ω and hzi= 0} (3.4) of H1(Ω) whose elements have zero mean and satisfy (3.1b). It is known that the H1-semi norm is a norm on W0,uBC which is equivalent to the H1- norm. Now consider the bilinear and linear forms aF : W ×W → R and lF :W →Rdefined by
aF(w, v) :=Z
Ω−hcofD2u∇wi· ∇v+ (b(u)· ∇w)vdx +Z
∂Ω
hcofD2u∇wi·nvds (3.5) lF(v) :=−
Z
ΩF(u)vdx=Z
Ω
−detD2u+ f g(∇u)
vdx, (3.6) with b(u) from (3.3). Now pick any fixed wg ∈H1(Ω) satisfying
ny(∇u)· ∇wg =−G(∇u) and hwgi=um0 − hui. Then, findw0 ∈W0,uBC such that
aF(w0, v) =lF(v)−a(wg, v) ∀v∈W0,uBC. (3.7) By the above comments, Lax–Milgram’s Theorem applies so that (3.7) has a unique solution. The Newton update is then given by
w=w0+wg. Remark 3.1
During the so far considered OT context, it is possible to rearrange the boundary conditions (3.1b) to describe natural Neumann boundary condi- tions fitting to classical FE methods. Unfortunately, we also want to consider PDE problems arising in optics which do not easily admit such reformula- tions. For the interested reader we show a derivation for natural boundary conditions in Appendix A, but in later numerical experiments we always enforce them weakly by Lagrangian multiplier as described in the remainder of this Section.
3.2.3 Embedding the boundary conditions
Discretising the lifting process in straight forward analogy to Dirichlet prob- lems is not clear. Instead of constructing first wg in (3.7) we describe next how to append the boundary condition (3.1b) weakly. We only require the Newton update wto satisfy
aG(w, q) =lG(q) ∀q ∈Q (3.8) with the spaceQ:=H1/2(∂Ω) and the bilinear formaG:V×Q→Rdefined by
aG(w, q) :=Z
∂Ω(ny(∇u)· ∇w)qds, (3.9) and the linear formlG:Q→Rdefined by
lG(q) :=− Z
∂ΩG(u)q ds=− Z
∂ΩH(∇u)q ds. (3.10) Now, we follow the strategy of [FB91; DK01] and consider the problem:
Find w∈W and the Lagrangian multipliers (p, λ)∈Q×R such that aF(w, v) +aG(v, p) +λhvi=lF(v) ∀v∈W, (3.11a)
aG(w, q) =lG(q) ∀q∈Q, (3.11b) hwi=um0 − h¯ui. (3.11c) We define the combined bilinear formaGM :W ×(Q×R)→R by
aGM(v,(p, λ)) =aG(v, p) +λhvi.
We know from the Ladyzhenskaya–Babuška–Brezzi theorem, see, e.g., [FB91], that the system (3.11) is uniquely solvable, if aF and aGM are continuous, aF is coercive on the kernel of aGM, and aGM satisfies the LBB condition.
Given w ∈ W which satisfies (3.11), we observe for any v0 ∈ W0,uBC that (3.11a) simplifies to
aF(w, v0) =lF(v0).
Thus, we obtain thatwalso satisfies (3.1a) weakly. Since (3.11b) and (3.11c) coincide with the weak version of the boundary conditions (3.1b) and the auxiliary mean fixing condition (3.1c), respectively, we obtain that the so- lutionw is our desired Newton update.
In (3.11) we have introduced Lagrangian multipliers for both the boundary conditions (3.8) and the mean condition (3.1c). However, the introduction
