Linearized inverse Schr\"odinger potential problem at a large wavenumber

TL;DR

This paper studies the linearized inverse Schrödinger potential problem at high wavenumbers, demonstrating improved stability and proposing an efficient Fourier mode reconstruction algorithm verified by numerical examples.

Contribution

It introduces a stability analysis at large wavenumbers and develops a new reconstruction algorithm for the potential function.

Findings

Hölder type stability at large wavenumbers

Exponential dependence of stability on attenuation constant

Successful numerical reconstruction of Fourier modes

Abstract

We investigate recovery of the (Schr\"odinger) potential function from many boundary measurements at a large wavenumber. By considering such a linearized form, we obtain a H\"older type stability which is a big improvement over a logarithmic stability in low wavenumbers. Furthermore we extend the discussion to the linearized inverse Schr\"odinger potential problem with attenuation, where an exponential dependence of the attenuation constant is traced in the stability estimate. Based on the linearized problem, a reconstruction algorithm is proposed aiming at the recovery of the Fourier modes of the potential function. By choosing the large wavenumber appropriately, we verify the efficiency of the proposed algorithm by several numerical examples.

Equations

$$\{\begin{array}{rlrl}-\Delta u-(k^2-c)u& =0& & \text{in }\Omega \subset {\mathbb{R}}^n,\\ u& =g_0& & \text{on }\partial \Omega ,\end{array}$$

$$\Lambda :g_0\mapsto g_1:={\partial }_{\nu }u\text{on }\partial \Omega$$

$$(I)\{\begin{array}{rlrl}-\Delta u-(k^2-c)u& =0& & \text{in }\Omega \subset {\mathbb{R}}^n,\\ u& =g_0& & \text{on }\partial \Omega ,\end{array}$$

$$\Lambda :g_0\mapsto g_1:={\partial }_{\nu }u\text{on }\partial \Omega .$$

$$(I_0)$$

$$(I_1)$$

$$u=u_0+u_1+\cdots$$

$${\Lambda }^{\prime }:g_0\mapsto {\partial }_{\nu }{u}_1\text{on }\partial \Omega .$$

$$-\Delta v-k^{2}v=0\text{in }\Omega ,$$

$${\int }_{\Omega }c{u}_{0}v={\int }_{\partial \Omega }\left({\partial }_{\nu }{u}_1\right)v.$$

$${\Vert c\Vert }_{L^2(\Omega )}^2⩽C(\Omega )\left(k^{n+4}+E^{n+4}\right){ϵ}^2+C(\Omega )E^{n+2}\left(ϵ+{ϵ}^3\right)+\frac{M_1^2}{1+E^2+3k^2}$$

$$u_0(x)=e^{\mathbf{i}\zeta \cdot x},v(x)=e^{\mathbf{i}{\zeta }^{\ast }\cdot x}.$$

$$\zeta :=\frac{\left|\xi \right|}{2}{e}_1+\sqrt{k^2-\frac{{\left|\xi \right|}^2}{4}}{e}_2,{\zeta }^{\ast }:=\frac{\left|\xi \right|}{2}{e}_1-\sqrt{k^2-\frac{{\left|\xi \right|}^2}{4}}{e}_2$$

$$\mathcal{F}[c](\xi ):={\int }_{\Omega }c(x)e^{\mathbf{i}\xi \cdot x}\mathrm{d}x={\int }_{\partial \Omega }\left({\partial }_{\nu }{u}_1\right)v.$$

$${\Vert u_0\Vert }_{H^1(\Omega )}^2={\Vert v\Vert }_{H^1(\Omega )}^2={\int }_{\Omega }\left(1+{\left|\zeta \right|}^2\right)\mathrm{d}x⩽\left(1+k^2\right){\mathrm{Vol}}_n\Omega ,$$

$${\Vert u_0\Vert }_{H^1(\Omega )}^2={\Vert v\Vert }_{H^1(\Omega )}^2$$

$$⩽\left(1+k^2\right){\mathrm{Vol}}_{n-1}\Omega \frac{e^{D\Xi }-e^{-D\Xi }}{2\Xi },$$

$${\Vert u_0\Vert }_{H^{\frac{1}{2}}(\partial \Omega )}⩽C(\Omega ){\Vert u_0\Vert }_{H^1(\Omega )},{\Vert v\Vert }_{H^{\frac{1}{2}}(\partial \Omega )}⩽C(\Omega ){\Vert v\Vert }_{H^1(\Omega )}.$$

$${\left|\mathcal{F}[c](\xi )\right|}^2$$

$$⩽{ϵ}^2{\Vert u_0\Vert }_{H^{\frac{1}{2}}(\partial \Omega )}^2{\Vert v\Vert }_{H^{\frac{1}{2}}(\partial \Omega )}^2$$

$$⩽{ϵ}^2{C}^4(\Omega ){\Vert u_0\Vert }_{H^1(\Omega )}^2{\Vert v\Vert }_{H^1(\Omega )}^2.$$

$${\left|\mathcal{F}[c](\xi )\right|}^2⩽C^4(\Omega ){\left(1+k^2\right)}^2{\left({\mathrm{Vol}}_n\Omega \right)}^2{ϵ}^2,\text{if }k⩾\frac{\left|\xi \right|}{2},$$

$${\left|\mathcal{F}[c](\xi )\right|}^2⩽C^4(\Omega ){\left(1+k^2\right)}^2{\left({\mathrm{Vol}}_{n-1}\Omega \right)}^2\frac{{\left(e^{D\Xi }-e^{-D\Xi }\right)}^2}{4{\Xi }^2}{ϵ}^2,\text{if }k<\frac{\left|\xi \right|}{2}.$$

$${\Vert c\Vert }_{L^2(\Omega )}^2$$

$$⩽C^4(\Omega ){\left(1+k^2\right)}^2{\left({\mathrm{Vol}}_n\Omega \right)}^2{\sigma }_n(2k{)}^n{ϵ}^2+\frac{M_1^2}{1+(2k{)}^2}$$

$$⩽C_1(\Omega )k^{n+4}{ϵ}^2+\frac{M_1^2}{1+4k^2}⩽C_1(\Omega )k^{n+4}{ϵ}^2+\frac{M_1^2}{1+E^2+3k^2}$$

$${\Vert c\Vert }_{L^2(\Omega )}^2$$

$$⩽C^4(\Omega ){\left(1+k^2\right)}^2{\left({\mathrm{Vol}}_n\Omega \right)}^2{\sigma }_n(2k{)}^n{ϵ}^2$$

$$+C^4(\Omega ){\left(1+k^2\right)}^2{\left({\mathrm{Vol}}_{n-1}\Omega \right)}^2\left({\int }_{k<\frac{|\xi |}{2}<\frac{\rho }{2}}\mathrm{d}\xi \right)\frac{D^2{\left(e^{\frac{E}{2}}-e^{-\frac{E}{2}}\right)}^2}{E^2}{ϵ}^2$$

$$+\frac{M_1^2}{1+{\rho }^2}.$$

$$\frac{e^y-e^{-y}}{y}=2\left(1+\frac{y^2}{3!}+\cdots +\frac{y^{2n}}{(2n+1)!}+\cdots \right)$$

$$\frac{e^{D\Xi }-e^{-D\Xi }}{2\Xi }⩽\frac{D\left(e^{\frac{E}{2}}-e^{-\frac{E}{2}}\right)}{E},\mathrm{since}k<\frac{\left|\xi \right|}{2}<\frac{\rho }{2}=\sqrt{\frac{E^2}{4D^2}+k^2}.$$

$${\int }_{k<\frac{|\xi |}{2}<\frac{\rho }{2}}\mathrm{d}\xi$$

$$={\sigma }_n\frac{E^n}{D^n}\left[{\left(1+{\left(2k\frac{D}{E}\right)}^2\right)}^{\frac{n}{2}}-{\left(2k\frac{D}{E}\right)}^n\right]$$

$$⩽{\sigma }_n\frac{E^n}{D^n}\left[{\left(1+{\left(2D\right)}^2\right)}^{\frac{n}{2}}-{\left(2D\right)}^n\right].$$

$${\Vert c\Vert }_{L^2(\Omega )}^2$$

Full text

{AMS}

35R30, 65N21

1 Introduction

We consider a linearized problem of recovering the potential function $c:=c(x)$ in the Schrödinger equation

[TABLE]

from many boundary measurements at a large wavenumber $k$ in the dimension $n\geqslant 2$. As the data, we use a linearized form of the standard Dirichlet-to-Neumann (DtN) map

[TABLE]

which will be specified below.

Reconstruction of the Schrödinger potential function $c(x)$ at a fixed wavenumber $k=0$ in (1) retrospects back to the original inverse conductivity problem, proposed by Calderón, arising in electrical impedance tomography where uniqueness for the linearization can be proven by using the complex exponential solutions [7]. Later on, the fundamental work by Sylvester and Uhlmann proved the global uniqueness for the inverse potential problem in three and higher dimensions by constructing almost complex exponential solutions, which also yields the global uniqueness of the inverse conductivity problem [17]. In two dimensions, recent work by Bukhgeim [6] demonstrated uniqueness, at any fixed $k$, by introducing stationary phase methods. At $k=0$, a logarithmic stability estimate of recovering the potential function is given by Alessandrini [1] and is further proven to be optimal by Mandache [16]. As for the numerical reconstruction of the potential function (or the conductivity function) at a fixed wavenumber, we refer to [10] and recent work [9]. We shall mention that the inverse medium problem in [2] also highly relates with the current work.

Recently it has been widely observed in different inverse boundary value problems that the stability estimate improves with growing wavenumbers (or energy) both analytically and numerically. The increasing stability in the inverse boundary value problem for the Schrödinger equation (1) was firstly observed within different ranges of the wavenumbers in [12]. In [15], the authors considered the increasing stability on inverse boundary value problems for the Schrödinger equation with attenuation where a linear exponential dependence on the attenuation constant is established. We note that all the above mentioned increasing stability results on inverse Schrödinger potential problems are considered in three and higher dimensions under different a priori regularity assumptions.

Another multifrequency set-up uses single observation for multiple wavenumbers. To have an overview of such results, we recommend a recent review paper [3] which nicely summarizes the theoretical and numerical evidences verifying the increasing stability in inverse medium and source problems for acoustic Helmholtz equations and Maxwell equations.

In this paper, we investigate a linearized form of the inverse Schrödinger potential problem (1) with a large wavenumber in two and higher dimensions. Such a linearization is carried out at a zero potential function which is reasonable if $c$ is sufficiently small compared with the squared wavenumber $k^{2}$. Noticing that the squared wavenumber $k^{2}$ may be large, we are allowed to reconstruct a potential function of moderate amplitude. In Section 2 we introduce the linearized problem and obtain an (increasing) Hölder-type stability estimate for the linearized inverse Schrödinger potential problem with a large wavenumber by using bounded complex exponential solutions. We extend, in Section 3, the discussion to the inverse Schrödinger potential problem with attenuation. An exponential dependence on the attenuation constant in the stability estimate is observed. A novel reconstruction algorithm is proposed to recover the Fourier modes of the potential function in Section 4 based on the Calderón’s method. By choosing the large wavenumber appropriately, we show various numerical examples confirming the efficiency of the proposed algorithm in Section 5. Finally a conclusion Section 6 ends the manuscript with further prospects.

2 Increasing stability in the linearized inverse Schrödinger potential problem at the large wavenumber

We recall that the original problem, which initializes the current work, is to find the Schrödinger potential function $c=c(x)$ defined in a bounded domain $\Omega$ in the following problem

[TABLE]

from the knowledge of the Dirichlet-to-Neumann (DtN) map

[TABLE]

We assume that $k^{2}$ is not a Dirichlet eigenvalue for the Laplace operator in $\Omega$.

If we assume that $c$ is small (or $k$ is large), we can justify the linearization of the Schrödinger equation. More precisely, let $u_{0}$, $u_{1}$ solve the following sub-problems

[TABLE]

then the solution $u$ of the original problem (2) is

[TABLE]

where the remaining ”$\cdots$” are ”higher” order terms. The linearzied DtN map $\Lambda^{\prime}$ of $\Lambda$ in (3) is defined accordingly as

[TABLE]

Multiplying both sides of the sub-problem (4b) with a test function $v\in H^{1}(\Omega)$ satisfying the equation

[TABLE]

we obtain

[TABLE]

Without loss of generality we may assume that $0\in\Omega$. Let $D=2\sup\limits_{x\in\Omega}\left|x\right|$ and $\epsilon$ be the operator norm of $\Lambda^{\prime}:H^{\frac{1}{2}}(\partial\Omega)\to H^{-\frac{1}{2}}(\partial\Omega)$. We present the main stability estimate below:

Theorem 2.1**.**

Let $D\leqslant 1$, $\left\|c\right\|_{H^{1}(\Omega)}\leqslant M_{1}$, and $k>1$, $\epsilon<1$, then the following estimate holds true

[TABLE]

*for the linearized system (4a)-(4b) with $E=-\ln\epsilon$ and the constant $C(\Omega)$ depending on the domain $\Omega$. *

Proof.

The proof is based on the complex exponential solutions suggested by Calderón [7] and Faddeev [11]. Let $\xi\in\mathbb{R}^{n}$ and $\zeta,\zeta^{*}\in\mathbb{C}^{n}$ with $\xi,\zeta,\zeta^{*}\neq 0$. We consider the following solutions

[TABLE]

To select $\zeta$ and $\zeta^{*}$, we introduce an orthonormal base $\left\{e_{1}:=\frac{\xi}{\left|\xi\right|},e_{2},\cdots,e_{n}\right\}$ of $\mathbb{R}^{n}$ and let

[TABLE]

where $\sqrt{k^{2}-\frac{\left|\xi\right|^{2}}{4}}=\mathbf{i}\sqrt{\frac{\left|\xi\right|^{2}}{4}-k^{2}}$ if $k<\frac{\left|\xi\right|}{2}$. For brevity, we denote $\Xi:=\sqrt{\frac{\left|\xi\right|^{2}}{4}-k^{2}}$. Then $u_{0}(x)v(x)=e^{\mathbf{i}\xi\cdot x}$ and from the equality (6) we derive

[TABLE]

Here $\mathcal{F}[\cdot]$ denotes the Fourier transform.

Observe that if $k\geqslant\frac{\left|\xi\right|}{2}$, then the norms of exponential solutions in (7) are bounded as follow

[TABLE]

where $\operatorname{Vol}_{n}\Omega$ is the volume of $\Omega\subset\mathbb{R}^{n}$. If $k<\frac{\left|\xi\right|}{2}$, then

[TABLE]

where $\operatorname{Vol}_{n-1}\Omega=\sup\limits_{\Omega^{\prime}}\{\operatorname{Vol}_{n-1}\Omega^{\prime}\}$ over all $(n-1)$-dimensional orthonormal projections $\Omega^{\prime}$ of $\Omega$.

By the trace theorem, there hold

[TABLE]

So from (8) and (11), we have

[TABLE]

Hence, due to (9) and (10),

[TABLE]

and

[TABLE]

Let $E:=-\ln\epsilon>0$ and $k>1$, $\epsilon<1$, we consider two cases

1. a)

$k>E$ (i.e. $\epsilon=e^{-E}>e^{-k}$), and 2. b)

$k\leqslant E$ (i.e. $\epsilon=e^{-E}\leqslant e^{-k}$).

In the case a), we have

[TABLE]

where $\sigma_{n}$ is the volume of an unit ball in $\mathbb{R}^{n}$, and the constant $C_{1}(\Omega)$ is defined as $C^{4}(\Omega)\left(\operatorname{Vol}_{n}\Omega\right)^{2}\sigma_{n}2^{n+2}$.

In the case b), we let $\rho^{2}:=\frac{E^{2}}{D^{2}}+4k^{2}$ and split

[TABLE]

The second term in the above inequality is derived by (12) and the fact that

[TABLE]

increases while $y>0$ and hence

[TABLE]

Observe that, since $k\leqslant E$,

[TABLE]

Thus

[TABLE]

where $C_{2}(\Omega):=C^{4}(\Omega)\left(\operatorname{Vol}_{n-1}\Omega\right)^{2}\sigma_{n}\frac{2^{2}}{D^{n-2}}\left[\left(1+\left(2D\right)^{2}\right)^{\frac{n}{2}}-\left(2D\right)^{n}\right]$.

Combining case a) $k>E$ with the bound (13) and case b) $k\leqslant E$ with the bound (14), we prove Theorem 2.1.

Remark 2.1**.**

*When $\Omega$ is a sphere, constants in the trace theorem are explicit and one can make $C(\Omega)$ in Theorem 2.1 explicitly. *

In view of (14) it is obvious that stability is increasing with growing wavenumber $k$, provided $k\leqslant E$. We comment on the increasing stability when $k>E$. In such a case, we have an upper bound estimate (13) of $c$ which is

[TABLE]

By an elementary calculus, the upper bound $\omega(k)$ decreases on $(0,k_{*})$ and increases on $(k_{*},\infty)$ where

[TABLE]

If $k_{*}\leqslant 1$ in (15), we have

[TABLE]

and the minimum of $\omega(k)$ is obtained at $k=1$ such that

[TABLE]

If $k_{*}>1$, then minimum of $\omega(k)$ is obtained at $k=k_{*}$ such that

[TABLE]

Now we consider the ”small” perturbation when $E\geqslant 1$ (i.e., $\epsilon\leqslant e^{-1}$). If $k_{*}\leqslant E$, then minimum of $\omega(k)$ (at $k=E$) is

[TABLE]

If $k_{*}>E$, then minimum of $\omega(k)$ is at $k=k_{*}$ which is the same as in the case $0<E<1$. So in all possible cases above, given a small error (or a perturbation) $\epsilon$ we can choose a large wavenumber $k$ producing a Hölder stability which is a big improvement over a logarithmic stability at low wavenumbers.

3 Linearized inverse Schrödinger potential problem at the large wavenumber with attenuation

Furthermore we provide some extended discussion on the linearized inverse Schrödinger potential problem with attenuation. The original problem in three (and higher) dimensions is studied in [15] where increasing stability estimates with a linearly exponential dependence of the attenuation constant are obtained for low/high frequencies separately by constructing almost complex exponential solutions and real solutions respectively.

In current section, we investigate the stability estimate recovering the potential function $c=c(x)$ below

[TABLE]

from the knowledge of linearized DtN map for a large wavenumber $k$. The constant $b>0$ is the attenuation constant. For the sake of simplicity, we use the notations in Section 2, for instance, that $u$ represents the solution of the Schrödinger potential equation with attenuation.

Referring to Section 2, we let $u_{0}$, $u_{1}$ solve the following sub-problems

[TABLE]

and the solution $u$ of (16) is

[TABLE]

with ”$\cdots$” denoting the remaining terms. Similarly the linearized DtN map $\Lambda^{\prime}$ is defined as

[TABLE]

Multiplying both sides of the sub-problem (17b) with a test function $v$ satisfying

[TABLE]

we obtain the equality

[TABLE]

which will again play important roles below.

Let $\epsilon$ be the operator norm of $\Lambda^{\prime}:H^{\frac{1}{2}}(\partial\Omega)\to H^{-\frac{1}{2}}(\partial\Omega)$, $0\in\Omega$ and $D:=2\sup\left|x\right|$, $x\in\Omega$. We present the main stability estimate in this section.

Theorem 3.1**.**

Let $D\leqslant 1$, $\left\|c\right\|_{H^{1}(\Omega)}\leqslant M_{1}$, and $k>1$, $\epsilon<1$, then the following estimate holds true

[TABLE]

*for the linearized system (17a)-(17b) with $E=-\ln\epsilon$ and the constant $C(\Omega)$ depending on the domain $\Omega$. *

Proof.

Similar to the proof of Theorem 2.1, we again use the exponential solutions.

Let $\xi\in\mathbb{R}^{n}$ and $\zeta,\zeta^{*}\in\mathbb{C}^{n}$ with $\xi,\zeta,\zeta^{*}\neq 0$. We consider the following solution

[TABLE]

The orthonormal base is chosen by $\left\{e_{1}=\frac{\xi}{\left|\xi\right|},e_{2},\cdots,e_{n}\right\}$ and we let

[TABLE]

Then $u_{0}(x)v(x)=e^{\mathbf{i}\xi\cdot x}$ and we derive, by (19)

[TABLE]

We can write $\sqrt{k^{2}-\frac{\left|\xi\right|^{2}}{4}-\mathbf{i}kb}=X+\mathbf{i}Y$, $X>0$. If $\left|\xi\right|^{2}\leqslant 3k^{2}$, then squaring the both sides and following the arguments in [15], we have

[TABLE]

If $3k^{2}<\left|\xi\right|^{2}\leqslant 4k^{2}$, we derive

[TABLE]

On the other hand, if $\left|\xi\right|^{2}>4k^{2}$, then

[TABLE]

We then derive that

• •

if $\left|\xi\right|^{2}\leqslant 3k^{2}$, then

[TABLE]

• •

If $3k^{2}<\left|\xi\right|^{2}\leqslant 4k^{2}$, then

[TABLE]

• •

If $\left|\xi\right|^{2}>4k^{2}$, there holds

[TABLE]

where a variable $y:=-Y>0$ due to (20).

Similarly to the proof of Theorem 2.1, we have

[TABLE]

which further yields that

• •

if $\left|\xi\right|^{2}\leqslant 3k^{2}$, then

[TABLE]

• •

If $3k^{2}<\left|\xi\right|^{2}\leqslant 4k^{2}$, then

[TABLE]

• •

If $\left|\xi\right|^{2}>4k^{2}$, then

[TABLE]

Let $E=-\ln\epsilon>0$, $k>1$ and $\epsilon<1$. We again consider the cases below

1. a)

$k>E$ (i.e. $\epsilon=e^{-E}>e^{-k}$), and 2. b)

$k\leqslant E$ (i.e. $\epsilon=e^{-E}\leqslant e^{-k}$).

In the case a), we have

[TABLE]

In the case b), we again choose $\rho^{2}:=\frac{E^{2}}{D^{2}}+4k^{2}$ and split

[TABLE]

As in Theorem 2.1, $\frac{e^{Dy}-e^{-Dy}}{y}$ increases while $y>0$ and thus it attains the maximum value if $\left|\xi\right|=\rho$, which yields

[TABLE]

Since

[TABLE]

and, by using $k\leqslant E$

[TABLE]

we then bound

[TABLE]

Meanwhile, by using $k\leqslant E$ and the proof of Theorem 2.1, we yield

[TABLE]

Combining above inequalities and using $k\leqslant E$ again, we conclude that

[TABLE]

We end the proof by combining the bounds (21) for $k>E$ and (22) for $k\leqslant E$.

Following the discussion in Section 2, we again observe the increasing stability with growing wavenumbers $k$. At the same time, because of the attenuation term, such an error estimate is deteriorated by the attenuation constant $b$ with a linearly exponential growth. Such a dependance is better than that for the multifrequency inverse source problem where a quadratically exponential growth is observed comparing the increasing stability estimates in [8, 14].

4 Reconstruction algorithm

In this section, we propose a reconstruction algorithm for the linearized inverse Schrödinger potential problem at a large wavenumber, observing that the linearized problems are (4a) and (4b). For the sake of simplicity, we fix the dimensionality $n=2$.

The linearized problem (4a)-(4b) represents a linearized DtN map $\Lambda^{\prime}:g_{0}\mapsto\partial_{\nu}u_{1}$ by varying $g_{0}$. In particular, we can choose the solutions $u_{0}$ to be the bounded complex exponential solutions (7) for $n=2$ such that

[TABLE]

by varying different vectors $\xi$ in the phase space and $g_{0}=\left.u_{0}\right|_{\partial\Omega}$, where $\xi^{\perp}$ satisfies $\xi\cdot\xi^{\perp}=0$ and $\left|\xi^{\perp}\right|=1$. At the same time, substituting each $u_{0}$ above into the second problem (4b), we obtain the solution $u_{1}$ and its Neumann trace $\partial_{\nu}u_{1}=\Lambda^{\prime}g_{0}$ which also varies with different choice of $\xi$. Similarly we can choose the test function $v$ to be the relevant exponential solutions

[TABLE]

To realize the reconstruction algorithm numerically, we shall approximate the linearized DtN map $\Lambda^{\prime}$ (5) accurately. In particular, we mention that the Neumann boundary data $\partial_{\nu}u_{1}$ depends on the unknown potential function $c$ referring to (4b). Assuming that the potential function $c$ is comparably small with respect to the squared wavenumber $k^{2}$ and following the standard linearization approach (for instance, in [10]), we define a linearized Neumann boundary data

[TABLE]

where $g_{1}$ is the Neumann boundary data in (3) of the original problem (2) and $\partial_{\nu}u_{0}$ is the Neumann boundary data of the sub-problem (4a). We thus approximate the linearized Neumann boundary data numerically by $\partial_{\nu}u_{1}\approx g_{1}^{\,\prime}$, where the ”higher” order terms are ignored in (25). If we substitute the above solutions (23), (24) into the equality (6), the recovery of $c$ is equivalent to solving the following integral equations

[TABLE]

where vectors $\xi$, $\xi^{\perp}$ vary in the phase space with a fixed wavenumber $k$. The left hand side of the above formula (26) is the Fourier coefficient $\mathcal{F}[c](\xi)$ of the potential function $c$ at a point $\xi$. Thus the reconstruction of $c$ can be achieved by varying vectors $\xi$, $\xi^{\perp}$ and the inverse Fourier transform.

Remark 4.1**.**

Notice that, in current section, the chosen orthonormal base of the phase space $\mathbb{R}^{2}$ is $\left\{e_{1}=\frac{\xi}{\left|\xi\right|},\,e_{2}=\xi^{\perp}\right\}$. Thus the formula of a vector $\zeta$ (or $\zeta^{*}$) is fixed by $k$. Noticing that $\left|\xi\right|>2k$ and $\sqrt{k^{2}-\frac{\left|\xi\right|^{2}}{4}}=\mathbf{i}\sqrt{\frac{\left|\xi\right|^{2}}{4}-k^{2}}$, we write the exponential solution $u_{0}$ (or $v$) explicitly by

[TABLE]

*Indeed, $u_{0}$ is a plane wave traveling along its wave-vector $\zeta$ with a length equal to the wavenumber $k$ if $\left|\xi\right|\leqslant 2k$. Meanwhile, if $\left|\xi\right|>2k$, the (high) oscillation of the exponential solution $u_{0}$ is observed in a direction $e_{1}=\frac{\xi}{\left|\xi\right|}$ with the spatial frequency $\frac{\left|\xi\right|}{2}$, and decays exponentially along the unit vector $e_{2}=\xi^{\perp}$. Then numerical calculation in (6) or (26) becomes unstable when $\left|\xi\right|>2k$. *

Now we describe a reconstruction algorithm based on discrete sets of length and angle of the vectors in phase space. We first choose a discrete and finite length set

[TABLE]

Here we choose $2\leqslant m\in\mathbb{N}_{+}$ and $mk$ is the maximum length of the vector $\xi$. Combining with two unit vector (or angle) sets

[TABLE]

which satisfy that $\hat{y}_{s}\cdot\hat{z}_{s}=0$, we denote

[TABLE]

which are the vectors (or points) chosen in the phase space, while $\ell=1,2,\cdots,M$ and $s=1,2,\cdots,N$. More precisely, the superscript notation $\cdot^{\langle\ell;s\rangle}$ will be referred to a vector $\xi^{\langle\ell;s\rangle}$ with the $\ell$th length $\kappa_{\ell}$ and the $s$th angle $\hat{y}_{s}$. Finally, for the inverse Fourier transform, a numerical quadrature rule can be constructed by a suitable choice of the weights $\sigma^{\langle\ell;s\rangle}$ according to these points $\xi^{\langle\ell;s\rangle}$.

We summarize our reconstruction algorithm below.

   Algorithm 1: Reconstruction Algorithm for the Linearized Schrödinger Potential Problem

  Input: $\{\kappa_{\ell}\}_{\ell=1}^{M}$, $\{\hat{y}_{s}\}_{s=1}^{N}$, $\{\hat{z}_{s}\}_{s=1}^{N}$ and $\sigma^{\langle\ell;s\rangle}$;

Output: Approximated Potential $c_{\rm inv}=c^{\langle M+1;1\rangle}$.

1. 1:

Set $c^{\langle 1;1\rangle}:=0$; 2. 2:

For $\ell=1,2,\cdots,M$ (length updating) 3. 3:

For $s=1,2,\cdots,N$ (angle updating) 4. 4:

Choose $u_{0}:=\exp\{\mathbf{i}\zeta^{\langle\ell;s\rangle}\cdot x\}$ and $g_{0}:=\left.u_{0}\right|_{\partial\Omega}$; 5. 5:

Measure the Neumann boundary data $\partial_{\nu}u$ of the forward problem (2) 6. given the Dirichlet boundary data $g_{0}$; 7. 6:

Calculate the approximated Neumann boundary data $g_{1}^{\,\prime}=\partial_{\nu}u-\partial_{\nu}u_{0}$ 8. on the boundary $\partial\Omega$; 9. 7:

Choose $v:=\exp\{\mathbf{i}\zeta^{\langle\ell;s\rangle}_{*}\cdot x\}$ and w:=\big{[}u_{0}v\big{]}^{-1}=\exp\{-\mathbf{i}\xi^{\langle\ell;s\rangle}\cdot x\}; 10. 8:

Compute $\mathcal{F}[c](\xi^{\langle\ell;s\rangle})\approx\int_{\partial\Omega}g_{1}^{\,\prime}\,v\,\mathrm{d}s_{x}$; 11. 9:

Update $c^{\langle\ell;s+1\rangle}=c^{\langle\ell;s\rangle}+\sigma^{\langle\ell;s\rangle}\mathcal{F}[c](\xi^{\langle\ell;s\rangle})\,w$; 12. 10:

End 13. 11:

Set $c^{\langle\ell+1;1\rangle}:=c^{\langle\ell;N+1\rangle}$; 14. 12:

End.

The above Algorithm 1 provides us with a reconstructed potential function $c_{\rm inv}$ to mimic the exact one. The main error between these two potential functions consists of two parts. The first one is the approximation error because of the discrete and finite length set and angle sets. The second one is the numerical error of approximated linearized DtN map $\Lambda^{\prime}$ in (26) where the elliptic equation solvers of large wavenumbers and the numerical differentiation of the Neumann boundary data may induce some additional ill-posedness. The error estimate of Algorithm 1 may be carried out by choosing the length and angle sets following the analysis in [4, 5]. We intend to report such results in a separate work.

5 Numerical examples

In this section, we provide some numerical examples verifying the efficiency of the above Algorithm 1. The main computational costs in Algorithm 1 is Step 5 for the forward problem (2) and Step 9 for the inverse Fourier transform. To avoid the inverse crime, we will use fine grids (for instance, $100\times 100$ or $200\times 200$ equal-distance points) for the forward problem and a coarse grid ($90\times 90$ equal-distance points) for the inversion both in the square domain $[-1,1]^{2}$.

5.1 CASE 1

In the first case, a circular domain $\Omega$ with a radius $r=0.7$ and a simple true potential function $c$ are chosen, which are shown in the upper left panel of Figure 1. The potential function contains one peak and one valley patterns in the bounded domain $\Omega$. The boundary measurement points on $\partial\Omega$ are also presented in the upper left panel of Figure 1, marked by ”$*$”. The forward problem grid is $100\times 100$ equal-distance points and the inversion grid is $90\times 90$ equal-distance points in the domain $[-1,1]^{2}$. In Subsection 5.1.3, we will show that one can improve the resolution of the reconstructed potential function if a finer grid of the forward problem is chosen.

5.1.1 CASE 1: General recovery

To realize the Fourier transform $\mathcal{F}[\cdot]$, we use the following sampling points $\xi$ in the phase space which are marked by ”$*$” in the upper right panel of Figure 1. The modulus of $\xi$, that is $\kappa=\left|\xi\right|$, are equally distributed in the interval $[1,50]$, and the step size is $0.2$. The degree of the angle between two adjacent slope lines is $\frac{2}{9}\pi$. For example, the red ”$\circ$” in the upper right panel of Figure 1 is a point in the phase space with an angle $\pi$ and a modulus $10.2$. Thus, the coordinate of this point $\xi$ in the phase space is $(-10.2,0)$, or $\kappa=10.2$ and $\hat{y}=(\cos(\pi),\sin(\pi))=(-1,0)$, i.e., $\xi=\kappa\hat{y}$. Fourier coefficients of the true potential function near the nine slope lines are presented in the bottom panel of Figure 1.

By varying different $\hat{y}$ (or corresponding $\hat{z}$) and different $\kappa=\left|\xi\right|$ in the given sets, we obtain different exponential solutions (23) (or (27)). In particular, we present two different $u_{0}$ (real part only) in Figures 2 for $\kappa\leqslant 2k$ and Figure 4 for $\kappa>2k$ while $k=15.2$ is fixed. The real part of the corresponding Dirichlet boundary data $g_{0}$ for the linearized problem are also presented in these two figures. For instance, in Figure 2 (left), the exponential solution $u_{0}$ with a small length $\kappa=8.4$ and $\hat{y}=(-0.17,0.98)$ is a plane wave indeed. In Figure 2 (middle), we present its Dirichlet boundary data $g_{0}$. By substituting the Dirichlet boundary data $g_{0}$ into the original problem (2), we obtain the solution $u$ and its Neumann boundary data $g_{1}=\partial_{\nu}u$, which are shown in Figure 2 (right). In order to obtain the measurable data $g_{1}^{\,\prime}$ of the linearized DtN map (5), we subtract the benchmark solution $u_{0}$ with respect to the solution $u$ in Figure 3 (left), and measure the linearized data $g_{1}^{\,\prime}=g_{1}-\partial_{\nu}u_{0}$ on the boundary $\partial\Omega$, which is shown in Figure 3 (right).

The second exponential solution $u_{0}$ with a large length $\kappa=32.6$ and $\hat{y}=(-0.77,-0.64)$ is shown in Figure 4 (left) where exponential decay in the direction $\hat{z}$ is observed. The Dirichlet boundary data $g_{0}$ and the Neumann boundary data $g_{1}=\partial_{\nu}u$ are shown in Figure 4 (middle and right) respectively. The difference between $u$ and $u_{0}$, the measurable data $g_{1}^{\,\prime}=g_{1}-\partial_{\nu}u_{0}$ of the linearized DtN map are shown in Figure 5. Compared with the moderate linearized difference for $\kappa=8.4$ as shown in Figure 3, we observe a large amplitude of the linearized difference in the right panel of Figure 5 for $\kappa=32.6$. Such phenomena further yields an unstable calculation of the Fourier coefficients when $\kappa>2k$.

Noticing that the exponential solution $v$ in (24) is similar to $u_{0}$, we use the same sampling points $\xi$ in the upper right panel of Figure 1 for $v$. However, to illustrate the difference, for any given wavenumber $k$ and vector $\xi=\kappa\hat{y}$, the directions of two wave-vectors $\zeta=\frac{\kappa}{2}\hat{y}+\sqrt{k^{2}-\frac{\kappa^{2}}{4}}\,\hat{z}$ and $\zeta^{*}=\frac{\kappa}{2}\hat{y}-\sqrt{k^{2}-\frac{\kappa^{2}}{4}}\,\hat{z}$ are not the same. Indeed, since these vectors obey $\hat{y}\cdot\hat{z}=0$, the wave-vectors $\zeta$, $\zeta^{*}$ of exponential solutions $u_{0}$ and $v$ are symmetric with respect to the vector $\xi$ (or $\hat{y}$). Thus, we could reflect the pattern of the exponential solution $u_{0}$ with respect to the vector $\xi$ to obtain the pattern of the solution $v$.

The key step of the reconstruction Algorithm 1 is to compute the value of each Fourier coefficient $\mathcal{F}[c](\xi)$ with different $\xi=\kappa\hat{y}$ for any fixed $k$ by employing the formula (26). The numerical coefficients are collected in Figure 6 when $k=15.2$. Each curve there represents the recovered value of Fourier coefficients on each slope of sampling points $\xi$ in the upper right panel of Figure 1, while the horizontal axis is the length interval $(0,50]$ for $\kappa=\left|\xi\right|$. We shall emphasize that the chosen wavenumber $k=15.2$ avoids the eigenvalues for the Laplacian operator in the domain $\Omega$. On the other hand, by choosing a wavenumber $k$ near the eigenvalues, we may recover a potential function with lower resolution which will be presented in next subsection.

Finally, we implement the inverse Fourier transform to reconstruct the potential function $c$. To obtain such a reconstruction, we need to choose a truncated value $K$ in the phase space such that all the Fourier coefficients with $\kappa\leq K$ shall be used to reconstruct an approximant. By choosing different truncated value $K=mk$ while $m=1,2,3$, we collect those reconstructed functions $c_{\rm inv}$ in Figure 7. More precisely, the dashed lines in Figure 6 highlight those threshold values as $k=15.2$, $2k=30.4$ and $3k=45.6$ when $k=15.2$. The Fourier coefficients $\mathcal{F}[c](\xi)$ with lengths $\kappa=\left|\xi\right|$ smaller than the threshold values, i.e. $\kappa\in(0,mk\,]$, will be used to generate $c_{\rm inv}$ respectively. As one can observe, we obtain a stable reconstructed potential function $c_{\rm inv}$ with threshold values $k$ and $2k$ in Figure 7 (left and middle), because all the Fourier coefficients $\mathcal{F}[c](\xi)$ are rather small for $\kappa\leqslant k$ or $\kappa\leqslant 2k$. The more Fourier coefficients when $\kappa\leqslant 2k$, the better resolution of the reconstructed potential $c_{\rm inv}$. On the other hand, if we increase the truncated value to $3k$, we observe some high frequency patterns in the right panel of Figure 7, which are also reflected in the (blowing-up) Fourier coefficients when $\kappa>2k$, see Figure 6, or referring to the estimate (12) in the proof of Theorem 2.1. In current work, we change different truncated value $K=mk$, $m=1,2,3$ but in real calculation we suggest to choose the truncated value $K=2k$.

We shall note that a large value of the wavenumber $k$ allows us to use more (stable) Fourier coefficients in the phase space when $\kappa\leqslant 2k$. To visualize such difference, we presents the reconstructed potential function $c_{\rm inv}$ in Figure 8 by choosing $k=5$ and the truncated values $K=5,10,15$ respectively.

5.1.2 CASE 1: Recovery at a large wavenumber near eigenvalues

In this subsection, we numerically investigate the consequence by choosing a large wavenumber near eigenvalues for a Laplacian operator in the circular domain $\Omega$.

A particular choice of such a large wavenumber is $k=12.3625$. In Figure 9, we present the reconstructed Fourier coefficients and the reconstructed potential function $c_{\rm inv}$ analogously as in the previous subsection by choosing truncated values $K=2k=24.7250$. Noticing that the high frequency patterns appear in Figure 9 (right), we also observe (blowing-up) Fourier coefficients compared with those in Figures 6.

5.1.3 CASE 1: Finer grids and less Fourier modes

We further focus on another two numerical aspects enhancing or weakening the resolution of the reconstructed potential function $c_{\rm inv}$.

The first aspect is the grid of the forward problem. In the above two subsections, we choose a $100\times 100$ equal-distance points grid in $[-1,1]^{2}$. The chosen grid may capture some solution patterns if the wavenumber $k$ is small. But when the wavenumber $k$ increases, for instance from $5$ to $15.2$, the high frequency patterns induced by large wavenumbers may not be accurately demonstrated in the grid. To weaken such difficulty, we consider a finer $200\times 200$ equal-distance points grid in $[-1,1]^{2}$. By implementing the same approach in the above subsection, we present the Fourier coefficients and their reconstructed potential function $c_{\rm inv}$ in Figure 10 for $k=15.2$ again. As one can observe, the Fourier coefficients of the finer grid become more accurate in the large length interval ($\kappa\in[20,30]$) than those of the original grid. The reconstructed potential functions $c_{\rm inv}$ also enhance the resolution (in the peak pattern) if we choose the truncated value $K=2k=30.4$, see Figure 10 (right).

We also investigate the influence of the observation angles, i.e. the limited-angle slope lines in the upper right panel of Figure 1. In Figure 11, we show the recovered potential functions with $2$, $3$, $7$ different slope lines by choosing the truncated value $K=2k$ and $k=15.2$. Compared with the result in Figure 10, we observe that the more slope lines, the better accurately recovered potential functions, because of the gained Fourier coefficients. In principle, one can include more angle slope lines to obtain better resolution.

5.2 Another numerical example

Another example considers a complicate potential function $c$ in Figure 12 (left) with the same circular domain $\Omega$. Similar to the previous example, the same set of the sampling points $\xi$ in the phase space is utilized referring to the upper right panel of Figure 1. For the sake of simplicity, we skip all the discussion on the reconstruction algorithm but provide the reconstructed potential function $c_{\rm inv}$ with the truncated value $K=2k$ and $k=5.0$, $k=20.0$ respectively in Figure 12 (middle & right). Similarly, we also have chosen a finer $200\times 200$ equal-distance points grid in the domain $[-1,1]^{2}$ for the forward problem. Comparing the reconstructed potential functions $c_{\rm inv}$ for $k=5$ and $k=20$, we clearly visualize the improved resolution at the large wavenumber.

6 Conclusion

The original inverse problem for the Schrödinger potential is ill-posed and nonlinear (moreover, non-convex). To study stability we prefer to focus on the most serious difficulty: ill-conditioning and we linearized inverse problem to avoid additional difficulties with multiple local minima. A numerical solution of the linearized problem is much faster and more reliable and is quite satisfactory in many applications.

We demonstrated that recovery of the Schrödinger potential from all boundary data at a fixed energy/wavenumber is dramatically improving for a larger $k$. There is no doubt that this improvement will be more significant when finding more complicated $c$ or using a complete nonlinear problem.

These results promise a better numerical reconstruction of the conductivity coefficient in the stationary Maxwell system at higher wavenumbers, as predicted analytically in [13]. We intend to work on the reconstruction at least in the linearized version in near future. There are good expectations that such numerical results will be a solid base for a serious improvement in the electrical impedance tomography which now suffers from a very low resolution. As known, this type of tomography has many geophysical and medical applications.