Numerical solution of inverse evolution problems via the nonlinear

Levitan equation

Frank Natterer

April 1999

1 Introduction

We consider a second order evolution problem

u_{tt}^{φ} = Lu^{φ} , u^{φ}(0) = 0 , u_{t}^{φ}(0) = φ  (1.1)

where L is a selfadjoint operator in a Hilbert space H. The problem is to recover L (i.e. all or some of its coefficients) from the observations

g_{φ,ψ}(t) = (ψ, u^{φ}(t)) , 0 ≤ t ≤ T ,  (1.2)

for sufficiently many φ, ψ. We shall show that we have the identity

(u^{ψ}(s), u^{φ}(t)) = G_{φ,ψ}(s, t) , 0 ≤ s ≤ t ,  (1.3)

G_{φ,ψ}(s, t) = \frac{1}{2} ∫_{0}^{s} (g_{φ,ψ}(t + s′) + g_{φ,ψ}(t − s′)) ds′ .

In the 1D case and L = \frac{d^{2}}{dx^{2}} −q(x), (1.3) will turn out to be Levitan’s nonlinear integral equation. We shall extend this to the case L = \frac{d}{dx}(c^{2} \frac{d}{dx}) − q and suggest a numerical method based on the Cholesky decomposition. Finally we indicate how (1.3) can be used for the numerical solution of nD problems.

2 Derivation of the identity

Let E_{0}(t), E_{1}(t) be solution operators for (1.1), i.e. v_{tt} = Lv with initial values v(0) = φ_{0}, v_{t}(0) = φ_{1} has the solution

v(t) = E_{0}(t)φ_{0} + E_{1}(t)φ_{1} .

We extend E_{0}, E_{1} to all of IR^{1} as even and odd functions, respectively. We have E′_{1} = E_{0}, and L^{∗} = L implies E_{0}^{∗} = E_{0} and E_{1}^{∗} = E_{1}.

Now let 0 ≤ s ≤ t. Exploiting the time invariance of (1.1) we get

g_{φ,ψ}(t + s) = (ψ, u^{φ}(t + s))
= (ψ, E_{0}(s)u^{φ}(t) + E_{1}(s)u_{t}^{φ}(t))
= (E_{0}(s)ψ, u^{φ}(t)) + (E_{1}(s)ψ, u_{t}^{φ}(t))
= (E_{0}(s)ψ, E_{1}(t)φ) + (E_{1}(s)ψ, E′_{1}(t)φ) .

Likewise, since E_{0} is even and E_{1} is odd,

g_{φ,ψ}(t − s) = (E_{0}(s)ψ, E_{1}(t)φ) − (E_{1}(s)ψ, E′_{1}(t)φ) .

Adding the last two equations yields

\frac{1}{2}(g_{φ,ψ}(t + s) + g_{φ,ψ}(t − s)) = (E_{0}(s)ψ, E_{1}(t)φ) .

The identity (1.3) is now obtained by integrating with respect to s.

We remark that for a non-selfadjoint operator L (1.3) reads

(v^{ψ}(s), u^{φ}(t)) = G_{φ,ψ}(s, t)

where v^{ψ} is the solution of

v_{tt}^{ψ} = L^{∗}v^{ψ} , v^{ψ}(0) = 0 , v_{t}^{ψ}(0) = ψ .

Thus the method of the following sections can be carried out also in the non-selfadjoint case, with the Cholesky decomposition replaced by the LUfactorization.

3 The nonlinear Levitan equation

Now consider L to be the differential operator

Lu = u′′ − q(x)u in (0, ∞)

with boundary conditions u′(0) = 0 and let φ = ψ = δ. In slightly simplified notation (1.3) reads

∫_{0}^{s} u(s, x)u(t, x)dx = G(s, t) , s ≤ t ,  (3.1)

where now g(t) = u(0, t). The restriction of the integration to [0, s] is due to causality, i.e. u(x, t) = 0 for x > t. (3.1) is clearly the well known nonlinear Levitan equation. It leads immediately to a numerical algorithm: Let

x_{i} = iΔx , t_{ℓ} = ℓΔt , u_{iℓ} = u(t_{ℓ}, x_{i}) , G_{kℓ} = \frac{1}{Δx}G(t_{ℓ}, t_{k}) .

Let U = (U_{iℓ})_{i,ℓ=0,...,n}, G = (G_{kℓ})_{k,ℓ=0,...,n}. Then, (3.1) reads approximately

U^{∗}U = G .  (3.2)

For Δt = Δx, U is upper triangular. From Kirchhoff’s formula we get

u(x, t) = \frac{1}{2}(H(x + t) − H(x − t)) − \frac{1}{2} ∫_{Δ(x,t)} qudx′dt′

with H the Heavyside function and Δ(x, t) the triangle with vertex (x, t) and base (x − t, x + t) on the x-axes. It follows that u = 1/2 along x = t. Hence the diagonal of U in positive. Thus, U can be computed from (3.2) simply by a Cholesky decomposition. This leads to an approximation to the function u, which in turn provides an approximation to q by the differential equation. Of course all this is well known, see e.g. [2].

This procedure can easily be extended to the differential operator

Lu = (c^{2}(x)u′)′ − q(x)u .

Assume we know an upper bound \overline{c} for c, i.e. 0 < c(x) ≤ \overline{c}. Then, by causality, u(x, t) = 0 for |x| > \overline{c}t. Choosing Δx = \overline{c}Δt leads again to an upper triangular matrix U, and (3.2) is still valid. However, the diagonal of U can now longer be assumed to be positive, and the Cholesky process may break down. Fortunately we have the following

Lemma 3.1 Let U be an upper triangular matrix with the property that u_{k,ℓ} = 0 implies u_{k+1,ℓ+1} = 0 and the non vanishing element in each column (if any) with largest row index is positive. Then U is uniquely determined by G = U^{∗}U.

Proof: Let m_{ℓ} = Max{k : u_{k,ℓ} ≠ 0}, ℓ = 1, …, n, and m_{ℓ} = 0 if column ℓ of U is zero. Since U is upper triangular, m_{ℓ} ≤ ℓ. Since u_{k+1,ℓ+1} = 0 when ever u_{k,ℓ} = 0, m_{ℓ+1} ≤ m_{ℓ} + 1. Moreover, u_{m_{ℓ},ℓ} > 0 for ℓ = 1, …, n. We show by induction on ℓ that the rows m_{ℓ}, ℓ = 1, …, n, of U are uniquely determined by G.

Let ℓ = 1. If m_{ℓ} = 0 there is nothing to prove. If m_{ℓ} = 1 then u_{1,1}^{2} = g_{1,1} and u_{1,1} is determined since u_{1,1} > 0. The rest of row m_{1} is now determined by u_{1,ℓ} = g_{1,ℓ}/u_{1,1}, ℓ = 2, …, n. Thus row m_{1} is determined.

Assume that rows m_{1}, …, m_{ℓ−1} are already determined for some ℓ ≤ n. The equation for g_{ℓ,ℓ} reads

∑_{k=1}^{m_{ℓ−1}} u_{k,ℓ}^{2} + ∑_{k=m_{ℓ−1}+1}^{m_{ℓ}} u_{k,ℓ}^{2} = g_{ℓ,ℓ} .

The first of these sums is known. Because of m_{ℓ} ≤ m_{ℓ−1} + 1 the second sum is either empty (this is the case m_{ℓ} < m_{ℓ−1} + 1) or consists of the single term u_{m_{ℓ},ℓ}^{2} (case m_{ℓ} = m_{ℓ−1} + 1)). In the former case, m_{ℓ} ≤ m_{ℓ−1}, hence row m_{ℓ} is already known. In the latter case, u_{m_{ℓ},ℓ} is uniquely determined since u_{m_{ℓ},ℓ} > 0. The remaining elements in row m_{ℓ} are computed from the elements g_{m_{ℓ},k} exactly as in the Cholesky algorithm.

□

The condition for the assumptions of the lemma to hold is essentially that the discretization be sufficiently fine. The solution u(x, t) is continuous for t ≥ τ(x) where τ(x) is the travel time of a signal travelling from 0 to x, and u(x, τ(x) + 0) > 0. Since τ is continuous (provided c is), U satisfies the assumptions of the lemma if the discretization is sufficiently fine.

4 Extension to higher dimensions

We start with the simple example

\begin{aligned}u_{tt} & = & Δu − q(x)u & in & Ω × [0, T] \\ \frac{∂u}{∂ν} & = & 0 & in & ∂Ω × [0, T] . \end{aligned}  (4.1)

The function q has to be reconstructed from the observations

g_{y,z}(t) = u^{y}(z, t) , y, z ∈ ∂Ω , 0 ≤ t ≤ T

where u^{y} is the solution of (4.1) with u^{y} = 0, u_{t}^{y} = δ_{y} for t = 0 where δ_{y} is the δ-function on ∂Ω concentred at y. It is clear that this problem is a special case of (1.1), (1.2). It is not a difficult problem since it can be reduced to the Radon problem by the progressive wave expansion, see e.g. [1]. We use it only for the explanation of our method in higher dimensions.

With a slight abuse of notation, (1.3) reads in this case

∫_{Ω} u^{z}(x, s)u^{y}(x, t)dx = G_{z,y}(s, t) , 0 ≤ s ≤ t , z, y ∈ ∂Ω .  (4.2)

We note in passing that (4.2) is essentially also used in [3], with the difference that in that paper z, y refer to arbitrary functions rather than δ functions on ∂Ω. G is determined in that paper from the Dirichlet-to-Neuman map by solving (\frac{∂^{2}}{∂t^{2}} − \frac{∂^{2}}{∂s^{2}} )G = 0, while we consider G as data.

To fix ideas let Ω be the disk with midpoint 0 and radius ρ in IR^{2}. Introduce the polar coordinate grid

x_{iν} = r_{i} ( \begin{matrix}cos θ_{ν} \\ sin θ_{ν} \end{matrix} ) , r_{i} = (n − i)Δr , θ_{ν} = νΔθ ,

Δr = \frac{ρ}{n} , Δθ = \frac{2π}{q} , i = 0, …, n − 1 , ν = 0, …, q − 1 .

Define the q × q matrices U_{iℓ} by

(U_{iℓ})_{μν} = u^{x_{0ν}}(x_{iμ}, t_{ℓ})\sqrt{r_{i}ΔrΔθ}

where t_{ℓ} = ℓΔt, ℓ = 0, …, n − 1 and put U = (U_{iℓ})_{i,ℓ=0,...,n−1}. U is a (nq, nq) matrix, arranged as a (n, n)-matrix of (q, q) blocks. The same applies to the matrix

G = (G_{kℓ})_{k,ℓ=0,...,n−1} , (G_{kℓ})_{νμ} = G_{x_{0ν},x_{0μ}}(t_{ℓ}, t_{k}) .

Applying a quadrature rule to (4.2) we obtain

ΔrΔθ ∑_{i=0}^{n−1} ∑_{λ=0}^{q−1} r_{i}u^{x_{0μ}}(x_{iλ}, t_{ℓ})u^{x_{0ν}}(x_{iλ}, t_{k}) = G_{x_{0μ},x_{0ν}}(t_{ℓ}, t_{k}) .

In terms of the quantities just introduced this reads

∑_{i=0}^{n−1} ∑_{λ=0}^{q−1}(U_{iℓ})_{λμ}(U_{ik})_{λν} = (G_{kℓ})_{νμ}

or

∑_{i=0}^{n−1} U_{iℓ}^{∗}U_{ik} = G_{kℓ}

or

U^{∗}U = G .  (4.3)

By causality, u^{y}(x, t) = 0 for |y − x| > t. Hence U_{iℓ} = 0 for i > ℓ provided Δr = Δt. In this case, U_{ℓℓ} is a diagonal matrix, again by causality. Hence (4.3) is simply a block Cholesky decomposition of G which determines U uniquely.

The extension to differential equations of the type

u_{tt} = div(c^{2}∇u) − q(x)u

is done exactly in the same way as in the 1D case. Putting Δr = \overline{c}Δt yields U_{iℓ} = 0 for i > ℓ, hence U is block upper triangular. However it is not clear what the analogue of the lemma is. Before this problem is dealt with it seems advisable to get numerical experience with the 1D case first.

We did not discuss the determination of c, q from the differential equation once u is determined. The exact extent to which this is possible is not yet clear.

5 First order evolutions

Of course the same can be done for first order evolutions

u_{t}^{φ} = Lu^{φ} , u^{φ}(0) = φ .  (5.1)

With g_{φ,ψ} as in (1.2) we readily see that the identity (1.3) now reads

(u^{ψ}(s), u^{φ}(t)) = g_{φ,ψ}(s + t) .  (5.2)

The identity (5.2) is now used exactly in the same way as (1.3) to solve inverse problems for (5.1). However, since inverse first order evolution problems are notoriously unstable, the practical usefulness of (5.2) is probably very limited.

References

[1] Rose, J.H. - Cheney, M. - DeFacio, B.: The connection between time - and frequency domain three-dimensional inverse scattering methods. J. Math. Phys. 25, 2995-3000 (1984).

[2] Bube, K.P. - Burridge, R.: The one-dimensional inverse problem of reflection seismology, SIAM Review 25, 497-559 (1983).

[3] Belishev, M.I.: On an approach to multidimensional inverse problems for the wave equation, Soviet Math. Dokl. 36, 481-484 (1988).

0LaTeX: berg f berger paperfne PapFNFeb97.tex