Numerik II - math.uni-bielefeld.demstephan/files/scripts/NumerikII.pdf · 1 Anfangswertprobleme...

Numerik II

Prof. C. Rohde∗

Version: 11. Oktober 2005

∗Eine Mitschrift von Matthias StephanEmail: [email protected]

Inhaltsverzeichnis

1 Anfangswertprobleme fur gewohnliche Differentialgleichungen 11.1 Grundlegende Definitionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Beispiele und Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.3 Analytische Grundlagen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.4 Einschrittverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.5 Explizite Runge-Kutta-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . 201.6 Implizite Runge-Kutta-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . 221.7 Runge-Kutta-Verfahren und Quadraturformeln . . . . . . . . . . . . . . . . . . . 251.8 Mehrschrittverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 291.9 Adaptive Schrittweitensteuerung . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

2 Randwertprobleme fur gewohnliche Differentialgleichungen 372.1 Beispiele und Definitionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 372.2 Lineare Randwertaufgaben/Sturm-Liouville Problem . . . . . . . . . . . . . . . . 402.3 Das klassische Finite Differenzen Verfahren . . . . . . . . . . . . . . . . . . . . . 412.4 Matrizenalgebra und Stabilitatssatz . . . . . . . . . . . . . . . . . . . . . . . . . 44

3 Die Finite-Elemente Methode in 1D 513.1 Das Sturm-Liouville Problem als Variationsproblem . . . . . . . . . . . . . . . . 513.2 Existenz (und Eindeutigkeit) einer schwachen Losung . . . . . . . . . . . . . . . 533.3 Sobolevraume in 1D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 563.4 Existenz einer schwachen Losung in (3.1),(3.2) . . . . . . . . . . . . . . . . . . . 593.5 Das Ritz-Galerkin-Verfahren als Motivation der FDM . . . . . . . . . . . . . . . 613.6 Der Teilraum Vh und das Finite-Elemente-Verfahren . . . . . . . . . . . . . . . . 633.7 Konvergenz der FEM in 1D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 653.8 Abschließende Bemerkungen zur FEM . . . . . . . . . . . . . . . . . . . . . . . . 68

4 Iterative Losungsverfahren fur lineare Gleichungssysteme 69

i

Inhaltsverzeichnis

ii

1 Anfangswertprobleme fur gewohnlicheDifferentialgleichungen

1.1 Grundlegende Definitionen

Fur m ∈ N sei der Zustandsraum U ⊆ Rm offen und f ∈ C(R × U,Rm). Weiter wahle t0 ∈R, u0 ∈ Rm und ein Intervall I ⊆ R mit t0 ∈ I.Finde ein u = (u1, . . . , um)T ∈ C1(I, U) mit

u′(t) = f(t, u(t)), (t ∈ I) (1.1)

u(t0) = u0. (1.2)

Gleichung (1.1) heisst (m-dimensionales) System gewohnlicher DGL (erster Ordnung).Die Gleichungen (1.1),(1.2) heissen Anfangswertproblem zur Anfangsbedingung (1.2).

Definition 1.1. Falls eine Funktion u ∈ C1(I, U) existiert, die das AWP (1.1),(1.2) erfullt,heisst u (klassische) Losung in I.

Bemerkung 1.2 (spezielle Typen).

i. Die Gleichung (1.1) heisst linear, falls sie von der Form

u′(t) = A(t)u(t) + b(t) (t ∈ R) (1.3)

ist. Dabei sei A ∈ C0(I,Rm×m), b ∈ C0(I,Rm).

ii. Falls in (1.3) b ≡ 0, heisst (1.3) homogen, sonst inhomogen oder affin.

iii. Die Gleichung (1.1) heisst autonom, falls f = f(u), d.h. nicht explizit von t abhangt.

1.2 Beispiele und Motivation

In diesem Unterkapitel soll gezeigt werden, weshalb numerische Methoden zu (approximativen)Losungen gewohnlicher DGL unerlasslich sind.

Beispiel 1.3 (Verbreitung eines Geruchts). Sei eine menschliche Population der Große N ∈ Ngegeben. Gesucht ist die Anzahl der Menschen I(t), die zum Zeitpunkt t ≥ 0 eine bestimmteNachricht kennen.Annahmen:

i. Zum Zeitpunkt t = 0 kennen I0 Menschen die Nachricht.

ii. Die Nachricht wird nur durch Gesprache unter vier Augen weitergegeben.(Nachricht ≈ Gerucht)

iii. Fur ∆t, k > 0 hat jeder Informierte in der Zeitspanne ∆t genau k∆t ”Vier-Augen-Kontakte“.

1

1 Anfangswertprobleme fur gewohnliche Differentialgleichungen

80

100

40

0

60

20

t

20151050

Abbildung 1.1: Verbreitungskurve fur N = 100 und k = 1

iv. I sei eine stetig differenzierbare Funktion.

Falls wir mit q(t) ∈ [0, 1] den Anteil der Nichtinformierten bezeichnen, gilt:

q(t)N = N − I(t) (t ≥ 0) (?)

Falls ∆t klein ist, kann man die Anzahl der Kontakte eines Informierten mit einem Nichtinfor-mierten in der Zeitspranne t bis t+ ∆t gut durch q(t)k∆t approximieren.Die Zunahme von I im Zeitraum ∆t ist dann gegeben durch

I(t+ ∆t)− I(t) = q(t)k∆tI(t).

Mit Gleichung (?) folgtI(t+ ∆t)− I(t)

∆t=N − I(t)

NkI(t)

mit (i),(iv) erhalten wir fur ∆t→ 0 das gewunschte AWP

I ′(t) = kN − I(t)

NI(t)

I(0) = I0

(??)

Mit Trennung der Variablen ergibt sich als Losung

I(t) =N

1 +(

NI0− 1)e−kt

.

Fur t→∞ erfahren alle Populationsmitglieder das Gerucht.Bemerkung: Die DGL (??) heisst auch logistische Gleichung.

In Beispiel 1.3 konnte eine Losung explizit bestimmt werden. Das nachste Beispiel zeigt, dassexakte Losungen nicht immer hilfreich sind.

Beispiel 1.4 (AWP mit komplizierter expliziten Losung). Betrachte das lineare homogene AWP

u′(t) = (1 + 2t)u(t) + 1 (t ∈ R)

2

1.2 Beispiele und Motivation

u(0) = 0

Die Funktion u ∈ C1(R) mit

u(t) = et(1+t)

∫ t

0e−s(1+s)ds (?)

ist eine (sogar die eindeutige) Losung des AWP.Die konkrete Bestimmung von Funktionswerten u(t) fuhrt auf ein Quadraturproblem, dessennumerische Durchfuhrung komplexer als die numerische Integration von (?) ist. Der Integrandvon (?) besitzt keine Stammfunktion, die sich durch Elementarfunktionen ausdrucken lasst.

Beispiel 1.5 (Bewegungsgleichungen). Fur N ∈ N betrachten wir die 3D-Bewegung von Teil-chen pi(i = 1, . . . , N) mit Massen mi > 0 (i = 1, . . . , N).Eine Anwendung ware zum Beispiel die Dynamik von Planeten, Molekul- und Atomwolken odergradulare Medien (Sand, Getreide, etc.).Wir schranken uns auf den Fall von Bewegungen ein, die unter dem Einfluss eines Potential-feldes V : R3N → R stattfinden. Falls das Potentialfeld durch ein Gravitationsfeld gegeben ist(Planetendynamik), gilt zum Beispiel

V (x1, . . . , xN ) = −∑i<j

i,j∈1,...,N

mimj

‖xi − xj‖2. (?)

Dabei sei xi ∈ R3 die Position des Teilchens pi, i = 1, . . . , N .Aus dem Newtonschen Kraftgesetz erhalt man die 3N Bewegungsgleichungen

mix′′i (t) = −∂V

∂xi(x1(t), . . . , xN (t)) (i = 1, . . . , N),

wobei der Vektor ∂V∂xi

durch

∂V

∂xi(x1, . . . , xN ) =

(∂V

∂xi1(x1, . . . , xN ), . . . ,

∂V

∂xi3(x1, . . . , xN )

)T

gegeben ist.

∂

∂xi1V (x1, . . . , xN )

(?)= −

N∑j=1j 6=i

∂

∂xi1

mimj√(xi1 − xj1)2 + . . .+ (xi3 − xj3)2

= −N∑

j=1j 6=i

mimj0− 2(xi1 − xj1)2−1

√(xi1 − xj1)2 + . . .+ (xi3 − xj3)2

−1√(xi1 − xj1)2 + . . .+ (xi3 − xj3)2

2

Wir erhalten im Fall (?) die Gleichungen

mix′′i (t) = −mi

N∑j=1j 6=i

mj

‖xi(t)− xj(t)‖32

(xi − xj) (i = 1, . . . , N)

Zum Zeitpunkt t = 0 sind die Anfangspositionen und Anfangsgeschwindigkeiten gegeben:

xi(0) = x0i ∈ Rn, x′i(0) = v0

i ∈ Rn (i = 1, . . . , N).

Wir schreiben dieses AWP 2.Ordnung auf ein AWP 1.Ordnung um:

u = u(t) : (u1(t), . . . , u2N (t)) := (x1(t), x′1(t), . . . , xN (t), x′N (t)) ∈ R6N .

3


Insgesamt erhalten wir das 6N-dimensionale AWP

u′(t) =

u2(t)

− 1m1

∂V∂x1

(u1(t), u3(t), . . . , u2N−1(t))...

u2N (t)− 1

mN

∂V∂xN

(u1(t), u3(t), . . . , u2N−1(t))

und

u(0) = (x01, v

01, . . . , x

0N , v

0N ).

In den Anwendungen ist in der Regel N 1. In diesem Fall ist nicht zu erwarten, dass eineexplizite Losung bestimmt werden kann.Bemerkung: Jedes AWP hoherer Ordnung kann auf ein AWP 1.Ordnung (aber hoherer Di-mension) umgeschrieben werden. Da fast alle Softwarepakete auf AWP 1.Ordnung zugeschnittensind, betrachtet man in der Numerik nur AWP 1.Ordnung.

Beispiel 1.6 (chemische Reaktionen). Wir betrachten zunachst fur Spezies A,B,C und Reak-tionsrate k > 0 die allgemeine bimolekulare Reaktion

A+Bk→ C.

Die zugehorigen Konzentrationen (in[mol/Vol]) seien mit a, b, c bezeichnet. Falls die Speziesinitial raumlich homogen verteilt sind, kann man folgenden Ansatz fur die zeitliche Entwicklungvon a = a(t), b = b(t), c = c(t) machen:

a′(t) = −ka(t)b(t) a(0) = a0

b′(t) = −ka(t)b(t) b(0) = b0c′(t) = ka(t)b(t) c(0) = c0

(?)

a0, b0, c0 ≥ 0 sind gegebene Anfangskonzentrationen.Dieser allgemeine Modellierungsansatz soll jetzt auf die komplexe CO-Oxidation an Pt ange-wandt werden (Katalysatorreaktion).

O2 + 2Pt k1→ 2PtO

CO + Ptk2

k−2

PtCO

PtCO + PtOk3→ 2Pt+ CO2

CO + Ptk4

k−4

[PtCO]

Dabei seien k1, k2, k−2, k3, k4, k−4 > 0 Reaktionsraten und [PtCO] bezeichne eine chemisch in-aktive Version von PtCO. Als Bezeichnungen fur die Konzentrationen fuhren wir ein:

PtO : x, P t : z, P tCO : y, [PtCO] : s

Die Konzentration von O2 und CO werden als konstant angesehen, da diese Stoffe standignachgefuhrt werden. O.B.d.A. seien diese ≡ 1. Wir erhalten mit (?) (unter Vernachlassigungder freien Variablen t)

x′ =2k1z2 − k3xy

z′ =− 2k1z2 − k2z + k−2y + 2k3xy − k4z + k−4s

y′ =k2z − k−2y − k3xy

s′ =k4z − k−4s

(1.4)

4

1.3 Analytische Grundlagen

y

20

10

0

-10

-20

x

1,510,50-0,5-1-1,5

1

0

0,5

10

-0,5

x

5-5 0-10

-1

Abbildung 1.2: Verlauf von tan(t) und tanh(t)

Dieses System kann nicht explizit gelost werden. Wir bemerken noch, dass aus (1.4) fur eineLosung (x, y, z, s) durch Addition der Gleichungen

(x′, y′, z′, s′)(t) = 0 ∀t ≥ 0

gilt. Und daraus folgt die Erhaltungseigenschaft

(x+ y + z + s)(t) = const., (t ≥ 0).


In diesem Kapitel werden einige Resultate zur Analysis gewohnlicher AWP zusammengetragen,soweit diese fur die Numerik von Bedeutung sind.

Beispiel 1.7. Sei m = 1, U = R, u0 = 0. Fur a ∈ R betrachten wir (1.1),(1.2) mit

f(t, u) = u2 + a

i. a < 0 :u(t) = −

√−a · tanh(

√−at)

ist eine Losung in I = R.

ii. a ≥ 0 :u(t) =

√a · tan(

√at)

ist eine Losung in(−π2√

a, π

2√

a

).

Falls (1.1),(1.2) eine Losung in I = R besitzt, sagt man, dass (1.1),(1.2) global losbar ist. ImAllgemeinen ist (1.1),(1.2) nicht global losbar (siehe Beispiel 1.7). Aber es gilt der folgende Satz:

5


Satz 1.8 (Picard-Lindelof). Sei a > 0 und b > 0 so gewahlt, dass die Bedingung Bb(u0) ⊂ U :=R

m gilt. Es existiere eine Konstante L > 0, so dass fur alle t ∈ [t0− a, t0 + a] und alle Vektorenw, w ∈ Bb(u0) die Abschatzung

|f(t, w)− f(t, w)| ≤ L|w − w|

gilt (Lipschitz-Stetigkeit in [t0 − a, t0 + a]×Bb(u0)). Dann gilt

i. Es existiert ein α ∈ (0, a], so dass (1.1),(1.2) eine Losung u ∈ C1((t0 − α, t0 + α),Rm)besitzt.

ii. Es existiert genau eine Losung von (1.1),(1.2) in (t0 − α, t0 + α).

Fur α in (i),(ii) gilt

α = min(a,b

A

)A := max

t∈[t0−α,t0+α]w∈Bb(u0)

|f(t, w)|

Fur den Beweis wird der Banachsche Fixpunktsatz benotigt:Sei D 6= ∅ eine abgeschlosssene Teilmenge eines Banachraumes (X, ‖.‖). Fur T : D → X gelte

i. T (D) ⊆ D

ii. ∃L ≥ 0 : ‖Tv − T v‖ ≤ L‖v − v‖ ∀v, v ∈ D.

Falls L ∈ (0, 1) gilt, besitzt die Fixpunktgleichung T (u) = u genau eine Losung u ∈ D.

Beweis: von Satz 1.8

i. Sei o.B.d.A. I := [t0, t0 + a] (die andere Seite analog).Wir betrachten das Fixpunktproblem

u(t) = (T (u))(t), (t ∈ I) (?)

fur die ungekannte Funktion u ∈ C0(I,Rm). Dabei sei der Operator T durch

T : C0(I) → C0(I),

u 7→ T (u) = (T (u))(t) := u0 +∫ t

t0

f(s, u(s))ds

definiert. Falls eine Losung von (?) existiert, gilt sogar u ∈ C1(I,Rm). Einsetzen der Losungin (1.1),(1.2) ergibt gerade, dass es auch eine Losung von (1.1),(1.2) ist. Es genugt also,(?) zu betrachten.

a) Zunachst nehmen wir an, dass f die globale Lipschitz-Bedingung

|f(t, w)− f(t, w)| ≤ L|w − w| (??)

fur t ∈ I und w, w ∈ Rm erfullt.Sei β := 2L. Der Raum (C1(I,Rm), ‖ · ‖β) mit

‖v‖β := maxt∈I

|v(t)|e−βt

ist ein Banachraum. Wir wenden das Kontraktionsprinzip an mit

X = (C0(I,Rm), ‖.‖β), D = X.

6


Es gilt fur v, v ∈ D und t ∈ I:

|T (v)(t)− T (v)(t)| =∣∣∣∣∫ t

t0

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

(??)

≤∫ t

t0

L|v(s)− v′(s)|e−βseβsds ≤ L‖v − v‖β

∫ t

t0

eβsds

= L‖v − v‖β ·1β

(eβt − eβt0

) β=2L≤ 1

2‖v − v‖βe

βt

⇒ ‖T (v)− T (v)‖β ≤12‖v − v‖β

Somit ist T eine Kontraktion in D. Also besitzt T einen Fixpunkt u ∈ C0(I,Rm) und(1.1),(1.2) eine Losung u ∈ C0(I,Rm).

b) Wir beweisen jetzt die Existenz ohne die Einschrankung (??). Dazu betrachten wirdas abgewandelte AWP:

u′(t) = f(t, u(t)) (t ∈ I)u′(t0) = u0

(? ? ?)

Dabei sei f ∈ C0(I ×Rm,Rm) eine beliebige Funktion mit

1. f(t, w) = f(t, w) fur t ∈ I und w ∈ Bb(u0)2. f erfulle (??) fur L > 0.3. maxt∈I,w∈Rmf(t, w) ≤ A

Nach (a) und (2) existiert eine Losung u von (? ? ?). Wegen (1) ist u eine Losung von(1.1),(1.2) in (t0, t0 +α) mit α > 0, falls u(t) ∈ Bb(u0)∀t ∈ [t0, t0 +α). Es ist fur t ∈ I

|u′(t)| ≤ |f(t, u(t))|(3)

≤ A.

Also kann man u(t) ∈ Bb(u0) fur α = min(a, b

A

)garantieren.

Teil (i) ist bewiesen. (Das Intervall [t0 − a, t0] wird analog behandelt.)

ii. Eindeutigkeit: ∫ t

t0

u′(t) =∫ t

t0

f(s, u(s))ds

ist das Fixpunktproblem. Dies hatte nur eine Losung. Somit folgt die Eindeutigkeit desAWP.

In einigen Fallen kann gezeigt werden, dass die durch Satz 1.8 gegebene lokale Losung auf I = Rfortgesetzt werden kann.

Lemma 1.9 (Apriori-Abschatzung). Es gelte folgende Implikation:Fur alle globale Losungen v ∈ C1(R, U) von (1.1),(1.2) gilt:

∃c > 0 : ‖v‖L∞(R) ≤ c.

Fur f ∈ C1(R× U,Rm) ist (1.1),(1.2) global losbar.

Beweis: Ubung.

7


Lemma 1.10. Die Funktion f sei global Lipschitz-stetig. Dann ist (1.1),(1.2) global losbar.

Beweis: Das ist in Beweis von (1.8) schon gezeigt worden.

Aus dem letzten Lemma (1.10) folgt fur affine AWP:

Korollar 1.11. Sei A ∈ C0(R,Rm×m) ∩ L∞(R,Rm×m). Dann besitzt das lineare AWP

u′(t) = A(t)u(t) + b(t) (t ∈ R)u(0) = u0

genau eine Losung u ∈ C1(I,Rm).

Die Losungen linearer Anfangswertprobleme lassen sich genauer charakterisieren.

Satz 1.12. Sei A ∈ C0(R,Rm×m) und b ∈ C0(R,Rm). Dann gilt:

i. Die Losungen der homogenen Differentialgleichung

u′(t) = A(t)u(t) (t ∈ R) (1.5)

bilden einen m - dimensionalen Unterraum V in C1(R,Rm). Eine Basis von V ist durchdie Losungen u1, . . . , um ∈ C1(R,Rm) mit der Eigenschaft

ui(0) = ei (i = 1, . . . ,m) (i-ter Einheitsvektor)

bestimmt. Die matrixwertige Funktion Y0 = Y0(t) := (u1| . . . |um) heisst normierte Funda-mentalmatrix von (1.5).

ii. Sei Y ∈ C1(R,Rm×m) so gegeben, dass die Spaltenvektoren von Y Losungen von (1.5)sind, dann gilt fur t ∈ R:

Y (t) = Y0(t)Y (0) (1.6)

Weiter ist det(Y (t)) 6= 0 ∀t ∈ R genau dann, wenn det(Y (0)) 6= 0. Y heisst Fundamen-talmatrix zu (1.5).

iii. Die Losungen der inhomogenen Differentialgleichung

u′(t) = A(t)u(t) + b(t) (t ∈ R) (1.7)

bilden fur u ∈ C1(R,Rm) einen affinen Unterraum u+V von C1(R,Rm). Dabei ist u eine(spezielle) Losung von (1.7).Fur die Losung des Anfangswertproblems fur (1.7) gilt

u(t) = Y0(t)u0 +∫ t

0Y0(t)(Y0(s))−1b(s)ds.

iv. Sei A ∈ Rm×m. In dem autonomen Fall

u′(t) = Au(t)

ist die normierte Fundamentalmatrix durch

Y0(t) = etA :=∞∑

n=1

tnAn

n!(t ∈ R)

gegeben.

8


Beweis:

i. Wir zeigen, dass u1(t), . . . , um(t) ∈ Rm fur alle t ∈ R linear unabhangig sind (Existenzder u1, . . . , um nach Korollar (1.11) klar).Annahme: Es existiert c = (c1, . . . , cm)T ∈ Rm\0, τ ∈ R,

∑mi=1 ci · ui(τ) = 0. Dann ist

v =∑m

i=1 ciui eine Losung des AWPes

v′(t) = A(t)v(t), v(τ) Ann.= 0

Wegen der Eindeutigkeit der Losung v ≡ 0 gilt

0 = v(τ) = v(0) =m∑

i=1

ciui(0) =m∑

i=1

ciei 6= 0

E

Jede andere Losung des AWPes fur (1.5) lasst sich als Linearkombination der u1, . . . , um

darstellen. Also ist V ein m-dimensionaler Unterraum mit Basis u1, . . . , um.

ii. Die Funktionen aus der i-ten Spalte von Y und Y0 ·Y (0) losen beide das homogene System(1.5) nach Definition. Weiter gilt

Y (0) = Y0(0)Y (0)

wegen Y0(0) = I.Eindeutigkeit der Losung des AWPes fur (1.5) liefert (1.6).Weiter ist fur t ∈ R

det(Y (t))(1.6)= det(Y0(t)Y (0)) = det(Y0(t))︸︷︷︸

6=0(nach Bew. von (i))

·det(Y (0))

Also ist (ii) bewiesen.

iii. Die Aussage im affinen Unterraum ist klar. Zum Beweis der Formel betrachten wir zunachstdie Funktion u ∈ C1(R,Rm) mit

u(t) :=∫ t

0Y0(t)(Y0(s))−1b(s)ds.

Es gilt

u′(t) = Y ′0(t)

∫ t

0(Y0(s))−1b(s)ds+ Y0(t)(Y0(t))−1b(t)

= A(t)Y0(t)∫ t

0(Y0(s))−1b(s)ds+ b(t) = A(t)u(t) + b(t)

u ist eine spezielle Losung von (1.7) mit u(0) = 0.Wegen u(0) = Y0(0)u0 + u(0) = Y0(0)︸︷︷︸

=Id

u0 = u0 erfullt u die Anfangsbedingung und, da

Y0(t)u0 eine Losung des homogenen Systems ist, das AWP fur (1.7).

9


iv. Zunachst zeigen wir, dass etA fur t ∈ R existiert. Fur eine bliebige Matrixnorm |.| gilt:

|etA| Def.=

∣∣∣∣∣∞∑

n=1

tnAn

n!

∣∣∣∣∣ ≤∞∑

n=1

|t|n

n!|A|n = e|t|·|A| <∞

Es gilt fur den Ausdruck∞∑

n=1

n

n!tn−1An =

∞∑n=1

tn−1

(n− 1)!An−1A = A

∞∑n=0

tn

n!An = AetA

Mit dem selben Argument wie fur etA kann gezeigt werden, dass die Potenzreihe∞∑

n=1

n

n!tn−1An,

das ist gerade die gliedweise Ableitung von etA, existiert. Damit ist Y ′0(t) = AetA und

offenbar lost Y0 die Dgl. Y ′0 = AY0.

Wegen Y0(0) = Id ist Y0 tatsachlich die normierte Fundamentalmatrix.

1.4 Einschrittverfahren

Wir betrachten (1.1),(1.2) und nehmen an, dass (1.1),(1.2) genau eine Losung u ∈ C1(I,Rm)besitzt. Sei o.B.d.A. t0 = 0 und fur T > 0 das Intervall I durch I := [0, T ] gegeben.Um ein numerisches Verfahren zur Approximation von u zu konstruieren, sei fur N ∈ N derVektor

h := (h0, . . . , hN−1)T ∈ ([0, T ])N

mitN−1∑j=0

hj = T

gegeben.Wir konstruieren das Gitter Ih zu I durch

Ih := t0 := 0, t1, t2, . . . , tN

wobei fur j = 1, . . . , N − 1 die Beziehung tj := tj−1 + hj−1 gelte. (Also insbesondere tN = T .)Falls fur das Gitter

h1 = h2 = . . . = hN−1

gilt, sprechen wir von einem aquidistanten Gitter. In diesem Fall betrachten wir h als Skalar mith = hj (j = 1, . . . , N − 1).Als Gitterweite bezeichnen wir die Zahl

|h| := maxhj |j = 0, . . . , N − 1.

Ziel aller numerischen Verfahren ist die Bestimmung einer Gitterfunktion

uh : Ih → Rm

fur gegebenes Gitter Ih. Wir benutzen auch die Schreibweise

uj := uh(tj) fur j = 0, . . . , N.

10


Beispiel 1.13 (Eulersches Einschritt-Verfahren).Sei zunachst m = 1 und sei ein aquidistantes Gitter Ih fest vorgegeben. Setze uh(t0) = u0.Fur die exakte Losung u von (1.1),(1.2) gilt fur ein τ01 ∈ [t0, t1] nach Satz von Taylor

u(t1) = u(t0) + (t1 − t0)u′(τ01) = u0 + hf(τ01, u(τ01))

Wir approximieren u(t1) durch

u1 = uh(t1) := u0 + hf(t0, u(t0))

Um u(t2) zu approximieren, wenden wir dieselbe Idee, aber mit Startpunkt (t1, u1) an:

u2 = uh(t2) := uh(t1) + hf(t1, uh(t1)) = u1 + hf(t1, u1)

Das durch die Iterationsvorschrift

uj := uj−1 + hj−1f(tj−1, uj−1) (j = 1, . . . , N)

auf einem allgemeinen Gitter Ih festgelegte numerische Verfahren heisst explizites Eulerverfah-ren oder Polygonzugverfahren.

Das explizite Eulerverfahren gehort zur Klasse der Einschrittverfahren, die folgendermaßen de-finiert sind.

Definition 1.14 (Einschrittverfahren (ESV)).Sei ein Gitter Ih gegeben und φ ∈ C([0, T ]2 ×Rm,Rm). Dann heisst das Verfahren

uj := uj−1 + hj−1φ(hj−1, tj−1, uj−1) (j = 1, . . . , N)

explizites Einschrittverfahren und φ die zugehorige Inkrementfunktion.

Beispiel 1.15.

i. Zur Motivation des Eulerverfahrens hatten wir in (Bsp 1.13) die Ableitung u(τ01) gemaß

u′(τ01) ≈ f(t0, u0)

approximiert. Eine andere Wahl ware die Mittelpunktwahl

u′(τ01) ≈ f

(t0 +

h

2, u

(t0 +

h

2

)).

Zur Approximation von u(t0 + h

2

)betrachte

u

(t0 +

h

2

)≈ u(t0) +

h

2u′(t0)

(1.1)= u(t0) +

h

2f(t0, u(t0))

Im Falle eines aquidistanten Gitters ergibt sich die Iterationsvorschrift

uj = uj−1 + hf

(tj−1 +

h

2, uj−1 +

h

2f(tj−1, uj−1)

)(j = 1, . . . , N)

Mit der Wahl

φ(k, t, w) = f

(t+

k

2, w +

k

2f(t, w)

)erhalten wir ein explizites ESV, das man als modifiziertes Eulerverfahren oder Collatz-verfahren bezeichnet.

11


ii. (Fast) ohne Motivation sei noch das Verfahren von Heun angegeben, dass sich aus Inkre-mentfunktion

φ(k, t, w) =12(f(t, w) + f(t+ k,w + kf(t, w)))

ergibt. Dazu:f(tj , uj) ≈ f(tj−1 + hj−1, uj−1 + hj−1f(tj−1, uj−1))

Um die Qualitat der verschiedenen Einschrittverfahren ESV zu beurteilen, werden wir Fehler-abschatzungen herleiten.

Definition 1.16 (Fehler). Fur h ∈ (0, T ]N sei uh eine gegebene Gitterfunktion auf Ih.

i. Die Funktion eh : Ih → Rm mit eh := u

∣∣Ih−uh heisst globale Fehlerfunktion. Als globalen

Diskretisierungsfehler bezeichnen wir den Ausdruck

eh = max|eh(tj)||j = 0, . . . , N.

ii. Sei uh durch ein ESV mit Inkrementfunktion φ ∈ C0([0, T ]2×Rm,Rm) gegeben. Die Funk-tion εh : Ih → R

m mit

εh(tj) =1hj

(u(tj+1)− u(tj)− hjφ(hj , tj , u(tj))) j = 0, . . . , N

heisst lokale Fehlerfunktion. Als lokalen Diskretisierungsfehler oder Konsistenzfehler be-zeichnen wir den Ausdruck

εh = max|εh(tj)|j = 0, . . . , N.

Der wirklich wichtige Fehler ist eh. Der Konsistenzfehler spielt aber eine wichtige technischeRolle und soll zunachst untersucht werden.

Definition 1.17 (Konsistenz). Fur h ∈ (0, T ]N sei uh eine gegebene Gitterfunktion auf Ih, diedurch ein ESV erzeugt sei. Das ESV heisst konsistent zu (1.1) genau dann, wenn

lim|h|→0

εh = 0

gilt. Das ESV heisst konsistent zur Ordnung p ∈ N zu (1.1) genau dann, wenn

εh = O(|h|p)

gilt.

Satz 1.18. Fur h ∈ (0, T ]N sei ein Gitter Ih und ein ESV mit Inkrementfunktion φ ∈ C0([0, T ]2×R

m,Rm) gegeben.

i. Das ESV ist genau dann konsistent zu (1.1), wenn

φ(0, t, u(t)) = f(t, u(t))

fur t ∈ [0, T ] gilt.

ii. Fur p ∈ N sei f ∈ Cp(I ×Rm,Rm) und φ ∈ Cp([0, T ]2×Rm,Rm). Das ESV ist konsistentzur Ordnung p zu (1.1) genau dann, wenn fur i = 0, . . . , p− 1

di

dti[f(t, u(t))] = (i+ 1)

∂i

∂kiφ(0, t, u(t)) (1.8)

fur t ∈ [0, T ] gilt.

12


Bemerkung 1.19. Um die hochste Ableitung dp−1

dtp−1 [f(t, u(t))] in (1.8) im klassischen Sinn auf-fassen zu konnen, muss u ∈ Cp−1(I ,Rm) sein. Das ist wegen f ∈ Cp−2(I × Rm,Rm) gegeben.Es ist sogar u ∈ Cp+1(I ,Rm), was im Beweis benotigt wird.

Beweis:

i. ”⇒“: Wahle t ∈ [0, T ] beliebig, aber fest. Fur jeden Vektor h ∈ (0, T ]N existiert eine Zahlj = j(h) ∈ N0 mit der Eigenschaft

tj(h) ≤ t ≤ tj(h)+1.

Fur |h| → 0 folgt tj(h) → t, tj(h)+1 → t. Nach dem Mittelwertsatz gilt fur |h| → 0

1hj

(u(tj+1)− u(tj)) → u′(t).

Andererseits gilt wegen der Stetigkeit der Inkrementfunktion φ fur |h| → 0

φ(hj , tj , u(tj)) → φ(0, t, u(t)).

Also gilt mit der Definition von εh

εh(tj)|h|→0→ u′(t)− φ(0, t, u(t) = f(t, u(t))− φ(0, t, u(t)).

Wegen εh(tj) → 0 fur |h| → 0 (Konsistenz) folgt

f(t, u(t)) = φ(0, t, u(t)).

”⇐“: Fur den Konsistenzfehler gilt

|εh(tj)|∞ ≤∣∣∣∣u(tj+1)− u(tj)

hj− u′(tj)

∣∣∣∣∞

+ |u′(tj)− φ(hj , tj , u(tj))∣∣∞

=: |I1(j)|∞ + |I2(j)|∞ (j = 0, . . . , N)

Fur eine beliebige Zahl δ > 0 sei Ih so gewahlt, dass |h| ≤ δ.Wir wenden den Mittelwertsatz komponentenweise auf I1 an. Dann existiert fur einenIndex l = 1, . . . ,m eine Zahl ξlj ∈ (tj , tj+1), so dass gilt

I1l = u′l(ξlj )− u′l(tj)

Wegen der gleichmassigen Stetigkeit von u′l auf dem abgeschlossenen Intervall [0, T ] exis-tiert eine von j unabhangige Funktion r1 ∈ C0(R≥0,R≥0) mit r1(0) = 0 und

|I1l| ≤ r1(δ) (l = 1, . . . ,m).

Um I2 abzuschatzen, bemerken wir zunachst, dass die Inkrementfunktion φ gleichmassigstetig auf dem Kompaktum (k, t, u(t))| |k| ≤ δ, t ∈ I ist.Somit existiert r2 ∈ C0(R≤0R≤0), r2(0) = 0 mit

|φ(0, t, u(t))− φ(k, t, u(t))| ≤ r2(δ) ∀|k| ≤ δ, t ∈ I .

⇒ |I2(j)|(1.1)= |f(tj , u(tj))− φ(hj , tj , u(tj))|

V orr.= |φ(0, tj , u(tj)− φ(hj , tj , u(tj))| ≤ r2(δ)

Also insgesamt fur alle j ∈ 0, . . . , N :

|I1(j)|∞ + |I2(j)|∞ ≤ r1(δ) + r2(δ)

Da δ beliebig war, muss |εh(tj)|∞ → 0 fur |h| → 0 gelten. Also ist das Verfahren konsistentzu (1.1).

13


ii. Wir diskutieren zunachst die Funktion

γ = γ(t, k) :=1k(u(t+ k)− u(t)− kφ(k, t, u(t))) (t, t+ k ∈ I)

Taylorentwicklung liefert

γ(t, k) =1k

(p∑

i=0

ki

i!u(i)(t)− u(t) +O(|h|p+1)

−p−1∑i=0

∂i

∂kiφ(0, t, u(t))

ki+1

i!+O(|h|p+1)

)

=p∑

i=1

ki−1

i!u(i)(t)−

p−1∑i=0

∂i

∂kiφ(0, t, u(t))

ki

i!+O(|h|p)

=p−1∑i=0

(1

(i+ 1)!u(i+1)(t)− ∂i

∂kiφ(0, t, u(t))

1i!

)ki +O(|h|p)

=∑(

di

dti(f(t, u(t))− (i+ 1)

∂i

∂kiφ(0, t, u(t))

)ki

(i+ 1)!+O(|h|p)

= O(|h|p)

(?)

Also gilt nach Definition von γ

|εk(tj)|∞ = |γ(tj , hj)|∞ = O(|h|p) (j = 0, . . . , N)

”⇒“: Das Verfahren sei konsistent zur Ordnung p. Definiere

zi(t) =di

dti[f(t, u(t))]− (i+ 1)

∂iφ(0, t, u(t))∂ki

(t ∈ [0, T ])

Annahme: Es existiert t ∈ I und ein Index i ∈ 0, . . . , p− 1 mit

a) zi(t) 6= 0

b) zν(t) = 0 ∀t ∈ I , ν = 0, . . . , i− 1.

Wir wahlen den Index j = j(h) wie im Beweis von (i)”⇒“. Dann existiert eine Konstanteα > 0, so dass fur |h| hinreichend klein gilt:

α(a)

≤ |zi(tj)|∞

⇒ 1(i+ 1)!

αhij ≤

∣∣∣∣∣zi(tj) hij

(i+ 1)!

∣∣∣∣∣∞

(b)=

∣∣∣∣∣i∑

ν=0

zν(tj)hν

j

(i+ 1)!

∣∣∣∣∣∞

(?)fur γ(tj ,hj)=

∣∣∣∣∣ 1hj

(u(tj + hj)− u(tj)− φ(hj , tj , u(tj))−p−1∑

ν=i+1

zν(tj)hν

i

(i+ 1)!+O(|h|p)

∣∣∣∣∣≤ |εj(tj)|∞ + c

p−1∑ν=i+1

hνj +O(|h|p) = O

(|h|minp,i+1

)⇒ 1

(i+ 1)!αhi

j ≤ O(|h|minp,i+1

)Wegen i ≤ p− 1 ergibt sich ein Widerspruch fur die Gitterweite |h| → 0. Damit muss 1.8gelten.

14


Beispiel 1.20 (Konsistenzordnung fur das Eulerverfahren). Fur das explizite Eulerverfahrengilt

φ(hj , tj , u(tj)) = f(tj , u(tj)).

Nach Satz (1.18)(i) ist das Eulerverfahren konsistent. Falls f ∈ C1(I × Rm,Rm) ist, ist dasEulerverfahren von Konsistenzordnung 1. Das Eulerverfahren ist nicht konsistent zur Ordnung2, denn:

d

dtf(t, u(t)) =

∂

∂tf(t, u(t)) +

∂

∂uf(t, u(t))u′(t) =

∂

∂tf(t, u(t)) +

∂

∂uf(t, u(t))f(u(t)) (?)

und(1 + 1)

∂

∂kφ(0, t, u(t)) = 0

Wir konnen nicht erwarten, dass die rechte Seite von (?) verschwindet. Fur z.B.

f(t, w) = Aw,A ∈ Rm×m

gilt etwa:d

dtf(t, u(t)) = Au′(t) = A2u(t) 6= 0 fur u0 6= 0

Also ist nach Satz (1.18)(ii) das Eulerverfahren hochstens von Ordnung 1 fur alle f ∈ C1(I ×R

m,Rm).

Um eine Fehlerabschatzung fur ESV zu beweisen, benotigen wir einige Notationen.

Definition 1.21 (Diskreter Funktionenraum). Fur h ∈ (0, T ]N sei Ih ein Gitter. Die Menge

Xh := uh : Ih\tN = T → Rm|∃c > 0 : |uh(tj)|∞ ≤ c ∀j ∈ 0, . . . , N − 1

heisst Raum der beschrankten Gitterfunktionen.Mit der Norm ‖.‖∞ gegeben durch

‖uj‖∞ := maxj=0,...,N−1

|uh(tj)|∞ (uh ∈ Xh)

ist Xh = (Xh, ‖.‖∞) ein Banachraum (da isomorph zu RmN ).

Definition 1.22 (Diskreter Operator). Sei fur ein Gitter Ih ein ESV mit Inkrementfunktionφ ∈ C0(I2 ×Rm,Rm) gegeben. Der Operator Th : Xh → Xh mit

(Th[uh])(tj) :=

uh(t0)− u0 fur j = 01hj

(uh(tj+1)− uh(tj)− hjφ(hj , tj , uh(tj))) fur j = 1, . . . , N − 1

heisst der dem ESV zugeordnete diskrete Operator.

uh ist genau dann die Gitterfunktion aus einem ESV, wenn Th[uh] = 0 fur den zugeordnetendiskreten Operator gilt.Fur ein konsistentes ESV (der Ordnung p ∈ N) gilt dann mit dem zugeordneten Operator Th

die Beziehunglim|h|→0

∥∥∥Th

(u∣∣Ih

)∥∥∥∞

= 0[∥∥∥Th

(u∣∣Ih

)∥∥∥∞

= O(|h|p)]

15


Definition 1.23 (Stabilitat). Fur alle h ∈ (0, T ]N sei ein Gitter Ih gegeben. Weiter sei einEinschrittverfahren mit Inkrementfunktion φ ∈ C0(I2 ×Rm,Rm). Das Einschrittverfahren mitzugeordnetem Operator Th : Xh → Xh heisst stabil genau dann, wenn gilt:

∃c, h > 0 : ‖uh1 − uh2‖∞ ≤ c‖Th(uh1)− Th(uh2)‖∞, (1.9)

fur alle uh1 , uh2 ∈ Xh und fur alle h ∈ (0, T ]N mit |h| ≤ h.

Eigenschaften von Stabilitat:

i. Wenn das Verfahren stabil ist, ist die Losung uh von Th(uh) = 0 eindeutig, denn:

‖uh − uh‖ ≤ c‖Th(uh)− Th(uh)‖ = 0

ii. Betrachte speziell uh1 = uh (Th(uh) = 0) und uh2 = 0

‖uh1 − 0‖∞ ≤ ‖Th(uh1)− Th(0)‖∞ = c‖Th(0)‖∞ ≤ cc0

‖uh‖∞ ≤ cc0 ∀|h| ≤ h

iii. Sei Th(w) = Aw − b fur A ∈ Rm×m, b ∈ Rm, so gilt

Aw − b = 0

‖w1 − w2‖∞ ≤ c‖Aw1 − b−Aw2 + b‖∞ = c‖A(w1 − w2)‖∞⇒ det(A) 6= 0

Satz 1.24. Sei ein ESV mit Inkrementfunktion φ ∈ C0(I2 ×Rm,Rm) gegeben. Weiter sei dasESV stabil. Dann gilt

i. Ist das ESV konsistent, dann giltlim|h|→0

eh = 0.

ii. Ist das ESV konsistent zur Ordnung p ∈ N, dann gilt

eh = O(|h|p)

fur |h| hinreichend klein.

Motto:Konsitenz und Stabilitat = Konvergenz

Beweis: Sei Th der dem ESV zugeordnete diskrete Operator. Dann ist

eh = ‖uh − u|Ih‖∞

stab.≤ c‖Th(u|Ih

)− Th(uh)︸︷︷︸=0

‖∞ = c‖Th(uh|Ih)‖∞

Konsist.=o(1) fur Fall (i)O(|h|p) fur Fall (ii)

Im Beweis wird die Tatsache, dass der Operator Th einem ESV fur eine gewohnliche Anfangs-wertaufgabe zugeordnet ist, nicht benutzt. Die Beweistechnik kann fur andere Verfahren zurLosung anderer Probleme eingesetzt werden. Das zentrale Resultat dieses Kapitels ist

16


Satz 1.25 (Konvergenz von ESV I). Sei ein Einschrittverfahren mit Inkrementfunktion φ ∈C0(I2 × Rm,Rm) und h > 0 gegeben. Weiterhin existiere eine Konstante M > 0, so dass furalle t ∈ I , w, w ∈ Rm und k ∈ [0, h]

|φ(k, t, w)− φ(k, t, w)|∞ ≤M |w − w|∞ (1.10)

gilt. Dann gelten folgende Aussagen:

i. Das ESV ist stabil.

ii. Falls das ESV konsistent ist, konvergiert das ESV.

iii. Falls das ESV konsistent zur Ordnung p ∈ N ist, existiert eine Konstante c > 0, so dassfur alle Gitter Ih mit |h| ≤ h die Abschatzung

eh ≤ ccs|h|p

gilt. Dabei ist cs explizit gegeben durch

cs = eMT · (T + 1).

Beweis:

i. Wir mussen die Stabilitatsungleichung (1.9) fur den dem ESV zugeordneten diskretenOperator Tn nachweisen. Dazu seien w1, w2 ∈ Xn. Setze

z := w1 − w2, v1 := Tn[w1], v2 := Tn[w2].

Fur j = 0, . . . , N − 1 gilt:

1hj

(z(tj+1)− z(tj)) =1hj

(w1(tj+1)− w1(tj)− hjφ(hj , tj , w1(tj))

− 1hj

(w2(tj+1)− w2(tj)− hjφ(hj , tj , w2(tj))

+ φ(hj , tj , w1(tj))− φ(hj , tj , w2(tj))=v1(tj)− v2(tj) + φ(hj , tj , w1(tj))− φ(hj , tj , w2(tj))

|z(tj+1)|∞ = |z(tj) + hjv1(tj)− hjv2(tj) + hjφ(hj , tj , w1(tj))− hjφ(hj , tj , w1(tj))|∞≤ |z(tj)|∞ + hj |v1(tj)− v2(tj)|+ hjM |w1(tj)− w2(tj)|∞≤ (1 + hjM)|z(tj)|∞ + hj max

j=0,...,N|v1(tj)− v2(tj)|∞

= (1 + hjM)|z(tj)|∞ + hj%

(?)

Sei a ∈ C1(I ,R≥0) a′(t) ≤ ca(t) + b. Mit Greenwell ist

a(t) ≤ ect(a(0) + bt)

Diskret:aj+1 ≤ (1 + hjM)aj + btjhj

aj+1 − aj

hj≤Maj + btj

17


mitaj+1 = a(tj+1)

Ein Induktionsbeweis (s.u.) fuhrt auf folgende Abschatzung

|z(tj)|∞ ≤j−1∏l=0

(1 + hlM)︸︷︷︸≤ehlM

|z(t0)|∞ + %

j−1∑l=0

j−1∏m=l+1

(1 + hmM)︸︷︷︸=ehmM

hl (??)

⇒ |z(tj)|∞ ≤ exp

(j−1∑l=0

hlM

)|z(t0)|∞ + %

j−1∑l=0

exp

(M

j−1∑m=l+1

hm

)hl

≤ exp (M(tj − t0)) |z(t0)|∞ + %

j−1∑l=0

exp ((tj − tl+1)M)hl

≤ exp(Mtj)|z(t0)|∞ + %eMtj tj

= exp(Mtj)(|z(t0)|∞ + %tj)

Wegen|z(t0)|∞ = |w1(t0)− u0 + u0 − w2(t0)|∞ = |v1(t0)− v2(t0)|∞ ≤ %

folgt‖z‖∞ ≤ eMT (1 + T )% = cs‖Th[w1]− Th[w2]‖∞ (? ? ?)

ii. ist eine Folgerung aus Satz 1.24.

iii. Aus (? ? ?) und dem Beweis von Satz 1.24 erhalten wir

eh ≤ eMT (1 + T )‖Th[u|Ih]− Th[uh]‖∞ = csεh ≤ csc|h|p fur |h| < 1

Wir mussen noch die Formel (??) durch Indunktion beweisen:j = 0 ist klar wegen

∏∅ = 1 und

∑∅ = 0, d.h.

|z(t0)|∞ ≤ |z(t0)|∞

j → j + 1:

|z(tj+1)|∞(?)

≤ (1 + hjM)|z(tj)|∞ + hj%

≤j∏

l=0

(1 + hjM)|z0(t0)|∞ + %

j−1∑l=0

j∏m=l+1

(1 + hmM)hl + hj%)

=j∏

l=0

(1 + hjM)|z0(t0)|∞ +j∑

l=0

(j∏

m=l+1

(1 + hmM))%hl

)

Bemerkung 1.26. Die Formel aus Satz 1.25.iii) fur cs zeigt, dass wir im Allgemeinen Kon-vergenz auf dem Intervall I = [0,∞) nicht erwarten konnen.(Die Formel ist scharf fur den Fall u′(t) = Au(t), A ∈ Rm×m\0.)

Satz 1.25 ist insofern unbefriedigend, als dass die Inkrementfunktion global Lipschitz-stetig seinmuss.

18


Satz 1.27 (Konvergenz von ESV II). Sei ein ESV mit Inkrementfunktion φ ∈ C2(I ×Rm,Rm)gegeben. Es existieren h > 0, M > 0 und ε > 0, so dass

|φ(k, t, w)− φ(k, t, w)|∞ ≤M |w − w|∞ (1.11)

fur k ∈ (0, h], t ∈ I und w, w ∈ v ∈ Rm|∃t ∈ I : |v − u(t)| ≤ ε gilt. Dann gelten alle Aussagenaus Satz 1.25:

i. Stabilitat

ii. Falls Konsistenz, dann Konvergenz

iii. Falls Konsistenz zur Ordnung p, dann Konvergenz zur Ordnung p

Beweis: Wir betrachten das ESV mit einer abgewandelten Inkrementfunktion

φ ∈ C0(I2 ×Rm,Rm)

definiert durchφ(k, t, w) = φ(k, t, S[w]).

Dabei sei die i-te Komponente von S[w] durch

(S[w])i = ψ(wi ;ui(t)− ε, ui(t) + ε)

gegeben durchψ(x; a, b) := minb,maxa, x

Offenbar ist ψ(., a, b) Lipschitz-stetig mit der Lipschitzkonstanten 1. Also gilt fur w, w ∈ Rnm

|S[w]− S[w]|∞ ≤ |w − w|∞

⇒ |φ(k, t, w)− φ(k, t, w)|∞ = |φ(k, t, S[w])− φ(k, t, s[w])|∞1.11≤ M · |S[w]− S[w]|∞ ≤M · 1 · |w − w|∞

(?)

fur k ∈ (0, h], t ∈ [0, T ].Somit erfullt also φ die Voraussetzung 1.10 aus Satz 1.25 und mit dem zugeordneten diskretenOperator Tn folgt uh1 , uh2 ∈ Xh mit h ∈ (0, h]

‖uh1 − uh2‖∞ ≤ Cs · ‖T [uh1 ]− T [uh2 ]‖∞ (??)

Weiterhin ist das ESV mit Inkrementfunktion φ konsistent (Konsistent zur Ordnung p), dennwir haben

Th[u|Ih] = Th[u|Ih

] (? ? ?)

Aus den Bedingungen (?),(??),(? ? ?) folgt mit Satz 1.25

eh = ‖uh − u|Jh‖∞

|h|→o−→ 0

bzw.‖uh − u|Jh

‖∞ = O(|h|p).Aus der Konvergenz der Funktionenfolge uh fur |h| → 0 ergibt sich aber, dass ein h > 0existiert, so dass

‖uh − u|Jh‖∞ ≤ ε ∀|h| ≤ h

gilt. Dann sind aber φ und φ identisch, d.h. Th = Th und

Th[uh] = Th[uh] = 0,

aber auch Th[uh] = 0.Da uh und uh offenbar eindeutig als Losungen eines ESV sind, folgt uh ≡ uh.

19


1.5 Explizite Runge-Kutta-Verfahren

Sei p ∈ N0 gegeben. Wir nehmen an, dass (1.1),(1.2) genau eine Losung u ∈ Cp+1(I,Rm) besitzt(etwa der Fall fur f ∈ Cp(I ×Rm,Rm)). Wie kann in einer systematischen Weise ein ESV mitKonsistenzordnung p konstruiert werden?

Beispiel 1.28 (Nochmal das Heunverfahren).In Beispiel 1.16 (iii) war das Heunverfahren

uj+1 = uj + hj ·12· [f(tj , uj) + f(tj + hj , uj + hjf(tj , uj))]

vorgestellt worden. Das Heunverfahren ist ein ESV mit Inkrementfunktion

φ = φ(k, t, w) =12(f(t, w) + f(t+ k,w + k · f(t, w))).

Das Heunverfahren erreicht die hohere Genauigkeit (Konsistenzordnung 2) durch mehrfacheinterative Anwendung der Funktion f .

Definition 1.29 (Explizite Runge-Kutta-Verfahren). Fur r ∈ N seien α2, . . . , αr ∈ R, γ1, . . . , γr ∈R, und βij ∈ R (i = 2, . . . , r, j = 1, . . . , r − 1, i > j) gegeben. Ein ESV mit Inkrementfunktion

φ = φ(k, t, w) :=r∑

i=1

γiki(k, t, w)

und

k1 = k1(k, t, w) = f(t, w)k2 = k2(k, t, w) = f(t+ α2k,w + k · β21 · k1(k, t, w))

...

kr = kr(k, t, w) = f(t+ αrk,w + k

r−1∑s=1

βrs · ks(k, t, w))

heisst allgemeines (explizites) Runger-Kutta-Verfahren (RK). Die Zahl r ∈ N heisst Stufe desRK-Verfahrens.Die Koeffizienten des allgemeinen RK-Verfahrens werden oft in der Form einer Tabelle (Butcher-Tableau) zusammengefasst:

α2 β21...

. . .αr βr1 . . . βr,r−1

γ1 . . . γr−1 γr

Das Eulerverfahren bzw. das Heunverfahren sind RK-Verfahren mit den Tabellen

φ(k, t, w) = f(t, w)

α1 := 01

bzw.φ(k, t, w) =

12(f(t, w) + f(t+ k,w + k · f(t, w)))

α1 := 01 1

1/2 1/2

20

1.5 Explizite Runge-Kutta-Verfahren

Bemerkung 1.30. Setzen wir noch α1 := 0, konnen die Koeffizientenfunktionen ki iterativgemaß

ki(k, t, w) = f

(t+ αik,w + k ·

i−1∑s=1

βisks(k, t, w)

)bestimmt werden. In dieser Form sollte die Implementierung realisiert sein.

Konnen die (freien) Koeffizienten des RK-Verfahren so bestimmt werden, dass ESV der Ordnungp resultieren?

Beispiel 1.31 (Das klassische RK-Verfahren). Ein RK-Verfahren der Ordnung p = 4 ergibt sichaus dem Butcher-Tableau:

12

12

12 0 1

21 0 0 1

16

13

13

16

Dieses wird als klassisches Runge-Kutta-Verfahren bezeichnet.Explizit ausgeschrieben ergibt sich das ESV

k1 = f(tj , uj)

k2 = f(tj +hj

2, uj +

12hjk1)

k3 = f(tj +hj

2, uj +

12hjk2)

k4 = f(tj + hj , uj + hjk3)

uj+1 = uj +hj

6(k1 + 2k2 + 2k3 + k4)

Wie viele RK-Verfahren 2-ter Stufe der Konsistenzordung 2 gibt es?Sei f ∈ C2(I ×Rm,Rm). Die Inkrementfunktion fur das RK-Verfahren hat die allgemeine Form

φ(k, t, w) = γ1 · f(t, w) + γ2 · f(t+ α2k,w + k · β21f(t, w))

Aus dem Konsistenzsatz 1.19 erhalten wir die beiden Bedingungen

f(t, w) != φ(0, t, w) = γ1f(t, w) + γ2f(t, w)

und

ft(t, w) + fw(t, w) · f(t, w)︸︷︷︸=u′

=d

dt[f(t, u(t))]

∣∣u(t)=w

==i+1︷︸︸︷

2 ·∂φ(k = 0, t, w)∂k

= 2 · γ2[ft(t, w) · α2 + fw(t, w) · β21 · f(t, w)]

Ein Koeffizientenvergleich fuhrt uns auf die algebraischen Gleichungen

1 = γ1 + γ2, 1 = 2γ2α2, 1 = 2γ2β21

Wir haben 3 Gleichungen mit 4 Unbekannten. Also nehmen wir zum Beispiel γ2 als Parameter,dann gilt

γ1 = 1− γ2, α2 =1

2γ2, β21 =

12γ2

oder als Tabelle geschrieben:01

2γ2

12γ2

1− γ2 γ2

Speziell erhalten wir fur

21


• γ2 = 12 das Heunverfahren

• γ2 = 1 das modifiziertes Eulerverfahren.

Bemerkung: Im Allgemeinen erhalt man ein nichtlineares Gleichungssystem.

Lemma 1.32. Sei r ∈ N. Die Konsistenzordnung eines r-stufigen RK-Verfahrens ist durch rnach oben beschrankt.

Beweis: Ubungsblatt 4.

Leider ist die Grenze r fur die Konsistenzordnung nicht scharf. Sei p = p(r) die maximaleKonsistenzordnung, die mit einem Verfahren r-ter Stufe erreicht werden kann. Folgendes istbekannt

r 1 2 3 4 5 6 7 8 r ≥ 9p(r) 1 2 3 4 4 5 6 6 ≤ r − 2

Fur Beweise und explizite Formeln entsprechender RK-Verfahren sei auf das Buch vonE. Hairer, S.P.Norse, C. Warner, Solving ordinary differential equations I, Springer(1993) verwiesen.Zum Abschluss des Kapitels zu expiziten RK-Verfahren wird noch ein konstruktives Resultatangegeben.

Proposition 1.33. Ein RK-Verfahren r-ter Stufe, r ∈ N, ist konsistent, falls

r∑i=1

γi = 1

gilt. Im Falle u′ 6≡ 0 ist diese Bedingung auch notwendig.

Beweis: Aus der Definition der ki (Def. 1.29) fur i = 1, . . . , r folgt

ki(0, t, w) = f(t, w)

⇒ φ(0, t, w) =r∑

i=1

γif(t, w) n.V.= f(t, w)

Satz 1.16 liefert die Behauptung.Da Satz 1.16 auch eine notwendige Bedingung liefert, die fur f 6= 0 auf

∑ri=1 γi = 1 fuhrt, folgt

der Rest der Aussage.

Bemerkung 1.34. Die Erhohung der Ordnung eines RK-Verfahrens wird durch eine hohereAnzahl von Elemtaroperationen erkauft. Daher ist nicht klar, ob Verfahren hoherer Ordnungauch effizienter sind.

1.6 Implizite Runge-Kutta-Verfahren

Die Ordnung eines r-stufigen expliziten RK-Verfahren ist durch r beschrankt. Um diese Schrankeaufzuheben, werden jetzt implizite ESV vom RK-Typ betrachtet.

Definition 1.35 (Implizite RK-Verfahren). Fur r ∈ N seien

α = (α1, . . . , αr)T ∈ Rr, (γ1, . . . , γr)T ∈ Rr

22

1.6 Implizite Runge-Kutta-Verfahren

und B = (bij) ∈ Rr×r gegeben.Ein ESV mit Inkrementfunktion

φ = φ(k, t, w) :=r∑

i=1

γiki(k, t, w) (1.12)

heisst allgemeines implizites RK-Verfahren, falls das nichtlineare Gleichungssystem

k1 =f

(t+ α1k,w + k

r∑s=1

b1sks

)...

kr =f

(t+ αrk,w + k

r∑s=1

brsks

) (1.13)

erfullt ist. Die Zahl r heisst Stufe des impliziten RK-Verfahrens.

In Analogie zu expliziten RK-Verfahren fasst man die Koeffizienten folgendermaßen zusammen(Implizites Butcher Tableau):

α1 b11 . . . b1r...

......

αr br1 . . . brr

γ1 . . . γr

Bemerkung 1.36. Das nichtlineare Gleichungssystem (1.13) ist von der Dimension rm ×m.Falls eine Losung des Systemes (1.13) existiert, erhalten wir durch Einsetzen in (1.12) einexplizites Verfahren. Dies ist aber oft analytisch nicht durchfuhrbar, so dass im Allgemeinen injedem Zeitschritt ein nichtlineares System numerisch zu losen ist.Die Konvergenztheorie aus 1.4 kann aber angewendet werden.

Beispiel 1.37 (Implizite RK-Verfahren der Stufe 1 und 2). Es sei m = 1 und r = 1. Dies fuhrtgemaß (1.13) auf die folgende nichtlineare Gleichung fur k1:

k1 = f(t+ α1k,w + kb11k1) (?)

Annahme: Es existiert genau ein k1 = k1(k, t, w) ∈ C1(I2 × R), so dass (?) gilt. Also habenwir die Inkrementfunktion

φ = φ(k, t, w) = γ1k1(k, t, w) (??)

Um Konsistenz zu erhalten, mussφ(0, t, w) = f(t, w)

gelten. Dies ist richtig, falls wir γ1 = 1 wahlen (k1(0, t, w) = f(t, w)). Differentiation von (??)nach k liefert (in k = 0)

∂

∂kφ(0, t, w) = ft(t, w)α1 + fw(t, w)b11k1.

Weiter istd

dtf(t, u(t)) = ft(t, u(t)) + fw(t, u(t))f(t, u(t)).

Nach dem Konsistenzsatz muss fur die Konsistenzordnung 2 gerade

(1 + 1)∂

∂kφ(0, t, u(t)) =

d

dtf(t, u(t))

23


gelten. Dies ist fur α1 = b11 = 12 gegeben.

(Das einzige explizite RK-Verfahren erster Stufe ist das Euler-Verfahren. Dies besitzt nur dieKonsistenzordnung 1.)Wie sieht das obige Verfahren konkret aus?

uj+1 = uj + hk1

⇒ k1 =uj+1 − uj

h

uj +12hk1 = uj +

12(uj+1 − uj) =

12(uj + uj+1)

⇒ 1h

(uj+1 − uj) = f

(12(tj + tj+1),

12(uj + uj+1)

)Wir erhalten also wirklich ein implizites Verfahren. Ein implizites RK-Verfahren der Stufe 2und der Ordnung 4 ist durch folgendes Tableau gegeben:

3−√

36

14

14 −

√3

63+

√3

614 +

√3

614

12

12

Beispiel 1.37:

r = 1 ⇒ p(r) = 2

r = 2 ⇒ p(r) = 4

Bemerkung 1.38. Das Beispiel 1.37 zeigt, dass implizite RK-Verfahren auf hohere Konsisten-zordnungen als explizite RK-Verfahren fuhren konnen.Als Nachteil muss erwahnt werden, dass ein hochdimensionales Gleichungssystem (approximativ)gelost werden muss. Dies fuhrt unter Umstanden zu einer weiteren Verringerung der Schrittwei-te.

Eine etwas einfachere Situation ergibt sich fur sogenannte halbimplizite RK-Verfahren, fur diedie Eintrage von B = (bis) die Bedingung

bis = 0 (i < s)

erfullt ist. Das System (1.13) vereinfacht sich dann zu

k1 =f(t+ α1k,w + kk1b11)...

kr =f

(t+ αrk,w +

r∑s=1

brskks

).

Man kann also die Koeffizientenfunktionen k1, . . . , kr sukzessive durch Losen von r Gleichungs-systemen der Dimension m×m bestimmen.Ein spezielles halbimplizites Verfahren wird in den Ubungen vorgestellt.

24

1.7 Runge-Kutta-Verfahren und Quadraturformeln


In diesem Kapitel soll ein Weg vorgestellt werden, der systematisch zu RK-Verfahren einergegebenen Ordnung p ∈W fuhrt.Allgemeines RK-Verfahren r-ter Stufe:

uj+1 = uj + hj

r∑i=1

γiki

k1 =f

(tj + α1hj , uj + hj

r∑s=1

b1sks

)...

ki =f

(tj + αihj , uj + hj

r∑s=1

bisks

)...

kr =f

(tj + αrhj , uj + hj

r∑s=1

brsks

)Vorbemerkungen:

i. Eine notwendige Bedingung fur Konsistenzordnung p ∈ N ergibt sich aus folgender Uber-legung. Betrachte das RK-Verfahren fur

u′(t) = g(t) (t ∈ I), u(0) = u0 ∈ Rm.

Dabei sei g : I → Rm. Wir erhalten

1hj

(uj+1 − uj) =r∑

i=1

γig(tj + αihj)

und damit fur den Konsistenzfehler

εj(tj) =1hj

(u(tj+1)− u(tj)− hj

r∑i=1

γig(tj + αihj)

)

=1hj

(∫ tj+1

tj

g(t)dt− hj

r∑i=1

γig(tj + αihj)

)

Um eine gegebene Konsistenzordnung p ∈ N zu erreichen, muss also gelten:∣∣∣∣∣∫ tj+1

tj

g(t)dt−r∑

i=1

hjγig(tj + αihj)

∣∣∣∣∣ = O(hp+1

j

)(1.14)

Dies ist ein reines Quadraturproblem!Die Koeffizienten α1, . . . , αr und γ1, . . . , γr sind damit festgelegt.

ii. Wir setzen die exakte Losung u = u(t) von (1.1), (1.2) in das RK-Verfahren ein:

u(tj+1)− u(tj)hj

≈r∑

i=1

γiki(hj , tj , u(tj))

25


”⇔ ”∫ tj+1

tj

u′(t)dt ≈ hj

r∑i=1

γiki(hj , tj , u(tj))

Damit die letzte Beziehung mit hoher Genauigkeit gilt, sollte fur die gegebenen Werteγ1, . . . , γr mit hoher Genauigkeit

ki ≈ u′(tj + αihj) (i = 1, . . . , r) (1.15)

gelten. Wegen der Definition der ki und der DGL (1.1) folgt

ki = f

(tj + αihj , u(tj) +

r∑s=1

hjβisks

)≈ f(tj + αihj , u(tj + αihj))

”⇔ ”u(tj + αihj) ≈ u(tj) + hj

r∑s=1

βisks

”⇔ ”∫ tj+αihi

tj

u′(t)dt ≈ hj

r∑s=1

βisu′(tj + αshj) (siehe (1.15) fur die rechte Seite.)

Die Zahlen βis (i, s = 1, . . . , r) sollten so gewahlt werden, dass fur eine Funktion h : I →R

m mit hoher Genauigkeit

1hj

∫ tj+αihj

tj

h(t)dt ≈r∑

s=1

βish(tj + αshj) (1.16)

gilt (i = 1, . . . , s). Dies ist wiederum ein reines Quadraturproblem!

Den Erfolg dieses Konzeptes belegt

Satz 1.39 (Butcher, Kunzmann ’69). Es sei ein RK-Verfahren fur α, γ ∈ Rr und B = (bis) ∈R

r×r gegeben. Fur p, q ∈ N seien die Koeffizienten so gewahlt, dass fur g ∈ Cp+1(I ,Rm) undh ∈ Cq+1(I ,Rm) gelte

i. ∣∣∣∣∣ 1hj

∫ tj+1

tj

g(t)dt−r∑

i=1

hjγig(tj + αihj)

∣∣∣∣∣ = O(hpj ) (j = 0, . . . , N − 1)

ii. ∣∣∣∣∣ 1hj

∫ tj+αihj

tj

h(t)dt−r∑

s=1

hjβish(tj + αihj)

∣∣∣∣∣ = O(hqj) (j = 0, . . . , N − 1, i = 1, . . . , r)

Dann ist das RK-Verfahren konsistent von der Ordnung minp, q + 1.

Beweis: Siehe Deuflhard & Bornemann, Numerische Mathematik II.

Um Quadraturprobleme aus Satz 1.39 moglichst gut zu losen, wird der Begriff der Exaktheiteingefuhrt.

Definition 1.40 (Exaktheit). Es sei h ∈ C0([0, 1],R) und τ ∈ (0, 1]. Fur γ1, . . . , γr ∈ R undα1, . . . , αr ∈ [0, 1] sei eine Quadraturformel

Q[h] :=r∑

i=1

γih(αi)

fur das Integral∫ τ0 h(t)dt gegeben. Q heisst exakt vom Grad l ∈ N

:⇔ Q[p]−∫ τ

0p(t)dt = 0 ∀p ∈ Pl

26


Satz 1.41. Fur l ∈ N sei h ∈ Cl+1([0, 1],R) und τ ∈ (0, 1]. Fur γ, α ∈ Rr mitα1, . . . , αr ∈ [0, 1] sei eine Quadraturformel

Q[h] :=r∑

i=1

γih(αi)

fur das Integral∫ τ0 h(t)dt gegeben, die exakt vom Grad l ∈ N sei. Dann gilt∫ τ

0h(t)dt = Q[h] +

∫ 1

0πl+1(t)h(l+1)(t)dt,

wobei πl+1 der Peano-Kern

πl+1(t) :=1

(l + 1)!

(((τ − t)+)l+1 − (l + 1)

r∑i=1

γi((αi − t)+)l

)(t ∈ [0, 1])

ist. Hier sei

α+(t) =a(t) fur a > 00 fur a ≤ 0

.

Beweis: Aus dem Satz von Taylor folgt fur t ∈ [0, 1].

h(t) =l∑

k=0

h(k)(t)k!

tk +∫ t

0

1l!h(l+1)(s)(t− s)lds

=l∑

k=0

h(k)(t)k!

tk +∫ 1

0

1l!h(l+1)(s)((t− s)+)lds

Mit I = I[h] =∫ τ0 h(t)dt folgt wegen der Exaktheit der Quadraturformel Q

I[h]−Q[h] =∫ τ

0h(t)dt−Q[h]

=∫ 1

0

1l!h(l+1)(s)(I −Q)[((· − s)+)l]ds

(?)

Wir berechnen

I[((· − s)+)l] =∫ τ

0((t− s)+)ldt

=

0 fur s ≥ τ∫ τs (t− s)ldt fur s < τ

=1

l + 1((τ − s)+)l+1

und

Q[((l − s)+)l] =r∑

i=1

γr((αi − s)+)l,

so dass mit (?) folgt

(I −Q)(h) =∫ 1

0

1l!h(l+1)(s)

(1

l + 1((τ − s)+)l+1 −

r∑i=1

γi((αi − s)+)l

)ds

=∫ 1

0

1(l + 1)!

h(l+1)(s)

(((τ − s)+)l+1 − (l + 1)

r∑i=1

γi((αi − s)+)l

)ds

27


Korollar 1.42. Sei g ∈ Cl+1([0, τ ],Rm) und Ih ein Gitter. Mit den Bezeichnungen aus Satz1.41 gilt dann fur n ∈ 1, . . . , r und j ∈ 0, . . . , N − 1∣∣∣∣∣ 1

hj

∫ tj+αnhj

tj

g(t)dt−r∑

i=1

γig(tj + αjhj)

∣∣∣∣∣ = O(|h|l+1

)Beweis: O.B.d.A. sei m = 1. Wir wahlen h = h(t) := g(tj + thj) in Satz 1.41 und erhalten furτ = αn ∫ αn

0h(t)dt =

∫ αn

0g(tj + thj)dt

s=tj+thj=∫ tj+αnhj

tj

g(s)1hjds

(1.41)⇒

∣∣∣∣∣ 1hj

∫ tj+αnhj

tj

g(t)dt−r∑

i=1

γig(tj + αihj)

∣∣∣∣∣=∣∣∣∣ ∫ 1

0πl+1(t) h(l+1)(t)︸︷︷︸

g(l+1)(tj+thj)hl+1j

dt

∣∣∣∣≤ max

t∈[0,1]

∣∣∣πl+1(t)g(l+1)(tj + thj)∣∣∣hl+1

j

= O(|h|l+1

)

Definition 1.43 (Legendre Polynome). Fur m ∈ N ist das m-te Legendre Polynom pm durch

pm(t) =m!

(2m)!dm

dtm(t2 − 1)m (t ∈ R)

definiert.

Offenbar gilt grad(pm) = m und weiter sind

p0 = 1, p1 = t, p2 =13(3t2 − 1), · · ·

Lemma 1.44 (Nullstellen & Orthogonalitat).

i. Die Legendrepolynome pm,m ∈ N, besitzen m paarweise verschiedene Nullstellen %1, . . . , %m

mit−1 < %1 < %2 < · · · < %m < 1.

ii. Fur n,m ∈ N und n 6= m gilt ∫ 1

−1pn(t)pm(t)dt = 0

Beweis: Ubungsblatt.

Satz 1.45. Sei h ∈ C([−1, 1],R). Betrachte die Gauß-Quadraturformel

Q = Q[h] =r∑

i=1

ωih(%i)

28

1.8 Mehrschrittverfahren

mit

ωi =∫ 1

−1

r∑j=1,j 6=i

t− %j

%i − %jdt (i = 1, . . . , r)

Dann gilt

Q[p]−∫ 1

−1p(t)dt = 0 ∀p ∈ P2r−1.

Beweis: Die Polynome p0, . . . , pr bilden eine Basis von Pr (wegen der Orthogonaleigenschaftaus Lemma (1.44)(ii)),

⇒∫ 1

−1pr(t)q(t)dt = 0 ∀q ∈ Pr−1. (?)

Fur p ∈ Pr−1 gilt bei Division mit Rest fur Polynome π, s ∈ Pr−1

p = pr · π + s

Wegen (?) folgt fur das beliebige Polynom p∫ 1

−1p(t)dt =

∫ 1

−1pr(t)π(t)︸︷︷︸

=0

+s(t)dt =∫ 1

−1s(t)dt = Q[s],

da Q insbesondere interpolatorische Quadraturformel mit r Stutzstellen ist, so dass s ∈ Pr−1

exakt integriert wird. Aus der Definition der Quadraturformel folgt weiter:∫ 1

−1p(t)dt = Q[s] =

r∑i=1

ωi · σ(%i) =r∑

i=1

ωi(s(%i) + pr(%i)π(si)) = Q[p],

da pr(%i)π(si) = 0 (%i Nullstellen der Legendrepolynome).


Definition 1.46 (Mehrschrittverfahren). Sei ein aquidistantes Gitter Ih gegeben und sei

ψ ∈ C0(Ik+2 ×Rm(k+1),Rm)

fur k ∈ N. Weiter seien a0, . . . , ak ∈ R und u0 = u(t0), u1, . . . , uk−1 ∈ Rm gegeben. DasVerfahren

1h

(a0uj + a1uj+1 + . . .+ akuj+k) = ψ(h, tj , . . . , tj+k, uj , . . . , uj+k) (1.17)

heisst K-Mehrschrittverfahren (K-MSV) mit Verfahrensfunktion ψ.

Im Gegensatz zu RK-Verfahren soll bei MSV die Ordnung durch Benutzung alter Approximati-onswerte und nicht durch vielfache Auswertung von f erhoht werden.

Bemerkung 1.47.

i. Falls ψ in (1.17) nicht von uj+k abhangt und ak 6= 0 gilt, heisst das K-MSV explizit.

ii. Fur k = 1 ergibt sich im expliziten Fall ein explizites ESV.

29


iii. Um (1.17) ausfuhren zu konnen, mussen u0, . . . , uk−1 bekannt sein. Diese sollten mit einemESV derselben Konsistenzordnung (siehe unten) berechnet werden.

iv. Falls die Verfahrensfunktion ψ von der Form

ψ(h, tj , . . . , tj+k, uj , . . . , uj+k) =k∑

i=0

bif(ti+j , uj+i)

fur b0, . . . , bk ist, heisst das zugehorige MSV linear.

Definition 1.48 (Fehlerdefinition fur lineare MSV). Sei uh : Ih → Rm durch ein lineares MSV

erzeugt.

i. Der globale Diskretisierungsfehler eh ist durch

eh := maxj=0,...,N

|eh(tj)|, eh := u∣∣Ih− uh

gegeben.

ii. Die Funktion εh : Ih\t0, . . . , tk−1 → Rm mit

εh(tj+k) =1h

k∑i=0

aiu(tj+i)−k∑

i=0

bif(tj+i, u(tj+i)) (j = 0, . . . , N − k)

heisst lokale Fehlerfunktion.Der lokale Diskretisierungsfehler εh ist durch

εh := maxj=0,...,N−k

|εh(tj+k)|

gegeben.

Wie mussen wir die Koeffizienten a0, . . . , ak und b0, . . . , bk eines linearen K-MSVs wahlen, damitdas Verfahren die Konsistenzordnung p ∈ N besitzt, d.h. wann gilt

εh = O(hp)?

Fur den Konsistenzfehler berechnen wir (j = 0, . . . , N − k)

εh(tj+k) =1h

k∑i=0

aiu(ti+j)−k∑

i=0

bif(tj+i, u(tj+i))

=1h

k∑i=0

aiu(tj + ih)−k∑

i=0

bi · u′(tj + ih)

Mit der Definition pl = pl(i) := (l!)−1il, l ∈ N0 ergibt sich aus dem Satz von Taylor

|u(t+ ih)| =

∣∣∣∣∣p∑

l=0

hlu(l)(t)pl(i)

∣∣∣∣∣+O(hp+1)

|u′(t+ ih)| =

∣∣∣∣∣p∑

l=0

hlu(l+1)(t)pl(i)

∣∣∣∣∣+O(hp+1)

=

∣∣∣∣∣p∑

l=1

hl−1u(l)(t)pl−1(i)

∣∣∣∣∣+O(hp)

=

∣∣∣∣∣p∑

l=1

hl−1u(l)(t)p′l(i)

∣∣∣∣∣+O(hp)

30


|εh(tj+k)| =

∣∣∣∣∣1hk∑

i=0

ai

(p∑

l=0

hlu(l)(tj)pl(i)

)−

k∑i=0

bi

(p∑

l=1

hl−1u(l)(t)p′l(i)

)∣∣∣∣∣+O(hp)

=

∣∣∣∣∣p∑

l=0

hl−1u(l)(tj+k)

(k∑

i=0

aipl(i)−k∑

i=0

bip′l(i)

)∣∣∣∣∣+O(hp)

Proposition 1.49. Falls die Koeffizienten a0, . . . , ak ∈ R, b0, . . . , bk ∈ R eines K-MSV dieBedingung

k∑i=0

aipl(i) =k∑

i=0

bi p′l(i)︸︷︷︸=pl−1(i)

(1.18)

fur l = 0, . . . , p erfullen, besitzt das K-MSV die Konsistenzordnung p.

(1.18) ist ein lineares System mit p+ 1 Gleichungen und 2k + 2 Unbekannten.Da a0 = . . . = ak = b0 = . . . bk = 0 eine unsinnige Losung ist, fugen wir noch die Normierungs-bedingung

k∑i=0

bi = 1 (1.19)

hinzu.

Beispiel 1.50 (Trapezverfahren). Wir betrachten den Fall k = 1 mit den 4 Unbekannten

a0, a1, b0, b1.

Fur p = 2 ergeben sich aus (1.18) und (1.19) die 4 Gleichungen l = 0:

a0 · 1 + a1 · 1 = 0 (p0 = 1)

l = 1:a0 · 0 + a1 · 1 = b0 + b1 (p1 = i)

l = 2:

a0 · 0 + a1 ·12

= b0 · 0 + b1 · 1(p2 =

i2

2

)1 = b0 + b1

⇒ a1 = b0 + b1 = 1, ⇒ a0 = −1 ⇒ 12

= b1 ⇒ b0 = 1− b1 =12

Das zugehorige 1-MSV ist

1h

(−uj + uj+1) =12(f(tj , uj) + f(tj+1, uj+1))

⇒ uj+1 = uj +h

2(f(tj , uj) + f(tj+1, uj+1)) (Trapezverfahren)

Zur Motivation betrachte

1h

((u(tj+1)− u(tj)) =1h

∫ tj+1

tj

u′(t)dt =1h

∫ tj+1

tj

f(t, u(t))dt

Trapezregel liefert

≈ 12(f(tj , u(tj)) + f(tj+1, u(tj+1)))

31


K-MSV kann man systematisch herleiten, indem man die zu (1.1),(1.2) aquivalente Integralglei-chung

u(t) = u(s) +∫ t

sf(r, u(r))dr

betrachtet, wobei t > s, t, s ∈ I sei.Wir wahlen speziell fur ein aquidistantes Gitter Ih und j ∈ 0, . . . , N − k t = tj+k, s = tj+k−1

u(tj+k) = u(tj+k−1) +∫ tj+k

tj+k−1

f(r, u(r))dr

= u(tj+k−1) +∫ tj+k

tj+k−1

u′(r)dr

Fur j = 0, . . . , N − k sei pj ∈ Pk−1 das eindeutig bestimmte Interpolationspolynom, das sichaus den k Bedingungen

pjk(tj+i) = u′(tj+i) = f(tj+i, u(tj+i)) (i = 0, . . . , k − 1)

ergibt.Einsetzen in die Integralgleichung liefert

u(tj+k) ≈ u(tj+k−1) +∫ tj+k

tj+k−1

pjk(r)dr︸︷︷︸=:I

(?)

Da pjk ein Polynom ist, kann I explizit bestimmt werden. Zur Berechnung werden nur Wertevon f(·, u(·)) in den Stutzstellen tj , . . . , tj+k−1 benotigt. Daher motiviert (?) ein K-MSV:

1h

(uj+k − uj+k−1) = ψ(h, tj , . . . , tj+k, uj , . . . , uj+k−1) :=1h

∫ tj+k

tj+k−1

pjk(r)dr

Dabei ist pjk ∈ Pn−1 das eindeutig bestimmte Interpolationspolynom, dass sich aus den Bedin-gungen

pjk(tj+i) = f(tj+i, uj+i) (i = 0, . . . , k − 1)

ergibt.Betrachtet man den Fall k = 4, erhalt man das explizite Adams-Bashforth-Verfahren der Stufe4:

1h

(uj+4 − uj+3) =124

(55f(tj , uj)− 59f(tj+1, uj+1) + 37f(tj+2, uj+2)− 9f(tj+3, uj+3))

Das Adams-Bashforth-Verfahren besitzt die Konsistenzordnung 4.

Bemerkung 1.51.

i. In jedem Zeitschritt eines Mehrschrittverfahrens benotigt man nur eine neue Auswertungvon f , namlich f(tj+k, uj+k). Man vergleiche dies mit der Situation fur RK-Verfahren.

ii. Mit der Vorgehensweise, die auf die expliziten MSV (Adams-Bashforth) fuhrt, kann manauch implizite Verfahren konstruieren. Dann ist zusatzlich (tj+k, uj+k) ein Stutzpunkt furdas Interpolationspolynom.Dann ergeben sich die Verfahren von Adams-Moulton. Fur k = 4 erhalt man

1h

(uj+4 − uj+3) =1

720(251f(tj+4, uj+4) + 646f(tj+3, uj+3)− 269f(tj+2, uj+2)

+ 106f(tj+1, uj+1)− 19f(tj , uj))

Das Verfahren besitzt die Konsistenzordnung 5.

32

1.9 Adaptive Schrittweitensteuerung

iii. Implizite und explizite Mehrschrittverfahren werden oft zu Pradiktor-Korrektor-Verfahrenkombiniert:Sei das explizite MSV

1h

k∑i=0

aiui+j = ψexp(h, tj , . . . , hj+k−1, uj , . . . , uj+k−1)

und das implizite MSV

1h

k∑i=0

αiui+j = ψimp(h, tj , . . . , tj+k, uj , . . . , uj+k)

gegeben. Fur das Pradiktor-Korrektor-Verfahren definiert man zunachst den Prediktor upk+j

vermoge1h

k−1∑i=0

aiui+j +1hup

k+j = ψexp(h, tj , . . . , tj+k−1, uj , . . . , uj+k−1)

und dann uk+j durch

1h

k∑i=0

αiui+j = ψimp(h, tj , . . . , tj+k, uj , . . . , uj+k−1, upj+k)

Insgesamt erhalten wir ein voll explizites Verfahren.


Sei ein ESV zur approximativen Losung von (1.1),(1.2) gegeben. Fur gegebenes Gitter Ih seiT (Ih) die Rechenzeit, die man benotigt, um uh zu bestimmen (auf einem Rechner). Weiterbezeichne d > 0 eine gegebene Fehlertoleranz.Aufgabe: Finde ein Gitter Iopt

h , so dass gilt:

i. eh ≤ d

ii. T (Iopth ) ≤ T (Ih) fur alle Gitter Ih mit eh = ‖u− uh‖∞ ≤ d.

Es ist unklar, ob ein solches ”ideales“ Gitter uberhaupt existent (und eindeutig ist).Die adaptive Gitterweitensteuerung soll nun eine moglichst gute Approximation von Iopt

h liefern.Grundlage der adaptiven Schrittweitensteuerung ist

Satz 1.52 (Fehlerentwicklung). Fur p ∈ N sei u ∈ Cp+2(I ,Rm) die Losung von (1.1),(1.2).Betrachte fur ein aquidistantes Gitter Ih ein stabiles Einschrittverfahren mit Inkrementfunktionφ ∈ Cp+1(I2 × Rm,Rm) und Konsistenzordnung p. Es sei uh(0) = u0. Dann existiert eineFunktion e0 ∈ C2(I ,Rm) mit e0(0) = 0, so dass die Beziehung

‖uh − (u− hpe0)‖∞ = O(hp+1)

gilt.

Beweis: Sei Th der dem Einschrittverfahren zugeordnete diskrete Operator, d.h.

Th[wh](tj) :=wh(0)− u(0) fur j = 01h(wh(tj+1)− wh(tj)− hφ(h, tj , wh(tj))) fur j = 1, . . . , N − 1

33


mit wh : Ih → Rm.

Es genugt zu zeigen, dass eine Funktion e0 ∈ C2(I ,Rn), e0(0) = 0 mit

‖Th(u+ hpe0)∣∣Ih‖∞ = O(hp+1) (?)

existiert.Wegen der Stabilitat des ESVs folgt namlich

‖uh − (u− hpe0)|Ih‖∞ ≤ C‖Th(uh)− Th(u− hpe0)

∣∣Ih

)‖∞= C‖Th(uh − hpe0)|Ih

‖∞= O(hp+1).

Zum Beweis von (?) definieren wir e0 als Losung des Anfangswertproblems

e′0(t) =∂

∂wφ(0, t, u(t))e0(t)−

1(p+ 1)!

u(p+1)(t) +1p!∂p

∂kpφ(0, t, u(t)),

e0(0) = 0.

Fur j = 0 gilt‖Th((u+ hpe0)

∣∣Ih

)‖∞ = uh(0)− u(0) = 0(= O(hp+1))

Fur j = 1, . . . , N − 1 ergibt sich nach dem Satz von Taylor

Th((u+ hpe0)∣∣Ih

)(tj) =1h

(u(tj+1)− u(tj) + hpe0(tj+1)− hpe0(tj)

− hφ(h, tj , u(tj) + hpe0(tj)))

=p+1∑l=1

hl−1

l!u(l)(tj) +O(hp+1) +

hp

h(he′0(tj) +O(h2))

−p∑

l=0

hl

l!∂lφ

∂kl(0, tj , u(tj) + hpe0(tj))

=p−1∑l=0

hl

(u(l+1)(tj)(l + 1)!

− 1l!∂lφ

∂kl(0, tj , u(tj))

)+O(hp+1)

+ hp 1(p+ 1)!

u(p+1)(tj) + h0φ(0, tj , u(tj))−1p!∂pφ

∂kp(0, tj , u(tj))hp

+ hpe′0(tj)− φ(0, tj , u(tj) + hpe0(tj))1.18(ii)

= hp 1(p+ 1)!

u(p+1)(tj) + φ(0, tj , u(tj))−1p!∂pφ

∂kp(0, tj , u(tj))hp

+ hpe′0(tj)− φ(0, tj , u(tj))−∂

∂wφ(0, tj , u(tj))hpe0(tj) +O(h2p)

+O(hp+1)

p≥1= hp

(u(p+1)(tj)(p+ 1)!

− 1p!∂pφ

∂kp(0, tju(tj)) + e′0(tj)

− ∂

∂wφ(0, tj , u(tj))e0(tj)

)+O(hp+1)

Wahl von e0= O(hp+1),

Daher folgt der Satz.

Satz 1.52 kann verfeinert werden, indem man weitere Fehlerfunktionen e1, e2, e3, . . . bestimmt.Wir benotigen

34


Satz 1.53. Fur p ∈ N sei u ∈ C(p+2)(I ,Rm) die Losung von (1.1),(1.2). Weiter gelten alleVoraussetzungen aus Satz 1.52 und e0 ∈ C2(I ,Rm) sei wie in Satz 1.52. Dann existiert eineFunktion e1 ∈ C3(I ,Rm) mit e1(0) = 0, so dass die Beziehung

‖uh − (u− hpe0 − hp+1e1)∣∣Ih‖∞ = O(hp+2)

gilt.

Beweis: Ubung 8.

Um mit den Satzen 1.52/1.53 eine adaptive Schrittweitensteuerung zu konstruieren, bemerkenwir zunachst, dass fur alle t ∈ Ih/2(⊇ Ih) gilt:

|uh(t+ h)− u(t+ h)| = hp|e0(t+ h)|+ hp+1|e1(t+ h)|+O(hp+2)

bzw.

|uh2(t+ h)− u(t+ h)| =

(h

2

)p

|e0(t+ h)|+(h

2

)p+1

|e1(t+ h)|+O(hp+2).

Im Folgenden beschranken wir uns auf den Fall m = 1. (Der Fall m > 1 ist komponentenweisezu behandeln.)Sei nun fur t ∈ Ih die Funktion v die exakte Losung zu

v′(s) = f(t, v(s)) (s > t),

v(t) = uj (= Approximation zum Zeitpunkt t = tj)

Wir erhalten (wegen der Schittinvarianz der Losung v) die Fehlerdarstellungen

uh(t+ h)− v(t+ h) = hpe0(t+ h) + hp+1e1(t+ h) +O(hp+2)

= hp(e0(t)︸︷︷︸=0

+he′0(t) + he1(t)︸︷︷︸=0

+h2e′1(t)) +O(hp+2)

= hp+1(e′0(t)) +O(hp+2)

bzw.

uh2(t+ h)− v(t+ h) =

(h

2

)p

he′0(t) +O(hp+2)

Aus den beiden Fehlerdarstellungen ergibt sich durch Auflosen nach v(t+ h)

uh(t+ h)− uh2(t+ h) =

(hp+1 −

(h

2

)p

h

)e′0(t) +O(hp+2)

⇒ he′0(t) =1

2p − 1

(h

2

)−p

(uh(t+ h)− uh2(t+ h)) +O(h2)

Einsetzen in die Fehlerdarstellung fur uh2

liefert

uh2(t+ h)− v(t+ h) =

(h

2

)p(h2

)−p 12p − 1

(uh(t+ h)− uh2(t+ h)) +O(hp+2)

⇒ uh2(t+ h)− v(t+ h) =

12p − 1

(uh(t+ h)− uh2(t+ h))︸︷︷︸

=:∆h

+O(hp+2)

35


Der Term ∆h ist eine Schatzung fur den Fehler uh2(t + h) − v(t + h), der nur berechenbare

Grossen benutzt. Es sei ∆h die relative Schatzung

∆h :=∆h

max1, |uh|.

Auf der Basis von ∆h geben wir (in Pseudocode) einen selbstadaptiven Algorithmus an ((h, h2 )-

Gittersteuerung):Initialisierung:

Startschrittweite: h0 ∈ [0, T ]Maximale/Minimale Schrittweite: hmin < h0 < hmax

Fehlertoleranz: d > 0Vergrosserungs-/Verkleinerungsfaktoren: κ−1

min, κmax > 1Verfahrensordnung: p ∈ N

t := 0, u := u(t0), h := h0.∆h := d+ 1while (t < T ) dowhile(∆h > d) dov := u+ hφ(h, t, u);z := u+ h

2φ(

h2 , t, u

);

w := z + h2φ(

h2 , t+ h

2 , u);

∆h := 12p−1

1max1,|uk|(|v − w|);

h := maxhmin, κmink; if (h = hminexit(0); %Schleife fur Schatzung%u := w;t := t+ h;hi = minhmax, κmaxh;if(h = hmax)exit(0); %Zeitschleife%

36

2 Randwertprobleme fur gewohnlicheDifferentialgleichungen

2.1 Beispiele und Definitionen

Beispiel 2.1 (Transport durch Diffusion und Reaktion). In dem Beispiel 1.6 haben wir einSystem von Differentialgleichungen zur Beschreibung der Reaktion

A+B → C

hergeleitet. Fur die Konzentration a = a(t) und b = b(t) gilt:

a′(t) = g(a(t), b(t)) = −ka(t)b(t) (t ∈ (0,∞))b′(t) = g(a(t), b(t))

Wir hatten vorausgesetzt, dass die Spezies A,B raumlich homogen verteilt sind.Jetzt soll der raumliche Transport berucksichtigt werden.Annahmen:

• A und B liegen in einer geraden Rohre der Lange L vor.

• Es kann nur ein Transport in Langsrichtung der Rohre stattfinden (x-Richtung).

• F sei die Querschnittflache der Rohre.

Die gesuchten Konstanten a, b sind Funktionen von x und t

a = a(x, t), b = b(x, t).

Die Massenbilanz der Spezies A in einem beliebigen, aber festen Volumen ∆V = F ·∆ξ und ineinem beliebigen, aber festen Zeitraum ∆t(= t+ ∆t− t) ist

F

∫ ξ+∆ξ

ξa(x, t+ ∆t)dx ≈ F

∫ ξ+∆ξ

ξa(x, t)dx+ F∆t

∫ ξ+∆ξ

ξg(a(x, t), b(x, t))dx

+ FJa(ξ, t)∆t− FJa(ξ + ∆ξ, t)∆t

Dabei sei Ja = Ja(ξ, t) der Materialdurchfluss von A in (positive) x-Richtung pro Zeitintervallder Lange 1.Wir ersetzen jetzt ≈ durch = und dividieren durch ∆t und ∆ξ:

F

∆ξ

∫ ξ+∆ξ

ξ

a(x, t+ ∆t)− a(x, t)∆t

dx =F

∆ξ

∫ ξ+∆ξ

ξg(a(x, t), b(x, t))dx

+F

∆ξ(Ja(ξ, t)− Ja(ξ + ∆ξ, t))

Wir nehmen an, dass a, b ∈ C2([0, L]× [0,∞]; [0,∞]) gilt. Der Grenzubergang ∆t→ 0 liefert

1∆ξ

∫ ξ+∆ξ

ξat(x, t)dx =

1∆ξ

∫ ξ+∆ξ

ξg(a(x, t), b(x, t))dx+

1∆ξ

(Ja(ξ, t)− Ja(ξ + ∆ξ, t))

37

2 Randwertprobleme fur gewohnliche Differentialgleichungen

Der Grenzubergang ∆ξ → 0 ergibt

at(ξ, t) = g(a(ξ, t), b(ξ, t))− Jax (ξ, t).

Mit derselben Argumentation fur die Spezies B erhalten wir insgesamt

at = g(a, b)− Jax

bt = g(a, b)− Jbx

in (0, L) × (0,∞). Dabei ist Jb der Materialfluss der Spezies B in (positiver) x-Richtung. Umein in a, b geschlossenes System zu erhalten, benotigen wir einen Zusammenhang zwischen Ja,b

und a, b.Falls raumlicher Transport nur durch Diffusion stattfindet, gilt das Fick’sche Gesetz(Proportionalitat zum Konzentrationsgradienten) (da, db > 0):

Ja(x, t) = −daax(x, t), Jb(x, t) = −dbbx(x, t)

Damit ergibt sichat = daaxx + g(a, b)

bt = dbbxx + g(a, b)

in (0, L)× (0,∞). Dies ist ein System partieller Differentialgleichungen (genauer vom paraboli-schen Typ).Was sind sinnvolle Anfangs- und Randbedingungen?Anfangsbedingungen:

a(·, 0) = a0, b(·, 0) = b0 (a0, b0 : [0, L] → [0,∞))

Randbedingungen in x = 0, x = L:

a(0, t) = a0(t), a(L, t) = aL(t)

b(0, t) = b0(t), b(L, t) = bL(t),

wobei a0, b0, aL, bL : (0,∞) → [0,∞) sind (Dirichletrandbedingungen).Falls kein Materialfluss uber den Rand x = 0, L stattfindet (geschlossener Behalter), gilt

ax(0, t) = ax(L, t) = bx(0, t) = bx(L, t) = 0 (t ≥ 0) (Neumannrandbedingugen).

In einem abgeschlossenen Behalter ergibt sich also folgendes Anfangsrandwertproblem:Finde a, b ∈ C2([0, L]× [0,∞); [0,∞)) mit

at = daaxx + g(a, b), bt = dbbxx + g(a, b) in (0, L)× (0,∞),

a(·, 0) = a0, b(·, 0) = b0 in (0, L),

ax(0, t) = ax(L, t) = bx(0, t) = bx(L, t) = 0.

Um ein Randwertproblem zu erhalten, betrachten wir das System fur t → ∞. Dann sollten dieSpezies im Gleichgewicht sein, d.h.

at(·,∞) = 0, bt(·,∞) = 0 in (0, L).

Aus dem Differentialgleichungssystem und den Randbedingungen folgt (unter Vernachlassigungder Zeitvariable)

−daaxx = g(a, b), −dbbxx = g(a, b) in (0, L)ax(0) = ax(L) = bx(0) = bx(L) = 0.

(?)

Dieses Problem heisst (elliptische) Randwertaufgabe. Hier erhalten wir ein System gewohnlicherDifferentialgleichungen.

38

2.1 Beispiele und Definitionen

Beispiel 2.2 (Warmeleitungsgleichung). Wir betrachten einen Stab der Lange L mit Tempera-turverteilung T = T (x, t). Falls die Temperatur nur in x-Richtung variabel ist, erhalten wir mitder Argumentation wie in Beispiel 2.1 fur die Temperatur T das Anfangswertproblem

Tt = Txx in (0, L)× (0,∞),

T (·, 0) = T0 ∈ C2([0, L]) in (0, L),

T (0, t) = T 0(t), T (L, t) = TL(t) (t ∈ [0,∞]).

Dabei seien auch T 0, TL ∈ C2([0,∞]). In diesem Fall sind die Dirichletbedingungen realistischer,da Temperaturen in einzelnen Punkten gut experimentell realisiert werden konnen.Fur den asymptotischen Fall t→∞ gilt das Randwertproblem

Txx = 0 in (0, L)

T (0) = T 0 T (L) = TL

fur die Temperatur T = T (x). Dabei seien T 0, TL ≥ 0 vorgegeben. Die Losung ist

T = T (x) =TL − T 0

Lx+ T 0.

Randwertaufgaben:Diffusionstransport in einer Raumdimension: Finde Konzentrationen

a : I = (0, 1)× [0,∞) → [0,∞)

mitat − ((d(a)ax)x = g(a) in I × (0,∞),

wobei a(·, 0) = a0 : I → [0,∞] in I.Diffusionstransport in meheren Raumdimensionen (Ω ⊂ Ra offen und beschrankt):Finde a : Ω× [0,∞) → [0,∞) mit

at − div(d(a)∇a) = g(a) in Ω× (0,∞),

wobei a(·, 0) = a0 : Ω → [0,∞] in Ω.∂

∂na(x, t) = 0 x ∈ ∂Ω

Annahme: Es existiert genau eine Losung a ∈ C2(I×[0,∞), [0,∞)) mit a = a(x) := limt→∞ a(·, t)exisitert, limt→∞ at(·, t) = 0.Daraus folgt:

(d(a)ax)x = g(a) in I,

ax(0) = ax(1) = 0

Allgemeines Transportproblem:

at − Jx = g(x) in I

i. Diffusion: J = d(a)ax

ii. Transport: J = va (v ∈ R sei die Transportgeschwindigkeit)

iii. Diffusion und Transport: J = va+ d(a)ax

Draus ergibt sich

i. at − (d(a)ax)x = 0

ii. at + vax = 0Eine Losung hierzu ist durch a0 = a0(x), a = a(x, t) = a0(x− vt) gegeben.

iii. at − (d(a)ax)x + vax = 0

39


2.2 Lineare Randwertaufgaben/Sturm-Liouville Problem

Es existiert keine allgemeine Losungstheorie fur nichtlineare Randwertprobleme. Wir betrachten(scheinbar spezielle) lineare Probleme.Seien a0, a1, h ∈ C1(I ,R) fur I := (a, b). Finde u ∈ C2(I) mit

u′′ + a1u′ + a0u = h in I, (2.1)

R1u := r11u(a) + r12u′(a) = s1

R2u := r21u(b) + r22u′(b) = s2

(2.2)

Dabei seien r11, r12, r21, r22 ∈ R mit (r11, r12), (r21, r22) 6= (0, 0), s1, s2 ∈ R.Die Gleichung 2.1 lasst sich als System erster Ordnung umschreiben:

w′ = −(a1w + a0u) + h, u′ = w

Seien u1, u2 die der u-Komponente entsprechenden Eintrage von zwei linear unabhangigenLosungen des zugehorigen homogenen Systemes (h ≡ 0, Existenz ist nach Satz 1.12 klar).

Satz 2.3. Das Randwertproblem (2.1),(2.2) ist genau dann eindeutig losbar, wenn

det(r11 r12r21 r22

)6= 0

gilt.

Beweis: Sei up eine spezielle Losung von (2.1). Dann existieren nach (1.12) Konstanten c1, c2,so dass fur jede Losung u von (2.1) die Darstellung

u(x) = up + c1u1(x) + c2u2(x) (x ∈ I)

gilt. Um die Randbedingungen (2.2) zu erfullen, muss die Bedingung

Rku = Rkup + c1Rku1 + c2Rku2 (k = 1, 2)

gelten. Das ist aquivalent zu(r11 r12r21 r22

)·(c1c2

)=(s1 −R1up

s2 −R2up

).

Daraus folgt die Aussage.

Das Problem (2.1),(2.2) heisst Sturm-Liouville Problem, falls Funktionen

p ∈ C1(I), q ∈ C0(I)

mit den Eigenschaftenp(x) > 0, q(x) ≥ 0

existieren, so dass (2.1) in der Form

−(pu′)′ + qu = h in I (2.3)

40

2.3 Das klassische Finite Differenzen Verfahren

geschrieben werden kann.Die DGL (2.3) ist deshalb von besonderer Bedeutung, weil sie Variationsstruktur besitzt:Wir definieren das Funktional F : C1(I) → R vermoge

F [v] :=12

∫ b

ap(x)(v′(x))2dx+

12

∫ b

aq(x)(v(x))2dx−

∫ b

ah(x)v(x)dx

und betrachten folgende Variationsaufgabe:Finde u ∈ C2(I) mit

F [u] ≤ F [v] ∀v ∈ C2(I) (2.4)

Wodurch ist nun die Losung u von (2.4) ausgezeichnet (wenn sie uberhaupt existiert)?

Lemma 2.4 (Euler-Lagrange Gleichung). Sei L = L(p, z, x) ∈ C1(R3,R). Es existiere eineFunktion u ∈ C2(I) mit ∫

IL(u′(x), u(x), x)dx ≤

∫IL(v′(x), v(x), x)dx

fur alle v ∈ C1(I). Dann gilt fur u

−[Lp(u′(x), u(x), x)]′ + Lz(u′(x), u(x), x) = 0 (x ∈ I).

Beweis: Ubung.

Wenden wir Lemma 2.4 auf das Funktional F an, erhalten wir

−[p(x) · v′(x)]′ + q(x) · v(x) = h(x) (x ∈ I).

Das ist das Sturm-Liouville Problem.Die Variationsstruktur werden wir ausnutzen, um die Finite-Elemente Methode (FEM) zu mo-tivieren.In der Analysis und Numerik von Randwertaufgaben ist es oft wichtig, dass nur Ableitungenhochster Ordnung in linearer Weise auftreten. Etwas allgemeiner als die Sturm-Niouville Glei-chung ist die Sturm Gleichung

−u′′ = f(x, u, u′) in I

fur f ∈ C0(R3,R). Diese losen wir jetzt mit der Finite-Differenzen-Methode (FDM).


Wir betrachten das Sturm’sche Problem

−u′′ = f(x, u, u′) in I := (a, b) (2.5)

mit Dirchlet-Randbedingungen

u(a) = ua, u(b) = ub fur ua, ub ∈ R. (2.6)

Wir nehmen an, dass (2.5),(2.6) eine Losung u ∈ C2(I) besitzt.Es sei Ih := x0, . . . , xN ein aquidistantes Gitter zu I mit

a =: x0 < x1 < x2 < . . . < xN−1 < xN := b.

Das heisst:xi = a+ ih, h =

b− a

N.

41


Auf Ih betrachten wir die Differenzenquotienten

u′(xi) ≈12h

(u(xi+1)− u(xi−1)), i = 1, . . . , N − 1,

u′′(xi) ≈1h2

(u(xi−1)− 2u(xi) + u(xi+1)), i = 1, . . . , N − 1.

Einsetzen in (2.5) liefert fur die Unbekannten u1, u2, . . . , uN−1 die Gleichungen

− 1h2

(ui−1 − 2ui + ui+1) = f

(xi, ui,

12h

(ui+1 − ui−1))

(i = 1, . . . , N − 1). (2.7)

Das Verfahren (2.7) zur Bestimmung von u0, . . . , uN mit u0 = ua, uN = ub heisst Finite-Differenzen-Methode (FDM) fur (2.5),(2.6).Die FD-Approximation uh auf Ih zu der Losung u von (2.5),(2.6) ist durch uh(xi) = ui definiert,falls eine Losung (u0, . . . , uN ) von (2.7) existiert.

Bemerkung 2.5. Um uh gemaß (2.7) zu bestimmen, ist ein (N − 1)× (N − 1)-dimensionales(nichtlineares) Gleichungssystem zu losen. Es ist fur (2.5),(2.6) nicht klar, ob eine Losung exis-tiert und eindeutig ist.

Wir konnen (2.7) in Funktionalschreibweise umformulieren, indem der zugehorige diskrete Ope-rator Th eingefuhrt wird:

Th :Xh → Xh

wh 7→ Th[wh]

wobei Xh := wh : Ih → R und

Th[wh] :=

wh(x0)− ua i = 0− 1

h2 [wh(xi−1)− 2wh(xi) + wh(xi+1)]−f(xi, wh(xi), 1

2h(wh(xi+1)− wh(xi−1)) i = 1, . . . , N − 1wh(xN )− ub i = N

ist. Dann ist (2.7) aquivalent zu Th[uh] = 0.

Definition 2.6 (Konsistenz, Stabilitat, Konvergenz).

i. Die FDM (2.7) heisst konsistent zur Ordnung p ∈ N genau dann, wenn gilt:

‖Th[u∣∣Ih

]‖∞ = O(hp)

ii. Die FDM (2.7) heisst stabil genau dann, wenn eine Konstante C > 0 existiert mit

‖wh − wh‖∞ ≤ C‖Th[wh]− Th[wh]‖∞ ∀wh, wh ∈ Xh

iii. Die FDM (2.7) heisst konvergent zur Ordnung p ∈ N genau dann, wenn fur den globalenFehler

eh = maxi=0,...,N

|uh(xi)− u(xi)|

die Beziehung eh = O(hp) gilt.

Satz 2.7 (Konsistenzsatz). Sei f = f(x, v, w) ∈ C(I ×R2,R) mit fw ∈ C(I ×R2,R). Dann istdie FDM (2.7) fur u ∈ C(4)(I) konsistent von der Ordnung 2.

42


Beweis: Wir mussen ‖Th[u∣∣Ih

]‖∞ = O(h2) zeigen.Fur i = 1, . . . , N − 1 gilt:

|Th[u∣∣Ih

](xi)| ≤∣∣∣∣− 1h2

(u(xi+1)− 2u(xi) + u(xi−1)) + u′′(xi)∣∣∣∣

+∣∣∣∣−u′′(xi)− f

(xi, u(xi),

12h

(u(xi+1)− u(xi−1))∣∣∣∣

≤ c1h2‖u(4)‖L∞(I) +

∣∣∣∣f(xi, u(xi), u′(xi))− f

(xi, u(xi),

12h

(u(xi+1)− u(xi−1)))∣∣∣∣

≤ c1h2‖u(4)‖L∞(I) +

∣∣f(xi, u(xi), u′(xi))− f(xi, u(xi), u′(xi) +Di(h)

)∣∣Dabei gilt Di(h) = O(h2) und es folgt wieder mit Taylorentwicklung

|Th[u∣∣Ih

](xi)| ≤ c1h2‖u(4)‖L∞(I) + c2c3h

2.

Dabei ist c3 so gewahlt, dass |Di(h)| < c3h2 gilt und c2 von f abhangt. Somit ist

|Th[u∣∣Ih

](xi)| = O(h2)

Im nachsten Schritt soll die Konvergenz der FDM (2.7) fur das einfachere Sturm-Liouville Pro-blem bewiesen werden:

−(pu′)′ + qu = h in I

⇔ −pu′′ − p′u′ + qu = h in I

p>0⇔ −u′′ − p′

pu′ +

q

pu = h in I

⇔ −u′′ + pu′ + qu = h in I

Als FDM (2.7) ergibt sich:

ui+1 − 2ui + ui−1

h2+ p(xi)

ui+1 − ui−1

2h+ q(xi)ui = h(xi) (i = 1, . . . , N − 1).

Der Operator Th lasst sich fur wh ∈ Xh dann in der Form

Th[wh] = Ahwh + rh

schreiben. Dabei ist:

Ah =

2h2 + q(x1) − 1

h2 + p(x1)2h 0 · · · 0

− 1h2 − p(x2)

2h2h2 + q(x2) − 1

h2 + p(x2)2h

...

0. . . . . . . . . 0

... − 1h2 − p(xN−2)

2h2h2 + q(xN−2) − 1

h2 + p(xN−2)2h

0 · · · 0 − 1h2 + p(xN−1)

2h2h2 + q(xN−1)

und

rh =(−h(x1)−

(1h2

+p(x1)2h

)ua, h(x2), . . . ,−h(xN−2), h(xN−1)−

(1h2− p(xN−1)

2h

)ub

)T

Th ist ein affiner Operator.Annahme:

43


i. detA 6= 0 fur h > 0

ii. ∃C > 0 : ‖A−1h ‖∞ < C fur h > 0

In (ii) sei ‖.‖∞ eine mit dem Vektornorm ‖.‖∞ vertragliche Matrixnorm.

2.4 Matrizenalgebra und Stabilitatssatz

Um die Stabilitat der FDM nachzuweisen, benotigen wir einige (abstrakte) Resultate der Ma-trizenalgebra.

Definition 2.8 (Halbordnungen auf Rm und Rm×m). Seien u, v ∈ Rm und A = (aij), B =(bij) ∈ Rm×m.

i. Wir schreiben u ≤ v (u < v) genau dann, wenn

ui ≤ vi (ui < vi) ∀i ∈ 1, . . . ,m

gilt.

ii. Wir schreiben A ≤ B (A < B) genau dann, wenn

aij ≤ bij , (aij < bij) ∀i, j ∈ 1, . . . ,m

gilt.

iii. die Matrix heisst nichtnegativ (oder monoton) genau dann, wenn 0 ≤ A gilt.

Lemma 2.9 (Charakterisierung der Monotonie). Fur A = (aij) ∈ Rm×m sind aquivalent:

i. A ist nichtnegativ,

ii. u, v ∈ Rm, u ≤ v ⇒ Au ≤ Av,

iii. v ∈ Rm, 0 ≤ v ⇒ 0 ≤ Av.

Beweis: (i) ⇒ (ii): Seien u, v ∈ Rm mit u ≤ v. Fur i = 1, . . . ,m gilt

(Au)i =m∑

j=1

aijuj

u≤v,aij≥0

≤m∑

j=1

aijvj ≤ (Av)i

⇒ Au ≤ Av

(ii) ⇒ (iii): Wahle u = 0.(iii) ⇒ (i): Sei ej ∈ Rm der j-te Einheitsvektor. Wegen (ii) und 0 ≤ ej gilt

A · 0 = 0 ≤ A · ej .

Aej ist der j-te Spaltenvektor von A. Da j beliebig war, folgt aus Definition 2.8, dass 0 ≤ aij .Das ist gerade (i).

Definition 2.10. Sei e ∈ Rm mit 0 < e. Die Norm ‖ · ‖e : Rm → R≥0 mit

‖u‖e := maxj=1,...,m

|uj |ej

heisst gewichtete Maximumnorm.

44


Sei zum Beispiel e = (1, . . . , 1)T . Dann ist ‖ · ‖e die Maximumnorm ‖ · ‖∞.

Lemma 2.11. Sei e ∈ Rm, 0 < e und A ∈ Rm×m. Fur die assoziierte Matrixnorm

‖ · ‖e : Rm×m → [0,∞)

mit‖A‖e = sup‖Au‖e

∣∣u ∈ Rm, ‖u‖e ≤ 1

gilt‖A‖e = ‖|A| · e‖e

mit(|A|)ij = |aij | (i, j = 1, . . . ,m).

Falls 0 ≤ A gilt, ist ‖A‖e = ‖Ae‖e.

Beweis: Sei

E = diag(e−1i ) =

1e1

0. . .

0 1em

.

Dann gilt die Darstellung‖u‖e = ‖Eu‖∞ (u ∈ Rm)

also

‖A‖e = sup‖Au‖e

∣∣u ∈ Rm, ‖u‖e ≤ 1= sup‖EAu‖∞

∣∣u ∈ Rm, ‖ Eu︸︷︷︸=:v

‖∞ ≤ 1

= sup‖EAE−1v‖∞∣∣v ∈ Rm, ‖v‖∞ ≤ 1

= ‖EAE−1‖∞

= maxi=1,...,m

m∑j=1

|aij |ejei

= ‖|A|e‖e.

Der Zusatz gilt wegen aij ≥ 0 fur 0 ≤ A.

Definition 2.12 (inversmonotone Matrizen). Eine Matrix A ∈ Rm×m heisst inversmonotongenau dann, wenn gilt

i. detA 6= 0

ii. A−1 ist monoton.

Lemma 2.13 (Charakterisierung der Inversmonotonie). Fur A = (aij) ∈ Rm×m sind aquivalent:

i. A ist inversmonoton.

ii. u, v ∈ Rm, Au ≤ Av ⇒ u ≤ v

iii. v ∈ Rm, 0 ≤ Av ⇒ 0 ≤ v.

Beweis: Analog zum Beweis von 2.9.

Der Begriff der Inversmonotonie ermoglicht jetzt eine Abschatzung der Norm einer Matrixinver-sen!

45


Satz 2.14 (Normabschatzung). Sei A ∈ Rm×m inversmonoton und e ∈ Rm mit 0 < e. Weitergelte fur ein c ∈ (0,∞) die Abschatzung

ce ≤ Ae.

Dann gilt

‖A−1‖e ≤1c.

Beweis: ce ≤ Ae. Da A−1 monoton ist, folgt aus Lemma 2.9(ii)

cA−1e = A−1(ce) ≤ A−1(Ae) = e (?)

Anwendung von Lemma 2.11 (Zusatz) liefert (0 ≤ A−1 monoton)

‖A−1‖e = ‖A−1e‖e ≤∥∥∥∥1ce

∥∥∥∥e

=1c‖e‖e =

1c

Die Inversmonotonie einer Matrix ist schwer zu beweisen. Wir fuhren noch ein Kriterium ein.

Definition 2.15 (M-Matrix). Eine Matrix A = (aij) ∈ Rm×m heisst M-Matrix genau dann,wenn gilt:

i. A ist inversmonoton

ii. aij ≤ 0 fur i, j ∈ 1, . . . ,m, i 6= j.

Satz 2.16 (M-Kriterium). Sei A = (aij) ∈ Rm×m mit aij ≤ 0 fur i, j ∈ 1, . . . ,m, i 6= jgegeben. Dann gilt: Existiert ein e ∈ Rm mit 0 < e und 0 < Ae. Dann ist A M-Matrix.

Beweis: Zerlege A = D −B mit D = diag(aii). Nach Voraussetzung gilt

0 ≤ B. (?)

Daher ist 0 ≤ Be wegen 0 < e und Lemma 2.9(iii). Nach Definition von e ist 0 < Ae = De−Be,also

0 ≤ Be < De.

Wir erhalten 0 < De und mit 0 < e folgt

aii > 0 (i = 1, . . . ,m).

Setze nun P = D−1︸︷︷︸≥0

· B︸︷︷︸≥0

≥ 0. Es gilt

A = D · (I − P )

undPe < e.

Aus Lemma 2.11 (Zusatz) ergibt sich

‖P‖e = ‖Pe‖e

Pe<e≤ ‖e‖e = 1.

46


Nach dem Storungslemma existiert dann (I − P )−1 und erlaubt die Darstellung

(I − P )−1 =∞∑

j=0

P j .

Offenbar ist mit 0 ≤ P auch 0 ≤∑∞

j=0 Pj , also

0 ≤ (I − P )−1

⇒ 0 ≤ (I − P )−1 ·D−1 = A−1

Betrachte nun wieder (2.3) und (2.6). Diese waren aquivalent zu

−u′′ + pu′ + q = h in I

u(a) = ua, u(b) = ub

Die FDM war gegeben durch

−ui+1 − 2ui + ui−1

h2+ p(xi)

ui+1 − ui−1

2h+ q(xi)ui = h(xi), i = 1, . . . , N − 1

u0 = ua, uN = ub

Dies war aquivalent zu

Th[uh] = Ahuh + rh!= 0, uh = (u0, . . . , uN )T

Wir wissen fur den Konsistenzfehler (Satz 2.7):

‖Th[u∣∣Ih

]‖∞ = O(h2)

Fragen:

i. Existiert die Losung uh, d.h. gilt det(Ah) 6= 0?

ii. Ist das Verfahren stabil?

Satz 2.16 besagte:Sei A = (aij) ∈ Rm×m, aij ≤ 0 fur i 6= j gegeben. Dann gilt:

∃e ∈ Rm, 0 < e, 0 < Ae⇒ A ist M-Matrix.

Ist A eine Matrix, dann ist

i. det(A) 6= 0

ii. ∃c > 0 : (0 ≤ Ae− ce⇒ ‖A‖−1e ≤ 1

c )

Satz 2.17 (Konvergenz FDM). Wir betrachten die Sturm-Liouville Gleichung (2.3) mit Dirich-letrandbedingungen (2.6). Weiter gelte p, q > 0 auf I und u ∈ C4(I) sei die eindeutige Losungvon (2.3),(2.6). Dann gilt

i. Es existiert ein h0 > 0, so dass die FDM (2.7), d.h.

Th[uh] = Ahuh + rh = 0

fur h ∈ (0, h0) genau eine Losung uh ∈ Xh besitzt.

47


ii.‖u∣∣Ih− uh‖∞ = O(h2).

Beweis:

i. Zum Beweis mussen wir zeigen, dass Ah fur h hinreichend klein invertierbar ist. Fur hhinreichend klein gilt wegen p > 0 und p ∈ C1(I) die Abschatzung

− 1h2± p(xi)

2h≤ − 1

h2+‖p‖L∞

2h= − 1

h2+

12h

∥∥∥∥−p′p∥∥∥∥

L∞

≤ − 1h2

+12h

‖p′‖L∞

infx∈I(p(x))< 0 (fur h hinreichend klein)

Die Außerdiagonalelemente von Ah sind negativ.Fur e = (1, . . . , 1)T gilt 0 < e und

Ahe =(

1h2

+ q(x1) +p(x1)2h

, q(x2), . . . , q(xN−2),1h2

+ q(xN−1)−p(xN−1)

2h

)T

Wegen q > 0 (q = qp) folgt (mit h hinreichend klein), dass 0 < Ahe gilt. Insgesamt ist dann

nach Satz 2.16 Ah eine M-Matrix. Damit folgt (i), da M-Matrizen invertierbar sind.

ii. Fur den Beweis von (ii) betrachte fur h hinreichend klein

‖eh‖∞ = ‖uh − u∣∣Ih‖∞

= ‖A−1h Ah(uh − u

∣∣Ih

)‖∞≤ ‖A−1

h ‖∞‖Ah(uh − u∣∣Ih

)‖∞≤ ‖A−1

h ‖∞‖Ahu∣∣Ih

+ rh − (Ahuh + rh)‖∞= ‖A−1

h ‖∞‖Th[u∣∣Ih

]− Th[uh]︸︷︷︸=0 (nach (i))

‖∞

(2.7)= ‖A−1

h ‖∞O(h2)

(??)

Es gilt fur B ∈ Rm×m

‖B‖∞ = max‖x‖∞≤1

‖Bx‖∞

= maxi=1,...,m

m∑j=1

|bij |

(Zeilensummennorm)

und fur e = (1, . . . , 1)T

‖B‖e = max‖x‖e≤1

‖Bx‖e = maxmaxi

|xi|ei≤1

(max

i

|(Bx)i|ei

)= max

‖x‖∞≤1‖Bx‖∞

Also ist hier ‖ · ‖e = ‖ · ‖∞.Nach Konstruktion existiert ein c > 0, so dass c · e ≤ Ahe fur e = (1, . . . , 1)T gilt. NachSatz 2.14 ist dann ‖A−1

h ‖∞ = ‖A−1h ‖e ≤ 1

c . Also gilt insgesamt

‖eh‖∞ = O(h2).

48


Exkurs - Tsumami:

Gleichung der Hydrodynamik:Sei % = %(x, y, t) > 0 die Dichte, u die Geschwindigkeit in x-Richtung, v die Geschwindigkeit iny-Richtung und g die Graviationskonstante. Dann gilt in Ω(t)× (0, T ):

%t + (%u)x + (%v)y = 0

(%u)t + (%u2)x + (%uv)y + p(%)x = ∆u

(%v)t + (%uv)x + (%u2)y + p(%)x = ∆v − g%

Flachwassergleichung fur die Hohe h = h(x, t) ≥ 0 und die Geschwindigkeit in x-Richtungu = u(x, t) in I × (0, T ).

ht + ux = 0

ut +(u2

h+

12gh2

)= −ghzx(x)

Hierbei beschreibt die Funktion zx das Bodenprofil. Mit ω = (h, u)T und f(u) = (u, h2

u + 12jh

2)T

folgtωt + f(ω)x = 0 in I × (0, T ) (Erhaltungssatze)

49


50

3 Die Finite-Elemente Methode in 1D

3.1 Das Sturm-Liouville Problem als Variationsproblem

Wir betrachten das Sturm-Liouville Problem

−(pu′)′ + qu = h in I = (0, 1) (3.1)

mitu(0) = 0, u′(1) = 0. (3.2)

Dabei moge fur p, q, h : I → R gelten

p, q ∈ L∞(I), h ∈ L2(I), p = p(x) ≥ p0 > 0, q = q(x) ≥ 0 (x ∈ I) (3.3)

Bemerkung 3.1.

i. (3.2) heissen naturliche Randbedingungen. Wir behandeln spater auch allgemeinere Falle.

ii. Die Bedingungen (3.3) sind wesentlich schwacher als die Bedingungen an p, q, h in Kapitel2.

Definition 3.2 (Klassische Losung). Sei p ∈ C1(I), q, h ∈ C0(I). Eine Funktion u ∈ C2(I), die(3.1) punktweise in I und (3.2) erfullt, heisst klassische Losung von (3.1),(3.2).

Dieses klassische Losungskonzept, das Grundlage der Finite-Differenzen Methode in Kapitel 2war, ist unter den Bedingungen (3.3) nicht anwendbar.Zur Einfuhrung eines ”schwachen“ Losungskonzepts definieren wir

V := v ∈ C1(I)|v(0) = 0.

Das Funktional F : V → R mit

F (v) :=12

∫ 1

0p(x)(v′(x))2dx+

12

∫ 1

0q(x)(v(x))2dx−

∫ 1

0h(x)v(x)dx (3.4)

heisst das zu (3.1) gehorige Funktional. Dem Funktional ist folgendes Variationsproblem zuge-ordnet:Finde u ∈ V mit

F (u) = infv∈V

F (v) (3.5)

Lemma 3.3.

i. Sei u ∈ V eine Losung von (3.5). Dann gilt fur alle φ ∈ V∫ 1

0pu′φ′ dx+

∫ 1

0quφ dx =

∫ 1

0hφ dx. (3.6)

ii. Sei u ∈ V ∩ C2(I) eine Losung von (3.5). Fur p ∈ C1(I), q, h ∈ C0(I) ist u eine klassischeLosung von (3.1),(3.2).

51


Definition 3.4 (schwache Losung).

i. (3.6) heisst schwache Formulierung von (3.1).

ii. Eine Funktion u ∈ V , die (3.6) fur alle φ ∈ V erfullt, heisst schwache Losung von(3.1),(3.2).

Beweis: von Lemma 3.3:

i. Da u eine Losung von (3.5) ist, gilt F (u) ≤ F (v) fur alle v ∈ V , insbesondere gilt

F (u) ≤ F (u+ εφ) =: G(ε) (?)

fur alle ε ∈ R und φ ∈ V . Fur G gilt

G(ε) =12

∫ 1

0p(u′ + εφ′)2 + q(u+ εφ)2dx−

∫ 1

0h(u+ εφ)dx

=12

∫ 1

0p(u′)2 + qu2 − 2hu dx+ εA(u, φ) +

12

∫ 1

0ε2((φ′)2p+ φ2q)dx

A(u, φ) :=∫ 1

0(pu′φ′ + quφ− hφ)dx

Es ist G ∈ C1(R) (als Polynom) und daher

G′(ε) = A(u, φ) +∫ 1

0ε((φ′)2p+ φ2q)dx

⇒ G′(0) = A(u, φ)

Wegen (?) muss G′(0) = 0 gelten, also mit A(u, φ) = 0 fur alle φ ∈ V . Das ist gerade (3.6).

ii. u erfullt die schwache Formulierung (3.6). Wegen u ∈ C2(I) und p ∈ C1(I) ergibt partielleIntegration [

pu′φ]10−∫ 1

0(pu′)′φdx+

∫ 1

0quφdx =

∫ 1

0hφdx

φ(0)=0,φ∈V⇒ p(1)u′(1)φ(1)−∫ 1

0((pu′)′ − qu+ h)φdx = 0 ∀φ ∈ V (??)

Also gilt fur φ ∈ C10(I) ⊂ V ∫ 1

0((pu′)′ − qu+ h)φdx = 0

Nach dem Variationslemma folgt dann, dass (pu′)′ − qu+ h = 0 in I gilt.Es is noch (3.2) zu zeigen. Wegen u(0) = 0 mit u ∈ V muss noch u′(1) = 0 gezeigt werden.Aus (??) und (3.1) ergibt sich

p(1)φ(1)u′(1) = 0 ∀φ ∈ V.

Wegen p(3.3)

≥ p0 > 0 muss entweder φ(1) = 0 oder u′(1) = 0 gelten. Da in V Funktionen ψmit ψ(1) 6= 0 existieren, muss u′(1) = 0 gelten.

Lemma 3.2 zeigt, dass wir statt dem RWP (3.1),(3.2) auch das Variationsproblem (3.5) losenkonnen.

52

3.2 Existenz (und Eindeutigkeit) einer schwachen Losung


Programm:

i. Zeige: Es gilt d := infv∈V F (v) > −∞. Somit existiert eine Folge (vn) ⊂ V mit

F (vn) n→∞−→ d := infv∈V

F (v).

ii. Zeige: Fur (vn) gilt: vn ist eine Cauchyfolge in einem normierten Raum (X, ‖ · ‖) mitX ⊆ V .

iii. Zeige: (X, ‖ · ‖) ist vollstandig. Dann existiert ein x ∈ X mit vnn→∞−→ x in X.

iv. Zeige: F ist ”stetig“, d.h. F (xn) n→∞−→ F (x), falls xnn→∞−→ x in X. Dann ist namlich

F (v) = limn→∞

F (vn) = d.

Hierbei stellt sich gerade der zweite Schritt als relativ problematisch heraus.

Lemma 3.5 (Poincare - Ungleichung). Fur jedes v ∈ V gilt∫ 1

0(v(x))2dx ≤ 1

2

∫ 1

0(v′(x))2dx.

Beweis: Wegen v(0) = 0 gilt

v(x) = v(x)− v(0) =∫ x

0v′(y) dy (x ∈ I)

Sei nun x > 0. Dann folgt mit der Jensen-Ungleichung

(v(x))2 =(∫ x

0v′(y)dy

)2x>0= x2

(1x

∫ x

0v′(y)dy

)2

≤ x2 1x

∫ x

0|v′(y)|2dy

Schließlich folgt∫ 1

0(v(x))2dx ≤

∫ 1

0

(x

∫ x

0(v′(y))2dy

)dx ≤

∫ 1

0

(x

∫ 1

0(v′(y))2dy

)dx =

12

∫ 1

0(v′(y))2dy

Sei nun v ∈ V . Man benotigt im folgenden die Young’sche Ungleichung (a, b ≥ 0, ε > 0):

a · b ≤ εa2 +14εb2 (?)

Beweis:4εab ≤ 4ε2a2 + b2

⇔ 0 ≤ 4ε2a2 − 4εab+ b2

⇔ 0 ≤ (2εa− b)2

53


Dann gilt fur F

F (v) =12

∫ 1

0(p(v′)2 + qv2)dx−

∫ 1

0hvdx

≥ p0

2

∫ 1

0(v′)2dx−

∫ 1

0hvdx

Holder≥ p0

2

∫ 1

0(v′)2dx−

(∫ 1

0h2dx

)1/2(∫ 1

0v2dx

)1/2

(3.5)

≥ p0

2

∫ 1

0(v′)2dx− 1√

2

(∫ 1

0h2dx

)1/2(∫ 1

0(v′)2dx

)1/2

(?)

≥ p0

2

∫ 1

0(v′)2dx− 1√

2

(ε

∫ 1

0nh2dx+

14ε

∫ 1

0(v′)2dx

)ε > 0

=(p0

2− 1√

2 · 4ε

)∫ 1

0(v′(x))2dx− ε√

2

∫ 1

0h2dx

= 0− (p0

√2 · 2)−1

√2

∫ 1

0h2dx (ε = (p0

√2 · 2)−1)

= − 14p0

‖h‖2L2(I)

> −∞

Definition 3.6 (Minimalfolge). Sei (vn) ⊂ V . Dann heisst (vn) Minimalfolge zu F genau dann,wenn

limn→∞

F (vn) = infv∈V

F (v)(=: d).

Um das Programm weiter durchzufuhren, muss jetzt ein normierter Raum gewahlt werden. Wirwahlen (V, ‖ · ‖V ), wobei

‖v‖V :=(∫ 1

0(v′(x))2dx+

∫ 1

0(v(x))2dx

)1/2

Die Frage, warum man nicht (V, ‖ · ‖C1) wahlt, bleibt erstmal offen.

Lemma 3.7. Jede Minimalfolge ist Cauchyfolge in (V, ‖ · ‖V ).

Beweis: Sei (vn)n∈N eine Minimalfolge. Wir definieren

B(v, w) :=∫ 1

0(pv′w′ + qvw)dx (v, w ∈ V )

B ist offensichtlich eine Bilinearform. Wegen (3.3) folgt fur v ∈ V

B(v, v) =∫ 1

0(p(v′)2 + qv2)dx ≥ p0

∫ 1

0(v′)2dx

und

B(v, v) ≥ p0

(23

∫ 1

0(v′)2dx+

13212

∫ 1

0(v′)2dx

)≥ p0

(23

∫ 1

0(v′)2dx+

23

∫ 1

0v2dx

)(Poincare)

=23p0‖v‖2

V

54


Fur zwei beliebige Folgenglieder vr, vs ∈ (vn)n∈N folgt dann

23p0‖vr − vs‖V ≤ B(vr − vs, vr − vs)

= 2B(vr, vr) + 2B(vs, vs)− 4B(vr + vs

2,vr + vs

2

)= 2 · 2

(12

∫ 1

0((v′r)

2p+ qv2r − 2hvr)dx+

12

∫ 1

0((v′s)

2p+ qv2s − 2hvs)dx

)− 4

(22

∫ 1

0

(p

(v′r + v′s

2

)2

+ q

(vr + vs

2

)2

− 2hvr + vs

3

)dx

)

+ 4∫ 1

0(hvr + hvs − h(vr + vs))︸︷︷︸

=0

dx

= 4(

(F (vr) + F (vs)− 2F(vr + vs

2

))Sei d := infv∈V F (v). Aus der Definition der Minimalfolge folgt

23p0‖vr − vs‖V ≤ 4 ((F (vr) + F (vs)− 2d)

r,s→∞−→ 4(d+ d− 2d) = 0

Somit ist (vn)n∈N Cauchyfolge in (V, ‖ · ‖V ).

Beispiel 3.8 ((V, ‖ · ‖V ) ist nicht vollstandig). Betrachte die Funktionenfolge

vn(x) :=

|x| fur x ≤ − 1

n ,12nx

2 + 12n fur − 1

n < x ≤ 1n ,

x fur x > 1n .

Fur vn(x) = vn(x− 12)− 1

2 gilt vn ∈ V .Die Funktionenfolge vnn∈N ist Cauchyfolge bezuglich (V, ‖ · ‖V ):∫ 1

0|vn+1 − vn|2dx =

∫ 1n+1

− 1n+1

∣∣∣∣12(n+ 1− n)x2 +1

2(n+ 1)− 1

2n

∣∣∣∣ ≤ ∣∣∣∣12 · 12

+1n

∣∣∣∣≤ 1

2dx+ 2

∫ − 1n+1

− 1n

∣∣∣∣12nx2 +12n

− |x|∣∣∣∣︸︷︷︸

| 12n

+ 12n− 1

n |≤2

dxn→∞−→ 0

v′n(x) :=

−1 fur x ≤ − 1

n ,nx fur − 1

n < x ≤ 1n ,

1 fur x > 1n .

Eine analoge Rechnung ergibt jetzt ∫ 1

0|v′n+1 − v′n|dx

n→∞−→ 0

Wir beobachten weiter

vn(x) n→∞−→∣∣∣∣x− 1

2

∣∣∣∣− 12

(α ∈ I)

v′n(x) n→∞−→ sign

(x− 1

2

)(x ∈ I\1/2).

55


Falls eine Funktion v ∈ V existiert mit limn→∞

‖v − vn‖V = 0, muss auch

limn→∞

‖v − vn‖L2 = 0

gelten (wegen der Definition von ‖ · ‖V , die ‖ · ‖L2 als Baustein enthalt). Der Limes in L2 isteindeutig, namlich der punktweise Limes von v(α) =

∣∣α− 12

∣∣ − 12 . Diese Funktion ist nicht in

v ∈ C1(I). E

Somit ist (V, ‖ · ‖V ) nicht vollstandig.

3.3 Sobolevraume in 1D

Das Vervollstandigungsproblem aus 3.1 kann durch den Prozess der Vervollstandigung von Vbezuglich ‖ · ‖V zumindest abstrakt gelost werden. Leider erkennt man dann die wesentlichenEigenschaften des resultierenden Banachraumes nur schwer. Wir gehen daher einen etwas ande-ren Weg.

Definition 3.9 (schwache Ableitung). Seien u, v ∈ L1loc(I) mit

L1loc(I) := w : I → R

∣∣∀K ⊂ I : w∣∣K∈ L1(K).

v ∈ L1loc(I) heisst schwache Ableitung der Ordnung k ∈ N von u genau dann, wenn fur alle

φ ∈ C∞0 (I) gilt ∫Iu(x) · φ(k)(x)dx =

∫I(−1)kv(x)φ(x)dx.

Lemma 3.10. Sei u ∈ L1loc(I). Dann gilt:

i. v, v schwache Ableitungen von u sind, dann gilt v = v f.u.

ii. Sei u ∈ L1loc(I) ∩ C1(I). Dann existiert eine schwache Ableitung von u und ist identisch

mit der (klassischen) Ableitung u′.

Beweis: Ubung.

Im Folgenden benutzen wir die selbe Notation fur schwache und klassische Ableitungen, also u′

bezeichnet auch die schwache Ableitung von u.

Beispiel 3.11 (Existenz nicht trivialer schwacher Ableitungen und Nichtexistenz).

i. Betrachte u ∈ L1loc(I) mit

u(x) =

x fur 0 < x ≤ 1/21/2 fur 1/2 < x < 1.

Behauptung: Es gilt u′(x) = v(x) fur fast alle x ∈ I mit

v(x) =

1 fur 0 < x ≤ 1/20 fur 1/2 < x < 1.

56

3.3 Sobolevraume in 1D

Zum Beweis der Behauptung sei φ ∈ C∞0 (I) gegeben. Dann ist∫ 1

0u(x)φ′(x)dx =

∫ 1/2

0u(x)φ′(x)dx+

∫ 1

1/2

12φ′(x)dx

=[xφ(x)

]1/2

0−∫ 1/2

01φ(x)dx+

12φ(1)− 1

2φ

(12

)=

12φ

(12

)− 0 · φ(0)−

∫ 1/2

0φ(x)dx+ 0− 1

2φ

(12

)= (−1)

∫ 1/2

0v(x)φ(x)dx = (−1)1

∫ 1

0v(x)φ(x)dx

ii. Die Funktion v aus (i) ist nicht schwach differenzierbar. Wir nehmen an, dass eine Funk-tion w ∈ L1

loc(I) mit ∫ 1

0v(x)φ′(x)dx = −

∫ 1

0w(x)φ(x)dx

fur alle φ ∈ C∞0 (I) existiert. Wir haben

−∫ 1

0w(x)φ(x)dx =

∫ 1

0v(x)φ′(x)dx =

∫ 1/2

01 · φ′(x)dx = φ(1/2) (?)

Wahle eine Folge φnn∈N ⊂ C∞0 (I) mit 0 ≤ φn(x) ≤ 1, φn(1/2) = 1, limn→∞

φn(x) = 0 fur

x ∈ I\1/2, supp (φn) ⊆[

14 ,

34

]fur n ∈ N.

Sei δ > 0 beliebig. Wir zeigen, dass ein n0 = n0(δ) existiert, so dass fur alle n > n0∣∣∣∣−∫ 1

0w(x)φn(x)dx

∣∣∣∣ ≤ δ (??)

gilt. Betrachte dazu fur ε ∈ (0, 1/n)∣∣∣∣−∫ 1

0w(x)φn(x)dx

∣∣∣∣ ≤ ∫ 1

0|w(x)| · |φn(x)|dx

=∫ 1/2+ε

1/2−ε|w(x)| · |φn(x)|dx+

∫ 1/2−ε

0|w(x)| · |φn(x)|dx

+∫ 1

1/2+ε|w(x)| · |φn(x)|dx

≤∫ 1/2+ε

1/2−ε|w(x)|dx+ ‖w‖L1([1/4,3/4]) · max

x∈I\[1/2−ε,1/2+ε]φn(x)

= Iε + Jε,n.

Wahle jetzt ε so klein, dass Iε ≤ δ2 gilt. Wahle zu diesem ε die Zahl n0 so groß, dass

Jε,n ≤ δ2 gilt. Also haben wir den Widerspruch zwischen (?) mit φ = φn0 und (??).

Definition 3.12 (Sobolevraum W k,p(I)).Sei p ∈ [1,∞) ∪ ∞ und k ∈ N0. Der Funktionenraum

W k,p(I) := u ∈ L1loc(I)

∣∣u(l) existiert fur l = 1, . . . , k und u(l) ist in Lp(I) fur l = 0, . . . , k

heisst Sobolevraum.

Bemerkung 3.13.

57


i. Der Sobolevraum W 0,p(I) = Lp(I).

ii. Fur p = 2 schreibt man auch Hk(I) := W k,2(I).

Lemma 3.14 (Eigenschaften von Sobolevraumen).Seien u, v ∈W k,p(I) und λ, µ ∈ R. Dann gilt:

i. u(l) ∈W k−l,p(I) fur l ≤ k

ii. ∂l

∂xl

(∂r

∂xr u)

= ∂r

∂xr

(∂l

∂xlu)

fur l, r ∈ N0, l + r ≤ k.

iii. (λu+ µv) ∈W k,p (λu+ µv)′ = λu′ + µv′

iv. J ⊆ I, J offen ⇒ u ∈W k,p(J)

v. Sei φ ∈ C∞0 (I), dann ist φu ∈W k,p(I) und es gilt die Leibnitzregel

(φ · u)(l) =l∑

r=0

(l

r

)φ(r) · u(l−r) fur l ≤ k

Definition 3.15 (Sobolevnorm). Sei u ∈W k,p(I). Dann heisst

‖u‖W k,p(I) :=

(∑k

l=0

∫I |u

(l)|pdx)1/p

fur p = [1,∞)∑kl=0 esssup|u(l)(x)| fur p = ∞

die Sobolevnorm. Hierbei ist

esssup f := infµ ∈ R∣∣ |f > µ| = 0

fur eine messbare, reelwertige Funktion.

Satz 3.16. Fur k ∈ N0 p ∈ [1,∞) ∪ ∞ ist (W k,p(I), ‖.‖W k,p(I)) ein Banachraum. Statt(W k,p(I), ‖.‖W k,p(I)) schreiben wir W k,p(I).

Beweis: Zunachst ist zu zeigen, dass ‖.‖W k,p(I) eine Norm ist.Es gilt offenbar fur u ∈W k,p(I) und λ ∈ R

‖λu‖W k,p(I) = |λ|‖u‖W k,p(I)

und‖u‖W k,p(I) = 0 ⇔ u = 0 fast uberall

Aus der Minkowski-Ungleichung folgt fur u, v ∈W k,p(I), p ∈ [1,∞)

‖u+ v‖W k,p(I) =

(k∑

l=0

∫ 1

0

(u(l) + v(l)

)pdx

)1/p

=

(k∑

l=0

‖u(l) + v(l)‖pLp

)1/p

≤

(k∑

l=0

(‖u(l)‖Lp + ‖v(l)‖Lp)p

)1/p

Mink. fur Vek.≤

(k∑

l=0

‖u(l)‖pLp

)1/p

+

(k∑

l=0

‖v(l)‖pLp

)1/p

def.= ‖u‖W k,p(I) + ‖v‖W k,p(I)

58

3.4 Existenz einer schwachen Losung in (3.1),(3.2)

Die Aussage fur die ∞-Norm folgt direkt mit der Dreiecksungleichung.Im nachsten Teil wird die Vollstandigkeit von W k,p(I) nachgewiesen.Sei vnn∈N eine Cauchyfolge in W k,p(I). Dann ist v(l)

n n∈N fur l = 0, . . . , k eine Cauchyfolgein Lp(I). Da Lp(I) vollstandig ist, existieren Funktionen vl ∈ Lp(I) mit

v(l)n → vl ∈ Lp(I) (l = 0, . . . , k).

Dabei sei v := v0.Behauptung: v ∈W k,p(I), v(l) = vl l = 0, . . . , kZum Beweis sei φ ∈ C∞

0 (I), dann gilt fur l = 0, . . . , k∫Ivφ(l)dx

Leb.Konv.= limn→∞

∫Ivnφ

(l)dx

= limn→∞

∫Iv(l)n φdx(−1)l

= (−1)l

∫Ivlφdx

Damit ist die Behauptung gezeigt und der Satz bewiesen.

Satz 3.17 (Approximation von Sobolev-Raumen).

i. Sei u ∈W k,p(I), dann existiert eine Funktion u ∈ Ck−1(I) mit den Eigenschaften

a) u = u fast uberall in I

b) ∃C > 0, so dass ‖u‖Ck−1(I) ≤ C‖u‖W k,p(I)

ii. Zu u ∈W k,p(I) existiert eine Folge unn∈N ⊂ Ck(I) mit der Eigenschaft

limn→∞

‖u− un‖W k,p(I) = 0

Falls u(0) = 0 gilt mit u aus (i)(a), kann die Folge unn∈N so gewahlt werden, dass un(0) = 0fur n ∈ N.

3.4 Existenz einer schwachen Losung in (3.1),(3.2)

Wir betrachten die Menge

V := u ∈W 1,2(I)∣∣u(0) = 0 (V = u ∈ C1(I)

∣∣u(0) = 0)

Wegen Satz 3.17(i) ist V wohldefiniert (falls man den stetigen Reprasentanten u nimmt, wasimplizit vorausgesetzt sei).V ist ein abgeschlossener Unterraum von W 1,2(I) (und als solcher vollstandig):Sei vnn∈N Cauchyfolge in V und sei v das Grenzelement in W 1,2(I) (existiert nach Satz 3.16).Es gilt

|v(0)| = | vn(0)︸︷︷︸=0

−v(0)| ≤ ‖vn − v‖C0(I)

3.17(i)

≤ C‖vn − v‖W 1,2(I)n→∞−→ 0

Wir kehren zum Variationsproblem (3.5) zuruck, dass wir jetzt nicht in V , sondern in V be-trachten. Finde u ∈ V mit

F [u] =12

∫ 1

0p(u′)2 + qu2dx−

∫ 1

0hudx = inf

v∈VF [v] (3.7)

59


Satz 3.18 (Existenz einer schwachen Losung). Es existiert u ∈ V von (3.7).

Beweis: Nach dem 10. Ubungsblatt wissen wir, dass die Poincare-Ungleichung auch fur Sobo-levfunktionen gilt. Also existiert (wie fur den Fall v ∈ V ) der Ausdruck

infv∈V

F [v] =: d.

Ebenso wie in Lemma 3.7 fur v ∈ V gilt hier: Jede Minimalfolge vnn∈N ⊂ V ist eine Cauchy-Folge in V . Da W 1,2(I) vollstandig ist, und V abgeschlossener Unterraum von W 1,2(I) ist,existiert v ∈ V mit

limn→∞

‖vn − v‖W 1,2(I) = 0,

wobei vnn∈N eine Minimalfolge ist.Es bleibt F [v] = d zu zeigen (Programmpunkt (iv)):

|F [vn]− F [v]| =∣∣∣∣12∫ 1

0((p(vn)′)2 + q(vn)2 − p(v′)2 − q(v2)dx−

∫ 1

0h(vn − v)dx

∣∣∣∣≤ 1

2‖p‖L∞(I)

∫ 1

0

∣∣(v′n)2 − v′2∣∣ dx+

12‖q‖L∞(I)

∫ 1

0|v2

n − v2|dx−∫ 1

0|h||vn − v|dx

Holder≤ 1

2‖p‖L∞(I)

∫ 1

0|v′n − v′| · |v′n + v′|dx+

12‖q‖L∞(I)

∫ 1

0|vn − v| · |vn + v|dx

− ‖h‖L2(I)‖v − vn‖L2(I)

Holder≤ 1

2‖p‖L∞(I)‖v′n + v′‖L2(I) · ‖v′n − v′‖L2(I)

+12‖q‖L∞(I)‖vn + v‖L2(I) · ‖vn − v‖L2(I) − ‖h‖L2(I) · ‖v − vn‖L2(I)

≤ C(p, q, h, ‖v‖W 1,2(I)

)· ‖vn − v‖W 1,2(I)

n→∞−→ 0, da vnn→∞−→ v in W 1,2(I)

Also gilt insgesamt F [vn] → F [v] fur n→∞ und damit ist der Satz bewiesen.

Auf der Basis von Satz 3.18 wird jetzt eine neue Definition der schwachen Losung gegeben.

Definition 3.19 (Schwache Losung von W 1,2(I)). Eine Funktion u ∈ V heisst schwache Losungdes Sturm-Liouville Problems (3.1),(3.2) genau dann, wenn∫ 1

0pu′φ′ + quφdx =

∫ 1

0hφdx

fur alle φ ∈ V .

Satz 3.20 (Existenz und Eindeutigkeit fur schwache Losungen in V ). Es existiert genau eineschwache Losung von (3.1),(3.2) in V .

Beweis:Existenz: Der Beweis ist identisch zum Beweis von Lemma 3.3 nur mit V statt V .Eindeutigkeit: Seien u1, u2 ∈ V zwei schwache Losungen von (3.1),(3.2). Fur w = u1 − u2 giltdann ∫ 1

0pw′φ′ + qwφdx = 0 ∀φ ∈ V

Setze φ = w ∈ V . Dann folgt ∫ 1

0p(w′)2 + qw2dx = 0

60

3.5 Das Ritz-Galerkin-Verfahren als Motivation der FDM

q≥0⇒∫ 1

0p(w′)2dx = 0

p≥p0⇒ p0

∫ 1

0(w′)2dx = 0

Somit ist w′ = u′1 − u′2 = 0 fast uberall in I. Nach Poincare gilt dann ‖w‖L2(I) ≤ ‖w′‖L2(I),woraus w = u1 − u2 = 0 fast uberall folgt.

3.5 Das Ritz-Galerkin-Verfahren als Motivation der FDM

Wir betrachten das Variationsproblem (3.7) jetzt fur einen endlichdimensionalen Teilraum Vh ⊂V mit den Eigenschaften

dim Vhh→0−→∞ (vh ∈ Vh ⇒ vh(0) = 0)

(h wird spater die Gitterweite sein.)Das Ritz-Galerkin Verfahren lasst sich nun wie folgt beschreiben:Finde uh ∈ Vh mit

F [uh] ≤ infvh∈Vh

F [vh] (3.8)

Wie fur die Variationsprobleme (3.5),(3.7) konnen wir nachweisen, dass eine Losung uh von (3.8)(eindeutig bestimmt) eine Losung des schwach formulierten Problems∫ 1

0(pu′hφ

′h + quhφh)dx =

∫ 1

0h · φhdx ∀φh ∈ Vh (3.9)

ist.

Satz 3.21 (Existenz und Eindeutigkeit einer diskreten Losung). Es gibt genau eine Funktionuh ∈ Vh die (3.9) erfullt. uh heisst diskrete Losung.

Beweis: Der Beweis wird wie fur Satz 3.20 gefuhrt (Die Minimalfolge ist eine Cauchy-Folge,aus der Abgeschlossenheit folgt die Konvergenz). Man beachte, dass Vh als endlichdimensionalerTeilraum abgeschlossen ist.

Bevor Beispiele fur Vh diskutiert werden, wird eine abstrakte, aber einfache Fehlerabschatzungbewiesen.

Satz 3.22 (Fehlerabschatzung I). Sei u ∈ V die schwache Losung von (3.1),(3.2). Weiter seiuh ∈ Vh die diskrete Losung. Dann existiert eine von dem Parameter h unabhangige KonstanteC > 0 mit

‖u− uh‖W 1,2(I) ≤ C · infvh∈Vh

‖u− vh‖W 1,2(I).

Das Problem der Fehlerabschatzung ist jetzt auf ein reines Approximationsproblem zuruck-gefuhrt.

Beweis: Sei wieder B : V × V → R durch

B(v, w) =∫ 1

0(pv′w′ + qvw)dx

gegeben. Wir zeigen zunachst

i. ∃C0 > 0, so dass

B(v, v) ≥ C0‖v‖2W 1,2(I) ∀v ∈ V (Koerzitivitat)

61


ii. ∃C1 > 0, so dass

B(v, w) ≤ C1‖v‖W 1,2(I)‖w‖W 1,2(I) ∀v, w ∈ V (Stetigkeit)

iii. B(u− uh, vh) = 0 vh ∈ Vh

(Der Fehler u− uh steht senkrecht auf dem Testraum bezuglich B.)

Abschatzung (i) ist schon im Beweis von Lemma 3.7 fur v ∈ V gezeigt worden. Genauer gilt

B(v, v) ≥ 23p0‖v‖2

W 1,2(I)

Allgemein fur v ∈ V folgt i) dann mit Satz 3.17(ii).Zum Beweis von (ii) betrachte

|B(v, w)| ≤∫ 1

0|pv′w′|+ |qvw|dx

≤ max‖p‖L∞(I), ‖q‖L∞(I)∫ 1

0(|v′w′|+ |vw|)dx

Holder≤ C(p, q) · ‖v‖W 1,2(I) · ‖w‖W 1,2(I)

Zu (iii): Es gilt fur alle vn ∈ Vn ⊂ V , also ist

B(u, vh) =∫ 1

0hvhdx,

da u eine schwache Losung ist.uh ist die diskrete Losung:

B(uh, vh) =∫ 1

0hvhdx

⇒ B(u− uh, vn) = B(u, vh)−B(uh, vh) = 0

Ausgehend von (i) erhalten wir fur beliebiges vh ∈ Vh:

‖u− uh‖2W 1,2(I)

(i)

≤ 1C0B(u− uh, u− uh) =

1C0

(B(u, u− uh)− B(uh, u− uh)︸︷︷︸

=0 wegen (iii),uh∈Vh

)=

1C0B(u, u− uh) =

1C0B(u, u− uh)− 1

C0B(vn, u− uh)︸︷︷︸=0 wegen (iii)

=1C0B(u− vh, u− uh)

(ii)

≤ C1

C0‖u− vh‖W 1,2(I)‖u− uh‖W 1,2(I)

Da vh beliebig war, folgt die Fehlerabschatzung.

62

3.6 Der Teilraum Vh und das Finite-Elemente-Verfahren

3.6 Der Teilraum Vh und das Finite-Elemente-Verfahren

Um ein durchfuhrbares, numerisches Verfahren zu konstruieren, muss noch der Teilraum Vh

bestimmt werden.

Beispiel 3.23 (Globale Ansatzfunktion). Als ersten Vorschlag fur Vh wahlen wir fur h > 0

vh := x, x2, x3, . . . , xn n =⌊

1h

⌋Durch die Polynome x, x2, x3, . . . , xn ist eine Basis von Vh gegeben. Wir stellen die gesuchteFunktion uh ∈ Vh als Linearkombination

uh(x) =n∑

j=1

b− jxj

dar. Dabei sind bj ∈ R die gesuchten Koeffizienten. (Beachte vh ∈ V wegen vh(0) = 0 fur allevh ∈ Vh.)Aus der diskreten schwachen Formulierung (3.9) ergibt sich dann:Finde fur j = 1, . . . , n Koeffizienten bj ∈ R mit∫ 1

0

p n∑j=1

bj · j · xj−1v′h(x) + q

n∑j=1

bjxjvh(x)

dx =∫ 1

0h(x)vh(x)dx ∀vh ∈ Vh

Dies ist aquivalent zu der Beziehung∫ 1

0

p n∑j=1

bjjxj−1kxk−1

+

q n∑j=1

bjxjxk

dx =∫ 1

0h(x)xkdx fur k = 1, . . . , n (?)

Mit p ≡ 1 und q ≡ 0 fuhrt (?) auf∫ 1

0

n∑j=1

(bj · j · xj−1xk−1 · k

)dx =

∫ 1

0h(x)xkdx fur k = 1, . . . , n

Dies fuhrt auf ein lineares Gleichunssystem der Form

Ab = g

fur b = (b1, . . . , bn)T . Dabei sind die Matrix A = (akj)kj ∈ Rn×n und g ∈ Rn gegeben durch

akj =∫ 1

0(k · j · xj−1xk−1)dx = kj

∫ 1

0xj+k−2dx = k · j · 1j+k−1 · 1

j + k − 1

⇒ akj =k · j

j + k − 1(j, k = 1, . . . , n)

und

gk =∫ 1

0h(x)xkdx (k = 1, . . . , n).

A ist (leider) eine vollbesetzte Matrix. Fur die Kondition (cond(A) = ‖A‖ · ‖A−1‖) von A giltfur n = 5:

cond2(A) ≈ 3 · 105

fur n = 10:cond2(A) ≈ 8, 6 · 1012

Mit diesem Ansatz ist das Ritz-Galerkin-Verfahren numerisch unbrauchbar.

63


Beispiel 3.23 zeigt, dass auf ganz I definierte Basisfunktionen auf vollbesetzte Matrizen fuhren.Zur Umgehung dieses Problems sei ein Gitter

Ih = 0 = x0, x1, . . . , xN+1 = 1

gegeben. Dabei seihj := xj+1 − xj und Ij = [xj , xj+1].

Fur k ∈ N sei nun

V kh = φh ∈ C0([0, 1])

∣∣ φh

∣∣Ij∈ Pk ∀j ∈ 0, . . . , N, φn(0) = 0 (3.10)

Dabei ist Pk die Menge der Polynome vom Grad ≤ k.Im nachsten Schritt geben wir eine Basis von V k

h fur den einfachsten Fall k = 1 an.Fur j = 1, . . . , N sei φj ∈ V 1

h durch

φj(x) =

x−xj−1

hj−1x ∈ Ij−1

xj+1−xhj

x ∈ Ij0 sonst

definiert.Weiter sei φN+1 ∈ V 1

h durch

φN+1(x) = x−xN

hNx ∈ IN

0 sonst

definiert.Die Menge φ1, . . . , φN+1 ist linear unabhangig und es gilt

V 1h = spanφ1, . . . , , φN+1 (3.11)

Die Funktionen φ1, . . . , φN+1 heissen P1-Basisfunktionen. Die lineare Finite-Elemente-Methode(FEM) fur (3.1),(3.2) ergibt sich folgendermaßen:Finde Koeffizienten b1, . . . , bN+1 ∈ R, so dass die Funktion

uh(x) =N+1∑j=1

bjφj (3.12)

die Gleichung ∫ 1

0(pu′hφ

′k + quhφk)dx =

∫ 1

0(hφk)dx (3.13)

fur k = 1, . . . , N + 1 erfullt.

Bemerkung 3.24.

i. Die Approximationsfunktion uh aus (3.12) heisst auch Ansatzfunktion. Im Gegensatz zurFDM gilt uh : I → R (und nicht uh : Ih → R).

ii. Weiterhin ist das Gitter Ih nicht als aquidistant vorausgesetzt. Die FEM lasst sich aufbeliebigen Gittern anwenden.

iii. Im Gegensatz zu den Basisfunktionen aus Beispiel 3.23 ist der Trager der Basisfunktionφ1, . . . , φN+1 ”klein“ (genau zwei Zellen Ij , j ∈ 0, . . . , N)

iv. Im Kontext der FEM wird das Intervall Ij auch Finites Element genannt.

64

3.7 Konvergenz der FEM in 1D

Aus (3.13) ergibt sich das Gleichungssystem

Ab = g fur b = (b1, . . . , bN+1)T ∈ RN+1.

Dabei ist im Falle p = const und q ≡ 0 die Matrix A = (ajk) ∈ R(N+1)×(N+1) durch

akj =∫ 1

0φ′jφ

′kdx (j, k = 1, . . . , N + 1)

gegeben. Fur g = (g1, . . . , gN+1)T ∈ RN+1 ist

gk =∫ 1

0hφkdx.

Fur die Koeffizienten akj ergibt sich wegen

supp (φj) ∩ supp (φk) = ∅ fur |j − k| ≥ 2

die Beziehungakj = 0 fur |j − k| ≥ 2.

Außerdem gilt akj = ajk fur j, k ∈ 1, . . . , N + 1. Die Matrix A ist also eine Bandmatrix mit 3Bandern und symmetrisch.


In Satz 3.22 haben wir die abstrakte Fehlerabschatzung

‖u− uh‖W 1,2(I) ≤ C infvhVh

‖u− vh‖W 1,2(I). (3.14)

Konnen wir die rechte Seite in (3.14) fur Vh = V 1h aus (3.10) weiter gegen |h| = maxj=1,...,Nhj

abschatzen?

Definition 3.25 (Interpolierende). Die Abbildung Ih : V → V heisst die lineare Interpolati-onsabbildung auf dem Gitter Ih genau dann, wenn gilt

i. v ∈ V ⇒ (Ihv)(xj) = v(xj) ∀j ∈ 0, . . . , N + 1

ii. Ihv ∈ V 1h

Wir nutzen jetzt die Ungleichnung

infvh∈V 1

h

‖u− vh‖W 1,2(I) ≤ ‖u− Ihu‖W 1,2(I)

Lemma 3.26 (Eine Art Poincare-Ungleichung). Sei v ∈W 1,2(I) und∫ 10 v(x)dx = 0 dann ist

‖v‖2L2(I) ≤ ‖v′‖2

L2(I)

Beweis: Zunachst sei v ∈W 1,2(I)∩C1(I). Nach Voraussetzung existiert ein ξ ∈ I mit v(ξ) = 0.Also ist fur x ∈ I

v(x) = v(x)− v(ξ) =∫ x

ξv′(y)dy

65


⇒ v(x)2 = (x− ξ)2(

1x− ξ

∫ x

ξv′(y)dy

)2

Jensen≤ (x− ξ)2

1|x− ξ|

∫ x

ξ|v′(y)|2dy

= |x− ξ|∫ x

ξ|v′(y)|2dy

⇒∫ 1

0|v(x)|2 dx ≤

∫ 1

0

(|x− ξ|︸︷︷︸≤1

∫ x

ξ|v′(y)|2dy

)dx ≤

∫ 1

0

(∫ 1

0|v′(y)|2dy

)dx =

∫ 1

0|v′(y)|2dy

Anwendung von Satz 3.17(ii) ergibt die Behauptung fur v ∈W 1,2(I).

Lemma 3.27 (Interpolation auf dem Referenzelement). Sei Ih = x0, x1 = 0, 1 und diezugehorige Interpolationsabbildung I = I1 gegeben.

i. Sei v ∈ V . Dann existiert eine Konstante C1 > 0 mit

‖v − Iv‖L2(I) ≤ C1‖v′‖L2(I)

‖(Iv)′‖L2(I) ≤ C1‖v′‖L2(I)

ii. Sei v ∈W 2,2(I) ∩ V . Dann existiert eine Konstante C2 > 0 mit

‖v − Iv‖L2(I) ≤ C2‖v′′‖L2(I)

‖(v − Iv)′‖L2(I) ≤ C2‖v′′‖L2(I)

Beweis:

i. Nach der Poincare-Ungleichung (in der Form von Lemma 3.5) gilt wegenv(0) = 0 = (Iv)(0) = 0:

‖v − Iv‖L2(I) ≤12‖v′ − (Iv)′‖L2 =

12

∥∥∥∥v′ − v(1)− v(0)1− 0

∥∥∥∥L2

=12‖v′ − v(1)‖L2

‖v − Iv‖2L2(I) ≤

14

∫ 1

0(v′(x)− v(1))2dx

=14

∫ 1

0(v′(x))2 − 2v(1)v′(x) + (v(1))2dx

≤ 14‖v′‖2

L2 −14(v(1))2

=14‖v′‖2

L2 −14‖(Iv)′‖2

L2

⇒ ‖v − Iv‖2L2(I) +

14‖(Iv)′‖2

L2(I) ≤14‖v′‖2

L2(I)

und somit ist (i) bewiesen.

ii. Wir haben

‖v − Iv‖L2

Lemma 3.5≤ 1

2‖v′ − (Iv)′‖L2

Lemma 3.26≤ 1

2‖(v − Iv)′′‖L2

v∈W 2,2(I)

≤ 12‖v′′‖L2

66


Die Abschatzungen aus Lemma 3.27 werden jetzt genutzt, um eine Abschatzung ‖u−Ihu‖W 1,2(I)

herzuleiten. Dazu ist

‖u− Ihu‖2L2(I) =

N∑j=0

‖u− Ihu‖2L2(Ij)

mit Ij = [xj , xj+1), hj = xj+1 − xj . Fur j = 0, . . . , N beliebig, aber fest, betrachten wir


=∫ xj+1

xj

(u(x)− (Ihu)(x))2dx

x=hj x+xj= hj

∫ 1

0(u(hj x+ xj)− (Ihu)(hj x+ xj))2dx

Wegen

(Ihu)(x) = u(xj) +1hj

(x− xj)(u(xj+1)− u(xj))

folgt


= hj

∫ 1

0

(u(hj x+ xj)−

(u(xj) +

1hj

(hj x+ xj − xj)(u(xj+1)− u(xj))))2

dx

= hj

∫ 1

0(u(hj x+ xj)︸︷︷︸

=:u(x)

− (u(xj) + x(u(xj+1)− u(xj)))︸︷︷︸I1u(x)

)2dx

= hj

∫ 1

0(u− I1u)2dx

Lemma 3.27(i)

≤ hjC1‖u′′‖2L2(I) = hjC1

∫ 1

0|u′′(x)|2dx

= hjC1

∫ 1

0|h2

ju′′(hj x+ xj)|2dx

= h5jC1

∫ 1

0|u′′(hj x+ xj)|2dx

= h5jC1

∫ xj+1

xj

|u′′(x)|2 1hjdx

= C1h4j‖u′′‖2

L2(Ij)

Insgesamt haben wir mit h := maxj=0,...,Nhj

‖u− Ihu‖2L2(I) ≤ C1

N∑j=0

h4j‖u′′‖2

L2(Ij)≤ C1h

4‖u′′‖2L2(I) (3.15)

Weiter gilt mit Aussage (ii) aus Lemma 3.27 und ahnlicher Rechnung

‖(u− Iju)′‖2L2(I) ≤ C2|h|2‖u′′‖2

L2(I) (3.16)

Aus (3.15) und (3.16) ergibt sich das Konvergenztheorem

Satz 3.28 (Konvergenz der FEM). Sei u ∈ V die schwache Losung des Sturm-Liouville Pro-blems (3.1),(3.2) und uh ∈ V 1

h die lineare FE-Approximation aus (3.12). Falls u ∈W 2,2(I) gilt,folgt

(‖u− uh‖L∞(I) ≤)‖u− uh‖W 1,2(I) ≤ C|h|‖u′′‖L2(I)(‖u− uh‖L2(I) ≤ C|h|2‖u′′‖L2(I)

)

67


3.8 Abschließende Bemerkungen zur FEM

Bemerkung 3.29 (Schwache Formulierung und FEM). Sei a ∈ L∞(I × R,R) und f, h ∈L∞(I ×R,R)× L2(I ×R,R). Betrachte die (nichtlineare) Differentialgleichung

−(a(x, u)u′)′f(x, u)′ = h(x, u) in I (?)

mit geeigneten Randbedingungen. Eine schwache Formulierung fur (?) ist∫Ia(x, u)u′φ′ − f(x, u)φ′)dx =

∫Ih(x, u)φdx ∀φ ∈ V (??)

Dabei sei V ein geeigneter Funktionenraum, etwa V = W 1,2(I). Man erhalt also eine FE-Approximation uh fur u, falls man uh und Testfunktionen φh aus einem endlichdimensionalenTeilraum Vh ⊂ V wahlt.

Bemerkung 3.30 (Dirichlet-Problem). Bisher haben wir die speziellen Randbedingungen (3.2)betrachtet. Um diese zu berucksichtigen, war der Raum

V := v ∈W 1,2(I)∣∣v(0) = 0

definiert worden.Im allgemeinen Fall schreibt man die Randbedingungen vollstandig in den Ansatzraum, etwa imFalle der Dirichlet-Bedingungen u(0) = u(1) = 0

V = v ∈W 1,2(I)∣∣v(0) = v(1) = 0.

Bemerkung 3.31 (Die Stetigkeitsmatrix A). Die Matrix A aus der FEM heisst Stetigkeitsma-trix. Es gilt fur vh =

∑N+1j=1 bjφj , bj ∈ R, vn 6= 0 im Spezialfall q ≡ 0, p ≡ 1

0 < B(vh, vh) =12

∫ 1

0((vh)′)2dx =

12

∫ 1

0

∑j

bjφ′h

∑k

bkφ′kdx

=12

N+1∑j,k=0

bjbk

∫ 1

0φ′jφ

′kdx︸︷︷︸

akj

=12bTAb

Also ist A (symmetrisch und) positiv definit, da b ∈ RN+1 beliebig.

kontinuierlich diskret

Variationsproblem u : F [u] ≤ F [v]∀v ∈ V uh : F [uh] ≤ F [v]∀v ∈ V 1

h

b ∈ RN+1 : f(b) ≤ f(c)∀c ∈ RN+1

⇓ ⇓

schwache Losung u : B(u, v) =∫hvdx∀v ∈ V

uh : B(uh, vh) =∫ 10 hvhdx∀vh ∈ V

b ∈ RN+1 : Ab = g

68

4 Iterative Losungsverfahren fur lineareGleichungssysteme

Die FDM und die FEM zur Losung z.B. des Sturm-Liouville Problems fuhren auf lineare Glei-chungssysteme vom Typ

Ab = g, (4.1)

wobei A ∈ Rm×m, g ∈ Rm fur m ∈ N sei. Dabei sei A charakterisiert durch folgende Eigenschaf-ten:

• A = AT , A positiv definit

• m 1

• A ist dunn besetzt, z.B. eine Bandmatrix.

• cond2(A) ≈ h−p p ≥ 1

Bemerkung 4.1 (direkte Verfahren). Direkte Verfahren fur derartige Matrizen (Cholesky)benotigen O(m3/2) Elementaroperationen. Die Zerlegungsmatrizen L ∈ Rm×m aus A = LTLsind (unterhalb der Hauptdiagonale) i.A. voll besetzt. Konnen die Anzahl der Elementaropera-tionen und der Speicheraufwand reduziert werden?Bei der FEM fuhrt die Diskretisierung auf das Variationsproblem:Finde b ∈ Rm : f(b) := 1

2bTAb−gb ≤ f(c) ∀c ∈ Rm. Dazu aquivalent ist die Losung von (4.1).

Definition 4.2 (Iterationsverfahren). Sei der Startwert b0 ∈ Rm gegeben. Weiter seien furk ∈ N0, αk ∈ R, dk ∈ Rm. Die Vorschrift

bk+1 = bk − αkdk (k ∈ N0)

fur die Folge bkk∈N0 heisst Iterationsverfahren zur Losung von (4.1). Fur k ∈ N0 heisst dk

Abstiegsrichtung und αk Abstiegsparameter.

Beispiel 4.3 (2× 2-Problem). Sei A ∈ R2×2 mit

A =(λ1 00 λ2

)(λ1, λ2 > 0)

Wir betrachten die Hohenlinien von f(b) = 12b

TAb− 0 · b (g = 0).(Also ist b = 0 die Losung von Ab = g)Wahle d0 so, dass gilt

∂

∂d0f(b0) = d0∇bf(b0) ≤ 0

Wahle α0 so, giltf(b0 + α0d

0) ≤ f(b0) b1 = b0 + α0d0

Wahle d1 so, dass gilt∂

∂d1f(b1) ≤ 0

Wahle α1 so, dass giltf(b1 + α1d

1) ≤ f(b1) b2 = b1 + α1d1

69

4 Iterative Losungsverfahren fur lineare Gleichungssysteme

Definition 4.4 (Gradientenverfahren). Sei A positiv definit und symmetrisch. Ein Iterations-verfahren heisst Gradientenverfahren fur (4.1) genau dann, wenn gilt

i. dh = −∇bf(bk) = −Abk + g (k ∈ N0) (steilster Abstieg)

ii. αk = α (k ∈ N0)

Satz 4.5 (Konvergenz). Fur das Gradientenverfahren gilt

‖bk − b‖2 ≤ γk‖b0 − b‖2 mit γ < 1.

Ein wesentlich besseres Iterationsverfahren ergibt sich aus

Definition 4.6 (cg-Verfahren, conjugate gradient). Ein Iterationsverfahren heisst cg-Verfahrenfur (4.1) genau dann, wenn gilt

i. dk = −rk + βk−1dk−1

βk−1 =rTk Ad

k−1

(dk−1)TAdk−1, rk = Abk − g

ii.

αk :=rTk d

k

(dk−1)TAdk−1

Bemerkung 4.7. Betrachte das cg-Verfahren fur (4.1). Fur α mit

f(bk + αdk) ≤ minα∈R

f(bk + αdk)

gilt∂

∂αf(bk + αdk) = 0

⇔ dk∇bf(bk + αdk) = 0

⇔ dk(A(bk + αdk)− g) = 0

α(dk)TAdk + dk(Abk − g) = 0

α = − dkrk

(dk)TAdk

Satz 4.8 (Konvergenzeigenschaften). Fur ein Iterationsverfahren gelte fur alle k ≥ k0, k0 ∈ N

‖b− bk‖2 ≤ ε‖b− b0‖2

Dabei sei ε > 0 gegeben:

i. Fur das Gradientenverfahren gilt

k0 ≥ C1cond2(A) ln(

2ε

)ii. Fur das cg-Verfahren gilt

k0 ≥ C2

√cond2(A) ln

(2ε

).

70

Numerik II - math.uni-bielefeld.demstephan/files/scripts/NumerikII.pdf · 1 Anfangswertprobleme...

Documents

Transcript of Numerik II - math.uni-bielefeld.demstephan/files/scripts/NumerikII.pdf · 1 Anfangswertprobleme...