Increasing stability in the linearized inverse Schrödinger potential problem with power type nonlinearities

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY-NC-ND 3.0 https://creativecommons.org/licenses/by-nc-nd/3.0/ Increasing stability in the linearized inverse Schrödinger potential problem with power type nonlinearities © 2022 IOP Publishing Ltd. Accepted version (Final draft) Lu, Shuai; Salo, Mikko; Xu, Boxi Lu, S., Salo, M., & Xu, B. (2022). Increasing stability in the linearized inverse Schrödinger potential problem with power type nonlinearities. Inverse Problems, 38(6), Article 065009. https://doi.org/10.1088/1361-6420/ac637a 2022 Inverse Problems ACCEPTED MANUSCRIPT Increasing stability in the linearized inverse Schrödinger potential problem with power type nonlinearities To cite this article before publication: Shuai Lu et al 2022 Inverse Problems in press https://doi.org/10.1088/1361-6420/ac637a Manuscript version: Accepted Manuscript Accepted Manuscript is “the version of the article accepted for publication including all changes made as a result of the peer review process, and which may also include the addition to the article by IOP Publishing of a header, an article ID, a cover sheet and/or an ‘Accepted Manuscript’ watermark, but excluding any other editing, typesetting or other changes made by IOP Publishing and/or its licensors” This Accepted Manuscript is © 2022 IOP Publishing Ltd. During the embargo period (the 12 month period from the publication of the Version of Record of this article), the Accepted Manuscript is fully protected by copyright and cannot be reused or reposted elsewhere. As the Version of Record of this article is going to be / has been published on a subscription basis, this Accepted Manuscript is available for reuse under a CC BY-NC-ND 3.0 licence after the 12 month embargo period. After the embargo period, everyone is permitted to use copy and redistribute this article for non-commercial purposes only, provided that they adhere to all the terms of the licence https://creativecommons.org/licences/by-nc-nd/3.0 Although reasonable endeavours have been taken to obtain all necessary permissions from third parties to include their copyrighted content within this article, their full citation and copyright line may not be present in this Accepted Manuscript version. Before using any content from this article, please refer to the Version of Record on IOPscience once published for full citation and copyright details, as permissions will likely be required. All third party content is fully copyright protected, unless specifically stated otherwise in the figure caption in the Version of Record. View the article online for updates and enhancements. This content was downloaded from IP address 130.234.90.39 on 04/04/2022 at 07:23 Increasing stability in the linearized inverse Schrödinger potential problem with power type nonlinearities‡ Shuai Lu1, Mikko Salo2and Boxi Xu3,4 1Shanghai Key Laboratory for Contemporary Applied Mathematics, Key Laboratory of Mathematics for Nonlinear Sciences and School of Mathematical Sciences, Fudan University, Shanghai, China. 2Department of Mathematics and Statistics, University of Jyväskylä, Jyväskylä, Finland. 3School of Mathematics, Shanghai University of Finance and Economics, Shanghai, China. 4Author to whom any correspondence should be addressed. E-mail: [email protected], [email protected], [email protected] Abstract. We consider increasing stability in the inverse Schrödinger potential problem with power type nonlinearities at a large wavenumber. Two linearization approaches, with respect to small boundary data and small potential function, are proposed and their performance on the inverse Schrödinger potential problem is investigated. It can be observed that higher order linearization for small boundary data can provide an increasing stability for an arbitrary power type nonlinearity term if the wavenumber is chosen large. Meanwhile, linearization with respect to the potential function leads to increasing stability for a quadratic nonlinearity term, which highlights the advantage of nonlinearity in solving the inverse Schrödinger potential problem. Noticing that both linearization approaches can be numerically approximated, we provide several reconstruction algorithms for the quadratic and general power type nonlinearity terms, where one of these algorithms is designed based on boundary measurements of multiple wavenumbers. Several numerical examples shed light on the efficiency of our proposed algorithms. Keywords: increasing stability, inverse Schrödinger potential problem, power type nonlinearities, reconstruction algorithms. Submitted to: Inverse Problems ‡S. Lu is supported by NSFC (No.11925104), Science and Technology Commission of Shanghai Municipality (19XD1420500, 21JC1400500). M. Salo is supported by the Academy of Finland (Finnish Centre of Excellence in Inverse Modelling and Imaging, grant 284715) and by the European Research Council under Horizon 2020 (ERC CoG 770924). B. Xu is supported by NSFC (No.12171301 and No.11801351). Page 1 of 28 AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript 1. Introduction 1.1. Background The inverse Schrödinger potential problem arises from electrical impedance tomography (EIT) [9] and has attracted much attention both theoretically and computationally. In a general setting, we can formulate the following Schrödinger equation (∆u+k2u−c(x)u= 0 in Ω⊂Rn, u=g0on ∂Ω,(1.1) where, throughout the article, Ω⊂Rnis assumed to be a bounded open domain with smooth boundary ∂Ωand dimension n≥2. The inverse Schrödinger potential problem is to identify the unknown potential function c(x)from many boundary measurements or the Dirichlet-to- Neumann map defined below. A classical result in [1] shows that if the wavenumber k= 0 in (1.1) the stability of the inverse Schrödinger potential problem is logarithmic. When the wavenumber is sufficiently large, increasing stability with respect to the wavenumber khas been observed and well documented, starting with [15] and with many further results given in [19,17,18] for (1.1) or its linearized form. These results are often stated as stability estimates involving a Hölder term and a logarithmic term which goes to zero as the wavenumber goes to infinity. An alternative way to observe increasing stability is to note that one can determine the Fourier transform of the unknown coefficient in a stable way for a range of frequencies, and that this range increases with the wavenumber. We note that these increasing stability results have also been verified both theoretically and numerically in other inverse source, obstacle or medium problems where we refer to [5,7,27,2,3,4,11,31,16,20,6,8] and references therein. There have also been several recent works on inverse problems for nonlinear elliptic equations. In such problems, it has been observed that higher order linearizations of the nonlinear Dirichlet-to-Neumann map carry information about the unknown coefficients. This method allows one to exploit nonlinear effects in order to obtain better results than those that are currently known for corresponding linear equations. The higher order linearization method goes back to [22] in the hyperbolic case and to [13,24] in the elliptic case. The method has been further applied to more general equations and partial data problems. See [23,21,26,10,25] for a selection of recent results. This article studies possible improvements in stability properties of inverse problems for nonlinear Schrödinger type equations with a large wavenumber. More specifically, we study the inverse Schrödinger potential problem with an arbitrary power type nonlinearity term and discuss its unique determination, increasing stability and numerical reconstruction algorithms. In particular, we consider the problem of recovering the potential function c(x), defined in Ω⊂Rn, in the following nonlinear Schrödinger equation, with an integer m≥2 2 Page 2 of 28AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript denoting the nonlinearity index, (∆u+k2u−c(x)um= 0 in Ω, u=g0on ∂Ω,(1.2) from many boundary measurements. Here, we assume that the squared wavenumber k2is sufficiently large, 0is not a Dirichlet eigenvalue of ∆ + k2in Ωand the Dirichlet boundary data g0is sufficiently small. Meanwhile, by assuming that c(x)is compactly supported in Ω, the well-posedness of the forward problem (1.2) can be verified following the variational framework developed in [12, Theorem 1]. Thus, the boundary measurements can be given by the nonlinear Dirichlet-to-Neumann (DtN) map Λc:g07→ ∂νuon ∂Ω.(1.3) The precise definition of Λcand its two linearized forms Dm 0Λc,Λ′ cwill be specified later. When the wavenumber k= 0 in (1.2), unique identification of the potential function c(x) has been provided in [25] by measurement of the DtN map in (1.3) and its linearized form Dm 0Λc. In current work, we particularly focus on the increasing stability estimate for (1.2) under two different linearized forms of Λc. 1.2. Linearization approaches To solve the nonlinear inverse Schrödinger potential problem stably, we implement linearization approaches and discuss recovery of the potential function by the linearized DtN maps accordingly. In this subsection, we briefly overview two linearization approaches with respect to small boundary data and small potential function, which have been studied in linear and nonlinear elliptic inverse problems, for instance in [9,24], when k= 0 in (1.1) and (1.2). To treat elliptic equations with power type nonlinearities, a novel linearization approach with respect to small boundary data has recently been discussed in [13,24] and been extended to a fractional nonlinearity index in [26]. We briefly introduce its extension to the nonlinear Schrödinger potential problem (1.2) below. Assume that c∈Cα(Ω) for some αwith 0< α < 1, and 0is not a Dirichlet eigenvalue of ∆ + k2in Ω. By [26, Proposition 2.1], we can find a constant τ > 0such that for any Dirichlet boundary value fin Uτ:= {f∈C2,α(∂Ω) : kfkC2,α(∂Ω) ≤τ}, there is a unique small solution u∈C2,α(Ω) and u|∂Ω=f. The nonlinear DtN map in the Hölder spaces is defined by Λc:Uτ⊂C2,α(∂Ω) →C1,α(∂Ω), f 7→ ∂νu|∂Ω. Let ε= (ε1, . . . , εm)where each εj>0is small, and consider the solution uεcorresponding to the Dirichlet boundary value fε=ε1f1+. . . +εmfm. 3 Page 3 of 28 AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript By [26, Proposition 2.1] the solution uεdepends smoothly on the parameters εj. We may thus differentiate the equation ∆uε+k2uε−c(x)um ε= 0 in Ω, uε|∂Ω=fε(1.4) with respect to the parameters εj. Writing vj=∂εjuε|ε=0, we observe that vjis the unique solution of ∆vj+k2vj= 0 in Ω, vj|∂Ω=fj. Similarly, applying ∂ε1· · · ∂εmto the equation (1.4) and setting ε= 0, we can define w=∂ε1· · · ∂εmuε|ε=0 which solves the equation ∆w+k2w= (m!)c(x)v1· · · vmin Ω, w|∂Ω= 0.(1.5) Moreover, the Neumann boundary data can be obtained in form of ∂νw|∂Ω=∂ε1· · · ∂εm(∂νuε|∂Ω)|ε=0 =∂ε1· · · ∂εmΛc(fε)|ε=0 =Dm 0Λc(f1, . . . , fm)(1.6) where Dm 0denotes the mth Fréchet derivative at 0considered as an m-linear form. If we integrate the equation (1.5) against another function vm+1 solving ∆vm+1 +k2vm+1 = 0 in Ω, vm+1|∂Ω=fm+1, we obtain a Calderón or Alessandrini type identity (m!) ZΩ c(x)v1· · · vmvm+1 dx=Z∂Ω Dm 0Λc(f1, . . . , fm)fm+1 dS(1.7) which will be revisited later. Noticing that the mth Fréchet derivative Dm 0Λc(f1, . . . , fm)is numerically hard to obtain, we further consider the case when c(x)is small compared to the wavenumber and study the linearization approach with respect to the potential function as investigated in the linear Schrödinger potential problem in [18]. Taking the asymptotic expansion with respect to the potential function c(x), we have u=u0+u1+u2+. . . (1.8) where the remaining “. . .” denotes the “higher” order term and following subproblems are satisfied such that ∆u0+k2u0= 0, ∆u1+k2u1=c(x)um 0, ∆u2+k2u2=mc(x)um−1 0u1. This shows that u0satisfies the Helmholtz equation ∆u0+k2u0= 0 in Ωand the first-order expansion term u1satisfies ∆u1+k2u1=c(x)um 0in Ω.(1.9) 4 Page 4 of 28AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript When u0|∂Ω=g0and u1|∂Ω=g1≡0, the linearized DtN map Λ′ cis formally defined by Λ′ c:g07→ ∂νu1on ∂Ω.(1.10) Note that Λ′ cis actually a nonlinear map, since it corresponds to linearization with respect to the potential. Multiplying the above equation (1.9) from both sides with another φsolving ∆φ+k2φ= 0 in Ω, we thus obtain another Calderón or Alessandrini type identity ZΩ c(x)um 0φdx=Z∂Ω ∂νu1φdS(1.11) which will also be revisited later. In current article, we consider the following problem: Recover the potential function c(x)from the linearized DtN map Dm 0Λc or Λ′ c. The first main result shows that from the knowledge of the mth Fréchet derivative Dm 0Λc, one can determine the Fourier transform F[c](ξ)of cin a stable way for frequencies |ξ| ≤ (m+ 1)k. Thus the range of frequencies that can be determined stably increases both with respect to the wavenumber kand the nonlinearity index m. However, determining Dm 0Λcfrom Λcbecomes numerically very difficult when mincreases. The second main result considers the case where the potential function is small compared to the wavenumber. In this case we consider the linearization Λ′ c. We show that in the quadratic case where m= 2, from the knowledge of Λ′ cone can stably determine F[c](ξ)for frequencies |ξ| ≤ 3k. This is in contrast with the linear case where one can only determine frequencies |ξ| ≤ 2kstably [18]. Thus in both main results above, the nonlinearity leads to improved stability properties in a certain sense. The theoretical stability results are confirmed by numerical results given in the end of the article. The article is organized as follows. In Section 2we show that the linearized DtN map Dm 0Λcprovides a uniform increasing stability where the range of frequencies that can be determined stably increases with respect to kand m. On the other hand, Λ′ conly yields the uniqueness of the potential function c(x)in the general setting m≥2. In Section 3 we further explore the linearized DtN map Λ′ cfor the inverse Schrödinger potential problem with a quadratic nonlinearity term. By calibrating the identity (1.11) carefully, we verify an improved increasing stability for the specific inverse Schrödinger potential problem with a quadratic nonlinearity term. Noticing that both linearized DtN maps Dm 0Λcand Λ′ ccan be numerically approximated, we extend the reconstruction algorithm in [18] to the inverse Schrödinger potential problem with quadratic and general nonlinearity terms in Section 4, respectively. We note that one of these reconstruction algorithms is realized by the linearized DtN map Λ′ cwith multiple wavenumbers. In the same Section 4we provide some numerical examples and extended discussion verifying the efficiency of our proposed algorithms. 5 Page 5 of 28 AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript 2. Linearized inverse Schrödinger potential problem with an arbitrary power type nonlinearity term In this section, we investigate the linearized inverse Schrödinger potential problem with an arbitrary power type nonlinearity term provided with the linearized DtN map Dm 0Λcor Λ′ c. The analysis is based on the Calderón type identities (1.7) and (1.11). For the map Dm 0Λcwith k= 0, [24] has verified that by linearizing the small boundary data, the stability estimate for the inverse potential problem is logarithmic, which is consistent with the classical result in EIT [1]. In the current section, we verify that by constructing an appropriate set of complex exponential solutions, there will be improved stability when the wavenumber is large. Recall the identity (1.7), (m!) ZΩ c(x)v1· · · vmvm+1 dx=Z∂Ω Dm 0Λc(f1, . . . , fm)fm+1 dS. Here vjsolve ∆vj+k2vj= 0 in Ωwith vj|∂Ω=fj. To derive the stability estimate, we rely on the above identity and observe that ZΩ c(x)v1· · · vmvm+1 dx ≤ϵ m! m Y j=1 kfjkC2,α(∂Ω)!kfm+1kL2(∂Ω) (2.1) where we define ϵ:= sup∥fj∥C2,α(∂Ω)≤1kDm 0Λc(f1, . . . , fm)kL2(∂Ω). Thus (2.1) further yields the inequality ZΩ c(x)v1· · · vm+1 dx ≤ϵ m! m+1 Y j=1 kvjkC2,α(Ω).(2.2) Let F[c](ξ)denote the Fourier transform of c(extended by zero outside Ω) at a frequency ξ∈Rn. The following result shows that frequencies |ξ| ≤ (m+ 1)kcan be recovered in a Lipschitz stable way from the knowledge of the linearized map Dm 0Λc. Theorem 2.1. Let m≥2be an integer, k≥1, and assume that |ξ| ≤ (m+ 1)k. Then |F[c](ξ)| ≤ ϵ m!3(1 + k6)m+1 2. Proof. We first claim that if ℓ≥2is an integer, then for any η∈Rnwith |η| ≤ ℓthere are unit vectors ω1, . . . , ωℓ∈Rnsuch that ℓ X j=1 ωj=η. This can be proved by induction. When ℓ= 2 and |η| ≤ 2, we may choose ω1=η 2+r1−|η|2 4ω, ω2=η 2−r1−|η|2 4ω, 6 Page 6 of 28AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript where ωis any unit vector orthogonal to η. We make the induction hypothesis that the claim holds for some ℓ≥2. Let ηbe a vector with |η| ≤ ℓ+ 1. We can write η=η0+ ˜ω where η0and ˜ωare parallel to ηand |η0| ≤ ℓ,|˜ω|= 1. Applying the induction hypothesis to η0gives unit vectors ω1, . . . , ωℓthat add up to η0. The induction step is completed by setting ωℓ+1 = ˜ω. To prove the theorem we choose special solutions of ∆vj+k2vj= 0 in Ω(j= 1, . . . , m+1) having the form vj= eiζj·x where ζj∈Cnsatisfy ζj·ζj=k2. Since |ξ k| ≤ m+ 1, the claim above shows that we can find unit vectors ω1, . . . , ωm+1 such that m+1 X j=1 ωj=ξ k. Thus, choosing ζj=kωj, we have m+1 X j=1 ζj=ξ. It follows that v1· · · vm+1 = eiξ·x. Now (2.2) shows that |F[c](ξ)| ≤ ϵ m! m+1 Y j=1 kvjkC3(Ω). The proof is completed upon observing that kvjk2 C3(Ω) ≤3(1 + k6)when k≥1. The assumption |ξ| ≤ (m+ 1)kensured that we could choose solutions vj= eiζj·xwith ζjpurely real in the proof. When |ξ|>(m+ 1)kthis will no longer be possible, and there will be a logarithmic component in the increasing stability estimate. We will next prove such an estimate for the linearized DtN map Dm 0Λcby making a more careful choice of the vectors ζj. Without loss of generality we assume that 0∈Ωand denote D:= 2 supx∈Ω|x|. Theorem 2.2. Let D≤1,kckC1(Ω) ≤M1, and k > 1,ϵ < 1, then the following estimate holds true kck2 L2(Ω) ≤Ckn+6(m+1)ϵ2+CEn+6(m+1)ϵ+M2 1 1 + m2k2+E2 for the linearized system (1.5) with E=−ln ϵand the constant Cdepending on the domain Ω, the nonlinearity index mand the dimensionality n. Proof. To prove the stability estimate, we shall choose the complex exponential solutions vjin (2.2) carefully. Let ξ∈Rnwith ξ6= 0 and choose an orthonormal base ne1:= ξ |ξ|, e2, . . . , eno 7 Page 7 of 28 AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript and the identity (3.3) yields 2F[c](ξ) = 2 ZΩ c(x)eiξ·xdx=Z∂Ω ∂νw1φdS−Z∂Ω ∂νu1φdS+Z∂Ω ∂νv1φdS. Hence we obtain, since w0=u0+v0in Ω, |F[c](ξ)|2≤1 4k∂νw1k2 L3 2(∂Ω) +k∂νu1k2 L3 2(∂Ω) +k∂νv1k2 L3 2(∂Ω)kφk2 L3(∂Ω) ≤1 4ϵ2kw0|∂Ωk4 C2(∂Ω) +ku0|∂Ωk4 C2(∂Ω) +kv0|∂Ωk4 C2(∂Ω)kφk2 L3(∂Ω) ≤Cϵ2kw0k4 C2(Ω) +ku0k4 C2(Ω) +kv0k4 C2(Ω)kφk2 L∞(Ω) ≤Cϵ2ku0k4 C2(Ω) +kv0k4 C2(Ω) with a generic constant Cdepending on the domain Ω, and kφk2 L∞(Ω) =keike1·xk2 L∞(Ω) ≤1. Noticing the fact that |ζℓ|2=k2,ℓ= 1,2,3, we thus obtain, if k≥|ξ| 3, ku0k4 C2(Ω) =kv0k4 C2(Ω) ≤C1 + k8. Then, there holds |F[c](ξ)|2≤Cϵ21 + k8,for k≥|ξ| 3. If k < |ξ| 3, by denoting Ξ := p|ξ|2−2k|ξ| − 3k2we then derive the following bounds ku0k4 C2(Ω) =kv0k4 C2(Ω) ≤C1 + k8keiζ1·xk4 L∞(Ω) ≤C1 + k8eDΞ. Consequently, we derive |F[c](ξ)|2≤Cϵ21 + k8eDΞ,for k < |ξ| 3. Let E:= −ln ϵ > 0and k > 1,ϵ < 1, we again consider two cases a) k > E (i.e. ϵ= e−E>e−k), and b) k≤E(i.e. ϵ= e−E≤e−k). In the case a), we have kck2 L2(Ω) =Z|F[c](ξ)|2dξ=Zk≥|ξ| 3 |F[c](ξ)|2dξ+Zk< |ξ| 3 |F[c](ξ)|2dξ ≤Cϵ21 + k8σn(3k)n+M2 1 1 + (3k)2 ≤Ckn+8ϵ2+M2 1 1 + 8k2+E2 where σnis the volume of an unit ball in Rn, and the constant Cdepends on the domain Ω and the dimensionality n. 14 Page 14 of 28AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript In the case b), we let ρ:= k+q4k2+E D2such that pρ2−2kρ −3k2=E Dand split kck2 L2(Ω) =Zk≥|ξ| 3 |F[c](ξ)|2dξ+Zk< |ξ| 3<ρ 3 |F[c](ξ)|2dξ +Zρ≤|ξ| |F[c](ξ)|2dξ. (3.5) The first term in the right-hand side of (3.5) can be bounded by Zk≥|ξ| 3 |F[c](ξ)|2dξ≤Ckn+8ϵ2≤CEn+8ϵ2, noticing k≤E. We focus on the second term in the right-hand side of (3.5) and estimate Zk< |ξ| 3<ρ 3 |F[c](ξ)|2dξ≤Cϵ2k8 Zk< |ξ| 3<ρ 3 eDΞdξ!≤Cϵk8 Zk< |ξ| 3<ρ 3 dξ! since eDΞ≤eE=ϵ−1when k < |ξ| 3<ρ 3. Meanwhile, noticing ρ≤3k+E Dand k≤E, we bound Zk< |ξ| 3<ρ 3 dξ=σn(ρn−(3k)n) ≤σn En Dn1 + 3kD En −3kD En ≤σn En Dn[(1 + 3D)n−(3D)n] where σnis the volume of an unit ball in Rn. The above inequalities yields Zk< |ξ| 3<ρ 3 |F[c](ξ)|2dξ≤Cϵk8σn En Dn[(1 + 3D)n−(3D)n]≤CEn+8ϵ. Furthermore, since ρ > q4k2+E D2, we bound the third term in the right-hand side of (3.5) by Zρ≤|ξ| |F[c](ξ)|2dξ≤M2 1 1 + ρ2≤M2 1 1 + 4k2+E2 D2 . We thus prove for both cases the proposed bound. Remark 3.2. The stability estimate in above Theorem 3.1, if kis sufficiently large, is similar as in [18, Theorem 2.1] where a linear elliptic equation is investigated ibid, i.e. (∆u+k2u−c(x)u= 0 in Ω, u=g0on ∂Ω. 15 Page 15 of 28 AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript A clear numerical evidence will be provided in Section 4and one can stably recover the Fourier coefficients with |ξ| ≤ 3kwhereas in [18] one can only recover those with |ξ| ≤ 2k. Such gain highly depends on the sophisticatedly selected complex exponential functions and the modified identity (3.3) considered above. It can be viewed as the advantage of the quadratic nonlinearity term when we solve the linearized inverse problems (1.2) with m= 2. 4. Reconstruction algorithm and numerical examples In this section, we provide two reconstruction algorithms stably recovering the unknown potential function by the linearized DtN map Λ′ cand a vanilla reconstruction algorithm by the linearized DtN map Dm 0Λc. In view of the quadratic nonlinearity term, we rely on the theoretical discussion in Section 3and deliver the first algorithm where boundary measurements of a single (large) wavenumber could offer a high resolution. Meanwhile, the second algorithm focuses on the high-order nonlinearity term discussed in Section 2and the linearized DtN map Λ′ cof multiple wavenumbers is included to recover sufficiently many Fourier coefficients of the unknown potential function. Finally a vanilla reconstruction algorithm by the (approximated) linearized DtN map Dm 0Λcis presented to verify the feasibility of the proposed linearization which, to the best of our knowledge, is the first attempt to realize the linearized DtN map Dm 0Λcnumerically. We shall emphasize that by implementing the linearized DtN maps Λ′ cand Dm 0Λcthe range of stably reconstructed Fourier mode F[c](ξ)is expanded to |ξ| ≤ (m+1)kas shown in Theorem 2.2 with an arbitrary finite integer mand Theorem 3.1 with m= 2. This truncated value (m+1)kcan be viewed as a regularization for stably recovering the unknown potential function c(x). 4.1. Reconstruction algorithm by Λ′ cfor a quadratic nonlinearity term Noticing that the linearized DtN map Λ′ ccan be numerically approximated, see [18, Eq.(4.3)], we present the first reconstruction algorithm based on the identity (3.3). As an illustration, we focus on the two-dimensional space n= 2. By selecting the complex exponential solutions (3.4) in the proof of Theorem 3.1, we know that the left-hand side of (3.3) reflects a Fourier coefficient of the potential function c(x). Then by choosing ξ∈Rnand recalling Remark 3.2, we aim to recovering all the Fourier coefficients F[c](ξ)of the potential function c(x)satisfying |ξ| ≤ 3k. The larger wavenumber k, the more Fourier coefficients can be recovered. To further address the reconstruction algorithm, we need the following discrete sets of lengths and angles of the vectors in the phase space. The discrete and finite length set is defined by {κi}I i=1 ⊂(0, Lk ]for any fixed k. 16 Page 16 of 28AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript Here we choose L≥3and Lk is the maximum length of the vector ξ. Two angle sets are defined by {ˆys}S s=1 ⊂Sn−1and {ˆzs}S s=1 ⊂Sn−1, which satisfy ˆys·ˆzs= 0. The vector ξ⟨i;s⟩:= κiˆysand following vectors ζ⟨i;s⟩ ℓ∈Cn,ℓ= 1,2,3are chosen            ζ⟨i;s⟩ 1:= 1 2(−k+κi)ˆys−1 2q3k2+ 2kκi−κ2 iˆzs, ζ⟨i;s⟩ 2:= 1 2(−k+κi)ˆys+1 2q3k2+ 2kκi−κ2 iˆzs, ζ⟨i;s⟩ 3:= kˆys, which further assign to the complex exponential solution as in (3.4) below u0(x) = eiζ⟨i;s⟩ 1·x, v0(x) = eiζ⟨i;s⟩ 2·x, φ(x) = eiζ⟨i;s⟩ 3·x for i= 1,2,· · · , I and s= 1,2,· · · , S. Here, the superscript notation ·⟨i;s⟩will be referred to a vector ξ⟨i;s⟩with the ith length κiand the sth angle ˆys. Finally, for the inverse Fourier transform, a numerical quadrature rule can be constructed by a suitable choice of the weights σ⟨i;s⟩according to these points ξ⟨i;s⟩. We summarize our reconstruction algorithm below, which is similar to that in [18] but one has to solve the nonlinear Schrödinger potential problem three times at each iteration because of the quadratic nonlinearity term. Algorithm 1: Reconstruction Algorithm for the Linearized Schrödinger Potential Problem, the quadratic nonlinearity term Input: k,{κi}I i=1,{ˆys}S s=1,{ˆzs}S s=1 and {σ⟨i;s⟩}; Output: Approximated Potential c⟨I+1;1⟩. 1: Set c⟨1;1⟩:= 0; 2: For i= 1,2, . . . , I (length updating) 3: For s= 1,2, . . . , S (angle updating) 4: Choose u0:= exp{iζ⟨i;s⟩ 1·x},v0:= exp{iζ⟨i;s⟩ 2·x}and w0:= u0+v0; 5: Measure the Neumann boundary data ∂νu,∂νv,∂νwof the forward problem (1.2) while the Dirichlet boundary data u0|∂Ω,v0|∂Ω,w0|∂Ωare given; 6: Calculate the approximated linearized Neumann boundary data g′ u:= (∂νu−∂νu0)|∂Ω,g′ v:= (∂νv−∂νv0)|∂Ω,g′ w:= (∂νw−∂νw0)|∂Ω; 7: Choose φ:= exp{iζ⟨i;s⟩ 3·x}and γ:= [u0v0φ]−1= exp{−iξ⟨i;s⟩·x}; 8: Compute F[c](ξ⟨i;s⟩)≈1 2R∂Ω(g′ w−g′ u−g′ v)φdS; 17 Page 17 of 28 AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript 9: Update c⟨i;s+1⟩:= c⟨i;s⟩+F[c](ξ⟨i;s⟩)γσ⟨i;s⟩, if κi≤3k; 10: End; 11: Set c⟨i+1;1⟩:= c⟨i;S+1⟩; 12: End. In fact, the linearized Neumann boundary data ∂νw1depends on the unknown potential function c(x)referring to (3.2). As mentioned in [18, Eq.(4.3)], we utilize g′ w:= (∂νw−∂νw0)|∂Ω(4.1) to approximate the non-measurable data ∂νw1|∂Ω. The similar approximation g′ uand g′ vare employed for the linearized Neumann data ∂νu1|∂Ωand ∂νv1|∂Ω, respectively. As one can observe, the computational cost of Algorithm 1 is quite high because of the nonlinearity term in the forward problem, see e.g. [14,28,29,30]. In particular in Steps 5-6 of Algorithm 1, we must solve the nonlinear elliptic equation (1.2) three times in order to derive their Neumann traces which are necessary to compute the Fourier coefficient in Step 8. Figure 1. The sampling points ξ= (ξ1, ξ2)in frequency domain. To numerically test Algorithm 1, we consider the domain Ω = B0.5(0) in a square [0.5,0.5]2. To avoid the inverse crime, we use a fine grids (200 ×200 equal-distance points) for the forward problem and a coarse grid (90 ×90 equal-distance points) for the inversion. The sampling points ξ= (ξ1, ξ2)in frequency domain are shown in Figure 1, marked by blue “∗” near which all the Fourier coefficients will be recovered. In Figure 2, the horizontal axis shows the length |ξ|of all ξ, and the vertical axis shows the absolute value |F[c](ξ)|of Fourier coefficients near the sampling points. By comparing the exact (top) and recovered (bottom) Fourier coefficients in each sub-figure of Figure 2: (a) k= 5 and (b) k= 10, we 18 Page 18 of 28AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript conclude that, while kis larger, the more Fourier modes can be recovered stably, i.e. F[c](ξ) with |ξ| ≤ 3k. Quadratic case: (a) k= 5 (b) k= 10 Figure 2. The exact (Top) and recovered (Bottom) Fourier coefficients F[c](ξ)in each sub-figure: (a) k= 5 and (b) k= 10. Here, the horizontal axis shows the length |ξ|of ξ; the vertical axis shows the absolute value |F[c](ξ)|of Fourier coefficients. Then, by using all the recovered Fourier coefficients F[c](ξ)with |ξ| ≤ 3k, we implement the inverse Fourier transform in Step 9 to reconstruct the potential function c(x). In Figure 3, we present the exact and reconstructed potential functions c(x)with different wavenumbers: (a) k= 5 and (b) k= 10, respectively. These results numerically verify the increasing stability in Theorem 3.1 while kbecomes large. As one can observe, the point-wise absolute errors between the exact (left) and recovered (middle) potential functions are shown in Figure 3, and Algorithm 1 reduces the maximum absolute error from 0.5to 0.08 when k increases from 5to 10. 19 Page 19 of 28 AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript Quadratic case: (a) k= 5 (b) k= 10 Figure 3. The exact (Left) and recovered (Middle) potential functions c(x)together with the point-wise absolute error (Right) when (a) k= 5 and (b) k= 10. Here, we use the Fourier coefficients F[c](ξ)in the range |ξ| ≤ 3k. 20 Page 20 of 28AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript 4.2. Reconstruction algorithm by Λ′ cfor high-order nonlinearity terms with multiple wavenumbers In this subsection, we show that the uniqueness result Theorem 2.4 in Section 2indeed could provide a stable reconstruction algorithm for the nonlinear inverse Schrödinger potential problem whose nonlinearity index is an arbitrary finite integer m≥2, if the linearized DtN map Λ′ cof multiple wavenumbers is provided. As highlighted in Remark 2.5, the complex exponential solutions constructed in the proof of Theorem 2.4 has a stable interval [(m−1)k, (m+1)k]for any fixed k. Suppose that the same discrete (phase space) length and angle sets of the vectors in Section 4.1 could be used, i.e. {κi}I i=1,{ˆys}S s=1,{ˆzs}S s=1 and the following vectors µ⟨i;s⟩ ℓ∈Cn,ℓ= 1,2are chosen          µ⟨i;s⟩ 1:= +(m2−1)k2+κ2 i 2mκi ˆys−p−(m2−1)2k4+ 2(m2+ 1)k2κ2 i−κ4 i 2mκi ˆzs, µ⟨i;s⟩ 2:= −(m2−1)k2+κ2 i 2κi ˆ ys+p−(m2−1)2k4+ 2(m2+ 1)k2κ2 i−κ4 i 2κi ˆ zs, similar to Algorithm 1, we summarize a plain reconstruction algorithm for high-order nonlinearity terms with a fixed wavenumber k, according to the identity (1.11). Algorithm 2: Reconstruction Algorithm for the Linearized Schrödinger Potential Problem, the high-order nonlinearity term Input: m,k,{κi}I i=1,{ˆys}S s=1,{ˆzs}S s=1 and {σ⟨i;s⟩}; Output: Approximated Potential c⟨I+1;1⟩. 1: Set c⟨1;1⟩:= 0; 2: For i= 1,2, . . . , I (length updating) 3: For s= 1,2, . . . , S (angle updating) 4: Choose u0:= exp{iµ⟨i;s⟩ 1·x}; 5: Measure the Neumann boundary data ∂νuof the forward problem (1.2) while the Dirichlet boundary data u0|∂Ωare given; 6: Calculate the approximated linearized Neumann boundary data g′ u:= (∂νu−∂νu0)|∂Ω; 7: Choose φ:= exp{iµ⟨i;s⟩ 2·x}and γ:= [um 0φ]−1= exp{−iξ⟨i;s⟩·x}; 8: Compute F[c](ξ⟨i;s⟩)≈R∂Ωg′ uφdS; 9: Update c⟨i;s+1⟩:= c⟨i;s⟩+F[c](ξ⟨i;s⟩)γσ⟨i;s⟩, if κi∈[(m−1)k, (m+ 1)k]; 10: End; 11: Set c⟨i+1;1⟩:= c⟨i;S+1⟩; 21 Page 21 of 28 AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript 12: End. Furthermore, if we could measure the boundary data by appropriate multiple wavenumbers, we could reconstruct sufficiently many Fourier coefficients of the unknown potential function. By choosing k1small and a threshold value Kas the maximum wavenumber, we choose a discrete set of multiple wavenumbers, namely {kj}J j=1 ⊂(0, K ],(4.2) which satisfies kj+1 =m+1 m−1kj. Below we present an updated reconstruction algorithm of Algorithm 2 for the linearized Schrödinger potential problem with a high-order nonlinearity term, i.e. the nonlinearity index m≥2, if the linearized DtN map Λ′ cof multiple wavenumbers can be obtained. Algorithm 2*: Reconstruction Algorithm for the Linearized Schrödinger Potential Problem with a high-order nonlinearity term (Multiple wavenumbers) Input: m,{kj}J j=1,{κi}I i=1,{ˆys}S s=1,{ˆzs}S s=1 and {σ⟨i;s⟩}; Output: Approximated Potential cinv := J P j=1 c⟨I+1;1⟩ j. 1: For j= 1,2, . . . , J (wavenumber updating) 2: Compute the approximated potential c⟨I+1;1⟩ jby using Algorithm 2 and a fixed kj; 3: End. As an illustration, we consider the linearized Schrödinger potential problem with a cubic nonlinear term (m= 3). The wavenumber set in (4.2) is set with k1= 1.25 and K= 10 where we recover the Fourier coefficients F[c](ξ)with 4wavenumbers k∈ {1.25,2.5,5,10}. In Figure 4, the red region indicates the Fourier coefficients within [2k1,4k1) = [2.5,5), the green region indicates the Fourier modes within [2k2,4k2) = [5,10), the blue region indicates the Fourier modes within [2k3,4k3) = [10,20), and the cyan region indicates the Fourier coefficients within [2k4,4k4) = [20,40). By using Fourier coefficients F[c](ξ)within |ξ| ∈ SJ j=1 [(m−1)kj,(m+ 1)kj) = [(m−1)k1,(m+ 1)kJ), we implement the inverse Fourier transform to reconstruct the potential function c(x). In Figure 5, we present the exact (left) and reconstructed (right) potential functions c(x)with 4wavenumbers k∈ {1.25,2.5,5,10}. It can be seen that, by including the boundary measurements of four wavenumbers, we have obtained a good approximation of the unknown potential function in (1.2) with a cubic nonlinear term m= 3. 22 Page 22 of 28AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript Cubic case: multiple wavenumbers k∈ {1.25,2.5,5,10} Figure 4. (Cubic case, m= 3) The exact (Top) and recovered (Bottom) Fourier coefficients F[c](ξ)with multiple wavenumbers k∈ {1.25,2.5,5,10}. Here, the horizontal axis shows the length |ξ|of ξ; the vertical axis shows the absolute value |F[c](ξ)|of Fourier coefficients. Cubic case: multiple wavenumbers k∈ {1.25,2.5,5,10} Figure 5. (Cubic case, m= 3) The exact (Left) and recovered (Middle) potential functions c(x)together with the point-wise absolute error (Right) when multiple wavenumbers k∈ {1.25,2.5,5,10}are considered. 23 Page 23 of 28 AUTHOR SUBMITTED MANUSCRIPT - IP-103351.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript