Diskussion eines Ansatzes zur Lösung des inversen Streuproblems

Diplomarbeit

nach einer Themenstellung von Prof. Dr. F. Natterer

Institut für Numerische und instrumentelle Mathematik

der Westfälischen Wilhelms-Universität Münster Einsteinstraße 62 D-4400 Münster

Frank Wübbeling Niemergs Weide 17 4420 Coesfeld

Münster, den 28. November 1990

Kapitel 1

Einleitung

Thema dieser Diplomarbeit ist die Analyse eines Ansatzes von Kristensson und Vogel zur Rekonstruktion des Umrisses eines Körpers aus der Kenntnis von vom Körper reflektierten Schallwellen bei gegebenen einfallenden Wellen.

Hier zunächst die Geometrie des Problems:

Das physikalische Problem stellt sich wie folgt dar: Gegeben sei ein Objekt (d.h. Teilmenge des IR^{n},n=2,3). Auf das Objekt fallen Wellen der Form w(r, t) = e^{ikr·α̂}e^{−iωt}. Hierbei ist r der Ortsvektor eines Punktes im IR^{n}, ω ist die Frequenz der Welle, t die Zeitkomponente, k = ω/c die Wellenzahl, α̂ der auf 1 normierte Richtungsvektor der einfallenden Welle, c die Wellengeschwindigkeit. Nach einiger Zeit schwingt sich die ausfallende Welle in eine Funktion der Form u^{s}(r)e^{−iωt} ein. Da die einfallende und ausfallende Welle offenbar beide den gemeinsamen Term e^{−iωt} besitzen, wird er im folgenden nicht berücksichtigt, d.h. das Problem ist immer zeitunabhängig.

Das direkte Problem besteht in der Berechnung der ausfallenden Welle. Diese wird eindeutig charakterisiert durch ihr Verhalten auf einem Kreis mit (theoretisch unendlich) grossem Radius, beschrieben durch die Fernfeldamplitude F (ω, k, α), ω ∈ S^{n−1}. Daher wird die ausfallende Welle u^{s} im Rahmen dieser Arbeit oft mit ihrer Fernfeldamplitude identifiziert, insbesondere reicht es aus, nur diese Fernfeldamplitude zur Lösung des direkten Problems anzugeben. Damit diese Lösung eindeutig ist, müssen der Umriß des Objekts und Bedingungen auf seinem Rand bekannt sein.

Das inverse Problem, das im Rahmen dieser Arbeit behandelt wird, besteht in der Berechnung von Objekt bzw. Randbedingungen aus gegebenen Daten, die durch die einfallenden und ausfallenden Wellen für einige Einfallsrichtungen charakterisiert werden.

Mathematisch können die Probleme wie folgt formuliert werden:

Sei Ω ∈ IR^{n},n=2,3. Seien u, u^{i}, u^{s} : IR^{n} ↦ C Funktionen mit:

(r = rω, r ∈ IR^{+}, ω ∈ S^{n−1}, α̂ ∈ IR^{n}, ∥α̂∥ = 1)
u = u^{i} + u^{s}  (1.1)
u^{i}(r) = e^{ikr·α̂}  (1.2)
Δu^{s} + k^{2}u^{s} = 0  (1.3)
lim_{r→∞} r^{(n−1)/2} \frac{∂u^{s}(r)}{∂r} − iku^{s} = 0  (1.4)

gleichmäßig für alle ω ∈ S^{n−1}, sowie Randbedingungen auf Ω:

entweder homogene von-Neumann-Randbedingungen:

u(r) = 0  (1.5)

oder homogene Dirichlet-Randbedingungen

\frac{∂u}{∂ν}(r) = 0  (1.6)

oder die gemischte homogene Randbedingung

\frac{∂u}{∂ν}(r) + w(r) · u(r) = 0  (1.7)

Hierbei ist r ein Punkt aus ∂Ω und ν die nach außen gerichtete Normale an ∂Ω im Punkt r. w ist eine Funktion ∂Ω ↦ C.

Das direkte Problem der Bestimmung von u^{s} ist, wenn alle anderen Daten gegeben sind, eindeutig und numerisch stabil lösbar. Dies wird durch die Ausstrahlungsbedingung sichergestellt. Die inversen Probleme der Bestimmung von w oder Ω sind bei einigen Datensätzen ebenfalls eindeutig lösbar (siehe Kapitel 3).

Das vorgestellte Verfahren zur Rekonstruktion von ∂Ω aus einem Satz von Daten u_{ι}^{i} und u_{ι}^{s} benutzt die Tatsache, daß das direkte Problem gut lösbar ist. Wenn die Eindeutigkeit von Ω bei einem gegebenen Datensatz geklärt ist, kann das Rekonstruktionsproblem für Objekte mit C^{2}-stetigem Rand innerhalb eines Kreises um 0 mit Radius R umgeschrieben werden zu:

Sei M die Menge aller Teilmengen des IR^{n} mit C^{2}-stetigem Rand innerhalb eines Kreises um 0 mit Radius R, Ω ∈ M gesucht. Es sei ein Satz von auf das Gebiet Ω einfallenden Wellen u_{ι}^{i} mit zugehörigen ausfallenden Wellen u_{ι}^{s} vorgegeben, ι ∈ I, I eine Indexmenge. Sei weiter F_{ι}, ι ∈ I, eine Schar von Funktionen auf M. F_{ι} ordne jedem Ω ∈ M die der einfallenden Welle u_{ι}^{i}, bezüglich des Gebiets Ω entsprechende, ausfallende Welle u_{ι}^{s} zu.

Dann ist Ω Lösung der Aufgabe

∑_{ι∈I} ||F_{ι}(Ω) − u_{ι}^{s}||_{ι} = 0
= min_{Ω̃∈M} ∑_{ι∈I} ||F_{ι}(Ω̃) − u_{ι}^{s}||_{ι}

|| · ||_{ι} ist hierbei eine geeignete Norm auf L^{2}(C_{R}), C_{R} das Äußere des Kreises um 0 mit Radius R.

Für die numerische Berechnung muß nun nur noch geeignet diskretisiert werden, d.h. wähle Ω aus einer Menge M′ ⊂ M und I endlich, also:

Bestimme ein Ω ∈ M′, so daß gilt:

∑_{ι∈I} ||F_{ι}(Ω) − u_{ι}^{s}||_{ι} = min_{Ω̃∈M′} ∑_{ι∈I} ||F_{ι}(Ω̃) − u_{ι}^{s}||_{ι}  (1.8)

Hierbei ergeben sich folgende Schwierigkeiten:

1. Das Problem ist schlecht gestellt, deshalb muß 1.8 regularisiert werden.

2. F ist nichtlinear, deshalb muß ein angemessenes Minimierungsverfahren zur Lösung von 1.8 benutzt werden.

3. Der Aufwand zur Berechnung von F ist groß. Da wegen 2 aber die Hessematrix von F benötigt wird, muß der Aufwand zur Berechnung von F verringert, die Hessematrix grob approximiert werden.

Das Streuproblem

2.1 Definitionen

Die in diesem Kapitel gegebenen Definitionen und ohne Beweis angegebenen Sätze werden während der gesamten Arbeit benutzt. Der Grundraum ist hierbei entweder der IR^{2} oder der IR^{3}. (Die Definitionen sind möglichst allgemein gehalten, genauere Erklärungen werden beim ersten Auftauchen des Satzes/der Definition in der Arbeit gegeben)

Definition 2.1.1 Sei D eine offene Teilmenge des Grundraums mit C^{2}- stetigem Rand ohne Löcher. Dann besteht das Helmholtz- Problem in der Lösung der (Helmholtz-) Gleichung

Δu + k^{2}u = 0  (2.1)

auf D für eine Funktion u : IR^{n} ↦ C, n = 2, 3, (in irgendeinem Funktionenraum, in dem die Bildung einer zweiten Ableitung sinnvoll ist) bei gegebenem k ∈ C mit Im k ≥ 0. In dieser Arbeit wird das Helmholtzproblem auch als Vorwärtsproblem bezeichnet.

Im Rahmen dieser Arbeit werden nur Lösungen im Sobolevraum H^{2}(D) betrachtet, bei denen die Bildung der Normalenableitung auf dem Rand und die Fortsetzung der Funktion auf den Rand von D möglich ist. Im wesentlichen werden zwei Probleme unterschieden:

Definition 2.1.2 Sei ein Helmholtz-Problem auf dem Gebiet D gegeben. Falls D das Äußere eines Kreises um 0 enthält, heißt das Problem ein äußeres Helmholtz-Problem.

Definition 2.1.3 Sei ein Helmholtz-Problem auf dem Gebiet D gegeben. Falls D in einem Kreis um 0 mit hinreichend großem Radius R enthalten ist, heißt das Problem ein inneres Helmholtz-Problem.

Um die Lösungen einzuschränken, werden meist Randbedingungen auf ∂D gefordert. Hier werden unterschieden:

Definition 2.1.4 Sei u eine Funktion auf D. u erfüllt Dirichlet-Randbedingungen, falls sich u auf den Rand von D fortsetzen läßt und dort gilt u = f. Hierbei sei f eine vorgegebene Funktion auf dem Rand von D. Wenn f nicht besonders spezifiziert ist, wird f = 0 angenommen, also homogene Randbedingungen.

Definition 2.1.5 Sei u eine Funktion auf D. u erfüllt von Neumann-Randbedingungen, falls sich \frac{∂u}{∂ν}, ν die äußere Normale von D, auf den Rand von D fortsetzen läßt und dort \frac{∂u}{∂ν} = f gilt. Hierbei sei f wieder eine vorgegebene Funktion auf dem Rand von D. Wenn f nicht besonders spezifiziert ist, wird f = 0 angenommen, also homogene Randbedingungen.

Definition 2.1.6 Sei u eine Funktion auf D. Seien weiter λ und χ komplexwertige Funktionen auf dem Rand von D. u erfüllt gemischte Randbedingungen, falls sich \frac{∂u}{∂ν}, ν die äußere Normale von D, und u auf den Rand von D fortsetzen lassen und dort λ\frac{∂u}{∂ν} + χu = f gilt.

In den physikalischen Problemen werden für die Gesamtwelle immer homogene Randbedingungen gefordert.

Definition 2.1.7 Um die Eindeutigkeit der Lösung des Helmholtzproblems beim äußeren Problem zu sichern, muß eine Glattheitsbedingung im Unendlichen eingeführt werden. Dies ist die Sommerfeldsche Ausstrahlungsbedingung: (Sei n die Dimension des Grundraums)

lim_{r↦∞} [r^{(n−1)/2} (\frac{∂u}{∂r} − iku)] = 0  (2.2)

Wichtige Grundlösung für die Helmholtzgleichung ist die zugehörige Greensche Funktion. Im IR^{2} hat sie die Form (R = |x − y|):

G(x, y) = \frac{1}{4}iH_{0}^{(1)}(kR)  (2.3)

Hierbei ist H_{0}^{(1)} die Hankelfunktion erster Art zum Index 0. Im IR^{3} ist G(x, y) rechentechnisch besonders angenehm. Es gilt:

G(x, y) = \frac{e^{ikR}}{4πkR}  (2.4)

2.2 Greensche Formeln

Im Verlauf der Beweise werden mehrere Male die Greenschen Formeln benutzt. Sie seien daher hier in der Form angegeben, in der sie benutzt werden:

Sei G ein Gebiet mit ∂G ∈ C^{2}. Sei R(G) im folgenden definiert als Menge aller komplexwertigen Funktionen u für die gilt u ∈ C^{2}(G) ∩ C(Ḡ), und die Normalenableitung

\frac{∂u}{∂ν}(x) = lim_{h↓0} ν(x) · grad u(x − hν(x)), x ∈ ∂G,

existiert gleichmäßig auf ∂G. Sei u ∈ R(G). Dann gelten die Formeln von Green,

∫_{G} uΔvdx = ∫_{∂G} u\frac{∂v}{∂ν} ds − ∫_{G} grad u · grad vdx  (2.5)

und

∫_{G} (uΔv − vΔu)dx = ∫_{∂G} (u\frac{∂v}{∂ν} − v \frac{∂u}{∂ν})ds  (2.6)

Beweis: Für u ∈ C^{2}(G) folgt der Beweis sofort aus dem Gaußschen Integralsatz. Seien nun u, v ∈ R(G).

Nach [CK83], Anmerkungen zu Abschnitt 2.1, sind die Mengen G_{ϵ} mittels

∂G_{ϵ} := {x : x = y − ϵν(y), y ∈ ∂G}  (2.7)

für 0 < ϵ < ϵ_{0} sinnvoll definiert, und das Oberflächenmaß der Kurve ∂G_{ϵ} im Punkt y − ϵν(y), y ∈ ∂G, konvergiert gleichmäßig in y gegen das Oberflächenmaß der Kurve G im Punkt y.

∫_{∂G} uΔv − vΔudx = lim_{ϵ→0} ∫_{G_{ϵ}} uΔv − vΔudx

= lim_{ϵ→0} ∫_{∂G_{ϵ}} (u\frac{∂v}{∂ν} − v \frac{∂u}{∂ν}) ds

= ∫_{∂G} (u\frac{∂v}{∂ν} − v \frac{∂u}{∂ν}) ds.  (2.8)

2.8 gilt, da sowohl die Integrandenfunktion (nach Voraussetzung) als auch das zugehörige Oberflächenmaß gleichmäßig konvergieren.

2.3 Darstellungssätze  für  Lösungen  der

Helmholtz-Gleichung

In diesem Abschnitt sei D immer ein beschränktes offenes Gebiet im IR^{3} mit C^{2}-stetigem Rand (Alle Sätze lassen sich direkt auch im IR^{2} beweisen, die Einfachheit der Greenschen Funktion im IR^{3} macht die Rechnungen aber wesentlich übersichtlicher). B_{x,r} sei immer die Kugel um x mit Radius r, B_{r} = B_{0,r}, Ω_{x,r} = ∂B_{x,r}, Ω_{r} = Ω_{0,r}. R sei immer so groß, daß D ganz im Kreis um 0 mit Radius R liegt.

Satz 2.3.1 Sei u ∈ R(D) eine Lösung der Helmholtzgleichung Δu+k^{2}u = 0 in D, und sei Φ(x, y) die in 2.4 eingeführte Fundamentallösung Φ(x, y) = \frac{1}{4π} \frac{e^{ik|x−y|}}{|x−y|}. Dann gilt:

∫_{∂D} \{u(y)\frac{∂Φ(x, y)}{∂ν(y)} − \frac{∂u}{∂ν}(y)Φ(x, y)\} ds(y) = \{ \begin{matrix}−u(x), & x ∈ D \\ 0, & x ∈ IR^{3} \ D̄ \end{matrix}

(2.9)

Beweis:

Sei x ∉ D̄. Dann ist

0 = ∫_{D} u(y)k^{2}Φ(x, y) − Φ(x, y)k^{2}u(y)dy
= ∫_{D} u(y)ΔΦ(x, y) − Φ(x, y)Δu(y)dy

= ∫_{∂D} u(y)\frac{∂Φ(x, y)}{∂ν} − \frac{∂u}{∂ν}(y)Φ(x, y)ds(y)  (2.10)

Hierbei beziehen sich alle Differentiationen von Φ natürlich immer auf die zweite Variable.

Sei x ∈ D und r so klein daß Ω_{x,r} ganz in D liegt (D ist offen). Mit der zweiten Greenschen Formel 2.6, angewandt auf D \ B_{x,r} gilt:

0 = ∫_{D\B_{x,r}} u(y)k^{2}Φ(x, y) − Φ(x, y)k^{2}u(y)dy

= ∫_{D\B_{x,r}} u(y)ΔΦ(x, y) − v(x, y)Δu(y)dy

= ∫_{∂D+Ω_{x,r}} u(y)\frac{∂Φ(x, y)}{∂ν} − \frac{∂u}{∂ν}(y)Φ(x, y)ds(y)  (2.11)

Die Behauptung folgt nun aus der Definition der Greenschen Funktion Φ mit Hilfe des Zwischenwertsatzes:

Φ(x, y) = \frac{1}{4π} \frac{e^{ikr}}{r} , r = |x − y|

\frac{∂Φ}{∂r}(x, y) = \frac{ikre^{ikr} − e^{ikr}}{4πr^{2}}
= Φ(x, y)(ik − \frac{1}{r})  (2.12)

∫_{|x−y|=r} \frac{∂u}{∂ν}(y)Φ(x, y) − \frac{∂Φ(x, y)}{∂ν}u(y)ds(y) = 4πr^{2} \frac{∂u}{∂ν}(ξ_{r})Φ(x, ξ_{r}) − 4πr^{2} \frac{∂Φ(x, μ_{r})}{∂ν}u(μ_{r})
⟶_{r→0} u(x)  (2.13)

(μ_{r}, ξ_{r} ∈ B_{x,r}), denn Φ(x, y) = \frac{1}{4πr} + O(1) für r gegen 0. (Bemerkung: Es gilt Im k ≥ 0, deshalb ist e^{ikr} von der Ordnung O(1).)

Ein entsprechender Satz läßt sich auch für die Lösungen der Helmholtzgleichung außerhalb des Gebiets D beweisen. Hierzu benötige ich das folgende Lemma:

Lemma 2.3.2 Sei u ∈ R(IR^{3} \ D̄) eine Lösung der Helmholtzgleichung, die die Sommerfeldsche Ausstrahlungsbedingung 2.1.7 erfüllt. Dann gilt:

lim sup_{R→∞} ∫_{Ω_{R}} |u|^{2}ds < ∞.  (2.14)

Beweis:

0 = lim_{R→∞} ∫_{Ω_{1}} R^{2} |\frac{∂u(Rω)}{∂ν} − iku(Rω)|^{2} ds(ω)  (2.15)

= lim_{R→∞} ∫_{Ω_{R}} |\frac{∂u}{∂ν} − iku|^{2} ds

= lim_{R→∞} ∫_{Ω_{R}} |\frac{∂u}{∂ν}|^{2} + |k|^{2}|u|^{2} − 2 Re(iku\frac{∂ū}{∂ν})ds

= lim_{R→∞} ∫_{Ω_{R}} |\frac{∂u}{∂ν}|^{2} + |k|^{2}|u|^{2} + 2 Im(ku\frac{∂ū}{∂ν})ds  (2.16)

2.15 gilt, da die Konvergenz in der Sommerfeldschen Ausstrahlungsbedingung gleichmäßig ist.

Sei nun R so groß, daß ∥x∥ < R ∀x ∈ D, und sei D_{R} := {x ∉ D|∥x∥ < R}. Aus der Greenschen Formel 2.5 für dies Gebiet ergibt sich:

∫_{D_{R}} −k̄^{2}|u|^{2}dx = ∫_{D_{R}} uΔūdx

= ∫_{∂D_{R}} u\frac{∂ū}{∂ν}ds − ∫_{D_{R}} | grad u|^{2}dx

= ∫_{Ω_{R}} u\frac{∂ū}{∂ν}ds − ∫_{∂D} u\frac{∂ū}{∂ν}ds − ∫_{D_{R}} | grad u|^{2}dx

und damit

∫_{Ω_{R}} u\frac{∂ū}{∂ν}ds = − ∫_{D_{R}} k̄^{2}|u|^{2}dx + ∫_{D_{R}} | grad u|^{2}dx + ∫_{∂D} u\frac{∂ū}{∂ν}ds.

Mit k multipliziert:

k ∫_{Ω_{R}} u\frac{∂ū}{∂ν}ds = k ∫_{∂D} u\frac{∂ū}{∂ν} − k̄|k|^{2} ∫_{D_{R}} |u|^{2}dy + k ∫_{D_{R}} | grad u|^{2}dy. (2.17)

Eingesetzt in 2.16:

lim_{R→∞} \{∫_{Ω_{R}} |\frac{∂u}{∂ν}|^{2} + |k|^{2}|u|^{2}ds + 2 Im(k) ∫_{D_{R}} |k|^{2}|u|^{2} + | grad u|^{2}dy\}

= −2 Im (k ∫_{∂D} u\frac{∂ū}{∂ν}ds) . (2.18)

Offensichtlich hängt die rechte Seite der Gleichung nicht mehr von R ab. Sei nun Im(k) ≥ 0. Dann sind die Terme auf der linken Seite alle ≥ 0, d.h. jeder einzelne ist nach oben beschränkt, denn sonst wäre der Limes ∞.

Insbesondere ist lim sup ∫_{Ω_{R}} |u|^{2}ds ≤ −2 Im (k ∫_{∂D} u\frac{∂ū}{∂ν}ds). Aus 2.18 werde

ich in 2.5.1 die Eindeutigkeit des Vorwärtsproblems folgern.

Satz 2.3.3 Sei u ∈ R(IR^{3} \ D̄) eine Lösung der Helmholtzgleichung, die die Sommerfeldsche Ausstrahlungsbedingung 2.1.7 erfüllt. Dann gilt:

∫_{∂D} u(y)\frac{∂Φ(x, y)}{∂ν} − \frac{∂u}{∂ν}(y)Φ(x, y)ds(y) = \{ \begin{matrix}0, & x ∈ D \\ u(x), & x ∈ IR^{3} \ D̄ \end{matrix} (2.19)

Beweis:

∫_{Ω_{R}} u(y)\frac{∂Φ(x, y)}{∂ν} − \frac{∂u}{∂ν}(y)Φ(x, y)ds(y) = ∫_{Ω_{R}} u(y) \{\frac{∂Φ(x, y)}{∂ν} − ikΦ(x, y)\} ds(y)

− ∫_{Ω_{R}} Φ(x, y) \{\frac{∂u}{∂ν}(y) − iku(y)\} ds(y)

=: I_{1} − I_{2}  (2.20)

Anwendung der Cauchy-Schwarzschen Ungleichung:

|I_{1}| ≤ (∫_{Ω_{R}} |u(y)|^{2} ds(y))^{1/2} (∫_{Ω_{R}} |\frac{∂Φ(x, y)}{∂ν} − ikΦ(x, y)|^{2} ds(y))^{1/2}

= O(\frac{1}{R}),  (2.21)

denn nach 2.3.2 ist der erste Term beschränkt, und wegen 2.12 ist (\frac{∂Φ}{∂ν} − ikΦ)^{2} = O(\frac{1}{R^{4}} ), also das Integral über D_{R} noch von der Ordnung O(\frac{1}{R^{2}} ).

|I_{2}| ≤ \{∫_{Ω_{R}} |Φ(x, y)|^{2} ds(y)\}^{1/2} \{∫_{Ω_{R}} |\frac{∂u}{∂ν}(y) − iku(y)|^{2} ds(y)\}^{1/2}
= o(1),  (2.22)

denn nach Definition der Greenschen Funktion gilt Φ(x, y) = O(1/R) , also ist der erste Term von der Ordnung O(1), und die Ausstrahlungsbedingung 2.1.7 liefert für den zweiten Term die Ordnung o(1), also lim_{R→∞} I_{2} = 0. Insbesondere gilt damit

lim_{r→∞} ∫_{Ω_{R}} u(y)\frac{∂Φ(x, y)}{∂ν} − \frac{∂u}{∂ν}(y)Φ(x, y)ds(y) = 0.

Sei nun x ∈ D, und sei wieder D_{R} = B_{R} \ D.

0 = ∫_{D_{R}} u(y)ΔΦ(x, y) − Δu(y)Φ(x, y)ds(y)

= ∫_{Ω_{R}∪∂D} u(y)\frac{∂Φ(x, y)}{∂ν} − \frac{∂u}{∂ν}(y)Φ(x, y)ds(y)

⟶_{R→∞} ∫_{∂D} u(y)\frac{∂Φ(x, y)}{∂ν} − \frac{∂u}{∂ν}(y)Φ(x, y)ds(y).  (2.23)

Sei nun x ∉ D̄, r so klein daß B_{x,r} nicht in D̄ liegt, und E_{R} = D_{R} \ B_{x,r}. Dann gilt:

0 = ∫_{E_{R}} u(y)ΔΦ(x, y) − Δu(y)Φ(x, y)ds(y)

= ∫_{Ω_{R}∪∂D∪Ω_{x,r}} u(y)\frac{∂Φ(x, y)}{∂ν} − \frac{∂u}{∂ν}(y)Φ(x, y)ds(y)

⟶_{R→∞} ∫_{∂D∪Ω_{x,r}} u(y)\frac{∂Φ(x, y)}{∂ν} − \frac{∂u}{∂ν}(y)Φ(x, y)ds(y)

⟶_{r→0} −u(x) + ∫_{∂D} u(y)\frac{∂Φ(x, y)}{∂ν} − \frac{∂u}{∂ν}(y)Φ(x, y)ds(y). (2.24)

Die letzte Identität folgt wie in 2.13.

2.4 Entwicklung von Lösungen der Helmholtzgleichung

Die Darstellungssätze des letzten Abschnitts werden benutzt um zu zeigen, daß Lösungen der Helmholtzgleichung im Äußeren eines beschränkten Gebiets, die die Ausstrahlungsbedingung erfüllen, einfach entwickelt werden können. Insbesondere wird dadurch die Definition der Fernfeldamplitude gegeben. Zunächst aber ein einführendes Lemma:

Lemma 2.4.1 Sei u Lösung der Helmholtzgleichung 2.1.1, u zweimal stetig differenzierbar. Dann ist u analytisch.

Beweis:

Sei y ∈ IR^{3} fest. Betrachte dann die Greensche Funktion Φ der Helmholtzgleichung in den drei Variablen x_{1}, x_{2}, x_{3} im IR^{3}. Sei x ∈ IR^{3} ebenfalls fest, x ≠ y und ϵ < |x − y|. Auf B_{x,ϵ} gibt es dann nach [Els86] einen holomorphen Logarithmus und damit eine holomorphe Wurzelfunktion von f(z) = ∑_{n=1}^{3}(z_{n} − y_{n})^{2}, die für positive reelle Werte von f mit der normalen Wurzelfunktion übereinstimmt. Mit dieser Wurzelfunktion ist die Funktion Φ(·, y) holomorph erweiterbar auf die komplexen Zahlen z mit |z − x| < ϵ. Insbesondere ist die (reelle) Funktion Φ(·, y) damit analytisch im Punkt x.

Sei nun u eine Lösung der Helmholtzgleichung in einem Gebiet D, x ∈ D. Sei ϵ < d(∂D, x)/2. Dann läßt sich u mit den Darstellungssätzen 2.3.1 bzw. 2.3.3 und der oben angegebenen Wurzelfunktion auf B_{x,ϵ} erweitern.

2.3.1 und 2.3.3 erfüllen für x ∈ B_{x,ϵ} die Voraussetzungen für das Differenzieren unter dem Integralzeichen, denn der Rand von D ist kompakt und der Integrand ist gleichmäßig stetig auf dem Kompaktum B̄_{x,ϵ}. Deshalb werden die Ableitungen in den Cauchy-Riemann Differentialgleichungen zu Ableitungen von Φ. Φ erfüllt diese Gleichungen, denn Φ ist analytisch. Also erfüllt u die Cauchy-Riemann Differentialgleichungen, also ist u komplex differenzierbar (die Ableitungen von u waren als stetig vorausgesetzt), also ist u (lokal) holomorph, also (lokal) analytisch.

Satz 2.4.2 Sei u ∈ C^{(2)}(IR^{3} \ D̄) eine Lösung der Helmholtzgleichung, die die Ausstrahlungsbedingung 2.1.7 erfüllt. Sei R_{0} so groß, daß D̄ ganz im Kreis um 0 mit Radius R_{0} liegt. Weiter sei x = rω_{x}, r ∈ IR^{+}, ω_{x} ∈ S^{2}, r > R_{0}. Dann gilt:

u(x) = \frac{e^{ikr}}{r} ∑_{n=0}^{∞} \frac{F_{n}(ω_{x})}{r^{n}}  (2.25)

mit analytischen Funktionen F_{n} auf dem Einheitskreis. Dabei ist die Konvergenz gleichmäßig absolut auf jedem Kompaktum mit r > R_{0} für jedes Element in dem Kompaktum, und es darf beliebig oft unter dem Summenzeichen differenziert werden.

Um die Bedeutung dieses Satzes hervorzuheben, hier zunächst die zwei wichtigsten Folgerungen:

Korollar 2.4.3 Sei u wie in Satz 2.4.2. Dann hat u die Form

u(x) = \frac{e^{ikr}}{r} F_{0}(ω_{x}) + O(\frac{1}{r^{2}} )  (2.26)

F_{0} heißt dann die Fernfeldamplitude von u. (Bemerkung: Es gilt immer Im k ≥ 0, deshalb läßt sich der Exponentialterm mit 1 abschätzen.)

Korollar 2.4.4 Sei u wie in Satz 2.4.2. Dann wird u durch die Fernfeldamplitude eindeutig bestimmt.

Beweis:

Nach [BS81] gilt:

Δu = \frac{1}{r^{2}} (\frac{∂}{∂r} (r^{2} \frac{∂u}{∂r}) + \frac{1}{sin θ} \frac{∂}{∂θ} (sin θ \frac{∂u}{∂θ}) + \frac{1}{sin^{2} θ} \frac{∂^{2}u}{∂ϕ^{2}} )  (2.27)

\frac{∂}{∂r} \frac{e^{ikr}}{r^{n}} = \frac{e^{ikr}}{r^{n}} (ik − \frac{n}{r})

\frac{∂}{∂r} (r^{2} \frac{∂ \frac{e^{ikr}}{r^{n+1}}}{∂r} ) = \frac{∂}{∂r} (r^{2} \frac{e^{ikr}}{r^{n+1}} (ik − \frac{n + 1}{r}))

= \frac{∂}{∂r} (\frac{e^{ikr}}{r^{n−1}} (ik − \frac{n + 1}{r}))
= −k^{2} \frac{e^{ikr}}{r^{n−1}} − \frac{e^{ikr}}{r^{n}} ((n − 1)ik + (n + 1)ik)
+n(n + 1)\frac{e^{ikr}}{r^{n+1}}  (2.28)

Durch Umstellung der Reihe, die ja absolut konvergiert, ergibt sich:

\frac{1}{r^{2}} \frac{∂}{∂r} (r^{2} \frac{∂u}{∂r}) = \frac{1}{r^{2}} \frac{∂}{∂r} (r^{2} \frac{∂}{∂r} ∑_{n=0}^{∞} (F_{n}(ω_{x})\frac{e^{ikr}}{r^{n+1}} ))

= ∑_{n=0}^{∞} F_{n}(ω_{x}) (−k^{2} \frac{e^{ikr}}{r^{n+1}} − \frac{e^{ikr}}{r^{n+2}} 2nik + n(n + 1)\frac{e^{ikr}}{r^{n+3}} )
= −k^{2}u + \frac{e^{ikr}}{r} ∑_{n=1}^{∞} \frac{1}{r^{n+1}} (2nikF_{n} + n(n − 1)F_{n−1})

Nach 2.4.2 darf unter dem Summenzeichen differenziert werden. Damit erhält man:

0 = Δu + k^{2}u
= \frac{e^{ikr}}{r} ∑_{n=1}^{∞} \frac{1}{r^{n+1}} (−2iknF_{n} + n(n − 1)F_{n−1} + BF_{n−1}) (2.29)

Dabei enthält B die nicht von r abhängigen Teile der Differentiation, also

B := \frac{1}{sin θ} \frac{∂}{∂θ} (sin θ \frac{∂}{∂θ} ) + \frac{1}{sin^{2} θ} \frac{∂^{2}}{∂ϕ^{2}}

Da 2.29 eine Potenzreihe in 1/r ist, müssen alle Koeffizienten der Reihe verschwinden, wenn diese konstant 0 sein soll. Damit erhält man die Rekursionsformel

F_{n} = \frac{1}{2ikn} (n(n − 1)F_{n−1} + BF_{n−1}) ,  (2.30)

aus der sofort die Behauptung des Satzes folgt. Insbesondere gilt:

F_{0} = 0 ⇒ u = 0  (2.31)

Beweis: (zu 2.4.2)

u erfüllt die Voraussetzungen von 2.3.3. Damit gilt für |x| > R, R > R_{0}: (S^{2} die Oberfläche der Einheitskugel im IR^{3},y = Rω):

u(x) = ∫_{D_{R}} u(y)\frac{∂Φ(x, y)}{∂ν} − \frac{∂u}{∂ν}(y)Φ(x, y)ds(y)

= R^{2} ∫_{S^{2}} u(y)\frac{∂Φ(x, y)}{∂ω} − \frac{∂u}{∂ω} (y)Φ(x, y)ds(ω)

= \frac{R^{2}}{4π} ∫_{S^{2}} u(y)(2r − 2R cos γ(ω))\frac{1}{ρ} \frac{d}{dρ} \frac{e^{ikρ}}{ρ} − \frac{∂u}{∂ω} (y)\frac{e^{ikρ}}{ρ} ds(ω)

= \frac{R^{2}}{4π} ∫_{S^{2}} u(y)(2r − 2R cos γ(ω))\frac{1}{ρ} \frac{e^{ikρ}}{ρ} (ik − \frac{1}{ρ}) − \frac{∂u}{∂ω} (y)\frac{e^{ikρ}}{ρ} ds(ω)

(2.32)

mit der Definition ρ = |x − y| = \sqrt{R^{2} − 2rR cos γ(ω) + r^{2}}, γ(ω) der Winkel zwischen ω_{x} und ω. Zum Beweis des Satzes ist nun nur noch zu zeigen, daß die von r abhängigen Teile des Integrals multipliziert mit \frac{r}{e^{ikr}} analytisch in 1/r sind.

Mit der Definition z := R/r, r > R, also z < 1, gilt:

ρ = \frac{R}{z}(z^{2} − 2z cos γ + 1)^{1/2} = \frac{R}{z} + ∑_{k=0}^{∞} a_{k}(γ)z^{k}
= r + ∑_{k=0}^{∞} a_{k}(γ)z^{k}  (2.33)

r − ρ = − ∑_{k=0}^{∞} a_{k}(γ)z^{k}

Die Nullstellen des Polynoms unter der Wurzel sind wegen

z^{2} − 2z cos γ + 1 = z^{2} − z {e^{iγ} + e^{−iγ}} + 1
= (z + e^{iγ})(z − e^{iγ})

beide vom Betrag 1, die zugehörige Wurzelfunktion ist also analytisch um 0 mit Konvergenzradius 1.

Für den zweiten Teil des Integrals 2.32 gilt damit:

\frac{re^{ik(ρ−r)}}{ρ} = \frac{e^{ik ∑_{k=0}^{∞} a_{k}(γ)z^{k}}}{(z^{2} − 2z cos γ + 1)^{1/2}} ,

damit ist der zweite Term entwickelbar in z in einem Kreis um 0 mit Radius 1, also nach 1/r für r > R.

Der erste Term in 2.32 enthält ebenfalls \frac{e^{ikρ}}{ρ} sowie \frac{r}{ρ} und \frac{1}{ρ}, diese sind wie gezeigt entwickelbar, also ist der Gesamtterm unter dem Integralzeichen von 2.32 entwickelbar, also auch der Integralterm selbst, denn Potenzreihen konvergieren gleichmäßig absolut auf jedem Kompaktum innerhalb ihres Konvergenzbereichs. Aus den Sätzen für Potenzreihen folgt auch die Differenzierbarkeit unter dem Summenzeichen, denn der Rand ist kompakt und alle Integrandenfunktionen sind stetig. Die entstehenden Koeffizientenfunktionen F_{n}(ω) sind analytisch, denn sie sind nach den Definitionen analytisch in γ, das wiederum analytisch von ω abhängt.

2.5 Die Eindeutigkeit von Lösungen der Helmholtzgleichung

Die Eindeutigkeit ist eine direkte Folgerung aus den Sätzen der letzten beiden Kapitel.

Lemma 2.5.1 Sei k ∈ IR, u ∈ C^{2}(IR^{3} \ D̄) eine Lösung der Helmholtzgleichung, die die Ausstrahlungsbedingung 2.1.7 erfüllt. Außerdem gelte

∫_{Ω_{R}} |u|^{2}ds ⟶_{R→0} 0  (2.34)

Dann ist u = 0.

Beweis: Sei R so groß, daß D ganz im Kreis um 0 mit Radius R liegt. Dann gilt:

∫_{Ω_{R}} |u|^{2}ds = R^{2} ∫_{Ω_{1}} |\frac{e^{ikR}}{R} ∑_{n=0}^{∞} \frac{F_{n}(Rω)}{R^{n}}|^{2}dω

= ∫_{Ω_{1}} |F_{0}|^{2}dω + O(\frac{1}{R})  (2.35)

Damit ist nach 2.34 F_{0} = 0, also gilt mit 2.4.4 u = 0.

Satz 2.5.2 Sei u eine Lösung der Helmholtzgleichung im Äußeren des Gebiets D, die die Ausstrahlungsbedingung 2.1.7 erfüllt. Außerdem gelte

Im (k ∫_{∂D} u\frac{∂ū}{∂ν}ds) = 0.  (2.36)

Dann ist u = 0.

Die wichtigste Folgerung lautet:

Korollar 2.5.3 Sei u eine Lösung des äußeren Helmholtzproblems 2.1.2, die die Ausstrahlungsbedingung 2.1.7 und die Randbedingungen 2.1.4 oder 2.1.5 erfüllt. Dann ist u eindeutig bestimmt.

Beweis des Korollars:

Seien u_{1} und u_{2} zwei Lösungen des äußeren Helmholtzproblems, die dieselben Dirichlet- oder von Neumann-Randbedingungen und die Ausstrahlungsbedingung erfüllen. Dann erfüllt u_{1} − u_{2} ebenfalls die Ausstrahlungsbedingung und homogene Randbedingungen. Nach 2.36 gilt dann u_{1} = u_{2}.

Bemerkung 2.5.4 Die Bedingung 2.36 ist notwendig, d.h. die Helmholtzgleichung ist für allgemeinere Randbedingungen u.U. nicht eindeutig lösbar.

Beweis: (zu 2.5.2) Nach 2.18 gilt für Im k > 0

∫_{D_{R}} |u|^{2}dx ⟶_{R→∞} 0,

also u = 0. Falls k reell, so folgt aus 2.18:

∫_{Ω_{R}} |u|^{2}ds ⟶_{R→∞} 0,

also wiederum u = 0 nach 2.5.1.

2.6 Existenz von Lösungen - Ansatz über Integralgleichungen

Für die Existenz von Lösungen der Helmholtzgleichung werden einige Sätze aus der Potentialtheorie benutzt. Die Beweise finden sich z.B. bei Leis, “Vorlesungen über partielle Differentialgleichungen zweiter Ordnung”([Lei67]).

Lemma 2.6.1 Sei ψ ∈ C(∂D). Definiere v für x ∉ ∂D durch

v(x) := ∫_{∂D} \frac{∂Φ(x, y)}{∂ν(y)}ψ(y)ds(y)  (2.37)

Dann besitzt v folgende Eigenschaften:

1. v kann von D auf D̄ durch die Funktion v_{−} und von IR^{3} \D̄ auf IR^{3} \D durch v_{+} stetig fortgesetzt werden mittels

v_{±} = ∫_{∂D} \frac{∂Φ(x, y)}{∂ν(y)}ψ(y)ds(y) ± \frac{1}{2}ψ(x) =: \frac{1}{2}(±I + K)ψ  (2.38)

für x ∈ ∂D.

2. K ist kompakt.

3. Es gilt für x ∈ ∂D

v_{+} − v_{−} = ψ  (2.39)

4. sowie

\frac{∂ν_{+}}{∂ν} = \frac{∂ν_{−}}{∂ν}  (2.40)

2.40 und 2.39 heißen Sprungrelationen. Die Normalen sind bei v_{+} nach innen und v_{−} nach außen zu nehmen.

5. Sei

(I + K)ψ = 2f  (2.41)

auf ∂D. Dann ist v Lösung des äußeren Dirichletproblems, bei dem auf ∂D die Randbedingung u = f vorgegeben ist.

Beweisidee:

Für konstantes ψ folgt die Behauptung 2.38 aus dem Gaußschen Integralsatz. 2.39 folgt durch Nachrechnen. 2.40 folgt aus 2.38. Zu 2.41: v löst die Helmholtzgleichung und erfüllt die Ausstrahlungsbedingung 2.1.7, denn ψ und Φ sind stetig und Φ ist stetig differenzierbar nach der zweiten Variablen auf jedem Kompaktum in IR^{3} \ D̄, also darf unter dem Integralzeichen differenziert werden. Wegen 2.38 ist die Randbedingung erfüllt.

Satz 2.6.2 Sei (I + K)ψ = 0. Dann ist ψ Einschränkung auf den Rand einer Eigenfunktion u des inneren Neumann-Problems der Potentialtheorie, also

ψ = u|_{∂D}, u ∈ R(D), Δu + k^{2}u = 0 auf D, \frac{∂u}{∂ν} = 0 auf ∂D  (2.42)

Beweis:

Definiere u := \frac{1}{2}Kψ wie im Lemma. Da Φ und ψ stetig sind und Φ innerhalb von D stetig differenzierbar ist, ist u Lösung der Helmholtzgleichung im Inneren von D. Nach 2.38 ist u_{+} = 0 auf ∂D. Da die Lösung des äußeren Problems eindeutig ist, gilt u = 0 auf IR^{3}\D, also insbesondere auch \frac{∂u_{+}}{∂ν} = 0 auf ∂D. Mit 2.40 gilt damit aber auch \frac{∂u_{−}}{∂ν} = 0 auf ∂D. Damit ist u Eigenfunktion des inneren Neumann-Problems der Potentialtheorie zum Eigenwert k^{2}. Da u_{+} = 0 auf ∂D, ist wegen 2.40 ψ = u_{−} auf ∂D. Daraus folgt die Behauptung des Satzes.

Satz 2.6.3 Sei k kein Eigenwert des inneren Neumannproblems. Dann ist das inhomogene äußere Dirichletproblem mit der Randbedingung u = f lösbar.

Nach 2.6.2 gilt N(I + K) = {0}. Da K kompakt ist, folgt aus der Fredholm-Alternative: (I + K)ψ = 2f ist eindeutig lösbar, also besitzt das Problem eine Lösung, die auch eindeutig ist, denn die Eindeutigkeit wurde bereits gezeigt.

Bemerkung 2.6.4 Falls k Eigenwert ist, findet man eine Lösung des äußeren Dirichletproblems durch Aufaddieren von Funktionen der Form

u_{j}(x) := ∫_{∂D} Φ(x, y)ϕ_{j}(y)ds(y)  (2.43)

auf eine Funktion u wie im Lemma.

Bemerkung 2.6.5 Der entsprechende Satz der Existenz einer Lösung für das Neumann-Problem wird genauso bewiesen, in diesem Fall ist jede Lösung des homogenen Problems (I + K)ψ = 0 mit der Definition

Kψ := 2 ∫_{∂D} Φ(x, y)ψ(y)ds(y)  (2.44)

Einschränkung einer Eigenfunktion des inneren Dirichletproblems der Potentialtheorie. Die zugehörigen Sprungrelationen lauten dann:

v_{+} = v_{−}  (2.45)
\frac{∂v_{+}}{∂ν} − \frac{∂v_{−}}{∂ν} = ψ  (2.46)

2.7 Vollständigkeit von Funktionensystemen

Die Idee der vereinfachten Lösungsmethoden für das Vorwärtsproblem lautet:

Es sei eine Lösung u der Helmholtzgleichung mit vorgegebenen Randbedingungen Ru = ψ gesucht. Sei {ϕ_{n}} ein vollständiges System von Lösungen der Helmholtzgleichung. Sei weiter {ϕ_{n}|_{∂D}} ein vollständiges System in L^{2}(∂D). Dann ist bei geeigneter Wahl der Koeffizienten a_{k} die Funktion ∑_{k=0}^{N} a_{k}ϕ_{k} eine Näherung für die Lösung. Die Koeffizienten a_{k} können durch eine normale Vektorraumapproximation berechnet werden.

Das Problem, die Vollständigkeit der Funktionensysteme zu zeigen, wurde von Müller und Kersten gelöst. Es folgt hier eine vereinfachte Darstellung für die äußeren Probleme:

Satz 2.7.1 Sei H_{0}(k, G) die Menge aller Funktionen in R(G) die die äußere Helmholtzgleichung und die Ausstrahlungsbedingung 2.1.7 erfüllen. Weiter sei S ein offenes Gebiet mit S̄ ⊂ D, und ∂D sei C^{1} -stetig. Dann sind die Einschränkungen der Funktionen in H_{0}(k, S) auf ∂D dicht in L^{2}(∂D).

Beweis:

Sei ∂D C^{2}-parametrisierbar über der kompakten Menge Ω ⊂ IR^{2} mittels der Karte u. Sei f in L^{2}(∂D). Sei |Du| die Wurzel der Gramschen Determinante der Funktion u, also

∫_{∂D} f(x)ds(x) = ∫_{Ω} f ∘ u(x)|Du|dx.

Dann ist f ∘ u in L^{2}(Ω), denn aus

∫_{Ω} (f ∘ u)^{2}|Du|dx < ∞

folgt

∫_{Ω} (f ∘ u)^{2}dx < ∞

mit der Cauchy-Schwarzschen Ungleichung (u ist stetig auf Ω, also nach oben beschränkt; ebenso ist |Du| nach unten beschränkt durch eine Konstante > 0). Auf Ω liegen die stetigen Funktionen dicht in L^{2}, es gibt also eine Folge von Funktionen f_{n} aus C^{0}(Ω), die in der L^{2}-Norm gegen f ∘ u konvergiert.

Wieder mit der Cauchy-Schwarzschen Ungleichung folgt wie oben aus der Konvergenz

∫_{Ω} |f_{n} − f ∘ u|^{2}dx → 0

die Konvergenz

∫_{Ω} |f_{n} − f ∘ u|^{2}|Du|dx → 0,

also

∫_{∂D} |f_{n} ∘ u^{−1} − f|^{2}ds → 0.

f_{n} ∘ u^{−1} ist eine Folge von stetigen Funktionen auf ∂D, d.h. die stetigen Funktionen auf ∂D liegen dicht in den L^{2}-Funktionen auf ∂D.

Es reicht daher zu zeigen: Die Einschränkungen von H_{0}(k, S) auf ∂D sind dicht in C^{0}(∂D) bezüglich L^{2}(∂D).

Sei dazu g ∈ C^{0}(∂D). Sei (g, f|_{∂D}) = 0 für alle f ∈ H_{0}(k, S). Zu zeigen ist: g = 0.

Definiere u(x) durch u(x) := ∫_{∂D} \overline{g(y)}Φ(x, y)ds(y) für x ∉ ∂D und h_{x} durch

h_{x}(y) := Φ(x, y) für y ≠ x. Sei nun x ∈ S. Dann ist h_{x} ∈ H_{0}(k, S), also h_{x}|_{∂D} ∈ H_{0}(k, S)|_{∂D}, also u(x) = (g, h_{x}|_{∂D}) = 0 für x ∉ S. u ist analytisch in D \ S, D \ S̄ hat innere Punkte, also ist u = ∇u = 0 auf D \ S, also u_{+} = \frac{∂u_{+}}{∂ν} = 0.

Mit der Sprungrelation 2.45,2.46 gilt damit: u_{−} = 0 auf ∂D, \frac{∂u_{−}}{∂ν} = ḡ auf ∂D. Da die Lösung der homogenen äußeren Dirichletaufgabe eindeutig ist, gilt schon u = 0 auf IR^{3} \ D, also g = 0.

Satz 2.7.2 Unter den gleichen Voraussetzungen gilt auch: Die Einschränkungen der Normalenableitungen der Funktionen aus H_{0}(k, S) auf ∂D liegen dicht in L^{2}(∂D).

Der Beweis verläuft analog zum Beweis des letzten Sates. Definiere u(x)

durch u(x) := ∫_{∂D} \overline{g(y)}\frac{∂Φ(x,y)}{∂ν}ds(y). Dann ist für x ∈ D \ S̄ u(x) =

(g, \frac{∂h_{x}}{∂ν} |_{∂D}) = 0. Es folgt wieder u_{+} = \frac{∂u_{+}}{∂ν} = 0. Mit den Sprungrelationen 2.39,2.40 erhält man diesmal u_{−} = ḡ, \frac{∂u_{−}}{∂ν} = 0. Aus der Eindeutigkeit des homogenen äußeren von-Neumann-Problems folgt u = 0, also g = 0.

Korollar 2.7.3 Das äußere Dirichlet-Problem und das äußere von Neumann-Problem sind eindeutig lösbar.

Das inverse Streuproblem

3.1 Eindeutigkeit des inversen Problems

Sei in diesem Abschnitt immer k > 0.

Nach den entsprechenden Sätzen über die Lösbarkeit und Eindeutigkeit des Vorwärtsproblems stellen sich die Fragen für das inversen Problem. Es wird sich zeigen, daß bei genügend großem Datenraum die Lösungen des inversen Problems zwar eindeutig bestimmt sind, daß das Problem aber andererseits nicht für jeden beliebigen Datensatz lösbar ist.

Die erste Frage lautet natürlich: Aus welchen Daten sollte die Rekonstruktion erstellt werden? Nach Satz 2.4.3 bestimmt die Fernfeldamplitude die reflektierte Welle u bereits eindeutig. Es reicht daher, nur diese zur Berechnung heranzuziehen.

Der wichtigste Satz lautet daher:

Satz 3.1.1 Sei D ganz in einer Kugel vom Radius R um den Nullpunkt enthalten. Dann gibt es eine Zahl N ∈ IN, so daß gilt: ∂D wird eindeutig bestimmt durch die Fernfeldamplituden F (x, k, α_{j}) für x ∈ Ω_{0}, j ∈ {1..n}. Ω_{0} ist eine Menge auf dem Einheitskreis mit Häufungspunkt.

Beweis:

Da F (x, k, α_{j}) in x analytisch ist für x auf dem Einheitskreis, wird durch die Werte auf Ω_{0} schon der Wert auf dem Einheitskreis bestimmt. Sei also F (x, k, α_{j}) auf dem Einheitskreis gegeben und damit u^{s}, die reflektierte Welle, außerhalb eines Kreises mit Radius R, der D enthält. Mit analytischer Fortsetzung kann u^{s} auf IR^{3} \ D fortgesetzt werden ([Sle82]).

Seien nun D_{1} und D_{2} zwei unterschiedliche Objekte mit derselben Fernfeldamplitude und den reflektierten Wellen u_{1}^{s} und u_{2}^{s}. Es gibt zwei Fälle:

a) D̄_{1}∩D̄_{2} = ∅: u^{s} läßt sich offenbar fortsetzen auf ganz IR^{3}, denn für x ∈ D_{1} definiere u^{s}(x) durch u_{2}^{s}(x), für x ∈ D_{2} durch u_{1}^{s}(x). Für x ∉ D_{1} und x ∉ D_{2} müssen die Werte übereinstimmen (Beide Funktionen sind holomorph und stimmen außerhalb des Kreises mit Radius R überein), setze daher hier u^{s}(x) = u_{1}^{s}(x) = u_{2}^{s}(x). Dann ist u^{s} ganz, Lösung der Helmholtzgleichung und erfüllt die Ausstrahlungsbedingung. Dann ist u^{s} = 0 nach dem folgenden Lemma:

Lemma 3.1.2 Sei u Lösung der Helmholtzgleichung im gesamten Raum und erfülle die Ausstrahlungsbedingung 2.1.7. Dann ist u = 0 konstant.

Beweis des Lemmas:

Sei x ∈ IR^{2} oder x ∈ IR^{3}. Sei D ein Kreis von Radius R um den Nullpunkt, R so groß, daß x ∈ D. Nach 2.3.1 gilt, da u Lösung der Helmholtzgleichung im Inneren von D ist,

∫_{∂D} \{u(y)\frac{∂Φ(x, y)}{∂ν} − \frac{∂u}{∂ν}(y)Φ(x, y)\} ds(y) = −u(x)  (3.1)

Andererseit gilt nach 2.19, da u Lösung der Helmholtzgleichung im Äußeren von D ist,

∫_{∂D} \{u(y)\frac{∂Φ(x, y)}{∂ν} − \frac{∂u}{∂ν}(y)Φ(x, y)\} ds(y) = 0  (3.2)

Also ist u = 0.

Da u^{s} auf dem Rand von D die Randbedingung u^{s}+u^{i} = 0 oder \frac{∂(u^{i}+u^{s})}{∂ν} = 0 erfüllen muß, ist dies ein Widerspruch.

b) Sei nun D_{3} := D_{1} ∩ D_{2} ≠ ∅. Sei ohne Einschränkung D_{1} eine echte Teilmenge von D_{3} und D_{4} eine Zusammenhangskomponente von D_{1} \ D_{3}.

(Falls sich D_{1} und D_{2} nur auf den Rändern schneiden, wähle D_{4} als D_{1}). Dann ist u^{i} +u_{2}^{s} auf D_{4} eine Funktion, die die Dirichlet- bzw. von Neumann-Randbedingung und die Helmholtzgleichung erfüllt. Damit ist u^{i} + u_{2}^{s} eine Eigenfunktion des zugehörigen Laplace-Operators zum Eigenwert k^{2}.

Sei nun ein Datensatz von einfallenden Wellen u_{n}^{i} und ausfallenden Wellen u_{n}^{s} auf D_{4} mit unterschiedlichen Einstrahlungsrichtungen α_{n}, n = 1..N gegeben.

Lemma 3.1.3 Sei D ein Gebiet, R so groß daß D ⊂ B_{R}. Sei u_{n}^{i}(x) = exp(ikx · α_{n}) ein Satz von unterschiedlichen einfallenden Wellen mit zugehörigen ausfallenden Wellen u_{n}^{s}, n ∈ {1..N}, die die Sommerfeldsche Ausstrahlungsbedingung erfüllen. Dann sind die Funktionen {u_{n}^{i} + u_{n}^{s}}, n ∈ {1..N}, linear unabhängig außerhalb von B_{R}.

Beweis: (Nach [Ram86], Kapitel II, Theorem 1, und [CS83])

Sei ∑_{n=1}^{N} c_{n}(u_{n}^{i} + u_{n}^{s}) = 0. Die ausfallenden Wellen u_{n}^{s} sind von der Form O(1/r). Damit gilt:

∑_{n=1}^{N} c_{n}e^{ikx·α_{n}} = O(|x|^{−1})

und insbesondere für beliebiges m ∈ {1..N}:

c_{m} + ∑_{n≠m} c_{n}e^{ikx·(α_{n}−α_{m})} = O(|x|^{−1})  (3.3)

Nach [Ram86] gilt für β ∈ IR^{3}, β ≠ 0:

∫_{|x−x_{0}|≤a} e^{ikβ·x}dx = O(a).  (3.4)

Dabei hängen die Konstanten im O nur noch von β ab.

Sei nun zunächst a fest. Sei x_{0} hinreichend groß. Wenn man 3.3 wie in 3.4 über die Kugel um x_{0} mit Radius a integriert und x_{0} gegen ∞ gehen läßt, konvergiert die rechte Seite gegen 0, denn O(|x|^{−1}) konvergiert natürlich gleichmäßig gegen 0.

Mit der Abschätzung 3.4 gilt damit:

\frac{4}{3}πc_{m}a^{3} + O(a) = 0, a → ∞  (3.5)

Diese Gleichung kann nur richtig sein für c_{m} = 0, also sind die u_{n}^{i} +u_{n}^{s} linear unabhängig.

Der Eigenraum des (kompakten!) Laplace-Operators zu einem von 0 verschiedenen Eigenwert im Inneren des kompakten Gebiets D ist endlichdimensional. Sei N_{0} die Dimension dieses Eigenraums. Dann gilt offensichtlich N ≤ N_{0}.

Wenn es nun gelingt, die Größe N_{0} (die noch von D_{4} abhängt) nach oben abzuschätzen, so ist der Beweis vollzogen, denn für N > N_{0} ist dann auch der Fall sich schneidender Gebiete ausgeschlossen.

Nach Voraussetzung gilt, daß D_{1} und D_{2} in einer Kugel mit Radius R liegen. Nach (Courant/Hilbert I, Kap. 6.2, Theorem 2 und 4) gilt: Die Eigenwerte eines Operators mit Dirichlet- oder von Neumann-Randbedingungen auf dem Rand eines im Inneren eines Gebiets G liegenden Gebiets G_{1} sind, wenn sie nach Größe geordnet und mit ihrer Vielfachheit gezählt werden, größer gleich den entsprechenden Eigenwerten des Gebiets G.

Insbesondere: Seien λ_{1}..λ_{p} die Eigenwerte des Laplaceoperators auf dem Kreis um 0 mit Radius R, die kleiner oder gleich k^{2} sind, ihrer Vielfachheit nach gezählt. Seien μ_{1}..μ_{q} die Eigenwerte kleiner oder kleich k^{2} auf D_{4}. Dann ist q ≤ p, und insbesondere ist auch die Dimension des Eigenraums zum Eigenwert k^{2} bezüglich des Laplaceoperators auf dem Gebiet D_{4} durch p beschränkt. Damit ist der Satz bewiesen.

Der Wert p kann sogar angegeben werden. Anhand des Fundamentalsystems 4.5 von Lösungen der Helmholtzgleichung zu festem k sieht man, daß es genau dann eine Eigenfunktion u zu einem Eigenwert k gibt, falls kR eine Nullstelle der sphärischen bzw. zylindrischen Besselfunktionen ist. p ist damit einfach die Anzahl der Nullstellen der zugehörigen Besselfunktion im IR^{2} bzw. IR^{3} im Bereich 0..kR, multipliziert mit der Vielfachheit der zugehörigen Eigenräume.

Satz 3.1.4 Sei D ganz in einer Kugel vom Radius R um den Nullpunkt enthalten. Sei {k_{n}} eine Folge von Zahlen mit Häufungspunkt H ≠ 0. Sei α_{n} eine Folge von Richtungen. Dann wird ∂D eindeutig bestimmt durch die Fernfeldamplituden F (x, k_{n}, α_{n}), n ∈ IN.

Beweis:

Fall a) wird ausgeschlossen wie oben. Bei Fall b) erhält man wie oben zu jedem k_{n}^{2} eine Eigenfunktion des Laplaceoperators. Da der Laplaceoperator aber kompakt ist, hat sein Spektrum keinen Häufungspunkt ≠ 0. Also ist Fall b) ebenfalls ausgeschlossen, und D wird eindeutig bestimmt.

3.2 Schlechtgestelltheit des Inversen Streuproblems

In der Numerik werden drei Forderungen an ein gut gestelltes Problem gestellt:

1. Das Problem muß für jeden gegebenen Datensatz lösbar sein.

2. Das Problem muß für jeden gegebenen Datensatz eindeutig lösbar sein.

3. Das Ergebnis muß stetig von den Daten abhängen.

In diesem Sinne ist das Vorwärtsproblem gut gestellt. Injektivität und Surjektivität sind gezeigt. [CK83] zeigen die stetige Abhängigkeit der Fernfeldamplitude von den gegebenen Randbedingungen.

Für das inverse Problem wurde die Eindeutigkeit gezeigt. Die Lösbarkeit kann dagegen nicht garantiert werden. Die Fernfeldamplituden sind analytische Funktionen, aber nicht jede L^{2}-Funktion ist analytisch, daher ist das inverse Problem nicht für jede L^{2}-Funktion lösbar.

Die instabilen numerischen Ergebnisse legen ausserdem nahe, daß die Ränder der streuenden Mengen nicht stetig von den Fernfeldamplituden abhängen, wobei als Maß auf den Rändern die L^{2}-Norm der Darstellung des Randes gewählt wird (im Sinne dieser Arbeit also die Radiusfunktion r(φ) bei einem Sterngebiet).

Lösung des Streuproblems

4.1 Berechnung des reflektierten Feldes I

In diesem Abschnitt wird die klassische Methode der Lösung des Vorwärtsproblems 2.1.2 vorgestellt. Sie stellte sich im Lauf der Implementation als vom rein rechentechnischen her den Verfahren II und III unterlegen heraus. Im Gegensatz zu diesen wird die Gleichung hier analytisch gelöst, und dann eine finite Approximation berechnet.

Hierzu werden zunächst zwei Fundamentalsysteme von Lösungen der Helmholtzgleichung definiert, eins (v_{n}^{i}) für die Lösungen im Inneren eines Gebiets, eins (v_{n}^{e}) für die Lösungen im Äußeren, die die Ausstrahlungsbedingung 2.1.7 erfüllen. Sei wie üblich r ∈ IR^{2} bzw. r ∈ IR^{3}, und sei r = (R, φ) bzw. r = (R, φ, ϑ) in Polar- bzw. sphärischen Koordinaten, und sei r̂ = r/R.

Im IR^{2}:

v_{m,σ}^{e}(r) = H_{m}^{(1)}(kR)Y_{m,σ}(r̂)  (4.1)
v_{m,σ}^{i}(r) = J_{m}(kR)Y_{m,σ}(r̂)  (4.2)
Y_{m,σ}(r̂) = \frac{1}{2}\sqrt{2 − δ_{m,0}} \{ \begin{matrix}cos(mϕ) & σ = e, \\ sin(mϕ) & σ = o \end{matrix}  (4.3)

(4.4)

Hierbei sind H_{m} die zylindrischen Hankelfunktionen, J_{m} die zylindrischen Besselfunktionen. σ kann die Werte e und o annehmen.

Im IR^{3}:

v_{l,m,σ}^{e}(r) = h_{l}^{(1)}(kR)Y_{l,m,σ}(r̂)  (4.5)
v_{l,m,σ}^{i}(r) = j_{l}(kR)Y_{l,m,σ}(r̂)  (4.6)
Y_{l,m,σ}(r̂) = (\frac{2 − δ_{m,0}}{2π} \frac{2l + 1}{2} \frac{(l − m)!}{(l + m)!})^{\frac{1}{2}} P_{l}^{m}(cos ϕ) \{ \begin{matrix}cos(mϕ) & σ = e, \\ sin(mϕ) & σ = o \end{matrix}

(4.7)

Hierbei sind die h_{m} die sphärischen Hankelfunktionen, die j_{m} die sphärischen Besselfunktionen. Die P_{l}^{m} sind die zugeordneten Legendre-Funktionen. Der Index σ kann wieder die Werte e und o annehmen, m läuft von 0 bis l, l ist frei.

Die Konstanten sind gerade so gewählt, daß Gleichung 4.20 gilt. In den weiteren Berechnungen wird nun immer für die Indexkombinationen der Gesamtindex n benutzt.

Da die v_{n}^{i} bzw. v_{n}^{e} Fundamentalsysteme sind, lassen sich die Funktionen u^{i} und u^{s} nach diesen entwickeln, d.h. es gibt Entwicklungskoeffizienten a_{n} und f_{n}, so daß gilt:

u^{i}(r) = ∑_{n} a_{n}v_{n}^{i}(r)  (4.8)
u^{s}(r) = ∑_{n} f_{n}v_{n}^{e}(r)  (4.9)

Der T-Matrix-Ansatz bestimmt nun eine (unendlichdimensionale) Matrix T mit:

f = T a  (4.10)

Im folgenden nun zunächst der klassische T-Ansatz, wie er von Waterman zuerst benutzt wurde:

u^{s} ist eine Lösung der Helmholtzgleichung im äußeren des Gebiets D, die die Ausstrahlungsbedingung erfüllt. u^{i} ist eine Lösung im Inneren des Gebiets D. Sei Φ wieder die zugehörige Greensche Funktion. Nach den Sätzen 2.3.1 und 2.3.3 gelten dann die Gleichungen:

Sei x ∉ D̄.

u^{s}(x) = ∫_{∂D} Φ(x, y)\frac{∂u^{s}}{∂ν} (y) − \frac{∂Φ(x, y)}{∂ν}u^{s}(y)ds(y)  (4.11)

0 = ∫_{∂D} Φ(x, y)\frac{∂u^{i}}{∂ν} (y) − \frac{∂Φ(x, y)}{∂ν}u^{i}(y)ds(y)  (4.12)

Für Dirichlet-Randbedingungen gilt damit durch Addition der Gleichungen, u = u^{i} + u^{s}:

u^{s}(x) = ∫_{∂D} Φ(x, y)\frac{∂u}{∂ν}(y) − \frac{∂Φ(x, y)}{∂ν}u(y)ds(y)

= ∫_{∂D} Φ(x, y)\frac{∂u}{∂ν}(y)ds(y)  (4.13)

Für von Neumann-Randbedingungen:

u^{s}(x) = ∫_{∂D} Φ(x, y)\frac{∂u}{∂ν}(y) − \frac{∂Φ(x, y)}{∂ν}u(y)ds(y)  (4.14)

= − ∫_{∂D} \frac{∂Φ(x, y)}{∂ν}u(y)ds(y)  (4.15)

Sei x ∈ D.

0 = ∫_{∂D} Φ(x, y)\frac{∂u^{s}}{∂ν} (y) − \frac{∂Φ(x, y)}{∂ν}u^{s}(y)ds(y)  (4.16)

u^{i}(x) = ∫_{∂D} Φ(x, y)\frac{∂u^{i}}{∂ν} (y) − \frac{∂Φ(x, y)}{∂ν}u^{i}(y)ds(y)  (4.17)

Wie oben erhält man für Dirichlet bzw. von Neumann:

u^{i}(x) = ∫_{∂D} Φ(x, y)\frac{∂u}{∂ν}(y)ds(y)  (4.18)

u^{i}(x) = − ∫_{∂D} \frac{∂Φ(x, y)}{∂ν}u(y)ds(y)  (4.19)

Nach [CH53a] besitzt die Greensche Funktion Φ die Entwicklung

G(y, x) = G(x, y) = i ∑_{n} v_{n}^{i}(x)v_{n}^{e}(y)  (4.20)

mit ∥ x ∥<∥ y ∥. Eingesetzt in 4.13, x ∈ D, und 4.18, x ∉ D̄, ergibt sich:

u^{s}(x) = ∫_{∂D} i ∑_{n} v_{n}^{i}(y)v_{n}^{e}(x)\frac{∂u}{∂ν}ds(y)

= ∑_{n} (i ∫_{∂D} v_{n}^{i}(y)\frac{∂u}{∂ν}(y)ds(y)) v_{n}^{e}(x)  (4.21)

u^{i}(x) = ∑_{n} (i ∫_{∂D} v_{n}^{e} \frac{∂u}{∂ν}ds(y)) v_{n}^{i}(x)  (4.22)

Sei nun χ_{n} ein vollständiges Funktionensystem auf dem Rand von D. Dann läßt sich die (unbekannte) Funktion \frac{∂u}{∂ν} auf dem Rand nach den χ_{n} entwikkeln mit Entwicklungskoeffizienten β_{n}, d.h.

\frac{∂u}{∂ν} = ∑_{n} β_{n}χ_{n}  (4.23)

Dann gilt:

u_{s}(x) = ∑_{n} (∑_{n′} i (∫_{∂D} v_{n}^{i}(y)χ_{n′}(y)ds(y)) β_{n′} ) v_{n}^{e}(x)  (4.24)

u_{i}(x) = ∑_{n} −i (∑_{n′} (∫_{∂D} v_{n}^{e}(y)χ_{n′}(y)ds(y)) β_{n′} ) v_{n}^{i}(x) (4.25)

Definiere nun die Matrizen Q^{i} und Q^{e} durch

Q_{n,n′}^{i} = ∫_{∂D} v_{n}^{i}(y)χ_{n′}(y)ds(y)

bzw.

Q_{n,n′}^{e} = ∫_{∂D} v_{n}^{e}(y)χ_{n′}(y)ds(y).

Dann gilt a = Q^{e}β, f = Q^{i}β (formal für unendlichdimensionale Vektoren) und damit

f = Q^{i}(Q^{e})^{−1},  (4.26)

falls die Matrix Q^{e} invertierbar ist.

Für von-Neumann-Randbedingungen lassen sich alle Überlegungen analog durchführen. In diesem Fall wird die Funktion u auf dem Rand in die Funktionen χ_{n} entwickelt, und man erhält die Definitionen

Q_{n,n′}^{i} = ∫_{∂D} \frac{∂v_{n}^{i}}{∂ν} (y)χ_{n′}(y)ds(y)

bzw.

Q_{n,n′}^{e} = ∫_{∂D} \frac{∂v_{n}^{i}}{∂ν} (y)χ_{n′}(y)ds(y).

Die Matrix Q^{e} ist tatsächlich invertierbar. Ihre Kondition hängt von den verwendeten Funktionen χ_{n} ab.

Der vorgestellte Algorithmus wurde zuerst von Waterman entwickelt. Eine Zusammenfassung aller Aspekte der T-Matrix findet sich in Varadan/Varadan. Eine alternative Herleitung wurde 1979 von Visscher entwickelt. Da das Konzept des Beweises auch in Algorithmus II benutzt wird, sei dieser hier kurz für Dirichlet-Randbedingungen skizziert:

Es gilt u = u^{i} + u^{s}, also

∑_{n} a_{n}v_{n}^{i} + f_{n}v_{n}^{e} = 0

auf ∂D. Sei {χ_{n}}, n ∈ IN, ein vollständiges Funktionensystem des L^{2} auf ∂D. Definiere (u, v) := ∫_{∂D} u(x)v̄(x)ds(x). Durch Multiplikation von 4.1 mit

χ̄_{n} und anschließende Integration erhält man

∑_{n} (a_{n}(v_{n}^{i}, χ_{j}) + f_{n}(v_{n}^{e}, χ_{j})) = 0, j ∈ IN  (4.27)

und diese Bedingung ist sogar äquivalent zu 4.1, denn das System χ_{j} ist vollständig. Also gilt Q^{i}a + Q^{e}f = 0 mit Q_{n,n′}^{i} = (v_{n}^{i}, χ_{n′}) und Q_{n,n′}^{e} = (v_{n}^{e}, χ_{n′}). Hierdurch wird eine Matrixmethode zur Bestimmung der ausfallenden Welle definiert. Visscher zeigt, daß bei geeigneter Wahl der Funktionen χ_{j} das dadurch bestimmte Verfahren äquivalent zum T-Matrix-Verfahren ist.

In der numerischen Implementation gibt es folgende wählbare Parameter:

1. Anzahl der Stützstellen für die Integration.

2. Anzahl der berücksichtigten Entwicklungskoeffizienten für die einfallende Welle.

3. Anzahl der berücksichtigten Entwicklungskoeffizienten für die ausfallende Welle.

4. Anzahl der berücksichtigten Entwicklungskoeffizienten für die Welle auf dem Rand des Gebiets D(∑_{n} β_{n}χ_{n}).

4.2 Berechnung des reflektierten Feldes II

Im vorhergenden Kapitel wurde die klassische Berechnung der T-Matrix gezeigt. Die Berechnung ist numerisch nicht problemlos (stark oszillierende Funktionen, gegen die integriert werden muß, verlangen eine feine Diskretisierung des Integrals) und aufwendig, da in großem Maß spezielle Funktionen berechnet werden müssen. Andererseits läßt sich das Problem auf eine sehr einfache Formel zusammenfassen (4.27). Es bietet sich daher an, das Problem abstrahiert von der physikalischen Seite rein numerisch zu betrachten. Gesucht ist nämlich eine Bestapproximation an die Nullfunktion von

u^{i} + ∑_{k=1}^{N} c_{k}u_{k}(Dirichlet)

bzw.

\frac{δ}{δν} u^{i} + ∑_{k=1}^{N} c_{k} \frac{δ}{δν} u_{k}(vonNeumann)

bezüglich der durch das Standardskalarprodukt auf dem Rand des Objekts D definierten Norm. Die u_{k} seien dabei Grundlösungen der Helmholtzgleichung im Äußeren des Gebiets D, die die Ausstrahlungsbedingung erfüllen.

Satz: Es gilt (∑_{i=0}^{N} c_{i}^{∗}u_{i} + u^{i}, u_{k}) = 0 (Dirichlet) bzw. (∑_{i=0}^{N} c_{i}^{∗} \frac{δ}{δν} u_{i} + \frac{δ}{δν} u^{i}, \frac{δ}{δν} u_{k}) = 0 (von Neumann) , c_{i}^{∗} die Koordinaten der Bestapproximation.

Beweis: (Für Dirichlet, von Neumann entsprechend)

Sei f = ∑_{n=1}^{N} c_{n}u_{n} + u^{s}, f^{∗} = ∑_{i=0}^{N} c_{i}^{∗}u_{i} + u^{s} so gewählt daß gilt: f^{∗} steht senkrecht auf der linearen Hülle der u_{k}, also auch auf f − f^{∗}. Damit gilt:

∥f∥^{2} = ∥f^{∗} + f − f^{∗}∥^{2}
= ∥f^{∗}∥^{2} + ∥f − f^{∗}∥^{2}
≥ ∥f^{∗}∥^{2}.  (4.28)

Damit ist f^{∗} eine Bestapproximation (siehe PraMa1 [Nat85a]). (Visscher kommt zu demselben Ergebnis unter Benutzung der allgemeineren Formel 4.27, indem er zeigt, daß die Wahl der χ_{k} als u_{k} zu einer Fehlerminimierung in der Approximationsfunktion führt. Visscher nennt dies Verfahren daher MOOT, Method of optimal truncation.) Sei nun {u_{k}} eine abzählbar unendliche Menge von linear unabhängigen Funktionen. Dann wird für jedes N ∈ IN eine Näherung f_{N}^{∗} definiert.

Offensichtlich gilt ∥u^{i} − f_{N}^{∗}∥ ≥ ∥u^{i} − f_{N+1}∥ ≥ 0, daher bilden die ∥u^{i} − f_{N}^{∗}∥ eine monoton fallende konvergente Folge. Falls die Einschränkungen von u_{n} auf den Rand von D dicht liegen im Raum der Einschränkungen von u^{i} auf den Rand von D, konvergiert das Verfahren gegen die Lösung des Helmholtzproblems.

4.3 Berechnung des reflektierten Feldes III

Die in diesem Abschnitt vorgestellte Methode zur Berechnung des reflektierten Feldes ist äquivalent zur Diskretisierung der Methode von Visscher, bietet jedoch Vorteile bei der Implementation. Sie geht von derselben Idee aus, benutzt jedoch nicht die Approximation sondern die überbestimmte Interpolation der ausfallenden Welle auf dem Rand des Objekts D.

Definition 4.3.1 Es seien eine einfallende Welle u^{i} und eine linear unabhängige Menge M = {u_{i}} von Lösungen der Helmholtzgleichung im äußeren eines Gebiets D gegeben, die die Ausstrahlungsbedingung 2.1.1 erfüllen. Außerdem sei {x_{k}} eine (endliche) Menge von Punkten auf dem Rand von D, k ∈ {1..N}. Sei dann f_{N}^{∗} die Funktion aus der linearen Hülle von M, die die Funktion

∑_{n=1}^{N}(f_{n}^{∗}(x_{k}) + u^{i}(x_{k}))^{2}

(Dirichlet-Randbedingungen) bzw.

∑_{n=1}^{N}(\frac{∂f_{n}^{∗}}{∂ν} (x_{k}) + \frac{∂u^{i}}{∂ν} (x_{k}))^{2}

(von Neumann-Randbedingungen) minimiert. Die so definierte Funktion ist eine Approximation an die zugehörige ausfallende Welle u^{s}.

Mit dieser Definition gilt: eine diskretisierte Version von Methode II und Methode III sind gleichwertig. Die überraschend genau übereinstimmenden Ergebnisse der Verfahren II und III bei einfachen Objekten und kleinen Frequenzen legen folgenden Satz nahe:

Satz 4.3.1 Die Methode der Berechnung des Vorwärtsproblems durch Approximation bezüglich der auf dem Rand der Kurve gegebenen Norm und die Methode der überbestimmten Interpolation im IR^{2} sind äquivalent in folgendem Sinne: Die bei der ersten Methode auftretenden Integrale seien numerisch ausgewertet durch n Stützstellen von gleichem Gewicht, die in gleichem Abstand auf der Kurve liegen. Dann liefert diese Methode dasselbe Ergebnis wie die zweite, interpoliert an den Stützstellen.

Beweis: (für Dirichlet, von Neumann entsprechend)

Methode I verlangt die Lösung der Gleichungen

∫_{∂D} (∑_{k=0}^{n} α_{k}u_{k}(x))\overline{u_{j}(x)}ds(x) = ∫_{∂D} −u^{i}(x)\overline{u_{j}(x)}ds(x)

für einen Vektor α und für j ∈ {0..n}. Nach 6.3 entspricht dies der Berechnung von

∫_{0}^{2π}(∑_{k=0}^{n} α_{k}u_{k}(x))\overline{u_{j}(x)}\sqrt{r′(φ)^{2} + r(φ)^{2}}dφ = ∫_{0}^{2π}−\overline{u^{i}(x)}u_{j}(x)\sqrt{r′(φ)^{2} + r(φ)^{2}}dφ,

x = ( \begin{matrix}r(φ) cos(φ) \\ r(φ) sin(φ) \end{matrix} ). Die Stützstellen seien nun äquidistant auf der Kurve

gewählt. Dann haben sie den Abstand \frac{c}{\sqrt{r′(φ)^{2}+r(φ)^{2}}} . Damit fällt ein Term

aus der Berechnung heraus, und die numerische Integration liefert das Gleichungssystem:

∑_{k=0}^{n} ∑_{m=0}^{N} u_{k}(x_{m})\overline{u_{j}(x_{m})}α_{k} = ∑_{m=0}^{N} −u^{i}(x_{m})\overline{u_{j}(x_{m})}, j ∈ \{0..n\}  (4.29)

Methode III verlangt die kleinste-Quadrate-Lösung der überbestimmten Interpolation, also des Gleichungssystems

∑_{k=0}^{n} u_{k}(x_{j})α_{k} = −u^{i}(x_{j}), j ∈ {0..N}.  (4.30)

Diese ist nach PraMa I gegeben durch die Lösung der Gleichung

A^{⋆}Aα = A^{⋆}b

mit den Definitionen

A_{mj} = u_{j}(x_{m})

und

b_{j} = −u^{i}(x_{j}).

Ausgerechnet ergibt sich:

(A^{∗}b)_{j} = ∑_{m=0}^{N} A^{∗}_{jm}b_{m}

= ∑_{m=0}^{N} −\overline{u_{j}(x_{m})}u^{i}(x_{m})  (4.31)

(A^{∗}Aα)_{j} = ∑_{m=0}^{N} (A^{∗})_{jm} (Aα)_{m}

= ∑_{m=0}^{N} \overline{(A)_{mj}} ∑_{k=0}^{n} (A)_{mk} α_{k}

= ∑_{m=0}^{N} ∑_{k=0}^{n} \overline{u_{j}(x_{m})}u_{k}(x_{m})α_{k}  (4.32)

und damit dasselbe Gleichungssystem wie oben.

Vom rein numerischen Standpunkt aus bietet eine direkte Implementation von Methode III Vorteile gegenüber einer Implementation von Methode II, da bei der Approximation die auftretende Matrix A^{⋆}A ausgewertet und das zugehörige Gleichungssystem tatsächlich gelöst wird. Nach PraMa 1 gilt allerdings, daß diese Matrix die quadrierte Kondition der Matrix A besitzt. Die Interpolation wird dagegen mit dem QR-Verfahren zur direkten Bestimmung der kleinsten Quadrate- Lösung aus der Matrix A gelöst, d.h. ohne explizite Berechnung von A^{⋆}A und damit ohne Konditionsvergrößerung.

In der Praxis zeigte sich, daß die auftretenden Matrizen relativ gut konditioniert sind. Daher wirkte sich dieser Vorteil nur gering aus. Da die in den anderen Verfahren auftretenden Integrationen hier aber nicht durchgeführt werden, ist es numerisch bei stark oszillierenden Integranden stabiler. Dies konnte in der Implementation gut beobachtet werden, bei Benutzung der Verfahren I und II traten häufig “Underflow Errors” auf, die bei Verfahren III nicht beobachtet wurden. Die etwas andere Berechnungsweise führte außerdem in bestimmten Situationen zu einem deutlichen Geschwindigkeitsvorteil von Methode III gegenüber Methode II.

Bemerkung: Im Beweis wurde immer von Dirichlet-Randbedingungen ausgegangen. Durch Übergang zu den Normalenableitungen in den einzelnen Sätzen gilt der Satz auch für von-Neumann- und gemischte Randbedingungen.

Da die numerischen Ergebnisse von Methode I und Methode II für die einfachste Kurve, den Kreis, übereinstimmen, vermutet man hier den Satz:

Satz 4.3.2 Sei D ein Kreis vom Radius R. Dann sind im IR^{2} alle drei Methoden äquivalent im Sinne von 4.3.1, wenn als Funktionenraum {χ_{n}} in Methode I die trigonometrischen Funktionen gewählt werden.

Beweis (für Dirichlet-Randbedingungen, von Neumann entsprechend):

Es ist nur zu zeigen, daß Methode I und II äquivalent sind.

Sei D ein Kreis vom Radius R. Bei Methode I ergibt sich für das Matrixelement Q_{n,n′}^{i}, n = (m, σ), n′ = (m′, σ′):

Q_{n,n′}^{i} = ∫_{∂D} v_{n}^{i}(y)χ_{n′}(y)ds(y)

= ∫_{0}^{2π}Rv_{n}^{i}(y) \{ \begin{matrix}cos(m′φ) \\ sin(m′φ) \end{matrix} \} dϕ

= ∫_{0}^{2π}RJ_{m}(kR) \{ \begin{matrix}cos(mφ) \\ sin(mφ) \end{matrix} \} \{ \begin{matrix}cos(m′φ) \\ sin(m′φ) \end{matrix} \} dϕ

= \{ \begin{matrix}J_{m}(kR)R ∫_{0}^{2π}\{ \begin{matrix}cos(mφ)^{2} \\ sin(mφ)^{2} \end{matrix} \} dϕ = c_{n}RJ_{m}(kR) & n = n′ \\ 0 & sonst \end{matrix}

(4.33)

Die letzte Gleichung folgt aus der Orthogonalität der trigonometrischen Funktion. Mit gleicher Rechnung erhält man für die Matrix Q^{e}:

Q_{n,n′}^{e} = \{ \begin{matrix}c_{n}H_{m}^{(1)}(kR) & n = n′ \\ 0 & sonst \end{matrix}  (4.34)

Die Matrizen Q^{i} und Q^{e} sind also Diagonalmatrizen. Damit lassen sich Matrixmultiplikation und -inversion einfach durchführen, und aus der Beziehung f = Q^{i}(Q^{e})^{−1}a folgt

f_{n} = \frac{J_{m}(kR)}{H_{m}^{(1)}(kR)}a_{n},  (4.35)

a_{n} wie in 6.2 berechnet.

Methode II verlangt die Gültigkeit der Beziehung

(∑_{i} f_{i}v_{i}^{e}, v_{n}^{e}) = (u^{i}, v_{n}^{e})  (4.36)

für n ∈ {1..N}, N die Dimension des Funktionenraumes der {χ_{n}}.

Aus der Orthogonalität der trigonometrischen Funktionen bezüglich des Integrals über einen Kreis folgt wie oben:

(f_{n}v_{n}^{e}, v_{n}^{e}) = (u^{i}, v_{n}^{e}) ,  (4.37)

und damit

f_{n} = \frac{(u^{i}, v_{n}^{e})}{(v_{n}^{e}, v_{n}^{e})}

= \frac{∫_{0}^{2π}Re^{ikR cos(φ−α)}H_{m}^{(1)}(kR) \{ \begin{matrix}cos(mφ) \\ sin(mφ) \end{matrix} \} dφ}{∫_{0}^{2π}RH_{m}^{(1)}(kR)^{2} \{ \begin{matrix}cos(m^{2}φ) \\ sin(m^{2}φ) \end{matrix} \} dφ}

= \frac{1}{H_{m}^{(1)}(kR)} \frac{∫_{0}^{2π}e^{ikR cos(φ−α)} \{ \begin{matrix}cos(mφ) \\ sin(mφ) \end{matrix} \} dφ}{∫_{0}^{2π}\{ \begin{matrix}cos(m^{2}φ) \\ sin(m^{2}φ) \end{matrix} \} dφ}

(4.38)

Der zweite Bruch ist der zu \{ \begin{matrix}cos(mφ) \\ sin(mφ) \end{matrix} \} gehörende Koeffizient der Fourier-

reihe zu e^{ikR cos(φ−α)}, denn durch den Nenner werden die trigonometrischen Funktionen im Zähler normiert. Dieser Koeffizient ist aber in 6.2 berechnet, er hat den Wert

a_{n}J_{m}(kR),  (4.39)

es gilt auch hier

f_{n} = \frac{J_{m}(kR)}{H_{m}^{(1)}(kR)}a_{n},  (4.40)

und damit sind die Enewicklungskeffizienten f_{n} in Methode I und II unter den genannten Bedingungen gleich.

Damit tritt der unangenehme Fall ein, daß ausgerechnet unter den einzigen analytisch durchrechenbaren Bedingungen die Algorithmen dasselbe Ergebnis liefern, so daß keine analytische Aussage über die Vorteile der Algorithmen gemacht werden kann.

Lösung des inversen Streuproblems

5.1 Beschreibung des Verfahrens

Im letzten Kapitel wurde die einfache Lösbarkeit des Vorwärtsproblems gezeigt. Das vorgestellte Verfahren benutzt dies zur Lösung des inversen Problems:

Sei D ein Sterngebiet. Ohne Einschränkung sei 0 ein Sternpunkt von D. Sei der Rand von D parametrisiert über den Raumwinkel, D = r(ω), ∥ω∥ = 1. Dann besteht das inverse Streuproblem in der Bestimmung der (eindeutig festgelegten) Funktion r. Sei A der Operator, der zu einer gegebenen Radiusfunktion r die gemessenen Daten (Fernfeldamplitude, cross-section...) liefert. Seien weiter b die gemessenen Daten in einem Datensatz, der Eindeutigkeit garantiert. Dann ist r die Lösung von:

min_{r∈C^{0}^{+}(0,2π)} ∥Ar − b∥  (5.1)

Natürlich sind im allgemeinen keine kontinuierlichen Daten vorgegeben, wie es die Eindeutigkeitssätze verlangen. Es muß daher eine diskretisierte Version benutzt werden. Diese garantiert nun natürlich keine Eindeutigkeit mehr. Zusätzlich wird nicht im gesamten Raum C^{0} approximiert, sondern nur in einem Teilraum. r_{disc} wird daher gesucht als (u.u. nicht eindeutige) Lösung der Aufgabe

min_{r∈V_{n}} ∥Ar − b∥_{disc}  (5.2)

Die Räume V_{n} werden dabei in dieser Arbeit als trigonometrische Polynome vom Grad ≤ n gewählt.

5.2 Regularisierung der Lösung

Da das inverse Streuproblem schlecht gestellt ist, muß die Lösung regularisiert werden. Unterläßt man dies, muß bei Störungen in den Daten mit großen Fehlern in den Lösungen gerechnet werden, wie bei den praktischen Ergebnissen gezeigt.

Ich habe für den gegebenen Algorithmus zwei Arten der Regularisierung ausprobiert. Die erste ist die von Kristensson/Vogel vorgeschlagene Regularisierung von Tychonoff-Phillips, die zweite ist eine Regularisierung, die den Rechenablauf erheblich vereinfacht und beschleunigt.

I) Regularisierung nach Tychonoff-Phillips: Statt, wie im Algorithmus angegeben, die Funktion |F (x)−d| zu minimieren, wird die Funktion |F (x)−d|+ α|x| minimiert. Hierbei sollte die Norm so gewählt werden, daß die Norm hochfrequenter Funktionen größer ist als die von niederfrequenten Funktionen. Man wählt daher die durch das folgende Skalarprodukt definierte Norm in H^{1}(0, 2π):

(x, y) := u(0)v(0) + ∫_{0}^{2π}x′(t)v′(t)dt  (5.3)

(Kristensson/Vogel benutzen hier eine kompliziertere, über die gebrochenen Sobolevräume definierte Norm. In der Praxis wirkte sich die Wahl der Norm aber nicht aus.) In den Versuchen wird jeweils zuerst ein hoher Regularisierungsparameter gewählt, mit diesem die Optimierung durchgeführt, der Parameter anschließend verkleinert und mit diesem neuen Wert weitergerechnet.

II) Einschränkung des Grundfunktionenraums:

Offensichtlich verhindert die Regularisierung I) das zu starke Anwachsen der hochfrequenten Funktionen zum Beginn der Optimierung, wenn die Näherung noch sehr weit vom Ergebnis entfernt liegt. Da die Anzahl der Funktionen ein wesentlicher Faktor in der Komplexität des Programms ist, lautet die Regularisierungsvorschrift hier:

a) Wähle eine kleine Zahl n als Dimension des Interpolationsraums V_{n}. In den Versuchen wurde n = 5 gewählt.

b) Führe eine Optimierung durch.

c) Erhöhe n. Falls noch kein Maximum erreicht ist, gehe zu b) zurück.

Da der verwendete Optimierungsalgorithmus dazu neigt, bei sehr guter erreichter Näherung nur noch kleine Schritte zu machen, liegt hier tatsächlich eine Regularisierung vor. Ein großer Vorteil dieses Verfahrens liegt in der erheblichen Zeitersparnis, sie lag in der Regel bei ca. 3/4 der Zeit von Verfahren I. Zusätzlich arbeitete der Optimierungsalgorithmus bei kleinen Dimensionen erheblich effektiver als bei großen, so daß zusätzlich eine Beschleunigung des Konvergenzverhaltens des Optimierungsverfahrens zu bemerken war.

III) In der Praxis zeigte auch eine Kombination der beiden Verfahren Vorteile. Mit Verfahren II) kann sehr schnell eine relativ gute Näherung erreicht werden. Diese Näherung wird dann als Startpunkt für Verfahren I gewählt. Als Norm wurde dann dort eine Halbnorm gewählt, die die durch Verfahren II gut bestimmten Koeffizienten kleiner Ordnung nicht bestrafte.

Implementation

In diesem Kapitel wird die Implementation des Verfahrens auf einem konkreten Rechner untersucht. Für alle Abschnitte gilt: Es wurde das Verfahren im IR^{2} implementiert mit Dirichlet- Randbedingungen.

Im letzten Kapitel ist wohl klar geworden, daß eine schnelle Lösung des Vorwärtsproblems Voraussetzung für eine effiziente Lösung des inversen Problems mit dem hier vorgestellten Verfahren ist. Aus diesem Grund wird dies hier besonders genau diskutiert.

Es wurden beide Lösungsansätze implementiert. Beide seien hier noch einmal für 2 Dimensionen definiert:

Verfahren a): Klassisch wie bei Kristensson und Vogel. Lösung durch das T-Matrix-Konzept: Seien (a_{n}) die Entwicklungskoefffizienten der einfallenden Welle u_{i} bezüglich eines Systems von Grundlösungen v_{n}^{i} der Helmholtzgleichung, und f_{n} die Entwicklungskoeffizienten der Lösung des Problems bezüglich eines Systems von Grundlösungen v_{n}^{e}, die die Ausstrahlungsbedingung 2.1.7 erfüllen. Dann gilt die Beziehung f = T a. Wir benutzen die in 4.26 hergeleitete Beziehung f = Q^{i}β und a = Q^{e}β, also f = Q^{i}(Q^{e})^{−1}a.

Verfahren b): In der modifizierten Version von Visscher. Es wurden im wesentlichen zwei Abwandlungen benutzt. Die erste geht über das Konzept der Approximation der Randbedingungen bezüglich der hergeleiteten Norm. Die zweite geht über das Konzept der überbestimmten Interpolation der Randbedingungen. Die Zeitvergleiche und die numerischen Ergebnisse zeigen, daß das letzte Verfahren von der Effizienz her den anderen überlegen ist.

Phantome: Zum Testen des Verfahrens wurden vier Phantome untersucht. Im einzelnen:

1. Der Erdnuß-Umriß (nach Kristensson und Vogel), dieser enthält alle Ellipsen. Das Phantom stellte sich als nicht sehr hilfreich heraus, denn der Effekt der Regularisierung läßt sich nur unzureichend beobachten, weil das Objekt zu glatt ist. (3,4)

2. Pinocchio, ein von Visscher vorgeschlagener Umriß. (1)

3. Quadrate, zur Untersuchung des Verhaltens an nicht C^{2}-stetigen Objekten. (6)

4. Kreisteil, wie oben. (5)

5. Teddy, wie Pinocchio, aber mit zwei Ausbuchtungen. (2)

Zur Erklärung der Zeichnungen: Innen liegt jeweis das Objekt (rekonstruiert oder original). Außen liegt ein Kreis. Auf diesem Kreis werden nun die einfallende Welle (der Übersichtlichkeit wegen nur in einigen Bildern, hier in 3-6), und die resultierende Welle (Real- und Imaginärteil) abgetragen. Also: Abgetragener Wert (φ) = Re u_{ges}(r(φ)), entsprechend mit Im u_{ges} und Re u_{i}.

Bei optimaler Erfüllung der Dirichlet-Randbedingungen sollten nun der Kreis und die Welle u_{ges} zusammenfallen. Dies ist in obigen Bildern auch tasächlich der Fall.

6.1 Entwicklung der einfallenden Welle bei Methode I

Die einfallende Welle u^{i} wird dargestellt durch eine Reihe der Form

∑_{n=0}^{∞} a_{n}J_{m}(kr)Y_{n}(r̂).

Die Funktionen Y_{n} sind von der Form cos(mφ) bzw. sin(mφ). Damit ist die Reihe eine Fourierreihe in φ. Sei u_{i}(r) = e^{ikr·α̂}, und sei r = R(cos φ, sin φ), α̂ der Einheitsvektor in Richtung α. Dann gilt für die Fourierkoeffizienten c_{n} der Fourierreihe ∑_{n=−∞}^{∞} c_{n}e^{inφ}:

c_{n} = \frac{1}{2π} ∫_{0}^{2π}e^{ikr·α̂}e^{inφ}dφ

= \frac{1}{2π} ∫_{0}^{2π}e^{ikR cos(φ−α)}e^{inφ}dφ

= \frac{1}{2π} ∫_{0}^{2π}e^{ikR cos(φ)}e^{in(φ−α)}

= \frac{1}{2π} ∫_{0}^{2π}e^{i(kR cos(φ)+nφ)}dφe^{−inα}

= J_{n}(kR)\frac{1}{2π} e^{−inα}  (6.1)

Für die Koeffizienten a_{n} in der obigen Reihe ergibt sich damit:

a_{0} = \frac{1}{2π}

a_{2n−1} = −\frac{1}{π} sin(nα)

a_{2n} = \frac{1}{π} cos(nα)  (6.2)

6.2 Auswertung der Integrale

Es sind Integrale der Form

∫_{∂D} u(x)ds(x)

auszuwerten. Da D ein Sterngebiet ist, wird der Rand von D über den Winkel parametrisiert. Sei also r(φ) der zum Winkel φ gehörige Abstand der Kurve vom Nullpunkt. Dann gilt:

∫_{∂D} u(x)ds(x) = ∫_{0}^{2π}u( \begin{matrix}r(φ) cos(φ) \\ r(φ) sin(φ) \end{matrix} ) ∥ \begin{matrix}r′(φ) cos(φ) − r(φ) sin(φ) \\ r′(φ) sin(φ) − r(φ) cos(φ) \end{matrix} ∥ dφ

= ∫_{0}^{2π}u( \begin{matrix}r(φ) cos(φ) \\ r(φ) sin(φ) \end{matrix} )\sqrt{r′(φ)^{2} + r(φ)^{2})}dφ  (6.3)

6.3 Optimierung

Um die Optimierungsaufgabe zu lösen, wurde die double-dogleg- Strategie gewählt. Sie ist eine Kombination aus gedämpftem Newton- und steilster Abstieg-Verfahren. Für die Implementation wurde das Verfahren dann leicht geändert, da die z.B. in [DS80] beschriebene Version nur mit Funktionen arbeitet, deren Hessematrix in einem gegebenen Gebiet positiv definit ist. Außerdem wurden die Hessematrizen durch grobe Näherungen ersetzt. Zur Beschreibung benötige ich zunächst einige Definitionen.

Definition 6.3.1 Sei f : IR^{n} ↦ IR zweimal stetig differenzierbar. Sei x_{0} ∈ IR^{n} fest. Falls die Hessematrizen ∇^{2}f(x_{k}) positiv definit sind, so definiert

∇^{2}f(x_{k})s_{k}^{N} = −∇^{2}f(x_{k})
x_{(k+1)} = x_{k} + s_{k}^{N}  (6.4)

eine Zahlenfolge, das Newtonsche Minimierungsverfahren. Falls die Folge gegen einen Wert x ∈ IR^{n} konvergiert und ∇^{2}f(x) positiv definit ist, so ist x ein lokales Minimum von f.

Beweis: Die Definition entspricht dem normalen Newton-Verfahren zum Finden von Nullstellen, angewandt auf die Funktion ∇f.

Zur Idee: f wird lokal durch eine Approximation zweiten Grades, nämlich

m_{c}(x_{n} + p) = f(x_{n}) + ∇f(x_{n})^{T} p + \frac{1}{2}p^{T} Hp,  (6.5)

beschrieben (H die Hessematrix). Das Minimum dieser Funktion wird dann gerade an der Stelle x_{n+1} angenommen.

Da die gegebenen Hessematrizen oft nicht positiv definit sind, wird meist eine modifizierte Version von 6.4 angewandt:

(∇^{2}f(x_{k}) + λ_{k}I)s_{k}^{N} = −∇^{2}f(x_{k})  (6.6)
x_{(k+1)} = x_{k} + s_{k}^{N}  (6.7)

Hierbei sei λ_{k} klein, d.h. λ_{k} >∼ min{λ : ∇^{2}f(x_{k})+λ_{k}I positiv definit)}. Für dieses Verfahren kann allerdings nur noch garantiert werden, daß es (wenn überhaupt) gegen eine Nullstelle x_{N} von ∇f konvergiert, dies muß nicht unbedingt ein lokales Minimum sein. x_{M} ist dann vielmehr ein Minimum der Funktion f + \frac{λ}{2}∥x − x_{n}∥^{2}, λ(x_{k}) = λ_{k}.

In diesem Fall lautet die Approximationsfunktion

m_{c}(x_{n} + p) = f(x_{n}) + ∇f(x_{n})^{T} p + \frac{1}{2}p^{T} H_{n}p,  (6.8)

(H_{n} = ∇^{2}f(x_{k}) + λ_{k}I).

Das Verfahren hat alle Eigenschaften des Newton-Verfahrens, d.h. in der Nähe eines Minimums konvergiert das Verfahren quadratisch (bei genügend glatter Zielfunktion). Leider erbt es auch die schlechte Eigenschaft, nicht global konvergent zu sein. In der Praxis schlägt das Newton-Verfahren oft Schritte vor, die sich weit vom Minimum entfernen.

6.4 Trust-region

Wie oben erwähnt, sind die im Newton-Verfahren gewählten Schritte oft sehr groß. Dadurch verläßt das Newton-Verfahren sehr schnell den Bereich, in dem die Approximation f = m_{c} (siehe 6.8) noch gültig ist. Bei den vorliegenden Problemen war dieses Phänomen sehr deutlich. Es wurde daher ein modifiziertes Newton-Verfahren benutzt.

Idee des Verfahrens: Es sei ein Minimierungsproblem mit der lokalen quadratischen Approximationsfunktion m_{c} am Punkt x_{0} gegeben. Es sei bekannt, daß m_{c} die zu minimierende Funktion f_{min} in einem Kreis um x_{0} mit Radius R in einem gewissen Sinne gut approximiert. Wähle dann als nächste Näherung für das Minimum das Minimum von m_{c} innerhalb des Kreises.

Aufbauend hierauf wurde das double-dogleg-Verfahren implementiert. Das Programm berechnet zuerst den Newtonschritt x_{k}^{N} und den Cauchypunkt x_{n}^{P}. x_{n}^{P} ist der nächste Schritt des Newtonschen Minimierungsverfahrens, angewandt auf die (eindimensionale) Funktion

f(x_{k} + λν_{k}), λ ∈ IR.

Hierbei ist ν_{k} die Richtung des steilsten Abstiegs.

Es werden nun drei Fälle unterschieden:

1. Der Newtonschritt liegt innerhalb der Trust-region. Dann wird für x_{n+1} der Newtonschritt gewählt.

2. Der Cauchypunkt liegt außerhalb der Trust-region. Dann wird für x_{n+1} der Schnitt der Trust-region mit der Richtung des steilsten Abstiegs gewählt.

3. Der Cauchypunkt liegt innerhalb, der Newtonschritt außerhalb der Trust-region. Dann wird für x_{n+1} der Schnitt der Trust-region mit der Verbindungsgeraden von Cauchypunkt und Newtonschritt gewählt.

Das Verfahren wurde nach [DS80] (leicht modifiziert) implementiert. Die Zeichnungen zeigen am besten, wie das Verfahren funktioniert.

Bemerkung: Falls der Newton-Schritt innerhalb des Kreises (der Trust-Region) liegt, fallen der double-dogleg-step und das Newtonverfahren zusammen. In der Nähe des Minimums bleibt also die quadratische Konvergenz erhalten.

6.5 Finite Differenzen zur Berechnung der Hessematrix

Da die zu minimierende Funktion in diesem Verfahren zur Lösung des inversen Helmholtz-Problems sehr schwer auszuwerten ist, wird eine Approximation der Hessematrix für die Optimierungsverfahren benutzt:

Zu minimieren sei die Funktion m : IR^{p} ↦ IR, m(x) = \frac{1}{2}|f(x) − d|^{2}, f : IR^{n} ↦ C^{n}, d ∈ C^{n}. Dann gilt:

\frac{∂m}{∂x_{j}} (x) = Re (∑_{i=0}^{p} ¯ \frac{∂F_{i}}{∂x_{j}} (x)(F_{i}(x) − d_{i}))  (6.9)

und damit für ∇m:

∇m(x) = Re (F′(x)^{∗}(F (x) − d))  (6.10)

Für die Hessematrix ergibt sich:

\frac{∂m}{∂x_{l}∂x_{j}} = Re (∑_{i=0}^{p} (\frac{∂F̄_{i}}{∂x_{l}∂x_{j}} (x)(F_{i} − d_{i}) + \frac{∂F̄_{i}}{∂x_{j}} (x)\frac{∂F_{i}}{∂x_{l}} (x)))  (6.11)

Da der erste Teil der Formel nur sehr langsam auszuwerten ist (die Berechnung einer normalen Approximation für die zweiten Ableitungen der Funktion F benötigt O(n^{2}) Auswertungen der Funktion), wird dieser in der Implementation weggelassen. Es gilt allerdings, daß in der Nähe des absoluten Minimums das Residuum F_{i}(x)−d_{i} klein ist und dort die Näherung und die Hessematrix (fast) übereinstimmen. Es sollte daher ein Minimierungsalgorithmus benutzt werden, der bei großer Entfernung vom Extrempunkt nicht einen Newtonschritt, sondern einen Schritt durchführt, der die zweiten Ableitungen nicht benötigt (Gradientenverfahren).

Der benutzte Algorithmus zur Berechnung der Hessematrix lautet damit:

H̃_{l,j} = Re (∑_{i=0}^{p} \frac{∂F̄_{i}}{∂x_{j}} (x)\frac{∂F_{i}}{∂x_{l}} (x))  (6.12)

oder

H̃ = Re(F^{′∗}F′)  (6.13)

Bemerkung 6.5.1 Die Näherungsmatrix H̃ ist immer positiv semidefinit.

Beweis:

Sei x ∈ IR^{p}. Dann ist

0 ≤ (F′x)^{∗}F′x = x^{t}(F^{′∗}F)x = Re (x^{t}(F^{′∗}F)x) = x^{t} Re (F^{′∗}F′) x (6.14)

Da H̃ positiv semidefinit ist, muß bei der Berechnung der Cholesky-Zerlegung von H̃ nur ein kleines Vielfaches der Einheitsmatrix addiert werden um die Existenz der Zerlegung zu sichern.

In der Praxis zeigte sich jedoch, daß H̃ sogar positiv definit ist, mit kleiner Kondition, und nicht sehr kleinen Eigenwerten.

6.6 Benutzte Rechner

Die Programme wurden auf SUN-Workstations der Versionen SUN-3/60 mit Prozessor 68020 und mathematischem Koprozessor 68881, SUN-4 mit SPARC-Prozessor und mathematischem Koprozessor sowie auf einem Großrechner IBM 3270 getestet.

Die mathematischen Ergebnisse unterschieden sich (im stabilen Bereich), wie nicht anders zu erwarten, nicht voneinander.

6.7 Parallele Verarbeitung

Aufgrund der hohen Rechenleistung, die fuer die Berechnung der Hessematrix des Problems benötigt wird, bietet sich eine Parallelisierung des Algorithmus auf einer Rechenanlage an.

Hierzu gibt es zwei Konzepte:

a) Bearbeitung des Problems auf einer Anlage mit Vektorprozessor, z.B. CYBER 205. Bei dieser Art der Parallelisierung können gleichartige Rechnungen mit sehr hoher Geschwindigkeit durchgeführt werden. Der Prozessor des Computers ist in der Lage, nicht nur skalare Größen, sondern direkt Vektoren zu verarbeiten. Typische Beispiele sind hier die Bildung eines Skalarprodukts oder die Addition zweier Vektoren, die jeweils mit maximalem Durchsatz gerechnet werden können.

b) Bearbeitung des Problems auf einer Multiprozessoranlage. Bei dieser Art der Parallelisierung wird die Gesamtaufgabe in viele Einzelaufgaben aufgeteilt und jeweils einem Prozessor zur Bearbeitung übergeben. Eine Kommunikation zwischen Steuerprozessor und den einzelnen ”Slave”-Prozessoren findet erst wieder statt wenn von diesen signalisiert wird, daß die Ergebnisse zum Abholen bereitstehen.

Analyse des Problems mit Vorwärtsmethode I:

Wie schon in Kapitel 4 gezeigt, besteht das Hauptproblem in der Berechnung der Matrizen Q_{e} und Q_{i} bei der Lösung des Vorwärtsproblems. Zur Berechnung des zugehörigen Integrals muß in jedem Interpolationsschritt eine Hankelfunktion berechnet werden. Diese Berechnung läßt sich mit Rechnern des Typs a) nur schlecht parallelisieren, da der Berechnungsvorgang nicht ”linear” ist, d.h. der Berechnungsvorgang selbst (und nicht nur die Werte, die in einen Algorithmus eingesetzt werden) hängt stark vom Argument der Hankelfunktion ab.

Hierzu ein Beispiel: Die Auswertung der Sinusfunktion läßt sich einfach parallelisieren, denn (mit Einschränkungen) sie läßt sich durch ein Polynom approximieren. Die Besselfunktion ist schlecht zu parallelisieren, denn bei ihrer Berechnung tritt eine Rekursion auf, deren Länge vom Wert des Arguments abhängig ist, d.h. die Art der einzelnen Rechenvorgänge hängt vom Wert des Arguments der Besselfunktion ab und eine gleichförmige Berechnung der Werte ist nicht möglich.

Somit lassen sich mit den Typ a) - Rechnern nur nebengeordnete Probleme wie das Invertieren der T-Matrix beschleunigen.

Andererseits sind die Rechner vom Typ b) geradezu auf das Problem zugeschnitten. Bei der Berechnung der Hessematrix ist für jedes Matrixelement jeweils ein Vorwärtsproblem zu lösen, daß von den anderen völlig unabhängig ist. Das legt nahe, jede dieser Rechnungen einem Prozessor zu übertragen und anschliessend nur die Ergebnisse wieder abzuholen. Bei n^{2} Prozessoren (n^{2} die Anzahl der Matrixelemente) wird die Komplexität somit von O(n^{2} ∗ V) auf O(V) herabgesetzt, O(V) die Komplexität des Vorwärtsproblems. Die Verwaltung und Verteilung der Aufgaben bereitet hier keine zusätzlichen Probleme, denn die Eingangsdaten bzw. Ergebnisse sind vom Umfang her sehr gering (als Parameter wird nur die Form des Sterngebiets benötigt, Ergebnis sind die berechneten Werte der Fernfeldamplitude).

Da kein solcher Rechner zur Verfügung stand, wurde eine Simulation auf einem ETHERNET-Netz implementiert. Hierzu wurde zunächst eine unterstützende Bibliothek geschaffen und dann das vorhandene Programm mit dieser Bibliothek auf verteilte Rechnungen angepaßt.

Konzept der Parallel-Simulation:

Die Bibliothek geht von folgender Architektur aus: In einem Netz von Unixrechnern gibt es auf einem Fileserver ein Verzeichnis, daß direkt von allen Rechnern erreicht werden kann. Über dieses Verzeichnis kann also jedes auf einem Rechner laufende Programm mit jedem auf einem anderen Rechner laufenden Programm in Verbindung treten.

Auf jeder Maschine, die für Rechnungen zur Verfügung steht, wird ein Hintergrundprogramm gestartet, das darauf wartet, daß ihm Eingangsdaten zur Verfügung gestellt werden. Sobald diese eingegangen sind, wird mit diesen Daten das Problem gelöst, die Ergebnisse werden dem Steuerprogramm zur Verfügung gestellt und gleichzeitig gemeldet, daß die Einheit wieder frei für neue Rechnungen ist.

Der Benutzer hat damit zwei Programme zu schreiben, einmal das Steuerprogramm, das die Aufgaben verteilt, und ein Slaveprogramm, das die Berechnungen auf den Daten, die ihm vom Steuerprogramm überreicht werden, ausführt.

Ablauf des Slaveprogramms: (Da die Konventionen in Fortran und C teilweise voneinander abweichen, sind die Eigenheiten für eine bestimmte Sprache jeweils gekennzeichnet.)

1. Initialisierung der Bibliotheksroutinen mit InitTasker().

2. Warten mit WaitTask(), bis eine zu rechnende Aufgabe vorliegt.

3. Einlesen und Rechnen der Aufgabe, dieser Teil (und nur dieser) ist natürlich von der Aufgabe abhängig.

(C) Der Name des Files mit den Eingabedaten steht in TaskerInputFile . (FORTRAN) Der Name des Files ist Ergebnis der Funktion (CHA-RACTER ∗ 80) TaskerInputFile.

4. Schreiben der Ergebnisse in das Ausgabefile.

(C) Der Name des Files mit den Ausgabedaten steht in der externen Variablen TaskerOutputFile . (FORTRAN) Der Name des Files ist Ergebnis der Funktion (CHA-RACTER ∗ 80) TaskerOutputFile.

5. Übergabe der Ergebnisse an das Steuerprogramm durch FinishTask().

6. Wiederholung ab Schritt 2.

Ablauf des Steuerprogramms:

1. Die Namen der Rechner, auf denen das Slaveprogramm läuft und die auf die Übergabe von Daten warten, werden mit addmachine(MachineName) bekanntgemacht. Beispiel: call addmachine(’polya’) call addmachine(’calderon’)

(C) Vereinbarung aller zur Verfügung stehenden Rechner auf einmal durch den Aufruf InitTasker(”calderon”,”polya”,0).

Die Namen der Files, in die das Steuerprogramm die Daten für die Slaveprogramme schreibt bzw. aus denen es die Resultate zurückholt, werden mit SetMainInputFile bzw. SetMainOutputFile gesetzt.

Für jede zu bearbeitende Aufgabe:

2. Schreibe die aktuellen Parameter der Aufgabe in das InputFile.

3. Verteile die Aufgabe auf einen freien Rechner mit SpreadTask(task), task ist die Tasknummer unter der später die Ergebnisse abgeholt werden können. Wenn kein Rechner frei ist: Warte bis ein Rechner frei wird, sichere die erhaltenen Ergebnisse und übergib diesem die neue Aufgabe. (Dies wird von den Bibliotheksroutinen erledigt.)

4. Wiederhole ab Schritt 2 bis alle parallel bearbeitbaren Aufgaben verteilt sind. Für jede bearbeitete Aufgabe:

5. Fordere die Ergebnisse einer Task an durch FetchTaskResults(task) und lies sie aus dem angegebenen OutputFile.

Diese Festlegung des Ablaufs hat zwei Vorteile:

1. Der Aufwand ein bestehendes Programm zu ändern, damit es auch mit dieser parallelen Architektur läuft, ist sehr gering.

2. Auch der in der Bibliothek versteckte Verwaltungsaufwand ist sehr gering solange die übertragenen Datenmengen relativ klein sind gemessen an dem Rechenaufwand, der mit ihnen getrieben wird.

Dank dieser Vorteile konnte die natürliche Grenze für die Beschleunigung des Algorithmus fast erreicht werden, bei einer Verteilung auf n Prozessoren ergibt sich eine Gesamtzeit von fast T/n, T die Zeit für die Bearbeitung auf einem Prozessor, für die Berechnung der Hessematrix.

Hierbei soll natürlich nicht verschwiegen werden, daß bei großen Datenmengen oder insbesondere kurzen Zeiten für die Berechnung einer einzelnen Aufgabe das vorgestellte Verfahren absolut ineffizient ist. In diesem Fall benötigt die Kommunikation zwischen den Programmen viel mehr Zeit als durch die Verteilung der Aufgaben eingespart wird.

Hier noch einmal schematisch die drei besprochenen Implementationen bzw.

Rechnerarchitekturen:

Bild 1: Parallelverarbeitung auf einem Vektorrechner. Alle Vorarbeiten werden auf dem skalaren Prozessor durchgeführt. Erst wenn viele gleiche Operationen (Addition, Multiplikation) an einem Stück durchgeführt werden, wird der Vektorprozessor aktiv. Nachteil: Das Programm muß speziell für die Vektorarchitektur umgeschrieben werden, viele Programme sind gar nicht vektorisierbar, obwohl sie unabhängig voneinander berechenbare Teile enthalten. Vorteil: Der Transfer vom Speicher in den Vektorprozessor ist sehr schnell, die Kommunikation zwischen Skalarprozessor und Vektorprozessor wird direkt von der Hardware gesteuert.

Bild 2: Parallelverarbeitung auf einem Multiprozessorsystem. Alle Vorarbeiten werden auf einem Hauptprozessor durchgeführt. Sobald unabhängige Programmteile auftreten, können diese erkannt und parallel berechnet werden.

Bild 3: Pseudo-Parallelverarbeitung in separaten Programmen. Wie (2), aber die Programme müssen für die Verarbeitung speziell aufbereitet werden. Nachteil: Die Kommunikation zwischen den Prozessen wird nicht vom Betriebssystem, sondern von den Programmen selbst gesteuert. Das macht diese langsam und für große auszutauschende Datenmengen ungeeignet.

6.8 Komplexität des Problems, Laufzeitverhalten

In 6.1 und 6.2 sind zwei Laufzeitprofile für die Verfahren I und III dargestellt. Es ergibt sich:

Bei Algorithmus I wird der größte Teil der Zeit zum Berechnen der speziellen Funktionen benutzt (Bessel-, Hankel-, trigonometrische Funktionen). Die Lösung der Gleichungssystem fällt dann nicht mehr ins Gewicht.

Bei Algorithmus III wird der größte Teil der Zeit verbracht in der kQl-Routine, die die QR-Zerlegung berechnet und das entsprechende Gleichungssystem löst. Algorithmus III wäre deshalb auf einem Vektorrechner eher parallelisierbar als Algorithmus I, denn für diese Standardaufgabe stehen fertige parallelisierte Prozeduren bereits zur Verfügung.

Numerische Ergebnisse

Als wesentliches Problem des Verfahrens stellte sich schnell die hohe Komplexität des Vorwärts-Algorithmus heraus. Es wurde deshalb besonderer Wert auf eine korrekte und effiziente Implementierung gelegt. Die drei vorgestellten Algorithmen werden zunächst von den Genauigkeiten und Zeiten her verglichen. Hierbei stellt sich eine deutliche Überlegenheit des Ansatzes von Visscher gegenüber dem klassischen T-Matrix Ansatz von Waterman heraus.

7.1 Ergebnisse der Vorwärtsalgorithmen

Alle Ergebnisse in diesem Kapitel beziehen sich auf Dirichlet-Randbedin-gungen. Zu den Zeichnungen: Im Zentrum des Bildes befindet sich jeweils das Objekt. In Abhängigkeit vom Winkel sind über einem Kreis Real- und Imaginärteil der resultierenden Welle auf dem Rand des Objekts abgetragen, d.h. wenn die Randbedingung optimal erfüllt ist, fallen der Kreis und die Wellen zusammen. Bei einigen Bildern ist zusätzlich die einfallende Welle abgetragen. Die einzelnen Diagramme sind zeilenweise durchnumeriert zu denken. Als relative Frequenz wird die Größe “Frequenz · Durchmesser des Objekts” bezeichnet.

Nach den theoretischen Ergebnissen müssen die Ergebnisse des Objekts “Kreis” bei den drei Verfahren übereinstimmen (Bilder 1.1 und 1.2).

In 1.1 wird das Verhalten am Kreis bei kleiner Frequenz demonstriert. Die Bilder 1 und 4, 2 und 3, sowie 5 und 6 wurden mit einer unterschiedlichen Anzahl von Entwicklungsfunktionen jeweils mit den Algorithmen I und III berechnet. Die Ergebnisse sind exakt gleich.

Für alle anderen Bilder gilt:

Die obere Reihe des Diagramms zeigt Bilder der ausfallenden Welle, berechnet mit den Verfahren II, I und III, resp. Die Anzahl der Entwicklungsfunktionen war hier hoch, es wurden viele Interpolationsschritte gewählt. Die untere Reihe zeigt jeweils das Verhalten bei wenigen Entwicklungsfunktionen und wenigen Interpolationsschritten.

Es zeigt sich: Verfahren I und II liefern oft numerische Instabilitäten, Verfahren III tut dies nicht. Wenn I nicht instabil wird, sind die Ergebnisse von I und III gleich. Dies deckt sich mit den theoretischen Resultaten.

7.2 Ergebnisse des Algorithmus bei gestörten und ungestörten Daten

Die Diagrammseiten 2.x demonstrieren Ergebnisse des Algorithmus bei einer kleinen relativen Frequenz (≤ 10).

7.3 Wahl eines geeigneten Regularisierungsparameters

Für das Programm wurden, wie oben beschrieben, zwei Methoden der Regularisierung angewandt. Die praktischen Ergebnisse von Methode II waren besser als die von Methode I. Außerdem war die Regularisierung II typischerweise in einem Viertel der Zeit der Regularisierung I durchführbar. Hier einige numerische Ergebnisse, es wird jeweils der L^{2}-Fehler ∗15 gegen die Anzahl der Iterationen abgetragen:

Pinocchio-Shape, relative Frequenz 5, 8 ∗ 8 gleichverteilte Meßdaten, Regularisierung 1:

Exakte Daten: Bei exakten Daten darf die Anzahl der berücksichtigten Funktionen zur Approximation des Randes hoch gewählt werden. Im Rahmen der Versuche war kein Ansteigen des Fehlers zu beobachten. Der Durchlauf wurde nach der 23. Iteration abgebrochen, weil der Gradient auch bei hoher Funktionenzahl fast 0 war und deshalb das Minimierungsverfahren keine Schritte mehr durchführte.

Schwach gestörte Daten, Standardabweichung 1 % des typischen Werts: Bei schwach gestörten Daten ist ein leichter Anstieg des Fehlers zum Ende der Iteration zu bemerken. Dieser wirkte sich jedoch in den Bildern nicht stark aus. Interessant ist, daß die gestörten Daten fast dasselbe Minimum wie die exakten Daten erreichen.

Stark gestörte Daten, Standardabweichung 10 % des typischen Werts: Bei stark gestörten Daten wird nach kurzer Zeit ein Minimum erreicht, anschließend steigt der Fehler stark an.

Hier einige typische Beispiele, wie sich die Wahl einer zu großen Funktionenschar auswirkt. Oben stehen jeweils alle Zwischenobjekte eines Iterationsdurchlaufs, unten ein Vergleich von Original und letzter Approximation vor dem Abbruch des Programms.

7.4 Genauigkeit der Approximation der Hessematrix

In den Diagrammen sind die relativen Fehler der Hessematrix bei genauer Ermittlung durch finite Differenzen, und bei Approximation wie im Teil Optimierung beschrieben, gegenübergestellt. 3.1 und 3.2 sind die Ergebnisse bei niedriger Frequenz (relative Frequenz 7) ohne und mit gestörten Daten. 3.3 zeigt das Ergebnis bei ungestörten Daten und der relativen Frequenz 15.

Obwohl die Hessematrix bei hoher Frequenz eine große Abweichung zeigt, ist der Einfluß auf das Optimierungsverfahren nur gering. Die verwandten Newtonrichtungen unterscheiden sich nur wenig, die erhaltenen Bilder sind gleich. Es gilt daher, daß das verwendete Optimierungsverfahren für hohe Frequenzen in jedem Fall in Nebenminima endet und dieses Verhalten nicht der schlechten Approximation zugeschrieben werden kann.

7.5 Verhalten bei reduzierter Anzahl von Einstrahlungsrichtungen

Bei reduzierter Zahl von Einstrahlungsrichtungen sind zwei Phänomene zu beobachten:

1. Die Optimierung dauert länger, es entstehen viele Zwischenschritte.

2. Die Objekte haben schlecht bestimmte Kanten, falls die Einstrahlrichtungen nicht gleichmäßig um sie verteilt sind.

Hier ein Bild bei Einstrahlung von unten. Die obere Kante ist schlecht bestimmt.

Eine Beobachtung, die sowohl von Kristensson/Vogel, als auch von Boström gemacht wurde, läßt sich am Pinocchio-Shape beobachten. Wenn die Einstrahlrichtung der einzigen einfallenden Welle auf einer Symmetrieachse des Objekts liegt, ist dies weiterhin rekonstruierbar (Bild 4.1). Ist dies nicht der Fall, wird die Berechnung sehr instabil (Bild 4.2). Isolinienplots der einfallenden Wellen geben hier die Begründung: Es entstehen mehrere Wellenfronten, die nur durch eine größere Anzahl von Daten eindeutig aufgelöst werden.

Zu den Isolinienplots 4.3 und 4.4: 4.3 zeigt das Verhalten von einfallender, ausfallender und Gesamtwelle am Objekt “Kreis”.

4.4 zeigt das Verhalten beim Objekt “Pinocchio”. In der oberen Reihe fällt die Welle parallel zur Symmetrieachse, in der unteren senkrecht zur Symmetrieachse ein. Die Bilder ganz rechts zeigen jeweils das Verhalten der ausfallenden Welle.

7.6 Verhalten bei hohen Frequenzen

Das Verfahren ist für höhere relative Frequenzen als bei Kristensson/Vogel vorgeschlagen nicht geeignet. Das Minimierungsverfahren landet sehr schnell in Nebenminima, die noch sehr weit vom Objekt entfernt liegen. Selbst bei sehr starker Regularisierung ist keine Annäherung an das Objekt erreichbar. Eine Konvergenz der Näherung gegen dem Umriß des Objekts wird nur durch Zufall erreicht.

Dies ist hier an einigen Beispielen gezeigt:

(5.1): Versuch der Rekonstruktion eines Pinocchio-Shapes. Der Algorithmus erkennt zwar, daß eine Störung auf der rechten Seite vorliegt, wählt die Ausbuchtung aber nach innen statt nach außen. Bei Erhöhung der Anzahl der Funktionen divergiert das Verfahren.

(5.2): Hier wurde eine Rekonstruktion eines Kreises aus einen Funktionenraum, in dem sich nur Kreise befanden, durchgerechnet. Selbst in diesem einfachsten Fall war es nicht möglich bei weit entferntem Originalobjekt eine korrekte Rekonstruktion durchzuführen.

(5.3): Rekonstruktion eines Peanut-Shape. Da das Testobjekt äußerst glatt ist, ist bei guter Anfangsnäherung hier immerhin eine einigermaßen gute Approximation erreichbar.

7.7 Überprüfung der Korrektheit der Programme zur Lösung des Vorwärtsproblems

Die Programme sind so konzipiert, daß die Methoden I,II und III in ihren wesentlichen Teilen unabhängig voneinander arbeiten. Da alle drei bei genügend groß gewählten Parametern dieselben Ergebnisse liefern, kann davon ausgegangen werden, daß die Implementationen korrekt sind.

Außerdem wurden zwei Sätze aus der Physik zur Überprüfung der Ergebnisse herangezogen:

Bei Drehen des Versuchsaufbaus um 180 Grad werden dieselben Ergebnisse erzielt, d.h. bei Einstrahlungsrichtung α und Meßrichtung β erhält man dieselben Werte wie bei Einstrahlungsrichtung α+π und Meßrichtung β +π.

Nach [Ram86]: Für die Fernfeldamplitude F gilt

Im F (α, k, α) = \frac{k}{4π} ∫_{S^{2}} |F (ω, k, α)|^{2}dω  (7.1)

Die Überprüfung dieser Gleichung wurde von Visscher zur Programmverifikation vorgeschlagen. Die vorhandene Implementation bestand diesen Test.

Zusammenfassung und Ausblick

8.1 Erweiterung  des  Verfahrens  auf

Nicht-Sterngebiete

Die Einschränkung des Verfahrens auf Sterngebiete ist natürlich nicht durch das Verfahren, sondern nur durch die Implementation vorgeschrieben. Zur Vereinfachung wurde immer von einem über einen Winkel ϕ parametrisierbaren Rand ausgegangen. Natürlich ist dieser Rand ebenso über dem Intervall [0, 2π] parametrisierbar durch den Ansatz:

x(ϕ) = ∑_{n=0}^{∞} (a_{n}cos(nϕ) + b_{n}sin(nϕ))  (8.1)
y(ϕ) = ∑_{n=0}^{∞} (c_{n}cos(nϕ) + d_{n}sin(nϕ))  (8.2)

Dabei wird natürlich die Komplexität des Verfahrens quadriert.

8.2 Anwendung auf das inverse Streuproblem bei mehreren Körpern

In mehreren Artikeln (z.B. Varadan/Varadan) wird auch auf das Streuproblem bei mehreren Körpern eingegangen. Die Methode von Visscher liefert hier eine Lösung des Vorwärtsproblems.

Deshalb kann die vorgestellte Methode natürlich auch auf das inverse Streuproblem bei mehreren Körpern angewandt werden. (Eindeutigkeitssätze zum Beispiel bei (Varadan/Varadan)).

Dazu muß nur ein innerer Punkt der Körper bekannt sein. Jeder einzelne der Körper wird dann parametrisiert über den Raumwinkel. Dadurch ist das Problem zurückgeführt auf die Approximation der entsprechenden Radiusfunktionen r_{1}..r_{N}, N die Anzahl der Körper. Unter diesen Voraussetzungen kann das Problem, wie in dieser Arbeit beschrieben, gelöst werden.

8.3 Implementation des Algorithmus im IR^{3}

Die Komplexität des IR^{3}-Problems ist ungefähr das Quadrat der Komplexität des IR^{2}-Problems. Mit allen aufgezeigten Optimierungen, die in der Arbeit besprochen wurden, und den neuen Rechnerentwicklungen gilt daher m.E. der Satz von Kristensson/Vogel, vor einer Implementation des Verfahrens für den IR^{3} müssen neue Verfahren zur Lösung des Vorwärtsproblems entwickelt werden, nicht mehr. (Es sei denn, man bezeichnet die Methode von Visscher als neu, aber die war damals schon entwickelt (s. Varadan/Varadan)).

8.4 Zusammenfassung

• Das von Kristensson/Vogel vorgeschlagene Verfahren liefert unter geeigneten Bedingungen eine Lösung des inversen Problems. Das Verfahren ist nicht geeignet bei hohen Frequenzen oder wenig glatten Rändern. In diesen Fällen führt das Verfahren auf falsche Ergebnisse.

• Die MOOT-Methode von Visscher und die Methode der überbestimmten Interpolation sind (numerisch) äquivalent. Bei Kreisen sind die MOOT und die T-Matrix-Methode von Waterman ebenfalls äquivalent.

• Die (modifizierte) MOOT-Methode von Visscher zur Lösung des Vorwärtsproblems stellte sich als effizienter heraus als Waterman’s T-Matrix Methode. MOOT verkürzte den Aufwand zur Lösung des Problems erheblich. Die Methode der überbestimmten Interpolation lieferte eine Idee zur Implementation der MOOT. Die dadurch erzielten Ergebnisse waren noch besser als die der Direktimplementation der Methode von Visscher.

• Die Methode der Regularisierung der Lösung durch kleine Approximationsräume ist effizienter als die Tychonoff-Philips-Regularisierung und dieser vom Ergebnis her nur knapp unterlegen.

[DS80] [Vis80] [VO80] [KV86] [Vog85] [VV80] [CH53a] [CH53b] [Lei67]

[Col84] [CS83] [KRS83] [Lav70] [LP67] [New80] [Sle82] [Wat65] [Wat69]

[CK83] [KKMZ88] [Ram86] [KK83] [Nat85a] [Nat85b] [Zin89] [Pre86]

Literaturverzeichnis

[BS81] I. N. Bronstein and K. A. Semendjajew. Taschenbuch der Mathematik. B. G. Teubner Verlagsgesellschaft, 1981.

[CH53a] R. Courant and D. Hilbert. Methods of Mathematical Physics, volume 1. Interscience Publishers Ltd., 1953.

[CH53b] R. Courant and D. Hilbert. Methods of Mathematical Physics, volume 2. Interscience Publishers Ltd., 1953.

[CK83] D. Colton and R. Kress. Integral equation methods in scattering theory. John Wiley, 1983.

[Col84] D. Colton. The inverse scattering problem for time-harmonic acoustic waves. SIAM Review, 26, 1984.

[CS83] D. Colton and B. Sleeman. Uniqueness theorems for the inverse problem of acoustic scattering. IMA Journal of Applied Mathematics, 31, 1983.

[DS80] Dennis and Schnabel. Numerical methods for unconstrained optimization and nonlinear equations. Prentice-Hall, 1980.

[Els86] J. Elstrodt. Skript zur Vorlesung Funktionentheorie I/II 1985/86. Fachbereich Mathematik der WWU, 1986.

[KK83] R. E. Kleinman and R. Kress. On the condition number of integral equations in acoustics using modified fundamental solutions. IMA Journal of Applied Mathematics, 31, 1983.

[KKMZ88] A. Kirsch, R. Kress, P. Monk, and A. Zinn. Two methods for solving the inverse acoustic scattering problem. Inverse Problems, 4, 1988.

[KRS83] G. Kristensson, A. G. Ramm, and S. Strom. Convergence of the t-matrix approach in scattering theory. Journal of Mathematical Physics, 24, 1983.

[KV86] G. Kristensson and C. R. Vogel. Inverse problems for acoustic waves using the penalised likelihood method. Inverse Problems, 1986.

[Lav70] M. Lavrientiev. Multidimensional inverse problems for diffeential equations. Lecture notes in Mathematics, 167, 1970.

[Lei67] Rolf Leis. Vorlesungen über partielle Differentialgleichungen zweiter Ordnung. BI Hochschultaschenbücher, 1967.

[LP67] P. Lax and R. Phillips. Scattering theory. Academic Press, 1967.

[Nat85a] F. Natterer. Skript zur Vorlesung Praktische Mathematik I WS 1984/85. Fachbereich Mathematik der WWU, 1985.

[Nat85b] F. Natterer. Skript zur Vorlesung Praktische Mathematik II SS 1985. Fachbereich Mathematik der WWU, 1985.

[New80] R. Newton. Inverse scattering. Journal of Mathematical Physics, 21, 1980.

[Pre86] William H. Press. Numerical recipes. Cambridge Unviersity Press, 1986.

[Ram86] A. G. Ramm. Scattering by Obstacles. D. Reidel Publishing, 1986.

[Sle82] B. Sleeman. The inverse problem of acoustic scattering. IMA Journal of Applied Mathematics, 29, 1982.

[Vis80] William M. Visscher. A new way to calculate scattering of acoustic and elastic waves 1, theory illustrated for scalar waves. Journal of Applied Physics, 51(2), 1980.

[VO80] William M. Visscher and Jon L. Opsal. Theory of elastic wave scattering: Applications of the method of optimal truncation. Journal of Applied Physics, 51(2), 1980.

[Vog85] C. R. Vogel. Numerical solution of a non-linear ill-posed problem arising in inverse scattering. Inverse Problems, 1985.

[VV80] V. K. Varadan and V. V. Varadan. Acoustic, Electromagnetic, and Elastic Wave Scattering – Focus on the T-matrix Approach. Pergamon, 1980.

[Wat65] P. Waterman. Matrix formulation of electromagnetic scattering. Proceedings of the IEEE, 53, 1965.

[Wat69] P. Waterman. New formulation of acoustic scattering. Journal of the Acoustical Society of America, 45, 1969.

[Zin89] A. Zinn. On an optimisation method for the full- and the limited-aperture problem in inverse acoustic scattering for a sound-soft obstacle. Inverse Problems, 5(2), 1989.

Inhaltsverzeichnis

1 Einleitung 2

2 Das Streuproblem 6 2.1 Definitionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Greensche Formeln . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Darstellungssätze für Lösungen der Helmholtz-Gleichung … 9 2.4 Entwicklung von Lösungen der Helmholtzgleichung . . . . . . 14 2.5 Die Eindeutigkeit von Lösungen der Helmholtzgleichung … 19 2.6 Existenz von Lösungen - Ansatz über Integralgleichungen . . 20 2.7 Vollständigkeit von Funktionensystemen . . . . . . . . . . . . 23

3 Das inverse Streuproblem 26 3.1 Eindeutigkeit des inversen Problems . . . . . . . . . . . . . . 26 3.2 Schlechtgestelltheit des Inversen Streuproblems . . . . . . . . 30

4 Lösung des Streuproblems 31 4.1 Berechnung des reflektierten Feldes I . . . . . . . . . . . . . . 31 4.2 Berechnung des reflektierten Feldes II . . . . . . . . . . . . . 36

4.3 Berechnung des reflektierten Feldes III . . . . . . . . . . . . . 37

5 Lösung des inversen Streuproblems 43 5.1 Beschreibung des Verfahrens . . . . . . . . . . . . . . . . . . . 43 5.2 Regularisierung der Lösung . . . . . . . . . . . . . . . . . . . 44

6 Implementation 46 6.1 Entwicklung der einfallenden Welle bei Methode I . . . . . . 48 6.2 Auswertung der Integrale . . . . . . . . . . . . . . . . . . . . 49 6.3 Optimierung . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 6.4 Trust-region . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 6.5 Finite Differenzen zur Berechnung der Hessematrix . . . . . . 52 6.6 Benutzte Rechner . . . . . . . . . . . . . . . . . . . . . . . . . 53 6.7 Parallele Verarbeitung . . . . . . . . . . . . . . . . . . . . . . 53 6.8 Komplexität des Problems, Laufzeitverhalten . . . . . . . . . 60

7 Numerische Ergebnisse 62 7.1 Ergebnisse der Vorwärtsalgorithmen . . . . . . . . . . . . . . 62 7.2 Ergebnisse des Algorithmus bei gestörten und ungestörten

Daten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 7.3 Wahl eines geeigneten Regularisierungsparameters . . . . . . 63 7.4 Genauigkeit der Approximation der Hessematrix . . . . . . . 67 7.5 Verhalten bei reduzierter Anzahl von Einstrahlungsrichtungen 67 7.6 Verhalten bei hohen Frequenzen . . . . . . . . . . . . . . . . . 68 7.7 Überprüfung der Korrektheit der Programme zur Lösung des

Vorwärtsproblems . . . . . . . . . . . . . . . . . . . . . . . . 69

8 Zusammenfassung und Ausblick 71 8.1 Erweiterung des Verfahrens auf Nicht-Sterngebiete . . . . . . 71 8.2 Anwendung auf das inverse Streuproblem bei mehreren Körpern 72 8.3 Implementation des Algorithmus im IR^{3} . . . . . . . . . . . . 72 8.4 Zusammenfassung . . . . . . . . . . . . . . . . . . . . . . . . 72