Numerische Verfahren

zur

Impedanz-Tomographie

Diplomarbeit

vorgelegt von

Stefan Biedenstein

Thema gestellt von

Prof. Dr. F. Natterer

22. April 1996

�

Westfälische Wilhelms-Universität Münster

Institut für Numerische Mathematik

Numerische Verfahren

zur

Impedanz-Tomographie

Diplomarbeit

vorgelegt von

Stefan Biedenstein

Thema gestellt von

Prof. Dr. F. Natterer

22. April 1996

�

Westfälische Wilhelms-Universität Münster

Institut für Numerische Mathematik

Inhaltsverzeichnis

1 Das Problem 4 1.1 Inverse Probleme . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Physikalische Gesetzmäßigkeiten . . . . . . . . . . . . . . . . . . . 5 1.3 Eindeutige Lösbarkeit der Differentialgleichung . . . . . . . . . . . 5 1.4 Elektrische-Impedanz-Tomographie . . . . . . . . . . . . . . . . . 7

2 Ein analytisches Verfahren 9

2.1 Inverse Scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.1 Vereinbarungen . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.2 Umformulierung des Problems . . . . . . . . . . . . . . . . 10

2.2 Hilfsprobleme und Fundamentallösungen . . . . . . . . . . . . . . 13 2.3 Beweis von Satz 2.1.9 . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.4 Berechnung der Scattering-Amplitude . . . . . . . . . . . . . . . . 19 2.5 Die Randwerte von σ . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.6 Die eindeutige Lösung des inversen Problems . . . . . . . . . . . . 24

3 Numerische Verfahren 25 3.1 Übersicht . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Lineare Rückprojektion . . . . . . . . . . . . . . . . . . . . . . . . 26

3.2.1 Die verallgemeinerte Radontransformation . . . . . . . . . 26 3.2.2 Linearisierung . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2.3 Die Poincaré-Abbildung . . . . . . . . . . . . . . . . . . . 29

3.2.4 Anwendung der Radontransformation . . . . . . . . . . . . 29 3.3 Variationsrechnung . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.3.1 Formulierung der Minimierungsaufgabe . . . . . . . . . . . 33 3.3.2 Numerische Methoden . . . . . . . . . . . . . . . . . . . . 33

3.4 Layer Stripping . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.4.1 Rekonstruktion auf dem Rand . . . . . . . . . . . . . . . . 37 3.4.2 Unteroberflächensynthese von Daten im Innern . . . . . . 40

3.5 Ermittlung unstetiger Leitfähigkeiten . . . . . . . . . . . . . . . . 45 3.5.1 Das Problem . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.5.2 Eine Approximation von u . . . . . . . . . . . . . . . . . . 45 3.5.3 Numerische Experimente . . . . . . . . . . . . . . . . . . . 47

4 Ein spezielles iteratives Verfahren 48

4.1 Das Iterationsverfahren . . . . . . . . . . . . . . . . . . . . . . . . 48 4.1.1 Der nichtlineare Fall . . . . . . . . . . . . . . . . . . . . . 48 4.1.2 Der lineare Fall . . . . . . . . . . . . . . . . . . . . . . . . 49 4.1.3 Die optimale Wahl von c_{j} . . . . . . . . . . . . . . . . . . 51

4.2 Anwendung auf die EIT . . . . . . . . . . . . . . . . . . . . . . . 51 4.2.1 Die Fréchet-Ableitung . . . . . . . . . . . . . . . . . . . . 51 4.2.2 Der adjungierte Operator . . . . . . . . . . . . . . . . . . 53

4.3 Zusammenfassung . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

5 Numerische Experimente 55

5.1 Approximation der Differentialgleichung . . . . . . . . . . . . . . 55 5.1.1 Der Fall Ω = (0, 1)^{2} . . . . . . . . . . . . . . . . . . . . . . 55 5.1.2 Der Fall Ω = {x ∈ IR^{2} : |x| < 1} . . . . . . . . . . . . . . . 57 5.1.3 Randdaten . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.1.4 Meßfehler . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.1.5 Numerische Durchführung einer Berechnung . . . . . . . . 58 5.2 Ergebnisse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

5.2.1 Erste Ergebnisse . . . . . . . . . . . . . . . . . . . . . . . 58 5.2.2 Rauschen . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.2.3 Trennung mehrerer Objekte . . . . . . . . . . . . . . . . . 62 5.2.4 Kleine Objekte . . . . . . . . . . . . . . . . . . . . . . . . 65 5.2.5 Objekte geringer Höhe . . . . . . . . . . . . . . . . . . . . 67 5.2.6 Geschwindigkeit . . . . . . . . . . . . . . . . . . . . . . . . 69

A Das C-Programm 71

B Mathematische Grundlagen 76

Kapitel 1

Das Problem

1.1 Inverse Probleme

Die Aufgabe in der Impedanz-Tomographie (EIT) ist es, die elektrische Leitfähigkeit im Innern eines Körpers nur durch Randmessungen, d.h. ohne selber in den Körper einzudringen, zu ermitteln.

Dieses inverse Problem präsentiert sich folgendermaßen: Sei σ eine passende Funktion, und davon abhängend sei L(σ) ein partieller Differentialoperator, der im zu behandelnden Fall ein streng elliptischer Operator zweiter Ordnung ist. Es bezeichnen N(σ) und D(σ) Lösungsoperatoren, die zu zwei unterschiedlichen Klassen von Randwertproblemen Lösungen u der Gleichung L(σ)u = 0 ergeben. Man kann annehmen, daß diese Lösungen jeweils eindeutig sind. Zu einem Paar (f, g) von Randwerten ist bei vorgegebenem σ das Vorwärtsproblem also überbestimmt.

Beim inversen Problem versucht man ganz allgemein aus dem Verhalten der Randwerte der Lösungen u von L(σ)u = 0 Rückschlüsse auf die Funktion σ zu ziehen. Man sucht ein σ, so daß für alle möglichen Paare (f, g) von Randwerten die entsprechenden Lösungen übereinstimmen, d.h.:

N(σ)f = D(σ)g für alle zulässigen Paare (f, g).

Die Randwertpaare ermittelt man aus Messungen am tatsächlichen Objekt.

Eine mögliche Anwendung dieses Verfahrens findet man in der Medizin. Barber und Brown haben in [B-B] eine Liste von Leitfähigkeiten unterschiedlicher Organe zusammengestellt. Mit Hilfe solcher Strukturinformationen kann man also Rückschlüsse auf den inneren Aufbau des Körpers ziehen.

1.2 Physikalische Gesetzmäßigkeiten

Sei Ω ⊂ IR^{n} ein beschränktes Gebiet. Für die Stromdichte J(x) und das elektrische Feld E(x) gilt das Kirchhoffsche Gesetz:

∇ · J = 0,
∇ × E = 0,

d.h. die Divergenz von J und die Rotation von E sind Null. Außerdem sind J und E durch das Ohmsche Gesetz miteinander verknüpft:

J(x) = σ(x)E(x),

dabei ist σ(x) die Leitfähigkeit im Punkt x. In diesem Modell und nimmt man an, daß σ eine reellwertige Funktion mit positiver unterer Schranke ist.

Aus ∇ × E = 0 folgt, daß es eine Funktion u gibt mit

E = −∇u.

Hier ist u die Spannung. So erhält man eine einzige Differentialgleichung zweiter Ordnung für u

∇ · σ∇u = 0, in Ω.  (1.1)

Wegen σ(x) > 0 ∀x ∈ Ω ist diese Differentialgleichung streng elliptisch.

In diesem Modell berechnet man die Spannung u in Ω, indem man Randwerte auf zwei unterschiedliche Arten vorgibt:

u = g auf ∂Ω  (1.2)
oder σ \frac{∂u}{∂ν} = f auf ∂Ω.  (1.3)

Dabei bezeichne ν die äußere Normale an Ω. Es sind f die Spannung und g die Stromdichte auf dem Rand des Körpers. Als Nebenbedingung hat man

∫_{∂Ω} u dS = 0

und ∫_{∂Ω} σ \frac{∂u}{∂ν} dS = 0.

1.3 Eindeutige Lösbarkeit der Differentialgleichung

In diesem Abschnitt möchte ich die eindeutige Lösbarkeit der Randwertprobleme (1.1) mit Randwerten (1.2) oder (1.1) mit Randwerten (1.3) zeigen. Dazu müssen vorher entsprechende Funktionenräume und schwache Formulierungen der Probleme angegeben werden.

Vereinbarungen

Sei Ω ⊂ IR^{n} ein Gebiet mit hinreichend glattem Rand. Als Funktionenraum für die Lösungen kann man folgenden linearen Teilraum von H^{1}(Ω) nehmen:

H := {u ∈ H^{1}(Ω), ∫_{∂Ω} u dS = 0}.  (1.4)

Da die Abbildung u ↦ u|_{∂Ω} von H^{1}(Ω) nach H^{1/2}(∂Ω) stetig ist, ist H abgeschlossen in H^{1}(Ω). Damit ist H ein Hilbertraum. Sei σ ∈ L^{∞}(Ω) mit positiver unterer Schranke σ_{o}. Dann kann man auf H die Bilinearform

a(u, v) := ∫_{Ω} σ∇u · ∇v dx  (1.5)

definieren. Da Null die einzige konstante Funktion in H ist, folgt aus

a(u, u) = ∫_{Ω} σ|∇u|^{2} dx
≥ σ_{0} ∫_{Ω} |∇u|^{2} dx, ∀ u ∈ H

und

|a(u, v)| ≤ ∥σ∥_{L^{∞}(Ω)}∥u∥_{H_{0}^{1}(Ω)}∥v∥_{H_{0}^{1}(Ω)}
≤ ∥σ∥_{L^{∞}(Ω)}∥u∥_{H^{1}(Ω)}∥v∥_{H^{1}(Ω)}, ∀ u, v ∈ H,

daß a eine koerzive beschränkte Bilinearform auf H ist.

Das Dirichlet-Problem

Um schwache Lösungen von (1.1) mit Randwerten (1.2) zu erklären, definiere ich das lineare Funktional D_{g} auf H durch

D_{g}(v) := ∫_{∂Ω} σ \frac{∂v}{∂ν} g dS.

Hier ist ν wieder die äußere Normale an Ω und dS das übliche Oberflächenmaß

auf ∂Ω. Da \frac{∂v}{∂ν} |_{∂Ω} ∈ H^{−1/2}(∂Ω) für v ∈ H^{1}(Ω), ist D_{g} für alle g ∈ H^{1/2}(∂Ω)

beschränkt, also stetig (H^{−1/2}(∂Ω) ist der Dualraum von H^{1/2}(∂Ω) bzgl. des Skalarprodukts auf L^{2}(∂Ω)). Man kann nun u ∈ H als schwache Lösung von (1.1) mit Randwerten (1.2) erklären, wenn u die Gleichung

a(u, v) = D_{g}(v) ∀ v ∈ H  (1.6)

erfüllt. Aus dem Gauß’schen Integralsatz folgt, daß jede klassische Lösung auch schwache Lösung ist.

Das von-Neumann-Problem

Analog zum Dirichlet-Problem erklärt man schwache Lösungen u ∈ H von (1.1) mit Randwerten (1.3) als Lösungen von

a(u, v) = N_{f}(v) ∀ v ∈ H  (1.7)

mit

N_{f}(v) := ∫_{∂Ω} fv dS

für f ∈ H^{−1/2}(∂Ω).

Eindeutige Lösbarkeit

Nun folgt aus dem Satz von Lax Milgram, daß die Probleme (1.6) und (1.7) in H eindeutig lösbar sind.

1.4 Elektrische-Impedanz-Tomographie

In der Mathematik hat das Problem zunächst folgende Gestalt: Es sei die Abbildung, die einer vorgegebenen Randableitung f die zugehörigen Randwerte g zuordnet, bekannt. Diese Abbildung will ich an dieser Stelle mit Λ_{σ} bezeichnen. Aufgrund der eindeutigen Lösbarkeit des Vorwärtsproblems kann man natürlich auch die Umkehrabbildung betrachten.

Die gesuchte Leitfähigkeit σ kann man nun aus der vollständigen Kenntnis von Λ_{σ} auf gewissen Funktionenräumen (z.B. H^{s±1}(∂Ω) für ein Gebiet Ω ⊂ IR^{n}) eindeutig bestimmen.

In der Praxis kennt man diese Abbildung aber nur auf endlich-dimensionalen linearen Teilräumen, weil man natürlich nur endlich viele Messungen durchführen kann.

Das inverse Problem stellt sich also folgendermaßen:

Sei Ω ⊂ IR^{n} ein Gebiet mit hinreichend glattem Rand, und seien {f_{i}}_{i=1}^{N} ⊂ H^{−1/2}(∂Ω) und {g_{i}}_{i=1}^{N} ⊂ H^{1/2}(∂Ω) Folgen von Randdaten. Sei {u_{i}}_{i=1}^{N} eine Folge von Lösungen des von-Neumann-Problems

\begin{aligned}∇ · σ∇u_{i} & = & 0 & in & Ω, \\ σ \frac{∂u_{i}}{∂ν} & = & f_{i} & auf & ∂Ω \end{aligned}  (1.8)

und {v_{i}}_{i=1}^{N} eine Folge von Lösungen des Dirichlet-Problems

∇ · σ∇v_{i} = 0 in Ω,
v_{i}|_{∂Ω} = g_{i} auf ∂Ω. (1.9)

Dabei bezeichnet ν die äußere Normale an Ω. Man sucht nun ein σ, so daß diese Lösungen u_{i} und v_{i} übereinstimmen für i = 1, …, N. Aus der eindeutigen Lösbarkeit der Randwertaufgaben folgt, daß Übereinstimmung auf dem Rand von Ω ausreicht.

Beim Vorwärtsproblem bewirken große Änderungen bei σ nur sehr geringe Änderungen in den Randdaten. Deshalb ist das inverse Problem schlecht gestellt im Sinne von Hadamard, da das Ergebnis nicht stetig von den Anfangsdaten abhängt.

In der Regel sind die {f_{i}}_{i=1}^{N} vorgegeben und die {g_{i}}_{i=1}^{N} Meßwerte. Im idealen Modell, d.h., falls die g_{i} nicht mit Meßfehlern behaftet sind, ist die Existenz einer Lösung des inversen Problems klar. Hier bleibt nur zu prüfen, ob diese Lösung auch eindeutig ist, was sich als ein Ergebnis der Nachman-Sylvester-Uhlman-Methode ergeben wird.

Ein analytisches Verfahren

In diesem Abschnitt wird ein analytisches Verfahren zur Bestimmung von σ vorgestellt. Dieses Verfahren kann man Inverse Scattering oder Nachman-Sylvester-Uhlman-Methode nennen. Leider ist eine Voraussetzung, daß die Dimension des Raumes mindestens drei ist. Aber diese konstruktive Methode beinhaltet gleichzeitig einen Eindeutigkeitsbeweis in gewissen Funktionenräumen.

2.1 Inverse Scattering

2.1.1 Vereinbarungen

Definition 2.1.1 Sei Ω ⊂ IR^{n} ein beschränktes Gebiet. C^{1,1}(Ω) ist der Raum aller stetig differenzierbaren Funktionen f, deren erste Ableitungen gleichmäßig Lipschitz-stetig sind, d.h.:

∃K > 0 derart, daß |\frac{∂f}{∂v} (x) − \frac{∂f}{∂v} (y)| ≤ K|x − y|, ∀x, y ∈ Ω

für eine beliebige Richtung v ∈ IR^{n}.

Es gilt: f ∈ C^{1,1}(\overline{Ω}) :⇔ es existiert ein beschränktes Gebiet U mit \overline{Ω} ⊂ U ⊂ IR^{n} und f ∈ C^{1,1}(U).

In diesem Abschnitt gilt:

• n ≥ 3, und Ω ⊂ IR^{n} ist ein beschränktes Gebiet.

• Der Rand von Ω ist C^{1,1}.

• Die Leitfähigkeit σ ist eine reellwertige Funktion aus C^{1,1}(\overline{Ω}) mit positiver unterer Schranke.

• Die Funktion q bezeichnet eine Funktion in L^{∞}(IR^{n}) mit kompaktem Träger (in der Regel im Abschluß von Ω).

• ν bezeichnet die äußere Normale an Ω (ebenso ν_{+} und ν_{−})

Definition 2.1.2 L_{δ}^{2}(IR^{n}) ist der Raum aller Funktionen, für die (1 + |x|^{2})^{δ/2}f in L^{2}(IR^{n}) liegt:

L_{δ}^{2}(IR^{n}) = \{f : ∥f∥_{δ} =: (\begin{matrix}∫ \\ IR & ^{n} \end{matrix} (1 + |x|^{2})^{δ}|f(x)|^{2} dx)^{\frac{1}{2}} < ∞\} .

Ähnlich wie in Abschnitt 1.4 definiert man für das Vorwärtsproblem und für die Helmholtz-Gleichung die Dirichlet-zu-Neumann-Abbildungen.

Definition 2.1.3 Seien f und g in H^{1/2}(∂Ω), σ und q wie oben und u und w Lösungen folgender Gleichungen:

−Δw + qw = 0, in Ω,
w|_{∂Ω} = f,  (2.1)
∇ · σ∇u = 0, in Ω,
u|_{∂Ω} = g.  (2.2)

Definiere nun dazu die linearen Abbildungen Λ_{q} und Λ_{σ} von H^{1/2}(∂Ω) nach H^{−1/2}(∂Ω) durch

Λ_{q}f = \frac{∂w}{∂ν}|_{∂Ω} ,  (2.3)
Λ_{σ}g = σ \frac{∂u}{∂ν}|_{∂Ω} .  (2.4)

2.1.2 Umformulierung des Problems

Zunächst transformiert man das Problem in ein inverses Problem für die Helmholtz-Gleichung. Auf dieses wendet man dann ein Verfahren aus der Akustik an: die inverse Scattering-Methode. Bei diesem Verfahren ersetzt man das Problem auf dem Rand von Ω durch ein Problem im Unendlichen.

Satz 2.1.4 Es gelte u = σ^{−1/2}w und q = σ^{−1/2}Δσ^{1/2}. Dann sind die partiellen Differentialgleichungen

∇ · σ∇u = 0
und  − Δw + qw = 0

äquivalent.

Bemerkung 2.1.5 Da σ ∈ C^{1,1}(\overline{Ω}) ist, folgt, daß q ∈ L^{∞}(IR^{n}) mit Träger im Abschluß von Ω ist.

Beweis:

Δw = Δ (σ^{1/2}u)
= ∑ \frac{∂^{2}}{∂x_{i}^{2}} (σ^{1/2}u)
= ∑ [\frac{∂^{2}(σ^{1/2})}{∂x_{i}^{2}}u + 2\frac{∂(σ^{1/2})}{∂x_{i}} \frac{∂u}{∂x_{i}} + σ^{1/2} \frac{∂^{2}u}{∂x_{i}^{2}}]
= ∑ [\frac{∂^{2}(σ^{1/2})}{∂x_{i}^{2}}u + 2\frac{1}{2}σ^{−1/2} \frac{∂σ}{∂x_{i}} \frac{∂u}{∂x_{i}} + σ^{1/2} \frac{∂^{2}u}{∂x_{i}^{2}}]
= (σ^{−1/2}Δ(σ^{1/2})) σ^{1/2}u + σ^{−1/2} (∇ · σ∇u)
= qw + σ^{−1/2}(∇ · σ∇u).
□

Um auch die Dirichlet-zu-Neumann-Abbildungen zu transformieren, benötigt man die Randwerte von σ. Deshalb zunächst ein Lemma. Der etwas umfangreichere Beweis wird in Abschnitt 2.5 nachgeliefert.

Lemma 2.1.6 Die Randwerte σ|_{∂Ω} und \frac{∂σ}{∂ν} |_{∂Ω} lassen sich aus Λ_{σ} berechnen.

Satz 2.1.7 Mit den Bezeichnungen aus Definition 2.1.3 und den Bedingungen aus Satz 2.1.4 gilt:

Λ_{q}f = σ^{−1/2}|_{∂Ω}Λ_{σ}σ^{−1/2}|_{∂Ω}f + \frac{1}{2}σ^{−1}|_{∂Ω} \frac{∂σ}{∂ν} |_{∂Ω} f.  (2.5)

Beweis: Es gilt w = σ^{1/2}u und deshalb g = σ^{−1/2}|_{∂Ω}f.

\frac{∂w}{∂v} = \frac{∂(σ^{1/2}u)}{∂v}
= σ^{1/2} \frac{∂u}{∂v} + \frac{1}{2}σ^{−1/2} \frac{∂σ}{∂v} u
= σ^{−1/2} (σ \frac{∂u}{∂v}) + \frac{1}{2}σ^{−1} \frac{∂σ}{∂v} w

für eine beliebige Richtung v ∈ IR^{n} \ {0}. Es folgt mit Lemma 2.1.6:

Λ_{q}f = \frac{∂w}{∂ν}|_{∂Ω}
= σ^{−1/2}|_{∂Ω}Λ_{σ}g + \frac{1}{2}σ^{−1}|_{∂Ω} \frac{∂σ}{∂ν} |_{∂Ω} f
= σ^{−1/2}|_{∂Ω}Λ_{σ}(σ^{−1/2}|_{∂Ω}f) + \frac{1}{2}σ^{−1}|_{∂Ω} \frac{∂σ}{∂ν} |_{∂Ω} f.
□

Eine wichtige Rolle in dieser Methode spielt die Scattering-Amplitude:

Definition 2.1.8 Sei q ∈ L^{∞}(IR^{n}) mit kompaktem Träger und ψ(x, ζ) eine Lösung von

−Δψ + qψ = 0 in IR^{n},  (2.6)

die sich asymptotisch wie e^{ix·ζ}, ζ · ζ = 0, verhält. Definiere nun die (modifizierte) Scattering-Amplitude t(ξ, ζ) durch

t(ξ, ζ) :=\begin{matrix}∫ \\ IR & ^{n} \end{matrix} e^{−ix·(ξ+ζ)}q(x)ψ(x, ζ) dx.  (2.7)

In Abschnitt 2.4 wird gezeigt, wie man t(ξ, ζ) aus Λ_{q} berechnen kann. Schließlich kann man aus t(ξ, ζ) nun q berechnen.

Satz 2.1.9 Seien q ∈ L^{∞}(IR^{n}) mit kompaktem Träger und δ = \frac{1}{2}. Dann gilt für die Fouriertransformation q̂ von q für |ζ| ≥ a

|t(ξ, ζ) − q̂(ξ)| ≤ \frac{c̃(δ, a)}{|ζ|}∥q∥_{δ}^{2}.  (2.8)

Der Beweis wird in Abschnitt 2.3 nachgereicht.

Bemerkung 2.1.10 Eine alternative Methode q̂ aus t(ξ, ζ) zu berechnen nennt Nachman in [N]: q̂ und damit auch q sind gegeben durch

q̂(ξ) = t(ξ, ξ_{0}) + ∫_{{|ζ|=M \atop ζ∈V_{ξ}}} t(ξ, ζ)K(ζ, ζ_{0})  (2.9)

−\frac{1}{(2π)^{n−1}} ∫ ∫_{{|ζ|>M \atop ζ∈V_{ξ}}} δ(η^{2} + 2ζ · η)t(ξ − η, ζ + η)t(η, ζ)η d\overline{ζ} dη ∧ K(ζ, ζ_{0})

mit

K(ζ, z) = (−1)^{((n−2)(n−3)/2)+1} \frac{(n − 3)!}{(2πi)^{n−2}} \frac{α(ζ, z) ∧ β(z)}{|ζ − z|^{2(n−2)}},

α(ζ, z) = 2 ∑_{1≤j<k<l≤n}(−1)^{j+k+l} | \begin{matrix}\overline{ζ}_{j} − \overline{z}_{j} & \overline{ζ}_{k} − \overline{z}_{k} & \overline{ζ}_{l} − \overline{z}_{l} \\ ζ_{j} + z_{j} & ζ_{k} + z_{k} & ζ_{l} + z_{l} \\ ξ_{j} & ξ_{k} & ξ_{l} \end{matrix} | ⋀_{m≠j,k,l} d\overline{ζ}_{m},

β(ζ) = \frac{1}{4 ∑_{1≤j<k≤n} |ζ_{j}ξ_{k} − ζ_{k}ξ_{j}|^{2}} ∑_{1≤j<k≤n}(−1)^{j+k}(\overline{ζ}_{j}ξ_{k} − \overline{ζ}_{k}ξ_{j}) ⋀_{l≠j,k} dζ_{l},
V_{ξ} = {ζ ∈ l C^{n}, ζ^{2} = 0, ξ^{2} + 2ζ · ξ = 0} .

2.2 Hilfsprobleme und Fundamentallösungen

Um Lösungen von −Δψ + qψ = 0 zu konstruieren, die sich für ζ ∈ l C^{n}, n ≥ 3, |ζ| groß, ζ · ζ = 0, wie e^{ix·ζ} verhalten, betrachtet man Lösungen der Form

ψ = e^{ix·ζ}(1 + v_{ζ}).

Dann gilt

e^{−ix·ζ}Δψ = ζ^{2}(1 + v_{ζ}) + 2iζ · ∇v_{ζ} + Δv_{ζ}
= (2iζ · ∇ + Δ)v_{ζ}.

Mittels Faltung und Fouriertransformation der Dirac’schen-Delta-Distribution δ erhält man eine Fundamentallösung zum Differentialoperator −2iζ − Δ:

g(x, ζ) = g_{ζ}(x) = \frac{1}{(2π)^{n}} \begin{matrix}∫ \\ IR & ^{n} \end{matrix} \frac{e^{ix·ξ}}{ξ^{2} + 2ζ · ξ} dξ
⇒ −Δg_{ζ}(x) − 2iζ · ∇g_{ζ}(x) = δ.

Dieses ist möglich, da (ξ^{2} + 2ζ · ξ)^{−1} als Funktion von ξ ∈ IR^{n} wegen n ≥ 3 lokal integrierbar ist. So bekommt man für v_{ζ} die Integralgleichung

v_{ζ} = g_{ζ} ∗ (e^{i⟨·,ζ⟩}q(1 + v_{ζ}))

und damit für ψ

ψ = e^{ix·ζ} −\begin{matrix}∫ \\ IR & ^{n} \end{matrix} G_{ζ}(x − y)q(y)ψ(y) dy  (2.10)

mit G_{ζ}(x) = e^{ix·ζ}g_{ζ}(x).

G_{0}(x) = \frac{|x|^{2−n}}{(n − 2)ω_{n}} , ω_{n} = \frac{2π^{n/2}}{Γ(\frac{n}{2})} ,

ist die klassische Fundamentallösung der Laplace-Gleichung, und daher gibt es eine harmonische Funktion H_{ζ} mit

G_{ζ} = G_{0} + H_{ζ}, ΔH_{ζ} = 0 in IR^{n}.  (2.11)

Deshalb lassen sich viele Ergebnisse, die für G_{0} gelten, auch auf G_{ζ} übertragen.

Für die Familie G_{ζ} Green’scher Funktionen für den Laplaceoperator definiert man das single-layer-Potential S_{ζ}, das double-layer-Potential D_{ζ} und das tracelayer-Potential B_{ζ} wie folgt:

G_{ζ}(x, y) := G_{ζ}(x − y),
S_{ζ}f(x) = ∫_{∂Ω} G_{ζ}(x, y)f(y) dS_{y}, x ∈ IR^{n} \ ∂Ω,

D_{ζ}f(x) = ∫_{∂Ω} \frac{∂G_{ζ}}{∂ν_{y}} (x, y)f(y) dS_{y}, x ∈ IR^{n} \ ∂Ω,

B_{ζ}f(x) = p.v. ∫_{∂Ω} \frac{∂G_{ζ}}{∂ν_{y}} (x, y)f(y) dS_{y}, x ∈ ∂Ω.

Für ζ = 0 haben Fabes, Jodeit und Rivère in [F-J-R] eine Reihe von Eigenschaften dieser Potentiale bewiesen, die im weiteren Verlauf der Arbeit noch benötigt werden. Hier fließt auch die Einschränkung n ≥ 3 ein, da G_{0} für n = 2 eine vollkommen andere Gestalt annimmt. Da sich G_{ζ} von G_{0} nur durch Addition einer harmonischen Funktion unterscheidet, behalten diese Aussagen auch noch für ζ ≠ 0 ihre Gültigkeit (siehe auch Nachman [N]).

Sei im folgenden Ω ⊂ IR^{n}, n ≥ 3, ein beschränktes Gebiet mit C^{1,1} Rand.

a) Für 0 ≤ s ≤ 1 gilt

∥S_{ζ}f∥_{H^{s+1}(∂Ω)} ≤ c(ζ, s)∥f∥_{H^{s}(∂Ω)}.  (2.12)

b) Für 0 ≤ s ≤ \frac{3}{2} ist B_{ζ} beschränkt auf H^{s}(Ω).

Sei nun ρ_{0} so groß, daß \overline{Ω} ⊂ \{x, |x| < ρ_{0}\}; für beliebige ρ > ρ_{0} sei Ω′_{ρ} := \{x : x ∉ \overline{Ω}, |x| < ρ\}. Für f ∈ H^{1/2}(∂Ω) besitzen u = S_{ζ}f und v = D_{ζ}f die Eigenschaften

c) u erfüllt

Δu = 0 in IR^{n} \ Ω.  (2.13)

d) u ist in H^{2}(Ω) und in H^{2}(Ω′_{ρ}).

e) Die Grenzwerte u_{+} und u_{−} von u gegen ∂Ω von außerhalb bzw. innerhalb von Ω stimmen überein, und für die entsprechenden Grenzwerte der Nor-

malenableitung \frac{∂u}{∂ν} gilt

\frac{∂u}{∂ν_{−}} − \frac{∂u}{∂ν_{+}} = f.  (2.14)

f) Der nichttangtiale Grenzwert v_{+} (v_{−}) von v gegen ∂Ω von außerhalb (innerhalb) von Ω erfüllt

v_{±}(x) = ±\frac{1}{2}f(x) + B_{ζ}f(x) für fast alle x ∈ ∂Ω.  (2.15)

Für die Beweise dieser Behauptungen sei auf Fabes, Jodeit und Rivère [F-J-R] zusammen mit Anhang A zu Nachman [N] verwiesen.

Sylvester und Uhlmann [S-U] waren die ersten, die die gewichteten Räume L_{δ}^{2} einführten (siehe Definition 2.1.2). In diesen Räumen kann man die Lösbarkeit von Gleichung (2.10) zeigen. Dazu benötigt man u.a. folgendes Lemma, welches sogar für n ≥ 2 gilt:

Lemma 2.2.1 Sei n ≥ 2, Ω ⊂ IR^{n} ein beschränktes Gebiet mit C^{1,1} Rand. Dann erfüllt ∀ζ ∈ l C^{n} mit ζ · ζ = 0 und |ζ| > a die Faltung mit g_{ζ}

∥g_{ζ} ∗ f∥_{δ−1} ≤ \frac{c(δ, a)}{|ζ|}∥f∥_{δ}, 0 < δ < 1.  (2.16)

Beweis: Für f ∈ L_{δ}^{2}(IR^{2}) definiert man den Operator Z durch

Zf(u) = f ∗ \frac{1}{u_{1} + iu_{2}} .  (2.17)

Für g ∈ L_{1−δ}^{2}(IR^{2}), dem Dualraum von L_{δ−1}^{2}(IR^{2}), gilt

⟨g, Zf⟩ =\begin{matrix}∫ \\ IR & ^{2} \end{matrix} \begin{matrix}∫ \\ IR & ^{2} \end{matrix} \frac{g(u)f(v)}{(u_{1} − v_{1}) + i(u_{2} − v_{2})} dv du.

Anwenden der Cauchy-Schwarz’schen-Ungleichung ergibt

|⟨g, Zf⟩|^{2} ≤ \begin{matrix}∫ \\ IR & ^{2} \end{matrix} \begin{matrix}∫ \\ IR & ^{2} \end{matrix} \frac{(1 + |u|^{2})^{δ}|f|^{2}}{(1 + |u|^{2})^{(δ−1)/2}|u − v|(1 + |v|^{2})^{(2−δ)/2}} dv du

·\begin{matrix}∫ \\ IR & ^{2} \end{matrix} \begin{matrix}∫ \\ IR & ^{2} \end{matrix} \frac{(1 + |v|^{2})^{1−δ}|g|^{2}}{(1 + |u|^{2})^{(δ+1)/2}|u − v|(1 + |v|^{2})^{−δ/2}} du dv

≤ M∥g∥_{L_{1−δ}^{2}}^{2}∥f∥_{L_{δ}^{2}}^{2} ∀ g ∈ L_{1−δ}^{2}(IR^{2})

mit M^{2} = ( sup_{u∈IR^{2}}\begin{matrix}∫ \\ IR & ^{2} \end{matrix} \frac{(1 + |u|^{2})^{(1−δ)/2}}{|u − v|(1 + |v|^{2})^{(2−δ)/2}} dv)

· ( sup_{v∈IR^{2}}\begin{matrix}∫ \\ IR & ^{2} \end{matrix} \frac{(1 + |v|^{2})^{δ/2}}{(1 + |u|^{2})^{(δ+1)/2}|u − v|} du) .

Wegen 0 < δ < 1 gilt

0 < \frac{1 − δ}{2} < \frac{1}{2} < \frac{2 − δ}{2} < 1
und 0 < \frac{δ}{2} < \frac{1}{2} < \frac{δ + 1}{2} < 1.

Deshalb ist M endlich.

Für x ∈ IR^{n}, n > 2, x_{1} = u_{1}, x_{2} = u_{2} gilt nun wegen δ > 0 und δ − 1 < 0

(1 + |u|^{2})^{δ−1} ≥ (1 + |x|^{2})^{δ−1}
und (1 + |u|^{2})^{δ} ≥ (1 + |x|^{2})^{δ}.

Also ist Z ein beschränkter Operator von L_{δ}^{2}(IR^{n}) nach L_{δ−1}^{2}(IR^{n}), n ≥ 2:

∥Zf∥_{L_{δ−1}^{2}} ≤ M∥f∥_{L_{δ}^{2}}.  (2.18)

Man betrachte nun das Problem L_{ζ}w = Δw + 2iζ · ∇w = f mit ζ · ζ = 0 und |ζ| > a > 0. Es sei ohne Einschränkung 2iζ = s(e_{1} + ie_{2}). Dann führt Fouriertransformation zu

l(η) = −|η|^{2} − sη_{2} + isη_{1},
l(η)ŵ = f̂,
M := {η | l(η) = 0}.

Definiere nun Überdeckungen von M und IR^{n} durch

U_{2} = M ∩ \{η | (η_{2} − \frac{s}{2})^{2} > \frac{s^{2}}{8n}\},
U_{j} = M ∩ \{η | η_{j}^{2} > \frac{s^{2}}{8n}\} j = 3, …, n,
V_{1} = IR^{n} \ N_{ε/2}(M),
V_{j} = N_{ε}(U_{j}), j = 2, …, n, \frac{a}{4} < ε < \frac{s}{4},

wobei N_{ε}(A) die ε-Umgebung einer Menge A bezeichne.

Sei ρ_{j} eine Zerlegung der Eins mit Träger in V_{j} und

f̂_{j} = ρ_{j}f̂.

In V_{1} gilt

|l(η)| > \frac{ηs}{2} > \frac{as}{8}.

Hier ist also nichts zu zeigen.

In V_{j}, j = 2, …, n führt man neue Koordinaten β^{1}, …, β^{n} ein, so daß

l(β) = s(β^{j} + iβ^{1}),

und man definiert

û_{j} = \frac{f̂_{j}}{l} .

Sylvester und Uhlmann haben in [S-U] gezeigt, daß die Behauptung auch unter den notwendigen Koordinatentransformationen ϕ_{j} erhalten bleibt. Also folgt mit Plancherel

∥u∥_{δ−1} ≤ ∥ϕ_{j}^{−1} (Z (ϕ_{j}(f_{j}))) ∥_{δ−1}
≤ \frac{C}{s} ∥f∥_{δ}.  (2.19)

Die Konstante C hängt nur von δ und a ab. Daraus ergibt sich die Behauptung für g_{ζ} ∗ f, da −g_{ζ} Fundamentallösung zu L_{ζ} ist. □

Lemma 2.2.2 Seien Ω und Ω′_{ρ} wie oben und u = S_{ζ}f (u = D_{ζ}f) für f ∈ H^{1/2}(∂Ω) (f ∈ H^{3/2}(∂Ω)). Dann gilt:

lim_{ρ→0} ∫_{|y|=ρ} [G_{ζ}(x, y)\frac{∂u}{∂ν_{y}} − u(y)\frac{∂G_{ζ}(x, y)}{∂ν_{y}}] dS_{y} = 0  (2.20)

und sogar

∫_{|y|=ρ} [G_{ζ}(x, y)\frac{∂u}{∂ν_{y}} − u(y)\frac{∂G_{ζ}(x, y)}{∂ν_{y}}] dS_{y} = 0.  (2.21)

Beweis: Man regularisiert G_{ζ} durch

G_{ζ}^{ε}(x, y) = \frac{(|x − y|^{2} + ε^{2})^{(2−n)/2}}{(n − 2)ω_{n}} + H_{ζ}. (vgl. (2.11))

Dann wendet man den Gauß’schen Integralsatz auf G_{ζ}^{ε}(x, y) und G_{ζ}^{ε′}(y, z) an. Der Grenzübergang ε → 0, die Vertauschung der Integrationsreihenfolge und anschließend der Grenzübergang ε′ → 0 liefern dann das gewünschte Ergebnis. □

2.3 Beweis von Satz 2.1.9

Korollar 2.3.1 Sei δ = \frac{1}{2} und |ζ| > c(δ, a)∥q(x)(1 + x^{2})^{δ}∥_{∞}. Dann existiert eine eindeutige Lösung ψ(x, ζ) von (2.10), und diese erfüllt

∥e^{−ix·ζ}ψ(·, ζ) − 1∥_{−δ} ≤ \frac{c̃(δ, a)}{|ζ|}∥q∥_{δ}.  (2.22)

Beweis: In diesem Beweis verwende ich zur besseren Verdeutlichung e^{i⟨x,ζ⟩} anstelle von e^{ix·ζ}.

Der Operator T : L_{−δ}^{2} → L_{δ}^{2}, ϕ ↦ q · ϕ ist durch ∥q(x)(1 + x^{2})^{δ}∥_{∞} beschränkt. Also ist der Operator K_{ζ} : L_{−δ}^{2} → L_{−δ}^{2}, definiert durch

K_{ζ}ϕ = g_{ζ} ∗ (T ϕ),

beschränkt durch c(δ, a)|ζ|^{−1}∥q(x)(1 + x^{2})^{δ}∥_{∞} < 1. Deshalb ist (I + K_{ζ}) invertierbar. Definiere ψ(x, ζ) durch

ψ(x, ζ) = e^{i⟨x,ζ⟩} (1 − (I + K_{ζ})^{−1}(g_{ζ} ∗ q)) .

Für eine beliebige Lösung ϕ von (2.10) gilt:

ϕ(x) = e^{i⟨x,ζ⟩} −\begin{matrix}∫ \\ IR & ^{n} \end{matrix} e^{i⟨x−y,ζ⟩}g_{ζ}(x − y)q(y)ϕ(y) dy
⇔  e^{−i⟨x,ζ⟩}ϕ(x) − 1 = −\begin{matrix}∫ \\ IR & ^{n} \end{matrix} g_{ζ}(x − y)q(y)e^{−i⟨y,ζ⟩}ϕ(y) dy
= −K_{ζ}(e^{−i⟨·,ζ⟩}ϕ).
⇒ (I + K_{ζ})(e^{−i⟨·,ζ⟩}ϕ(·) − 1)(x) = e^{−i⟨x,ζ⟩}ϕ(x) − 1 − (e^{−i⟨x,ζ⟩}ϕ(x) − 1) − K_{ζ}(1)(x)
= −(g_{ζ} ∗ q)(x).

Für zwei Lösungen ψ_{1} und ψ_{2} gilt

(I + K_{ζ}) (e^{−i⟨x,ζ⟩}(ψ_{1}(x, ζ) − ψ_{2}(x, ζ))) = 0.

Da (I + K_{ζ}) invertierbar ist, folgt ψ_{1} = ψ_{2}. Die Abschätzung (2.22) ergibt sich aus

e^{−i⟨x,ζ⟩}ψ(x, ζ) − 1 = −K_{ζ}(e^{−i⟨·,ζ⟩}ψ(·, ζ))

mit dem vorher Erwähnten und Lemma 2.2.1.

□

Damit kann man nun Satz 2.1.9 beweisen: Beweis:

t(ξ, ζ) = \begin{matrix}∫ \\ IR & ^{n} \end{matrix} e^{−ix·(ξ+ζ)}q(x)ψ(x, ζ) dx

= q̂(ξ) +\begin{matrix}∫ \\ IR & ^{n} \end{matrix} e^{−ix·ξ}q(x) (ψ(x, ζ)e^{−ix·ζ} − 1) dx

und

∥q(x) (ψ(x, ζ)e^{−ix·ζ} − 1)∥_{L^{1}(IR^{n})} ≤ ∥q∥_{δ}∥ψ(x, ζ)e^{−ix·ζ} − 1∥_{−δ}
≤ \frac{c̃(δ, a)}{|ζ|}∥q∥_{δ}^{2}.
□

2.4 Berechnung der Scattering-Amplitude

Im Verlauf dieses Abschnitts möchte man den Gauß’schen Integralsatz auf G_{ζ} anwenden. Da dieses nicht geht, verwendet man zunächst wie im Beweis von Lemma 2.2.2

G_{ζ}^{ε}(x, y) = \frac{(|x − y|^{2} + ε^{2})^{(2−n)/2}}{(n − 2)ω_{n}} + H_{ζ}.  (2.23)

In Definition 2.1.8 betrachtet man das Problem (2.6). Es sei Ω′ = IR \ \overline{Ω}. Man teilt dieses Problem in ein inneres und ein äußeres auf:

Inneres Problem: Für Lösungen u von

−Δu + qu = 0 in Ω, u|_{∂Ω} = f,

definiert man den Lösungsoperator

P_{q} : H^{3/2}(∂Ω) → H^{2}(Ω), P_{q}f := u.

Äußeres Problem: Es sei h ∈ C^{2}(IR^{n}) harmonisch.

Δψ = 0 in Ω′_{ρ},  (2.24)
ψ ∈ H^{2}(Ω′_{ρ})  für alle ρ > ρ_{0},  (2.25)
ψ(x, ζ) − h(x)  erfüllt (2.20),  (2.26)
\frac{∂ψ}{∂ν_{+}} = Λ_{q}ψ auf ∂Ω.  (2.27)

Es sei ψ eine Lösung von (2.6), die sich wie, in Definition 2.1.8 beschrieben, wie e^{ix·ζ} verhält. Mit dem Gauß’schen Integralsatz in Ω gilt für (ξ + ζ)^{2} = 0 und f_{ζ} := ψ|_{∂Ω}

t(ξ, ζ) = ∫_{Ω} e^{−ix·(ξ+ζ)}q(x)ψ(x) dx
= ∫_{Ω} e^{−ix·(ξ+ζ)}Δψ(x) dx
= ∫_{∂Ω} e^{−ix·(ξ+ζ)} (Λ_{q}f_{ζ} + i(ξ + ζ) · νf_{ζ}) dS.  (2.28)

Da Λ_{q} bekannt ist, benötigt man nur die Randwerte f_{ζ}. Das Problem ist aber, wie sieht man den Randwerten an, daß die korrespondierende Lösung sich wie e^{ix·ζ}, ζ · ζ = 0, verhält.

Es sei

T_{ζ} := −S_{ζ}Λ_{q} + B_{ζ} + \frac{1}{2}I,
R_{ζ} := −S_{ζ}Λ_{q} + D_{ζ},

und ψ löse (2.24)-(2.27). Da h harmonisch ist, gilt

h(x) = ∫_{|y|=ρ} G_{ζ}(x, y)\frac{∂h(y)}{∂ν_{y}} − \frac{∂G_{ζ}(x, y)}{∂ν_{y}}h(y) dS_{y}.

Der Gauß’sche Integralsatz angewandt auf ψ und G_{ζ}^{ε} in Ω′_{ρ} mit anschließendem Grenzübergang ε → 0 liefert mit f = ψ|_{∂Ω}

ψ(x) = ∫_{|y|=ρ} G_{ζ}(x, y)\frac{∂(ψ(y) − h(y))}{∂ν_{y}} − \frac{∂G_{ζ}(x, y)}{∂ν_{y}}(ψ(y) − h(y)) dS_{y}

+ ∫_{|y|=ρ} G_{ζ}(x, y)\frac{∂h(y)}{∂ν_{y}} − \frac{∂G_{ζ}(x, y)}{∂ν_{y}}h(y) dS_{y} + R_{ζ}f

→_{ρ→0} h(x) + R_{ζ}f

wegen (2.26). Also erfüllen die Randwerte f von ψ aufgrund von (2.15) für den nichttangentialen Grenzwert von x ∈ Ω′_{ρ} gegen ∂Ω die Integralgleichung

f = h + T_{ζ}f.  (2.29)

Sei andersherum f ∈ H^{3/2}(∂Ω) eine Lösung von (2.29). Nun definiert man ψ durch

ψ := h + R_{ζ}f.  (2.30)

Wegen der Eigenschaft e) von S_{ζ}, der eindeutigen Lösbarkeit der Laplacegleichung in Ω′_{ρ} unter den gegebenen Voraussetzungen und schließlich (2.14) für f = 0 ist S_{ζ} injektiv. Nun folgt mit den gleichen Argumenten wie oben

S_{ζ} (Λ_{q}f − \frac{∂ψ}{∂ν_{+}} ) = 0

und damit

Λ_{q}f = \frac{∂ψ}{∂ν_{+}} .

Also ist ψ aufgrund der Eigenschaften von S_{ζ} und D_{ζ} eine Lösung von (2.24)- (2.27).

Lösen des äußeren Problems ist also gleichbedeutend mit dem Lösen der Gleichung (2.29).

Man betrachte nun die Integralgleichung

ψ(x) = h(x) −\begin{matrix}∫ \\ IR & ^{n} \end{matrix} G_{ζ}(x, y)q(y)ψ(y) dy,  (2.31)

h wie oben. Betrachtet man die rechte Seite von (2.31) für eine Lösung ψ ∈ L_{loc}^{2}(IR^{n}), dann ist aufgrund der Eigenschaften von G_{ζ} die Lösung ψ ∈ H_{loc}^{2}(IR^{n}). In der Kugel |y| < ρ gilt

∫_{|y|=ρ} G_{ζ}(x, y)\frac{∂ψ(y)}{∂ν_{y}} − \frac{∂G_{ζ}(x, y)}{∂ν_{y}}ψ(y) dS_{y}

= ∫_{|y|<ρ} G_{ζ}(x, y)q(y)ψ(y) dy + ψ(x) fast überall

→_{ρ→∞} h(x).

Da h harmonisch ist, löst ψ also (2.24)-(2.27).

Andersherum kann man Lösungen ψ ∈ H_{loc}^{2}(IR^{n} \ \overline{Ω}) von (2.24)-(2.27) durch

ψ̃ = P_{q}(ψ|_{∂}Ω)

auf ganz IR^{n} fortsetzen. Aus den Eigenschaften von ψ und G_{ζ}, und da der Träger von q im Abschluß von Ω enthalten ist, ist ψ̃ eine Lösung von (2.31). Aufgrund

der Eigenschaften der Grenzwerte von ψ(x) und \frac{∂ψ(x)}{∂ν} für x gegen ∂Ω ist die

Fortsetzung ψ̃ eindeutig in H_{loc}^{2}(IR^{n}).

Lösen des äußeren Problems ist also gleichbedeutend mit dem Lösen der Gleichung (2.31).

Es sei die Abbildung γ definiert durch

γ : H^{2}(Ω) → H^{3/2}(Ω), γu := u|_{∂}Ω.

Man rechnet leicht nach, daß man T_{ζ} durch

T_{ζ}f = (∫_{Ω} G_{ζ}(x − y)q(y)(P_{q}f)(y) dy)|_{∂Ω}
= γG_{ζ}qIP_{q}f

mit der kompakten Inklusion I : H^{2}(Ω) → L^{2}(Ω) faktorisieren kann. Deshalb ist auch T_{ζ} kompakt. Hier bezeichnet q die Multiplikation mit q(x) als Operator und

G_{ζ}w(x) := ∫_{Ω} G_{ζ}(x − y)w(y) dy.

Nun kann man folgenden Satz beweisen:

Satz 2.4.1 Sei Ω ⊂ IR^{n}, n ≥ 3, ein beschränktes Gebiet mit C^{1,1} Rand, q ∈ L^{∞} reellwertig. Sei 0 kein Eigenwert von −Δ + q in Ω. Dann existiert ein M > 0 derart, daß

f_{ζ} = e^{ix·ζ} − (S_{ζ}Λ_{q} − B_{ζ} − \frac{1}{2}I) f_{ζ}  (2.32)

eine eindeutige Lösung f_{ζ} ∈ H^{3/2}(∂Ω) ∀ζ ∈ l C^{n} mit ζ · ζ = 0 und |ζ| ≥ M existiert. Zusätzlich erfüllt diese Lösung

t(ξ, ζ) = ∫_{∂Ω} e^{−ix·(ξ+ζ)} [Λ_{q} + i(ξ + ζ) · ν] f_{ζ} dS  (2.33)
∀ξ ∈ IR^{n} mit (ξ + ζ)^{2} = 0

Beweis: Man wählt mit den Bezeichnungen aus Korollar 2.3.1 M > c(δ, a)∥q(x)(1 + x^{2})∥_{∞}. Sei g im Nullraum von I − T_{ζ}. Man kann durch ϕ = R_{ζ}g die eindeutige Lösung von (2.24)-(2.27) für h ≡ 0 bestimmen.

Diese kann aber zu der eindeutigen Lösung ϕ̃ von (2.31), h ≡ 0, fortgesetzt werden. Dann gilt

e^{−ix·ζ} ϕ̃ = −K_{ζ} (e^{−ix·ζ} ϕ̃) .

Wegen |ζ| ≥ M folgt also aus Korollar 2.3.1, daß ϕ̃ ≡ 0 und damit g ≡ 0 ist. Also ist I − T_{ζ} injektiv, und da T_{ζ} kompakt ist, sogar invertierbar. Man kann f_{ζ} also durch

f_{ζ} = (\frac{1}{2}I + S_{ζ}Λ_{q} − B_{ζ})^{−1} (e^{ix·ζ})  (2.34)

definieren. Für ζ · ζ = 0 ist e^{ix·ζ} in ganz IR^{n} harmonisch. Deshalb ist f_{ζ} ein Randwert zu einer Lösung von (2.6), die sich für |ζ| groß wie e^{ix·ζ} verhält. □

2.5 Die Randwerte von σ

Satz 2.5.1 Seien Ω, Ω′ und σ wie vorhin. Es gilt:

a)

lim_{|η|→∞} ∥2e^{−ix·η}S_{0}Λ_{σ}e^{ix·η} − σ(x)∥_{L^{2}(∂Ω)} = 0.

b)

lim_{|η|→∞} ∥e^{−ix·η} (σΛ_{1} + Λ_{1}σ − 2Λ_{σ}) e^{ix·η} − \frac{∂σ}{∂ν} (x)∥_{L^{2}(∂Ω)} = 0.

Damit sind σ und \frac{∂σ}{∂ν} auf ∂Ω bekannt.

Beweis:

a) Sei u ∈ C^{2}(Ω) ∩ C(\overline{Ω}) Lösung von ∇ · σ∇u = 0, u|_{∂Ω} = f und v ∈ C^{2}(Ω) ∩ C(\overline{Ω}) beliebig. Der Gauß’sche Integralsatz in Ω liefert

0 = ∫_{Ω} (∇ · σ∇u) v dx

= ∫_{Ω} u (∇ · σ∇v) dx + ∫_{∂Ω} σ (\frac{∂u}{∂ν}v − u\frac{∂v}{∂ν} ) dS.

Nun setzt man G_{ζ}^{ε} aus dem Beweis zu Lemma 2.2.2 für v ein um

S_{0}Λ_{σ}f − \frac{1}{2}σf = B_{0}(σf) − (∫_{Ω} u(y)∇_{y}G_{0}(x, y) · ∇σ(y) dy)|_{∂Ω} (2.35)

zu zeigen. Die rechte Seite der Gleichung (2.35) läßt sich, wie in Nachman [N], Kapitel 4 unter zu Hilfenahme von Fabes, Jodeit und Rivère [F-J-R] gezeigt, gegen C|η|^{−1/2} abschätzen.

b) Man zeigt leicht, daß für u_{η} und v_{η} mit

Δv_{η} = 0 in Ω, v_{η}|_{∂Ω} = e^{ix·η},
∇ · σ∇u = 0 in Ω, u_{η}|_{∂Ω} = e^{ix·η},
∥u_{η} − v_{η}∥_{H^{1}(Ω)} ≤ c|η|^{−1/2}

gilt. Außerdem ist S_{0}, wie in Abschnitt 2.4 gezeigt, injektiv und, da Ω beschränkt ist, sogar invertierbar. Die Behauptung folgt ähnlich wie im Beweis zu Teil a) (siehe Nachman [N], Kapitel 4, Satz 4.3, Lemma 4.4 und Anhang A).

□

2.6 Die eindeutige Lösung des inversen Problems

Satz 2.6.1 Sei Ω ⊂ IR^{n}, n ≥ 3, ein beschränktes Gebiet mit C^{1,1} Rand, und sei Ω′ = IR^{n} \ Ω zusammenhängend. Sei σ ∈ C^{1,1}(\overline{Ω}) reellwertig mit positiver unterer Schranke. Dann ist σ eindeutig durch Λ_{σ} bestimmt und kann folgendermaßen rekonstruiert werden:

1. Berechne die Randwerte σ|_{∂Ω} und \frac{∂σ}{∂ν}|_{∂Ω} aus Λ_{σ}.

2. Löse (2.32) um t(ξ, ζ) aus Λ_{q} für |ζ| ≥ M, wie in Satz 2.4.1 gezeigt, zu berechnen.

3. Benutze Formel (2.8) aus Satz 2.1.9 oder Formel (2.9) um q aus t(ξ, ζ) zu ermitteln.

4. Berechne σ aus q durch Lösen der partiellen Differentialgleichung

−Δσ^{1/2} + qσ^{1/2} = 0 in Ω,

dabei ist σ^{1/2}|_{∂Ω} aus Satz 2.5.1 bekannt.

Numerische Verfahren

3.1 Übersicht

In diesem Kapitel stelle ich einige numerische Verfahren zur Rekonstruktion der Leitfähigkeit σ vor. Es gibt noch andere Verfahren, die ich an dieser Stelle nicht nennen werde, weil sie nicht stabil sind oder nur aufgrund sehr großer Einschränkungen funktionieren. Trotzdem habe ich die entsprechenden Originalarbeiten im Literaturverzeichnis erwähnt.

Das erste Verfahren ist die gefilterte Rückprojektion nach Barber und Brown. Dieses klassische Verfahren stammt aus der Computer-Tomographie (CT). Die Methode beruht auf der Arbeit von Radon aus dem Jahre 1917.

Danach betrachte ich einen Ansatz aus der Variationsrechnung. Hier minimiert man ein nichtlineares Funktional. Dazu wird eine quadratische Approximation des Funktionals hergeleitet. Die numerische Implementierung erfolgt mit finiten Elementen.

Beim Layer-Stripping rekonstruiert man die Leitfähigkeit nur auf dem Rand eines Gebietes. Um aber Daten aus dem Innern zu erhalten, stellt man das Problem in Untergebieten. Numerisch wird dieses Verfahren mit Fourierreihen implementiert.

In dem letzen Verfahren in diesem Kapitel will man stückweise konstante Leitfähigkeiten rekonstruieren. Man sucht also nicht die Leitfähigkeit σ an sich, sondern das Gebiet, in dem σ einen bestimmten Wert annimmt. Dieses Verfahren stelle ich vor, weil ich diesen komplett anderen Ansatz interessant finde.

Schließlich wird in den anschließenden beiden Kapiteln ein iteratives Verfahren hergeleitet. Dazu werde ich auch eigene numerische Ergebnisse präsentieren.

3.2 Lineare Rückprojektion

3.2.1 Die verallgemeinerte Radontransformation

Als Analogon zum Skalarpordukt in der klassischen Radontransformation führt man für die verallgemeinerte Radontransformation eine Funktion ϕ(x, θ) definiert auf X × (IR^{n} \ {0}) mit X ⊂ IR^{n} ein. Dabei soll ϕ folgende Bedingungen erfüllen:

(i) ϕ(x, θ) ist eine reellwertige C^{∞}-Funktion auf X × (IR^{n} \ {0}).

(ii) ϕ(x, θ) ist homogen vom Grad eins bzgl. θ:

ϕ(x, λθ) = λϕ(x, θ) ∀λ ∈ IR.

(iii) ∇_{x}ϕ verschwindet nirgendwo in X × (IR^{n} \ {0}).

(iv) Für

h(x, θ) := det (\frac{∂^{2}ϕ(x, θ)}{∂x^{j}∂θ^{k}})

gilt: h(x, θ) > 0 in X × (IR^{n} \ {0}).

Zu gegebenem ϕ konstruiert man nun eine Familie H von (n − 1)-dimensionalen Oberflächen H_{s,ω} im Gebiet X mit

H_{s,ω} = {x ∈ X|s = ϕ(x, ω), s ∈ IR^{1}, ω ∈ S^{n−1}} .

Aus (ii) folgt, daß H_{s,ω} = H_{−s,−ω} gilt. Also ist H isomorph zu (IR^{1} × S^{n−1})/ZZ_{2}. Zu gegbenem x ∈ X sei Y_{x} die Menge aller Hyperoberflächen der Familie H, die durch den Punkt x verlaufen,

Y_{x} = {[ϕ(x, ω), ω] ∈ H|ω ∈ S^{n−1}} .

Die Funktion ϕ muß noch zwei zusätzliche Bedingungen erfüllen:

(v) Aus H_{s,ω} = H_{s′,ω′} folgt schon [s, ω] = [s′, ω′].

(vi) Ist Y_{x} = Y_{x′}, dann folgt x = x′.

Nun definiert man die verallgemeinerte Radontransformation durch

(Ru)([s, ω]) = ∫_{H_{s,ω}} u(x)a(x, ω)Ω  (3.1)

mit einer positiven Gewichtsfunktion a(x, ω). Dabei ist die Differentialform

Ω = ∑_{j=1}^{n}(−1)^{j−1} \frac{∂ϕ(x, ω)/∂x_{j}}{|∇_{x}ϕ(x, ω)|^{2}} dx_{1} ∧ … ∧ dx_{j−1} ∧ dx_{j+1} ∧ … ∧ dx_{n}  (3.2)

so gewählt, daß

d_{x}ϕ(x, ω) ∧ Ω = dx_{1} ∧ … ∧ dx_{n}  (3.3)

gilt.

Der zu R duale Operator R^{∗} lautet:

R^{∗}v(y) = ∫_{|ω|=1} b(y, ω)v(s, ω)|_{s=ϕ(y,ω)} dS_{ω}  (3.4)

mit

b(x, ω) = \frac{h(x, ω)}{a(x, ω)}.

In Beylkin [B1] wird nun gezeigt, daß folgender Satz gilt:

Satz 3.2.1 Sei K der Integraloperator definiert als Faltung mit dem Kern

K(s) = \frac{1}{2(2π)^{n}} ∫_{−∞}^{∞} |r|^{n−1}e^{irs} dr.  (3.5)

Dann gilt für F = R^{∗}KR:

F = I + T  (3.6)

mit einem kompakten Operator T.

Bemerkung 3.2.2 K(s) existiert nicht im klassischen Sinne, sondern als Distribution, d.h., man integriert K(s) gegen eine schnell fallende C^{∞}-Funktion (Schwartz’scher Raum S, siehe Anhang B) und vertauscht die Integrationsreihenfolge. Die Existenz folgt nun aus der Glattheit der Testfunktion.

Mit Hilfe dieses Satzes kann man nun die Radontransformation approximativ invertieren, indem man ganz einfach den kompakten Operator T vernachlässigt.

Einen anderen Ansatz die Radontransformation zu invertieren wird in Berenstein [B-C-T] gegeben. In diesem Artikel wird gezeigt, daß für den Laplace-Beltrami-Operator Δ in gewissen Koordinaten

−\frac{1}{4π} ΔSR^{∗}R = I

gilt. Dabei sei der Operator S als Faltung mit dem Kern

S̃(r) = coth(r) − 1

definiert.

3.2.2 Linearisierung

Man betrachte das in Kapitel 1 hergeleitete mathematische Modell:

∇ · (σ∇u) = 0 in Ω,
σ \frac{∂u}{∂ν} = f auf ∂Ω.  (3.7)

Zur Vereinfachung linearisiert man dieses Problem um u, σ: Für eine kleine Störung η in σ erfüllt nach Vernachlässigung von Termen zweiter Ordnung die zugehörige Störung w in u die Differentialgleichung

∇ · (σ∇w) = −∇ · (η∇u) in Ω,
σ \frac{∂w}{∂ν} = −η \frac{∂u}{∂ν} auf ∂Ω.  (3.8)

Im weiteren Verlauf betrachtet man das Problem noch unter folgenden Voraussetzungen:

1. Ω = {x ∈ IR^{2} : |x| < 1}.

2. σ ≡ 1. u sei eine Dipollösung bzgl. ω, d.h. u löst

Δu = 0 in Ω
\frac{∂u}{∂ν} = −πδ′_{ω} auf ∂Ω.  (3.9)

Dieses entspricht in der Praxis zwei aufeinander folgenden, entgegengesetzt gepolten Elektroden.

3. Da man Effekte aus dem Innern sucht, nimmt man an, daß η = 0 in einer Umgebung von ω ist.

Unter diesen Annahmen vereinfacht sich (3.8) zu

Δw = −∇η · ∇u in Ω,
\frac{∂w}{∂ν} = 0 auf ∂Ω.  (3.10)

Die derart vereinfachte Aufgabe besteht nun darin, zu gegebenen Störungen w von verschiedenen Dipollösungen u ein η zu suchen, so daß (3.10) für alle Dipol-Positionen ω erfüllt wird.

3.2.3 Die Poincaré-Abbildung

Sei Ω die Einheitskugel in IR^{2}, ω ∈ ∂Ω. Zu gegebenem ω definiert man U und V als

x′_{1} = ω^{⊥} · x,
x′_{2} = 1 − ω · x,
U = \frac{x′_{1}}{x′_{1}^{2} + x′_{2}^{2}} ,  (3.11)
V = \frac{x′_{2}}{x′_{1}^{2} + x′_{2}^{2}} .

Santosa nennt die Abbildung x ↦ (−U, V) die Poincaré-Abbildung. Wie man leicht nachrechnet, bildet sie Ω konform auf P = \{(x, y), x > \frac{1}{2}\} ab.

Außerdem ist U die Dipollösung von (3.9) in der Einheitskugel.

Betrachtet man nun das Problem (3.10) in (U, V)-Koordinaten, so erhält man folgende Gleichung:

ΔW = −\frac{∂η}{∂U} in P,
\frac{∂W}{∂V} = 0 auf ∂P.  (3.12)

W und η hängen in diesen Koordinaten noch von der Dipol-Position ω ab.

3.2.4 Anwendung der Radontransformation

In der EIT wählt man die in Abschnitt 3.2.1 eingeführten Funktionen ϕ und a folgendermaßen:

ϕ(x, ξ) = |ξ|U(x, ξ/|ξ|), x ∈ Ω, ξ ∈ IR^{2} \ {0},  (3.13)
a(x, ω) = |∇_{x}ϕ(x, ω)|^{2}.  (3.14)

Man rechnet leicht nach, daß ϕ die Eigenschaften (i)-(vi) erfüllt und sich die Gewichtsfunktion b(x, ω) von R^{∗} zu

b(x, ω) = \frac{h(x, ω)}{a(x, ω)} = 2V (x, ω) − 1, x ∈ Ω, |ω| = 1,  (3.15)

errechnet. Der Beweis kann in Santosa [S] nachgelesen werden.

Satz 3.2.3 Seien W, U, V und ω wie oben. Dann ist in der Notation von Abschnitt 3.2.1

KR(η) ≈ −\frac{1}{2π} ( \frac{∂W}{∂U} |_{V =1/2}) (s, ω).

Beweis: In IR^{2} ist

γ_{2}(x, y) = \frac{1}{2π} log \frac{1}{\sqrt{x^{2} + y^{2}}}

eine Fundamentallösung der Gleichung −Δu = f. Nun betrachte man

G̃(x, y) = −\frac{∂γ_{2}}{∂x} = \frac{x}{2π(x^{2} + y^{2})}.

Mittels der Theorie über distributive Ableitungen erhält man so eine Fundamentallösung zu −Δu = \frac{∂f}{∂x}. Die Greensche Funktion in P muß nun noch

\frac{∂G(u, v, U, V)}{∂V}|_{V = \frac{1}{2}} = 0 ∀ u, U ∈ IR, v ∈ [\frac{1}{2}; ∞[ ,

erfüllen und ergibt sich zu

G(u, v, U, V) = \frac{1}{2} (G̃(U − u, V − v) + G̃(U − u, 1 − V − v))
= −\frac{U − u}{4π} (\frac{1}{(U − u)^{2} + (V − v)^{2}} + \frac{1}{(U − u)^{2} + (1 − V − v)^{2}} ) .

Nun konstruiert man Lösungen von (3.12) durch

W (U, V, ω) = ∫_{P} G(u, v, U, V)η(u, v, ω) du dv.

Differenziert man W nach U, so ergibt Einsetzen von G

− ( \frac{∂W}{∂U} |_{V =1/2}) (s, ω) = \frac{1}{2π} ∫_{1/2}^{∞} ∫_{−∞}^{∞} \frac{(1/2 − v)^{2} − (s − u)^{2}}{((s − u)^{2} + (1/2 − v)^{2})^{2}} η(u, v, ω) du dv.

(3.16)

Man betrachte nun für q > 0

g(p) = \frac{q^{2} − p^{2}}{(q^{2} + p^{2})^{2}} ,
h(r) = \{ \begin{matrix}r & , & r ≥ 0 \\ 0 & , & r < 0 \end{matrix}
und f(r) = |r|e^{−q|r|} = h(r)e^{−qr} + h(−r)e^{qr}.

KAPITEL 3. NUMERISCHE VERFAHREN 31 Die Fouriertransformation von g bestimmt man durch partielle Integration:

∫_{−∞}^{∞} h(r)e^{−qr+irp} dr = lim_{R→0} ∫_{0}^{R} re^{(ip−q)r} dr
= lim_{R→0} (\frac{Re^{(ip−q)R}}{ip − q} − \frac{e^{(ip−q)R} − 1}{(ip − q)^{2}})
= \frac{1}{(ip − q)^{2}} .

Analog folgt

∫_{−∞}^{∞} h(−r)e^{qr+irp} dr = \frac{1}{(ip + q)^{2}} .

Damit gilt

f̃(p) = \frac{2}{\sqrt{2π}} \frac{q^{2} − p^{2}}{(q^{2} + p^{2})^{2}}
oder ĝ(r) = \sqrt{2π} \frac{|r|}{2}e^{−q|r|}.

Mit q = v − \frac{1}{2} kann man nun (3.16) als

− ( \frac{∂W}{∂U} |_{V =1/2}) (s, ω) = \frac{1}{2π} ∫_{1/2}^{∞} g ∗ η(·, v, ω) dv

schreiben und Fouriertransformation bezüglich s führt zu

− ( \frac{∂W}{∂U} |_{V =1/2})^{∧} (r, ω) = \frac{|r|}{2} ∫_{1/2}^{∞} e^{−q|r|}η̂(r, v, ω) dv.

Vernachlässigt man nun den Term e^{−q|r|} unter dem Integral, dann folgt nach Vertauschung der Integrationsreihenfolge

− ( \frac{∂W}{∂U} |_{V =1/2})^{∧} (r, ω) ≈ \frac{|r|}{2} ( ∫_{1/2}^{∞} η(·, v, ω) dv)^{∧} (r)
= 2πK ( ∫_{1/2}^{∞} η(·, v, ω) dv) (s).

Der Operator K sei wie in Satz 3.2.1 für n = 2 gegeben.

KAPITEL 3. NUMERISCHE VERFAHREN 32 Variablentransformation ergibt nun

∫_{1/2}^{∞} η(s, v, ω) dv = ∫_{H_{s,ω}} η|∇_{x}U(x, ω)| dS
= R(η)(s, ω).
□

Korollar 3.2.4 Zusammenfassend gilt mit den Bezeichnungen aus Satz 3.2.1:

B(x) := \frac{1}{2π} ∫_{|ω|=1} ( \frac{∂W}{∂U} |_{V =1/2}) (U(x, ω), ω)(1 − 2V (x, ω)) dS_{ω}
≈ R^{∗}KR (η) (x)
= (I + T) (η) (x).

Damit ist B(x) eine Approximation von η(x).

3.3 Variationsrechnung

3.3.1 Formulierung der Minimierungsaufgabe

Satz 3.3.1 Minimieren des Funktionals

I(σ, u_{1}, …, u_{N}, v_{1}, …, v_{N}) = ∑_{i=1}^{N} \frac{1}{2} ∫_{Ω} |σ^{1/2}∇v_{i} − σ^{1/2}∇u_{i}|^{2} dx  (3.17)

mit den Randbedingungen

σ \frac{∂u_{i}}{∂ν} = f_{i},  v_{i}|_{∂Ω} = g_{i}  (3.18)

über alle Argumente löst (1.6), (1.7) und das inverse Problem.

Beweis: Es gilt:

I = ∑_{i=1}^{N} [\frac{1}{2} ∫_{Ω} σ |∇v_{i}|^{2} dx + \frac{1}{2} ∫_{Ω} σ |∇u_{i}|^{2} dx + ∫_{∂Ω} f_{i}g_{i} dS] .  (3.19)

Lösen der Variationsgleichung a(u, w) = F (w) ∀w (vgl. Abschnitt 1.3) ist äquivalent zur Bestimmung von inf_{u} (a(u, v) − 2F (u)). Falls u_{i} und v_{i} schwache Lösungen von (1.6) und (1.7) sind und übereinstimmen, ist I(σ, u_{1}, …, u_{N}, v_{1}, …, v_{N}) = 0, also minimal. Ist I aber minimal, so müssen u_{i} und v_{i} Lösungen von (1.6) und (1.7) sein und übereinstimmen. Also ist σ Lösung des inversen Problems. □

3.3.2 Numerische Methoden

Bezeichnungen:

J_{i} := −σ∇u_{i},
E_{i} := −∇v_{i},
E := −(∇v_{1}, …, ∇v_{N}),  (3.20)
J := (J_{1}, …, J_{N}),
|E|^{2} = ∑_{i=1}^{N} |∇v_{i}|^{2},

|J|^{2} = ∑_{i=1}^{N} |J_{i}|^{2},  (3.21)

E · J = ∑_{i=1}^{N} ∇v_{i} · J_{i}.

(3.17) läßt sich dadurch in folgender Form schreiben:

I(σ; E; J) = \frac{1}{2} ∫_{Ω} |σ^{1/2}E − σ^{−1/2}J|^{2} dx.  (3.22)

Die Gutman-Methode

Die erste in Kohn und McKenny [K-M] vorgestellte Methode kann man die Gutman-Methode nennen. Kohn und McKenny bezeichnen sie als die alternating direction implicit Methode, was aber eventuell im Widerspruch zur Literatur steht. Das Verfahren ist in zwei Schritte unterteilt:

1. Für ein festes vorgegebenes σ_{n} minimiere I über v_{1}, …, v_{N} und J_{1}, …, J_{N}. Man erhält also E(σ_{n}) und J(σ_{n}).

2. Mit E(σ_{n}) und J(σ_{n}) aus dem 1. Schritt minimiere I über σ, um eine neue Approximation σ_{n+1} zu berechnen: Wie oben gezeigt, läßt sich I(σ; E; J) auch als

I(σ; E; J) = ∫_{Ω} σ|E|^{2} + \frac{1}{σ} |J|^{2} dx + Randintegral unabhängig von σ

schreiben. Das Integral wird minimal, wenn der Integrand minimal wird. Also:

σ_{n+1} = (∑_{i=1}^{N} |J_{i}|^{2})^{\frac{1}{2}} (∑_{i=1}^{N} |∇v_{i}|^{2})^{− \frac{1}{2}} = \frac{|J(σ_{n})|}{|E(σ_{n})|}.  (3.23)

Man erhält so eine Folge I(σ_{n+1}; E(σ_{n}); J(σ_{n})), die gegen das Minimum von I konvergiert. Gutman [G] zeigt, daß sich diese Folge folgendermaßen einschachteln läßt:

I(σ_{n+1}; E(σ_{n+1}); J(σ_{n+1})) ≤ I(σ_{n+1}; E(σ_{n}); J(σ_{n})) ≤ I(σ_{n}; E(σ_{n}); J(σ_{n})).

Die modifizierte Newtonmethode

Um die Konvergenzgeschwindigkeit zu erhöhen, schlagen Kohn und McKenny in [K-M] für den zweiten Schritt einen quadratischen Algorithmus vor. Dabei ist jedoch nicht immer die positive Definitheit der 2. Ableitung des betrachteten Funktionals gegeben, so daß dort die Newtonmethode modifiziert werden muß.

Streng genommen wird (3.22) nicht über alle Variablen minimiert, sondern man behandelt σ als bestimmt durch E und J. I wird also unter Zuhilfenahme von (3.23) minimiert.

Die quadratische Approximation von I ermittelt man durch Auswerten von (3.22) in σ = σ_{0} + δσ, E = E_{0} + δE und J = J_{0} + δJ bis zur 2. Ordnung. Dabei müssen für i = 1, …, N folgende Randbedingungen erfüllt sein:

δv_{i}|_{∂Ω} = 0,
δJ_{i} · ν|_{∂Ω} = 0,  (3.24)
∇ · δJ_{i} = 0.

Dann ist

I(σ, E, J) = I(σ_{0}, E_{0}, J_{0}) + δI + \frac{1}{2}δ^{2}I + O (|δE|^{2} + |δE|^{2} + (δσ)^{2}) (3.25)

mit

\frac{1}{σ} = \frac{1}{σ_{0} + δσ} = \frac{1}{σ_{0}} − \frac{δσ}{σ_{0}^{2}} + \frac{(δσ)^{2}}{σ_{0}^{3}} + O((δσ)^{3}),  (3.26)
δI = \frac{1}{2} ∫_{Ω} [(|E_{0}|^{2} − \frac{1}{σ_{0}^{2}} |J_{0}|^{2}) δσ + 2σ_{0}E_{0} · δE + \frac{2}{σ_{0}} J_{0} · δJ] dx,  (3.27)

δ^{2}I = ∫_{Ω} [\frac{1}{σ_{0}^{3}} |J_{0}|^{2}(δσ)^{2} + 2 (E_{0} · δE − \frac{1}{σ_{0}^{2}} J_{0} · δJ) δσ + σ_{0}|δE|^{2} + \frac{1}{σ_{0}} |δJ|^{2}] dx.

Aus (3.23) ergibt sich nun

σ_{0} = \frac{|J_{0}|}{|E_{0}|},  (3.28)
δσ = \frac{J_{0} · δJ − σ_{0}^{2}E_{0} · δE}{|E_{0}||J_{0}|} + O(|δE|^{2} + |δJ|^{2}).  (3.29)

Einsetzen von (3.28) und (3.29) in (3.27) führt zu

δI′ = ∫_{Ω} [σ_{0}E_{0} · δE + \frac{1}{σ_{0}} J_{0} · δJ] dx,  (3.30)

δ^{2}I′ = ∫_{Ω} σ_{0}|δE|^{2} + \frac{1}{σ_{0}} |δJ|^{2} − \frac{σ_{0}}{|E_{0}|^{2}} (E_{0} · δE − \frac{1}{σ_{0}^{2}} J_{0} · δJ)^{2} dx. (3.31)

Da (3.31) im Allgemeinen nicht positiv definit ist, ersetzen Kohn und Mc Kenny δ^{2}I′ durch

δ^{2}I′′ = ∫_{Ω} σ_{0}|δE|^{2} + \frac{1}{σ_{0}} |δJ|^{2} − \frac{1}{1 + ε} \frac{σ_{0}}{|E_{0}|^{2}} (E_{0} · [δE − \frac{1}{σ_{0}} δJ])^{2} dx. (3.32)

Es gilt für |J − σ_{0}E| → 0 und ε → 0, daß

|J − \frac{σ}{\sqrt{1 + ε}}E| → 0  (3.33)

und deshalb

δ^{2}I′′ − δ^{2}I′ = ∫_{Ω} \frac{σ_{0}}{|E_{0}|^{2}} [(E_{0} · δE − \frac{1}{σ_{0}^{2}} J_{0} · δJ)^{2} − \frac{1}{1 + ε} (E_{0} · [δE − \frac{1}{σ_{0}} δJ])^{2}] dx

= ∫_{Ω} \frac{σ_{0}}{|E_{0}|^{2}} [\frac{1}{σ_{0}^{2}} ((J_{0} − \frac{σ_{0}}{\sqrt{1 + ε}}E_{0}) · δJ)

((1 + \frac{1}{\sqrt{1 + ε}}) E_{0} · δE − \frac{1}{σ_{0}^{2}} (\frac{σ_{0}}{\sqrt{1 + ε}}E + J_{0}) · δJ)] dx

→ 0.

Nun folgt

\frac{1}{1 + ε} \frac{σ_{0}}{|E_{0}|^{2}} (E_{0} · [δE − \frac{1}{σ_{0}} δJ])^{2} ≤ \frac{σ_{0}}{1 + ε} |δE_{0} − \frac{1}{σ_{0}} δJ|^{2}
= \frac{1}{1 + ε} (σ_{0}|δE|^{2} + \frac{1}{σ_{0}} |δJ|^{2} − 2δE · δJ) ,

und partielles Integrieren unter Verwendung von (3.20) und (3.24) ergibt

∫_{Ω} δE · δJ dx = 0.

Deshalb ist δ^{2}I′′ positiv definit:

δ^{2}I′′ ≥ \frac{ε}{1 + ε} ∫_{Ω} σ_{0}|δE|^{2} + \frac{1}{σ_{0}} |δJ|^{2} dx > 0.  (3.34)

Man kann also

I ≈ I_{0} + δI′ + \frac{1}{2}δ^{2}I′′

als Approximation der nach dem dritten Glied abgebrochenen Taylorentwicklung von I betrachten. Da I_{0} unabhängig von δE und δJ ist, genügt es lediglich

δI′ + \frac{1}{2}δ^{2}I′′

zu minimieren. Nach Diskretisierung mittels finiter Elemente ergibt sich δI′ zu einem linearen Funktional für die Knotenvektoren, und aus δ^{2}I′′ wird eine Bilinearform. So gelangt man zu einem linearen Gleichungssystem für die Knotenvektoren.

3.4 Layer Stripping

Grundidee der folgenden Methode ist es, das Problem in zwei Schritte aufzuteilen:

1. Berechnung der Leitfähigkeit σ auf dem Rand eines beschränkten Gebietes Ω.

2. Umrechnung des Problems auf ein kleineres Gebiet Ω̃, das ganz in Ω enthalten ist, also \overline{Ω̃} ⊂ Ω.

Auf diese Art und Weise kann man die Leitfähigkeit σ auf gewissen Kurven bestimmen.

3.4.1 Rekonstruktion auf dem Rand

Man betrachtet die etwas allgemeinere partielle Differentialgleichung

\begin{aligned}∇ · σ∇u & = & f & in & Ω, \\ σ \frac{∂u}{∂ν} & = & j & auf & ∂Ω \end{aligned}  (3.35)

zusammen mit den Nebenbedingungen

\begin{aligned}∫_{∂Ω} u dS & = & 0, \\ ∫_{Ω} f dx & = & ∫_{∂Ω} j dS. \end{aligned}  (3.36)

Zur Formulierung des Problems und des Lösungsweges verwendet man einige Operatoren:

Sei s ≥ 0. Für eine Lösung u von (3.35) definiert man

\begin{aligned}γ & : & H^{s+2}(Ω) → H^{s+3/2}(∂Ω), \\ γ & : & u ↦ u|_{∂Ω}, \\ T_{σ} & : & H^{s}(Ω) × H^{s+1/2}(∂Ω) → H^{s+2}(Ω), \\ T_{σ}(f, j) & := & u, \\ R & = & γT_{σ}|_{\{0\}×H^{s+1/2}(∂Ω)}. \end{aligned}  (3.37)

Schließlich bezeichnet man mit T die Abbildung T_{σ} für σ ≡ 1 und entsprechend R_{0} = γT |_{{0}×H^{s+1/2}(∂Ω)}.

KAPITEL 3. NUMERISCHE VERFAHREN 38 Satz 3.4.1 Sei σ hinreichend glatt (z.B. σ ∈ H^{s+2}(Ω)). Dann läßt sich der Operator T_{σ} in

T_{σ} = \frac{1}{σ} T + T_{R}  (3.38)

aufteilen mit T_{R} : H^{s}(Ω) × H^{s+1/2}(∂Ω) → H^{s+3}(Ω), s ≥ 0. Beweis: Seien u = T_{σ}(f, j) und u_{0} = T (f, j). Es ist zu zeigen, daß

u_{R} = u − \frac{1}{σ} u_{0} ∈ H^{s+3}(Ω).

u_{R} erfüllt wegen ∇ · σ∇u = 0 und Δu_{0} = 0 die Gleichungen

∇ · σ∇u_{R} = ∇ · u_{0}∇η,
σ \frac{∂u_{R}}{∂ν} = u_{0} \frac{∂η}{∂ν}

mit η = log σ. Es folgt

u_{R} = T_{σ}Bu_{0} = T_{σ}BT (f, j)  (3.39)

mit

Bψ = (∇ · (ψ∇η), ψ \frac{∂η}{∂ν} |_{∂Ω}).

Wenn nun η hinreichend glatt ist, so gilt

B : H^{s+2}(Ω) → H^{s+1}(Ω) × H^{s+3/2}(∂Ω).  (3.40)

Zusammen mit (3.37) ist nun u_{R} ∈ H^{s+3}(Ω). Setze T_{R} := T_{σ}BT. □ Aus dem gerade gezeigten Satz folgt nun für die Randwerte

R = γT_{σ}|_{{0}×H^{s+1/2}(∂Ω)}
= γ \frac{1}{σ} T |_{\{0\}×H^{s+1/2}(∂Ω)} + γT_{R}|_{\{0\}×H^{s+1/2}(∂Ω)}  (3.41)
= \frac{1}{σ} |_{∂Ω} R_{0} + R_{R}
und R_{R} : H^{s+1/2}(∂Ω) → H^{s+5/2}(∂Ω).

Die Bedeutung dieser Formel sieht man am besten an einem kleinen Beispiel.

Beispiel

In dem Gebiet Ω := {x ∈ IR^{2}, |x| < 1} kann man den Operator R_{0} folgendermaßen berechnen:

Man betrachtet Lösungen von

Δu_{0} = 0 in Ω,
\frac{∂u}{∂r} = j auf ∂Ω

in Polarkoordinaten, die sich in der Form

u_{0}(r, θ) = ∑_{{m=−∞ \atop m≠0}}^{∞} c_{m}r^{|m|}e^{imθ}

schreiben lassen. Die Randbedingung wird so zu

∑_{{m=−∞ \atop m≠0}}^{∞} |m|c_{m}e^{imθ} = j(θ).

Schreibt man j als Fourierreihe und sei j_{m} der m-te Fourierkoeffizient von j, dann folgt daraus

c_{m} = \frac{1}{|m|}j_{m}.

Zusammenfassend gilt

R_{0} (e^{in·}) (θ) = \frac{1}{|n|}e^{inθ}, ∀ n ∈ ZZ, n ≠ 0.  (3.42)

Den Fall n = 0 kann man wegen (3.36) ausschließen.

Aufbauend auf Satz 3.4.1 kann man damit folgenden Satz über die Fourierkoef-

fizienten von \frac{1}{σ(1, θ)} beweisen:

Satz 3.4.2 Seien Ω ⊂ IR^{2} der Einheitskreis und σ hinreichend glatt mit positiver unterer Schranke. Dann ist

lim_{|n|→∞} |n|R_{n+k,n} = ρ_{k}, n, k ∈ ZZ, n ≠ 0  (3.43)

mit

R_{m,n} = ⟨e^{imθ}, R (e^{in·}) (θ)⟩ ,  (3.44)
ρ_{n} = \frac{1}{2π} ∫_{0}^{2π}e^{−inθ} \frac{1}{σ(1, θ)} dθ.

Beweis: Aus (3.42) und (3.44) folgt

R_{m,n} = \frac{1}{|n|}ρ_{m−n} + ⟨e^{imθ}, R_{R} (e^{in·}) (θ)⟩ .  (3.45)

Es gilt:

m^{2} ⟨e^{imθ}, R_{R} (e^{in·}) (θ)⟩ = m^{2} \frac{1}{2π} ∫_{0}^{2π}e^{−imθ}R_{R} (e^{in·}) (θ) dθ

= −\frac{1}{2π} ∫_{0}^{2π} \frac{∂^{2}e^{−imθ}}{∂θ^{2}} R_{R} (e^{in·}) (θ) dθ

Durch Anwenden der Cauchy-Schwartz’schen-Ungleichung erhält man

|m^{2} ⟨e^{imθ}, R_{R} (e^{in·}) (θ)⟩| ≤ \frac{1}{2π} ∥e^{im·}∥_{H^{−(s+5/2)+2}(S^{1})} ∥R_{R} (e^{in·})∥_{H^{s+5/2}(S^{1})} .

mit s ≥ 0. Nach Satz 3.4.1 ist R_{R} : H^{s+1/2}(S^{1}) ⟶ H^{s+5/2}(S^{1}) beschränkt, deshalb gilt

|m^{2} ⟨e^{imθ}, R_{R} (e^{in·}) (θ)⟩| ≤ C ∥e^{im·}∥_{H^{−(s+1/2)}(S^{1})} ∥e^{in·}∥_{H^{s+1/2}(S^{1})}
≤ C \frac{(1 + n^{2})^{(2s+1)/4}}{(1 + m^{2})^{(2s+1)/4}} .

Für m = n + k folgt dann

|⟨e^{imθ}, R_{R} (e^{in·}) (θ)⟩| ∈ O(n^{−2}).
□

3.4.2 Unteroberflächensynthese von Daten im Innern

Somersalo, Cheney, Isaacson und Isaacson [S-C-I-I] entwickeln in sphärischen Gebieten eine Differentialgleichung für R. Diese nennen sie die Riccatigleichung. Sei Ω_{t} = {x ∈ IR^{n}, |x| < t}, 0 < t ≤ 1, und Ω = Ω_{1}. Für jedes t kann man nun das Randwertproblem lösen:

\begin{aligned}∇ · σ∇u & = & 0 & für |x| < t, \\ σ \frac{∂u}{∂ν} & = & j & für |x| = t. \end{aligned}  (3.46)

Damit definiert man eine Familie von Operatoren

\begin{aligned}R(t) & : & H^{s}(S^{n−1}) & → & H^{s+1}(S^{n−1}), s ≥ \frac{1}{2}, \\ R(t) & : & j & ↦ & u|_{∂Ω_{t}}. \end{aligned}  (3.47)

Satz 3.4.3 Die Operatoren R(t), 0 < t ≤ 1, bilden eine linksstetige Familie bzgl. t in der starken Operatortopologie von L(H^{s}(S^{n−1}), H^{s+1}(S^{n−1})), s ≥ \frac{3}{2}. Zusätzlich existiert der starke Limes

lim_{ε→0+} \frac{1}{ε} (R(t) − R(t − ε)) =: \frac{d}{dt_{−}} R(t)

in L(H^{s}(S^{n−1}), H^{s+1}(S^{n−1})), und es gilt

\frac{d}{dt_{−}} R(t) = \frac{1}{σ} + \frac{n − 1}{t}R(t) + R(t)∇_{τ} · σ∇_{τ}R(t).  (3.48)

Dabei ist ∇_{τ} ·σ∇_{τ} der tangentiale Teil von ∇·σ∇ entlang der Oberfläche |x| = t.

Beweis:

1. Stetigkeit: Sei u eine Lösung von (3.46).

(R(t) − R(t − ε))j = (R(t)j − u|_{|x|=t−ε}) + (u|_{|x|=t−ε} − R(t − ε)j)
= (u|_{|x|=t} − u|_{|x|=t−ε}) + R(t − ε) (σ \frac{∂u}{∂r}|_{|x|=t−ε} − σ \frac{∂u}{∂r}|_{|x|=t}) .

Die Abbildung u ↦ u|_{|x|=r} ist für r ≤ t stetig von H^{s+3/2}(Ω_{t}) nach H^{s+1}(S^{1}). Deshalb gilt

∥u|_{|x|=t} − u|_{|x|=t−ε}∥_{H^{s+1}(S^{n−1})} →_{ε→0} 0

und entsprechend

∥ \frac{∂u}{∂r}|_{|x|=t} − \frac{∂u}{∂r}|_{|x|=t−ε}∥_{H^{s}(S^{n−1})} →_{ε→0} 0.

Mit R(t) ist auch R(t−ε) für kleine ε beschränkt. Daraus folgt die Stetigkeit.

2. Differenzierbarkeit: Es ist zu zeigen

\frac{1}{ε}(u|_{|x|=t} − u|_{|x|=t−ε}) ⟶_{ε→0} \frac{∂u}{∂r}|_{|x|=t} in H^{s}(S^{n−1}) und

\frac{1}{ε} (σ \frac{∂u}{∂r}|_{|x|=t} − σ \frac{∂u}{∂r}|_{|x|=t−ε}) ⟶_{ε→0} \frac{∂}{∂r}σ \frac{∂u}{∂r}|_{|x|=t} in H^{s−1}(S^{n−1}).

Für t = 1 gilt

∥\frac{1}{ε}(u|_{|x|=1} − u|_{|x|=1−ε}) − \frac{∂u}{∂r}|_{|x|=1}∥_{H^{s}(S^{n−1})}

= ∥\frac{1}{ε} ∫_{1−ε}^{1} \frac{∂u}{∂r}|_{|x|=t} dt − \frac{∂u}{∂r}|_{|x|=1}∥_{H^{s}(S^{n−1})}

≤ \frac{1}{ε} ∫_{1−ε}^{1} ∥ \frac{∂u}{∂r}|_{|x|=t} − \frac{∂u}{∂r}|_{|x|=1}∥_{H^{s}(S^{n−1})} dt ⟶_{ε→0} 0.

Der zweite Limes folgt analog mit σ \frac{∂u}{∂r} anstatt u.

3. Die Differentialgleichung für R(t): Es gilt

(\frac{d}{dt_{−}} R(t)) j = \frac{∂u}{∂r}|_{|x|=t} − R(t) (\frac{∂}{∂r}σ \frac{∂u}{∂r})|_{|x|=t}

= \frac{1}{σ} j − R(t) (\frac{∂}{∂r}σ \frac{∂u}{∂r})|_{|x|=t} .  (3.49)

Außerdem ist Δu ∈ H^{s+1/2}({|x| < 1}), so daß die Einschränkung auf |x| = t in H^{s}(S^{n−1}) existiert und

0 = ∇ · σ∇u|_{|x|=t}
= (\frac{∂}{∂r}σ \frac{∂u}{∂r} + \frac{n − 1}{r}σ \frac{∂u}{∂r} + ∇_{τ} · σ∇_{τ}u)|_{|x|=t}

= (\frac{∂}{∂r}σ \frac{∂u}{∂r})|_{|x|=t} + \frac{n − 1}{t}j + ∇_{τ} · σ∇_{τ}R(t)j.  (3.50)

Einsetzen von (3.50) in (3.49) ergibt die Behauptung.

□

Beispiel

Im folgenden beachte man, daß die Funktionen e^{inθ} ∈ C^{∞}(] − π; π[), n ∈ ZZ, ein vollständiges Orthonormalsystem des Raumes L^{2}(] − π; π[) bilden.

Zurückkehrend zum vorherigen Beispiel nimmt die Gleichung (3.48) in Ω_{t} = {x ∈ IR^{2}, |x| < t} folgende Gestalt in Polarkoordinaten (r, θ) an:

\frac{d}{dt_{−}} R(t) = \frac{1}{σ} + \frac{1}{t}R(t) + \frac{1}{t^{2}} R(t)\frac{∂}{∂θ}σ \frac{∂}{∂θ}R(t).

KAPITEL 3. NUMERISCHE VERFAHREN 43 Dabei ist mit \frac{1}{σ} die Multiplikation mit dieser Funktion als Operator gemeint. Ersetzt man R(t) durch tW (t), so ergibt dieses

t\frac{d}{dt_{−}} W (t) = \frac{1}{σ} + W (t)\frac{∂}{∂θ}σ \frac{∂}{∂θ}W (t).  (3.51)

Mit Hilfe von

w_{m,n}(t) := \frac{1}{2π} ∫_{−π}^{π} e^{−imθ} (W (t)e^{in·}) (θ) dθ  (3.52)

kann man W (t) durch

(W (t)f) (θ) = ∑_{m,n=−∞}^{∞} w_{m,n}(t)f_{n}e^{imθ}

ausdrücken, falls f die Fourierreihe

f(θ) = ∑_{n=−∞}^{∞} f_{n}e^{inθ}

besitzt. Es seien σ_{n}(t) und ρ_{n}(t) die Fourierkoeffizienten von σ(t, θ) bzw. \frac{1}{σ(t, θ)}. Um aus (3.51) eine Gleichung für w_{m,n}(t) abzuleiten, betrachte man

(W (t)e^{in·}) (θ) = ∑_{l=−∞}^{∞} w_{l,n}(t)e^{ilθ}
⇒  \frac{∂}{∂θ} (W (t)e^{in·}) (θ) = ∑_{l=−∞}^{∞} w_{l,n}(t)ile^{ilθ}
⇒  σ(t, θ)\frac{∂}{∂θ} (W (t)e^{in·}) (θ) = ∑_{k,l=−∞}^{∞} σ_{k−l}(t)w_{l,n}(t)ile^{ikθ}
⇒  \frac{∂}{∂θ}σ(t, θ)\frac{∂}{∂θ} (W (t)e^{in·}) (θ) = − ∑_{k,l=−∞}^{∞} kσ_{k−l}(t)w_{l,n}(t)le^{ikθ}
⇒ W (t) (\frac{∂}{∂θ}σ(t, θ)\frac{∂}{∂θ} (W (t)e^{in·}) (θ)) = − ∑_{k,l=−∞}^{∞} kσ_{k−l}(t)w_{l,n}(t)l (W (t)e^{ik·}) (θ).

Setzt man dieses in (3.52) ein, dann folgt

t\frac{dw_{m,n}}{dt_{−}} (t) = ρ_{m−n}(t) − ∑_{k,l=−∞}^{∞} w_{m,k}(t)kσ_{k−l}(t)w_{l,n}(t).

Zusammenfassend hat man nun KAPITEL 3. NUMERISCHE VERFAHREN 44

1. w_{n+k,n}(t) = \frac{t}{|n|}ρ_{k}(t) + O(\frac{1}{n^{2}} ), d.h., man kann \frac{1}{σ} aus w_{m,n}(t) und damit aus R(t) berechnen in Ω_{t} = {x ∈ IR^{2}, |x| < t}.

2. w_{m,n}(t − ε) = w_{m,n}(t) − ε\frac{dw_{m,n}}{dt_{−}} (t) + O(ε^{2}), ε > 0, d.h., man kann aus dem Problem in dem Gebiet Ω_{t} eines in Ω_{t−ε} ableiten.

3.5 Ermittlung unstetiger Leitfähigkeiten

3.5.1 Das Problem

Sei Ω ⊂ IR^{N}, N ≥ 2, beschränkt und einfach zusammenhängend mit C^{2}-Rand, und seien D ⊂ Ω offen und K > 0. Definiere σ(x) als

σ(x) = \{ \begin{matrix}1, & x & ∈ & Ω \ D, \\ k, & x & ∈ & D. \end{matrix}

Sei u die Lösung von

L_{σ}u := ∇ · σ∇u = 0 in Ω,  (3.53)
\frac{∂u}{∂ν}|_{∂Ω} = g,
∫_{∂Ω} u dS = 0,
mit ∫_{∂Ω} g dS = 0.

Gesucht ist bei dieser Formulierung des inversen Problems das Gebiet D, wobei k à priori bekannt ist. Zur Vereinfachung werde angenommen, daß eine Abbildung von einer kompakten Menge Q ⊂ IR^{M} in die Menge der Gebiete D existiert: q ↦ D(q). Sei u^{q} die Lösung von (3.53) mit D = D(q) und u die Lösung von (3.53) zum gesuchten Gebiet D.

Definition 3.5.1 Seien N_{i} offene Teilmengen von ∂Ω, und sei F : L^{2}(∂Ω) → IR^{M} definiert durch

F_{i}(ϕ) = ∫_{N_{i}} f_{i}ϕ dS, i = 1, …, M

mit Gewichtsfunktionen f_{i}. Definiere nun die fit-to-data Funktion durch

J(q) = \frac{1}{2} ∑_{i=1}^{M} |F_{i}(u^{q}) − F_{i}(u)|^{2}.  (3.54)

3.5.2 Eine Approximation von u

Im Einheitskreis ist

G(x, y) = −\frac{1}{2π} (log(|x − y|) + log(|\frac{x}{|x|^{2}} − y|))

eine Green’sche Funktion für −Δu = f, d.h., G(x, y) ist Lösung von

Δ_{y}G(x, y) = −δ(x) für y ∈ Ω, ∀x ∈ Ω,
\frac{∂}{∂ν_{y}} G(x, y)|_{y∈∂Ω} = −\frac{1}{|∂Ω|},
∫_{∂Ω} G(x, y) dS_{y} = 0.

Um Lösungen von

Δu_{0} = f in Ω,
\frac{∂u_{o}}{∂ν} = g auf ∂Ω,
mit ∫_{∂Ω} u_{0} dS = 0

zu ermitteln, kann man nun G(x, y) verwenden:

u_{0}(x) = ∫_{∂Ω} G(x, y)g(y) dS_{y} − ∫_{Ω} G(x, y)f(y) dy, ∀x ∈ Ω.  (3.55)

Zu gegebenem D sei u(k, x) die Lösung von (3.53). Definiere ψ als

ψ(x) = ∫_{∂D} u_{0}(y)\frac{∂}{∂ν_{y}} G(x, y) dS_{y}, ∀x ∈ ∂Ω.  (3.56)

Damit kann man u(k, x) folgendermaßen approximieren:

Satz 3.5.2 Für festes x ∈ ∂Ω ist die Abbildung k ↦ (u(k, x)) in k = 1 differenzierbar bzgl. k, und die Ableitung ist durch (3.56) gegeben. Es gilt sogar

|u(k, x) − u_{0}(x) − (k − 1)ψ(x)| ≤ C(k − 1)^{2},

für eine Konstante C, die lediglich von D und u_{0} abhängt.

Beweis: Aufgrund der Eigenschaften von G(x, y) und des Gauß’schen Integralsatzes kann man für die Lösung u(k, x) von (3.53) folgende zwei Integralgleichungen herleiten:

Für x ∈ Ω̄ \ D gilt

u(k, x) = u_{0}(x) + (1 − k) ∫_{∂D} u(k, y)\frac{∂}{∂ν_{y}} G(x, y) dS_{y}.  (3.57)

Für x ∈ ∂D genügt u(k, x) der Gleichung

u(k, x) = \frac{2}{1 + k} u_{0}(x) + \frac{2(1 − k)}{1 + k} ∫_{∂D} u(k, y)\frac{∂}{∂ν_{y}} G(x, y) dS_{y}.  (3.58)

Die Behauptung beweist man nun, indem man (3.58) in der Supremumsnorm abschätzt und das Ergebnis dann in (3.57) einsetzt. Den genauen Beweis kann man in Bryan [B3] nachlesen. □

Eine lineare Approximation ũ(x) der Lösung u(k, x) von (3.53) ist also durch

ũ(x) = u_{0}(x) + (k − 1)ψ(x).

gegeben.

3.5.3 Numerische Experimente

Nun kann man das originale inverse Problem umformulieren, indem man u durch ũ aus Definition 3.5.2 ersetzt. Für numerische Berechnungen müssen noch Q und D(q) definiert werden. Als einfachstes Beispiel für Ω = {x ∈ IR^{2}, |x| < 1} nennt Bryan in [B3]

Q := \{(a, b, r) ∈ IR^{3}, \sqrt{a^{2} + b^{2}} + r ≤ ρ < 1\} und
D(q) := {x ∈ Ω, (x_{1} − a)^{2} + (x_{2} − b)^{2} ≤ r^{2}} .

Damit führt Bryan zunächst eine erste Approximation durch. Um aber allgemeinere Gebiete zu ermitteln, parametrisiert er ∂D(q) folgendermaßen:

x(t) = a + r(t) cos(t),
y(t) = b + r(t) sin(t),

für 0 ≤ t ≤ 2π und a und b aus der ersten Approximation. r(t) sei ein kubischer Spline mit r(0) = r(2π), r′(0) = r′(2π) und r(t_{i}) = q_{i}, i = 1 … m mit t_{i} = 2πi/m. Nun wendet man den Algorithmus auf q_{i} … q_{m} an.

Ein spezielles iteratives Verfahren

4.1 Das Iterationsverfahren

4.1.1 Der nichtlineare Fall

Sei R_{j}(f) eine Familie nichtlinearer Fréchet-differenzierbarer Operatoren, und gesucht sei eine Lösung des Systems

R_{j}(f) = 0, für j = 1, …, p.  (4.1)

Da R_{j}(f) Fréchet-differenzierbar ist, gilt

R_{j}(f + h) − R_{j}(f) = R′_{j}(f)h + o(∥h∥^{2}).  (4.2)

Nun soll f + h eine Approximation der Lösung f^{∗} von (4.1) sein. Man setzt deshalb R_{j}(f + h) = 0. Vernachlässigt man den o(∥h∥^{2}) Term und löst nach h auf, so erhält man als kleinste-Quadrate-Lösung

h = −R′_{j}(f)^{∗} (R′_{j}(f)R′_{j}(f)^{∗})^{−1} R_{j}(f).  (4.3)

Ersetzt man R′_{j}(f)R′_{j}(f)^{∗} zur weiteren Vereinfachung durch eine Konstante c_{j}, dann führt dieses zu folgendem Verfahren:

f_{0} = f^{0},
f_{j} = f_{j−1} − \frac{1}{c_{j}} R′_{j}(f)^{∗}R_{j}(f), j = 1, …, p,  (4.4)
f^{1} = f_{p}.

KAPITEL 4. EIN SPEZIELLES ITERATIVES VERFAHREN 49

4.1.2 Der lineare Fall

Satz 4.1.1 Für ein konsistentes System linearer Gleichungen

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

konvergiert das Iterationsverfahren

f_{0} = f^{0},
f_{j} = f_{j−1} − \frac{1}{c_{j}} R_{j}^{∗}R_{j}f_{j−1},  (4.6)
f^{1} = f_{p}

für c_{j} > \frac{1}{2}ρ (R_{j}R_{j}^{∗}). Gilt zusätzlich

f^{0} ∈ ∑_{j=1}^{p} Bild (R_{j}^{∗}) ,  (4.7)

dann konvergiert (4.6) gegen die Lösung von (4.5) mit minimaler Norm. Bemerkung 4.1.2 Dieses Iterationsverfahren ist als algebraic reconstruction technique (ART) oder Kaczmarzmethode bekannt. Beweis: Sei

f_{j} = f^{0} + ∑_{k=1}^{j} R_{k}^{∗}u_{k}.

Dann gilt:

f^{0} + ∑_{k=1}^{j} R_{k}^{∗}u_{k} = (I − \frac{1}{c_{j}} R_{j}^{∗}R_{j}) (f^{0} + ∑_{k=1}^{j−1} R_{k}^{∗}u_{k})
⇒  R_{j}^{∗}u_{j} = −\frac{1}{c_{j}} R_{j}^{∗}R_{j}f^{0} − \frac{1}{c_{j}} R_{j}^{∗} ∑_{k=1}^{j−1} R_{j}R_{k}^{∗}u_{k}
⇒  c_{j}u_{j} = −R_{j}f^{0} − ∑_{k=1}^{j−1} R_{j}R_{k}^{∗}u_{k}.

Setze nun

R = ( \begin{matrix}R_{1} & & & & 0 \\ & . & . & . \\ 0 & & & & R_{p} \end{matrix} ) ,  D = ( \begin{matrix}R_{1}R_{1}^{∗} & & & & 0 \\ & . & . & . \\ 0 & & & & R_{p}R_{p}^{∗} \end{matrix} ) ,
C = ( \begin{matrix}c_{1}I & & & & 0 \\ & . & . & . \\ 0 & & & & c_{p}I \end{matrix} ) ,
A = RR^{∗} = L + L^{∗} + D.

Daraus folgt

Cu = −Rf^{0} − Lu
⇒  u = −(C + L)^{−1}Rf^{0}.

Man nehme jetzt f^{k} = ∑_{j=1}^{p} R_{j}^{∗}u_{j}^{k} = R^{∗}u^{k} und damit

R^{∗}u^{1} = f^{1}
= f^{0} + R^{∗}u = R^{∗}u^{0} + R^{∗}u
= R^{∗}u^{0} − R^{∗}(C + L)^{−1}RR^{∗}u^{0}
⇒  u^{1} = u^{0} − (C + L)^{−1}RR^{∗}u^{0}
= (C + L)^{−1}(C + L − L − L^{∗} − D)u^{0}
= (C + L)^{−1}(C − L^{∗} − D)u^{0}.

Betrachte nun das Eigenpaar (λ, v) von B = (C + L)^{−1}(C − L^{∗} − D):

(C + L)^{−1}(C − L^{∗} − D)v = λv
⇒  λ = \frac{⟨(C − L^{∗} − D)v, v⟩}{⟨(C + L)v, v⟩}.

Seien nun o.E. ⟨Dv, v⟩ = 1, ⟨Cv, v⟩ = γ > \frac{1}{2} und ⟨Lv, v⟩ = α + iβ. Da 0 ≤ ⟨RR^{∗}v, v⟩ = ⟨(L + L^{∗} + D)v, v⟩, folgt α ≥ −\frac{1}{2} und deshalb

γ − α − 1 ≤ γ + \frac{1}{2} − 1
≤ γ + α
und α + 1 − γ < α + 1 − \frac{1}{2}
< α + γ
⇒ |γ − α − 1| ≤ |γ + α| = γ + α.

Daraus folgt nun

|λ|^{2} = |\frac{γ − α − 1 + iβ}{γ + α + iβ} |^{2}
= \frac{(γ − α − 1)^{2} + β^{2}}{(γ + α)^{2} + β^{2}} ≤ 1,

und für ⟨RR^{∗}v, v⟩ ≠ 0 gilt sogar

α > −\frac{1}{2} und deshalb |λ| < 1.  (4.8)

Also ist der Spektralradius ρ(B) < 1, und das Iterationsverfahren konvergiert nach dem Banachschem Fixpunktsatz. □

4.1.3 Die optimale Wahl von c_{j}

Zurückkehrend zum nichtlinearen Fall sei s = \frac{1}{c_{j}}. Dann gilt

∥f^{∗} − f_{j}∥^{2} = ∥f^{∗} − f_{j−1} + sR_{j}^{∗}R_{j}f_{j−1}∥^{2}
= ∥f^{∗} − f_{j−1}∥^{2} + s^{2}∥R_{j}^{∗}R_{j}f_{j−1}∥^{2} + 2s⟨f^{∗} − f_{j}, R_{j}^{∗}R_{j}f_{j−1}⟩
= ∥f^{∗} − f_{j−1}∥^{2} + s^{2}∥R_{j}^{∗}R_{j}f_{j−1}∥^{2} + 2s⟨R_{j}(f^{∗} − f_{j−1}), R_{j}f_{j−1}⟩
= ∥f^{∗} − f_{j−1}∥^{2} + s^{2}∥R_{j}^{∗}R_{j}f_{j−1}∥^{2} − 2s⟨R_{j}f_{j−1}, R_{j}f_{j−1}⟩,

da R_{j}f^{∗} = 0. Differenziert man nun diese Gleichung nach s, um das Minimum zu finden, und setzt die Ableitung gleich Null, so erhält man als Bedingung

0 = \frac{d}{ds}∥f^{∗} − f_{j}∥^{2} = 2s∥R_{j}^{∗}R_{j}f_{j−1}∥^{2} − 2∥R_{j}f_{j−1}∥^{2}.

Die optimal Wahl von c_{j} erfüllt also

\frac{1}{c_{j}} = \frac{∥R_{j}f_{j−1}∥^{2}}{∥R_{j}^{∗}R_{j}f_{j−1}∥^{2}}  (4.9)

oder

c_{j} = \frac{∥R_{j}^{∗}R_{j}f_{j−1}∥^{2}}{∥R_{j}f_{j−1}∥^{2}}
= \frac{⟨R_{j}R_{j}^{∗}v, v⟩}{∥v∥^{2}} mit v = R_{j}f_{j−1}.  (4.10)

Da ρ(R_{j}R_{j}^{∗}) = sup_{v≠0} \frac{⟨R_{j}R_{j}^{∗}v, v⟩}{∥v∥^{2}}, ist also zu hoffen, daß v “geeignet” ist, so daß

c_{j} > \frac{1}{2}ρ(R_{j}R_{j}^{∗}), und damit die Konvergenz im linearen Fall gewährleistet ist.

4.2 Anwendung auf die EIT

4.2.1 Die Fréchet-Ableitung

Seien Ω ⊂ IR^{n}, n ≥ 2, ein beschränktes, hinreichend glattes Gebiet, σ ∈ L^{∞}(Ω) mit positiver unterer Schranke (σ(x) ≥ σ_{0} > 0 ∀ x ∈ Ω) und f_{j} ∈ H^{−1/2}(∂Ω). Dann besitzt die Randwertaufgabe

\begin{aligned}∇ · σ∇u_{j} & = & 0 & in & Ω, \\ σ \frac{∂u_{j}}{∂ν} & = & f_{j} & auf & ∂Ω \end{aligned}  (4.11)

eine eindeutige (schwache) Lösung u_{j} ∈ H^{1}(Ω).

Da Ω ein beschränktes Gebiet ist, also für das Lebesguemaß λ gilt: λ(Ω) < ∞, ist L^{2}(Ω) ⊂ L^{∞}(Ω) wegen ∥f∥_{L^{2}} ≤ \sqrt{λ(Ω)}∥f∥_{L^{∞}}. Seien nun g_{j} ∈ H^{1/2}(∂Ω) Meßwerte, die zu der angelegten Stromdichte f_{j} gehören. Für allgemeine σ (ohne positive untere Schranke) definiert man den Operator R_{j} : L^{2}(Ω) → H^{1/2}(∂Ω) durch

R_{j}(σ) = u_{j}|_{∂Ω} − g_{j}.  (4.12)

Da u_{j} ∈ H^{1}(Ω) ist, gilt u_{j}|_{∂Ω} ∈ H^{1/2}(Ω) (siehe Anhang B). Für allgemeine σ in L^{2}(Ω) geht natürlich die Eindeutigkeit in (4.11) verloren, aber darauf komme ich im nächsten Abschnitt zurück.

Um die Fréchet-Ableitung zu bilden, linearisiert man (4.11) und (4.12), indem man zunächst u_{j} durch u_{j} + w_{j} und σ durch σ + η ersetzt:

0 = ∇((σ + η)∇(u_{j} + w_{j}))
= ∇(σ∇u_{j}) + ∇(η∇u_{j}) + ∇(σ∇w_{j}) + ∇(η∇w_{j})
in Ω.
f_{j} = (σ + η)\frac{∂(u_{j} + w_{j})}{∂ν}
= σ \frac{∂u_{j}}{∂ν} + η \frac{∂u_{j}}{∂ν} + σ \frac{∂w_{j}}{∂ν} + η \frac{∂w_{j}}{∂ν}
auf ∂Ω.
R_{j}(σ + η) = (u_{j} + w_{j})|_{∂Ω} − g_{j}
= R_{j}(σ) + w_{j}|_{∂Ω}.

Benutzt man, daß u_{j} (4.11) löst, und vernachlässigt man die Terme, in denen sowohl w_{j}, als auch η auftauchen (dieses ist möglich, da sie von zweiter Ordnung sind), so erhält man die Fréchet-Ableitung

∇(σ∇w_{j}) = −∇(η∇u_{j}) in Ω,  (4.13)
σ \frac{∂w_{j}}{∂ν} = −η \frac{∂u_{j}}{∂ν} auf ∂Ω,
R′_{j}(σ)η = w_{j}|_{∂Ω},  (4.14)

und für festes σ ist die Richtungsableitung R′_{j}(σ) ein Operator von L^{2}(Ω) nach H^{1/2}(∂Ω). Da u_{j} bekannt ist und die Richtung η vorgegeben ist, sind die rechten Seiten in (4.13) bekannt.

4.2.2 Der adjungierte Operator

Satz 4.2.1 Der zu R′_{j}(σ) adjungierte Operator R′_{j}(σ)^{∗} : H^{−1/2}(∂Ω) → L^{2}(Ω) ist gegeben durch

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

wobei z ∈ H^{1}(Ω) die (eindeutige) Lösung von

∇(σ∇z) = 0 in  Ω
σ \frac{∂z}{∂ν} = −g auf ∂Ω, g ∈ H^{−1/2}(∂Ω), (4.16)

und u_{j} ∈ H^{1}(Ω) die Lösung von (4.11) ist.

Beweis: ∀z, w ∈ C^{2}(Ω) gilt

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

Durch Einsetzen von (4.14) und (4.16) erhält man

− ∫_{Ω} ∇(η∇u_{j})z dx = ∫_{∂Ω} wg − η \frac{∂u_{j}}{∂ν} z dS

= ∫_{∂Ω} (R′_{j}(σ)η) g − η \frac{∂u_{j}}{∂ν} z dS.

Nach dem Gauß’schen Integralsatz gilt:

∫_{Ω} ∇(η∇u_{j})z + η∇u_{j} · ∇z dx = ∫_{Ω} ∇ · (zη∇u_{j}) dx

= ∫_{∂Ω} η \frac{∂u_{j}}{∂ν} z dS

und deshalb

∫_{Ω} η∇u_{j} · ∇z dx = ∫_{∂Ω} (R′_{j}(σ)η) g dS.

Da C^{2}(Ω) dicht in H^{1}(Ω) liegt und H^{1/2}(∂Ω) bzgl. ⟨·, ·⟩_{L^{2}(∂Ω)} der Dualraum von H^{−1/2}(∂Ω) ist, gilt:

⟨η, ∇u_{j} · ∇z⟩_{L^{2}(Ω)} = ⟨R′_{j}(σ)η, g⟩_{L^{2}(∂Ω)},

also

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

4.3 Zusammenfassung

In der Impedanz-Tomographie hat der i-te Schritt des Iterationsverfahrens folgende Gestalt:

• σ_{0} = σ^{i}.

• Für j = 1, …, p und σ_{j−1} fest löse

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

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

h_{j} = \frac{1}{c_{j}} ∇u_{j} · ∇z_{j},

c_{j} = \frac{∥h_{j}∥_{L^{2}(Ω)}^{2}}{∥u_{j} − g_{j}∥_{L^{2}(∂Ω)}^{2}} ,

σ_{j} = σ_{j−1} + ωh_{j},

für ein “geeignetes” ω > 0 (z.B. ω = 1).

• σ^{i+1} = σ_{p}.

Wie im vorherigen Abschnitt schon erwähnt, benötigt man für Eindeutigkeit in (4.11), daß σ_{j} ∈ L^{2}(Ω) eine positive untere Schranke besitzt. In der Anwendung kann man sicherlich annehmen, daß für die gesuchte Leitfähigkeit σ^{∗} eine positive untere Schranke σ_{0}^{∗} bekannt ist. Den Startwert σ^{0} wählt man so, daß er dieses auch erfüllt.

Im linearen Fall ist die einzige Bedingung an ω, daß

\frac{c_{j}}{ω} > \frac{1}{2}ρ(R_{j}R_{j}^{∗})

gilt. Deshalb kann man ω in jedem Schritt so klein wählen, daß z.B. immer noch

σ_{j}(x) > \frac{σ_{0}^{∗}}{2} , ∀ x ∈ Ω ∀ j = 1, …, p,

erfüllt wird. Also ist die Positivität der σ_{j} in jedem einzelnen Schritt garantiert. Damit sind (4.11) und (4.16) eindeutig (schwach) lösbar, da die partielle Differentialgleichung zweiter Ordnung streng elliptisch ist.

Numerische Experimente

Die Fähigkeiten des Algorithmus will ich anhand einiger Ergebnisse erläutern, die ich mit einem kleinen, selbstgeschriebenen Programm erzielt habe. Für einen Schritt des Algorithmus muß zweimal das von-Neumann-Problem

∇ · σ∇u = 0 in Ω,
σ \frac{∂u}{∂ν} = f auf ∂Ω  (5.1)

mit verschiedenen Randdaten f für festes σ gelöst werden. Dieses Problem löse ich numerisch mit finiten Differenzen und einem SOR-Iterationsverfahren. Als Grundgebiet sind sowohl Ω = (0, 1)^{2} in euklidischen Koordinaten, als auch Ω = {x ∈ IR^{2} : |x| < 1} in Polarkoordinaten implementiert.

5.1 Approximation der Differentialgleichung

Für u ∈ C^{4}(Ω) werden Ableitungen erster und zweiter Ordnung durch folgende Approximationen ersetzt:

\frac{u(x + h, y) − u(x − h, y)}{2h} = \frac{∂u}{∂x}(x, y) + h^{2} \frac{∂^{3}u}{∂x^{3}}(ξ, y), (5.2)
\frac{u(x + h, y) − 2u(x, y) + u(x − h, y)}{h^{2}} = \frac{∂^{2}u}{∂x^{2}}(x, y) + h^{2} \frac{∂^{4}u}{∂x^{4}}(ξ̃, y). (5.3)

5.1.1 Der Fall Ω = (0, 1)^{2}

Die Differentialgleichung läßt sich auch in folgender Form schreiben:

∇σ · ∇u + σΔu = 0.  (5.4)

Legt man über Ω ein äquidistantes Gitter mit Schrittweite h = 1/M und ersetzt nun die Ableitungen erster und zweiter Ordung durch (5.2) und (5.3), so erhält man:

16σ_{i,j}u_{i,j} = (σ_{i+1,j} − σ_{i−1,j} + 4σ_{i,j})u_{i+1,j}
+(σ_{i−1,j} − σ_{i+1,j} + 4σ_{i,j})u_{i−1,j}
+(σ_{i,j+1} − σ_{i,j−1} + 4σ_{i,j})u_{i,j+1}
+(σ_{i,j−1} − σ_{i,j+1} + 4σ_{i,j})u_{i,j−1}  (5.5)

mit einem lokalen Fehler O(h^{2}). Für die Normalableitung auf ∂Ω gilt folgende Beziehung:

σ \frac{u(x + hν) − u(x − hν)}{2h} = σ \frac{∂u}{∂ν}(x) + O(h^{2}) = f(x) + O(h^{2}).  (5.6)

Dadurch lassen sich die “Randterme” u_{−1,j}, u_{M+1,j}, u_{i,−1} und u_{i,M+1} bestimmen:

u_{−1,j} = \frac{2h}{σ_{0,j}} f(0, y_{j}) + u_{1,j}.  (5.7)

Die “Randterme” σ_{−1,j}, σ_{M+1,j}, σ_{i,−1} und σ_{i,M+1} werden durch lineare Interpolation ermittelt. Dieses ist gleichbedeutend mit Verwendung von Vorwärts- bzw. Rückwärtsdifferenzen anstatt zentrierter Differenzen:

σ_{M+1,j} := 2σ_{M,j} − σ_{M−1,j}
⟹ σ_{M+1,j} − σ_{M−1,j} = 2σ_{M,j} − σ_{M−1,j} − σ_{M−1,j}
= 2(σ_{M,j} − σ_{M−1,j}).

Dadurch ist allerdings der lokale Fehler am Rand O(h). Im Programm kann man die Interpolation leichter implementieren, indem man den Test, ob ein Parameterpaar auf dem Rand liegt, in einer externen Prozedur durchführt und sich so in der eigentlichen Prozedur nicht mehr um Randterme kümmern muß.

Im SOR-Verfahren habe ich ω = \frac{2}{1 + sin(πh)} gesetzt. Dieser Wert ist für das

Dirichlet Problem und σ ≡ 1 optimal. Es gilt nämlich:

u_{i,j}^{(k,l)} = sin(πhik) sin(πhjl), 1 ≤ k, l < M,

ist Eigenvektor von

(Au)_{i,j} = u_{i+1,j} + u_{i−1,j} + u_{i,j+1} + u_{i,j−1} − 4u_{i,j}

zum Eigenwert 1 − sin^{2}(πhk) − sin^{2}(πhl). Da die u^{(k,l)} linear unabhängig sind, bilden sie eine Basis des IR^{2(M−1)}. Daraus folgt, daß ρ(A) = cos(πh). Die optimale Wahl von ω für dieses SOR-Verfahren lautet:

ω = \frac{2}{1 + \sqrt{1 − ρ(A)^{2}}} = \frac{2}{1 + sin(πh)}.

Beweis: Siehe Skript zu Cryer, Numerik Partieller Differentialgleichungen II. □

Da die tatsächlich verwendete Matrix A(σ) diesem Problem sehr ähnlich ist, hoffe ich, damit eine bessere Konvergenzgeschwindigkeit zu erzielen.

5.1.2 Der Fall Ω = {x ∈ IR^{2} : |x| < 1}

Für Ω = {x ∈ IR^{2} : |x| < 1} führe ich zunächst Polarkoordinaten ein:

x(r, ϕ) = r cos(ϕ),
y(r, ϕ) = r sin(ϕ).

In (r, ϕ)-Koordinaten lautet (5.1):

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

Für jede Funktion v in Polarkoordinaten muß v(0, ϕ) = v(0, 0), also konstant sein. Deshalb folgt für r → 0 nach de l’Hôspital

\frac{∂σ}{∂r} (0, ϕ)\frac{∂u}{∂r}(0, ϕ) + 2σ(0, ϕ)\frac{∂^{2}u}{∂r^{2}}(0, ϕ) = 0.  (5.9)

Die Ableitungen und Randterme werden wie oben approximiert. Hier habe ich für

das SOR-Verfahren ω = \frac{2}{1 + sin(k)} gewählt, dabei entspricht k der Gitterweite

in Winkelrichtung. Dieses begründe ich wie oben, indem ich die Eigenfunktionen und Eigenwerte in Radiusrichtung vernachlässige.

5.1.3 Randdaten

Zunächst fällt auf, daß die Lösung u von (5.1) linear von den Randdaten f abhängt. Im weiteren Verlauf dieses Abschnitts werde ich im Gegensatz zum vorherigen Kapitel die Randwerte mit zwei Indizes (k, l) numerieren. Mit E bezeichne ich die Anzahl Elektroden auf ∂Ω. Seien (x_{k}, y_{k}) und (x_{l}, y_{l}) Randgitterpunkte, fortlaufend numeriert, dann ist f_{k,l} für k ≠ l gegeben durch

f_{k,l}(x, y) = \{ \begin{matrix}1 & für & (x, y) = (x_{k}, y_{k}) & ∈ & ∂Ω, \\ −1 & für & (x, y) = (x_{l}, y_{l}) & ∈ & ∂Ω, \\ 0 & sonst. \end{matrix}

Hier genügt es aus Symmetriegründen, daß l ≤ E/2 ist (E sei ohne Einschränkung gerade). Um effizient eine (exakte) Messung von Randdaten g_{k,l} := u_{k,l}|_{∂Ω} zu simulieren, nutze ich die lineare Abhängigkeit der Lösungen von den Randdaten aus. Eine Basis des Raumes der f_{k,l} sind die Funktionen f_{k,k+1}. Durch Addition kann man alle anderen Randdaten errechnen. Rauschen wird natürlich erst nach dieser Addition zu den g_{k,l} hinzugefügt.

5.1.4 Meßfehler

Um Meßfehler zu simulieren, werden die Meßwerte g_{k,l} mit dem Faktor

1 + p · r_{k,l}(x)

multipliziert. Dabei sind die r_{k,l}(x) ∈ [−1, 1] Zufallszahlen (für jedes g_{k,l} für jeden Gitterpunkt x unabhängig voneinander bestimmt), und p ∈ [0, 1] gibt an, um wieviel die so erzeugten “fehlerhaften” Daten maximal von den korrekten Werten abweichen. Im weiteren Verlauf werde ich, falls p > 0 ist, p in Prozent Rauschen angeben.

5.1.5 Numerische Durchführung einer Berechnung

Die beiden von-Neumann-Probleme eines einzelnen Schrittes werden, wie oben beschrieben, berechnet. Das Skalarprodukt der Gradienten wird mit zentrierten Differenzen approximiert. Dabei nimmt man an, daß σ auf dem Rand von Ω bekannt ist. So tauchen bei den zentrierten Differenzen auch keine undefinierten Randwerte auf.

5.2 Ergebnisse

5.2.1 Erste Ergebnisse

Zuerst beschäftige ich mich mit einem einfachen Problem (siehe Abb. 5.1). Hieran kann man gut erkennen, wie der Algorithmus arbeitet, und was er überhaupt leistet. Wie bei allen weiteren Experimenten hat die Leitfähigkeit σ normalerweise den Wert Eins. Die Aufgabe besteht darin, eine Kreisscheibe zu indentifizieren, in der σ den Wert 2 hat. Abbildung 5.1 wurde in Polarkoordinaten mit 16 Gitterpunkten in Radiusrichtung und 64 Gitterpunkten in Winkelrichtung erstellt. Das Verhältnis von viermal so vielen Gitterpunkten in Winkelrichtung erwies sich bei ersten Versuchen als besonders günstig, außerdem kann man so die Ergebnisse besser mit den Ergebnissen in orthographischen Koordinaten vergleichen:

Abbildung 5.1: Kreis der Höhe 2 in einem Gebiet der Höhe 1

Bei gleicher Gesamtanzahl an Gitterpunkten bleibt so auch die Anzahl Gitterpunkte in horizontaler bzw. vertikaler Richtung über die Mitte hinweg konstant. Es wurde immer ω = 3 gesetzt, was sich als der beste Kompromiß zwischen Konvergenzgeschwindigkeit und Stabilität erwies.

Man sieht deutlich, wie in der Approximation die Inhomogenität verschwimmt, besonders zur Mitte hin. Außerdem bildet sich am Rand ein kleiner “Graben”, und die Spitze wandert leicht nach außen. Ferner zeigte sich, daß der Algorithmus nur so lange schnell konvergiert, bis die Approximation ungefähr die gleiche l_{2}- Norm wie das Original besitzt. Deshalb ist auch nicht verwunderlich, daß die Höhe der rekonstruierten Inhomogenität nur ungefähr 1,6 beträgt, da sie sich über ein viel größeres Gebiet verteilt. (l_{2}-Norm des Originals: 46.8188; l_{2}-Norm der Rekonstruktion: 42.0642) Dieser Effekt wird sich im weiteren besonders bei der Identifizierung von “kleinen” Objekten zeigen. Ein ähnliches Ergebnis erhält man in karthesischen Koordinaten im Quadrat mit einer rechteckigen Inhomogenität, in der σ = 2 gilt, mit einem 32x32 Gitter (Abb. 5.2).

Abbildung 5.2: Rechteck der Höhe 2 in einem Gebiet der Höhe 1

5.2.2 Rauschen

Hier habe ich die gleichen Berechnungen wie oben durchgeführt, mit dem einzigen Unterschied, daß die simulierten Meßdaten mit 10 Prozent Rauschen versehen wurden (vergleiche Abschnitt 5.1.4). Das Ergebnis wird noch mehr verwischt, und die Artefakte am Rand nehmen deutlich zu. Trotzdem lassen sich noch qualitative Schlußfolgerungen aus der Grafik ziehen (siehe Abb. 5.3).

KAPITEL 5. NUMERISCHE EXPERIMENTE 61 Abbildung 5.3: Wie oben mit 10 Prozent Rauschen

5.2.3 Trennung mehrerer Objekte

Abbildung 5.4: Zwei Kreise der Höhe 2, gleiches Gitter

Der Algorithmus trennt größere, gleichhohe Objekte mit einem relativ großen Abstand recht gut (Abb. 5.4 und Abb. 5.6). In Abbildung 5.5 sieht man den Einfluß von 10 Prozent Rauschen in den Anfangsdaten auf die Rekonstruktion. Problematisch werden erst Objekte verschiedener Größe. Abbildung 5.7 zeigt, wie ein niedriges Objekt (Höhe 1,5) im “Schatten” eines hohen Objekts (Höhe 3,0) verschwindet. Das Objekt der Höhe 2,5 wird recht gut rekonstruiert, während das Objekt der Höhe 2 nur noch gerade eben getrennt werden kann.

KAPITEL 5. NUMERISCHE EXPERIMENTE 63

Abbildung 5.5: Das gleiche mit 10 Prozent Rauschen Abbildung 5.6: 3 Objekte der Höhe 2 KAPITEL 5. NUMERISCHE EXPERIMENTE 64 Abbildung 5.7: Vier Objekte mit Höhen 1,5 bis 3,0 auf einem Gebiet mit Höhe 1

5.2.4 Kleine Objekte

Abbildung 5.8: Ein kleiner Kreis der Höhe 2

Das Problem bei kleinen Inhomogenitäten ist, daß der Algorithmus zwar das Objekt findet, aber die ermittelte Höhe ist so gering, daß das Ergebnis im Rauschen untergeht. Das Ergebnis in Abbildung 5.8 stimmt mit der Beobachtung überein, daß der Algorithmus nur solange konvergiert, bis die Rekonstruktion ungefähr die gleiche l_{2}-Norm wie das Original hat (l_{2}-Norm des Originals: 34.322; l_{2}-Norm der Rekonstruktion: 33.783). Dieses ist auch nicht weiter verwunderlich, da der Algorithmus die kleinste-Quadrate-Lösung von (4.12) approximieren soll. Allerdings verschwinden so kleine Inhomogenitäten im Rauschen, wie Abbildung 5.9 zeigt.

KAPITEL 5. NUMERISCHE EXPERIMENTE 66

Abbildung 5.9: Ein kleiner Kreis der Höhe 2 bei 5 Prozent Rauschen Abbildung 5.10: Zwei kleine Kreise der Höhe 2

5.2.5 Objekte geringer Höhe

Abbildung 5.11: Ein Kreis der Höhe 1.1

Bei Objekten geringer Höhe taucht prinzipiell das gleiche Problem wie bei kleinen Objekten auf (In Abbildung 5.11: l_{2}-Norm des Originals: 33.6248; l_{2}-Norm der Rekonstruktion: 33.6743).

KAPITEL 5. NUMERISCHE EXPERIMENTE 68 Abbildung 5.12: Ein Kreis der Höhe 1.1 bei 5 Prozent Rauschen

5.2.6 Geschwindigkeit

Abbildung 5.13: Wie Abb. 5.1, aber nur ein Durchlauf

Abbildung 5.13 zeigt die Rekonstruktion des gleichen Objekts wie in Abbildung 5.1, aber mit nur einem statt 20 Durchläufen. (l_{2}-Norm des Originals: 46.8188; l_{2}-Norm der Rekonstruktion: 43.4529) Ein Durchlauf benötigt bei 250 Iterationen pro von-Neumann-Problem ungefähr 46 Sekunden Systemzeit auf einer SUN SPARCstation 20 (arminius). Abb. 5.14 zeigt das gleiche Problem, rekonstruiert mit 100 Durchläufen. Es gibt zwar noch eine quantitative, aber keine qualitative Verbesserung.

KAPITEL 5. NUMERISCHE EXPERIMENTE 70 Abbildung 5.14: Wie Abb. 5.1, aber 100 Durchläufe

Anhang A

Das C-Programm

Die X-Oberfläche

Abbildung A.1: Screendump

ANHANG A. DAS C-PROGRAMM 72

Zum Ausprobieren und interaktiven Verändern der vielen Parameter habe ich eine kleine X-Oberfläche für mein Programm geschrieben. So kann man dem Algorithmus direkt bei der Arbeit zuschauen. Dies war vor allem bei der Erkennung und Beseitigung von Fehlern nützlich. Außerdem konnte ich nur so sehen, welche Werte von ω noch sinnvoll sind.

Teile des Quellcodes

Auszüge aus der Datei citlib.c.

extern float rauschen, Mult;
extern void checkAndProcess();
extern int RAND, CJSW, TRY, POLAR, Anzahl, dAnz, M, N, E;
extern float u[E_MAX][E2_MAX][I][J], f[E_MAX][E2_MAX][JJ];
extern float v[E_MAX][E2_MAX][I][J], g[E_MAX][E2_MAX][JJ];
extern float s_org[I][J], sigma[I][J], n_u_g, Cj, min_sigma, max_sigma;
extern char coordChoiceText[2][6];
inline float rauschFaktor() {
return 1.+rauschen*(2.*drand48()-1.);
}
inline int phi( int j, int M ) {
if ((j>=0)&&(j<M))
      return j;
if (j<0)
      return M-1;
else
      return 0;
}

void lift_polar( float u[I][J], int M, int N ) {
int i, j;
float h=0;
for (j=0; j<N; j++)
  h+=u[M][j];
h/=N;
for (i=0; i<=M; i++) for (j=0; j<N; j++)
  u[i][j]-=h;
}
void neumann_polar( float u[I][J], float f[JJ], float sigma[I][J], int M, int N, int anz ){
int i, j, k, jj=0, NN=N/2;
float tmp=N*N/(4*M_PI*M_PI), tmp2, delta=eps+1, max_u_r;
float srp[I][J], srm[I][J], sphip[I][J], sphim[I][J];
float w=2./(1.+sin(M_PI/N));
for (j=0; j<N; j++) {
      tmp2=8.0*(M*M+tmp)*sigma[M][j];
      srp[M][j]=(M*M*2*(sigma[M][j]-sigma[M-1][j])+(4.*M*M+2.*M)*sigma[M][j])
        /tmp2;
      srm[M][j]=(M*M*2*(sigma[M-1][j]-sigma[M][j])+(4.*M*M-2.*M)*sigma[M][j])
        /tmp2;
      sphip[M][j]=tmp*(sigma[M][phi(j+1,N)]-sigma[M][phi(j-1,N)]+4.*sigma[M][j])
        /tmp2;
      sphim[M][j]=tmp*(sigma[M][phi(j-1,N)]-sigma[M][phi(j+1,N)]+4.*sigma[M][j])
        /tmp2;
}
for (i=1; i<M; i++) for (j=0; j<N; j++) {

  tmp2=8.0*(i*i+tmp)*sigma[i][j];
  srp[i][j]=(i*i*(sigma[i+1][j]-sigma[i-1][j])+(4.*i*i+2.*i)*sigma[i][j])
    /tmp2;
  srm[i][j]=(i*i*(sigma[i-1][j]-sigma[i+1][j])+(4.*i*i-2.*i)*sigma[i][j])
    /tmp2;
  sphip[i][j]=tmp*(sigma[i][phi(j+1,N)]-sigma[i][phi(j-1,N)]+4.*sigma[i][j])
    /tmp2;
  sphim[i][j]=tmp*(sigma[i][phi(j-1,N)]-sigma[i][phi(j+1,N)]+4.*sigma[i][j])
    /tmp2;
}
for (k=0; k<anz; k++) {
delta=0.0;
max_u_r=0.0;
for (j=0; j<N; j++)
u[M][j]+=w*(srp[M][j]*(2*f[j]/(sigma[M][j]*M)+u[M-1][j])
                            +srm[M][j]*u[M-1][j]
                            +sphip[M][j]*u[M][phi(j+1,N)]+sphim[M][j]*u[M][phi(j-1,N)]
                            -u[M][j]);
for (i=M-1; i>0; i--) {
u[i][0]+=w*(srp[i][0]*u[i+1][0]+srm[i][0]*u[i-1][0]
                            +sphip[i][0]*u[i][1]+sphim[i][0]*u[i][N-1]
                            -u[i][0]);
for (j=1; j<N-1; j++)
          u[i][j]+=w*(srp[i][j]*u[i+1][j]+srm[i][j]*u[i-1][j]
                                  +sphip[i][j]*u[i][j+1]+sphim[i][j]*u[i][j-1]
                                  -u[i][j]);
j=N-1;
u[i][j]+=w*(srp[i][j]*u[i+1][j]+srm[i][j]*u[i-1][j]
                            +sphip[i][j]*u[i][0]+sphim[i][j]*u[i][N-2]
                            -u[i][j]);
}
/* Mitte */

u[0][0]+=w*( ( (sigma[1][jj]-sigma[1][jj+NN]+8.*sigma[0][0])*u[1][jj]
             +(sigma[1][jj+NN]-sigma[1][jj]+8.*sigma[0][0])*u[1][jj+NN])
           /(16.*sigma[0][0])-u[0][0]);
for (j=1; j<N; j++)
u[0][j]=u[0][0];
jj++;
if (jj>=NN)
jj=0;
/*
u[0][0]=u[1][0];
for (j=1; j<N; j++)
u[0][0]+=u[1][j];
u[0][0]/=N;
for (j=1; j<N; j++)
u[0][j]=u[0][0];
*/
}
lift_polar(u,M,N);
}
void cit_start_polar() {
int i, j, l=1;
float n;
int k, ll, E2=E/2;
n=((float)N)/E;
for (k=0; k<E; k++) {
ll=k+1;
if (ll>=E)
ll-=E;

for (j=0; j<N; j++)
  f[k][1][j]=0;
f[k][1][(int)(n*k)]=1.;
f[k][1][(int)(n*ll)]=-1.;
for (i=0; i<=M; i++) for (j=0; j<N; j++)
  u[k][1][i][j]=v[k][1][i][j]=0.;
neumann_polar(u[k][1],f[k][1],s_org,M,N,Anzahl);
for (j=0; j<N; j++)
  g[k][1][j]=u[k][1][M][j];
}
for (k=0; k<E; k++) for (l=1; l<E2; l++) {
ll=l+1;
for (j=0; j<N; j++) {
  f[k][ll][j]=f[k][1][j]+f[k+1][l][j];
  g[k][ll][j]=g[k][1][j]+g[k+1][l][j];
}
for (i=0; i<=M; i++) for (j=0; j<N; j++) {
  u[k][ll][i][j]=u[k][1][i][j]+u[k+1][l][i][j];
  v[k][ll][i][j]=0.;
}
}
for (k=0; k<E; k++) for (l=1; l<=E2; l++) for (j=0; j<N; j++)
g[k][l][j]*=rauschFaktor();
}
void cit_turn_polar( int k, int l ) {
int i, j;
float M2=M*M, N2=N*N, Pi4_2=4*M_PI*M_PI, tmp, c;
float gg[JJ], h[I][J], ng2, nh2;
neumann_polar(u[k][l],f[k][l],sigma,M,N,dAnz);
checkAndProcess();
ng2=0;
for (j=0; j<N; j++) {
    gg[j]=u[k][l][M][j]-g[k][l][j];
    ng2+=gg[j]*gg[j];
}
n_u_g=sqrt(ng2);
if (ng2>0) {
neumann_polar(v[k][l],gg,sigma,M,N,dAnz);
checkAndProcess();
nh2=0;
for (j=0; j<N; j++) {
  h[M][j]=f[k][l][j]*gg[j]
            +0.25*N2/Pi4_2*(u[k][l][M][phi(j+1,N)]-u[k][l][M][phi(j-1,N)])
              *(v[k][l][M][phi(j+1,N)]-v[k][l][M][phi(j-1,N)]);
  nh2+=h[M][j]*h[M][j];
}
for (i=1; i<M; i++) for (j=0; j<N; j++) {
  h[i][j]=0.25*M2*( (u[k][l][i+1][j]-u[k][l][i-1][j])
                                       *(v[k][l][i+1][j]-v[k][l][i-1][j])
                                       +N2/(Pi4_2*i*i)*(u[k][l][i][phi(j+1,N)]-u[k][l][i][phi(j-1,N)])
                                       *(v[k][l][i][phi(j+1,N)]-v[k][l][i][phi(j-1,N)]));
  nh2+=h[i][j]*h[i][j];
}
c=Mult*M*ng2/nh2;
if (RAND)
  for (j=0; j<N; j++) {
            sigma[M][j]+=c*h[M][j];
            if (sigma[M][j]<min_sigma)
              sigma[M][j]=min_sigma;
  }
for (i=1; i<M; i++) for (j=0; j<N; j++) {
  sigma[i][j]+=c*h[i][j];
  if (sigma[i][j]<min_sigma)

    sigma[i][j]=min_sigma;
  if (sigma[i][j]>max_sigma)
    sigma[i][j]=max_sigma;
}
tmp=sigma[1][0];
for (j=1; j<N; j++)
  tmp+=sigma[1][j];
sigma[0][0]=tmp/N;
for (j=1; j<N; j++)
  sigma[0][j]=sigma[0][0];
}
}

Anhang B

Mathematische Grundlagen

Definition B.1 (Topologische Räume) Sei Ω ⊂ IR^{d}, d ∈ IN ein Gebiet.

• C^{k}(Ω) ist der Raum aller k-mal stetig differenzierbaren Funktionen f : Ω → l C.

• C_{0}^{∞}(Ω) ist der Raum aller unendlich oft stetig differenzierbaren Funktionen mit kompaktem Träger in Ω.

• D(Ω) ist der Raum C_{0}^{∞}(Ω) versehen mit folgendem Konvergenzbegriff (Topologie): Eine Funktionenfolge {ϕ_{n}}_{n=1}^{∞} in D(Ω) konvergiert genau dann gegen Null, wenn

– ∃ K ⊂ Ω kompakt, d.d. supp(ϕ_{n}) ⊂ K ∀n

– \frac{∂^{α}}{∂x^{α_{1}} … ∂x^{α_{d}}} ϕ_{n} konvergiert gleichmäßig gegen Null ∀α ∈ IN^{d}.

• D′(Ω) ist der Dualraum von D(Ω), d.h. der Raum aller stetigen linearen Funktionale l : D(Ω) → l C. Man nennt l Distribution.

• S(IR^{d}) ist der Schwartzsche Raum der schnell fallenden C^{∞}-Funktionen. S′(IR^{d}) bezeichnet den zugehörigen Dualraum.

• L^{2}(Ω) ist der Raum aller Funktionen f : Ω → l C mit

∥f∥_{L^{2}(Ω)} := (∫_{Ω} |f|^{2} dx)^{1/2} < ∞.

Versehen mit dem Skalarprodukt

⟨f, g⟩_{L^{2}(Ω)} := ∫_{Ω} f\overline{g} dx

ist L^{2}(Ω) ein Hilbertraum.

Definition B.2 (Fouriertransformation) Die Fouriertransformation f̂ einer Funktion f ∈ S(IR^{d}) ist gegeben durch

f̂(k) := Ff(k) := (2π)^{−d/2}\begin{matrix}∫ \\ IR & ^{d} \end{matrix} e^{−ix·k}f(x) dx.

Die inverse Fouriertransformation lautet

f̃(k) := F^{−1}f(k) := (2π)^{−d/2}\begin{matrix}∫ \\ IR & ^{d} \end{matrix} e^{ix·k}f(x) dx.

Satz B.3 Da S′(IR^{d}) der Dualraum von S(IR^{d}) ist, läßt sich die Fouriertransformation auch in diesem Raum erklären. Der zu F duale Operator ist wieder F selber. Sei nun T ∈ S′(IR^{d}) ein lineares Funktional, dann ist FT gegeben durch

⟨FT, ϕ⟩ = ⟨T, Fϕ⟩, ∀ϕ ∈ S(IR^{d}).

Beweis: Siehe Rauch [R2], Satz 2.4.5. □

Korollar B.4 Wegen L^{2}(IR^{d}) ⊂ S′(IR^{d}) ist die Fouriertransformation auch hier erklärt: Sei f ∈ L^{2}(IR^{d}), dann ist Ff gegeben durch

Ff(k) = \frac{1}{(2π)^{d/2}} lim_{R→∞} ∫_{|x|<R} e^{−ikx}f(x) dx.

Hier ist Konvergenz im Sinne von L^{2} gemeint.

Bemerkung B.5 Es gibt keine einheitliche Position der Vorfaktoren (2π)^{−d/2} in der Fouriertransformation.

Definition B.6 (Sobolevräume) Der Sobolevraum H^{k}(Ω), k ∈ IN ist definiert als Vervollständigung von C^{k}(Ω) versehen mit der Norm

∥f∥_{H^{k}(Ω)}^{2} := ∑_{|α|≤k} ∥D^{α}f∥_{L^{2}(Ω)}^{2}.

Satz B.7 (Plancherel) Für f ∈ L^{2}(IR) gilt:

∥Ff∥_{L^{2}(Ω)} = ∥f∥_{L^{2}(Ω)}.

Definition B.8 (Sobolevräume gebrochener Ordnung) Für s ∈ IR ist H^{s}(IR^{d}) der Raum aller Distributionen f ∈ S′(IR^{d}) mit

∥f∥_{s} := (\begin{matrix}∫ \\ IR & ^{d} \end{matrix} (1 + |k|^{2})^{s}|Ff(k)|^{2} dk)^{1/2} < ∞.

Bemerkung B.9 Aus dem Satz von Plancherel folgt, daß beide Definitionen von H^{k} für k ∈ IN übereinstimmen.

Bemerkung B.10 Versehen mit dem Skalarprodukt

⟨f, g⟩_{H^{s}} :=\begin{matrix}∫ \\ IR & ^{d} \end{matrix} (1 + |k|^{2})^{s}f̂(k)¯ĝ(k) dk

ist H^{s} ein Hilbertraum. Aber mit dem L^{2}-Skalarprodukt ist H^{−s} der Dualraum von H^{s}.

Satz B.11 (Randwerte) Für Ω ⊂ IR^{d} mit C^{m}-Rand und f ∈ H^{s}(Ω) gilt:

D^{α}f|_{∂Ω} ∈ H^{s−|α|−1/2}(∂Ω) für |α| < m.

Beweis: Siehe Adams [A]. □

Satz B.12 (Gauß’scher Integralsatz) Seien Ω ⊂ IR^{d} ein hinreichend glattes beschränktes Gebiet und F ein stetig differenzierbares Vektorfeld. Dann gilt

∫_{Ω} ∇ · F dx = ∫_{∂Ω} F · ν dS.

Dabei bezeichnet ν die äußere Normale an Ω und dS das übliche Oberflächenmaß.

Beweis: Siehe Forster, Analysis 3, Satz 15.3. □

Korollar B.13 Seien u, v ∈ C^{2}(Ω) und σ ∈ C^{1}(Ω). Dann folgt aus dem Gauß’- schen Integralsatz

a_{σ}(u, v) := ∫_{Ω} v∇ · σ∇u dx = ∫_{∂Ω} vσ \frac{∂u}{∂ν} dS − ∫_{Ω} σ∇v · ∇u dx.

Deshalb läßt sich a_{σ}(u, v) auf H^{1}(Ω) fortsetzen mit σ ∈ L^{∞}(Ω).

Literaturverzeichnis

[A] Robert A. Adams Sobolev Spaces Academic Press, New York San Fransisco London 1975

[B-K] H.T. Banks, K. Kunisch Estimations Techniques for Distributed Parameter Systems Birkhäuser Boston, 1989

[B-B] D.C. Barber, B.H. Brown Applied potential tomography J.Phys.E.Sci.Instrum., Vol.17, 1984

[B-C-T] Carlos A. Berenstein, Enrico Casadio Tarabusi Inversion Formulas for the k-Dimensional Radon Transformation in Real Hyperbolic Spaces to appear in Duke Math. J.

[B1] Gregory Beylkin The Inversion Problem and Applications of the Generalized Radon Transform Cumm.Pure Appl.Math. 37 (1984), pp 579-599

[B2] William Breckon A Preliminary Report on the Reconstruction Problem Januar 1986

[B3] Krut Bryan Numerical recovery of certain discontinuous electrical conductivities Inverse Problems 7 (1991), pp 827-840

[C-I] Margaret Cheney, David Isaacson, Eli L. Isaacson Exact solutions to a linearized inverse boundary value problem Inverse Problems 6 (1990), pp 923-934 [C] Edward B. Curtis, James A. Morrow

Determining the resistors in a network SIAM J.Appl.Math. Vol.50,3 (1990), pp 918-941

[F-J-R] E.B. Fabes, M. Jodeit, N.M. Rivière, Potential techniques for boundary value problems on C^{1}-domains Acta Math. 141 (1978), pp 165-186

[G] Semion Gutman Alternating directions method for impedance computed tomography Inverse Problems 5 (1989), pp 1027-1047

[K-M] Robert V. Kohn und Alan McKenny Numerical implementation of a variational method for electrical impedance tomography Inverse Problems 6 (1990), pp 389-414

[K-V1] Robert Kohn, Michael Vogelius Determining Conductivity by Boundary Measurements Comm.Pure Appl.Math. 37 (1984), pp 289-298

[K-V2] Robert Kohn, Michael Vogelius Determining Conductivity by Boundary Measurements II. Interior Results Comm.Pure Appl.Math. 38 (1985), pp 643-667

[L] A.K. Louis Medical imaging: state of the art and future development Inverse Problems 8 (1992), pp 709-738

[M-G] Stephen McCormick, J. Gordon Multigrid solution of a linearized, regularized least-squares problem in electrical impedance tomography Inverse Problems 9 (1993), pp 637-713

[N] Adrian I. Nachman Reconstructions from boundary measurements Annals of Math. 128 (1988), pp 513-576

[FN1] Frank Natterer A constructive solution to the inverse scattering problem of the wave equation at one frequency and plane wave irradiation August 1989

[FN2] Frank Natterer Numerical Solution of Bilinear Inverse Problems Juli 1995 [N-G-I] J.C. Newell, David G. Gisser, David Isaacson

An Electrical Current Tomograph IEEE Transactions on biomedical engineering, Vol.35, No.10, 1988

[P] Kevin Paulson, William Breckon, Michael Pidock Electrode Modelling in Electrical Impedance Tomography

[R1] A.G. Ramm Recovery of the potential from fixed-energy scattering data Inverse Problems 4 (1988), pp 877-996

[R2] Jeffrey Rauch Partial Differential Equations Springer-Verlag, 1991

[S] Fadil Santosa und Michael Vogelius A backprojection algorithm for electrical impedance imaging SIAM Journal of Applied Mathematics 50 (1990), pp 216-243

[S-C-I-I] Erkki Somersalo, Margaret Cheney, David Isaacson, Eli Isaacson Layer stripping: a direct numerical method for impedance imaging Inverse Problems 7 (1991), pp 899-926

[S-U] John Sylvester, Gunther Uhlmann A global uniqueness theorem for an inverse boundary value problem Annals of Mathematics 125 (1987), pp 153-169

[Y] Thomas J. Yorkey, John G. Webster, Willis J. Tompkins Comparing reconstruction algorithms for electrical impedance tomography IEEE Transactions on biomedical engineering, Vol.34, No.11, 1987 Hiermit versichere ich, diese Diplomarbeit eigenständig und ohne unzulässige Hilfsmittel erstellt zu haben.