Numerische Verfahren f¨ur Differentialgleichungen

Numerische Verfahren f¨ ur Differentialgleichungen

Steffen B¨ orm

Stand 1. Juli 2014

Alle Rechte beim Autor.

Inhaltsverzeichnis

1 Einleitung 7

1.1 Federpendel . . . . 7

1.2 Mehrk¨ orperprobleme . . . . 10

1.3 Schwingende Saite . . . . 11

1.4 W¨ armeleitung . . . . 14

2 Einige theoretische Aussagen ¨ uber gew¨ ohnliche Differentialgleichungen 17 2.1 Allgemeine Problemstellung . . . . 17

2.2 Existenz und Eindeutigkeit . . . . 18

2.3 St¨ orungen der Daten . . . . 22

3 Einschrittverfahren 25 3.1 Euler-Verfahren . . . . 25

3.2 Konvergenz . . . . 29

3.3 Konsistenz . . . . 37

3.4 Lokalisierte Konvergenzaussagen . . . . 41

3.5 Konsistenzkriterium . . . . 45

3.6 Runge-Kutta-Verfahren . . . . 50

4 Verfeinerte Techniken f¨ ur gew¨ ohnliche Differentialgleichungen 57 4.1 Extrapolation . . . . 57

4.2 Schrittweitensteuerung . . . . 61

4.3 Steife Differentialgleichungen . . . . 67

4.4 Differential-algebraische Gleichungen . . . . 73

5 Beispiele f¨ ur partielle Differentialgleichungen 79 5.1 Hyperbolische Gleichungen und das Verfahren der Charakteristiken . . . . 80

5.2 Elliptische Differentialgleichungen und das Finite-Differenzen-Verfahren . 84 5.3 Parabolische Differentialgleichungen und die Linienmethode . . . . 91

6 Variationsformulierungen und das Finite-Elemente-Verfahren 97 6.1 Variationsformulierung . . . . 97

6.2 Sobolew-R¨ aume . . . 101

6.3 Existenz und Eindeutigkeit . . . 106

6.4 Galerkin-Verfahren . . . 115

6.5 Interpretation als Minimierungsproblem . . . 118

6.6 Eindimensionale finite Elemente . . . 119

6.7 Mehrdimensionale finite Elemente . . . 124 6.8 Analyse des Approximationsfehlers . . . 130 7 L¨ osungsverfahren f¨ ur schwachbesetzte Matrizen 137 7.1 Gradientenverfahren . . . 137 7.2 Verfahren der konjugierten Gradienten . . . 144

Index 149

Literaturverzeichnis 151

Vorwort

Differentialgleichungen sind ein wichtiges Werkzeug bei der Behandlung vieler mathema- tischer, physikalischer, biologischer, chemischer oder wirtschaftswissenschaftlicher Fra- gestellungen.

• In der Mathematik lassen sich beispielsweise mit ihrer Hilfe Geod¨ aten beschreiben, die k¨ urzesten Verbindungen zwischen zwei Punkten einer gekr¨ ummten Oberfl¨ ache.

• In der Physik werden beispielsweise die Bewegung von K¨ orpern, die Ausbreitung elektromagnetischer Wellen und die Ausbreitung von W¨ arme mit Hilfe von Diffe- rentialgleichungen beschrieben.

• In der Biologie werden Differentialgleichungen beispielsweise eingesetzt, um die zeitliche Entwicklung von Bakterienpopulationen zu modellieren.

• In der Chemie l¨ asst sich mit Hilfe von Differentialgleichungen die Dynamik von Molek¨ ulen beschreiben.

• In den Wirtschaftswissenschaften kommen Differentialgleichungen bei der Simula- tion der Entwicklung von Kursen am Aktien- oder Geldmarkt zum Einsatz.

Den meisten dieser Gleichungen ist gemeinsam, dass sie sich in der Regel nicht

” mit Pa- pier und Bleistift“ l¨ osen lassen, mit den Methoden der reinen Mathematik k¨ onnen h¨ aufig nur R¨ uckschl¨ usse auf das Verhalten der L¨ osung ziehen, aber nur in seltenen Spezialf¨ allen die L¨ osung angeben.

Deshalb kommen numerische L¨ osungsverfahren zum Einsatz. Mit Hilfe dieser Ver- fahren kann die exakte L¨ osung zwar auch nicht berechnet werden, aber sie l¨ asst sich beliebig genau approximieren, und f¨ ur die meisten praktischen Anwendungen gen¨ ugt eine hinreichend gute N¨ aherung der L¨ osung.

Diese Vorlesung gibt einen ¨ Uberblick ¨ uber einige der wichtigsten Verfahren f¨ ur die Approximation der L¨ osungen von Differentialgleichungen. Besonders eingehend behan- delt werden dabei gew¨ ohnliche Differentialgleichungen, die vor allem bei der Simulation zeitabh¨ angiger Ph¨ anomene eine wichtige Rolle spielen, und elliptische partielle Diffe- rentialgleichungen, mit denen sich Kraft- und Spannungsfelder beschreiben lassen, bei- spielsweise bei der Untersuchung elektromagnetischer Felder oder strukturmechanischer Fragestellungen.

Zwei weitere Typen partieller Differentialgleichungen, hyperbolische und parabolische

Gleichungen, werden nur am Rande behandelt. Parabolische Gleichungen lassen sich

allerdings als Kombination gew¨ ohnlicher und elliptischer Gleichungen interpretieren und

deshalb mit den hier beschriebenen Verfahren behandeln.

Die Behandlung hyperbolischer Differentialgleichungen erfordert dagegen h¨ aufig spe- zialisierte L¨ osungsverfahren, die den Rahmen dieser Vorlesung sprengen w¨ urden. Deshalb wird lediglich einer der einfachsten L¨ osungsans¨ atze anhand eines einfachen Beispielpro- blems demonstriert und die Diskussion allgemeinerer Verfahren spezialisierten Vorle- sungen ¨ uberlassen, beispielsweise aus dem Bereich der Erhaltungsgleichungen oder der Str¨ omungsmechanik.

Die Vorlesung setzt Kenntnisse der Analysis und der linearen Algebra voraus: Aus der Analysis sollten neben Grundbegriffen (Konvergenz von Folgen und Reihen, Ste- tigkeit) nat¨ urlich grundlegende S¨ atze der Differential- und Integralrechnung bekannt sein (Hauptsatz, Mittelwerts¨ atze). Aus der linearen Algebra sollten Kenntnisse ¨ uber das Rechnen in Vektorr¨ aumen und den Umgang mit Skalarprodukten mitgebracht werden.

Danksagung

Ich bedanke mich bei Janina Gnutzmann, Hendrik Felix Pohl und Sven Christophersen

f¨ ur Korrekturen und Verbesserungsvorschl¨ age.

1 Einleitung

Bevor wir uns der Analyse und numerische Behandlung von gew¨ ohnlichen Differential- gleichungen, insbesondere von Anfangswertproblemen, zuwenden, sollen zun¨ achst einige mehr oder weniger einfache Probleme vorgestellt werden, die sich mit Hilfe derartiger Gleichungen beschreiben lassen.

1.1 Federpendel

Ein sehr einfaches Beispiel f¨ ur ein Anfangswertproblem ist das abstrakte Federpendel:

Es besteht aus einer Masse m, die mittels einer Feder mit einem festen Punkt verbunden ist und nach oben oder unten ausgelenkt werden kann. Die Auslenkung aus der Ruhelage zu einem bestimmten Zeitpunkt t bezeichnen wir mit u(t).

Falls die Masse sich nicht im Nullpunkt befindet, ist die Feder angespannt und ¨ ubt eine Kraft aus, die die Masse in den Nullpunkt zur¨ uckzieht. Im einfachten Fall ist diese Kraft F (t) durch das Hookesche Gesetz [7]

F (t) = − c ˆ

` u(t) (1.1)

gegeben, wobei ˆ c eine vom Material der Feder abh¨ angende Konstante und ` die L¨ ange der Feder im Ruhezustand ist.

u(t)

Abbildung 1.1: Modell eines Federpendels

Gem¨ aß der Newtonschen Axiome [8] bewirkt die Kraft F eine Beschleunigung a(t) der Masse, die proportional zum Kehrwert von m ist, es gilt also

F (t) = ma(t), a(t) = 1

m F(t).

Die Beschleunigung ist die Ableitung der Geschwindigkeit v(t) der Masse, und die Ge- schwindigkeit ist die Ableitung der Auslenkung u(t), so dass wir die Gleichungen

u

(t) = v(t), v

(t) = a(t) = 1

m F (t) = − ˆ c m` u(t) erhalten. Indem wir λ := ˆ c/(m`) einf¨ uhren und u(t) und v(t) zu einem Vektor

y(t) = u(t)

v(t)

zusammenfassen, erhalten wir die kompakte Schreibweise y

(t) =

0 1

−λ 0

y(t),

mit der die Bewegungen des abstrakten Pendels vollst¨ andig beschrieben werden k¨ onnen.

Wenn die Auslenkung u(0) und die Geschwindigkeit v(0) zum Startzeitpunkt bekannt sind, k¨ onnen wir Auslenkung und Geschwindigkeit zu jedem sp¨ ateren Zeitpunkt t ≥ 0 als L¨ osung des Systems

y(0) = u(0)

v(0)

, y

(t) =

0 1

−λ 0

y(t) f¨ ur alle t ∈ R

(1.2) bestimmen. Ein derartiges Gleichungssystem, bei dem die L¨ osung zu einem Anfangs- zeitpunkt und die Ableitung der L¨ osung zu jedem Zeitpunkt t bekannt sind, nennt man Anfangswertproblem.

In unserem Fall haben wir es mit einem besonders einfachen System zu tun, dass sich analytisch l¨ osen l¨ asst: Wir f¨ uhren die Matrix

A :=

0 1

−λ 0

ein und erhalten

y

(t) = Ay f¨ ur alle t ∈ R

,

ein lineares Anfangswertproblem. Wenn A eine Zahl w¨ are, k¨ onnten wir t 7→ α exp(At) als L¨ osungsansatz verwenden, wobei exp(x) = e

die Exponentialfunktion bezeichnet.

Dann m¨ ussen wir nur α so bestimmen, dass die Anfangsbedingung erf¨ ullt ist.

Da A in unserem Fall eine Matrix ist, bietet es sich an, nach einer Verallgemeinerung der Exponentialfunktion zu suchen. Die naheliegende Definition

exp(C) :=

∞

X

j=0

C

j! (1.3)

1.1 Federpendel f¨ ur eine Matrix C ∈ R

beruht auf der Exponentialreihe. F¨ ur eine beliebige induzierte Matrixnorm gilt

k exp(C)k =

∞

X

j=0

C

j!

≤

∞

X

j=0

C

j!

≤

∞

X

j=0

kCk

j! = exp(kCk),

also ist die Reihe absolut konvergent, also insbesondere konvergent, und die Matrix- Exponentialfunktion damit durch (1.3) wohldefiniert.

Im Kontext der Anfangswertprobleme sind wir an der Ableitung der Funktion t 7→

exp(tA) interessiert, die durch

∂

∂t exp(tA) = ∂

∂t

∞

X

j=0

t

A

j! =

∞

X

j=0

jt

A

j! =

∞

X

j=1

t

A

(j − 1)!

= A

∞

X

j=1

t

A

(j − 1)! = A exp(tA) f¨ ur alle t ∈ R gegeben ist. Da aus (1.3) auch exp(0) = I folgt, erf¨ ullt die durch

y(t) := exp(tA)y(0) f¨ ur alle t ∈ R

definierte Funktion gerade

y(t

) = exp(0A)y(0) = y(0),

y

(t) = A exp(tA)y(0) = Ay(t) f¨ ur alle t ∈ R

, ist also eine L¨ osung des linearen Anfangswertproblems.

Statt durch direkte Auswertung der Exponentialsumme k¨ onnen wir die Exponential- funktion auch berechnen, indem wir A mit einer ¨ Ahnlichkeitstransformation diagonali- sieren: Mit der Matrix

T :=

1 1 i √

λ −i √ λ

, T

= 1

2

1 −i/ √ λ 1 i/ √

λ

, erhalten wir

T

AT = i √

λ

−i √ λ

=: D und k¨ onnen die Exponentialfunktion durch

exp(tA) =

∞

X

i=0

t

A

i! =

∞

X

i=0

t

(TDT

)

i! = T

∞

X

i=0

t

D

i!

! T

= T P

i=0 (ti√

λ)ⁱ i!

P

(−ti√ λ)ⁱ i!

!

T

= T

exp(ti √ λ)

exp(−ti √ λ)

T

auf die Berechnung der Exponentialfunktion f¨ ur skalare Werte zur¨ uckf¨ uhren.

Die Exponentialfunktion eines rein imagin¨ aren Werts steht in enger Beziehung zu Sinus- und Cosinus-Funktionen, deshalb ¨ uberrascht es nicht, dass wir auch direkt den Ansatz

u(t) = α sin(ωt) + β cos(ωt),

v(t) = αω cos(ωt) − βω sin(ωt) f¨ ur alle t ∈ R

verwenden k¨ onnen. Der Parameter ω = √

λ = p ˆ

c/(m`) beschreibt die Frequenz der Schwingung in Abh¨ angigkeit von Masse, Materialeigenschaften und Federl¨ ange, w¨ ahrend die Parameter α und β verwendet werden k¨ onnen, um sicherzustellen, dass die Anfangs- bedingungen erf¨ ullt sind.

1.2 Mehrk¨ orperprobleme

Das abstrakte Federpendel ist ein relativ einfaches Beispiel, weil die Ableitung y

und die Funktion y lediglich durch eine Matrix, also eine lineare Abbildung, gekoppelt sind und sich deshalb die L¨ osung analytisch angeben l¨ asst.

Die in der Praxis auftretenden Probleme sind in der Regel nicht so einfach zu behandeln. Ein Beispiel ist das Mehrk¨ orperproblem, bei dem n Massen m

, . . . , m

zu einem Zeitpunkt t an n verschiedenen Positionen x

(t), . . . , x

(t) im zwei- oder h¨ oherdimensionalen Raum liegen und mittels der Gravitation aufeinander einwirken.

In diesem Fall ¨ ubt die Masse m

auf die Masse m

eine Kraft von F

(t) = %m

m

x

(t) − x

(t) kx

(t) − x

(t)k

aus, wobei % die Gravitationskonstante ist. Insgesamt wirkt also eine Kraft von F

(t) = X

j=1 j6=i

F

(t) = %m

X

j=1 j6=i

m

x

(t) − x

(t) kx

(t) − x

(t)k

auf die Masse m

, und entsprechend der Newton-Axiome entsteht dadurch eine Beschleu- nigung von

a

(t) = %

n

X

j=1 j6=i

m

x

(t) − x

(t)

kx

(t) − x

(t)k

. (1.4) Wie im Falle des Federpendels benutzen wir die Newtonschen Axiome um festzustellen, dass a

die Ableitung der Geschwindigkeit v

und v

die Ableitung des Ortes x

ist, also

x

(t) = v

(t), v

(t) = a

(t) = %

n

X

j=1 j6=i

m

x

(t) − x

(t)

kx

(t) − x

(t)k

f¨ ur alle t ∈ R

1.3 Schwingende Saite gilt. Wir fassen die Orte x

, die Geschwindigkeiten v

und die Beschleuigungen a

zu Vektoren

x(t) :=





 x

(t)

.. . x

(t)





 , v(t) :=





 v

(t)

.. . v

(t)





 , a(t) :=





 a

(t)

.. . a

(t)







zusammen und schreiben die Differentialgleichung in der Form x

(t)

v

(t)

= v(t)

a(t)

f¨ ur alle t ∈ R

. Gem¨ aß (1.4) k¨ onnen wir a(t) als Funktion A von x schreiben, erhalten also

a(t) = A(x(t)) f¨ ur alle t ∈ R

.

Zur Vereinheitlichung der Darstellung fassen wir x(t) und v(t) zu einem Vektor y(t) :=

x(t) v(t)

f¨ ur alle t ∈ R

zusammen und f¨ uhren die Funktion f (y) :=

y

A(y

)

ein, um die kompakte Darstellung

y

(t) = f (y(t)) f¨ ur alle t ∈ R

zu erhalten.

Im Falle des Federpendels war f lediglich eine lineare Abbildung, im Falle des Mehrk¨ orperproblems ist f nicht linear, und die einfachen analytischen L¨ osungsans¨ atze f¨ ur lineare Probleme lassen sich nicht mehr verwenden.

F¨ ur n ≤ 3 ist es noch m¨ oglich, die L¨ osung y wenigstens formal (etwa durch spezielle Reihenentwicklungen) darzustellen, f¨ ur n > 3 dagegen sind numerische Approximations- verfahren das Mittel der Wahl.

1.3 Schwingende Saite

Wir k¨ onnen Mehrk¨ orpersysteme auch verwenden, um N¨ aherungen im Wesentlichen kon- tinuierlicher Ph¨ anomene zu gewinnen. Als Beispiel untersuchen wir eine Saite, die zwi- schen zwei Punkten eingespannt ist und die zwischen diesen Punkten frei schwingen kann.

Als Modell verwenden wir Punktmassen, die durch Federn aneinander gekoppelt sind.

Sei dazu n ∈ N . Wir verwenden n + 2 Punktmassen, die wir mit {0, . . . , n + 1} durchnu-

merieren. Dabei sollen die Punktmassen mit den Nummern 0 und n + 1 den Endpunkten

der Saite entsprechen. Die Punktmasse mit der Nummer i ∈ {1, . . . , n} soll mit den Punktmassen mit den Nummern i − 1 und i + 1 durch Federn verbunden sein.

Wenn wir die L¨ ange der Saite im Ruhezustand mit ` ∈ R

bezeichnen und wir sie in gleich lange St¨ ucke einteilen, hat jede Feder die Ruhel¨ ange

h := ` n + 1 .

Um die von den Federn ausge¨ ubten Kr¨ afte in Beschleunigungen umrechnen zu k¨ onnen, m¨ ussen wir den Punktmassen Massen zuordnen. Dazu nehmen wir an, dass die Masse einer Saite der Ruhel¨ ange 1 durch ˆ m gegeben ist. F¨ ur i ∈ {1, . . . , n} ersetzt die i-te Punktmasse ein St¨ uck der Saite der L¨ ange h, f¨ ur die Endpunkte i ∈ {0, n + 1} dagegen nur ein Restst¨ uck der L¨ ange h/2. Wenn wir davon ausgehen, dass die Masse gleichm¨ aßig

¨ uber die Saite verteilt ist, ergibt sich f¨ ur die i-te Punktmasse die Masse m

:=

( mh/2 ˆ falls i ∈ {0, n + 1}, ˆ

mh ansonsten.

Da wir uns f¨ ur Schwingungen der Saite interessieren, bietet es sich an, die Positionen der einzelnen Punktmassen zu untersuchen. Wir bezeichnen mit x

(t) die Position der i-ten Masse zu einem Zeitpunkt t ∈ R.

Die von der Feder zwischen der i-ten und der (i + 1)-ten Punktmasse auf erstere ausge¨ ubte Kraft ist nach dem Hooke’schen Gesetz (1.1) durch

ˆ c

h (x

(t) − x

(t))

gegeben. Die von der Feder zwischen der (i −1)-ten und der i-ten Punktmasse ausge¨ ubte Kraft addiert sich hinzu, so dass wir insgesamt die Kraft

F

(t) = ˆ c

h (x

(t) − x

(t)) ˆ c

h (x

(t) − x

(t))

= ˆ c

h (x

(t) − 2x

(t) + x

(t)),

erhalten. Die Beschleunigung ergibt sich wie zuvor, indem wir durch die Masse m

divi- dieren:

a

(t) = 1 m

ˆ c

h (x

(t) − 2x

(t) + x

(t))

= c ˆ ˆ m

x

(t) − 2x

(t) + x

(t)

h

.

Anders als im Fall der Gravitationskraft sind die einzelnen Komponenten der Vektoren x

(t), F

(t) und a

(t) voneinander v¨ ollig unabh¨ angig, so dass wir sie getrennt voneinander untersuchen k¨ onnen.

Im Folgenden beschr¨ anken wir uns auf eine der Komponenten, verwenden aber weiter- hin die bisher eingesetzte Notation. Mit den Newton-Axiomen erhalten wir das System

x

(t) = v

(t), v

(t) = c ˆ ˆ m

x

(t) − 2x

(t) + x

(t)

h

f¨ ur alle t ∈ R . (1.5)

1.3 Schwingende Saite

Wie zuvor k¨ onnen wir das System kompakt schreiben, indem wir die Gr¨ oßen zu Vektoren

x(t) :=





 x

(t)

.. . x

(t)





 , v(t) :=





 v

(t)

.. . v

(t)







zusammenfassen. Die Berechnung der Beschleunigung aus x(t) l¨ asst sich durch die Matrix

L := c ˆ ˆ mh













2 −1

−1 . .. ...

. .. ... −1

−1 2













kompakt in die Form

v

(t) = −Lx(t) bringen, so dass sich insgesamt

x

(t) v

(t)

=

I

−L

x(t) v(t)

ergibt. Man erkennt eine gewisse ¨ Ahnlichkeit zu der Gleichung (1.2), durch die wir das Federpendel beschrieben haben: Die Matrix L tritt an die Stelle des Parameters λ.

Die durch Federn verbundenen Punktmassen sind lediglich als N¨ aherung einer kontinu- ierlichen Saite gedacht. Um zu einer entsprechenden Gleichung zu gelangen, untersuchen wir den Quotienten in der Gleichung (1.5) etwas n¨ aher.

Lemma 1.1 (Differenzenquotient f¨ ur die zweite Ableitung) Sei h ∈ R

, und sei g : [−h, h] → R viermal stetig differenzierbar. Dann existiert ein η ∈ [−h, h] mit

g(h) − 2g(h) + g(−h)

h

= g

(0) + h

12 g

(η).

Beweis. Mit dem Satz von Taylor erhalten wir g(h) = g(0) + hg

(0) + h

2 g

(0) + h

6 g

(0) + h

24 g

(η

), g(−h) = g(0) − hg

(0) + h

2 g

(0) − h

6 g

(0) + h

24 g

(η

)

f¨ ur geeignete Zwischenpunkte η

∈ [0, h] und η

∈ [−h, 0]. Indem wir beide Gleichungen addieren ergibt sich

g(h) + g(−h) = 2g(0) + h

g

(0) + h

12

g

(η

) + g

(η

)

2 .

Wir bringen 2g(0) auf die linke Seite und dividieren durch h

, um zu g(h) − 2g(0) + g(−h)

h

= g

(0) + h

12

g

(η

) + g

(η

) 2

zu gelangen. Da g

stetig ist, finden wir mit dem Zwischenwertsatz ein η ∈ [η

, η

] ⊆ [−h, h], das

g

(η

) + g

(η

)

2 = g

(η) erf¨ ullt. Damit ist unsere Gleichung bewiesen.

Um Lemma 1.1 auf die Gleichung (1.5) anwenden zu k¨ onnen, legen wir f¨ ur jeden Zeitpunkt t ∈ R eine viermal stetig differenzierbare Kurve s 7→ x(t, s) derart durch die Punkte x

(t), das

x(t, s) = x

(t) f¨ ur alle i ∈ {0, . . . , n + 1}, t ∈ R , s = ih gilt. Dann folgt mit dem Lemma

x

(t) − 2x

(t) + x

(t)

h

= x(t, ih + h) − 2x(t, ih) + x(t, ih − h) h

= ∂

x

∂s

(t, s) + h

12

∂

x

∂s

(t, η

)

mit einem Parameter η

. Wenn wir h gegen null streben lassen und annehmen, dass die vierten Ableitungen beschr¨ ankt bleiben, wird so aus der Gleichung (1.5) die partielle Differentialgleichung

∂x

∂t (t, s) = v(t, s), ∂v

∂t (t, s) = c ˆ ˆ m

∂

x

∂s

(t, s) f¨ ur alle t ∈ R, s ∈ (0, `), aus der sich durch Elimination der Geschwindigkeit v(t, s) die eindimensionale Wellen- gleichung

∂

x

∂t

(t, s) = ˆ c ˆ m

∂

x

∂s

(t, s) f¨ ur alle t ∈ R , s ∈ (0, `) ergibt. Sie beschreibt die Schwingung einer kontinuierlichen Saite.

1.4 W¨ armeleitung

Ein weiteres Beispiel f¨ ur eine partielle Differentialgleichung ist die eindimensionale W¨ armeleitungsgleichung

∂u

∂t (x, t) = κ ∂

u

∂x

(x, t) f¨ ur alle t ∈ R

, x ∈ [0, 1]. (1.6)

1.4 W¨ armeleitung Sie beschreibt die Erw¨ armung oder Abk¨ uhlung eines Drahtes der L¨ ange 1: x ∈ [0, 1] gibt die Position auf dem Draht an, t den Zeitpunkt, und u(x, t) ist die Temperatur im Punkt x zum Zeitpunkt t. Wir nehmen zur Vereinfachung an, dass die Randbedingungen

u(0, t) = u(1, t) = 0 f¨ ur alle t ∈ R

gelten, dass also die Temperatur an den beiden Endpunkten des Drahts fixiert ist.

Um diese Gleichung numerisch behandeln zu k¨ onnen, kehren wir den Weg um, den wir bei der Wellengleichung gegangen sind: Mit Lemma 1.1 gilt

u(x − h, t) − 2u(x, t) + u(x + h, t)

h

= ∂

u

∂x

(x, t) + h

12

∂

u

∂x

(η, t) (1.7) f¨ ur ein geeignetes η ∈ [0, 1]. Falls die vierte Ableitung gleichm¨ aßig beschr¨ ankt ist, k¨ onnen wir also die zweite Ableitung durch den Differenzenquotienten auf der linken Seite appro- ximieren, und die Approximation wird wie h

gegen die korrekte Ableitung konvergieren.

F¨ ur ein n ∈ N w¨ ahlen wir eine Schrittweite h := 1

n + 1

und ersetzen das kontinuierliche Intervall [0, 1] durch (n + 2) ∈ N diskrete Punkte 0 = x

< x

< . . . < x

< x

= 1, die durch

x

:= ih f¨ ur alle i ∈ {0, . . . , n + 1}

gegeben sind, und wir beschreiben entsprechend die Funktion u(x, t) durch den Vektor y(t) = (y

(t))

mit

y

(t) = u(x

, t) f¨ ur alle t ∈ R

, i ∈ {1, . . . , n}.

Indem wir die zweite Ableitung durch (1.7) approximieren, stellen wir fest, dass die W¨ armeleitungsgleichung (1.6) durch die Gleichungen

y

(t) ≈



 

 

κ

h²

(y

(t) − 2y

(t) + y

(t)) falls 1 < i < n,

κ

h²

(−2y

(t) + y

(t)) falls i = 1,

κ

h²

(y

(t) − 2y

(t)) falls i = n,

f¨ ur alle i ∈ {1, . . . , n}

und t ∈ R

approximiert wird. Wenn wir also das System

y

(t) = Ay(t) f¨ ur alle t ∈ R

mit der Matrix A ∈ R

, gegeben durch

A := − κ h















2 −1

−1 2 −1 . .. ... ...

−1 2 −1

−1 2















,

l¨ osen, d¨ urfen wir darauf hoffen, dass die Komponenten y

(t) des Vektors y(t) gute

N¨ aherungen f¨ ur die Werte u(x

, t) der tats¨ achlichen L¨ osung sind. Es handelt sich also

wieder um ein lineares Anfangswertproblem, allerdings steht diesmal die Approximati-

on eines Differentialoperators hinter der Matrix A. Eine genauere Analyse zeigt, dass

sich bestimmte Eigenschaften dieses Operators auf die Matrix A ¨ ubertragen und dazu

f¨ uhren, dass sich einfache Verfahren zur Behandlung von Anfangswertproblemen f¨ ur die-

ses Problem nicht gut eignen. Ein wichtiges Ziel wird deshalb darin bestehen, Techniken

zu entwickeln, mit denen sich auch dieses Problem effizient behandeln l¨ asst.

2 Einige theoretische Aussagen ¨ uber gew¨ ohnliche Differentialgleichungen

Bevor wir numerische L¨ osungsverfahren f¨ ur gew¨ ohnliche Differentialgleichungen untersu- chen k¨ onnen, m¨ ussen wir zun¨ achst kl¨ aren, unter welchen Bedingungen diese Gleichungen uberhaupt eine L¨ ¨ osung besitzen. In Hinblick auf die numerische Behandlung ist eben- falls wichtig, wie empfindlich die L¨ osung auf St¨ orungen der Parameter, insbesondere des Startwerts, reagiert.

2.1 Allgemeine Problemstellung

Wir konzentrieren uns auf die Analyse des Anfangswertproblems

y(a) = y

, y

(t) = f (t, y(t)) f¨ ur alle t ∈ [a, b] (2.1) auf einem kompakten Intervall [a, b] mit einem Startwert y

in einem Banachraum V und einer Funktion f : [a, b] × V → V . Gesucht ist eine mindestens einmal stetig diffe- renzierbare Funktion y : [a, b] → V .

Das allgemeinere Problem

y(a) = y

, y

(a) = y

, . . . , y

(a) = y

, y

(t) = f(t, y(t), y

(t), . . . , y

(t)) f¨ ur alle t ∈ [a, b]

l¨ asst sich auf die Form (2.1) zur¨ uckf¨ uhren, indem wir den Hilfsvektor

w(t) :=









 y(t) y

(t) .. . y

(t)











f¨ ur alle t ∈ [a, b]

einf¨ uhren und das erweiterte System

w(a) =









 y

y

.. . y











, w

(t) =















w

(t) w

(t)

.. . w

(t)

f (t, w

(t), w

(t), . . . , w

(t))















f¨ ur alle t ∈ [a, b] l¨ osen.

2.2 Existenz und Eindeutigkeit

Bei der Untersuchung der Eigenschaften eines Anfangswertproblems hat es sich als sehr n¨ utzlich erwiesen, anstelle der differentiellen Formulierung (2.1) eine Integralformulie- rung zu verwenden, die ohne die Forderung nach Differenzierbarkeit auskommt.

Lemma 2.1 (Integralformulierung) Sei eine stetige Funktion f ∈ C([a, b] × V, V ) gegeben. Falls eine Funktion y ∈ C

([a, b], V ) das Anfangswertproblem (2.1) l¨ ost, gilt

y(t) = y

+ Z

a

f (s, y(s)) ds f¨ ur alle t ∈ [a, b]. (2.2) Falls umgekehrt eine stetige Funktion y ∈ C([a, b], V ) die Integralgleichung (2.2) erf¨ ullt, ist sie auch stetig differenzierbar und l¨ ost das Anfangswertproblem (2.1).

Beweis. Im ersten Schritt gehen wir davon aus, dass y das Anfangswertproblem l¨ ost.

Nach Hauptsatz der Differential- und Integralrechnung gilt dann y

+

Z

f (s, y(s)) ds = y

+ Z

a

y

(s) ds = y(a) + y(t) − y(a) = y(t), f¨ ur alle t ∈ [a, b], also die Integralgleichung (2.2).

Im zweiten Schritt gehen wir davon aus, dass y ∈ C([a, b], V ) die Gleichung (2.2) erf¨ ullt. F¨ ur t = a folgt aus ihr unmittelbar y(a) = y

. Wir m¨ ussen nachweisen, dass y differenzierbar ist und seine Ableitung die gew¨ unschte Eigenschaft aufweist.

Sei t ∈ [a, b). Die Ableitung y

(t) ist als Grenzwert des Quotienten y(t + h) − y(t)

h = 1

h

Z

f(s, y(s)) ds − Z

a

f (s, y(s)) ds

= 1 h

Z

f (s, y(s)) ds f¨ ur h → 0 definiert. Wir m¨ ussen zeigen, dass dieser Grenzwert existiert und mit f(t, y(t))

¨ ubereinstimmt, also sollte

f (t, y(t)) − y(t + h) − y(t) h

= 1 h

Z

f (t, y(t)) − f (s, y(s)) ds

≤ 1 h

Z

kf (t, y(t)) − f (s, y(s))k ds f¨ ur h → 0 gegen null konvergieren.

Mit dem Mittelwertsatz der Integralrechnung finden wir einen Zwischenpunkt η ∈ [t, t + h] derart, dass

Z

kf (t, y(t)) − f (s, y(s))k ds = hkf (t, y(t)) − f (η, y(η))k gilt, also folgt

f (t, y(t)) − y(t + h) − y(t) h

≤ kf (t, y(t)) − f (η, y(η))k.

2.2 Existenz und Eindeutigkeit Wegen η ∈ [t, t + h] impliziert h → 0 auch η → t, und aufgrund der Stetigkeit der Funktionen f und y d¨ urfen wir auf f (η, y(η)) → f (t, y(t)) schließen, so dass wir insgesamt

h→0

lim

f (t, y(t)) − y(t + h) − y(t) h

= 0 bewiesen haben, also y

(t) = f (t, y(t)).

F¨ ur t = b folgt die Aussage, indem wir entsprechend den linksseitigen Differenzenquo- tienten zur Approximation der Ableitung einsetzen.

Bei genauerer Betrachtung stellt sich heraus, dass die Integralgleichung (2.2) eine Fixpunktgleichung ist: Wir definieren auf dem Raum U := C([a, b], V ) den Operator Ψ : U → U durch

Ψ[y](t) = y

+ Z

a

f (s, y(s)) ds f¨ ur alle y ∈ U, t ∈ [a, b]

und halten fest, dass (2.2) sich dann kurz als

y = Ψ[y] (2.3)

schreiben l¨ asst. Nach Lemma 2.1 ist also das L¨ osen des Anfangswertproblems ¨ aquivalent dazu, einen Fixpunkt des Operators Ψ zu finden.

Ein zentrales Hilfsmittel f¨ ur den Beweis von Existenz und Eindeutigkeit von Fixpunk- ten ist der folgende Fixpunktsatz von Banach.

Satz 2.2 (Banach) Sei X eine vollst¨ andige Teilmenge eines normierten Raumes. Sei Ψ : X → X eine Abbildung, und sei L ∈ [0, 1) eine Zahl mit

kΨ(u) − Ψ(v)k ≤ Lku − vk f¨ ur alle u, v ∈ X. (2.4) Dann besitzt Ψ einen Fixpunkt in X, es existiert also ein u

∈ X mit

Ψ(u

) = u

. Dieser Fixpunkt ist eindeutig bestimmt.

Beweis. (vgl. [1, Th´ eor` eme 6] und [4]) Sei u

∈ X. Wir definieren die Folge (u

)

durch

u

= Ψ(u

) f¨ ur alle m ∈ N

. Unser Ziel ist es, nachzuweisen, dass (u

)

eine Cauchy-Folge ist.

Zun¨ achst beweisen wir

ku

− u

k ≤ L

ku

− u

k (2.5)

durch Induktion f¨ ur alle m ∈ N

. F¨ ur m = 0 ist (2.5) trivial.

Gelte nun (2.5) f¨ ur ein m ∈ N

. Nach Voraussetzung gilt dann

ku

− u

k = kΨ(u

) − Ψ(u

)k

≤ Lku

− u

k

≤ LL

ku

− u

k = L

ku

− u

k, und der Induktionsschritt ist bewiesen.

Seien nun m ∈ N

und n ∈ N

gegeben. Dann gilt ku

− u

k =

n−m−1

X

j=0

u

− u

≤

m−n−1

X

j=0

ku

− u

k

(2.5)

≤

m−n−1

X

j=0

L

ku

− u

k = ku

− u

kL

m−n−1

X

j=0

L

≤ ku

− u

kL

∞

X

j=0

L

= ku

− u

k L

1 − L

dank der geometrischen Summenformel. Mit dieser Absch¨ atzung k¨ onnen wir nachweisen, dass (u

)

eine Cauchy-Folge ist: Sei ∈ R

. Wir w¨ ahlen m

∈ N

so, dass

ku

− u

k L

1 − L ≤ gilt. F¨ ur alle m, n ∈ N

mit m

≤ m ≤ n gilt dann

ku

− u

k ≤ ku

− u

k L

1 − L ≤ ku

− u

k L

1 − L ≤ , also ist (u

)

eine Cauchy-Folge.

Da X vollst¨ andig ist, muss es ein u

∈ X mit

m→∞

lim ku

− u

k = 0

geben, und wir m¨ ussen nur noch nachpr¨ ufen, dass u

auch ein Fixpunkt von Ψ ist.

Sei dazu ∈ R

. Da (u

)

gegen u

konvergiert, gibt es ein m ∈ N

so, dass ku

− u

k ≤ /2, ku

− u

k ≤ /2

gelten, und wir erhalten wegen u

= Ψ(u

) die Absch¨ atzung ku

− Ψ(u

)k = ku

− u

+ Ψ(u

) − Ψ(u

)k

≤ ku

− u

k + kΨ(u

) − Ψ(u

)k

≤ ku

− u

k + Lku

− u

k

≤ /2 + L/2 < .

2.2 Existenz und Eindeutigkeit Da beliebig gew¨ ahlt werden kann, folgt u

= Ψ(u

), also ist u

in der Tat ein Fixpunkt.

Zum Nachweis der Eindeutigkeit w¨ ahlen wir einen zweiten Fixpunkt u

∈ X und erhalten

ku

− u

k = kΨ(u

) − Ψ(u

)k ≤ Lku

− u

k, also folgt aus L < 1 bereits u

= u

.

Indem wir diesen Satz auf die alternative Formulierung (2.3) anwenden, erhalten wir die folgende f¨ ur uns zentrale Aussage ¨ uber Existenz und Eindeutigkeit der L¨ osung eines Anfangswertproblems.

Satz 2.3 (Picard-Lindel¨ of ) Die Funktion f ∈ C([a, b] × V, V ) erf¨ ulle die globale Lipschitz-Bedingung

kf(t, x) − f (t, y)k ≤ L

kx − yk f¨ ur alle t ∈ [a, b] und x, y ∈ V. (2.6) Dann besitzt das Anfangswertproblem (2.1) eine eindeutige L¨ osung y ∈ C

([a, b], V ).

Beweis. (vgl. [10, Abschnitt 1.6]) Wir f¨ uhren den Beweis mit Hilfe des Banachschen Fixpunktsatzes 2.2: Dazu f¨ uhren wir den Operator Ψ durch

Ψ[u](t) := y

+ Z

a

f (s, u(s)) ds f¨ ur alle t ∈ [a, b], u ∈ C([a, b], V )

ein, der den Banachraum X := C([a, b], V ) in sich abbildet, und untersuchen die von ihm induzierte Fixpunktiteration. Nach Lemma 2.1 wissen wir n¨ amlich, dass ein Fixpunkt des Operators Ψ gerade eine L¨ osung des Anfangswertproblems (2.1) ist.

Damit wir Satz 2.2 anwenden k¨ onnen, m¨ ussen wir eine geeignete Norm auf dem Raum C([a, b], V ) einf¨ uhren. Wir verwenden die gewichtete Supremumsnorm

kuk

:= sup

e

ku(x)k : x ∈ [a, b] , f¨ ur alle u ∈ C([a, b], V ),

die wegen 0 < e

≤ e

¨ aquivalent zu der ¨ ublichen Supremumsnorm ist, so dass C([a, b], V ) auch mit dieser Norm vollst¨ andig ist. Bez¨ uglich dieser Norm gilt

e

kΨ[u](t) − Ψ[v](t)k = e

Z

f (s, u(s)) − f(s, v(s)) ds

≤ e

Z

a

kf(s, u(s)) − f(s, v(s))k ds

≤ L

e

Z

a

ku(s) − v(s)k ds

= L

e

Z

a

e

e

ku(s) − v(s)k ds

≤ L

e

Z

a

e

ku − vk

ds

= 1

2 e

ku − vk

Z

a

2L

e

ds

= 1

2 e

ku − vk

(e

− e

) ≤ 1

2 ku − vk

f¨ ur alle t ∈ [a, b] und alle u, v ∈ C([a, b], V ), wobei wir im vorletzten Schritt ausgenutzt haben, dass s 7→ 2L

e

die Ableitung der Funktion s 7→ e

ist, so dass sich das Integral mit dem Hauptsatz berechnen l¨ asst. Indem wir zu dem Maximum ¨ uber alle t ∈ [a, b] ¨ ubergehen folgt

kΨ[u] − Ψ[v]k

≤ 1

2 ku − vk

f¨ ur alle u, v ∈ C([a, b], V ),

so dass wir Satz 2.2 anwenden k¨ onnen, um zu folgern, dass ein eindeutig bestimmter Fixpunkt y ∈ C([a, b], V ) mit Ψ[y] = y existiert.

Nach Lemma 2.1 ist diese Funktion y auch die eindeutig bestimmte L¨ osung des An- fangswertproblems.

Satz 2.3 ist nicht nur ein Existenz- und Eindeutigkeitsresultat, er bietet uns auch ein Konstruktionsverfahren f¨ ur die L¨ osung des Anfangswertproblems:

Bemerkung 2.4 (Picard-Iteration) Ausgehend von einer beliebigen Funktion u

k¨ onnen wir, wie im Satz 2.2, die Folge u

:= Ψ(u

) konstruieren, und Satz 2.3 im- pliziert, dass diese Folge gegen die L¨ osung des Anfangswertproblems (2.1) konvergieren wird. Diese Konstruktion tr¨ agt den Namen Picard-Iteration.

F¨ ur die Praxis ist diese Konstruktion nur dann anwendbar, wenn sich die einzelnen Iterierten u

geeignet im Rechner darstellen lassen, etwa mit Hilfe einer Diskretisierung.

2.3 St¨ orungen der Daten

F¨ ur die numerische Behandlung des Anfangswertproblems (2.1) ist neben der prinzipiel- len L¨ osbarkeit auch der Einfluss von St¨ orungen relevant, schließlich wird im praktischen Algorithmus in der Regel mit Gleitpunktarithmetik beschr¨ ankter Genauigkeit gearbeitet.

Ein wichtiges Hilfsmittel f¨ ur die Analyse ist die Gr¨ onwallsche Ungleichung, von der wir hier nur die folgende vereinfachte Variante ben¨ otigen:

Lemma 2.5 (Gr¨ onwall) Seien [a, b] ⊆ R ein Intervall, sei α ∈ C[a, b] eine monoton wachsende Funktion, sei β ∈ R

. Falls eine Funktion u ∈ C[a, b] die Absch¨ atzung

u(t) ≤ α(t) + β Z

a

u(s) ds f¨ ur alle t ∈ [a, b] (2.7) erf¨ ullt, gilt die Ungleichung

u(t) ≤ α(t)e

f¨ ur alle t ∈ [a, b].

2.3 St¨ orungen der Daten Beweis. (vgl. [2] und [5]) Sei u ∈ C[a, b] eine Funktion, die (2.7) erf¨ ullt.

Wir f¨ uhren die Hilfsfunktion v ∈ C([a, b)) mit v(t) := e

Z

βu(s) ds f¨ ur alle t ∈ [a, b]

ein und erhalten mit der Produktregel und (2.7) die Absch¨ atzung v

(t) = −βe

Z

βu(s) ds + e

βu(t)

= βe

u(t) − Z

a

βu(s) ds

≤ βe

α(t) f¨ ur alle t ∈ [a, b]. Aus v(a) = 0 folgt

e

Z

a

βu(s) ds = v(t) = v(t) − v(a) = Z

a

v

(s) ds ≤ β Z

a

α(s)e

ds

≤ βα(t) Z

a

e

ds

= βα(t)

− 1 β

(e

− e

) = α(t) − α(t)e

und indem wir mit e

multiplizieren

β Z

a

u(s) ds ≤ α(t)e

− α(t).

Durch Einsetzen in (2.7) gelangen wir zu u(t) ≤ α(t) +

Z

βu(s) ds ≤ α(t) + α(t)e

− α(t) = α(t)e

, und das ist die zu beweisende Ungleichung.

Mit Hilfe dieses Korollars und des Lemmas 2.1 k¨ onnen wir nun den Einfluss von St¨ orungen der Anfangsdaten untersuchen:

Satz 2.6 (St¨ orungen) Sei U ⊆ V . Die Funktion f ∈ C([a, b] × U, U) erf¨ ulle die bereits aus Satz 2.3 bekannte globale Lipschitz-Bedingung (2.6). Sei g ∈ C([a, b] × U, U ) eine weitere Funktion.

Seien y

, z

∈ U , und seien y, z ∈ C

([a, b], U ) L¨ osungen der Anfangswertprobleme y(a) = y

, y

(t) = f(t, y(t)),

z(a) = z

, z

(t) = g(t, z(t)) f¨ ur alle t ∈ [a, b].

Dann gilt die Absch¨ atzung ky(t) − z(t)k ≤ e

ky

− z

k + Z

a

kf (s, z(s)) − g(s, z(s))k ds

f¨ ur alle t ∈ [a, b],

kleine St¨ orungen der Anfangsdaten und der rechten Seite f¨ uhren also auch nur zu kleinen

St¨ orungen der L¨ osung.

Beweis. Mit Lemma 2.1 erhalten wir y(t) − z(t) = y

− z

+

Z

f (s, y(s)) − g(s, z(s)) ds, ky(t) − z(t)k ≤ ky

− z

k +

Z

kf (s, y(s)) − g(s, z(s))k ds

≤ ky

− z

k + Z

a

kf (s, y(s)) − f (s, z(s))k + kf (s, z(s)) − g(s, z(s))k ds

(2.6)

≤ ky

− z

k + Z

a

L

ky(s) − z(s)k ds + Z

a

kf (s, z(s)) − g(s, z(s))k ds. (2.8) Wir definieren

β := L

, α(t) := ky

− z

k + Z

a

kf(s, z(s)) − g(s, z(s))k ds, u(t) := ky(t) − z(t)k und stellen fest, dass (2.8) gerade

u(t) ≤ α(t) + β Z

a

u(s) ds f¨ ur alle t ∈ [a, b]

entspricht. Da β nicht-negativ und α monoton wachsen ist, k¨ onnen wir die Gr¨ onwall- Ungleichung aus Lemma 2.5 anwenden und erhalten

ky(t) − z(t)k = u(t) ≤ α(t)e

= e

ky

− z

k + Z

a

kf (s, z(s)) − g(s, z(s))k ds

, also die gew¨ unschte Absch¨ atzung.

Beispiel 2.7 (Entfernung vom Anfangswert) Neben der offensichtlichen Anwen- dung auf gest¨ orte Anfangswerte und rechte Seiten l¨ asst sich Satz 2.6 auch anderweitig verwenden.

Beispielsweise k¨ onnen wir g = 0 und z

= y

einsetzen. Dann gilt offenbar z(t) = y

und wir erhalten

ky(t) − y

k ≤ e

Z

a

kf(s, y

)k ds,

k¨ onnen also absch¨ atzen, wie schnell sich die L¨ osung des Anfangswertproblems vom An-

fangswert entfernt.

3 Einschrittverfahren

Die exakte L¨ osung eines Anfangswertproblems der Form (2.1) wird sich im allgemei- nen Fall nicht exakt berechnen lassen. Stattdessen m¨ ussen wir auf eine Approximati- on zur¨ uckgreifen: Statt nach einer geschlossenen Formel f¨ ur die L¨ osung zu suchen, be- schr¨ anken wir uns darauf, sie nur in einzelnen Punkten t

, . . . , t

∈ [a, b] n¨ aherungsweise zu berechnen.

Wir sind nat¨ urlich an Verfahren interessiert, die uns eine m¨ oglichst genaue N¨ aherung zur Verf¨ ugung stellen, und das bei m¨ oglichst geringem Rechen- und Speicheraufwand.

Ein m¨ oglicher Zugang w¨ are etwa die Picard-Iteration (vgl. Bemerkung 2.4): Wir k¨ onnten die Iterierten durch ihre Werte in Punkten des Intervalls approximieren und zur Berechnung der Integrale eine Quadraturformel verwenden, die nur diese Punktwer- te ben¨ otigt. Der Nachteil dieses Zugangs besteht darin, dass alle Punktwerte gleichzeitig im Speicher gehalten werden m¨ ussen.

Wir suchen stattdessen nach einem Verfahren, bei dem wir die Werte zu den verschie- denen Zeitpunkten der Reihe nach berechnen k¨ onnen. Ein Einschrittverfahren versucht, den Wert zu einem Zeitpunkt t

nur auf Grundlage des Wertes zum unmittelbar vor- angehenden Zeipunkt t

zu approximieren.

Um die Diskussion der L¨ osbarkeit zu vermeiden, setzen wir, sofern nicht gesondert erw¨ ahnt, im folgenden Kapitel voraus, dass die rechte Seite f des Anfangswertproblems (2.1) im zweiten Argument Lipschitz-stetig ist (vgl. Bedingung (2.6)). Satz 2.3 impliziert dann die eindeutige L¨ osbarkeit f¨ ur beliebige Startwerte in V und Startpunkte in [a, b].

3.1 Euler-Verfahren

Wir untersuchen zun¨ achst ein besonders einfaches Einschrittverfahren: Das Euler- Verfahren l¨ asst sich aus der in Lemma 2.1 eingef¨ uhrten Integralformulierung gewinnen.

Wir beschr¨ anken uns f¨ ur den Moment auf den Fall t = b, f¨ ur den y(b) = y

+

Z

f (s, y(s)) ds

gilt. Es bietet sich an, das Integral mit einer Quadraturformel mit Quadraturgewichten w

, . . . , w

∈ R und Quadraturpunkten s

, . . . , s

∈ [a, b] zu approximieren:

y(b) ≈ y

+

m

X

j=0

w

f(s

, y(s

)).

Leider k¨ onnen wir die rechte Seite dieser Gleichung im Allgemeinen nicht auswerten, da

uns die Werte y(s

), . . . , y(s

) nicht zur Verf¨ ugung stehen.

Allerdings kennen wir y(a) = y

, so dass wir immerhin eine Quadraturformel mit m = 0 und s

= a verwenden k¨ onnten. Damit wenigstens konstante Funktionen von dieser Quadraturformel exakt integriert werden, m¨ ussen wir das Gewicht w

= b − a verwenden und erhalten

y(b) ≈ y(b) := ˜ y

+ (b − a)f (a, y(a)). (3.1) Falls y zweimal stetig differenzierbar ist, k¨ onnen wir den Fehler wie folgt absch¨ atzen:

Lemma 3.1 (Genauigkeit) Sei y ∈ C

([a, b], V ), und sei y(b) ˜ wie in (3.1) definiert.

Dann existiert ein η ∈ [a, b] mit

ky(b) − y(b)k ≤ ˜ (b − a)

2 ky

(η)k.

Beweis. Um die Aufgabe etwas besser zug¨ anglich zu machen, f¨ uhren wir sie auf das Einheitsintervall [0, 1] zur¨ uck, indem wir die Funktion

ˆ

y : [0, 1] → V, t 7→ y(a + (b − a)t),

untersuchen. Sie erf¨ ullt offenbar y(a) = ˆ y(0), y(b) = ˆ y(1) sowie nach Definition des Anfangswertproblems und Kettenregel

˜

y(b) = y

+ (b − a)f (a, y(a)) = y(a) + (b − a)y

(a) = ˆ y(0) + ˆ y

(0).

Mit dem Hauptsatz der Integral- und Differentialrechnung erhalten wir y(b) − y(b) = ˆ ˜ y(1) − y(0) ˆ − y ˆ

(0) =

Z

0

ˆ

y

(t) dt − y ˆ

(0) = Z

0

ˆ

y

(t) − y ˆ

(0) dt.

Da mit y auch ˆ y zweimal stetig differenzierbar ist, k¨ onnen wir den Hauptsatz erneut anwenden, um zu

y(b) − y(b) = ˜ Z

0

Z

ˆ

y

(s) ds dt

zu gelangen. Um die Integrationsgrenzen des inneren Integrals von t unabh¨ angig zu machen, substituieren wir s = tr und erhalten

y(b) − y(b) = ˜ Z

0

t Z

0

ˆ

y

(tr) dr dt.

Wir wollen den Beweis mit dem Mittelwertsatz der Integralrechnung abschließen, der nur f¨ ur reellwertige Funktionen gilt. Also gehen wir zu der Norm ¨ uber und erhalten

ky(b) − y(b)k ≤ ˜ Z

0

t Z

0

kˆ y

(tr)k dr dt.

Nun k¨ onnen wir den Mittelwertsatz der Integralrechnung erst auf das ¨ außere und dann auf das innere Integral anwenden, um η

, η

∈ [0, 1] mit

ky(b) − y(b)k ˜ = Z

0

t dt Z

0

kˆ y

(η

r)k dr = 1

2 k y ˆ

(η

η

)k

3.1 Euler-Verfahren

zu finden. Per Kettenregel folgt ky(b) − y(b)k ˜ = k 1

2 y ˆ

(η

η

)k = k (b − a)

2 y

(a + (b − a)η

η

)k, also mit η := a + (b − a)η

η

die Behauptung.

Bemerkung 3.2 (Optimale Absch¨ atzung) Aus Lemma 3.1 folgt die Absch¨ atzung ky(b) − y(b)k ≤ ˜ (b − a)

2 ky

k

.

F¨ ur den Fall V = R k¨ onnen wir den Beweis des Lemmas 3.1 so modifizieren, dass wir y(b) − y(b) = ˜ (b − a)

2 y

(η)

f¨ ur ein η ∈ [a, b] erhalten. Wenn wir die Absch¨ atzung auf y(t) = (t − a)

/2 anwenden, gelten y

(t) = t − a, y

(t) = 1 und y(b) = 0, so dass wir ˜

y(b) − y(b) = ˜ (b − a)

2 − 0 = (b − a)

2 ky

k

erhalten. Wir haben also ein Beispiel gefunden, in dem sich unsere Absch¨ atzung nicht verbessern l¨ asst.

Offenbar ist der Fehler um so kleiner, je k¨ urzer das Intervall ist, auf dem die N¨ aherung verwendet wird. Deshalb zerlegen wir das Intervall [a, b] in Teilintervalle: Wir w¨ ahlen n ∈ N sowie t

, . . . , t

∈ [a, b] mit a = t

< t

< . . . < t

= b. Unser Ziel ist es, N¨ aherungswerte der L¨ osung y in diesen Punkten zu berechnen.

Mit dem Hauptsatz der Integral- und Differentialrechnung finden wir y(t

) = y(t

) +

Z

ti−1

y

(s) ds = y(t

) + Z

ti−1

f (s, y(s)) ds f¨ ur alle i ∈ {1, . . . , n}, und indem wir (3.1) auf die Intervalle [t

, t

] anwenden, folgt

y(t

) ≈ y(t

) + (t

− t

)f (t

, y(t

)) f¨ ur alle i ∈ {1, . . . , n}.

Um die Formel etwas zu verk¨ urzen definieren wir die Schrittweiten h

:= t

− t

f¨ ur alle i ∈ {1, . . . , n}

und schreiben die Gleichung in der Form

y(t

) ≈ y(t

) + h

f (t

, y(t

)) f¨ ur alle i ∈ {1, . . . , n}.

Da wir y(t

), . . . , y(t

) nicht kennen, k¨ onnen wir diese Approximation nicht direkt einsetzen, wir k¨ onnen allerdings der Reihe nach N¨ aherungsl¨ osungen berechnen, die die Stelle der exakten Werte annehmen. Damit erhalten wir die Rechenvorschrift

˜

y(t

) := y

, y(t ˜

) := ˜ y(t

) + h

f (t

, y(t ˜

)) f¨ ur alle i ∈ {1, . . . , n}.

Offenbar ist dieses Verfahren sehr effizient durchf¨ uhrbar: In jedem Schritt muss f einmal ausgewertet und eine Linearkombination berechnet werden, und es brauchen nur jeweils

˜

y(t

) und ˜ y(t

) gleichzeitig im Speicher gehalten zu werden. Da der Wert ˜ y(t

) jeweils direkt berechnet werden kann, spricht man von einem expliziten Verfahren, n¨ amlich von dem expliziten Euler-Verfahren.

Bemerkung 3.3 (Diskretisierung) Das Euler-Verfahren ist das erste Diskretisie- rungsverfahren, das wir behandeln. Der Name stammt daher, dass das kontinuierliche Intervall [a, b] durch die diskrete Punktmenge {t

, . . . , t

} ersetzt wird. Im allgemeinen Fall spricht man schon von einer Diskretisierung, wenn ein unendlich-dimensionaler Funktionenraum durch einen endlich-dimensionalen Raum ersetzt wird. Ein derartiger Schritt ist fast immer erforderlich, wenn Differentialgleichungen mit Hilfe eines Com- puters gel¨ ost werden sollen, da einem Computer nur endlich viel Speicher und seinem Benutzer nur endlich viel Zeit zur Verf¨ ugung steht.

Wir k¨ onnen uns bei der Approximation der Gleichung y(b) = y(a) +

Z

f (s, y(s)) ds

auch auf eine Quadraturformel st¨ utzen, die y(b) statt y(a) verwendet. So erhalten wir y(b) ≈ y(a) + (b − a)f(b, y(b)),

also die N¨ aherung einer Fixpunktgleichung. Dementsprechend k¨ onnen wir eine N¨ aherung

˜

y(b) des Werts y(b) durch

˜

y(b) = y(a) + (b − a)f (b, y(b)) ˜ (3.2) definieren, falls sich dieses, im allgemeinen nichtlineare, Gleichungssystem l¨ osen l¨ asst.

Ein brauchbarer Ansatz hierzu ist eine Fixpunkt-Iteration mit dem Operator Ψ(x) := y(a) + (b − a)f(b, x) f¨ ur alle x ∈ V.

Falls f Lipschitz-stetig im zweiten Argument ist, also

kf (b, x

) − f (b, x

)k ≤ L

kx

− x

k f¨ ur alle x

, x

∈ V gilt, erhalten wir

kΨ(x

) − Ψ(x

)k = k(b − a)f(b, x

) − (b − a)f (b, x

)k

3.2 Konvergenz

≤ (b − a)L

kx

− x

k f¨ ur alle x

, x

∈ V,

und der Satz 2.2 von Banach garantiert Konvergenz gegen einen eindeutig bestimmten Fixpunkt x

, falls wir b − a < 1/L

sicherstellen k¨ onnen. Dieser Fixpunkt erf¨ ullt

x

= Ψ(x

) = y(a) + (b − a)f (b, x

),

ist also die gesuchte L¨ osung ˜ y(b) der Gleichung (3.2). Je kleiner b − a wird, desto schnel- ler konvergiert die Fixpunktiteration. Falls f hinreichend oft differenzierbar ist und gute Startwerte bekannt sind, kann man nat¨ urlich statt der Fixpunkt-Iteration auch alterna- tive Ans¨ atze wie beispielsweise das Newton-Verfahren verwenden, um ˜ y(b) zu berechnen.

Entsprechend der Vorgehensweise f¨ ur das explizite Euler-Verfahren k¨ onnen wir die N¨ aherungswerte wieder der Reihe nach berechnen und erhalten die Vorschrift

˜

y(t

) := y

, y(t ˜

) = ˜ y(t

) + hf (t

, y(t ˜

)) f¨ ur alle i ∈ {1, . . . , n}

des impliziten Euler-Verfahrens.

Bemerkung 3.4 (Genauigkeit) Die Analyse des Approximationsfehlers f¨ ur das im- plizite Euler-Verfahren gestaltet sich etwas schwieriger als f¨ ur die explizite Variante. Wir k¨ onnen von

y(b) − y(b) = ˜ y(b) − y(a) − (b − a)f (b, y(b)) ˜

= y(b) − y(a) − (b − a)f (b, y(b)) + (b − a)(f (b, y(b)) − f (b, y(b))) ˜ ausgehen, den ersten Term wie in Lemma 3.1 behandeln und den zweiten mit Hilfe der Lipschitz-Stetigkeit der Funktion f absch¨ atzen, um

ky(b) − y(b)k ≤ ˜ (b − a)

2(1 − L

(b − a)) ky

(η)k (3.3) f¨ ur ein η ∈ [a, b] zu erhalten.

3.2 Konvergenz

Nat¨ urlich ist das Euler-Verfahren nur dann n¨ utzlich, wenn es auch eine hinreichend gute Approximation der tats¨ achlichen L¨ osung berechnet. Wir m¨ ussen also untersuchen, ob und, falls ja, wie schnell die approximative L¨ osung gegen die echte L¨ osung konvergiert.

Dazu gehen wir davon aus, dass die Voraussetzungen des Existenzsatzes 2.3 erf¨ ullt sind.

Die Theorie basiert auf Vergleichen zwischen exakten und approximativen L¨ osungen zu verschiedenen Startwerten: Wir wollen auch L¨ osungen untersuchen k¨ onnen, die zu anderen Startzeitpunkten t

∈ [a, b] von anderen Startwerten y

∈ V ausgehen.

Gl¨ ucklicherweise ist auch deren Existenz durch den erw¨ ahnten Satz gesichert:

Lemma 3.5 (Partielle L¨ osungen) Sei t

∈ [a, b], und sei y

∈ V . Dann existiert eine Funktion y(·; t

, y

) ∈ C

([t

, b], V ), die die Gleichungen

y(t

; t

, y

) = y

, ∂

∂t y(t; t

, y

) = f (t, y(t; t

, y

)) f¨ ur alle t ∈ [t

, b] (3.4) erf¨ ullt. Die Funktion y(·; t

, y

) ist durch diese Gleichungen eindeutig bestimmt.

Beweis. Satz 2.3 angewendet auf dem Teilintervall [t

, b].

Eine wichtige Konsequenz der Eindeutigkeit der L¨ osungen besteht darin, dass zwei partielle L¨ osungen, die in einem Punkt ¨ ubereinstimmen, bereits auf dem gesamten De- finitionsbereich identisch sein m¨ ussen:

Lemma 3.6 (Fortsetzung) Seien t

, s

∈ [a, b] mit t

≤ s

gegeben, und sei y

∈ V . Dann gilt

y(t; t

, y

) = y(t; s

, y(s

; t

, y

)) f¨ ur alle t ∈ [s

, b]. (3.5) Beweis. Sei x

= y(s

; t

, y

) der Wert der Funktion y(·; t

, y

) in dem Zwischenpunkt s

. Nach Definition (3.4) gilt

y(s

; s

, x

) = x

= y(s

; t

, y

),

also stimmen die Funktionen y(·; s

, x

) und y(·; t

, y

) im Punkt s

¨ uberein. Dieselbe Definition beinhaltet auch

∂

∂t y(t; s

, x

) = f (t, y(t; s

, x

)) f¨ ur alle t ∈ [s

, b].

Wegen t

≤ s

gilt außerdem

∂

∂t y(t; t

, y

) = f (t, y(t; t

, y

)) f¨ ur alle t ∈ [s

, b],

also erf¨ ullen y(·; t

, y

) und y(·; s

, x

) auf [s

, b] dieselbe Differentialgleichung mit dem- selben Anfangswert. Nach Satz 2.3 m¨ ussen sie deshalb identisch sein.

Nun ben¨ otigen wir eine ¨ ahnliche Aussage f¨ ur die approximativen L¨ osungen. Wir un- tersuchen ein allgemeines Einschrittverfahren:

Definition 3.7 (Einschrittverfahren) Sei h

∈ R

∪ {∞}. F¨ ur die Menge

∆ := {(t, h) : t ∈ [a, b], h ∈ [0, b − t] ∩ [0, h

]} (3.6) sei eine Funktion

Φ : ∆ × V → V

fixiert. Diese Funktion definiert ein Einschrittverfahren durch

˜

y(t

) = ˜ y(t

) + h

Φ(t

, h

, y(t ˜

)) f¨ ur alle i ∈ {1, . . . , n}, (3.7)

3.2 Konvergenz

falls

(t

, h

) ∈ ∆ f¨ ur alle i ∈ {1, . . . , n}

gilt. Die Funktion Φ bezeichnen wir in diesem Kontext als die Verfahrensfunktion des Einschrittverfahrens, w¨ ahrend wir h

als die maximale Schrittweite bezeichnen.

In der Definition verwenden wir die Menge ∆, um sicher zu stellen, dass f¨ ur jedes Paar (t, h) ∈ ∆ auch der n¨ achste zu berechnende Punkt t + h im Intervall [a, b] enthalten ist.

Damit ist sicher gestellt, dass unsere Algorithmen wohldefiniert sind. Die Schranke h

f¨ ur die Schrittweite ist erforderlich, um beispielsweise bei dem impliziten Euler-Verfahren die L¨ osbarkeit der definierenden Fixpunktgleichung sicher zu stellen. Mit Hilfe einer Verfahrensfunktion k¨ onnen wir, ¨ ahnlich wie in Lemma 3.5, partielle diskrete L¨ osungen definieren:

Definition 3.8 (Diskrete partielle L¨ osungen) F¨ ur beliebige j ∈ {0, . . . , n} und y

∈ V definieren wir analog zu (3.7) die Werte

˜

y(t

; t

, y

) :=



 

 

y

falls i = j,

˜

y(t

; t

, y

)

+h

Φ(t

, h

, y(t ˜

; t

, y

)) ansonsten

f¨ ur i ∈ {j, . . . , n}.

Aus der Definition folgt insbesondere

˜

y(t

; t

, x) − x = ˜ y(t

; t

, x) + h

Φ(t

, h

, y(t ˜

; t

, x)) − x = h

Φ(t

, h

, x), so dass es sich anbietet, die Gleichung

Φ(t, h, x) = y(t ˜ + h; t, x) − x

h f¨ ur alle (t, h) ∈ ∆, x ∈ V (3.8) zu verwenden, um den bisher definierten N¨ aherungsverfahren eine Verfahrensfunktion zuzuordnen.

F¨ ur das explizite Euler-Verfahren erhalten wir Φ(t, h, x) = x + hf (t, x) − x

h = f (t, x) f¨ ur alle (t, h) ∈ ∆, x ∈ V

und d¨ urfen h

= ∞ w¨ ahlen. F¨ ur das implizite Euler-Verfahren ist es etwas schwieriger, Φ zu definieren: F¨ ur (t, h) ∈ ∆ und x ∈ V gilt nach Definition

˜

y(t + h; t, x) = x + hf (t + h, y(t ˜ + h; t, x)), also ergibt sich

Φ(t, h, x) = x + hf (t + h, y(t ˜ + h; t, x)) − x h

= f (t + h, y(t ˜ + h; t, x)) = f (t + h, x + hΦ(t, h, x)).