MADRID LECTURE # 5

1. MADRID LECTURE # 5 Numerical solution of Eikonal equations

2. 1. Historical Background. Motivated by the analysis of nonlinear models from micro-Mechanics and Material Science, B. Dacorogna (EPFL-Math) was lead to investigate the properties of the solutions (including the non-existence of smooth solutions) of the following system of Eikonal-like equations: Find u H10(Ω) such that (EIK-L) |∂u/∂xi| = 1a.e. on Ω,  i = 1, …, d, with Ωa bounded domain of Rdof boundary Γ.

3. Since problem (EIK-L) has a infinity of solutions we are going to focus on those solutions which maximize (or “nearly” maximize) the functional v →∫Ωvdx. This leads us to consider the following problem from Calculus of Variations: u E, (EIK-L-MAX) ∫Ω udx≥∫Ωvdx, v E, with E = {v|v H10(Ω), |∂v/∂xi| = 1 a.e. on Ω,  i = 1, …, d }.

4. Some exact solutions (useful for numerical validation): Suppose that Ω= {x| x = {x1,x2}, |x1x2| < 1}; the unique solution of (EIK-L-MAX) is given on Ω by: u(x) = 1 – |x1| – |x2|. In general, (EIK-L-MAX) has no solution. Assuming howe er that such u exists, we can easily show that u ≥ 0onΩΓ. Morever since |v|2=d inE, problem (EIK-L-MAX) is equivalent to:

5. (UP)u  E+; J(u)  J(v), v E+, with E+ = {v| v E, v ≥ 0 in Ω} and J(v) = ½ ∫Ω|v|2dx – C ∫Ωvdx, Cbeing an arbitrary positive constant. The iterative solution of (UP) will be discussed in the following sections. Our methodology will be elliptic in nature, unlike the one developed by PA Gremaud for the

6. solution of (EIK-L) . Gremaud methology is hyperbolic in nature and relies on TVD schemes. PAG considers (EIK-L) as a kind of Hamilton-Jacobi equation. 2. Penalty/Regularization of (UP). Motivated by Ginzburg-Landau equation, we are going to treat,  i = 1, …, d, the relations |∂v/∂xi| = 1 by exterior penalty. Moreover, in order to control mesh-related oscillations, we are going to bound||2v||L2(Ω) (without this additional constraint our method did not work; PAG did something like that, too). This leads to approximate the functional J in (UP) by the following Jε:

7. Jε (v) = ½ε1∫Ω|Δv|2dx + J(v) + ¼ ε2– 1Σdi= 1∫Ω(|∂v/∂xi|2–1)2dx; above ε={ε1,ε2 }withε1,ε2both positiveandsmall. Then, we approximate (UP) by (UP)εdefined as follows:

8. The problem (UP)ε uε K+ H2(Ω), (UP)ε Jε(uε)  Jε (v),   K+H2(Ω), withK+ = {v| v  H10(Ω), v ≥ 0 in Ω}. (UP)εis a ‘beautiful’ obstacle problem. Proving the existence of solutions by compactness methods is easy, the real difficulty being the actual computation of the solutions. r

9. 3. An equivalent formulation of (UP)ε Let us denote (L2(Ω))d by Λanduεby pε; problem (UP)ε is clearly equivalent to: pε Λ , (UP-E)ε jε(pε)  jε(q),  q Λ , with q = {qi}dI = 1and,  q Λ ,

10. jε(q) = ½ ∫Ω|q|2dx– C ∫Ω1.qdx + ¼ ε2– 1Σdi= 1∫Ω(|qi|2–1)2dx + I+(q); here:  1 is the unique solution in H10(Ω) of – Δ1= 1 inΩ, 1= 0 on Γ.  I+(.)is defined by ½ ε1∫Ω |.q|2dx if q  (K+H2(Ω)), I+(q) = +∞ elsewhere.

11. I+(.) is convex, proper, l.s.c. overthe spaceΛ. The Euler-Lagrangeequation of (UP-E)εreads as follows (after dropping most εs): ∫Ω p.q dx+ ε2– 1Σdi = 1∫Ω(|pi|2–1)pi qidx + <∂I+(p),q> = (E-L) C ∫Ω1.qdx, q  Λ ; p  Λ .

12. To (E-L) we associate the following flow to capture asymptotically solutions of (EL) p(0) = p0, ∫Ω(∂p/∂t).q dx+ ∫Ω p.q dx+ ε2– 1Σdi = 1∫Ω(|pi|2–1)pi qidx (E-L-F) + <∂I+(p),q> = C ∫Ω1.qdx, q  Λ ; p(t)  Λ , t  (0,+∞). The structure of (E-L-F) strongly calls for an Operator – Splitting solution; motivated by its simplicity*, we will apply theMarchuk-Yanenko scheme to the solution of (E-L-F).

13. Operator-splitting solution of (E-L-F)(I) (0) p0 = p0( = 0, or1, for exemple); then, forn ≥ 0, pn→pn + ½ →pn + 1 as follows: (1/Δt)∫Ω(pn + ½ – pn).q dx+ ∫Ω pn + ½ .q dx+ (1)ε2– 1Σdi = 1∫Ω(|pin + ½ |2–1)pin + ½qidx = C ∫Ω1.qdx,  q  Λ; pn + ½  Λ ,

14. Operator-splitting solution of (E-L-F)(II) (1/Δt)∫Ω(pn +1 – pn + ½).q dx+ + <∂I+(pn +1 ),q> = 0, (2)  q  Λ , pn +1. What about the solution of (1) and (2)? (i) Problem (1) has a unique solution “as soon” as: Δtε2 Problem (1) can be solvedpointwise; indeed,

15. pin+ ½(x) is then the unique solution of a single variablecubic equation of the following type: (1 – Δtε2– 1 + Δt)z + Δtε2– 1z3 = RHS whoseNewton’ssolution is trivial. (ii) The solution of (ii) is given by pn+1 = un+1 where pn+1 is the solution of the following obstacle type variational inequality:

16. (EVI)un +1K+H2(Ω), ∫Ωun +1.(v – un +1)dx + Δtε1 ∫Ω2un +12(v – un +1)dx ≥ ∫Ωpn+ ½.(v – un +1)dx ,  v K+H2(Ω). In order to facilitate the numerical solution of (EVI)we perform a ‘variational crime’ by ‘approximately’ factoring (EVI) as follows:

17. ωn +1K+, (P1) ∫Ωωn +1.(v – ωn +1)dx ≥ ∫Ωpn+ ½.(v – ωn +1)dx ,  v K+ followed by (P2) un +1 – Δtε12un +1 =ωn +1 in Ω, un +1 = 0 on Γ. The maximum principle implies the ≥ 0 of un +1( H2(Ω)).

18. The finite element implementation is easy, the main requirement being that the mesh preserves the maximum principle at the discrete level. Of course, everything we say applies to: (i) Non-homogeneous boundary conditions. (ii) The ‘true’Eikonal equation |u| = f (≥ 0).

19. 4. Numerical experiments Test problems 1 & 2: Ω = (0,1) × (0,1) |∂u/∂x1| = |∂u/∂x2| = 1inΩ, u = gon Γ, with g(x1,0) =g(x1,1) = min(x1,1 – x1), 0  x1  1, g(0, x2) =g(1, x2) = min(x2,1 – x2), 0  x 2 1. We have then: umax(x) = min(x1,1 – x1) + min(x2,1 – x2), umin(x) = |x1 – x2| or |1 – (x1 – x2)| depending wherex is in Ω. C = 10 or – 10, h = 1/128, Δtε1= h2/36, ε2 = 10– 3, Δt =10– 4

24. ||uh– u||2≈ 10–3, ||uh– u||∞≈ 10–2 Test problem 3: Ω = (0,1) × (0,1) |∂u/∂x1| = |∂u/∂x2| = 1inΩ, u = 0 on Γ, with C = 10, h = 1/512 and 1/1024, Δtε1= h2/9, ε2 = 10– 3, Δt =10– 4 This problem has no solution; the weak limit should be the distance function. The approximate solutions exhibit self- -similarmulti-scale structures (fractals)?

29. Solution of the true Eikonal equation Test problem 4:Ω = (0,1) × (0,1) |u| = 1inΩ, u = 0 on Γ. The maximal solution is thedistance to the boundary function Parameters:C = 10, h = 1/256, Δtε1= h2/49, ε2 = 10– 3, Δt =10– 4

32. Solution of the true Eikonal equation Test problem 5 :Ω = (0,1) × (0,1) |u| = 1inΩ, u = gon Γ, with g(x1,0) = 0, g(x1,1) = min(x1,1 – x1), 0  x1  1, g(0, x2) =g(1, x2) = 0, 0  x 2 1. Parameters:C = 10, h = 1/256, Δtε1= h2/49, ε2 = 10– 3, Δt =10– 4