Schlecht gestellte Probleme bei der Wärmeleitungsgleichung

Diplomarbeit von Dieter Trautwein

Angefertigt nach einer Aufgabenstellung von Prof. Dr. F. Natterer, Institut für Numerische und Instrumentelle Mathematik der Westfälischen Wilhelms–Universität. Münster, im März 1991

Inhaltsverzeichnis

1 Grundlagen 1 1.1 Einleitung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Maximumprinzip . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Fundamentallösung . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4 Hilfsmittel aus der Funktionalanalysis . . . . . . . . . . . . . . 14

2 Allgemeine Theorie der Probleme 15

2.1 Ein Wärmeleitungsvorgang mit umgekehrter

Zeitrichtung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Ein Cauchy Problem mit negativer räumlicher Ausdehnung . . 18 2.3 Ein nichtcharakteristisches Cauchy Problem . . . . . . . . . . 25 2.4 Das inverse, eindimensionale Ein–Phasen Stefan Problem … 31 2.5 Parameteridentifikation . . . . . . . . . . . . . . . . . . . . . . 39

3 Regularisierungsmethoden 42 3.1 Allgemeine Theorie . . . . . . . . . . . . . . . . . . . . . . . . 42 3.2 Die Methode des ”quasi–reversibilité“ . . . . . . . . . . . . . . 45 3.3 Ein modifiziertes Cauchy Problem . . . . . . . . . . . . . . . . 50 3.4 Die Tikhonov–Regularisierung . . . . . . . . . . . . . . . . . . 54

4 Numerische Verfahren 61

4.1 Das diskrete Problem und die Kondition

der Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.2 Regularisierung des Diskreten Problems und

numerische Ergebnisse . . . . . . . . . . . . . . . . . . . . . . 64

5 Zusammenfassung 85

1 Grundlagen

1.1 Einleitung

Im Titel dieser Arbeit werden zwei Aspekte der angewandten Mathematik angesprochen. Es sind dies die schlecht gestellten Probleme und die Wärmeleitungsgleichung. Der französische Mathematiker Hadamard (1885–1963) definierte ein Problem als gut gestellt, falls gilt:

1. eine Lösung des Problems existiert,

2. die Lösung ist eindeutig,

3. die Lösung des Problems hängt stetig von den Daten ab.

Erfüllt ein Problem einen der drei Punkte nicht, so wird es als schlecht gestellt bezeichnet. Hadamard vertrat die Ansicht, daß derartige Probleme nicht von größerem Interesse seien, da sie meistens keine wichtigen physikalischen Anwendungen besäßen, doch es zeigte sich, daß viele im Sinne von Hadamard schlecht gestellte Probleme einen realen physikalischen Hintergrund besitzen. Ein Beispiel dafür liefern die in dieser Arbeit angesprochenen Aufgaben bei der Wärmeleitungsgleichung. Um eine Differentialgleichung für den Wärmefluß herzuleiten, benötigt man zwei Gesetze aus der Physik. Bezeichnet man mit q die Wärmestromdichte, also die Größe des Wärmeflusses während einer Zeiteinheit durch einen Einheitskörper, so besagt das erste Gesetz

q = −k∇u ,  (1.1)

wobei k eine Proportionalitätskonstante und ∇u der Gradient der Temperatur u ist. In einer räumlichen Variablen ergibt sich so

q = −ku_{x} .  (1.2)

Das zweite Gesetz besagt, daß die Veränderung der enthaltenen Wärmeenergie ΔQ direkt proportional zur Temperaturveränderung Δu und dem Maß ρ des Mediums ist. Also

ΔQ = cρΔu  (1.3)

mit der Kapazität c als Proportionalitätskonstante. Gegeben sei nun das Rechteck

R = {(ξ, τ) : x − Δx < ξ < x + Δx und t − Δt ≤ τ ≤ t + Δt} . (1.4)

Der Zuwachs der enthaltenen Wärmeenergie im Zeitintervall 2Δt im Bereich [x − Δx, x + Δx] wird somit nach (1.3) gegeben durch

cρ ∫_{x−Δx}^{x+Δx}\{u(ξ, t + Δt) − u(ξ, t − Δt)\}dξ = cρ ∫∫_{R} \frac{∂u}{∂τ}dτ dξ .  (1.5)

Dieser Energiezuwachs kann aber nur über die Differenz des Energieflusses an den Rändern x − Δx und x + Δx erfolgen. Dies ergibt nach (1.2)

k ∫_{t−Δt}^{t+Δt} −\{\frac{∂u}{∂x}(x − Δx, τ) − \frac{∂u}{∂x}(x + Δx, τ)\}dτ = k ∫∫_{R} \frac{∂^{2}u}{∂ξ^{2}}dτ dξ . (1.6)

Mit dem Energieerhaltungsprinzip schließt man

∫∫_{R}{cρu_{τ} − ku_{ξξ}}dτ dξ = 0 ,  ∀ Δx, Δt .  (1.7)

Aus (1.7) folgt

cρu_{τ} − ku_{xx} = 0  (1.8)

oder mit k_{1} = c^{−1}ρ^{−1}k

u_{τ} − k_{1}u_{xx} = 0 .  (1.9)

Durch Substitution von t = k_{1}τ ergibt sich dann die klassische Form der sogenannten Wärmeleitungsgleichung

Lu = u_{xx} − u_{t} = 0 .  (1.10)

1.2 Maximumprinzip

Ein wichtiges Hilfsmittel bei der Betrachtung von Problemen, die mit der Wärmeleitungsgleichung zusammenhängen, ist das sogenannte Maximumprinzip.

1.2.1 Satz. (Schwaches Maximumprinzip) Sei D_{T} := {(x, t) : s_{1}(t) < x < s_{2}(t), 0 < t ≤ T } für T > 0 ein Gebiet in R I ^{2} mit stetigen Kurven s_{1}, s_{2}. B_{T} sei der Rand von D_{T}, also

B_{T} := {(s_{1}(t), t), 0 ≤ t ≤ T } ∪ {(s_{2}(t), t), 0 ≤ t ≤ T } ∪ {(x, 0), a ≤ x ≤ b}

mit s_{1}(0) = a und s_{2}(0) = b. Ist u nun eine auf D_{T} ∪ B_{T} stetige Lösung von Lu = 0, so gilt:

max_{D_{T} ∪B_{T}} u = max_{B_{T}} u .

Beweis. Setzt man v(x, t) = u(x, t) + εx^{2} für ε > 0, so nimmt v sein Maximum auf B_{T} an. Begründung: Falls nicht, so existiert ein Punkt (x_{0}, t_{0}) ∈ D_{T} mit

v(x_{0}, t_{0}) = max_{D_{T} ∪B_{T}} v .

In (x_{0}, t_{0}) gilt weiter:

Lv ≤ 0 , da v_{t}(x_{0}, t_{0}) ≥ 0 und v_{xx}(x_{0}, t_{0}) ≤ 0 .

Dieses ist aber ein Widerspruch zu

Lv = Lu + Lεx^{2} = 2ε > 0 in ganz D_{T} .

Da u ≤ v und v ≤ max_{B_{T}} v , ist

u ≤ max_{B_{T}} v ≤ max_{B_{T}} u + ε max_{B_{T}} x^{2} ,

und mit beliebig kleinem ε folgt dann

u ≤ max_{B_{T}} u .

□

Mit −u statt u beweist man unter Anwendung von Satz 1.2.1 folgenden Satz.

1.2.2 Satz. Für eine Lösung u von Lu = 0 in einem beschränkten Gebiet D_{T}, die stetig auf D_{T} ∪ B_{T} ist, gilt:

min_{D_{T} ∪B_{T}} u = min_{B_{T}} u .

Die Sätze 1.2.1 und 1.2.2 erlauben den Beweis des sogenannten Vergleichsprinzips.

1.2.3 Satz. Seien u und v auf D_{T} ∪ B_{T} stetige Lösungen von (1.10), für die gilt: u ≤ v auf ganz B_{T} . Dann gilt:

u ≤ v in D_{T} ∪ B_{T} .

Beweis. Man setze w = u − v, so folgt Lw = 0 in D_{T} und w ≤ 0 auf B_{T}. Mit dem Maximumprinzip ergibt sich: w ≤ 0 in D_{T} ∪ B_{T} . □

Mit Hilfe der obigen Sätze, ist man nun in der Lage, einen Eindeutigkeitssatz für Lösungen der Gleichung (1.10) zu beweisen.

1.2.4 Satz. Sind u und v auf D_{T} ∪ B_{T} stetige Lösungen von Lu = 0 bzw. Lv = 0 in D_{T}, und gilt weiter u = v auf B_{T}, so ist

u = v in D_{T} ∪ B_{T} .

Beweis. Mit w = v − u ist Lw = 0 in D_{T} und w = 0 auf B_{T}. Satz 1.2.3 angewandt auf w und 0 ergibt dann w = 0 auf D_{T} ∪ B_{T} . □

Die folgenden Sätze werden sich im Verlauf der Arbeit noch als hilfreich erweisen. Da ihre Beweise jedoch einiges an Vorarbeit benötigen und somit vom Thema abweichen, sollen die Sätze hier nur zitiert werden. Die Beweise sind nachzulesen in §15 von [2]. Die ersten beiden Sätze beschäftigen sich mit der Frage der Eindeutigkeit von Lösungen der Differentialgleichung bei unstetigen Rändern.

1.2.5 Satz. Seien u und v Lösungen von (1.10) in D_{T}, die stückweise stetig auf D_{T} ∪ B_{T} sind. Die Anzahl der Unstetigkeitsstellen sei endlich und die Funktionen dort beschränkt. Ist weiter, außer an den Unstetigkeitsstellen, u ≤ v auf B_{T}, so folgt:

u ≤ v in D_{T} .

1.2.6 Satz. Es seien u und v wieder Lösungen der Differentialgleichung (1.10) in D_{T} mit maximal endlich vielen, beschränkten Unstetigkeitsstellen auf B_{T} ∪ D_{T} . Sei, außer an diesen Stellen, u = v auf B_{T}, dann gilt:

u = v in D_{T} .

Der nächste Satz ist eine Verschärfung des schwachen Maximumprinzips.

1.2.7 Satz. (Starkes Maximumprinzip) Die Funktion u(x, t) erfülle die Wärmeleitungsgleichung (1.10) in D_{T} und sei stetig auf D_{T} ∪ B_{T}, dann gilt für jedes t, 0 < t ≤ T, entweder

inf_{B_{t}} u < u(x, t) < sup_{B_{t}} u , (x, t) ∈ D_{T} ,

oder

u ≡ konstant in D_{t} ∪ B_{t} .

Der letzte Satz beschäftigt sich mit dem Verhalten von u_{x} auf den Randpunkten, in denen u sein Maximum oder Minimum annimmt.

1.2.8 Satz. (Hopf) Seien die Randkurven von D_{T}, s_{1}(t) und s_{2}(t), lipschitzstetig. Weiter sei u eine Lösung von (1.10) in D_{T} und stetig auf D_{T} ∪ B_{T}. Ist u nicht in jedem D_{t}, 0 < t ≤ T, konstant und nimmt sein \left\{{Maximum \atop Minimum}\right\} in (s_{2}(t_{0}), t_{0}) an, so ist:

lim_{x↑s_{2}(t_{0})} \left\{{inf \atop sup}\right\} (\frac{u(x, t_{0}) − u(s_{2}(t_{0}), t_{0})}{x − s_{2}(t_{0})}) > < 0 .

Falls also u_{x} existiert, ist

u_{x}(s_{2}(t_{0}), t_{0}) > < 0 .

1.3 Fundamentallösung

Um zu einer ersten Lösung für (1.10) zu gelangen, betrachte man die folgende Aufgabe. Gesucht ist eine Funktion u(x, t), die der Differentialgleichung und der Anfangsbedingung

\begin{aligned}u_{xx}(x, t) & = & u_{t}(x, t) , & −∞ < x < ∞ , t > 0 , \\ u(x, 0) & = & f(x) , & −∞ < x < ∞ \end{aligned}  (1.11)

genügt. Die Bedingung u(x, 0) = f(x) ist hierbei als Grenzwertprozeß zu verstehen, also: lim_{t→0^{+}} u(x, t) = f(x) ∀x.

\overline{}^{x} u_{t} = u_{xx}^{u = f}

Eine Standardmethode, um zu einer Darstellung für eine Lösung zu kommen, ist die Anwendung der Fourier–Transformation. Sei

û(s, t) = (2π)^{−1/2} ∫_{−∞}^{∞} u(x, t) e^{−isx}dx  (1.12)

die Fourier–Transformierte von u. Die Inversion ist gegeben durch

u(x, t) = (2π)^{−1/2} ∫_{−∞}^{∞} û(s, t) e^{isx}ds .  (1.13)

Bei Anwendung der Transformation auf die Differentialgleichung u_{xx} = u_{t} erhält man im Bildraum

û_{t}(s, t) = −s^{2}û(s, t) .  (1.14)

Dieses Ergebnis ist nur unter einigen Annahmen an u zu erreichen; so sollten u, u_{x} und u_{xx} sowie u_{t} fouriertransformierbar sein. Speziell muß gelten:

u_{x}(x, t) ⟶^{|x|→∞} 0, u(x, t) ⟶^{|x|→∞} 0.

Nimmt man weiter an, daß der Grenzwertprozeß lim_{t→0^{+}} u(x, t) = f(x) mit der Fourier–Transformation vertauschbar ist, so erhält man

û(s, 0) = f̂(s) .  (1.15)

(1.14) und (1.15) bilden eine gewöhnliche Differentialgleichung 1.Ordnung in t, deren Lösung sich sofort ergibt als

û(s, t) = f̂(s)e^{−s^{2}t} .  (1.16)

Wendet man nun den Faltungssatz für Fourier–Transformationen an und benutzt, daß

(\frac{1}{\sqrt{2t}}e^{− \frac{x^{2}}{4t}} )^{∧} = e^{−s^{2}t} , t > 0 ,  (1.17)

(siehe [17]), so folgt

u(x, t) = ∫_{−∞}^{∞} K(x − ξ, t)f(ξ)dξ ,  (1.18)

mit

K(x, t) = \frac{1}{\sqrt{4πt}} exp(−\frac{x^{2}}{4t}) ,  für t > 0 .  (1.19)

K(x, t) wird auch als Fundamentallösung der Wärmeleitungsgleichung bezeichnet. Im folgenden sollen einige Eigenschaften dieser Lösung gezeigt werden.

1.3.1 Satz. (Eigenschaften der Fundamentallösung)

A. K(x, t) > 0 für t > 0.

B. Für jedes feste t > 0 konvergieren K und sämtliche Ableitungen von K exponentiell gegen Null, für |x| gegen unendlich.

C. Für jedes feste δ > 0 ist lim_{t→0^{+}} K(x, t) = 0, gleichmäßig für alle |x| ≥ δ.

D. Für jedes feste δ > 0 ist lim_{t→0^{+}} ∫_{|x|≥δ} K(x, t)dx = 0.

E. Für jedes t > 0 ist ∫_{−∞}^{∞} K(x, t)dx = 1.

F. lim_{x→0^{−}} − ∫_{0}^{t}(∂K/∂x)(x, t − τ)dτ = −1/2.

G. lim_{x→0^{+}} − ∫_{0}^{t}(∂K/∂x)(x, t − τ)dτ = +1/2.

H. Für t > 0 ist K eine analytische Funktion in x und t.

I. K und sämtliche Ableitungen von K erfüllen für t > 0 die Differentialgleichung (1.10).

Beweis. Eigenschaft A ist sofort ersichtlich. Für B zeigt man, daß sich für t > 0 alle Ableitungen von K darstellen lassen als

K_{x^{n}t^{m}} = R_{m,n}(x, t^{1/2})K ,

wobei R_{m,n} eine rationale Funktion in x und t^{1/2} ist. Für |x| → ∞ folgt die Behauptung. C: Für |x| ≥ δ ist

K(x, t) ≤ \frac{1}{\sqrt{4πt}} exp(−\frac{δ^{2}}{4t})

und somit

|K_{x^{n}t^{m}}| ≤ |R_{m,n}|\frac{1}{\sqrt{4πt}} exp(−\frac{δ^{2}}{4t}) .

Es ergibt sich, daß für |x| ≥ δ, K und K_{x^{n}t^{m}} (für t gegen unendlich) gleichmäßig gegen null konvergieren. Zum Beweis von E führt man eine Variablentransformation ρ = x/\sqrt{4t} durch und erhält

\frac{1}{\sqrt{4πt}} ∫_{−∞}^{∞} exp(−\frac{x^{2}}{4t})dx = \frac{1}{\sqrt{π}} ∫_{−∞}^{∞} exp(−ρ^{2})dρ = 1 .  (1.20)

Für Eigenschaft D folgt mit derselben Transformation

∫_{|x|≥δ} K(x, t)dx = \frac{1}{\sqrt{π}} ∫_{|ρ|≥δ/\sqrt{4t}} exp(−ρ^{2})dρ .

Da δ/\sqrt{4t} gegen unendlich geht für t gegen Null, folgt D aus der Konvergenz des zweiten Integrals in (1.20).

Um die Eigenschaften F und G zu beweisen, führt man die Transformation ρ = x/2\sqrt{t − τ} durch und erhält

− ∫_{0}^{t} \frac{∂K}{∂x} (x, t − τ)dτ = \{ \frac{1}{\sqrt{π}} ∫_{x/(2\sqrt{t})}^{∞} exp(−ρ^{2})dρ, x > 0 , ; −\frac{1}{\sqrt{π}} ∫_{−∞}^{x/(2\sqrt{t})} exp(−ρ^{2})dρ, x < 0 .

Mit x → 0^{+} bzw. x → 0^{−} folgen so die Behauptungen, da

∫_{−∞}^{∞} exp(−ρ^{2})dρ = \sqrt{π} .

Der Beweis von Eigenschaft H ist sofort ersichtlich, da K das Produkt einer analytischen Funktion mit einer Komposition von analytischen Funktionenen ist. I ist ebenso offensichtlich, da K die Differentialgleichung erfüllt, und da für eine beliebig oft stetig differenzierbare Funktion v(x, t), die der Gleichung (1.10) genügt, v_{x} und v_{t} ebenfalls dieser Gleichung genügen. □

Um unabhängig von den Vorraussetzungen für die Fourier–Transformation zu zeigen, daß (1.18) eine Lösung von (1.11) darstellt, soll als nächstes gezeigt werden, für welche f das Faltungsintegral existiert.

1.3.2 Lemma. Für alle stückweise stetigen Funktionen f mit |f(x)| ≤ c_{1} exp(c_{2}|x|^{1+α}) , 0 ≤ α < 1 (c_{1}, c_{2} konstant), konvergiert

u(x, t) = ∫_{−∞}^{∞} K(x − ξ, t)f(ξ)dξ , für alle x und t > 0 .

Verhält sich f asymptotisch wie c_{1} exp(c_{2}x^{2}) für |x| gegen unendlich, so konvergiert

u(x, t) = ∫_{−∞}^{∞} K(x − ξ, t)f(ξ)dξ , für 0 < t < (1/4c_{2}) .

Verhält sich f wie c_{1} exp(c_{2}x^{2+α}), α > 0, so divergiert das Integral.

Beweis. Unter Verwendung von 2ab ≤ ε^{−1}a^{2} + εb^{2}, ε > 0 folgt,

exp (−\frac{(1 + ε^{−1})x^{2}}{4t} ) exp (−\frac{(1 + ε)ξ^{2}}{4t} )

≤ exp (−\frac{(x − ξ)^{2}}{4t} ) ≤ exp (−\frac{(1 − ε^{−1})x^{2}}{4t} ) exp (−\frac{(1 − ε)ξ^{2}}{4t} ) .

D.h. es existieren zwei positive Funktionen F_{1}(x, t) und F_{2}(x, t) mit

F_{1}(x, t) exp (−\frac{(1 + ε)ξ^{2}}{4t} ) ≤ K(x − ξ, t) ≤ F_{2}(x, t) exp (−\frac{(1 − ε)ξ^{2}}{4t} ) .

Ist nun f(x) = c_{1} exp(c_{2}|x|^{2+α}), α > 0, so gilt:

K(x − ξ, t)f(ξ) ≥ F_{1}(x, t) c_{1} exp (c_{2}|ξ|^{2+α} − \frac{(1 + ε)ξ^{2}}{4t} ) .

Falls ξ groß genug ist, so ist c_{2}|ξ|^{α} ≥ (1 + ε)/(4t), und es folgt

∫_{−∞}^{∞} K(x − ξ, t)f(ξ)dξ = ∞

für alle Funktionen f, die asymtotisch zu c_{1} exp(c_{2}|x|^{2+α}) bei |x| = ∞ sind. Als nächstes sei f(x) = c_{1} exp(c_{2}x^{2}). Dann ist

K(x − ξ, t)f(ξ) ≤ F_{2}(x, t)c_{1} exp {(c_{2} − (1 + ε)/(4t))ξ^{2}} .

Da c_{2} − (1 − ε)/(4t) < 0 für 0 < ε < 1 und 0 < t < (1 − ε)/(4c_{2}), folgt mit 0 < t < (1/4c_{2})

∫_{−∞}^{∞} K(x − ξ, t)f(ξ)dξ < ∞,

für alle stückweise stetigen Funktionen f, die sich bei |x| = ∞ verhalten wie c_{1} exp(c_{2}x^{2}). Mit f(x) = c_{1} exp(c_{2}|x|^{1+α}), 0 ≤ α < 1, folgt als letztes

K(x − ξ, t)f(ξ) ≤ F_{2}(x, t)c_{1} exp (c_{2}|ξ|^{1+α} − \frac{(1 − ε)ξ^{2}}{4t} ) .

Nun ist aber für 0 < ε < 1 und |ξ| groß genug c_{2} − (4t)^{−1}(1 − ε)|ξ|^{1−α} < 0, und somit folgt für alle f mit |f(x)| ≤ c_{1} exp(c_{2}|x|^{1+α}), 0 ≤ α < 1,

∫_{−∞}^{∞} K(x − ξ, t)f(ξ)dξ < ∞.
□

Um zu zeigen, daß u(x, t) = ∫_{−∞}^{∞} K(x−ξ, t)f(ξ)dξ eine Lösung von (1.11) darstellt, ist es nötig, das Verhalten von u(x, t) für t → 0^{+} zu untersuchen.

1.3.3 Lemma. Für stückweise stetige Funktionen f mit |f(x)| ≤ c_{1} exp(c_{2}x^{2}) (c_{1}, c_{2} konstant), konvergiert u(x, t) für t → 0^{+} gleichmäßig gegen f(x) für alle x innerhalb einer kompakten Teilmenge der Stetigkeitsintervalle von f.

Beweis. Aus Eigenschaft E von Satz 1.3.1 folgt

f(x) = \frac{1}{\sqrt{4πt}} ∫_{−∞}^{∞} exp(−\frac{(x − ξ)^{2}}{4t})f(x)dξ .

Für festes x ergibt sich so

u(x, t) − f(x) = J_{1}(t) + J_{2}(t) + J_{3}(t)

mit

J_{1}(t) = \frac{1}{\sqrt{4πt}} ∫_{−∞}^{x−δ} exp(−\frac{(x − ξ)^{2}}{4t})[f(ξ) − f(x)]dξ ,
J_{2}(t) = \frac{1}{\sqrt{4πt}} ∫_{x−δ}^{x+δ} exp(−\frac{(x − ξ)^{2}}{4t})[f(ξ) − f(x)]dξ ,
J_{3}(t) = \frac{1}{\sqrt{4πt}} ∫_{x+δ}^{∞} exp(−\frac{(x − ξ)^{2}}{4t})[f(ξ) − f(x)]dξ .

Sei K eine kompakte Teilmenge eines Intervalles, in dem f stetig ist. Es folgt, daß f in K gleichmäßig stetig ist. Somit existiert für ein solches K und für jedes ε > 0 ein δ_{ε} > 0, so daß

|f(x) − f(ξ)| < \frac{ε}{3} , ∀x, ξ ∈ K mit |x − ξ| < δ_{ε} .

O.B.d.A. sei x − δ_{ε} ≤ ξ ≤ x + δ_{ε} in K. Man erhält mit Eigenschaft E aus Satz 1.3.1,

|J_{2}(t)| ≤ ∫_{x−δ_{ε}}^{x+δ_{ε}} K(x − ξ, t)|f(ξ) − f(x)|dξ ≤ \frac{ε}{3} ∫_{x−δ_{ε}}^{x+δ_{ε}} K(x − ξ, t)dξ ≤ \frac{ε}{3} .

(1.21)

Betrachtet man J_{3}(t) so ergibt sich

|J_{3}(t)| ≤ 2c_{1} \frac{1}{\sqrt{4πt}} ∫_{x+δ_{ε}}^{∞} exp(−\frac{(x − ξ)^{2}}{4t}) · exp(c_{2}ξ^{2}) dξ .

Substituiert man ρ = x − ξ/2\sqrt{t}, so folgt

|J_{3}(t)| ≤ \frac{2c_{1}}{\sqrt{π}} ∫_{−∞}^{−δ_{ε}/2\sqrt{t}} exp(−ρ^{2} + c_{2}(x − 2\sqrt{t}ρ)^{2})dρ .

Da 2ab ≤ a^{2} + b^{2}, folgt:

c_{2}(x − 2\sqrt{t}ρ)^{2} = c_{2}(x^{2} − 2x · 2\sqrt{t}ρ + 4tρ^{2}) ≤ 2c_{2}x^{2} + 8tρ^{2}

und somit

|J_{3}(t)| ≤ \frac{2c_{1}}{\sqrt{π}} exp(2c_{2}x^{2}) ∫_{−∞}^{−δ_{ε}/2\sqrt{t}} exp((8c_{2}t − 1)ρ^{2})dρ .

Damit ergibt sich für 0 < t < (16c_{2})^{−1}

|J_{3}(t)| ≤ c ∫_{−∞}^{δ_{ε}/2\sqrt{t}} exp(−\frac{1}{2}ρ^{2})dρ ,

mit c = (2c_{1}/\sqrt{π}) exp(2c_{2}A^{2}), wobei |x| ≤ A. Aus der Existenz von ∫_{−∞}^{0} exp(−1/2 ρ^{2})dρ folgt dann:

∃ T_{ε} so daß ∀ 0 < t < T_{ε} |J_{3}(t)| < \frac{ε}{3} ,  ∀ |x| < A .  (1.22)

Analog zu J_{3} beweist man für J_{1},

∃ T_{ε} so daß ∀ 0 < t < T_{ε} |J_{1}(t)| < \frac{ε}{3} ,  ∀ |x| < A .  (1.23)

Sei jetzt A so gewählt, daß K in |x| < A enthalten ist, so folgt für ein beliebiges ε > 0 mit (1.21), (1.22), (1.23):

∃ T_{ε} > 0 so daß ∀ 0 < t < T_{ε} |u(x, t) − f(x)| < ε , x ∈ K .

Erwähnt sei noch, daß durch die gleichmäßige Konvergenz von u(x, t) gegen f auf kompakten Teilmengen der Stetigkeitsintervalle und die stückweise Stetigkeit von f sich die zweidimensionale stückweise Stetigkeit von u bei t = 0 ergibt. Begr.: Für x_{1}, x_{2} ∈ K existiert ein δ_{ε}(x_{1}) > 0, so daß

|u(x_{2}, t) − f(x_{1})| ≤ |u(x_{2}, t) − f(x_{2})| + |f(x_{2}) − f(x_{1})| ≤ ε

für (x_{1} − x_{2})^{2} < (δ_{ε}(x_{1}))^{2} .

Somit ist u bei t = 0 stückweise zweidimensional stetig und nimmt dort die Werte von f an. □

Mit den beiden Lemmata ist es nun möglich, einen Existenzsatz für Lösungen der Aufgabe (1.11) zu formulieren.

1.3.4 Satz. (Existenz) Für alle stückweise stetigen Funktionen f mit |F (x)| ≤ c_{1} exp(c_{2}|x|^{1+α}) , 0 ≤ α < 1 , c_{1}, c_{2} positive Konstanten, ist

u(x, t) = ∫_{−∞}^{∞} K(x − ξ, t)f(ξ)dξ , t > 0 ,

eine Lösung von (1.11). Für alle stückweise stetigen Funktionen f, die bei |x| = ∞ asymptotisch zu c_{1} exp(c_{2}x^{2}) sind, ist u nur für 0 < t < (4c_{2})^{−1} definiert, und löst (1.11) nur für diese Werte.

Beweis. Aus Lemma 1.3.2 folgt die Konvergenz der Integrale. Da

\frac{∂^{m+n}K}{∂x^{n}∂t^{m}} (x, t) = R_{m,n}(x, t^{1/2})K(x, t) ,

mit einer rationalen Funktion R_{m,n} (siehe Bew. zu 1.3.1 B ), folgt aus der Leibniz–Regel für alle f, die den Bedingungen des Satzes genügen, daß u_{t} und u_{xx} für die entsprechenden t–Werte existieren und gleich sind. Die Grenzwerteigenschaften folgen aus Lemma 1.3.3 . □

Die eindeutige Lösbarkeit wird hier nicht mehr gezeigt, da sie für die eigentlichen Aufgaben nicht von Interesse ist. Sie ist nachzulesen bei [2] Seite 40 - 42.

1.4 Hilfsmittel aus der Funktionalanalysis

1.4.1 Definition. Eine lineare Abbildung A eines normierten Raumes E in einen normierten Raum F heißt kompakt oder vollstetig, falls das Bild Ax_{n} jeder beschränkten Folge x_{n} ∈ E eine konvergente Teilfolge in F enthält.

1.4.2 Satz. Sei (Af)(x) = ∫_{Ω_{1}} K(x, ξ)f(ξ)dξ , x ∈ Ω_{2}, Ω_{1}, Ω_{2} Teilmengen von R I ^{m_{1}}, R I ^{m_{2}}, ein Operator dessen Kern K in L_{2}(Ω_{1} × Ω_{2}) liegt. Dann gilt: A ist kompakt.

Beweis. Siehe [1] Seite 50–52.

1.4.3 Korollar. Ein kompakter Operator A : E → F hat nur dann eine stetige Inverse, wenn sein Definitionsraum von endlicher Dimension ist.

Begr. Sei E ein normierter Raum mit dim(E) = ∞. Dann existiert eine beschränkte Folge x_{n} in E die keine konvergente Teilfolge hat. Da A kompakt, besitzt y_{n} = Ax_{n} eine konvergente Teilfolge y_{n_{k}} → y. Es gilt aber, x_{n_{k}} = A^{−1}y_{n_{k}} ↛ A^{−1}y = x. □

2 Allgemeine Theorie der Probleme

2.1 Ein Wärmeleitungsvorgang mit umgekehrter Zeitrichtung

Unter den zahlreichen Problemen, die zur Formulierung einer Randwertaufgabe für die Wärmeleitungsgleichung führen, ist das Problem der Temperaturverteilung u(x, t) in einem Stab mit festen Randtemperaturen bei x = 0 und x = 1 und einer vorgegebenen Anfangstemperatur zur Zeit t = 0 das bekannteste. Gesucht wird hier die Temperaturverteilung u(x, T) zur Zeit t = T > 0. Mit Hilfe der in Abschnit 1.1 hergeleiteten Differentialgleichung für die Wärmeleitung (1.10) läßt sich die Aufgabe wie folgt formulieren: Gesucht wird bei gegebenem f(x) die Funktion g(x) = u(x, T) unter den Bedingungen

\begin{aligned}u_{xx}(x, t) & = & u_{t}(x, t) , & 0 < x < 1 , 0 < t ≤ T , \\ u(0, t) & = & u(1, t) = 0 , & 0 < t ≤ T , \\ u(x, 0) & = & f(x) , & 0 < x < 1 . \end{aligned}  (2.1)

Die Rand– und Anfangswertbedingungen sind hierbei wieder als Grenzprozesse zu verstehen.

[figure: x; u = f; t; 1; T; u = g; (0,0); u_{t} = u_{xx}; u = 0; u = 0]

Für viele Anwendungen ist es nun interessant, nicht das direkte Problem (2.1) zu betrachten, sondern eines mit umgekehrter Zeitrichtung. Man möchte aus Messungen der Temperatur des Stabes zur Zeit t = T die Anfangstemperatur u(x, 0) = f(x) bestimmen.

Gegeben ist also das Problem der Bestimmung von f(x) bei gegebenem g(x) unter den Bedingungen

\begin{aligned}u_{xx}(x, t) & = & u_{t}(x, t) , & 0 < x < 1 , 0 < t ≤ T , \\ u(0, t) & = & u(1, t) = 0 , & 0 < t ≤ T , \\ u(x, T) & = & g(x) , & 0 < x < 1 , T > 0 . \end{aligned}  (2.2)

Um eine Lösung für (2.2) herzuleiten, betrachte man zunächst das direkte Problem. Durch Fortsetzen von f(x) auf (−1, 0) durch −f(−x) und Definieren von: f(0) := 0, f(1) := 0 und

f̃(2k + x) := f(x) , für x ∈ (−1, 1], k = 0, ±1, ±2, ⋯,

ist

u(x, t) = ∫_{−∞}^{∞} K(x − ξ, t)f̃(ξ)dξ ,  (2.3)

für t > 0, eine Lösung der Anfangswertaufgabe (1.11) für dieses f̃. Durch Umformen der Lösung erhält man

u(x, t) = ∫_{−∞}^{∞} K(x − ξ, t)f̃(ξ)dξ  t > 0
= ∑_{n=−∞}^{∞} ∫_{n}^{n+1} K(x − ξ, t)f̃(ξ)dξ
= ∑_{n=−∞}^{∞} ∫_{0}^{1} K(x − (ξ + n), t)f̃(ξ + n)dξ
= ∫_{0}^{1} ∑_{n=−∞}^{∞} K(x − (ξ + 2n), t)f̃(ξ + 2n)dξ

+ ∫_{0}^{1} ∑_{n=−∞}^{∞} K(x − (ξ + 2n − 1), t)f̃(ξ + 2n − 1)dξ

= ∫_{0}^{1} ∑_{n=−∞}^{∞} K(x − (ξ + 2n), t)f̃(ξ)dξ

+ ∫_{0}^{1} ∑_{n=−∞}^{∞} K(x − (ξ − 1 + 2n), t)f̃(ξ − 1)dξ

= ∫_{0}^{1} ∑_{n=−∞}^{∞} K(x − (ξ + 2n), t)f(ξ)dξ

− ∫_{0}^{1} ∑_{n=−∞}^{∞} K(x + (1 − ξ) − 2n, t)f(1 − ξ)dξ

= ∫_{0}^{1} ∑_{n=−∞}^{∞} K(x − (ξ + 2n), t)f(ξ)dξ

− ∫_{0}^{1} ∑_{n=−∞}^{∞} K(x + ξ − 2n), t)f(ξ)dξ

u(x, t) = ∫_{0}^{1} ∑_{n=−∞}^{∞} {K(x − (ξ + 2n), t) − K(x + (ξ − 2n), t)}f(ξ)dξ .

Mit den Eigenschaften von K aus Satz 1.3.1 ist sofort ersichtlich, daß die Reihe für t > 0 gleichmäßig konvergent ist. Durch Aufteilen der Reihe erhält man, daß die Randbedingungen u(0, t) = u(1, t) = 0 erfüllt sind. Führt man den Grenzprozeß t → 0^{+} durch, so bleibt zu zeigen:

lim_{t→0^{+}} ∫_{0}^{1} K(x − ξ, t)f(ξ)dξ = f(x) , 0 < x < 1 .

Setzt man f außerhalb von [0, 1] mit Null fort, so ergibt sich diese Behauptung mit Lemma 1.3.3 . Die Eindeutigkeit der Lösung folgt aus dem Maximumprinzip 1.2.6 . Um nun zu einer Lösung des Problems (2.2) zu gelangen, betrachtet man u an der Stelle T > 0 und erhält

u(x, T) = ∫_{0}^{1} ∑_{n=−∞}^{∞} {K(x − (ξ + 2n), T) − K(x + (ξ − 2n), T)}f(ξ)dξ
= g(x) .  (2.4)

(2.4) ist eine Fredholm Integralgleichung erster Art mit glattem Kern. Betrachtet man die Gleichung (2.4) als Operatorengleichung (Af)(x) = g(x) , A : L_{2}[0, 1] → L_{2}[0, 1] , mit

Af(x) = ∫_{0}^{1} ∑_{n=−∞}^{∞} {K(x − (ξ + 2n), T) − K(x + (ξ − 2n), T)}f(ξ)dξ,

so ist A nach Satz 1.4.2 ein kompakter Operator, da der Kern von A aus L_{2}([0, 1] × [0, 1]) stammt. (2.2) ist somit ein schlecht gestelltes Problem, da nach Korollar 1.4.3 kompakte Operatoren auf unendlichdimensionalen Räumen keine stetige Inverse besitzen. Die Tatsache, daß es sich um einen Fredholm-Integraloperator erster Art handelt, ergibt sogar, daß das Problem sehr schlecht gestellt ist.

2.2 Ein Cauchy Problem mit negativer räumlicher Ausdehnung

Eine Vielzahl von Anwendungen der Wärmeleitungstheorie befaßt sich mit der Bestimmung von Temperaturen in einem Körper, der im Inneren nicht zugänglich ist. Man stelle sich eine sehr dicke Wand vor, deren eine Seite nicht erreichbar ist, da dort, wie z.B. bei einem Schmelztiegel, sehr hohe Temperaturen herrschen. Um die Temperatur im Inneren zu bestimmen, bohrt man nun, so weit es geht, ein Loch in die Wand. Mit Temperraturdaten aus diesem Loch versucht man nun, die Temperatur im Reststück der Wand, vom Lochende bis zur Innenseite des Tiegels, zu bestimmen. Auf eine räumliche Koordinate beschränkt, läßt sich das Problem wie folgt formulieren:

\begin{aligned}u_{t} & = & u_{xx} , & 0 < x < ∞ , 0 < t < ∞ , \\ u(x, 0) & = & 0 , & 0 ≤ x < ∞ , \\ u(1, t) & = & g(t) , & 0 ≤ t < ∞ . \end{aligned}  (2.5)

Weiter sei an die Temperaturausbreitung noch die zusätzliche Bedingung gestellt, daß u(x, t) < c, für eine Konstante c und x → ∞. Bei gegebenem g(t) ist nun das Temperaturverhalten an der Stelle x = 0, u(0, t) = f(t) zu bestimmen.

[figure: x; u = 0; t; 1; (0,0); u_{t} = u_{xx}; u = f; u = g]

Um eine Integralgleichung für die Lösung herzuleiten, betrachte man zunächst das direkte Problem der Bestimmung von u(1, t) = g(t) bei gegebenem f(t). Als Hilfsmittel benötigt man die Laplace–Transformation.

Ist

ĥ(s) = ∫_{0}^{∞} e^{−st}h(t)dt

die Laplace–Transformierte von h, so ergibt sich bei Anwendung der Transformation auf die Differentialgleichung (2.5)

û_{xx}(x, s) = sû(x, s) − u(x, 0) .  (2.6)

Da u(x, 0) = 0, erhält man im Bildraum die gewöhnliche Differentialgleichung

û_{xx}(x, s) = sû(x, s) .  (2.7)

Aus dem Ansatz û(x) = e^{αx} erhält man, mit α = ±\sqrt{s} , die beiden Partikulärlösungen e^{x\sqrt{s}} und e^{−x\sqrt{s}}, wobei \sqrt{s} der Hauptzweig der komplexen Wurzel ist. Eine allgemeine Lösung von (2.7) läßt sich also darstellen als

û(x, s) = c_{1}(s) exp(x\sqrt{s}) + c_{2}(s) exp(−x\sqrt{s}) .  (2.8)

Führt man einen rechten Rand bei x = l ein und sucht nach Lösungen û_{0}(x, s), û_{1}(x, s), für die gelten:

\begin{matrix}û_{0}(0, s) = 1 , & û_{0}(l, s) = 0 , \\ û_{1}(0, s) = 0 , & û_{1}(l, s) = 1 , \end{matrix}  (2.9)

so erhält man:

û_{0}(x, s) = \frac{exp ((l − x)\sqrt{s}) − exp (− (l − x)\sqrt{s})}{exp(l\sqrt{s}) − exp(−l\sqrt{s})} = \frac{sinh ((l − x)\sqrt{s})}{sinh(l\sqrt{s})}

û_{1}(x, s) =  \frac{exp(x\sqrt{s}) − exp(−x\sqrt{s})}{exp(l\sqrt{s}) − exp(−l\sqrt{s})}  = \frac{sinh(x\sqrt{s})}{sinh(l\sqrt{s})} .

Schreibt man in (2.9) allgemeine Datenfunktionen f̂(s) b.z.w. f̂_{l}(s) vor, so ist die zugehörige Lösung

û(x, s) = û_{0}(x, s)f̂(s) + û_{1}(x, s)f̂_{l}(s) .  (2.10)

Läßt man nun l → ∞ streben, so folgt:

û_{1}(x, s) ⟶^{l→∞} 0 und û_{0}(x, s) ⟶^{l→∞} exp(−x\sqrt{s})

und im Bildraum ergibt sich die Lösung als

û(x, s) = exp(−x\sqrt{s})f̂(s) .  (2.11)

Unter Ausnutzung der Beziehung

(\frac{x}{2\sqrt{π}t^{3/2}} exp(−\frac{x^{2}}{4t}))^{∧} = exp(−x\sqrt{s})  (2.12)

(siehe [7]) folgt mit dem Faltungssatz für Laplace–Transformationen

u(x, t) = (\frac{x}{2\sqrt{π}t^{3/2}} exp(−\frac{x^{2}}{4t})) ∗ f(t)
= \frac{x}{2\sqrt{π}} ∫_{0}^{t} \frac{f(τ)}{(t − τ)^{3/2}} exp (\frac{−x^{2}}{4(t − τ)}) dτ .  (2.13)

Vergleicht man das so erhaltene Ergebnis mit der Fundamentallösung aus Abschnitt 1.3, so erhält man

u(x, t) = −2 ∫_{0}^{t} \frac{∂K}{∂x} (x, t − τ)f(τ)dτ .  (2.14)

Da Formel (2.14) nur unter einer Reihe implizit vorgenommener Annahmen zustande gekommen ist, muß man noch unabhängig von diesen Annahmen zeigen, daß (2.14) eine Lösung des Problems darstellt.

2.2.1 Satz. Für stückweise stetige Funktionen f mit |f(t)| ≤ ct^{−α}, 0 < t < ε , 0 < α < 1 löst

u(x, t) = −2 ∫_{0}^{t} \frac{∂K}{∂x} (x, t − τ)f(τ)dτ ,  x > 0, t < 0 ,  (2.15)

das initialisierte Randwertproblem

\begin{aligned}u_{t} & = & u_{xx} , & 0 < x < ∞ , 0 < t , \\ u(0, t) & = & f(t) , & 0 < t , \\ u(x, 0) & = & 0 , & 0 < x < ∞ . \end{aligned}  (2.16)

Beweis. Aus dem Beweis zu Eigenschaft B von Satz 1.3.1 ersieht man, daß

\frac{∂^{m+n}K}{∂x^{m}∂t^{n}} (x, t − τ) = R_{m,n}(x, (t − τ)^{1/2}) · exp ( − \frac{x^{2}}{4(t − τ)}) ,  (2.17)

mit einer rationalen Funktion R_{n,m}. Es folgt nun für x > 0

lim_{τ ↑ t} \frac{∂^{m+n}K}{∂x^{m}∂t^{n}} (x, t − τ) = 0  (2.18)

und somit, daß (2.17) beschränkt ist durch Konstanten, die nur von m, n und x abhängen. Da für f(t) ein Wachstum nicht stärker als ct^{−α} , 0 < α < 1 , bei t = 0 gefordert wurde, existiert das Faltungsintegral (2.15), und u(x, t) ist wohldefiniert. Mit der Leibniz–Regel und den obigen Ergebnissen ist unter Anwendung von Eigenschaft I aus Satz 1.3.1 sofort ersichtlich, daß u die Differentialgleichung für x > 0, t > 0 erfüllt. Für x > 0 erhält man leicht, daß lim_{t→0^{+}} u(x, t) = 0 ist. Wenn sogar x ≥ δ > 0 gilt, so ist die Konvergenz gleichmäßig, denn aus e^{−x} ≤ p!x^{−p} x > 0, p = 1, 2, 3 ⋯, folgt für p = 2

2 |\frac{∂K}{∂x} (x, t − τ)| ≤ \frac{x · 2!}{2\sqrt{π}(t − τ)^{3/2}} (\frac{x^{2}}{4(t − τ)})^{−2}
= \frac{16(t − τ)^{1/2}}{\sqrt{π} x^{3}} ≤ 16π^{−1/2}δ^{3}(t − τ)^{1/2} ,

und somit

|u(x, t)| ≤ c_{1} ∫_{0}^{t} τ^{−α}(t − τ)^{1/2}dτ
≤ c_{1}t^{1/2} ∫_{0}^{t} τ^{−α}dτ
= c_{2}t^{(3−2α)/2} ,  0 < α < 1 , c_{1}, c_{2} konstant.

Als letztes bleibt die Konvergenz von u(x, t) gegen f(t) für x → 0 zu zeigen. Dazu betrachte man zunächst einen Stetigkeitspunkt t > 0 von f. Aus Satz 1.3.1 G folgt

lim_{x→0^{+}} f(t) · −2 ∫_{0}^{t} \frac{∂K}{∂x} (x, t − τ)dτ = f(t) .  (2.19)

Aus dem Beweis der Eigenschaft G von Satz 1.3.1 folgt mit t ≥ δ > 0, x > 0 :

1 + 2 ∫_{0}^{t} \frac{∂K}{∂x} (x, t − τ)dτ = 1 − \frac{2}{\sqrt{π}} ∫_{x/(2\sqrt{t})}^{∞} exp(−ρ^{2})dρ

= \frac{2}{\sqrt{π}} ∫_{0}^{x/(2\sqrt{t})} exp(−ρ^{2})dρ
≤ \frac{2}{\sqrt{π}} ∫_{0}^{x/(2\sqrt{δ})} exp(−ρ^{2})dρ
≤ π^{−1/2}δ^{−1/2}x .

Somit ist für t ≥ δ > 0 die Konvergenz in (2.19) gleichmäßig. Dies hat zur Folge, daß es im weiteren reicht zu zeigen, daß

lim_{x→0^{+}} \underbrace{−2 ∫_{0}^{t} \frac{∂K}{∂x} (x, t − τ)(f(τ) − f(t))dτ}_{=: J} = 0 .  (2.20)

Sei nun ε > 0, dann gibt es ein δ_{ε} > 0, so daß

|f(t) − f(τ)| < \frac{ε}{2} für |t − τ| ≤ δ_{ε},

da f stetig in t. Definiert man

J_{1} := −2 ∫_{0}^{t−δ_{ε}} \frac{∂K}{∂x} (x, t − τ)(f(τ) − f(t))dτ  und
J_{2} := −2 ∫_{t−δ_{ε}}^{t} \frac{∂K}{∂x} (x, t − τ)(f(τ) − f(t))dτ ,

so ist J = J_{1} + J_{2}. Mit −2(∂K/∂x) > 0 für x > 0 und der Substitution ρ = x/(2\sqrt{t − τ}) folgt

|J_{2}| ≤ \frac{ε}{2} \{−2 ∫_{t−δ_{ε}}^{t} \frac{∂K}{∂x} (x, t − τ)dτ \} = \frac{ε}{\sqrt{π}} ∫_{x/(2\sqrt{δ_{ε}})}^{∞} exp(−ρ^{2})dρ < \frac{ε}{2} ,

falls x klein genug ist. Für J_{1} gilt:

|J_{1}| ≤ \frac{x}{2\sqrt{π}} ∫_{0}^{t−δ_{ε}} |(t − τ)^{−3/2}| exp (\frac{−x^{2}}{4(t − τ)})|f(τ) − f(t)|dτ
≤ \frac{x}{2\sqrt{π}} δ_{ε}^{−3/2} ∫_{0}^{t−δ_{ε}} |f(τ) − f(t)|dτ
≤ \frac{x}{2\sqrt{π}} δ_{ε}^{−3/2} ∫_{0}^{T} |f(τ) − f(t)|dτ ,  0 < t < T .

Es ist also für

0 < x < \sqrt{π} δ_{ε}^{3/2} \{∫_{0}^{T} |f(τ) − f(t)|dτ \}^{−1} ε

(o.B.d.A. ∫_{0}^{T} |f(τ) − f(t)|dτ ≠ 0 )

|J_{1}| < \frac{ε}{2} .

Zusammen ergibt sich also, daß zu jedem ε > 0 ein γ_{ε} existiert, mit dem

|J| ≤ |J_{1}| + |J_{2}| < ε ,

für alle 0 < x < γ_{ε}. Daß die Konvergenz von u(x, t) → f(t) sogar gleichmäßig ist für t aus einer kompakten Teilmenge eines Stetigkeitsintervalles von f, ergibt sich, da man für solch ein Intervall δ_{ε} unabhängig von t wählen kann, da f lokal gleichmäßig stetig ist. □

Eine Eindeutigkeit der Lösung läßt sich im allgemeinen nicht beweisen, da z.B. u = ∂K/∂x eine Lösung mit u(0, t) = 0 und u(x, 0) = 0 ist. Fordert man jedoch, daß die Unstetigkeitsstellen auf dem Rand beschränkt sein sollen, so läßt sich für |u(x, t)| ≤ c_{1} exp(c_{2}x^{2}) ein Eindeutigkeitsbeweis führen. Der Beweis benutzt im wesentlichen das Vergleichsprinzip 1.2.6 und ist nachzulesen auf Seite 52 von [2].

Um das eigentliche Problem (2.5) zu Lösen, betrachte man u an der Stelle x = 1. Es ergibt sich

u(1, t) = ∫_{0}^{t} \frac{f(τ)}{2\sqrt{π}(t − τ)^{3/2}} exp (\frac{−1}{4(t − τ)})dτ = g(t).

Das Problem kann also wie folgt formuliert werden: bestimme f(t) von g(t) durch lösen der Volterraschen Integralgleichung 1.Art

(Af)(t) = g(t)

mit

(Af)(t) = \frac{1}{2\sqrt{π}} ∫_{0}^{t} \frac{f(τ)}{(t − τ)^{3/2}} exp (\frac{−1}{4(t − τ)})dτ .  (2.21)

Genaueres über den Operator (2.21) verrät der folgende Satz.

2.2.2 Satz. Für T > 0 ist der durch (2.21) definierte Operator A : L_{2}[0, T] → L_{2}[0, T] auf seinem Bildraum invertierbar. Die Inverse ist jedoch unstetig, d.h. A^{−1} ist unbeschränkt.

Beweis. Sei h(t) definiert als

h(t) = \{ \begin{matrix}\frac{1}{2\sqrt{π}t^{3/2}} exp(\frac{−1}{4t}) & t > 0 \\ 0 & t ≤ 0 . \end{matrix}

h(t) ist eine beschränkte stetige Funktion auf −∞ < t < ∞ und h(t) ∈ L_{2}(−∞, ∞). Mit diesem h läßt sich nun (2.21) als Fredholm– Integralgleichung umschreiben.

(Af)(t) = ∫_{0}^{T} h(t − τ)f(τ)dτ, 0 < t ≤ T.

Aus Abschnitt 1.4 folgt sofort, daß A kompakt ist und somit keine stetige Inverse besitzt. Um zu zeigen, daß A invertierbar ist, sei (Af)(t) = 0 f.ü. auf [0, T]. Mit der Laplace-Transformation von h(t) (siehe(2.12)) erhält man

(∫_{0}^{t} h(t − τ)f(τ)dτ )^{∧} = (0)^{∧}
⇒  (h(t))^{∧} · (f(t))^{∧}  = 0
⇒  e^{−\sqrt{s}} · f̂(s)  = 0 .

Es folgt f̂(s) = 0 und somit f(t) = 0 a.e. für 0 ≤ t < ∞ . □

Der Beweis hat gezeigt, daß A zwar invertierbar ist, aber die Inverse unstetig ist. Die Tatsache, daß sich der Operator in einen Fredholm–Operator 1.Art umschreiben läßt, verrät, wie schon beim ersten Problem, daß es sich hierbei um ein ernsthaft schlecht gestelltes Problem handelt.

2.3 Ein nichtcharakteristisches Cauchy Problem

Man betrachte einen Stab dessen Temperaturverhalten auf der einen Seite durch die Funktion g(t) = u(0, t) beschrieben wird. Weiter sei der Stab so isoliert, daß an dieser Seite keine Wärme hinzukommt oder abfließt, also u_{x}(0, t) = 0. Die Temperaturverteilung zur Zeit t = 0 sei vorgegeben, z.B. durch u(x, 0) = 0. Ist man nun am Temperaturverhalten auf der anderen Seite des Stabes interressiert, so stellt sich folgendes Problem: Bestimme f(t) = u(1, t) bei gegebenem g(t) unter den Bedingungen:

\begin{aligned}u_{t} & = & α u_{xx} , & 0 ≤ x < 1 ,  0 < t ≤ 1 , \\ u(0, t) & = & g(t) , & 0 < t ≤ 1 , \\ u_{x}(0, t) & = & 0 , & 0 < t ≤ 1 , \\ u(x, 0) & = & 0 , & 0 ≤ x < 1 , \end{aligned}  (2.22)

mit 0 < c_{1} ≤ α ≤ c_{2} < ∞.

[figure: x; u = 0; t; 1; 1; (0,0); u_{t} = αu_{xx}; \begin{matrix}u = g \\ u_{x} = 0 \end{matrix}; u = f]

(2.22) wird als nicht charakteristisches Cauchy Problem bezeichnet, da die Cauchy Daten u(0, t) = g(t), u_{x}(0, t) an einer Kurve (x = 0) gegeben sind, die keine Charakteristik der Differentialgleichung ist. (Die Charakteristiken der Diffgl. sind t = konstant.) Aufgaben mit allgemeinen Daten u(x, 0) = φ(x) und u_{x}(0, t) = h(t) lassen sich durch die Linearität der Probleme leicht im wesentlichen auf solche der Form (2.22) zurückführen. Mittels Laplace–Transformation soll auch hier versucht werden, formal eine Lösung herzuleiten. Dazu betrachte man zunächst folgendes Randwertproblem,

\begin{aligned}u_{t} & = & α u_{xx} , & 0 ≤ x < 1 , 0 < t ≤ 1 , \\ u_{x}(0, t) & = & 0 , & 0 < t ≤ 1 , \\ u(1, t) & = & f(t) , & 0 < t ≤ 1 , \\ u(x, 0) & = & 0 , & 0 ≤ x < 1 . \end{aligned}  (2.23)

Setzt man u(x, t) und f(t) mit 0 fort für t > 1, und transformiert die Differentialgleichung, so erhält man in Analogie zu (2.7)

αû_{xx}(x, s) = s · û(x, s) ,  (2.24)

und somit

û(x, s) = a_{1}(s) exp(\sqrt{\frac{s}{α}} x) + a_{2}(s) exp(−\sqrt{\frac{s}{α}} x) ,

wobei \sqrt{s/α} wieder der Hauptteil der komplexen Wurzel ist. Unter Ausnutzung der Randbedingungen ergibt sich weiter

û_{x}(0, s) = a_{1}(s)\sqrt{\frac{s}{α}} − a_{2}(s)\sqrt{\frac{s}{α}} = 0 , ∀ s .

Es folgt sofort, a_{1}(s) = a_{2}(s) =: a(s) . Die Lösung ergibt sich also zu

û(x, s) = a(s) [exp(\sqrt{\frac{s}{α}} x) + exp(−\sqrt{\frac{s}{α}} x)] = a(s) · 2 cosh(\sqrt{\frac{s}{α}} x) .

Die zweite Randbedingung hilft nun bei der Bestimmung von a(s), denn es ist

û(1, s) = a(s) · 2 cosh(\sqrt{\frac{s}{α}} ) = f̂(s) ,

also

a(s) = \frac{f̂(s)}{2 cosh(\sqrt{s/α})}

und somit

û(x, s) = \frac{cosh(\sqrt{s/α} x)}{cosh(\sqrt{s/α})} · f̂(s) .

2 ALLGEMEINE THEORIE DER PROBLEME 27 Wendet man die inverse Laplace–Transformation auf die Funktion

ĥ(x, s) = (\frac{cosh(\sqrt{s/α} x)}{cosh(\sqrt{s/α})}) .

an, so erhält man mit dem Faltungssatz

u(x, t) = f(t) ∗ h(x, t) .  (2.25)

h(x, t) ergibt sich aus,

ĥ(x, s) = \frac{exp(\sqrt{s/α} x) + exp(−\sqrt{s/α} x)}{exp(\sqrt{s/α}) + exp(−\sqrt{s/α})}
= \frac{exp (\sqrt{s/α} (x − 1)) + exp ( − \sqrt{s/α} (x + 1))}{1 + exp(−2\sqrt{s/α})}
= [exp (\sqrt{s/α} (x − 1)) + exp ( − \sqrt{s/α} (x + 1))]

· ∑_{n=0}^{∞}(−1)^{n} (e^{−2\sqrt{s/α}})^{n}

= ∑_{n=0}^{∞}(−1)^{n} exp (\sqrt{s/α} (x − 2n − 1))

+ ∑_{n=0}^{∞}(−1)^{n} exp ( − \sqrt{s/α} (x + 2n + 1)) .

Unter der Annahme, daß die inverse Laplace–Transformation mit der Summation vertauschbar ist, bekommt man mit Beziehung (2.12)

h(x, t) = ∑_{n=0}^{∞}(−1)^{n} \frac{(2n + 1 − x)}{2\sqrt{πα} t^{3/2}} exp ( − \frac{(2n + 1 − x)^{2}}{4α t} )

+ ∑_{n=0}^{∞}(−1)^{n} \frac{(2n + 1 + x)}{2\sqrt{πα} t^{3/2}} exp ( − \frac{(2n + 1 + x)^{2}}{4α t} )

= ∑_{n=0}^{∞}(−1)^{n} · −2α [\frac{∂K}{∂x} (2n + 1 − x, αt) + \frac{∂K}{∂x} (2n + 1 + x, αt)] .

Für die Lösung u ergibt sich somit aus (2.25)

u(x, t) = ∫_{0}^{t} \{∑_{n=0}^{∞}(−1)^{n} · −2α [\frac{∂K}{∂x} (2n + 1 − x, α(t − τ))  (2.26)

+\frac{∂K}{∂x} (2n + 1 + x, α(t − τ))]\} · f(τ) dτ .

Der nun folgende Satz soll zeigen, daß (2.26) auch unabhängig von den für die Herleitung benutzten Einschränkungen eine Lösung von (2.23) ist.

2.3.1 Satz. Für stückweise stetiges f(t) mit |f(t)| < t^{−β} , 0 < β < 1 , 0 < t < ε ist u(x, t), definiert durch (2.26), eine Lösung der Differentialgleichung (2.23).

Beweis. Für 0 ≤ x < 1 , 0 < t ≤ 1 konvergiert die Summe im Integral gleichmäßig, und das Integral existiert für alle f, die die Vorraussetzungen des Satzes erfüllen. Aus Satz 1.3.1 ist ersichtlich, daß die Summe auch für ∂^{3}K/∂x^{3} und ∂^{2}K/(∂x∂t) gleichmäßig konvergiert. Somit ergibt sich, daß u die Differentialgleichung erfüllt, da ∂K/∂x sie erfüllt. Beim Grenzübergang x → 1 erhält man

lim_{x→1} u(x, t) = lim_{x→1} ∫_{0}^{t} −2α\frac{∂K}{∂x} (1 − x, α(t − τ))f(τ)dτ
= f(t) (nach den Ergebnissen aus Abschnit 2.2),

da für die restlichen Summmanden gilt,

lim_{x→1} [∑_{n=1}^{∞}(−1)^{n} · −2α\frac{∂K}{∂x} (2n + 1 − x, α(t − τ))f(τ)

+ ∑_{n=0}^{∞}(−1)^{n} · −2α\frac{∂K}{∂x} (2n + 1 + x, α(t − τ))f(τ)]

= 0 .

Daß lim_{t→0} u(x, t) = 0 ist, läßt sich ganz in Analogie zu dem Beweis in Abschnit 2.2 zeigen und soll hier nicht wiederholt werden.

Unter Verwendung von

lim_{x→0} \frac{∂}{∂x} (\frac{∂K}{∂x} (2n + 1 + x, α(t − τ)))

= − lim_{x→0} \frac{∂}{∂x} (\frac{∂K}{∂x} (2n + 1 − x, α(t − τ)))

ist auch die letzte Randbedingung u_{x}(0, t) = 0 erfüllt, da sich die Summanden paarweise aufheben. □

Um zu einer Lösung des eigentlichen Problems (2.22) zu gelangen, betrachte man u(x, t) an der Stelle x = 0.

u(0, t) = ∫_{0}^{t} [\frac{1}{\sqrt{πα} (t − τ)^{3/2}} ∑_{n=0}^{∞}(2n + 1)(−1)^{n} exp (\frac{−(2n + 1)^{2}}{4α(t − τ)} )] f(τ)dτ
= g(t) .

Definiert man k_{α} als

k_{α}(t) := \frac{1}{\sqrt{πα} t^{3/2}} ∑_{n=0}^{∞}(2n + 1)(−1)^{n} exp (\frac{−(2n + 1)^{2}}{4αt} ) ,  (2.27)

so erhält man die Volterrasche Integralgleichung

(Af)(t) = ∫_{0}^{t} k_{α}(t − τ)f(τ)dτ = g(t) .  (2.28)

Um einige Eigenschaften des Integralkerns k_{α} zu erarbeiten, betrachte man für 0 < t ≤ 1

∑_{n=0}^{∞} |(2n + 1)(−1)^{n} exp (\frac{−(2n + 1)^{2}}{4αt} )| ≤ ∑_{n=0}^{∞}(2n + 1) exp (\frac{−(2n + 1)^{2}}{4α · 1} ) .

Da die Summme auf der rechten Seite unabhängig von t konvergiert, ist k_{α} beliebig oft differenzierbar, denn

k_{α}^{(ν)}(t) = ∑_{n=0}^{∞} R_{n,ν}(1/t) · exp (\frac{−(2n + 1)^{2}}{4αt} ) ,

wobei R_{n,ν} eine rationale Funktion in 1/t ist. Aus der Darstellung wird auch klar, daß

lim_{t→0^{+}} k_{α}^{(ν)}(t) = 0 für ν = 0, 1, 2 · · · .

k_{α} läßt sich also samt seiner Ableitungen bei t = 0 mit Null fortsetzen. Versucht man nun die Integralgleichung 1.Art mittels Differentiation in eine 2.Art umzuformen, so ergibt sich

(Af)′(t) = k_{α}(0) + ∫_{0}^{t} kα^{′}(t − τ)f(τ)dτ = g′(t) .

Da aber k_{α}^{(ν)} = 0, ist die Umformulierung in eine Intrgralgleichung 2.Art nicht möglich. Setzt man k_{α}(t) = 0 für t ≤ 0 , so läßt sich mit dem erweiterten Kern Gleichung (2.28) schreiben als

(Af)(t) = ∫_{0}^{1} k_{α}(t − τ)f(τ)dτ .  (2.29)

(2.29) ist wieder eine Fredholm–Integralgleichung 1.Art mit glattem Kern, die zu lösen, wie schon in den vorherigen Abschnitten, ein sehr schlecht gestelltes Problem ist.

2.4 Das inverse, eindimensionale Ein–Phasen Stefan Problem

Das Schmelzverhalten eines fiktiven, schmalen, im Wasser liegenden Eisblokkes, der sich auf das Intervall [b, ∞) erstreckt, soll untersucht werden. Nimmt man an, daß der Block überall die konstante Temperatur 0^{o} C besitzt und keine Zirkulationen im Wasser stattfinden, so soll die Temperaturverteilung im Wasser (u(x, t)) und der Verlauf der Schmelzkurve (s(t)) im Intervall 0 ≤ t ≤ T berechnet werden. Eine mathematische Beschreibung des Problems ergibt sich durch: Bestimme die Temperatur im Wasser u(x, t) und die Berührungskurve zwischen Wasser und Eis s(t) unter den Bedingungen,

\begin{aligned}u_{xx} & = & u_{t} , & 0 < x < s(t) , 0 < t ≤ T , \\ u_{x}(0, t) & = & g(t) ≤ 0 , & 0 < t ≤ T , \\ u(x, 0) & = & f(x) ≥ 0 , & 0 ≤ x ≤ b = s(0) , \\ u(s(t), t) & = & 0 , & 0 < t ≤ T . \end{aligned}  (2.30)

Die Vorzeichenbedingungen für g(t) und f(x) besagen dabei, daß bei x = 0 keine Wärme abfließt (u_{x}(0, t) ≤ 0) und daß zur Zeit t = 0 die Temperatur des Wassers nicht unterhalb des Gefrierpunktes liegt und somit ein Schmelz– statt eines Verfestigungsprozesses entsteht.

[figure: x; u = f; t; T; b; (0,0); u_{t} = u_{xx}; \begin{matrix}u = 0 \\ u_{x} = −ṡ \end{matrix}; u_{x} = g]

Da für vorgegebene, vernünftige Funktionen s(t) (mit s(0) = b) und glatte Funktionen g(t) und f(x) (2.30) immer eine Lösung u(x, t) besitzt, muß man, um Eindeutigkeit zu erhalten, eine zusätzliche Bedingung an u und s stellen. Aus der Annahme, daß das Eis überall 0^{o} C habe, folgt, daß die gesamte Wärmeenergie, die an der Schmelzkurve abfließt, in den Schmelzprozeß übergehen muß. Somit gelangt man zu der Zusatzbedingung

ṡ(t) = −u_{x}(s(t), t) ,  0 ≤ t ≤ T .  (2.31)

Das direkte, eindimensionale Ein–Phasen Stefan Problem ergibt sich so aus der Bestimmung von u(x, t) und s(t) unter den Bedingungen (2.30) und (2.31). Das Inverse Stefan Problem (ISP) beschäftigt sich nun mit der Frage, wie man die Randbedingung u_{x}(0, t) = g(t) verändern, damit die Kurve s(t) einen vorgegebenen Verlauf einnimmt. Es ist also g(t) so zu bestimmen, daß ∥s(t) − s̃(t)∥ minimal wird, wobei s̃(t) eine vorgegebene Datenfunktion ist.

Um Eindeutigkeit und Existenz von Lösungen des direkten Problems beweisen zu können, benötigt man folgende Bedingung:

Bedingung A. f und g seien stückweise stetige Funktionen mit höchstens endlich vielen beschränkten Unstetigkeitsstellen. Weiter sei g nicht positiv und f genüge der Ungleichung

0 ≤ f(x) ≤ N(b − x) für eine positive konstante N .

Als Hilfsmittel ist noch ein Lemma vonnöten, daß jedoch hier nicht bewiesen werden soll.

2.4.1 Lemma. Sei Bedingung A erfüllt und s(t) eine auf [0, T] lipschitzstetige Kurve. Ist u(x, t) eine Lösung von (2.30) mit diesem s, so gilt,

u_{x}(s(t), t) existiert und ist stetig für 0 < t ≤ T .

Der Beweis des Lemmas ist zu finden in [4], Seite 14–19.

Unter Verwendung von Lemma 2.4.1 läßt sich nun die Gleichung (2.31) in eine äquivalente Integralgleichung umschreiben. Dazu sei, wie auch im folgenden immer, die Kurve s(t) als lipschitzstetig vorrausgesetzt.

2.4.2 Satz. Unter Bedingung A ist Gleichung (2.31) äquivalent zu

s(t) = b + ∫_{0}^{b} f(x)dx − ∫_{0}^{t} g(τ)dτ − ∫_{0}^{s(t)} u(x, t)dx .  (2.32)

2 ALLGEMEINE THEORIE DER PROBLEME 33

Beweis. Sei (u, s) eine Lösung von (2.30), (2.31), dann ist

0 = ∫_{0}^{t} ∫_{0}^{s(τ)} u_{xx}(x, τ) − u_{τ}(x, τ)dx dτ
= ∫_{0}^{t} \underbrace{u_{x}(s(τ), τ)}_{−ṡ(τ)} dτ − ∫_{0}^{t} \underbrace{u_{x}(0, τ)}_{g(τ)} dτ − ∫_{0}^{t} ∫_{0}^{s(τ)} u_{τ}(x, τ)dx dτ .

Da

\frac{∂}{∂τ} (∫_{0}^{s(τ)} u(x, τ)dx) = u(s(τ), τ) · ṡ(τ) + ∫_{0}^{s(τ)} u_{τ}(x, τ)dx ,

folgt

0 = ∫_{0}^{t} −ṡ(τ)dτ − ∫_{0}^{t} g(τ)dτ

− ∫_{0}^{s(t)} u(x, t)dx + ∫_{0}^{s(0)=b} \underbrace{u(x, 0)}_{f(x)} dx + ∫_{0}^{t} \underbrace{u(s(τ), τ)}_{=0} ·ṡ(τ)dτ

= −s(t) + b − ∫_{0}^{t} g(τ)dτ − ∫_{0}^{s(t)} u(x, t)dx + ∫_{0}^{b} f(x)dx .
⟹  s(t) = b + ∫_{0}^{b} f(x)dx − ∫_{0}^{t} g(τ)dτ − ∫_{0}^{s(t)} u(x, t)dx .

Erfüllt das Paar (u, s) nun die Bedingungen (2.30) und (2.32), so existiert nach Lemma 2.4.1 u_{x}(s(t), t), und man erhält aus Differentiation von (2.32)

ṡ(t) = −g(t) − [\underbrace{u(s(t), t)}_{=0} ·ṡ(t) + ∫_{0}^{s(t)} u_{t}(x, t)dx]
= −u_{x}(0, t) − ∫_{0}^{s(t)} u_{xx}(x, t)dx
= −u_{x}(0, t) − u_{x}(s(t), t) + u_{x}(0, t)
= −u_{x}(s(t), t) .

□

Für den Beweis eines Eindeutigkeitssatzes benötigt man folgende Abwandlung des schwachen Maximumprinzips.

2.4.3 Lemma. Für D_{T} = {(x, t) | 0 < x < s(t), 0 < t ≤ T } sei v eine Lösung von

\begin{aligned}v_{xx} & = & v_{t} , & in D_{T} , \\ v_{x}(0, t) & ≤ & 0 , & 0 < t ≤ T , \\ v(x, 0) & ≥ & 0 , & 0 ≤ x ≤ s(0) = b , \\ v(s(t), t) & ≥ & 0 , & 0 < t ≤ T . \end{aligned}  (2.33)

Dann gilt: v(x, t) ≥ 0 in D_{T} .

Beweis. Sei p > sup_{0≤t≤T} s(t) und

w^{ε}(x, t) := v(x, t) + ε(p − x) ,

so ist

\begin{aligned}w_{xx}^{ε} & = & w_{t}^{ε} , & inD_{T} , \\ w_{x}^{ε}(0, t) & = & v_{x}(0, t) − ε < 0 , & 0 < t ≤ T , \\ w^{ε}(x, 0) & > & 0 , & 0 ≤ x ≤ s(0) = b , \\ w^{ε}(s(t), t) & > & 0 , & 0 < t ≤ T . \end{aligned}  (2.34)

Es ist nun w^{ε} > 0 in D_{T}, da anderenfalls mit Satz 1.2.2 folgt, daß in einem Punkte (0, t_{0}), t_{0} > 0 w^{ε} sein Minimum M ≤ 0 annehmen würde. Aber w_{x}^{ε}(0, t_{0}) ≤ −ε < 0, und somit hätte w^{ε} ein Minimum M′ < M ≤ 0 in D_{T}. Also ist w^{ε} > 0 in D_{T}. Mit ε → 0 folgt die Behauptung. □

Als Korollar kann man eine Aussage über die Monotonie von s(t) tätigen.

2.4.4 Korollar. Jede Lösung von (2.30) ist nicht negativ, und für jede Lösung von ((2.30), (2.31)) gilt: ṡ(t) = −u_{x}(s(t), t) ≥ 0.

Beweis. Die Aussage über das Vorzeichen von u ergibt sich aus Lemma 2.4.3 . Zum Beweis der zweiten Aussage nehmen wir nun an, daß −ṡ(t_{0}) = u_{x}(s(t_{0}), t_{0}) > 0 sei, für ein t_{0} > 0. Dann existiert ein Punkt (ξ_{0}, t_{0}) in D_{T} mit u(ξ_{0}, t_{0}) < u(s(t_{0}), t_{0}) = 0. Dies ist aber ein Widerspruch zur ersten Aussage. □

Als nächstes soll ein Satz über die monotone Abhängigkeit der Lösung von den Daten bewiesen werden.

2.4.5 Satz. Seien (s_{1}, u_{1}) und (s_{2}, u_{2}) Lösungen des Stefan Problems ((2.30), (2.31)) zu den Daten (b_{1}, f_{1}, g_{1}) b.z.w. (b_{2}, f_{2}, g_{2}). Im Falle 0 < b_{1} ≤ b_{2}, f_{1} ≤ f_{2}, g_{1} ≥ g_{2} gilt, s_{1} ≤ s_{2}. Ist zusätzlich b_{1} < b_{2}, so ist s_{1} < s_{2}.

Beweis. Man betrachte zunächst den Fall 0 < b_{1} < b_{2}. Nimmt man an, daß s_{1} < s_{2} nicht gilt, so folgt die Existenz eines t_{0} > 0 mit s_{1}(t_{0}) = s_{2}(t_{0}) und s_{1}(t) < s_{2}(t) ∀ 0 ≤ t < t_{0}. Bei der Untersuchung von w(x, t) = u_{2}(x, t) − u_{1}(x, t) ergibt sich nach Lemma 2.4.3, daß w(x, t) ≥ 0 ist in D_{t_{0}} = {(x, t) | 0 < x < s_{2}(t), 0 < t ≤ t_{0}}. Da w(s_{2}(t), t) = 0 und w(x, 0) ≥ 0 folgt mit dem starken Maximumprinzip, daß sogar w > 0 in D_{t_{0}} ist, denn mit w_{x}(0, t) ≤ 0 ist auch w(0, t) ≥ 0. Nach Satz 1.2.8 (Hopf) folgt somit, da w(s_{2}(t_{0}), t_{0}) = 0,

w_{x}(s_{2}(t_{0}), t_{0}) < 0
⟹ u_{2}_{x}(s_{2}(t_{0}), t_{0}) − u_{1}_{x}(s_{1}(t_{0}), t_{0}) < 0
⟹ ṡ_{1}(t_{0}) − ṡ_{2}(t_{0}) < 0 oder ṡ_{1}(t_{0}) < ṡ_{2}(t_{0}) ,

was einen Widerspruch zur Annahme darstellt. Für den Fall b_{1} = b_{2} > 0 definiert man b_{2}^{δ} = b_{2} + δ. Sei (s_{2}^{δ}, u_{2}^{δ}) die Lösung des Stefan Problems ((2.30), (2.31)) mit den Daten (b_{2}^{δ}, f_{2}^{δ}, g), wobei

f_{2}^{δ} = \{{f_{2} 0 ≤ x ≤ b_{1} = b_{2} \atop 0 b_{2} < x ≤ b_{2} + δ .}

Wie eben folgt s_{1} < s_{2}^{δ}, und s_{2} < s_{2}^{δ}. Unter Verwendung von (2.32) erhält man

s_{2}^{δ}(t) − s_{2}(t) = δ − ∫_{0}^{s_{2}(t)} (u_{2}^{δ}(ξ, t) − u_{2}(ξ, t)) dξ − ∫_{s_{2}(t)}^{s_{2}^{δ}(t)} u_{2}^{δ}(ξ, t)dξ .

Da nach Lemma 2.4.3 u_{2}^{δ} ≥ 0 und u_{2}^{δ} − u_{2} ≥ 0, folgt s_{2}^{δ}(t) ≤ s_{2}(t) + δ. Daraus ergibt sich aber

s_{1}(t) < s_{2}(t) + δ .

Da δ beliebig klein sein darf, folgt s_{1}(t) ≤ s_{2}(t). □

Es ist nun einfach, den Eindeutigkeitssatz zu beweisen.

2.4.6 Satz. Es existiert höchstens eine Lösung des Stefan Problems ((2.30), (2.31)).

Beweis. Seien (s_{1}, u_{1}), (s_{2}, u_{2}) Lösungen des Problems mit den Daten (b, f, g). Aus Satz 2.4.5 ergibt sich, s_{1} ≤ s_{2} und s_{2} ≤ s_{1}, also s_{1} ≡ s_{2}. Mit Lemma 2.4.3 folgt u_{1} ≤ u_{2} und u_{2} ≤ u_{1}, was u_{1} ≡ u_{2} ergibt. □

Bevor die Existenz einer Lösung gezeigt wird, läßt sich für monoton steigendes s(t) unter Bedingung A eine Abschätzung für u_{x}(s(t), t) beweisen, die im weiteren von Nutzen sein wird.

2.4.7 Lemma. Sei u eine Lösung von (2.30) mit monoton steigendem Rand s. Unter Bedingung A existiert eine positive Konstante c = max(N, M), mit M = sup_{0≤t≤T} |g(t)| und N wie in Bedingung A, so daß

0 ≤ ρ^{−1}u(s(t) − ρ, t) ≤ c , ∀ 0 ≤ t ≤ T , 0 < ρ < s(t) .

Beweis. Definiert man die Hilfsfunktion w(x, t_{0}) := c · (s(t_{0}) − x) (w unabhängig von t), so ergibt sich mit Lemma 2.4.3

w − u ≥ 0 ,
⟹ 0 ≤ u(x, t) ≤ c · (s(t_{0}) − x) ,  0 ≤ t ≤ t_{0} , 0 < x ≤ s(t) .

Setzt man nun x = (s(t_{0}) − ρ) , 0 < ρ < s(t_{0}), so folgt

0 ≤ u(s(t_{0}) − ρ, t) ≤ c · ρ .

Da t_{0} beliebig, folgt somit

0 ≤ ρ^{−1}u(s(t) − ρ, t) ≤ c .
□

2.4.8 Korollar. Falls u_{x}(s(t), t) existiert, gilt

−c ≤ u_{x}(s(t), t) ≤ 0 .

Begründung. ρ → 0 in Lemma 2.4.7 □

Mit diesem Ergebnis ist man nun in der Lage, eine Existenzaussage zu beweisen.

2.4.9 Satz. Unter Bedingung A und b > 0 existiert eine Lösung (s, u) des Stefan Problems ((2.30), (2.31)), dessen freier Rand s die Ungleichung 0 ≤ ṡ(t) ≤ c, c wie in Lemma 2.4.7, erfüllt.

Beweis. Zum Beweis konstruiert man eine Familie von Lösungen (s^{τ}, u^{τ}) , 0 < τ ≤ τ_{0}, von (2.30), mit dem Rand s^{τ}, gegeben durch

s^{τ}(t) = b , 0 ≤ t ≤ τ

und

ṡ^{τ}(t) = −u_{x}^{τ}(s^{τ}(t − τ), t − τ ) , t > τ .

Das zugehörige f habe dabei die Form

f^{τ}(x) = \{{f(x) 0 ≤ x < b − τ \atop 0 b − τ ≤ x ≤ b ,}.

Das so konstruierte s^{τ}(t) erfüllt nach Korollar 2.4.8 die Bedingung

0 ≤ ṡ^{τ}(t) ≤ c ,  0 ≤ t ≤ T .  (2.35)

Aus (2.35) folgt, daß s^{τ}(t) eine gleichmäßig beschränkte, gleichgradig stetige Familie von Funktionen ist. Nach dem Satz von Arzela–Ascoli existiert somit eine Teilfolge s^{τ_{ν}}, die gleichmäßig gegen eine Funktion s konvergiert. Sei o.B.d.A. s^{τ} = s^{τ_{ν}}. Bezeichnet man mit u die Lösung des Problems (2.30) bezüglich diesem s, und definiert w := u^{τ} − u , 0 ≤ x ≤ α^{τ}(t), mit α^{τ}(t) = min(s^{τ}(t), s(t)), so folgt

|w(α^{τ}(t), t)| = \{ \begin{matrix}u^{τ}(s(t), t) , & s(t) < s^{τ}(t) \\ u(s^{τ}(t), t) & sonst \end{matrix}  (2.36)
≤ c · |s^{τ}(t) − s(t)| ,

unter Verwendung von ρ = |s^{τ}(t) − s(t)| in Lemma 2.4.7. Weiter gilt, w_{x}(0, t) = 0, w(x, 0) = 0 für 0 ≤ x ≤ b − τ und |w(x, 0)| ≤ c · τ , b − τ ≤ x ≤ b. Wendet man nun Lemma 2.4.3 auf

z^{τ}(x, t) := c · ∥s^{τ} − s∥_{[0,T]} + c · τ ± w(x, t)

an, so folgt für 0 ≤ x ≤ α^{τ}(t), 0 ≤ t ≤ T

|u^{τ}(x, t) − u(x, t)| ≤ c · (∥s^{τ} − s∥_{[0,T]} + τ) ,  (2.37)

wobei ∥s^{τ} − s∥_{[0,T]} = sup_{0≤t≤T} |s^{τ}(t) − s(t)|.

Betrachtet man nun w(x, t) im Bereich α^{τ}(t) ≤ x ≤ β^{τ}(t), 0 ≤ t ≤ T , β^{τ}(t) = max(s^{τ}(t), s(t)) , so folgt nach (2.36)

|w(α^{τ}(t), t)| ≤ c · |s^{τ}(t) − s(t)| .

Da weiter w(β^{τ}(t), t) = 0 und α^{τ}(0) = β^{τ}(0) = b folgt mit dem Maximumprinzip angewandt auf z^{τ}(x, t) = c · ∥s^{τ} − s∥_{[0,T]} ± w(x, t)

∥u^{τ}(x, t) − u(x, t)∥ ≤ c · ∥s^{τ} − s∥_{[0,T]} ,  (2.38)

für α^{τ}(t) ≤ x ≤ β^{τ}(t), 0 ≤ t ≤ T. Aus der gleichmäßigen Konvergenz der Folge s^{τ} gegen s, folgt somit aus (2.37) und (2.38) die Konvergenz der zugehörigen u^{τ} gegen u, für τ → 0 . Da ṡ^{τ}(t + τ) = −u_{x}^{τ}(s^{τ}(t), t), t > 0 folgt analog zum Beweis von Satz 2.4.2

s^{τ}(t + τ) = b + ∫_{0}^{b−τ} f(ξ)dξ − ∫_{0}^{t} g(η)dη − ∫_{0}^{s^{τ} (t)} u^{τ}(ξ, t)dξ .

Für τ → 0 ergibt sich so, daß die oben konstruierte Lösung (s, u) die Bedingung (2.32) erfüllt, und somit nach Satz 2.4.2 (s, u) eine Lösung des Stefan Problems ((2.30), (2.31)) ist. □

Kehrt man nun zum ISP zurück und stellt sich die Aufgabe, eine Funktion g(t) zu finden, so daß die vorgegebene Kurve s(t) genau erreicht wird, so erhält man folgendes Problem. Finde g(t) = u_{x}(0, t) unter den Bedingungen

\begin{aligned}u_{xx} & = & u_{t} , & in D_{T} , \\ u(s(t), t) & = & 0 , & 0 < t ≤ T , \\ u_{x}(s(t), t) & = & ṡ(t) , & 0 < t ≤ T , \\ u(x, 0) & = & f(x) , & 0 < x < s(0) = b . \end{aligned}  (2.39)

(2.39) ist ein nichtcharakteristisches Cauchy Problem, da die Cauchy Daten auf einer Kurve s(t) gegeben sind, die keine Charakteristik ist. Das ISP ist somit ebenso schlecht gestellt wie das nichtcharakteristische Cauchy Problem in Abschnitt 2.3 .

2.5 Parameteridentifikation

In vielen partiellen Differentialgleichungen tauchen neben der zu bestimmenden Funktion noch orts– oder zeitabhängige Parameter auf. Diese beschreiben oft physikalische Eigenschaften des zu untersuchenden Mediums. Um diese Parameter zu bestimmen, ist es in einigen Fällen einfacher, die Lösung der Differentialgleichung an einigen Punkten zu messen, als den Parameter direkt zu bestimmen. Probleme dieser Art werden als ”inverse Probleme“ bezeichnet. Betrachtet man nun die in Abschnitt 1.1 gegebene Gleichung (1.2) mit dem ortsabhängigem Parameter a(x) statt der Konstanten k so ergibt sich

q = −a(x)u_{x} .  (2.40)

Kombiniert man (2.40) mit einem ”Erhaltungsgesetz“ der Form

u_{t} + q_{x} = f ,  (2.41)

so erhält man die inhomogene partielle Differentialgleichung

u_{t} = (a(x)u_{x})_{x} + f(x, t) .  (2.42)

Aufgabe ist es nun, a(x) aus Meßwerten der Funktion u(x, t) im Bereich 0 ≤ x ≤ 1, 0 ≤ t ≤ T zu bestimmen.

[figure: x; t; T; 1; (0,0); u_{t} = (a(x)u_{x})_{x} + f(x, t)]

Als erstes stellt sich dabei die Frage nach der Eindeutigkeit der Lösung a(x). Diese wird durch folgenden Satz beantwortet.

2.5.1 Satz. Sei u(x, t) definiert und stetig auf [0, 1] × [0, T]. Weiter existiere ein t_{0} ∈ [0, T], so daß u_{x}(x, t_{0}) ≠ 0 f.ü. . Die Menge A sei gegeben durch

A = {a(x) | a(x) ist stetig und 0 < a_{1} ≤ a(x) ≤ a_{2}} .

Dann gilt :

1. Falls zu jedem t ∈ [0, T] ein x ∈ [0, 1] existiert, so daß u_{x}(x, t) = 0 (z.B. falls u(0, t) = u(1, t)), dann gibt es höchstens ein a(x) ∈ A, mit dem u eine Lösung von (2.42) ist.

2. Existiert ein t_{1} ∈ [0, T] mit u_{x}(x, t_{1}) ≠ 0, ∀ x ∈ [0, 1], und sei a gegeben an einem Punkt x̃, d.h. a(x̃) = ã, so existiert höchstens ein a(x) ∈ A, mit dem u eine Lösung von (2.42) darstellt.

Beweis. Annahme: Es gebe zwei Lösungen a, a^{∗} ∈ A. Es folgt,

u_{t}(x, t) = (a(x)u_{x}(x, t))_{x} + f(x, t) = (a^{∗}(x)u_{x}(x, t))_{x} + f(x, t)

⟹  ([a(x) − a^{∗}(x)]u_{x}(x, t))_{x} = 0
⟹ [a(x) − a^{∗}(x)]u_{x}(x, t) = h(t) ,  (2.43)

mit einer von x unabhängigen Funktion h : [0, t] → R I . Im 1. Fall gilt nun h ≡ 0, da zu jedem t ein x existiert, so daß u_{x}(x, t) und somit die linke Seite von (2.43) null ist. Im 2. Fall ist ebenfalls h ≡ 0, da unabhängig von t, a(x̃) = a^{∗}(x̃) = ã. Insgesamt folgt in beiden Fällen

[a(x) − a^{∗}(x)]u_{x}(x, t) = 0 .

Da nun aber u_{x}(x, t_{0}) ≠ 0 f.ü. , folgt a(x) − a^{∗}(x) = 0 a.e. und somit aus der Stetigkeit von a und a^{∗} die Behauptung. □

Der erste Ansatz zur Bestimmung einer Lösung von (2.42) ist die einfache Integration der Differentialgleichung. Man erhält hierbei für u_{x}(x, t) ≠ 0

a(x) = \frac{∫ u_{t}(x, t) − f(x, t) dx + c(t)}{u_{x}(x, t)} ,  (2.44)

wobei c(t) durch die Vorgaben in Satz 2.5.1 eindeutig bestimmt ist. Bei der numerischen Berechnung von a(x) ergibt sich nun als erstes Problem, daß a(x) nicht direkt von u, sondern von u_{x} und u_{t} abhängt. Die Bestimmung von u_{t} und u_{x} aus den Daten über u ist an sich schon ein schlecht gestelltes Problem, welches aber nicht so schwer wiegt, da es Verfahren gibt, die akzeptable Lösungen liefern. Das zweite Problem entsteht bei der Bestimmung einer Näherung ã(x) von a(x) mit den eben erwähnten fehlerhaften Funktionen u_{x} und u_{t}. Für ã erhält man,

ã(x) = \frac{∫ ũ_{t}(x, t) − f(x, t) dx + c_{1}(t)}{ũ_{x}(x, t)} .  (2.45)

Besitzt nun u_{x} bei x_{0} eine Nullstelle, so müssen aufgrund der Beschränktheit von a(x) sowohl Nenner als auch Zähler in Gleichung (2.44) eine Nullstelle besitzen. Da dies für Nenner und Zähler in Gleichung (2.45) nicht der Fall sein muß, kann ã(x) in einer Umgebung von x_{0} Werte von −∞ bis ∞ annehmen. Die Bestimmung einer Lösung des Problems ist also als äußerst schlecht gestellt zu bezeichnen, da sich im Nullstellenbereich von u_{x} das Verhalten der Näherung ã(x) nicht erfassen läßt.

3 Regularisierungsmethoden

3.1 Allgemeine Theorie

Um schlecht gestellte Probleme der Form Af = g, wie sie in Abschnitt 2 beschrieben wurden, zu lösen, ist es nötig, sowohl die Eindeutigkeit und Existenz einer Lösung zu garantieren als auch die stetige Abhängigkeit dieser regularisierten Lösung von den Daten. Die ersten beiden Forderungen lassen sich dadurch erfüllen, daß man eine sogenannte verallgemeinerte Lösung betrachtet.

3.1.1 Definition. Seien X, Y Hilberträume, A : X → Y ein stetiger, linearer Operator. Man bezeichnet mit A^{+}g die verallgemeinerte Lösung von Af = g, falls gilt:

f = A^{+}g ⟺ \{ \begin{matrix}∥Af − g∥_{Y} ≤ ∥Au − g∥_{Y} ∀u ∈ X , \\ ∥f∥_{X} ≤ ∥v∥_{X} ∀v ∈ X mit ∥Af − g∥_{Y} = ∥Av − g∥_{Y} . \end{matrix}

A^{+} wird auch als Moore–Penrose–Inverse bezeichnet.

Man zeigt leicht, daß diese Definition äquivalent ist zu

3.1.2 Satz. Für g ∈ R(A) + R(A)^{⊥} gilt:

f = A^{+}g ⟺ \{{A^{∗}Af = A^{∗}g , \atop f ∈ \overline{R(A^{∗})} .}

Beweis. Siehe [16], Seite I – 2.5 .

Das Problem der stetigen Abhängigkeit des Ergebnisses von den Daten ist damit allerdings noch nicht gelöst,da die verallgemeinerte Inverse A^{+} im allgemeinen unstetig und A^{+}g somit nicht direkt berechenbar ist. Um auch dieses Hindernis zu überwinden, ersetzt man A^{+} durch einen stetigen Operator R_{γ}, dessen Lösung f̃ = R_{γ}g nicht zu weit von A^{+}g abweicht.

3.1.3 Definition. Sei A : X → Y ein linearer, stetiger Operator zwischen Hilberträumen. Für γ > 0 sei R_{γ} : X → Y linear und stetig. R_{γ} heißt Regularisierung von A, falls gilt:

R_{γ}g ⟶^{γ→0} A^{+}g ∀g ∈ R(A) + R(A)^{⊥} .

γ wird als Regularisierungsparameter des Operators bezeichnet.

Um die Wirkungsweise der einzelnen Regularisierungen besser untersuchen zu können, benötigt man den Begriff der Singulärwertzerlegung.

3.1.4 Definition. Der Operator A besitzt eine Singulärwertzerlegung (SVD) genau dann, wenn gilt:

Af = ∑_{k=1}^{∞} σ_{k}(f, f_{k})_{X}g_{k} ,

wobei {f_{k}} ein Orthonormalsystem (ONS) in X und {g_{k}} ein ONS in Y ist. Ferner muß für die Singulärwerte σ_{k} von A gelten:

σ_{k} > 0 ∀k (insbesondere also σ_{k} ∈ R I ).

Aus (Af, g)_{Y} = ∑_{k=1}^{∞} σ_{k}(f, f_{k})_{X}(g, g_{k})_{Y} = (f, ∑_{k=1}^{∞} σ_{k}(g, g_{k})_{Y} f_{k})_{X} folgt

A^{∗}g = ∑_{k=1}^{∞} σ_{k}(g, g_{k})_{Y} f_{k} .

Insbesondere gelten die Relationen:

Af_{k} = σ_{k}g_{k}  sowie A^{∗}g_{k} = σ_{k}f_{k}  ∀k ∈ N I .

Durch leichtes Nachrechnen ergibt sich, daß f_{k} bzw. g_{k} die Eigenfunktionen von A^{∗}A bzw. AA^{∗} zum Eigenwert σ_{k}^{2} sind. Für die verallgemeinerte Inverse A^{+}g folgt

3.1.5 Satz. Sei g ∈ R(A) + R(A)^{⊥}, dann gilt:

A^{+}g = ∑_{k=1}^{∞} σ_{k}^{−1}(g, g_{k})_{Y} f_{k} .

Beweis. Siehe [16], Seite I – 2.6 .

Betrachtet man nun die Norm von A^{+}g so ergibt sich

∥A^{+}g∥^{2} = ∑_{k=1}^{∞} |σ_{k}^{−1}(g, g_{k})_{Y} |^{2} ≤ sup_{k}(1/σ_{k}^{2}) ∑_{k=1}^{∞} |(g, g_{k})_{Y} |^{2} ≤ sup_{k}(1/σ_{k}^{2}) ∥g∥^{2} .

Es ist also A^{+} stetig ⇔ die σ_{k}^{−1} sind beschränkt. Ordnet man die Singulärwerte von A der Größe nach an, also σ_{k} ≥ σ_{k+1} ∀k und gilt σ_{k} ⟶^{k→∞} 0, so folgt 1/σ_{k} ⟶^{k→∞} ∞ und damit die Unstetigkeit von A^{+}.

Eine erste einfache Regularisierungsmethode ist in diesem Fall dadurch gegeben, daß man die kleinen Singulärwerte von A bei der Berechnung von A^{+}g nicht beachtet. Man setzt also

R_{γ}g = ∑_{σ_{k}>γ} σ_{k}^{−1}(g, g_{k})f_{k} .

Offensichtlich gilt R_{γ}g ⟶ A^{+}g für γ → 0. Die Stetigkeit von R_{γ} folgt aus

∥R_{γ}g∥^{2} = ∑_{σ_{k}>γ} |σ_{k}^{−1}(g, g_{k})_{Y} |^{2} ≤ 1/γ^{2} ∑_{σ_{k}>γ}(g, g_{k})_{Y}^{2} ≤ 1/γ^{2}∥g∥^{2} ,

also ∥R_{γ}∥ ≤ 1/γ < ∞. Diese Methode wird auch die abgeschnittene SVD genannt.

3.2 Die Methode des ”quasi–reversibilité“

Um die Wirkungsweise der Methode auf das Problem (2.2) aus Abschnitt 2.1

\begin{aligned}u_{xx}(x, t) & = & u_{t}(x, t) , & 0 < x < 1 , 0 < t ≤ T , \\ u(0, t) & = & u(1, t) = 0 , & 0 < t ≤ T , \\ u(x, T) & = & g(x) , & 0 < x < 1 , T > 0 , \\ u(x, 0) & = & f(x) = ? , & 0 < x < 1 . \end{aligned}

zu untersuchen, ist es sinnvoll, eine SVD für den zugehörigen Operator herzuleiten. Dazu betrachte man den Ansatz:

u(x, t) = ∑_{k=1}^{∞} a_{k}(t) sin(kπx) .

Für die partiellen Ableitungen folgt dann

u_{xx}(x, t) = ∑_{k=1}^{∞} −a_{k}(t)k^{2}π^{2} sin(kπx)
u_{t}(x, t) = ∑_{k=1}^{∞} a′_{k}(t) sin(kπx) ,

und nach Anwendung der Differentialgleichung ergibt sich für die Koeffizienten

a′_{k}(t) = −a_{k}(t)k^{2}π^{2}

und somit

a_{k}(t) = c_{k}e^{−k^{2}π^{2}t} ,

mit Konstanten c_{k} für k > 0 . Um c_{k} zu bestimmen, betrachte man nun die Funktion u(x, t) an der Stelle t = 0 . Es ergibt sich aus

u(x, 0) = f(x) = ∑_{k=1}^{∞} c_{k} sin(kπx) ,

daß die Werte für c_{k} die Fouriersinuskoeffizienten von f sind, also

c_{k} = 2 ∫_{0}^{1} f(ξ) sin(kπξ)dξ .

Insgesamt erhält man für u(x, t):

u(x, t) = ∑_{k=1}^{∞} e^{−k^{2}π^{2}t} 2 ∫_{0}^{1} f(ξ) sin(kπξ)dξ sin(kπx)

und für t = T

\begin{aligned}u(x, T) = Af(x) & = & ∑_{k=1}^{∞} e^{−k^{2}π^{2}T} 2 ∫_{0}^{1} f(ξ) sin(kπξ)dξ sin(kπx) \\ & = & g(x) . \end{aligned} (3.1)

Da aufgrund des Maximumprinzips die Lösung des Problems eindeutig ist, ist die Lösung (3.1) die SVD des Operators A aus Abschnitt 2.1. Betrachtet man nun die verallgemeinerte Inverse, so erhält man

A^{+}g(x) = ∑_{k=1}^{∞} e^{k^{2}π^{2}T} 2 ∫_{0}^{1} g(ξ) sin(kπξ)dξ sin(kπx) .  (3.2)

Es ist sofort ersichtlich, daß die Singulärwerte von A^{+}, 1/σ_{k} = e^{k^{2}π^{2}T} für beliebig große Werte von k unbeschränkt sind. Die verallgemeinerte Inverse läßt sich also nicht direkt berechnen.

Um eine Regularisierung durchzuführen, betrachte man anstatt des ursprünglichen Problems die Differentialgleichung

\begin{aligned}v_{t}(x, t) − γv_{xxt}(x, t) & = & v_{xx}(x, t) , & 0 < x < 1 , 0 < t ≤ T , \\ v(0, t) = v(1, t) & = & 0 , & 0 < t ≤ T , \\ v(x, T) & = & g(x) , & 0 < x < 1 , \\ v(x, 0) & = & f(x) = ? , & 0 < x < 1 . \end{aligned} (3.3)

Mit dem gleichen Ansatz wie oben erhält man

v_{t}(x, t) = ∑_{k=1}^{∞} a′_{k}(t) sin(kπx) ,
v_{xx}(x, t) = ∑_{k=1}^{∞} −a_{k}(t)k^{2}π^{2} sin(kπx) ,
v_{xxt}(x, t) = ∑_{k=1}^{∞} −a′_{k}(t)k^{2}π^{2} sin(kπx) .

Es folgt dann

a′_{k}(t)(1 + γk^{2}π^{2}) = −a_{k}(t)k^{2}π^{2}

oder

a_{k}(t) = c_{k} exp (−\frac{k^{2}π^{2}}{1 + γk^{2}π^{2}} t) .

Aus der Identität

v(x, 0) = f(x) = ∑_{k=1}^{∞} c_{k} sin(kπx)

folgt

c_{k} = 2 ∫_{0}^{1} f(ξ) sin(kπξ)dξ ,

und man erhält das Ergebnis:

v(x, t) = ∑_{k=1}^{∞} exp (−\frac{k^{2}π^{2}}{1 + γk^{2}π^{2}} t) 2 ∫_{0}^{1} f(ξ) sin(kπξ)dξ sin(kπx) .

Für t = T ist dann

A_{γ}f(x) := ∑_{k=1}^{∞} exp (−\frac{k^{2}π^{2}}{1 + γk^{2}π^{2}} T ) 2 ∫_{0}^{1} f(ξ) sin(kπξ)dξ sin(kπx)
= g(x) .

Betrachtet man auch hier die verallgemeinerte Inverse, so ist:

R_{γ}g(x) = ∑_{k=1}^{∞} exp (\frac{k^{2}π^{2}}{1 + γk^{2}π^{2}} T ) 2 ∫_{0}^{1} g(ξ) sin(kπξ)dξ sin(kπx) . (3.4)

Es ist sofort ersichtlich, daß ∥R_{γ}∥ ≤ e^{(1/γ)T} ist. Um nun zu zeigen, daß R_{γ}g eine Regularisierung darstellt, ist noch die Konvergenz von R_{γ}g → A^{+}g zu beweisen. Hierzu betrachte man:

∥(R_{γ} − A^{+})g∥^{2} = ∥ ∑_{k=1}^{∞} [exp (\frac{k^{2}π^{2}}{1 + γk^{2}π^{2}} T ) − exp (k^{2}π^{2}T )] (g, g_{k})f_{k}∥^{2}
= ∑_{k=1}^{∞} \underbrace{| [exp (\frac{k^{2}π^{2}}{1 + γk^{2}π^{2}} T ) − exp (k^{2}π^{2}T )] (g, g_{k}) |^{2}}_{:=c_{k}(γ)} ,

wobei f_{k}(x) = g_{k}(x) = sin(kπx) die Singulärfunktionen von A sind. Da mit γ → 0 auch c_{k}(γ) → 0 geht, ist für die Konvergenz von R_{γ}g gegen A^{+}g nur noch die gleichmäßige Konvergenz der Summe über die c_{k}(γ) zu zeigen. Dies geschieht mit dem Weierstraßschen Majorantenkriterium. Es ist

| exp (\frac{k^{2}π^{2}}{1 + γk^{2}π^{2}} T ) − exp (k^{2}π^{2}T ) |
= exp (k^{2}π^{2}T ) | exp ((\frac{k^{2}π^{2}}{1 + γk^{2}π^{2}} − k^{2}π^{2}) T ) − 1|
= exp (k^{2}π^{2}T ) \underbrace{| exp (\frac{−γk^{4}π^{4}}{1 + γk^{2}π^{2}} T ) − 1|}_{≤1}

≤ exp (k^{2}π^{2}T ) = 1/σ_{k} ,

und somit

sup_{γ∈ R I ^{+}} |c_{k}(γ)| ≤ |1/σ_{k} (g, g_{k})|^{2} .

Daraus ergibt sich

∑_{k=1}^{∞} c_{k}(γ) ≤ ∑_{k=1}^{∞} |1/σ_{k} (g, g_{k})|^{2} = ∥A^{+}g∥^{2} < ∞ ,

und man erhält abschließend

lim_{γ→0} ∥(R_{γ} − A^{+})g∥^{2} = 0 .

Vergleicht man nun die Regularisierung (3.4) mit der verallgemeinerten Inversen (3.2), so ergibt sich, daß die Singulärwerte von A^{+}

ρ_{k} = exp (k^{2}π^{2}T )

ersetzt wurden durch Werte der Form

ρ̃_{k} = exp (\frac{k^{2}π^{2}}{1 + γk^{2}π^{2}} T ) .

Das heißt, für relativ kleine Singulärwerte, bei denen k^{2} ∼ 1/γ gilt, ist ρ_{k} ∼ ρ̃_{k}. Ist jedoch k^{2} ≫ 1/γ oder auch ρ_{k} ≫ 1 , so wird ρ_{k} durch einen kleineren Wert ersetzt, für den gilt: ρ̃_{k} < e^{1/γ}. Die Regularisierung läßt sich also auf eine Filterung der Singulärwerte von A^{+} zurückführen, bei der relativ kleine Singulärwerte kaum verändert und große Singulärwerte stark gefiltert werden.

3.3 Ein modifiziertes Cauchy Problem

Da es nicht immer möglich ist, die SVD eines Operators explizit anzugeben, ist es oft nötig, um die Wirkungsweise einer Regularisierung zu untersuchen, eine andere Darstellung der Lösung zu betrachten.

So wurde zum Beispiel im Abschnitt 2.2 gezeigt, daß sich die Lösung des Problems (2.5) ergibt aus dem Lösen der Integralgleichung

Af(t) = \frac{1}{2\sqrt{π}} ∫_{0}^{t} \frac{f(τ)}{(t − τ)^{3/2}} exp (\frac{−1}{4(t − τ)})dτ
= g(t) .

Wendet man auf beiden Seiten der Gleichung eine Laplace–Transformation an, so erhält man mit dem Faltungssatz und der Identität (2.12) für x = 1:

e^{−\sqrt{s}} · f̃(s) = g̃(s) ,  (3.5)

wobei \sqrt{s} der Hauptzweig der komplexen Wurzel und f̃(s) die Laplace– Transformierte der Funktion f(t) ist.

Bezeichnet man mit f̂(ξ) die Fourier–Transformierte von f(t), so gilt für Funktionen, die auf der negativen t–Achse verschwinden:

f̃(iξ) = \sqrt{2π} · f̂(ξ) .

Setzt man s = iξ, so ergibt sich damit in (3.5)

e^{−\sqrt{iξ}} f̂(ξ) = ĝ(ξ) oder f̂(ξ) = e^{\sqrt{iξ}} ĝ(ξ) .

Ist nun f̂(ξ) ∈ L_{1}( R I ), so erhält man mit der inversen Fourier–Transformation

f(t) = \frac{1}{\sqrt{2π}} ∫_{−∞}^{+∞} e^{iξt}e^{\sqrt{iξ}} ĝ(ξ)dξ ,  (3.6)

wobei \sqrt{iξ} der Hauptzweig der komplexen Wurzel, also

\sqrt{iξ} = \{{(1 + i)\sqrt{\frac{|ξ|}{2}}  ξ ≥ 0 , \atop (1 − i)\sqrt{\frac{|ξ|}{2}}  ξ < 0}

ist.

Man hat somit eine Darstellung der Lösung f(t) in Abhängigkeit von der Fourier–Transformierten der Datenfunktion ĝ(ξ) gefunden.

Um eine Regularisierung durchzuführen, betrachte man die Aufgabe

\begin{aligned}v_{t}(x, t) + γ^{2}v_{ttxx}(x, t) & = & v_{xx}(x, t) , & 0 < x < ∞ , 0 < t < ∞ , \\ v(x, 0) & = & 0 , & 0 ≤ x < ∞ , \\ v_{t}(x, 0) & = & 0 , & 0 ≤ x < ∞ , \\ v(1, t) & = & g(t) , & 0 ≤ t < ∞ , \\ v(0, t) & = & f_{γ}(t) = ? , & 0 ≤ t < ∞ . \end{aligned} (3.7)

mit der zusätzlichen Bedingung, daß v(x, t) für x → ∞ beschränkt ist.

In Analogie zu Abschnitt 2.2 erhält man mit Hilfe der Laplace–Transformation angewandt auf (3.7) die gewöhnliche Differentialgleichung

ṽ_{xx} = sṽ + s^{2}γ^{2}ṽ_{xx}

oder auch

ṽ_{xx} = \frac{s}{1 − γ^{2}s^{2}} ṽ

und mit dem gleichem Verfahren wie dort die Lösung

ṽ(x, s) = exp (−x\sqrt{\frac{s}{1 − γ^{2}s^{2}}} ) f̃_{γ}(s) .

Für x = 1 ergibt sich

g̃(s) = exp (−\sqrt{\frac{s}{1 − γ^{2}s^{2}}} ) f̃_{γ}(s) .  (3.8)

Geht man auch in (3.8) zur Fourier–Transformierten über, so erhält man mit s = iξ:

f̂_{γ}(ξ) = exp (\sqrt{\frac{iξ}{1 + ξ^{2}γ^{2}}} ) · ĝ(ξ) .  (3.9)

Mit einer abschließenden inversen Fourier–Transformation ist dann

\begin{aligned}f_{γ}(t) & = & \frac{1}{\sqrt{2π}} ∫_{−∞}^{+∞} exp(iξt) exp (\sqrt{\frac{iξ}{1 + ξ^{2}γ^{2}}} ) · ĝ(ξ)dξ \\ & = & R_{γ}g . \end{aligned}  (3.10)

3 REGULARISIERUNGSMETHODEN 53 Es folgt nun sofort mit (3.9) und der L_{2}–Norm:

∥R_{γ}g∥^{2} = ∥f̂_{γ}∥^{2} ≤ sup_{ξ∈ R I} | exp (\sqrt{\frac{iξ}{1 + ξ^{2}γ^{2}}} ) |^{2} · ∥ĝ∥^{2}
≤ sup_{ξ∈ R I} \{ exp (\sqrt{\frac{|ξ|}{2}} · \sqrt{\frac{1}{1 + ξ^{2}γ^{2}}} ) \}^{2} · ∥ĝ∥^{2}
= sup_{ξ∈ R I ^{+}} \{ \underbrace{exp (\sqrt{\frac{ξ}{2 + 2ξ^{2}γ^{2}}} )}_{max. bei ξ=1/γ} \}^{2} · ∥ĝ∥^{2}
≤ (exp (\sqrt{\frac{1}{4γ}} ))^{2} · ∥ĝ∥^{2}
⇒  ∥R_{γ}∥ ≤ exp (\sqrt{\frac{1}{4γ}} ) < ∞ .

Um die Konvergenz von R_{γ}g ⟶^{γ→0} f zu zeigen, betrachte man

| exp(iξt) exp (\sqrt{\frac{iξ}{1 + γ^{2}ξ^{2}}} ) · ĝ(ξ)| = exp (\sqrt{\frac{1}{1 + γ^{2}ξ^{2}}} ) · | exp (\sqrt{iξ} ) · ĝ(ξ)|
≤ e^{1} · | exp (\sqrt{iξ} ) · ĝ(ξ) |
= e^{1} · |f̂(ξ) | .

Da f̂ ∈ L_{1}( R I ) ist, folgt nun mit dem Satz von der majorisierten Konvergenz:

lim_{γ→0} R_{γ}g(t) = lim_{γ→0} \frac{1}{\sqrt{2π}} ∫_{−∞}^{+∞} exp(iξt) exp (\sqrt{\frac{iξ}{1 + γ^{2}ξ^{2}}} )ĝ(ξ)dξ
= \frac{1}{\sqrt{2π}} ∫_{−∞}^{+∞} exp(iξt) lim_{γ→0} (exp (\sqrt{\frac{iξ}{1 + γ^{2}ξ^{2}}} )) ĝ(ξ)dξ
= \frac{1}{\sqrt{2π}} ∫_{−∞}^{+∞} exp(iξt) exp (\sqrt{iξ} )ĝ(ξ)dξ
= f(t) .

Vergleicht man die Lösung (3.6) mit der Regularisierung (3.10), so ist sofort ersichtlich, daß lediglich der Faktor exp (\sqrt{iξ} ) ersetzt wurde durch exp (\sqrt{\frac{iξ}{1+γ^{2}ξ^{2}}} ).

Da | exp (\sqrt{iξ})| = exp (\sqrt{|ξ|/2} ) für |ξ| → ∞ beliebig groß wird, werden die hochfrequenten Teile der Datenfunktion durch diesen Faktor verstärkt. Dies geschieht bei der Regularisierung nicht, da hier gilt:

| exp (\sqrt{\frac{iξ}{1 + γ^{2}ξ^{2}}} ) | = exp (\sqrt{\frac{|ξ|}{2(1 + γ^{2}ξ^{2})}} ) ⟶^{|ξ|→∞} 1 .

Für kleine Werte von ξ hingegen ist exp (\sqrt{iξ}) ∼ exp (\sqrt{\frac{iξ}{1 + γ^{2}ξ^{2}}} ).

Die Regularisierung beläßt also die niederfrequenten Teile fast unverändert, während ein Aufblähen der hochfrequenten Teile verhindert wird. Es handelt sich bei der Regularisierung also um eine Art Filter für die hochfrequenten Teile der Datenfunktion.

3.4 Die Tikhonov–Regularisierung

Es sei der Operator A : X → Y gegeben, und gesucht wird ein f ∈ X, so daß die Gleichung Af = g erfüllt wird. Ist anstelle der Funktion g nur eine Näherung g_{δ} bekannt, so ist es nötig (siehe Def. 3.1.1), um eine Lösung zu bekommen, die Aufgabe zu ersetzen durch

min_{f∈X} (∥Af − g_{δ}∥_{Y}^{2} ) .  (3.11)

Ersetzt man nun (3.11) durch

min_{f∈X} (∥Af − g_{δ}∥_{Y}^{2} + γ^{2}∥f∥_{X}^{2}) ,  (3.12)

so erhält man die sogenannte Tikhonov–Regularisierung. Um zu zeigen, daß es sich tatsächlich um eine Regularisierung handelt, betrachte man die Gleichung

\binom{A}{γI}f = \binom{g_{δ}}{0} ,

wobei \binom{A}{γI} ein Operator von X nach X × Y darstellt. Es ergibt sich

∥\binom{A}{γI}f − \binom{g_{δ}}{0}∥_{X×Y}^{2} = ∥Af − g_{δ}∥_{Y}^{2} + γ^{2}∥f∥_{X}^{2} .

Die Lösung der Minimierungsaufgabe (3.12) ist nach Satz 3.1.2 also charakterisiert durch

\binom{A}{γI}^{∗}\binom{A}{γI}f = \binom{A}{γI}^{∗}\binom{g_{δ}}{0}

oder, wie sich leicht nachrechnen läßt, durch

(A^{∗}A + γ^{2}I)f = A^{∗}g_{δ} .  (3.13)

(3.13) ist die regularisierte Normalgleichung. Um zu zeigen, daß f als Lösung der Gleichung (3.13) eindeutig bestimmt ist und stetig von den Daten g_{δ} abhängt, benutzt man

3.4.1 Lemma. B : X → X sei linear und stetig, ferner sei B selbstadjungiert (B = B^{∗}) und positiv definit, d.h. ∃γ > 0 s.d.:(Bf, f) ≥ γ(f, f) ∀f ∈ X, dann besitzt B eine auf ganz X definierte Inverse B^{−1}, und es gilt ∥B^{−1}∥ ≤ 1/γ.

Beweis. Siehe [16], Seite I – 4.4 .

Setzt man B = (A^{∗}A + γ^{2}I), so ist sofort ersichtlich, daß B selbstadjungiert ist. Die positive Definitheit ergibt sich aus

(A^{∗}Af + γ^{2}f, f) = \underbrace{(Af, Af)}_{≥0} +γ^{2}(f, f) ≥ γ^{2}(f, f) .

Mit Lemma 3.4.1 folgt nun also (A^{∗}A + γ^{2}I) ist invertierbar und

∥(A^{∗}A + γ^{2}I)^{−1}∥ ≤ 1/γ^{2} .

Definiert man R_{γ} := (A^{∗}A + γ^{2}I)^{−1}A^{∗}, so folgt

R_{γ}g_{δ} = f und ∥R_{γ}∥ ≤ 1/γ^{2}∥A^{∗}∥ .

Nach der Stetigkeit und der Wohldefiniertheit ist nun noch die punktweise Konvergenz von R_{γ} gegen A^{+} zu zeigen.

3.4.2 Satz. Der Operator A : X → Y besitze eine Singulärwertzerlegung, dann gilt:

R_{γ}g ⟶^{γ→0} A^{+}g , falls g ∈ R(A) + R(A)^{⊥} .

Beweis. Betrachtet man die Singulärwertdarstellungen von A , A^{∗} und A^{∗}A

Af = ∑_{k=1}^{∞} σ_{k}(f, f_{k})g_{k} , A^{∗}g = ∑_{k=1}^{∞} σ_{k}(g, g_{k})f_{k}

A^{∗}Af = ∑_{k=1}^{∞} σ_{k}^{2}(f, f_{k})f_{k} ,

wobei {f_{k}} ein ONS in X und {g_{k}} ein ONS in Y ist, so folgt für die regu-3 REGULARISIERUNGSMETHODEN 57 larisierte Normalgleichung:

A^{∗}Af + γ^{2}f = ∑_{k=1}^{∞} σ_{k}^{2}(f, f_{k})f_{k} + γ^{2}f = ∑_{k=1}^{∞} σ_{k}(g, g_{k})f_{k} = A^{∗}g  (3.14)
⇒  f ist Linearkombination der f_{k}
⇒  f = ∑_{k=1}^{∞}(f, f_{k})f_{k} .

Durch Einsetzen der Linearkombination für f in (3.14) und anschließenden Koeffizientenvergleich folgt:

σ_{k}^{2}(f, f_{k}) + γ^{2}(f, f_{k}) = σ_{k}(g, g_{k}) , ∀k

⇒  (f, f_{k}) = \frac{σ_{k}}{σ_{k}^{2} + γ^{2}} (g, g_{k}) = \frac{1}{σ_{k}} · \frac{1}{(γ^{2}/σ_{k}^{2}) + 1} (g, g_{k}) ,  ∀k
⇒  R_{γ}g = f = ∑_{k=1}^{∞} \frac{1}{(γ^{2}/σ_{k}^{2}) + 1} · \frac{1}{σ_{k}} (g, g_{k})f_{k} .  (3.15)

Mit A^{+}g = ∑_{k=1}^{∞} \frac{1}{σ_{k}} (g, g_{k})f_{k} folgt also:

∥(R_{γ} − A^{+})g∥^{2} = ∑_{k=1}^{∞} \underbrace{|(\frac{1}{(γ^{2}/σ_{k}^{2}) + 1} − 1) \frac{1}{σ_{k}} (g, g_{k})|^{2}}_{:=c_{k}(γ)} .  (3.16)

Die Summe in (3.16) konvergiert gleichmäßig bezüglich γ, da

c_{k}(γ) = |(\frac{1}{(γ^{2}/σ_{k}^{2}) + 1} − 1) \frac{1}{σ_{k}} (g, g_{k})|^{2}
= \underbrace{|\frac{−γ^{2}/σ_{k}^{2}}{(γ^{2}/σ_{k}^{2}) + 1}|^{2}}_{≤1} |\frac{1}{σ_{k}} (g, g_{k})|^{2}
≤ |\frac{1}{σ_{k}} (g, g_{k})|^{2} =: c_{k}

und ∑_{k=1}^{∞} c_{k} = ∥A^{+}g∥^{2} < ∞.

Mit dem Weierstraßschen Majorantenkriterium und lim_{γ→0} c_{k}(γ) = 0 ∀k folgt dann

lim_{γ→0} ∥(R_{γ} − A^{+})g∥^{2} = 0 .

□

Der Beweis zeigt nicht nur, daß R_{γ}g → A^{+}g konvergiert, sondern gibt auch Aufschluß über die Wirkungsweise der Tikhonov–Regularisierung. So ist anhand von Gleichung (3.16) ersichtlich, daß für die Regularisierung lediglich

die Singulärwerte von A^{+} ersetzt wurden durch \frac{1}{(γ^{2}/σ_{k}^{2}) + 1} · \frac{1}{σ_{k}} .

Es wurde also in diesem Fall ein Filter F_{γ}(σ) verwendet, für den gilt:

F_{γ}(σ) = \frac{1}{(γ^{2}/σ_{k}^{2}) + 1} = \{ ∼ 1 , für große σ , ∼ 0 , für kleine σ .

Die explizite Durchführung einer Tikhonov–Regularisierung wird im 4. Abschnitt anhand eines nichtcharakteristischen Cauchy Problems dargestellt. Dort ist es möglich, den zu regularisierenden Operator A direkt anzugeben. Dies ist nicht immer der Fall. Betrachtet man zum Beispiel das in Abschnitt 2.5 vorgestellte Problem der Parameteridentifikation, so läßt sich durch Einführen von ”Basisfunktionen“ e_{1}(x) … e_{n}(x) und Definieren von

P_{n}(x) = ∑_{j=1}^{n} α_{j}e_{j}(x) das Problem zurückführen auf die Berechnung von

α = (α_{1} … α_{n}) derart, daß die Folge {P_{n}(x)} gegen den Koeffizienten a(x) konvergiert. Bezeichnet man hierzu mit u_{n}(x, t) die Lösung der Differentialgleichung

\begin{aligned}\frac{∂u_{n}}{∂t} (x, t) & = & \frac{∂}{∂x} (P_{n}(x)\frac{∂u_{n}}{∂x} (x, t)) + f(x, t) , & 0 < x < 1 , 0 < t ≤ T , \\ u_{n}(x, 0) & = & g(x) , & 0 < x < 1 , \\ u_{n}(0, t) & = & ϕ(t) , u_{n}(1, t) = ψ(t) , & 0 < t ≤ T \end{aligned}

(3.17)

und mit ũ(x, t) die Meßwerte für die Funktion u(x, t), die die Differentialgleichung (3.17) mit dem Koeffizienten a(x) erfüllt, so ist folgende Minimierungsaufgabe zu lösen:

min_{α=(α_{1}...α_{n})} \frac{1}{2} ∫_{0}^{T} ∫_{0}^{1} [u_{n}(x, t) − ũ(x, t)]^{2} dx dt .  (3.18)

Oft ist es nötig, weitere a–priori–Informationen in die Aufgabe einzubringen, um neben einer Näherung u_{n}(x, t) für ũ(x, t) auch tatsächlich die korrekte Näherung für a(x) zu bekommen, da für u_{n}(x, t) → u(x, t) nicht P_{n}(x) → a(x) gelten muß. Um die a–priori–Informationen einzubringen, ersetzt man das Integral in (3.18) durch

\begin{aligned}J(α) & = & \frac{1}{2} ∫_{0}^{T} ∫_{0}^{1} [u_{n}(x, t) − ũ(x, t)]^{2} dx dt \\ & & +\frac{λ}{2} ∫_{0}^{1} [u_{n}(x, T) − ũ(x, T)]^{2} dx \\ & & +\frac{1}{2} ∑_{j=1}^{l} μ_{j} [P_{n}(x_{j}) − ã_{j}]^{2} , \end{aligned}  (3.19)

wobei ã_{j} eventuell vorher bekannte Näherungswerte für a(x_{j}) an einigen wenigen Stellen x_{j} sind. Das 2. Integral berücksichtigt die physikalische Begebenheit, daß die Messungen für t = T oft bessere Werte liefern als für kleinere Werte von t. Sollte dies nicht der Fall sein, oder keine a–priori–Informationen über das Verhalten von a(x) vorliegen, so sind λ , μ_{1} … μ_{l} Null zu setzen.

Eine Tikhonov–Regularisierung ergibt sich nun aus

min_{α=(α_{1}...α_{n})} \{J(α) + \frac{γ^{2}}{2} ∥P_{n}∥^{2}\} = min_{α=(α_{1}...α_{n})} \{I(α)\} .  (3.20)

Die Aufgabe (3.20) ist endlichdimensional und läßt sich mit Hilfe von Gradientenverfahren lösen. Ein Problem bei diesen Verfahren ist oft die effiziente Berechnung des Gradienten. Hierzu setze man

I_{k} := \frac{∂I}{∂α_{k}} , sowie w_{k} := \frac{∂u_{n}}{∂α_{k}} , k = 1 … n ,

und es ergibt sich

∂I/∂α_{k} = I_{k} = ∫_{0}^{T} ∫_{0}^{1} [u_{n}(x, t) − ũ(x, t)] w_{k}(x, t)dx dt (3.21)
+λ ∫_{0}^{1} [u_{n}(x, T) − ũ(x, T)] w_{k}(x, T)dx (3.22)
+ ∑_{j=1}^{l} μ_{j} [P_{n}(x_{j}) − ã_{j}] e_{k}(x_{j})  (3.23)
+γ^{2} ∫_{0}^{1} P_{n}(x) · e_{k}(x)dx .  (3.24)

Das Problem bei der Berechnung des Gradienten besteht nun darin, die Integrale (3.21) und (3.22) zu bestimmen. Differenziert man die Gleichungen in (3.17), so erhält man für die Werte von w_{k}

\frac{∂w_{k}}{∂t} (x, t) = \frac{∂}{∂x} (P_{n}(x)\frac{∂w_{k}}{∂x} (x, t)) + \frac{∂}{∂x} (e_{k}(x)\frac{∂u_{n}}{∂x} ) ,

w_{k}(x, 0) = 0 ,  (3.25)
w_{k}(0, t) = w_{k}(1, t) = 0 .

Durch Lösen einer partiellen Differentialgleichung der Art (3.17) für u_{n} und n partiellen Differentialgleichung der Art (3.25) für die w_{k} ist es nun möglich, die Integrale (3.21) und (3.22) zu bestimmen und die Gradienten zu berechnen.

Um das aufwendige Berechnen der w_{k} zu umgehen, betrachte man folgenden Trick: Sei v(x, t) die Lösung der Differentialgleichung

\begin{matrix}v_{t}(x, t) + [P_{n}(x)v_{x}(x, t)]_{x} = u_{n}(x, t) − ũ(x, t) , & 0 < x < 1 , 0 < t ≤ T , \\ v(x, T) = −λ (u_{n}(x, T) − ũ(x, T)) , & 0 < x < 1 , \\ v(0, t) = 0 , v(1, t) = 0 , & 0 < t ≤ T . \end{matrix}

(3.26)

Mit diesem v folgt dann:

∫_{0}^{T} ∫_{0}^{1}(u_{n}(x, t) − ũ(x, t)) w_{k}(x, t)dx dt + \underbrace{λ ∫_{0}^{1} (u_{n}(x, T) − ũ(x, T)) w_{k}(x, T)dx}_{=:I_{2}}

=^{(3.26)} ∫_{0}^{T} ∫_{0}^{1} {v_{t}(x, t) + [P_{n}(x)v_{x}(x, t)]_{x}} w_{k}(x, t)dx dt + I_{2}
=^{part. Int.} ∫_{0}^{1} [v(x, t)w_{k}(x, t)]_{0}^{T} dx − ∫_{0}^{1} ∫_{0}^{T} v(x, t)(w_{k})_{t}(x, t)dt dx

+ ∫_{0}^{T} [P_{n}(x)v_{x}(x, t)w_{k}(x, t)]_{0}^{1} dt − ∫_{0}^{T} ∫_{0}^{1} P_{n}(x)v_{x}(x, t)(w_{k})_{x}(x, t)dx dt + I_{2}

=^{(3.25)}  \underbrace{∫_{0}^{1} v(x, T)w_{k}(x, T)dx}_{=−I_{2} siehe (3.26)} − ∫_{0}^{1} ∫_{0}^{T} v(x, t)(w_{k})_{t}(x, t)dt dx

− ∫_{0}^{T} ∫_{0}^{1} P_{n}(x)v_{x}(x, t)(w_{k})_{x}(x, t)dx dt + I_{2}
=^{(3.25)} − ∫_{0}^{1} ∫_{0}^{T} v(x, t) {[P_{n}(x)(w_{k})_{x}(x, t)]_{x} + [e_{k}(x)(u_{n})_{x}(x, t)]_{x}} dx dt
− ∫_{0}^{T} ∫_{0}^{1} P_{n}(x)v_{x}(x, t)(w_{k})_{x}(x, t)dx dt
=^{part. Int.} − ∫_{0}^{T} [v(x, t)P_{n}(x)(w_{k})_{x}(x, t)]_{0}^{1} dt + ∫_{0}^{T} ∫_{0}^{1} v_{x}(x, t)P_{n}(x)(w_{k})_{x}(x, t)dx dt

− ∫_{0}^{T} [v(x, t)e_{k}(x)(u_{n})_{x}(x, t)]_{0}^{1} dt + ∫_{0}^{T} ∫_{0}^{1} v_{x}(x, t)e_{k}(x)(u_{n})_{x}(x, t)dx dt

− ∫_{0}^{T} ∫_{0}^{1} v_{x}(x, t)P_{n}(x)(w_{k})_{x}(x, t)dx dt
=^{(3.26)} ∫_{0}^{T} ∫_{0}^{1} v_{x}(x, t)e_{k}(x)(u_{n})_{x}(x, t)dx dt
=: L_{k} .

Der Algorithmus läßt sich nun wie folgt beschreiben:

- Initialisiere α^{(0)}

- Berechne u_{n} = u_{n}(α^{ν}, (x, t)) durch Lösen der Differentialgleichung (3.17)

- Berechne v(x, t) durch Lösen der Differentialgleichung (3.26)

- Werte das Integral L_{k} für k = 1 … n aus

- Berechne I_{k} mit Hilfe von L_{k}, (3.23) und (3.24)

- Berechne α^{ν+1} unter Verwendung von I_{k} so, daß I(α^{ν+1}) ≤ I(α^{ν})

Abschließend sei noch bemerkt, daß die Differentialgleichung (3.26), obwohl sie rückwärts in der Zeit verläuft, keine Instabilität in den Algorithmus einbringt. Führt man eine Transformation τ = T − t durch, so ergibt sich für (3.26)

[P_{n}(x)v_{x}(x, τ)]_{x} = v_{τ}(x, τ) + u_{n}(x, τ) − ũ(x, τ) ,
v(x, 0) = −λ (u_{n}(x, 0) − ũ(x, 0)) ,
v(0, τ) = v(1, τ) = 0 .

Die Lösung v(x, τ) dieser Differentialgleichung läßt sich nun stabil berechnen.

4 Numerische Verfahren

4.1 Das diskrete Problem und die Kondition der Matrix

In Abschnitt 2.3 wurde das nichtcharakteristische Cauchy Problem (2.22) umformuliert in eine Volterrasche Integralgleichung 2. Art der Form

Af = ∫_{0}^{t} k_{α}(t − τ)f(τ)dτ = g(t)  (4.1)

mit

k_{α}(t) = \frac{1}{\sqrt{πα}t^{3/2}} ∑_{n=0}^{∞}(2n + 1)(−1)^{n} exp (\frac{−(2n + 1)^{2}}{4αt} ) .  (4.2)

Berechnet man für einige Werte von α den Verlauf des Integralkerns k_{α}(t) (siehe Bild 1.1 ), so ist sofort zu erkennen, daß die Kurven für relativ kleine α ”glatter“ verlaufen, als die für große Werte von α. Dieses Verhalten von k_{α} wird sich, wie im Folgenden noch gezeigt wird, auf die Kondition des Problems auswirken.

Bild 1.1: Kernfunktion für α = 5 (gepunktete Linie), α = 1 (gestrichelte Linie) und α = 0.5 (durchgezogene Linie)

Um das Problem zu diskretisieren betrachte man die Mittelpunktsregel. Das Intervall [0, 1] wird in n gleiche Teile der Länge h = 1/n aufgeteilt und die Integralgleichung (4.1) ersetzt durch das Gleichungssystem

h · ∑_{i=1}^{j} k_{α}(t_{j} − t_{i−1/2}) · f(t_{i−1/2}) = g(t_{j}) ,  j = 1 … n ,  (4.3)

wobei t_{i−1/2} := (i − \frac{1}{2}) · h und t_{j} = j · h ist. Setzt man für feste Werte von α , g_{j} = g(t_{j}) , f_{i} = f(t_{i−1/2}) sowie k_{j} = k_{α}(t_{j−1/2}), so erhält man das diskrete Problem

Kf = g

mit

f = ( f_{1} ; ⋮ ; f_{n} ) , g = ( g_{1} ; ⋮ ; g_{n} ) , K = ( \begin{matrix}k_{1} \\ k_{2} & k_{1} & 0 \\ ⋮ & ⋮ & ... \\ k_{n} & k_{n−1} & … & k_{1} \end{matrix} ) .  (4.4)

Berechnet man die Singulärwertzerlegung der so gewonnenen Systemmatrix für n = 50 und verschiedene Werte von α, so wird für kleinere Werte von α ein immer stärkeres Abfallen der Singulärwerte gegen 0 ersichtlich.

Bild 1.2: Singulärwerte für α = 5 (durchgezogene Linie), α = 1 (gestrichelte Linie) und α = 0.5 (gepunktete Linie)

Für die Kondition von K, K(K) = σ_{1}/σ_{n} , σ_{1} … σ_{n} Singulärwerte von K erhält man

\begin{matrix}α =5 : K(K) ∼ & 5.0 · 10^{17} \\ α =1 : K(K) ∼ & 7.1 · 10^{25} \\ α =0.5 : K(K) ∼ & 9.8 · 10^{36} . \end{matrix}

Es ist nun zu erkennen, daß, je kleiner α gewählt wird, d.h. um so ”glatter“ die Kernfunktion der Integralgleichung verläuft, um so schlechter ist das Problem konditioniert. Anhand der hohen Kondition der Matrix ist sofort ersichtlich,(Bild 1.3 ) daß ein direktes Auflösen des Gleichungssystems nicht möglich ist.

Bild 1.3: α = 5, Orginalfunktion f (durchgezogene Linie), Rekonstruktion ohne Regularisierung f̃ (gestrichelte Linie), ungestörte Datenfunktion g (gepunktete Linie). (n = 50)

4.2 Regularisierung des Diskreten Problems und numerische Ergebnisse

Betrachtet man in Bild 1.2 den Verlauf der Singulärwerte, so ist sofort ersichtlich, daß die kleinen Singulärwerte für die großen Fehler bei der Rekonstruktion einer Lösung verantwortlich sind. Verzichtet man nun bei der Konstruktion einer Lösung auf diese Werte, so ergibt sich die abgeschnittene Singulärwertzerlegung

A^{+}g = ∑_{σ_{k}>γ} 1/σ_{k}(g, g_{k})f_{k} .  (4.5)

Wählt man γ im Bereich von 10^{−2}, so erhält man (siehe Bild 2.1 ) durchaus akzeptable Ergebnisse. Der Nachteil dieser Methode liegt darin, daß sie keinerlei Gebrauch von der speziellen Struktur der Matrix in (4.4) macht. Da es sich hier um eine untere Dreiecks–Toeplitz–Matrix handelt, wäre eine direkte Inversion ohne weiteres in O(n^{2}) R.O. möglich, während der Rechenaufwand für die abgeschnittene Singulärwertzerlegung bei O(n^{3}) R.O. anzusiedeln ist. Unter Ausnutzung der speziellen Struktur möchte man nun auch für die Regularisierungen einen Rechenaufwand von O(n^{2}) nicht überschreiten.

Bild 2.1: Rekonstruktion mit einer abgeschnittenen SVD. Orginalfunktion f (durchgezogene Linie), Näherung f̃ (gestrichelte Linie) und Datenfunktion g (gepunktete Linie). (n = 50, α = 1)

Um eine Tikhonov–Regularisierung durchzuführen, betrachte man anstelle von Kf = g die Aufgabe

min_{f} {∥Kf − g∥_{2}^{2} + γ^{2}∥f∥_{2}^{2}}

oder die verallgemeinerte Aufgabe

min_{f} {∥Kf − g∥_{2}^{2} + γ^{2}∥Lf∥_{2}^{2}} ,  (4.6)

wobei L eine Matrix ist, die z.B. eine Diskretisierung der 1.Ableitung, also eine (n − 1) × n Matrix der Form

L = \frac{1}{h} ( \begin{matrix}1 & −1 \\ & ... & ... \\ & & 1 & −1 \end{matrix} )

sein kann. Das Problem (4.6) ist nun äquivalent zu

min_{f} ∥\binom{K}{γL}f − \binom{g}{0}∥_{2}^{2} .

Für die Beschreibung des Algorithmus und zur einfacheren Berechnung ist es sinnvoll, die untere Toeplitz–Matrix durch eine obere Toeplitz–Matrix

K = ( \begin{matrix}k_{1} & k_{2} & … & k_{n} \\ & k_{1} & … & k_{n−1} \\ & 0 & ... & ⋮ \\ & & & k_{1} \end{matrix} )

zu ersetzen. Dies ist durch einfaches Umsortieren der Komponenten von f =

( f_{1} ; ⋮ ; f_{n} ) und g = ( g_{1} ; ⋮ ; g_{n} ) zu erreichen. Für den Algorithmus seien also K

und L obere Dreiecks–Toeplitz–Matrizen. Zu lösen ist dann die Aufgabe

min_{f} ∥\binom{K}{γL}f − \binom{g}{0}∥_{2}^{2}  (4.7)

mit der Matrix

\binom{K}{γL} = ( \begin{matrix}k_{1} & k_{2} & … & k_{n} \\ & k_{1} & … & k_{n−1} \\ & 0 & ... & ⋮ \\ & & & k_{1} \\ l_{1} & l_{2} & … & l_{n} \\ & l_{1} & … & l_{n−1} \\ & 0 & ... & ⋮ \\ & & & l_{1} \end{matrix} ) ,

wobei l_{i} das γ fache der entsprechenden Komponente von L darstellt. Die Minimierungsaufgabe (4.7) ist nach Satz 3.1.2 äquivalent zum Lösen der Normalgleichung

\binom{K}{γL}^{∗}\binom{K}{γL}f = \binom{K}{γL}^{∗}\binom{g}{0} .  (4.8)

Führt man eine Q–R Zerlegung von \binom{K}{γL} durch, also

\binom{K}{γL} = Q · \binom{K_{γ}}{0}

mit einer orthogonalen Matrix Q und einer oberen Dreiecksmatrix K_{γ}, so ergibt sich für (4.8)

\binom{K_{γ}}{0}^{T} Q^{T} Q\binom{K_{γ}}{0}f = \binom{K_{γ}}{0}^{T} Q^{T} \binom{g}{0}

oder

\binom{K_{γ}}{0}^{T} \binom{K_{γ}}{0}f = \binom{K_{γ}}{0}^{T} Q^{T} \binom{g}{0} .

Setzt man g̃ = Q^{T} \binom{g}{0} = (g̃_{1} ; g̃_{2}) und f̃ = \binom{K_{γ}}{0}f = \binom{f̃_{1}}{0}, so erhält man

(K_{γ}^{T} , 0) · \binom{f̃_{1}}{0} = (K_{γ}^{T} , 0) · (g̃_{1} ; g̃_{2})

⇔ K_{γ}^{T} f̃_{1} = K_{γ}^{T} g̃_{1} ,

und da mit \binom{K}{γL} auch K_{γ} und K_{γ}^{T} vollen Rang haben gilt

f̃_{1} = g̃_{1} oder K_{γ}f = g̃_{1} .  (4.9)

Eine Lösung von (4.7) ergibt sich nun aus dem einfachen Lösen des Gleichungssystems (4.9) mit der oberen Dreiecksmatrix K_{γ}. Da sich die Matrix Q aus Produkten von Givens–Rotationen zusammensetzt, wird sie nicht explizit berechnet, sondern es werden die einzelnen Rotationen direkt auf \binom{K}{γL} und \binom{g}{0} angewandt.

Im erten Schritt des Algorithmus werden die Elemente der Hauptdiagonalen der Matrix L eleminiert. Hierzu werden nacheinander Givens–Rotationen auf die Zeilenpaare (1, n+1), (2, n+2), …, (n, 2n) durchgeführt. Anhand der Toeplitzstruktur der Matrizen K und L ist sofort ersichtlich, daß bei allen Rotationen der Drehwinkel derselbe ist und sich somit die Rotationen der Zeilenpaare (2, n + 2), …, (n, 2n) aus der Rotation von (1, n + 1) ergeben. Das System

[{K \atop γL} ⋮ ⋮ {g \atop 0}] = [ \begin{matrix}k_{1} & k_{2} & … & k_{n} ⋮ & g_{1} \\ k_{1} & … & k_{n−1}⋮ & g_{2} \\ 0 & ... & ⋮ ⋮ & ⋮ \\ k_{1} ⋮ & g_{n} \\ l_{1} & l_{2} & … & l_{n} ⋮ & 0 \\ l_{1} & … & l_{n−1} ⋮ & 0 \\ 0 & ... & ⋮ ⋮ & ⋮ \\ l_{1} ⋮ & 0 \end{matrix} ]

wird so überführt in das System

[ \begin{matrix}k′_{1} & k′_{2} & k′_{3} & … & k′_{n} ⋮ & g′_{1} \\ k′_{1} & k′_{2} & … & k′_{n−1}⋮ & g′_{2} \\ ... & ... & ⋮ ⋮ & ⋮ \\ 0 & k′_{1} & k′_{2} ⋮ & g′_{n−1} \\ k′_{1} ⋮ & g′_{n} \\ 0 & l′_{2} & l′_{3} & … & l′_{n} ⋮ & g′_{n+1} \\ 0 & l′_{2} & … & l_{n−1} ⋮ & g′_{n+2} \\ ... & ... & ⋮ ⋮ & ⋮ \\ 0 & 0 & l′_{2} ⋮ & g′_{2n−1} \\ 0 ⋮ & g′_{2n} \end{matrix} ]

Im zweiten Schritt werden nun die Elemente der 1.Nebendiagonalen der (transformierten) Matrix L eliminiert. Hierzu werden wieder Givens–Rotationen, diesmal auf den Zeilenpaaren (2, n + 1), (3, n + 2), …, (n, 2n − 1), durchgeführt. Da auch hier die Drehwinkel wieder dieselben sind, lassen sich die Rotationen aus der ersten Rotation berechnen. Nach dem 2.Schritt hat das System die Gestallt:

[ \begin{matrix}k′_{1} & k′_{2} & k′_{3} & k′_{4} & … & k′_{n} ⋮ & g′_{1} \\ k1^{′′} & k2^{′′} & k3^{′′} & … & kn−1^{′′}⋮ & g2^{′′} \\ ... & ... & ... & ⋮ ⋮ & ⋮ \\ 0 & k1^{′′} & k2^{′′} & k3^{′′} ⋮ & gn−2^{′′} \\ k1^{′′} & k2^{′′} ⋮ & gn−1^{′′} \\ k1^{′′} ⋮ & gn^{′′} \\ 0 & 0 & l3^{′′} & l4^{′′} & … & ln^{′′} ⋮ & gn+1^{′′} \\ 0 & 0 & l3^{′′} & … & ln−1^{′′} ⋮ & gn+2^{′′} \\ ... & ... & ... & ⋮ ⋮ & ⋮ \\ 0 & 0 & 0 & l3^{′′} ⋮ & g2n−2^{′′} \\ 0 & 0 ⋮ & g2n−1^{′′} \\ 0 ⋮ & g2n^{′} \end{matrix} ]

Man beachte, daß die Toeplitzstruktur der Matrix nur noch in den Zeilen 2 bis n besteht. Im dritten Schritt werden analog zu den ersten beiden Schritten Givens– Rotationen auf den Zeilen (3, n + 1), (4, n + 2), …, (n, 2n − 2) durchgeführt. Nach n Schritten hat das Gleichungssystem die gewünschte Form

[{K_{γ} \atop 0} ⋮ ⋮ g̃_{1} ; g̃_{2} ] = [ \begin{matrix}k_{1}^{1} & k_{2}^{1} & … & k_{n}^{1} ⋮ & g_{1}^{1} \\ k_{1}^{2} & … & k_{n−1}^{2}⋮ & g_{2}^{2} \\ 0 & ... & ⋮ ⋮ & ⋮ \\ k_{1}^{n} ⋮ & g_{n}^{n} \\ ⋮ & g_{n+1}^{n} \\ 0 & ⋮ & ⋮ \\ ⋮ & g_{2n}^{1} \end{matrix} ] ,

wobei k_{1}^{1} = k′_{1}, k_{1}^{2} = k′′_{1} u.s.w. sind. Es ist nun leicht zu sehen, daß im i.Schritt 4(n − i + 1) Elemente neu berechnet werden müssen. Bei Verwendung von Givens–Rotationen ist jede Transformation eines Elementes mit zwei Multiplikationen durchführbar. Für die Transformation des Gleichungssystems ergibt sich also ein Rechenaufwand von

2 · ∑_{i=1}^{n} 4(n − i + 1) ≈ 4n^{2} R.O.

Für den gesamten Algorithmus kommen dann noch O(n^{2}/2) R.O. für die Lösung des Gleichungssystems K_{γ}f = g̃_{1} hinzu.

Falls die Matrix L nicht quadratisch, also L eine (p×n) Matrix mit p < n ist, läßt sich der Algorithmus analog durchführen. In diesem Fall werden die Rotationen nur auf die Zeilen 1 … p, b.z.w. n + 1 … n + p durchgeführt. Die restlichen Zeilen werden getrennt betrachtet.

Um den Algorithmus zu testen, wurde die Datenfunktion g(t) dadurch erzeugt, daß für eine vorgegebene Funktion f(t) das direkte, gut gestellte Problem gelöst wurde. Die so erhaltene Datenfunktion wurde dann noch mit einem normalverteilten Fehler mit einer Standardabweichung im Bereich von 10^{−3} versehen. In Bild 2.2 und Bild 2.3 sind die Rekonstruktionen f̃ (gestrichelte Linie) die Orginalfunktionen f (durchgezogene Linie) und die gestörten Datenfunktionen g (gepunktete linie) aufgezeichnet. Es ist zu sehen, daß sowohl für

L = D = \frac{1}{h} ( \begin{matrix}1 & −1 \\ 0 & ... & ... & 0 \\ & & 1 & −1 \end{matrix} ) als auch für L = I = ( \begin{matrix}1 \\ 0 & ... & 0 \\ & & 1 \end{matrix} )

akzeptable Näherungen berechnet wurden. Dabei ist zu erkennen, daß erwartungsgemäß die Näherung für α = 1 besser ist als die für α = 0.5. Dies ist mit der höheren Kondition der Matrix für α = 0.5 zu erklären.

Bild 2.2: Rekonstruktionen für α = 1. Linke Seite: L = I, γ = 0.012. Rechte Seite: L = D, γ = 0.021. (n = 50)

Bild 2.3: Rekonstruktionen für α = 0.5. Linke Seite: L = I, γ = 0.0063. Rechte Seite: L = D, γ = 0.01. (n = 50)

Um den Algorithmus an realen Problemen zu testen, betrachte man folgende Situation: Gegeben ist ein Schlauch, dessen Enden verschlossen sind. Dieser wird mit einer konstanten Temperatur von c = 30^{∘} in einen Heißluftofen gelegt, in dem eine Temperatur von c.a. 110^{∘} vorhanden ist. Während der Aufheizphase wird nun an einer Stelle des Schlauches der Temperaturverlauf an der Innenseite und an der Außenseite gemessen.

[figure] [figure]

[figure]

[figure]

[figure]

Meßpunkt innen

Meßpunkt außen

Transformiert man das Problem in ein räumlich eindimensionales und betrachtet die Temperaturverteilung an der eingezeichneten Linie, so liegt folgende Situation vor:

\begin{aligned}u_{t} & = & α u_{xx} , & 0 ≤ x < l ,  0 < t ≤ T , \\ u(x, 0) & = & 30 , & 0 ≤ x < l , \\ u(0, t) & = & g(t) , & 0 < t ≤ T , \\ u_{x}(0, t) & = & 0 , & 0 < t ≤ T , \\ u(l, t) & = & f(t) , & 0 < t ≤ T . \end{aligned}  (4.10)

[figure: x; u = 30^{∘}; t; T; l; (0,0); u_{t} = αu_{xx}; \begin{matrix}u = g \\ u_{x} = 0 \end{matrix}; u = f]

Zu gegebenen Meßwerten an der Innenseite (g(t)) möchte man nun eine Näherung f̃(t) für die Meßwerte f(t) an der Außenseite berechnen. Der Wärmeleitkoeffizient des verwendeten Materials hat den Wert α = 2.0 · 10^{−7} \frac{m^{2}}{sec} . Transformiert man das Problem auf das Einheitsquadrat [0, 1] × [0, 1], so erhält man für l = 10mm und T = 960 sec einen Wert für α von α = 1.92.

Führt man die Berechnungen mit dem beschriebenen Algorithmus für L = I, n = 80 und dem Regularisierungsparameter γ = 0.08 durch, so zeigt sich (Bild 2.4 ), daß die Funktion f̃ bis auf die letzten Komponenten eine akzeptable Näherung für f darstellt. Die letzten Werte von f̃ jedoch fallen plötzlich stark gegen 0 ab.

Bild 2.4: Rekonstruktion mit realen Daten. Gemessene Werte (durchgezogeneLinie), berechnete Rekonstruktion (gestrichelte Linie). (n = 80, α = 1.92)

Um dieses Phänomen zu erklären, benötigt man folgenden Satz über das Verhalten der Singulärvektoren in den letzten Komponenten.

4.2.1 Satz. Es sei K = ( \begin{matrix}k_{1} \\ ⋮ & ... & 0 \\ k_{n} & … & k_{1} \end{matrix} ) eine n × n untere Toeplitz–

Matrix. K habe die Singulärwertzerlegung Kf = ∑_{k=1}^{n} σ_{k} (f, f_{k}) g_{k} ,

f, f_{k}, g_{k} ∈ R I ^{n}. Bezeichnet man mit f_{k}(i) die i–te Komponente von f_{k} (für g_{k} analog) so gilt:

a) K hat den n–fachen Eigenwert k_{1} zum Eigenvektor e_{n}, wobei e_{n} der n–te Einheitsvektor, also e_{n} = (0, …, 0, 1)^{T} ist.

b) Für die Singulärvektoren von K gilt: f_{k}(j) = g_{k}(n + 1 − j) , k, j = 1 … n.

4 NUMERISCHE VERFAHREN 75

c) |f_{k}(n)| ≤ \frac{k_{1}}{σ_{k}} , k = 1 … n.

Beweis. a) trivial.

b) Mit der Matrix P = ( \begin{matrix}& & 1 \\ 0 & .. & ˙ 0 \\ 1 \end{matrix} ) folgt

K^{#} = P K = ( \begin{matrix}k_{n} & … & k_{1} \\ ⋮ & .. & ˙ 0 \\ k_{1} \end{matrix} ) .

Da alle k_{i} ∈ R I und (K^{#})^{T} = K^{#} ist K^{#} symetrisch. Also gibt es ein ONS (w_{k})_{k=1...n} aus Eigenvektoren von K^{#}, so daß

K^{#}w_{k} = σ_{k}w_{k} , k = 1 … n .

Definiert man W := (w_{1}, …, w_{n}) als Matrix der Eigenvektoren und Σ := ( \begin{matrix}σ_{1} \\ 0 & ... & 0 \\ & & σ_{n} \end{matrix} ) als Matrix der Eigenwerte, so erhält man

K^{#}W = W Σ

⇒ P KW = W Σ
⇒ KW = P W Σ
⇒ K = P W ΣW^{T} , da W^{T} W = I .

Setzt man f_{k} = w_{k} = k–Zeile von W^{T}, sowie g_{k} = P w_{k} = k–Spalte von P W,

so ist Kf = ∑_{k=1}^{n} σ_{k}(f, f_{k})g_{k} und f_{k}(j) = g_{k}(n + 1 − j).

c) Für K gilt nach a) und b):

Ke_{n} = ∑_{k=1}^{n} σ_{k}f_{k}(n) · g_{k} = k_{1}e_{n}

⇒ ∑_{k=1}^{n} σ_{k}f_{k}(n)g_{k}(n) = k_{1}
⇒^{b)}  ∑_{k=1}^{n} σ_{k}f_{k}(n)f_{k}(1) = k_{1} .  (4.11)

Da Kf_{k} = σ_{k}g_{k} folgt in der ersten Komponente

k_{1}f_{k}(1) = σ_{k}g_{k}(1)
⇒^{b)} k_{1}f_{k}(1) = σ_{k}f_{k}(n) .  (4.12)

Eingesetzt in (4.11) folgt:

∑_{k=1}^{n} σ_{k} \frac{k_{1}}{σ_{k}} f_{k}(1)f_{k}(1) = k_{1}

⇒ ∑_{k=1}^{n} (f_{k}(1))^{2} = 1
⇒ |f_{k}(1)| ≤ 1 k = 1 … n .

Mit (4.12) folgt dann schließlich

|f_{k}(n)| ≤ \frac{k_{1}}{σ_{k}} .

□

Betrachtet man nun die Singulärvektoren zu den Singulärwerten, für die gilt σ_{k} ≥ γ, also zu den relativ großen Singulärwerten, so folgt mit Satz 4.2.1 b):

|f_{k}(n)| ≤ \frac{k_{1}}{γ} .  (4.13)

Spaltet man die verallgemeinerte Inverse in zwei Teile auf

f(n) = (A^{+}g) (n) = ∑_{σ_{k}≥γ} \frac{1}{σ_{k}} (g, g_{k})f_{k}(n) + ∑_{σ_{k}<γ} \frac{1}{σ_{k}} (g, g_{k})f_{k}(n) ,

so erhält man für den ersten Teil mit (4.13)

∑_{σ_{k}≥γ} \frac{1}{σ_{k}} (g, g_{k})f_{k}(n) ≤ \frac{1}{γ} ∥g∥\frac{k_{1}}{γ} = \frac{k_{1}}{γ^{2}} ∥g∥ =: w .

Da bei dem betrachteten Problem γ > 10^{−2}, k_{1} < 10^{−8} und ∥g∥ < 10^{4} sind, folgt in diesem Fall w < 1. Der Wert von f(n) wird also in der Hauptsache durch den zweiten Teil der Summe bestimmt. Dieser wird aber bei den Regularisierungen stark gefiltert, oder sogar ganz weggelassen. Für die Tikhonov–Regularisierung erhält man z.b.:

|f̃(n)| ≤ ∑_{k=1}^{n} | \frac{1}{\frac{γ^{2}}{σ_{k}} + σ_{k}} (g, g_{k})f_{k}(n)|

≤ ∑_{k=1}^{n} | \frac{1}{\frac{γ^{2}}{σ_{k}} + σ_{k}} (g, g_{k})\frac{k_{1}}{σ_{k}} |

= ∑_{k=1}^{n} \frac{k_{1}}{γ^{2} + σ_{k}^{2}} |(g, g_{k})|
≤ \frac{k_{1}}{γ^{2}} ∥g∥ .

Mit den Werten von oben folgt somit |f̃(n)| < 1. Die theoretischen Betrachtungen lassen sich auch numerisch nachvollziehen. So ist in den Bildern 4.5 bis 4.7 sofort zu erkennen, daß die Singulärvektoren f_{1} bis f_{79} in der letzten Komponente nahezu 0 sind. Nur f_{80} hat an dieser Stelle einen deutlich von 0 verschiedenen Wert.

Bild 2.5: Singulärvektoren f_{1} und f_{2}. (n = 80, α = 1.92)

Bild 2.6: Singulärvektoren f_{77} und f_{78}. (n = 80, α = 1.92)

Bild 2.7: Singulärvektoren f_{79} und f_{80}. (n = 80, α = 1.92)

Um das Problem zu lösen betrachte man folgende Methoden.

1. Es wird als a–priori–Information angenommen, daß sich die Werte für f in den letzten Komponenten nicht sehr stark verändern. Dann verwendet man für die Berechnung der letzten Komponenten die Werte der vorherigen Komponenten. Dies läßt sich z.B. dadurch erreichen, daß

anstelle der Matrix L = I die Matrix L = D = ( \begin{matrix}1 & −1 \\ 0 & ... & ... & 0 \\ & & 1 & −1 \end{matrix} )

verwendet wird. (Siehe Bild 2.8 )

2. Man erweitert die Dimension n des Problems und hängt an die Datenfunktion g eine Funktion g_{0} an, die relativ glatt gegen 0 abfällt. Die mit den neuen Daten berechnete Näherung wird dann nur in den ersten n Komponenten betrachtet. (Siehe Bild 2.9 )

Bild 2.8: Regularisierung mit L = D. (n = 80, α = 1.92)

Bild 2.9: Regularisierung mit L = I und erhöhter Dimension. (α = 1.92)

Bei der Berechnung von Funktionen, von denen bekannnt ist daß sie am Ende gegen 0 abfallen, oder bei der Verwendung der zweiten Methode, ist es naheliegend, die Informationen über das Verhalten der Näherungen in den letzten Komponenten mit in die Berechnungen einzubeziehen. Man betrachte also das Problem:

min_{f∈B} ∥Af − g∥_{2}^{2} ; B = {f : f_{j} = 0 , j = p + 1 … n} .

Dies ist äquivalent zu

min_{f̄} ∥K̄f̄− g∥_{2}^{2} mit K̄ = ( \begin{matrix}k_{1} \\ ⋮ & ... & 0 \\ k_{p} & … & k_{1} \\ ⋮ & & ⋮ \\ k_{n} & … & k_{n−p+1} \end{matrix} ) , f̄ = ( f_{1} ; ⋮ ; f_{p} ) .

Eine Lösung der Aufgabe ergibt sich wieder aus einer Q–R Zerlegung der Matrix K̄. Diesmal wird die Matrix Q so bestimmt, daß

K̄ = Q^{T} · \binom{R}{0} , g̃ = Q^{T} g = (g̃_{1} ; g̃_{2})

ist und dann das Gleichungssysthem Rf̄ = g̃_{1} gelöst. Betrachtet man die Ergebnisse des Verfahrens für p = 48, so ist ersichtlich, daß für α = 5 eine gute Rekonstruktion berechnet wurde, während bei α = 1 die Lösung stark oszilliert. (Siehe Bild 2.10 )

Bild 2.10: Rekonstruktion für α = 5 (links) und für α = 1 (rechts). (n = 50, p = 48)

Als letztes soll noch die in [10] vorgeschlagene Möglichkeit der Verwendung von gemischten Verfahren untersucht werden. Es soll also sowohl eine Tikhonov–Regularisierung, als auch die Projektion auf den Unterraum durchgeführt werden. Bei der Betrachtung der Ergebnisse dieses gemischten Verfahrens ist sofort ersichtlich, daß die Ergebnisse kaum von denen bei einer einfachen Tikhonov– Regularisierung erzielten Werten abweichen. (Siehe Bild 2.11 )

Bild 2.11: Rekonstruktion mit γ = 0.0063 und p = 48. (n = 50, α = 0.5)

Betrachtet man das Problem min_{f∈B} ∥Kf −g∥_{2}^{2} genauer, so ist zu erkennen, daß

sich die Lösung aus der Berechnung der kleinsten Quadrate–Lösung minimaler Norm von

min_{f} ∥K̃f − g∥ mit K̃ = ( \begin{matrix}k_{1} \\ ⋮ & ... & & 0 \\ k_{p} & … & k_{1} \\ k_{p+1} & … & k_{2} & 0 \\ ⋮ & & ⋮ & ⋮ & ... \\ k_{n} & … & k_{n−p+1} & 0 & … & 0 \end{matrix} )

ergibt. Bei der Berechnung der Singulärwertzerlegung von K̃ fällt auf, daß sich die Singulärvektoren f̃_{1} … f̃_{p} bzw. g̃_{1} … g̃_{p} (p = 48) von K̃ kaum von den entsprechenden Singulärvektoren von K unterscheiden.(Siehe Bild 2.12 und Bild 2.13 ) Nur in den Komponenten p+1 … n ist ein sichtlicher Unterschied zwischen den Werten von f̃_{k} und f_{k} zu erkennen, da hier die Singulärvektoren f̃_{k} erwartungsgemäß den Wert 0 haben. Bei den Singulärwerten verhält es sich ähnlich. Hier unterscheiden sich die ersten p Singulärwerte von K̃ und K kaum, während die Singulärwerte σ̃_{p+1} … σ̃_{n} von K̃ den Wert 0 besitzen.

4 NUMERISCHE VERFAHREN 83 Bild 2.12: Singulärvektoren f_{1}, f_{30} und f_{48} für K (linke Seite) und für K̃ (rechte Seite). (n = 50, α = 1, p = 48)

4 NUMERISCHE VERFAHREN 84 Bild 2.13: Singulärvektoren g_{1}, g_{30} und g_{48} für K (linke Seite) und für K̃ (rechte Seite). (n = 50, α = 1, p = 48)

Bezeichnet man mit f̃(i) = K̃^{+}g(i) die i–te Komponente der verallgemeinerten Lösung, so folgt:

f̃(i) = K̃^{+}g (i) = ∑_{k=1}^{p} \frac{1}{σ̃_{k}} (g, g̃_{k})f̃_{k}(i)
∼ ∑_{k=1}^{p} \frac{1}{σ_{k}} (g, g_{k})f_{k}(i) ,  i = 1 … p .

Da f̃ in den Komponenten p+1 … n nach Vorraussetzung den Wert 0 hat, ist die Projektion auf den Unterraum für eine geeignete Wahl von p vergleichbar mit einer abgeschnittenen Singulärwertzerlegung. In diesem Lichte betrachtet ist letztlich auch sofort ersichtlich, daß ein gemischtes Verfahren (siehe Bild 2.11 ) keinerlei Vorteile gegenüber einer einfachen Tikhonov–Regularisierung (siehe Bild 2.3 ) besitzt.

5 Zusammenfassung

Die in dieser Arbeit vorgenommene Aufstellung von schlecht gestellten Problemen bei der Wärmeleitungsgleichung kann aufgrund der hohen Anzahl derartiger Probleme nicht vollständig sein. Es ist jedoch klar zu erkennen, daß immer dann, wenn die ”Berechnungsrichtung“ entgegengesetzt zur Ausbreitungsrichtung der Wärme verläuft, ein schlecht gestelltes Problem vorliegt. Dieses läßt sich auch physikalisch nachvollziehen, da bei einem Wärmeleitungsprozess Schwankungen in der Temperaturverteilung ausgeglichen werden.

Um dennoch brauchbare Berechnungen durchführen zu können, wurden einige Regularisierungsmethoden untersucht. Dabei war zu erkennen, daß sich alle Regularisierungen auf verschiedene Arten von Filtern zurückführen lassen. Die Wirkungsweise dieser Filter kann jedoch verschieden sein. Bei Berechnungen mit realen Daten zeigt sich, daß die benutzten Regularisierungen gute Ergebnisse liefern. Es ist jedoch auch zu erkennen, daß durch die spezielle Wirkungsweise der einzelnen Methoden ungewollte Effekte auftreten können. Da diese Effekte von der Art der verwendeten Regularisierung abhängen, ist es sinnvoll, bei der Wahl der Regularisierung die spezielle Strucktur des Problems zu beachten, um so ein möglichst gutes Ergebnis zu erziehlen.

Literatur

[1] Baumeiser, J. : Stable solution of inverse problems. Vieweg, 1987.

[2] Cannon, J.R. : The one–dimensional heat equation. Adisson–Wesley, Menlo Park, 1984.

[3] Cannon, J.R., & Ewing R.E. : A direct numerical procedure for the heat equation. J. Math. Anal. and Appl., 56 (1976), 7–17.

[4] Cannon, J.R., & Hill, C.D. : Existence, uniqueness, stability, and monotone dependence in a Stefan problem for the heat equation. J. Math. Mech., 17 (1967), 1–20.

[5] Cannon, J.R., & Primicerio, M. : Remarks on the one–phase Stefan problem for the heat equation with flux prescribed on the fixed boundary. J. Math. Anal. and Appl., 35 (1971), 361–373.

[6] Carasso, A. : Determining surface temperatures from interior observations. SIAM J. Appl. Math., 42 (1982), 558–574.

[7] Doetsch, G. : Anleitung zum praktischem Gebrauch der Laplace– Transformation. R. Oldenbourg, München, 1961.

[8] Eldén, L. : An efficient algorithm for the regularization of ill– conditioned least squares problems with triangular Toeplitz matrix. SI-AM J. Sci. Comp., 5 (1984), 229–236.

[9] Eldén, L. : Approximations for a Cauchy problem for the heat equation. Inverse Probl., 3 (1987) 263–213.

[10] Eldén, L. : The numerical solution of a non–characteristic Cauchy problem for a parabolic equation. Numerical Treatment of Inverse Problems in Differential and Integral Equations, P. Deuflhard, H. Hairer (eds.), Birkhäuser, Bosten, 1983.

[11] Ewing, R. E. : The approximation of certain parabolic equations backward in time by Sobolev equations. SIAM J. Math. Anal., 6 (1975) 283–294.

[12] Ewing, R. E., Falk, R. S. : Numerical approximation of a Cauchy problem for a parabolic partial differential equation. Math. Comp., 33 (1979), 1125–1144.

[13] Gajewski, H., Zacharias, K. : Zur Regularisierung einer Klasse nichtkorrekter Probleme bei Evolutionsgleichungen. J. Math. Anal. Appl., 38 (1972) 784–789.

[14] Guenther, R., Hudspeth, R., McDougal, W., Gerlach, J. : Remarks on Parameter Identification I. Numer. Math., 47 (1985), 355–361.

[15] Kravaris, C., & Steinfeld, J. : Identification of parameters in distributed parameter systems by regularization. SIAM J. Control Optimization, 23 (1985), 217–241.

[16] Natterer, F. : Schlecht gestellte Probleme. Vorlesungsskript, Münster, Sommersemester 1986.

[17] Sneddon, I. : Fouriertransforms. Mc Graw–Hill, New York, 1951.

[18] Widder, D.V. : The heat equation. Academic Press, New York, 1975.

Ich versichere, daß ich die vorliegende Diplomarbeit selbständig verfaßt und keine anderen als die angegebenen Hilfsmittel verwendet habe.

Münster, 7. Februar 1994