Numerical solution of some geometric inverse problems

Numerical solution of some geometric inverse problems Anna DOUBOVA Dpto. E.D.A.N. - Univ. of Sevilla joint work with E. FERN ´ ANDEZ-CARA - Dpto. E.D.A.N. - Univ. of Sevilla Numerical Resolution for Inverse Problems BCAM, Bilbao, 8-9 January 2015 A. Doubova Numerical solution of some geometric inverse problemsm Outline 1Motivation: Elastography 2Wave equation 3Elasticity systems 4Reconstruction and numerical algorithms 5Numerical results 2D wave equation 2D Lam´ e system 3D wave equation 3D Lam´ e system 6Work in progress A. Doubova Numerical solution of some geometric inverse problemsm Geometric Inverse Problems governed by PDEs Motivation: Elastography We consider: Geometric inverse problems Wave equation and Lam´ e systems Motivation: Elastrography A non-invasive method of tumor detection: when a mechanical compression or vibration is applied, the tumor deforms less than the surrounding tissue A technique to detect elastic properties of tissue from acoustic wave generators (applications in Medicine) A. Doubova Numerical solution of some geometric inverse problemsm Motivation Elastography Classical detection methods in mammography: Figure: Palpation Figure: x-rays A. Doubova Numerical solution of some geometric inverse problemsm Motivation Elastography Elastography (“imaging palpation”) is better suited than palpation and x-rays techniques: — Tumors can be far from the surface —orsmall — or may have properties indistinguishable through palpation or x-rays Figure: Stiffness is represented by a color spectrum, ranging from dark red (very stiff) through orange, yellow, and green, to blue (very soft). A. Doubova Numerical solution of some geometric inverse problemsm Wave equation N-dimensional wave equation (N=2 or 3) (a) Direct problem: Data: Ω,T>0, ϕ,Dand γ⊂∂Ω Result: the solution u (1)       utt −∆u=0 in (Ω \D)×(0,T) u=ϕon (∂Ω) ×(0,T) u=0 on (∂D)×(0,T) u(x,0) = u0,ut(x,0) = u1in Ω Information: (2)α=∂u ∂non γ×(0,T) (b) Inverse problem: (Partial) data: Ω,T,ϕand γ⊂∂Ω (Additional) information: α Goal: Find Dsuch that the solution to (1)satisfies (2) A. Doubova Numerical solution of some geometric inverse problemsm N-dimensional wave equation Uniqueness        ui tt −∆ui=0 in Ω\Di×(0,T),i=0,1 ui=ϕin ∂Ω×(0,T) ui=0 in ∂Di×(0,T) ui(x,0) = 0,ui t(x,0) = 0 in Ω\Di Theorem T>T∗(Ω, γ),D0,D1are convex, ϕ6=0 ∂u0 ∂n=∂u1 ∂non γ×(0,T)    =⇒D0=D1 Fundamental results: H¨ ormander, Lions Attention: Weaker than the geometric condition (Only uniqueness, not observability!) A. Doubova Numerical solution of some geometric inverse problemsm Elasticity systems Isotropic case        utt − ∇ · (µ(x)(∇u+∇ut) + λ(x)(∇ · u)Id.) = 0 in Ω\D×(0,T) u=ϕon ∂Ω×(0,T) u=0 on ∂D×(0,T) u(0) = u0,ut(0) = u1in Ω\D Observation: σ(u)·n:= µ(x)(∇u+∇ut) + λ(x)(∇ · u)Id.·non γ×(0,T) Explanations: u= (u1,u2,u3)is the displacement vector Small displacements. Hence, linear elasticity Isotropy assumptions. The tissue is described by λand µ A. Doubova Numerical solution of some geometric inverse problemsm N-dimensional elasticity system Uniqueness        ui tt − ∇ · (µ(x)(∇ui+∇)(ui)t) + λ(x)(∇ · ui)Id.) = 0 in Ω\Di×(0,T) ui=ϕon ∂Ω×(0,T) ui=0 on ∂Di×(0,T) ui(0) = 0,ui t(0) = 0 in Ω\Di Theorem (Constant coefficients) T>T∗(Ω, γ),D0,D1are convex, ϕ6=0 σ(u0)·n=σ(u1)·n on γ×(0,T))=⇒D0=D1 For uniqueness, the key point is Unique continuation property (Imanuvilov–Yamamoto, 2008, complex conditions on µ, λ) AD, E. Fern´ andez-Cara, work in progress (∃Other unique continuation results for stationary problems: Lin–Wang, 2005; Escauriaza, 2005; Alessandrini–Morasi, 2001; Nakamura–Wang, 2006; Imanuvilov–Yamamoto, 2012) A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 2-D wave equation I Case of a ball Test 1:T=5, u0=10x,u1=0, ϕ=10x x0des = -3,y0des = 0,rdes = 0.4 x0ini = 0,y0ini = 0,rini = 0.6 NLopt (AUGLAG + CRS2), NoIter = 1007, FreeFem++: x0cal = -2.998645439,y0cal = 0.000425214708 rcal = 0.4001667063 Figure: Initial mesh: triangles 992, vertices 526 Figure: The desired center and radius of the ball A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 2-D wave equation II Case of a ball Figure: Computed center and radius: AUGALG + CRS2 x0cal= -2.998645439,y0cal = 0.000425214708,rcal=0.4001667063 A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 2-D wave equation III Case of a ball 0 50 100 150 200 250 300 350 400 450 500 0 500 1000 1500 2000 2500 3000 3500 Iterations Current Function Values Function value Figure: Evolution of Jduring the first 500 iterations of CRS2 A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 2-D wave equation IV Case of a ball 0.2 0.3 0.4 0.5 0.6 0.7 0.8 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 r J(x0,r) Figure: The functional Jwith respect to the variable r A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 2-D wave equation V Case of a ball −5 −4 −3 −2 −1 −2 −1 0 1 2 0 20 40 60 80 100 120 Figure: The functional Jwith respect to the variables x0and y0 A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 2-D wave equation VI Case of a ball Test 2:T=5, u0=10x,u1=0, ϕ=10x x0des = -3,y0des = 0,rdes = 0.4 x0ini = 0,y0ini = 0,rini = 0.6 NLopt (AUGLAG + DIRECT), NoIter = 1001, FreeFem++: x0cal = -2.962962963 y0cal = -0.01219326322 rcal = 0.4220164609 A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 2-D wave equation VII Case of a ball 20 ANNA DOUBOVA AND ENRIQUE FERN ´ ANDEZ-CARA Figure 14. Some experiences with AUGLAG and DIRECTNoScal as the subsidiary optimization algorithm. Remark 3. As in the one-dimensional case, it is possible to introduce another method for the partial identification on B(x0;r) relying on di↵erentiation with respect to the domain. Thus, for each sufficiently small m2R2⇥R, let us assume that ↵=@u @n|⇥(0,T )and ↵m=@um @n|⇥(0,T )are known, with m=(d, s),D+m=B(x0+d;r+s). Here, u(reps. um) is the solution to (35) corresponding to B=B(x0;r)(resp. B=B(x0+d;r+s)). Using domain variation techniques, we obtain the following formula: ↵m↵=L(d, s)+1 2Q((d, s),(d, s)) + o(|d|2+|s)|2), where L(d, s)= @z @n⇥(0,T ),Q((d, s),(d, s)) = @w @n⇥(0,T ), Figure: Some experiences with AUGLAG and DIRECT A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 2-D wave equation VIII Case of a ball Figure: Desired and computed radius and centers of the ball A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 2-D wave equation I Case of an ellipse Optimization problem: case of an ellipse Given: eα=eα(x,t). Find x0,y0and θand a,bsuch that (x0,y0, θ, a,b)∈Xeand J(x0,y0, θ, a,b)≤J(x0 0,y0 0, θ0,a0,b0)∀(x0 0,y0 0, θ0,a0,b0)∈Xe,(2) the function J:Xe7→ Ris defined by J(x0,y0, θ, a,b) := 1 2ZZγ×(0,T) |α[x0,y0, θ, a,b]−eα|2ds dt with α[x0,y0, θ, a,b] = ∂u ∂non γ×(0,T) Xe:= {(x0,y0, θ, a,b)∈R5:E(x0,y0, θ, a,b)⊂Ω} A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 2-D wave equation II Case of an ellipse Resolution of an optimization problem. Now, J=J(x0,y0, θ, a,b) Test 3:T=5, u0=10x,u1=0, ϕ=10x x0des=-3,y0des=-3,sin(thetades)=0,ades=0.8,bdes=0.4 x0ini=-1,y0ini=-1,sin(thetaini)=0,aini=0.5,bini=0.5 NLopt (AUGLAG + DIRECTNoScal), NoIter = 2001, FreeFem++ x0cal = -2.963301665 y0cal = -3.035106437 sin(thetacal) = 0.112178021 acal = 0.8446502058 bcal = 0.4166666667 A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 2-D Lam´ e system III Case of an ellipse Figure: Computed solution at the final time A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 3-D wave equation I Case of a sphere Test 6:T=5, u0=10x,u1=0, ϕ=10x x0des = -2,y0des = -2,z0des = -2,rdes = 1 x0ini = 0,y0ini = 0,z0des = 0,rini = 0.6 NLopt (AUGLAG + DIRECTNoScal), NoIter = 438, FreeFem++ x0cal = -1.975308642 y0cal = -2.232383275 z0cal = -2.305542854 rcal = 1.05 A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 3-D wave equation II Case of a sphere Figure: Initial mesh. Points: 829, tetrahedra: 4023, faces: 8406, edges: 5210, boundary faces: 720, boundary edges: 1080 A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 3-D wave equation III Case of a sphere Figure: Desired configuration A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 3-D wave equation IV Case of a sphere Figure: Computed observation and mesh A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 3-D wave equation V Case of a sphere Figure: Solution of the wave equation corresponding to computed data A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 3-D Lam´ e system I Case of a sphere Test 7:T=5, u01 =10x,u02 =10y,u03 =10z u11 =0, u12 =0, u13 =0 ϕ1=10x,ϕ2=10y,ϕ3=10z x0des = -2,y0des = -2,z0des = -2,rdes = 1 x0ini = 0,y0ini = 0,z0des = 0,rini = 0.6 NLopt (AUGLAG + DIRECTNoScal), NoIter = 444 , FreeFem++: x0cal = -1.981405274 y0cal = -2.225232904 z0cal = -2.148084171 rcal = 0.9504115226 A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 3-D Lam´ e system II Case of a sphere Figure: Initial mesh. Points: 829, tetrahedra: 4023, faces: 8406, edges: 5210, boundary faces: 720, boundary edges: 1080 A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 3-D Lam´ e system III Case of a sphere Figure: Desired configuration A. Doubova Numerical solution of some geometric inverse problemsm Numerical results: 3-D Lam´ e system IV Case of a sphere A. Doubova Numerical solution of some geometric inverse problemsm