Numerical Solution of a Non-Smooth Eigenvalue Problem

1. Numerical Solution of a Non-Smooth Eigenvalue Problem An Operator-Splitting Approach A. Caboussat & R. Glowinski

2. 1. Formulation. Motivation Our main objective is the numerical solution of the following problem from Calculus of Variations Compute γ = inf v Σ∫Ω|v|dx(NSEVP) where: Ω is a bounded domain of R2and Σ = {v| v  H01(Ω), ∫Ω|v|2dx = 1}.

3. Actually, γ = 2√ π, independently of the shape and size of Ω(holds even for non-simply connectedΩand in fact for unboundedΩ) (G. Talenti). A natural question is then: Why solve numerically a problem whose exact solution is known ? (i) If I claim that it is a new method to compute πnobody will believe me. (ii) (NSEVP) is a fun problem to test solution methods for non-smooth & non-convexoptimization problems.

4. (iii) ∫Ω|v|dxarises in a variety of problems from Image Processing and Plasticity. Actually, our motivation for investigating (NSEVP)arises from the following problem from visco-plasticity : u L2(0,T; H01(Ω)) C0([0,T ]; L2(Ω)); u(0) = u0, (BFP) ρ∫Ω(∂u/∂t)(t)(v – u(t))dx + μ∫Ωu(t).(v – u(t))dx + g[ ∫Ω|v|dx– ∫Ω|u(t)|dx] ≥ C(t)∫Ω(v – u(t))dx, v  H01(Ω), a.e. t (0, T), withρ > 0, μ> 0, g > 0, Ωa bounded domain ofR2and u0  L2(Ω).

5. (BFP) models the flow of a Bingham visco-plastic fluid in an infinitely long cylinder of cross sectionΩ, Cbeing the pressure drop per unit length. Suppose that C = 0 and that T = +∞; we can show that (C-O.PR) u(t) = 0,  t ≥ Tc, with Tc = (ρ/μλ0)ln[1 + (μλ0/γg)||u0||L2(Ω)], λ0being the smallest eigenvalue of – 2 in H01(Ω). A similar cut-off property holds if after space discretization we use the backward Euler scheme for the time discretization of (BFP), with λ0 and γreplaced by their discrete analoguesλ0hand γh.

6. Suppose that the space discretization is achieved via C0-piecewise linear finite element approximations, we have then |λ0h– λ0| = O(h2). But what can we say about |γh– γ| ? The main goal of this lecture is to look for answers to the above question !

7. 2. Some regularization procedures There are several ways to approximate (NSEVP) – at the continuous level – by a better posed and/or smoother variational problem. The most obvious candidate is clearly γε = infv Σ∫Ω(|v|2 + ε2)½dx, (NSEVP.1)ε a regularization quite popular in Image Processing. Assuming that the above problem has a minimizeruε, this minimizer verifies the following Euler-Lagrange equation (reminiscent of the mean curvature equation):

8. First regularized problem:

9. (RP.1)is clearly a nonlinear eigenvalue problem for a close variant of the mean curvature operator, the eigenva lue being γε. Another regularization, more sophisticated in some sense, since this time the regularized problem has minimizers, is provided (with ε > 0) by γε = min v Σ [ ½ ε∫Ω|v|2dx +∫Ω|v|dx ]. (NSEVP.2)ε An associated Euler-Lagrange (multivalued) equation reads as follows, also of the nonlinear (in fact, non- smooth) eigenvalue type (as above the eigenvalue is γε):

10. – ε2uε+∂j(uε)γεuεin Ω, (RP.2)uε= 0 on ∂Ω, ∫Ω|uε|2dx = 1; in (RP.2), ∂j(uε) is thesubgradientatuεof the functional j : H01(Ω) → R defined by j(v) = ∫Ω|v|dx. The solution of problems such as (RP.2) is discussed in GKM (2007); the method used in the above reference is of the operator-splitting/inverse powermethod type.

11. In order to avoid handling simultaneously two small parameters, namely ε and h, we will address the solution of γ = inf v Σ∫Ω|v|dx without using any regularization (unless we consider the space approximation as a kind of regularization, that it is indeed).

12. 3. Finite Element Approximation (i) First, we introduce a family {Ωh}hof polygonalapproxi-mations of Ω, such that limh→0Ωh= Ω. (ii) With each Ωhwe associate a triangulationTh verifying the usual assumptions of: (a) compatibility between triangles, and (b) regularity. (iii) With each Th we associate the finite dimensional space V0hdefined (classically) as follows:

13. V0h ={v| v C0(Ωh∂Ωh), v|T P1, T  Th, v = 0 on ∂Ωh}. (iv) We approximate γ = inf v Σ∫Ω|v|dx(NSEVP) by γh = min v Σh∫Ωh|v|dx(NSEVP)h

14. with Σh= {v| v V0h, ||v||L2(Ωh) = 1}. It is easy to prove that: (i) Problem (NSEVP)h has a solution, that is there exists uh Σhsuch that ∫Ωh|uh|dx= γh. (ii) limh→0γh = γ( = 2√π). We would like to investigate (computationally) the order of the convergence of γh to γ. From the non-smoothness of the problem, we do not expect O(h2).

15. 4. An iterative method for the solutionof (NSEVP)h We are going to look for robustness, modularity and simplicity of programming instead of performance measured in number of elementary operations (this is not image processing and/or real time). At ADI 50 ( December 2005, at Rice University), we showed that the inverse power method for eigenvalue computations has an operator-splitting interpretation; we also showed the equivalence between some augmented Lagrangian algorithms and ADI methods such as Douglas- Rachford’s and Peaceman-Rachford’s. For problem (NSEVP)h we think that it is simpler to take the ALapproa- ch, keeping in mind that it will lead to a ‘disguised’ ADI method.

16. For formalism simplicity, we will use the continuous problem notation. We observe that there is equivalence between γ = inf v Σ∫Ω|v|dx and γ = inf {v, q, z}  E∫Ω|q|dx, where E = {{v, q, z}| v  H01(Ω), q  (L2(Ω))2, z L2(Ω), v – q = 0, v – z= 0, ||z||L2(Ω) = 1}.

17. The above equivalence suggests introducing the following augmented Lagrangian functional Lr : (H01(Ω)×Q×L2(Ω))×(Q×L2(Ω)) → R defined as follows, with Q = (L2(Ω))2 and r = {r1, r2}, ri > 0, Lr(v, q, z; μ1, μ2) = ∫Ω|q|dx + ½ r1 ∫Ω|v – q|2dx +½ r2 ∫Ω|v – z|2dx + ∫Ω(v – q).μ1dx + ∫Ω(v – z)μ2dx

18. We consider then, the following saddle-point problem Find{{u, p, y}, {λ1, λ2}}  (H01(Ω)×Q×S)×(Q×L2(Ω)) such that Lr(u, p, y; μ1, μ2) ≤Lr(u, p, y; λ1, λ 2) ≤Lr(v, q, z; λ1, λ 2), (SDP-P)  {{v, q, z}, {μ1, μ2}} (H01(Ω)×Q×S)×(Q×L2(Ω)), with S = {z| z L2(Ω), ||z||L2(Ω) = 1}. Suppose that the above saddle-point problem has a solut ion. We have then p = u, y = u, ubeing a minimizer for the original mimimization problem (the primal one).

19. To solve the above saddle-point problem, we advocate the algorithm below (called ALG 2 by some practitioners (BB)): (1) {u –1, {λ10, λ20}} is given in H01(Ω)×(Q×L2(Ω)); for n ≥ 0, assuming that{un – 1, {λ1n, λ20}} is known, solve: (2) {pn, yn} = arg min{q, z} Q×SLr(un – 1, q, z; λ1n, λ 2n), then (3)un =arg minvLr(v, pn, yn; λ1n, λ 2n), v  H01(Ω), (4) λ1n+1 = λ1n+ r1(un– pn), λ2n+1 = λ2n+ r2(un– yn).

20. The above algorithm is easy to implement since: (i) Problem (3) is equivalent to the following linearvariational problem in H01(Ω) un H01(Ω), r1∫Ωun.v dx + r2 ∫Ωunv dx = ∫Ω(r1pn – λ1n ).v dx + ∫Ω(r2yn – λ2n )v dx,  vH01(Ω). The solution of the discrete analogue of the above problem is a simple task nowadays.

21. (ii) Problem (2) decouples as (a) pn = arg minq Q[½ r1∫Ω|q|2dx+ ∫Ω|q|dx –∫Ω(r1un+ λ1n).qdx ]. (b) yn = arg minz S[½ r2∫Ω|z|2dx–∫Ω(r2un+ λ2n)zdx ]. Both problems have closed form solutions; indeed, since ||z||L2(Ω) = 1,  z S, one has yn = (r2un+ λ2n) / ||r2un+ λ2n ||L2(Ω).

22. Similarly, the minimization problem in (a)can be solved point-wise (one such elementary problem for each triangle of Th, in practice). We obtain then, a.e. on Ω, pn(x) = (1/r1) (1 – 1/|Xn(x)|)+ Xn(x), where Xn(x) = r1un(x)+ λ1n(x).

23. 5. Numerical experiments First Test Problem:Ωis the unit disk

25. Unit Disk Test Problem Variation ofγhversush

26. Unit Disk Test Problem Variation ofγh– γversush

27. Unit Disk Test Problem Visualisation of the coarse mesh solution

28. Unit Disk Test Problem Visualisation of the fine mesh solution

29. Unit Disk Test Problem Coarse mesh solution contours

30. Unit Disk Test Problem Fine mesh solution contours

31. Unit Disk Test Problem Fine mesh solution contours (details)

32. Second Test Problem:Ωis the unit square Coarse mesh

33. Unit Square Test Problem Fine mesh

34. Unit Square Test Problem Variation ofγhversush

35. Unit Square Test Problem Variation ofγh– γversush

36. Unit Square Test Problem Visualisation of the coarse mesh solution

37. Unit Square Test Problem Visualisation of the fine mesh solution

38. Unit Square Test Problem Contours of the coarse mesh solution

39. Unit Square Test Problem Contours of the fine mesh solution

40. Unit Square Test Problem Contours of the fine mesh solution (details)

41. Circular Ring Test Problem (coarse mesh)

42. Circular Ring Test Problem (fine mesh)

43. A GENERALIZATION Compute for Ω R2 γ* = infv ∫Ω|v|dx with = {v| v (H10(Ω))2, ||v||(L2(Ω))2 = 1}.

44. Conjecture (unless it is a classical result):

45. Square(coarse mesh)

46. Square (fine mesh)

47. Disk(coarse mesh)

48. Disk(fine mesh)

49. The results of our numerical computations suggest very strongly that the value we conjectu- red for γ* is the good one.

50. APPLICATION to a SEDIMENTATION PROBLEM The following problem has been considered by C.Evans & L. Prigozhin u/t +IK(u)  fin Ω × (0, T), (SP) u(0) = u0, with Ω R2 and K = {v | v  H1(Ω), |v|  C, v = g on Γ0 (  Ω)}.