Theoretische und numerische Untersuchung einiger …botkin/all_kaust_files/papers/...November 2010...

Technische Universitat Munchen

Fakultat fur Mathematik

Theoretische und numerischeUntersuchung einiger Optimalitats-

und Erhaltungsgleichungen

Diplomarbeit von Natalie Mayer

Aufgabensteller: Prof. Dr. Dr. h.c. mult. Karl-Heinz Hoffmann

Betreuer: Dr. Nikolai Botkin,

Abgabetermin: 20. November 2010

Ich erklare hiermit, dass ich die Diplomarbeit selbststandig und nur mit den angegebenenHilfsmitteln angefertigt habe.

Munchen, den 20 Dezember 2010

(Natalie Mayer)

Vorwort

Die vorliegende Diplomarbeit entstand wahrend meiner Tatigkeit am Lehrstuhl M6 furmathematische Modellbildung und Simulation an der Technische Universitat Munchen.

Mein Dank gilt Herrn Prof. Dr. Dr. h.c. mult. Karl-Heinz Hoffmann und Dr. Nikolai Botkinfur die Betreuung der Arbeit und die Schaffung der notwendigen Rahmenbedingungen, ohnedie diese Diplomarbeit nicht moglich gewesen ware.

Munchen, den 20 Dezember 2010

(Natalie Mayer)

Verwendete Formelzeichen und Abkurzungen

• R ist die Menge der reellen Zahlen.

• RN ist der euklidische Raum.

• 〈·, ·〉 ist das Skalarprodukt in RN .

• | · | ist der Betrag des Skalars, die Euklidische Norm des Vektors in RN , die Opera-

tornorm der N ×N -Matrix oder die Lange des großten Intervalls in der Partition.

• Ck(Q) ist die Menge aller k-mal stetig differenzierbaren Funktionen mit der Defini-tionsmenge Q.

• Ckb (Q) ist die Menge aller k-mal stetig differenzierbaren Funktionen mit der Defini-

tionsmenge Q, die Ableitungen bis zur Ordnung k sind beschrankt.

• Ck0 (Q) ist die Menge aller k-mal stetig differenzierbaren Funktionen mit der Defini-

tionsmenge Q und kompaktem Trager.

• C0,1b (Q) ist die Menge aller beschrankten reellen Lipschitz-stetigen Funktionen mit

der Definitionsmenge Q.

• ‖v‖ = supx∈Q

|v(x)| fur v ∈ C0,1b (Q).

• ‖Dv‖ ist die Lipschitz-Konstante fur v ∈ C0,1b (Q).

• Dv(x) = ( ∂v∂x1

(x), ∂v∂x2

(x), ..., ∂v∂xN

(x)) ist die Frechet-Ableitung der differenzierbarenFunktion.

• ‖Dkv‖ = supx∈Q

|Dkv(x)| fur v ∈ Ckb (Q), k = 1, 2 .

• supp(v) ist der Trager von v.

• H(τ, x, p) = −maxβ∈B

minα∈A

〈p, f(T − τ, x, α, β)〉 ist die Hamilton-Funktion fur das Diffe-

rentialspiel mit der Zustandsgleichung x = f(x, α, β).

• BN (x0, r) = x ∈ RN : |x− x0| ≤ r .

Inhaltsverzeichnis

1 Einleitung, Hamilton-Jacobi Gleichung 21.1 Mechanik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2 Geometrische Optik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.3 Steuerungstheorie, Spieltheorie . . . . . . . . . . . . . . . . . . . . . . . . . 51.4 Transportgleichungen, Erhaltungsgleichungen . . . . . . . . . . . . . . . . . 71.5 Zielsetzung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2 Hamilton-Jacobi-Gleichung aus der Konflikt-Steuerungstheorie mit Zustandsbe-schrankungen 92.1 Kostenfunktional und Wertefunktion . . . . . . . . . . . . . . . . . . . . . . 92.2 Formulierung der Bedingungen fur die viskose Losungen . . . . . . . . . . . 12

3 Ein stabiles Verfahren zur Losung von Hamilton-Jacobi-Gleichung aus der Konflikt-Steuerungstheorie mit Zustandsbeschrankungen 133.1 Abstrakter Zeit-Schritt-Operator . . . . . . . . . . . . . . . . . . . . . . . . 133.2 Zwei Beispiele des Zeit-Schritt-Operators . . . . . . . . . . . . . . . . . . . 143.3 Approximation der Wertefunktion . . . . . . . . . . . . . . . . . . . . . . . 153.4 Punktweise Konvergenz der Approximation . . . . . . . . . . . . . . . . . . 163.5 Optimale Strategien fur die Spieler . . . . . . . . . . . . . . . . . . . . . . . 273.6 Modellbeispiele . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

3.6.1 Nichtlineares Pendel (Nonlinear pendulum) . . . . . . . . . . . . . . 283.6.2 Materialpunkt (Material point) . . . . . . . . . . . . . . . . . . . . . 293.6.3 Das Spiel der zwei Autos (Game of two cars) . . . . . . . . . . . . . 30

4 Anwendung: Optimierung von Kuhlprotokollen bei der Kryopreservation von le-benden Zellen 324.1 Mathematisches Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334.2 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

5 Anwendung: Mehrphasenstromungen 415.1 Mathematisches Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

5.1.1 Phasen-Gleichungen . . . . . . . . . . . . . . . . . . . . . . . . . . . 435.1.2 Globale Gleichungen . . . . . . . . . . . . . . . . . . . . . . . . . . . 445.1.3 Gemeinsames System . . . . . . . . . . . . . . . . . . . . . . . . . . 45

5.2 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 465.2.1 Linienmethode zur Berechnung des Druckes . . . . . . . . . . . . . . 465.2.2 Relaxation-Algorithmus zur Losung des gemeinsamen Systems, Upwind-

Approximationen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

6 Zusammenfassung 58

1

1 Einleitung, Hamilton-Jacobi Gleichung

Der Typ der in dieser Diplomarbeit betrachteten partiellen Differentialgleichungen ist dieHamilton-Jacobi Gleichung erster Ordnung. Die Hamilton-Jacobi Gleichungen

∂u

∂t+H(t, x,Dxu) = 0 (1.1)

sind typische nichtlineare Differentialgleichungen, die in verschiedenen Gebieten der Wis-senschaft auftreten (Mathematik, Physik).

Das Cauchy-Problem fur Hamilton-Jacobi Gleichungen lautet

∂u

∂t+H(t, x,Dxu) = 0 in (0,∞) × R

N , (1.2)

u(0, x) = u0(x) auf RN , (1.3)

wobei H die Hamilton-Funktion ist und u ist die gesuchte Funktion. Die Koordinaten xsind Ortskoordinaten, Dxu Impulskoordinaten und t ist die Zeit.

Das Cauchy-Problem besitzt klassische Losungen aus C1, wenn H und u0 ausreichend glattauf einem kleinen Zeitintervall sind. Die Losungen kann man mit der klassischen Methodeder Charakteristiken bestimmen. Die Idee ist dabei folgende: Um die partielle Differential-gleichung in ein System von gewohnlichen Differentialgleichungen (charakteristische Glei-chungen) zu uberfuhren, werden die Koordinaten uber neue Koordinaten parametrisiert.

Wir betrachten einen speziellen Fall

ut +H(x,Dxu) = 0. (1.4)

Die charakteristischen Gleichungen fur nichtlineare partielle Differentialgleichung ersterOrdnung sind

x = DpH(x, p), (1.5)

p = −DxH(x, p), (1.6)

wobei p := Dxu ist.Die Gleichungen (5),(6) heißen das Hamilton’sche System.

Die charakteristischen Gleichungen ermoglichen uns die gegebenen Anfangsdaten entlangeiner (N-1) dimensionalen Mannigfaltigkeit in eine N dimensionale Losungsmannigfaltigkeitzu entwickeln.

2


Wir mussen dabei sicherstellen, dass die gewohnlichen Differentialgleichungen ohne Pro-blem gelost werden konnen. Zum Beispiel kann die Losung u entlang der charakteristischenGrundkurven existieren, allerdings konnen sich die charakteristischen Grundkurven schnei-den. Oder die globale Losung kann gar nicht existieren.Ein alternativer Zugang ist bei den physikalischen Systemen die gesamte Energie zu mini-mieren (das Prinzip der kleinsten Wirkung).Wir werden im nachsten Unterkapitel dieses Prinzip genauer betrachten. In den weiterenUnterkapiteln werden wir andere typische Anwendungen fur die Hamilton-Jacobi Gleichungprasentieren.

1.1 Mechanik

Das Hamilton’sche Prinzip oder das Prinzip der kleinsten Wirkung: Physikalischen Teilchenund Felder verhalten sich so, dass eine Große, die die Teilchenbahnen und Felder bewertet,kleiner ist als bei allen anderen denkbaren Teilchenbahnen oder Feldwerten.Die Bewertung nennen Physiker die Wirkung, mathematisch ist die Wirkung ein Wert einesFunktionals. Genauer betrachtet erweist sich in vielen Fallen die Wirkung nicht als minimal,sondern nur als stationar. Dann ist das Hamilton’sche Prinzip ein Prinzip der stationarenWirkung.

Die Bahn eines Korpers ist die Bewegung eines mechanischen Systems von einem gegebenenAnfangszustand y0 zu einem gegebenen Endzustand y1 derart, dass die Wirkung S minimalist. S ist ein Funktional, das allen moglichen Bahnen eine reelle Zahl zuordnet. Man versuchtdurch variieren der Bahn jene zu finden, die S minimiert.Sei S das Funktional

S[x(·)] :=

∫ t

0L(x(s), x(s))ds, (1.7)

wobei L : RN × R

N → R die Lagrangefunktion und x(·) eine Trajektorie mit x(0) =y0, x(t) = y1 sind. Die Lagrangefunktion lautet L = T − V , wobei T die kinetische und Vdie potenzielle Energie des betrachteten Systems bezeichnen.

Die Trajektorie x(·), die das Integral (7) minimiert, lost das System von Euler-LagrangeGleichungen

d

ds(∂L

∂q(x(s), x(s))) − ∂L

∂x(x(s), x(s)) = 0, 0 ≤ s ≤ t, q := x(s). (1.8)

Die mit der Lagrangefunktion assoziierte Hamilton-Funktion ist

H(x, p) := p · q(x, p) − L(x,q(x, p)), (1.9)

wobei p := DqL(x, x) und q(x, p) ist die glatte Funktion von x und p, die man erhalt, wennman die Definition der Impulse p = DqL(x, x) nach den Geschwindigkeiten auflost.Dabei x(·) und p(·) genugen dem Hamilton’schen System (1.5), (1.6).

Im nachsten Unterkapitel betrachten wir die Anwendung in der geometrischen Optik.

3


1.2 Geometrische Optik

Strahlenoptik ist eine Naherung der Optik, in der die Welleneigenschaften des Lichts ver-nachlassigt werden, weil die mit dem Licht wechselwirkenden Strukturen (Spiegel, Linsen,Blenden) und die abgebildeten Objektdetails im Verhaltnis zur Wellenlange des Lichts großsind. Sichtbares Licht hat eine Wellenlange von ca. 400 bis 750 nm. Es ist wesentlich kleinerals die Maße von Alltagsgegenstanden und optischen Bauteilen. Man kann die geometrischeOptik mathematisch aus der Wellenoptik als Grenzfall fur verschwindende Wellenlange her-leiten.

Das Ziel der geometrische Optik ist die Untersuchung der Ausbreitung von elektrischen undmagnetischen Feldern im Fall kurzer Wellenlangen.

Analog zur ebenen Welle wird die Phase S durch S(x, t) = ω · t− ~k · x definiert, wobei

ω = St(t, x), (1.10)

~k = −∇S(t, x). (1.11)

Die Phase S hangt also von den zwei Parametern: Wellenvektor ~k und Kreisfrequenz ω ab.Nach der Anwendung der WKB-Methode1 wird unser Gleichungsystem aus der Eikonal-und Transportgleichungen bestehen

(St)2 − c2(∇S)2 = 0, (1.12)

Sttv0 + 2Stv

0t − c2∆Sv0 − 2c2∇ · S · ∇v0 = 0, (1.13)

wobei v0 die Richtung des Energietransportes ist, c ist die Phasengeschwindigkeit. DasSystem beschreibt mittels des Fermatschen Prinzips den kurzesten Weg zwischen zwei durchoptische Medien getrennte Punkte ((1.12) fur die Phase und (1.13) fur den Transport derEnergie) (siehe [4]).Die Losung S genugt der Gleichung

St − c |∇S| = 0. (1.14)

Die Gleichung (1.14) ist die Hamilton-Jacobi Gleichung erster Ordnung.

Als nachstes werden wir die Ausbreitung der Wellen in einem anisotropen Medium betrach-ten (wie in [22]).Die Hamilton-Jacobi Gleichung der Form

H(x, p) := |p|c(x, p|p|) = 1,

wobei p = Sx und c ist die Phasengeschwindigkeit, wird nach der Anwendung der WKB-Methode fur die Elastizitatsgleichungen

ρuitt −∂

∂xj(Cijkl(x)

∂ul

∂xk) = 0, i, j, k, l = 1, 2, 3. (1.15)

1Die WKB(Wentzel–Kramers–Brillouin)-Methode aus der Physik liefert eine Naherung an die Losung. DieNaherung basiert auf der Annahme, dass sich das Potential langsam mit der Position andert und sichdaher eine Losung aus dem konstanten Potential finden lasst. Im Fall der Wellen wird angenommen, dassentweder die Amplitude oder die Phase sich langsam verandert.

4


auftreten.Hier sind ui die Komponenten des Verschiebungsvektors, ρ die Dichte und Ci,j,k,l ein Elas-tizitatstensor.Der WKB-Ansatz fur die Losung ist

uεj(t, x) = εeiS(t,x)/εvε

j (t, x), (1.16)

vεj (t, x) = v0

j (t, x) + εv1j (t, x) + ... (1.17)

mit Anfangsbedingungen

uεj(0, x) = εeiS(0,x)/εφj(0, x), x ∈ R

3, (1.18)

uεjt(0, x) = ψj(0, x), x ∈ R

3. (1.19)

Nach dem Einsetzen von (1.16) in (1.15) erhalt man die nichtstationaren Eikonal- undTransportgleichungen

det

∣∣∣∣∣1

ρCijkl(x)

∂S

∂xj

∂S

∂xk−

(∂S

∂t

)2

I

∣∣∣∣∣ = 0, (1.20)

ρSttv0i +2ρStv

0it−Cijkl(x)

∂2S

∂xj∂xkv0l −2Cijkl(x)

∂S

∂xj

∂v0l

∂xk−

[∂

∂xjCijkl(x)

]∂S

∂xkv0l = 0, (1.21)

fur i = 1, 2, 3.Sei die Phasengeschwindigkeit cα(x, n), α = 1, 2, 3, die Losung des Eigenwertsproblems

det

∣∣∣∣1

ρCijkl(x)njnk − c2αI

∣∣∣∣ = 0 (1.22)

mit der Ausbreitungsrichtung der Welle |n| = 1. Die nichtstationare Eikonal-Gleichung istaquivalent zu den folgenden drei Gleichungen

Sαt − |∇Sα|cα(x,

∇Sα

|∇Sα|

)= 0, α = 1, 2, 3. (1.23)

Die Substitution Sα(t, x) = t+Tα(x) mit der Zeit Tα(x) fuhrt zu den folgenden Gleichungen

|∇Tα(x)|cα(x,

∇Tα(x)

|∇Tα(x)|

)= 1, α = 1, 2, 3. (1.24)

Ein weiteres großes Anwendungsgebiet ist die Steuerungstheorie.

1.3 Steuerungstheorie, Spieltheorie

Wir beginnen mit dem Konzept der Wertefunktion und der dynamischen Programmierung.Dabei beschranken wir uns auf die Spieltheorie mit einem Spieler.Wir betrachten das Zustandssystem mit der Dynamik x = f(x, α), α ist die Steuerungsva-riable und Trajektorien ξ(·; t, x, α(·)), die zur Zeit t im Punkt x starten. Das Kostenfunk-tional(Zielfunktional) fur die Trajektorien ist

C(t, x, α(·)) =

∫ 0

tl(ξ(τ ; t, x, α(·))) dτ + g(ξ(0; t, x, α(·))), (1.25)

5


wobei l, g : RN → R die Kosten charakterisieren, t ∈ [−T, 0], T > 0.

Der Spieler versucht das Funktional zu minimieren.

In der Steuerungstheorie oder in der Spieltheorie wird die Wertefunktion V : [−T, 0]×RN →

R betrachtetV (t, x) = inf

α(·)∈A(t)C(t, x, α(·)). (1.26)

Die Wertefunktion definiert die Kosten fur die optimale Trajektorie.

Die dynamische Programmierung definiert die Vorgehensweise fur die Berechnung der Wer-tefunktion. Die Wertefunktion wird durch die folgende Gleichung

V (t, x) := minα(·)∈A(t)

[ ∫ t+h

tl(ξ(τ ; t, x, α(·))) dτ + V (ξ(t+ h; t, x, α(·)), t + h)

](1.27)

fur −T ≤ t ≤ t+ h ≤ 0 und V (0, x) = g(x) definiert.Wenn die Wertefunktion in einem Punkt differenzierbar ist, dann erfullt sie in diesem Punktdie Isaacs-Bellman Gleichung

DtV (t, x) + minα∈A

DxV (t, x)T f(x, α)︸︷︷︸

H(x,DxV (x,t))

= −l(x). (1.28)

Die Gleichung (1.28) entsteht aus (1.27) mittels Division durch h, Taylor-Entwicklung derFunktion V (ξ, t+ h) und h→ 0.

Die Wertefunktion ist aber nicht uberall differenzierbar und deshalb existiert keine klassischeglobale Losung. Aufgrund dessen wird das Konzept der klassischen Losungen verallgemei-nert. Um die Existenz einer Losung nachzuweisen, benutzt man die Viskositatsmethode

∂uε

∂t− ε · ∇2uε +H(t, x, uε,Dxu

ε) = 0. (1.29)

Als viskose Losung u betrachtet man den Limes von uε gegen Null.M.G. Crandall und P.L. Lions [9, 10, 11] haben das Maximumprinzip verwendet, um dieExistenz der eindeutigen Losung zu beweisen.Spater haben sie das Konzept der viskosen Losungen, das auf der Differentialunglechungenbasiert, vorgeschlagen.Die viskose Losung ist folgendermassen definiert: Sei t = T − τ , dann ist eine beschrankteund stetige Funktion u eine viskose Losung von (2),(3) wenn

1. u = u0 auf (t = 0 × RN )

2. fur jede Testfunktionen φ,ψ ∈ C∞([0, T ] × RN ) gilt

a) wenn (u− φ) das lokale Maximum in (t0, x0) annimmt, dann ist

φt(t0, x0) +H(x0,Dφ(t0, x0)) ≤ 0 (1.30)

b) wenn (u− ψ) das lokale Minimum in (t0, x0) annimmt, dann ist

ψt(t0, x0) +H(x0,Dψ(t0, x0)) ≥ 0 (1.31)

6


Die Wertefunktion fur das Steuerungsproblem ist die eindeutige viskose Losung fur dieHamilton-Jacobi(Bellman-Isaacs) Gleichung.

Als nachstes betrachten wir das Differentialspiel mit zwei Spieler. Sei α die Steuerungsvaria-ble fur den ersten Spieler und β fur den zweiten Spieler, dann ist die Dynamik des Systemsfolgendermaßen darstellbar:

x = f(x, α, β). (1.32)

Beim Konfliktproblem aus der Spieltheorie mit der Zustandsgleichung (1.32) und Neben-bedingungen treten Schwierigkeiten auf. Die Struktur der Wertefunktion ist dann nochkomplizierter und enthalt mehrere Singularitaten, bei dennen ihr Gradient nicht stetig ist.

Weiter wird den Zusammenhang zwischen der Hamilton-Jacobi Gleichungen und Erhal-tungsgleichungen betrachtet.

1.4 Transportgleichungen, Erhaltungsgleichungen

Wir betrachten die lineare Transportgleichung

ut + V (x) · ∇u = 0 (1.33)

und die skalare Erhaltungsgleichung

ut + div(V (x) · u) = 0. (1.34)

Die Gleichungen beschreiben die Verbindung zwischen einer Dichtefunktion u einer physi-kalischen Große und einer zugehorigen Flussfunktion V (RN -wertiges Vektorfeld). Sie re-prasentieren den Spezialfall der Klasse von semiliniaren und nichtlinearen Systemen (hy-perbolische Systeme), die keine klassische Losung besitzen.Wir beschranken uns auf die skalaren linearen Gleichungen, weil (1.33) und (1.34) in dieserDiplomarbeit Mehrphasenstromungen beschreiben.Die Gleichung (1.33) ist exakt die Hamilton-Jacobi Gleichung erster Ordnung.

Zur Losung der Gleichungen verwendet man haufig sogenannte Upwind-Approximationen.Die prinzipielle Idee dieser Methoden ist, die Interpolation abhangig von der Richtungdes Geschwindigkeitsvektors zu machen. Damit nutzt man die Transporteigenschaft vonKonvektions-Diffusions-Problemen aus. Die wichtigsten Verfahren sind das Lax-FriedrichsSchema (fur Konvektionsgleichung, Level-Set Methode) und das Godunov Schema.Das Lax-Friedrichs Schema ist eine numerische Methode zur Losung hyperbolischer partiel-ler Differentialgleichungen, die auf der Finite-Differenzen Methode basiert (FTCS (forwardin time, centered in space) scheme)

H(x, pR, pL) = H(x,pL + pR

2) − 1

2αT (pR − pL), (1.35)

wobei pR und pL die rechte und die linke Approximationen fur die Ableitung p der gesuchtenFunktion u sind:

pLi =

u(x1, ..., xi, ..., xN ) − u(x1, ..., xi − ∆xi, ..., xN )

∆xi,

7


pRi =

u(x1, ..., xi + ∆xi, ..., xN ) − u(x1, ..., xi, ..., xN )

∆xi

und αi = maxp∈I

∣∣∣∂H∂pi

∣∣∣.

Das Schema von Godunov ist in der numerischen Analyse ein konservatives numerischesSchema, das von S. K. Godunov 1959 eingefuhrt wurde, um die partiellen Differentialglei-chungen zu losen. In dieser Methode werden die konservativen Variablen als stuckweise Kon-stanten betrachtet und dem Verfahren wird die Finite-Volumen-Methode zugrunde gelegt.Das Verfahren benutzt in seiner Herleitung eine integrale Form der Erhaltungsgleichungenund erlaubt damit auch unstetige Losungen, die fur solche Gleichungen typisch sind. Fernerwerden nur geringe Anforderungen an die Gitterzellen gestellt, was unstrukturierte und fle-xible Geometrien erlaubt. Daruber hinaus werden die konservativen Großen der Gleichungtatsachlich erhalten.

Ein anderes Upwind-Verfahren wird im Kapitel 3 praziser betrachtet.

1.5 Zielsetzung

Die Hamilton-Jacobi Gleichungen wurden bis zum heutigen Tag sehr intensiv untersucht.Die verschiedenen Methoden wurden fur spezielle Steuerungsprobleme entwickelt (Sougani-dis [25, 26], Krasovskii [20], Subbotin [27, 28, 29], Tarasyev [30], Isaacs [19], Mitchell [24],Falcone [12], Bardi [1, 2], Soravia [3], Botkin [5, 6, 22]).Schwierigkeiten treten meistens auf, wenn Nebenbedingungen betrachtet werden. In derSpieltheorie wird das Kostenfunktional komplexer und die Bedingungen, die die viskoseLosungen festlegen, fur die Probleme mit Zustandsbeschrankungen schwer zu formulieren.

Ziel der vorliegenden Arbeit ist eine Formulierung der Bedingungen fur die viskose Losungeneines Konflikt-Steuerunsproblems mit Zustandsbeschrankungen, im Kapitel 2 werden wirdiese Bedingungen prasentieren.Der Beweis der aufgestellten Behauptung und die Entwicklung eines stabilen numerischenVerfahrens zur Losung der Hamilton-Jacobi Gleichungen wird im Kapitel 3 betrachtet.Außerdem wird im Kapitel 4 und im Kapitel 5 die Anwendung des Vefahrens auf die Proble-me der aktuellen forschungsbezogenen Gebiete, wie die Optimierung von Kuhlprotokollenbei der Kryopreservation von lebenden Zellen sowie Simulation von Mehrphasenstromungenprasentiert.

8

2 Hamilton-Jacobi-Gleichung aus derKonflikt-Steuerungstheorie mitZustandsbeschrankungen

2.1 Kostenfunktional und Wertefunktion

Die Differentialspieltheorie benutzt Methoden der klassischen Spieltheorie sowie der op-timalen Steuerungstheorie. Die entwickelte Theorie der viskosen Losungen und die ent-sprechenden numerischen Methoden konnen zur Losung des Differentialspiels angewendetwerden. Bisher wurden in vielen wissenschaftlichen Arbeiten nur die numerischen Methodenfur die Losung der Hamilton-Jacobi(Bellman-Isaacs) Gleichung vorgeschlagen, die aus demDifferentialspiel mit dem folgenden Funktional entsteht

γ0(x(·)) := mint∈[t0,T ]

χ(t, x(t)), (2.1)

wobei t0 die Anfangszeit ist und T ist die Endzeit, x(·) Trajektorie des Kontrollsystems undχ die gegebene Funktion.Die Hamilton-Jacobi (Bellman-Isaacs) Gleichung kommt ins Spiel im Zusammenhang mitder viskosen Losungen. Es wurde in [30] gezeigt, dass die “minmax“-Losungen der Bellman-Isaacs Gleichung mit der viskosen Losungen der Hamilton-Jacobi Gleichung ubereinstimmen.Allerdings bereitet die Erweiterung der Formulierung des Kostenfunktionals auf die Proble-me mit Zustandsbeschrankungen die Schwierigkeiten bei der Losung. Wir werden hier einspezielles Kostenfunktional definieren. Dieses Funktional ermoglicht uns solche Problemezu losen.

In dieser Diplomarbeit wird die Bellman-Isaacs Gleichung betrachtet, die beim Differenti-alspiel mit dem folgenden Funktional auftritt

γ(x(·)) := max mint∈[t0,T ]

χ(t, x(t)), maxt∈[t0 ,T ]

θ(t, x(t)), (2.2)

wobei χ : [0, T ] × RN → R, θ : [0, T ] × R

N → R gegebene Funktionen sind und

θ(t, x(t)) < χ(t, x(t), t ∈ [0, T ].

Das Differentialspiel startet zur Zeit t0 ∈ [0, T ] und wird beendet zur Zeit T . Das Kosten-funktional γ ist auf den Trajektorien x(·) definiert mit dem dynamischen System

x = f(t, x, α, β), t ∈ [0, T ], x ∈ RN , α ∈ A ⊂ R

a, β ∈ B ⊂ Rb, (2.3)

wobei x ein Zustandsvektor ist, α, β sind die entsprechenden Steuerungsvariablen fur denersten und den zweiten Spieler, A und B sind gegebene kompakte Mengen.

9

2 Hamilton-Jacobi-Gleichung aus der Konflikt-Steuerungstheorie mit ZB

Der erste Teil des Funktionals (2.2) mint∈[t0,T ] χ(t, x(t)) wird verantwortlich fur die Qualitatder Steuerung und der zweite Teil maxt∈[t0,T ] θ(t, x(t)) fur die Zustandsbeschrankungen. Daskann man sehen, wenn man ein Differentialspiel mit der Zielmenge

M := (t, x) : t ∈ [t0, T ], χ(t, x) ≤ 1

und der Zustandsbeschrankungsmenge

N := (t, x) : t ∈ [t0, T ], θ(t, x) ≤ 1

betrachtet. Offensichtlich ist es, dass M in N liegt (M ⊂ N ).

Wenn wir die Wertefunktion fur das Differentialspiel (die genaue Definition der Wertefunk-tion wird auf der nachsten Seite eingefuhrt) kleiner oder gleich 1 in der Startposition ist,dann existiert eine Strategie fur den ersten Spieler, dass fur alle Strategien fur den zweitenSpieler und alle Trajektorien x(·) die folgenden Behauptungen gelten:

• mint∈[t0,T ] χ(t, x(t)) ≤ 1, d.h. die Trajektorie x(t) erreicht die Zielmenge zu einemZeitpunkt t ≤ T ,

• maxt∈[t0,T ] θ(t, x(t)) ≤ 1, d.h. die Trajektorie x(t) bleibt in der Zustandsbeschrankungs-menge fur alle t ∈ [t0, T ].

Mit anderen Worten erlaubt uns das modifizierte Funktional (2.2) die Konflikt-Steuerungs-probleme mit Zustandsbeschrankungen zu behandeln.

Der erste Spieler versucht dabei mit der Steuerungsvariable α das Kostenfunktional (2.2)zu minimieren. Das Ziel des zweiten Spielers ist das Gegenteil mit der Steuerungsvariableβ zu erreichen und wir konnen dann die Resultate aus [20, 29] benutzen.Die Spieler des Differentialspiels verwenden entgegengesetzte Strategien (Ruckkopplungs-strategien), die die beliebigen Funktionen sind

A : [0, T ] × RN → A, B : [0, T ] × R

N → B.

Fur einen beliebigen Anfangspunkt (t0, x0) ∈ [0, T ]×RN und beliebige Strategien A,B sind

zwei Funktional-Mengen X1(t0, x0,A) and X2(t0, x0,B) definiert. Die Mengen enthalten dieGrenzwerte der Schritt-fur-Schritt Losungen von (2.3), die mit Strategien A und B generiertwerden (siehe [20, 29]).

Unsere Voraussetzungen an die Funktionen f, χ, θ sind:

(f1) Die Funktion f ist gleichmaßig stetig in [0, T ] × RN ×A×B.

(f2) f ist beschrankt1, d.h.|f(t, x, α, β)| ≤M,

fur alle (t, x, α, β) ∈ [0, T ] × RN ×A×B.

1Die Beschrankheit vereinfacht unsere Funktion. Im Allgemein wird nur die Monotonie-Eigenschaft ver-langt.

10


(f3) f ist Lipschitz-stetig in t, x, d.h.

|f(t1, x1, α, β) − f(t2, x2, α, β)| ≤ N(|t1 − t2| + |x1 − x2|),

fur alle (ti, xi, α, β) ∈ [0, T ] × RN ×A×B, i = 1, 2.

(f4) χ ist beschrankt und Lipschitz-stetig in t, x, d.h.

|χ(t, x)| ≤ C00,

und|χ(t1, x1) − χ(t2, x2)| ≤ L00(|t1 − t2| + |x1 − x2|),

fur alle (t, x), (ti, xi) ∈ [0, T ] × RN , i = 1, 2.

θ ist beschrankt und Lipschitz-stetig in t, x, d.h.

|θ(t, x)| ≤ C01,

und|θ(t1, x1) − θ(t2, x2)| ≤ L01(|t1 − t2| + |x1 − x2|),

fur alle (t, x), (ti, xi) ∈ [0, T ] × RN , i = 1, 2.

(f5) Die Funktion f genugt der Sattelpunkt Bedingung:

minα∈A

maxβ∈B

〈p, f(t, x, α, β)〉 = maxβ∈B

minα∈A

〈p, f(t, x, α, β)〉,

fur jede p ∈ RN , (t, x) ∈ [0, T ] × R

N .

Definition ([20, 29]). Unter den Vorraussetzungen (f1)-(f5), fur das Differentialspiel (2.2)-(2.3) ist die Wertefunktion c : (t, x) → c(t, x) folgendermaßen definiert

c(t, x) := minA

maxx(·)∈X1(t,x,A)

γ(x(·)) = maxB

minx(·)∈X2(t,x,B)

γ(x(·)). (2.4)

Die Wertefunktion ist beschrankt und Lipschitz-stetig in t, x, d.h.

|c(t, x)| ≤ C,

|c(t1, x1) − c(t2, x2)| ≤ L(|t1 − t2| + |x1 − x2|),fur jede (t, x), (ti, xi) ∈ [0, T ] × R

N , i = 1, 2 (siehe [20]).

Bemerkung. Es gilt

θ(t, x) ≤ c(t, x) ≤ χ(t, x), ∀(t, x) ∈ [0, T ] × Rn. (2.5)

Die erste Ungleichung folgt aus der Definition des Funktionals γ. Die zweite Ungleichungist nicht offensichtlich. Diese Ungleichung wird spater bewiesen, wenn wir die Approxima-tionfunktion fur die Wertefunktion c konstruieren, die der Ungleichung genugt. Nach derKonstruktion zeigen wir, dass die Approximationsfunktion punktweise gegen der Funktionc konvergiert.

11


2.2 Formulierung der Bedingungen fur die viskose Losungen

In [29] wurde gezeigt, dass einige Differentialungleichungen unter den gewissen Vorausset-zungen die Wertefunktion c fur das Differentialspiel (2.2) mit dem Kostenfunktional (2.1) be-stimmen. Dabei wurde die Gateaux-Differenzierbarkeit der Wertefunktion verlangt. In [28]wurde diese Forderung relaxiert und die Differentialungleichungen wurden durch die obereund untere Richtungsableitungen definiert. Spater wurde in [30] gezeigt, dass die Unglei-chungen fur die obere und untere Richtungsableitungen aquivalent zu den Ungleichungen,die die viskosen Losungen definieren, sind (vgl. die Differentialungleichungen (1.30),(1.31)).Die Aquivalenz ist lokal: Wenn die Ungleichungen fur die Richtungsableitungen in einer klei-nen Umgebung eines Punktes erfullt sind, dann sind die entsprechenden Ungleichungen furdie viskose Losung im gleichen Punkt erfullt. Insbesondere gilt, dass alle diese Argumenteauf das Differentialspiel (2.2) mit dem Zielfunktional (2.3) leicht verallgemeinern lassen. Zudiesem Zweck ist es genugend zu zeigen, dass die v- Stabilitatseigenschaft der Wertefunk-tion sollte nur in Positionen (t, x) uberpruft werden, wo die Ungleichung c(t, x) > θ(t, x)gilt. Beweis ist dabei dieselbe wie in Lemma 3.6.2 von [?]. Entsprechend gilt die nachsteBehauptung.

Satz 2.2.1 ([9, 28, 30])Eine stetige Funktion c : [0, T ] × R

N → R ist die Wertefunktion des Differentialspiels (2.2)- (2.3), wenn die Funktionen

u(τ, x) := c(T − τ, x), u0(τ, x) := χ(T − τ, x), u1(τ, x) := θ(T − τ, x)

den folgenden Ungleichungen genugen:

a) Fur jede (τ, x) ∈ [0, T ] × RN

u(τ, x) ≥ u1(τ, x), u(0, x) = u0(0, x). (2.6)

b) Fur jeden Punkt (s0, y0) ∈ (0, T ) × RN mit u(s0, y0) > u1(s0, y0) und jede Testfunktion

ϕ ∈ C1((0, T ) × RN fur die Differenz u − ϕ das lokale Maximum in (s0, y0) annimmt gilt

die folgende Ungleichung:

∂ϕ

∂τ(s0, y0) +H(s0, y0,Dϕ(s0, y0)) ≤ 0. (2.7)

c) Fur jeden Punkt (s0, y0) ∈ (0, T ) × RN) mit u(s0, y0) < u0(s0, y0) und jede Testfunktion

ϕ ∈ C1((0, T ) × RN ) fur die Differenz u − ϕ das lokale Minimum in (s0, y0) annimmt gilt

die folgende Ungleichung:

∂ϕ

∂τ(s0, y0) +H(s0, y0,Dϕ(s0, y0)) ≥ 0. (2.8)

Hier istH(s, y, p) = −max

β∈Bminα∈A

〈p, f(T − s, y, α, β)〉

die Hamilton-Funktion des Differentialspiels.

Als nachstes betrachten wir die Entwicklung von stabilen Verfahren zur Berechnung derFunktion u, die der Ungleichungen des Satzes 2.2.1 genugt.

12

3 Ein stabiles Verfahren zur Losung vonHamilton-Jacobi-Gleichung aus derKonflikt-Steuerungstheorie mitZustandsbeschrankungen

Die Differentialungleichungen aus dem Satz 2.2.1 fur das Differentialspiel (2.2) - (2.3) de-finieren die viskosen Losungen der Hamilton-Jacobi (Bellman-Isaacs) Gleichung, d.h. dieWertefunktion (2.4). Wir werden jetzt einen abstrakten Operator betrachten, der die Ap-proximationen fur die viskosen Losungen bzw. Wertefunktion liefert.

3.1 Abstrakter Zeit-Schritt-Operator

Sei u(τ, x) := c(T − τ, x), u0(τ, x) := χ(T − τ, x), u1(τ, x) := θ(T − τ, x). Dann geltenfur u, u0, u1 die Behauptungen a) - c) des Satzes 2.2.1 und

|u0(t1, x1) − u0(t2, x2)| ≤ L00(|t1 − t2| + |x1 − x2|), |u0(t, x)| ≤ C00,

|u1(t1, x1) − u1(t2, x2)| ≤ L01(|t1 − t2| + |x1 − x2|), |u1(t, x)| ≤ C01,

|u(t1, x1) − u(t2, x2)| ≤ L(|t1 − t2| + |x1 − x2|), |u(t, x)| ≤ C

fur jede (t, x), (ti, xi) ∈ [0, T ] × RN , i = 1, 2.

Wir fuhren einen Operator (siehe [25])

F (τ, ρ, ·) : C0,1b (RN ) → C0,1

b (RN ), 0 ≤ τ ≤ T, ρ ≥ 0

ein, der die folgenden Eigenschaften hat:

(F1) F (τ, 0, v) = v, v ∈ C0,1b .

(F2) Die Abbildung (τ, ρ) → F (τ, ρ, v) ist stetig bezuglich der Maximumsnorm ‖ · ‖.

(F3) F (τ, ρ, v + k) = F (τ, ρ, v) + k fur jede k ∈ R.

(F4) ‖F (τ, ρ, v)− v‖ ≤ C1ρ, wobei C1 kann abhangig von ‖v‖ und ‖Dv‖ sein, ‖Dv‖ ist dieLipschitz’sche Konstante.

(F5) Es existieren ein r > 0 und ein L1 > 0, so:wenn

v1(x) ≤ v2(x) fur jede x ∈ RN

ist, dann gilt fur jede y ∈ RN mit

|v1(y + w) − v1(y + w)|, |v2(y +w) − v2(y + w)| ≤ L|w − w|, ∀w, w ∈ Bn(0, ρr)

13

3 Ein stabiles Verfahren zur Losung von Hamilton-Jacobi-Gleichung

F (τ, ρ, v1)(y) ≤ F (τ, ρ, v2)(y).

Hier ist L := max2L+ 1, L1.

(F6) Es existiert eine Konstante C2 > 0, so dass

‖F (τ, ρ, v)‖ ≤ eρC2(‖v‖ + ρC2),

d.h. ‖Dv‖ ≤ L.

(F7) Es existieren Konstanten C3, C4 > 0, so dass

eT (C3+C4)(L00 + TC4) ≤ L

und‖DF (τ, ρ, v)‖ ≤ eρ(C3+C4)(‖Dv‖ + ρC4),

d.h. ‖v‖ ≤ eTC2(C00 + TC2) und ‖Dv‖ ≤ L.

(F8) Fur jede φ ∈ C2b (RN ) und x ∈ R

N , so dass |Dφ(x)| < 2L+ 1, gilt

∣∣∣∣F (τ, ρ, φ)(x) − φ(x)

ρ+H(τ, x,Dφ(x))

∣∣∣∣ ≤ C5 · (1 + ‖Dφ‖ + ‖D2φ‖)ρ.

3.2 Zwei Beispiele des Zeit-Schritt-Operators

Eine Klasse der Schritt-Operatoren, die die Eigenschaften (F1) - (F8) besitzen, ist in [11, 25]beschrieben. Seien ai positive Konstanten und ∆xi = aiρ raumliche Schrittweiten, dieproportional zu den zeitlichen Schrittweiten ρ, i = 1, n sind.Der Operator

F (τ, ρ, v)(x) =

v(x) − ρH(τ, x,v(x1, ..., xn) − v(x1 − ∆x1, ..., xn)

∆x1, ...,

v(x1, ..., xn) − v(x1, ..., xn − ∆xn)

∆xn),

erfullt alle Eigenschaften, außer der Monotonie (F5).Wenn die Hamilton-Funktion H(τ, x, p1, p2, ..., pn) in den Impuls-Variablen p1, p2, ..., pn mo-noton wachst und ∆xi/ρ ≡ ai ≥ √

nΩ, wobei Ω die Lipschitz-Konstante fur H in p ist,dann gilt die Monotonie-Eigenschaft (F5).Die geforderte Monotonie der Hamilton-Funktion kann mit der Transformation

x = x− C(T − t)

in Kontrollsystem erreicht werden, wenn C ≥ Ω konstant ist. Eine solche Transformationandert die Lipschitz-Konstante fur die Hamilton-Funktion in den Impuls-Variablen.Die Bedingung C ≥ Ω′, wobei Ω′ die Lipschitz-Konstante ist, kann verletzt werden. Das istaber vermeidbar mit der gewunschten Monotonie.Wir sollten auch erwahnen, dass die Technik nur anwendbar ist, wenn T klein ist. Das istfast unmoglich fur realistische Probleme.

14


Wir betrachten deshalb einen anderen Operator, der in [21] prasentiert wurde. Er erfulltdie Monotonie-Eigenschaft (F5) und verlangt keine Variablen-Transformationen.Der Operator ist durch

F (τ, ρ, v)(x) = v(x) + ρmaxβ∈B

minα∈A

n∑

i=1

(pLi · f+

i + pRi · f−i ) (3.1)

definiert, wobeifi = fi(T − τ, x1, ..., xn, α, β),

pLi =

v(x1, ..., xi, ..., xn) − v(x1, ..., xi − ∆xi, ..., xn)

∆xi,

pRi =

v(x1, ..., xi + ∆xi, ..., xn) − v(x1, ..., xi, ..., xn)

∆xi,

∆xi = aiρ, f+i = max (fi, 0) und f−i = min (fi, 0), i = 1, n sind.

Bemerkung. 1. Sei M die Beschrankung fur die rechte Seite des Kontroll-Systems. Wenn∆xi/ρ ≡ ai ≥

√n·M, i = 1, n, dann ist der durch (3.1) gegebene Operator F (τ, ρ) monoton

(siehe [6]).

2. Wir halten fest: Wenn ρ fixiert ist, dann konnen die beiden Schritt-Operatoren auf diespeziellen Funktionen beschrankt sein. Die Funktionen sind auf einem rechteckigen Gittermit der Schrittweite ∆xi in der ite Koordinate definiert, i = 1, n. Diese Operatoren erge-ben sich auch aus diskreten Differenzen Schematas, wenn die betrachtende Approximationangewandt wird.

3. Wenn die Funktion f linear in α und β fur alle feste (t, x) ist, dann kann die Operationmaxβ∈B

minα∈A

in der Definition des Operators (3.1) durch maxβ∈ext B

minα∈ext A

ersetzt werden, wobei

“ext” die Menge der extremalen Punkten ist. Die Operation “maxmin” kann uber die Men-gen der Ecken berechnet werden, wenn A und B Polyeder sind (siehe [6]). Diese Bemerkungist sehr wichtig fur die numerische Implementierung des Operators (3.1), weil die Operation“max min” auf die Funktionen, die nicht linear und nicht konvex (konkave) in α und β sind,angewendet wird.

3.3 Approximation der Wertefunktion

Fur die Zeitzerlegung P = 0 = t0 < t1 < ... < tm(P ) = T sei die Funktion

uP : [0, T ] × RN → R

mit

uP (0, x) := maxu0(0, x), u1(0, x) = u0(0, x),

uP (τ, x) := max

min u0(τ, x), F (τ, τ − ti−1, uP (ti−1, ·))(x), u1(τ, x)

(3.2)

definiert, wobei τ ∈ (ti−1, ti] fur jedes i = 1, 2, ...,m(P ) ist.

15


Lemma 3.3.1

(a) Fur alle τ ∈ [0, T ] gilt‖uP (τ, ·)‖ ≤ eτC2(C00 + τC2),

und uP (τ, ·) ∈ C0,1b (RN ) mit

‖DuP (τ, ·)‖ ≤ eτ(C3+C4)(L00 + τC4).

Und wenn ti > τ ≥ ti−1 fur jedes i = 1, ...,m(P ) ist, dann gilt

‖uP (τ, ·) − uP (ti−1, ·)‖ ≤ maxL00, C1(τ − ti−1),

wobei C1 = C1(eTC2(C00 + TC2), L) ist.

(b) Die Funktion uP is beschrankt und gleichmaßig stetig auf [0, T ] × RN .

(c) Fur jede (τ, x) ∈ [0, T ] × RN genugt die Funktion uP (τ, x) der folgenden Ungleichung

u1(τ, x) ≤ uP (τ, x) ≤ u0(τ, x). (3.3)

Der Beweis basiert auf den Operator-Eigenschaften und der Definition der Funktion uP .

3.4 Punktweise Konvergenz der Approximation

Satz 3.4.1

Unter den gemachten Vorraussetzungen (f1)-(f5) und (F1)-(F8) gilt die folgende Behaup-tung

|u(τ, x) − uP (τ, x)| ≤ K|P |1/2, (τ, x) ∈ [0, T ] × RN .

Die Konstante K ist nur von Konstanten M, N, C00, L00, C01, L01 abhangig, die die Funk-tionen f, u0 und u1 charakterisieren, und von Konstanten, die in den Eigenschaften (F1) -(F8) enthalten sind. |P | ist die großte Schrittweite in der Zeitzerlegung.

Beweis. Wir verwenden die Resultate aus [5] und [25]. Die meisten Argumentationen sindim Beweis des Satzes 3.4.1 parallel zu den Argumentationen in [5] und [25].

Es genugt zu zeigen dass eine Konstante K1 existiert, so dass

MP = sup(τ,x)∈[0,T ]×RN

(uP (τ, x) − u(τ, x))± ≤ K1|P |1/2,

wobei die Operatoren “+”,“–” bedeuten: r+ = max(r, 0), r− = max(−r, 0). Wir haltenfest: Wenn wir “±” schreiben, meinen wir, dass die Falle “+” und “–” separat betrachtetwerden sollten.Ohne Beschrankung der Allgemeinheit (o.B.d.A.) setzen wir voraus, dass

MP > 0

16


und mit dem Lemma 3.3.1 gilt:Es existiert eine Konstante ℜ1, die unabhangig von der Zerlegung P ist, so dass gilt

|uP (τ, x)| ≤ ℜ1, ∀(τ, x) ∈ [0, T ] × RN ,ℜ1 ist die Konstante.

Mit der Funktion ℜ := max(ℜ1, C) und ε := |P |1/4 betrachten wir die Funktion Φ, diefolgendermaßen definiert ist

Φ : RN × R

N × [0, T ] × [0, T ] → R

Φ(x, y, τ, s) = (uP (τ, x) − u(s, y))± + 4(ℜ + 1)βε(x− y) + 4(ℜ + 1)γε(τ − s) − τ + s

4TMP .

Hier sind βε(·) = β(·/ε) und γε(·) = γ(·/ε), wobei β, γ Funktionen mit folgenden Eigen-schaften sind:

β ∈ C20 (RN ), β(0) = 1, 0 ≤ β ≤ 1,

‖Dβ‖ ≤ 2, ‖D2β‖ ≤ 15, β(w) = 0 fur |w| > 1,

β(w) = 1 − |w|2 fur |w| ≤√

3

4,

β(w) <1

4fur |w| >

√3

4,

γ ∈ C20 (R), γ(0) = 1, 0 ≤ γ ≤ 1,

|Dγ| ≤ 2, |D2γ| ≤ 15, γ(t) = 0 fur |t| > 1,

γ(t) = 1 − t2 fur |t| ≤√

3

4,

γ(t) <1

4fur |t| >

√3

4.

Wir berucksichtigen dabei, dass die Beschrankungen fur die erste und zweite Ableitungenvon β und γ kompatibel mit den anderen Eigenschaften der Funktionen sind.Weiter betrachten wir eine

”Permutation-“ Funktion

Ψ(x, y, τ, s) = Φ(x, y, τ, s) + 2δξδ(x, y),

wobei δ > 0 ein kleiner Parameter ist. Dabei sind ξδ ∈ C20(RN ×R

N ) Funktionen, so dass dieFunktion Ψ ihren Maximalwert auf R

N ×RN × [0, T ]× [0, T ] annimmt. Eine solche Funktion

existiert in der Tat. Da Φ beschrankt auf RN × R

N × [0, T ] × [0, T ] fur jede δ > 0 ist, dannexistiert einen Punkt (x1, y1, τ1, s1) ∈ R

N × RN × [0, T ] × [0, T ], so dass

Φ(x1, y1, τ1, s1) > supΦ(x, y, τ, s) : (x, y, τ, s) ∈ RN × R

N × [0, T ] × [0, T ] − δ.

Wir wahlen ξδ ∈ C∞0 (RN × R

N ) mit:

0 ≤ ξδ ≤ 1, ξδ(x1, y1) = 1, ‖Dξδ‖ ≤ 1, ‖D2ξδ‖ ≤ 1.

17


Offensichtlich ist Ψ = Φ mit dem Trager außerhalb ξδ und

Ψ(x1, y1, τ1, s1) = Φ(x1, y1, τ1, s1) + 2δ

> supΦ(x, y, τ, s) : (x, y, τ, s) ∈ RN × R

N × [0, T ] × [0, T ] + δ

> supΨ(x, y, τ, s) : (x, y, τ, s) ∈ (RN × RN \ supp(ξδ)) × [0, T ] × [0, T ].

Es existiert (x0, y0, τ0, s0) ∈ RN × R

N × [0, T ] × [0, T ], so dass

Ψ(x0, y0, τ0, s0) ≥ Ψ(x, y, τ, s) (3.4)

fur alle (x, y, τ, s) ∈ RN × R

N × [0, T ] × [0, T ] gilt.Obwohl der Maximalpunkt abhangig von δ ist, werden die nachsten Betrachtungen dieseTatsache nicht benotigen. Und der Index δ wird weggelassen.

Die Idee ist es, eine Funktion mit einer doppelten Menge von Variablen zu betrachten. DieFunktion Ψ in unserem Fall wurde in [9] vorgeschlagen. Mit der Funktion wird den Unter-schied zwischen zwei viskosen Losungen abgeschatzt. Diese Idee wurde in [25] verwendet,um den Fehler fur das finite Differenzen Schema abzuschatzen. Mit der Funktion Ψ wird eineFunktion konstruiert, die der Ungleichung (2.8) aus dem Satz (2.2.1) in dem Punkt(s0, y0)genugt, und eine andere Funktion, die der Ungleichung (2.7) aus dem Satz (2.2.1) approxi-mativ in (τ0, x0) genugt. Die zwei Funktionen geben uns die gewunschte Abschatzung fur(uP − u)+. Die Abschatzung (uP − u)− wird mit der gleichen Vorgehensweise erreicht.

Um den Beweis weiter fortzusetzen sollten wir den Abstand zwischen den Punkten (τ0, x0)und (s0, y0) abschatzen.

Lemma 3.4.2

Fur δ < min34 ,

18MP gelten die folgenden Ungleichungen:

|x0 − y0| ≤√

3/4ε, |τ0 − s0| ≤√

3/4ε,

|x0 − y0| ≤L+ 2δ

4(ℜ + 1)ε2, |τ0 − s0| ≤

L+ ℜ/T4(ℜ + 1)

ε2, (3.5)

0 <1

4MP ≤ (uP (τ0, x0) − u(s0, y0))

±. (3.6)

Dieses Lemma erlaubt und den Beweis fortzusetzen. Der Beweis vom Lemma 3.4.2 wirdnach dem Beweis des Satzes 3.4.1 erklart.

1. Zuerst betrachten wir den Fall “+”, so dass gilt


(uP (τ, x) − u(τ, x))+.

Wir fangen mit dem Fall τ0 = 0, s0 ≥ 0 an.

Mit dem Lemma 3.4.2 haben wir

0 <1

4MP ≤ (uP (0, x0) − u(s0, y0))

+

= uP (0, x0) − u(s0, y0) = u0(0, x0) − u(0, y0) + u(0, y0) − u(s0, y0)

18


≤ L00|x0 − y0| + L|τ0 − s0| ≤(L00

L+ 4δ

4(ℜ + 1)+ L

L+ ℜ/T4(ℜ + 1)

)|P |1/2.

Die Konstanten C,L sind Funktionen, die von M,N,C00, C01, L00, L01, T abhangig sind(siehe [29]). ℜ ist die Funktion mit M,N,C00, C01, L00, L01, T und die Behauptung desSatzes 3.4.1 ist wahr in dem Fall.

Jetzt betrachten wir den Hauptfall τ0 > 0, s0 ≥ 0.

(a) Wenn u(s0, y0) > u0(s0, y0) ist, dann erhalten wir

0 <1

4MP ≤ uP (τ0, x0) − u(s0, y0) ≤ uP (τ0, x0) − u0(s0, y0)

= uP (τ0, x0) − u0(τ0, x0) + u0(τ0, x0) − u0(s0, y0).

Das Lemma 3.3.1 gibt uns die Ungleichung uP (τ0, x0) ≤ u0(τ0, x0), dann ergibt sich

0 < uP (τ0, x0) − u(s0, y0) ≤ u0(τ0, x0) − u0(s0, y0),

und mit dem Lemma 3.4.2 erhalten wir die gewunschte Abschatzung

MP ≤ 4L00

(L+ 4δ

4(ℜ + 1)+L+ ℜ/T4(ℜ + 1)

)|P |1/2.

Die Behauptung des Satzes 3.4.1 ist damit fur den Fall (a) bewiesen. Dabei werden dieBemerkungen fur C,L und ℜ berucksichtigt.

(b) Wenn uP (τ0, x0) = u1(τ0, x0) ist, erhalten wir von (3.6) und (3.2)

u1(τ0, x0) = uP (τ0, x0) > u(s0, y0) ≥ u1(s0, y0).

Dann gilt0 > u(s0, y0) − uP (τ0, x0) ≥ u1(s0, y0) − u1(τ0, x0)

und mit dem Lemma 3.4.2 bekommen wir die gewunschte Abschatzung

MP ≤ 4L00

(L+ 4δ

4(ℜ + 1)+L+ ℜ/T4(ℜ + 1)

)|P |1/2.

Die Behauptung des Satzes 3.4.1 ist damit fur den Fall (b) bewiesen. Dabei werden dieBemerkungen fur C,L und ℜ berucksichtigt.

(c) Wenn u(s0, y0) = u0(s0, y0) ist, dann erhalten wir von (3.6) und (3.2)

u0(s0, y0) = u(s0, y0) < uP (τ0, x0) ≤ u0(τ0, x0).

Dann gilt0 < uP (τ0, x0) − u(s0, y0) ≤ u0(τ0, x0) − u0(s0, y0)

und mit dem Lemma 3.4.2 bekommen wir die gewunschte Abschatzung

MP ≤ 4L00

(L+ 4δ

4(ℜ + 1)+L+ ℜ/T4(ℜ + 1)

)|P |1/2.

19


Die Behauptung des Satzes 3.4.1 ist damit fur den Fall (c) bewiesen. Dabei werden dieBemerkungen fur C,L und ℜ berucksichtigt.

(d) Wir nehmen jetzt an: u(s0, y0) < u0(s0, y0) und uP (τ0, x0) > u1(τ0, x0), die insbesonderes0 > 0 bedeutet.

Schritt 1. Wir betrachten eine Testfunktion ϕ

ϕ(s, y) = uP (τ0, x0) + 4(ℜ + 1)βε(x0 − y) + 4(ℜ + 1)γε(τ0 − s) + 2δξδ(x0, y) −τ0 + s

4TMP .

Es ist klar, dass u(s, y)−ϕ(s, y) ≥ −Ψ(x0, y, τ0, s) und dann nimmt u−ϕ das lokale Minimumim Punkt (s0, y0) an. Die Substitution fur ϕ in (2.8) bedeutet dass die Ungleichung

0 ≤ −4(ℜ + 1)γ′ε(τ0 − s0) −1

4TMP +H(s0, y0,Dϕ(s0, y0)).

Daraus folgt1

4TMP ≤ −4(ℜ + 1)γ′ε(τ0 − s0) +H(s0, y0,Dϕ(s0, y0)). (3.7)

Schritt 2. Wir betrachten jetzt eine andere Testfunktion ψ,

ψ(τ, x) = u(s0, y0) − 4(ℜ + 1)βε(x− y0) − 4(ℜ + 1)γε(τ − s0) − 2δξδ(x, y0) +τ + s0

4TMP .

Wir brauchen das folgende Lemma, um die Eigenschaft (F5) des Operators F mit denFunktionen v1(·) = uP (ti−1, ·) und v2(·) = ψ(ti−1, ·) + k (k ist eine Konstante) im Punkt x0

anzuwenden.

Lemma 3.4.3

Sei r die Konstante in (F5) und ti−1 ist so, dass τ0 ∈ (ti−1, ti]. Dann gelten fur δ < 113 und

|P |1/2 < min1, 1/(240(ℜ + 1)r) die Ungleichungen

|ψ(ti−1, x0 + w) − ψ(ti−1, x0 + w)| < L|w − w|,

|uP (ti−1, x0 + w) − uP (ti−1, x0 + w)| < L|w − w|,fur jede w, w ∈ Bn(0, (τ0 − ti−1)r).

Der Beweis vom Lemma 3.4.3 wird nach dem Beweis des Satzes prasentiert.

Bevor wir das Lemma 3.4.3 anwenden, bemerken wir, dass uP (τ, x)−ψ(τ, x) ≤ Ψ(x, y0, τ, s0)ist und uP − ψ nimmt das lokale Maximum im Punkt (τ0, x0) an. Dann gilt fur ti−1 mitτ0 ∈ (ti−1, ti] die Ungleichung

(uP − ψ)|(ti−1,x) ≤ (uP − ψ)|(τ0 ,x0)

fur jede x ∈ RN . Mit anderen Wortern

uP (ti−1, x) ≤ ψ(ti−1, x) + (uP (τ0, x0) − ψ(τ0, x0)). (3.8)

20


Mit der Definition (3.2) und der Voraussetzung uP (τ0, x0) > u1(τ0, x0) gilt

uP (τ0, x0) = minu0(τ0, x0), F (τ0, τ0 − ti−1, uP (ti−1, ·))(x0).

Wir beachten jetzt das Lemma 3.4.3, Eigenschaften (F3) und (F5) des Operators F undUngleichung (3.8) und damit erhalten wir

uP (τ0, x0) ≤ F (τ0, τ0 − ti−1, uP (ti−1, ·))(x0)

≤ F (τ0, τ0 − ti−1, ψ(ti−1, ·))(x0) + (uP (τ0, x0) − ψ(τ0, x0)). (3.9)

Nach der Addition des Terms ψ(ti−1, x0) auf beiden Seiten bekommen wir

0 ≤ F (τ0, τ0 − ti−1, ψ(ti−1, ·))(x0) − ψ(ti−1, x0)

τ0 − ti−1+ψ(ti−1, x0) − ψ(τ0, x0)

τ0 − ti−1. (3.10)

Und mit der Definition fur ψ erhalten wir

ψ(ti−1, x0) − ψ(τ0, x0)

τ0 − ti−1= 4(ℜ + 1)

γε(τ0 − s0) − γε(ti−1 − s0)

τ0 − ti−1− 1

4TMP .

Dann bekommen wir mit (3.10), (F8) und der letzten Ungleichung

1

4TMP ≤ 4(ℜ + 1)

γε(τ0 − s0) − γε(ti−1 − s0)

τ0 − ti−1

−H(τ0, x0,Dψ(ti−1, x0)) + C5 · (1 + ‖Dψ(ti−1, ·)‖ + ‖D2ψ(ti−1, ·)‖)|P |, (3.11)

die die obere Grenze fur δ: δ < 113 ( siehe (3.21) und (F8)) garantiert.

Jetzt addieren wir (3.7) und (3.11) und bekommen1

2TMP ≤ 4(ℜ + 1)(

γε(τ0 − s0) − γε(ti−1 − s0)

τ0 − ti−1− γ′ε(τ0 − s0))

+H(s0, y0,Dϕ(s0, y0)) −H(τ0, x0,Dψ(ti−1, x0))

+C5 · (1 + ‖Dψ(ti−1, ·)‖ + ‖D2ψ(ti−1, ·)‖)|P |.

Da|Dϕ(s0, y0) −Dψ(ti−1, x0)| ≤ 4δ,

|H(t1, x1, p1) −H(t2, x2, p2)| ≤ |p2|N(|t1 − t2| + |x1 − x2|) +M |p1 − p2|(die Konstanten N und M sind in (f2) und (f3) definiert) gilt, ergibt sich mit (3.19), (3.20),(3.21) die Abschatzung

|D2γε(τ)| ≤15

ε2(3.12)

1

2TMP ≤ 60(ℜ + 1)

|P |ε2

+ (2L+ 6δ)N

(L+ 2δ

4(ℜ + 1)+L+ ℜ/T4(ℜ + 1)

)ε2 + 4δM

+C5 ·(1 + 8(ℜ + 1)(

1

ε+ δ) + 60(ℜ + 1)(

1

ε2+ δ)

)|P |.

Mit der Tatsache, dass ε = |P |1/4, |P | < 1 und δ geht gegen 0 erhalten wir

MP ≤ K1|P |1/2,

21


wobei

K1 = 2T

(60(ℜ + 1) +

2LN(2L + ℜ/T )

4(ℜ + 1)+ C5 · (1 + 68(ℜ + 1))

)

garantiert, dass

|P |1/2 < min1, 1

240(ℜ + 1)r.

2. Wir betrachten jetzt den Fall “–”, d.h.


(uP (τ, x) − u(τ, x))− .

Wir starten mit dem Fall τ0 ≥ 0, s0 = 0.

Mit Lemma 3.3.1 und Lemma 3.4.2 haben wir die folgende Kette von Ungleichungen

0 <1

4MP ≤ (uP (τ0, x0) − u(0, y0))

−

= u(0, y0) − uP (τ0, x0) = u0(0, y0) − u(0, x0) + u(0, x0) − uP (τ0, x0)

≤ L00|x0−y0|+maxL00, C1|τ0−s0| ≤(L00

L+ 4δ

4(ℜ + 1)+ maxL00, C1

L+ ℜ/T4(ℜ + 1)

)|P |1/2.

Jetzt betrachten wir den Hauptfall τ0 ≥ 0, s0 > 0.

(a) Wenn u(s0, y0) = u1(s0, y0) ist, dann ergibt sich

u1(s0, y0) = u(s0, y0) > uP (τ0, x0) ≥ u1(τ0, x0),

und deshalb gilt

0 > uP (τ0, x0) − u(s0, y0) ≥ u1(τ0, x0) − u1(s0, y0).

Aus (3.5), (3.6) des Lemmas 3.4.2 erhalten wir die Abschatzung

MP ≤ 4L00

(L+ 4δ

4(ℜ + 1)+L+ ℜ/T4(ℜ + 1)

)|P |1/2.

Daraus folgt die Behauptung des Satzes 3.4.1 fur den Fall (a).

(b) Wenn uP (τ0, x0) = u0(τ0, x0) ist, dann ergibt sich

u0(τ0, x0) = uP (τ0, x0) < u(s0, y0) ≤ u0(s0, y0),

und deshalb gilt

0 < u(s0, y0) − uP (τ0, x0) ≤ u0(s0, y0) − u0(τ0, x0).

Aus (3.5), (3.6) des Lemmas 3.4.2 erhalten wir die Abschatzung

MP ≤ 4L00

(L+ 4δ

4(ℜ + 1)+L+ ℜ/T4(ℜ + 1)

)|P |1/2.

22


Daraus folgt die Behauptung des Satzes 3.4.1 fur den Fall (b).

(c) Wir nehmen jetzt an: uP (τ0, x0) < u0(τ0, x0) und u(s0, y0) > u1(s0, y0) , die insbesondereτ0 > 0 bedeutet.

Schritt 1. Wir betrachten die Testfunktion ϕ,

ϕ(s, y) = uP (τ0, x0) − 4(ℜ + 1)βε(x0 − y) − 4(ℜ + 1)γε(τ0 − s) − 2δξδ(x0, y) +τ0 + s

4TMP .

Offensichtlich ist u(s, y) − ϕ(s, y) ≤ Ψ(x0, y, τ0, s) und deshalb nimmt u − ϕ das lokaleMaximum im Punkt (s0, y0) an. Aus (2.7) erhalten wir die Ungleichung

0 ≥ 4(ℜ + 1)γ′ε(τ0 − s0) +1

4TMP +H(s0, y0,Dϕ(s0, y0)),

und d.h.1

4TMP ≤ −4(ℜ + 1)γ′ε(τ0 − s0) −H(s0, y0,Dϕ(s0, y0)). (3.13)

Schritt 2. Wir betrachten jetzt eine andere Testfunktion ψ,

ψ(τ, x) = u(s0, y0) + 4(ℜ + 1)βε(x− y0) + 4(ℜ + 1)γε(τ − s0) + 2δξδ(x, y0) −τ + s0

4TMP .

Beachte, dass uP (τ, x)−ψ(τ, x) ≥ −Ψ(x, y0, τ, s0) und dann nimmt uP −ψ das lokale Mini-mum im Punkt (τ0, x0) an. Fur die Zeit ti−1, so dass τ0 ∈ (ti−1, ti] ist, gilt die Ungleichung

(uP − ψ)|(ti−1,x) ≥ (uP − ψ)|(τ0 ,x0)

fur jede x ∈ RN und

uP (ti−1, x) ≥ ψ(ti−1, x) + (uP (τ0, x0) − ψ(τ0, x0)). (3.14)

Die Behauptung des Lemmas 3.4.3 gilt fur die Funktion ψ. Mit der VoraussetzunguP (τ0, x0) < u0(τ0, x0) gilt

uP (τ0, x0) = F (τ0, τ0 − ti−1, uP (ti−1, ·))(x0).

Und mit dem Lemma 3.4.3 und Eigenschaften (F3), (F5) des Operators F erhalten wir von(3.14):

uP (τ0, x0) = F (τ0, τ0 − ti−1, uP (ti−1, ·))(x0)

≥ F (τ0, τ0 − ti−1, ψ(ti−1, ·))(x0) + (uP (τ0, x0) − ψ(τ0, x0)). (3.15)

Wir addiern den Term ψ(ti−1, x0) zu beiden Seiten in (3.15), dann bekommen wir

0 ≥ F (τ0, τ0 − ti−1, ψ(ti−1, ·))(x0) − ψ(ti−1, x0)

τ0 − ti−1+ψ(ti−1, x0) − ψ(τ0, x0)

τ0 − ti−1. (3.16)

Aus der Definition fur ψ gilt

ψ(ti−1, x0) − ψ(τ0, x0)

τ0 − ti−1= 4(ℜ + 1)

γε(ti−1 − s0) − γε(τ0 − s0)

τ0 − ti−1+

1

4TMP .

23


Aus (3.16),(F8) und der letzten Ungleichung bekommen wir

1

4TMP ≤ 4(ℜ + 1)

γε(τ0 − s0) − γε(ti−1 − s0)

τ0 − ti−1

+H(τ0, x0,Dψ(ti−1, x0)) + C5 · (1 + ‖Dψ(ti−1, ·)‖ + ‖D2ψ(ti−1, ·)‖)|P |, (3.17)

was die obere Grenze fur δ: δ < 16 ( siehe (3.21) und (F8) ) garantiert. Addition von (3.13)

und (3.17) gibt

1

2TMP ≤ 4(ℜ + 1)

(γε(τ0 − s0) − γε(ti−1 − s0)

τ0 − ti−1− γ′ε(τ0 − s0)

)

+H(τ0, x0,Dψ(ti−1, x0)) −H(s0, y0,Dϕ(s0, y0))

+C5 · (1 + ‖Dψ(ti−1, ·)‖ + ‖D2ψ(ti−1, ·)‖)|P |.Da

|Dψ(ti−1, x0) −Dϕ(s0, y0)| ≤ 4δ,

|H(t1, x1, p1) −H(t2, x2, p2)| ≤ |p1|N(|t1 − t2)| + |x1 − x2|) +M |p1 − p2|(die Konstanten N , M sind in (f2) and (f3) definiert) gilt, ergibt sich mit (3.19), (3.20),(3.21) und (3.12) die Abschatzung

1

2TMP ≤ 60(ℜ + 1)

|P |ε2

+ (2L+ 6δ)N

(L+ 2δ

4(ℜ + 1)+L+ ℜ/T4(ℜ + 1)

)ε2 + 4δM

+C5 ·(

1 + 8(ℜ + 1)(1

ε+ δ) + 60(ℜ + 1)(

1

ε2+ δ)

)|P |.

Mit der Tatsache, dass ε = |P |1/4, |P | < 1 und δ geht gegen 0 erhalten wir

MP ≤ K1|P |1/2,

wobei

K1 = 2T

(60(ℜ + 1) +

2LN(2L+ ℜ/T )

4(ℜ + 1)+C5 · (1 + 68(ℜ + 1))

),

garantiert, dass

|P |1/2 < min

1,

1

240(ℜ + 1)r

.

Der Satz 3.4.1 ist damit bewiesen.

Bemerkung. Die Konvergenz der Approximationslosung zur exakten Losung ist mit demSatz 3.4.1 bewiesen und die Ungleichung (3.3) zeigt, dass

u0(τ, x) ≥ u(τ, x) ≥ u1(τ, x), ∀(τ, x) ∈ [0, T ] × RN (3.18)

gilt.

24


Im Beweis des Satzes haben wir ein paar Aussagen verwendet, die wir jetzt zeigen. Wirbeginnen mit dem Beweis vom Lemma 3.4.2

Beweis vom Lemma 3.3. Wir schatzen zuerst den Abstand |x0−y0| ab. Man betrachtetdazu

Ψ(x0, y0, τ0, s0) − Ψ(x0, x0, τ0, s0)

= (uP (τ0, x0) − u(s0, y0))± − (uP (τ0, x0) − u(s0, x0))

±+

+2δξδ(x0, y0) − 2δξδ(x0, x0) + 4(ℜ + 1)βε(x0 − y0) − 4(ℜ + 1)

≤ |u(s0, y0) − u(s0, x0)| + 2δ(ξδ(x0, y0) − ξδ(x0, x0)) + 4(ℜ + 1)βε(x0 − y0) − 4(ℜ + 1).

Wenn |x0 − y0| >√

3/4ε ist, dann ist βε(x0 − y0) < 1/4 und mit δ < 3/4 bekommen wir

Ψ(x0, y0, τ0, s0) < Ψ(x0, x0, τ0, s0)

mit der bestimmten Wahl der ℜ. Diese Ungleichung widerspricht (3.4).D.h wir haben |x0 − y0| ≤

√3/4ε und

Ψ(x0, y0, τ0, s0) − Ψ(x0, x0, τ0, s0)

≤ (L+ 2δ)|x0 − y0| − 4(ℜ + 1)|x0 − y0|2

ε2= |x0 − y0|

4(ℜ + 1)

ε2

(L+ 2δ

4(ℜ + 1)ε2 − |x0 − y0|

).

Die letzte Ungleichung zeigt uns, dass

|x0 − y0| ≤L+ 2δ

4(ℜ + 1)ε2.

Andererseits gilt Ψ(x0, y0, τ0, s0) < Ψ(x0, x0, τ0, s0) und das widerspricht (3.4).

Als nachstes schatzen wir |τ0 − s0| ab. Man betrachtet

Ψ(x0, y0, τ0, s0) − Ψ(x0, y0, τ0, τ0)

= (uP (τ0, x0) − u(s0, y0))± − (uP (τ0, x0) − u(τ0, y0))

±−

−τ0 + s04T

MP +τ0 + τ0

4TMP + 4(ℜ + 1)γε(τ0 − s0) − 4(ℜ + 1)

≤ |u(τ0, y0) − u(s0, y0)| +|τ0 − s0|

4TMP + 4(ℜ + 1)γε(τ0 − s0) − 4(ℜ + 1).

Wenn |τ0 − s0| >√

3/4ε ist, dann ist γε(τ0 − s0) < 1/4, und mit der bestimmten Wahl derFunktion ℜ erhalten wir

Ψ(x0, y0, τ0, s0) < Ψ(x0, y0, τ0, τ0).

Das ist wieder ein Widerspruch. Wir haben dann |τ0 − s0| ≤√

3/4ε und

Ψ(x0, y0, τ0, s0) − Ψ(x0, y0, τ0, τ0)

≤ (L+ ℜ/T )|τ0 − s0| − 4(ℜ + 1)|τ0 − s0|2

ε2= |τ0 − s0|

4(ℜ + 1)

ε2

(L+ ℜ/T4(ℜ + 1)

ε2 − |τ0 − s0|).

25


Wir konnen auch aus der letzen Ungleichung sehen, dass

|τ0 − s0| ≤L+ ℜ/T4(ℜ + 1)

ε2

gilt. Auf der anderen Seite gilt Ψ(x0, y0, τ0, s0) < Ψ(x0, y0, τ0, τ0) und das ist ein Wider-spruch.

Aus der Definition fur Ψ und (3.4) haben wir die Ungleichung

(uP (τ0, x0) − u(s0, y0))± + 8(ℜ + 1) + 2δ

≥ Ψ(x0, y0, τ0, s0) ≥ Ψ(x, x, τ, τ) ≥ (uP (τ, x) − u(τ, x))± + 8(ℜ + 1) − 1

2MP

fur jede (τ, x) ∈ [0, T ] × RN .

Allerdings

(uP (τ0, x0) − u(s0, y0))± + 2δ ≥ 1

2MP ,

und fur δ < 18MP haben wir

(uP (τ0, x0) − u(s0, y0))± ≥ 1

4MP > 0.

Das beweist das Lemma (3.4.2) und ermoglicht uns den Beweis des Satzes 3.4.1.

Beweis vom Lemma 3.4. Aus der Definition fur ψ haben wir

Dψ(τ, x) = −4(ℜ + 1)Dβε(x− y0) − 2δDξδ(x, y0).

Die rechte Seite ist nicht abhangig von τ . Aus der Wahl der Funktionen βε, ξδ erhalten wir

|Dψ(τ, x)| ≤ 8(ℜ + 1)(1

ε+ δ), (3.19)

|D2ψ(τ, x)| ≤ 60(ℜ + 1)(1

ε2+ δ). (3.20)

Und mit dem Lemma 3.4.2 bekommen wir

|Dψ(τ, x0)| = | − 4(ℜ + 1)Dβε(x0 − y0) − 2δDxξδ(x0, y0)|

≤ 8(ℜ + 1)|x0 − y0|ε2

+ 2δ ≤ 2L+ 6δ. (3.21)

Als nachstes beachten wir, dass

|ψ(ti−1, x0 + w) − ψ(ti−1, x0 + w)|≤ |ψ(ti−1, x0 + w) − ψ(ti−1, x0 + w) −Dψ(ti−1, x0) · (w − w)|

+|(Dψ(ti−1, x0) −Dψ(τ0, x0)) · (w − w)| + |Dψ(τ0, x0) · (w − w)|.

Mit (3.20), (3.21), den Voraussetzungen fur |P | und fur w, w: w, w ∈ Bn(0, (τ0 − ti−1)r)erhalten wir

|ψ(ti−1, x0 + w) − ψ(ti−1, x0 + w)|

26


≤ (120(ℜ + 1)(τ0 − ti−1)r(1

ε2+ δ) + 2L+ 6δ)|w − w|

≤ (2L+1

2+δ

2+ 6δ)|w − w|,

wobei ε = |P |1/4 ist. Die Wahl der Konstante δ liefert die erste Behauptung des Lemmas.Die zweite Behauptung ist die Folge des Lemmas 3.3.1(a). Dies beweist das Lemma 3.4.3und erlaubt uns den Beweis des Satzes 3.4.1 fortzusetzen.

Das dargestellte stabile Verfahren wurde in ein paar Modellbeispielen getestet. Bevor dieBeispiele prasentiert werden betrachten wir im nachsten Unterkapitel eine Methode, dieman benutzt um die optimalen Strategien fur die beiden Spieler zu bestimmen.

3.5 Optimale Strategien fur die Spieler

Die Wahl der optimalen Steuerungsstrategien fur die Spieler basiert auf einem Verfahrender extremalen Ziele, das in [20] beschrieben ist.Seien ε eine kleine positive Zahl, tn die gegenwartige Zeit. Man betrachtet eine ε-Umgebung

Uε = x ∈ Rn : |x− x(tn)| ≤ ε

fur den gegenwartigen Zustand x(tn) des Kontrollsystems (2.3).

Alle Gitterpunkten x# ∈ Uε werden durchsucht um einen Punkt x#∗ mit

cP (tn, x#∗ ) = min

x#∈Uε

cp(tn, x#)

zu finden, wobei cp ist die Gitterapproximation der Wertefunktion ist. Die gegenwartigeSteuerung α(tn)( β(tn)), die in nachsten Zeitintervall [tn, tn+1) verwendet wird, ist ausder Bedingung fur die maximale(minimalen) Projektion der Geschwindigskeit f auf die

Richtung des Vektors s = x#∗ − x(tn) zu bestimmen. Damit ergibt sich

α(tn) : maxβ∈Q

〈s, f(tn, x(tn), α(tn), β)〉 = minα∈P

maxβ∈Q

〈s, f(tn, x(tn), α, β)〉 , (3.22)

β(tn) : minα∈P

〈s, f(tn, x(tn), α, β(tn))〉 = maxβ∈Q

minα∈P

〈s, f(tn, x(tn), α, β)〉 . (3.23)

3.6 Modellbeispiele

Die Beispiele, die wir hier betrachten, demonstrieren die Probleme, die man mit demFunktional (2.2) beschreiben kann. Die Beispiele zeigen auch die effektive Anwendung desUpwind-Operators (3.1).

27


3.6.1 Nichtlineares Pendel (Nonlinear pendulum)

Das erste Beispiel ist das nichtlineare Pendel mit dem dynamischen Kontrollsystem

x1 = x2 + β, x2 = −gl sinx1 −α

lcos x1, |α| ≤ 1, |β| ≤ 1. (3.24)

Hier ist x1 die Winkelauslenkung, x2 die Winkelgeschwindigkeit, α die waagerechte Kraftdes ersten Spielers, β die Storunginformation (die wird vom zweiten Spieler beherrscht), ldie Lange und g die Gravitationskraft. Die Masse ist gleich 1 angenommem. Sei

χ(x1, x2) := 10√x2

1 + x22 und θ(x1, x2) := 2.5|x1|.

Wir betrachten das Differentialspiel mit der Zielmenge M = (x1, x2) : χ(x1, x2) ≤ 1 undder Zustandsbeschrankungsmenge N = (x1, x2) : θ(x1, x2) ≤ 1. Die Anfangszeit ist t0 = 0und die Endzeit ist T = 5. Wir betrachten das Differentialspiel mit dem dynamischen System(3.24) und dem Kostenfunktional der Form (2.2). Wir nehmen an, dass der Anfangszustandx0 = (x10, x20) in der Losungsmenge W = (x1, x2) : c(t0, x1, x2) ≤ 1 liegt, wobei c dieWertefunktion ist. Eine optimale Trajektorie ist in der Abbildung 3.1 dargestellt, die dieZielmenge erreicht und bleibt immer in der Zustandsbeschrankungsmenge.

-1 -0.5 0 0.5 1

-1.5

-1

-0.5

0

0.5

1

1.5

x0)

x2

x1

W

N

MAA

Abbildung 3.1: Die Losungsmenge W (rot), die Zielmenge M (blau), die Zustandsbe-schrankungsmenge N (purpurrot), und die optimale Trajektorie mit derStartposition x0 (grun).

28


3.6.2 Materialpunkt (Material point)

Das nachste bekannte Beispiel ist das Differentialspiel der zwei Materialpunkte auf derGerade. In relativen Koordinaten wird das folgende dynamische System betrachtet

x1 = x2 + β, x2 = α, |α| ≤ 1, |β| ≤ 0.2. (3.25)

Hier ist x1 die relative Position der zwei Punkten, x2 die relative Geschwindigkeit. Der ersteSpieler verwendet die Steuerungsvariable α um die Relation

√x1(τ)2 + x2(τ)2 ≤ 0.2 in der

Zeit τ zu garantieren. Die Zustandsbeschrankungsmenge ist durch |x2(t)| ≤ 0.5 fur allet ≤ τ gegeben. Das Ziel des zweiten Spielers ist das Gegenteil mit der Steuerungsvariableβ zu erreichen. Wir betrachtnen das Differentialspiel mit dem dynamischen System (3.25),dem Kostenfunktional der Form (2.2) und mit

χ(x1, x2) := 5√x2

1 + x22 und θ(x1, x2) := 2|x2|.

Die Anfangszeit ist t0 = 0 und die Endzeit ist T = 5.Die Abbildung 3.2 zeigt optimale Trajektorien fur die beiden Spieler. Die Trajektorie auf derlinke Seite beginnt in einem Punkt der Losungsmenge W = (x1, x2) : c(t0, x1, x2) ≤ 1 underreicht die Zielmenge. Wenn die Trajektorie nahe der Grenze, die Zustandsbeschrankungdefiniert, ist, dann wird das Zick-Zack-Verhalten beobachtet (siehe die Abbildung 3.3). DieTrajektorie auf der linken Seite beginnt in einem Punkt außerhalb der Losungsmenge undin diesem Fall gewinnt das Differentialspiel der zweite Spieler.

-1.5 -1 -0.5 0 0.5 1 1.5

-1

-0.5

0

0.5

1

x2

x1

W N

M

Abbildung 3.2: Die Losungsmenge W (rot), die Zielmege M (blau), die Zustandsbe-schrankungsmenge N (purpurrot), die optimale Trajektorie fur den erstenSpieler (grun) und die Gewinntrajektorie fur den zweiten Spieler (schwarz).

29


-1.5 -1.4 -1.3 -1.2

0.4

0.5

0.6

x2

x1

Abbildung 3.3: Das Zick-Zack-Verhalten.

3.6.3 Das Spiel der zwei Autos (Game of two cars)

Das Differentialspiel der zwei Autos wurde in [19] prasentiert. In diesem Spiel hat das ersteAuto einen Verfolger-Charakter und das zweite Auto hat einen Fluchtling-Charakter. Genauwie in der Arbeit [23] werden wir einen Fall betrachten, wobei die zwei Autos die gleichenlinearen Geschwindigkeiten haben.In einem beweglichen Referenz-Koordinatensystem (sieheAbbildung 3.4) konnte die Dynamik der relativen Bewegungen durch

x = −yα+ sinφ, y = xα+ cosφ− 1, φ = −α+ β, |α| ≤ 1, |β| ≤ 1 (3.26)

beschrieben werden. Die terminalen Bedingungen sind√x2 + y2 ≤ r, cos (φ− φf (x, y)) − cosφf (x, y) ≤ 0 at

√x2 + y2 = r,

wobei φf (x, y) =

arccos ( x/

√x2 + y2), y ≥ 0,

π + arccos (−x/√x2 + y2), y < 0,

D.h. die Distanz zwischen zwei Autos ist kleiner oder gleich dem gegebenen Radius r, unddie relative radiale Geschwindigkeit ist nicht positiv. Die Endzeit ist tf = 10. Wir habendieses Spiel mit der Zustandsbeschrankung |φ + 2| ≤ s fur die dritte Variable erweitert.Dabei haben wir das Kostenfunktional der Form (2.2) mit

χ(x, y, φ) = max√

x2 + y2 − r, cos (φ− φf (x, y)) − cosφf (x, y),

θ(x, y, φ) = |φ+ 2| − s

betrachtet.Sei cP die Approximation der Wertefunktion. Die Abildung 3.5 zeigt die Losungsmenge, diedurch (x, y, φ) : cP (x, y, φ) ≤ 0 definiert ist und den Querschnitt bei φ = −2 hat. DieBerechnungen wurden auf einem Linux SMP-Computer mit dem Gitter 300×300×300 undder Zeitschritweite 10−3 durchgefuhrt und die Laufzeit betrag dabei circa 10 Minuten.

30


VF

VV

x

y

(x, y)

φ

F

V

Abbildung 3.4: Bewegliches Koordinatensystem. Hier ist y die Achse, die parallel zumVerfolger-Geschwindigkeitsvektor ist. φ ist die Winkel zwischen der y-Achseund dem Fluchtling-Geschwindigkeitsvektor, x ist so gewahlt, dass x, y, φdas rechte Koordinatensystem ist.

x

y

φ

φ = −2

6

-

y

x

Abbildung 3.5: Das rechte Bild: Die Losungsmenge fur den Fall ohne Zustandsbe-schrankungen (blau) und der Querschnitt bei φ = −2 (rot). Das linke Bild:Der Querschnitt der drei Losungsmengen: a) ohne Zustandsbeschrankungen(der großte); b) mit der Zustandsbeschrankungsmenge, die mit s = 1 de-finiert ist (der mittlere); c) mit der Zustandsbeschrankungsmenge, die mits = 0.8 definiert ist (der kleinste).

Im nachsten Kapitel betrachten wir eine Anwendung der Hamilton-Jacobi Gleichungen beider Kryopreservation von lebenden Zellen.

31

4 Anwendung: Optimierung vonKuhlprotokollen bei der Kryopreservationvon lebenden Zellen

In diesem Kapitel werden wir auf die Probleme eingehen, die das Bestandteil des Scwer-punktprogramms der Deutsche Forschungsgemeinschaft SPP1253: Optimierung mit partiel-len Differentialgleichungen. Unter dem Projekt Optimale Steuerung der Kryopreservationvon lebenden Zellen und Gewebe (Optimal control in cryopreservation of cells andtissues, siehe [16]), das beim Professor Dr. Dr.h.c.mult Karl-Heinz Hoffmann geleitet wird,wird die Optimierung von Kuhlprotokollen bei der Kryopreservation von lebenden Zellenbetrachtet. Innerhalb der Forschungsgruppe mit Dr. Nikolay Botkin und Dr. Varvara Tu-rova wurde die Entwicklung von Kuhlprotokollen vorgenommen.

In der Entwicklung von Kuhlprotokollen liegt der Schwerpunkt bei der Minimierung derSchadigungen der lebenden Zellen wahrend der Kuhlung. Die Schadigungen sind die Folgender zwei folgenden physikalischen Effekte:

1. Das Wasser gefriert nicht gleichzeitig innerhalb und außerhalb der Zelle, was Span-nungen verursacht.

2. Das Eis bildet sich wahrend der Kuhlung außerhalb der Zelle schneller als innerhalb.

Die großten Schaden werden durch die entstehenden Spannungen auf der Zellmembranwahrend der langsamen Abkuhlung verursacht. Die schnelle Abkuhlung hilft bei dem Pro-blem leider nicht. In dem Fall bilden sich in der Zelle kleine Eiskristalle. Die Eiskristallesind in der Regel instabil mit einer unregularen Form.

Die entsprechenden Modelle sind in der Regel durch die partiellen Differentialgleichungendarstellbar. Die Gleichungen beschreiben die Dynamik der Phasen im jeden raumlichenPunkt [7, 14, 18]. Die raumliche Verteilung ist nicht sehr wichtig, wenn die kleinen Objektewie die lebenden Zellen untersucht werden. Die praktischen Beobachtungen zeigen, dassdie Mittelungsmethoden das Modell zum System der gewohnlichen Differentialgleichungenvereinfachen. Die gewohnlichen Differentialgleichungen enthalten in der Regel Steuerungs-parametern und nichtlineare Abhangigkeiten, die als Daten in der Tabellen gegeben sind.Dies verkompliziert die Anwendung der traditionellen Methoden, die auf dem Maximum-prinzip von Pontrjagin basieren.

Die Anwendung der Kontrolltheorie wird ein gleichzeitiges Gefrieren des Fluids außerhalbund innerhalb der Zelle fur die langsame Abkuhlung ermoglichen. Es wird die Theorieder Hamilton-Jacobi-Bellman-Isaacs Gleichungen angewandt. Der numerische Ansatz zurBestimmung eines optimalen Kuhlprotokoll, wurde in [17, 18] entwickelt, der auf dieserTheorie grundet.

32

4 Optimierung von Kuhlprotokollen bei der Kryopreservation von lebenden Zellen

4.1 Mathematisches Model

Wir betrachten eine Zelle, die in einer Pore liegt. Die Pore ist mit einer extrazellularenFlussigkeit gefullt. Das Fluid außerhalb der Zelle gefriert schneller als innerhalb und dasVolumen des entstehenden Eises lost den Druck auf die Zelle auf. Der Effekt konnte appro-ximativ erfasst werden

p ≈ Ceis · a · (1 − βl), (4.1)

wobei p der Druck ist und Ceis ist die Elastizitat des Eises, a die Ausdehnungkoeffizient furdas Eis, βl der Anteil der nicht gefrorenen Wasser in der Zelle. Die grobe Abschatzung furden Druck liefert: p ≈ 1Bar, was die totalen Schaden an der Zellmembran verursacht. Umdies zu vermeiden, wird eine bestimmte Methode angewandt. Die Methode ermoglicht denniedriegen Gefrierpunkt des Wassers. Dabei konnte das Wasser außerhalb und innerhalbder Zelle gleichzeitig gefrieren.

Wir betrachten das folgende Schema fur eine Zelle:

Γ1 Ω1 θ1 θ1s

Γ2 Ω2 θ2 θ2s E2

λ1 θext

λ2

@@Iextrazellulares Fluid

@@

@@@R

Zelle

Abbildung 4.1: Das Schema der Zelle in der Pore.

Das Gebiet Ω1 ⊂ R3 mit einem C0,1-Rand Γ1, d.h. Γ1 ist Lipschitz-stetig, bezeichnet das Ge-

biet der Pore. Das Gebiet Ω2 ⊂ R3 mit einem C0,1-Rand Γ2 bezeichnet das Gebiet der Zelle.

θ1 und θ2 sind die Temperaturen innerhalb und außerhlab der Zelle. θext ist die Temperaturaußerhalb der Pore. Als nachstes fuhren wir die Koeffizienten λ1 und λ2 ein. Sie bezeichnendie dunnschichtbildenden Warmeubergangskoeffizienten des Randes der Pore und der Zell-membran. θ1s und θ2s sind die Gefrierpunkte des Wassers innerhalb und außerhalb der Zelle.

Außerhalb der Zelle werden wir das Phasen-Feld Modell nach G.Caginalp [7, 8] betrachten:

ρθt + ρLφt −K∆θ = 0, (4.2)

τφt − ε∆φ+W ′(φ) − 2θ = 0, (4.3)

−K∂θ

∂ν|∂Ω = λ(θ − θe)|∂Ω, (4.4)

33


∂φ

∂ν= 0, (4.5)

θ|t=0 = θ0 = const > 0, φ|t=0 = φ0 = 1, (4.6)

wobei Ω ein Gebiet aus R3 (Ω ⊂ R

3) mit einem C0,1-Rand ist und θ ist die Celsius Tempera-tur, θe die Temperatur außerhalb des Gebietes (z.B. außerhalb der Pore), K Warmeleitfahig-keitkoeffizient, λ Warmeubergangskoeffizient, C die spezifische Warmekapazitat, L die la-tente Warme, ε eine beliebig kleine Zahl, ε > 0, τ ein Relaxationsparameter, ν die außereNormale, ρ die Dichte. Wir werden annehmen, dass die Dichten des Fluids und des Eisesgleich sind. Die Funktion φ, die unabhangig von der Temperatur ist, ist die Phasen-Funktion.φ beschreibt den Anteil des nicht gefrorenen Fluids im Gebiet Ω,W Doppelmuldenpotential-Funktion (double well function).Die Funktion φ = −1 fur den Eis-Zustand und φ = 1 fur den flussigen Zustand.

Im unseren Fall ist Ω = Ω1, Γ1 ∈ C1, θ = θ1, θe = θext, λ = λ1, W (φ) = W (2φ − 1) mitW (φ) = 1

8(1− φ2)2 (siehe die Abb. 4.2) und die Funktion φ beschreibt den Anteil des nichtgefrorenen Fluids außerhalb der Zelle, φ = 0 fur die Eis-Phase und φ = 1 fur die flussigePhase.

-

6

-1 0 1 φ

W (φ)

W (φ)@@R

Abbildung 4.2: Doppelmuldenpotential-Funktionen als Funktionen der Phasenfeldvariab-len.

Der Anteil des nicht gefrorenen Fluids in der Zelle kann mit Hilfe eines Phasen-Feld Modellsnach M. Fremond berechnet werden[14]. Der Hauptunterschied dieses Modells liegt beider Phasen-Funktion im Vergleich zum Modell vom G.Caginalp. Die Phasen-Funktion istabhangig von der Temperatur.Wir betrachten die Phasen-Funktion βl. Sie bezeichnet den Anteil des nicht gefrorenenFluids innerhalb eines Gebietes Ω. Die Funktion βl wird nicht durch eine Gleichung definiert,sondern ist sie eine gegebene Funktion von der Temperatur. Dann gilt:

∂E(θ)

∂t−K∆θ = 0, (4.7)

34


−K∂θ

∂ν|∂Ω = λ(θ − θe)|∂Ω, (4.8)

E(θ) = ρCθ + ρLβl(θ), (4.9)

βl(θ) = ψ(L(θ − θs)

(T0 + θs)(T0 + θ)), (4.10)

wobei θ, θe die Celsius Temperaturen sind, θs der Gefrierpunkt, T0 der Celsius Nullpunkt(273K), K Warmeleitfahigkeitkoeffizient, λ Warmeubergangskoeffizient, C die spezifischeWarmekapazitat, L die latente Warme, ρ die Dichte. Die Funktion E(θ) ist die innere Ener-gie. Die Funktion ψ ergibt sich experimentell.

Im unseren Fall ist Ω = Ω2, Γ2 ∈ C1, θ = θ2, θe = θ1, θs = θ2s, λ = λ2, E(θ) = E2(θ2) unddie Funktion βl beschreibt den Anteil des nicht gefrorenen Fluids innerhalb der Zelle.

Um zu garantieren, dass der Gefrierprozess außerhalb der Zelle langsamer verlauft, wird diefolgende Bedingung gestellt:

θ1s − θ2s = −13.

Das bedeutet, dass die Gefrierpunkte in der Pore und in der Zelle verschieden sind. Damitwird das gleichzeitige Gefrieren ermoglicht.

Dann ist unser System folgendermaßen darstellbar:

ρθ1 + ρLφ−K∆θ1 = 0,

τ φ− ε∆φ+W ′(φ) − 2(θ1 + 13) = 0, (4.11)

ρθ2 + ρLβlt(θ2) −K∆θ2 = 0.

Weiter integrieren wir uber die Gebiete Ω1 und Ω2. Mit der folgenden Notation

E2 =1

|Ω2|

∫

Ω2

E2dV , (4.12)

θi =1

|Γi|

∫

Γi

θidS, (4.13)

θext =1

|Γ1|

∫

Γ1

θEdS, (4.14)

αi =Γ2

|Ωi|λ2, (4.15)

λ =Γ1

|Ω1|λ1, (4.16)

und der Vorraussetzung |Γ1|−1 ∫Γ1θ1dS ≈ |Γ2|−1 ∫

Γ1θ1dS (da θ1 fast uberall konstant auf

einem kleinen Gebiet ist) erhalten wir das System der gewohnlichen Differentialgleichungen:

ρ˙θ1 + ρL

˙φ = −α1(θ1 − θ2) − λ(θ1 − θext),

35


τ˙φ+W ′(φ) − 2(θ1 + 13) = 0, (4.17)

ρ˙θ2 + ρLβlt(

˙θ2) = −α2(θ2 − θ1).

Mit dem Zusammenhang zwischen E2 und θ2:

E2 = ρθ2 + ρLβl(θ2), (4.18)

wobei βl durch (4.10) mit θ2s definiert ist, ergibt sich

ρ˙θ1 + ρL

˙φ = −α1(θ1 − θ2) − λ(θ1 − θext),

τ˙φ+W ′(φ) − 2(θ1 + 13) = 0, (4.19)

˙E2 = −α2(θ2 − θ1).

Wir merken an, dass die Ableitung ∂E2/∂θ2 nicht stetig ist. Deshalb werden wir die Tem-peratur θ2 durch die Energie E2 ausdrucken. Dabei wird (4.18) invertiert: θ2 = Θ2(E2). AlsFolge erhalten wir das System der gewohnlichen Differentialgleichungen:

ρ˙θ1 + ρL

˙φ = −α1(θ1 − Θ2(E2)) − λ(θ1 − θext),

τ˙φ+W ′(φ) − 2(θ1 + 13) = 0, (4.20)

˙E2 = −α2(Θ2(E2) − θ1).

Eine weitere Aufgabe dieser Diplomarbeit war es, sowohl ein Funktional als auch geeigneteZustandbeschrankungen einzufuhren, um das gewunschte System zu losen.

Wir betrachten das dynamische Kontrollsystem mit x1 := θ1, x2 := φ, x3 := E2, x4 := θext,αi := αi, λ = λ:

ρx1 + ρLx2 = −α1(x1 − Θ2(x3)) − λ(x1 − x4) + β1,

τ x2 +W ′(x2) − 2(x1 + 13) = 0, (4.21)

x3 = −α2(Θ2(x3) − x1) + β2,

x4 = α.

Hier ist x4 die Temperatur außerhalb des Gebietes, α die Kontrollvariable fur die Abkuhlungmit |α| ≤ µ, β1 und β2 die Storungen, |β1| ≤ ν, |β2| ≤ ν.Laut der Definition der Funktionen βl und φ kann das gleichzeitige Gefrieren außerhalb undinnerhalb der Zelle mit dem Verschwinden des folgenden Funktionals interpretiert werden:

36


J = maxt∈[0,tf ]

γ(x2(t), x3(t)), (4.22)

mit γ(x2, x3) = |x2 − βl(Θ2(x3))| .Die Funktion γ charakterisiert die Abweichung des Anteils vom Eis in der Region außerhalbder Zelle vom Anteil innerhalb der Zelle.Wir betrachten als nachstes das Differentialspiel (4.21) und (4.22). Das Ziel des erstenSpielers ist das Funktional (4.22) mit der Kontrollvariable α zu minimieren. Das Ziel deszweiten Spielers ist mit den Storungen das Gegenteil zu erreichen. Die folgenden Zustand-beschrankungen werden sinnvoll:

1. Eine Beschrankung auf die Phasen-Funktion

ζ1 < x2 < ζ2, ζ1, ζ2 < 0, x2 = φ. (4.23)

In dem Fall wird der Gefrierprozess außerhalb der Zelle langsamer verlaufen.

2. Eine Beschrankung auf die Geschwindigkeit

x3 < ζ3, ζ3 < 0, x3 =˙E2. (4.24)

Dieser Fall stimuliert einen schnellen Verlauf des Gefrierprozesses innerhalb der Zelle.

Das Differentialspiel wird numerisch mit einem Approximationschema fur die viskose Losungender Hamilton-Jacobi Gleichungen gelost.

4.2 Simulation

Wir betrachten das folgende Kontrollsystem:

x1 = l1

τ(W ′(x2) − 2(x1 + 13)) − α1(x1 − Θ2(x3)) − λ(x1 − α) + β1,

x2 = −1

τ(W ′(x2) − 2(x1 + 13)) = 0, (4.25)

x3 = −α2(Θ2(x3) − x1) + β2,

wobei l = L/ρ ist, |α+ 20| ≤ µ, |β1| ≤ ν, |β2| ≤ ν.

Nach der Vereinfachung sieht unser System folgendermaßen aus:

x1 = l1

τ(W

′(2x2 − 1) − x1 − 13) − α1(x1 − Θ2(x3)) − λ(x1 − α) + β1,

x2 = −1

τ(W

′(2x2 − 1) − x1 − 13), (4.26)

x3 = −α2(Θ2(x3) − x1) + β2,

Weiter betrachten wir das Differentialspiel dem dynamischen System (4.26) und dem Kos-tenfunktional der Form (4.22). Die Anfangszeit ist t = 0 und die Endzeit T = 900. Die

37


Funktion W wird nochmal modifiziert, um ihre Steigung auf den Intervallen [−∞, 0]×[1,∞]zu versteilen. Die folgenden Koeffizienten und die Beschrankungen fur die Kontrollvariablenwurden gewahlt: α1 = α2 = 0.01, λ = 0.2, µ = 4, ν = 0.2.

Die Notation βi = 1 − βl(Θ(x3)) und φi = 1 − φ = 1 − x2 wurde verwendet um denEis-Anteil darzustellen. Die Abbilding 4.3 zeigt wie die kontrollierten Temperaturen sichinnerhalb und außerhalb der Zelle ohne Zustandbeschrankungen verlaufen und wie denEis-Anteil sich bildet. Auf der Abbildung 4.4 kann man den Verlauf der Steuerung sehen.

Abbildung 4.3: Links: Die Temperaturen innerhalb und außerhalb der Zelle, θ1 = x1 ist rot,θ2 = Θ2(E2) = Θ2(x3) ist grun. Rechts: Eis-Anteile innerhalb und außerhalbder Zelle, φi ist rot, βi ist grun.

Abbildung 4.4: Der Verlauf der Steuerung α.

38


Danach legen wir die Beschrankung auf die Geschwindigkeit x3 an.Sei

x3 + 0.15 < 0.

Die Anfangszeit ist t = 0 und die Endzeit T = 580. Die Abbilding 4.5 zeigt wie die kon-trollierten Temperaturen sich innerhalb und außerhalb der Zelle verlaufen und wie denEis-Anteil sich bildet. Auf der Abbildung 4.6 kann man den Verlauf der Steuerung und derEnergie x3 sehen.


Abbildung 4.6: Links: Der Verlauf der Steuerung α. Rechts: Der Verlauf der Energie x3.

39


Weiter betrachten wir die Beschrankung auf die Phasen-Funktion x2.Sei

1 − 1

700t < x2 < 1 − 1

900t.

Die Anfangszeit ist t = 0 und die Endzeit T = 1000. Die Abbilding 4.7 zeigt wie diekontrollierten Temperaturen sich innerhalb und außerhalb der Zelle verlaufen und wie denEis-Anteil sich bildet. Auf der Abbildung 4.8 kann man den Verlauf der Steuerung und derPhasen-Funktion, die zwischen zwei Beschrankungsgeraden liegt, sehen.


Abbildung 4.8: Links: Der Verlauf der Steuerung α. Rechts: Der Verlauf der Phasen-Funktion x2 (rot).

40

5 Anwendung: Mehrphasenstromungen

Im diesen Kapitel werden wir eine Anwendung der Hamilton-Jacobi-Gleichung in einenModell fur Mehrphasenstromungen in einem porosen Medium (zum Beispiel Sand) be-sprechen. Die Modelle fur Mehrphasenstromungen sind aktuell in der Industrie, die sichmit der Frage CO2-Sequestration (CO2-Abscheidung und -Speicherung) beschaftigt. CO2-Abscheidung und -Speicherung ist die Abscheidung von Kohlendioxid (CO2) insbesondereaus Verbrennungs-Abgasen sowie dessen Injektion und behalterlose Lagerung in tiefen un-terirdischen Gesteinsschichten auf unbegrenzte Zeit (Sequestrierung, Lagerung). Durch dieSpeicherung soll weniger CO2 in die Atmosphare gelangen, denn dort wirkt es als Treib-hausgas.Das Hauptanwendungsgebiet soll die klimaschonende Nutzung fossiler Rohstoffe bei derStromerzeugung in Kraftwerken werden. Momentan befindet sich die CO2-Abscheidungund -Speicherung noch im Entwicklungsstadium.Bis auf Weiteres existieren nur kleine Pilotanlagen mit geringer Kapazitat. Angestrebtwird eine Langzeitsicherheit, also ein Zustand, der gewahrleistet, dass das gespeicherteKohlendioxid unter Berucksichtigung der erforderlichen Vorsorge gegen Beeintrachtigungenvon Mensch und Umwelt vollstandig und auf unbegrenzte Zeit in dem Kohlendioxidspei-cher zuruck gehalten werden kann. Ein großtechnischer Einsatz in Kraftwerken erscheintfruhestens in 10–20 Jahren moglich. Die Abtrennung vom Rauchgas kann mit unterschied-lichen Verfahren erfolgen, z.B. nach einer Kohlevergasung, Verbrennung in Sauerstoffatmo-sphare, oder CO2-Wasche aus dem Rauchgas. Als mogliche CO2-Lager gelten insbesonderegeologische Formationen wie ausgeraumte Erdol– und Erdgaslagerstatten sowie salzhaltigetiefe Grundwasserleiter (Aquifere). Auch eine Lagerung in der Tiefsee ist nicht ausgeschlos-sen (siehe Abbildung 5.1).

Mit diesem Thema beschaftigen zahlreiche Projekte. An der Technische Universitat Munchenwurde im Jahr 2009 ein ProjektKAUST−TUM, Simulating C02 Seguestration begonnen.Dieses Projekt ist ein Partnerschaft-Projekt mit der King Abdullah Technische Universitat(King Abdullah University of Science and Technology) in Saudi Arabia. Das Ziel dieses Pro-jekts ist, neuartige Untersuchung und neue Ansatze zum Modellieren und zur Simulation derCO2-Sequestration Prozessen, insbesondere im Rahmen der erhohten Olwiederaufnahme.Es ist eine Zusammenarbeit der Forschungsgruppen mit Sachkenntnise aus Bereichen, wieFlussphysik, mathematisches Modellieren, numerische Algorithmen, Optimierung und Da-tenverarbeitung. Das UnterprojektModelling of C02−Sequestration Including ParameterIdentification and Numerical Simulations wird beim Professor Dr. Dr.h.c.mult Karl-Heinz Hoffmann geleitet.In diesem Unterprojekt wird mathematisches Modellieren der CO2-Sequestration betrach-tet. Das mathematische Modellieren fuhrt zu Berucksichtigung von Mehrphasenstromungenin einem porosen Medium einschließlich spezielle Eigenschaften wie chemische Reaktionen,Hysteresisverhalten und Phasenubergange.

41


Speicherung

CO2

Abgase

Kohle Luft,Sauerstoff

BioGas

Oil

CO2

Sequestrierung

CO2unterirdische

CO2Transport

Pipeline, Schiff Verbrauch

Chemische Anlagen

der Miniralen

Karbonisation Industriell

AscheWärme Energieanlagen

Brenngebiete

geologische

Abbildung 5.1: Das Schema des Prozesses.

Im Modell, das in Rahmen der Diplomarbeit betrachtet wurde, existieren drei Phasen desKohlenstoffdioxides (siehe [15]):

• Das Wasser mit den gelosten Fremdstoffe (die mittlere Phase),

• Das geloste CO2 im Wasser (die schwere Phase),

• Das CO2 in einem kritischen Zustand, als das Gas und das Fluid gleichzeitig(die leichte Phase).

Es ist enorm wichtig den Prozess der Umwandlung einer Phase in die andere Phase zu ver-stehen. Das Verstehen hilft bei der Bestimmung eines Ortes, wohin den CO2 gelagert wird.

Allgemeiner Zugang zur solchen Probleme ist uber die Navier-Stokes Gleichungen. Wir wer-den aber im unserem Fall das Darcy-Gesetz betrachten. Das Darcy-Gesetz gibt die in einerZeiteinheit durch einen bestimmten Querschnitt eines Porengrundwasserleiters fliessendeWassermenge an. Fur den Anwendungsbereich des Darcy-Gesetzes muss die gemesseneTragheitskraft des stromenden Mediums gegenuber der Kraft der inneren Reibung ver-nachlassigbar klein sein. Sie gilt demnach nicht mehr, wenn die Stromungsgeschwindigkeitin den Porenraumen so gross wird, dass die Tragheitskraft merklich wird. Ein Mass fur dieobere Gultigkeitsgrenze stellt die Reynolds-Zahl dar. Das Darcy-Gesetz trifft dann zu, wennsich die Reynolds-Zahl im Bereich zwischen 1 und 10 befindet. In naturlichen Grundwas-serstromen werden Reynolds-Zahl gleich 10 in der Regel nicht uberschritten. Stromungen

42


mit hoheren Reynolds-Zahlen werden als turbulente Stromungen bezeichnet.

Wir werden annehmen, dass der Kapillardruck in unserem Modell nicht vorhanden ist undfur den allgemeinen Druck auf dem Rand eines Gebietes wird die Dirichlet-Bedingung be-tractet. Wir werden die Finite Methoden und Relaxation-Methoden fur die Losung dernichtlinearen Gleichungen anwenden. Die großte Rolle werden die Upwind-Approximationenbei der Losung der Erhaltung-Gleichungen spielen, d.h. die Interpolation von der Richtungdes Geschwindigkeitsvektor abhangig sein wird.

5.1 Mathematisches Model

In unserem Modell haben wir die folgenden Bezeichnungen eingefuhrt:

• φ ist die Porositat. Die Porositat ist als ein Verhaltniss zwischen des Volumens desFluides und dem gesamten Volumen definiert: φ = WF

W .

• Sα definiert die Saturation fur die Phase α. Die Saturation ist als ein Verhaltnisszwischen dem Volumen der Phase α und dem gesamten Volumen des Fluides definiert:Sα = Wα

WF.

• ρ0α(pα) ist die physikalische Dichte der Phase α. Die physikalische Dichte ist in der

Regel eine Funktion, die vom Druck abhangig ist.

• ρα ist die Dichte der Phase α und ρα=Sαρ0αφ.

• ~Vα ist die Geschwindigkeit der Phase α.

• p bezeichnet den Druck, der als gleichwertig fur jede Phase angenommen wird: pα = p.

Das Modell basiert auf den physikalischen Gesetzen: dem Prinzip der Massenerhaltungund dem Darcy-Gesetz. Das Gebiet, in dem ein Fluid betrachtet wird ist im Allgemeinenzeitabhangig. Das Gleichungssystem besteht aus der Phasen-Gleichungen und globalen Glei-chungen. Das Geschwindigkeitsfelde ~Vα(x, t) als vektorielle Große, ρα(x, t) und der Druckp(x, t) als skalare Großen sind dabei die Unbekannten des Systems. Die Raumkoordinate xliegt im R

3, also x = (x1, x2, x3). Um die Schreibweise zu erleichten, werden werden wir beider Großen des Systems die Abhangigkeit von der Variablen ab jetzt nicht mehr explizitschreiben.

5.1.1 Phasen-Gleichungen

Prinzip der Massenerhaltung

Wir betrachten zuerst das Prinzip der Massenerhaltung fur jede Phase:

∂ρα

∂t+ ∇ · (ρα

~Vα) = 0, (5.1)

wobei ρα die Dichte und Vα das Geschwindigkeitsfeld der Phase α sind.Diese Gleichung beschreibt, dass die gesamte Masse der Phase α konstant ist. Es gibt nureine Moglichkeit die gesamte Masse zu verandern, nahmlich mit einem physikalischen Flußdurch den Rand des Gebietes. Die Gleichung (5.1) kann auch modifiziert werden, um denFall der Formation einer Phase aus der andere zu berucksichtigen.

43


Darcy-Gesetz

Die nachste Gleichung stellt das Darcy-Gesetz fur jede Phase dar:

~Vα = −kα(∇p− ρ0α · ~g). (5.2)

Hier sind ~Vα das Geschwindigkeitsfeld mit der Diffusionsfahigkeit, ∇p der Gradient desDruckes, ρ0

α die physikalische Dichte der Phase α, ~g die Schwerkraft und kα der Koeffizientfur die Porositat der Phase α.Ursprunglich wurde das Darcy-Gesetz experimentell erworben, aber zur Zeit kann man dasGesetz aus dem Navier-Stokes-Gesetz mit der Homogenisierung herleiten. Um die Auftriebs-kraft in unserem Modell zu berucksichtigen, sollte die Darcy-Gleichung (5.2) umgeschriebenwerden:

~Vα = −kα(∇p− (ρ0α −

∑β ρ

0βkβρβ∑

γ kγργ) · ~g).

Das Vorzeichen vom Term (ρ0α −

P

β ρ0βkβρβ

P

γ kγργ) · ~g berucksichtigt das Verhalten der Phase α,

ob sie nach unten oder nach oben geht. Dieser Term beschreibt die Auftriebskraft.

Es sollte noch erwahnt werden, dass wir mit dem kalibrierten Druck p die nummerischenProbleme, die wegen der Eliminierung des Gravitationsterm in der nichtlinearen Gleichun-gen fur den Druck erscheinen, vermeiden werden.

Als nachstes diskutieren wir die globale Gleichungen des Systems.

5.1.2 Globale Gleichungen

Prinzip der Massenerhaltung

Zuerst wird das Prinzip der Massenerhaltung betrachtet. Die Gleichung ist fur die globaleDichte ρ sehr ahnlich wie die Gleichung (5.1):

∂ρ

∂t+ ∇ · (ρ~V ) = 0, (5.3)

wobei ρ =∑

α ρα ist und ~V ist die baryzentrische Geschwindigkeit, ~V = 1ρ

∑α ραVα.

Im Gegensatz zur Gleichung fur jede Phase, andert die gesamte Masse sich nicht beimPhasentransfer.

Gleichung fur den globalen Druck

Um die Gleichung fur den globalen Druck herzuleiten, werden wir mit den bereits erwahntenGesetzen manipulieren: Darcy-Gesetz und Massenerhaltung. Die Vorgehnsweise ist folgende:

• Zuerst multiplizieren wir die Gleichung (5.2) fur das Darcy-Gesetzt mit der Dichteρα:

ρα~Vα = −ραkα(∇p− (ρ0

α −∑

β ρ0βkβρβ∑

γ kγργ) · ~g),

44


oderρα~Vα = −ραkα(∇p− (ρ0

α −∑

β

ρ0βωβ) · ~g), (5.4)

mit ωβ :=kβρβ∑γ kγργ

,∑

β

ωβ = 1.

• Weiter summieren wir uber alle Phasen auf, d.h mit der Gleichung (5.4) bekommenwir:

∑

α

ρα~Vα =

∑

α

−ραkα∇p+ (∑

α

−ραkαρ0α) · ~g + (

∑

α

ραkα

∑

β

ρ0βωβ) · ~g,

oder

∑

α

ρα~Vα =

∑

α

−ραkα∇p+ (∑

α

−ραkαρ0α) · ~g + (

∑

β

ρ0βkβρβ) · ~g.

Nach der Vereinfachung erhalten wir:

ρ~V = −∑

α

ραkα∇p. (5.5)

• Im nachsten Schritt setzen wir die Gleichung (5.5) in die Gleichung fur die Massener-haltung (5.3) ein. Sei K :=

∑α ραkα, dann gilt:

∂ρ

∂t−∇ · (K∇p) = 0.

5.1.3 Gemeinsames System

Jetzt sind wir in der Lage das gemeinsame System der relevanten Gleichungen fur die Losungaufzustellen:

∂ρ

∂t− div(K∇p) = 0, (5.6)

~Vα = −kα(∇p− (ρ0α −

∑

β

ρ0βωβ) · ~g), (5.7)

∂ρα

∂t+ div(ρα

~Vα) = 0. (5.8)

In jedem Zeitpunkt werden wir aus der Gleichung (5.6) den globalen Druck des Systemsbestimmen und in die Gleichung (5.7) einsetzen. Auf der Weise erhalten wir das Geschwin-digkeitsfeld ~Vα fur jede Phase und setzen das in die Gleichung (5.8) ein, um die Dichte derPhase in jedem Zeitpunkt zu berechnen.

In nachstem Unterkapitel wird die genaue Simulationsmethode fur die Berechnung des Sys-tems prasentiert, die auf Finite-Differenz-Methode und Upwind-Approximationen basiert.

45


5.2 Simulation

Das matematische Modell wurde in Zusammenarbeit mit dem Austausch-Student Luiz Fa-ria aus der Universitat KAUST untersucht und simuliert (siehe [13]). Die Aufgabe dieserDiplomarbeit war die Simulation zu erweitern, in dem wir fur die Massenerhaltungsglei-chung (5.6) neuen Ansatz zur Losung vorgenommen haben. Der neue Ansatz basiert aufFinite-Volumen-Methode und Upwind-Approximationen-Methode.Wir werden jetzt ganzes Schema zur Losung des Systems (5.6),(5.6),(5.6) durchgehen.

5.2.1 Linienmethode zur Berechnung des Druckes

Zu Begin wird eine Methode zur Berechnung des Druckes betrachtet. Um die Gleichung(5.6) zu losen, werden die Zeitdiskretisierung und raumliche Diskretisierung hyperbolischerGleichung im Zusammenhang mit der Linienmethode vorgenommen. Darunter versteht mandie Diskretisierung von zeitabhangigen partiellen Differentialgleichungen in zwei Schritten,die klar voneinander getrennt werden: Diskretisierung der Ortsachse und Diskretisierung derZeitachse. Falls zuerst der Ort diskretisiert wird und dann die Zeit, spricht man von ver-tikaler Linienmethode, andernfalls von horizontaler Linienmethode (auch Rothe-Methodegenannt). Die vertikale Methode fuhrt auf ein System gewohnlicher Differentialgleichungenin der Zeit, die horizontale auf eine Folge stationarer Probleme im Ort. In unserem Fallwerden wir die horizontale Linienmethode anwenden.

Zeitdiskretisierung

Der globale Druck ist gegeben durch die Gleichung:

∂ρ

∂t+ div(K∇p) = 0, K =

∑

α

ραkα.

Nach der Diskretisierung der Ortsachse ergibt sich folgende Gleichung:

ρn+1 − ρn

τ− div(Kn+1∇pn+1) = 0, τ > 0,

oder

ρn+1 − τdiv(Kn+1∇pn+1) = ρn, τ > 0.

Sei ρ =∑

α ρα, dann erhalten wir:

∑

α

ρn+1α − τdiv(Kn+1∇pn+1) =

∑

α

ρnα

Danach wird die Definition fur die Dichte ρα angewendet:

ρα = Sαρ0α ≈ Sαρ

00α (1 + γp),

wobei ρ00 und γ konstant sind. Die letze Gleichung zeigt, dass wir die physikalische Dichteρ0

α als eine Funktion betrachten, die vom Druck abhangig ist. Nach der Substitution wirddie Gleichung fur den Druck folgendermaßen umgeschrieben:

46


∑

α

(Sn+1α ρ00

α γ)pn+1 − τdiv(Kn+1∇pn+1) =

∑

α

Snαρ

00α −

∑

α

Sn+1α ρ00

α . (5.9)

Die erworbene Gleichung ist nichtlinear in pn+1 und enthalt die unbekante Große Sn+1α .

Um die Gleichung (5.9) zu losen, wird ein Relaxation-Algorithmus angewandt. Bevor wirdas Relaxation-Algorithmus nahe betrachten, werden wir auf die Diskretisierung im Orteingehen.

Diskretisierung im Ort, Finite-Differenzen-Methode

Nach dem wir die Gleichung (5.9) erhalten haben, wird die Diskretisierung im Ort vorgenom-men. Wir werden uns dabei auf einen Fall mit zwei Dimensionen und einem Rechteckgebietbeschranken. Der dreidimensionale Fall ist auf die gleiche Weise erzielbar.

Der erste Schritt zur numerischen Losung auf einem Rechteckgebiet ist die Einfuhrung einesGitters, das wir hier nicht aquidistant in beiden Raumdimensionen wahlen, hx und hy sinddie Schrittweiten der Diskretisierung in x- und y-Richtung. Das Gitter besteht dabei aus(Nx + 2) × (Ny + 2) Punkten.Fur die Diskretisirung des Terms div(Kn+1∇pn+1) wird die Differenzen-Methode ange-wandt. Dabei setzen wir voraus, dass ein Koordinatensystem mit dem Ursprung links untendes rechteckigen Gebites eingefuhrt wurde. In jedem Gitterpunkt (xi, yj) ist eine Approxi-mation p(i, j) der exakten Losung p(xi, yj) gesucht. Wir gehen von Dirichlet-Randbedingungp = g auf ∂Ω aus. In diesem Fall sind die diskreten Werte p(i, j) auf dem Rand durch ein-fache Auswertung der Funktion g gegeben, p(i, j) = g(xi, yj) fur (xi, yj). Am Rand, wodie Phasen eingepresst werden, ist p(i, j) 6= 0. Das Koeffizient K nehmen wir abhangigvon (x, y), d.h. K(x, y). Dadurch konnten wir das Geschwindigkeitfeld in der Nahe vomRand beeinflußen. In unserem Fall wird das Geschwindigkeitsfeld an zwei Rande: j = 0 undj = Ny + 1 als Null angenommen. Sei nun (xi, yj) ein Punkt im Innern des Gebietes. Wirapproximieren div(K(x, y)∇pn+1) durch den diskreten Differenzenoperator.

Die Methode wird in drei Schritte zusammengefasst:

• Im ersten Schritt werden die Werte px und py fur den Gradient ∇pmittels der Zentral-

Differenzen approximiert. Sei hx

2 undhy

2 die Schrittweiten, x := xi und y := yj danngilt:

div(K(x, y)∇p) ≈

div

(K(x, y)

p(x+ hx

2 , y) + p(x− hx

2 , y)

hx,K(x, y)

p(x, y +hy

2 ) + p(x, y − hy

2 )

hy

),

• Im zweiten Schritt wird den Divergenzoperator mittels der Zentral-Differenzen appro-ximiert:

div(K(x, y)∇p) ≈K(x+ hx

2 , y) · (p(x+ hx, y) + p(x, y)) −K(x− hx

2 , y) · (p(x, y) − p(x− hx, y))

h2x

+

47


K(x, y +hy

2 ) · (p(x, y + hy) + p(x, y)) −K(x, y − hy

2 ) · (p(x, y) − p(x, y − hy))

h2x

,

• Im dritten Schritt fuhren wir ein paar Hilfskoeffizienten ein um den diskreten Diffe-renzenoperator zu prasentieren:

K1(i, j) :=K(i+ hx, j) +K(i, j)

2h2x

≈ K(i+ hx

2 , j)

h2x

,

K2(i, j) :=K(i, j + hy) +K(i, j)

2h2y

≈ K(i, j +hy

2 )

h2y

,

KK(i, j) := K1(i− 1, j) +K1(i, j) +K2(i, j − 1) +K2(i, j).

Dann ergibt sich den diskreten Differenzenoperator durch:

div(K(i, j)∇p) =

K1(i− 1, j)p(i− 1, j)+K1(i, j)p(i+1, j)+K2 (i, j− 1)p(i, j − 1)+K2(i, j)p(i, j +1)−−KK(i, j)p(i, j). (5.10)

p(i, j)

p(i, j + 1)

p(i, j − 1)

p(i− 1, j) p(i+ 1, j)

−KK(i, j)

K2(i, j)

K2(i, j − 1)

K1(i, j)K1(i− 1, j)

Abbildung 5.2: Modifizierte Funf-Punkte-Stern.

Fur jeden Gitterpunkt des Gebietes wird die lineare Gleichung gelost, um den Druckzu bestimmen. Deshalb ergibt sich das System der Nx × Ny-Gleichungen mit Nx × Ny-Unbekannten. Wenn wir die Unbekanten in einen Vektor ph lexikographish einordnen, lasstsich das ganze System schreiben als

48


Ahph = RechteSeite.

In der lexikographischen Anordnung der Unbekannten hat die Matrix Ah Blocktridiagonalgestalt.Zum Beispiel fur den Fall N := Nx = Ny wird die Matrix folgendermaßen darstellbar:

Ah =

B1 I1 0I1 B2 I2

I2. . .

. . .. . . BN−1 IN−1

0 IN−1 BN

mit

Bl =

−KK(l, 1) K2(l, 1) 0K2(l, 1) −KK(l, 2) K2(l, 2)

K2(l, 2). . .

. . .. . . −KK(l,N − 1) K2(l,N − 1)

0 K2(l,N − 1) −KK(l,N)

und

Il =

K1(l, 1) 0K1(l, 2)

. . .

K1(l,N − 1)0 K1(l,N)

, l = 1..N.

Dieses System der linearen Gleichungen wird mit einem Programm spoolesa gelost. Die-ses Programm ist ein Teil der FeliCs-Software. Der Solver FeliCs wurde speziell fur dieTechnische Universitat Munchen entwickelt. Das Prinzip des Programms ist die Zerlegungder Matrix (QR-Zerlegung fur die quadratischen Matrizen). Diese Methode ermoglicht dieBerechnung des Systems fur die große Anzahl der Unbekannten in optimaler Zeit. Die Be-rechnung wird fur die dunnbesetzten Matrizen vorgenommen. Es ist leicht zu sehen, dassdie Matrix in unserem Fall dunnbesetzt ist.

Als nachstes wird das Relaxation-Algorithmus zur Losung des gemeinsamen Systems be-trachtet.

5.2.2 Relaxation-Algorithmus zur Losung des gemeinsamen Systems,Upwind-Approximationen

Die Relaxation-Methode ist eine Iteration-Methode, die durch Vergleich alter Funktions-werte in jedem Schritt zur gesuchten numerischen Losung kommt. Es wird angenommen,dass die Losung zur Zeit t = n bekannt ist. In unserem Fall wird angenommen, dass die

49


Sattigungswerte(Saturation) und den Druck zur Zeit t = n bekannt sind und die Werte zurZeit t = n+ 1 werden mit einem Relaxationparameter m gesucht. Das bedeutet, dass zumAnfang des Verfahrens die Initialwerte bekannt sind:

Sn+1,mα = Sn

α, pn+1,m = pn, m = 0.

Druckgleichung

In jedem Iterationsschritt wird die Gleichung (5.6) mit dem Relaxationsparameter m gelost:

∑

α

(Sn+1,mα ρ00

α γ)pn+1,m+1 − τdiv(Kn+1,m∇pn+1,m+1) =

∑

α

Sn,mα ρ0n

α −∑

α

Sn+1,mα ρ00

α ,

(5.11)wobei die Werten fur die Sattigung(Saturation) Sn+1,m

α und den Druck pn+1,m bekannt sind.Mit den neuen Bezeichnungen

C1(x, y) :=∑

α

(Sn+1,mα ρ00

α γ)

undC2(x, y) :=

∑

α

Sn,mα ρ0n

α −∑

α

Sn+1,mα ρ00

α

wird die Gleichung (5.11) zu:

C1(x, y)pn+1,m+1 − τdiv(K(x, y)∇pn+1,m+1) = C2(x, y).

Als nachstes wird die Finite-Differenzen-Methode mit dem diskreten Differenzenoperator(5.10) angewandt. Es ergibt sich das System der linearen Gleichungen:

C1(i, j)pn+1,m+1 − τ(K1(i− 1, j)p(i − 1, j) +K1(i, j)p(i + 1, j)

+K2(i, j − 1)p(i, j − 1) +K2(i, j)p(i, j + 1) −KK(i, j)p(i, j)) = C2(i, j),

Der Index (i, j) lauft uber alle Gitterpunkte des recheckigen Gebietes Ω.Die Losung des Matrix-Systems ergibt den Wert fur den Druck mit dem Relsaxationspara-meter m+ 1: pn+1,m+1.

Geschwindigkeitsfeldgleichung

Nachdem der Wert pn+1,m+1 bekannt ist, wird das Geschwindigkeitsfeld fur jede Phase mitdem Darcy-Gesetz bestimmt. Es wird im m+ 1-Relaxationschritt die folgende Gleichunggelost:

~Vα = −kα(∇p− (ρ0α −

∑β ρ

0βkβρβ∑

γ kγργ) · ~g),

wobei p = pn+1,m+1 ist.

50


Massenerhaltungsgleichung, Upwind-Approximation-Verfahren

Nachdem die Geschwindigkeit fur jede Phase bestimmt wurde, konnte die Gleichung derMassenerhaltung nach der Dichte fur jede Phase aufgelost werden. Dabei wird die Finite-Volumen-Methode angewandt. Das Finite-Volumen-Verfahren ist ein numerisches Verfahrenzur Losung von Erhaltungsgleichungen, also von speziellen, haufig hyperbolischen, partiel-len Differentialgleichungen, denen ein Erhaltungssatz zugrunde liegt. Das Verfahren benutztin seiner Herleitung eine integrale Form der Erhaltungsgleichungen und erlaubt damit auchunstetige Losungen, die fur solche Gleichungen typisch sind. Ferner werden nur geringeAnforderungen an die Gitterzellen gestellt, was unstrukturierte und flexible Geometrienerlaubt. Daruber hinaus werden die konservativen Großen der Gleichung tatsachlich er-halten. Die Dichten mit dem Relaxationsparameter m+ 1 werden danach durch Upwind-Approximationen bestimmt.Also, wir integrieren die Gleichung (5.1) uber das ganze Gebiet Ω:

∂

∂t

∫

Ωρα +

∫

Ωdiv(ρα

~Vα) = 0,

wobei ρα = ρn+1,m+1α ist.

Mit dem Gauß’schen Integrasatz ergibt sich:

∂

∂t

∫

Ωρα +

∫

∂Ω

~Vαn · ρα = 0. (5.12)

Das letzte Integrall wird in vier Integrallen auf dem Gitter zerfallen. Nach der Zeit- undOrtdiskretisierung betrachten wir die Gleichung (5.12) in einem Gitterpunkt (i, j) mit derUmgebung [i− 1/2, i+ 1/2]× [j − 1/2, j + 1/2] (siehe Abbildung 5.3), dann bekommen wirfolgende Gleichung:

hxhy ·pn+1(i, j) − pn(i, j)

τ+

hx ·(V x

α (i+1

2, j)ρn

α(i+1

2, j) − V x

α (i− 1

2, j)ρn

α(i− 1

2, j)

)+ (5.13)

hy ·(V y

α (i, j +1

2)ρn

α(i, j +1

2) − V y

α (i, j − 1

2)ρn

α(i, j − 1

2)

)= 0.

Um die Geschwindigkeit V xα (i + 1

2 , j) zu bestimmen, wird die folgende Approximation ver-wendet:

V xα (i+

1

2, j) :=

V xα (i+ 1, j) + V x

α (i, j)

2≈ O(h2

x).

Mit der Taylor-Entwicklung kann man die Ordnung O(h2x) nachweisen.

Um die Geschwindigkeit V xα (i − 1

2 , j) zu bestimmen, wird die andere Approximation ver-wendet:

V xα (i− 1

2, j) :=

V xα (i− 1, j) + V x

α (i, j)

2≈ O(h2

x).

Die Geschwindigkeitswerte V yα (i, j + 1

2), V yα (i, j − 1

2) werden analog berechnen.Mit diesen berechneten Werten fur das Geschwindigkeitsfeld wird die entsprechende Dichtein den Punkten des Gitters i − 1 oder i oder i + 1 ausgewahlt, um die Berechnung fort-zusetzen. Die Methode, die wir hier anwenden, ist die Upwind-Approximation-Methode.

51


j + 1

j

j − 1

i− 1 i i+ 1

Abbildung 5.3: Gitter fur Upwind-Approximationen.

Abhanging von der Richtung des Geschwindigkeitsfeldes (vom Vorzeichen der Approxima-tionswerte) V x

α (i + 12 , j), V

xα (i − 1

2 , j), Vyα (i, j + 1

2 ), V yα (i, j − 1

2 ) wird die Dichte in bereitsbekannten Punkten des Gitter verwendet um die Werte in neuen Punkten zu berechnen.

Wir werden zuerst die x-Richtung auf den Gitter betrachten.Sei V x

α (i+ 12 , j) bekannt, dann

ρnα(i+

1

2, j) =

ρn

α(i, j), falls V xα (i+ 1

2 , j) ≥ 0ρn

α(i+ 1, j), falls V xα (i+ 1

2 , j) < 0(5.14)

Sei V xα (i− 1

2 , j) bekannt, dann

ρnα(i− 1

2, j) =

ρn

α(i− 1, j), falls V xα (i− 1

2 , j) ≥ 0ρn

α(i, j), falls V xα (i− 1

2 , j) < 0(5.15)

Die y-Richtung auf den Gitter wird analog approximiert.Sei V y

α (i, j + 12) bekannt, dann

ρnα(i, j +

1

2) =

ρn

α(i, j), falls V xα (i, j + 1

2) ≥ 0ρn

α(i, j + 1), falls V xα (i, j + 1

2 ) < 0(5.16)

Sei V xα (i, j − 1

2) bekannt, dann

ρnα(i, j − 1

2) =

ρn

α(i, j − 1), falls V xα (i, j − 1

2 ) ≥ 0ρn

α(i, j), falls V xα (i, j − 1

2) < 0(5.17)

Mit den erworbenen Upwind-Approximationen wird die Dichte ρα = ρn+1,m+1α fur jede

Phase aus der Gleichung (5.13) besimmt.Es bleibt noch zu zeigen, dass das vorgeschlagene Schema die Monotonie-Eigenschaft besitzt.

Behauptung 5.2.1 Das Upwind-Schema (5.13)-(5.17) ist monoton.

52


Beweis. Wir betrachten das vereinfachte eindimensionale Upwind-Schema:

ρn+1α (i) = ρn

α(i) − τ

h

(V 1

α (i)ρnα(i+

1

2) − V 2

α (i)ρnα(i− 1

2)

), (5.18)

wobei

V 1α (i) :=

Vα(i+ 1) + Vα(i)

2und

V 2α (i) :=

Vα(i− 1) + Vα(i)

2sind.

Sei ρα,1 und ρα,2 zwei Losungsfunktionen mit der Eigenschaft

ρnα,1 ≤ ρn

α,2.

Es ist zu zeigen dass aus der Monotonie zur Zeit n die Monotonie zur Zeit n+ 1 folgt, d.h.

ρnα,1 ≤ ρn

α,2 =⇒ ρn+1α,1 ≤ ρn+1

α,2 .

Um dies zu beweisen, betrachten wir die Differenz:

ρn+1α,1 (i) − ρn+1

α,2 (i) = ρnα,1(i) − ρn

α,2(i)−

τ

h

(V 1

α (i)(ρnα,1(i+

1

2) − ρn

α,2(i+1

2)) − V 2

α (i)(ρnα,1(i−

1

2) − ρn

α,2(i−1

2))

). (5.19)

In nachsten Schriit wird die Fallunterscheidung bezuglich der Approximationen V 1α und V 2

α

durchgefuhrt:

• Sei V 1α (i) ≤ 0 und V 2

α (i) ≥ 0, dann ist die rechte Seite der Gleichung (5.19) immernegativ:

ρn+1α,1 (i) − ρn+1

α,2 (i) = ρnα,1(i) − ρn

α,2(i)−τ

h

(V 1

α (i)(ρnα,1(i+ 1) − ρn

α,2(i+ 1)) − V 2α (i)(ρn

α,1(i− 1) − ρnα,2(i− 1))

)≤ 0

Es folgt aus:(V 1

α (i)(ρnα,1(i+ 1) − ρn

α,2(i+ 1)) − V 2α (i)(ρn

α,1(i− 1) − ρnα,2(i− 1))

)≥ 0

mit ρnα,1(i) ≤ ρn

α,2(i), ρnα,1(i− 1) ≤ ρn

α,2(i− 1), ρnα,1(i+ 1) ≤ ρn

α,2(i+ 1).

• Sei V 1α (i) ≥ 0 und V 2

α (i) ≤ 0, dann ist die rechte Seite der Gleichung (5.19) auchnegativ, falls 1 − τ

h ≥ 0, τ ist klein genug.

ρn+1α,1 (i) − ρn+1

α,2 (i) = ρnα,1(i) − ρn

α,2(i)−τ

h

(V 1

α (i)(ρnα,1(i) − ρn

α,2(i)) − V 2α (i)(ρn

α,1(i) − ρnα,2(i))

)=

(1 − τ

h)(ρn

α,1(i) − ρnα,2(i)) |Vα(i)| ≤ 0,

da ρnα,1(i) ≤ ρn

α,2(i) ist.

53


Die Falle V 1α (i) ≤ 0 und V 2

α (i) ≤ 0 oder V 1α (i) ≥ 0 und V 2

α (i) ≥ 0 folgen aus bewiesenenFallen automatisch. Also, wir haben gezeigt:

ρn+1α,1 ≤ ρn+1

α,2 .

Die Monotonie des prasentierten Upwind-Schemas wurde bewiesen. Um die Stabilitat desSchemas zu garantieren, d.h.

τ < min

(min (hx, hy)

max (Vα)√

2,

1

maxdiv(Vα)

),

wird die Zeitvariable τ klein genug ausgewahlt.

Bestimmung der Sattigung

Wenn die Dichte fur jede Phase bestimmt wurde, wird die Sattigung mit dem Relaxations-parameter m+ 1 fur jede Phase bestimmt

Sn+1,m+1α =

ρn+1,m+1α

ρ0,n+1α

.

Danach werden die berechneten Sattigungswerte skaliert um der Bedingung zu genugen:∑Sα = 1. Wir berechnen deshalb das folgende:

Sn+1,m+1α =

Sn+1,m+1α∑

β Sn+1,m+1β

Zum Schluß wird das Maximum der Differenzen zwischen alter Funktionswerte Sn+1,mα und

den neuen Sn+1,m+1α betrachtet. Die Relaxationiteration wird solange wiederholt, d.h m =

m+ 1 , bis die Sattigungswerte konvergieren. Nach dem Konvergenz wird die Berechnungzum nachsten Zeitpunkt starten: n = n+ 1.

Die Ergebnise der Berechnung nach unserer Methode ist auf den Abbildungen 5.4, 5.5, 5.6zu sehen. Die Schwerkraft ist dabei entlang der y-Achse, die z-Achse zeigt die Werte derSattigung fur jede Phase von 0 bis 1. Die leichte Phase ist in grun, die schwere Phase ist in rotund die mittlere Phase ist in blau dargestellt. Die Anfangswerte fur die leichte und schwerePhase sind ungleich Null. Die mittlere Phase ist zum Anfanf ist fast Null angenommen.Wir pressen dabei die schwere Phase am Rand mit x = 0 ein. Die Simulationen wurden aneinem LRZ-Rechner durch paralleles Rechnen berechnet. Die Parallelitat wurde mit einemComputer mit 30 CPU’s realisiert.

54


X Y

Z

X Y

Z

X Y

Z

X Y

Z

Abbildung 5.4: Zeitentwicklung der Sattigungen S1, S2, S3, S1 ist die leichte Phase, S2 istdie mittlere Phase, S3 ist die schwere Phase

55


X Y

Z

X Y

Z

X Y

Z

X Y

Z


56


X Y

Z

X Y

Z

X Y

Z


57

6 Zusammenfassung

58

Literaturverzeichnis

[1] M. Bardi. Some applications of viscosity solutions to optimal control and differentialgames. Viscosity Solutions and Applications, Lecture Notes in Mathematic, Springer,1660:44–97, 1997.

[2] M. Bardi and I. Capuzzo-Dolcetta. Optimal Control and Viscosity Solutions ofHamilton-Jacobi-Bellman equations. Birkhauser, Boston, 1997.

[3] M. Bardi, M. Falcone, and P. Soravia. Fully-discrete schemes for the value function ofpursuit-evasion games. In T. Basar and A. Hauri (eds.), Advances in Dynamic Gamesand Applications, Annals. of the Int. Society of Dynamic Games I,Birkhauser, Boston,pages 89–105, 1994.

[4] G. Barles. Remarks on a flame propagation model. Rapports de Recherche, 451, No-vember 1985.

[5] N.D. Botkin. Approximation schemes for finding the value functions for differentialgames with nonterminal payoff functional. Analysis, 14:203–220, 1994.

[6] N.D. Botkin, K.-H. Hoffmann, and V.L. Turova. Stable solutions of Hamilton-Jacobiequations. application to control of freezing processes. Preprint-Number SPP1253-080,http://www.am.uni-erlangen.de/home/spp1253, September 2009.

[7] G. Caginalp. An analysis of a phase field model of a free boundary. Arch Rat. Mech.Anal., 92:205–245, 1986.

[8] G. Caginalp and X. Chen. Convergence of the phase field model to its sharp interfacelimits. Euro Jnl of Applied Mathematics, 9:417–445, 1998.

[9] M.G. Crandall, L.C. Evans, and P.L. Lions. Some properties of viscosity solutions ofHamilton-Jacobi equations. Trans. Amer. Math. Soc., 282:487–502, 1984.

[10] M.G. Crandall and P.L. Lions. Viscosity solutions of Hamilton-Jacobi equations.Trans. Amer. Math. Soc., 277:1–47, 1983.

[11] M.G. Crandall and P.L. Lions. Two approximations of solutions of Hamilton-Jacobiequations. Math. Comp., 43:1–19, 1984.

[12] E. Cristiani and M. Falcone. Fully-discrete schemes for value function of pursuit-evasion games with state constraints. In P. Bernhard, V. Gaitsgory, and O. Pourtallier(eds.), Advances in Dynamic games and Their Applications, Annals. of the Int. Societyof Dynamic Games X, Birkhauser, Boston, pages 177–206, 2009.

[13] Luiz Faria. Simplified model for multiphase flow. Technical report, TUM, M6, 2010.

[14] M. Fremond. Non-smooth thermomechanics. Springer-Verlag, Berlin, 2002.

59

6 Zusammenfassung

[15] K.-H. Hoffmann and N. Botkin. Multiphase flows with phase transitions and chemicalreactions: Models and numerics. Preprint submited to Math. and Comp. in Simulation(IMACS).

[16] K.-H. Hoffmann, N. Botkin, and V. Turova. Projekt: Optimal control in cryopreserva-tion of cells and tissues. http://www-m6.ma.tum.de/ botkin/spp1253.html.

[17] K.-H. Hoffmann and N.D Botkin. Optimal control in cryopreservation of cells andtissues. Advances in Math. Sciences and Appl., 29:177–200, 2008.

[18] K.-H. Hoffmann and Jiang Lishang. Optimal control of a phase field model for solidi-fication. Numer. Funct. Anal. Optimiz., 13:11–27, 1992.

[19] R. Isaacs. Differential Games. New York, John Wiley, 1965.

[20] N.N. Krasovskii and A.I. Subbotin. Game-Theoretical Control Problems. New York:Springer-Verlag, 1988.

[21] O.A. Malafeyev and M.S. Troeva. A weak solution of Hamilton-Jacobi equation for adifferential two-person zero-sum game. in Preprints of the Eighth International Sympo-sium on Differential Games and Applications, Maastricht, Netherlands, pages 366–369,1998.

[22] A. Melikyan, N.D. Botkin, and V.L. Turova. Propagation of disturbances in inhomoge-neous anisotopic media. Preprint N181, SFB611, Singulare Phanomene und Skalierungin mathematischen Modellen, Bonn, pages 1–25, 2004.

[23] A.W. Merz. The game of two identical cars. JOTA, 9(5):324–343, 1972.

[24] I.M. Mitchell. Application of level set methods to control and reachability problems incontinuous and hybrid systems. PhD thesis, Standford University, Scientific Computingand Computational Mathematics Program, 2002.

[25] P.E. Souganidis. Approximation schemes for viscosity solutions of Hamilton-Jacobiequations. Journal of Differential Equations, 59:1–43, 1985.

[26] P.E. Souganidis. Max - min representation and product formulas for the viscosity solu-tions of Hamilton-Jacobi equations with applications to differential games. NonlinearAnalysis, Theory, Methods and Applications, 9:217–257, 1985.

[27] A.I. Subbotin. A generalization of the basic equation of the theory of differential games.Dokl. Akad. Nauk SSSR, 254:293–297, 1980.

[28] A.I. Subbotin. Generalization of the main equation of differential game theory. J.Optimiz. Theory and Appl., 43:103–133, 1984.

[29] A.I. Subbotin and A.G. Chentsov. Optimization of Guaranteed Result in Control Pro-blems. Moscow, Nauka (in Russian), 1981.

[30] A.I. Subbotin and A.M. Taras’yev. Stability properties of the value function of adifferential game and viscosity solutions of Hamilton-Jacobi equations. Problems ofControl and Information Theory, 15:451–463, 1986.

60

Theoretische und numerische Untersuchung einiger …botkin/all_kaust_files/papers/...November 2010...

Documents

Transcript of Theoretische und numerische Untersuchung einiger …botkin/all_kaust_files/papers/...November 2010...