Vorlesungsskript Inverse Probleme partieller

Differentialgleichungen

F. Natterer Institut für Numerische

und instrumentelle Mathematik

SS 1997, Stand: September 2000

Inhaltsverzeichnis

1 Einleitung 3

2 Die Gelfand-Levitan-Theorie 8

2.1 Hyperbolische Differentialgleichungen in einer räumlichen Va-

riablen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 Das inverse Problem einer hyperbolischen Gleichung in einer

räumlichen Variablen . . . . . . . . . . . . . . . . . . . . . . . 14 2.3 Ein Beispiel aus der

Reflektions - Seismologie . . . . . . . . . . . . . . . . . . . . . 20 2.4 Näherungslösung für kleine q . . . . . . . . . . . . . . . . . . . 22 2.5 Diskrete inverse Evolutionsprobleme . . . . . . . . . . . . . . . 23

3 Näherungen für das inverse Problem der Wellengleichung 29 3.1 Klassische Migration (Hagedorn 1954) . . . . . . . . . . . . . 30 3.2 Migration durch die Wellengleichung . . . . . . . . . . . . . . 33 3.3 Fourier-Transformation und

Distributionen . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.4 Asymptotische Lösung der

Wellengleichung . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.5 Die Bornsche Näherung . . . . . . . . . . . . . . . . . . . . . . 46 3.6 Die “distorted” Born-Approximation . . . . . . . . . . . . . . 52

4 Integral-Geometrie 55 4.1 Die Radon-Transformation . . . . . . . . . . . . . . . . . . . . 55 4.2 Rekonstruktion einer Funktion

aus sphärischen Mitteln . . . . . . . . . . . . . . . . . . . . . 60 4.3 Lösung integralgeometrischer Aufgaben mit Symmetrien … 65

4.4 Näherungsweise Inversion verallgemeinerter Radon-Transformation 69 5 Numerische Methoden 70

5.1 Das Kaczmarz-Verfahren . . . . . . . . . . . . . . . . . . . . . 70 5.2 Das Kaczmarz-Verfahren in der

Impedanz-Tomographie . . . . . . . . . . . . . . . . . . . . . . 76 5.3 Layer - Stripping . . . . . . . . . . . . . . . . . . . . . . . . . 79

Literaturverzeichnis 85

Kapitel 1

Einleitung

Gegenstand der Vorlesung sind inverse Probleme bei gewöhnlichen und partiellen Differentialgleichungen. Während beim direkten Problem die Lösung aus der gegebenen Differentialgleichung zu berechnen ist, versucht man beim inversen Problem die Differentialgleichung aus Informationen über ihre Lösung zu bestimmen. Wir erläutern dies an Hand einiger Beispiele.

1. Das inverse Streuproblem der Akustik

Beim direkten Streuproblem sucht man eine Lösung u_{θ} der Differentialgleichung

Δu_{θ} + k^{2}u_{θ} = qu_{θ}  ,  q = q(x) gegeben,

der Form

u_{θ}(x) = e^{ikx·θ} + v_{θ}(x) ,

wobei v_{θ} der Sommerfeldschen Ausstrahlungsbedingung

\frac{∂}{∂r}v_{θ}(x) − ikv_{θ}(x) → 0 , |x| = r → ∞

(in geeignetem Sinne) genügt. Dabei beschreibt das “Potential” q die akustischen Eigenschaften eines beschränkten Objekts. Außerhalb des Objekts sei q = 0. Wird dieses Objekt aus der Richtung θ mit der Frequenz k beschallt, so entsteht die gestreute Welle v_{θ}.

Beim inversen Problem möchte man q innerhalb des Objekts aus den Werten von v_{θ} für alle Richtungen θ außerhalb des Objektes bestimmen. Dies ist z.B. das mathematische Problem der Ultraschall-Tomographie.

2. Impedanz-Tomographie

Ein beschränktes Objekt Ω besitze die Leitfähigkeit σ. Dann gilt für das Potential u einer stationären Stromverteilung in Ω

div(σ∇u) = 0 in Ω .  (0.1)

Beim direkten Problem berechnet man u (bis auf eine additive Konstante) aus den Randwerten

j = σ \frac{∂u}{∂ν}  auf ∂Ω .  (0.2)

Man hat also eine Neumannsche Randwertaufgabe zu lösen.

Beim inversen Problem nimmt man an, man kenne für jedes j auf ∂Ω das Potential u, welches (0.1), (0.2) genügt, auf ∂Ω. Mit anderen Worten, man kennt die Neumann-zu-Dirichlet-Abbildung Λ_{σ}, welche folgendermaßen erklärt ist:

Schritt 1 : Löse (0.1), (0.2) für vorgegebenes j .
Schritt 2 :  Λ_{σ}j = u|_{∂Ω}  .

Technisch wird das so gemacht: Man appliziert mit Hilfe von Elektroden auf ∂Ω alle möglichen Stromdichten j auf ∂Ω und mißt die zugehörigen Spannungen Λ_{σ}j. Aus Λ_{σ} ist dann σ zu berechnen.

3. Inverse Probleme bei chemischen Reaktionen

In einem Gefäß der Länge L seien zwei Substanzen mit den Konzentrationen u(x, t), v(x, t) an der Stelle x ∈ [0, L] und der Zeit t enthalten, welche gemäß

\frac{∂u}{∂t} = D\frac{∂^{2}u}{∂x^{2}} − kuv , \frac{∂v}{∂t} = −kuv ,  0 ≤ x ≤ L , t ≥ 0 (0.3)

miteinander reagieren. Dabei bedeutet D den Diffusionskoeffizienten, k die Reaktionsgeschwindigkeit. Das Gefäß sei bei x = L für die erste Substanz undurchdringlich, bei x = 0 sei die erste Substanz nicht vorhanden, also

u(0, t) = 0 , \frac{∂u}{∂x}(L, t) = 0 , t ≥ 0 .  (0.4)

Zur Zeit t = 0 seien beide Konzentrationen 1, also

u(x, 0) = 1 , v(x, 0) = 1 , 0 ≤ x ≤ 1 .  (0.5)

Das direkte Problem besteht darin, bei bekanntem D, k das (nichtlineare, parabolische) Anfangsrandwertproblem (0.3 - 0.5) zu lösen. Beim inversen Problem bestimmt man aus Information über u, v, z.B. der Gesamtmasse von u, also

h(t) = ∫_{0}^{L} u(x, t)dx ,

die Parameter k, D.

4. Ein inverses Eigenwertproblem der Seismologie

Wir betrachten die Erde als Kugel vom Radius R und nehmen an, daß ihre Dichte ρ und Steifigkeit (rigidity) μ im Mantel r = |x| > R_{c} nur von r abhängen. Dann genügen die Eigenschwingungen der Erde der Differentialgleichung

− (r^{4}μu′)′ + \frac{(ℓ + 2)(ℓ − 1)}{r^{2}}r^{4}μu = λr^{4}ρu , R_{c} ≤ r ≤ R ,  (0.6)
u(R_{c}) = u(R) = 0 .

Hier ist ℓ = 0, 1, . . ..

Man kann (unter geeigneten Voraussetzungen an μ, ρ) zeigen, daß dieses Problem nur für diskrete Werte λ_{ℓj}, j = 0, 1, … Lösungen u hat, die nicht identisch verschwinden. λ_{ℓj} heißen Eigenwerte von (0.6), und \sqrt{λ_{ℓj}/2π} sind die Frequenzen der Eigenschwingungen.

Das direkte Problem besteht darin, bei bekannten μ, ρ die λ_{ℓj} zu berechnen. Praktisch interessanter ist das inverse Problem: Aus den Eigenfrequenzen \sqrt{λ_{ℓj}/2π}, die man bei jedem Erdbeben messen kann, die Funktionen μ, ρ zu berechnen.

5. Das inverse Streuproblem der Quantenmechanik

Die Schrödinger-Gleichung für das Potential V lautet

− Δψ + V ψ = Eψ .  (0.7)

Unter gewissen Voraussetzungen an V kann man zeigen, daß es für jede Richtung θ und jede Frequenz k mit k^{2} = E eine Lösung ψ_{θ}(k, x) gibt mit

ψ_{θ}(k, x) = e^{ikθ·x} + \frac{e^{ik|x|}}{|x|} A(k, \frac{x}{|x|}, θ) + O (\frac{1}{|x|^{2}} ) .

A ist das Fernfeld, das bei Einstrahlung einer Welle aus Richtung θ mit Frequenz k entsteht. Das direkte Problem besteht in der Berechnung von A aus V. Beim inversen Problem möchte man V aus A berechnen.

6. Das inverse kinematische Problem

Sei Ω ein Objekt mit dem Brechungsindex n. Für x_{0}, x_{1} ∈ ∂Ω sei T (x_{0}, x_{1}) die Zeit, die ein Signal von x_{0} nach x_{1} benötigt. Dieses Signal läuft entlang der Geodätischen Γ_{n}(x_{0}, x_{1}) bezüglich der Metrik

ds = n\sqrt{∑_{i=1}^{3}(dx_{i})^{2}} ;

diese ist Lösung des Systems

ẍ = \frac{n′}{n^{3}} \frac{x}{|x|} − \frac{2n′}{n} \frac{x · ẋ}{|x|}ẋ , n′ = ∇n · \frac{x}{|x|} .  (0.8)

Das direkte Problem besteht darin, für vorgegebenes n die Laufzeit

T (x_{0}, x_{1}) = ∫_{Γ_{n}(x_{0},x_{1})} n\sqrt{∑_{i=1}^{3}(dx_{i})^{2}}

für jedes Paar x_{0}, x_{1} ∈ ∂Ω zu berechnen. Das inverse Problem verlangt die Berechnung von n aus der Funktion T (x_{0}, x_{1}).

Die Aufgaben des Mathematikers sind die folgenden:

1) Man beweise, daß das inverse Problem eindeutig lösbar ist.

2) Man untersuche, wie die Lösung von den Daten abhängt (Stabilität).

3) Man entwickle numerische Verfahren.

Die Geschichte der inversen Probleme beginnt mit der Lösung des inversen kinematischen Problems für die Erde durch Herglotz und Wiechert in den Jahren 1905-07. Der Aufbau der mathematischen Theorie beginnt mit der Arbeit von Borg zum inversen Sturm-Liouville-Problem 1946. Ab 1950 bemühten sich vor allem russische Mathematiker (I.M. Gelfand, B.M. Levitan, V.A. Marchenko) um eine Lösung des quantenmechanischen inversen Streuproblems. Etwa ab 1970 setzten weltweit intensive Bemühungen zum mathematischen Verständnis, aber auch zur numerischen Lösung inverser Probleme ein. Treibende Kraft waren dabei die Anforderungen aus der Praxis, insbesondere aus Medizin, Materialprüfung, Seismologie, Fernerkundung, Exploration.

Die Gelfand-Levitan-Theorie

2.1 Hyperbolische Differentialgleichungen in einer räumlichen Variablen

Wir betrachten das Anfangswertproblem

\frac{∂^{2}u}{∂t^{2}} = \frac{∂^{2}u}{∂x^{2}} + q(x)u , x ∈ IR^{1} , t ≥ 0 ,  (1.1)
u(x, 0) = f(x) , u_{t}(x, 0) = g(x) , x ∈ IR^{1} .  (1.2)

Dabei seien q, f, g zunächst Funktionen aus C_{0}^{∞}(IR^{1}).

Satz 2.1.1 Das Problem (1.1)-(1.2) besitzt eine eindeutig bestimmte Lösung u ∈ C^{1}(IR × [0, ∞)) ∩ C^{2}(IR × (0, ∞)).

Beweis: Jede Lösung der genannten Art kann mit Hilfe der Kirchhoffschen Formel (siehe PDE, Kap. 4.5) in der Form

u(x, t) = \frac{1}{2}(f(x + t) + f(x − t)) + \frac{1}{2} ∫_{x−t}^{x+t}g(y)dy

(1.3)

+\frac{1}{2} ∫_{0}^{t} ∫_{x−(t−τ)}^{x+t−τ} q(y)u(y, τ)dydτ

geschrieben werden. Sei X der Banachraum C(IR×[0, t_{0}]) mit t_{0} > 0 und der Norm des maximalen Betrags, und sei T : X → X der Operator

(T v)(x, t) = \frac{1}{2}(f(x + t) + f(x − t)) + \frac{1}{2} ∫_{x−t}^{x+t} g(y)dy

+\frac{1}{2} ∫_{0}^{t} ∫_{x−(t−τ)}^{x+t−τ} q(y)v(y, τ)dydτ .

Für u, v ∈ X gilt dann

|(T v − T u)(x, t)| ≤ \frac{1}{2} ∫_{0}^{t} ∫_{x−(t−τ)}^{x+t−τ} |q(y)||v(y, τ) − u(y, τ)|dydτ ,

also

||T v − T u|| ≤ Q||v − u|| , Q = \frac{1}{2} t_{0}^{2} max_{IR^{1}} |q(x)| .

Ist t_{0} so klein, daß Q < 1 gilt, so ist T also kontrahierend. u = T u besitzt dann eine eindeutig bestimmte Lösung in X. Dieses u löst dann (1.3) und ist die Lösung der behaupteten Art zunächst in dem Streifen 0 < t < t_{0}. Entsprechend schließt man in den anschließenden Streifen.

□

Wir interessieren uns für Anfangswerte, welche einen Impuls darstellen, also z.B. f = 0 und g = δ mit der Dirac-Funktion δ. Dies ist eine Funktion mit den Eigenschaften δ(x) ≥ 0, δ(x) = 0 für x ≠ 0, und

∫_{IR^{1}} δ(x)dx = 1 .

Dies widerspricht natürlich dem üblichen Begriff der Funktion. Man kann δ aber im Rahmen der Theorie der Distribution befriedigend definieren. Wir wollen diese Theorie weitgehend vermeiden und definieren δ durch einen Grenzprozeß. Sei δ_{1} ∈ C_{0}^{∞}(IR^{1}), δ_{1}(x) ≥ 0, δ_{1}(x) = 0 für x ∉ [−1, +1] und

∫_{IR^{1}} δ_{1}(x)dx = 1 .

Für n = 1, 2, … setzen wir δ_{n}(x) = nδ(nx). Dann ist δ_{n}(x) ≥ 0, δ_{n}(x) = 0 für x ∉ [−1/n, 1/n], und

∫_{IR^{1}} δ_{n}(x)dx = 1 .

Für jede stetige Funktion f gilt nun

lim_{n→∞} ∫_{IR^{1}} δ_{n}(x)f(x)dx = f(0) .

Hierfür schreiben wir

∫_{IR^{1}} δ(x)f(x)dx = f(0) ,

setzen also formal δ = lim_{n} δ_{n} und vertauschen Limesbildung und Integration. δ erscheint hier also nicht als Funktion, sondern als lineares Funktional auf dem linearen Raum der stetigen Funktionen.

Nach Satz 1 ist die Lösung u_{n} von (1.1) mit Anfangswerten u_{n}(x, 0) = 0, \frac{∂u_{n}}{∂t} (x, 0) = δ_{n}(x) eindeutig bestimmt.

Satz 2.1.2 Für |x| ≠ t konvergiert u_{n} gegen eine Lösung u von (1.1), und zwar gleichmäßig auf jedem Kompaktum. Es ist u = 0 für |x| > t und u(|x|, |x| + 0) = 1/2 für x ≠ 0.

Beweis: Nach (1.3) erfüllt u_{n}

u_{n}(x, t) = \frac{1}{2} ∫_{x−t}^{x+t}δ_{n}(y)dy + \frac{1}{2} ∫_{0}^{t} ∫_{x−(t−τ)}^{x+(t−τ)}q(y)u_{n}(y, τ)dydτ .  (1.4)

Sei nun K = \{(x, t) : 0 ≤ t < |x| − ϵ, t ≤ t_{0}, |x| ≤ \frac{1}{ε}\}, ϵ > 0. Für hinreichend großes n ist dann

∫_{x−t}^{x+t}δ_{n}(y)dy = 0 , (x, t) ∈ K .

Sei weiter t_{0} so klein, daß Q (siehe Beweis zu Satz 1) < 1. Sei X = C(K). Mit der Norm in X gilt dann ∥u_{n}∥ ≤ Q∥u_{n}∥, also u_{n} = 0 in K. Wie im Beweis zu Satz 1 muß auch diese Argumentation für die “Streifen” mit t > t_{0} wiederholt werden. In jedem Kompaktum in t < |x| gilt also von einem gewissen n an u_{n} = 0. Also haben wir u_{n} → 0 in t < |x| mit gleichmäßiger Konvergenz auf kompakten Teilmengen.

Setzen wir hingegen K = {(x, t) : t_{0} ≥ t > |x| + ε}, ε > 0, so gilt für hinreichend großes n

∫_{x−t}^{x+t}δ_{n}(y)dy = 1 , (x, t) ∈ K .

Damit ist für hinreichend große n in K

u_{n}(x, t) = \frac{1}{2} + \frac{1}{2} ∫_{0}^{t} ∫_{x−(t−τ)}^{x+t−τ} q(y)u_{n}(y, τ)dydτ .  (1.5)

Sei nun wieder t_{0} so klein, daß Q < 1. Dann folgt aus (1.4), daß die Folge u_{n} in X = C(IR^{1} × [0, t_{0}]) beschränkt ist. Nach Arzela-Ascoli enthält der Integralterm in (1.5) eine in K konvergente Teilfolge. Wie üblich schließt man, daß die gesamte Folge in K gleichmäßig konvergieren muß. Die Grenzfunktion u erfüllt in K natürlich ebenfalls (1.5). Also ist u Lösung von (1.1). Die Beziehung u(|x|, |x| + 0) = \frac{1}{2} ergibt sich aus u = 0 für |x| > t.

□

Definition: Die Funktion u aus Satz 2 nennen wir die Lösung der Anfangswertaufgabe

\frac{∂^{2}u}{∂t^{2}} = \frac{∂^{2}u}{∂x^{2}} + q(x)u , x ∈ IR^{1} , t > 0

(1.6)

u(x, 0) = 0 , \frac{∂u}{∂t}(x, 0) = δ(x) , x ∈ IR^{1} .

Mit Satz 2 ist das direkte Problem (1.6) erledigt: Zu vorgegebenem q gibt es eine wohlbestimmte Lösung u von (1.6).

Die Lösung aus Satz 2 ist in t > |x| nichts anderes als die Lösung der charakteristischen Anfangswertaufgabe, vgl. PDE, Kap. 4.6, mit den Anfangswerten u = 1/2 auf den Charakteristiken t = |x|. u läßt sich auch als schwache Lösung von (1.1) in t > 0 auffassen. Damit meint man folgendes. Sei L ein linearer Differentialoperator in einem Gebiet Ω ⊆ IR^{n}, z.B.

Lu = ∑_{i,j=1}^{n} a_{ij} \frac{∂^{2}u}{∂x_{i}∂x_{j}} + ∑_{i=1}^{n} b_{i} \frac{∂u}{∂x_{i}} + cu .

Sei L^{∗} der zu L adjungierte Operator, also

L^{∗}u = ∑_{i,j=1}^{n} \frac{∂^{2}}{∂x_{i}∂x_{j}} (a_{ij}u) − ∑_{i=1}^{n} \frac{∂}{∂x_{i}} (b_{i}u) + cu .

Dann ist für φ ∈ C_{0}^{∞}(Ω)

∫_{Ω} φLudx = ∫_{Ω} (L^{∗}φ)udx ,

vgl. PDE Kap. 4.6. Ist also u ∈ C^{2}(Ω) eine Lösung von Lu = f in Ω, so gilt

∫_{Ω} (L^{∗}φ)udx = ∫_{Ω} φfdx ,  ∀ φ ∈ C_{0}^{∞}(Ω) .  (1.7)

Diese Beziehung macht Sinn für u ∈ L_{loc}(Ω), die Menge der in Ω lokal integrierbaren Funktionen. Wir nennen daher ein u ∈ L_{loc}(Ω) mit (1.7) schwache Lösung von Lu = f.

Satz 2.1.3 Die Lösung u aus Satz 2 ist schwache Lösung von (1.1) in t > 0.

Beweis: Sei a_{ij} = a_{ji} in Ω konstant und

Lu = ∑_{i,j=1}^{n} a_{ij} \frac{∂^{2}u}{∂x_{i}∂x_{j}} + c(x)u .

Dann ist

∫_{Ω} (Lu · φ − uLφ)dx = ∫_{∂Ω} (\frac{∂u}{∂ν}φ − u\frac{∂φ}{∂ν}) dσ

mit der (verallgemeinerten) Normalen ν_{i} = ∑_{j=1}^{n} a_{ij}n_{j}, wo n die äußere Normale

auf ∂Ω ist. Wir wenden dies an auf Lu = u_{tt} − u_{xx} − q(x)u und Ω = {(x, t) : t > |x|}. (x_{1} = x, x_{2} = t). Es ist auf ∂Ω = {(x, t) : t = |x|}

ν = −\frac{1}{\sqrt{2}} ( \begin{matrix}sgn(x) \\ 1 \end{matrix} ) .

Ist u die Lösung aus Satz 2, so ist Lu = 0 in Ω und u = 1/2 auf ∂Ω; außerhalb \overline{Ω} ist u = 0. Also ist für φ ∈ C_{0}^{∞}(t > 0)

∫_{0}^{+∞} ∫_{−∞}^{∞} uLφdxdt = ∫_{Ω} uLφdxdt = ∫_{Ω} Lu · φdx

−\frac{1}{\sqrt{2}} ∫_{−∞}^{+∞}((sgn(x)u_{x} + u_{t})φ − u(sgn(x)φ_{x} + φ_{t})) (x, |x|)dx

= − ∫_{∂Ω} (\frac{∂u}{∂τ} · φ − u\frac{∂φ}{∂τ}) dσ ,

wo τ einen Tangenteneinheitsvektor auf ∂Ω bedeutet. Wegen Lu = 0 in Ω verschwindet das Integral über Ω. Wegen u = 1/2 auf ∂Ω und φ ∈ C_{0}^{∞}(Ω) verschwindet auch das Integral über ∂Ω. u ist also in der Tat schwache Lösung von Lu = 0 in t > 0.

□

2.2 Das inverse Problem einer hyperbolischen Gleichung in einer räumlichen Variablen

Das inverse Problem zu (1.6) lautet nun: Gegeben sei g(t) = u(0, t), t ≥ 0. Man bestimme q.

In den Anwendungen wird das inverse Problem gewöhnlich für den Viertelraum 0 < x < ∞, 0 < t < ∞ formuliert:

u_{tt} = u_{xx} + q(x)u , x, t > 0
u_{x}(0, t) = 0 , t > 0 ,
u(x, 0) = 0 , u_{t}(x, 0) = δ(x) , x > 0 .

Man bestimme q aus g(t) = u(0, t), t > 0. Wir führen dieses Problem auf unseres zurück, indem wir q gerade auf IR^{1} fortsetzen.

Satz 2.2.1 (Nichtlineare Gelfand-Levitan-Gleichung): Sei

G(s, t) = \frac{1}{2} ∫_{0}^{s}(g(t + s′) + g(t − s′))ds′ .

Dabei sei g als ungerade Funktion auf IR^{1} fortgesetzt. Dann gilt

∫_{−s}^{s} u(x, s)u(x, t)dx = G(s, t) .

Wegen u(x, t) = 0 für x > t genügt es im Falle t ≤ s bereits, von −t bis t zu integrieren.

Beweis: Wir setzen Lu = u_{xx} + q(x)u. Für φ_{0}, φ_{1} ∈ C_{0}^{∞}(IR^{1}) können wir die Lösung von

u_{tt} = Lu in IR^{2} , u(x, 0) = φ_{0} , u_{t}(x, 0) = φ_{1} in IR^{1}

für −∞ < t < ∞ in der Form

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

schreiben. Dabei ist E_{0} gerade und E_{1} ungerade in t. Sei für φ, ψ ∈ C_{0}^{∞}(IR^{1})

w(s, t) = (E_{0}(s)ψ, E_{0}(t)φ)

mit dem inneren Produkt von L_{2}(IR^{1}), und sei v(s) = E_{0}(s)ψ, u(t) = E_{0}(t)φ. Dann ist

w_{ss}(s, t) = ∫_{IR^{1}} v_{ss}(x, s)u(x, t)dx

= ∫_{IR^{1}} (Lv)(x, s)u(x, t)dx

= ∫_{IR^{1}} v(x, s)(Lu)(x, t)dx

= ∫_{IR^{1}} v(x, s)u_{tt}(x, t)dx
= w_{tt}(s, t) .

Außerdem ist

w_{t}(s, 0) = (E_{0}(s)ψ , \frac{d}{dt}E_{0}(t)φ|_{t=0})) = 0 .

Nach der Kirchhoffschen Formel gilt also

w(s, t) = \frac{1}{2}(w(s − t, 0) + w(s + t, 0)) .

Da E_{0} gerade ist, ist auch w(s, 0) gerade, nach der letzten Formel also w(s, t) = w(t, s). Für s = 0 ergibt dies

(ψ, E_{0}(t)φ) = (E_{0}(t)ψ, φ) .

E_{0} ist also ein symmetrischer Operator, und das gleiche gilt für E_{1}.

Wir zeigen weiter, daß \frac{d}{dt}E_{1}(t)φ = E_{0}(t)φ für φ ∈ C_{0}^{∞}(IR^{1}), wofür wir kurz E′_{1} = E_{0} schreiben. Dazu setzen wir u = E_{1}(t)φ und v = u_{t}. Dann ist

v_{tt} = v_{xx} + q(x)v ,

\begin{aligned}v(x, 0) & = & u_{t}(x, 0) & = & φ(x) , \\ v_{t}(x, 0) & = & u_{tt}(x, 0) & = & u_{xx}(x, 0) + q(x)u(x, 0) = 0 , \end{aligned}

weil u = 0 für t = 0. Diese drei Beziehungen sagen aber gerade v = E_{0}φ. Als weitere Eigenschaft der E_{i} haben wir die Halbgruppeneigenschaft

u(t + s) = E_{0}(s)u(t) + E_{1}(s)u_{t}(t)

für jede Lösung u von u_{tt} = Lu. Dies ergibt sich einfach aus der Zeitinvarianz der Differentialgleichung. Seien nun φ, ψ ∈ C_{0}^{∞}(IR^{1}) und u(t) = E_{1}(t)φ. Dann gilt

(ψ, 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)φ) .

Weil E_{0} gerade und E_{1} ungerade ist, folgt

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

Addition der letzten beiden Beziehungen ergibt

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

also, wegen E_{1}(s) = ∫_{0}^{s} E_{0}(s′)ds′,

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

Mit φ = δ_{n}, ψ = δ_{m} und u_{n} = E_{1}δ_{n} lautet dies

\frac{1}{2} ∫_{0}^{s} ∫_{−∞}^{+∞} δ_{m}(x)(u_{n}(x, t + s′) + u_{n}(x, t − s′))dxds′ = ∫_{−∞}^{+∞}u_{m}(x, s)u_{n}(x, t)dx .

Lassen wir hier erst n → ∞, dann m → ∞ streben, so erhalten wir

\frac{1}{2} ∫_{0}^{s}(u(0, t + s′) + u(0, t − s′))ds′ = ∫_{−∞}^{+∞}u(x, s)u(x, t)dx

mit der Lösung u von (1.6). Wegen u(x, t) = 0 für |x| > t ist dies die behauptete Gleichung.

□

Satz 2.2.2 (Lineare Gelfand-Levitan-Gleichung): Es gibt eine in |t| ≤ x stetig differenzierbare Funktion w(x, t) mit q(x) = 4\frac{d}{dx}w(x, x) und

∫_{−x}^{x} g(t − s)w(x, s)ds = \frac{1}{2}(g(t − x) + g(t + x)) , |t| < x .

Dabei ist wie üblich g als ungerade Funktion auf IR^{1} fortgesetzt.

Beweis: Wir nehmen δ_{n} gerade an und setzen

g_{n}(t) = ∫_{IR^{1}} δ_{n}(t − s)g(s)ds .

Weiter definieren wir

u_{n}(x, t) = ∫_{IR^{1}} δ_{n}(t − s)u(x, s)ds .

Dabei ist u wie üblich als ungerade Funktion auf −∞ < t < ∞ fortgesetzt. Die Funktion u_{n} erfüllt die Differentialgleichung und die Anfangs- und Randbedingungen

u_{n}(x, 0) = ∫_{IR^{1}} δ_{n}(−s)u(x, s)ds = 0  (2.1)

u_{n}(0, t) = ∫_{IR^{1}} δ_{n}(t − s)g(s)ds = g_{n}(t)  (2.2)

\frac{∂u_{n}}{∂x} (0, t) = 0 .  (2.3)

Wir zeigen nun: Es gibt eine Funktion v_{n}, so daß

u_{n}(x, t) = ∫_{IR^{1}} g(t − s)v_{n}(x, s)ds .  (2.4)

Wir definieren v_{n} als Lösung von

\frac{∂^{2}v_{n}}{∂t^{2}} = \frac{∂^{2}v_{n}}{∂x^{2}} + q(x)v_{n} in IR^{2}

v_{n}(0, t) = δ_{n}(t) , \frac{∂v_{n}}{∂x} (0, t) = 0 , t ∈ IR^{1} .

v_{n} ist also die (als gerade Funktion in x fortgesetzte) Lösung der Cauchyschen Anfangswertaufgabe, wobei als Anfangsmannigfaltigkeit jetzt die t-Achse auftritt. Bei dieser Wahl von v_{n} erfüllt offenbar u_{n} die Bedingungen (2.1) - (2.3), welche u_{n} als Lösung der Differentialgleichung eindeutig festlegen.

Ein Blick auf die Charakteristiken zeigt, daß

u_{n}(x, t) = 0 , |t| < |x| − \frac{1}{n} ,
v_{n}(x, t) = 0 , |t| > |x| + \frac{1}{n} .

Also folgt aus (2.4) für x > 0

∫_{−x−1/n}^{x+1/n} g(t − s)v_{n}(x, s)ds = 0 , |t| < x − 1/n .

Für v_{n} ergibt sich aus der Kirchhoffschen Formel

v_{n}(x, t) = \frac{1}{2}(δ_{n}(t + x) + δ_{n}(t − x)) − w_{n}(x, t)
w_{n}(x, t) = +\frac{1}{2} ∫_{Δ} ∫_{(x,t)} q(x′)v_{n}(x′, t′)dx′dt′ ,

wobei nun Δ(x, t) das Dreieck ist, das von den Charakteristiken durch (x, t) und der t-Achse gebildet wird. Damit bekommen wir schließlich

∫_{−x−1/n}^{x+1/n} g(t−s)w_{n}(x, s)ds = \frac{1}{2} ∫_{−x−1/n}^{x+1/n} g(t−s)(δ_{n}(s+x)+δ_{n}(s−x))ds , |t| < −1/n .

Für n → ∞ ergibt sich daraus die behauptete Integralgleichung. Für die rechte Seite ist dies klar. Für die linke schreiben wir w_{n} in der Form

w_{n}(x, t) = r_{n}(x, t) − \frac{1}{2} ∫_{Δ} ∫_{(x,t)} q(x′)w_{n}(x′, t′)dx′dt′ ,

r_{n}(x, t) = \frac{1}{4} ∫_{Δ} ∫_{(x,t)} q(x′)(δ_{n}(t′ + x′) + δ_{n}(t′ − x′))dx′dt′ .

Die δ_{n}-Funktionen verschwinden außerhalb von Streifen der Breite 2/n um die Charakteristiken. Diese treffen den Rand von Δ(x, t) für |t| ≤ x bei (x + t)/2 bzw. (x − t)/2. Also gilt

lim_{n→∞} r_{n}(x, t) = \frac{1}{4} ( ∫_{0}^{(x−t)/2}q(x′)dx′ + ∫_{0}^{(x+t)/2}q(x′)dx′) = r(x, t)

gleichmäßig in |t| ≤ x. Es folgt, daß w_{n} für |t| ≠ x gegen eine Funktion w konvergiert, und zwar gleichmäßig auf kompakten Teilmengen. Diese Funktion erfüllt

w(x, t) = r(x, t) − \frac{1}{2} ∫_{Δ} ∫_{(x,t)} q(x′)w(x′, t′)dx′dt′ .

Wegen w(x, t) = 0 für |t| > x folgt w(x, x) = r(x, x), also \frac{d}{dx}w(x, x) = 4q(x).

□

Die Integralgleichung aus Satz 2 ist von erster Art. Man kann eine Gleichung zweiter Art erhalten, indem man nach t differenziert und dabei beachtet, daß g(+0) = \frac{1}{2}, g(−0) = −\frac{1}{2} ist. Bezeichnen wir mit g′ die gerade Fortsetzung von g′ auf IR^{1}, so entsteht

w(x, x) + ∫_{−x}^{x} g′(t − s)w(x, s)ds = \frac{1}{2}(g′(t − x) + g′(t + x)) , |t| ≤ x .

2.3 Ein Beispiel aus der Reflektions - Seismologie

Sei t die Zeit und z die Tiefe in der Erde, die wir uns als isotropes elastisches Medium mit Dichte ρ(z) und Lame’-Parametern λ(z), μ(z) vorstellen. Sei w(z, t) die z-Komponente der Teilchengeschwindigkeit und p(z, t) der Druck. Dann gilt

ρw_{t} + p_{z} = 0 , p_{t} + (λ + 2μ)w_{z} = 0 .  (3.1)

Wir nehmen λ, μ und ρ als stetig differenzierbar an.

Wir wollen das System (3.1) in eine Gleichung zweiter Ordnung verwandeln. Dazu differenzieren wir die zweite Gleichung nach t, die erste nach z und eliminieren w_{tz}. Es entsteht

p_{tt} = c^{2} (p_{zz} − \frac{ρ_{z}}{ρ} p_{z}) , c = (\frac{λ + 2μ}{ρ})^{1/2} .

c ist also die Ausbreitungsgeschwindigkeit von Wellen in der Erde. Wir wollen diese auf 1 normieren und führen dazu die neue unabhängige Tiefen-Variable

x = ∫_{0}^{z} \frac{dz′}{c(z′)} , also p_{z} = \frac{1}{c}p_{x} , p_{zz} = −\frac{c_{z}}{c^{2}} p_{x} + \frac{1}{c^{2}} p_{xx}

ein. Es entsteht

p_{tt} = p_{xx} − (c_{z} + \frac{ρ_{z}}{ρ} c) p_{x} .

Um hier eine rechte Seite der Form u_{xx} + qu zu erhalten, setzen wir p = Au.

Es entsteht

u_{tt} = u_{xx} + (2\frac{A_{x}}{A} − c_{z} − \frac{ρ_{z}}{ρ} c) u_{x} + (\frac{A_{xx}}{A} − (c_{z} + \frac{ρ_{z}}{ρ} c) \frac{A_{x}}{A} ) u .

Wir wählen nun A so, daß der Koeffizient von u_{x} verschwindet, also

\frac{A_{x}}{A} = \frac{1}{2} (c_{z} + \frac{ρ_{z}}{ρ} c)

oder

\frac{A_{z}}{A} = \frac{1}{2} (\frac{c_{z}}{c} + \frac{ρ_{z}}{ρ} )

also

A = (ρc)^{1/2} .

Dann lautet die Differentialgleichung

u_{tt} = u_{xx} + qu ,
q = \frac{A_{xx}}{A} − (c_{z} + \frac{ρ_{z}}{ρ} c) \frac{A_{x}}{A} = (\frac{A_{x}}{A} )_{x} − 2c (\frac{A_{x}}{A} )^{2} .

Aus q kann man also höchstens die Funktion A bestimmen, nicht ρ, λ, μ separat.

Wir nehmen nun an, daß man für den Impuls

u(x, 0) = 0 , u_{t}(x, 0) = δ(x)

die Funktion g(t) = u(0, t) messen kann. Dann ist q nach den Methoden von §2 zu bestimmen.

2.4 Näherungslösung für kleine q

Ist q klein, so kann man leicht eine Näherungslösung für das inverse Problem angeben. Die Integralgleichung für u können wir mit Hilfe der Heaviside-Funktion

H(t) = \{ \begin{matrix}1 & , & t > 0 , \\ 0 & , & sonst \end{matrix}

als

u(x, t) = \frac{1}{2}H(t − |x|) + \frac{1}{2} ∫ ∫_{D(x,t)} q(x′)u(x′, t′)dx′dt′  (4.1)

schreiben. Es folgt

u(x, t) = \frac{1}{2}H(t − |x|) + O(q) .

Dies setzen wir nun auf der rechten Seite von (4.1) ein. Es folgt

u(x, t) = \frac{1}{2}H(t − |x|) + \frac{1}{2} ∫ ∫_{D(x,t)} q(x′)H(t′ − |x′|)dx′dt′ + O(q^{2}) .

Für x = 0 bekommen wir also bis auf O(q^{2})

g(t) = \frac{1}{2} + \frac{1}{4} ∫_{−t}^{t} q(x′) ∫_{0}^{t−|x′|}H(t′ − |x′|)dt′dx′

= \frac{1}{2} + \frac{1}{4} ∫_{−t/2}^{t/2} q(x′)(t − 2|x′|)dx′ .

Durch Differenzieren bekommt man näherungsweise

g′(t) = \frac{1}{4} ∫_{−t/2}^{t/2} q(x′)dx′ ,

g′′(t) = \frac{1}{8} (q(\frac{t}{2}) + q(−\frac{t}{2})) = \frac{1}{4}q(\frac{t}{2}) ,

falls q gerade ist.

Für kleine q bekommt man q also aus g durch zweimalige Differentiation. Differentiation ist ein instabiler Prozeß. Zum ersten Mal sehen wir an dieser Stelle, daß die Lösung inverser Probleme zu Instabilitäten neigt.

2.5 Diskrete inverse Evolutionsprobleme

Sei T Tridiagonalmatrix

T = ( \begin{matrix}a_{0} & b_{0} \\ c_{0} & ... & ... \\ & ... & & & b_{n−1} \\ & & & c_{n−1} & a_{n} \end{matrix} ) , a_{i}, b_{i}, c_{i} ∈ IR

und seien u^{ℓ} ∈ IR^{n+1} eine Lösung von

u^{ℓ+1} = T u^{ℓ} .  (5.1)

Wir stellen das inverse Problem: Gegeben u^{0} und g^{ℓ} = u_{0}^{ℓ}, ℓ = 0, 1, 2, . . .. Man bestimme T!

Wir werden in diesem Abschnitt verschiedene Aufgaben dieser Art mit Hilfe eines diskreten Analogons der nichtlinearen Gelfand-Levitan-Gleichung lösen.

1. Fokussierter Anfangszustand u^{0} = e^{0} = (1, 0, …, 0)^{t}.

Wir führen die Vektoren

v^{ℓ+1} = T^{t}v^{ℓ} , ℓ = 0, 1, …, v^{0} = e^{0}  (5.2)

ein und haben dann

g^{ℓ+k} = e^{0t}u^{ℓ+k} = e^{0t}T^{ℓ+k}e^{0} = (T^{ℓt}e^{0})^{t}(T^{k}e^{0})
= v^{ℓt}u^{k} .

Mit den (n + 1, n + 1)-Matrizen

G = ( g^{0} g^{1} … g^{n} ; ⋮ ; g^{n} g^{n+1} … g^{2n} ) , U = (u^{0}, …, u^{n}) , V = (v^{0}, …, v^{n})

lautet dies

G = V^{t}U .  (5.3)

U, V sind obere Dreiecksmatrizen, weil T tridiagonal ist.

Wir betrachten nun verschiedene Spezialfälle, je nachdem welche Information über T vorhanden ist.

(a) c_{i} = b_{i} > 0, i = 0, …, n − 1, d.h. T ist symmetrisch. Dann ist U = V, und die Diagonale ist positiv. Also ist U durch (5.2) eindeutig bestimmt und durch eine Cholesky-Zerlegung berechenbar.

Ist U bestimmt, so berechnen wir T aus Gleichung ℓ, ℓ + 1 von (5.1):

\begin{aligned}u_{ℓ}^{ℓ+1} & = & c_{ℓ−1}u_{ℓ−1}^{ℓ} + a_{ℓ}u_{ℓ}^{ℓ} \\ & & & , & ℓ = 0, …, n − 1 , \\ u_{ℓ+1}^{ℓ+1} & = & c_{ℓ}u_{ℓ}^{ℓ} \end{aligned}  (5.4)

wo c_{−1} = u_{−1}^{ℓ} = 0 gesetzt wurde. Da u_{ℓ}^{ℓ} ≠ 0 ist, sind dadurch c_{ℓ} = b_{ℓ} und a_{ℓ} für ℓ = 0, …, n − 1 bestimmt. T ist also bestimmt bis auf das (n, n)-Element a_{n}.

(b) c_{i}, i = 0, …, n − 1 bekannt, c_{i}, b_{i} ≠ 0. U und V sind nach wie vor rechte Dreiecksmatrizen, und die Diagonale (1, c_{0}, c_{0}c_{1}, …, c_{0} ⋯ c_{n−1}) von U ist bekannt und enthält keine Null. (5.3) ist jetzt eine LR-Zerlegung von G. Diese ist eindeutig bestimmt, weil die Diagonale von U bekannt ist und keine Nullen enthält. Ist V bestimmt, so berechnen wir T aus den Zeilen ℓ und ℓ + 1 von (5.2), also

\begin{aligned}v_{ℓ}^{ℓ+1} & = & b_{ℓ−1}v_{ℓ−1}^{ℓ} + a_{ℓ}v_{ℓ}^{ℓ} \\ & & & , & ℓ = 0, …, n − 1 \\ v_{ℓ+1}^{ℓ+1} & = & b_{ℓ}v_{ℓ}^{ℓ} . \end{aligned}  (5.5)

Wegen b_{ℓ} ≠ 0 sind alle v_{ℓ}^{ℓ} = b_{0} ⋯ b_{ℓ−1} ≠ 0, also sind b_{0}, …, b_{n−1} durch die zweite Gleichung eindeutig bestimmt. Danach bestimmt man die a_{ℓ} bis auf a_{n} aus der ersten Gleichung.

(c) c_{i}, b_{i} ≠ 0, i = 0, …, n − 1. Jetzt sind die Diagonalen von U und V nach wie vor frei von Nullen, sind aber beide unbekannt. Sei D die Diagonale von V. (5.3) lautet

G = (D^{−1}V)^{t}DU .

DV^{t} ist eine linke Dreiecksmatrix mit Diagonale 1. Also ist die rechte Dreiecksmatrix W = DU eindeutig als rechter Faktor einer LR-Zerlegung von G bestimmt. Sei W = (w^{0}, …, w^{n}). Dann lautet (5.1)

w^{ℓ+1} = DT D^{−1}w^{ℓ} , ℓ = 0, …, n − 1 .  (5.6)

Da die w^{ℓ} bekannt und die Außerdiagonalelemente der Tridiagonalmatrix DT D^{−1} ungleich 0 sind, kann man hieraus DT D^{−1} wie oben bis auf das (n, n)-Element berechnen. Damit sind a_{0}, …, a_{n−1} bestimmt. Über b_{i}, c_{i} kann man nichts sagen.

2. Nicht-fokussierter Anfangszustand

Sei u^{0} = e und e ein Vektor mit von Null verschiedenen Komponenten. (5.3) gilt nach wie vor, aber U ist keine rechte Dreiecksmatrix mehr. Also machen wir eine LR-Zerlegung von U, U = LR (falls es diese gibt). L wählen wir so, daß die Diagonale gerade die von V^{−1} ist (möglich, falls die b_{j} ≠ 0). Dann ist

G = (V^{t}L)R

ein LR-Zerlegung von G, deren linker Faktor die Diagonale 1 hat. R ist dadurch eindeutig festgelegt und durch eine LR-Zerlegung von G zu berechnen. Aus (5.1) folgt für die Spalten r^{ℓ} von R

r^{ℓ+1} = Sr^{ℓ} , S = L^{−1}T L .

Aus der ersten Beziehung folgt, daß S eine rechte Hessenberg-Matrix ist (d.h. s_{ℓj} = 0, ℓ > j + 1), aus der zweiten, daß S auch eine linke Hessenberg-Matrix ist (d.h. s_{ℓj} = 0, j > ℓ + 1). Damit ist S wie T eine Tridiagonalmatrix. Wie oben sieht man, daß S durch R bis auf das (n, n)-Element bestimmt ist.

Wir nehmen nun an, die Außerdiagonalelemente b_{i}, c_{i}, i = 0, …, n − 1 seien bekannt, ebenso der Vektor e. Seien ℓ_{j} die Zeilen von L. Wegen e = Lr^{0} sind dann die Elemente in der Position 0 der ℓ_{j} bekannt. Insbesondere ist ℓ_{0} bekannt. LS = T L lautet durch die ℓ_{j} ausgedrückt (ℓ_{−1} = ℓ_{n+1} = 0)

ℓ_{j}S = c_{j−1}ℓ_{j−1} + a_{j}ℓ_{j} + b_{j}ℓ_{j+1} ,  j = 0, …, n .  (5.7)

Hieraus lassen sich die a_{j} rekursiv bestimmen: Für j = 0 ist

ℓ_{0}S = a_{0}ℓ_{0} + b_{0}ℓ_{1} .

ℓ_{0} ist bekannt, also (für n > 0) auch ℓ_{0}S. Die Komponente 1 von ℓ_{1} kann also bestimmt werden. a_{0} ergibt sich dann aus der Komponente 0 dieser Gleichung. Somit sind a_{0}, ℓ_{0}, ℓ_{1} bekannt. Nehmen wir an, es seien bereits a_{0}, …, a_{j−1}, ℓ_{0}, …, , ℓ_{j} für ein j < n bekannt. Dann bestimmen wir aus (5.7) zunächst einmal die Komponente j + 1 von ℓ_{j+1}. Dies ist möglich, weil ℓ_{j}S (trotz des fehlenden (n, n)-Elements von S) berechnet werden kann und b_{j} bekannt und ≠ 0 ist. Dann bestimmen wir aus der Komponenten 0 dieser Gleichung a_{j}. Dies ist möglich, weil die c_{j} bekannt und die Komponente 0 der ℓ_{j0} bekannt und ≠ 0 sind. Damit sind a_{0}, …, a_{j}, ℓ_{0}, …, ℓ_{j+1} bekannt. Also lassen sich a_{0}, …, a_{n−1} bestimmen.

Wenden wir uns nun Evolutionsgleichungen zweiter Ordnung zu. Im Hinblick auf die Anwendungen beschränken wir uns auf Gleichungen der Form

u^{ℓ+1} + u^{ℓ−1} = T u^{ℓ} ,  ℓ = 0, 1, ….  (5.8)

Dabei seien nun die Anfangszustände u^{0}, u^{−1} = 0 gegeben. Wie oben wollen wir T aus den Zahlen

g^{ℓ} = u_{0}^{ℓ} , ℓ = 0, 1, …

bestimmen.

Lemma 2.5.1 Seien die (n + 1, n + 1)-Matrizen T_{ℓ} bestimmt durch die Rekursion

T_{ℓ+1} + T_{ℓ−1} = T T_{ℓ} , ℓ = 0, 1, …  (5.9)
T_{0} = I , T_{−1} = 0 .

Dann gilt für die Vektoren aus (5.8)

u^{ℓ+k} = T_{k}u^{ℓ} − T_{k−1}u^{ℓ−1} , k, ℓ = 0, 1, ….

Insbesondere ist also u^{k} = T_{k}u^{0}.

Beweis: Wir führen den Beweis durch Induktion nach k. Für k = 0, 1 ist die Behauptung richtig. Sei die Behauptung richtig bis zu einem k ≥ 1. Dann ist nach (5.8)

u^{ℓ+k+1} + u^{ℓ+k−1} = T u^{ℓ+k}

und nach Induktionsannahme

u^{ℓ+k+1} + (T_{k−1}u^{ℓ} − T_{k−2}u^{ℓ−1}) = T (T_{k}u^{ℓ} − T_{k−1}u^{ℓ−1})

oder

u^{ℓ+k+1} = (T T_{k} − T_{k−1})u^{ℓ} − (T T_{k−1} − T_{k−2})u^{ℓ−1} .

Wegen (5.9) ist dies die Behauptung für k + 1.

□

Wir betrachten nun den fokussierten Anfangszustand

u^{0} = e_{0} = (1, 0, …, 0)^{t} , u^{−1} = 0 .

Dann haben wir

g^{k+ℓ} = e^{0t}u^{ℓ+k} = e^{0t}(T_{k}u^{ℓ} − T_{k−1}u^{ℓ−1})
= (T_{k}^{t}e^{0})^{t}u^{ℓ} − (T_{k−1}^{t}e^{0})^{t}u^{ℓ−1}
= v^{kt}u^{ℓ} − v^{(k−1)t}u^{ℓ−1} ,  (5.10)

wo wir

v^{k} = T_{k}^{t}e^{0}

gesetzt haben. Die v^{k} sind Lösung von

v^{k+1} + v^{k−1} = T^{t}v^{k} , k = 0, 1, …
v^{0} = e^{0} , v^{−1} = 0 .

Dies folgt aus dem Lemma mit T^{t} an Stelle von T. Aus (5.10) folgt durch Addition für ν ≤ Min (k, ℓ)

g^{k+ℓ} + g^{k+ℓ−2} + ⋯ + g^{k+ℓ−2ν} = v^{kt}u^{ℓ} − v^{(k−ν−1)t}u^{(ℓ−ν−1)} .

Für ℓ ≥ k wählen wir ν = k und erhalten wegen v^{−1} = 0

g^{k+ℓ} + g^{k+ℓ−2} + ⋯ + g^{ℓ−k} = v^{kt}u^{ℓ} .

Entsprechend erhalten wir für ℓ ≤ k und ν = ℓ wegen u^{−1} = 0

g^{k+ℓ} + g^{k+ℓ−2} + ⋯ + g^{k−ℓ} = v^{kt}u^{ℓ} .

Definieren wir also die (n + 1, n + 1)-Matrix G = (g_{ℓk}) durch

g_{ℓk} = g^{k+ℓ} + g^{k+ℓ−2} + ⋯ + g^{|k−ℓ|} , ℓ, k = 0, …, n und U = (u^{0}, …, u^{n}), V = (v^{0}, …, v^{n}), so haben wir

G = V^{t}U .  (5.11)

V, U sind wieder rechte Dreiecksmatrizen. Wir können also wie im Falle von Evolutionsgleichungen 1. Ordnung verfahren.

(a) Seien b_{i} = c_{i} > 0, i = 0, …, n − 1. Dann ist T = T^{t} und damit U = V mit positiver Diagonale. U kann also durch Cholesky-Zerlegung von G bestimmt werden. Die Beziehungen (5.4) gelten nach wie vor und bestimmen die b_{i}, c_{i}, a_{i} mit Ausnahmen von a_{n}.

(b) c_{i}, b_{i} ≠ 0 und c_{i} bekannt, i = 0, …, n − 1. Dann ist die Diagonale von U bekannt und enthält keine Nullen, so daß U, V aus G durch LR-Zerlegung eindeutig berechnet werden können. Es gilt nach wie vor (5.5), woraus die b_{i}, a_{i} mit Ausnahme von a_{n} berechnet werden können.

(c) b_{i}, c_{i} ≠ 0, aber sonst unbekannt, i = 0, …, n − 1. Wieder argumentiert man wie im Falle der Evolution erster Ordnung. Die Beziehung (5.6) lautet jetzt

w^{ℓ+1} + w^{ℓ−1} = DT D^{−1}w^{ℓ} , ℓ = 0, …, n − 1 ,

erfüllt aber den gleichen Zweck und bestimmt DT D^{−1} bis auf das (n, n)- Element. Damit sind die a_{i} bis auf a_{n} bestimmt.

Näherungen für das inverse Problem der Wellengleichung

Sei Ω ein Gebiet in IR^{n}, und sei c = c(x) die lokale Schallgeschwindigkeit eines Mediums in Ω. Sei s ∈ ∂Ω eine Quelle und r ∈ ∂Ω ein receiver. Zum Zeitpunkt t = 0 werde an der Stelle s ein Impuls mit dem Zeitverlauf f(t) gegeben. An der Stelle r werde die Funktion g(t) gemessen. Zu bestimmen ist c aus den Messungen von g(t) für viele Paare r, s. Eine wichtige Anwendung ist die Seismologie. Hier ist Ω der Halbraum (Erde) {(x_{1}, x_{2}, x_{3}) : x_{3} > 0}, wo x_{3} die Tiefe bedeutet. Die Menge der gemessenen Funktionen g heißt Seismogramm. Für die Wahl von r, s gibt es verschiedene Möglichkeiten:

1. Zero-offset. Es ist r = s, und beide laufen gemeinsam über die Erdoberfläche x_{3} = 0. Meistens sind sie auf eine Gerade, z.B. x_{2} = 0, beschränkt.

2. Fixed source. Hier ist s fest, z.B. (0, 0, 0), und r läuft über die Oberfläche oder nur über eine Gerade, z.B. x_{2} = x_{3} = 0.

Wir werden in diesem Kapitel eine Reihe von Näherungslösungen zunehmender Raffinesse kennenlernen. In der Seismologie spricht man von Migration. Am Ende geben wir eine einheitliche Darstellung dieser Techniken durch die Born-Approximation.

3.1 Klassische Migration (Hagedorn 1954)

Die hierzu zu besprechenden Migrationstechniken machen keinen Gebrauch von der Wellengleichung. Sie beruhen auf einfachen Laufzeitüberlegungen und können nur auf sehr einfache Erdmodelle angewendet werden. Betrachten wir etwa ein zweidimensionales Erdmodell, bei dem c oberhalb einer Kurve Γ den konstanten Wert c_{0}, unterhalb den Wert c_{1} ≠ c_{0} hat. Wir stellen uns vor, daß Γ aus Reflektoren unterschiedlicher Stärke aufgebaut ist. Sei g(x_{1}, t) das zero-offset Seismogramm, also die Zeitreihe des Seismometers an dem gemeinsamen Ort x_{1} von Quelle und receiver. Ist Γ die Gerade x_{3} = z ∈ IR, so peaked g(x_{1}, t) für t = 2z/c_{0}, und der Funktionswert im peak ist ein Maß für die Stärke des Reflektors. Bei allgemeinen Kurven peaked g(x_{1}, t) entlang einer gewissen, Γ ähnlichen Kurve Γ′.

Aufgabe der Migration ist es, aus Γ′ die ursprüngliche Kurve Γ wiederzugewinnen. Wir lösen dieses Problem durch drei Techniken. In jedem Fall konvertieren wir im Seismogramm die Zeitangaben in Tiefenangaben, so daß also die Zeit t der Tiefe \frac{c_{0}}{2} t entspricht. Ebensogut können wir c_{0} = 2 annehmen.

1. Geometrische Migration. Sei a ein Reflektor in der Erde. Die von ihm erzeugten peaks im Seismogramm sind auf der Kurve R_{a} = {(x_{1}, t) : t = \frac{2}{c_{0}} \sqrt{(x_{1} − a_{1})^{2} + a_{3}^{2}}\} (vgl. Abb. 3.1) zu finden. Dies ist eine Hyperbel, deren Scheitel direkt über a liegt. Γ′ ist nun gerade die Einhüllende all dieser Kurven, wenn a Γ durchläuft. Dies sieht man wie folgt ein. Um die Einhüllende der Kurvenschar R_{a} : a ∈ Γ zu berechnen, stellt man Γ als Kurve a = a(λ) mit dem Parameter λ dar. Die Einhüllende bekommt man dann dadurch, daß man die Gleichung von R_{a} nach λ differenziert, nach λ auflöst und in die Gleichung für R_{a} einsetzt. Differentiation der Gleichung nach λ bedeutet aber, |r − a| für a ∈ Γ zu minimieren. Für das minimierende a ergibt sich aus der Gleichung für R_{a} aber gerade ein Punkt auf Γ′!

Um nun Γ aus Γ′ zu gewinnen, berechnet man für jeden Punkt von Γ′ die berührende Hyperbel R_{a} und addiert den Funktionswert von g in den Scheitel von R_{a}. Dies ergibt das migrierte Seismogramm.

2. Wave front migration. Für jeden Punkt (x_{1}, t) im Seismogramm sei I_{x_{1},t} die Kurve \{a = (a_{1}, a_{2}) : t = \frac{2}{c_{0}} \sqrt{(x_{1} − a_{1})^{2} + a_{3}^{2}}\} (vgl. Abb. 3.1) der Abbildung 3.1: Darstellung von Seismogramm und Schallgeschwindigkeit im Boden. Die Unstetigkeit von c(x) im Boden zeigt sich als Kurve im Seismogramm.

Reflektoren, die zu einem Signal bei (x_{1}, t) beitragen können. Man ordnet den Punkten dieser Kurve den Wert g(x_{1}, t) zu. Um Γ aus Γ′ zu gewinnen, mittelt man für jeden Punkt a über alle I_{x_{1},t}, welche durch a gehen.

3. Maximal convexity migration. Für jeden Punkt a in der Erde betrachtet man die Kurve R_{a} im Seismogramm. Jeder Punkt (x_{1}, t) dieser Kurve entspricht einem receiver, der von einem Reflektor in a ein Signal bekommt. Man mittelt den Wert g(x_{1}, t) über alle diese Punkt und weist diese Mittel dem migrierten Seismogramm im Punkt a zu.

Das migrierte Seismogramm ist also im letzten Fall

∫_{(x_{1},t)∈R_{a}} g(x_{1}, t) ,

im zweiten Fall

∫_{a∈I_{(x_{1},t)}} g(x_{1}, t) .

Hier wird offenbar über die gleichen Mengen integriert. Die beiden Methoden liefern also das gleiche Resultat.

Die letzten beiden Methoden können weitgehend verallgemeinert werden. Sei τ(x, y) die Laufzeit eines Signals von x nach y. Sei d = (r, s, t) ein Punkt im Seismogramm und g(d) der zugehörige Funktionswert. Sei

\begin{aligned}R_{a} & = & \{d = (r, s, t) : τ(s, a) + τ(a, r) = t\} , \\ I_{d} & = & \{a : τ(s, a) + τ(a, r) = t\} . \end{aligned}

R_{a} heißt Fläche maximaler Konvexität (maximal convexity curve, reflectiontime surface), I_{d} die isochrone Fläche (maximal convexity surface). Die migrierten Seismogramme sind dann

∫_{d∈R_{a}} g = ∫_{I_{d}∋a} g .

3.2 Migration durch die Wellengleichung

Wir betrachten wieder den zweidimensionalen Fall (x_{2} = 0) mit zero-offset-Daten, welche von einer reflektierenden Kurve Γ herrühren, oberhalb derer die Schallgeschwindigkeit c konstant ist. Im exploding reflector model nimmt man an, die Quellen seien entlang Γ verteilt, und man halbiert alle Laufzeiten. Dann hat man näherungsweise folgendes Modell: Es ist

u_{tt} = c^{2}Δu c = c_{0} ∈ IR_{+}  (2.1)

zwischen Erdoberfläche und Γ. An der Erdoberfläche x_{3} = 0 ist

u(x_{1}, 0, t) = g(x_{1}, t)  (2.2)

gegeben, und für t = 0 hat man

u(x_{1}, x_{3}, 0) = δ_{Γ}(x_{1}, x_{3}) .  (2.3)

Dabei ist δ_{Γ} eine Funktion, deren Träger in etwa mit Γ übereinstimmt. Man versucht nun, δ_{Γ} (und damit Γ) aus g zu bestimmen. Selbst durch diese groben Vereinfachungen ist noch kein vernünftiges inverses Problem entstanden. Wir werden aber trotzdem so eine Art Lösung angeben. Dazu bemerken wir zunächst, daß für jedes f ∈ C^{2} die Funktion

u(x_{1}, x_{3}, t) = f(ωt + ξ_{1}x_{1} + ξ_{3}x_{3})  (2.4)

die Wellengleichung erfüllt, wenn nur die Dispersionsrelation

\frac{ω^{2}}{c^{2}} = ξ_{1}^{2} + ξ_{3}^{2}  (2.5)

erfüllt ist. Diese hat die beiden Lösungen ω = ±c\sqrt{ξ_{1}^{2} + ξ_{3}^{2}}. u ist entlang ωt + ξ_{1}x_{1} + ξ_{3}x_{3} = const. konstant. Man nennt u daher eine ebene Welle und die Ebenen ξ_{1}x_{1} + ξ_{2}x_{2} = const. ihre Wellenfronten. Haben ω und ξ_{3} das gleiche Vorzeichen, so läuft die Welle mit wachsendem t in Richtung fallender x_{3}-Werte, also nach oben, andernfalls nach unten. Da im exploding reflector model die Quellen in der Erde sitzen, sind wir nur an aufsteigenden Wellen interessiert, wählen also

ω = cξ_{3}\sqrt{1 + ξ_{1}^{2}/ξ_{3}^{2}} .  (2.6)

Spezielle ebene Wellen sind

e^{i(ωt+ξ_{1}x_{1}+ξ_{3}x_{3})} .

Linearkombinationen von Lösungen der Wellengleichung sind wieder Lösungen. Also ist auch

u(x_{1}, x_{3}, t) = ∫ ∫ Â(ξ_{1}, ξ_{3})e^{i(ωt+ξ_{1}x_{1}+ξ_{3}x_{3})}dξ_{1}dξ_{3}

eine Lösung. Wir wollen Â so bestimmen, daß u die Lösung unseres Problems wird. Dazu muß zunächst einmal die Bedingung an der Oberfläche, also

g(x_{1}, t) = ∫ ∫ Â(ξ_{1}, ξ_{3})e^{i(ωt+ξ_{1}x_{1})}dξ_{1}dξ_{3}  (2.7)

erfüllt sein. Wir werden von jetzt an von der Fourier-Transformation Gebrauch machen (siehe §3). Danach ist

g(x_{1}, t) = \frac{1}{2π} ∫ ∫ e^{i(ξ_{1}x_{1}+ωt)}ĝ(ξ_{1}, ω)dξ_{1}dω ,  (2.8)
ĝ(ξ_{1}, ω) = \frac{1}{2π} ∫ ∫ e^{−i(ξ_{1}x_{1}+ωt)}g(x_{1}, t)dx_{1}dt .

Mit ω aus (2.6) ergibt sich

\frac{dω}{dξ_{3}} = c\frac{|ξ_{3}|}{\sqrt{ξ_{1}^{2} + ξ_{3}^{2}}} = \frac{c}{\sqrt{1 + ξ_{1}^{2}/ξ_{3}^{2}}}

und damit aus (2.8) nach der Transformationsformel

g(x_{1}, t) = \frac{c}{2π} ∫ ∫ e^{i(ξ_{1}x_{1}+ωt)}ĝ(ξ_{1}, ω)\frac{dξ_{1}dξ_{3}}{\sqrt{1 + ξ_{1}^{2}/ξ_{3}^{2}}} .  (2.9)

Vergleich von (2.7), (2.9) führt zu

Â(ξ_{1}, ξ_{3}) = \frac{c}{\sqrt{1 + ξ_{1}^{2}/ξ_{3}^{2}}} ĝ (ξ_{1}, cξ_{3}\sqrt{1 + ξ_{1}^{2}/ξ_{3}^{2}}) .

Aus (2.3) bekommt man jetzt

δ_{Γ}(x_{1}, x_{3}) = ∫ ∫ Â(ξ_{1}, ξ_{3})e^{i(ξ_{1}x_{1}+ξ_{3}x_{3})}dξ_{1}dξ_{3}
= 2πA(x_{1}, x_{3}) .

Da Γ im wesentlichen der Träger von δ_{Γ} ist, ist Γ damit näherungsweise bestimmt.

Wir wollen nachweisen, daß die Methoden für lineares Γ, also x · θ = s, θ = (θ_{1}, 0, θ_{3}), θ_{1}^{2} + θ_{3}^{2} = 1 das exakte Resultat liefert. Die Daten sind in diesem Fall g(x_{1}, t) = δ(t+x_{1}θ_{1} −s) für c = 1. Für θ_{3} > 0 sind dies die Werte der aufsteigenden Welle δ(t + x · θ − s) an der Oberfläche. Damit wird mit Hilfe von III.3.4

ĝ(ξ_{1}, ω) = \frac{1}{2π} ∫ ∫ e^{−i(ξ_{1}x_{1}+ωt)}δ(t + x_{1}θ_{1} − s)dx_{1}dt
= \frac{1}{2π} ∫ e^{−i(ξ_{1}x_{1}−ω(x_{1}θ_{1}−s))}dx_{1}
= e^{−iωs}δ(ξ_{1} − ωθ_{1}) .

Also ist

Â(ξ_{1}, ξ_{3}) = \frac{1}{\sqrt{1 + ξ_{1}^{2}/ξ_{3}^{2}}} ĝ (ξ_{1}, ξ_{3}\sqrt{1 + ξ_{1}^{2}/ξ_{3}^{2}})

= \frac{1}{√}e^{−isξ_{3}√} δ(ξ_{1} − ξ_{3}θ_{1}√ ) .

Inverse Fouriertransformation ergibt

A(x_{1}, x_{3}) = \frac{1}{2π} ∫ ∫ \frac{1}{√} e^{−isξ_{3}√} δ(ξ_{1} − ξ_{3}θ_{1}√ )e^{i(x_{1}ξ_{1}+x_{3}ξ_{3})}dξ_{1}dξ_{3} .

An Stelle von ξ_{1} führen wir die neue Variable

y = ξ_{1} − ξ_{3}θ_{1}√

ein. Man berechnet

\frac{dy}{dξ_{1}} = 1 − \frac{θ_{1}ξ_{1}/ξ_{3}}{√} .

Für y = 0 ist für θ_{3} > 0

√ = \frac{1}{θ_{3}} , ξ_{1} = \frac{θ_{1}}{θ_{3}} ξ_{3} , \frac{dy}{dξ_{1}} = θ_{3}^{2} .

Damit erhält man mittels III.3.4

A(x_{1}, x_{3}) = \frac{1}{2π} ∫ ∫ \frac{1}{√} e^{−isξ_{3}√} δ(y)e^{i(x_{1}ξ_{1}+x_{3}ξ_{3})} \frac{dydξ_{3}}{|\frac{dy}{dξ_{1}} |}

= \frac{1}{2πθ_{3}} ∫ e^{i(−sξ_{3}/θ_{3}+x_{1}ξ_{3}θ_{1}/θ_{3}+x_{3}ξ_{3})}dξ_{3}
= \frac{1}{θ_{3}} δ(\frac{1}{θ_{3}} (−s + x_{1}θ_{1} + x_{3}θ_{3}))
=  δ(x · θ − s) .

3.3 Fourier-Transformation und

Distributionen

Für f ∈ L_{1}(IR^{n}) definieren wir die Fourier-Transformation

f̂(ξ) = (2π)^{−n/2} ∫ e^{−ix·ξ}f(x)dx .

Offenbar ist f̂ stetig und beschränkt,

∥f̂∥_{∞} ≤ (2π)^{−n/2}∥f∥_{1} , ∥f∥_{∞} = lim_{p→∞} (∫ |f|^{p}dx)^{1/p} .

Die inverse Fourier-Transformation ist f̃(ξ) = (2π)^{−n/2} ∫ e^{ix·ξ}f(x)dx .

Satz 3.3.1 Sei f ∈ L_{1}(IR^{n}) ∩ C^{1}(IR^{n}) und g = f̂. Dann existiert g̃ im Sinne von

g̃(x) = lim_{k→∞}(2π)^{−n/2} ∫_{|ξ_{j}|≤k} e^{iξ·x}g(ξ)dξ ,

und es gilt g̃ = f.

Beweis: Es ist

g̃(x) = lim_{k→∞}(2π)^{−n/2} ∫_{|ξ_{j}|<k} e^{iξ·x}(2π)^{−n/2} ∫ e^{−iξ·y}f(y)dydξ
= lim_{k→∞} ∫ f(y)(2π)^{−n} ∫_{|ξ_{j}|≤k} e^{iξ·(x−y)}dξdy
= lim_{k→∞} ∫ f(y)δ_{k}(x − y)dy ,
δ_{k}(x) = (2π)^{−n} ∫_{|ξ_{j}|≤k} e^{iξ·x}dξ .  (3.1)

Für n = 1 ist

δ_{k}(x) = \frac{k}{π} sinc(kx) , sinc(x) = \frac{sin x}{x} .

Man berechnet leicht

∫ δ_{k}(x)dx = 1 .

Damit wird für f ∈ C^{1}(IR^{1}) ∩ L_{1}(IR^{1})

∫ f(x)δ_{k}(x)dx = ∫ (f(x) − f(0))δ_{k}(x)dx + f(0)
= \frac{1}{π} ∫ \frac{f(x) − f(0)}{x} sin kx dx + f(0) .

Mit k → ∞ geht das Integral gegen 0. Also gilt

lim_{k→∞} ∫ f(x)δ_{k}(x)dx = f(0) .

Für n > 1 ist δ_{k} das Produkt der eindimensionalen δ_{k}-Funktion. Es folgt g̃(x) = f(x).

□

Die Funktion δ_{k} aus (3.1) ist also eine approximative δ- Funktion. Ist

χ_{k}(ξ) = (2π)^{−n/2} \{ \begin{matrix}1 & , & |ξ_{j}| ≤ k \\ 0 & , & sonst \end{matrix} ,

so gilt χ̃_{k} = δ_{k} und δ̂_{k} = χ_{k}. Wir schreiben deshalb

((2π)^{−n/2})^{∼} = δ , δ̂ = (2π)^{−n/2} .

Als Beispiel berechnen wir das (im klassischen Sinn nicht existierende) Integral

u(x, t) = −\frac{1}{2(2π)^{3}} ∫_{IR^{3}} e^{ix·ξ}(e^{iωt} − e^{−iωt})\frac{dξ}{|ξ|} ,  (3.2)

wo ω die Dispersionsrelation ω = c|ξ| erfüllt. Wir führen für ξ Polarkoordinaten ρ, φ, ψ um x ein, so daß also |ξ| = ρ, x · ξ = rρ cos ψ mit r = |x| wird. Es wird dann

u(x, t) = −\frac{1}{2(2π)^{3}} ∫_{0}^{∞}∫_{0}^{2π}∫_{0}^{π} ρ sin ψe^{irρ cos ψ}(e^{icρt} − e^{−icρt})dψdφdρ

= −\frac{1}{2(2π)^{2}} ∫_{0}^{∞}(e^{icρt} − e^{−icρt}) ∫_{0}^{π} ρ sin ψ e^{irρ cos ψ}dψdρ

= \frac{1}{2(2π)^{2}r} ∫_{0}^{∞}(e^{icρt} − e^{−icρt})(e^{−irρ} − e^{irρ})dρ

= \frac{1}{2(2π)^{2}r} ∫_{−∞}^{+∞}(e^{−icρt} − e^{icρt})e^{irρ}dρ

= \frac{1}{4πr}(δ(r − ct) − δ(r + ct)) .  (3.3)

Damit haben wir Kugelwellen als Linearkombination von ebenen Wellen dargestellt.

Sei ϕ ∈ C^{1}(Ω) und sei Γ = {x ∈ Ω : ϕ(x) = 0}. Sei |∇ϕ| ≠ 0 entlang Ω. Wir zeigen

∫_{Ω} f(x)δ(ϕ(x))dx = ∫_{Γ} f \frac{dσ}{|∇ϕ|} .  (3.4)

Dazu nehmen wir an, daß ϕ(x) = z nach einem x_{i}, etwa x_{1} auflösbar ist. Dies ist zumindest lokal möglich. Dann ist mit x′ = (x_{2}, …, x_{n})

∫ f(x)δ(ϕ(x))dx = ∫ ∫ f(x_{1}, x′)δ(z)dx_{1}dx′
= ∫ ∫ f(ψ(z, x′), x′)δ(z) \frac{dzdx′}{|ϕ_{x_{1}}(x_{1}, x′)|}

= ∫ f(ψ(0, x′), x′) \frac{dx′}{|ϕ_{x_{1}}(x_{1}, x′)|} .

x_{1} = ψ(0, x′) ist eine Darstellung von Γ. Das Oberflächenelement auf Γ ist

dσ(x) = (1 + ∑_{i=2}^{n} (\frac{∂ψ}{∂x_{i}} )^{2})^{1/2} dx′ .

Wegen ϕ(ψ(0, x′), x′) = 0 entlang Γ ist

\frac{∂ϕ}{∂x_{1}} \frac{∂ψ}{∂x_{i}} + \frac{∂ϕ}{∂x_{i}} = 0 , i = 2, …, n .

Also folgt

dσ(x) = \frac{1}{|\frac{∂ϕ}{∂x_{1}} |} |∇ϕ|dx′

und damit (3.4). Durch eine ähnliche Rechnung bestätigt man

δ(rx) = r^{−n}δ(x)

für die n-dimensionale δ-Funktion und r > 0. Die Fourier-Transformation hat folgende Eigenschaften:

1. Sei f_{r}(x) = f(rx), r > 0. Dann ist f̂_{r}(ξ) = r^{−n}f̂(ξ/r).

2. Sei f_{y}(x) = f(x + y), y ∈ IR^{n}. Dann ist f̂_{y}(ξ) = e^{iξ·y}f̂(ξ).

3. Ist f, \frac{∂f}{∂x_{j}} ∈ L_{1}(IR^{n}), so ist (\frac{∂f}{∂x_{j}})^{∧}(ξ) = iξ_{j}f̂(ξ).

4. Sei f ⋆g(x) = ∫ f(x−y)g(y)dy, f, g ∈ L_{1}(IR^{n}). Dann ist f ⋆g ∈ L_{1}(IR^{n}),

und

(f ⋆ g)^{∧} = (2π)^{n/2}f̂ĝ .

Sei ϕ ∈ L_{1}(IR^{n}). ϕ definiert einen Filter gemäß

f̂ → ϕ̂f̂ , f ∈ L_{1}(IR^{n}) .

Nach Eigenschaft 4 ist hierzu äquivalent

f → (2π)^{−n/2}ϕ ⋆ f .

Einer Faltung im Orts- (oder Zeit-) raum entspricht also einer Multiplikation im Fourier-Raum. Ein Filter ϕ mit der Eigenschaft ϕ̂(ξ) = 0, |ξ| > Ω heißt Tiefpassfilter, Ω seine Bandbreite. Wichtigstes Beispiel ist der ideale Tiefpass

ϕ̂(ξ) = \{ \begin{matrix}1 & , & |ξ| ≤ Ω \\ 0 & , & sonst . \end{matrix}

Für diesen ist (vgl. Aufgabe 16)

ϕ(x) = (2π)^{−n/2}Ω^{n} \frac{J_{n/2}(Ω|x|)}{(Ω|x|)^{n/2}} = δ_{Ω}(x) .

ϕ⋆f ist eine geglättete Version von f, bei der Details ≤ π/Ω geglättet wurden. Ein Filter ψ mit ψ̂(ξ) = 0 für |ξ| ≤ Ω heißt Hochpass. Der ideale Hochpass

ψ̂(ξ) = \{ \begin{matrix}1 & , & |ξ| ≥ Ω \\ 0 & , & sonst \end{matrix}

ist zwar nicht in L_{1}. Man kann seine Wirkung aber beschreiben durch

f → f − δ_{Ω} ⋆ f = ψ ⋆ f .

Er hat folgende Eigenschaft: Für eine beliebig offene Menge U ⊆ IR^{n} ist f ∈ C^{∞}(U) genau dann, wenn ψ ⋆ f ∈ C^{∞}(U). Sucht man also nur die Stellen von f, wo f nicht C^{∞} ist (der singuläre Träger s sup(f)), so kann man genausogut ψ ⋆ f wie f studieren.

3.4 Asymptotische Lösung der Wellengleichung

Wieder betrachten wir die Wellengleichung

\frac{∂^{2}u}{∂t^{2}} = c^{2}(x)Δu

mit einer Funktion c(x) > 0, und wir suchen Lösungen der Form

u(x, t) = e^{−iωt}v(x) ,
Δv + \frac{ω^{2}}{c^{2}} v = 0 .  (4.1)

Wir nehmen nun an, daß ω groß ist, und suchen eine Lösung von (4.1) in der Form

v(x) = e^{iωϕ(x)} ∑_{j=0}^{∞} \frac{v_{j}(x)}{(iω)^{j}} .  (4.2)

Setzen wir dies in (4.1) ein, so entsteht

(iωΔϕ − ω^{2}|∇ϕ|^{2}) ∑_{j=0}^{∞} \frac{v_{j}}{(iω)^{j}} + 2iω∇ϕ · ∑_{j=0}^{∞} \frac{∇v_{j}}{(iω)^{j}} + ∑_{j=0}^{∞} \frac{Δv_{j}}{(iω)^{j}}

+\frac{ω^{2}}{c^{2}} ∑_{j=0}^{∞} \frac{v_{j}}{(iω)^{j}} = 0 .

Diese Beziehung soll für alle ω gelten. Vergleichen wir die Faktoren der Potenzen ω^{2}, ω, 1, ω^{−1}, …, so entsteht der Reihe nach

\begin{aligned}ω^{2} & : & |∇ϕ|^{2} = \frac{1}{c^{2}} \\ ω & : & 2∇ϕ · ∇v_{0} + v_{0}Δϕ = 0 \\ 1 & : & 2∇ϕ · ∇v_{1} + v_{1}Δϕ + Δv_{0} = 0 \\ ω^{−1} & : & 2∇ϕ · ∇v_{2} + v_{2}Δϕ + Δv_{1} = 0 \end{aligned}

usw. Die erste Gleichung ist die uns wohlbekannte Eikonal-Gleichung. Wir schreiben sie in der Form

F (x, p, ϕ) = |p| − \frac{1}{c(x)} = 0 .  (4.3)

Die weiteren Gleichungen sind sogenannte Transportgleichungen. Diese Gleichungen werden nun sukzessive nach ϕ, v_{0}, …, v_{j} aufgelöst und dadurch eine Funktion v bestimmt, welche (4.1) bis auf Terme der Ordnung ω^{−j} erfüllt. Z.B. wird (4.1) für j = 0 bis auf den Term Δv_{0} erfüllt.

Betrachten wir zunächst die Funktion

v(x) = e^{iωϕ(x)} ,

wobei ϕ die Eikonal-Gleichung erfüllt. Sie führt zu der Näherung

u(x, t) = e^{iω(ϕ(x)−t)}

für die Wellengleichung. Sie erfüllt (4.1) bis auf O(ω). u beschreibt eine Welle, deren Phase in allen Punkten der Fläche ϕ(x) = t die gleiche ist. Diese Fläche heißt daher Wellenfront. Für konstantes c > 0 ist z.B. ϕ(x) = p · x mit |p| = \frac{1}{c} eine Lösung der Eikonal-Gleichung, die Wellenfront also die Ebene p · x = t. Diese bewegt sich mit der Geschwindigkeit c. Ein anderes Beispiel ist ϕ(x) = \frac{1}{c}|x|. Die Wellenfront hat jetzt die Gleichung |x| = ct und stellt ebenfalls eine sich mit der Geschwindigkeit c fortbewegende Fläche dar. Dies stimmt auch für nicht konstantes c. Ist nämlich x(t) ein Punkt einer Wellenfront, der sich senkrecht zu ihr fortbewegt, so gilt

ϕ(x(t)) = t ,  ẋ(t) = α(t)∇ϕ(x(t))

mit einem Skalar α(t). Differenzieren der ersten Beziehung nach t ergibt zusammen mit der zweiten

α(t)|∇(ϕ(x(t))|^{2} = 1

und damit wegen der Eikonal-Gleichung

α(t) = c^{2}(x(t)) , |ẋ(t)| = c^{2}(x(t))|∇ϕ(x(t))| = c(x(t)) .

Also bewegt sich unser Punkt x(t) mit der Geschwindigkeit c(x(t)) und damit die ganze Wellenfront in jedem ihrer Punkte x mit der Geschwindigkeit c(x) senkrecht zu sich selbst.

Das charakteristische System der Eikonal-Gleichung ist nach PDE 1

ẋ = \frac{p}{|p|} ,  ṗ = ∇(\frac{1}{c(x)}) ,  ϕ̇ = \frac{1}{c(x)} .  (4.4)

Wir wollen zeigen, daß die Lösungen x(t) von (4.4), also die Kurven, welche charakteristische Streifen tragen, gerade die Kurven sind, entlang denen sich Signale in dem Geschwindigkeitsfeld c(x) ausbreiten. Nach dem Fermat’schen Prinzip erfolgt diese Ausbreitung nämlich so, daß die Laufzeit zwischen zwei Punkten minimal ist. Mit anderen Worten: Ist K eine beliebige Kurve, welche die Punkte x_{0}, x_{1} verbindet, so läuft ein Signal, das bei x_{0} startet und bei x_{1} beobachtet wird entlang der Kurve K^{∗}, für welche

∫_{K} \frac{ds}{c(x)} ,  s = Bogenlänge  (4.5)

möglichst klein ist.

Wir brauchen nun folgendes Hilfsmittel aus der Variationsrechnung.

Satz 3.4.1 Sei L ∈ C^{1}(IR^{2n}), und sei x_{0}, x_{1} ∈ IR^{n}. Unter allen Kurven x, welche x_{0} mit x_{1} verbinden, sei x^{∗} so, daß

∫_{0}^{1} L(x^{∗}, ẋ^{∗})dt ≤ ∫_{0}^{1} L(x, ẋ)dt .

Ist x^{∗} ∈ C^{2}[0, 1], so gelten die Euler’schen Gleichungen

\frac{d}{dt} \frac{∂L}{∂ẋ_{i}} (x^{∗}, ẋ^{∗}) − \frac{∂L}{∂x_{i}} (x^{∗}, ẋ^{∗}) = 0 , i = 1, …, n .

Wenden wir dies nun an auf das Problem, (4.5) zu minimieren. Hier haben wir

L(x, ẋ) = \frac{|ẋ|}{c(x)} .

Die Euler’schen Gleichungen lauten

\frac{d}{dt} \frac{ẋ}{c(x)|ẋ|} − |ẋ|∇\frac{1}{c(x)} = 0

Diese Gleichungen sind invariant gegenüber Parametertransformationen. Also können wir z.B. die Bogenlänge als Parameter verwenden. Dann ist |ẋ| = 1, und wir erhalten

\frac{d}{dt} \frac{ẋ}{c(x)} = ∇\frac{1}{c(x)} .

Mit p = ẋ/c(x), also |p| = c(x), wird hieraus gerade (4.4).

Damit steht fest: Die (Träger der) Charakteristiken der Eikonal-Gleichung genügen dem Fermatschen Prinzip (und werden daher als “Strahlen” angesehen).

3.5 Die Bornsche Näherung

Wir nehmen an, die gesuchte Schallgeschwindigkeit c unterscheide sich nur in einem Gebiet Ω von der Hintergrundschallgeschwindigkeit c_{0}:

\frac{1}{c^{2}} = \frac{1}{c_{0}^{2}} (1 + f) ,  supp(f) ⊆ Ω .  (5.1)

Die Quellen s und die receiver r liegen außerhalb von Ω. Zur Zeit t = 0 werde an der Stelle s ein Impuls gegeben, dessen Zeitverlauf durch das wavelet w(t) mit w(t) = 0, t < 0 beschrieben wird. Typische wavelets sind

w(t) = δ(t) , w(t) = e^{iωt}χ(t) , χ ∈ C_{0}^{∞}(0, ∞) .

Wir beschreiben das Experiment durch

\frac{∂^{2}u}{∂t^{2}} = c^{2}(Δu − δ(x − s)w(t)) in IR^{3} × IR^{1} ,

(5.2)

u = 0 für t < 0 .

Sei u^{0} die Lösung für c = c_{0}. Wir setzen u = u^{0} + v und finden für v

\frac{∂^{2}v}{∂t^{2}} = c_{0}^{2}Δv − f \frac{∂^{2}u}{∂t^{2}} in IR^{3} × IR^{1} ,
v = 0 für t < 0 .

Die Born-Approximation besteht nun darin, daß man hier u durch u^{0} ersetzt:

\frac{∂^{2}v}{∂t^{2}} = c_{0}^{2}Δv − f \frac{∂^{2}u^{0}}{∂t^{2}} in IR^{3} × IR^{1}

(5.3)

v = 0 für t < 0 .

Wir führen die weiteren Rechnungen zunächst für den Fall einer konstanten Hintergrundschallgeschwindigkeit c_{0} durch.

Wir benutzen die Kirchhoffsche Lösung

v(x, t) = \frac{1}{4πc_{0}^{2}} ∫_{|x−x′|≤c_{0}t} \frac{h(x′, t − |x − x′|/c_{0})}{|x − x′|}dx′  (5.4)

der Anfangswertaufgabe

\frac{∂^{2}v}{∂t^{2}} = c_{0}^{2}Δv + h(x, t) , v = v_{t} = 0 bei t = 0 .

Mit h = −fu_{tt}^{0} erhalten wir

v(x, t) = −\frac{1}{4πc_{0}^{2}} ∫_{|x−x′|≤c_{0}t} f(x′) \frac{u_{tt}^{0}(x′, t − |x − x′|/c_{0})}{|x − x′|} dx′ .

Da u^{0}(x, t) = 0 für t < 0 können wir die Einschränkung des Integrationsbereiches weglassen und dann die zweite Ableitung bei u^{0} vor das Integral ziehen. Da f außerhalb Ω verschwindet, erhalten wir

v(x, t) = −\frac{1}{4πc_{0}^{2}} \frac{∂^{2}}{∂t^{2}} ∫_{Ω} f(x′) \frac{u^{0}(x′, t − |x − x′|/c_{0})}{|x − x′|} dx′ .  (5.5)

Ebenfalls aus der Kirchoffschen Formel (5.4) erhalten wir

u^{0}(x, t) = −\frac{1}{4π} ∫_{|x−x′|≤c_{0}t} δ(x′ − s)w(t − |x − x′|/c_{0}) \frac{dx′}{|x − x′|}

= −\frac{1}{4π} \frac{w(t − |x − s|/c_{0})}{|x − s|} .

Setzen wir dies in (5.5) ein, so entsteht schließlich

v(x, t) = \frac{1}{16π^{2}c_{0}^{2}} \frac{∂^{2}}{∂t^{2}} ∫_{Ω} f(x′) \frac{w(t − (|x − x′| + |x′ − s|)/c_{0})}{|x − x′||x′ − s|} dx′ . (5.6)

Dies ist die Lösung des direkten Problems in der Bornschen Näherung. Sie wird zur Lösung des inversen Problems in folgender Weise benutzt. Sei g(r, s, t) die Daten-Funktion, also g(r, s, t) = v(r, t), mit der Lösung u = u^{0} + v von (5.3). Dann lautet (5.6)

g(r, s, t) = \frac{1}{16π^{2}c_{0}^{2}} \frac{∂^{2}}{∂t^{2}} ∫_{Ω} f(x)\frac{w(t − (|r − x| + |x − s|)/c_{0})}{|r − x||x − s|} dx .

Dies ist eine lineare Integralgleichung für f. Für w = δ lautet sie (vgl. (3.4))

g(r, s, t) = \frac{1}{16π^{2}c_{0}} \frac{∂^{2}}{∂t^{2}} ∫_{I_{d}} f(x) \frac{dσ(x)}{||r − x|(s − x) + |s − x|(r − x)|}  (5.7)

mit der Isochronen I_{d} = {x ∈ Ω : c_{0}t = |r − x| + |x − s|}, d = (r, s, t) aus §1.

Die Born-Approximation reduziert also inverse Streuprobleme auf die Bestimmung von f aus Integralen über Mannigfaltigkeiten niedriger Dimension. Dies ist das Problem der Integralgeometrie, das wir in Kap. IV behandeln werden.

An Stelle der wavelets w kann man auch die zeitharmonische Einstrahlung mit w(t) = e^{−iωt} bei fester Frequenz ω verwenden. Dann kann man u(x, t) = e^{−iωt}u(x) annehmen und hat dann an Stelle von (5.3)

Δu + \frac{ω^{2}}{c^{2}} u = δ(x − s) in IR^{3}

oder, mit k = ω/c_{0}

Δu + k^{2}u = −k^{2}fu + δ(x − s) in IR^{3} .  (5.8)

Aus PDE I wissen wir, daß

Δu + k^{2}u = h

zusammen mit der Ausstrahlungsbedingung

u = O (\frac{1}{r}) , \frac{∂u}{∂r} − iku = O (\frac{1}{r^{2}} )

mit r = |x| für h ∈ C_{0}^{1}(IR^{3}) eindeutig lösbar ist, und zwar ist

u(x) = ∫ γ(x − y)h(y)dy

mit der Fundamentallösung

γ(x) = \frac{e^{ik|x|}}{4π|x|} .  (5.9)

Die zugehörige Lösung

e^{−iωt}γ(x) = \frac{e^{iω(|x|/c_{0}−t)}}{4π|x|}

entspricht dann einer Welle, deren Wellenfronten |x| = tc_{0} + const. sich mit wachsender Zeit ausbreitenden Kugeloberflächen entsprechen. Das eben besagt das Wort “Ausstrahlungsbedingung”.

Wir suchen wieder Lösungen in der Form u^{0} +v, wo u^{0} die Lösung mit c = c_{0} darstellt und alle Funktionen die Ausstrahlungsbedingung erfüllen. Offenbar ist

u^{0}(x) = γ(x − s) .

Aus (5.9) folgt dann

v(x) = −k^{2} ∫ γ(x − x′)f(x′)u(x′)dx′ .

Wieder erhalten wir die Born-Approximation, indem wir hier u durch u^{0} ersetzen:

v(x) = −k^{2} ∫_{Ω} γ(x − x′)f(x′)γ(x′ − s)dx′

Mit der Daten-Funktion g(r, s, ω) = v(r) erhalten wir

g(r, s, ω) = −\frac{k^{2}}{16π^{2}} ∫_{Ω} f(x) \frac{e^{ik(|r−x|+|x−s|)}}{|r − x||x − s|} dx .  (5.10)

Hier machen wir eine weitere Approximation. Wir nehmen an, daß r und s auf einer großen Kugel vom Radius R liegen und lassen R → ∞ streben. Mit r = Rα, s = Rβ, α, β ∈ S^{2} ist dann

\begin{aligned}|r − x| & = & R − α · x & + & O(\frac{1}{R}) \\ |s − x| & = & R − β · x & + & O(\frac{1}{R}) \end{aligned}

und damit bis auf O (\frac{1}{R^{3}} ), d.h. bis auf einen relativen Fehler von O (\frac{1}{R})

g(r, s, ω) = −\frac{k^{2}e^{ikR}}{16π^{2}R^{2}} ∫_{Ω} f(x)e^{−ik(α+β)·x}dx

= −\frac{k^{2}e^{ikR}}{16π^{2}R^{2}} (2π)^{3/2}f̂(k(α + β)) .

Abbildung 3.2: Ewald-Sphäre für eine feste receiver-Position α. Der gestrichelte Anteil entspricht dem backscatter, der durchgezogene der Transmissionsmessung.

Für eine feste Quelle s liefert g also die Werte der Fourier-Transformation auf einer Sphäre um kβ mit Radius k. Man nennt diese Sphären Ewald-Sphären. Entsprechendes gilt für eine feste receiver-Position.

Mißt man nur den backscatter (d.h. α · β ≥ 0), so kann man f̂(ξ) nur für |ξ| ≥ k\sqrt{2} bestimmen. Dies bedeutet, daß man nur eine hochpassgefilterte Version von f bestimmen kann. Dies reicht aus, um z.B. Unstetigkeiten von f bestimmen zu können. Macht man hingegen nur Transmissionsmessungen (α · β ≤ 0), so kann man f̂(ξ) nur für |ξ| ≤ \sqrt{2}k bestimmen. Dadurch wird die räumliche Auflösung beschränkt.

Zur Gültigkeit der Born-Approximation wollen wir noch eine heuristische Bedingung angeben. Dazu beschränken wir uns auf eine konstante Störung f in [−ρ, ρ]. Die ungestörte Helmholtz-Gleichung lautet dann

u′′_{0} + k^{2}u_{0} = 0 ,

die gestörte

u′′ + k^{2}(1 + f)u = 0 .

Die Lösungen sind

u_{0}(x) = e^{±ikx} , u(x) = e^{±ik\sqrt{1+f}x} .

Die Born-Approximation verlangt u ∼ u_{0}. Dazu ist (hochgradig) notwendig, daß die Phasen sich um deutlich weniger als π unterscheiden, d.h.

|k\sqrt{1 + f} − k|2ρ ≪ π .

Für kleines f bedeutet dies

k|f|ρ ≪ π .

Ist D der Durchmesser eines Objekts mit konstantem f, so erhalten wir

D|f| ≪ \frac{2π}{k} .

3.6 Die “distorted” Born-Approximation

Ist der Hintegrund nicht konstant, so kann man nicht mehr die Fundamentallösung (5.9) benutzen. Statt dessen nehmen wir die in §4 abgeleitete Näherung

u(x) = A(x)e^{iωτ}  (6.1)
|∇τ| = \frac{1}{c_{0}(x)} , 2∇τ ∇A + AΔτ = 0 .  (6.2)

Sei s eine Quelle außerhalb supp (f) und Γ_{s} eine kleine Kugelsphäre um s. Dann legen wir A, τ durch die Anfangsbedingungen

τ(x) = \frac{|x − s|}{c_{0}(s)} ,  A(x) = \frac{1}{4π|x − s|}  (6.3)

fest. Die Anfangswertaufgaben (6.2), (6.3) sind - zumindest in einer Umgebung von Γ_{s} - eindeutig lösbar. Dazu erinnern wir an die Theorie der quasilinearen Differentialgleichung F (x, u, p) = 0 mit p = ∇u in IR^{n}. Die Anfangswertaufgabe lautet: Sei eine Streifenmannigfaltigkeit

S : x = φ(t) ,  u = χ(t) , p = p(t) ,  t ∈ IR^{n−1}

\frac{dχ}{dt_{j}} = p · \frac{∂φ}{∂t_{j}} j = ℓ, …, n − 1 ,

F (φ(t), χ(t), p(t)) = 0

gegeben. Gesucht ist eine Lösung u der Differentialgleichung, welche S enthält, also

\begin{aligned}u(φ(t)) & = & χ(t) , \\ ∇u(φ(t)) & = & p(t) . \end{aligned}

Nach II.5 aus PDE I ist u eindeutig bestimmt, wenn S nicht charakteristisch ist, d.h.

Rg (∇_{p}F , \frac{∂φ}{∂t_{1}} , …, \frac{∂φ}{∂t_{n−1}} ) = n .

Wir wenden dies an auf |∇τ| = 1/c_{0} mit der Anfangsbedingung τ = ε/c_{0}(y) für |x − y| = ε. Wir ergänzen die Anfangsmannigfaltigkeit zu einer Streifenmannigfaltigkeit

φ(t) = s + εθ(t) , χ(t) = \frac{ε}{c_{0}(s)} , p = p(t) ,

wobei θ(t) eine Parametrisierung von S^{n−1} ist und p die Streifenbedingung

\frac{∂χ}{∂t_{j}} = 0 = p · \frac{∂φ}{∂t_{j}} , j = 1, …, n

und |p(t)| = 1/c_{0}(φ(t)) erfüllt. Die Rangbedingung lautet

Rg (\frac{p}{|p|} , \frac{∂φ}{∂t_{1}} , …, \frac{∂φ}{∂t_{n−1}} ) = n

und ist offenbar erfüllt.

Für die Differentialgleichung 2∇τ ·∇A+AΔτ = 0 lautet die Rangbedingung

Rg (2∇τ , \frac{∂φ}{∂t_{1}} , …, \frac{∂φ}{∂t_{n−1}} ) = n .

Wegen ∇τ = p ist auch diese Bedingung erfüllt.

Die Lösungen von (6.2), (6.3) bezeichnen wir auch mit τ(x, s), A(x, s). (5.10) nimmt dann die Form

g(r, s, ω) = −ω^{2} ∫_{Ω} \frac{f(x)}{c_{0}(x)}A(r, x)A(x, s)e^{iω(τ(r,x)+τ(x,s))}dx  (6.4)

an. Dies ist eine Näherung für großes ω.

Um eine entsprechende Approximation im Zeitbereich zu bekommen, also für das Problem (5.1), (5.3) mit w(t) = δ(t) nehmen wir an, daß (6.4) auch für kleine ω näherungsweise gilt. Fourier-Transformation von (6.4) bezüglich ω ergibt dann

ĝ(r, s, t) = ∫_{Ω} f(x)A(r, x)A(x, s)(2π)^{−1/2} ∫ ω^{2}e^{−iω(t−τ(r,x)−τ(x,s))}dωdx

= (2π)^{1/2} ∫_{Ω} f(x)A(r, x)A(x, s)δ′′(t − τ(r, x) − τ(x, s))dx

= \frac{∂^{2}}{∂t^{2}} ∫ f(x)A(r, x, s)δ(t − τ(r, x, s))dx ,
A(r, x, s) = (2π)^{1/2}A(r, x)A(x, s) , τ(r, x, s) = τ(r, x) + τ(x, s) .

Nach (3.4) erhalten wir schließlich als Analogon zu (5.7)

g(r, s, t) = \frac{∂^{2}}{∂t^{2}} ∫_{I_{d}} f(x)A(r, x, s)\frac{dσ(x)}{c_{0}(x)|∇_{x}τ(r, x, s)|}  (6.5)

mit der Isochronen I_{d} = {x : t = τ(r, x, s)} und d = (r, s, t). Die Größen A, τ sind durch das bekannte c_{0} bestimmt.

In Kap. IV werden wir integralgeometrische Aufgaben dieser Art durch Inversionsformeln lösen.

Kapitel 4 Integral-Geometrie

4.1 Die Radon-Transformation

Prototyp der integralgeometrischen Transformationen ist die Radon-Transformation

Rf(θ, s) = ∫_{x·θ=s} f(x)dσ(x) , θ ∈ S^{n−1} , s ∈ IR^{1}.

Dies hat Sinn z.B. für f ∈ C_{0}(IR^{1}). Andere Schreibweisen sind

(Rf)(θ, s) = ∫_{y∈θ^{⊥}} f(sθ + y)dy
= ∫_{IR^{n}} f(x)δ(s − x · θ)dx

Satz 4.1.1 (Projektionssatz): Sei f ∈ C_{0}(IR^{n}). Dann ist

(Rf(θ, ·))^{∧}(σ) = (2π)^{(n−1)/2}f̂(σθ) .

Dabei steht links die Fourier-Transformation in IR^{1}, rechts in IR^{n}. Beweis: Es ist

((Rf)(θ, ·))^{∧}(σ) = (2π)^{−1/2} ∫_{IR^{1}} e^{−isσ}(Rf)(θ, s)ds

= (2π)^{−1/2} ∫_{IR^{1}} e^{−isσ} ∫_{θ^{⊥}} f(sθ + y)dyds .

Mit x = sθ + y lautet dies

((Rf)(θ, ·))^{∧}(σ) = (2π)^{−1/2} ∫_{IR^{n}} e^{−iσx·θ}f(x)dx
= (2π)^{(n−1)/2}f̂(σθ) .
□

Für eine Funktion g ∈ C(S^{n−1} × IR^{1}) definieren wir

(R^{∗}g)(x) = ∫_{S^{n−1}} g(θ, x · θ)dσ(θ) .

Betrachten wir R als Operator von L_{2}(IR^{n}) in L_{2}(S^{n−1}×IR^{1}), so ist R^{∗} gerade der adjungierte zu R:

(Rf, g) = ∫_{S^{n−1}} ∫_{IR^{1}} (Rf)(θ, s)\overline{g}(θ, s)dsdσ(θ)

= ∫_{S^{n−1}} ∫_{IR^{1}} ∫_{θ^{⊥}} f(sθ + y)dy\overline{g}(θ, s)dσ(θ)ds

= ∫ f(x) ∫_{S^{n−1}} \overline{g}(θ, x · θ)dσ(θ)dx
= (f, R^{∗}g) .

Satz 4.1.2 Für f ∈ C_{0}(IR^{n}) gilt

R^{∗}Rf = |S^{n−2}| |x|^{−1} ⋆ f .

Beweis: Es ist

(R^{∗}Rf)(x) = ∫_{S^{n−1}} (Rf)(θ, x · θ)dσ(θ)

= ∫_{S^{n−1}} ∫_{θ^{⊥}} f((x · θ)θ + y)dydσ(θ)

= ∫_{S^{n−1}} ∫_{θ^{⊥}} f(x + y)dydσ(θ) .

Nun benutzen wir die Formel

∫_{S^{n−1}} ∫_{θ^{⊥}} f(y)dydσ(θ) = |S^{n−2}| ∫_{IR^{n}} |y|^{−1}f(y)dy ,

welche etwa für f ∈ C_{0}(IR^{n}) richtig ist. Damit erhalten wir

(R^{∗}Rf)(x) = |S^{n−2}| ∫_{IR^{n}} |y|^{−1}f(x + y)dy

= |S^{n−2}| ∫_{IR^{n}} |y − x|^{−1}f(y)dy .
□

Satz 4.1.3 (Radonsche Inversionsformel): Sei f ∈ C_{0}(IR^{n}) und g = Rf. Dann ist

f = R^{∗}Kg ,

(Kg(θ, ·))^{∧}(σ) = c|σ|^{n−1}ĝ(θ, σ) , c = \frac{1}{2}(2π)^{1−n} .

Beweis: Es ist

f(x) = (2π)^{−n/2} ∫_{IR^{n}} e^{+ix·ξ}f̂(ξ)dξ

= (2π)^{−n/2} ∫_{0}^{∞}|σ|^{n−1} ∫_{S^{n−1}} e^{iσx·θ}f̂(σθ)dσ(θ)dσ

= (2π)^{−n/2}(2π)^{(1−n)/2} ∫_{0}^{∞}|σ|^{n−1} ∫_{S^{n−1}} e^{iσx·θ}ĝ(θ, σ)dσ(θ)dσ .

Ersetzen wir θ durch −θ und σ durch −σ, so ändert sich das Integral nicht. Also können wir das so veränderte Integral addieren und die Summe durch 2 dividieren:

f(x) = \frac{1}{2}(2π)^{−n/2}(2π)^{(1−n)/2} ∫_{−∞}^{+∞}|σ|^{n−1} ∫_{S^{n−1}} e^{iσx·θ}ĝ(θ, σ)dσ(θ)dσ

= \frac{1}{2}(2π)^{−n/2}(2π)^{(1−n)/2} \frac{1}{c} ∫_{S^{n−1}} ∫_{−∞}^{+∞}(Kg)^{∧}(θ, σ)e^{iσx·θ}dσdσ(x)

= \frac{1}{2}(2π)^{−n/2}(2π)^{(1−n)/2}(2π)^{1/2} \frac{1}{c} ∫_{S^{n−1}} (Kg)(θ, x · θ)dσ(x)
= (R^{∗}Kg)(x) .
□

Wir betrachten zunächst den Fall n = 3. In diesem Fall ist

(Kg)^{∧}(σ) = cσ^{2}ĝ(σ) = −c(g′′)^{∧}(σ) ,

also Kg = −cg′′. Damit nimmt die Radonsche Inversionsformel die Gestalt

f(x) = −\frac{1}{8π^{2}} ∫_{S^{2}} g′′(θ, x · θ)dσ(θ)

an. Wir können sie auch in der Form

f(x) = −\frac{1}{8π^{2}} Δ ∫_{S^{2}} g(θ, x · θ)dσ(θ)

mit dem Laplaceschen Operator Δ schreiben.

Der Fall n = 2 ist nicht so einfach. Es ist jetzt

(Kg)^{∧}(σ) = \frac{1}{4π} |σ|ĝ(σ)
= \frac{1}{4π} i sgn(σ)(g′)^{∧}(σ) .

In Aufgabe 14 haben wir gesehen, daß die Fourier-Transformation von \frac{1}{s} im Sinne des Hauptwertes existiert und −i(\frac{π}{2})^{1/2}sgn(σ) ist. Formal ergibt also der Faltungssatz

(\frac{1}{s} ⋆ g)^{∧} = (2π)^{1/2}(\frac{1}{s})^{∧}ĝ
= −i π sgn(σ)ĝ

oder, mit der Hilbert-Transformation

(Hg)(s) = \frac{1}{π} ∫ \frac{g(t)}{s − t}dt ,
(Hg)^{∧}(σ) = −i sgn(σ)ĝ(σ) .

Damit erhalten wir

(Kg)^{∧} = −\frac{1}{4π} (Hg′)^{∧}

oder K = −\frac{1}{4π}H \frac{d}{ds}. Es folgt

f(x) = −\frac{1}{4π} ∫_{S^{1}} (Hg′)(θ, x · θ)dσ(θ)

= −\frac{1}{4π^{2}} ∫_{S^{1}} ∫ \frac{g′(θ, s)}{x · θ − s}dsdσ(θ) .

Dies ist die von Radon 1917 angegebene Formel.

4.2 Rekonstruktion einer Funktion aus sphärischen Mitteln

Als nächstes Beispiel eines integralgeometrischen Problems betrachten wir die Rekonstruktion einer Funktion f in IR^{n} aus ihren Integralen über Sphären, deren Mittelpunkte auf einer Hyperebene liegen. Diese Hyperebene nehmen wir als x_{n} = 0 an und haben dann die Integraltransformation

(Rf)(x′, r) = r^{1−n} ∫_{|y|=r} f (( \begin{matrix}x′ \\ 0 \end{matrix} ) + y) dσ(y), x′ ∈ IR^{n−1}, r > 0. (2.1)

Offenbar ist Rf = 0 für in x_{n} ungerade Funktionen f. Wir setzen daher f als gerade in x_{n} voraus. Überdies nehmen wir f ∈ C_{0}(IR^{n}) an, um technische Schwierigkeiten zu vermeiden. Eine äquivalente Darstellung von R ist

(Rf)(x′, r) = 2 ∫_{|y′|<1} f(x′ + ry′, r\sqrt{1 − |y′|^{2}}) \frac{dy′}{\sqrt{1 − |y′|^{2}}} .

Wir machen eine (n − 1)-dimensionale Fourier-Transformation bezüglich x′:

(Rf)^{∧}(ξ′, r) = (2π)^{(n−1)/2} ∫_{IR^{n−1}} e^{−ix′·ξ′}(Rf)(x′, r)dx′

= 2(2π)^{(n−1)/2} ∫_{IR^{n−1}} e^{−ix′·ξ′} ∫_{|y′|<1} f(x′ + ry′, r\sqrt{1 − |y′|^{2}}) \frac{dy′}{\sqrt{1 − |y′|^{2}}} dx′ .

Wir vertauschen die Integrationsreihenfolge und führen für x′ die neue Variable z′ = x′ + ry′ ein. Dies ergibt

(Rf)^{∧}(ξ′, r) = 2(2π)^{(n−1)/2} ∫_{|y′|<1} e^{iry′·ξ′} ∫ f(z′, r\sqrt{1 − |y′|^{2}})e^{−iz′·ξ′}dz′ \frac{dy′}{\sqrt{1 − |y′|^{2}}}

= 2 ∫_{|y′|<1} e^{iry′} f̌(ξ′, r\sqrt{1 − |y′|^{2}})\frac{dy′}{\sqrt{1 − |y′|^{2}}}

mit der (n − 1)-dimensionalen Fourier-Transformation f̌ von f bezüglich der Variablen x′. Hier führen wir für y′ Polarkoordinaten y′ = sθ, θ ∈ S^{n−2} ein.

Dann wird

(Rf)^{∧}(ξ′, r) = 2 ∫_{0}^{1} \frac{s^{n−2}}{\sqrt{1 − s^{2}}} f̌(ξ′, r\sqrt{1 − s^{2}}) ∫_{S^{n−2}} e^{irsθ·ξ′}dθds .  (2.2)

Dies ist für jedes ξ′ ∈ IR^{n−1} eine eindimensionale Integralgleichung für f̌(ξ′, ·). Die Fourier-Transformation hat also das n-dimensionale Problem auf eine (n − 1)-dimensionale Familie eindimensionaler Probleme reduziert.

Wir machen nun Gebrauch von der Hankel-Transformation H_{n}. Diese ist erklärt als die Fourier-Transformation einer radialsymmetrischen Funktion f(|x|) in IR^{n}:

(H_{n}f)(|ξ|) = (2π)^{−n/2} ∫_{IR^{n}} e^{−ix·ξ}f(x)dx

= (2π)^{−n/2} ∫_{0}^{∞}r^{n−1}f(r) ∫_{S^{n−1}} e^{−irθ·ξ}dθdr

= ∫_{0}^{∞}r^{n−1}f(r)(|ξ|r)^{1−n/2}J_{\frac{n−2}{2}} (r|ξ|)dr

mit der Bessel-Funktion J_{k}. Diese Formel haben wir für n = 2 in Aufgabe 15 bewiesen. Es folgt

(H_{n}f)(ρ) = ρ^{1−n/2} ∫_{0}^{∞}r^{n/2}J_{\frac{n−2}{2}} (rρ)f(r)dr .

Für die inverse Fourier-Transformation hätten wir das gleiche Resultat erhalten. Also ist H_{n}^{−1} = H_{n}.

Satz 4.2.1 Sei f ∈ C_{0}(IR^{n}) gerade in x_{n}. Dann gilt

(H_{n}(Rf)^{∧})(ξ′, ρ) = \{ \begin{matrix}2(2π)^{(n−1)/2}f̂(ξ′, \sqrt{ρ^{2} − |ξ′|^{2}})\frac{ρ^{2−n}}{\sqrt{ρ^{2}−|ξ′|^{2}}} , ρ > |ξ′| , \\ 0 , sonst . \end{matrix}

Dabei ist (Rf)^{∧} die (n − 1)-dimensionale Fourier-Transformation von Rf bezüglich des ersten Arguments, und H_{n} wirkt auf das zweite Argument.

Beweis: Wir beschränken uns auf n ≥ 3. Mit Hilfe der Formel

∫_{S^{n−1}} e^{−irθ·ξ}dσ(θ) = (2π)^{n/2}(r|ξ|)^{1−n/2}J_{\frac{n−2}{2}} (r|ξ|)

(vgl. Aufgabe 15) lautet (2.2)

(Rf)^{∧}(ξ′, r) = 2(2π)^{(n−1)/2} ∫_{0}^{1} \frac{s^{n−2}}{\sqrt{1 − s^{2}}} f̌(ξ′, r\sqrt{1 − s^{2}})(rs|ξ′|)^{1−(n−1)/2}J_{\frac{n−3}{2}} (rs|ξ′|)ds .

Die n-dimensionale Hankel-Transformation bezüglich r ergibt (H_{n}(Rf)^{∧})(ξ′, ρ) = 2(2π)^{(n−1)/2}ρ^{1−n/2}

∫_{0}^{1} \frac{s^{n−2}}{\sqrt{1 − s^{2}}} ∫_{0}^{∞}r^{n/2}J_{\frac{n−2}{2}} (rρ)J_{\frac{n−3}{2}} (rs|ξ′|)

(rs|ξ′|)^{1−(n−1)/2}f̌(ξ′, r\sqrt{1 − s^{2}})drds
= 2(2π)^{(n−1)/2}ρ^{1−n/2}|ξ′|^{(3−n)/2} ∫_{0}^{1} \frac{s^{(n−1)/2}}{\sqrt{1 − s^{2}}} ∫_{0}^{∞}r^{3/2}J_{\frac{n−2}{2}} (rρ)J_{\frac{n−3}{2}} (rs|ξ′|)
f̌(ξ′, r\sqrt{1 − s^{2}})drds .

Wir substituieren für r die Variable t = r\sqrt{1 − s^{2}}:

= 2(2π)^{(n−1)/2}ρ^{1−n/2}|ξ′|^{(3−n)/2} ∫_{0}^{1} \frac{s^{(n−1)/2}}{\sqrt{1 − s^{2}}} ∫_{0}^{∞}(\frac{t}{\sqrt{1 − s^{2}}} )^{3/2} J_{\frac{n−2}{2}} (\frac{ρt}{\sqrt{1 − s^{2}}} )

J_{\frac{n−3}{2}} (\frac{st|ξ′|}{\sqrt{1 − s^{2}}} ) ∨ f (ξ′, t)\frac{dtds}{\sqrt{1 − s^{2}}}

= 2(2π)^{(n−1)/2}ρ^{1−n/2}|ξ′|^{(3−n)/2} ∫_{0}^{∞} ∨ f (ξ′, t)t^{3/2}I_{n}(t, ρ, |ξ′|)dt ,

I_{n}(t, ρ, |ξ′|) = ∫_{0}^{1} s^{(n−1)/2}(1 − s^{2})^{−7/4}J_{\frac{n−2}{2}} (\frac{ρt}{\sqrt{1 − s^{2}}} ) J_{\frac{n−3}{2}} (\frac{st|ξ′|}{\sqrt{1 − s^{2}}} ) ds .

In I_{n} substituieren wir für s die Variable x = \frac{s}{\sqrt{1−s^{2}}} , also

x^{2} + 1 = \frac{1}{1 − s^{2}} , s = \frac{x}{\sqrt{x^{2} + 1}} , \frac{dx}{ds} = (1 − s^{2})^{−3/2} .

Das ergibt

I_{n}(t, ρ, |ξ′|) = ∫_{0}^{∞}(\frac{x}{\sqrt{x^{2} + 1}})^{(n−1)/2} (x^{2} + 1)^{7/4}J_{\frac{n−2}{2}} (ρt\sqrt{1 + x^{2}})
J_{\frac{n−3}{2}} (xt|ξ′|)(1 + x^{2})^{−3/2}dx
= ∫_{0}^{∞}x^{(n−1)/2}(x^{2} + 1)^{\frac{2−n}{4}} J_{\frac{n−2}{2}} (ρt\sqrt{1 + x^{2}})J_{\frac{n−3}{2}} (xt|ξ′|)dx .

Dies ist ein Sonine-Gegenbauer-Integral. In: Magnus-Oberkettenger-Soni: Formulas and Theorems for the Special Functions of Mathematical Physics, Springer 1966, S. 104, findet man die Formel (ν > μ > −1)

∫_{0}^{∞}J_{μ}(bx)J_{ν}(a(x^{2} + z^{2})^{1/2})(x^{2} + z^{2})^{−ν/2}x^{μ+1}dx

= \{ \begin{matrix}0 , a < b , \\ b^{μ}a^{−ν}z^{1+μ−ν}(a^{2} − b^{2})^{\frac{ν}{2} − \frac{μ}{2} − \frac{1}{2}} J_{ν−μ−1}(z(a^{2} − b^{2})^{\frac{1}{2}} ) ,  a > b . \end{matrix}

Setzen wir hier z = 1, ν = \frac{n−2}{2}, μ = \frac{n−3}{2}, so ergibt sich

∫_{0}^{∞}J_{\frac{n−3}{2}} (bx)J_{\frac{n−2}{2}} (a\sqrt{x^{2} + 1})x^{\frac{n−1}{2}} (\sqrt{x^{2} + 1})^{\frac{2−n}{2}} dx

= \{ \begin{matrix}0 , a < b \\ b^{(n−3)/2}a^{(2−n)/2}(a^{2} − b^{2})^{−1/4}J_{−1/2}((a^{2} − b^{2})^{1/2}) , a > b , \end{matrix}

= \{ \begin{matrix}0 , a < b \\ \sqrt{\frac{2}{π}}b^{(n−3)/2}a^{(2−n)/2}(a^{2} − b^{2})^{−1/2} cos((a^{2} − b^{2})^{1/2}) , a > b , \end{matrix}

wobei wir zuletzt

J_{−1/2}(z) = \sqrt{\frac{2}{πz}} cos z

benutzt haben. Mit a = ρt, b = t|ξ′| folgt

I_{n}(t, ρ, |ξ′|) = \{ \begin{matrix}0 , |ξ′| > ρ \\ \sqrt{\frac{2}{π}}|ξ′|^{(n−3)/2}ρ^{(2−n)/2}(ρ^{2} − |ξ′|^{2})^{−1/2}t^{−3/2} cos(t(ρ^{2} − |ξ′|^{2})^{1/2}), |ξ′| < ρ . \end{matrix}

Damit erhalten wir schließlich

(H_{n}(Rf)^{∧})(ξ′, ρ) = 2(2π)^{(n−1)/2}\sqrt{\frac{2}{π}} ρ^{2−n} ∫_{0}^{∞}f̌(ξ′, t)\frac{cos(t(ρ^{2} − |ξ′|^{2})^{1/2})}{(ρ^{2} − |ξ′|^{2})^{1/2}}dt

für |ξ′| < ρ und 0 sonst. Da f und damit f̌ im letzten Argument gerade sind, ist das Integral bis auf den Faktor \frac{1}{2}(2π)^{1/2} die Fourier-Transformation von f̌ bezüglich des letzten Argumentes, also f̂. Dieses ergibt die Formel des Satzes.

□

4.3 Lösung integralgeometrischer Aufgaben

mit Symmetrien

Sei G eine abelsche Gruppe, z.B. G = IR^{n}, G = IR_{+}, G = IR mod 2π (Torus-Gruppe). Eine Abbildung χ : G → C\{0} mit χ(s)χ(t) = χ(st) heißt Charakter. Mit χ_{1}, χ_{2} ist auch χ_{1}χ_{2} ein Charakter. Die Charaktere von G bilden wieder eine abelsche Gruppe, die Charaktergruppe oder den Dual Ĝ.

Beispiele:

1) G = IR^{n}, χ_{ξ}(x) = e^{−iξ·x} , Ĝ = IR^{n}.

2) G = IR_{+}, χ_{t}(x) = x^{t} , Ĝ = IR_{+} .

3) G = IR mod 2π , χ_{ℓ}(x) = e^{−iℓx} , Ĝ = Z . Ein Maß μ auf G heißt invariant, wenn μ(sH) = μ(H) für s ∈ G, H ⊂ G. Beispiele:

1) G = IR^{n} , μ = Lebesguesche Maß.

2) G = IR_{+} , dμ(x) = \frac{dx}{x} .

3) G = IR mod 2π, μ = Restriktion des Lebesgueschen Maßes auf [0, 2π). Die Fourier-Transformation einer Funktion f auf G ist nun erklärt als Funktion f̂ auf Ĝ, und zwar ist

f̂(χ) = ∫ χ(s)f(s)dμ(s)

mit dem geeignet normierten invarianten Maß auf G. Beispiele:

1) G = IR^{n} , f̂(ξ) = (2π)^{−n/2} ∫ f(x)e^{−iξ·x}dx.

Dies ist die übliche Fourier-Transformation.

2) G = IR_{+} , f̂(t) = ∫_{0}^{∞}f(x)x^{t} \frac{dx}{x} .

Dies ist die Mellin-Transformation.

3) G = IR mod 2π , f̂(ℓ) = \frac{1}{2π} ∫_{0}^{2π}f(φ)e^{−iℓφ}dφ .

Dies sind die Fourier-Koeffizienten.

Sei nun Rf = g ein n-dimensionales integralgeometrisches Problem. Dieses sei invariant gegenüber einer Gruppe G von Transformationen. Es gilt folgende Regel: Ist G abelsch und von der Dimension d, so kann man Rf = g auf eine d-parametrische Familie (n − d)-dimensionaler Probleme reduzieren. Im Falle d = n bedeutet dies eine Inversionsformel f = R^{−1}g.

Beispiele:

1) Sei R die Radon-Transformation in IR^{n}. Mit f_{a}(x) = f(x + a) hat man dann

(Rf_{a})(θ, s) = (Rf)(θ, s + a · θ) .

Wir haben also Invarianz bezüglich der n-dimensionalen abelschen Gruppe IR^{n} der Translationen x → x + a. Wir erwarten also eine Inversionsformel. In der Tat gilt nach Satz 1.1

ĝ(θ, σ) = (2π)^{(n−1)/2}f̂(σθ)

und dies führt über die Fouriersche Inversionsformel direkt zur Radonschen Inversionsformel, vgl. Satz 1.2.

2) Die Radon-Transformation besitzt noch weitere Invarianzen. Ist U ∈ SO(n) und f_{U}(x) = f(Ux), so ist

(Rf_{U})(θ, s) = ∫_{x·θ=s} f(Ux)dσ(x) = ∫_{y·Uθ=s} f(y)dσ(y)
= (Rf)(Uθ, s) .

Setzt man für t ∈ IR_{+} f_{t}(x) = f(tx), so ist

(Rf_{t})(θ, s) = ∫_{x·θ=s} f(tx)dσ(x)

= t^{1−n} ∫_{y·θ=ts} f(y)dσ(y)
= t^{1−n}(Rf)(θ, ts) .

Wir haben also das direkte Produkt von SO(n) und IR_{+} als Invarianzgruppe. Für n = 2 ist diese Gruppe abelsch und führt mittels Fourier-Reihe und Mellin-Transformation zur Cormackschen Inversionsformel (Aufgabe 30).

3) Sei

(Cf)(x) = ∫_{|y−x|=|x|} f(y)dσ(y) .

Dieses Problem hat die gleiche Invarianzgruppe wie die Radon-Transformation:

(Cf_{U})(x) = ∫_{|y−x|=|x|} f(Uy)dσ(y) = (Cf)(Ux) , U ∈ SO(n)

(Cf_{t})(x) = ∫_{|y−x|=|x|} f(ty)dσ(y) = t^{1−n}(Cf)(tx) .

Für n = 2 muß also eine Entwicklung in eine Fourier-Reihe, gefolgt von der Mellin-Transformation, zu einer Inversionsformel führen. Sei

f(rθ) = ∑_{ℓ} f_{ℓ}(r)e^{iℓφ} , θ = ( \begin{matrix}cos φ \\ sin φ \end{matrix} ) , r > 0 ,

f_{ℓ}(r) = \frac{1}{2π} ∫_{0}^{2π}f(rθ)e^{−iℓφ}dφ .

Mit y = r(θ + ω), ω = ( \begin{matrix}cos ψ \\ sin ψ \end{matrix} ) haben wir

(Cf)(rθ) = ∫_{0}^{2π}f(r(θ + ω))dψ

= ∫_{0}^{2π}f (r ( \begin{matrix}cos φ & + & cos ψ \\ sin φ & + & sin ψ \end{matrix} )) dψ

= ∫_{0}^{2π}f (2r ( \begin{matrix}cos(φ + ψ)/2 \\ sin(φ + ψ)/2 \end{matrix} ) cos \frac{φ − ψ}{2} ) dψ .

Da der Integrand in ψ die Periode 2π hat, können wir anstatt über [0, 2π] ebensogut über [φ − π, φ + π] integrieren. Dann ist cos \frac{φ−ψ}{2} > 0 im Innern des Integrationsintervalls, und wir können die Fourier-Reihe von f benutzen:

(Cf)(r, θ) = ∑_{ℓ} ∫_{φ−π}^{φ+π}f_{ℓ} (2r cos \frac{φ − ψ}{2} ) e^{iℓ(φ+ψ)/2}dψ

= ∑_{ℓ} e^{iℓφ} ∫_{−π}^{π} f_{ℓ} (2r cos \frac{ψ′}{2} ) e^{−iℓψ′/2}dψ′

= ∑_{ℓ} 2e^{iℓφ} ∫_{0}^{π} f_{ℓ}(2r cos \frac{ψ′}{2} ) cos(ℓψ′/2)dψ′

= ∑_{ℓ} e^{iℓφ} \frac{2}{r} ∫_{0}^{2r}f_{ℓ}(s) cos(ℓ arccos(\frac{s}{2r})) \frac{ds}{\sqrt{1 − (\frac{s}{2r})^{2}}} .

Ist also

(Cf)(r, θ) = ∑_{ℓ} g_{ℓ}(r)e^{iℓφ}

die Fourier-Reihe von Cf, so hat man

g_{ℓ}(r) = \frac{2}{r} ∫_{0}^{2r}f_{ℓ}(s) cos(ℓ arccos(\frac{1}{2r})) \frac{ds}{\sqrt{1 − (\frac{s}{2r})^{2}}} .

Dies ist für jedes ℓ eine Integralgleichung, die man durch die Mellin-Transformation lösen kann.

4.4 Näherungsweise Inversion verallgemeinerter Radon-Transformation

Die Radonsche Inversionsformel für n = 3 ist der Prototyp integralgeometrischer Inversionsformeln. Sie lautet

f = −\frac{1}{8π^{2}} R^{∗}g , g = \frac{∂^{2}}{∂s^{2}} Rf

vgl. IV.4.1. Diese Inversionsformel wollen wir näherungsweise auf (III.6.5) übertragen. R^{∗} entspricht dabei einer Integration über alle d mit x ∈ I_{d}, und der Ableitung nach s entspricht die nach t. Also kommt man zu

f(x) = ∫_{R_{x}} g(r, s, t)dσ(r, s, t) ,

wo R_{x} = {(r, s, t) : τ(r, x) + τ(x, s)) = t} die Fläche maximierender Konvexität (vgl. §1) und σ ein geeignetes Maß auf R_{x} ist.

Numerische Methoden

5.1 Das Kaczmarz-Verfahren

Sei A eine (komplexe) n × m-Matrix. Wir interessieren uns für die Lösung von Af = g. Wir nehmen an, es gäbe eine Lösung, die aber nicht eindeutig bestimmt zu sein braucht.

Wir spalten Af = g auf in p Systeme

A_{j}f = g_{j} , j = 1, …, p  (1.1)

mit den (q, m)-Matrizen A_{j} und g_{j} ∈ C^{q}. Es sei n = pq. Die Verallgemeinerung des Verfahrens auf Matrizen A_{j} mit unterschiedlicher Anzahl q_{j} von Zeilen ist trivial. Sei P_{j} die Orthogonalprojektion von C^{m} auf den affinen Unterraum A_{j}f = g_{j}.

Satz 5.1.1 A_{j} habe den Rang q. Dann ist

P_{j}f = f + A_{j}^{∗}S_{j}^{−1}(g_{j} − A_{j}f) S_{j} = A_{j}A_{j}^{∗} .

Beweis: P_{j}f ist charakterisiert durch

A_{j}P_{j}f = g_{j} , P_{j}f − f ⊥ ker (A_{j}) .

Die letzte Bedingung bedeutet wegen C^{m} = ker(A_{j}) ⊕ range (A_{j}^{∗})

P_{j}f = f + A_{j}^{∗}h

mit h ∈ C^{q}. Damit ergibt die erste Bedingung

A_{j}(f + A_{j}^{∗}h) = g_{j} .

Auflösen nach h ergibt die behauptete Formel.

□

Sei nun P_{j}(ω) = (1 − ω)I + ωP_{j}, also P_{j}(1) = P_{j}. Nach Satz 5.1 ist

P_{j}(ω)f = (1 − ω)f + ω(f + A_{j}^{∗}S_{j}^{−1}(g_{j} − A_{j}f))

(1.2)

= f + ωA_{j}^{∗}S_{j}^{−1}(g_{j} − A_{j}f) .

Das Kaczmarz-Verfahren ist ein Iterationsverfahren. Wir erklären den Schritt f^{0} → f^{1}:

f_{0} = f^{0}
f_{j} = P_{j}(ω)f_{j−1} ,  j = 1, …, p  (1.3)
f^{1} = f_{p} .

Für ω = 1 bedeutet dies also die sukzessive orthogonale Projektion auf A_{j}f = g_{j}:

A_{2}f = g_{2} A_{3}f = g_{3}  f_{1} f_{0} = f^{0} A_{1}f = g_{1} f_{2} f^{1} = f_{3}

Die Berechnung von S_{j} ist sehr zeitaufwendig. Wir werden das Verfahren (1.3) daher mit vereinfachten Matrizen S_{j} ausführen, die also nicht notwendig mit A_{j}A_{j}^{∗} übereinstimmen.

Satz 5.1.2 Sei S_{j} positiv definit und S_{j} ≥ A_{j}A_{j}^{∗}. A_{j} habe den Rang q. Ist (1.1) lösbar und erfüllt die Ausgangsnäherung

f^{0} ∈ range (A_{1}^{∗}) ⊕ ⋯ ⊕ range (A_{p}^{∗}) ,

so konvergiert (1.3) gegen die Lösung von (1.1) mit kleinster Norm.

Beweis: Die Iteration lautet explizit

f_{j} = f_{j−1} + ωA_{j}^{∗}S_{j}^{−1}(g_{j} − A_{j}f_{j−1}) .

f_{j} ist also von der Form

f_{j} = f_{0} + A_{1}^{∗}u_{1} + ⋯ + A_{j}^{∗}u_{j}  (1.4)

mit u_{j} ∈ C^{q}. Damit lautet die Rekursion

f_{0} + ∑_{k=1}^{j} A_{k}^{∗}u_{k} = f_{0} + ∑_{k=1}^{j−1} A_{k}^{∗}u_{k}

+ ωA_{j}^{∗}S_{j}^{−1} (g_{j} − A_{j}f_{0} − ∑_{k=1}^{j−1} A_{j}A_{k}^{∗}u_{k})

oder, vereinfacht,

A_{j}^{∗}u_{j} = ωA_{j}^{∗}S_{j}^{−1} (g_{j} − A_{j}f_{0} − ∑_{k=1}^{j−1} A_{j}A_{k}^{∗}u_{k}) .

Da A_{j} vollen Rang hat, ist A_{j}^{∗} injektiv, und man kann den Faktor A_{j}^{∗} kürzen. Nachfolgende Multiplikation mit S_{j} ergibt

S_{j}u_{j} = ω (g_{j} − A_{j}f_{0} − ∑_{k=1}^{j−1} A_{j}A_{k}^{∗}u_{k}) .  (1.5)

Mit

A = ( A_{1} ; ⋮ ; A_{p} ) , S = ( \begin{matrix}S_{1} \\ & . & . & . \\ & & & & S_{p} \end{matrix} ) , g = ( g_{1} ; ⋮ ; g_{p} ) , u = ( u_{1} ; ⋮ ; u_{p} )

haben wir

AA^{∗} = (A_{j}A_{k}^{∗})_{j,k=1,...,p} = L + L^{∗} + D ,

wobei L, D aus den q × q-Blöcken A_{j}A_{k}^{∗} bestehen und L die linke Hälfte, L^{∗} die rechte Hälfte und D die Diagonale von AA^{∗} bedeutet. (1.5) lautet damit

Su = ω(g − Af_{0} − Lu) .

Auflösen nach u ergibt

u = ω(S + ωL)^{−1}(g − Af_{0}) .

Aus (1.4) mit j = p entnehmen wir

f^{1} = f^{0} + A^{∗}u .

Erfüllt f^{0} die Bedingung des Satzes, so ist f^{0} = A^{∗}u^{0} und damit f^{1} = A^{∗}u^{1} mit

u^{1} = u^{0} + u = u^{0} + ω(S + ωL)^{−1}(g − AA^{∗}u^{0}) .

Dies ist von der Form

u^{1} = C(ω)u^{0} + c(ω) ,
C(ω) = I − ω(S + ωL)^{−1}AA^{∗} ,  (1.6)
c(ω) = ω(S + ωL)^{−1}g .

Wir zeigen nun, daß die Eigenwerte von C entweder 1 oder aber betragsmäßig < 1 sind, falls die S_{j} die Voraussetzung des Satzes erfüllen. Sei λ ein Eigenwert mit Eigenvektor u von C(ω), d.h.

(I − ω(S + ωL)^{−1}AA^{∗})u = λu

oder, nach Multiplikation mit S + ωL,

(S + ωL − ωAA^{∗})u = λ(S + ωL)u

Mit AA^{∗} = D + L + L^{∗} lautet dies

(S − ωD − ωL^{∗})u = λ(S + ωL)u .

Wir normieren u so, daß (u, Du) = 1 und setzen (Lu, u) = α + iβ, γ = (Su, u), wobei α, β, γ ∈ IR. Skalare Multiplikation der letzten Gleichung mit u ergibt dann

γ − ω − ω(α − iβ) = λ(γ + ω(α + iβ))
λ = \frac{γ − ω − ω(α − iβ)}{γ + ω(α + iβ)} .

Wegen AA^{∗} ≥ 0 ist

α = \frac{1}{2}((α + iβ) + (α − iβ))
= \frac{1}{2}((Lu, u) + (L^{∗}u, u))
= \frac{1}{2}((AA^{∗}u, u) − (Du, u))
≥ −\frac{1}{2} ,

und wegen S_{j} ≥ A_{j}A_{j}^{∗} ist

γ = (Su, u) ≥ (Du, u) = 1 .

Für α = −\frac{1}{2} ist λ = 1. Für α > −\frac{1}{2} ist |γ − ω − ωα| < |γ + ωα|, wenn nur 0 < ω < 2γ, also erst recht für 0 < ω < 2. Also ist |λ| < 1 für α > −\frac{1}{2} und 0 < ω < 2. In der Tat hat also C(ω) für 0 < ω < 1 nur Eigenwerte 1 oder vom Betrage < 1.

Kehren wir nun zurück zu (1.6). Wir zerlegen C^{n} = H_{0} + H_{1}, wo H_{0} der invariante Unterraum zu den Eigenwerten λ mit |λ| < 1 und H_{1} derjenige zu λ = 1 ist. Wir zeigen, daß H_{1} ⊆ ker (A^{∗}). Ist u_{1} ∈ H_{1} Eigenvektor von C(ω), so ist

(C(ω) − I)u_{1} = −ω(S + ωL)^{−1}AA^{∗}u_{1} = 0

und damit A^{∗}u_{1} = 0. Ist u_{2} ∈ H_{1} Hauptvektor der Stufe 2 von C(ω), so gibt es einen Hauptvektor der Stufe 1, also einen Eigenvektor u_{1} mit

(C(ω) − I)u_{2} = −ω(S + ωL)^{−1}AA^{∗}u_{2} = u_{1} .

Anwendung von A^{∗} ergibt

−ωA^{∗}(S + ωL)^{−1}AA^{∗}u_{2} = 0 .

Da A^{∗} injektiv ist, folgt A^{∗}u_{2} = 0. Ebenso zeigt man A^{∗}u = 0 für Hauptvektoren höherer Stufe und damit H_{1} ⊆ ker (A^{∗}).

Der Spektralradius von C(ω) in H_{0} ist < 1. Die Anteile von u^{k} in H_{0} sind also konvergent. Wegen H_{1} ⊆ ker (A^{∗}) konvergieren damit die f^{k} = A^{∗}u^{k}. Der Grenzwert f erfüllt natürlich alle Gleichungen A_{j}f = g_{j} und liegt in range (A^{∗}). Also ist f Lösung kleinster Norm.

□

Wir wollen nun das Kaczmarz-Verfahren auf das nichtlineare Problem

R_{j}(f) = g_{j} , j = 1, …, p

anwenden. Dabei ist R_{j} ein stetig differenzierbarer Operator von einem Hilbertraum F in einen Hilbertraum G. Damit ist die Fréchet-Ableitung R′(f) ein linearer Operator von F nach G, und wir können seine Adjungierte bilden. Als Verallgemeinerung von (1.3) bietet sich dann an:

f_{j} = f_{j−1} + ωR′_{j}(f_{j−1})^{∗}S_{j}^{−1}(g_{j} − R_{j}(f_{j−1})) .

Dabei ist S_{j} ein linearer Operator, von dem man

S_{j} ≥ R′_{j}(f_{j−1})R_{j}^{′^{∗}}(f_{j−1})

verlangen wird. Dies werden wir im nächsten Abschnitt auf ein inverses Problem anwenden.

5.2 Das Kaczmarz-Verfahren in der Impedanz-Tomographie

Sei Ω ⊆ IR^{n} ein Normalgebiet und σ ∈ C^{1}(\overline{Ω}) sei positiv und σ = 1 auf ∂Ω. Seien f_{j}, g_{j} ∈ C(∂Ω), j = 1, …, p, gegeben. Wir betrachten die Aufgabe

div(σ∇u) = 0 in Ω ,  σ \frac{∂u}{∂ν} = f_{j} auf ∂Ω .  (2.1)

Dabei sei ν die äußere Normale auf ∂Ω. Die Lösung u_{j} ist eindeutig bestimmt durch die Forderung, daß ihr Mittelwert über Ω verschwinden soll. Wir wollen σ näherungsweise aus

u_{j} = g_{j} auf ∂Ω , j = 1, …, p

bestimmen.

Dazu definieren wir den nichtlinearen Operator R_{j} : L_{2}(Ω) → L_{2}(∂Ω) durch

R_{j}(σ) = u_{j}|_{∂Ω}

und haben dann das System

R_{j}(σ) = g_{j} , j = 1, …, p

zu lösen. Zunächst berechnen wir die Fréchet-Ableitung R′_{j}(σ). Dazu berechnen wir für ein h ∈ C_{0}^{1}(Ω) den Wert R_{j}(σ + h). Die zugehörige Lösung von (2.1) sei u_{j} + w_{j}. Es ist also

div ((σ + h)∇(u_{j} + w_{j})) = 0 in Ω
σ \frac{∂(u_{j} + w_{j})}{∂ν} = s_{j} auf ∂Ω .

Vernachlässigen wir die Terme mit hw_{j}, so lautet dies

div(σ∇w_{j}) = −div(h∇u_{j}) in Ω ,

(2.2)

σ \frac{∂w_{j}}{∂ν} = 0 auf ∂Ω .

Wieder erzwingen wir Eindeutigkeit von w_{j} durch die Forderung, daß der Mittelwert verschwindet. Es ist nun

R′_{j}(σ)h = w_{j}|_{∂Ω} .

Wir betrachten R′_{j}(σ) als Operator von L_{2}(Ω) in L_{2}(∂Ω). Dann ist R′_{j}(σ)^{∗} : L_{2}(∂Ω) → L_{2}(Ω). Zu seiner Berechnung gehen wir aus von der 2. Greenschen Formel in der Form

∫_{Ω} (div(σ∇w)z − w div(σ∇z))dx

(2.3)

= ∫_{∂Ω} σ (\frac{∂w}{∂ν}z − w\frac{∂z}{∂ν} ) ds .

Sei nun z die Lösung mit Mittelwert 0 von

div(σ∇z) = 0 in Ω ,

(2.4)

\frac{∂z}{∂ν} = −g auf ∂Ω .

Mit diesem z und mit w = w_{j} lautet (2.3)

− ∫_{Ω} div(h∇u_{j})zdx = ∫_{∂Ω} w_{j}gds

oder, da h_{|∂Ω} = 0,

∫_{Ω} h∇u_{j} · ∇zdx = ∫_{∂Ω} w_{j}gds .

Dies bedeutet aber

(h, ∇u_{j} · ∇z)_{L_{2}(Ω)} = (R′_{j}(σ)h, g)_{L_{2}(∂Ω)}

und damit

R′_{j}(σ)^{∗}g = ∇u_{j} · ∇z .

Das Kaczmarz-Verfahren mit S_{j} = γI lautet damit:

Für j = 1, …, p:

1) Bereitstellung von R_{j}(σ_{j−1}) durch Lösen von

div(σ_{j−1}∇u_{j}) = 0 in Ω , \frac{∂u_{j}}{∂ν} = f_{j} auf ∂Ω .

Es ist dann R_{j}(σ_{j−1}) = u_{j}|_{∂Ω}.

2) Anwendung des Adjungierten auf R_{j}(σ_{j−1}) − g_{j} durch Lösen von

div(σ_{j−1}∇z_{j}) = 0 in Ω , \frac{∂z_{j}}{∂ν} = u_{j} − g_{j} auf ∂Ω .

3) Update

σ_{j} = σ_{j−1} + \frac{w}{γ}∇u_{j} · ∇z_{j} .

5.3 Layer - Stripping

Sei f die bei einem inversen Problem gesuchte Funktion. Mit Layer - Stripping bezeichnet man ein Verfahren, bei dem f zunächst in der Nähe des Randes von supp (f) rekonstruiert wird. Sind die randnahen Teile von f bekannt, rekonstruiert man eine etwas weiter im Innern von supp (f) gelegene Schicht. So fährt man fort, bis das ganze f rekonstruiert ist.

Wir geben drei Beispiele.

1. Radon - Transformation

Sei supp (f) ⊆ {x ∈ IR^{2} : |x| < ρ}, und sei g = Rf die Radon-Transformation, vgl. 4.1. Wir machen Gebrauch von der Cormackschen Inversionsformel aus Aufgabe 30. Sei

g(θ, s) = ∑_{ℓ} g_{ℓ}(s)e^{iℓφ} , θ = ( \begin{matrix}cos φ \\ sin φ \end{matrix} )
f(rθ) = ∑_{ℓ} f_{ℓ}(r)e^{iℓφ} .

Satz 5.3.1 Für ℓ ∈ IR gilt

f_{ℓ}(r) = −\frac{1}{π} ∫_{r}^{∞}(s^{2} − r^{2})^{−1/2}T_{|ℓ|}(\frac{s}{r})g′_{ℓ}(s)ds .  (3.1)

Da hier nur über s ≥ r integriert wird, kann man f(x) für |x| ≥ r aus g(θ, s) für s ≥ r berechnen. Dies ist also ein besonders einfacher Fall von Layer-Stripping. Die randnahen Teile von f werden zur Bestimmung der tiefergelegenen Schichten gar nicht benötigt. Wir können hier aber eine mögliche Gefahr des Layer - Stripping erkennen, nämlich Instabilität. Es gilt nämlich für u ≥ 1 und ℓ ≥ 0

T_{ℓ}(u) = cosh(ℓ arccos u) ≥ \frac{1}{2}e^{ℓ arccosh u}
≥ \frac{1}{2}(u + \sqrt{u^{2} − 1})^{ℓ} ,

so daß also T_{ℓ}(u) für u > 1 mit ℓ → ∞ sehr schnell anwächst. Dies zeigt, daß die Auswertung von (3.1) hochgradig instabil ist.

2. Impedanz - Tomographie

Wir betrachten wieder - wie schon in 5.1 und 1.1 - das Problem, die Funktion σ aus

div(σ∇u) = 0 in Ω

(3.2)

σ \frac{∂u}{∂ν} = f , u = g auf ∂Ω

bei Kenntnis “vieler” Paare f, g zu bestimmen. Wir nehmen jetzt Ω = {x ∈ IR^{2} : |x| < 1}. “Viele” Paare bedeutet, daß die Neumann-zu-Dirichlet - Abbildung S : L_{2}(∂Ω) → L_{2}(∂Ω), g = Sf bestimmt ist. Dies ist z.B. der Fall, wenn die g_{n} = Sf_{n} für ein vollständiges Orthonormalsystem {f_{n}} bekannt

sind, etwa f_{n}(θ) = e^{inφ} für θ = ( \begin{matrix}cos φ \\ sin φ \end{matrix} ), n ∈ Z.

Beim Layer - Stripping wird (3.2) in eine Schar von Aufgaben in den Gebieten |x| < ρ, 0 < ρ ≤ 1 eingebettet:

div(σ∇u_{ρ}) = 0 in |x| < ρ

(3.3)

σ \frac{∂u_{ρ}}{∂ν} = f , u_{ρ} = g auf |x| = ρ .

Hier ist ν (wie oben) die äußere Normale. In Polarkoordinaten (r, φ) ist also ∂/∂ν = ∂/∂r. Wir betrachten jetzt die Funktionen f, g als Funktionen auf S^{1}. Die Neuman-zu-Dirichlet - Abbildung S(ρ) wird dann zu einer Abbildung in L_{2}(S^{1}).

Das Hilfsmittel zum Layer - Stripping ist eine Differentialgleichung für S(ρ). Wir leiten sie her, ohne auf mathematische Exaktheit zu achten.

Satz 5.3.2 Der Operator S′(ρ) = lim_{h→0}(S(ρ + h) − S(ρ))/h erfüllt die Riccati - Differentialgleichung

S′(ρ) = \frac{1}{σ_{ρ}} I + \frac{1}{ρ}S(ρ) + \frac{1}{ρ^{2}} S(ρ)\frac{∂}{∂φ} (σ_{ρ} \frac{∂}{∂φ}S(ρ)) .

Dabei ist σ_{ρ}(θ) = σ(ρθ), θ = ( \begin{matrix}cos φ \\ sin φ \end{matrix} ) und I die Identität auf S^{1}. Beweis: Es ist

(S(ρ)f)(θ) = u_{ρ}(ρθ)

mit der Lösung u_{ρ} von (3.3). Wir berechnen zunächst die Ableitung u′_{ρ}(x) = \frac{∂}{∂ρ}u_{ρ}(x). Sei 0 < ρ_{2} < ρ_{1}. Dann ist

div(σ∇(u_{ρ_{1}} − u_{ρ_{2}})) = 0 in |x| < ρ_{2} ,

und für |x| = ρ_{2} ist bis auf O ((ρ_{2} − ρ_{1})^{2})

(σ \frac{∂(u_{ρ_{1}} − u_{ρ_{2}})}{∂ν}) (ρ_{2}θ) = (σ \frac{∂u_{ρ_{1}}}{∂r} ) (ρ_{2}θ) − f(θ)
= (σ \frac{∂u_{ρ_{1}}}{∂r} ) (ρ_{1}θ) + \frac{∂}{∂r} (σ \frac{∂u_{ρ_{1}}}{∂r} ) (ρ_{2}θ)(ρ_{2} − ρ_{1}) − f(θ)
= \frac{∂}{∂r} (σ \frac{∂u_{ρ_{1}}}{∂r} ) (ρ_{2}θ)(ρ_{2} − ρ_{1}) ,

weil u_{ρ_{1}} wie u_{ρ_{2}} die Randbedingung aus (3.3) erfüllt. Für ρ_{2} → ρ_{1} = ρ folgt

div(σ∇u′_{ρ}) = 0 in |x| < ρ ,

(3.4)

σ \frac{∂u′_{ρ}}{∂ν} = −\frac{∂}{∂r} (σ \frac{∂u_{ρ}}{∂r} ) auf |x| = ρ .

Die Differentialgleichung in (3.3) lautet in Polarkoordinaten

\frac{∂}{∂r} (σ \frac{∂u}{∂r}) + \frac{σ}{r} \frac{∂u}{∂r} + \frac{1}{r^{2}} (σ \frac{∂u}{∂φ}) = 0 .

Damit gilt auf |x| = ρ

σ \frac{∂u′_{ρ}}{∂ν} = \frac{σ_{ρ}}{ρ} \frac{∂u_{ρ}}{∂r} + \frac{1}{ρ^{2}} \frac{∂}{∂φ} (σ_{ρ} \frac{∂u_{ρ}}{∂φ} )
= \frac{1}{ρ}f + \frac{1}{ρ^{2}} \frac{∂}{∂φ} (σ_{ρ} \frac{∂}{∂φ}S(ρ)f) .

(3.4) bedeutet nun nichts anderes als

u′_{ρ} = \frac{1}{ρ}S(ρ)f + \frac{1}{ρ^{2}} S(ρ)\frac{∂}{∂φ} (σ_{ρ} \frac{∂}{∂φ}S(ρ)f) .  (3.5)

Aus

(S(ρ)f)(θ) = u_{ρ}(ρθ)

bekommen wir durch Differentiation

(S′(ρ)f)(θ) = \frac{d}{dρ}u_{ρ}(ρθ)
= u′_{ρ}(ρθ) + \frac{∂u_{ρ}}{∂r} (ρθ)
= u′_{ρ}(ρθ) + \frac{1}{σ_{ρ}} f(θ) .

Zusammen mit (3.5) ist dies die behauptete Differentialgleichung.

□

Wir illustrieren Satz 2 an dem Beispiel σ = 1. Für f(θ) = e^{inφ}, θ = ( \begin{matrix}cos φ \\ sin φ \end{matrix} ) kann man (3.3) leicht lösen, und zwar ist für x = rθ

u_{ρ}(x) = \frac{ρ}{n}(\frac{r}{ρ})^{n}e^{inφ} , n ≠ 0 .

Für r = ρ erhalten wir

S(ρ)e^{inφ} = \frac{ρ}{n}e^{inφ} .

Man berechnet sofort

\frac{1}{ρ^{2}} S(ρ)\frac{∂}{∂φ} (σ_{ρ} \frac{∂}{∂φ}S(ρ)) e^{inφ} = −e^{inφ} ,

und die behauptete Differentialgleichung für S ist offenbar erfüllt. Der Layer - Stripping - Algorithmus arbeitet nun wie folgt. Wir schreiben F (ρ, S) für die rechte Seite der Differentialgleichung aus Satz 2, also

F (ρ, S) = \frac{1}{σ_{ρ}} I + \frac{1}{ρ}S + \frac{1}{ρ^{2}} S \frac{∂}{∂φ} (σ_{ρ} \frac{∂}{∂φ}S) ,

Γ_{l} Γ_{+}  Γ_{−} V Γ_{l+1} φ_{0} ϵ

Abbildung 5.1: Skizze zur Bestimmung von σ_{l+1} aus σ_{l}

so daß die Differentialgleichung einfach

S′ = F (ρ, S)

lautet. S(1) ist bekannt. Wir setzen ρ_{ℓ} = 1 − hℓ, ℓ = 0, …, p − 1, h = \frac{1}{p}, und ersetzen die Differentialgleichung durch das Euler - Verfahren

S_{ℓ+1} = S_{ℓ} − hF (ρ_{ℓ}, S_{ℓ}) ,  ℓ = 0, 1, …, p − 1 ,

(3.6)

S_{0} = S(1) .

Wir benötigen noch eine Formel zur Berechnung von σ_{ℓ}(θ) = σ(ρ_{ℓ}θ). Seien σ_{ℓ}, S_{ℓ} bekannt. Dann können wir zunächst S_{ℓ+1} nach (3.6) berechnen und danach σ_{ℓ+1} in folgender Weise. Sei V = {(r, φ) : ρ_{ℓ−1} < r < ρ_{ℓ}, |φ−φ_{0}| < ε} mit beliebigem φ_{0}, vgl. Abb. 5.1. Sei g ∈ L_{2}(S^{1}) beliebig. Sei u_{ℓ} = u_{ρ_{ℓ}} die Lösung von (3.3) mit u_{ℓ}(ρ_{ℓ}θ) = g(θ), also f = S_{ℓ}^{−1}g = σ_{ℓ}∂u_{ℓ}/∂r. Dann ist nach dem Gaußschen Integralsatz

∫_{∂V} σ \frac{∂u_{ℓ}}{∂ν} ds = 0 .  (3.7)

Auf |x| = ρ_{ℓ+1} ist bis auf Terme der Ordnung h^{2}

u_{ℓ}(ρ_{ℓ+1}θ) = u_{ℓ}(ρ_{ℓ}θ) − h\frac{∂u_{ℓ}}{∂r} (ρ_{ℓ}θ)

(3.8)

= g(θ) − \frac{h}{σ_{ℓ}(θ)}f(θ) .

Damit ist auf |x| = ρ_{ℓ+1}, ebenfalls bis auf O(h^{2})

σ_{ℓ+1} \frac{∂u_{ℓ}}{∂ν} = S_{ℓ+1}^{−1}(g − \frac{h}{σ_{ℓ}} f) .  (3.9)

Sei nun Γ_{ℓ} der Teil von ∂V mit |x| = ρ_{ℓ}, Γ_{±} derjenige mit φ = φ_{0} ± ε. Da auf Γ_{±} ∂/∂ν = ±\frac{1}{r}∂/∂φ gilt, lautet (3.7) wegen (3.9)

∫_{Γ_{ℓ}} fds − ∫_{Γ_{ℓ+1}} S_{ℓ+1}^{−1}(g − \frac{h}{σ_{ℓ}} f)ds + ∫_{Γ_{−}} \frac{σ}{r} \frac{∂u_{ℓ}}{∂φ} ds − ∫_{Γ_{+}} \frac{σ}{r} \frac{∂u_{ℓ}}{∂φ} ds = 0 .

Die Integrale über Γ_{±} approximieren wir durch die Trapezregel und benutzen (3.9). Das ergibt schließlich

∫_{Γ_{ℓ}} fds − ∫_{Γ_{ℓ+1}} S_{ℓ+1}^{−1}(g − \frac{h}{σ_{ℓ}} f)ds + \frac{h}{2} \frac{σ_{ℓ}}{ρ_{ℓ}} (g′(φ_{0} − ε) − g′(φ_{0} + ε))

= \frac{h}{2} \frac{σ_{ℓ+1}}{ρ_{ℓ+1}} ((g − \frac{h}{σ_{ℓ}} f)′(φ_{0} + ε) − (g − \frac{h}{σ_{ℓ}} f)′(φ_{0} − ε)) .

Wendet man dies auf geeignete Funktionen g an, kann man σ_{ℓ+1} bestimmen.

Literaturverzeichnis

Kirsch, A.: An Introduction to the Mathematical Theory of Inverse Pro-

blems. Springer 1996. Colton, D. - Kress, R.: Inverse Acoustic and Electromagnetic Scattering

Theory. Springer 1992. Louis, A.K.: Inverse und schlecht gestellte Probleme. Teubner 1989. Ramm, A.R.: Multidimensional Inverse Scattering Problems. Longman 1992. Chadan, K. - Sabatier, P.C.: Inverse Problems in Quantum Scattering

Theory. Second Edition, Springer 1989. Romanov, V.G.: Integral Geometry and Inverse Problems for Hyperbolic

Equations. Springer 1974. Lavrentev, M.M. - Romanov, V.G. - Shishatskii, S.P.: Ill-Posed Pro-

blems of Mathematical Physics and Analysis. AMS 1986. Bui, H.D.: Inverse Problems in the Mechanics of Materials: An Introduc-

tion. CRC Press 1994. Romanov, V.G. - Kabanikhin, S.I.: Inverse Problems for Maxwell’s Equa-

tion. VSP 1994. Levitan, B.M.: Inverse Sturm-Liouville Problems. VNU 1987. Marchenko, V.A.: Sturm-Liouville Operators and Applications. Birkhäuser

1986. Newton, R.G.: Scattering Theory of Waves and Particles. Second Edition,

Springer 1982. Isakov, V.: Inverse Source Problems. AMS 1990. Romanov, V.G.: Inverse Problems of Mathematical Physics. VNU 1987. Lavrentiev, M.M.: Some Improperly Posed Problems of Mathematical Phy-

sics. Springer 1967.

Lax, P.D. - Phillips, R.S.: Scattering Theory. Academic Press 1967. Gladwell, G.M.L.: Inverse Problems in Scattering. Kluwer 1993.