Numerische Mathematik 2 - Universität Stuttgartharrach/lehre/Numerik2.pdf · Numerische Mathematik...

Numerische Mathematik 2

Prof. Dr. Bastian von Harrach

Universität Stuttgart,

Fachbereich Mathematik - IMNG

Lehrstuhl für Optimierung und inverse Probleme

Sommersemester 2013

http://www.mathematik.uni-stuttgart.de/oip

Inhaltsverzeichnis

1 Gewöhnliche Differentialgleichungen 1

1.1 Einführung und Beispiele . . . . . . . . . . . . . . . . . . . . . 1

1.1.1 Mathematische Beispiele gewöhnlicher DGL . . . . . . 1

1.1.2 Anwendungsbeispiele für gewöhnliche DGL . . . . . . . 3

1.1.3 Anfangswertprobleme . . . . . . . . . . . . . . . . . . . 6

1.1.4 Elementare Lösungsmethoden . . . . . . . . . . . . . . 6

1.2 Theorie gewöhnlicher DGL . . . . . . . . . . . . . . . . . . . . 8

1.2.1 Eine allgemeine Form . . . . . . . . . . . . . . . . . . . 8

1.2.2 Existenz, Eindeutigkeit und Stabilität . . . . . . . . . . 9

1.3 Erste Lösungsmethoden . . . . . . . . . . . . . . . . . . . . . 16

1.3.1 Das Richtungsfeld . . . . . . . . . . . . . . . . . . . . . 16

1.3.2 Explizites Euler-Verfahren . . . . . . . . . . . . . . . . 16

1.3.3 Implizites Euler-Verfahren . . . . . . . . . . . . . . . . 17

1.3.4 Weitere explizite und implizite Methoden . . . . . . . . 18

1.4 Einschrittmethoden höherer Ordnung . . . . . . . . . . . . . . 19

1.4.1 Konsistenz und Konvergenz . . . . . . . . . . . . . . . 20

1.4.2 Runge Kutta Methoden . . . . . . . . . . . . . . . . . 24

1.4.3 Wohldefiniertheit impliziter Methoden . . . . . . . . . 26

1.4.4 Runge-Kutta Ordnungsbedingungen . . . . . . . . . . . 28

1.5 Numerik steifer Differentialgleichungen . . . . . . . . . . . . . 31

1.5.1 Steife Differentialgleichungen . . . . . . . . . . . . . . . 31

1.5.2 Die Testgleichung . . . . . . . . . . . . . . . . . . . . . 33

i

INHALTSVERZEICHNIS

1.5.3 Die Stabilitätsfunktion . . . . . . . . . . . . . . . . . . 34

1.5.4 Stabilität . . . . . . . . . . . . . . . . . . . . . . . . . 35

1.5.5 Nachteile expliziter Verfahren . . . . . . . . . . . . . . 38

1.6 Linear implizite Methoden . . . . . . . . . . . . . . . . . . . . 39

1.7 Mehrschrittverfahren . . . . . . . . . . . . . . . . . . . . . . . 45

1.7.1 Adams-Bashforth Methoden . . . . . . . . . . . . . . . 45

1.7.2 Weitere auf Integration basierende Methoden . . . . . . 46

1.7.3 Auf Differentiation basierende Methoden . . . . . . . . 47

1.7.4 Konvergenz linearer Mehrschrittmethoden . . . . . . . 49

1.8 Randwertprobleme . . . . . . . . . . . . . . . . . . . . . . . . 50

1.8.1 Motivation: Diffusionsprozesse . . . . . . . . . . . . . . 50

1.8.2 Differenzenverfahren . . . . . . . . . . . . . . . . . . . 52

1.8.3 Konsistenz, Stabilität und Konvergenz . . . . . . . . . 54

1.9 Mehrdimensionale Randwertprobleme . . . . . . . . . . . . . . 58

1.9.1 Das Maximumsprinzip . . . . . . . . . . . . . . . . . . 59

1.9.2 Finite Differenzen im Mehrdimensionalen . . . . . . . . 62

1.9.3 Ein diskretes Maximumsprinzip . . . . . . . . . . . . . 65

1.9.4 Konsistenz, Stabilität und Konvergenz . . . . . . . . . 67

2 Eigenwertprobleme 71

2.1 Einleitung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

2.2 Zwei Einschließungssätze . . . . . . . . . . . . . . . . . . . . . 72

2.3 Stabilität des Eigenwertproblems . . . . . . . . . . . . . . . . 76

2.4 Die Potenzmethode . . . . . . . . . . . . . . . . . . . . . . . . 78

2.5 Die Rayleigh-Quotient-Iteration . . . . . . . . . . . . . . . . . 80

2.6 Das QR-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . 84

2.6.1 Motivation: simultane Potenzmethode . . . . . . . . . 84

2.6.2 Das QR-Verfahren . . . . . . . . . . . . . . . . . . . . 87

A Hilfsmittel aus der Analysis 93

A.1 Der mehrdimensionale Mittelwertsatz . . . . . . . . . . . . . . 93

ii

Kapitel 1

Gewöhnliche

Differentialgleichungen

In diesem Kapitel beschäftigen wir uns mit der Lösung gewöhnlicher Differen-tialgleichungen (engl.: ordinary differential equations, ODE) bzw. Anfangs-wertproblemen. Wir setzen dabei keinerlei Vorkenntnisse über die Theoriegewöhnlicher Differentialgleichungen voraus.

Wir beginnen mit einigen Beispielen.

1.1 Einführung und Beispiele

1.1.1 Mathematische Beispiele gewöhnlicher DGL

Differentialgleichung: Gleichung die eine unbekannte Funktion zusammenmit ihren Ableitungen enthält.

Beispiel 1.1Einfache mathematische Beispiele für Differentialgleichungen:

(a) Betrachte die Differentialgleichung

y′(x) = 0 (kurz: y′ = 0),

d.h. gesucht ist eine differenzierbare Funktion

y : R→ R, y : x 7→ y(x) mit y′(x) = 0 ∀x ∈ R.

1

KAPITEL 1. GEWÖHNLICHE DIFFERENTIALGLEICHUNGEN

• Spezielle Lösungen sind z.B.:

y(x) = 0, y(x) = 1, y(x) = −37, . . .

• Die allgemeine Lösung ist y(x) = C, C ∈ R, d.h. man kannzeigen, dass jede solche Funktion die ODE löst und jede Lösungin dieser Form geschrieben werden kann.

(b) Betrachte die Differentialgleichung

y′(x) = ry(x), r ∈ R (kurz: y′ = ry).


y(x) = 0, y(x) = erx, y(x) = −37erx, . . .

• Die allgemeine Lösung ist y(x) = Cerx, C ∈ R.

(c) Betrachte die Differentialgleichung höherer Ordnung

y′′(x) = −y(x) (kurz: y′′ = −y).


y(x) = 0, y(x) = sin(x), y(x) = 37 cos(x), . . .

• Die allgemeine Lösung ist C1 sin(x) + C2 cos(x), C1, C2 ∈ R.

(d) Betrachte das Differentialgleichungssystem

y′′1(x) = −y1(x)

y′2(x) = y3(x)

y′3(x) = y3(x)

• Eine spezielle Lösung ist z.B.:

y1(x) = sin(x)

y2(x) = 37

y3(x) = 0

• Die allgemeine Lösung ist

y1(x) = C1 sin(x) + C2 cos(x)

y2(x) = C3ex + C4

y3(x) = C3ex

mit C1, . . . , C4 ∈ R.

2

1.1. EINFÜHRUNG UND BEISPIELE

(e) Betrachte die Differentialgleichung für eine vektorwertige Funktion:

Gesucht ist (eine differenzierbare Funktion)

y : R→ R

2, y : x 7→ y(x) :=

(y1(x)y2(x)

)

.

mit y′ = y, also

y′(x) = y(x) ⇐⇒(

y′1(x)y′2(x)

)

=

(y1(x)y2(x)

)

⇐⇒

y′1(x)= y1(x)y′2(x)= y2(x)

Gewöhnliche Differentialgleichung: Alle Ableitung beziehen sich auf dieselbe eindimensionale Variable.

Beispiel 1.2Betrachte y : R2 → R, (x1, x2) 7→ y(x) = y(x1, x2)

Eine nicht gewöhnliche, partielle Differentialgleichung (engl.: Partial Diffe-rential Equation, PDE) ist

∂2y

∂x21

+∂2y

∂x22

= 0

Spezielle Lösungen sind z.B.: y(x) = 0, y(x) = x1 + x2, . . .

PDEs sind Gegenstand der im nächsten Semester angebotenen VorlesungNumerik partieller Differentialgleichungen.

1.1.2 Anwendungsbeispiele für gewöhnliche DGL

Je nach Anwendung verwenden wir in diesem Abschnitt unterschiedliche Be-zeichner für die Funktion und die Variable. Für Ableitungen bezüglich derZeit schreiben wir auch y(t) = y′(x).

Continuous compounding – stetige Verzinsung

y(t): Sparguthaben zum Zeitpunkt t.

Bank zahlt Zinsen proportional zu Guthaben und Zeitspanne:

∆y(t) = y(t+∆t)− y(t) = ry(t)∆t

Gutgeschriebene (kapitalisierte) Zinsen generieren Zinseszinsen.

Stetige Verzinsung: instantane Kapitalisierung, d.h. ∆t → 0:

y(t) = ry(t).

3


Populationsdynamik

y(t): Anzahl von Individuen einer Population

• Model 1 (Malthus): Population wächst mit konstanter Rate r ∈ R. Ineinem kurzen Zeitintervall ∆t, wächst die Population um ∆y ≈ ry∆t.

y = ry.

• Model 2 (Verhulst): Es existiert Maximalbevölkerung M > 0, Popula-tion wächst mit von y abhängiger Rate

y << M : Wachstum mit Rate r

y ≈ M : Wachstum mit Rate 0

In Zeitintervall ∆t, wächst die Population um ∆y ≈ r(1− y/M)y∆t.

y = r(

1− y

M

)

y = ry − r

My2.

rMy2: Todesrate aufgrund zu hoher Population.

• Model 3 (Lotka-Volterra / Räuber-Beute Modell):

Zwei Populationen y1(t), y2(t), Räuber und Beutetiere:

y1 = r1y1 − f1y1y2

y2 = −r2y2 + f2y1y2

Chemische Reaktionen

A(t), B(t), C(t), . . . : Konzentration der Chemikalien A,B,C,. . .

• Ak−→ B

In einem kurzen Zeitintervall ∆t wandeln sich kA(t)∆t Moleküle vonA in B um.

A(t) = −kA(t), B(t) = kA(t).

• A+Bk−→ C + 2D

A = −kAB, B = −kAB, C = kAB, D = 2kAB.

4


Newtonsche Gesetze

x(t) = (x1(t), x2(t), x3(t))T : Position eines Körpers zum Zeitpunkt t

v(t) = x(t) = (x1(t), x2(t), x3(t))T : Geschwindigkeit

a(t) = v(t) = x(t) = (x1(t), x2(t), x3(t))T : Beschleunigung

Newtonsches Gesetz: Wirkt eine Kraft F (t) auf einen Körper der Masse m,so ist

F (t) = ma(t) = mx(t)

Elektrische Schaltkreise Betrachte eine RLC-Reihenschaltung, d.h. eineReihenschaltung von Widerstand (R), Induktivität (L), und Kapazität (C).

• Stromstärke IR(t) durch einen Widerstand R bei Anlegen einer Span-nung UR(t):

IR(t) = UR(t)/R.

• Stromstärke IC(t) durch einen Kondensator mit Kapazität C bei An-legen einer Spannung UC(t):

IC(t) = UC(t)C

• Stromstärke IL(t) durch eine Spule mit Induktivität L bei Anlegeneiner Spannung UL(t):

IL(t) = UL(t)/L.

Aus den Kirchhoffschen Gesetzen folgt

U(t) = UR(t) + UC(t) + UL(t) (Spannungsbilanz)

IL(t) = IC(t) = IR(t) (Strombilanz)

Es folgt, dass UL(t) = CLUC(t) und UR(t) = RCUC(t), und damit

U(t) = UC(t) +RCUC(t) + LCUC(t).

5


1.1.3 Anfangswertprobleme

Die Lösung einer gewöhnlichen DGL ist üblicherweise nicht eindeutig. Dieallgemeine Lösung von y(t) = ry(t) ist y(t) = Cert mit einem ParameterC ∈ R. In vielen Anwendungen sind die Parameter eindeutig bestimmt durchdie Anfangswerte von y. Bei der stetigen Verzinsung ist z.B. C = y(0) dasanfängliche Sparguthaben.

Intuitiv erwarten wir in den anderen Beispielen, dass die Lösung durch fol-gende Informationen eindeutig bestimmt wird:

• Populationsdynamik: anfängliche Population y(0), bzw., y1(0), y2(0)

• Chemische Reaktionen: anfängliche Konzentrationen A(0), B(0), . . .

• Newtonsche Gesetze: anfängliche Position x1(0), x2(0), x3(0) und Ge-schwindigkeit v1(0), v2(0), v3(0)

• Elektrische Schaltkreise: Anfängliche Strom- und Spannungswerte amKondensator UC(0) and UC(0)

Eine gewöhnliche DGL zusammen mit Anfangsbedingungen heißt auch An-fangswertproblem (AWP).

1.1.4 Elementare Lösungsmethoden

Ähnlich wie Integrale lassen sich manche (aber nicht alle) gewöhnliche Dif-ferentialgleichungen analytisch, d.h. in geschlossener Form lösen. Wir gebenhier nur Beispiele besonders einfacher Lösungsmethoden an. Es existierennoch einige weitere wichtige Lösungsmethoden, aber im Allgemeinen lassensich gewöhnliche Differentialgleichungen nur numerisch lösen.

Raten/Wissen der Lösung siehe Beispiel 1.1.

Separation der Variablen Gewöhnliche DGL der Form

y′(x) = g(x)h(y(x))

lassen sich formal(!) schreiben als

dy

h(y)= g(x) dx (1.1)

6


und durch Integration lösen∫

1

h(y)dy =

∫

g(x) dx. (1.2)

Formal bedeutet dabei, dass dies keine mathematisch rigorosen Umformulie-rungen sind. Die Ausdrücke dx sind nicht definiert!

Man kann diese Methode rigoros formulieren und rechtfertigen. Aber auchohne rigorose Rechtfertigung haben formale Methoden einen großen Nutzen(nicht nur im Bereich gewöhnlicher DGL). Oft lässt sich nämlich durch for-males Vorgehen ein Lösungskandidat bestimmen und für diesen dann rigorosüberprüfen, ob er tatsächlich eine Lösung ist.

Beispiel 1.3Betrachte Modell 2 aus der Populationsdynamik (M = r = 1):

dy

dt= (1− y) y

mit Anfangwert y(0) > 1.

Obiges formales Vorgehen liefert∫

1

(1− y) ydy =

∫

1 dt.

Wir erwarten anschaulich, dass die Population von oben gegen ihren Ma-ximalwert M = 1 konvergieren diesen aber nicht unterschreiten wird. Wirlösen also die Integrale auf beiden Seiten unter der Annahme y > 1:

∫

1 dt = t+ const.∫

1

(1− y) ydy = −

∫1

y − 1dy +

∫1

ydy

= − ln(y − 1) + ln(y) + const. = lny

y − 1+ const.

Wir erhalten

lny

y − 1= t + const. =⇒ y

y − 1= Cet

=⇒ y =Cet

Cet − 1.

7


C kann aus dem Anfangswert y(0) = y0 bestimmt werden:

y0 =C

C − 1=⇒ C =

y0y0 − 1

.

Man prüft leicht (und mathematisch rigoros) nach, dass der so bestimmteLösungskandidat für y0 > 1 tatsächlich das AWP

dy

dt= (1− y) y, y(0) = y0

löst.

Variation der Konstanten Betrachte y′(x) = ry(x) + z(x).

Allgemeine Lösung der homogenen Gleichung y′(x) = ry(x) ist y(x) = Cerx

mit einer Konstante C ∈ R.

Ansatz: Ersetze zur Lösung der inhomogenen Gleichung die Konstante durcheine Funktion C(x):

C ′(x)erx + C(x)rerx = y′(x) = ry(x) + z(x) = rC(x)erx + z(x)

=⇒ C ′(x) = z(x)e−rx

Durch Integration erhalten wir C(x) und damit die Lösung y(x) = C(x)erx.

1.2 Theorie gewöhnlicher DGL

1.2.1 Eine allgemeine Form

Von nun an betrachten wir stets Anfangswertprobleme in der folgenden all-gemeinen Form

y′(x) = f(x, y(x)), y(x0) = y0

wobei x ∈ [x0,∞) ⊂ R, y(x) = (y1(x), . . . , yd(x))T ∈ R

d vektorwertig ist,d ∈ N, und

f : Rd+1 → R

d, f(x, y) =

f1(x, y1, . . . , yd)...

fd(x, y1, . . . , yd)

.

8

1.2. THEORIE GEWÖHNLICHER DGL

Die Differentialgleichung y′ = f(x, y(x)) ist äquivalent zum Differentialglei-chungssystem

y′1(x) = f1(x, y1(x), . . . , yd(x))...

y′d(x) = fd(x, y1(x), . . . , yd(x))

Gleichungen höherer Ordnung (d.h. solche, die höhere Ableitungen von yenthalten) können in diese Form transformiert werden, indem y und seineAbleitungen (bis zur zweithöchsten) in einer vektorwertigen Hilfsfunktionu = (u1, u2, . . .) zusammengefasst werden

u1(x) := y(x), u2(x) := y′(x), u3(x) = y′′(x), usw.

Beispiel 1.4y′′ = −y kann in obige Form transformiert werden durch

u =

(u1(x)u2(x)

)

:=

(y(x)y′(x)

)

.

Damit ist y′′ = −y äquivalent zu

u′ =

(u′1(x)

u′2(x)

)

=

(y′(x)y′′(x)

)

=

(y′(x)−y(x)

)

=

(u2(x)−u1(x)

)

=: f(x, u(x)).

1.2.2 Existenz, Eindeutigkeit und Stabilität

Ein Problem heißt wohlgestellt (nach Hadamard) wenn

(a) eine Lösung existiert (Existenz ),

(b) die Lösung eindeutig ist (Eindeutigkeit),

(c) die Lösung stetig von den Eingabeparametern abhängt (Stabilität).

Für das Anfangswertproblem

y′(x) = f(x, y(x)), y(x0) = y0.

bedeutet Wohlgestelltheit, dass eine eindeutige Lösung y(x) existiert und die-se Lösung stetig von den Anfangswerten y0 (und ggf. weiteren in f vorhande-nen Parametern) abhängt. Dies werden wir in diesem Abschnitt untersuchen.

9


Existenz und Eindeutigkeit. In diesem Abschnitt beweisen wir den Exis-tenz- und Eindeutigkeitssatz von Picard-Lindelöf in seiner einfachsten Form

Satz 1.5 (Picard-Lindelöf)Seien x0, xend ∈ R, x0 < xend, y0 ∈ Rd. Für

f : [x0, xend]×Rd → R

d, f : (x, y) 7→ f(x, y)

gelte

(a) f ist stetig

(b) f ist (global und bzgl. x gleichmäßig) Lipschitz-stetig in y, d.h. es existiertein L > 0 so dass

‖f(x, y)− f(x, z)‖ ≤ L ‖y − z‖ ∀x ∈ [x0, xend], y, z ∈ Rd.

Dann existiert genau eine stetig differenzierbare Funktion y : [x0, xend] → R

d

mity′(x) = f(x, y(x)) ∀x ∈ (x0, xend) und y(x0) = y0.

Wir werden Satz 1.5 mit Hilfe des Banachschen Fixpunktsatz und einigenHilfssätzen beweisen.

Satz 1.6 (Banachscher Fixpunktsatz)Sei (X, d) ein vollständiger metrischer Raum und Φ eine kontrahierendeSelbstabbildung von X, d.h. es existiere ein q < 1 mit

Φ : X → X und d(Φ(x),Φ(y)) ≤ q d(x, y) ∀x, y ∈ X.

Dann besitzt Φ genau einen Fixpunkt x, d.h. genau ein x ∈ K mit Φ(x) = x.

Für jeden Startwert x(0) ∈ X konvergiert die durch Fixpunktiteration defi-nierte Folge

(x(k))k∈N0, x(k+1) := Φ(x(k)) ∀k ∈ N0.

(in der Metrik d) gegen x, d.h.

d(x(k), x) → 0.

Beweis: aus der Numerik I bekannt.

Als vollständigen metrischen Raum werden wir den Raum stetiger Funktio-nen betrachten:

10


Lemma 1.7Die Menge der auf [x0, xend] stetigen Funktionen

C([x0, xend])d := y : [x0, xend] → R

d stetigist bezüglich der Supremums- (auch: Maximumsnorm)

‖y‖∞ := maxx∈[x0,xend]

‖y(x)‖

ein Banachraum, d.h. ein vollständiger normierter Vektorraum (und damitinsbesondere bzgl. der induzierten Metrik d(y1, y2) := ‖y1 − y2‖∞ ein voll-ständiger normierter Raum).

Beweis: Die Vektorraum- und Normeigenschaften lassen sich einfach nach-rechnen (siehe Übungsaufgabe 2.2).

Zum Beweis der Vollständigkeit müssen wir zeigen, dass (bzgl. ‖·‖∞) jedeCauchy-Folge konvergiert. Sei also (y(n))n∈N ⊂ C([x0, xend]) eine Cauchy-Folge bzgl. ‖·‖∞, d.h. für jedes ǫ > 0 existiert N ∈ N, so dass

ǫ >∥∥y(n) − y(m)

∥∥ = max

x∈[x0,xend]

∥∥y(n)(x)− y(m)(x)

∥∥ ∀n,m ≥ N.

Für jedes feste x ∈ [x0, xend] ist dann

ǫ >∥∥y(n)(x)− y(m)(x)

∥∥ ∀n,m ≥ N,

also (y(n)(x))n∈N ⊂ Rd eine Cauchy-Folge. Diese besitzt einen Grenzwert, sodass wir

y : [x0, xend] → R

d durch y(x) := limn→∞

y(n)(x).

definieren können.

Ausǫ >

∥∥y(n)(x)− y(m)(x)

∥∥ ∀n,m ≥ N

folgt dassǫ ≥

∥∥y(n)(x)− y(x)

∥∥ ∀n ≥ N

und so gibt es also zu jedem ǫ > 0 ein N ∈ N, so dass

ǫ ≥ maxx∈[x0,xend]

∥∥y(n)(x)− y(x)

∥∥ .

Die Funktionenfolge y(n) konvergiert also gleichmäßig gegen y, womit dieStetigkeit von y folgt, also y ∈ C([x0, xend])

d. Damit gilt also∥∥y(n) − y

∥∥∞

= maxx∈[x0,xend]

∥∥y(n)(x)− y(x)

∥∥→ 0,

d.h. y(n) → y bzgl. der Supremumsnorm.

11


Wir benötigen aber noch eine etwas abgewandelte Norm:

Folgerung 1.8Für eine Funktion w ∈ C([x0, xend]) gelte

∃c, C > 0 : c ≤ w(x) ≤ C ∀x ∈ [x0, xend].

Dann ist C([x0, xend])d auch bzgl. der gewichteten Supremumsnorm

‖y‖w := maxx∈[x0,xend]

(w(x) ‖y(x)‖)

ein Banachraum.

Beweis: Offenbar ist ‖·‖w tatsächlich eine Norm und es gilt

c ‖y‖∞ ≤ ‖y‖w ≤ C ‖y‖∞ ∀y ∈ C([x0, xend])d.

Insbesondere ist jede Cauchy-Folge bzgl. ‖·‖w auch eine bzgl. ‖·‖∞ und derGrenzwert bzgl. ‖·‖∞ ist auch der Grenzwert bzgl. ‖·‖w.

Wir formulieren noch die Anfangswertaufgabe in eine Fixpunktgleichung um.

Lemma 1.9Es gelten die Voraussetzungen von Satz 1.5. Eine stetig differenzierbare Funk-tion y : [x0, xend] → R

d erfüllt genau dann das AWP

y′(x) = f(x, y(x)) ∀x ∈ (x0, xend) und y(x0) = y0,

wenn y stetig ist und die folgende Fixpunktgleichung löst:

y(x) = y0 +

∫ x

x0

f(t, y(t)) dt ∀x ∈ [x0, xend].

Beweis: Ist y stetig differenzierbar und erfüllt das AWP, so ist

y(x) = y(x0) +

∫ x

x0

y′(t) dt = y0 +

∫ x

x0

f(t, y(t)) dt ∀x ∈ [x0, xend].

Ist umgekehrt y stetig und erfüllt die Fixpunktgleichung, so ist auch t 7→f(t, y(t)) stetig und

y(x) = y0 +

∫ x

x0

f(t, y(t)) dt

ist differenzierbar, y′(x) = f(x, y(x)) und y(x0) = y0.

12


Jetzt können wir Satz 1.5 beweisen:

Beweis von Satz 1.5: Nach unserer Vorarbeit genügt es zu zeigen, dassdurch

y 7→ Φ(y), Φ(y) : x 7→ y0 +

∫ x

x0

f(t, y(t)) dt

eine bzgl. einer gewichteten Supremumsnorm kontrahierende Selbstabbildungvon C([x0, xend])

d ist.

Die Selbstabbildungseigenschaft ist klar. Für die Kontraktionseigenschaft be-trachten wir für zwei Funktionen y(1), y(2) ∈ C([x0, xend])

d

∥∥Φ(y(1))(x)− Φ(y(2))(x)

∥∥

=

∥∥∥∥

∫ x

x0

(f(t, y(1)(t))− f(t, y(2)(t))

)dt

∥∥∥∥

≤∫ x

x0

∥∥f(t, y(1)(t))− f(t, y(2)(t))

∥∥ dt ≤

∫ x

x0

L∥∥y(1)(t)− y(2)(t)

∥∥ dt

Hieraus folgt

∥∥Φ(y(1))− Φ(y(2))

∥∥∞

≤ L(xend − x0)∥∥y(1) − y(2)

∥∥∞,

so dass Φ nur für kleine L oder nah an x0 liegendes xend eine Kontraktionist.

Durch Einfügen einer Gewichtsfunktion w(x) (mit den in Folgerung 1.8 ge-nannten Eigenschaften) erhalten wir jedoch

w(x)∥∥Φ(y(1))(x)− Φ(y(2))(x)

∥∥

≤ w(x)

∫ x

x0

L1

w(t)w(t)

∥∥y(1)(t)− y(2)(t)

∥∥ dt

≤ w(x)L∥∥y(1) − y(2)

∥∥w

∫ x

x0

1

w(t)dt

und damit

∥∥Φ(y(1))− Φ(y(2))

∥∥w≤∥∥y(1) − y(2)

∥∥wL max

x∈[x0,xend]

(

w(x)

∫ x

x0

1

w(t)dt

)

.

Φ ist also eine Kontraktion bzgl. ‖·‖w wenn wir eine Gewichtsfunktion wfinden mit

L

(

w(x)

∫ x

x0

1

w(t)dt

)

< 1 ∀x ∈ [x0, xend].

13


Offenbar gilt dass

e−a(x−x0)

∫ x

x0

1

e−a(t−x0)dt ≤ 1

a.

Mit w(x) = e−L/2(x−x0) ist also Φ eine Kontraktion (mit Kontraktionskon-stante 1/2) bzgl. ‖·‖w, womit Satz 1.5 bewiesen ist.

Beispiel 1.10Auf die Voraussetzung der Lipschitz-Stetigkeit kann nicht verzichtet werden,wie die folgenden Beispiele zeigen:

(a) Betrachte das AWPy′ =

√y, y(0) = 0.

Die rechte Seite f(x, y) :=√y ist nicht Lipschitz-stetig. Tatsächlich ist

die Lösung des AWP nicht eindeutig. Zwei Lösungen sind z.B.

y(x) = 0 und y(x) =

0 für 0 ≤ x ≤ 1,14(x− 1)2 für x > 1.

(b) Betrachte das AWPy′ = y2, y(0) = 1.

Die rechte Seite f(x, y) = y2 ist nur lokal aber nicht global Lipschitzstetig. Man kann zeigen, dass nur

y(x) =1

1− x

das AWP lösen kann. Die Lösung existiert also nur auf dem Intervall[0, 1).

Bemerkung 1.11Beispiel 1.10(b) zeigt in gewisser Weise das Schlimmste, was für eine nurlokal Lipschitz-stetige rechte Seite passieren kann. f erfülle die Lipschitz-Bedingung∥∥f(x, y(1))− f(x, y(2))

∥∥ ≤ L

∥∥y(1) − y(2)

∥∥ ∀x ∈ [x0, xend], y(1), y(2) ∈ Q.

nur in einem Rechteck Q := (a1, b1)× · · · × (ad, bd). Dann kann man zeigen,dass für jeden Anfangswert y0 ∈ Q, ein nichtleeres Teilintervall [x0, x1] exi-stiert auf dem genau eine Lösung existiert. Darüber hinaus kann man zeigen,dass sich die Lösungskurve (x, y(x)) bis zum Rand des Rechtecks eindeutigfortsetzen lässt.

14


Stabilität Nun untersuchen wir wie sich eine Störung der Anfangswerteauf die Lösung auswirkt.

Satz 1.12Es gelten weiterhin die Voraussetzungen von Satz 1.5. y und z seien zweiLösungen der gleichen DGL, aber mit verschiedenen Anfangswerten, also

y′(x) = f(x, y(x)), y(x0) = y0,

z′(x) = f(x, z(x)), z(x0) = z0.

Dann gilt

‖y(x)− z(x)‖ ≤ eL(x−x0) ‖y0 − z0‖ ∀x ∈ [x0, xend].

Beweis: Betrachte die Differenz

s(x) := ‖y(x)− z(x)‖2 = (y(x)− z(x))T (y(x)− z(x)).

Es ist

s′(x) = 2(y′(x)− z′(x))T (y(x)− z(x))

= 2(f(x, y(x))− f(x, z(x)))T (y(x)− z(x))

≤ 2 ‖f(x, y(x))− f(x, z(x))‖ ‖y(x)− z(x)‖≤ 2L ‖y(x)− z(x)‖2 = 2Ls(x).

Im Fall y0 = z0 gilt nach Satz 1.5 y(x) = z(x) für alle x ≥ x0 und dieBehauptung ist bewiesen. Ansonsten existiert ein größtmögliches Intervall[x0, x0+ δ) auf dem s(x) 6= 0 gilt und im Falle x0+ δ < xend ist s(x0+ δ) = 0(siehe Übungsaufgabe 2.5). Es genügt die Behauptung in diesem Intervall zuzeigen, da für s(x0 + δ) = 0 nach Satz 1.5 y(x) und z(x) ab x ≥ x0 + δübereinstimmen.

Tatsächlich gilt für alle x ∈ [x0, x0 + δ)

d

dxln s(x) =

s′(x)

s(x)≤ 2L

also

ln s(x) =

∫ x

x0

d

dtln s(t) dt+ ln s(x0) ≤ 2L(x− x0) + ln s(x0)

d.h.

‖y(x)− z(x)‖2 = s(x) ≤ e2L(x−x0)s(x0) = e2L(x−x0) ‖y0 − z0‖2 ,

womit die Behauptung gezeigt ist.

15


Bemerkung 1.13(a) Auch wenn f nur lokal Lipschitz-stetig, gilt das Stabilitätsresultat in Satz

1.12 noch dort, wo die eindeutige Existenz der Lösungen gesichert ist.

(b) Wir werden in dieser Vorlesung im Folgenden annehmen, dass f belie-big oft stetig differenzierbar (und damit inbesondere in jeder kompaktenMenge Lipschitz-stetig) ist, so dass das AWP

y′(x) = f(x, y(x)), y(x0) = y0.

zumindest in einem hinreichend kleinen Intervall [x0, xend] wohlgestelltist, also eine eindeutige Lösung existiert und diese (im Sinne von Satz1.12) stabil vom Anfangswert abhängt.

Außerdem kann man zeigen, dass die Lösung des AWP dann ebenfallsunendlich oft stetig differenzierbar ist.

1.3 Erste Lösungsmethoden

1.3.1 Das Richtungsfeld

Zur anschaulichen Herleitung einfacher erster Lösungsmethoden betrachtenwir ein skalares AWP, in dem wir die Lösung y : [x0,∞) → R von

y′(x) = f(x, y(x)), y(x0) = y0 ∈ Rsuchen.

Wir können uns die DGL y′(x) = f(x, y(x)) durch das dazugehörige Rich-tungsfeld veranschaulichen: Zu jedem Punkt (x, y) ∈ R2 zeichnen wir einenRichtungspfeil mit Steigung f(x, y), z.B. den Vektor (1, f(x, y))T (vgl. die inder Vorlesung gezeichnete Skizze). Eine Funktion löst die DGL genau dann,wenn an jedem Punkt durch den die Funktion geht, die Steigung der Funktionund die Steigung des Richtungspfeils übereinstimmen. Wir können die DGLzeichnerisch lösen, indem wir ausgehend vom Startwert (x, y0) die Funktionpassend zu den Richtungspfeilen zeichnen.

1.3.2 Explizites Euler-Verfahren

Basierend auf dieser zeichnerischen Idee gehen wir nun systematischer vorund entwickeln ein erstes numerisches Lösungsverfahren. Betrachte das all-gemeine (vektorwertige) AWP

y′(x) = f(x, y(x)), y(x0) = y0 ∈ Rd

16

1.3. ERSTE LÖSUNGSMETHODEN

auf dem Intervall x ∈ [a, b], a = x0. Wir diskretisieren das Intervall durchn+ 1 Punkte

a = x0 < x1 < . . . < xn = b,

also unter Verwendung der Schrittweite hi := xi+1 − xi, i = 0, . . . , n − 1.Beginnend mit x0 (wo wir die Lösung kennen, y(x0) = y0) berechnen wir nunsukzessiv Approximationen yi ≈ y(xi).

Der einfachste Möglichkeit ist die explizite Eulermethode, bei der wie dieSteigung im aktuellen Punkt verwenden, um die Approximation im nächstenPunkt zu berechnen:1

y1 := y0 + h0f(x0, y0)

y2 := y1 + h1f(x1, y1)

...

yi+1 := yi + hif(xi, yi)

...

yn := yn−1 + hn−1f(xn−1, yn−1)

Dies entspricht dem zeichnerischen Lösen der DGL im Richtungsfeld durcheine stückweise lineare Funktion, bei der die Steigung der Liniensegmente imlinken Punkt mit dem Richtungsfeld übereinstimmt.

Das gleiche Verfahren erhalten wir auch durch Diskretisierung der DGL mitfiniten Differenzen (Vorwärtsdifferenzenquotient):

yi+1 − yixi+1 − xi

≈ y(xi+1)− y(xi)

xi+1 − xi≈ y′(xi) = f(xi, y(xi)) ≈ f(xi, yi).

Das explizite Eulerverfahren heißt auch Vorwärts-Euler-Verfahren (engl.: fo-ward Euler).

1.3.3 Implizites Euler-Verfahren

Bei der zeichnerischen Lösung der DGL im Richtungsfeld durch eine stückwei-se lineare Funktion, könnten wir auch versuchen die Funktion so zu wählen,

1Hier und im Folgenden bezeichen wir mit y0, y1, . . . , yn ∈ Rd d-dimensionale Vektorenund nicht die Einträge eines Vektors.

17


dass die Steigung der Liniensegmente im rechten Punkt mit dem Richtungs-feld übereinstimmt.

Dann haben wir keine explizite Formel für die Berechnung von yi+1 aus yi.Statt dessen ist yi+1 implizit gegeben als Lösung von

yi+1 = yi + hif(xi+1, yi+1).

Bei diesem impliziten Euler-Verfahren müssen wir also für jedes yi ∈ Rd, i =2, . . . , n, ein (üblicherweise nicht-lineares) d-dimensionales Gleichungssytemlösen.

Wieder erhalten wir das Verfahren auch durch Diskretisierung der DGL mitfiniten Differenzen, diesmal mit dem Rückwärtsdifferenzenquotienten:

yi+1 − yixi+1 − xi

≈ y(xi+1)− y(xi)

xi+1 − xi≈ y′(xi+1) = f(xi+1, y(xi+1)) ≈ f(xi+1, yi+1).

Enstprechend heißt das implizite Euler-Verfahren auch Rückwärts-Euler-Ver-fahren (engl.: backward Euler).

Bemerkung 1.14Wie wir noch sehen werden, besitzen implizite Verfahren für manche (diesogenannten steifen) Differentialgleichungen so große Vorteile, dass sie denerhöhten Aufwand wert sind.

1.3.4 Weitere explizite und implizite Methoden

Die anschauliche Idee, stückweise lineare Funktion ins Richtungsfeld zu zeich-nen, führt auf viele weitere explizite und implizite Methoden:

Implizite Mittelpunktsregel. Zuerst zeichnen wir die Liniensegmente so,dass die Steigung der Segmente in ihrem Mittelpunkt mit dem Richtungsfeldübereinstimmt.

Wir bestimmen also erst yi+1/2 durch Lösung der Gleichung

yi+1/2 = yi +hi

2f

(xi+1 + xi

2, yi+1/2

)

und setzen dann

yi+1 := yi + hif

(xi+1 + xi

2, yi+1/2

)

Dies ist die sogenannte implizite Mittelpunktsregel, die auch als Kombinationeines halben impliziten Eulerschrittes mit einem halben expliziten Euler-schritt interpretiert (und implementiert) werden kann.

18

1.4. EINSCHRITTMETHODEN HÖHERER ORDNUNG

Verfahren von Runge. Wir können auch versuchen, eine explizite Formelfür yi+1/2 zu finden. Dafür setzen wir

yi+1/2 := yi +hi

2f(xi, yi)

und wählen dann

yi+1 := yi + hif

(xi+1 + xi

2, yi+1/2

)

.

Dies ist das Verfahren von Runge. Beachte dass dies nicht nur die Kombina-tion zweier Halbschritte des expliziten Eulerverfahrens ist.

Crank-Nicolson Methode. Die nächste Methode erhalten wir durch dieForderung, dass die Steigung der Liniensegmente mit dem Mittelwert derSteigungen im Richtungsfeld im linken und rechten Randpunkt des Segmentsübereinstimmen soll:

yi+1 = yi + hif(xi, yi) + f(xi+1, yi+1)

2.

Dies ist die Crank-Nicolson Methode.

Verfahren von Heun. Eine explizite Alternative zur Crank-Nicolson Me-thode erhalten wir, indem wir yi+1 auf der rechten Seite der Crank-Nicolson-Formel durch einen Schritt mit dem expliziten Euler-Verfahren approximie-ren:

η := yi + hif(xi, yi)

yi+1 := yi + hif(xi, yi) + f(xi+1, η)

2.

1.4 Einschrittmethoden höherer Ordnung

Generalvoraussetzung: In diesem Abschnitt sei [a, b] ⊂ R, b > a stetsein festes Intervall. Wir werden im Folgenden nur Differentialgleichungen

y′(x) = f(x, y(x))

mit unendlich oft differenzierbarer rechter Seite f betrachten. Insbesondereist f also lokal Lipschitz stetig.

19


Der Einfachheit halber fordern wir zusätzlich dass alle partiellen Ablei-tungen (jeder Ordnung) von f in [a, b]×Rd beschränkt sind. Insbesondereexistiert ein L > 0 so dass

‖fy(x, y)‖ ≤ L ∀(x, y) ∈ [a, b]×Rd,

wobei

fy(x, y) :=∂f

∂y(x, y) =

(∂fi∂yj

(x, y)

)d

i,j=1

∈ Rd×d.

Damit gilt also für alle x ∈ [a, b] und y, z ∈ Rd

‖f(x, y)− f(x, z)‖ =

∥∥∥∥

∫ 1

0

fy(x, y + t(z − y))(z − y) dt

∥∥∥∥≤ L ‖y − z‖

und nach Abschnitt 1.2.2 ist somit für jede Wahl der Anfangswerte x0 ∈[a, b] und y0 ∈ R

d die eindeutige Existenz von Lösungen und ihre stetigeAbhängigkeit von den Anfangswerten garantiert.

Außerdem kann man zeigen, dass dann alle Ableitungen jeder Lösung einesAWP mit rechter Seite f beschränkt sind. Präzise ausgedrückt: Zu jedersolchen rechten Seite f existieren Konstanten Ck > 0 (k ∈ N0), so dass fürjede Wahl der Anfangswerte x0 ∈ [a, b] und y0 ∈ Rd die zugehörige Lösungvon

y′(x) = f(x, y(x)), y(x0) = y0 ∈ Rd

erfüllt, dasssup

x∈[x0,b]

∥∥y(k)(x)

∥∥ ≤ Ck.

Auf diese globale Zusatzvoraussetzung kann verzichtet werden, wenn die Lös-barkeit sichergestellt ist und die Approximationen an die Lösung beschränktbleiben.

1.4.1 Konsistenz und Konvergenz

Zur Lösung des AWP

y′(x) = f(x, y(x)), y(x0) = y0 ∈ Rd

diskretisieren wir das Intervall in n+ 1 Punkte

a = x0 < x1 < . . . < xn = b.

Beginnend mit x0 (wo wir die Lösung y(x0) = y0 kennen) versuchen wirsukzessive Approximationen yi ≈ y(xi) ∈ Rd zu bestimmen.

20


Einschrittmethoden: Nur der vorherige Punkt yi (und xi, xi+1, hi :=xi+1 − xi) wird zur Berechnung von yi+1 verwendet.

Mehrschrittmethoden: Mehrere vorherige Punkte yi, yi−1, . . . yi−m (undxi+1, xi, . . . , xi−m ) werden zur Berechnung von yi+1 verwendet. (m+1-Schritt Methode benötigt Startprozedur zur Berechnung der ersten mWerte y1,. . . , ym).

Die Methoden in Abschnitt 1.3 sind allesamt Einschrittmethoden.Definition 1.15Eine Einschrittmethode heißt konsistent, falls für jede (unsere Generalvor-aussetzung erfüllende) rechte Seite f gilt, dass

limh→0

supxi∈[a,b],yi∈Rd

|yi+1 − y(xi + h)|h

= 0

wobei y die Lösung des AWP

y′(x) = f(x, y(x)), y(xi) = yi

ist, und yi+1 durch Anwendung der Methode auf yi mit Schrittweite h erzeugtwurde.

Die Methode besitzt Konsistenzordnung p ∈ N, falls für jede (unsere Gene-ralvoraussetzung erfüllende) rechte Seite f

lim suph→0


|yi+1 − y(xi + h)|hp+1

< ∞,

d.h. Konstanten C > 0, h0 > 0 existieren, so dass


|yi+1 − y(xi + h)| ≤ Chp+1 ∀0 < h < h0.

Eine Methode der Konsistenzordnung p macht in jedem Intervall [xi, xi+1]einen lokalen Fehler der Größenordnung hp+1. Da es n = b−a

hsolche Intervalle

gibt, erwarten wir dass der globale Fehler in der Größenordnung hp liegenwird. Der folgende Satz zeigt, dass dies tatsächlich der Fall ist.

Satz 1.16Sei y : [a, b] → R

d die Lösung des AWP

y′(x) = f(x, y(x)), y(x0) = y0 ∈ Rd

Wir betrachten die Anwendung einer Einschrittmethode mit äquidistanterDiskretisierung, d.h. mit Schrittweite h := b−a

n> 0,

a = x0 < x1 < . . . < xn = b, xi = x0 + ih (i = 1, . . . , n).

21


(a) Ist die Methode konsistent, so gilt

maxi=1,...,n

‖yi − y(xi)‖ → 0 für n → ∞.

(b) Besitzt die Methode Konsistenzordnung p, so gilt

maxi=1,...,n

‖yi − y(xi)‖ ≤ eL(b−a) − 1

LChp ∀0 < h < h0.

wobei C, h0 > 0 die Konstanten aus der Definition der Konsistenzord-nung und L die Lipschitz-Konstante aus der Generalvoraussetzung ist.

Beweis: Wir beginnen mit (b). Für die Schrittweite gelte h = b−an

< h0. Dawir mit dem korrekten Startwert y0 = y(x0) beginnen, gilt für den Fehlernach dem ersten Schritt

‖y1 − y(x1)‖ ≤ Chp+1.

Im nächsten Schritt, bei dem wir y2 aus y1 berechnen, gibt es zwei Fehler-quellen:

(i) Die Berechnung von y2 aus y1 mit dem Einschrittverfahren entsprichtder Anwendung des Verfahrens auf das (wegen unserer Generalvoraus-setzung eindeutig lösbare) AWP

z′(x) = f(x, z(x)), z(x1) = y1 ∈ Rd (1.3)

Dabei macht das Verfahren den Fehler

‖z(x2)− y2‖ ≤ Chp+1.

(ii) Da y1 nur eine Approximation an y(x1) ist, stimmt die Lösung z(x) von(1.3) nicht mit y(x) überein. Nach Satz 1.12 gilt aber

‖y(x2)− z(x2)‖ ≤ eL(x2−x1) ‖y(x1)− y1‖ .

Insgesamt erhalten wir also

‖y2 − y(x2)‖ ≤ ‖y2 − z(x2)‖+ ‖z(x2)− y(x2)‖≤ Chp+1 + eLh ‖y1 − y(x1)‖≤ (1 + eLh)Chp+1.

22


Mit trivialer Induktion erhalten wir für alle i = 1, . . . , n:

‖yi − y(xi)‖ ≤ Chp+1 + eLh ‖|yi−1 − y(xi−1)|‖≤ Chp+1 + eLh

(Chp+1 + eLh ‖yi−2 − y(xi−2)‖

)

≤ . . . ≤(1 + eLh + e2Lh + e(i−1)Lh

)Chp+1

≤n−1∑

j=0

(eLh)jChp+1 =

(eLh)n − 1

eLh − 1Chp+1.

und mit eLh ≥ 1 + Lh folgt

maxi=1,...,n

‖yi − y(xi)‖ ≤ enhL − 1

LhChp+1 =

eL(b−a) − 1

LChp.

Der Beweis von (a) geht analog.

BemerkungIm Folgenden verwenden wir im Zusammenhang mit Anfangswertproblemendie Landau-Notation O(hp), o(h), etc., mit der Konvention, dass die darinvorkommenden Konstanten von der rechten Seite f nicht jedoch von demAnfangswert x0 ∈ [a, b], y0 ∈ Rd abhängen dürfen. Mit dieser Konvention istein Einzelschrittverfahren

• konsistent, falls aus yi = y(xi) folgt, dass

y(xi+1)− yi+1 = o(h).

• konsistent mit Ordnung p, falls aus yi = y(xi) folgt, dass

y(xi+1)− yi+1 = O(hp+1).

Beispiele 1.17(a) Explizites Euler-Verfahren: Für y(xi) = yi erhalten wir durch Tay-

lorentwicklung

‖y(xi+1)− y(xi) + hy′(xi)‖ ≤ h2

2max

ξ∈[xi,xi+1]‖y′′(ξ)‖ ≤ h2

2C2

mit einer (aufgrund unserer Generalvoraussetzung nur von der rechtenSeite f abhängigen) Konstante C2. Es ist also mit obiger Konvention

y(xi+1) = y(xi) + hy′(xi) +O(h2)

= y(xi) + hf(xi, yi) +O(h2) = yi+1 +O(h2)

und das explizite Euler-Verfahren besitzt somit Konsistenzordnung 1.

23


(b) Implizites Euler-Verfahren: Für y(xi) = yi erhalten wir wiederumdurch Taylorentwicklung

y(xi) = y(xi+1)− hy′(xi+1) +O(h2)

= y(xi+1)− hf(xi+1, y(xi+1)) +O(h2)

Zusammen mit yi+1 = yi + hf(xi+1, yi+1) erhalten wir

‖yi+1 − y(xi+1)‖ = ‖yi + hf(xi+1, yi+1)− y(xi+1)‖= h ‖f(xi+1, yi+1)− f(xi+1, y(xi+1))‖+O(h2)

≤ hL ‖yi+1 − y(xi+1)‖+O(h2).

Für hinreichend kleine h gilt also

‖yi+1 − y(xi+1)‖ =1

1− hLO(h2) = O(h2).

Das implizite Euler-Verfahren besitzt also Konsistenzordnung 1.

1.4.2 Runge Kutta Methoden

Wir betrachten nun einen allgemeinen Ansatz um Einschrittmethoden hoherKonsistenzordnung zu konstruieren. Im ersten Schritt (mit h = x1 − x0) sollgelten

y1 ≈ y(x1) = y0 +

∫ x1

x0

y′(t) dt = y0 +

∫ x1

x0

f(t, y(t)) dt

Durch Approximation des Integrals auf der rechten Seite durch eine Quadra-turformel erhalten wir

∫ x1

x0

f(t, y(t)) dt ≈ h

s∑

j=1

bjf(x0 + cjh, ηj).

Die Quadraturformel sollte zumindest konstante Funktion exakt integrieren,deshalb fordern wir

s∑

j=1

bj = 1.

Die cj heißen Knoten, bj Gewichte und s heißt Stufenzahl des Verfahrens.

Dieser Ansatz verallgemeinert die Ideen aus Abschnitt 1.3.4, indem nun einegewichtetes Mittel aus s unterschiedlichen Steigungnen im Richtungsfeld anden Punkten (x0 + cjh, ηj), j = 1, . . . , s verwendet werden. Für das expliziteEulerverfahren ist s = 1, c1 = 0, and η1 = y0.

24


bj und cj sollten aus den Knoten und Gewichten eines möglichst guten Qua-draturverfahren bestimmt werden. Zur Wahl der ηj fordern wir dass

ηj ≈ y(x0 + cjh) = y0 +

∫ x0+cjh

x0

y′(t) dt = y0 +

∫ x0+cjh

x0

f(t, y(t)) dt.

Wir wenden wiederum ein Quadraturverfahren für die Integrale auf der rech-ten Seite an und verwenden dabei die gleichen Quadraturpunkte wie für daserste Integral, d.h. für j = 1, . . . , s verwenden wir

∫ x0+cjh

x0

f(t, y(t)) dt. ≈ h

s∑

l=1

ajlf(x0 + clh, ηl).

Wiederum sollten die Quadraturformeln zumindest konstante Funktionen ex-akt integrieren, deshalb fordern wir

s∑

l=1

ajl = cj.

So erhalten wir die allgemeinen Runge-Kutta Methoden:

Gegeben ajl, bj , cj ∈ R, j = 1, . . . , s, l = 1, . . . , s mit

s∑

j=1

bj = 1 unds∑

l=1

ajl = cj.

• Bestimme ηj ∈ Rd, j = 1, . . . , s aus

ηj = yi + hi

s∑

l=1

ajlf(xi + clhi, ηl), j = 1, . . . , s. (1.4)

• Setze

yi+1 := yi + hi

s∑

j=1

bjf(xi + cjhi, ηj).

Runge-Kutta Methoden können explizit oder implizit sein. Ist

ajl = 0 für l ≥ j

dann ist η1 = y1, η2 kann aus η1 berechnet werden, usw. (explizite Runge-Kutta Methoden). Ansonsten ist (1.4) ein System aus sd Gleichungen für die

25


sd unbekannten Einträge der d-dimensionalen Vektoren ηj ∈ Rd, j = 1, . . . , s(implizite Runge-Kutta Methoden).

Eine äquivalente Formulierung erhalten wir durch

kj := f(xi + cjhi, ηj)

Gegeben ajl, bj , cj ∈ R, j = 1, . . . , s, l = 1, . . . , s mit

s∑

j=1

bj = 1 unds∑

l=1

ajl = cj .

• Bestimme kj ∈ Rd, j = 1, . . . , s aus

kj = f(xi + cjh, yi + hi

s∑

l=1

ajlkl), j = 1, . . . , s.

• Setze

yi+1 := yi + hi

s∑

j=1

bjkj.

Die Koeffizienten A = (ajl) ∈ Rs×s, b = (bj) ∈ Rs und c = (cj) ∈ Rs einerRunge-Kutta Methode lassen sich im sogenannten Butcher Tableau zusam-menfassen:

c AbT

c1 a11 a12 . . . a1sc2 a21 a22 . . . a2s...

......

...cs as1 as2 . . . ass

b1 b2 . . . bs

Mit dieser Notation erhalten wir für das explizite und implizite Euler Ver-fahren

0 01

1 11

1.4.3 Wohldefiniertheit impliziter Methoden

Die Vorteile impliziter Runge-Kutta Methoden werden wir erst in Abschnitt1.5 kennenlernen. Wir zeigen aber an dieser Stelle schon, dass das nicht-lineare Gleichungssystem (1.4) für hinreichend kleine Schrittweiten eindeutiglösbar ist.

26


Satz 1.18Zu jeder rechten Seite f (die die Generalvoraussetzung erfüllt) und jedemRunge-Kutta Verfahren (A, b, c) existiert eine Schrittweite h0 > 0, so dassfür jedes x ∈ [a, b], y ∈ Rd und 0 < h ≤ h0 das Gleichungsystem für die ηj

ηj = y + hs∑

l=1

ajlf(x+ clh, ηl), j = 1, . . . , s.

eindeutig lösbar ist.

Beweis: Wir schreiben das Gleichungssystem als Fixpunktgleichung

η = Φ(η)

mit

η :=

η1...ηs

∈ Rds, Φ(η) :=

Φ1(η)...

Φs(η)

∈ Rds

und

Φj(η) := y + h

s∑

l=1

ajlf(x+ clh, ηl) ∈ Rd.

Es ist

Φ′(η) =

Φ′1(η)...

Φ′s(η)

∈ Rds×ds

Φ′j(η) =

(∂Φj

∂η1. . .

∂Φj

∂ηs

)

∈ Rd×ds

∂Φj

∂ηl= hajlfy(x+ clh, ηl) ∈ Rd×d.

Aufgrund unserer Generalvoraussetung und der Äquivalenz aller Normen aufdem R

d×d existiert ein C > 0, so dass für alle j, l = 1, . . . , d jeder Eintrag derMatrix ∂Φj

∂ηldurch Ch beschränkt ist. Damit ist jeder Eintrag von Φ′(η) durch

Ch beschränkt und (wegen der Äquivalenz aller Normen auf dem R

ds×ds)existiert ein C ′ > 0 mit

‖Φ′(η)‖ ≤ Ch.

Für hinreichend kleine h ist Φ also eine Kontraktion und die Behauptung wiein der Numerik I aus dem Banachschen Fixpunktsatz.

27


1.4.4 Runge-Kutta Ordnungsbedingungen

Im letzten Abschnitt haben wir gesehen, dass jede Wahl der Runge-KuttaKoeffizienten (A, b, c) auf lösbare implizite (oder sogar explizite) Gleichungführt, das Verfahren also für jede Wahl der Koeffizienten durchführbar ist.Jetzt wenden wir uns der Frage zu, wie die Koeffizienten gewählt werdenmüssen, um ein Verfahren möglichst hoher Ordnung zu erhalten.

Satz 1.19Seien A = (aij)i,j=1,...,s, b = (bi)i=1,...,s und c = (ci)i=1,...,s die Koeffizienteneines Runge-Kutta-Verfahrens.

(a) Auss∑

j=1

bj = 1 unds∑

k=1

ajk = cj

folgt dass das Verfahren (mindestens) Konsistenzordnung 1 besitzt.

(b) Das Verfahren hat genau dann mindestens Konsistenzordnung 2, wennzusätzlich gilt, dass

s∑

j=1

bjcj =1

2.

(c) Das Verfahren hat genau dann mindestens Konsistenzordnung 3, wennzusätzlich gilt, dass

s∑

j=1

bjc2j =

1

3und

s∑

j=1

bj

s∑

k=1

ajkck =1

6.

Beweis: Wegen Übungsaufgabe 4.2 genügt es, die Behauptung für autonomeDifferentialgleichungen

y′ = f(y), f : Rd → R

d

zu beweisen.

Sei y eine Lösung der DGL, xi ∈ [a, b], yi = y(xi), h = xi+1 − xi und yi+1 diedurch das Runge-Kutta Verfahren erzielte Näherung. Nach Übungsaufgabe3.2 gilt

y(xi+1) = yi + hf + 1/2h2f ′f + 1/6h3(fTf ′′f + (f ′)2f) +O(h4)

wobei wir das Argument (yi) bei f und seinen Ableitungen zur Vereinfachungder Schreibweise weglassen.

28


Genauso entwickeln wir die Näherung yi+1. Dabei verwenden wir immer wie-der

ηj = yi + h

s∑

k=1

ajkf(ηk),

s∑

k=1

ajk = cj (1.5)

Zuerst erhalten wir damit

ηj = yi +O(h) =⇒ f(ηj) = f(yi) +O(h) ∀j = 1, . . . , s.

Nochmalige Anwendung von (1.5) führt zu

ηj − yi = hs∑

k=1

ajkf(ηk) = hs∑

k=1

ajkf(yi) +O(h2) = hcjf(yi) +O(h2).

Ein weiteres Mal verwenden wir (1.5) und kombinieren es mit

f(ηk) = f(yi) + f ′(yi)(ηk − yi) +O(h2).

So erhalten wir

ηj = yi + hs∑

k=1

ajkf(ηk) = yi + hs∑

k=1

ajk(f(yi) + f ′(yi)(ηk − yi) +O(h2))

= yi + hs∑

k=1

ajk(f(yi) + f ′(yi)(hckf(yi) +O(h2)

)+O(h2))

= yi + hcjf(yi) + h2f ′(yi)f(yi)

s∑

k=1

ajkck +O(h3)

Mit

yi+1 := yi + h

s∑

j=1

bjf(ηj),

s∑

j=1

bj = 1

29


folgt (wobei wir das Argument (yi) bei f und seinen Ableitungen weglassen)

yi+1 = yi + h

s∑

j=1

bjf(ηj)

= yi + hs∑

j=1

bj

(

f + f ′(ηj − yi) +1

2(ηj − yi)

Tf ′′(ηj − yi) +O(h3)

)

= yi + hf + hf ′s∑

j=1

bj(ηj − yi) +1

2h

s∑

j=1

bj(ηj − yi)Tf ′′(ηj − yi) +O(h4)

= yi + hf + hf ′s∑

j=1

bj

(

hcjf + h2f ′fs∑

k=1

ajkck +O(h3)

)

+1

2h

s∑

j=1

bj(hcjf +O(h2)

)Tf ′′(hcjf +O(h2)

)+O(h4)

= yi + hf + h2f ′fs∑

j=1

bjcj + h3(f ′)2fs∑

j=1

bj

s∑

k=1

ajkck

+1

2h3fTf ′′f

s∑

j=1

bjc2j +O(h4)

Die Behauptung folgt nun aus dem Vergleich der Entwicklungen von y(xi+1)und yi+1.

Bemerkung 1.20(a) Mit Satz 1.19 lässt sich die Ordnung der Verfahren in 1.3.4 bestimmen.

(b) Mit einem systematischeren symbolischen Ansatz können Methoden be-liebig hoher Ordnung konstruiert werden.

(c) Man kann zeigen, dass die Ordnung eines s-stufigen Runge-Kutta Verfah-rens höchstens 2s ist. Explizite Runge-Kutta Verfahren können höchstensOrdnung s haben (siehe Abschnitt 1.5.5).

(d) Die in der Praxis wohl am häufigsten verwendete explizite Runge-KuttaMethode ist eine Methode 5. Ordnung von Dormand und Prince, die

30

1.5. NUMERIK STEIFER DIFFERENTIALGLEICHUNGEN

durch folgendes Tableau gegeben ist:

015

15

310

340

940

45

4445

−5615

329

89

193726561

−253602187

644486561

−212729

1 90173168

−35533

467325247

49176

− 510318656

35384

0 5001113

125192

−21876784

1184

Dieses Verfahren ist (in Kombination mit einer zur adaptiven Schritt-weitensteuerung verwendeten Methode 4. Ordnung) unter dem Namendopri5 oder ode45 in vielen Programmpaketen der Standardlöser fürAnfangswertprobleme.

1.5 Numerik steifer Differentialgleichungen

1.5.1 Steife Differentialgleichungen

Bei dem Pendel aus Übungsaufgabe 3.4 waren implizite Verfahren (trotzgleicher Konsistenzordnung) den expliziten weit überlegen. Differentialglei-chungen, in denen dieser Effekt auftritt, werden steif genannt. Steif ist dabeikein mathematisch präzise definierter Begriff, sondern wird anschaulich fürsolche Differentialgleichungen verwendet, bei denen naheliegende Standard-verfahren (z.B. explizite Runge-Kutta-Verfahren) keine (oder nur für extremkleine Schrittweiten) zufriedenstellenden Ergebnisse liefern.

Betrachten wir das einfache Beispiel

y′(x) = λy, y(0) = 1, λ < 0.

Offenbar ist die Lösung y(x) = eλx. Aufgrund der Annahme λ < 0 konvergiertdie Lösung mit exponentiell Geschwindigkeit gegen Null.

Durch Anwendung des expliziten Euler-Verfahrens mit Schrittweite h auf

31


dieses AWP erhalten wir

y1 = y0 + hλy0 = (1 + hλ)

y2 = y1 + hλy1 = (1 + hλ)y1 = (1 + hλ)2

...

yn = (1 + hλ)n

Da λ < 0 folgt für n → ∞

yn > 0 und yn → 0 für 1 + hλ ≥ 0yn alterniert im Vorzeichen, yn → 0 für 0 > 1 + hλ > −1yn alterniert zwischen +1 und −1 für 1 + hλ = −1yn alterniert im Vorzeichen, |yn| → ∞ für 1 + hλ < −1

Nur für 1 + hλ ≥ −1 (d.h. h ≤ −2/λ) zeigen die Approximationen alsodas korrekte Langzeitverhalten und konvergieren gegen Null, und für h >−1/λ oszillieren die Approximationen. Für dieses AWP mit stark negativemλ liefert die explizite Euler Methode also nur für extrem kleine Schrittweitenbrauchbare Ergebnisse.

Betrachten wir zum Vergleich die implizite Euler-Methode:

y1 = y0 + hλy1 =⇒ y1 = (1− hλ)−1

y2 = y1 + hλy2 =⇒ y2 = (1− hλ)−1y1 = (1− hλ)−2

...

yn = (1− hλ)−n

Für jede Schrittweite h, ist 1 − hλ > 1. yn bleibt also stets positiv undkonvergiert gegen Null für n → ∞.

Implizites und explizites Euler-Verfahren besitzen die gleiche Konsistenzord-nung. Für h → 0 konvergieren sie gleich schnell gegen die wahre Lösung.Jedoch gibt es zwei Eigenschaften der wahren Lösung, Positivität und Ab-fallverhalten, die nur die Iterierten des impliziten Euler-Verfahren für alleSchrittweiten zeigen. Die Iterierten des expliziten Euler-Verfahrens habendiese Eigenschaften erst für extrem kleine Schrittweiten.

Das betrachtete AWP ist also steif in dem Sinne, dass die wahre Lösunggewisse Eigenschaften besitzt, die so wichtig sind, dass man nur solche nu-merischen Approximationen akzeptieren wird, die diese Eigenschaften auchbesitzen.

32


1.5.2 Die Testgleichung

Wir motivieren nun heuristisch, dass das im letzten Abschnitt beobachteteVerhalten auch in allgemeine Anfangswertprobleme wiederfinden lässt.

Betrachten wir das allgemeine AWP

y′(x) = f(x, y), y(x0) = y0 ∈ Rd

Gemäß Übungsaufgabe 4.2 können wir es in eine autonome Gleichung um-formen. Außerdem können wir durch Verschiebung annehmen, dass x0 = 0.

y′(x) = f(y), y(0) = y0 ∈ Rd

Für kleine x wird sich die Lösung nur wenig verändern. Wir erwarten also,dass sich y lokal wie die Lösung linearisierten Gleichung

y′(x) = f(y) ≈ f(y0) + f ′(y0)(y − y0), y(0) = y0 ∈ Rd

verhält.

Wir nehmen noch an, dass sich die Shifts f(y0) und y0 durch geeignete Trans-formationen eliminieren lassen. Lokal lässt sich das AWP dann durch dieLösung der Gleichung

y′(x) = My

mit einer Matrix M ∈ Rd×d approximieren. Ist M diagonalisierbar mit Ei-genwerten λ1, . . . , λd, dann ist dies äquivalent zu d skalaren Testgleichungen

y′j = λjyj, j = 1, . . . , d.

Die Eigenwerte λj werden im Allgemeinen komplex sein. Offenbar gelten aberalle Ergebnisse dieses Kapitels auch genauso für komplexwertige Gleichungen.

Insgesamt scheint es also erstrebenswert, Methoden zu konstruieren, die nichtnur möglichst schnell konvergieren, sondern auch qualitativ richtiges Verhal-ten zeigen für die komplexe skalare Testgleichung

y′ = λy, λ ∈ C.

Aufgrund der Linearität der Gleichung können wir dabei den Anfangswertauf y(0) = 1 setzen.

33


1.5.3 Die Stabilitätsfunktion

Nach Abschnitt 1.5.1 gilt für die Iterierten des expliziten und implizitenEuler-Verfahrens bei Anwendung auf die Testgleichung (mit λ < 0)

y(expl)i = (1 + hλ)iy0 = R(expl)(hλ)iy0,

y(impl)i = (1− hλ)−iy0 = R(impl)(hλ)iy0

wobeiR(expl)(ζ) := (1 + ζ), and R(impl)(ζ) = (1− ζ)−1

Offenbar gilt das auch für λ ∈ C. Wie gut die Verfahren für die Testgleichungfunktionieren, lässt sich also vollständig mit der Funktion R(ζ) beschreiben.Gleiches gilt für allgemeine Runge Kutta Methoden.

Definition 1.21Seien

A = (aij)i,j=1,...,s ∈ Rs×s, b = (bi)i=1,...,s ∈ Rs, und c = (ci)i=1,...,s ∈ Rs

die Koeffizienten einer Runge-Kutta-Methode. Sei 1 := (1, . . . , 1)T ∈ Rs undI ∈ Rs×s sei die Einheitsmatrix. Für ζ ∈ C definieren wir

R(ζ) := 1 + ζbT (I − ζA)−11 ∈ C

falls I − ζA invertierbar ist. Ansonsten schreiben wir formal R(ζ) = ∞.(Offenbar ist dies genau dann der Fall, wenn 1

ζein Eigenwert von A ist, also

für höchstens s komplexe Zahlen).

Satz 1.22Betrachte die Anwendung des Runge-Kutta Verfahrens mit Koeffizienten A ∈R

s×s, b, c ∈ Rs auf die Testgleichung

y′(x) = λy(x), y0 = 1.

mit λ ∈ C und Schrittweite h > 0.

Ist die Matrix I − hλA ∈ Cs×s invertierbar, so ist die Runge Kutta Methodeanwendbar (d.h. die möglicherweise impliziten Gleichungen lösbar) und ihreIterierten sind gegeben durch

yi = (R(hλ))i.

34


Beweis: Anwendung der Runge-Kutta Methode liefert das (möglicherweiseimplizite) Gleichungssystem

ηj = yi + hs∑

l=1

ajlληl, j = 1, . . . , s.

Mit η := (η1, . . . , ηs) ∈ Cs ist das äquivalent zu

η = yi1+ hλAη ⇐⇒ (I − hλA)η = yi1

Ist I−hλA invertierbar, so existiert eine eindeutige Lösung η und wir erhalten

yi := yi−1 + hs∑

j=1

bjληj = yi−1 + hλbT η

= yi−1 + hλbT (I − hλA)−1yi−11 = (1 + hλbT (I − hλA)−11)yi−1

= (1 + ζbT (I − ζA)−11)iy0 = (R(ζ))i, ζ := hλ

Beispiel 1.23(a) Die Stabilitätsfunktion des expliziten Eulerverfahrens ist R(ζ) := 1 + ζ.

(b) Die Stabilitätsfunktion des impliziten Eulerverfahrens ist R(ζ) = (1 −ζ)−1.

(c) Das Butcher Tableau für die implizite Mittelpunktsformel aus Abschnitt1.3.4 ist (vgl. Übungsaufgabe 4.1)

1/2 1/21

Die Stabilitätsfunktion ist also

R(ζ) = 1 + ζbT (I − ζA)−11 = 1 + ζ1(1− ζ 1/2)−11 =

1 + ζ/2

1− ζ/2.

1.5.4 Stabilität

Die exakte Lösung der Testgleichung y′ = λy, y(0) = 1 ist

y(x) = eλx

35


Es gilt also

|y(x)|

→ ∞ für x → ∞ wenn Re(λ) > 0→ 0 für x → ∞ wenn Re(λ) < 0= 1 für alle x ≥ 0 wenn Re(λ) = 0

und außerdem ist, für alle x > 0,

|y(x)| → 0, wenn Re(λ) → −∞.

Dies motiviert die folgende Definition.

Definition und Satz 1.24yi ≈ y(hi) seien die Approximationen einer Einschrittmethode auf die Test-gleichung y′ = λy, y(0) = 1. Die Methode heißt

• A-stabil falls für Re(λ) ≤ 0 stets gilt, dass

|yi+1| ≤ |yi| für alle i und alle Schrittweiten h

• Isometrie-erhaltend wenn für Re(λ) = 0 stets gilt, dass

|yi+1| = |yi| für alle i und alle Schrittweiten h

• L-stabil, wenn sie A-stabil ist und (für alle h > 0)

|y1| → 0 für |λ| → ∞.

Eine Runge-Kutta Methode ist genau dann

• A-stabil, wenn |R(ζ)| ≤ 1 für all ζ ∈ C mit Re(ζ) ≤ 0,

• Isometrie-erhaltend, wenn |R(ζ)| = 1 für alle ζ ∈ C mit Re(ζ) = 0,

• L-stabil, wenn A-stabil und |R(ζ)| → 0 für |ζ | → ∞.

Beweis: Die Äquivalenzen folgen aus yi = (R(hλ))i.

Man kann zeigen, dass die Stabilitätsfunktion eines Runge-Kutta Verfah-rens stets eine rationale Funktion ist, so dass das (bei der Definition derL-Stabilität) das Verhalten für |ζ | → ∞ mit dem für Re(ζ) → −∞ überein-stimmt.

36


Beispiel 1.25(a) Explizites Euler-Verfahren:

|R(i)| = |1 + i| =√2 > 1.

Das explizite Euler-Verfahren ist also weder A-stabil noch Isometrie-erhaltend.

(b) Implizites Euler-Verfahren:

Für alle ζ ∈ C mit Re(ζ) ≤ 0 ist

|R(ζ)| = 1

|1− ζ | =1

(1− Re(ζ))2 + Im(ζ)2≤ 1.

Das implizite Euler-Verfahren ist also A-stabil. Außerdem ist |R(ζ)| → 0für |ζ | → ∞, das Verfahren ist also auch L-stabil.

Es ist jedoch R(i) = 1/|1 − i| = 1/√2 < 1, das Verfahren ist also nicht

Isometrie-erhaltende.

(c) Die implizite Mittelpunktsformel ist A-stabil (jedoch nicht L-stabil) undIsometrie-erhaltend (siehe Übungsaufgabe ??).

Betrachten wir noch einmal die Testgleichung mit Re(λ) < 0. A-Stabilität be-deutet, dass die Approximationen das korrekte qualitative Verhalten |yi+1| ≤|yi| für jede Schrittweite h > 0 zeigen. Auch wenn eine Methode nicht A-stabil ist, kann sie dennoch dieses korrekte Verhalten zeigen, wenn nur dieSchrittweite klein genug gewählt ist (so dass |R(hλ)| ≤ 1). Dies motiviert diefolgende Definition.

Definition 1.26Zu einer Runge-Kutta Methode mit Stabilitätsfunktion R(ζ) definieren wirdas Stabilitätsgebiet durch

S := ζ ∈ C : |R(ζ)| ≤ 1 ⊆ C.

Beispielsweise besteht das Stabilitätsgebiet des expliziten Euler-Verfahrensaus allen ζ ∈ C mit

1 ≥ |R(ζ)|2 = |1 + ζ |2 = (1 + Re(ζ))2 + Im(ζ)2,

d.h. dem abgeschlossenen Kreis mit Radius 1 um z = −1 in der komplexenEbene.

37


1.5.5 Nachteile expliziter Verfahren

In unseren Beispielen waren nur implizite Verfahren A-stabil oder Isometrie-erhaltend. Tatsächlich gibt es expliziten Verfahren diesen Eigenschaften, wiewir in diesem Abschnitt zeigen.

Satz 1.27Die Stabilitätsfunktion einer expliziten Runge-Kutta Methode mit s Stufenist ein Polynom der Ordnung s.

Beweis: Seien A ∈ R

s×s, b ∈ R

s, c ∈ R

s die Koeffizienten der Methode.Da die Methode explizit ist, ist A eine strikte untere Dreiecksmatrix. Manzeigt leicht, dass in höheren Potenzen von A immer mehr Diagonalen durchNullen aufgefüllt werden, und schließlich As = 0 gilt:

A =

0 0 0 0 . . . 0∗ 0 0 0 . . . 0∗ ∗ 0 0 . . . 0∗ ∗ ∗ 0 . . . 0...

......

.... . .

...∗ ∗ ∗ ∗ . . . 0

, A2 =

0 0 0 0 . . . 00 0 0 0 . . . 0∗ 0 0 0 . . . 0∗ ∗ 0 0 . . . 0...

......

.... . .

...∗ ∗ ∗ ∗ . . . 0

,

A3 =

0 0 0 0 . . . 00 0 0 0 . . . 00 0 0 0 . . . 0∗ 0 0 0 . . . 0...

......

.... . .

...∗ ∗ ∗ ∗ . . . 0

, As =

0 0 0 0 . . . 00 0 0 0 . . . 00 0 0 0 . . . 00 0 0 0 . . . 0...

......

.... . .

...0 0 0 0 . . . 0

Aus As = 0 folgt dass

(I − ζA)(I + ζA+ . . . ζs−1As−1) = I

also (I − ζA)−1 = I + ζA+ . . . ζs−1As−1.

R(ζ) := 1 + ζbT (I − ζA)−11 ist also ein Polynom der Ordnung s.

Satz 1.28Explizites Runge-Kutta Methode sind weder A-stabil noch Isometrie-erhaltend.

Beweis: Für jedes Polynom R(ζ) gilt |R(ζ)| → ∞ für |ζ | → ∞.

Außerdem erhalten wir noch die schon in Bemerkung 1.20 angesprocheneHöchstgrenze in der Ordnung expliziter Verfahren:

38

1.6. LINEAR IMPLIZITE METHODEN

Satz 1.29Die Konsistenzordnung einer expliziten Runge-Kutta Methode mit s Stufenist höchstens s.

Beweis: Nach Satz 1.27 ist die Stabilitätsfunktion eine Polynom der Ord-nung s:

R(ζ) = r0 + r1ζ + . . .+ rsζs, r0, . . . , rs ∈ R.

Betrachte die Anwendung der Methode auf die Testgleichung mit λ = 1, alsoy′ = y, y(0) = 1:

y1 = R(h) = r0 + r1h + . . .+ rshs.

y(x1) = eh = 1 + h +1

2h2 + . . .+

1

s!hs +

1

(s+ 1)!hs+1 +O(hs+2)

Höchstens die ersten s Terme der Entwicklungen können übereinstimmen, sodass der lokale Fehler einer expliziten Methode höchstens O(hs+1), also dieOrdnung höchstens s sein kann.

1.6 Linear implizite Methoden

Wir haben gesehen, dass steife Differentialgleichungen implizite Methodenerfordern. Im Allgemeinen erfordert die Anwendung eines impliziten Runge-Kutta Verfahrens aber die Lösung von s gekoppelten d-dimensionalen nicht-linearen Gleichungen

kj = f(xi + cjh, yi + h

s∑

l=1

ajlkl), j = 1, . . . , s,

nach den sd unbekannten Einträgen der kj, j = 1, . . . , s. Ziel dieses Ab-schnitts ist die Herleitung von einfacheren und weniger Rechenaufwand er-fordernden, aber dennoch stabilen Methoden.

Wir beschränken uns dabei auf autonome AWP

y′ = f(y), y(x0) = y0

(nach Aufgabe 4.2 kann jedes AWP in diese Form gebracht werden).

Die erste Vereinfachung ist, dass wir eine Runge-Kutta Methode verwenden,für die A eine linke untere Dreiecksmatrix ist, also ajl = 0 für l > j. Dannkönnen die Gleichungen für die kj,

kj = f(yi + h

j−1∑

l=1

ajlkl + ajjhkj), j = 1, . . . , s,

39


beginnend mit k1 eine nach der anderen gelöst werden. Statt eines sd-dimensionalennicht-linearen Gleichungssystems müssen wir so nur s mal ein d-dimensionalesnicht-lineares Gleichungssystem lösen. Diese bringen wir auf Nullstellenform,also gegeben k1, . . . , kj−1 ist kj so zu bestimmen, dass

0 = Fj(kj) := kj − f(yi + h

j−1∑

l=1

ajlkl + ajjhkj).

Anwendung des Newton-Verfahrens ergibt ausgehend von einer Startnähe-rung k

(0)j die Iterationen

k(n+1)j := k

(n)j − F ′

j(k(n)j )−1Fj(k

(n)j ),

wobei

F ′j(kj) = I − f ′(yi + h

j−1∑

l=1

ajlkl + ajjhkj)ajjh.

Als weitere Vereinfachung ersetzen wir für alle j die wahre Jacobi-MatrixF ′j(kj) durch

F ′j(kj) ≈ I − ajjhJ, J := f ′(yi).

Außerdem führen wir nur einen einzelnen Newton-Schritt durch, d.h. für allej = 1, . . . , s setzen wir

kj := k(0)j − (I − ajjhJ)

−1Fj(k(0)j )

= k(0)j − (I − ajjhJ)

−1

(

k(0)j − f(yi + h

j−1∑

l=1

ajlkl + ajjhk(0)j )

)

= (I − ajjhJ)−1

(

f(yi + h

j−1∑

l=1

ajlkl + ajjhk(0)j )− ajjhJk

(0)j

)

Es bleibt noch die Wahl des Startwerte k(0)j zu klären. Hierzu verwenden wir

eine lineare Kombination der bereits berechneten kl, l = 1, . . . , j − 1:

k(0)j :=

j−1∑

l=1

djl/ajjkl

mit noch zu bestimmenden Koeffizienten djl. Insgesamt erhalten wir so dielinear impliziten (auch: Rosenbrock -) Runge Kutta Methoden.

40


Gegeben ajl, djl, bj , cj ∈ R, j = 1, . . . , s, l = 1, . . . , s.

• Setze J := f ′(yi) und bestimme kj, j = 1, . . . , s nacheinanderaus

kj := (I − ajjhJ)−1

(

f(yi + h

j−1∑

l=1

(ajl + djl)kl)− hJ

j−1∑

l=1

djlkl

)

• Setze

yi+1 := yi + h

s∑

j=1

bjkj .

Bemerkung 1.30I − ajjhJ ist offenbar invertierbar für 0 < h < 1

|ajj | ‖J‖.

Satz 1.31Seien (A, b, c) die Koeffizienten einer Runge-Kutta Methode, wobei A eine lin-ke untere Dreiecksmatrix und ajj 6= 0 sei. Dann hat die dazugehörige linearimplizite Runge-Kutta Methode im folgenden Sinne die gleichen Stabilitäts-eigenschaften wie die ursprüngliche Methode:

Ist R(ζ) die Stabilitätsfunktion der ursprünglichen Methode, dann ergibt sichfür jede Wahl der djl bei Anwendung der linear impliziten Methode auf dieTestgleichung

y′ = λy, y(0) = 1

die Approximationenyi = R(hλ)i,

wenn I − hλA invertierbar ist (also 1hλ

6= ajj für alle j).

Beweis: Wir wenden die linear implizite Methode auf die Testgleichung an

y′ = λy =: f(y), y(0) = 1.

Für alle y ist J = f ′(y) = λ und damit

kj := (1− ajjhλ)−1

(

λ(yi + h

j−1∑

l=1

(ajl + djl)kl)− hλ

j−1∑

l=1

djlkl

)

= (1− ajjhλ)−1

(

λyi + hλ

j−1∑

l=1

ajlkl

)

.

41


Mit k := (k1, . . . , ks)T ist das äquivalent zu

1− a11hλ 0 . . . 0−a21hλ 1− a22hλ . . . 0

......

. . .−as1hλ −as2hλ . . . 1− asshλ

k1k2...ks

=

λyiλyi...

λyi

.

Wenn I − hλA invertierbar ist, dann ist also

k = λyi(I − hλA)−11

und damit

yi+1 = yi + hbTk = (1 + hλbT (I − hλA)−11)yi = R(hλ)yi.

Beispiel 1.32(a) Linear-implizites Euler-Verfahren

Für das implizite Euler-Verfahren

1 11

ist A eine linke untere Dreiecksmatrix und keine d-Koeffizienten nötig.Das dazugehörige linear-implizite Euler Verfahren lautet

yi+1 := yi + hk, mit k := (I − hf ′(yi))−1f(yi).

(b) Linear-implizites Mittelpunktsverfahren

Genauso erhalten wir das linear-implizite Mittelpunktsverfahren:

yi+1 := yi + hk, with k := (I − h/2f ′(yi))−1f(yi).

(c) ode23s

Das wohl am häufigsten verwendete linear-implizite Verfahren besteht ausder folgenden Kombination einer zweistufigen (y) und einer dreistufigen(y) Methode:

k1 := (I − ahJ)−1f(yi)

k2 := (I − ahJ)−1 (f(yi + h/2 k1)− ahJk1)

k3 := (I − ahJ)−1 (f(yi + hk2)− d31hJk1 − d32hJk2)

yi+1 := yi + hk2

yi+1 := yi +h

6(k1 + 4k2 + k3),

42


mit

J := f ′(yi), a :=1

2 +√2, d31 := −4 +

√2

2 +√2, d32 :=

6 +√2

2 +√2.

y und y werden wir in Übungsaufgabe ?? zur adaptiven Schrittweiten-steuerung kombiniert. Das Verfahren ist in Matlab unter dem Namenode23s eines der zur Lösung steifer DGL empfohlenen Verfahren.

Lemma 1.33Das linear implizite Euler Verfahren ist L-stabil, das linear implizite Mittel-punktsverfahren ist A-stabil und Isometrie-erhaltend.

Beweis: Dies folgt aus Satz 1.31 und den Stabilitätseigenschaften des im-pliziten Eulerverfahrens und des impliziten Mittelpunktsverfahrens.

BemerkungGemäß Satz 1.31 definieren wir die Stabilitätsfunktion eines linear impli-ziten Verfahrens durch die des zugrundeliegenden Runge-Kutta-Verfahrens.Entsprechend nennen wir (wie in Lemma 1.33 schon praktiziert) ein line-ar implizites Verfahren A-stabil, L-stabil oder Isometrie-erhaltend, wenn daszugrundeliegende Runge-Kutta-Verfahren diese Eigenschaften hat.

Eine linear implizite Method besitzt die gleichen Stabilitätseigenschaften wiedie ursprüngliche Methode aber (je nach Wahl der djl) kann sich die Kon-sistenzordnung unterscheiden. Wie in Satz 1.19, lassen sich Ordnungsbedin-gungen für die Koeffizienten von linear impliziten Verfahren herleiten. Wirzeigen nur exemplarisch am Beispiel ode23s die Berechnung der Ordnungeines linear impliziten Verfahrens.

Satz 1.34Die in Beispiel 1.32 beschriebene zweistufige Methode zur Berechnung von yin ode23s besitzt Konsistenzordnung 2.

Beweis: Für jedes k ∈ Rd ist

‖(I − ahJ)k‖ ≥ ‖k‖ − ah ‖J‖ ‖k‖

und J = f ′(yi) ist aufgrund unserer Generalvoraussetzung unabhängig vomAnfangswert xi, yi beschränkt.

Für hinreichend kleine h > 0 ist die Matrix I − ahJ also invertierbar und esgilt (mit unserer Konvention bzgl. der Landau-Notation aus Abschnitt 1.4.1)

∥∥(I − ahJ)−1

∥∥ ≤ 1

1− ah ‖J‖ =1

1 +O(h)= O(1).

43


Wir gehen nun wie im Beweis von Satz 1.19 vor. Nach Übungsaufabe 3.2 giltfür die Lösung von y′ = f(y), y(xi) = yi

y(xi+1) = yi + hf + 1/2h2f ′f +O(h3)

wobei wir wieder das Argument (yi) von f und f ′ weglassen.

Wir wollen dies mit

yi+1 = yi + hk2,

vergleichen und müssen dazu also k2 bis zur Ordnung O(h2) entwickeln. Dazubenötigen wie die Entwicklung von k1. Aus

k1 = (I − ahJ)−1f und∥∥(I − ahJ)−1

∥∥ = O(1)

folgt k1 = O(1). Wir verwenden die Definition von k1 nocheinmal und erhal-ten

k1 = f + ahJk1 = f +O(h).

Für k2 folgt zuerst

k2 = (I − ahJ)−1 (f(y0 + h/2 k1)− ahJk1) = O(1).

und dann durch nochmalige Anwendung der Definition von k2

k2 = f(y0 + h/2 k1)− ahJk1 + ahJk2

= f + h/2k1f′ +O(h2)− ahJk1 + ahJk2 = f +O(h).

Noch ein weiteres Mal verwenden wir die Definition von k2 und erhaltenzusammen mit k1 = f +O(h), dass

k2 = f + h/2 k1f′ +O(h2)− ahJk1 + ahJk2

= f + h/2 ff ′ − ahJf + ahJf +O(h2) = f + h/2ff ′ +O(h2).

Insgesamt ist also

yi+1 = yi + hk2 = yi + hf + h2/2 ff ′ +O(h3) = y(xi) +O(h3),

die Methode besitzt also Konsistenzordnung 2.

44

1.7. MEHRSCHRITTVERFAHREN

1.7 Mehrschrittverfahren

Wir beschreiben nun noch kurz die wesentliche Idee der Mehrschrittverfah-ren. Dabei beschränken wir uns in diesem Abschnitt auf äquidistant gewählteGitterpunkte

xi = x0 + ih, h > 0.

In einem m-Schritt Verfahren verwenden wir die letzten m Approximationen

yi−m+1 ≈ y(xi−m+1), . . . , yi ≈ y(xi)

zur Bestimmung der nächsten Approximation yi+1 ≈ y(xi+1). Für die Be-stimmung der dafür nötigen ersten Werte y1, . . . , ym−1 können dabei Einzel-schrittverfahrens oder Mehrschrittverahren mit weniger Schritten verwendetwerden

1.7.1 Adams-Bashforth Methoden

Zur Bestimmung von yi+1 ≈ y(xi+1) aus yi−m+1, . . . , yi verwenden wir zuerstwie bei der Herleitung der Runge Kutta Methoden

yi+1 − yi ≈ y(xi+1)− y(xi) =

∫ xi+1

xi

y′(x) dx =

∫ xi+1

xi

f(x, y(x)) dx

Die Funktionx 7→ f(x, y(x))

ist (zumindest näherungsweise) an den Stellen

fj := f(xj , yj) ≈ f(xj, y(xj)), j = i−m+ 1, . . . , i

bekannt. Es liegt daher nahe, die unbekannte Funktion x 7→ f(x, y(x)) durchihr Interpolationspolynom f(x, y(x)) ≈ p(x), p ∈ Πm−1 durch die Stützstellen(xj , fj), j = i − m + 1, . . . , i zu ersetzen. Mit Hilfe der (aus der Numerik Ibekannten) Lagrange-Grundpolynome

lk(x) =∏

l=i−m+1,...,il 6=j

x− xl

xk − xl, k = i−m+ 1, . . . , i

können wir das Interpolationspolynom schreiben als

p(x) =

i∑

k=i−m+1

fklk(x).

45


So erhalten wir

yi+1 − yi =

∫ xi+1

xi

p(t) dx =

i∑

k=i−m+1

fk

∫ xi+1

xi

lk(x) dx

= h

i∑

k=i−m+1

fk

∫ 1

0

lk(xi + th) dt

= h

i∑

k=i−m+1

fk

∫ 1

0

∏

l=i−m+1,...,il 6=k

xi + th− xl

xk − xldt

= h

i∑

k=i−m+1

fk

∫ 1

0

∏

l=i−m+1,...,il 6=k

i− l + t

k − ldt

Mit der Umnummerierung k = i−m+ j und l = i−m+ j′ können wir dasschreiben als

yi+1 − yi = hm∑

j=1

fi−m+j

∫ 1

0

∏

j′=1,...,mj′ 6=j

m− j′ + t

j − j′dt

︸︷︷︸

=:βj

= hm∑

j=1

βjfi−m+j ,

mit (von h und i unabhängigen!) Konstanten βj ∈ R.

Die so erhaltenen Methoden heißen explizite Adams Methoden oder Adams-Bashforth Methoden.

Beispiel 1.35Für den Spezialfall m = 1 ergibt sich die explizite Euler-Methode. Für m = 2ist

β1 :=

∫ 1

0

2− 2 + t

1− 2dt = −

∫ 1

0

t dt = −1

2

β2 :=

∫ 1

0

2− 1 + t

1− 2dt =

∫ 1

0

(t+ 1) dt =3

2,

also yi+1 := yi + h(32fi − 1

2fi−1).

1.7.2 Weitere auf Integration basierende Methoden

Analog lassen sich implizite Adams Methoden (die Adams-Moulton-Methoden)aufstellen, indem das Interpolationspolynom durch die Stützstellen (xj , fj)

46


für j = i−m+ 1, . . . , i+ 1, also inklusive der noch zu bestimmenden Stütz-stelle gewählt wird. Dies führt auf Formeln der Form

yi+1 = yi + hm+1∑

j=1

βjfi−m+j = yi + hm+1∑

j=1

βjf(xi−m+j , yi−m+j) (1.6)

Eine verbreitete Methode diese impliziten Gleichungen zu lösen, ist es zu-erst eine Näherung an yi+1 (und damit an fi+1) durch die explizite AdamsMethode zu bestimmen. Diese Näherung wird dann als Startwert für eineFixpunktiteration der Gleichung (1.6) verwendet (üblich sind ein oder zweiIterationsschritte). Dieses Vorgehen heißt Predictor-Corrector-Verfahren.

Das Integrationsintervall bei der Herleitung der Methoden könnte auch vorxi liegende Bereiche umfassen, z.B.

yi+1 − yi−1 ≈ y(xi+1)− y(xi−1) =

∫ xi+1

xi−1

y′(x) dx =

∫ xi+1

xi−1

f(x, y(x)) dx

Analag zum Adams-Verfahren können wir f durch sein Interpolationspo-lynom (mit oder ohne Verwendung der unbekannten Stützstelle xi+1, fi+1)annähern und erhalten (implizite bzw. explizite) bzw. Formeln der Form

yi+1 = yi−1 + hm+1∑

j=1

βjfi−m+j bzw. yi+1 := yi−1 + hm∑

j=1

βjfi−m+j .

Diese Formeln heißen Nyström-Methoden (die explizite Variante) oder Milne-Simpson (die implizite Variante).

1.7.3 Auf Differentiation basierende Methoden

Die bisher betrachteten Mehrschrittverfahren beruhten auf der Idee die Funk-tion x 7→ f(x, y(x)) = y′(x) durch ein Interpolationspolynom zu approximie-ren und dieses zu integrieren. Statt dessen können wir auch die (Approxima-tionen an die) Funktionswerte yi−m+1, . . . , yi+1 durch ein Polynom interpo-lieren. Wie bei der Herleitung der Adams-Bashforth Methode lässt sich dasInterpolationspolynom q ∈ Πm schreiben als

q(x) =i+1∑

k=i+1−m

yklk(x), lk(x) =∏

l=i+1−m,...,i+1l 6=k

x− xl

xk − xl

47


also

q(xi + th) =i+1∑

k=i+1−m

yk∏

l=i+1−m,...,i+1l 6=k

xi + th− xl

xk − xl

=

i+1∑

k=i+1−m

yk∏

l=i+1−m,...,i+1l 6=k

i− l + t

k − l

=m+1∑

j=1

yi−m+j

∏

j′=1,...,m+1j′ 6=j

m− j′ + t

j − j′

Wir können nun versuchen, den unbekannten Wert yi+1 so zu bestimmen,dass das Interpolationspolynom q im aktuellen Gitterpunkt xi die Differen-tialgleichung erfüllt, also

q′(xi) = f(xi, yi).

Wegen

q′(xi) =1

h

∂

∂tq(xi + th)|t=0

=1

h

m+1∑

j=1

yi−m+j

∂

∂t

∏

j′=1,...,m+1j′ 6=j

m− j′ + t

j − j′

∣∣∣t=0

︸︷︷︸

=:αj

führt dies auf Formeln der Form

m+1∑

j=1

αjyi−m+j = hf(xi, yi),

die sich (für αm+1 6= 0) explizit nach yi+1 auflösen lassen.

Genauso führt die Forderung, dass das Interpolationspolynom q im nächstenGitterpunkt xi+1 die Differentialgleichung erfüllt, mittels

q′(xi+1) =1

h

∂

∂tq(xi+1 + th)|t=0

=1

h

m+1∑

j=1

yi−m+j

∂

∂t

∏

j′=1,...,m+1j′ 6=j

m+ 1− j′ + t

j − j′

∣∣∣t=0

︸︷︷︸

=:αj

48


auf implizite Methoden der Formm+1∑

j=1

αjyi−m+j = hf(xi+1, yi+1).

Die so entstandenen impliziten Formeln heißen auch BDF-Methoden (Back-ward differentiation formula).

Beispiel 1.36Für die implizite BDF-Methode mit m = 1 ergibt sich

α1 =∂

∂t

1 + 1− 2 + t

1− 2

∣∣∣t=0

= −1,

α2 =∂

∂t

1 + 1− 1 + t

2− 1

∣∣∣t=0

= 1,

also−1yi + 1yi+1 = hf(xi+1, yi+1),

und damit gerade die implizite Euler-Methode.

1.7.4 Konvergenz linearer Mehrschrittmethoden

Alle bisher kennengelernten Mehrschrittmethoden können wir in der allge-meinen Form

m+1∑

j=1

αjyi−m+j = h

m+1∑

j=1

βjfi−m+j

mit Konstanten α1, . . . , αm+1, β1, . . . , βm+1 ∈ R schreiben.

Wir geben in dieser Vorlesung nur eine ganz kurze Zusammenfassung derTheorie dieser Methoden. Eine ausführlichere Darstellung findet sich z.B. in[HairerNorsettWanner].

Analog zu Einschrittmethoden wird die Konsistenzordnung eines Mehrschritt-verfahren durch Betrachtung des lokalen Fehlers ‖yi+1 − y(xi+1‖) definiert,der sich ergibt, wenn die Methode auf m exakte Werte

yi−m+1 = y(xi−m+1), . . . , yi = y(xi)

angewendet wird.

Satz 1.16 lässt sich jedoch nicht unmittelbar auf Mehrschrittverfahren über-tragen. Im Gegensatz zu Einschrittverfahren folgt aus der Konsistenz einesMehrschrittverfahrens nicht automatisch Konvergenz, sondern dies erforderteine zusätzliche Stabilitätseigenschaft. Die in Abschnitt 1.7.1 und 1.7.2 vor-gestellten Adam-Varianten erfüllen diese zusätzliche Eigenschaften, die BDF-Formeln jedoch nur bis m ≤ 6.

49


1.8 Randwertprobleme

1.8.1 Motivation: Diffusionsprozesse

Neben Anfangswertproblemen treten in der Praxis auch Randwertproblemefür gewöhnliche Differentialgleichungen auf. Die Theorie und Numerik dieserProbleme ist eng mit der für partielle Differentialgleichung verwandt, da (wiein der folgenden Motivation) Randwertprobleme für gewöhnliche Differenti-algleichungen oft als eindimensionale stationäre Spezialfälle von Randwert-problemen für PDGL auftreten. Die folgende Modellierung von Diffusions-prozessen folgt dem sehr lesenswerten Buch [FulfordBroadbridge].

Wir betrachten ein Rohr mit Querschnitt A, das von einer Lösung durchflos-sen wird. Wir bezeichnen mit

x: die Position innerhalb des Rohres, etwa x ∈ [0, 1]

C(x, t): die Konzentration des gelösten Stoffes am Ort x zur Zeit t

J(x, t): Flussdichte der Lösung, d.h. welche Masse des Stoffes einen Einheits-querschnitt pro Zeiteinheit durchquert.

Wir betrachten den Rohrabschnitt zwischen x und x+δx. Dabei nehmen wiran, dass δx so klein ist, dass die Konzentration in diesem Abschnitt räumlichkonstant ist. Die Gesamtmasse innerhalb des Abschnitts ist also

AδxC(x, t).

Nun nehmen wir an, dass δt so klein ist, dass der Fluss im Zeitabschnittzwischen t und t+ δt zeitlich konstant ist. Augrund des Flusses wird sich imbetrachteten Abschnitt des Rohres die Gesamtmasse in diesem Zeitabschnittändern um

J(x, t)Aδt− J(x+ δx, t)Aδt,

vgl. die in der Vorlesung gemalten Skizzen.

Wenn es keine anderen die Masse ändernden Phänomene gibt, so gilt also

AδxC(x, t + δt) = AδxC(x, t) + J(x, t)Aδt− J(x+ δx, t)Aδt

alsoC(x, t+ δt)− C(x, t)

δt= −J(x+ δx, t)− J(x, t)

δx

50

1.8. RANDWERTPROBLEME

und mit δx → 0, δt → 0 erhalten wir die Bilanzgleichung

∂C(x, t)

∂t= −∂J(x, t)

∂x.

Die einfachste Model für Diffusion ist Fick’s Gesetz, das besagt, dass dieFlussdichte proportional ist zum Konzentrationsgefälle

J(x, t) = −D(x)∂C(x, t)

∂x.

(D(x, t) heißt Diffusionskonstante).

So erhalten wir eine partielle Differentialgleichung, die sogenannte Diffusi-onsgleichung oder auch Wärmeleitungsgleichung

∂C(x, t)

∂t=

∂

∂x

(

D(x, t)∂C(x, t)

∂x

)

.

Konvektion und Absorption können ähnlich modelliert werden. Wenn dieFlüssigkeit sich mit der Geschwindigkeit v(x, t) bewegt, dann muss der Termv(x, t)C(x, t) zum Fluss addiert werden. Wenn pro Zeiteinheit und Raumein-heit die Masse M(x, t) hinzugegeben wird oder a(x, t)C(x, t) z.B. aufgrundeiner chemischen Reaktion verbraucht wird, so müssen diese Änderungen inder Massenbilanz berücksichtigt werden. Insgesamt erhalten wir

∂C

∂t(x, t) =

∂

∂x

(

D(x, t)∂

∂xC(x, t)

)

− ∂

∂x(v(x, t)C(x, t))−a(x, t)C(x, t)+M(x, t).

Es erscheint natürlich, dass diese partielle Differentialgleichung Anfangsbe-dingungen C(x, 0) für alle x ∈ (0, 1) und Randbedingungen für x = 0 undx = 1 benötigt. Als Randbedingungen können wir z.B. die KonzentrationC(0, t) and C(1, t) für alle t > 0 (Dirichlet Randbedingungen) oder den Fluss−D(0)∂C(0,t)

∂xund −D(1)∂C(1,t)

∂x(Neumann Randbedingungen) vorschreiben.

Sind alle Koeffizienten der Gleichung von der Zeit unabhängig, so stellt sichoft mit der Zeit ein Gleichgewichtszustand ein, d.h. die Konzentration ändertsich nicht mehr. Für diesen muss also gelten

∂

∂x

(

D(x)∂

∂xC(x)

)

− ∂

∂x(v(x)C(x))− a(x)C(x) +M(x) = 0

Dies ist wieder eine gewöhnliche Differentialgleichung, für die wir jedoch(Dirichlet- oder Neumann-)Randwerte anstelle von Anfangswerten kennen.

51


1.8.2 Differenzenverfahren

Motiviert durch den letzten Abschnitt betrachten wir nun die leicht verein-fachte Diffusionsgleichung

L[u] := −u′′(x) + b(x)u′(x) + c(x)u(x) = f(x) x ∈ (0, 1)

und zwar zuerst mit homogenen Dirichlet-Randbedingungen u(0) = u(1) = 0.

Es ist naheliegend, dass Randwertproblem zu lösen, indem wir die Funktionu diskretisieren durch ein äquidistantes Gitter xi = ih, i = 0, . . . n + 1,h := 1/(n+ 1). Es bezeichne

U := (u(x1), . . . , u(xn))T und F := (f(x1), . . . , f(xn))

T

die Auswertungen von u und f auf diesem Gitter.

Wir ersetzen die Ableitungen durch zentrale finite Differenzen

u′(x) ≈ Dh[U ](x) :=u(x+ h)− u(x− h)

2h

u′′(x) ≈ D2h[U ](x) :=

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

h2

(wobei wir am Rand u(0) = 0 = u(1) verwenden).

Aus der Gleichung L[u] = f ergibt sich so das LGS

LhUh = F

mit einer Matrix Lh ∈ R

n×n. Durch Lösung des LGS erhalten wir einenVektor

Uh := (u1, . . . , un)T ∈ Rn

von Approximationen an (u(x1), . . . , u(xn))T .

Finite Differenzen für ein einfaches Beispiel. Mit diesem Ansatz er-gibt sich für das einfache Beispiel −u′′ = f

f(x1)f(x2)

...f(xn)

︸︷︷︸

=:F

=

u′′(x1)u′′(x2)

...u′′(xn)

≈ h−2

2 −1 0−1 2 −1

. . . . . . −10 −1 2

︸︷︷︸

=:Lh

u(x1)u(x2)

...u(xn)

.

Wir können daher erwarten, dass wir durch Lösung von F = LhUh einenVektor Uh = (u1, . . . , un)

T aus Approximationen an (u(x1), . . . , u(xn))T er-

halten.

52


FD für die Diffusionsgleichung. Genauso diskretisieren wir

−u′′(x) + b(x)u′(x) + c(x)u(x) = f(x), u(0) = u(1) = 0

und erhalten

f(x1)f(x2)

...f(xn)

≈ h−2

d1 s1 0r2 d2 s2

. . . . . . sn−1

0 rn dn

u(x1)u(x2)

...u(xn)

mit

di = 2 + h2c(xi),

ri = −1 − hb(xi)/2,

si = −1 + hb(xi)/2.

Wiederum ergibt sich ein LGS F ≈ LhU , und wir können erwarten, dass dieLösung Uh := L−1

h F die wahren Lösungswerte in U approximiert.

Inhomogene Dirichlet-Bedingungen. Im Falle inhomogener Dirichlet-Bedingungen u(0) = α ∈ R, u(1) = β ∈ R ergibt sich

f(x1)f(x2)

...f(xn−1)

≈ h−2

d1 s1 0r2 d2 s2

. . . . . . sn−1

0 rn dn

u(x1)u(x2)

...u(xn)

+ h−2

r1α0...

snβ

.

Wir erhalten das LGS F ≈ LhU +Bh und können erwarten, dass

Uh := L−1h (F − Bh) ≈ U.

Neumann Randbedingungen. Neumann Randbedingungen u′(0) = α,u′(1) = β können behandelt werden, in dem wir die unbekannten Auswer-tungen an den Randwerten x0 und xn+1 zu den Vektoren hinzufügen

f(x0)f(x1)f(x2)

...f(xn)f(xn+1)

≈ h−2

r1 d1 s1 0r2 d2 s2

. . . . . . sn−1 00 rn dn sn

u(x0)u(x1)u(x2)

...u(xn)u(xn+1)

.

53


Um Gleichungen für u(x0) = u(0) und u(xn) = u(1) zu erhalten vewendenwir die Näherungen

u(−h) ≈ u(0)− hu′(0) = u(0)− αh

u(1 + h) ≈ u(1) + hu′(1) = u(1) + βh.

Damit ist

f(x0)f(x1)f(x2)

...f(xn−1)f(xn)

≈ h−2

d0 s0r1 d1 s1

. . . . . . . . .rn dn sn

rn+1 dn+1

u(x0)u(x1)

...u(xn+1)

+ h−2

r0(u(x0)− αh)0...

sn+1(u(xn+1) + βh)

= h−2

d0 + r0 s0r1 d1 s1

. . . . . . . . .rn dn sn

rn+1 dn+1 + sn+1

u(x0)u(x1)

...u(xn+1)

+ h−1

−r0α0...

sn+1β

Wiederum ergibt sich ein LGS f ≈ Lhu+ bh, und wir erwarten dass

uh := L−1h (f − bh) ≈ u.

1.8.3 Konsistenz, Stabilität und Konvergenz

Wir betrachten in diesem Abschnitt nur das spezielle Randwertproblem

L[u] := −u′′(x) + b(x)u′(x) + c(x)u(x) = f(x) x ∈ (0, 1)

mit homogenen Dirichletrandbedingungen und die dazugehörige Diskretisie-rung

LhUh = F

54


aus dem letzten Abschnitt. Außerdem nehmen wir an, dass b ∈ C2[0, 1],c ∈ C[0, 1] sowie c > 0 gilt2, und dass das Randwertproblem eine eindeutigeLösung u ∈ C4([0, 1]) besitzt.

Zuerst charakterisieren wir, wie gut die wahren Lösungswerte

U := (u(x1), . . . , u(xn))T

die diskretisierte Gleichung lösen.

Lemma 1.37Es existiert ein C > 0, so dass

‖LhU − F‖∞ ≤ Ch2.

Man sagt auch, das Differenzenverfahren hat Konsistenzordnung 2.

Beweis: Der i-te Eintrag von (i = 1, . . . , n, x0 = 0, xn+1 = 0) von LhU − Fist

D2h[U ](xi) + b(xi)Dh[U ](xi) + c(xi)u(xi)− f(xi).

Da u die DGL u′′(xi) + b(xi)u′(xi)+ c(xi)u(xi)− f(xi) = 0 löst, genügt es zu

zeigen, dass

Dh[u](xi) = u′(xi) +O(h2) und D2h[u](xi) = u′′(xi) +O(h2).

In der Tat erhalten wir durch Taylorentwichlung

u(xi + h) = u(xi) + hu′(xi) +1

2h2u′′(xi) +

1

3!h3u′′′(xi) +O(h4)

u(xi − h) = u(xi)− hu′(xi) +1

2h2u′′(xi)−

1

3!h3u′′′(xi) +O(h4)

und damit

Dh[u](xi) =u(xi + h)− u(xi − h)

2h=

2hu′(xi) +O(h3)

2h= u′(xi) +O(h2)

D2h[u](xi) =

u(xi + h)− 2u(xi) + u(xi − h)

h2=

h2u′′(xi) +O(h4)

h2

= u′′(xi) +O(h2),


2Da die stetige Funktion c auf dem Kompaktum [0, 1] ihr Minimum annimmt gilt damitsogar c(x) ≥ c0 := minx∈[0,1] c(x) > 0.

55


Aus Konsistenz (im Sinne von Lemma 1.37) folgt mit dem folgenden einfachenArgument Konvergenz

‖U − Uh‖∞ =∥∥L−1

h Lh(U − Uh)∥∥∞

≤∥∥L−1

h

∥∥∞‖LhU − F‖∞ ,

wenn wir zeigen können, dass Lh invertierbar ist und∥∥L−1

h

∥∥∞

(gleichmäßigin h) beschränkt ist. Die zweite Eigenschaft heißt auch Stabilität des Diffe-renzenverfahrens. Um die Stabilität zu zeigen, konstruieren wir eine Lösungw eines Randwertproblems, für den die zugehörigen Auswertungen W

LhW ≥ 1

erfüllen. Zusammen mit einer noch zu zeigenden Monotonieeigenschaft vonL−1h folgt dann

∥∥L−1

h

∥∥∞

=∥∥L−1

h 1

∥∥∞

≤∥∥L−1

h LhW∥∥∞

≤ maxx∈[0,1]

w(x).

Bemerkung 1.38Eine komponentenweise nicht-negative Matrix M = (mij)

ni,j=1 hat die Mono-

tonieeigenschaftx ≤ y =⇒ Mx ≤ My,

wobei die Ungleichheitszeichen für die Vektoren x, y,Mx,My ∈ Rn kompo-nentenweise zu verstehen sind.

Um die Monotonieeigenschaft von unserem L−1h zu zeigen, benötigen wir noch

ein Hilfsresultat:

Lemma 1.39 (Neumannsche Reihe)Ist R ∈ Rn×n und gilt ‖R‖ < 1 in einer submultiplikativen Matrixnorm, soist I −R invertierbar und es gilt (I − R)−1 =

∑∞k=0R

k.

Beweis: Betrachte die Folge der Partialsummen SN :=∑N

k=0Rk. Da

∥∥∥∥∥

N∑

k=0

Rk −M∑

k=0

Rk

∥∥∥∥∥≤

N∑

k=M+1

‖R‖k ∀N > M

und die geometrische Reihe∑∞

k=0 ‖R‖k konvergiert, folgt dass (SN)N∈N eineCauchy-Folge ist und damit konvergiert. Genauso folgt Rk → 0, so dasswegen

(I −R)

N∑

k=0

Rk =

N∑

k=0

Rk(I − R) = I − Rk+1 → I.

der Grenzwert S := limN→∞ SN =∑∞

k=0Rk die Inverse von (I −R) ist.

56


Lemma 1.40Ist A ∈ Rn×n eine strikt diagonaldominante Matrix mit positiven Diagonal-elementen und nicht-positiven Nichtdiagonalelementen, dann ist A invertier-bar und A−1 komponentenweise nicht-negativ.

Beweis: Wir zerlegen A = D−N in seinen Diagonal- und Nichtdiagonalan-teil. Nach Voraussetzung ist D ≥ 0 und R ≥ 0. Für R = D−1N gilt offenbar

A = D(I −R), R ≥ 0, ‖R‖∞ =∥∥D−1N

∥∥∞

< 1.

Aus Lemma 1.39 folgt die Invertierbarkeit von I − R und damit die von A.Außerdem folgt

A−1 = (I −R)−1D−1 =∞∑

k=0

RkD−1

Die Einträge von A−1 sind also Grenzwerte von Summen und Produktennicht-negativer Zahlen und damit nicht-negativ.

Lemma 1.41Es existieren h0 > 0 und C > 0 mit

∥∥L−1

h

∥∥∞

≤ C für alle 0 < h < h0.

Beweis: Nach Übungsaufgabe ?? existiert eine eindeutige Lösung w ∈ C4[0, 1]des Randwertproblems

−w′′(x) + b(x)w′(x) = 1 x ∈ (0, 1), w(0) = 0 = w(1).

Da w stetig ist, besitzt w sein globales Minimum in [0, 1]. Da in jedem innerenMinimum w′(x) = 0 ≤ w′′(x) gilt und damit die DGL nicht erfüllt sein kann,muss das Minimum auf dem Rand liegen und es folgt

w(x) ≥ 0 ∀x ∈ [0, 1].

w erfüllt

L[w] = −w′′(x) + b(x)w′(x) + c(x)w(x) = 1 + c(x)w(x).

Nach Lemma 1.37 existiert deshalb ein C ′ > 0, so dass für W = (w(x1), . . . , w(xn))T

und G = (1 + c(x1)w(x1), . . . , 1 + c(xn)w(xn))T gilt

‖LhW −G‖∞ ≤ C ′h2,

57


und damit insbesondere

LhW ≥ G− C ′h21.

Da c und w nicht-negativ sind, ist G ≥ 1 und es folgt

LhW ≥ 1− C ′h21.

Für hinreichend kleine h ist 1 − C ′h2 > 12

und die Matrix Lh erfüllt dieVoraussetzungen von Lemma 1.40. Mit Bemerkung 1.38 folgt dann

LhW ≥ 1

21 =⇒ W ≥ 1

2L−1h 1

und damit∥∥L−1

h

∥∥∞

=∥∥L−1

h 1

∥∥∞

≤ 2 ‖W‖∞ ≤ maxx∈[0,1]

w(x),


Folgerung 1.42Es existieren h0 > 0 und C > 0 mit

‖U − Uh‖∞ ≤ Ch2 für alle 0 < h < h0.

Beweis: Mit

‖U − Uh‖∞ =∥∥L−1

h Lh(U − Uh)∥∥ ≤

∥∥L−1

h

∥∥ ‖LhU − F‖

folgt die Behauptung aus Lemma 1.37 und Lemma 1.41.

1.9 Mehrdimensionale Randwertprobleme

Mit Finite Differenzen Verfahren lassen sich auch mehrdimensionale (partiel-le) Differentialgleichungen lösen. Wir betrachten dies zum Abschluss des Ka-pitels examplarisch an der mehrdimensionalen stationären Wärmeleitungs-gleichung

−∑ ∂2

∂x2j

u(x) = −∆u(x) = f(x). (1.7)

in einer beschränkten offenen Menge Ω ⊆ Rn. Dabei sei f ∈ C(Ω) eine steti-ge Funktion. Ähnlich wie in der Modellierung der eindimensionalen Diffusionkönnen wir uns f als die Verteilung angelegter Wärmequellen vorstellen und

58

1.9. MEHRDIMENSIONALE RANDWERTPROBLEME

die Lösung u beschreibt dann die sich im Gleichgewicht einstellende Tempe-ratur.

Es ist anschaulich klar, dass für Gleichgewichtstemperatur auch der Randdes Gebiets ∂Ω eine Rolle spielen wird, etwa wenn dieser Rand immer aufeiner konstanten Temperatur gehalten wird. Tatsächlich werden wir sehen,dass u durch (1.7) und Vorgabe von u|∂Ω eindeutig bestimmt ist.

Damit die Gleichung (1.8) und eventuelle Randwerte überhaupt einen Sinnergeben, betrachten wir als Lösungskandidaten nur Funktionen u ∈ C2(Ω)∩C(Ω) (sogenannte klassische Lösungen). Lösungen der Laplace-Gleichung∆u = 0 heißen auch harmonische Funktionen.

1.9.1 Das Maximumsprinzip

Satz 1.43 (Maximumsprinzip)Es sei f ∈ C(Ω) punktweise nicht-positiv und die Funktion u ∈ C2(Ω)∩C(Ω)erfülle

−∆u(x) = f(x) ≤ 0 ∀x ∈ Ω. (1.8)

Dann nimmt u sein Maximum auf dem Rand ∂Ω an (d.h. mindestens einglobales Maximum von u liegt auf ∂Ω).

Beweis: (i) Betrachte zunächst den Fall f(x) < 0 für alle x ∈ Ω.

Angenommen es gibt ein inneres Maximum, also y ∈ Ω mit

u(y) ≥ u(x) ∀x ∈ Ω.

Dann ist y insbesondere ein Maximum in jeder Koordinatenrichtung,also

∂2

∂x2j

u(y) ≤ 0 j = 1, . . . , n

und damit −∆u ≥ 0, was −∆u = f < 0 widerspricht. u kann also keininnereres Maximum haben. Da u als stetige Funktion auf dem Kom-paktum Ω aber mindestens ein Maximum besitzt, muss ein Maximumauf dem Rand liegen.

(ii) Nun sei f(x) ≤ 0 für alle x ∈ Ω. Angenommen es liegt kein Maximumauf dem Rand, dann gibt es ein inneres Maximum y ∈ Ω mit

u(y) ≥ u(x) ∀x ∈ Ω und u(y) > u(x) ∀x ∈ ∂Ω.

59


Mit diesem y definieren wir die Funktion

h(x) := ‖x− y‖2 =n∑

j=1

(xi − yi)2.

Da Ω beschränkt ist, ist ∂Ω kompakt. h ist auf also auf ∂Ω beschränktund es gilt h(y) = 0. Für hinreichend kleines δ > 0 nimmt deshalb

w(x) := u(x) + δh(x)

sein Maximum nicht auf dem Rand an. Aus ∆h(x) = 2n folgt aber

−∆w(x) = f − 2nδ < 0,

und wir erhalten den Widerspruch aus der in Teil (i) gezeigten Aussage.

Satz 1.44Sei f ∈ C(Ω).

(a) Ist f ≥ 0 und u ∈ C2(Ω) ∩ C(Ω) eine Lösung von −∆u = f ≥ 0 in Ω,so nimmt u sein Minimum auf dem Rand ∂Ω an. (Minimumsprinzip).

(b) Gilt für u, v ∈ C2(Ω) ∩ C(Ω)

−∆u ≤ −∆v in Ω und u ≤ v auf ∂Ω

so gilt u ≤ v in ganz Ω.

(c) Die Nullfunktion ist die einzige Lösung u ∈ C2(Ω) ∩ C(Ω) von

−∆u = 0, u|∂Ω.Eine Lösung u ∈ C2(Ω) ∩ C(Ω)

−∆u = f

ist also (wenn sie existiert) eindeutig durch f und u|∂Ω bestimmt.

(d) Erfüllen u1, u2 ∈ C2(Ω) ∩ C(Ω)

−∆u1 = f = −∆u2

so ist

‖u1 − u2‖∞ := maxx∈Ω

|u1(x)− u2(x)| = maxx∈∂Ω

|u1(x)− u2(x)|.

Die Lösung des Dirichlet-Problems hängen also (so sie denn existieren)stetig von den vorgegebenen Dirichlet-Randdaten ab.

60


(e) Es existiert ein (vom Gebiet Ω abhängiges) C > 0, so dass für alle u ∈C2(Ω) ∩ C(Ω)

‖u‖∞ = maxx∈Ω

|u(x)| ≤ maxx∈∂Ω

|u(x)|+ C supx∈Ω

|∆u|.

In diesem Sinne hängt eine Lösung von ∆u = f (so sie denn existiert)also auch stetig von der rechten Seite ab.

(f) Sind c, f ∈ C(Ω), c ≥ 0, f ≤ 0 und u ∈ C2(Ω) ∩ C(Ω) erfüllt

−∆u + cu = f ≤ 0,

so giltmaxx∈Ω

u(x) ≤ max0,maxx∈∂Ω

u(x)

Beweis: (a) folgt aus Anwendung des Maximumprinzips auf −u.

(b) folgt aus Anwendung des Maximumprinzips auf u− v.

(c) folgt aus Anwendung des Maximumsprinzips und des Minimumsprinzips.

(d) folgt aus Anwendung des Maximums- und Minimumprinzips auf u− v.

(e) Für supx∈Ω |∆u(x)| = ∞ ist die Aussage erfüllt. Sei also ∆u(x) be-schränkt.

Wähle R > 0 so groß, dass Ω ⊆ BR(0) := x ∈ Rn | ‖x‖ = R. Für

w(x) := R2 − 1

n‖x‖2

gilt −∆w = 2 und w(x) ≥ 0 für alle x ∈ Ω.

Definiere außerdem

v(x) := maxz∈∂Ω

|u(z)|+ w(x)

2supz∈Ω

|∆u(z)|

Dann ist

−∆v(x) = supz∈Ω

|∆u(z)| ≥ −∆u(x) ∀x ∈ Ω und v|∂Ω ≥ u|∂Ω

also nach (b) u ≤ v auf Ω. Genauso folgt u ≥ −v auf Ω und damit

|u(x)| ≤ |v(x)| = v(x) ≤ maxz∈∂Ω

|u(z)|+ w(x)

2supz∈Ω

|∆u(z)|

≤ maxz∈∂Ω

|u(z)|+ R2

2supz∈Ω

|∆u(z)|,

also folgt die Behauptung mit C := R2

2.

61


(f) Angenommen, es existiert ein y ∈ Ω mit

u(y) = maxx∈Ω

u(x) > maxx∈∂Ω

u(x) und u(y) > 0.

Dann gilt −∆u(y) = f(y)− c(y)u(y) ≤ 0 und wir erhalten den gleichenWiderspruch wie im Beweis von Satz 1.43.

1.9.2 Finite Differenzen im Mehrdimensionalen

Wie im Eindimensionalen betrachten wir finite Differenzen und definieren füreine Funktion u(x)

D+h,iu(x) :=

u(x+ hei)− u(x)

h

D−h,iu(x) :=

u(x)− u(x− hei)

h

Dh,iu(x) :=u(x+ hei)− u(x− hei)

2h

D2h,iu(x) := D+

h,iD−h,iu(x) = D−

h,iD+h,iu(x) =

u(x+ hei)− 2u(x) + u(x− hei)

h2

Wir verwenden für mehrfache partielle Ableitungen auch die Multiindex-Notation: Für α := (α1, . . . , αn) ∈ Nn

0 ist

Dαu =∂|α|u

∂xα11 · · ·∂xαn

n

,

und |α| = α1 + . . .+αn, α! = α1! · · ·αn!. Für α, β ∈ Nn0 , β ≤ α ist außerdem

(αβ

)= α!

β!(α−β)!wobei

α ≤ β :⇐⇒ αj ≤ βj für alle j = 1, . . . , n.

Mit dieser Notation definieren wir für u ∈ Ck(Ω) die Halbnormen

|u|Ck(Ω) = max|α|=k

maxx∈Ω

|Dαu(x)|.

und die Norm‖u‖Ck(Ω) = max

|α|≤kmaxx∈Ω

|Dαu(x)|.

Man rechnet leicht nach, dass letzteres tatsächlich eine Norm ist, die Ck(Ω)zu einem Banachraum (also einem vollständigen normierten Vektorraum)macht.

62


Lemma 1.45Sei x ∈ Ω, j ∈ 1, . . . , n und h > 0 hinreichend klein, so dass

x± thej ∈ Ω für alle t ∈ [−1, 1].

(a) Für u ∈ C2(Ω) ist

∣∣∣∣D+

h,iu(x)−∂

∂xi

u(x)

∣∣∣∣≤ h

2|u|C2(Ω).

(b) Für u ∈ C2(Ω) ist

∣∣∣∣D−

h,iu(x)−∂

∂xiu(x)

∣∣∣∣≤ h

2|u|C2(Ω).

(c) Für u ∈ C3(Ω) ist

∣∣∣∣Dh,iu(x)−

∂

∂xiu(x)

∣∣∣∣≤ h2

6|u|C3(Ω).

(d) Für u ∈ C4(Ω) ist

∣∣∣∣D2

h,iu(x)−∂2

∂x2i

u(x)

∣∣∣∣≤ h2

12|u|C4(Ω).

Beweis: wie im Eindimensionalen

Beispiel 1.46Wir betrachten die Gleichung

−∆u = f in Ω und u|∂Ω = 0,

zu gegebenen Quelltermen f : Ω → R, wobei Ω = (0, 1)2 ⊂ R

2 ein zweidi-mensionales Quadrat mit Seitenlänge 1 sei.

Wir diskretisieren Ω durch ein Punktegitter mit der Schrittweite

h = 1/(k + 1), k ∈ N,

und ordnen die inneren Punkte entsprechend der in der Vorlesung gemaltenSkizze an

x(i+(j−1)k) := (ih, jh) ∈ Ω, i, j = 1, . . . , k.

Wir setzen noch F = (f (i))i=1,...,k2 ∈ Rk2, f (i) = f(x(i)) und versuchen einenVektor von Approximationen

Uh = (u(i)h )i=1,...,k2 ≈ (u(i))i=1,...,k2 = U ∈ Rk2, u(i) = u(x(i))

zu finden. Für jeden inneren Gitterpunkt x(i) ∈ Ω erhalten wir durch dieobigen Differenzenverfahren 2. Ordnung die lineare Gleichung

f (i) = f(x(i)) = −∆u|x=x(i) ≈ −D2h,1u(x)|x=x(i) −D2

h,2u(x)|x=x(i)

= − 1

h2

(u(x(i) + he1)− 2u(x(i)) + u(x(i) − he1)

+ u(x(i) + he2)− 2u(x(i)) + u(x(i) − he2))

=1

h2

(4u(i) − u(i−1) − u(i+1) − u(i−k)) − u(i+k)

).

63


falls alle benachbarten Punkte ebenfalls innere Punkte sind. Ist ein benach-barter Punkt ein Randpunkt, so erhalten wir einen analogen Ausdruck, beidem der zum Randpunkt gehörige Summand (wegen u|∂Ω = 0) fehlt.

Insgesamt erhalten wir so ein lineares Gleichungssystem

LhUh = F

wobei die Matrix Lh die folgende Blocktridiagonalgestalt besitzt

Lh =1

h2

C −I−I C −I

−I. . . . . .. . . . . . −I

−I C

∈ Rk2×k2

mit

C :=

4 −1−1 4 −1

−1. . . . . .. . . . . . −1

−1 4

∈ Rk×k,

und der k × k-Einheitsmatrix I.

Bemerkung 1.47Für zwei Matrizen K = (ki,j) ∈ R

l×m und L ∈ R

r×s ist das KroneckerProdukt definiert durch

K ⊗ L =

k1,1L k1,2L . . . k1,mLk2,1L k2,2L . . . k1,mL

......

...kl,1L kl,2L . . . kl,mL

∈ Rlr×ms.

Damit lässt sich die Matrix Lh aus 1.46 schreiben als

Lh =1

h2(I ⊗ T + T ⊗ I)

mit

T :=

2 −1−1 2 −1

−1. . . . . .. . . . . . −1

−1 2

∈ Rk×k.

64


Inhomogene Dirichlet-Probleme

−∆u = f, u|∂Ω = g

lassen sich analog behandeln, indem entweder (wie im Eindimensionalen vor-geführt) der Effekt der Randpunkte auf die Differenzenverfahren der randna-hen Punkte in der rechten Seite berücksichtigt wird, oder indem die Rand-punkte als Unbekannte hinzugenommen werden und für jeden Randpunktx(j) die Gleichung

u(j) = g(x(j)).

aufgenommen wird.

1.9.3 Ein diskretes Maximumsprinzip

Offensichtlich lassen sich auch allgemeinere Gleichungen (mit veränderlichenDiffusionskoeeffizienten, Absorptions- und Konvektionstermen) analog be-handeln und führen wiederum auf lineare Gleichungssysteme der Form

LhUh = F.

Die so entstehenden Diskretisierungsmatrizen Lh sind üblicherweise von derGestalt, dass die Diagonaleinträge positiv sind, die Nebendiagonaleinträgenegativ sind, und zeilenweise der Diagonaleintrag die Zeile dominiert. Wirwerden sehen, dass solche Matrizen ein diskretes Maximumsprinzip erfüllenund benötigen dafür noch das folgende Lemma:

Lemma 1.48 (Sternlemma)Sei M ∈ N0 und αm, xm, m = 0, . . . ,M erfüllen

αm < 0 ∀m > 0,

M∑

m=0

αm ≥ 0,

M∑

m=0

αmxm ≤ 0 und x0 ≥ 0.

Ist x0 ≥ maxm=1,...,M xm, so gilt x0 = x1 = . . . = xM .

Beweis: Für M = 0 ist die Aussage trivial. Ansonsten ist

M∑

m=1

αm︸︷︷︸

<0

(xm − x0)︸︷︷︸

≤0

=M∑

m=0

αm(xm − x0) =M∑

m=0

αmxm

︸︷︷︸

≤0

−x0

M∑

m=0

αm

︸︷︷︸

≥0

≤ 0

und es folgt xm = x0 für alle m = 1, . . . ,M .

65


Offenbar kann für∑M

m=0 αm = 0 auf die Voraussetzung x0 ≥ 0 in Lemma 1.48verzichtet werden.

Satz 1.49Sei Lh = (ahij)

Ni,j=1 ∈ RN×N eine (nicht-notwendigerweise strikt) diagonaldo-

minante Matrix mit positiven Diagonal- und negativen Nebendiagonaleinträ-gen, also

ahii ≥∑

i 6=j

|ahij| ∀i = 1, . . . , N, ahij ≤ 0 ∀i 6= j.

Sei F ∈ RN und Uh ∈ RN sei eine Lösung von LhUh = F . Außerdem seiF ≤ 0, also komponentenweise nicht-positiv.

Betrachte die i-te Zeile des LGS LhUh = F ,

N∑

j=1

ahijuhj = fh

i ≤ 0.

Ist uhi ≥ max0,maxj: j 6=i,aij 6=0 u

hj , so ist uh

i = uhj für alle j mit aij 6= 0.

Beweis: Wir setzen α0 = ahii, x0 = uhi . Entsprechend seien αl und xl,

l = 1, . . . , k die anderen von Null verschiedenen Einträge ahij , j 6= i unddazugehörigen uh

j . Dann folgt die Behauptung aus Lemma 1.48.

Bemerkung 1.50(a) Im Kontext unserer Finite-Differenzen-Diskretisierungen kann Theorem

1.49 als diskretes Maximumsprinzip (in Analogie zu Theorem 1.44(f))interpretiert werden. Die diskrete Lösung Uh kann nicht in einem Gitter-punkt einen nicht-negativen Wert annehmen, der strikt maximal ist unterallen im Rahmen der Differenzenquotienten betrachten Nachbarwerten.Mehr noch: nimmt die diskrete Lösung einen in diesem Sinne maxima-len Wert an, so besitzen alle (in diesem Sinne) benachbarten Werte dengleichen maximalen Wert.

(b) Mit dem diskreten Maximumsprinzip lässt sich oft die Lösbarkeit der dis-kretisierten Gleichung beweisen. Betrachten wir exemplarisch die MatrixLh aus der Diskretisierung des Poisson-Problems auf dem Einheitsqua-drat. Lh ist genau dann invertierbar, wenn es injektiv ist, also nur derNullvektor das homogenen LGS LhUh = 0 löst. Sei also Uh eine solche Lö-sung. Dann besitzt Uh einen maximalen Eintrag uh

j . O.B.d.A. sei ujh ≥ 0,

ansonsten betrachten wir −Uh. Da wir jeden anderen Gitterpunkt übereinen Weg aus (bei den Differenzenquotienten vorkommenden) Nachbar-werten erreichen können (das Gitter ist diskret zusammenhängend), folgt

66


aus dem diskreten Maximumsprinzip, dass alle Einträge von uhj überein-

stimmen. Aus der ersten Zeile von Lh folgt dann, dass Uh = 0 geltenmuss.

Die folgende Verschärfung von Lemma 1.40 können wir als diskretes Analogonzur Monotonieeigenschaft aus Theorem 1.44(b) interpretieren:

Satz 1.51Sei A ∈ RN×N eine invertierbare, diagonaldominante Matrix mit positivenDiagonalelementen und nicht-positiven Nichtdiagonalelementen, wobei

aii ≥N∑

j=1

|aij| ∀i ∈ 1, . . . , n.

Dann ist A−1 komponentenweise nicht-negativ.

Insbesondere gilt (komponentenweise)

Au ≤ Av =⇒ u ≤ v.

Beweis: (a) Nach Lemma 1.40 ist jede strikt-diagonaldominante MatrixA ∈ RN×N (mit positiven Diagonalelementen und nicht-positiven Nicht-diagonalelementen) invertierbar und A−1 ist komponentenweise nicht-negativ.

(b) Für (nicht notwendigerweise strikt) diagonaldominantes und invertierba-res A (mit positiven Diagonalelementen und nicht-positiven Nichtdiago-nalelementen) erhalten wir aus Teil (a), dass (A + ǫI) invertierbar ist,und dass (A+ǫI)−1 komponentenweise nicht-negativ ist. Da A nach Vor-aussetzung invertierbar ist, konvergiert (für ǫ → 0) (A + ǫI)−1 → A−1

und es folgt die Behauptung.

1.9.4 Konsistenz, Stabilität und Konvergenz

Sei LhUh = F die entsprechend den letzten Abschnitten erstellte Diskretisie-rung des Poisson-Problems

−∆u = f, u|∂Ω = g

und U ∈ RN der Vektor der Auswertungen der Funktion u(x) auf den Git-terpunkten der Diskretisierung.

67


Wie im Eindimensionalen gilt der folgende Zusammenhang zwischen Konsi-stenz, Stabilität und Konvergenz:

‖U − Uh‖∞ =∥∥L−1

h Lh(U − Uh)∥∥∞

≤∥∥L−1

h

∥∥∞‖LhU − F‖∞ .

Erfüllt also (für h → 0) die wahre Lösung (genauer: ihre Auswertungen) im-mer besser die diskretisierte Gleichung (Konsistenz), d.h. limh→0 ‖LhU − F‖∞ =0, und bleibt

∥∥L−1

h

∥∥∞

beschränkt (Stabilität), so nähert die Lösung der dis-kreten Gleichung die wahre Lösung immer besser an.

Die Konsistenz von finiten Differenzenverfahren erhalten wir sofort aus Ab-schätzungen, wie sie in Lemma 1.45 vorkommen, die ihrerseits leicht ausAnwendungen der Taylor-Formel folgen. Auch im Mehrdimensionalen las-sen sich konsistente Differenzenverfahren meist sehr leicht für eine gegebenepartielle Differentialgleichung konstruieren.

Der Beweis von Stabilität ist wie im Eindimensionalen ungleich schwerer.Wir zeigen im Rahmen dieser Vorlesung die Stabilität nur für das hier exem-plarisch betrachtete Problem −∆u = f mit homogenen Dirichletrandwertenauf dem Einheitsquadrat.

Satz 1.52Sei Ω = (0, 1)2 ⊆ R2. Dann existiert ein C > 0, so dass für die entsprechendden letzten Abschnitten erstellte Diskretisierungsmatrix Lh zum Problem

−∆u = f, u|∂Ω = 0

gilt∥∥L−1

h

∥∥∞

< C für alle h > 0.

Beweis: Sei R > 0 hinreichend groß, so dass Ω ⊆ BR und betrachte dieFunktion w(x) = 1

4

(R2 − ‖x‖2

)+ 1. Offenbar gilt −∆w = 1. w ist ein Poly-

nom zweiten Grades und die verwendeten Differenzenquotienten sind dafürin allen nicht randnahen Gitterpunkten exakt, so dass an diesen Punkten(LhW )j = 1 gilt. (Dabei bezeichnet W wieder den Vektor der Auswertun-gen von w(x) in den Gittepunkten.) Für die randnahen Punkte bleiben eineroder mehrere positive Randwerte w(x)|∂Ω ≥ 0 mit einem negativen Gewichtunberücksichtigt, so dass in diesen Punkten (LhW )j ≥ 1 gilt. Insgesamt istalso

LhW ≥ 1.

Damit folgt

∥∥L−1

h

∥∥∞

=∥∥L−1

h 1

∥∥∞

≤∥∥L−1

h LhW∥∥∞

≤ maxx∈BR(0)

|w(x)| = R2

4.

und damit die Behauptung.

68


Folgerung 1.53Falls eine Lösung des homogenen Dirichletproblems

−∆u = f, u|∂Ω = 0, Ω = (0, 1)2 ⊆ R2

existiert, so konvergieren die durch Finite Differenzen erhaltenenen Appro-ximationen Uh gegen die Lösung.

69


70

Kapitel 2

Eigenwertprobleme

2.1 Einleitung

Wir untersuchen in diesem Kapitel, wie sich die Eigenwerte einer gegebenenMatrix numerisch bestimmen lassen. In diesem Kapitel sei A ∈ Cn×n stetseine quadratische Matrix mit Einträgen aij ∈ C. Oft werden wir uns dabei aufden Fall hermitescher A = A∗ (oder zumindest diagonalisierbarer) Matrizenbeschränken. Für die Praxis reicht dies oft aus, da für allgemeine Matrizendamit die (auf den Eigenwerten von A∗A beruhende) Singulärwertzerlegungbestimmt werden kann.

Wir erinnern zuerst an ein paar grundlegende Resultate aus der linearenAlgebra.

Definition 2.1Zu einer quadratischen Matrix A ∈ Cn×n heißt λ ∈ C Eigenwert, falls einVektor 0 6= x ∈ Cn existiert, so dass

Ax = λx.

x heißt Eigenvektor.

Die Menge aller Eigenwerte heißt Spektrum

σ(A) := λ ∈ C : ∃x ∈ Cn \ 0 : Ax = λx .

Der Betrag des betragsgrößten Eigenwertes heißt Spektralradius

ρ(A) = max|λ| : λ ∈ σ(A).

71

KAPITEL 2. EIGENWERTPROBLEME

Bemerkung 2.2(a) λ ist genau dann Eigenwert von A, wenn A − λI nicht injektiv ist. Da

A−λI quadratisch ist, ist das äquivalent dazu, dass A−λI nicht surjektivist. N (A− λI) ist der Eigenraum zum Eigenwert λ. Die Dimension desEigenraums heißt auch geometrische Vielfachheit des Eigenwerts.

(b) Eigenwerte sind genau die Nullstellen des komplexen charakteristischenPolynoms p(λ) = det(A − λI) vom Grad n. Aus dem Fundamentalsatzder Algebra folgt deshalb, dass jedes A ∈ C

n×n mindestens einen undhöchstens n Eigenwerte besitzt. Die Vielfachheit des Eigenwerts als Null-stelle von p heißt algebraische Vielfachheit.

(c) Eine Matrix ist genau dann diagonaliserbar, wenn eine Basis aus Eigen-vektoren besteht. Dies ist genau dann der Fall wenn für jeden Eigenwertdie geometrische und algebraische Vielfachheit übereinstimmen.

(d) ‖A‖ ist der größte Singulärwert von A, d.h. die Wurzel des größten Ei-genwertes der Matrix A∗A. Es gilt also

‖A‖ =√

ρ(A∗A).

(e) Für eine hermitesche Matrix A = A∗ sind alle Eigenwerte reell und esexistierte eine Orthogonalbasis aus Eigenvektoren. Offenbar folgt damitauch

‖A‖ =√

ρ(A2) = ρ(A).

2.2 Zwei Einschließungssätze

Wir beginnen mit zwei Resultaten, mit denen sich die Lage der Eigenwerteauf recht einfache Weise grob abschätzen lässt. Zuerst zeigen wir, dass alleEigenwerte in der Vereinigung der sogenannten Gerschgorin-Kreise liegenmüssen.

Satz 2.3 (Satz von Gerschgorin)Für jedes A ∈ Cn×n gilt

σ(A) ⊂n⋃

i=1

Ki, Ki :=

ζ ∈ C : |ζ − aii| ≤n∑

j=1j 6=i

|aij|

Die Ki heißen Gerschgorin-Kreise.

72

2.2. ZWEI EINSCHLIEßUNGSSÄTZE

Beweis: Sei λ ∈ C ein Eigenwert von A, also Ax = λx mit einem x 6= 0. Seii ein Index eines betragsmäßig größten Eintrags von x, also

|xi| ≥ |xj | ∀j 6= i.

Aus λxi = (Ax)i =∑n

j=1 aijxj folgt

|λ− aii| =∣∣∣∣∣

n∑

j=1

aijxj

xi− aii

∣∣∣∣∣=

∣∣∣∣∣∣∣∣

n∑

j=1j 6=i

aijxj

xi

∣∣∣∣∣∣∣∣

≤n∑

j=1j 6=i

|aij |∣∣∣∣

xj

xi

∣∣∣∣≤

n∑

j=1j 6=i

|aij|

und damit λ ∈ Ki.

Bemerkung 2.4Eine Matrix und ihre Transponierte AT = (aji)

ni,j=1 haben dieselbe Determi-

nante, detA = detAT . Aus det(A − λI) = det(AT − λI) folgt, dass A undAT dieselben Eigenwerte haben. Genauso folgt, dass die Eigenwerte von A∗

die komplex konjugierten Eigenwerte von A sind.

Aus Satz 2.3 folgt deshalb auch dass

σ(A) = σ(AT ) ⊂n⋃

i=1

Ki, Ki := ζ ∈ C : |ζ − aii| ≤n∑

j=1j 6=i

|aji|

Beispiel 2.5Betrachte

A =

2 1 −10 −1 11 1 1

.

Die Gerschgorin-Kreise für A sind die (in der komplexen Ebene befindlichen)Kreise um 2 mit Radius 2, um −1 mit Radius 1 und um 1 mit Radius 2. AlleEigenwerte von A liegen in diesen Kreisen.

Die Gerschgorin-Kreise für AT sind die Kreise um 2 mit Radius 1, um −1mit Radius 2 und um 1 mit Radius 2. Alle Eigenwerte von A liegen auch indiesen Kreisen.

Das zweite Einschließungsresultate verwendet das Konzept des Wertebereichseiner Matrix.

73


Definition 2.6Die Menge

W (A) =

ζ =x∗Ax

x∗x: 0 6= x ∈ Cn

= x∗Ax : x ∈ Cn, ‖x‖ = 1 ⊆ C

heißt Wertebereich von A. Ausdrücke der Form x∗Axx∗x

heißen auch Rayleigh-Quotienten.

Lemma 2.7Es gilt

σ(A) ⊂ W (A).

Beweis: Für einen Eigenwert λ ∈ C und zugehörigen Eigenvektor x ∈ Cn

ist x∗Axx∗x

= x∗λxx∗x

= λ.

Lemma 2.8(a) W (A) ist wegzusammenhängend, d.h. zu ζ0, ζ1 ∈ W (A) existiert stets

eine stetige Abbildung

α : [0, 1] 7→ W (A),

mit α(0) = ζ0 und α1 = ζ1.

(b) Ist A ∈ Cn×n hermitesch (d.h. A = A∗), dann ist W (A) ⊂ R die konve-xe Hülle der Eigenwerte, d.h. das reelle Intervall [λn, λ1] zwischen demgrößten Eigenwert λ1 und dem kleinsten Eigenwert λn.

(c) Ist A ∈ C

n×n schiefhermitesch (d.h. A = −A∗), dann ist W (A) ⊂ iRdie konvexe Hülle der Eigenwerte, d.h. das Intervall auf der imaginärenAchse i[a, b], wobei a der kleinste und b der größte Imaginärteil einesEigenwertes sind.

Beweis: (a) Seien ζ0, ζ1 ∈ W (A). O.B.d.A. ζ0 6= ζ1.

Aus der Definition von W (A) folgt dass

ζ0 =x∗0Ax0

x∗0x0

, ζ1 =x∗1Ax1

x∗1x1

mit x0 6= 0 6= x1.

Der naheliegende Kandidat für α ist

α(t) :=x∗tAxt

x∗txt

, xt := x0 + t(x1 − x0)

74

2.2. ZWEI EINSCHLIEßUNGSSÄTZE

Ist xt 6= 0 für alle t ∈ [0, 1], so erfüllt α offensichtlich das Gewünschte.Wir müssen also nur zeigen, dass xt 6= 0 für alle t ∈ [0, 1]. Angenommendies wäre nicht der Fall, also 0 = x0+t(x1−x0). Dann wäre aber offenbart 6= 0 und x1 =

t−1tx0, womit ζ0 = ζ1 und damit einen Widerspruch folgt.

(b) Ist A = A∗, so ist für alle x ∈ Cn

x∗Ax

x∗x=

(x∗Ax)∗

x∗x=

(Ax)∗x

x∗x=

x∗A∗x

x∗x=

x∗Ax

x∗x,

also W (A) ⊆ R.

Nach Lemma 2.7 gilt λn, λ1 ⊆ W (A) und zusammen mit Lemma 2.7folgt

[λn, λ1] ⊆ W (A) ⊆ R.

Wir müssen also nur noch zeigen, dass für alle ζ ∈ W (A) gilt λn ≤ ζ ≤ λ1.Für α > 0 hinreichend groß (z.B. α := 2 ‖A‖+ 1) ist A+ αI hermiteschund

x∗(A+ αI)x ≥ −‖A‖ ‖x‖2 + α ‖x‖2 > 0 ∀x 6= 0,

also A + αI positiv definit. A besitzt also eine Cholesky-Zerlegung A =LL∗ mit L ∈ Cn×n. Es ist

x∗Ax = x∗(A + αI)x− α ‖x‖2 = ‖L∗x‖2 − α ‖x‖2 ≤ (‖L∗‖2 − α) ‖x‖2

und mit Bemerkung 2.2 folgt

x∗Ax

x∗x≤ ‖L∗‖2 − α = ρ(LL∗)− α = ρ(A+ αI)− α.

Da A+αI positiv definit ist, sind alle Eigenwerte positiv und der größteist offenbar λ1 + α. Es gilt also ζ ≤ λ1 für alle ζ = x∗Ax

x∗x∈ W (A).

Der größte Eigenwert von −A ist −λn. Mit dem gleichen Argument er-halten wir deshalb

x∗(−A)x

x∗x≤ −λn

und damit ζ ≥ λn für alle ζ = x∗Axx∗x

∈ W (A).

(c) Ist A∗ = −A, dann ist (iA)∗ = −iA∗ = iA, also iA hermitesch. Daoffenbar W (iA) = iW (A) ist, folgt (c) aus (b).

Jede Matrix A ∈ Cn×n lässt sich in einen hermiteschen und einen schiefher-miteschen Teil zerlegen:

A =A+ A∗

2+

A− A∗

2.

75


Nach Lemma 2.8 ist

W

(A + A∗

2

)

+W

(A− A∗

2

)

=

ζ = ζ1 + ζ2 ∈ C : ζ1 ∈ W (A+ A∗

2), ζ2 ∈ W (

A−A∗

2)

=

ζ = ξ + iη ∈ C : ξ ∈ W (A+ A∗

2), iη ∈ W (

A− A∗

2)

.

ein Rechteck in der komplexen Ebene. Der folgende Satz zeigt, dass es W (A)und damit insbesondere alle Eigenwerte von A enthält.

Satz 2.9 (Satz von Bendixson)Für jedes A ∈ Cn×n gilt

W (A) ⊆ W (A+ A∗

2) +W (

A− A∗

2)

Insbesondere ist also σ(A) ⊂ W (A+A∗

2) +W (A−A∗

2).

Beweis: Für jedes x ∈ Cn mit ‖x‖ = 1 ist

x∗Ax = x∗

(A+ A∗

2+

A−A∗

2

)

x ∈ W (A+ A∗

2) +W (

A− A∗

2).

W (A+A∗

2) +W (A−A∗

2) lässt sich mit Gerschgorin-Kreisen abschätzen.

2.3 Stabilität des Eigenwertproblems

Der naheliegender Ansatz zur Berechnung der Eigenwerte von A ist, zuerstdas charakteristische Polynom aus A zu berechnen und dann dessen Nullstel-len zu bestimmen. In der Praxis wird man aber A und damit die Koeffizientendes char. Polynoms nur näherungsweise kennen. Die Berechnung der Null-stellen eines Polynoms aus seinen Koeffizienten ist jedoch im Allgemeinenschlecht konditioniert (d.h. stark fehlerverstärkend), vgl. die an der Tafelgemalte Skizze.

Das Eigenwertproblem selbst ist (zumindest für normale Matrizen) nichtschlecht konditioniert, wie wir im folgenden Satz zeigen. Zur Erinnerung:

76

2.3. STABILITÄT DES EIGENWERTPROBLEMS

A ∈ C

n×n heißt normal, falls AA∗ = A∗A. Normale Matrizen sind diago-nalisierbar bzgl. einer Orthogonalbasis aus Eigenvektoren. Es existieren alsounitäre Matrizen X ∈ Cn×n (d.h. X−1 = X∗) mit

A = XΛX−1, Λ =

λ1 0 0

0. . . 0

0 0 λn

Satz 2.10Ist A ∈ Cn×n normal und E ∈ Cn×n beliebig. Ist λ ein Eigenwert der gestör-ten Matrix A+ E, dann existiert ein λ ∈ σ(A) mit

|λ− λ| ≤ ‖E‖ .

Beweis: Sei λ ein Eigenwert von A+E und x 6= 0 ein dazugehöriger Eigen-vektor. Ist λ ∈ σ(A), dann stimmt die Behauptung. Sei also λ 6∈ σ(A). Dannexistiert (λI − A)−1 und aus

Ex = (A+ E)x− Ax = (λI − A)x

folgt x = (λI − A)−1Ex und damit

1 ≤∥∥∥(λI − A)−1E

∥∥∥ ≤

∥∥∥(λI − A)−1

∥∥∥ ‖E‖ =

∥∥∥(λI −XΛX−1)−1

∥∥∥ ‖E‖

≤ ‖X‖∥∥X−1

∥∥

∥∥∥(λI − Λ)−1

∥∥∥ ‖E‖

Da X und X−1 unitär sind, ist ‖X‖ = ‖X−1‖ = 1. Außerdem ist (λI −Λ)−1

eine Diagonalmatrix mit den Einträgen 1λ−λj

wobei λj die Eigenwerte von

A sind.∥∥∥λI − Λ

∥∥∥ ist der betragsgrößte Diagonaleintrag, d.h. es existiert ein

λ ∈ σ(A) mit∥∥∥(λI − Λ)−1

∥∥∥ =

1

|λ− λ|und damit

|λ− λ| ≤ ‖E‖ .

Ein noch stärkeres Stabilitätsresultat ist der folgende Satz, den wir ohneBeweis angeben.

77


Satz 2.11 (Satz von Wielandt-Hoffman)Seien A,E ∈ Cn×n hermitesch. λ1 ≥ . . . ≥ λn und λ1 ≥ . . . ≥ λn seien dieEigenwerte von A und A+ E.

Dann giltn∑

i=1

(λi − λi)2 ≤ ‖E‖2F .

2.4 Die Potenzmethode

Wir stellen nun die erste iterative Methode zur Berechnung von Eigenwer-ten vor, die Potenzmethode (auch: von Mises-Verfahren). Zur Motivationbetrachten wir eine Diagonalmatrix

A =

λ1 0 0

0. . . 0

0 0 λn

Wiederholtes Anwenden von A auf einen Vektor x = (xi)ni=1 ∈ Cn liefert

Akx =

λk1x1...

λknxn

=

n∑

i=1

λki xiei,

wobei ei der i-te Einheitsvektor ist. Gibt es genau einen betragsgrößten Ei-genwert λj und ist xj 6= 0, dann wird der dazugehörige Summand schnellerwachsen (bzw. langsamer fallen für |λj| < 1) als alle anderen Summanden.Er wird die Summe also immer mehr dominieren. Wir erwarten daher, dass

Akx ≈ λkjxjej .

Das gleiche gilt für beliebige diagonalisierbare Matrizen, wie der folgendeSatz zeigt:

Satz 2.12 (Potenzmethode)Sei A ∈ C

n×n diagonalisierbar. A besitze einen genau einen (aber nichtnotwendigerweise einfachen) betragsgrößten Eigenwert, d.h. es sei σ(A) =λ1, . . . , λm

|λ1| > |λ2| ≥ . . . ≥ |λm|.

Zu einem normierten Startvektor x(0) ∈ Cn definieren wir die Folge

x(k) := Ax(k−1), x(k) :=x(k)

‖x(k)‖ .

78

2.4. DIE POTENZMETHODE

An x(0) sei ein Eigenvektor zu λ1 beteiligt, d.h. es sei

x(0) = ξ1v1 +m∑

j=2

ξjvj, ξ1, . . . , ξm ∈ C, ξ1 6= 0.

wobei vj normierte Eigenvektoren zu λj seien. Dann gilt:

x(k) =λk1

|λ1|kξ1|ξ1|

v1 +O(|λ2/λ1|k),

(x(k−1))∗x(k) = λ1 +O(|λ2/λ1|k).

(x(k−1))∗x(k) konvergiert also gegen den betragsgrößten Eigenwert λ1 und x(k)

konvergiert bis auf skalare Vielfache gegen einen dazugehörigen Eigenvektor.Die Konvergenzgeschwindigkeit ist linear mit Konvergenzfaktor |λ2/λ1|.Beweis: Durch triviale Induktion sieht man, dass

x(k) =Akx(0)

‖Akx(0)‖ .

Es ist

Akx(0) = λk1ξ1v1 +

m∑

j=2

λkj ξjvj

und∣∣∥∥Akx(0)

∥∥− |λk

1ξ1|∣∣ ≤

∥∥∥∥∥

m∑

j=2

λkj ξjvj

∥∥∥∥∥≤ |λ2|k max |ξj|,

also∥∥Akx(0)

∥∥ = |λ1|k|ξ1|

(

1 +O(|λ2/λ1|k))

.

Damit erhalten wir

x(k) =Akx(0)

‖Akx(0)‖ =Akx(0)

|λ1|k|ξ1|1

1 +O(|λ2/λ1|k)=

Akx(0)

|λ1|k|ξ1|(

1 +O(|λ2/λ1|k))

=

(λk1

|λ1|kξ1|ξ1|

v1 +O(|λ2/λ1|k))(

1 +O(|λ2/λ1|k))

=λk1

|λ1|kξ1|ξ1|

v1 +O(|λ2/λ1|k).

und

(x(k−1))∗x(k) = (x(k−1))∗Ax(k−1) = v∗1Av1 +O(|λ2/λ1|k) = λ1 +O(|λ2/λ1|k).womit die Behauptung gezeigt ist.

79


Bemerkung 2.13Die Konvergenz der Potenzmethode kann auch ohne die Annahme der Dia-gonalisierbarkeit gezeigt werden. Auf die Voraussetzung, dass genau ein be-tragsgrößter Eigenwert existiert, kann jedoch nicht verzichtet werden.

Algorithm 1 Potenzmethode (von Mises-Verfahren)Gegeben Matrix A ∈ Cn×n und Startvektor x ∈ Cn.repeatx := Axλ := x∗xx := x/ ‖x‖

until STOPreturn λ ∈ C, x ∈ Cn approximiert betragsgrößten Eigenwert von A unddazugehörigen Eigenvektor (unter den Voraussetzungen von Satz 2.12).

Bemerkung 2.14Auch zur Bestimmung der anderer Eigenwerte kann die Potenzmethode ver-wendet werden, indem A zuerst geignet transformiert wird. Sei A ∈ C

n×n

diagonalisierbar.

(a) Inverse Iteration: Sei A invertierbar. Die betragsgrößten Eigenwerte vonA−1 sind dann 1/λj, wobei λj die betragskleinsten Eigenwerte von A sind.Existiert genau ein betragskleinster Eigenwert, so lässt sich dieser durchAnwendung der Potenzmethode auf A−1 bestimmen.

(b) Gebrochene Iteration von Wielandt: Ist µ 6∈ σ(A), dann ist A−µI inver-tierbar. Die betragsgrößten Eigenwerte von (A−µI)−1 sind 1

λj−µwobei λj

diejenigen Eigenwerte sind die am nächsten an µ liegen. Gibt es davongenau einen, so lässt er sich durch Anwendung der Potenzmethode auf(A− µI)−1 bestimmen.

2.5 Die Rayleigh-Quotient-Iteration

Das in Bemerkung 2.14(b) skizzierte Verfahren wird umso schneller konver-gieren, umso besser µ den gesuchten Eigenwert approximiert. Es liegt dahernahe, in jedem Schritt µ durch die aktuelle Näherung an den Eigenwert zuersetzen.

Für die Konvergenzuntersuchung dieses Verfahrens benötigen wir den Begriffdes Betrags einer normalen Matrix.

80

2.5. DIE RAYLEIGH-QUOTIENT-ITERATION

Algorithm 2 Rayleigh-Quotient-IterationGegeben Matrix A ∈ Cn×n und Startwerte x ∈ Cn, µ ∈ C.repeat

Löse (A− µI)x = xx := x/ ‖x‖µ := x∗Ax

until STOPreturn µ ∈ C approximiert Eigenwert, der am nächsten am Startwertliegt, x approximiert dazugehörigen Eigenvektor (vgl. Satz 2.17 für dierigorose Formulierung der Konvergenzaussage).

Definition 2.15Für eine Diagonalmatrix

Λ =

λ1 0 0

0. . . 0

0 0 λn

definieren wir |Λ| durch

Λ =

|λ1| 0 0

0. . . 0

0 0 |λn|

.

Sei A ∈ Cn×n normal. Dann existiert eine Diagonalmatrix Λ (mit den Ei-genwerten auf der Diagonale) und eine unitäre Matrix X ∈ C

n×n so dassA = XΛX−1. Damit definieren wir

|A| = X|Λ|X−1.

Offenbar ist |A| dadurch wohldefiniert.

Lemma 2.16Für jedes normale A ∈ Cn×n gilt

(a) |A| ist hermitesch, positiv semidefinit und es gilt |A|2 = A∗A.

(b) Für alle x ∈ Cn ist|x∗Ax| ≤ x∗|A|x.

81


Beweis: (a) Aus |Λ| = |Λ|∗ und X∗ = X−1 folgt die Hermitizität. Sei(vj)

nj=1 ⊂ C

n eine ONB aus Eigenvektoren von A und λ1, . . . , λn ∈ C

seien die zugehörigen Eigenwerte. Dann gilt für jedes x =∑n

j=1 ξjvj

A∗Ax = A∗n∑

j=1

λjξjvj =n∑

j=1

λjλjξjvj = |A|2x.

(b) Zum Beweis von (b) beachte, dass für jede Diagonalmatrix Λ ∈ C

n×n

und jedes y ∈ Cn offenbar gilt

|y∗Λy| =∣∣∣∣∣

n∑

j=1

λj |yj|2∣∣∣∣∣≤

n∑

j=1

|λj||yj|2 = y∗|Λ|y.

Damit folgt

|x∗Ax| = |x∗XΛX−1x| = |(X∗x)Λ(X∗x)| ≤ (X∗x)|Λ|(X∗x) = x∗|A|x.

Satz 2.17Sei A ∈ C

n×n normal. µ0 ∈ C und x(0) ∈ C

n seien hinreichend gute Nä-herungen an einen Eigenwert λ ∈ C und einen dazugehörigen Eigenvektorv ∈ Cn von A. µk und x(k) seien gemäß Algorithmus 2 definiert, also

x(k) := (A− µk−1I)−1x(k−1), x(k) :=

x(k)

‖x(k)‖ , µk := (x(k))∗Ax(k).

Dann erreicht µk entweder nach endlich vielen Schritten einen Eigenwertoder es existiert eine kubisch konvergente Folge ǫk → 0 mit

|λ− µk+1| ≤ ǫk ∀k ∈ N0.

Beweis: µk erreiche nicht nach endlich vielen Schritten einen Eigenwert, sodass die Rayleigh-Quotient-Iteration wohldefiniert ist. O.B.d.A. sei außerdemx(0) normiert.

(a) Eigenraumzerlegung und Konvergenz

Jedes x(k) kann eindeutig in seine Eigenraumanteile zerlegt werden

x(k) = v(k) + w(k), w(k) =∑

λ∈σ(A)\λ

w(k)

λ,

82

2.5. DIE RAYLEIGH-QUOTIENT-ITERATION

wobei v(k) ein Eigenvektor zu λ sei und w(k)

λEigenvektoren zu den an-

deren Eigenwerten λ sind. Da sie Eigenvektoren zu verschiedenen Eigen-werten sind, stehen all diese Vektoren paarweise orthogonal aufeinander.

Ähnlich wie im Beweis von Satz 2.12 kann man zeigen, dass (wennµ0 ∈ C und x(0) ∈ C

n hinreichend gute Näherungen an λ und einendazugehörigen Eigenvektor v sind) w(k) → 0 konvergiert, also x(k) immermehr zu einem Eigenvektor zum Eigenwert λ wird. Damit folgt auchµk = (x(k))∗Ax(k) → λ.

(b) Definition der ǫk:

Nach Lemma 2.16 gilt

|λ− µk| = |λ− (x(k))∗Ax(k)| = |(x(k))∗(λI − A)(x(k))∗|≤ (x(k))∗|λI − A|(x(k))∗ =: ǫk.

Wir werden zeigen, dass ein C > 0 existiert so dass (wenn µ0 ∈ C undx(0) ∈ Cn seien hinreichend gute Näherungen an λ und v sind)

ǫk+1 ≤ Cǫ3k (2.1)

Liegt x(0) hinreichend nahe an v, dann ist ǫ0 hinreichend klein, so dassaus (2.1) auch ǫk → 0 folgt.

(c) Wir zeigen nun, dass

ǫk+1 ≤ǫ3k

‖v(k)‖2max

λ∈σ(A)\λ|λ− µk|−2 (2.2)

Da∥∥v(k)

∥∥ =

∥∥x(k) − w(k)

∥∥→ 1 und µk → λ folgt daraus (2.1) und damit

die Behauptung.

Zum Beweis von (2.2) verwenden wir die Iterationsvorschriften und er-halten

ǫk+1 = (x(k+1))∗|λI − A|x(k+1) =(x(k+1))∗|λI − A|x(k+1)

‖x(k+1)‖2

=(x(k))∗(A− µkI)

−∗|λI − A|(A− µkI)−1x(k)

‖(A− µkI)−1x(k)‖2

Wir schätzen Zähler und Nenner ab.

Nenner: Wir beginnen mit dem Nenner und benutzen unsere orthogo-nale Zerlegung x(k) = v(k) + w(k). Offensichtlich ist

(A− µkI)−∗(A− µkI)

−1v(k) =1

|λ− µk|2v(k).

83


Aus v(k) ⊥ w(k) folgt deshalb

∥∥(A− µkI)

−1x(k)∥∥2

= (x(k))∗(A− µkI)−∗(A− µkI)

−1x(k)

= (v(k))∗|λ− µk|−2v(k) + (w(k))∗(A− µkI)−∗(A− µkI)

−1w(k)

≥ (v(k))∗|λ− µk|−2v(k) ≥ ǫ−2k

∥∥v(k)

∥∥2.

Zähler: Nun schätzen wir den Zähler ab. Mit der Zerlegung von w(k) inseine Eigenraumanteile w

(k)

λerhalten wir aufgrund der paarweisen Or-

thogonalität der Eigenraumanteile

(x(k))∗(A− µkI)−∗|λI −A|(A− µkI)

−1x(k)

= (x(k))∗∑

λ∈σ(A)\λ

|λ− λ||λ− µk|2

w(k)

λ=

∑

λ∈σ(A)\λ

|λ− λ||λ− µk|2

∥∥∥w

(k)

λ

∥∥∥

2

≤ maxλ∈σ(A)\λ

|λ− µk|−2∑

λ∈σ(A)\λ

|λ− λ|∥∥∥w

(k)

λ

∥∥∥

2

= maxλ∈σ(A)\λ

|λ− µk|−2 (x(k))∗|λI − A|x(k) = ǫk maxλ∈σ(A)\λ

|λ− µk|−2.

Insgesamt folgt (2.2).

2.6 Das QR-Verfahren

Das prominenteste Verfahren, um alle Eigenwerte und -vektoren einer Matrixzu bestimmen ist das QR-Verfahren. Wir können im Rahmen dieser Vorle-sung nur die grundlegende Idee vorstellen.

2.6.1 Motivation: simultane Potenzmethode

Betrachten wir zur Motivation noch einmal die Potenzmethode. Um alle Ei-genwerte und -vektoren einer Matrix zu bestimmen, können wir das Potenz-verfahren auf n Startvektoren x1, . . . , xn ∈ Cn gleichzeitig anwenden. In Ma-trixnotation beginnen wir also mit

(x1 . . . xn

)∈ Cn×n

84

2.6. DAS QR-VERFAHREN

und berechnen im Wesentlichen (d.h. bis auf Normierung)

Ak(x1 . . . xn

)=(Akx1 . . . Akxn

).

Im Allgemeinen wird jedoch an alle Startvektoren der Eigenvektor zum be-tragsgrößten Eigenwert λ1 beteiligt sein und alle Spalten Akxi werden gegeneinen Eigenvektor zu λ1 konvergieren.

Um das zu verhindern, orthonormalisieren wir die Vektoren in jedem Schritt.Wir nehmen also im ersten Schritt nicht

(Ax1 . . . Axn

)

sondern Vektoren (

x(1)1 . . . x

(1)n

)

mit den folgenden Eigenschaften:

• Für jedes k = 1, . . . , n die ersten k Vektoren jeweils den gleichen Raumaufspannen, also

span(x(1)1 , . . . , x

(1)k ) = span(Ax1, . . . , Axk)

• x(1)1 ,. . . , x(1)

n orthonormiert sind, also

(x(1)j )∗x

(1)k = δjk

Im nächsten Schritt bilden wir(

Ax(1)1 . . . Ax

(1)n

)

und orthonormalisieren diese zu Vektoren(

x(2)1 . . . x

(2)n

)

u.s.w.

Für den ersten Vektor entspricht dieses Verfahren gerade der Potenzmetho-de. Wir können also erwarten, dass er gegen einen Eigenvektor v1 zum be-tragsgrößten Eigenwert λ1 konvergiert. Beim zweiten Vektor wird bei diesemVerfahren die Potenzmethode angewandt und zusätzlich durch die Orthogo-nalisierung der Anteil in Richtung des erste Vektors eliminiert. Wir könnenalso erwarten, dass der zweite Vektor gegen einen Eigenvektor zum größtenEigenwert der Matrix A eingeschränkt auf den Unterraum

span(v1)⊥ = v : v∗1v = 0

85


konvergiert. Dies wird (zumindest wenn die Eigenvektoren senkrecht aufein-ander stehen) gerade der zweitbetragsgrößte Eigenwert sein. Enstprechenderwarten wir Konvergenz aller n-Vektoren gegen Eigenvektoren zu allen nEigenwerten. Die Eigenwerte erhalten wir dann als Diagonaleinträge von

(

x(k)1 . . . x

(k)n

)∗

A(

x(k)1 . . . x

(k)n

)

.

Für die Orhogonalisierung bietet sich das aus der Linearen Algebra bekann-te Gram-Schmidtsche Orthogonalisierungsverfahren an. Wir kennen jedochbereits eine Implementierung eines Orthogonalisierungsverfahren aus der Nu-merik I. Sei

A = QR

die QR Zerlegung von A. Da Q unitär ist, sind die Spalten von Q ortho-normiert. Da R obere rechte Dreiecksgestalt hat. ergibt sich die k-te Spaltevon A durch Kombination von k Spalten von Q, d.h. die jeweils ersten kSpalten von A und Q spannen den selben Raum auf. Q enthält also eineOrthonormalisierung der Spaltenvektoren von A.

Es erscheint natürlich, als Startvektoren die Einheitsvektoren zu wählen, also

(

x(0)1 . . . x

(0)n

)

=

1 0 0

0. . . 0

0 0 1

.

Dann ist im ersten Schritt keine Orthonormalisierung nötig und es ergibt(durch Anwendung von A) im den ersten Schritt immer A, so dass wir auchgleich mit

A0 :=(

x(0)1 . . . x

(0)n

)

= A.

starten können. Insgesamt erhalten wir also den folgenden Algorithmus

Algorithm 3 Simultane Potenzmethode

Gegeben Matrix A ∈ Cn×n. Setze A0 = A.for k = 0, 1, . . . do

Berechne QR-Zerlegung Ak = QkRk

Ak+1 := AQk

end forreturn Qk approximiert spaltenweise Eigenvektoren von A. Die Diago-naleinträge von Q∗

kAQk approximieren die Eigenwerte.

86


2.6.2 Das QR-Verfahren

Durch Umformulierung der simultanen Potenzmethode in Algorithmus 3 er-halten wir das QR-Verfahren in seiner einfachsten Form:

Algorithm 4 QR-VerfahrenGegeben Matrix A ∈ Cn×n. Setze A0 = A.for k = 0, 1, . . . do

Berechne QR-Zerlegung Ak = QkRk

Ak+1 := RkQk

end forreturn Die Diagonaleinträge von Ak approximieren die Eigenwerte.

Lemma 2.18Sei A ∈ Cn×n. Ak, Qk und Rk seien durch Algo. 4 definiert. Dann gilt:

(a) Für alle k ≥ 0 ist

Ak+1 = Q∗kAkQk = (Q0 . . . Qk)

∗A(Q0 . . . Qk).

(b) Für alle k ≥ 1 ist

Ak = (Q0 . . . Qk−1)(Rk−1 . . . R0).

(c) Die Folgen Qk := Q0 . . . Qk, Rk := Rk und Ak := Qk−1Ak erfüllen dieIterationsvorschrift aus Algorithmus 3

Ak = QkRk, Ak+1 = AQk ∀k ≥ 0.

Beweis: (a) ist trivial.

(b) ist offenbar richtig für k = 1. Gilt (b) für ein k ≥ 1, dann ist

Ak+1 = AAk = A(Q0 . . . Qk−1)(Rk−1 . . . R0).

Wegen (a) gilt A(Q0 . . . Qk−1) = (Q0 . . . Qk−1)Ak und zusammen mitAk = QkRk folgt

Ak+1 = (Q0 . . . Qk−1)QkRk(Rk−1 . . . R0).

87


(c) Aus der Definition von Ak und Qk sowie aus Ak = QkRk folgt

Ak = Qk−1Ak = Q0 . . . Qk−1QkRk = QkRk

Außerdem folgt aus (a)

Ak+1 = QkAk+1 = Q0 . . . Qk(Q0 . . . Qk)∗A(Q0 . . . Qk) = AQk.

Im Rahmen dieser Vorlesung beweisen wir die Konvergenz des QR-Verfahrensnur für den (durch die Motivation durch das simultane Potenzfahren nahe-liegenden) Spezialfall, dass A diagonalisierbar ist und paarweise betragsver-schiedene, einfache Eigenwerte besitzt. Außerdem benötigen eine Zusatzvor-aussetzung, die sicherstellt, dass die Eigenvektoren von A in geeigneter WeiseAnteile in Richtung der Einheitsvektoren besitzen.

Satz 2.19A ∈ C

n×n sei diagonalisierbar und besitze paarweise betragsverschiedene,einfache Eigenwerte, also

A = XΛX−1, Λ =

λ1 0 0

0. . . 0

0 0 λn

mit X ∈ Cn×n und dem Betrag nach sortierten

|λ1| > . . . > |λn| > 0.

Zusätzlich besitze X−1 eine LR-Zerlegung X−1 = LU . 1

Dann konvergieren die Diagonaleinträge der durch Algorithmus 4 erzeugtenMatrixfolge Ak mit (mindestens) linearer Geschwindigkeit gegen die Eigen-werte von A.

Beweis: (a) Aufgrund der Existenz der LR-Zerlegung ist

Ak = (XΛX−1)k = XΛkX−1 = XΛkLU = XΛkLΛ−k︸︷︷︸

=:Xk

ΛkU.

1Jede invertierbare Matrix besitzt nach Zeilen- oder Spaltenvertauschungen eine LR-Zerlegung. Der nicht-trivale Kern unserer Zusatzvoraussetzung besteht darin, dass eineLR-Zerlegung existiert für die Matrix X−1, in der die Eigenvektoren entsprechend derGröße der Eigenwerte sortiert sind. Zeilen-/Spaltenvertauschungen würden diese Sortie-rung zerstören!

88


Die Einträge von ΛkLΛ−k sind

(ΛkLΛ−k

)

ij= λk

i lijλ−kj =

0 für i < j,1 für i = j,lijλ

ki /λ

kj für i > j,

Aufgrund der Sortierung der Eigenwerte ist |λi/λj | < 1 für alle i > j undes folgt dass

ΛkLΛ−k = I +O(qk)

wobei I die n×n-Einheitsmatrix ist und q := max|λi/λj| : i > j < 1.

Damit folgt auch

Xk = XΛkLΛ−k = X + Ek mit Ek = O(qk). (2.3)

(b) Man zeigt leicht, dass die Inverse einer rechten oberen Dreiecksmatrixwieder eine rechte ober Dreiecksmatrix ist. Offenbar ist auch das Pro-dukt rechter oberer Dreiecksmatrizen wieder eine rechte obere Dreiecks-matrix, die Diagonaleinträge des Produkts sind das Produkt der Diago-naleinträge und die Diagonaleinträge der Inversen sind die Inversen derDiagonaleinträge.

Sind Q1R1 = Q2R2 zwei QR-Zerlegungen einer Matrix, dann existierteine unitäre Diagonalmatrix S mit Q1 = Q2S

∗, R1 = SR2.

(c) Nach Lemma 2.18 ist

Ak = (Q0 . . . Qk−1)(Rk−1 . . . R0) =: Qk−1Rk−1

eine QR-Zerlegung von Ak.

Sei Xk = PkUk eine QR-Zerlegung von Xk. Dann ist

Ak = XkΛkU = Pk(UkΛ

kU)

eine weitere QR-Zerlegung von Ak.

Es exisitert also eine unitäre Diagonalmatrix Sk mit

Qk−1 = PkS∗k und Rk−1 = SkUkΛ

kU.

Nun gilt

Qk = (Q0 . . . Qk−1)−1(Q0 . . . Qk) = Q−1

k−1Qk = (PkS∗k)

−1Pk+1S∗k+1

= SkP−1k Pk+1S

∗k+1.

Rk = (RkRk−1 . . . R0)(Rk−1 . . . R0)−1 = RkR−1

k−1

= Sk+1Uk+1Λk+1U

(SkUkΛ

kU)−1

= Sk+1Uk+1ΛU−1k S∗

k .

89


und damit

Ak = QkRk = SkP−1k Pk+1Uk+1ΛU

−1k S∗

k

= Sk UkU−1k P−1

k Pk+1Uk+1ΛU−1k S∗

k .

Wir verwenden noch Xk = PkUk und Xk+1 = Pk+1Uk+1 und erhalten

Ak = SkUkX−1k Xk+1ΛU

−1k S∗

k . (2.4)

(d) Nach (2.3) giltXk = X + Ek mit Ek = O(qk).

Damit folgt, dass

X−1k = (X + Ek)

−1 = X−1 +O(qk)

und damit

X−1k Xk+1 = (X + Ek)

−1(X + Ek+1) = I + Fk mit Fk = O(qk).

Einsetzen in (2.4) liefert

Ak = SkUkΛU−1k S∗

k + SkUkFkΛU−1k S∗

k .

Betrachte zuerst den ersten Summanden: SkUkΛU−1k S∗

k ist eine rechteobere Dreiecksmatrix, deren Diagonaleinträge jeweils das Produkt derentsprechenden Diagonaleinträge der Faktoren sind. Die Diagonaleinträ-ge von Sk und S∗

k = S−1k sowie die von Uk und U−1

k sind jeweils zueinanderinvers. Die Diagonaleinträge der ersten Summanden sind also gerade dievon Λ, also die Eigenwerte von A.

Wir müssen also nur noch zeigen, dass der zweite Summand gegen Nullkonvergiert. Es ist

∥∥SkUkFkΛU

−1k S∗

k

∥∥ ≤ ‖Sk‖ ‖Uk‖ ‖Fk‖ ‖Λ‖

∥∥U−1

k

∥∥ ‖S∗

k‖ .

Da Uk = P ∗kXk und Pk und Sk unitär sind folgt∥∥SkUkFkΛU

−1k S∗

k

∥∥ ≤ ‖Xk‖

∥∥X−

k 1∥∥ ‖Λ‖ ‖Fk‖ .

Da Xk → X und X−1k → X−1 ist ‖Xk‖

∥∥X−1

k

∥∥ beschränkt. Es folgt

∥∥SkUkFkΛU

−1k S∗

k

∥∥ = O(qk)


90


Bemerkung 2.20(a) Nach Lemma 2.18 und Satz 2.19 konvergiert

Ak = (Q0 . . . Qk−1)∗A(Q0 . . . Qk−1)

für große k gegen eine rechte obere Dreiecksmatrix. Die Eigenvektorenvon A ergeben sich aus den Eigenvektoren der Dreiecksmatrix durch dieunitäre Basistransformation Q0 . . . Qk−1.

(b) Ähnlich wie in Abschnitt 2.5 lässt sich die Konvergenzgeschwindigkeiterhöhen, indem A durch A− µkI mit geeigneten Shifts µk ersetzt wird.

(c) Algorithmus 4 benötigt in jedem Schritt eine QR-Zerlegung. Für eine ef-fiziente Implementierung bringt man A zuerst in eine spezielle Form, fürdie sich QR-Zerlegungen schnell berechnen lassen und die unter Algo-rithmus 4 erhalten bleibt, siehe z.B. [Hanke].

91


92

Anhang A

Hilfsmittel aus der Analysis

A.1 Der mehrdimensionale Mittelwertsatz

Alle im Folgenden vorkommenenden Funktionen werden als stetig differen-zierbar vorausgesetzt.

Für eine Funktion f : R → R und a, b ∈ R, a < b gilt der eindimensionaleMittelwertsatz

∃ξ ∈ [a, b] :f(b)− f(a)

b− a= f ′(ξ)

Die in der Numerik am häufigsten verwendete Folgerung daraus ist

|f(b)− f(a)| ≤ |b− a| supξ∈[a,b]

|f ′(ξ)|.

Die Änderung von f zwischen a und b lässt sich also durch Kenntnis derAbleitung an allen Mittelwerten zwischen a und b abschätzen.

Der eigentliche Mittelwertsatz lässt sich nicht ohne weiteres ins Mehrdimen-sionale übertragen, die Folgerung jedoch schon. Das folgende Lemma wirddeshalb oft als Mittelwertsatz im Mehrdimensionalen bezeichnet.

Lemma A.1Für F : Rn → R

m und x, y ∈ Rn gilt

F (x)−F (y) =

∫ 1

0

F ′(x+ t(y−x))(y−x) dt =

∫ 1

0

F ′(x+ t(y−x)) dt(y−x),

also insbesondere

‖F (x)− F (y)‖ ≤ ‖y − x‖ supt∈[0,1]

‖F ′(x+ t(y − x))‖ .

93

ANHANG A. HILFSMITTEL AUS DER ANALYSIS

Beweis: (a) Wir erinnern zunächst an die Kettenregel im Mehrdimensiona-len. Sind

g : Rs → Rn, h : Rn → R

m

zwei Funktionen und ist

f : Rs → R

m, f(x) := h(g(x))

ihre Komposition, so gilt für die zugehörigen Jacobi Matrizen

f ′(x) =

∂f1(x)∂x1

. . . ∂f1(x)∂xs

......

∂fm(x)∂x1

. . . ∂fm(x)∂xs

∈ Rm×s,

g′(x) =

∂g1(x)∂x1

. . . ∂g1(x)∂xs

......

∂gn(x)∂x1

. . . ∂gn(x)∂xs

∈ Rn×s

h′(y) =

∂h1(x)∂y1

. . . ∂h1(x)∂yn

......

∂hn(x)∂y1

. . . ∂hm(x)∂yn

∈ Rm×n

die Kettenregelf ′(x) = h′(g(x))g′(x).

(b) Außerdem benötigen wir den Hauptsatz der Differential- und Integral-rechnung. Für f : R→ R und x, y ∈ R gilt

f(y)− f(x) =

∫ y

x

f ′(t) dt.

Für eine vektorwertige Funktion

f : R→ R

m, f(x) =

f1(x)...

fm(x)

gilt diese Aussage ebenfalls, da sie komponentenweise gilt:

f(y)− f(x) =

f1(y)− f1(x)...

fm(y)− fm(x)

=

∫ y

xf ′1(t) dt...

∫ y

xf ′m(t) dt

=

∫ y

x

f ′(t) dt.

94

A.1. DER MEHRDIMENSIONALE MITTELWERTSATZ

Nun können wir den Mittelwertsatz im Mehrdimensionalen beweisen. SeiF : Rn → R

m und x, y ∈ Rn. Wir definieren

f : R→ R

m, f(t) := F (x+ t(y − x))

Offenbar ist f(t) = F (g(t)), wobei

g : R→ R

n, g(t) := x+ t(y − x)

und es ist

g′(t) =

∂g1∂t...

∂gn∂t

= y − x.

Aus der mehrdimensionalen Kettenregel (siehe (a)) folgt

f ′(t) = F ′(g(t))g′(t) = F ′(x+ t(y − x))(y − x)

und mit dem Hauptsatz der Differential- und Integralrechnung (siehe (b))erhalten wir

F (y)− F (x) = f(1)− f(0) =

∫ 1

0

f ′(t) dt =

∫ 1

0

F ′(x+ t(y − x))(y − x) dt


95

ANHANG A. HILFSMITTEL AUS DER ANALYSIS

96

Literaturverzeichnis

[FulfordBroadbridge] G. R. Fulford, P. Broadbridge: Industrial Mathematics:Case Studies in the Diffusion of Heat and Matter, Australian Mathe-matical Society Lecture Series 16, Cambridge University Press, Cam-bridge, 2002.

[Hanke] M. Hanke-Bourgeois: Grundlagen der Numerischen Mathematik unddes Wissenschaftlichen Rechnens, Teubner Verlag, Wiesbaden, 2009.

[HairerNorsettWanner] E Hairer, SP Nørsett, G Wanner: Solving ordinarydifferential equations I. Nonstiff problems, Springer, 1987.

97

Numerische Mathematik 2 - Universität Stuttgartharrach/lehre/Numerik2.pdf · Numerische Mathematik...

Documents

Transcript of Numerische Mathematik 2 - Universität Stuttgartharrach/lehre/Numerik2.pdf · Numerische Mathematik...