Matthias Vögeler

Ein

Propagations–Rückpropagations–Algorithmus

zur Lösung des inversen Streuproblems für die Maxwell–Gleichungen

Angewandte Mathematik

Ein

Propagations–Rückpropagations–Algorithmus zur Lösung des inversen Streuproblems

für die Maxwell–Gleichungen

Inaugural–Dissertation zur Erlangung des Doktorgrades

der Naturwissenschaften im Fachbereich Mathematik und Informatik der Mathematisch–Naturwissenschaftlichen Fakultät der Westfälischen Wilhelms–Universität Münster

vorgelegt von

Matthias Vögeler

aus Beelen

Münster, 2001

Dekan:  Prof. Dr. W. Lange
Erster Gutachter:  Prof. Dr. F. Natterer
Zweiter Gutachter:  Prof. Dr. C. Cryer

Tag der mündlichen Prüfungen: Tag der Promotion:

Zusammenfassung

Diese Arbeit hat das nichtlineare inverse Streuproblem für die Maxwell–Gleichungen im Bereich der Mikrowellen–Tomographie zum Thema. Dazu stelle man sich ein in einem homogenen Medium eingebettetes Objekt vor, das von einer planaren zeitharmonischen elektromagnetischen Welle bestrahlt wird. Das inverse Problem besteht nun darin, aus Messungen des gestreuten elektrischen Feldes in einiger Entfernung vom Objekt den komplexwertigen Brechungsindex des Objektes zu rekonstruieren. Das Lösen des Anfangswertproblems für die Maxwell–Gleichungen führt zu einem nichtlinearen Operator vom Raum der „Brechungsindizes“ in den Raum der „Messdaten“. Wird die Fréchet–Ableitung dieses Operators benutzt und stehen Messungen des gestreuten elektrischen Feldes für verschiedene Einfallsrichtungen der planaren zeitharmonischen elektromagnetischen Welle zur Verfügung, dann kann das inverse Problem iterativ gelöst werden. Ein Iterationsschritt besteht dabei aus einer Propagation (Lösung des direkten Problems) und einer Rückpropagation (Lösung des „adjungierten“ linearisierten direkten Problems). Der verwendete Algorithmus hat sich bereits auf dem Gebiet der Ultraschall–Tomographie etabliert. Hier wird zum ersten Mal seine Anwendbarkeit auf dem Gebiet der Mikrowellen–Tomographie untersucht. Unter anderem musste eine äquivalente Formulierung für die Maxwell–Gleichungen gefunden werden, die bedingte Stabilität der Anfangswertaufgabe und die Fréchet–Differenzierbarkeit des erwähnten nichtlinearen Operators wurde gezeigt. Numerische Simulationen belegen die analytischen Resultate für den Algorithmus und demonstrieren seine praktische Brauchbarkeit.

Inhaltsverzeichnis

1 Einleitung 1 1.1 Übersicht . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Praktische Relevanz . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Problembeschreibung . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.3.1 Mathematische Formulierung des Streuproblems . . . . . . . 3 1.3.2 Beschreibung des inversen Streuproblems . . . . . . . . . . . 6

1.4 Literaturüberblick . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.4.1 Rekonstruktion der Form und Oberflächenimpedanz eines Ob-

jektes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4.2 Rekonstruktion der Form von allgemeineren Objekten …. 9 1.4.3 Rekonstruktion von im Boden eingelassenen zylindrischen Ob-

jekten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.4.4 Dreidimensionale Mikrowellen–Tomographie . . . . . . . . . . 12 1.4.5 Bemerkung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2 Das direkte Problem 14 2.1 Regularität des elektrischen Feldes E . . . . . . . . . . . . . . . . . . 15 2.2 Elliptizität . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3 Äquivalente Formulierung der Maxwell–Gleichungen . . . . . . . . 19 2.4 Das Anfangswertproblem für das elektrische Feld . . . . . . . . . . . 24

2.4.1 Stabilitätsanalyse . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.4.2 Die Helmholtz–Gleichung als Tiefpassfilter . . . . . . . . . 33 2.4.3 Zusammenfassung . . . . . . . . . . . . . . . . . . . . . . . . 35

3 Das inverse Problem 37 3.1 Linearisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2 Der PBP–Algorithmus . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4 Numerische Simulationen 49

4.1 Numerische Implementierung . . . . . . . . . . . . . . . . . . . . . . 50 4.2 Simulation 1: Rekonstruktion einer homogenen Kugel . . . . . . . . . 51

4.2.1 Verifikation des Vorwärtslösers . . . . . . . . . . . . . . . . . 53 4.2.2 Rekonstruktionsergebnisse . . . . . . . . . . . . . . . . . . . . 56

4.3 Simulation 2: Rekonstruktion eines Herzmodells . . . . . . . . . . . . 59

4.3.1 Konvergenz des Vorwärtslösers . . . . . . . . . . . . . . . . . 64 4.4 Simulation 3: Der inhomogene Zylinder . . . . . . . . . . . . . . . . 65 4.5 Bornsche Näherung . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.6 Approximation durch die Helmholtz–Gleichung . . . . . . . . . . . 69 4.7 Abschließende Bemerkung . . . . . . . . . . . . . . . . . . . . . . . . 74

Literaturverzeichnis 75

Kapitel 1

Einleitung

Zunächst wird eine Übersicht über den Aufbau der Arbeit gegeben:

1.1 Übersicht

Das erste Kapitel dient der Einführung in die Problemstellung und der Skizzierung des Verfahrens zum Lösen des inversen Streuproblems für die Maxwell– Gleichungen. Dabei werden sowohl die geometrischen Rahmenbedingungen als auch wesentliche Begriffe und Bezeichnungen eingeführt. Desweiteren wird ein Überblick zu anderen aktuellen Arbeiten auf dem Gebiet des inversen Streuproblems für die Maxwell–Gleichungen gegeben.

Als wichtiger Zwischenschritt zur Lösung des inversen Streuproblems befasst sich das zweite Kapitel mit dem direkten Streuproblem. Zum einen wird eine geeignetere äquivalente Formulierung für die Maxwell–Gleichungen hergeleitet. Zum anderen ist die bedingte Stabilität des Anfangswertproblems für die Maxwell–Gleichungen von großem Interesse.

Im dritten Kapitel wird der Propagations–Rückpropagations–Algorithmus zur Lösung des inversen Streuproblems hergeleitet. Dabei ist die Fréchet–Ableitung eines nichtlinearen Operators und ihre Adjungierte von großer Bedeutung.

Im vierten und letzten Kapitel werden drei numerische Simulationen vorgestellt. Analytische Ergebnisse und die Brauchbarkeit des Propagations–Rückpropagations– Algorithmus werden gezeigt. Desweiteren wird untersucht, wie gut die Maxwell– Gleichungen durch die skalare Helmholtz–Gleichung approximiert werden können.

1.2 Praktische Relevanz

Auf dem Gebiet der bildgebenden Verfahren in der Medizin ist die Mikrowellen– Tomographie zunehmend Gegenstand intensiver Forschung. Sie hat gegenüber der Röntgen–, CT– und NMR–Tomographie entscheidende Vorteile, da sie ohne Verwendung ionisierender Strahlung oder invasiven Kontrastmitteln auskommt. Die Mikrowellen–Tomographie ist in der Lage, dielektrische Eigenschaften von biologischem Gewebe zu bestimmen. Dabei ist entscheidend, dass die dielektrischen Eigenschaften lebenden Gewebes von ihrem physiologischen Zustand abhängen. So ist bereits in [Semenov et al;1996] gezeigt worden, dass Unterschiede zwischen gesundem Herzgewebe und Infarkt gefährdetem Herzgewebe bzgl. ihrer dielektrischen Eigenschaften bestehen. Auch konnte nachgewiesen werden, dass von Krebs befallenes Brustgewebe im Bereich der Mikrowellen–Strahlung signifikante Änderungen der dielektrischen Leitfähigkeit gegenüber gesundem Gewebe zeigt [Hagnes et al;1998]. In [Joines et al;1980] und [Schepps,Foster;1980] wurde aufgeführt, dass es an Hand der Unterschiede dielektrischer Leitfähigkeit im Gewebe möglich ist verschiedene Tumorarten zu bestimmen.

1.3 Problembeschreibung

Diese Arbeit behandelt das nichtlineare inverse Streuproblem für die Maxwell– Gleichungen. Ziel ist die Rekonstruktion einer bestimmten physikalischen Eigenschaft eines Objektes aus noch näher zu definierenden Messdaten. Diese Eigenschaft wird mathematisch durch die komplexwertige Ortsfunktion n(x) beschrieben. Das Objekt befindet sich in einer Kugel D mit dem Radius ρ, was mathematisch durch die Forderung supp(n(x)−1) ⊂ D ausgedrückt wird. Die Funktion n(x) wird sich als Brechungsindex herausstellen. Üblicherweise wird durch n = 1 auf das Referenzmedium Vakuum verwiesen. Für alle anderen Brechungsindizes würde damit Re n > 1 gelten. Es ist natürlich möglich für ein anderes Referenzmedium n = 1 zu fordern, was für andere Brechungsindizes Re n > 0 bedeuten würde (Begrenzung durch die Lichtgeschwindigkeit). Weiter befindet sich die Kugel D in einem Würfel Q_{θ} der Seitenlänge 2ρ. Die Seiten des Würfels werden mit Γ_{θ}, Γ_{θ}^{−} und Γ_{θ}^{+}, wie in Figur 1.1 illustriert, beschrieben. Dabei beschreibt der Einheitsvektor θ die Orientierung des Würfels. Die Seiten Γ_{θ}^{−} und Γ_{θ}^{+} stehen senkrecht zum Vektor θ, wobei Γ_{θ}^{+} von Γ_{θ}^{−} aus gesehen in Richtung θ liegt.

Für die Rekonstruktion des Objektes in D, also der Funktion n(x), werden Messdaten benötigt. Dazu wird der Würfel Q_{θ} mit einem ebenen zeitharmonischen 1.3. PROBLEMBESCHREIBUNG 3

[figure: θ; θ_{1}^{⊥}; θ_{2}^{⊥}; Γ_{θ}; Γ_{θ}; θ; Γ_{θ}^{−}; Γ_{θ}; Γ_{θ}; Γ_{θ}^{+}; Q_{θ}]

Abbildung 1.1: Geometrie des Würfels Q_{θ}

elektromagnetischen Feld „bestrahlt“. Gemessen wird das resultierende elektrische Feld an allen Seiten des Würfels. Formal lässt sich der Streuprozess wie folgt beschreiben.

1.3.1 Mathematische Formulierung des Streuproblems

Die Ausbreitung elektromagnetischer Felder in einem inhomogenen, isotropen und nicht perfekten Leiter^{1} wird durch die Maxwell–Gleichungen der folgenden Form beschrieben:

rot E(x, t) + μ_{0} \frac{∂H(x, t)}{∂t} = 0,  rot H(x, t) − ε(x)\frac{∂E(x, t)}{∂t} = σ(x)E(x, t). (1.1)

Dabei sind E(x, t), H(x, t) das elektrische und magnetische Feld am Ort x ∈ ℝ^{3} zur Zeit t ∈ ℝ. Weiter ist die magnetische Permeabilität (Durchdringbarkeit) μ_{0} eine positive Konstante. Das ist keine wesentliche Einschränkung, da für die meisten Materialien die magnetische Permeabilität nahezu gleich ist, [Born,Wolf;1999] Seite 3. Ausgeschlossen werden durch die Einschränkung an μ_{0} zum Beispiel einige ferromagnetische Materialien. Die Funktion ε(x) ist eine reelle Funktion des Ortes und wird elektrische Permeabilität genannt. Außerhalb der Kugel D soll sie gleich der Konstanten ε_{0} sein. Die elektrische Permeabilität ist immer positiv. Als Kehrwert des elektrischen Widerstands wird die elektrische Leitfähigkeit durch die reelle positive Ortsfunktion σ(x) beschrieben. Das zu rekonstruierende Objekt in der Kugel D wird also beschrieben durch seine Materialeigenschaften ε(x) und σ(x).

1Ein perfekter Leiter besitzt keinen elektrischen Widerstand. Damit wäre die elektrische Leitfähigkeit σ = ∞.

Nimmt man einen ladungsfreien Raum an, so werden die beiden vektorwertigen Gleichungen (1.1) komplettiert durch zwei skalare Beziehungen:

div(ε(x)E(x, t)) = 0, div H(x, t) = 0.

Ist das elektromagnetische Feld zeitharmonisch mit der Frequenz ω, lässt es sich als

E(x, t) = \frac{1}{\sqrt{ε_{0}}} E(x) e^{−iωt}, H(x, t) = \frac{1}{\sqrt{μ_{0}}} H(x) e^{−iωt},

schreiben. Die Ortsfunktionen E(x) und H(x) sind damit Lösung der zeitharmonischen Maxwell–Gleichungen

rot E(x) − ikH(x) = 0,  rot H(x) + ikn(x)E(x) = 0.  (1.2)

Dabei sind die Konstante k und die Funktion n(x) definiert durch

k^{2} := ε_{0}μ_{0}ω^{2},  n(x) := \frac{1}{ε_{0}} (ε(x) + i\frac{σ(x)}{ω}) .  (1.3)

Man nennt k die Wellenzahl und n(x) den Brechungsindex, der die reellen Materialeigenschaften (elektrische Permeabilität und Leitfähigkeit) umkehrbar–eindeutig zusammenfasst. Setzen wir voraus, dass für das magnetische Feld H ∈ C^{2}(ℝ^{3}) gilt, liefert die Anwendung des div-Operators auf die zweite Gleichung von (1.2) wegen der dann gültigen Operator–Identität div rot = 0 die skalare Gleichung

div(n(x)E(x)) = 0.  (1.4)

Die Forderung div rot H = 0 ist nach [Born,Wolf;1999] Kapitel 1 übrigens gleichbedeutend mit der Forderung nach der Ladungserhaltung in der Umgebung jedes Punktes (Kontinuitätsgleichung). Das Streuproblem lässt sich nun wie folgt formulieren: Seien die so genannten einfallenden Wellen E^{i} und H^{i} Lösung des homogenen Streuproblems, das heißt Lösung von

rot E(x) − ikH(x) = 0,  rot H(x) + ikE(x) = 0,  (1.5)

dann ist eine Lösung E(x), H(x) von (1.2) gesucht, so dass

E(x) = E^{i}(x) + E^{s}(x),  H(x) = H^{i}(x) + H^{s}(x)  (1.6)

gilt und die gestreuten Felder E^{s}, H^{s} die Silver–Müller Ausstrahlungsbedingung

lim_{|x|→∞} H^{s} × x − |x|E^{s} = 0  (1.7)

gleichmäßig für alle Richtungen x/|x| erfüllen. In [Colton,Kress;1992] wird das Streuproblem auf die gleiche Weise definiert. Legen wir uns auf die Einfallsrichtung θ = (0, 0, 1)^{t} fest, so können die einfallenden ebenen Wellen als

E^{i}(x) = ( \begin{matrix}a_{1} e^{ikx_{3}+δ_{1}} \\ a_{2} e^{ikx_{3}+δ_{2}} \\ 0 \end{matrix} ) , H^{i}(x) = ( \begin{matrix}−a_{2} e^{ikx_{3}+δ_{2}} \\ a_{1} e^{ikx_{3}+δ_{1}} \\ 0 \end{matrix} )

mit a_{1}, a_{2}, δ_{1}, δ_{2} ∈ ℝ gewählt werden.^{2} Wegen div E^{i} = 0 dürfen nur die ersten beiden Komponenten von E^{i} von Null verschieden sein. Man sagt, dass die ebene einfallende Welle elliptisch polarisiert ist. Im Allgemeinen beschreibt E^{i} als Funktion von x_{3} eine Ellipse in der (x_{1}, x_{2})–Ebene. Das kleinstmögliche Rechteck in der (x_{1}, x_{2})–Ebene, das diese Ellipse enthält hat die Seitenlängen 2|a_{1}| und 2|a_{2}|. Es werden noch drei Spezialfälle gesondert bezeichnet:

• Ist δ_{2} − δ_{1} = mπ mit m ∈ ℤ, dann wird E^{i} als linear polarisiert bezeichnet. Die Ellipse entartet zu einer Geraden.

• Ist a_{1} = a_{2} und δ_{2} − δ_{1} = π/2, dann wird E^{i} als rechtszirkular bezeichnet.

• Ist a_{1} = a_{2} und δ_{2} − δ_{1} = −π/2, dann wird E^{i} als linkszirkular bezeichnet.

Weitere Details findet man in [Born,Wolf;1999] Abschnitt 1.4.2. Da sich eine allgemein elliptisch polarisierte Welle aus der Überlagerung von zwei linear polarisierten Wellen zusammensetzen lässt, werden wir die einfallende ebene Welle immer als linear polarisiert annehmen. Ihre Ausbreitungsrichtung sei so gewählt, dass sie mit der Ausrichtung des Würfels Q_{θ} übereinstimmt:

E^{i}(x) := E_{0}^{i} e^{ikθ·x}, E_{0}^{i} ∈ ℝ^{3} und θ ⊥ E_{0}^{i}.

Wegen div E^{i} = 0 ist die Orthogonalität von E_{0}^{i} und θ eine notwendige Bedingung.

Bemerkung

Es genügt entweder das gestreute elektrische Feld E^{s} oder das gestreute magnetische Feld H^{s} zu kennen. Denn durch (1.2) ist das jeweilige andere Feld sofort bestimmt.

2Legt man E^{i}(x) fest, ist H^{i}(x) nach (1.5) festgelegt.

1.3.2 Beschreibung des inversen Streuproblems

Während beim direkten Streuproblem bei bekannter Wellenzahl k und bekanntem Brechungsindex n(x) das elektrische Feld E(x) gesucht ist, ist bei dem inversen Streuproblem der Brechungsindex n(x) unbekannt. Gegeben ist dafür neben der Wellenzahl k das elektrische Feld E(x) auf der Oberfläche des Würfels Q_{θ} für endlich viele verschiedene Orientierungen θ ∈ S^{2}. Für die konkrete Lösung des inversen Streuproblems für die Maxwell–Gleichungen wird hier ein Propagations– Rückpropagations Algorithmus (PBP–Algorithmus) hergeleitet. Er wurde so schon für die Helmholtz–Gleichung im Fall der Ultraschall–Tomographie implementiert, siehe [Natterer,Wübbeling;1995]. Der Algorithmus ist iterativ. Das heißt: Für P verschiedene Orientierungen θ = θ_{j} ∈ S^{2}, j = 0, …, P − 1, ergeben sich P verschiedene Messanordnungen für das Objekt in dem Würfel Q_{j} ⊃ D. Werden nun Messdaten G_{j}(x) auf ∂Q_{j} für das elektrische Feld E genommen, erhalten wir ein Randwertproblem für das Innere des Würfels Q_{j}. Nun greifen wir zu einem Trick. Mit den Messdaten G_{j}(x) auf ∂Q_{j} ist nämlich auch das Außenraumproblem in ℝ^{3} \ \overline{Q_{j}} lösbar. Denn mit Hilfe von Satz 2.9 ist das elektrische Feld E in ℝ^{3} \ \overline{Q_{j}} berechenbar, so dass zusätzlich noch ∂G_{j}(x)/∂θ_{j} auf Γ_{j}^{−} zu Verfügung steht. Mit den Daten G_{j}(x) auf Γ_{j} ∪ Γ_{j}^{−} und ∂G_{j}(x)/∂θ_{j} können wir ein Anfangswertproblem formulieren. Sei hier vorausgesetzt, dass wir dieses numerisch lösen können ^{3}, dann erhalten wir bei gegebenem Brechungsindex eine nichtlineare Abbildung R_{j} aus der Menge der „Brechungsindizes“ auf die Menge der „Messdaten“ am Rand Γ_{j}^{+}. Gesucht ist also der Brechungsindex n(x) in D, für den

R_{j}(n) = G_{j}_{|Γ_{j}^{+}},  j = 0, …, P − 1,  (1.8)

gilt. Dabei ist G_{j}_{|Γ_{j}^{+}} die Einschränkung von G_{j} auf dem Rand Γ_{j}^{+}. Wenn nun anstatt des exakten Brechungsindex n nur eine Näherung ñ bekannt war, unterscheidet sich das gemessene elektrische Feld auf Γ_{j}^{+} von R_{j}(ñ). Die Idee ist nun, diese Differenz zu benutzen, um die Näherung ñ zu verbessern. Obwohl die Normalableitung von G_{j} auf Γ_{j}^{−} aus den Messdaten G_{j} auf Q_{j} berechenbar ist, wird im weiteren Text aus Gründen der einfacheren Schreibweise dennoch von Messdaten ∂G_{j}/∂θ_{j} auf Γ_{j}^{−} die Rede sein. In Grundzügen wird der Algorithmus zur Lösung des inversen Streuproblems nun skizziert: Weil keine explizite Inversionsformel von (1.8) existiert, wird diese Abbildung linearisiert:

R_{j}(n + d) = R_{j}(n) + A_{j}(n)d + O(∥d∥^{2}).

3Die numerische Lösung des Anfangswertproblems ist Thema in Abschnitt 2.4.

Hierbei ist A_{j}(n) die Fréchet–Ableitung von R_{j}. Ausgehend von einer Startnäherung n^{0} werden iterativ bessere Näherungen n^{r+1} = n^{r} + ωd^{r}, r = 0, 1, … durch Lösen von

A_{j}(n^{r})d^{r} = G_{j}^{+} − R_{j}(n^{r}), j = r mod P

gebildet. Da dieses Gleichungssystem stark unterbestimmt ist, wird seine Minimum– Norm–Lösung bestimmt. Mit der adjungierten Abbildung A_{j}(n)^{∗} von A_{j}(n) schreibt sich die Lösung als

d^{r} = A_{j}(n^{r})^{∗}(A_{j}(n^{r})A_{j}(n^{r})^{∗})^{−1}(G_{j}^{+} − R_{j}(n^{r})), j = r mod P .

Hierbei ist vorausgesetzt worden, dass (A_{j}(n^{r})A_{j}(n^{r})^{∗})^{−1} existiert. Offensichtlich muss für jedes Update d^{r} das direkte Problem gelöst werden. Aus numerischer Sicht ist das keine leichte Aufgabe. Deshalb wird die Lösung des direkten Problems Gegenstand des nächsten Kapitels sein.

1.4 Literaturüberblick

Wie durch die Literaturangaben in Kapitel 2 deutlich werden wird, ist die mathematische Beschreibung des direkten Streuproblems für die Maxwell–Gleichungen, intensiv in den Büchern [Colton,Kress;1983] und [Colton,Kress;1992] untersucht worden. Die folgenden Abschnitte geben einen Überblick zur aktuellen Literatur zum Thema des inversen Streuproblems für die Maxwell–Gleichungen.

1.4.1 Rekonstruktion der Form und Oberflächenimpedanz eines Objektes

Rainer Kress und William Rundell zeigen in [Kress,Rundell;2000] wie gleichzeitig die Form und Impedanz eines nicht durchdringbaren Objektes aus dem Fernfeld eines Satzes ebener einfallender akustischer oder elektromagnetischer Wellen u_{i} berechnet werden kann. Die Autoren beschränken sich dabei weiter auf unendlich lange Objekte mit Zylindergeometrie und beschränktem Querschnitt D ⊂ ℝ^{2}, so dass die mathematische Beschreibung und numerische Implementierung für den ℝ^{2} erfolgen kann. Weiter liegt in der Natur des Algorithmus, dass die Anzahl der zusammenhängenden Objekte a priori bekannt sein muss und die Objekte Sterngebiete sind. Das hier vorliegende Außenraumproblem wird durch die Helmholtz–Gleichung

Δu + k^{2}u = 0 in ℝ^{2} \ \overline{D},  (1.9)

k > 0, mit der Randbedingung

\frac{∂u}{∂ν} + iλku = 0 auf ∂D  (1.10)

für die Impedanz beschrieben. Für den akustischen Fall sind diese Gleichungen sofort eingängig. Für die Maxwell–Gleichungen kommen wir auf diese skalare Formulierung, wenn das magnetische Feld im so genannten TM–Mode (transvers magnetic) polarisiert wird. Damit ist E = (0, 0, u)^{t} und E^{i} = (0, 0, u^{i})^{t}, wenn (0, 0, 1)^{t} die Achse des Zylinders darstellt. Gleichung (1.9) folgt sofort aus den Maxwell–Gleichungen und Gleichung (1.10) aus der eigentlichen Randbedingung für die Impedanzbedingung elektromagnetischer Felder

ν × rot E − iλ(ν × E) × ν = 0

auf ∂D durch einfaches Nachrechnen. Die Randbedingung (1.10) impliziert für λ = 0 die Neumann–Randbedingung ∂u/∂ν = 0 auf ∂D für einen perfekten Leiter und für λ → ∞ die Dirichlet–Randbedingung u = 0 auf ∂D. Zusammen mit der Sommerfeldschen Ausstrahlungsbedingung für das gestreute Feld u^{s} = u − u^{i}

lim_{r→∞} r (\frac{∂u^{s}}{∂r} − iku^{s}) = 0,  r = |x|,  (1.11)

gleichmäßig für alle Richtungen x/|x|, geben (1.9) und (1.10) Anlass zu einer nichtlinearen Abbildung F des Randes ∂D und der Impedanz λ in das Fernfeld. Für den elektromagnetischen Fall geht die Silver–Müller Ausstrahlungsbedingung (1.7) direkt in die Sommerfeldsche Ausstrahlungsbedingung (1.11) über. Das zugehörige inverse Problem für den Operator F ist schlecht gestellt. Seine Lösung basiert auf einem regularisierten iterativen Newton–Verfahren, welches die Berechnung in zwei Schritte zerlegt: Zuerst wird das gestreute Feld u^{s} für eine Approximation von ∂D und λ durch das Lösen einer linearen Integralgleichung erster Ordnung auf ∂D berechnet, was ein sehr schlecht gestelltes Problem ist. Dann folgt die Invertierung einer nichtlinearen Gleichung durch Linearisierung. Der letzte Schritt führt iterativ auf ein einfaches „Update“ für ∂D und λ. Die Autoren bemerken dazu, dass ihre Aufteilung des Algorithmus in zwei Schritte im wesentlichen Vorteile für die numerische Implementierung beinhaltet. Die Anwendbarkeit dieses Verfahrens wird durch die Qualität einer vorgestellten Beispielrechnung dargelegt. Die Autoren bemerken weiter, dass die Ergebnisse mit zunehmendem Wert für die Impedanz besser und mit zunehmender Spanne des Wertebereichs der Impedanz schlechter werden.

Weiter wird bemerkt, dass die geforderte Menge an Messdaten mit zunehmender Komplexität von ∂D steigen muss, um keine Abstriche bei der Qualität der Rekonstruktion machen zu müssen. Anders gesagt liefert eine konstante Menge an Messdaten schlechtere Ergebnisse, wenn ∂D aus einer größeren Anzahl von „Basisfunktionen“ zusammengesetzt ist. Vorwissen über die Komplexität der zu rekonstruierenden Oberfläche ist somit von Vorteil. Diese Arbeit von Kress und Rundell [Kress,Rundell;2000] ist in dem Kontext von [Kress;2000] zu sehen, in der das inverse Problem darin besteht, lokale Störungen in der Gestalt einer unendlich großen, ideal leitenden planaren Platte durch Messungen des Fernfelds nach Einfall ebener zeitharmonischer Felder E^{i} zu rekonstruieren.

1.4.2 Rekonstruktion der Form von allgemeineren Objekten

David Colton und Andreas Kirsch stellen in ihrer Arbeit [Colton,Kirsch;1996] einen Algorithmus vor, der die Rekonstruktion der Form von endlich vielen Objekten erlaubt, egal ob sie für das elektrische Feld durchlässig sind oder nicht. Dabei versteht man hier unter einem Objekt eine zusammenhängende Menge, die disjunkt mit den anderen Objekten ist. Außerdem befinden sich dabei alle Objekte in einem bekannten Raumvolumen. Die Daten, die in diesem Algorithmus einfließen, bestehen aus dem Fernfeld für alle Richtungen. Die Autoren beschränken sich allerdings auch hier auf unendlich lange Objekte mit Zylindergeometrie, wobei die Symmetrieachse für alle Objekte dieselbe ist. Damit ist wieder eine skalare Beschreibung der Maxwell–Gleichungen in ℝ^{2} im TM–Mode wie im Abschnitt 1.4.1 möglich. Das Gesamtfeld u erfüllt somit

Δu(x) + k^{2}u(x) = 0 in ℝ^{2} \ \overline{D},
u(x) = e^{ikθ·x} +u^{s}(x)

und u^{s} die Sommerfeldsche Ausstrahlungsbedingung gleichmäßig für alle Richtungen x/|x|. Weiter ist die Wellenzahl k > 0 konstant und θ ∈ S^{2} die Richtung der einfallenden Welle. Ein Objekt D muss hier durch eine der folgenden physikalischen Eigenschaften beschrieben werden können:

1. Es erfüllt die Dirichlet–Randbedingung

u = 0 auf ∂D,

2. die Neumann–Randbedingung

\frac{∂u}{∂ν} = 0 auf ∂D,

mit der äußeren Normalen ν,

3. die Impedanz-Randbedingung

\frac{∂u}{∂ν} + λu = 0 auf ∂D,  (1.12)

wobei λ > 0 ist, oder

4. für den Brechungsindex n(x) (s. Gl. (1.3)) muss für mindestens ein x ∈ D Im(n(x)) > 0 gelten.

Wenn nun k^{2} für die ersten drei Bedingungen weder ein Dirichlet– noch ein Neumann–Eigenwert und für den vierten Fall kein Dirichlet–Eigenwert ist, zeigen die Autoren die Existenz einer Funktion g(·, y_{0}) ∈ L_{2}(−π, π) mit der Eigenschaft

lim_{y_{0}∈∂D, y→y_{0}} ∥g(·, y)∥_{L_{2}(−π,π)} = ∞.

Die Funktion g ist dabei Lösung der linearen Integralgleichung

∫_{−π}^{π} u_{∞}(ϕ, θ)g(θ, y_{0}) dθ = e^{−ikρ cos(ϕ−α)} für alle ϕ ∈ [−π, π],  (1.13)

wobei y_{0} = ρ(cos α, sin α)^{t}. Dieses Resultat legt folgenden einfachen Algorithmus nahe: Finde die Funktion g aus (1.13) in einem Raumbereich P, von dem a priori bekannt ist, dass er alle Objekte enthält. Suche dann die Punkte y_{0} ∈ P für die ∥g(·, y_{0})∥_{L_{2}(−π,π)} sein „Maximum“ annimmt. Die Menge dieser Punkte definiert dann die Ränder der gesuchten Objekte. Hier wird sofort ein Schwachpunkt des Verfahrens deutlich: Das Kriterium (1.12) lässt sich numerisch schwer fassen. Man kann vom gesetzten Maximalwert (Abbruchkriterium) für g(·, y_{0}) ∈ L_{2}(−π, π) nicht auf die Güte der Näherung schließen.

1.4.3 Rekonstruktion von im Boden eingelassenen zylindrischen Objekten

In dem Artikel [Hiraoka et al;1996] wird die Brauchbarkeit elektromagnetischer Wellen zur Detektion von im Boden eingelassenen Objekten untersucht. Auch hier wird das Problem, wie in Abschnitt 1.4.1 dargelegt, auf zwei Dimensionen reduziert. Die Halbebene y ≤ 0 entspricht dabei der Luft, die mit dem konstanten Brechungsindex n_{0} modelliert wird. Im Boden (y > 0) eingelassen, ist nun ein Objekt Ω, das durch den komplexen Brechungsindex n_{s}(x, y) ausgezeichnet ist, wobei der Boden selbst den konstanten Brechungsindex n_{1} haben soll. Die Geometrie lässt sich also wie folgt durch den Brechungsindex n(x, y) wiedergeben:

n(x, y) = \begin{cases}n_{0} & falls y ≤ 0, \\ n_{s}(x, y) & falls y ∈ Ω, \\ n_{1} & sonst. \end{cases}

Bei bekanntem n_{0}, n_{1} soll die Funktion n_{s}(x, y) aus Messungen des gestreuten elektrischen Feldes Ẽ^{s} an der Oberfläche (y < 0) berechnet werden. Genauer hat man für P verschiedene diskrete Einfallswinkel θ_{0} = θ_{1}, …, θ_{P} der einfallenden elektrischen Welle E^{i}(x, y) Q Messungen an den Stellen (x_{q}, y_{0}), q = 1, …, Q, y_{0} < 0. In Abwesenheit des Streuers Ω bzw. n_{s} = n_{1} schreibt sich das elektrische Feld als

E(θ_{0}, x, y) = \begin{cases}e^{−ik_{0}(x sin θ_{0}+y cos θ_{0})} +R e^{−ik_{0}(x sin θ_{0}+y cos θ_{0})} & falls y ≤ 0, \\ T e^{−ik(x sin θ+y cos θ)} & falls y > 0. \end{cases}

Dabei sind k_{0}, k = n_{1}k_{0} die Wellenzahlen in Luft und Boden, θ_{0}, θ = arcsin(n^{−1} sin θ_{0}) Einfalls- bzw. Brechungswinkel und die Konstanten R, T Reflexions- bzw. Transmissionskoeffizient. Im Allgemeinen gilt für das totale elektrische Feld

E^{t}(n_{s}, θ_{0}, x, y) = E(θ_{0}, x, y) + E^{s}(n_{s}, θ_{0}, x, y),

wobei das gestreute Feld E^{s} die lineare Integralgleichung

E^{s}(n_{s}, θ_{0}, x, y) = ∫_{Ω} k_{0}^{2}[n_{s}^{2}(x′, y′) − n^{2}]E^{t}(n_{s}, θ_{0}, x′, y′)G(x, y, x′, y′) dx′dy′ (1.14)

mit der zweidimensionalen Greenschen Funktion G(x, y, x′, y′) erfüllt. G ist dabei die Greensche Funktion für die Brechung an einer Halbebene. Das inverse Problem kann nun durch Minimierung des Kostenfunktionals

J(n_{s}) = ∑_{p=1}^{P} ∑_{q=1}^{Q} |E^{s}(n_{s}, θ_{p}, x_{q}, y_{0}) − Ẽ^{s}(θ_{p}, x_{q}, y_{0})|^{2}

gelöst werden. Dazu verwenden die Autoren das Verfahren der konjugierten Gradienten, wobei der Gradient von J(n_{s}) durch die Fréchet–Ableitung des Funktionals gegeben ist. Numerische Ergebnisse für konstantes n_{s} demonstrieren die prinzipielle Umsetzbarkeit des Verfahrens.

Bei der Diskretisierung wurde allerdings nur ein (16×16)-Gitter verwendet. Für eine allgemeine (n × n)-Diskretisierung werden alleine zum Lösen der Integralgleichung (1.14) mindestens O(n^{4}) Operationen benötigt, was den Algorithmus mit zunehmender feinerer Diskretisierung schnell unbrauchbar werden lässt. Würde Gleichung (1.14) zum Beispiel auf einem 32 × 32 Gitter gelöst, hätte sich die Rechenzeit wenigstens um das 16fache erhöht. Die Brauchbarkeit der in [Hiraoka et al;1996] vorgestellte Methode kann somit zumindest für Rechnungen im ℝ^{3} in Frage gestellt werden.

1.4.4 Dreidimensionale Mikrowellen–Tomographie

Die Arbeit [Bulyshev et al;2000] gehört zu den wenigen Arbeiten, die Mikrowellen– Tomographie im dreidimensionalen Raum behandeln. In ihrem mathematischem Modell approximieren sie die Maxwell–Gleichungen durch die skalare Helmholtz– Gleichung

Δu_{j}^{s} + k^{2}u_{j}^{s} = (k_{0}^{2} − k^{2})u_{j}^{i}.

Dabei ist u_{j}^{i} das einfallende Feld, welches vom j-ten Transmitter (Punktquelle) erzeugt worden ist. Befindet sich dieser an der Stelle r_{j} ∈ ℝ^{3}, so gilt

u_{j}^{i}(x) = \frac{1}{4π|x − r_{j}|} e^{ik_{0}|x−r_{j}|} .

Das gestreute elektrische Feld am l-ten Receiver wird durch

u_{j,l}^{s} = ∫ \frac{1}{4π|x − r_{l}|} e^{ik_{0}|x−r_{l}|}(k^{2} − k_{0}^{2})(u_{j}^{s}(x) + u_{j}^{i}(x)) dx

berechnet. Das inverse Problem kann nun als Variationsproblem formuliert werden:

J[k^{2}] = (∑_{j,l} |u_{j,l}^{theor} − u_{j,l}^{exper}|^{2}) => min.

Hierbei ist bei aktiviertem Transmitter j der Term u_{j,l}^{theor} die theoretische Voraussage für das Feld bei Receiver l und u_{j,l}^{exper} das bei Receiver l gemessene Feld. Um das Funktional J zu minimieren, wählen die Autoren die Gradienten–Methode. Dazu ist es nötig die Ableitung J′ in ℝ^{3} berechnen zu können. Werden dafür die für den zweidimensionalen Fall bestehenden Algorithmen einfach portiert, ist das Verfahren für die Praxis wohl nicht mehr geeignet, weil zu zeitaufwendig. Die Autoren stellen nun einen effizienten Algorithmus für die Berechnung von J′ in ℝ^{3} vor. Numerische Simulationen demonstrieren die Brauchbarkeit ihres Verfahrens. Sie verwenden zum Beispiel ein Herzmodell, das hier in Abschnitt 4.3 qualitativ als Vorlage für eine numerische Simulation diente.

1.4.5 Bemerkung

Würde sich auch diese Arbeit auf Objekte mit Zylindergeometrie beschränken, das heißt eine geometrische Formulierung wie in Abschnitt (1.4.1) heranziehen, oder die Kopplung der drei elektrischen Feldkomponenten E = (E_{1}, E_{2}, E_{3})^{t} vernachlässigen wie in [Bulyshev et al;2000], dann führten die Maxwell–Gleichungen direkt zur skalaren Helmholtz–Gleichung. Damit läge mathematisch der Fall der Ultraschall–Tomographie vor. Dies würde den in Abschnitt 1.3.2 skizzierten Algorithmus zur Lösung des inversen Streuproblems auf die Arbeiten von Frank Natterer und Frank Wübbeling zurückwerfen, siehe zum Beispiel [Natterer,Wübbeling;1995], [Natterer;1999]. Auch soll nochmal deutlich gesagt werden, dass es hier nicht nur um die Rekonstruktion der Geometrie von supp(n − 1) geht, sondern um die Rekonstruktion der komplexwertigen Funktion n selbst.

Das direkte Problem

Ist das elektrische Feld Lösung des direkten Streuproblems, wie in Abschnitt 1.3.1 beschrieben, dann hängt seine Regularität von der des Brechungsindex ab. Um diese Abhängigkeit darstellen zu können, wird der Begriff der Hölder–Stetigkeit gebraucht.

Definition 2.1 Sei G ⊂ ℝ^{3} offen. Man nennt eine Funktion f auf G gleichmäßig Hölder–stetig mit Exponenten 0 < α ≤ 1, wenn eine Konstante C > 0 existiert, so dass

|f(x) − f(y)| ≤ C|x − y|^{α}

für alle x, y ∈ G gilt. Für α = 1 nennt man f Lipschitz–stetig. Die Hölder– Räume C^{p,α}(G) sind nun definiert als die Menge der beschränkten Funktionen aus C^{p}(G), deren Ableitungen der Ordnung p gleichmäßig Hölder–stetig mit Exponent α sind.

Existenz und Eindeutigkeit des direkten Streuproblems und die stetige Abhängigkeit der Lösung vom einfallenden Feld können als Satz formuliert werden:

Satz 2.2 Sei n ∈ C^{1,α}(ℝ^{3}) mit dem Exponenten 0 < α < 1 und supp(n − 1) kompakt. Dann besitzt das Streuproblem (1.2), (1.5), (1.6) und (1.7) eine eindeutige Lösung E, H ∈ C^{1}(ℝ^{3}), die stetig von dem einfallenden Feld E^{i}, H^{i} im Sinne der Maximumsnorm abhängt.

Beweis: Siehe [Colton,Kress;1992] Theorem 9.4 und Theorem 9.5. □

Werden explizit Lösungen des direkten Streuproblems in C^{1,α}(ℝ^{3}) gesucht, hat man:

2.1. REGULARITÄT DES ELEKTRISCHEN FELDES E 15

Satz 2.3 Sei der Brechungsindex n wie in Satz 2.2. Dann ist das Lösen des Streuproblems (1.2), (1.5), (1.6) und (1.7) mit Lösungen in C^{1,α}(ℝ^{3}) äquivalent mit dem Lösen der Integralgleichung

E(x) = E^{i}(x) − k^{2} ∫_{ℝ^{3}} Φ(x, y)(1 − n(y))E(y) dy

+ ∇ ∫_{ℝ^{3}} \frac{∇n(y)}{n(y)} · E(y)Φ(x, y) dy (2.1)

in ℝ^{3}, wobei Φ(x, y) = e^{ik|x−y|} /4π|x − y| die Fundamentallösung der Helmholtz– Gleichung ist.

Beweis: Siehe [Colton,Kress;1992] Theorem 9.1 und 9.2. □

Wir werden das direkte Problem ausschließlich für das elektrische Feld angehen. Im Folgenden werden höhere Regularitätsvoraussetzungen an das elektrische Feld E gestellt. Wir wollen untersuchen, was die hinreichenden Bedingungen dafür sind.

2.1 Regularität des elektrischen Feldes E

Mit Blick auf die Integralgleichung (2.1) für E macht es Sinn, Aussagen über die Regularität des Volumenpotentials

u(x) := ∫_{ℝ^{3}} φ(y)\frac{1}{4π} \frac{e^{ik|x−y|}}{|x − y|} dy = ∫_{ℝ^{3}} φ(x − z)\frac{1}{4π} \frac{e^{ik|z|}}{|z|} dz  (2.2)

mit φ ∈ C_{0}(ℝ^{3}) zu treffen. Entscheidend ist nun der folgende Satz:

Satz 2.4 Sei 0 < α < 1. Das Volumenpotential u aus Gleichung (2.2) existiert als uneigentliches Integral für alle x ∈ ℝ^{3} und hat folgende Eigenschaften: Ist φ ∈ C_{0}(ℝ^{3}), dann ist u ∈ C^{1,α}(ℝ^{3}) und es darf unter dem Integralzeichen differenziert werden. Ist weiter φ ∈ C_{0}(ℝ^{3}) ∩ C^{p,α}(ℝ^{3}) für p ≥ 0, dann ist u ∈ C^{p+2,α}(ℝ^{3}).

Beweis: Die obigen Aussagen und die letzte für p = 0 sind [Gilbarg,Trudinger;1977] Kapitel 4 und [Colton,Kress;1992] Theorem 8.1 entnommen. Die Regularitätsaussage für p ≥ 1 ist dann ein Korollar. □

Mit Hilfe des obigen Satzes sieht man ein: Ist E ∈ C(ℝ^{3}) Lösung der Integralgleichung (2.1) mit n ∈ C^{1,α}(ℝ^{3}), so ist sogar E ∈ C^{1,α}(ℝ^{3}), denn n − 1 hat einen kompakten Träger. Wird nun die Regularität von n erhöht, nämlich n ∈ C^{2,α}(ℝ^{3}), dann folgt wieder aus obigen Satz und der Integralgleichung (2.1), dass E ∈ C^{2,α}(ℝ^{3}) gilt. Induktiv aus dem eben Gesagten und Satz 2.3 folgt:

Satz 2.5 Ist E ∈ C(ℝ^{3}) Lösung der Integralgleichung (2.1) und n ∈ C^{p,α}(ℝ^{3}) mit p ≥ 1, dann ist E ∈ C^{p,α}(ℝ^{3}) und Lösung des Streuproblems (1.2), (1.5), (1.6) und (1.7).

In der weiteren Darstellung wird die Annahme gemacht, dass für Lösungen der Maxwell–Gleichungen E ∈ C^{2}(ℝ^{3}) gilt. Beschränken wir uns auf den Teilraum C^{2,α}(ℝ^{3}) ⊂ C^{2}(ℝ^{3}), so ist es nach Satz 2.5 hinreichend, wenn für den Brechungsindex n ∈ C^{2,α}(ℝ^{3}) gilt. Damit kann die Operatorgleichung rot rot = −Δ + ∇ div auf E angewendet werden. Dies und die beiden Gleichungen in (1.2) ergeben eine Differentialgleichung zweiter Ordnung für das elektrische Feld E:

ΔE(x) − ∇ div E(x) + k^{2}n(x)E(x) = 0.  (2.3)

Für jedes n ∈ C^{2,α}(ℝ^{3}) erhalten wir also eine Lösung E ∈ C^{2,α}(ℝ^{3}) des direkten Streuproblems, die auch Gleichung (2.3) erfüllt. Andere Autoren (siehe zum Beispiel [Ammari,Bao;2001]) verzichten auf höhere Differenzierbarkeitseigenschaften und arbeiten mit der Differentialgleichung

rot rot E − k^{2}n(x)E(x) = 0.  (2.4)

Eine Lösung von (2.3) zusammen mit der Definition

H(x) := \frac{1}{ik} rot E(x)  (2.5)

ist andererseits wieder Lösung der Maxwell–Gleichungen (1.2). Denn (2.5) entspricht schon der ersten Gleichung aus (1.2) und durch Anwendung der Operatorgleichung rot rot = −Δ + ∇ div auf (2.3) erhält man die Darstellung (2.4). Nach der Definition (2.5) ist dies aber nichts anderes als

rot H + ikE(x) = 0,

was die Gültigkeit der zweiten Gleichung aus (1.2) zeigt. Unter der Annahme E ∈ C^{2}(ℝ^{3}) sind die Differentialgleichungen (2.3) und (1.2) also äquivalent.

2.2 Elliptizität

Die Differentialoperatoren rot rot aus (2.4) und Δ − ∇ div aus (2.3) sind nicht elliptisch. Um dies einzusehen werden einige Sprechweisen eingeführt: Sei P (ξ) ein Polynom mit konstanten Koeffizienten in den m Variablen ξ_{1}, …, ξ_{m}, dann verstehen wir unter der Bezeichnung P (D) den Differentialoperator, den man erhält, wenn die ξ_{j} durch ∂/∂x_{j} ersetzt werden. Die Ordnung von P bezeichnen wir mit ord P. Das charakteristische Polynom P^{0}(ξ) erhält man aus P (ξ) durch Streichung aller Terme kleinerer Ordnung als ord P. Das heißt, für ein konstantens Polynom ist P^{0}(ξ) = P (ξ) und für ein nicht konstantes Polynom ∂P^{0}(ξ)/∂ξ_{j} = ∂P (ξ)/∂ξ_{j} für alle j. Sei weiter

∑_{j=1}^{N} P_{i,j}(D)u_{j} = F_{i},  i = 1, …, N  (2.6)

ein System von N partiellen Differentialgleichungen mit konstanten Koeffizienten in N Unbekannten. In [Hörmander;1976] Kapitel 10 findet sich:

Definition 2.6 Das System (2.6) heißt elliptisch im Sinne von Douglis und Nirenberg, wenn es ganze Zahlen t_{j}, s_{i} ≥ 0, i, j = 1, …, N, gibt, so dass ord P_{i,j} ≤ t_{j} − s_{i} und

det(P̃_{i,j}^{0}(ξ)) ≠ 0,  falls 0 ≠ ξ ∈ ℝ^{m}.  (2.7)

Dabei gelte mit dem charakteristischen Polynom P_{i,j}^{0}(ξ) des Polynoms P_{i,j}(ξ):

P̃_{i,j}^{0}(ξ) := \begin{cases}P_{i,j}^{0}(ξ) & ord P_{i,j}^{0} = t_{j} − s_{i}, \\ 0 & sonst. \end{cases}

Für N = 1 erhält man die übliche Definition für einen elliptischen Operator. Es wird nun gezeigt, dass der Differentialoperator Δ − ∇ div nicht elliptisch ist. Angewendet auf E lautet er ausführlich

ΔE − ∇ div E = ( \begin{matrix}\frac{∂^{2}E_{1}}{∂x_{2}^{2}} + \frac{∂^{2}E_{1}}{∂x_{3}^{2}} − \frac{∂^{2}E_{2}}{∂x_{1}∂x_{2}} − \frac{∂^{2}E_{3}}{∂x_{1}∂x_{3}} \\ \frac{∂^{2}E_{2}}{∂x_{1}^{2}} + \frac{∂^{2}E_{2}}{∂x_{3}^{2}} − \frac{∂^{2}E_{1}}{∂x_{1}∂x_{2}} − \frac{∂^{2}E_{3}}{∂x_{2}∂x_{3}} \\ \frac{∂^{2}E_{3}}{∂x_{1}^{2}} + \frac{∂^{2}E_{3}}{∂x_{2}^{2}} − \frac{∂^{2}E_{1}}{∂x_{1}∂x_{3}} − \frac{∂^{2}E_{2}}{∂x_{2}∂x_{3}} \end{matrix} ) .

Seine zugehörigen Polynome P_{i,j}(ξ), i, j = 1, 2, 3, haben in einer Matrix zusammengefasst die Gestalt

P := ( \begin{matrix}ξ_{2}^{2} + ξ_{3}^{2} & −ξ_{1}ξ_{2} & −ξ_{1}ξ_{3} \\ −ξ_{1}ξ_{2} & ξ_{1}^{2} + ξ_{3}^{2} & −ξ_{2}ξ_{3} \\ −ξ_{1}ξ_{3} & −ξ_{2}ξ_{3} & ξ_{1}^{2} + ξ_{2}^{2} \end{matrix} ) .

Mit der Matrix

A = ( \begin{matrix}0 & −ξ_{3} & ξ_{2} \\ ξ_{3} & 0 & −ξ_{1} \\ −ξ_{2} & ξ_{1} & 0 \end{matrix} )

ist P = −A^{2}. Da ord P_{i,j} = 2, ist t_{j} − s_{i} ≥ 2, i, j = 1, 2, 3. Wählen wir zunächst t_{j} − s_{i} = 2 für alle i, j, so ist die Bedingung aus Gleichung (2.7) nicht erfüllt, denn es gilt det A = 0 und P̃_{i,j}^{0} = P_{i,j}. Für alle weiteren gültigen Konfigurationen von t_{j}, s_{i} wird mindestens immer eine Spalte von (P̃_{i,j}^{0})_{i,j} zu Null, so dass auch dann Bedingung (2.7) nicht gilt. Insgesamt ist somit gezeigt, dass der Differentialoperator Δ − ∇ div nicht elliptisch ist. Das macht eine Stabilitätsanalyse für die numerische Lösung des Anfangswertproblems mit den Methoden aus Abschnitt 2.4.1 unmöglich. Wir wollen deshalb eine andere Differentialgleichung für E herleiten:

Umstellen von Gleichung (1.4) liefert

div E(x) = −\frac{∇n(x)}{n(x)} · E(x).  (2.8)

Die rechte Seite von (2.8) ist aus physikalischen Gründen immer wohldefiniert, denn wegen der Beschränkung für die Lichtgeschwindigkeit ist nach [Born,Wolf;1999] Seite 11 Re n(x) > 0 für alle x ∈ ℝ^{3}. Die Gleichungen (2.3) und (2.8) lassen sich nun kombinieren:

ΔE(x) + ∇ (\frac{∇n(x)}{n(x)} · E(x)) + k^{2}n(x)E(x) = 0.  (2.9)

Der Differentialoperator obiger Gleichung ist von der Form

Δ + ∑_{|α|≤1} a_{α}(x)\frac{∂^{|α|}}{∂x_{1}^{α_{1}}∂x_{2}^{α_{2}}∂x_{3}^{α_{3}}} ,

wobei α = (α_{1}, α_{2}, α_{3}) und |α| = α_{1}+α_{2}+α_{3} mit α_{i} ∈ ℕ_{0}, i = 1, 2, 3. Er ist also eine gestörte Form des elliptischen Operators Δ und damit selbst elliptisch. Außerhalb von supp(n − 1) reduziert sich die Differentialgleichung (2.9) sogar auf die Helmholtz–Gleichung. Gleichung (2.9) ist die zentrale Gleichung für weitere Untersuchungen. Sie wurde in der Literatur über „Streuung elektromagnetischer Wellen“ schon früher verwendet, siehe zum Beispiel [Semenov et al;1999] oder [Born,Wolf;1999] Kapitel 13. Dort wird allerdings nicht darauf eingegangen, ob eine Lösung von (2.9) die Maxwell–Gleichungen (1.2) bzw. Gleichung (2.3) erfüllt, was keineswegs auf der Hand liegt. Denn für eine Lösung E der Gleichung (2.9) ist nicht einsichtig, dass auch div(nE) = 0 gilt. Dies zu beweisen, ist Gegenstand des nächsten Abschnitts. Außerdem wird eine Stabilitätsanalyse in Abschnitt 2.4.1 zeigen, dass Gleichung (2.9) für die numerische Lösung des Anfangswertproblems geeignet ist.

2.3. ÄQUIVALENTE FORMULIERUNG DER MAXWELL–GLEICHUNGEN 19

2.3 Äquivalente Formulierung der Maxwell–Gleichungen

Zunächst werden einige Lemmata zitiert:

Lemma 2.7 Erfüllen Lösungen der homogenen Maxwell–Gleichungen (1.5) die Silver–Müller Ausstrahlungsbedingung (1.7), so erfüllen die kartesischen Komponenten von E die Sommerfeldsche Ausstrahlungsbedingung (1.11). Ist umgekehrt E Lösung der homogenen Helmholtz–Gleichung und div E = 0, dann erfüllen E und H := rot E/ik die Silver–Müller Ausstrahlungsbedingung (1.7), wenn die kartesischen Komponenten von E die Sommerfeldsche Ausstrahlungsbedingung (1.11) erfüllen.

Beweis: Siehe [Colton,Kress;1983] Korollar 4.3 und 4.7. □

Lemma 2.8 Sei E Lösung der Helmholtz–Gleichung

ΔE + k^{2}E = 0.

Dann ist äquivalent:

• E erfüllt die Ausstrahlungsbedingung

rot E × \frac{x}{|x|} + \frac{x}{|x|} div E − ikE = o(|x|^{−1})  (2.10)

gleichmäßig für alle Richtungen x/|x|.

• Die kartesischen Komponenten von E erfüllen die Sommerfeldsche Ausstrahlungsbedingung (1.11).

Beweis: Siehe [Colton,Kress;1983] Korollar 4.14. □

Lemma 2.9 Sei E ∈ C^{2}(ℝ^{3} \ \overline{D}) ∩ C(\overline{ℝ^{3} \ D}) Lösung der Helmholtz–Gleichung

ΔE + k^{2}E = 0,

wobei zusätzlich gefordert wird, dass die Funktionen div E, rot E stetig auf den Rand von ℝ^{3} \ \overline{D} fortgesetzt werden können und E der Ausstrahlungsbedingung

rot E × \frac{x}{|x|} + \frac{x}{|x|} div E − ikE = o(|x|^{−1})

gleichmäßig für alle Richtungen x/|x| genügt. Dann ist

rot ∫_{∂D} ν(y) × E(y)Φ(x, y) ds(y) − ∇ ∫_{∂D} ν(y) · E(y)Φ(x, y) ds(y)

− ∫_{∂D}[rot E(y) × ν(y) + ν(y) div E(y)]Φ(x, y) ds(y) = \begin{cases}0 & falls x ∈ D, \\ E(x) & falls x ∈ ℝ^{3} \ \overline{D}. \end{cases}

Dabei ist ν(y) die äußere Normale von D an der Stelle y ∈ ∂D und Φ(x, y) = exp(ik|x − y|)/4π|x − y| die Fundamentallösung der Helmholtz–Gleichung.

Beweis: Siehe [Colton,Kress;1983] Theorem 4.13. □

Bemerkung

Das Lemma 2.7 bezieht sich eigentlich nur auf Lösungen der homogenen Maxwell– Gleichungen (1.5). Weil hier supp(n − 1) ⊂ D als kompakte Menge vorausgesetzt ist, erfüllt eine Lösung E, H der allgemeinen Maxwell–Gleichungen (1.2) im Außenraum ℝ^{3} \ \overline{D} die homogene Maxwell–Gleichung (1.5), so dass zum Beispiel durch Anwendung von Theorem 4.1 aus [Colton,Kress;1983] eine Lösung der homogenen Maxwell–Gleichung in ganz ℝ^{3} konstruiert werden kann, die im Außenraum ℝ^{3} \ \overline{D} mit der Lösung der allgemeinen Maxwell–Gleichung übereinstimmt. Somit gilt Lemma 2.7 hier auch für die allgemeine Maxwell–Gleichung und wegen der höheren Regularitätsvoraussetzung E ∈ C^{2}(ℝ^{3}) auch für die Differentialgleichung zweiter Ordnung (2.3). Analoge Argumente gelten auch für Lemma 2.8.

Das folgende Lemma wird den Äquivalenzbeweis zwischen den Maxwell–Gleichungen (1.2) und der Differentialgleichung zweiter Ordnung (2.9) für das elektrische Feld E auf einen Eindeutigkeitsbeweis für eine Helmholtz–Gleichung zurückführen. Dabei wird allerdings eine höhere Glattheit für die Komponenten des elektrischen Feldes E und damit des Brechungsindex n verlangt, nämlich E, n ∈ C^{3}(ℝ^{3}). Dass für n ∈ C^{3,α}(ℝ^{3}) eine Lösung E ∈ C^{3}(ℝ^{3}) der Differentialgleichung (2.9) existiert, ist schon gezeigt worden, denn wir könnten die Lösung der Integralgleichung (2.1) nehmen.

Lemma 2.10 Sei E ∈ C^{3}(ℝ^{3}) Lösung von

ΔE(x) + ∇ (\frac{∇n(x)}{n(x)} · E(x)) + k^{2}n(x)E(x) = 0  (2.11)

und n ∈ C^{3}(ℝ^{3}), dann erfüllt

u(x) := \frac{1}{n(x)} div(n(x)E(x))  (2.12)

die Helmholtz–Gleichung

Δu + k^{2}nu = 0.

Beweis: Der Beweis wird durch einfaches Nachrechnen geführt. Wegen E ∈ C^{3}(ℝ^{3}) ist der div-Operator auf (2.9) anwendbar und kann mit dem Δ-Operator vertauscht werden:

Δ div E(x) + Δ (\frac{∇n(x)}{n(x)} · E(x)) + k^{2} div(n(x)E(x)) = 0.

Wegen n ∈ C^{3}(ℝ^{3}) ist weiter

Δu(x) = Δ (div E(x) + \frac{∇n(x)}{n(x)} · E(x)) = Δ div E(x) + Δ (\frac{∇n(x)}{n(x)} · E(x)) .

Wendet man nun die Definition (2.12) für u auf den dritten Term in (2.11) an, folgt

Δu + k^{2}nu = 0,

also die Behauptung. □

Ein notwendiges weiteres Resultat mit einem längeren Beweis wird hier ebenfalls in Form eines Lemmas wiedergegeben.

Lemma 2.11 Sei E wie in Lemma 2.9, dann erfüllt div E die Sommerfeldsche Ausstrahlungsbedingung, nämlich

\frac{x}{|x|} · ∇ div E(x) − ik div E(x) = o(|x|^{−1})

gleichmäßig für alle Richtungen x/|x|.

Beweis: Für x ∉ \overline{D} gilt Lemma 2.9:

E(x) = rot ∫_{∂D} ν(y) × E(y)Φ(x, y) ds(y) − ∇ ∫_{∂D} ν(y) · E(y)Φ(x, y) ds(y)

− ∫_{∂D}[rot E(y) × ν(y) + ν(y) div E(y)]Φ(x, y) ds(y).

Man beachte, dass hier auf der rechten Seite beliebig oft unter dem Integral differenziert werden darf. Wird nun die Operatorgleichung div rot = 0 und Δ_{x}Φ = −k^{2}Φ benutzt, folgt

div E(x) = k^{2} ∫_{∂D} ν(y) · E(y)Φ(x, y) ds(y)

− ∫_{∂D}[rot E(y) × ν(y) + ν(y) div E(y)] · ∇_{x}Φ(x, y) ds(y).

Nach Anwendung der Richtungsableitung x/|x| · ∇, schreibt sich obige Gleichung als

\frac{x}{|x|} · ∇ div E^{s}(x) = k^{2} ∫_{∂D} ν(y) · E^{s}(y)\frac{x}{|x|} · ∇Φ(x, y) ds(y)

− ∫_{∂D}[rot E^{s}(y) × ν(y) + ν(y) div E^{s}(y)] · (\frac{x}{|x|} · ∇) ∇_{x}Φ(x, y) ds(y).

Um die Sommerfeldsche Ausstrahlungsbedingung für div E zu überprüfen, subtrahiert man die beiden letzten Gleichungen voneinander und schätzt die Differenz für große |x| ab. Für x ∉ \overline{D} folgt

\frac{x}{|x|} · ∇ div E(x) − ik div E(x) =

k^{2} ∫_{∂D} ν(y) · E(y) (\frac{x}{|x|} · ∇Φ(x, y) − ikΦ(x, y)) ds(y)

−∫_{∂D}[rot E(y)×ν(y)+ν(y) div E(y)]·[(\frac{x}{|x|} · ∇_{x}) ∇_{x}Φ(x, y) − ik∇_{x}Φ(x, y)] ds(y). (2.13)

Es genügt nun zu zeigen, dass der Betrag der Integranden für große |x| von der Ordnung O(|x|^{−2}) ist. Mit der Beziehung

∇_{x}Φ(x, y) = (ik \frac{x − y}{|x − y|} − \frac{x − y}{|x − y|^{2}} ) Φ(x, y)  (2.14)

gilt für folgenden Term aus dem ersten Integral

\frac{x}{|x|} · ∇Φ(x, y) − ikΦ(x, y) =

(ik [\frac{x · (x − y)}{|x||x − y|} − 1] − \frac{x · (x − y)}{|x||x − y|^{2}} ) \frac{e^{ik|x−y|}}{4π|x − y|}.

Der Betrag des Terms in den [ ]-Klammern lässt sich durch O(|x|^{−1}) abschätzen, denn sei o.B.d.A y ∈ S^{2}, dann ist

0 ≤ 1 − \frac{x · (x − y)}{|x||x − y|} ≤ 1 − \frac{x · (x − \frac{x}{|x|})}{|x||x − \frac{x}{|x|}|}

= \frac{|x||x − \frac{x}{|x|}| − |x|^{2} + |x|}{|x||x − \frac{x}{|x|}|}

≤ \frac{2|x|}{|x||x − \frac{x}{|x|}|} = \frac{2}{|x||1 − \frac{1}{|x|}|}.  (2.15)

Also ist der Betrag des Integranden im ersten Integral in (2.13) von der Ordnung O(|x|^{−2}). Wird diesmal Formel (2.14) auf den Term im zweiten Integral in (2.13) angewendet, folgt

(\frac{x}{|x|} · ∇_{x}) ∇_{x}Φ(x, y) − ik∇_{x}Φ(x, y) =  (2.16)

(\frac{x}{|x|} · ∇) [ik \frac{x − y}{|x − y|} − \frac{x − y}{|x − y|^{2}} ] Φ(x, y) + [k^{2} \frac{x − y}{|x − y|} − ik \frac{x − y}{|x − y|^{2}} ] Φ(x, y).

Der erste Term auf der rechten Seite von (2.16) berechnet sich zu

[ik \frac{x}{|x||x − y|} − ik|x|\frac{x − y}{|x − y|^{3}} − \frac{x}{|x||x − y|^{2}} + 2|x|\frac{x − y}{|x − y|^{4}} ] \frac{e^{ik|x−y|}}{4π|x − y|}

+ (ik \frac{x − y}{|x − y|} − \frac{x − y}{|x − y|^{2}} ) (ik \frac{x · (x − y)}{|x||x − y|} − \frac{x · (x − y)}{|x||x − y|^{2}} ) \frac{e^{ik|x−y|}}{4π|x − y|}.

Der Betrag des ersten Terms im obigen Ausdruck verhält sich für große |x| wie O(|x|^{−2}). Der zweite Term enthält allerdings den Ausdruck

−k^{2} \frac{x − y}{|x − y|} \frac{x · (x − y)}{|x||x − y|} \frac{e^{ik|x−y|}}{4π|x − y|}.

Sein Betrag ist für große |x| jedoch nur von der Ordnung O(|x|^{−1}). Jedoch zusammen mit dem zweiten Term auf der rechten Seite von (2.16), dessen Betrag sich ebenfalls wie O(|x|^{−1}) verhält, bekommen wir den Ausdruck

k^{2} \frac{x − y}{|x − y|} (1 − \frac{x · (x − y)}{|x||x − y|}) \frac{e^{ik|x−y|}}{4π|x − y|}.

Der Betrag dieses Ausdrucks (siehe (2.15)) sowie der Betrag aller verbleibenden Ausdrücke lässt sich durch O(|x|^{−2}) abschätzen. Insgesamt ist somit die Sommerfeldsche Ausstrahlungsbedingung für div E erfüllt. □

Durch die getroffenen Vorbereitungen lässt sich das zentrale Ergebnis dieses Abschnitts beweisen:

Satz 2.12 Sei E ∈ C^{3}(ℝ^{3}) Lösung von

ΔE(x) + ∇ (\frac{∇n(x)}{n(x)} · E(x)) + k^{2}n(x)E(x) = 0.  (2.17)

Ist weiter n ∈ C^{3}(ℝ^{3}) und erfüllt E^{s} die Sommerfeldsche Ausstrahlungsbedingung

\frac{x}{|x|} · ∇E_{j}^{s}(x) − ikE_{j}^{s}(x) = o(|x|^{−1}),  j = 1, 2, 3,  (2.18)

gleichmäßig für alle Richtungen x/|x|, wobei wieder E = E^{i} + E^{s} mit E^{i} Lösung von (1.5) ist, so sind E und H := rot E/ik Lösung der Maxwell–Gleichungen (1.2). Außerdem erfüllen E^{s} und H^{s} = H − rot E^{i}/ik die Silver–Müller Ausstrahlungsbedingung (1.7).

Beweis: Nach Lemma 2.10 erfüllt u := div(nE)/n die Helmholtz–Gleichung

Δu(x) + k^{2}n(x)u(x) = 0.  (2.19)

Nach Voraussetzung erfüllt E^{s} die Sommerfeldsche Ausstrahlungsbedingung, so dass nach Lemma 2.8 E^{s} die Ausstrahlungsbedingung

rot E^{s}(x) × \frac{x}{|x|} + \frac{x}{|x|} div E^{s}(x) = o(|x|^{−1})

gleichmäßig für alle Richtungen x/|x| erfüllt. Deshalb erfüllt u die Sommerfeldsche Ausstrahlungsbedingung nach Lemma 2.11, denn wegen supp(n − 1) ⊂ D gilt außerhalb der Kugel D

u = div E = div E^{s}.

Folglich ist u = 0 die eindeutige Lösung der Helmholtz–Gleichung (2.19). Mit der Definition für u ist also

div E(x) = −\frac{∇n(x)}{n(x)}.

Die Lösung E von (2.17) erfüllt also

ΔE(x) − ∇ div E(x) + k^{2}n(x)E(x) = 0.

Wegen E ∈ C^{3}(ℝ^{3}) kann die Operatorgleichung rot rot = −Δ + ∇ div angewendet werden. Es folgt

rot rot E(x) − k^{2}n(x)E(x) = 0.

Schlussendlich definieren wir H := rot E/ik. Damit ist gezeigt, dass E und H die Maxwell–Gleichungen (1.2) erfüllen. Dass E^{s} und H^{s} zudem die Silver– Müller Ausstrahlungsbedingung erfüllen, zeigt Lemma 2.7, denn außerhalb von D gilt div E^{s} = 0. □

2.4 Das Anfangswertproblem für das elektrische Feld in drei Dimensionen

Wesentlich für den PBP–Algorithmus ist das Lösen des Anfangswertproblems für das elektrische Feld, wie schon in Abschnitt 1.3.2 gezeigt wurde. Dies ist eine numerisch höchst unstabile Aufgabe. Für die skalare Helmholtz–Gleichung in ℝ^{3}

Δu + k^{2}u = 0

liefert eine Fourier–Transformation nach x_{1}, x_{2}

\frac{∂^{2}}{∂x_{3}^{2}} u(ξ_{1}, ξ_{2}, x_{3}) + (k^{2} − ξ_{1}^{2} − ξ_{2}^{2})û(ξ_{1}, ξ_{2}, x_{3}) = 0.

Dies ist eine gewöhnliche Differentialgleichung, ihre allgemeine Lösung lautet

û(ξ_{1}, ξ_{2}, x_{3}) = ∑_{j=1}^{2} c_{j} e^{(−1)^{j}i\sqrt{k^{2}−ξ_{1}^{2}−ξ_{2}^{2}} x_{3}} v̂(ξ_{1}, ξ_{2}).

Die Unstabilität besteht nun darin, dass Fehler in v für ξ_{1}^{2} + ξ_{2}^{2} > k^{2} exponentiell verstärkt werden. Ziel des nächsten Abschnitts ist zu zeigen, dass die Unstabilität verschwindet, wenn die „hohen“ Frequenzen von E vernachlässigt werden. So entsteht eine gefilterte Version von E, für die in Abschnitt 2.4.2 gezeigt wird, dass sie eine gute Näherung darstellt.

2.4.1 Stabilitätsanalyse

Sei das elektrische Feld E = E^{i} + E^{s} Lösung von

ΔE(x) + ∇ (\frac{∇n(x)}{n(x)} · E(x)) + k^{2}n(x)E(x) = 0.

E^{i} soll dabei die homogene Gleichung

ΔE^{i}(x) + k^{2}E^{i}(x) = 0

erfüllen und E^{s} die Sommerfeldsche Ausstrahlungsbedingung (2.18). Das gestreute elektrische Feld E^{s} erfüllt dann

\begin{matrix}ΔE^{s}(x) + ∇ (\frac{∇n(x)}{n(x)} · E^{s}(x)) + k^{2}n(x)E^{s}(x) = R(x), \\ R(x) := −∇ (\frac{∇n(x)}{n(x)} · E^{i}(x)) − k^{2}(n(x) − 1)E^{i}(x). \end{matrix}  (2.20)

Für die weitere Darstellung erweist sich folgende Definition als sinnvoll:

S_{x_{3}} := {(x_{1}, x_{2}, x_{3}) ∈ ℝ^{3} | (x_{1}, x_{2}) ∈ ℝ^{2}}.  (2.21)

Wir betrachten hier das Anfangswertproblem für E^{s} im Halbraum x_{3} > 0. Auf der Ebene S_{0} seien die Anfangswerte

E^{s}(x_{1}, x_{2}, 0) = G_{1}(x_{1}, x_{2}), \frac{∂E^{s}(x_{1}, x_{2}, 0)}{∂x_{3}} = G_{2}(x_{1}, x_{2})

vorgegeben. Ziel ist zu zeigen, dass die Unstabilität verschwindet, wenn anstatt der Lösung E^{s} eine Näherung für E^{s} betrachtet wird: Wir definieren

Ê_{κ}^{s}(ξ_{1}, ξ_{2}, x_{3}) := \begin{cases}Ê^{s}(ξ_{1}, ξ_{2}, x_{3}) & falls |(ξ_{1}, ξ_{2})| < κ, \\ 0 & sonst \end{cases}  (2.22)

mit einer noch festzulegenden Konstanten κ. Dabei ist Ê^{s} die Fourier–Transformierte von E^{s} in den Variablen x_{1}, x_{2}:

Ê^{s}(ξ_{1}, ξ_{2}, x_{3}) := \frac{1}{2π} ∫_{ℝ^{2}} E^{s}(x_{1}, x_{2}, x_{3}) e^{−i(x_{1}ξ_{1}+x_{2}ξ_{2})} dx_{1}dx_{2}.

Bemerkung 2.13 Die Bildung von Ê^{s}(ξ_{1}, ξ_{2}, x_{3}) ist nicht im Raum der quadratintegrierbaren Funktionen L_{2}(ℝ^{2}) möglich^{1}, da E^{s} die Sommerfeldsche Ausstrahlungsbedingung erfüllt und damit automatisch

E^{s}(x) = O(|x|^{−1})

gilt, [Colton,Kress;1992] Seite 71. Die Bildung der Fourier–Transformation von E^{s} ist daher nur im distributiven Sinne zu verstehen. Das heißt Ê^{s} ist ein bzgl. x_{3} parametrisiertes lineares Funktional auf dem Raum S aller beliebig oft differenzierbaren Funktionen ϕ(x_{1}, x_{2}), die nebst jeder Ableitung für |x| → ∞ stärker als jede Potenz von 1/|x| gegen Null konvergieren. Siehe dazu §11 in [Walter;1994]. Es ist auch nicht klar, ob die Distribution Ê^{s} eine reguläre Distribution ist, das heißt Ê_{κ}^{s}(ξ_{1}, ξ_{2}, x_{3}) kann eventuell durch (2.22) nicht definiert werden. In diesem Fall könnte (2.22) ersetzt werden durch

Ê_{κ}^{s}(x_{3}) := \begin{cases}Ê^{s}(x_{3}) & für alle offenen Mengen G ⊂ \{(ξ_{1}, ξ_{2}) | |(ξ_{1}, ξ_{2})| < κ\}, \\ 0 & sonst. \end{cases}

Unabhängig davon ob Ê_{κ}^{s}(ξ_{1}, ξ_{2}, x_{3}) eine reguläre oder irreguläre Distribution ist, ist E_{κ}(x_{1}, x_{2}, x_{3}) in den Variablen x_{1}, x_{2} eine ganze holomorphe Funktion. Denn Ê_{κ}^{s} ist in den ersten beiden Variablen eine Distribution mit kompaktem Träger, Satz von Paley–Wiener [Walter;1994].

Im Hinblick auf die oben erwähnten theoretischen Aspekte, werden im Vorgriff einige noch folgende Resultate genannt. Sei x_{3} ∈ ℝ so gewählt, dass S_{x_{3}} ∩ supp(n − 1) = ∅, dann wird im Abschnitt 2.4.2 gezeigt:

• Ê^{s}(ξ_{1}, ξ_{2}, x_{3}) bzw. Ê_{κ}^{s}(ξ_{1}, ξ_{2}, x_{3}) ist eine reguläre Distribution.

• Für κ < k gilt Ê_{κ}^{s}, E_{κ}^{s} ∈ L_{2}(ℝ^{2}) bzgl. der ersten beiden Variablen.

• |Ê^{s}(ξ_{1}, ξ_{2}, x_{3})| fällt exponentiell in ξ_{1}, ξ_{2} für |(ξ_{1}, ξ_{2})| > k.

Doch zunächst wird die Stabilität des Anfangswertproblems untersucht. Dazu benötigen wir eine Erweiterung von Proposition 1 aus [Natterer;1999]:

1Elemente aus L_{2}(ℝ^{2}) werden hier in einer laxen Sprechweise als Funktionen behandelt.

Lemma 2.14 Sei V ∈ C^{2}([0, ∞], H) eine Lösung von

V′′(t) + b(t)V′(t) + A(t)V (t) = R, t > 0,  (2.23)

wobei b(t) und A(t) = A_{1}(t) + A_{2}(t) lineare beschränkte Operatoren auf dem Hilbertraum H sind, die außerdem folgende Bedingungen erfüllen:

A_{1}^{∗} = A_{1},  (2.24)
0 ≤ (V, A_{1}V) ≤ β_{1}^{2}(V, V),  (2.25)
(V, A′_{1}V) ≤ γ_{1}(V, V),  (2.26)
∥A_{2}∥ ≤ β_{2},  (2.27)
∥b∥ ≤ β.  (2.28)

Dabei sind β, β_{1}, β_{2}, γ_{1} reelle Zahlen und (· , ·) das innere Produkt in H. Mit der Definition ϕ := (V′, V′) + (V, A_{1}V) gilt dann die Abschätzung

ϕ(t) ≤ (ϕ(0) + 2t(β + γ_{1})∥V (0)∥ + ∫_{0}^{t} ∥R(τ)∥^{2} dτ) e^{(1+β_{2}+2β+T^{2}(β_{2}+γ_{1}))} (2.29)

für 0 < t ≤ T.

Beweis: Linksmultiplikation von (2.23) mit V′ führt zu

(V′, V′′) + (V′, bV′) + (V′, AV) = (V′, R).  (2.30)

Wegen (2.24) ist

\frac{1}{2} \frac{d}{dt}(V′, V′) = Re(V′, V′′),
\frac{1}{2} \frac{d}{dt}(V, A_{1}V) = Re(V′, A_{1}V) + \frac{1}{2}(V, A′_{1}V),

was zusammen mit dem Realteil von (2.30) die Darstellung

\frac{1}{2} \frac{d}{dt}(V′, V′) + \frac{1}{2} \frac{d}{dt}(V, A_{1}V)

= Re(V′, R) − Re(V′, A_{2}V) − Re(V′, bV′) + \frac{1}{2}(V, A′_{1}V)

erlaubt. Integration nach t liefert

ϕ(t) = ϕ(0) + 2 ∫_{0}^{t} [Re(V′, R) − Re(V′, A_{2}V) − Re(V′, bV′) + \frac{1}{2}(V, A′_{1}V)] dτ.

Dies ergibt zusammen mit (2.26), (2.27), (2.28) und der Cauchy–Schwarzschen Ungleichung die Abschätzung

ϕ(t) ≤ ϕ(0) + ∫_{0}^{t}[2∥V′∥∥R∥ + 2β_{2}∥V′∥∥V ∥ + 2β∥V′∥^{2} + γ_{1}∥V ∥^{2}] dτ.

Wird nun die Ungleichung ab ≤ (a^{2} + b^{2})/2 für die ersten beiden Terme unter dem Integral auf der rechten Seite obiger Gleichung benutzt, folgt

ϕ(t) ≤ ϕ(0) + ∫_{0}^{t}[∥V′∥^{2} + ∥R∥^{2} + β_{2}∥V′∥^{2} + β_{2}∥V ∥^{2} + 2β∥V′∥^{2} + γ_{1}∥V ∥^{2}] dτ.

Mit Hilfe der Abschätzung

∫_{0}^{t} ∥V (τ)∥^{2} dτ ≤ 2t∥V (0)∥^{2} + t^{2} ∫_{0}^{t} ∥V′(τ)∥^{2} dτ

des folgenden Lemmas 2.15 und ∥V′∥^{2} ≤ ϕ(t) folgt weiter aus der letzten Abschätzung für ϕ(t)

ϕ(t) ≤ ϕ(0) + 2t(β_{2} + γ_{1})∥V (0)∥^{2} + ∫_{0}^{t} ∥R∥^{2} dτ
+ (1 + β_{2} + 2β + T^{2}(β_{2} + γ_{1})) ∫_{0}^{t} ∥V′∥^{2} dτ
≤ c_{1}(t) + c_{2} ∫_{0}^{t} ϕ(τ) dτ,

für 0 < t ≤ T. Hierbei ist

c_{1}(t) := ϕ(0) + 2t(β + γ_{1})∥V (0)∥^{2} + ∫_{0}^{t} ∥R∥^{2} dτ

und

c_{2} := 1 + β_{2} + 2β + T^{2}(β_{2} + γ_{1}).

Die Behauptung des Lemmas folgt nun aus Gronwalls Ungleichung [Walter;1972] Seite 214. □

Nachzutragen ist noch das Lemma

Lemma 2.15 Sei V ∈ C^{1}([0, ∞], H), wobei H ein Hilbertraum ist. Dann gilt die Abschätzung

∫_{0}^{t} ∥V (τ)∥^{2} dτ ≤ 2t∥V (0)∥^{2} + t^{2} ∫_{0}^{t} ∥^{2}V′∥ dτ.

Beweis: Nach dem Hauptsatz der Differential- und Integralrechnung ist

V (t) = V (0) + ∫_{0}^{t} V′(τ) dτ.

Dies führt zu der Abschätzung

∥V (t)∥ ≤ ∥V (0)∥ + ∫_{0}^{t} ∥V′(τ)∥ dτ,

bzw. mit der Abschätzung (a + b)^{2} ≤ 2a^{2} + 2b^{2} zu

∥V (t)∥^{2} ≤ 2∥V (0)∥^{2} + 2 (∫_{0}^{t} ∥V′(τ)∥ dτ)^{2}
≤ 2∥V (0)∥^{2} + 2t ∫_{0}^{t} ∥V′(τ)∥^{2} dτ,

wobei im letzten Schritt noch die Cauchy–Schwarzsche Ungleichung benutzt wurde. Integration beider Seiten nach t liefert

∫_{0}^{t} ∥V (τ)∥^{2} dτ ≤ 2t∥V (0)∥^{2} + 2 ∫_{0}^{t} τ ∫_{0}^{τ} ∥V′(s)∥^{2} dsdτ.

Der letzte Term auf der rechten Seite lässt sich noch weiter abschätzen:

∫_{0}^{t} τ ∫_{0}^{τ} ∥V′(s)∥^{2} dsdτ = ∫_{0}^{t} ∥V′(s)∥^{2} ∫_{s}^{t} τ dτds

= \frac{1}{2} ∫_{0}^{t}(t^{2} − s^{2})∥V′(s)∥^{2} ds

≤ \frac{t^{2}}{2} ∫_{0}^{t} ∥V′(s)∥^{2} ds.

Damit hat man die Behauptung. □

Mit Hilfe von Lemma 2.14 lässt sich nun die Stabilität der Anfangswertaufgabe (2.20) beweisen. Der folgende Satz findet sich entsprechend für die Helmholtz–Gleichung in zwei Raumdimensionen als Theorem 2 in [Natterer;1999]. Satz 2.16 Sei E^{s} für x_{3} > 0 Lösung der Differentialgleichung (2.20) unter Beachtung der Sommerfeldschen Ausstrahlungsbedingung und der Randwerte

E^{s}(x_{1}, x_{2}, 0) = G_{1}(x_{1}, x_{2}), \frac{∂E^{s}(x_{1}, x_{2}, 0)}{∂x_{3}} = G_{2}(x_{1}, x_{2}).

Weiter sei n ∈ C^{1}(ℝ^{3}) und n = n_{1} + \frac{i}{k}n_{2} eine Aufspaltung in Real- und Imaginärteil. Es gebe reelle Konstanten \underline{M}_{1}, \overline{M}_{1}, \overline{M}_{2}, \overline{M}_{3}, \overline{M}_{4}, so dass

0 < \underline{M}_{1} ≤ n_{1} ≤ \overline{M}_{1},  |\frac{∂n_{1}}{∂x_{3}} | ≤ \overline{M}_{1},  |n_{2}| ≤ \overline{M}_{2}
|\frac{∇n}{n}| ≤ \overline{M}_{3},  |\frac{∂}{∂x_{3}} \frac{∇n}{n}| ≤ \overline{M}_{4}

gilt. Mit κ ≤ k\sqrt{\underline{M}_{1}} sei außerdem Ê_{κ}^{s}(x_{3}) ∈ L_{2}(ℝ^{2}) für alle x_{3} mit S_{x_{3}} ∩ D ≠ ∅. Dann ist für 0 < x_{3} ≤ \overline{X}_{3}

∥E_{κ}^{s}′(x_{3})∥^{2} ≤ e^{αx_{3}} [(1 + 2x_{3}(\overline{M}_{3} + k^{2}\overline{M}_{1}))∥G_{2,κ}∥^{2} + k^{2}(1 + \overline{M}_{1})∥G_{1,κ}∥^{2}

+ 2 ∥∇ (\frac{∇n}{n} · (E^{s} − E_{κ}^{s}))∥^{2} + 2k^{4}∥(n − 1)(E^{s} − E_{κ}^{s})∥^{2}],

wobei

α = 1 + (1 + \overline{X}_{3}^{2})\sqrt{2k^{2}\overline{M}_{2}^{2} + 2(\overline{M}_{4}^{2} + 2κ^{2}\overline{M}_{3}^{2})} + 2\overline{M}_{4} + \overline{X}_{3}^{2}k^{2}\overline{M}_{1}.

Beweis: Die Fourier–Transformation der Gleichung (2.20) bzgl. x_{1} und x_{2} schreibt sich als

\frac{∂^{2}}{∂x_{3}^{2}} Ê^{s}(ξ_{1}, ξ_{2}, x_{3}) + (k^{2} − ξ_{1}^{2} − ξ_{2}^{2})Ê^{s}(ξ_{1}, ξ_{2}, x_{3}) + ( iξ_{1} ; iξ_{2} ; \frac{∂}{∂x_{3}} ) (\frac{∇n}{n} · E^{s})^{∧} (ξ_{1}, ξ_{2}, x_{3})

+ k^{2}((n − 1)E^{s})^{∧}(ξ_{1}, ξ_{2}, x_{3}) dη_{1}dη_{2} = R̂(ξ_{1}, ξ_{2}, x_{3}).

Wie am Anfang des Kapitels erläutert, ist der Ausdruck Ê^{s}(ξ_{1}, ξ_{2}, x_{3}) nur im distributiven Sinne zu verstehen. Hingegen sind die Fourier–Transformationen von (n − 1)E und R im klassischen Sinne möglich. Für |(ξ_{1}, ξ_{2})| ≤ κ schreibt sich obige Gleichung unter Verwendung der Definition (2.22) für E_{κ}^{s} als

\frac{∂^{2}}{∂x_{3}^{2}} Ê_{κ}^{s}(ξ_{1}, ξ_{2}, x_{3}) + (k^{2} − ξ_{1}^{2} − ξ_{2}^{2})Ê_{κ}^{s}(ξ_{1}, ξ_{2}, x_{3})

+ \frac{1}{2π} ( iξ_{1} ; iξ_{2} ; \frac{∂}{∂x_{3}} ) ∫_{ℝ^{2}} (\frac{∇n}{n})^{∧} (ξ_{1} − η_{1}, ξ_{2} − η_{2}, x_{3}) · Ê_{κ}^{s}(η_{1}, η_{2}, x_{3}) dη_{1}dη_{2}

+ \frac{k^{2}}{2π} ∫_{ℝ^{2}}(n − 1)^{∧}(ξ_{1} − η_{1}, ξ_{2} − η_{2}, x_{3})Ê_{κ}^{s}(η_{1}, η_{2}, x_{3}) dη_{1}dη_{2}

= ε̂(ξ_{1}, ξ_{2}, x_{3}) + R̂(ξ_{1}, ξ_{2}, x_{3}), (2.31)

wobei

ε̂(ξ_{1}, ξ_{2}, x_{3}) = − ( iξ_{1} ; iξ_{2} ; \frac{∂}{∂x_{3}} ) (\frac{∇n}{n} · (E^{s} − E_{κ}^{s}))^{∧} (ξ_{1}, ξ_{2}, x_{3})

− k^{2}((n − 1)(E^{s} − E_{κ}^{s})^{∧}(ξ_{1}, ξ_{2}, x_{3}).

In den letzten beiden Formeln ist die Fourier–Transformation ausschließlich im Raum L_{2}(ℝ^{2}) ausgeführt worden. Das Lemma 2.19 im nächsten Abschnitt und die Voraussetzungen in diesem Satz garantieren nämlich Ê_{κ}^{s}(x_{3}) ∈ L_{2}(ℝ^{2}) für alle x_{3}. Die inverse Fourier–Transformation bzgl. ξ_{1} und ξ_{2} von ε̂ ergibt

ε = −∇ (\frac{∇n}{n} · (E^{s} − E_{κ}^{s})) − k^{2}(n − 1)(E^{s} − E_{κ}^{s}).

Für festes ξ_{1}, ξ_{2} kann für (2.31)

Ê_{κ}^{s}′′(x_{3}) + b(x_{3})Ê_{κ}^{s}′(x_{3}) + A(x_{3})Ê_{κ}^{s}(x_{3}) = ε̂(x_{3}) + R̂(x_{3})

geschrieben werden. Mit der Faltung ∗ bzgl. der ersten beiden Variablen ξ_{1} und ξ_{2} ist dabei

b(x_{3})Ê_{κ}^{s}′(x_{3}) = ( \begin{matrix}0 \\ 0 \\ 1 \end{matrix} ) \frac{1}{2π} ∑_{j=1}^{3} ((\frac{∇n}{n})_{j}^{∧} ∗ \frac{∂}{∂x_{3}} Ê_{κ,j}^{s}) (x_{3})

und

A(x_{3})Ê_{κ}^{s}(x_{3}) = (k^{2} − ξ_{1}^{2} − ξ_{2}^{2})Ê_{κ}^{s}(ξ_{1}, ξ_{2}, x_{3}) + \frac{k^{2}}{2π} ((n − 1)^{∧} ∗ Ê_{κ}^{s}) (x_{3})

+ ( \begin{matrix}0 \\ 0 \\ 1 \end{matrix} ) \frac{1}{2π} ∑_{j=1}^{3} (\frac{∂}{∂x_{3}} (\frac{∇n}{n})_{j}^{∧} ∗ Ê_{κ,j}^{s}) (x_{3})

+ \frac{1}{2π} ( \begin{matrix}iξ_{1} \\ iξ_{2} \\ 0 \end{matrix} ) ∑_{j=1}^{3} ((\frac{∇n}{n})_{j}^{∧} ∗ Ê_{κ,j}^{s}) (x_{3}).

Um Lemma 2.14 anzuwenden, wird der Hilbertraum

H := {Ê = (Ê_{1}, Ê_{2}, Ê_{3})^{t} | Ê_{j} = Ê_{j}(ξ_{1}, ξ_{2}, x_{3}), Ê_{j}(x_{3}) ∈ L_{2}(ℝ^{2})}

mit dem Skalarprodukt

(Ê, F̂)_{|x_{3}} := ∫_{ℝ^{2}} Ê · \overline{F̂} dξ_{1}dξ_{2}

benutzt. Die Variable x_{3} übernimmt dabei die Rolle des Parameters t aus Lemma 2.14. Weiter benötigt man

A_{1}(x_{3})Ê_{κ}^{s}(x_{3}) := (k^{2} − ξ_{1}^{2} − ξ_{2}^{2})Ê_{κ}^{s}(ξ_{1}, ξ_{2}, x_{3}) + \frac{k^{2}}{2π} ((n_{1} − 1)^{∧} ∗ Ê_{κ}^{s})(x_{3})

und

A_{2}(x_{3})Ê_{κ}^{s}(x_{3}) := \frac{ik}{2π} ((n_{2} ∗ Ê_{κ}^{s})(x_{3}) + \frac{1}{2π} ( \begin{matrix}0 \\ 0 \\ 1 \end{matrix} ) ∑_{j=1}^{3} ((\frac{∂}{∂x_{3}} \frac{∇n}{n})_{j}^{∧} ∗ Ê_{κ,j}^{s}) (x_{3})

+ \frac{1}{2π} ( \begin{matrix}iξ_{1} \\ iξ_{2} \\ 0 \end{matrix} ) ∑_{j=1}^{3} ((\frac{∇n}{n})_{j}^{∧} ∗ Ê_{κ,j}^{s}) (x_{3}).

32 KAPITEL 2. DAS DIREKTE PROBLEM Wird die Isometrie der Fourier–Transformation in L_{2} benutzt, folgt

(A_{1}Ê, F̂)_{|x_{3}} = ((k^{2} − ξ_{1}^{2} − ξ_{2}^{2})Ê, F̂) + \frac{k^{2}}{2π} ∫_{ℝ^{2}}((n_{1} − 1)^{∧} ∗ Ê) · \overline{F̂} dξ_{1}dξ_{2}
= ((k^{2} − ξ_{1}^{2} − ξ_{2}^{2})Ê, F̂) + k^{2} ∫_{ℝ^{2}}(n_{1} − 1)E · \overline{F} dx_{1}dx_{2},

was

β_{1}^{2} = k^{2}(1 + \overline{M}_{1}), γ_{1} = k^{2}\overline{M}_{1}

impliziert. Die Beschränktheit des Operators A_{2} wird durch folgende Rechnung einsichtig:

∥A_{2}∥^{2} = sup_{∥E∥=1} ∥A_{2}Ê∥^{2}

= sup_{∥E∥=1} ∫_{ℝ^{2}} |ik(n_{2}E)^{∧} + ( \begin{matrix}0 \\ 0 \\ 1 \end{matrix} ) ((\frac{∂}{∂x_{3}} \frac{∇n}{n}) · E)^{∧} + ( \begin{matrix}iξ_{1} \\ iξ_{2} \\ 0 \end{matrix} ) (\frac{∇n}{n} · E)^{∧}|^{2} dξ_{1}dξ_{2}

≤ sup_{∥E∥=1} 2 ∫_{ℝ^{2}} [|ikn_{2}E|^{2} + |(\frac{∂}{∂x_{3}} \frac{∇n}{n}) · E|^{2} + κ^{2} |\frac{∇n}{n} · E|^{2}] dx_{1}dx_{2}

≤ 2k^{2}\overline{M}_{2}^{2} + 2(\overline{M}_{4}^{2} + κ^{2}\overline{M}_{3}^{2}) sup_{∥E∥=1} ∫_{ℝ^{2}} |E|^{2} dx_{1}dx_{2}

= 2k^{2}\overline{M}_{2}^{2} + 2(\overline{M}_{4}^{2} + κ^{2}\overline{M}_{3}^{2}). Damit ist

β_{2} = \sqrt{2k^{2}\overline{M}_{2}^{2} + 2(\overline{M}_{4}^{2} + κ^{2}\overline{M}_{3}^{2})}.

Es fehlt noch der Nachweis über die Beschränktheit des Operators b:

∥b∥^{2} = sup_{∥E∥=1} ∥bÊ∥^{2}

= sup_{∥E∥=1} ∫_{ℝ^{2}} |(\frac{∇n}{n} · E)^{∧}|^{2} dξ_{1}dξ_{2}

= sup_{∥E∥=1} ∫_{ℝ^{2}} |\frac{∇n}{n} · E|^{2} dx_{1}dx_{2}

≤ \overline{M}_{3}^{2} sup_{∥E∥=1} ∫_{ℝ^{2}} |E|^{2} dx_{1}dx_{2}

= \overline{M}_{3}^{2}.

Also ist β = \overline{M}_{4}. Wird nun Lemma 2.14 angewendet, folgt

∥E_{κ}^{s}′(x_{3})||^{2} ≤ ϕ(x_{3})

≤ e^{αx_{3}} [(1 + 2x_{3}(\overline{M}_{3} + k^{2}\overline{M}_{1}))∥G_{2,κ}∥^{2} + k^{2}(1 + \overline{M}_{1})∥G_{1,κ}∥^{2}

+ 2 ∥∇ (\frac{∇n}{n} · (E^{s} − E_{κ}^{s}))∥^{2} + 2k^{4}∥(n − 1)(E^{s} − E_{κ}^{s})∥^{2} + 2∥R∥^{2}],

also die Behauptung. □

Der eben gezeigte Satz zeigt die Beschränktheit von ∥E_{κ}^{s}′(x_{3})∥. Das Anfangswertproblem für E_{κ}^{s} ist also stabil lösbar wenn ∥(E^{s} − E_{κ}^{s})(x_{3})∥ klein ist. Außerdem wäre die Eigenschaft E^{s} ≈ E_{κ}^{s} eine sinnvolle Forderung für den Vorwärtslöser. Dies alles ist Gegenstand des nächsten Abschnitts.

2.4.2 Die Helmholtz–Gleichung als Tiefpassfilter

Außerhalb der kompakten Menge supp(n − 1) ⊂ D erfüllt das gestreute elektrische Feld die Helmholtz–Gleichung

ΔE^{s} + k^{2}E^{s} = 0.

Weil E^{s} die Silver–Müller Ausstrahlungsbedingung erfüllt, gilt nach Lemma 2.7, dass die kartesischen Komponenten von E^{s} der Sommerfeldschen Ausstrahlungsbedingung genügen. Mit Lemma 2.8 erlaubt Lemma 2.9 die Darstellung

E^{s}(x) = ∫_{∂D} [E(y)\frac{∂Φ(x, y)}{∂ν(y)} − \frac{∂E(y)}{∂ν(y)}Φ(x, y)] ds(y) für x ∈ ℝ^{3} \ \overline{D}. (2.32)

Hierbei ist

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

der Helmholtz–Kern. Zur Untersuchung des Frequenzverhaltens von Ê^{s}(ξ_{1}, ξ_{2}, x_{3}) in einer Ebene S_{x_{3}} := \{(x_{1}, x_{2}, x_{3}) ∉ \overline{D} | x_{3} = const\}, genügt es wegen Gleichung (2.32) das Frequenzverhalten von Φ(x, y) bzw. von ∇_{x}Φ(x, y) auf S zu untersuchen. Der Helmholtz–Kern Φ(x, y) ist für festes y ∈ ℝ^{3} und festes x_{3} ∈ ℝ nicht in L_{1}(ℝ^{2}) oder L_{2}(ℝ^{2}). Die Fourier–Transformation für Φ(x, y) in x_{1}, x_{2} ist deshalb im distributiven Sinne zu verstehen. Man siehe dazu zum Beispiel [Walter;1994] §11. In der folgenden Rechnung wird auf zwei Formeln aus [Oberhettinger;1973] zurückgegriffen. Es ist zu beachten, dass hier die eindimensionale Fourier–Transformation im Gegensatz zur Darstellung in [Oberhettinger;1973] mit dem Vorfaktor 1/\sqrt{2π} definiert wird:

f̂(ξ) := \frac{1}{\sqrt{2π}} ∫_{ℝ} f(x) e^{−iξx} dx.

In [Oberhettinger;1973] Tabelle 1 Zeile 95 findet sich die Formel

f(x) = \begin{matrix}e & −b\sqrt{a^{2}+x^{2}} \\ √ & \overline{a^{2} + x^{2}} \end{matrix} ,
f̂(ξ) = \sqrt{\frac{2}{π}} K_{0}(a\sqrt{b^{2} + ξ^{2}}),

wobei

K_{0}(x) = ∫_{0}^{∞} cos(x sinh t) dt

die modifizierte Bessel–Funktion der Ordnung 0 ist, siehe [Abramowitz,Stegun;1965] Kapitel 9.6. Mit a^{2} = (x_{2} − y_{2})^{2} + (x_{3} − y_{3})^{2} und b = −ik ist

Φ̂(ξ_{1}, x_{2}, x_{3}, y) = \frac{1}{\sqrt{2π}4π} ∫_{ℝ} \frac{e^{ik|x−y|}}{|x − y|} e^{ikξ_{1}x_{1}} dx_{1}

= \frac{1}{\sqrt{2π}4π} e^{iy_{1}ξ_{1}} ∫_{ℝ} \begin{matrix}e & ik\sqrt{x_{1}^{2}+(x_{2}−y_{2})^{2}+(x_{3}−y_{3})^{2}} \\ √ & \overline{x_{1}^{2} + (x_{2} − y_{2})^{2} + (x_{3} − y_{3})^{2}} \end{matrix} e^{ikξ_{1}x_{1}} dx_{1}

= \frac{1}{4π^{2}} e^{iy_{1}ξ_{1}} K_{0}(\sqrt{ξ_{1}^{2} − k^{2}}\sqrt{(x_{2} − y_{2})^{2} + (x_{3} − y_{3})^{2}}).

Hierbei wurde im ersten Schritt die Integrationsvariable x_{1} durch x_{1}+y_{1} substituiert. Ähnlich wird eine weitere Formel aus [Oberhettinger;1973] Tabelle 1 Zeile 439

f(x) = K_{0}(b\sqrt{a^{2} + x^{2}}),
f̂(ξ) = \sqrt{\frac{π}{2}} \frac{1}{\sqrt{b^{2} + ξ^{2}}} e^{−a\sqrt{b^{2}+ξ^{2}}}

benutzt, um

Φ̂(ξ_{1}, ξ_{2}, x_{3}, y) = \frac{1}{8π^{2}} e^{i(y_{1}ξ_{1}+y_{2}ξ_{2})} \frac{1}{\sqrt{ξ_{1}^{2} + ξ_{2}^{2} − k^{2}}} e^{−|x_{3}−y_{3}|\sqrt{ξ_{1}^{2}+ξ_{2}^{2}−k^{2}}}  (2.33)

zu erhalten. Für |(ξ_{1}, ξ_{2})| > k fällt |Φ̂(ξ_{1}, ξ_{2}, x_{3}, y)| exponentiell in ξ_{1}, ξ_{2} ab, wie der letzte Term auf der rechten Seite von (2.33) zeigt. Bei der Untersuchung von ∇_{x}Φ(x, y) benutzt man einfach, dass

(∇_{x}Φ(x, y))^{∧} = i ( ξ_{1} ; ξ_{2} ; \frac{∂}{∂x_{3}} ) Φ̂(ξ_{1}, ξ_{2}, x_{3}, y)

gilt. Das Resultat dieses Abschnitts kann so formuliert werden:

Satz 2.17 Sei das gestreute elektrische Feld E^{s} wie in Satz 2.17 Lösung des Streuproblems und x_{3} fest gewählt, so dass die Ebene S_{x_{3}} den Träger D von n − 1 nicht schneidet. Dann fällt die Fourier–Transformierte Ê^{s}(ξ_{1}, ξ_{2}, x_{3}) von E^{s}(x_{1}, x_{2}, x_{3}) für |(ξ_{1}, ξ_{2})| > k exponentiell in ξ_{1}, ξ_{2} ab.

Bemerkung 2.18 Wie üblich wird hier dieselbe Bezeichnung für eine reguläre Distribution und ihre sie erzeugende Funktion gewählt. Die jeweilige Bedeutung wird dabei durch den Kontext festgelegt. So ist zum Beispiel im obigen Satz mit Ê^{s}(ξ_{1}, ξ_{2}, x_{3}) die erzeugende Funktion der gleichnamigen regulären Distribution gemeint.

Der Satz zeigt, dass außerhalb des Streuers in D Frequenzen (ξ_{1}, ξ_{2}) mit |(ξ_{1}, ξ_{2})| > k für das gestreute elektrische Feld E^{s} eine vernachlässigbare Rolle spielen. Während in Satz 2.16 die Stabilität des Anfangswertproblems für κ etwas kleiner als k gezeigt worden ist, gilt E^{s} ≈ E_{κ}^{s}, wenn κ etwas größer als k ist. Dieser Konflikt spielt in der Praxis keine Rolle. Für κ etwas kleiner als k zeigten die numerischen Simulationen sowohl Stabilität als auch Genauigkeit. Innerhalb des Streuers in D hängen die niedrigen Frequenzen durchaus von den höheren ab. Diese Abhängigkeit wird zum Beispiel durch die Faltung (n − 1)^{∧} ∗ Ê_{κ}^{s} in Gleichung (2.31) deutlich.

Aus der expliziten Darstellung (2.33) von Φ̂(ξ_{1}, ξ_{2}, x_{3}) lässt sich noch folgendes Lemma ableiten:

Lemma 2.19 Seien E^{s} und x_{3} wie im Satz 2.17. Dann ist für κ < k die Funktion Ê_{κ}^{s} quadratintegrierbar bzgl. ξ_{1}, ξ_{2}:

Ê_{κ}^{s}(x_{3}) ∈ L_{2}(ℝ^{2}).

Zur Definition von Ê_{κ}^{s} siehe (2.22).

Beweis: Zum Beweis ist lediglich zu sagen, dass Ê_{κ}^{s} eine stückweise stetige Funktion mit kompaktem Träger ist, denn die Singularität (ξ_{1}^{2}+ξ_{2}^{2}−k^{2})^{−1/2} in der Darstellung für Φ̂(ξ_{1}, ξ_{2}, x_{3}) in Gleichung (2.33) wird nach Definition (2.22) für Ê_{κ}^{s} abgeschnitten. □

2.4.3 Zusammenfassung

Wie schon in der Einleitung gezeigt, ist das Lösen des direkten Problems Teil des Algorithmus für das Lösen des inversen Problems. In diesem Kapitel ist deutlich geworden, wie das direkte Problem gelöst werden kann. In Abschnitt 2.3 wurde gezeigt, dass die Maxwell–Gleichungen äquivalent zu einer gestörten Form des elliptischen Δ–Operators sind. Für das gestreute elektrische Feld E^{s} lautet die Differentialgleichung

\begin{matrix}ΔE^{s}(x) + ∇ (\frac{∇n(x)}{n(x)} · E^{s}(x)) + k^{2}n(x)E^{s}(x) = R(x), \\ R(x) := −∇ (\frac{∇n(x)}{n(x)} · E^{i}(x)) − k^{2}(n(x) − 1)E^{i}(x). \end{matrix}

Diese Gleichung kann nun diskretisiert werden. Ihr Anfangswertproblem wird dann auf einem Gitter gelöst. Auf der Grundlage des Stabilitätsresultats aus Satz 2.16 wird jede neu berechnete Gitterebene mit Hilfe der Fast–Fourier–Transformation (FFT) gefiltert. Dass durch dieses Vorgehen keine zu großen Fehler zu erwarten sind, hat die Analyse des Frequenzverhaltens des gestreuten elektrischen Feldes E^{s} in Abschnitt 2.4.2 gezeigt.

Das inverse Problem

Zunächst soll das inverse Problem nochmal definiert werden:^{1} Für endlich viele verschiedene Richtungen θ_{j} ∈ S^{2}, j = 0, …, P − 1, seien folgende Informationen gegeben:

• Gegeben sei die einfallende Welle E^{i}(x) = E_{0}^{i} e^{ikθ_{j}·x} mit Einfallsrichtung θ_{j}, Amplitude E_{0}^{i} ∈ ℝ^{3}, E_{0}^{i} ⊥ θ_{j} und Wellenzahl k.

• Durch die Streuung der einfallenden Welle E^{i}(x) gibt es Messdaten E_{|Γ_{θ_{j}}}, E_{|Γ_{θ_{j}}^{+}} und E_{|Γ_{θ_{j}}^{−}} auf den Rändern des Würfels Q_{θ_{j}}. Durch Lösen des Außenraumproblems erhält man zusätzlich ∂E_{|Γ_{θ_{j}}^{−}} /∂θ_{j}.

• Außerdem gelte supp(n − 1) ⊂ D ⊂ Q_{θ_{j}}.

Gesucht ist nun der komplexwertige Brechungsindex n(x) in der Kugel D ⊂ Q_{θ_{j}}.

Das Lösen von

ΔE + ∇ (\frac{∇n}{n} · E) + k^{2}nE = 0,  (3.1)

wobei

E(x) = E^{s}(x) + E_{0}^{i} e^{ikθ·x} θ, E_{0}^{i} ∈ S^{2} mit θ ⊥ E_{0}^{i}  (3.2)

und E^{s} die Sommerfeldsche Ausstrahlungsbedingung (2.18) gleichmäßig für alle Richtungen x/|x| erfüllt, gibt Anlass zu einer nichtlinearen Abbildung

R̃_{j}(n) := E in Q_{θ_{j}}.

Wegen der Existenz und Eindeutigkeit der Lösung E ist die Abbildung R̃_{j} wohldefiniert. In diesem Abschnitt wird ein Algorithmus (PBP–Algorithmus) zur Lösung des inversen Problems entwickelt. Eine wesentliche Idee dabei ist, das direkte Problem zu linearisieren. Dies wird im folgenden Abschnitt geschehen.

1Die folgenden Bezeichnungen wurden in der Einleitung Abschnitt 1.3 und Abbildung 1.1 eingeführt.

3.1 Linearisierung

Der zentrale Begriff in diesem Abschnitt wird durch folgende Definition aus der nichtlinearen Funktionalanalysis bekannt gemacht (siehe [Werner;2000]):

Definition 3.1 Seien X und Y normierte Räume, U ⊂ X offen und f : U → Y eine Abbildung.

1. f heißt Fréchet–differenzierbar bei x_{0} ∈ U, falls ein stetiger linearer Operator T : X → Y mit

lim_{h→0} \frac{f(x_{0} + hv) − f(x_{0})}{h} = T v für alle v ∈ X

existiert, wobei die Konvergenz gleichmäßig bezüglich v ∈ {x ∈ X | ∥x∥ ≤ 1} zu sein hat. Wegen der Offenheit von U ist der Limes wohldefiniert.

2. f heißt Fréchet–differenzierbar auf U, falls f für alle x_{0} ∈ U Fréchet– differenzierbar ist.

Die Definition für die im Folgenden auftretenden Sobolevräume findet man in [Werner;2000] Kapitel 5.

Eine heuristisch hergeleitete Linearisierung von R̃_{j} ist einfach zu erhalten. Sei dazu E = R̃_{j}(n) Lösung von (3.1) mit (3.2) und Ẽ = R̃_{j}(n + d) Lösung von

ΔẼ + ∇ (\frac{∇n + d}{n + d} · Ẽ) + k^{2}(n + d)Ẽ = 0,

wobei

Ẽ(x) = Ẽ^{s}(x) + E_{0}^{i} e^{ikθ·x}

und Ẽ^{s} die Sommerfeldsche Ausstrahlungsbedingung (2.18) gleichmäßig für alle Richtungen x/|x| erfüllt. Bilden wir nun die Differenz V := Ẽ − E, so erfüllt V die Differentialgleichung

ΔV + ∇ (\frac{∇(n + d)}{n + d} · Ẽ − \frac{∇n}{n} · E) + k^{2}nV + k^{2}d(V + E) = 0.

Vernachlässigen wir die nichtlinearen Terme in der obigen Differentialgleichung, folgt

ΔV + ∇ (\frac{∇n}{n} · V ) − ∇ (\frac{d∇n}{n} · E) + ∇ (\frac{∇d}{n} · E) + k^{2}nV + k^{2}dE = 0,

wobei wir die Taylorentwicklung

\frac{1}{n + d} = \frac{1}{n} − \frac{d}{n^{2}} + O(∥d∥^{2})

benutzt haben. Es ist natürlich, wenn wir annehmen, dass wir durch Lösen obiger Differentialgleichung unter Beachtung der Sommerfeldschen Ausstrahlungsbedingung für V die Fréchet–Ableitung Ã_{j}(n) von R̃_{j} durch

Ã(n)d = V_{|∂Q_{θ_{j}}}

erhalten. Um ein solideres mathematisches Fundament zu erhalten, wird hier auf ein Resultat aus der Arbeit von [Ammari,Bao;2001] zurückgegriffen: Um die Differenzierbarkeitsanforderungen gering zu halten, beschreiben die Autoren das direkte Problem durch die Maxwell–Gleichung

rot rot E − k^{2}nE = 0  (3.3)

mit der Separation (3.2) und der Ausstrahlungsbedingung

rot E × \frac{x}{|x|} − ikE = o(|x|^{−1}).  (3.4)

Der so entstandene Lösungsoperator in Q_{θ_{j}} wird im Folgenden mit S(n) bezeichnet. Der Brechungsindex n(x) entstammt dabei dem Raum X_{M_{0}}, der Abschluss des Raums

{(f − 1) ∈ C_{0}^{2}(Q_{θ_{j}}) | supp(f) ⊂ D, ∥f∥_{∞}^{2} + ∥∇f∥_{∞}^{2} ≤ M_{0}^{2}}

im Sobolevraum H^{1}(Q_{θ_{j}}) mit der üblichen Norm ∥ · ∥_{H^{1}} ist. Sei E = S(n), dann bezeichnen wir den Lösungsoperator in Q_{θ_{j}} von

rot rot V − k^{2}nV − k^{2}dE = 0,  (3.5)
rot V × \frac{x}{|x|} − ikV = o(|x|^{−1})  (3.6)

mit T (d). Im Folgenden gelte im Zusammenhang mit der Hölder–Stetigkeit für den Exponenten 0 < α < 1. Ähnlich zu Satz 2.3 hat man

Satz 3.2 Sei n, d ∈ C^{1,α}(ℝ^{3}) und sei der Träger von n−1 und d kompakt. Dann gilt für T (d) ∈ C^{1,α}(ℝ^{3}) und der Lösung der Integralgleichung

Ṽ = −k^{2} ∫_{ℝ^{3}} Φ(x, y)(1 − n(x))Ṽ (y) dy + ∇ ∫_{ℝ^{3}} Φ(x, y)\frac{∇n(y)}{n(y)} · Ṽ (y) dy

− k^{2} ∫_{ℝ^{3}} Φ(x, y)dE dy − ∇ ∫_{ℝ^{3}} div(dE)Φ(x, y) dy (3.7)

T (d) = Ṽ. Hierbei ist wieder Φ(x, y) die Fundamentallösung der Helmholtz– Gleichung.

Beweis: Siehe [Colton,Kress;1992] Theorem 9.1 und 9.2 und [Ammari,Bao;2001] Lemma 2.1. □

Wollen wir somit die Regularität von T (d) untersuchen, können wir wieder Satz 2.4 anwenden. Es folgt analog zu Satz 2.5

Satz 3.3 Ist Ṽ ∈ C(ℝ^{3}) Lösung der Integralgleichung (3.7) und n, d ∈ C^{p,α}(ℝ^{3}) mit p ≥ 1, dann ist Ṽ ∈ C^{p,α}(ℝ^{3}) und Lösung der Differentialgleichung (3.5) mit der Ausstrahlungsbedingung (3.6).

Das folgende Resultat zeigt, unter welchen Bedingungen T (d) die Fréchet– Ableitung von S(n) ist.

Satz 3.4 Sei n, (d + 1) ∈ X_{M_{0}}. Dann gibt es eine Konstante K, die nur von M_{0}, Q_{θ_{j}} und k abhängt, so dass die folgende Abschätzung gilt:

∥S(n + d) − S(n) − T (d)∥_{H^{1}(rot,Q_{θ_{j}})} ≤ K∥d∥_{H_{0}^{3/2}(Q_{θ_{j}})}^{2}.

Hierbei ist

H^{1}(rot, Q_{θ_{j}}) = {f ∈ H^{1}(Q_{θ_{j}}) | rot f ∈ H^{1}(Q_{θ_{j}})}

mit der Norm

∥f∥_{H^{1}(rot,Q_{θ_{j}})}^{2} = ∥f∥_{H^{1}(Q_{θ_{j}})}^{2} + ∥ rot f∥_{H^{1}(Q_{θ_{j}})}^{2}.

Beweis: Siehe [Ammari,Bao;2001] Lemma 4.3. □

Mit der Standardnorm ∥ · ∥_{H_{0}^{2}} wird der Raum C_{0}^{3,α}(ℝ^{3}) zu einem normierten Raum, in dem die Abschätzung ∥ · ∥_{H_{0}^{3/2}} ≤ ∥ · ∥_{H_{0}^{2}} gilt. Das folgende Korollar ergibt sich daher aus dem bereits Gesagtem.

Korollar 3.5 Ist sogar n, (d+1) ∈ C^{3,α}(ℝ^{3}) und haben n−1 und d kompakten Träger in Q_{θ_{j}}, dann folgt

∥S(n + d) − S(n) − T (d)∥_{H^{1}(Q_{θ_{j}})} ≤ K∥d∥_{H_{0}^{2}(Q_{θ_{j}})}^{2}.

Wir wissen bereits, dass für n ∈ C^{3,α}(ℝ^{3}) die Gleichheit R̃(n) = S(n) gilt. In diesem Fall haben wir mit T und den Funktionenräumen wie oben die Fréchet– Ableitung von R̃ gegeben. Allerdings ist T kein elliptischer Operator. Es wäre daher von Vorteil wenn dieser mit der heuristisch hergeleiteten Linearisierung Ã_{j}(n) übereinstimmt, weil ihm eine elliptische Differentialgleichung zu Grunde liegt. Unter welchen Bedingungen dies stimmt, zeigen die folgenden beiden Sätze.

Satz 3.6 Seien n, d ∈ C^{3,α}(ℝ^{3}) und sei der Träger von n−1 und d kompakt. Erfüllt V die Differentialgleichung

rot rot V − k^{2}nV − k^{2}dE = 0,  (3.8)

so erfüllt V ebenfalls

ΔV + ∇ (\frac{∇n}{n} · V ) − ∇ (\frac{d∇n}{n} · E) + ∇ (\frac{∇d}{n} · E) + k^{2}nV + k^{2}dE = 0 (3.9)

Beweis: Wegen Satz 3.3 können wir V ∈ C^{3,α}(ℝ^{3}) annehmen, so dass durch Anwendung des div–Operators auf (3.8) folgt

div(nV) + div(dE) = 0.

Damit erhalten wir

\begin{aligned}div V & = −\frac{∇n}{n} · V − \frac{div(dE)}{n} \\ & = −\frac{∇n}{n} · V + \frac{d∇n}{n^{2}} · E − \frac{∇d}{n} · E \\ & = −\frac{∇n}{n} · V − (∇\frac{d}{n}) · E. \end{aligned}  (3.10)

Dies zusammen mit der Operatoridentität rot rot = −Δ + ∇ div liefert gerade (3.9). □

Die Bedingungen für die Umkehrung zeigt der nächste Satz:

Satz 3.7 Seien n, d ∈ C^{3,α}(ℝ^{3}) und sei der Träger von n−1 und d kompakt. Erfüllt V ∈ C^{3,α}(ℝ^{3}) die Differentialgleichung

ΔV + ∇ (\frac{∇n}{n} · V ) − ∇ (\frac{d∇n}{n} · E) + ∇ (\frac{∇d}{n} · E) + k^{2}nV + k^{2}dE = 0 (3.11)

und die Sommerfeldsche Ausstrahlungsbedingung

\frac{x}{|x|} · ∇V_{j}(x) − ikV_{j}(x) = o(|x|^{−1}), j = 1, 2, 3,

so erfüllt V ebenfalls

rot rot V − k^{2}nV − k^{2}dE = 0.  (3.12)

und die Ausstrahlungsbedingung (3.4).

Beweis: Sei zunächst V_{1} Lösung von

rot rot V_{1} − k^{2}nV_{1} − k^{2}dE = 0

mit der Ausstrahlungsbedingung (3.4). Nach Satz 3.6 erfüllt V_{1} auch die Differentialgleichung (3.11). Die Differenz V_{2} := V_{1} − V erfüllt dann

ΔV_{2} + ∇ (\frac{∇n}{n} · V_{2}) + k^{2}nV_{2} = 0,

denn nach Voraussetzung erfüllt V ebenfalls (3.11). Wegen Lemma 2.8 erfüllen die kartesischen Komponenten von V_{1} die Sommerfeldsche Ausstrahlungsbedingung. Mit der Ausstrahlungsbedingung für V gilt dies somit auch für V_{2}. Nach dem Beweis von Satz 2.12 ist

0 = div(nV_{2}) = div(n(V_{1} − V)).

Weil nach dem Beweis von Satz 3.6 div(nV_{1}) = − div(nE) gilt, folgt damit

div(nV) = div(nV_{1}) = − div(dE).

Die Rechnung (3.10) und die Operatoridentität rot rot = −Δ + ∇ div zeigen, dass V die Differentialgleichung (3.12) erfüllt. Dass V zusätzlich die Ausstrahlungsbedingung (3.4) erfüllt, folgt aus div V = div E = 0 außerhalb von supp(n − 1) ∪ supp d und Lemma 2.7. □

Wir haben somit das Resultat

Korollar 3.8 Seien n, d ∈ C^{3,α}(ℝ^{3}) und sei der Träger von n−1 und d kompakt. Dann existiert R(n) ∈ C^{3,α}(ℝ^{3}), so dass mit einer Konstanten K

∥R̃_{j}(n + d) − R̃_{j}(n) − Ã_{j}(n)d∥_{H^{1}(Q_{θ_{j}})} ≤ K∥d∥_{H_{0}^{2}(Q_{θ_{j}})}^{2}

gilt.

Bemerkung

Die Differentialgleichung (3.11) hat die notwendige Form, um die Resultate der Stabilitätsanalyse in Abschnitt 2.4.1 anzuwenden. Geht man davon aus, dass der Brechungsindex höchstens innerhalb der Kugel D unbekannt ist, so ist supp d ⊂ D eine vernünftige Annahme. Außerhalb von D erfüllt V also die Helmholtz–Gleichung. Somit sind auch die Resultate des Abschnitts 2.4.2 übertragbar und eine Filterung wie in Definition (2.22) auch hier sinnvoll.

3.2 Der PBP–Algorithmus

Der im Abschnitt 1.3.2 skizzierte Algorithmus zur Lösung des inversen Streuproblems wird hier in aller Ausführlichkeit dargestellt.

Gesucht ist die Lösung von

ΔE + ∇ (\frac{∇n}{n} · E) + k^{2}nE = 0  (3.13)

in Q_{j}. Im Hinblick auf die numerische Implementierung wird folgende Skalierung verwendet:

E(x) = e^{ikθ_{j}·x}(E^{s}(x) + E_{0}^{i}).  (3.14)

Hierbei ist E^{i}(x) = E_{0}^{i} e^{ikθ_{j}·x} die einfallende Welle mit Einfallsrichtung θ_{j} und Schwingungsrichtung E_{0}^{i} und für die konstanten Vektoren θ, E_{0}^{i} ∈ ℝ^{3} muss wegen div E^{i} = 0 zwingend θ_{j} ⊥ E_{0}^{i} gelten. Die einfallende Welle propagiert in dem Würfel Q_{j} von der Seite Γ_{j}^{−} nach Γ_{j}^{+}. Wenn weiter E^{s} die Sommerfeldsche Ausstrahlungsbedingung (2.18) erfüllt, existiert nach Satz 2.2 für obiges Problem eine eindeutige Lösung, die stetig im Sinne der Maximumsnorm vom einfallenden Feld E^{i} abhängt. Aus der Differentialgleichung (3.13) für E und der Skalierung (3.14) folgt die partielle Differentialgleichung für E^{s}

ΔE^{s} + 2ik(θ · ∇E_{l}^{s})_{l=1,2,3} + ∇ (\frac{∇n}{n} · (E^{s} + E_{0}^{i}))

+ ikθ (\frac{∇n}{n} · (E^{s} + E_{0}^{i})) + k^{2}(n − 1)(E^{s} + E_{0}^{i}) = 0. (3.15)

Durch die Vorgabe von Randdaten für E^{s} durch

E^{s} = G_{j}  auf Γ_{j}^{−} ∪ Γ_{j},  (3.16)
\frac{∂E^{s}}{∂θ_{j}} = \frac{∂G_{j}}{∂θ_{j}}  auf Γ_{j}^{−},  (3.17)

erhalten wir ein Anfangswertproblem. Es gibt Anlass zu einer nichtlinearen Abbildung R_{j} : H_{0}^{2}(Q_{j}) → L_{2}(Γ_{j}^{+}), R_{j}(n) = E_{|Γ_{j}^{+}}^{s}. Hierbei ist H_{0}^{2}(Q_{j}) die Gesamtheit der regulären Distributionen f ∈ D′(Q_{j}), deren Distributionsableitungen D^{α}f für den Multiindex α mit |α| ≤ 2 ebenfalls regulär sind und zum Raum L_{2,0}(Q_{j}) gehören. Der Grund für die Wahl der Funktionenräume H_{0}^{2}(Q_{j}) bzw. L_{2}(Γ_{j}^{+}) wird im Folgenden klar. Das inverse Problem besteht nun darin, den Brechungsindex n(x) so zu bestimmen, dass

R_{j}(n) = G_{j}_{|Γ_{j}^{+}}  (3.18)

für alle j = 1, …, P gilt. Ausgehend von einer Anfangsnäherung n^{0}(x) vom wahren Brechungsindex n(x) soll d^{r}(x) so bestimmt werden, dass n^{r+1} = n^{r} + ωd^{r} eine bessere Näherung für n(x) darstellt als n^{r}(x), r = 0, 1, . . .. Dazu wird Gleichung (3.18) linearisiert

R_{j}(n + d) = R_{j}(n) + A_{j}(n)d + O(∥d∥^{2}),

wobei der Operator A_{j}(n) die Fréchet–Ableitung des nichtlinearen Operators R_{j} ist. Zu lösen ist also

A_{j}(n^{r})d^{r} = G_{j}_{|Γ_{j}^{+}} − R_{j}(n^{r}) j = r mod P.  (3.19)

Dieses Gleichungssystem ist stark unterbestimmt. Ein vernünftiger Ansatz wäre, die Minimum–Norm–Lösung von (3.19) zu bestimmen. Dafür wird neben dem linearen Operator A_{j}(n) noch sein adjungierter Operator A_{j}(n)^{∗} benötigt.

Der Operator A_{j}(n) : H_{0}^{2}(Q_{j}) → L_{2}(Γ_{j}^{+})

Der Operator A_{j}(n) angewendet auf d berechnet sich durch Beschränkung der Lösung auf Γ_{j}^{+} von

ΔW + 2ik(θ · ∇W_{l})_{l=1,2,3} + ∇ (\frac{∇n}{n} · W ) + ikθ (\frac{∇n}{n} · W ) + k^{2}(n − 1)W

+ ∇ (∇\frac{d}{n} · (E^{s} + E_{0}^{i})) + ikθ (∇\frac{d}{n} · (E^{s} + E_{0}^{i})) + k^{2}d(E^{s} + E_{0}^{i}) = 0

in Q_{j} mit Randdaten W = 0 auf Γ_{j}^{−} ∪ Γ_{j} und ∂W/∂θ_{j} = 0 auf Γ_{j}^{−}, wobei E^{s} die Lösung von (3.15) mit den Randdaten (3.16) und (3.17) ist.

Obiges Ergebnis erhält man durch Verwendung von (3.11) und den Skalierungen

V (x) = e^{ikθ_{j}·x} W (x)

bzw. (3.14) für E.

Der Operator A_{j}(n)^{∗} : L_{2}(Γ_{j}^{+}) → H_{0}^{2}(Q_{j})

Zur Berechnung des adjungierten Operators A_{j}(n)^{∗} erfolgt die Skalierung am Ende der Rechnung. Sei Y Lösung von

ΔY − \frac{∇n̄}{n̄} div Y + k^{2}n̄Y = 0  (3.20)

in Q_{j} mit Randdaten Y = 0 auf Γ_{j}^{+}∪Γ_{j} und ∂Y/θ_{j} = −∂G̃_{j}/θ_{j} := − e^{−ikθ_{j}·x} ∂G_{j}/∂θ_{j} auf Γ^{+}. Dann ist mit der äußeren Normalen ν auf ∂Q_{j}

I := ∫_{Q_{j}} \{[ΔV + ∇ (\frac{∇n}{n} · V ) + k^{2}nV ] · \overline{Y} − V · \overline{[ΔY − \frac{∇n̄}{n̄} div Y + k^{2}n̄Y ]}\} dx

= ∫_{∂Q_{j}} (\frac{∂V}{∂ν} · \overline{Y} − V · \frac{∂\overline{Y}}{∂ν} ) ds (nach der zweiten Greenschen Formel)

= − ∫_{Γ_{j}^{+}} V · \frac{∂\overline{Y}}{∂θ_{j}} ds = ∫_{Γ_{j}^{+}} V · \overline{G̃}_{j} ds

= (A(n)d, G̃_{j})_{L_{2}(Γ_{j}^{+})}.  (3.21)

Andererseits berechnet sich I durch Anwendung von (3.10), (3.11) und (3.20) als

I = ∫_{Q_{j}} [−k^{2}dE − ∇ (\frac{1}{n} div(dE))] · \overline{Y} dx

= ∫_{Q_{j}} (−k^{2}dE · \overline{Y} + \frac{1}{n} div(dE) div \overline{Y} ) dx

= ∫_{Q_{j}} d (−k^{2}E · \overline{Y} − E · ∇\frac{div \overline{Y}}{n} ) dx,  (3.22)

wobei zweimal partiell integriert worden ist. Die dabei auftretenden Randterme verschwinden wegen d ∈ H_{0}^{2}(Q_{j}). Für das innere Produkt zweier Funktionen f, g ∈ H_{0}^{2}(Q_{j}) wird hier

(f, g)_{H_{0}^{2}(Q_{j})} = ∫_{Q_{j}} ΔfΔ\overline{g} dx

gewählt. Damit ist

(A_{j}(n)d, G̃_{j})_{L_{2}(Γ_{j}^{+})} = (d, A_{j}(n)^{∗}G̃_{j})_{H_{0}^{2}(Q_{j})}
= (Δd, ΔA_{j}(n)^{∗}G̃_{j})_{L_{2}(Q_{j})}
= (d, Δ^{2}A_{j}(n)^{∗}G̃_{j})_{L_{2}(Q_{j})}.

Deshalb folgt mit (3.21) und (3.22)

Δ^{2}A_{j}(n)^{∗}G̃_{j} = −k^{2}\overline{E} · Y − \overline{E} · ∇\frac{div Y}{n̄} .

Mit den Skalierungen

Y (x) = e^{ikθ_{j}·x} Z(x)

und (3.14) für E schreibt sich Δ^{2}A_{j}(n)^{∗} als

Δ^{2}A_{j}(n)^{∗}G_{j} = −k^{2}(\overline{E}^{s} + E_{0}^{i}) · Z − ik(\overline{E}^{s} · θ_{j})\frac{div Z}{n̄}

− (\overline{E}^{s} + E_{0}^{i}) · ∇\frac{div Z}{n̄} ,

(3.23) wobei Z Lösung von

ΔZ + 2ik(θ_{j} · ∇Z_{l})_{l=1,2,3} − \frac{∇n̄}{n̄} div Z − ik \frac{∇n̄}{n̄}(θ_{j} · Z) + k^{2}(n̄ − 1)Z = 0 (3.24)

mit Randdaten Z = 0 auf Γ_{j}^{+} ∪ Γ_{j} und ∂Z/θ_{j} = −∂G_{j}/∂θ_{j} auf Γ^{+} ist.

Damit A_{j}(n)^{∗}G_{j} ∈ H_{0}^{2}(Q_{j}) gilt, werden zum Lösen von (3.23) Null–Randdaten auf ∂Q_{j} vorgegeben. Doch um diesen Rechenschritt zu vermeiden, hilft die Beobachtung, dass das Lösen von Gleichung (3.23) wie ein Tiefpassfilter wirkt. Denn wendet man Δ^{2} auf die Fourier–Reihe einer in Q_{j} periodischen Funktion f : Q_{j} → ℂ an, so ist

Δ^{2}f(ξ_{1}, ξ_{2}, ξ_{3}) = ∑_{k_{1},k_{2},k_{3}∈ℕ}c_{k_{1}k_{2}k_{3}}ω_{0}^{4}(k_{1}^{4} + k_{2}^{4} + k_{3}^{4}) e^{iω_{0}(k_{1}ξ_{1}+k_{2}ξ_{2}+k_{3}ξ_{3})}

mit ω_{0} = π/ρ und

c_{k_{1}k_{2}k_{3}} = \frac{1}{(2ρ)^{3}} ∫_{Q_{j}} f(x_{1}, x_{2}, x_{3}) e^{−iω_{0}(k_{1}x_{1}+k_{2}x_{2}+k_{3}x_{3})} dx.

Oder anders gesagt, wirkt „Δ^{−2}“ als Tiefpassfilter. Bei der numerischen Implementierung wird auf dem aktualisierten Brechungsindex n^{r+1} der Tiefpassfilter

n_{κ}(x) = 1 + \frac{1}{(2π)^{3/2}} ∫_{ℝ^{3}} n̂_{κ}(ξ) e^{iξ·x} dξ

mit

n̂_{κ}(ξ) := \begin{cases}(n − 1)^{∧}(ξ) & falls |ξ| ≤ κ, \\ 0 & sonst \end{cases}

angewendet. Daher kann der Differentialoperator Δ^{2} in (3.23) als diagonaler Operator vernachlässigt werden.

Die Moore–Penrose–Inverse

Um die Minimum–Norm–Lösung von (3.19) zu bestimmen, wird noch (in laxer Schreibweise) der Operator

C_{j} = (A_{j}(n^{r})A_{j}(n^{r})^{∗})^{−1}, j = r mod P,

benötigt. Da seine numerische Berechnung sehr zeitintensiv ist, ersetzen wir ihn durch seine Werte für großes k und kleines |n|, |∇n|, nämlich C_{j} ∼ 2/(ρk^{2})I, wobei der Operator I die Identität darstellt. Das sieht man so ein: Für die folgende Rechnung seien E_{l}^{s}, E_{0,l}^{i} ≠ 0, l = 1, 2, 3. Dann vereinfacht sich die Differentialgleichung (3.15) für die l-te Komponente des gestreuten elektrischen Feldes zu

2ikθ · ∇E_{l}^{s} = k^{2}(1 − n)E_{0,l}^{i}.  (3.25)

Dies ist eine einfache Form einer Transportgleichung für die l-te Komponente des gestreuten elektrischen Feldes. Der Einfachheit halber sei nun θ_{j} = (0, 0, 1)^{t} und E_{0,l}^{i} = 1. Für x ∈ Γ_{j}^{+} (dann ist x_{3} = ρ fest) lautet Gleichung (3.25)

\frac{d}{dx_{3}} E_{j}^{s}(x) = \frac{k}{2i}(1 − n).

Integration dieser Gleichung nach x_{3} liefert

E_{j}^{s}(x) = \frac{k}{2i} ∫_{−ρ}^{x_{3}}(1 − n(x_{1}, x_{2}, τ)) dτ

= \frac{k}{2i}(P (1 − n))(θ_{j}, x).

Dabei ist P die Röntgen–Transformation (Ray Transform), siehe zum Beispiel [Natterer,Wübbeling;2001] Abschnitt 2.2. Die Röntgen–Transformation ist linear in n. Deshalb schreibt sich der lineare Operator A_{j} als

A_{j}d = \frac{k}{2i}(P (d))(θ_{j}, x).

Der dazu adjungierte^{2} Operator ist

A_{j}^{∗}G_{j}(x_{1}, x_{2}, x_{3}) = −\frac{k}{2i}G_{j}(x_{1}, x_{2}),

wobei G_{j} auf Γ_{j}^{+} definiert ist. Damit ist

A_{j}(A_{j}^{∗}G_{j})(x) = \frac{k}{2i} ∫_{−ρ}^{ρ}(A_{j}^{∗}G_{j})(x_{1}, x_{2}, τ) dτ

= \frac{k^{2}}{4} ∫_{ℝ} G_{j}(x_{1}, x_{2})) dτ

= \frac{ρk^{2}}{2} G_{j}(x_{1}, x_{2}).

Asymptotisch ist also C_{j} ∼ 2/(ρk^{2})I.

Der Algorithmus

Ausgehend von r = 0 und einer Anfangsnäherung n^{0}(x) für den Brechungsindex n(x) lässt sich der Algorithmus wie folgt beschreiben:

2Der Begriff „adjungiert“ bezieht sich hier nur auf die Ortsvariable x und nicht zusätzlich auf die Richtungsvariable θ_{j} wie sonst üblich.

1. Berechne

d^{r} = A_{j}(n^{r})^{∗}C_{j} (R_{j}(n^{r}) − G_{j}_{|Γ_{j}^{+}}) , j = r mod P.

2. Mache das Update

n^{r+1} = n^{r} + ωd^{r}.  (3.26)

3. Wenn Güte der Näherung n^{r+1} noch zu schlecht, setze r := r + 1 und springe zu 1.

Punkt 1 bedarf noch einer genaueren Erläuterung: Zunächst wird das direkte Problem durch das Berechnen des Operators R_{j}(n^{r}) gelöst (Propagation), siehe (3.15). Danach wird das Residuum durch R_{j}(n^{r}) − G_{j}_{|Γ_{j}^{+}} gebildet, das anschließend mittels C_{j} skaliert wird. Schlussendlich wendet man darauf den Operator A_{j}(n^{r})^{∗} an, siehe (3.24) und (3.23), das ist die Rückpropagation^{3}. Damit steht das Update d^{r} fest. Der Parameter ω in (3.26) muss wegen der Approximation von C_{j} und A_{j}(n)^{∗} durch numerische Experimente bestimmt werden. Der Algorithmus zur Lösung des inversen Problems unterscheidet sich bei der Verwendung von Messdaten in folgendem Punkt von Algorithmen zum Beispiel aus [Bulyshev et al;2000] oder [Colton,Kirsch;1996]. Dort werden alle Messdaten, die aus Streuexperimenten für verschiedene Richtungen θ erhalten worden sind, für ein Update verwendet. Ein Zuviel an Daten führt dort eher zu Problemen bei der Implementierung als bei der Verwendung des Algorithmus, der hier vorgestellt worden ist. Hier werden lediglich die Messdaten aus Streuexperimenten für jeweils eine Richtung für ein Update verwendet. Man beachte die Analogie zum Kaczmarz–Verfahren bei linearen Gleichungssystemen. Siehe auch [Meyn;1983].

Einzelheiten zur Diskretisierung und Filterung sind in Kapitel 4 aufgeführt.

3Damit erklärt sich auch die Namensgebung Propagation–Backpropagation Algorithmus (PBP– Algorithmus)

Numerische Simulationen

In diesem Kapitel wird die praktische Brauchbarkeit des PBP–Algorithmus an Hand von drei numerischen Simulationen belegt. Dafür sind Messdaten auf dem Rand ∂Q_{θ} nötig. Zu realen Messdaten aus Streuexperimenten mit elektromagnetischen Wellen bestand für diese Arbeit leider kein Zugang. Für die erste numerische Simulation werden die Messdaten einer analytischen Lösung des Streuexperimentes entnommen, nämlich aus der Streuung einer homogenen ebenen Welle an einer Kugel mit konstantem komplexwertigen Brechungsindex. Diese Simulation soll zeigen, dass der PBP–Algorithmus grundsätzlich in der Lage ist mit Messdaten aus realen Streuexperimenten elektromagnetischer Wellen zu arbeiten. Weiterhin dient sie der Verifikation des hier verwendeten Vorwärtslösers. Die numerischen Simulationen 2 und 3 werden demonstrieren, dass die Rekonstruktion in dem Gebiet D auch für viel komplexere Objekte als eine homogene Kugel möglich ist. In Simulation 2 wird das einfache Modell eines Herzens genommen. Die hier verwendete Geometrie dieses Modells findet sich so qualitativ auch in [Bulyshev et al;2000]. In Simulation 3 dient ein inhomogener Zylinder als Vorlage. Sein Brechungsindex weist in x_{1}–Richtung einen konstanten von Null verschiedenen Gradienten auf.

In einem weiteren Abschnitt wird der Frage nachgegangen inwieweit die Maxwell– Gleichungen durch die Helmholtz–Gleichung ersetzt werden können.

Zunächst wird allerdings die numerische Umsetzung des in 3.2 beschriebenen PBP– Algorithmus dargestellt.

4.1 Numerische Implementierung

Das zu rekonstruierende Objekt wird in der Kugel D mit dem Radius ρ = 0.5 m vermutet. Die Kugel ist enthalten in dem Quader Q_{θ}, θ ∈ S^{2}. Seine Seiten Γ_{θ}^{−} und Γ_{θ}^{+} haben einen Abstand von 1 m. Die anderen Seiten haben von ihnen gegenüber liegenden Seiten einen Abstand von 4 m. Diese Unsymmetrie hat praktische Gründe, die weiter unten erläutert werden. Der Würfel Q_{θ} wird durch ein äquidistantes Gitter der Dimension 128 × 512 × 512 diskretisiert, was einer Gitterschrittweite von h = 1/128 m ≈ 7.8 mm entspricht. Die Wellenzahl der einfallenden Welle wird auf k = 100 m^{−1} festgelegt. Dies entspricht einer Frequenz von ν ≈ 4, 77 GHz bzw. einer Wellenlänge von λ ≈ 6.3 cm. Die Differentialoperatoren in der Differentialgleichung

ΔE^{s} + 2ik(θ · ∇E_{l}^{s})_{l=1,2,3} + ∇ (\frac{∇n}{n} · (E^{s} + E_{0}^{i}))

+ ikθ (\frac{∇n}{n} · (E^{s} + E_{0}^{i})) + k^{2}(n − 1)(E^{s} + E_{0}^{i}) = 0.

aus (3.15) zur Lösung des Vorwärtsproblems werden durch Standarddifferenzenquotienten ersetzt. Man erhält so ein Differenzenverfahren, für welches das gestreute elektrische Feld in der Schicht m + 1 von den Werten des Feldes in den Schichten m und m−1 abhängt. Genauer gesagt, ist E_{j}^{s,m+1} = E_{j}^{s,m+1}(E^{s,m}, E^{s,m−1}), j = 1, 2, 3. Die Anfangswerte für E^{s,0} und E^{s,1} werden den Randdaten für E^{s}

E^{s} = G_{θ}  auf Γ_{θ}^{−} ∪ Γ_{θ},
\frac{∂E^{s}}{∂θ} = \frac{∂G_{θ}}{∂θ}  auf Γ_{θ}^{−},

entnommen, nämlich E^{s,0} = G_{θ} und E^{s,1} = h\frac{∂G_{θ}}{∂θ} + G_{θ}. Die Schichten des Würfels Q_{θ} werden also sukzessive von Γ_{θ}^{−} nach Γ_{θ}^{+} berechnet. Wegen des Stabilitätsresultates aus Satz 2.16 wird jede berechnete Schicht m gefiltert. Das heißt, man geht über von E^{s,m} nach E_{κ}^{s,m}. Dazu wendet man auf die diskreten Daten E^{s,m} eine zweidimensionale schnelle Fourier–Transformation (FFT) an und bildet

Ê_{κ}^{s,m}(ξ_{1}, ξ_{2}, x_{3}) = \begin{cases}Ê^{s,m}(ξ_{1}, ξ_{2}, x_{3}) & falls |(ξ_{1}, ξ_{2})| ≤ κ, \\ 0 & sonst. \end{cases}

4.2. SIMULATION 1: REKONSTRUKTION EINER HOMOGENEN KUGEL 51

Beachte, dass hier ξ_{1}, ξ_{2} und x_{3} diskret sind. E_{κ}^{s,m} erhält man nun durch Anwenden der inversen zweidimensionalen FFT auf Ê_{κ}^{s,m}. Weiter ist hier κ = 0.93Ωk, wobei der Faktor Ω den 2π-Raum (FFT) auf die wirkliche Geometrie anpasst (Ω = 2/π). Die FFT ist auch der Grund für den großen Abstand der Seiten Γ_{θ} voneinander. Bei dieser Distanz und den hier präsentierten Beispielen ist der Betrag des gestreuten elektrischen Feldes |E^{s}| auf Γ_{θ} sehr klein, so dass Null–Randdaten auf Γ_{θ} verwendet werden können und die FFT gut anwendbar ist (Periodizität der Daten). Der Backpropagation–Schritt wird völlig analog implementiert. Für die numerischen Simulationen sind noch folgende Parameter fixiert:

• Die „Objekte“ haben einen Durchmesser von ≈ 52 Gitterpunkten, was 40 cm entspricht.

• Die einfallende Welle ist fixiert auf E^{i}(x) = E_{0}^{i} e^{ikx_{3}}, das heißt sie propagiert in x_{3}–Richtung. In den Simulationen 1 und 2 ist ihre Amplitude E_{0}^{i} = (1, 0, 0)^{t}, in Simulation 3 gilt hingegen E_{0}^{i} = 2^{−1/2}(1, 1, 0)^{t}.

• Für jede Simulation sind Messdaten aus 120 Streuexperimenten verwendet worden. Das Objekt wurde dabei zufällig mit einem Winkel zwischen 60^{∘} und 300^{∘} in der (x_{2}, x_{3})–Ebene rotiert. Eine Verteilung der Einfallsrichtungen zeigt zum Beispiel Abbildung 4.6 auf Seite 59.

• Für alle 3 Simulationen wurde als Anfangsnäherung n^{0}(x) = 1 gewählt.

4.2 Simulation 1: Rekonstruktion einer homogenen Kugel

Gegeben sei eine im Koordinatenursprung zentrierte Kugel mit dem Radius R und dem konstanten Brechungsindex n_{K}(x) = n_{K,1}+in_{K,2} mit n_{K,1}, n_{K,2} ∈ ℝ. Außerhalb der Kugel sei der Brechungsindex konstant und reell: n(x) = n_{1}. Betrachten wir nun die Streuung der ebenen Welle E^{i}(x) = (1, 0, 0)^{t} e^{ikx_{3}} an dieser Kugel, dann finden sich die Lösungen der Maxwell–Gleichungen (1.2) in [Born,Wolf;1999] Abschnitt 14.5. Mit der Trennung E = E^{i} + E^{s} in einfallendes und gestreutes Feld schreibt sich das gestreute Feld E^{s} außerhalb der Kugel in Kugelkoordinaten als

E_{r}^{s}(r, ϑ, φ) = −\frac{1}{k^{2}n_{1}} \frac{cos φ}{r^{2}} ∑_{l=1}^{∞} l(l + 1)B_{l}^{e}ζ_{l}^{(1)}(k\sqrt{n_{1}}r)P_{l}^{(1)}(cos ϑ),
E_{ϑ}^{s}(r, ϑ, φ) = −\frac{1}{k\sqrt{n_{1}}} \frac{cos φ}{r} ∑_{l=1}^{∞} [B_{l}^{e}ζ_{l}^{(1)}′(k\sqrt{n_{1}}r)P_{l}^{(1)}′(cos ϑ) sin ϑ

− iB_{l}^{m}ζ_{l}^{(1)}(k\sqrt{n_{1}}r)P_{l}^{(1)}(cos ϑ)\frac{1}{sin φ}],

E_{φ}^{s}(r, ϑ, φ) = −\frac{1}{k\sqrt{n_{1}}} \frac{sin φ}{r} ∑_{l=1}^{∞} [B_{l}^{e}ζ_{l}^{(1)}′(k\sqrt{n_{1}}r)P_{l}^{(1)}(cos ϑ)\frac{1}{sin ϑ}

− iB_{l}^{m}ζ_{l}^{(1)}(k\sqrt{n_{1}}r)P_{l}^{(1)}′(cos ϑ) sin ϑ].

Die einzelnen Summanden in den unendlichen Summen werden auch elektrische Partialwellen genannt. Die konstanten Koeffizienten B_{l}^{e}, B_{l}^{m} charakterisieren die gestreute Welle E^{s}. Sie hängen lediglich von den Umgebungsparametern ab, nämlich von den Materialkonstanten (Radius R und Brechungsindex n_{K}), sowie von dem Brechungsindex außerhalb der Kugel n_{1} und der Wellenzahl der einfallenden Welle k. Sie berechnen sich als

B_{l}^{e} = i^{l+1} \frac{2l + 1}{l(l + 1)} \frac{k\sqrt{n_{K}}ψ′_{l}(k\sqrt{n_{1}}R)ψ_{l}(k\sqrt{n_{K}}R) − kψ′_{l}(k\sqrt{n_{K}}R)ψ_{l}(kR)}{k\sqrt{n_{K}}ζ_{l}^{(1)}′(k\sqrt{n_{1}}R)ψ_{l}(k\sqrt{n_{K}}R) − kψ′_{l}(k\sqrt{n_{K}}R)ζ_{l}^{(1)}(kR)},

B_{l}^{m} = i^{l+1} \frac{2l + 1}{l(l + 1)} \frac{k\sqrt{n_{K}}ψ_{l}(k\sqrt{n_{1}}R)ψ′_{l}(k\sqrt{n_{K}}R) − kψ′_{l}(kR)ψ_{l}(k\sqrt{n_{K}}R)}{k\sqrt{n_{K}}ζ_{l}^{(1)}(k\sqrt{n_{1}}R)ψ′_{l}(k\sqrt{n_{K}}R) − kζ_{l}^{(1)}′(kR)ψ_{l}(k\sqrt{n_{K}}R)}.

Die auftretenden unbekannten Funktionen sind

ψ_{l}(ρ) = \sqrt{\frac{πρ}{2}}J_{l+ \frac{1}{2}} (ρ), ζ_{l}^{(1)}(ρ) = \sqrt{\frac{πρ}{2}}H_{l+ \frac{1}{2}} (ρ)

mit der Bessel–Funktion J_{l+ \frac{1}{2}} (ρ) und der Hankel–Funktion erster Art H_{l+ \frac{1}{2}}^{(1)} (ρ). Weiter sind die assoziierten Legendre–Polynome erster Art definiert als

P_{l}^{(1)}(cos ϑ) = \frac{1}{sin ϑ}[P_{l−1}(cos ϑ) − cos ϑP_{l}(cos ϑ)],

wobei

P_{l}(cos ϑ) = ∑_{m=0}^{[l/2]}(−1)^{m} \frac{(2l − 2m)!}{2^{l}m!(l − m)!(l − 2m)!}(cos ϑ)^{l−2m}

die Legendre–Polynome sind. Für einen im Koordinatenursprung zentrierten Würfel Q_{θ} mit Orientierung θ ∈ S^{2} und Seitenlänge ρ > R ist damit theoretisch das elektrische Feld E auf dem Rand ∂Q_{θ} berechenbar. Für seine praktische Berechnung benötigt man noch ein Resultat aus [Debye;1909], nachdem für kR < l + 1/2 die elektrischen Partialwellen sehr schnell gegen Null konvergieren. Die oben auftreten unendlichen Summen können also in guter Näherung durch die endliche Summe der ersten [kR] + 1 Terme ersetzt werden.

4.2.1 Verifikation des Vorwärtslösers

Damit der PBP–Algorithmus auch auf reale Daten G_{θ} funktioniert, muss zwingend R_{θ}(n) ≈ G_{θ}_{|Γ_{θ}^{+}} gelten. Mit Hilfe der analytischen Lösung für die homogene Kugel ist diese Verifikation leicht möglich. Die folgenden Grafiken zeigen das gestreute elektrische Feld in der Schicht Γ_{θ}^{+}. Dabei wird jeweils die analytische Lösung mit der gefilterten analytischen Lösung und der mit dem Vorwärtslöser berechneten Lösung gegenübergestellt. Die Anfangswerte auf Γ_{θ}^{−} für den Vorwärtslöser stammen aus der analytischen Lösung. Die Darstellung beschränkt sich ausschließlich auf die Gegenüberstellung des Realteils der drei Komponenten des gestreuten Feldes, weil die Gegenüberstellung des Imaginärteils qualitativ keine neuen Informationen liefern würde. Bei dem Vergleich des Realteils der ersten Komponente des gestreuten elektrischen Feldes in Abbildung 4.1 zeigt sich eine hervorragende Übereinstimmung der Lösungen. Hier kann der Vorwärtslöser vollständig überzeugen. Bei dem Vergleich des Realteils der zweiten Komponente des gestreuten elektrischen Feldes in Abbildung 4.2 fällt auf, dass besonders zum Rand hin die mit dem Vorwärtslöser berechnete Lösung die Vorgabe durch die analytische Lösung nicht besonders gut wiedergeben kann. Der Grund dafür ist vor allem darin zu sehen, dass wegen der Null–Randbedingung auf Γ_{θ} bei der Berechnung der einzelnen Schichten in dem Quader Q_{θ} die Wellen an den Rändern reflektiert werden. Deshalb entstehen durch Superposition diese Ungenauigkeiten. Abhilfe kann hier nur ein noch größeres Rechengebiet oder die Verwendung von exakten Randdaten auf Γ_{θ} bringen. Weil bei dem Realteil der ersten (und dritten) Komponente des gestreuten elektrischen Feldes der Abfall des Feldes zum Rand hin sehr viel stärker ist als bei der zweiten Komponente, spielen reflektierte Wellen dort fast keine Rolle. Man beachte außerdem, dass der Betrag des Realteils der zweiten Komponente ca. um den Faktor 100 kleiner ist als bei dem Realteil der ersten Komponente. Daher liegt die Vermutung nahe, dass die Informationen für die Rekonstruktion des Brechungsindex, die in der zweiten Komponente des gestreuten elektrischen Feldes stecken, kein großes Gewicht haben. Der Vergleich des Realteils der dritten Komponente des gestreuten elektrischen Fel-Abbildung 4.1: Von links nach rechts zeigen die oberen Grafiken die erste Komponente des Realteils der analytischen Lösung, der gefilterten analytischen Lösung und der durch den Vorwärtslöser berechneten Lösung. Die zweite Zeile zeigt einen Querschnitt durch die oberen Grafiken und den für die oberen Grafiken verwendeten Farbkeil.

des in Abbildung 4.2 fällt ebenso positiv aus wie bei dem Vergleich der Resultate für die erste Komponente. Für alle drei Komponenten soll noch festgehalten werden, dass die analytische Lösung mit ihrer gefilterten sehr gut übereinstimmt. Dies war nach Satz 2.17, nachdem Frequenzen oberhalb von k vernachlässigbar sind, zu erwarten. Insgesamt kann festgestellt werden, dass der hier verwendete Vorwärtslöser geeignet ist, für den PBP–Algorithmus verwendet zu werden.

Abbildung 4.2: Die oberen beiden Zeilen betreffen die zweite Komponente des Realteils des gestreuten elektrischen Feldes, die beiden unteren die dritte Komponente des Realteils des gestreuten elektrischen Feldes. Von links nach rechts zeigen die oberen Grafiken die zweite Komponente des Realteils der analytischen Lösung, der gefilterten analytischen Lösung und der durch den Vorwärtslöser berechneten Lösung. Die zweite Zeile zeigt einen Querschnitt durch die oberen Grafiken und den für die oberen Grafiken verwendeten Farbkeil.

4.2.2 Rekonstruktionsergebnisse

In diesem Abschnitt wird nun die Qualität des PBP–Algorithmus für die Rekonstruktion des Brechungsindex n_{K} = 1.1 + 0.03i der Kugel mit dem Radius 20 cm gezeigt. In Abbildungen 4.3 bzw. 4.4 werden Querschnitte des Real- und Imaginärteils durch die Kugel bei x_{3} = 0 cm bzw. x_{3} = 14.0625 cm gezeigt. Da das Update d für den approximierten Brechungsindex n^{r} ebenfalls gefiltert wird, kann ein fairer Vergleich der Rekonstruktionsergebnisse nur mit der gefilterten Version des exakten Brechungsindex erfolgen. Um also eine bessere Auflösung zu erhalten, muss die Wellenzahl k vergrößert werden. Alternativ kann man natürlich Löser der partiellen Differentialgleichungen verwenden, für die eine Filterung unnötig ist.

Die Rekonstruktionsergebnisse für die Schicht x_{3} = 0 cm und x_{3} = 14.0625 cm fallen ziemlich gut aus. Vergleicht man die Rekonstruktionsqualität des Realteils mit dem Imaginärteil, wird deutlich, dass der Imaginärteil etwas schlechter rekonstruiert werden kann. Diese Beobachtung gilt in etwa auch für die Simulationen 2 und 3. Bei allen drei Simulationen konnte beobachtet werden, dass die Konvergenzgeschwindigkeit des Imaginärteils gegen ihren exakten Wert langsamer war als die des Realteils gegen seinen exakten Wert.

Abbildung 4.3: Schicht x_{3} = 0 cm: Von links nach rechts zeigen die oberen Grafiken in der ersten bzw. dritten Zeile den Realteil bzw. Imaginärteil des exakten –, des gefilterten – und des rekonstruierten Brechungsindex. Links in der zweiten und vierten Zeile ist der in den oberen Grafiken eingezeichnete Querschnitt zu sehen. Daneben befindet sich der Querschnitt durch die Mitte des Objektes.

Abbildung 4.4: Schicht x_{3} = 14.0625 cm: Von links nach rechts zeigen die oberen Grafiken in der ersten bzw. dritten Zeile den Realteil bzw. Imaginärteil des exakten – , des gefilterten – und des rekonstruierten Brechungsindex. Links in der zweiten und vierten Zeile ist der in den oberen Grafiken eingezeichnete Querschnitt zu sehen. Daneben befindet sich der Querschnitt durch die Mitte des Objektes.

4.3 Simulation 2: Rekonstruktion eines Herzmodells

In dieser Simulation soll ein Herzmodell rekonstruiert werden. Seine Geometrie ist qualitativ der Arbeit [Bulyshev et al;2000] entnommen. Hier werden andere Materialeigenschaften angenommen, was aber schon alleine wegen der stark unterschiedlichen Frequenz der einfallenden ebenen Welle so sein muss (1 GHz zu 4.77 GHz). Doch ist hier der Kontrast mit ca. 10% geringer gewählt als ca. 30% wie in der Arbeit von [Bulyshev et al;2000]. Das Herzmodell hat folgende Geometrie und Materialeigenschaften: Es besteht aus einer Kugelschale mit einem äußeren und inneren Radius von 20 cm bzw. 12 cm. Darin eingelassen ist eine zentrierte planare Trennwand der Dicke von 8 cm. Die Schale und die Trennwand werden durch den Brechungsindex von n_{H} = 1.1 + 0.05i beschrieben. Die beiden Hohlräume sind hingegen charakterisiert durch den Brechungsindex von n_{I} = 1.04 + 0.02i. Die Grafik 4.5 zeigt einen Querschnitt senkrecht zur Trennwand durch die Mitte des Herzmodells. Die Richtungen θ_{j} der einfallenden

Abbildung 4.5: Querschnitt durch das Herzmodell in der (x_{1}, x_{2})–Ebene

elektrischen Felder E^{i} liegen alle senkrecht zur Ebene der Trennwand des Herzmodells. Sie werden dort zufällig variiert. Ihre Verteilung zeigt Grafik 4.6.

[figure: OBJEKT]

Abbildung 4.6: Verteilung der Einstrahlungsrichtungen

Weil hier kein Zugang zu analytischen- oder experimentellen Daten besteht, werden auf Γ_{θ}^{−} Null–Randdaten (keine Rückstreuung) und auf Γ_{θ}^{+} mit dem Vorwärtslöser berechnete Randdaten angenommen. Das Verwenden desselben Vorwärtslöser innerhalb des Rekonstruktionsalgorithmus und zur Erzeugung von Messdaten wird als Inverse-Crime bezeichnet. Diese numerische Simulation kann daher nicht alleine für die Verifikation des PBP–Algorithmus herangezogen werden. Auf die so erzeugten Randdaten wird deshalb ein relativer Fehler von 20% addiert. Das heißt, für jeden Messwert g wird eine Zufallszahl p_{g} ∈ [−1, 1] bestimmt, so dass ein modifizierter Messwert

g̃ = g(1 + 0.2p_{g})

entsteht. Die Zufallszahl p_{g} ist dabei gleichverteilt über [−1, 1].

Als nächstes werden die Rekonstruktionsergebnisse vorgestellt. Verglichen wer-

Abbildung 4.7: Schicht x_{3} = −12.5 cm: Von links nach rechts zeigen die oberen Grafiken den Realteil des exakten –, des gefilterten – und des rekonstruierten Brechungsindex. Links unten ist der in den oberen Grafiken eingezeichnete Querschnitt zu sehen. Daneben befindet sich der Querschnitt durch die Mitte des Objektes.

den die Schichten für x_{3} = −12.5 cm und x_{3} = 0 cm. Die x_{3}–Achse liegt in der Ebene der Trennwand und das Zentrum des Herzmodells ist im Ursprung platziert. Der Vergleich findet jeweils zwischen dem exakten Brechungsindex, seiner gefilterten Version und dem mit dem PBP–Algorithmus rekonstruierten Brechungsindex statt. Theoretisch kann der PBP–Algorithmus nur so gut sein wie der gefilterte exakte Brechungsindex. Die Abbildung 4.7 zeigt den Vergleich des Realteils des Brechungsindex in der Schicht x_{3} = −12.5 cm. Sie macht deutlich, dass der Realteil des Brechungsindex des Herzmodells in der Schicht x_{3} = −12.5 cm sehr gut durch den PBP–Algorithmus rekonstruiert werden kann. Die nächste Abbildung 4.8 zeigt den Vergleich des Imaginärteils des Brechungsindex in der Schicht x_{3} = −12.5 cm. Ihr ist zu entnehmen, dass das Niveau des Ima-

Abbildung 4.8: Schicht x_{3} = −12.5 cm: Von links nach rechts zeigen die oberen Grafiken den Imaginärteil des exakten –, des gefilterten – und des rekonstruierten Brechungsindex. Links unten ist der in den oberen Grafiken eingezeichnete Querschnitt zu sehen. Daneben befindet sich der Querschnitt durch die Mitte des Objektes.

ginärteils des exakten Brechungsindex zum Rand hin nicht erreicht werden kann. Hingegen ist das Niveau der gefilterten Version zu hoch. Der Querschnitt durch die Mitte zeigt jedoch eine sehr gute Approximation. Für den Imaginärteil hat der PBP–Algorithmus in der Schicht x_{3} = −12.5 cm insgesamt eine gute Rekonstruktion geliefert. Als nächstes wird in Abbildung 4.9 die Qualität der Rekonstruktion des Brechungsindex in der Schicht x_{3} = 0 cm gezeigt, wobei hier die Werte für seinen Realteil und Imaginärteil in einer Abbildung gegenübergestellt werden. Dort wird deutlich, dass der Realteil des Brechungsindex in der Schicht x_{3} = 0 cm sehr gut approximiert wird. Die Qualität der Rekonstruktion übertrifft sogar die Qualität der gefilterten Version des exakten Brechungsindex. Die Rekonstruktion des Imaginärteils des Brechungsindex ist hingegen nicht so gut gelungen. Der erste Querschnitt durch die untere Kammer des Herzmodells zeigt, dass das Niveau von 0.02i für die Kammer und 0.05i für die Wände nach 120 Iterationen des PBP–Algorithmus nicht erreicht werden konnte. Der zweite Querschnitt durch die Mitte des Objektes zeigt hingegen noch eine befriedigende Approximation.

Abbildung 4.9: Schicht x_{3} = 0 cm: Von links nach rechts zeigen die oberen Grafiken in der ersten bzw. dritten Zeile den Realteil bzw. Imaginärteil des exakten –, des gefilterten – und des rekonstruierten Brechungsindex. Links in der zweiten und vierten Zeile ist der in den oberen Grafiken eingezeichnete Querschnitt zu sehen. Daneben befindet sich der Querschnitt durch die Mitte des Objektes.

4.3.1 Konvergenz des Vorwärtslösers

Mit jeder besseren Näherung n^{r} für den Brechungsindex erreicht der Vorwärtslöser mit R_{θ}(n^{r}) eine bessere Übereinstimmung mit den vorgegebenen Randdaten bei Γ_{θ}^{+}. Dies soll hier beim Herzmodell verdeutlicht werden. Abbildung 4.10 zeigt diese Entwicklung für den Realteil der dritten Komponente des elektrischen Feldes.

Abbildung 4.10: Die erste Grafik links oben zeigt, die verwendeten Messwerte. Die anderen Grafiken zeigen ihre Approximation nach 2, 5, 10, 50 und 120 Iterationen des PBP–Algorithmus. Gezeigt ist jeweils der Realteil der dritten Komponente des gestreuten elektrischen Feldes.

4.4 Simulation 3: Der inhomogene Zylinder

In diesem Abschnitt werden die Rekonstruktionsergebnisse für einen Zylinder mit dem Brechungsindex

n_{Z}(x) = 1.1 − \frac{32}{75}|x_{1}| + 0.05i für |x_{1}|, |(x_{2}, x_{3})| ≤ 0.2

vorgestellt. Während der Imaginärteil des Brechungsindex konstant ist, besitzt der Realteil einen konstanten Gradienten in x_{1}-Richtung. Abhängig von x_{1} verändert sich der Realteil von 1.1 im Zentrum bis ca. 1.015 am Deckel des Zylinders. Grafik 4.11 zeigt einen Querschnitt in der (x_{1}, x_{2})–Ebene durch die Mitte des Zylinders.

Abbildung 4.11: Querschnitt durch die Rotationsachse des Zylinders. Links ist der Realteil und rechts der Imaginärteil zu sehen.

Im Gegensatz zu den Simulationen 1 und 2 wird hier die Amplitude der einfallenden Welle nicht parallel zu einer kartesischen Koordinatenachse gewählt: E_{0}^{i} = 2^{−1/2}(1, 1, 0)^{t}.

Abbildung 4.12 bzw. 4.13 zeigt die Rekonstruktionsergebnisse jeweils für den Realund Imaginärteil des Brechungsindex in den Schichten x_{3} = 0 cm bzw. x_{3} = 18.75 cm. Zu den Rekonstruktionsergebnissen ist noch anzumerken, dass die Rekonstruktion des Realteils des Brechungsindex in der Schicht 0 cm dem exakten Brechungsindex etwas näher kommt als die gefilterte Version des exakten Brechungsindex. In dieser Schicht konnte auch der Imaginärteil sehr gut wiedergegeben werden. Die Rekonstruktionsergebnisse in der Schicht x_{3} = 18.75 cm sind etwas schwächer, was bei einem Abstand von 1.25 cm vom Rand und einer Wellenlänge von ≈ 6.3 cm nicht verwundert.

Abbildung 4.12: Schicht x_{3} = 0 cm: Von links nach rechts zeigen die oberen Grafiken in der ersten bzw. dritten Zeile den Realteil bzw. Imaginärteil des exakten –, des gefilterten – und des rekonstruierten Brechungsindex. Links in der zweiten und vierten Zeile ist der in den oberen Grafiken eingezeichnete Querschnitt zu sehen. Daneben befindet sich der Querschnitt durch die Mitte des Objektes.

Abbildung 4.13: Schicht x_{3} = 18.75 cm: Von links nach rechts zeigen die oberen Grafiken in der ersten bzw. dritten Zeile den Realteil bzw. Imaginärteil des exakten – , des gefilterten – und des rekonstruierten Brechungsindex. Links in der zweiten und vierten Zeile ist der in den oberen Grafiken eingezeichnete Querschnitt zu sehen. Daneben befindet sich der Querschnitt durch die Mitte des Objektes.

4.5 Bornsche Näherung

Wir wollen noch ein Gefühl dafür bekommen, wie weit entfernt wir im Fall der homogenen Kugel von der Bornschen Näherung sind. Unter der Bornschen Näherung verstehen wir das elektrische Feld E, das wir erhalten, wenn in der Integralgleichung (2.1) auf der rechten Seite E durch das einfallende Feld E^{i} ersetzt wird. Also:

E(x) = E^{i}(x) − k^{2} ∫_{ℝ^{3}} Φ(x, y)(1 − n(x))E^{i}(y) dy

+ ∇ ∫_{ℝ^{3}} \frac{∇n(y)}{n(y)} · E^{i}(y)Φ(x, y) dy.

Wenn der Term ∇n/n klein ist, ist das gestreute elektrische Feld E^{s} = E − E^{i} also „fast“ linear in n − 1. Die Grafiken in Abbildung 4.14 zeigen, dass sich das gestreute elektrische Feld in

Abbildung 4.14: Die erste Zeile zeigt den Realteil der ersten Komponente des gestreuten elektrischen Feldes bei Γ_{θ}^{+} für die Streuung an einer Kugel nach Simulation 1 für den Brechungsindex n = 1.05, 1.1, 1.15, 1.2, 1.25.

Abhängigkeit vom Brechungsindex sehr früh nichtlinear verhält, so dass für n ≈ 1.1 der Gültigkeitsbereichs der Bornschen Näherung verlassen wird.

4.6 Approximation durch die Helmholtz–Gleichung

In diesem Abschnitt wird der Frage nachgegangen, inwieweit die Maxwell–Gleichungen für die hier vorgestellten Simulationen durch die Helmholtz–Gleichung approximiert werden kann. Das heißt, ob der Term (die Kopplung)

∇ (\frac{∇n}{n} · E)  (4.1)

in der Differentialgleichung

ΔE(x) + ∇ (\frac{∇n(x)}{n(x)} · E(x)) + k^{2}n(x)E(x) = 0  (4.2)

vernachlässigt werden kann. Nehmen wir dazu an, dass die einfallende ebene Welle linear polarisiert ist:

E_{0}^{i}(x) = ( \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ) e^{ikx_{3}} .

Mit den obigen Notationen ist also hier E_{0}^{i} = (1, 0, 0)^{t} und θ_{j} = (0, 0, 1)^{t}. Weil nun die zweite und dritte Komponente des gestreuten elektrischen Feldes E^{s} erst durch die Kopplung (4.1) entstehen, können diese nicht durch die Helmholtz–Gleichung approximiert werden. Wenn also die Approximation von (4.2) durch die Helmholtz–Gleichung gelingen soll, darf nur die erste Komponente von E^{s} berücksichtigt werden. Die Differentialgleichung (3.15) für E^{s} mit E = exp(ikx_{3})((1, 0, 0)^{t} + E^{s}) schreibt sich dann als

ΔE_{1}^{s} + 2ik \frac{∂}{∂x_{3}} E_{1}^{s} + k^{2}(n − 1)(E_{1}^{s} + 1) = 0.  (4.3)

Die Operatoren R_{j}(n), A_{j}(n) und A_{j}(n)^{∗} werden genauso wie in Kapitel 3 konstruiert. Im Fall der Helmholtz–Gleichung haben sie allerdings eine viel einfachere Struktur. Nach dem Artikel [Natterer,Wübbeling;1995] berechnet sich zum Beispiel A_{j}(n)^{∗} durch

A_{j}(n)^{∗}G_{j,1} = −k^{2}(\overline{E}_{1}^{s} + 1)Z,

wobei die skalare Funktion Z Lösung von

ΔZ + 2ik \frac{∂}{∂x_{3}} Z + k^{2}(n̄ − 1)Z = 0

mir Randdaten Z = 0 auf Γ_{j}^{+} ∪ Γ_{j} und ∂Z/∂x_{3} = −∂G_{j,1}/∂x_{3} auf Γ_{j}^{+} ist. Der PBP–Algorithmus vereinfacht sich also erheblich. Bei der Implementierung drittelt sich der Speicherbedarf und die Rechendauer verkürzt sich sogar um mehr als den Faktor 3.

Um das einfachere Helmholtz–Modell mit dem Maxwell–Modell zu vergleichen, werden die Simulationen 1, 2 und 3 herangezogen. Für beide Modelle sind natürlich dieselben Messdaten auf Γ_{j}^{−} und Γ_{j}^{+} verwendet worden. Bis auf den Austausch der Modelle wurden alle übrigen Parameter beibehalten. Als erstes Beispiel greifen wir die Simulation 1 nochmals auf. Bei der homogenen Kugel muss die Kopplung (4.1) bei der Implementierung nur am Rand der Kugel berücksichtigt werden. Vermutlich wird sie daher keinen großen Einfluss auf die Lösung haben. In Abbildung 4.15 sind die Rekonstruktionsergebnisse durch die Mitte der Kugel gegenübergestellt. Es kann festgestellt werden, dass die Approximation des elektromagnetischen Streuproblems durch die Helmholtz–Gleichung in (4.3) qualitativ zu den gleichen Rekonstruktionsergebnissen führt. Die Beschreibung des Streuproblems durch die Maxwell–Gleichungen bringt somit hier überhaupt keinen Vorteil. Als nächstes ziehen wir Simulation 2 heran. Bei dem Vergleich der Rekonstruktionsqualität für den Realteil in Abbildung 4.16 konnte beobachtet werden, dass die harten Übergänge (Kanten) in andere Bereiche des Modells durch das Maxwell– Modell etwas besser rekonstruiert wird als beim Helmholtz-Modell. Allerdings zeigt der Querschnitt durch die Mitte des Objektes, dass das Helmholtz–Modell das Niveau besser trifft. Die Qualität der Rekonstruktionen des Imaginärteils des Brechungsindex ist gleichwertig. Als drittes Beispiel benutzen wir die Simulation 3. Abbildung 4.17 zeigt, dass es kleine qualitative Unterschiede in der Rekonstruktion gibt. Der erste Querschnitt in Zeile zwei ist insoweit repräsentativ, als dass er zeigt, dass die Rekonstruktion mit dem Helmholtz–Modell fast immer über dem tatsächlichen Niveau der Vorlage liegt. Hier hat das Maxwell–Modell kleine Vorteile. Auch die Rekonstruktion des Imaginärteils gelingt bei Verwendung des Maxwell-Modells ein wenig besser. Hier kann also eine kleine Verbesserung durch die Verwendung der exakteren Maxwell–Gleichung erreicht werden. Inwieweit dies in der Praxis relevant ist, kann hier nicht entschieden werden.

Abbildung 4.15: Schicht x_{3} = 0 cm: Von links nach rechts zeigen die oberen Grafiken in der ersten bzw. dritten Zeile den Realteil bzw. Imaginärteil des exakten –, des mit dem Maxwell-Modell rekonstruierten – und des mit dem Helmholtz–Modell rekonstruierten Brechungsindex. Links in der zweiten und vierten Zeile ist der in den oberen Grafiken eingezeichnete Querschnitt zu sehen. Daneben befindet sich der Querschnitt durch die Mitte des Objektes.

Abbildung 4.16: Schicht x_{3} = 0 cm: Von links nach rechts zeigen die oberen Grafiken in der ersten bzw. dritten Zeile den Realteil bzw. Imaginärteil des exakten –, des mit dem Maxwell-Modell rekonstruierten – und des mit dem Helmholtz–Modell rekonstruierten Brechungsindex. Links in der zweiten und vierten Zeile ist der in den oberen Grafiken eingezeichnete Querschnitt zu sehen. Daneben befindet sich der Querschnitt durch die Mitte des Objektes.

Abbildung 4.17: Schicht x_{3} = 0 cm: Von links nach rechts zeigen die oberen Grafiken in der ersten bzw. dritten Zeile den Realteil bzw. Imaginärteil des exakten –, des mit dem Maxwell-Modell rekonstruierten – und des mit dem Helmholtz–Modell rekonstruierten Brechungsindex. Links in der zweiten und vierten Zeile ist der in den oberen Grafiken eingezeichnete Querschnitt zu sehen. Daneben befindet sich der Querschnitt durch die Mitte des Objektes.

4.7 Abschließende Bemerkung

Durch die Arbeiten [Natterer,Wübbeling;1995] und [Natterer;1999] wurde demonstriert, dass der PBP–Algorithmus exzellent im Bereich der Ultraschall–Tomographie funktioniert. In der hier vorliegenden Arbeit konnte gezeigt werden, dass sich sowohl die Implementierung als auch theoretische Resultate auf den Bereich der Mikrowellen– Tomographie übertragen lassen. Der PBP–Algorithmus ist in der Lage das nichtlineare inverse Problem für die Maxwell–Gleichungen zu lösen. Ein Vorteil der hier gewählten Implementierung ist, dass sie sich sehr leicht auf ein Multiprozessorsystem übertragen lässt, da Rechnungen im 3D–Raum auf aktuellen Einzelprozessorsystemen sehr zeitintensiv sind. Es sollte nicht unerwähnt bleiben, dass ein Konvergenzbeweis für den PBP–Algorithmus fehlt. Allerdings wurde für den diskreten Fall unter bestimmten Voraussetzungen lokale Konvergenz in [Meyn;1983] gezeigt. Der Abschnitt über die Approximation des Maxwell–Modells durch das Helmholtz–Modell zeigt an, dass es durchaus sinnvoll sein kann, das genauere Maxwell–Modell zu nehmen. Allerdings wird für viele Brechungsindizes das Helmholtz–Modell ein guter Vorkonditionierer sein, der für simple Fälle sogar ausreichend sein wird. Doch eine abschließende Antwort auf die Frage, für welche Objekte die Verwendung des Helmholtz–Modells im Bereich der Mikrowellen–Tomographie eine ausreichende Näherung darstellt, kann nur die Simulation mit realen Daten geben.

Literaturverzeichnis

[Abramowitz,Stegun;1965] Abramowitz, M.; Stegun, I.A.: Handbook of Mathematical Functions, Dover Publications, Inc., New York, 1965

[Ammari,Bao;2001] Ammari, H.; Bao, G.: Analysis of the scattering map of all linearized medium problem for electromagnetic waves; Inverse Problems 17 (2001) 219 − 234

[Born,Wolf;1999] Born, M.; Wolf, E.: Principles of Optics, Camebridge University Press, 1999, 7th edition

[Bulyshev et al;2000] Bulyshev, A.E. et al: Three–dimensional microwave tomography. Theory and computer experiments in scalar approximation, Inverse Problems, 16, 2000, 863-85

[Colton,Kress;1983] Colton, D.; Kress, R.: Integral Equation Methods in Scattering Theory, John Wiley & Sons, 1983

[Colton,Kress;1992] Colton, D.; Kress, R.: Inverse Acoustic and Electromagnetic Scattering Theory, Springer-Verlag, 1992

[Colton,Kirsch;1996] Colton, D., Kirsch, A.: A simple method for solving inverse scattering problems in the resonance region, Inverse Problems, 383-395, 1996

[Debye;1909] Debye, P.: Ann. d. Physik(4), 30, 1909; Mat. Ann., 67, 1909; Sitzungsb. Münch. Akad. Wiss, Math. Phys. Kl., 5 Abh., 1910

[Gilbarg,Trudinger;1977] Gilbarg, D.; Trudinger, N.S.: Elliptic Partial Differential Equations of Second Order, Springer–Verlag, 1977

[Hagnes et al;1998] Hagness, S.C.; Taflove, A.; Bridges, J.E.: Two–dimensional FDTD analysis of a pulsed microwave confocal system for breast cancer detection: fixed–focus and antenna–array sensors, IEEE Trans. Biomed. Eng. 45 1470-9, 1998 [Hiraoka et al;1996] Hiraoka, T.; Tanaka, M.; Takeneka, T.; Harada, H.; Wall, D.;

Furuoku, T.: Electromagnetic Detection of Cylindrical Buried Objects, Proceedings of The Fourth Japanese and Polish Joint Seminar on Electromagnetic Phenomena Applied to Technology, June 5-7, 1995, Oita, Japan

[Hörmander;1976] Hörmander, L.: Linear Partial Differential Operators, Springer– Verlag, 1976

[Joines et al;1980] Joines, W.T.; Jirtle, R.L.; Rafal, M.D.; Schaefer, D.J.:Microwave power absorption differences between normal and malignant tissue, Int. J. Radiat. Oncol. Biol. Phys. 6 681-7, 1980

[Kress;2000] Kress, R.: Iterative solution of the inverse scattering problem for a locally pertrubed half-plane, Kleinheubacher Berichte, 2000

[Kress,Rundell;2000] Kress, R.; Rundell, W.: Inverse Scattering for Shape and Impedance, erscheint in Inverse Problems, http://www.num.math.unigoettingen.de/kress/researchlist.html

[Meyn;1983] Meyn, K.-H.: Solution of Underdetermined Nonlinear Equations by Stationary Iteration Methods, Numer. Math. 42, 161-172, 1983

[Magnus et al;1966] Magnus, W.; Oberhettinger, F.; Soni, R.P.: Formulas and Theorems for Special Functions of Mathematical Physics, Springer–Verlag Berlin Heidelberg New York, 1966

[Natterer,Wübbeling;1995] Natterer, F.; Wübbeling, F.: A propagationbackpropagation method for ultrasound tomography, Inverse Problems 11 (1995), 1225-1232

[Natterer;1999] Natterer, F.: An initial value approach to the inverse Helmholtz problem at fixed frequency, Inverse Problems in Medical Imaging and Nondestructive Testing, Proceedings of the Congerence in Oberwolfach, Germany, 4.–10. Feb. 1997

[Natterer,Wübbeling;2001] Natterer, F.; Wübbeling, F.: Mathematical Methods in Image Reconstruction, Society for Industrial and Applied Mathematics (SIAM), 2001

[Oberhettinger;1973] Oberhettinger, F.: Fourier Transforms of Distributions and Their Inverses, A Collection of Tables; Academic Press, New York and London, 1973 [Schepps,Foster;1980] Schepps, J.L.; Foster, K.R.: The UHF and microwave dielec-

trical properties of normal and tumor tissues: variation in dielectrical properties with tissues water content, Phys. Med. Biol. 25 1149-59, 1980

[Semenov et al;1996] Semenov, S.Y.; Svenson R.H.; Dezern, K.R.; Quinn M.E.; Thompson M.; Tatsis G.P.: Myocardial ischemia and infarction can be detected by microwave spectroscopy, Experimental evidence, Proc. 18th Annual Int. Conf. of the IEEE Eng. in Medicine and Biology Soc., 1996

[Semenov et al;1999] Semenov, S.Y.; Svenson, R.H.; Bulyshev, A.E.; Souvorov, A.E.; et al: Three–Dimensional Microwave Tomography: Experimental Prototype of the System and Vector Born Reconstruction Method, IEEE Transactions on Biomedical Engineering, Vol. 46, 937–946, 1999

[Walter;1972] Walter, W.: Gewöhnliche Differentialgleichungen, Springer-Verlag, 1972

[Walter;1994] Walter, W.: Einführung in die Theorie der Distributionen; BI Wissenschaftsverlag, Mannheim, Leipzig, Wien, Zürich; 1994

[Werner;2000] Werner, D.: Funktionalanalysis, Springer, 2000

Lebenslauf

Matthias Vögeler geboren am 01.03.1974 in Beelen Familienstand: ledig Eltern: Helmut Vögeler, Hildegard Vögeler (geb. Kleine–Dütting)

1980–1984  Grundschule Beelen
1984–1986  Von–Galen–Hauptschule Beelen
1986–1990  Von–Galen–Realschule Warendorf

1990–1993 Höhere Berufsfachschule mit gymnasialer Oberstufe für Wirtschaft und Verwaltung

14.06.1993  Datum des Abiturzeugnisses
07.1993–10.1994  Zivildienst im Josephs–Hospital Warendorf

WS 1994/95 – SoSe 1999 Diplomstudiengang Mathematik (Nebenfach Physik) an der Westfälischen–Wilhelms Universität Münster

03.1997–03.1998 Studentische Hilfskraft bei Baetge & Partner, Technologiehof Münster

WS 1998/99 Studentische Hilfskraft am Institut für Numerische und instrumentelle Mathematik der WWU Münster

19.04.1999  Datum des Diplomzeugnisses

Seit Mai 1999 Beginn der Promotion: Wissenschaftlicher Mitarbeiter am Institut für Numerische und instrumentelle Mathematik der WWU Münster

Münster, Oktober 2001