Numerik Gewoehnlicher DGL
Numerik Gewoehnlicher DGL
Numerik Gewoehnlicher DGL
ohnlicher Differentialgleichungen
Bernd Simeon
4. Mehrschrittverfahren
6. Stochastische Differentialgleichungen
Literatur:
Deuflhard/Bornemann: Numerische Mathematik II,
3. Auflage; Walter de Gruyter, Berlin 2002
Hairer/Norsett/Wanner: Solving Ordinary Differential Equations I,
2. Auflage, Springer, Heidelberg 1993
Stoer/Bulirsch: Numerische Mathematik II, 4. Auflage, Springer, Berlin 2000
Strehmel/Weiner: Numerik gewohnlicher Differentialgleichungen,
Teubner, Stuttgart 1995
Vorwort
Inhalt dieser Vorlesung ist die numerische Losung eines Systems gewohnlicher
Dierentialgleichungen
y (x) = f (x, y(x)) (1.1)
bzw. ausgeschrieben in Komponenten
y1 (x) = f1 (x, y1 (x), . . . , yn (x)),
y2 (x) = f2 (x, y1 (x), . . . , yn (x)),
.. ..
. .
yn (x) = fn (x, y1 (x), . . . , yn (x)).
Man unterscheidet dabei zwischen Anfangswertproblemen und Randwert-
problemen. Bei ersteren fordert man
y(x0 ) = y0
an einer Stelle x0 , bei letzteren dagegen
r(y(a), y(b)) = 0
mit einer vorgegebenen Funktion r und einem Intervall [a, b].
Bevor wir die Problemklasse naher analysieren, wollen wir in diesem Ka-
pitel einige Beispiele aus Naturwissenschaft und Technik kennenlernen. Als
1
erstes wird die chemische Reaktionskinetik behandelt, danach folgt die Me-
chanik starrer Korper, und abschlieend werden elektrische Schaltkreise ein-
gefu
hrt. In allen drei Beispielen verwendet man mathematische Modelle, um
einen realen Vorgang zu beschreiben. Ohne Vereinfachungen und Model-
lannahmen ist dies im Allgemeinen nicht moglich, doch wird man von ei-
nem vernu nftigen Modell verlangen, dass es die interessierenden Phanomene
bzw. Eekte hinreichend gut widerspiegelt.
Betrachten wir dazu zwei Gase A und B unter den Annahmen konstan-
ter Druck, konstantes Volumen und konstante Temperatur. Im Fall einer
monomolekularen Reaktion
k
A B
mit Reaktionsgeschwindigkeit k gilt fu
r die zeitliche Anderung der Konzen-
trationen [A] und [B]
d = k[A], = k[A]
[A] = [A] [B]
dt
Die Exponentialfunktion als Losung beschreibt also den zeitlichen Verlauf
+ [B]
der Reaktion. Wegen [A] = 0 folgt [A] + [B] = konstant, was man als
Massenerhaltung bezeichnet.
2
mit Substanzen A, B, C, D. In diesem Fall ergeben sich die Dierentialglei-
chungen fu
r die Konzentrationen zu
= k[A][B],
[A] = k[A][B],
[B]
= k[A][B], = k[A][B]. (1.2)
[C] [D]
Im allgemeinen Fall liegen insgesamt n Substanzen sowie m Reaktionen
zwischen ihnen vor,
n
kj
n
ij Ai ij Ai , j = 1, . . . , m.
i=1 i=1
Dabei bezeichnen die ganzen positiven Zahlen ij , ij die Anteile der Sub-
stanzen (stochiometrische Konstanten) und die reellen Parameter kj die
unterschiedlichen Reaktionsgeschwindigkeiten. Als Dierentialgleichung fur
die Konzentration der Substanz Ai hat man dann
m
n
[A i ] = (ij ij )kj [Al ]lj , i = 1, . . . , n. (1.3)
j=1 l=1
A = k1 A, X = k1 A,
B = k2 BX, X = k2 BX,
Y = k2 BX, D = k2 BX,
X = k4 X, E = k4 X.
Die dritte Reaktion ist autokatalytisch trimolekular und liefert nach der allgemeinen Regel (1.3)
die Beitrage
X = k3 X 2 Y, Y = k3 X 2 Y.
3
Mit Sortieren nach Unbekannten und Aufsummieren der rechten Seiten folgen die gekoppelten
Dierentialgleichungen
A = k1 A,
B = k2 BX,
D = k2 BX, (1.4)
E = k4 X,
X = k1 A k2 BX + k3 X 2 Y k4 X,
Y = k2 BX k3 X 2 Y.
Als Vorbereitung fur die Diskussion als dynamisches System in Kapitel 2 werden die Gleichungen
(1.4) vereinfacht: Die Dierentialgleichungen fur D und E werden weggelassen (warum?), A und
B werden als konstant angenommen. Wir setzen ki = 1, schreiben x statt t f ur die unabhangige
Variable sowie y1 (x) := X(x), y2 (x) := Y (x).
(Notationsvereinbarung: y(t) mit Zeit t : dt d
y = y, d
y(x) mit unabh. Var. x : dx y = y)
y1 = A + y12 y2 (B + 1)y1 ,
(1.5)
y2 = By1 y12 y2 ,
Die Bewegung solcher Systeme wurde schon von Euler und Lagrange unter-
sucht, und auf diese beiden beru
hmten Mathematiker geht auch der Kalku l
zuru
ck, mit dem man die zugehorigen Dierentialgleichungen aufstellt. Sei
4
T (q, q)
die kinetische Energie des Systems in den Koordinaten q und Ge-
schwindigkeiten q, dann genu gt die Bewegung den Lagrange-Gleichungen
2. Art ( )
d
T T = Q(q, q, t) (1.6)
dt q q
mit eingepragten und angewandten Kraften Q (Gravitation, Feder-Dampfer,
etc). Im Spezialfall eines konservativen Systems folgt Q aus der potentiellen
Energie U (q) uber Q = q
U . (In diesem Fall konnte man auch das dazu
d
aquivalente Hamiltonprinzip c (T U )dt stationar heranziehen).
Das Doppelpendel ist ein System aus zwei starren Korpern mit konstantem Querschnitt. Die
Modellierung geht in folgenden Schritten vor:
5
Inertialsystem
Y
y2
Schwerkraft
2
x2
y1 x1
1 X
mit Gravitationskonstante .
Schritt 2: Einf
uhrung von Koordinaten zur Lokalisierung der materiellen Punkte:
( ) ( ) ( )
X1 x1 cos 1 sin 1
= A(1 ) , A(1 ) := Drehmatrix,
Y1 y1 sin 1 cos 1
und ( ) ( ) ( )
X2 cos 1 x2
= l1 + A(2 ) , l1 Lange Pendel 1.
Y2 sin 1 y2
und
U = 1/2 (l1 (m1 + 2m2 ) sin 1 + l2 m2 sin 2 )
mit Massen mi := i dVi und Tragheitsmomenten Ji := i (x2i + yi2 ) dVi , i = 1, 2.
6
mit der Massenmatrix
( )
J1 + m2 l12 1/2 l1 l2 m2 cos(1 2 )
M (1 , 2 ) = .
1/2 l1 l2 m2 cos(1 2 ) J2
= 2 ?
enthalten. Wie kommt man von dieser Gleichung auf die Schwingungsgleichung
7
wobei die Unbekannten y aus Spannungen und Stromen zusammengesetzt
sind. Die Eintrage in der Matrix C bestehen aus Kapazitaten. Falls keine
solchen auftreten (und keine Induktivitaten), ist C 0, und die Gleichungen
reduzieren sich auf den stationaren Fall.
In vielen Anwendungen ist die Matrix C singular, und es liegt somit keine
gewohnliche Dierentialgleichung vor, sondern ein sogenanntes dierentiell-
algebraisches System. Einzelheiten dazu in Kapitel 6. Das folgende Beispiel
dagegen enthalt neben einer Kapazitat noch eine Induktivitat und wird
durch eine Dierentialgleichung 2. Ordnung beschrieben.
Beispiel 1.3: Elektrischer Schwingkreis
Wir betrachten als einfaches Beispiel einen elektrischen Schwingkreis, der aus einer Spannungs-
quelle U (t), einer Kapazitat C, einer Induktivitat L und einem Widerstand R besteht, siehe
Abb. 2.
U(t)
C
F
ur die Spannungsabfalle UR , UC sowie UL gelten die Gesetze
UR = RI, UC = Q/C, UL = LI
wobei I den ieenden Strom und Q die Ladung der Kapazitat bezeichnen. Die Kirchhoschen
Regeln (hier die Maschenregel) ergeben
U = UR + UC + UL .
8
Kapitel 2
Exkurs in die Theorie
gewohnlicher Differentialgleichungen
Dieses Kapitel wiederholt zunachst die zentralen Ergebnisse aus der Ana-
lysis und beschaftigt sich dann vor allem mit der Kondition des Anfangs-
wertproblems. Neben dem Einuss der Anfangswerte spielt hier auch die
Abhangigkeit von Parametern eine wichtige Rolle. Schlielich folgt eine kur-
ze Einfu
hrung in die Dynamischen Systeme, um einige wesentliche Begrie
zur Charakterisierung von Losungen bereitzustellen.
Als Einstieg und Wiederholung beginnen wir mit einem Blick auf lineare
Dierentialgleichungen. Fu
r das skalare lineare Anfangswertproblem (AWP)
y = a y + b, y(x0 ) = y0 ,
9
Dieses Vorgehen verallgemeinert man schnell auf Systeme linearer Dieren-
tialgleichungen
y = A y + b, A Rnn , b Rn .
Wir nehmen an, dass A diagonalisierbar ist (d.h. n l.u. EV besitzt). Mit
z(x) = exp(A(x x0 ))z0 hat man eine Losung des homogenen Systems
z = Az. Dabei stellen die Spalten der Matrix
(A )k
exp(A ) = Matrizenexponentielle (2.1)
k!
k=0
Es ist vorteilhaft, die Gleichung (2.2) auf die neue Zeit x = t/ zu transformieren (die Periode
hangt dann nicht mehr von ab). Mit u(x) = z(x) = z(t(x)), u = z, folgt
1
u + (u2 1)u + u = 0.
2
Von besonderem Interesse ist der Fall 1, bei dem wir := 1/2 1 setzen (ein sogenanntes
singular gestortes System, siehe Kapitel 5 und 6). In einem letzten Schritt f
uhren wir noch die
Koordinaten nach Lienhard ein,
10
2.5
1.5
0.5
y2
0.5
1.5
2.5
1.5 1 0.5 0 0.5 1 1.5
y1
Bevor wir zum Existenz- und Eindeutigkeitssatz kommen, seien noch drei
oft hilfreiche Bemerkungen erwahnt:
y = 2 y y = v, v = 2 y.
11
Wir betrachten nun das allgemeine Anfangswertproblem
D := {(x, y) : a x b, y Rn },
falls
f (x, y) f (x, z) L y z (x, y), (x, z) D. (2.5)
Einfaches Kriterium: Falls f sowie f /y stetig auf dem Denitionsgebiet
D sind und f /y beschrankt ist, f /y L, folgt daraus die Lipschitz-
stetigkeit.
Beweisskizze: F
ur je zwei Punkte (x, y), (x, z) in D gilt (MWS der Dierentialrechnung)
Satz 2.1 Falls die rechte Seite f auf dem Definitionsgebiet D stetig ist
und der Lipschitzbedingung (2.5) genu gt, besitzt das Anfangswertproblem
(2.4) eine eindeutige, stetig differenzierbare Losung y auf dem ganzen In-
tervall [a, b].
Eine hohere Glattheit der Losung hangt von der Glattheit der rechten Seite
f bzw. ihrer Ableitungen p f /y p ab. Falls alle partiellen Ableitungen bis
zur Ordnung p stetig sind, ist auch y insgesamt p mal stetig dierenzierbar.
Lasst man dagegen die Lipschitzstetigkeit als Voraussetzung weg und geht
nur von einer stetigen rechten Seite f aus, kann man nur noch die Existenz
einer Losung nachweisen, die Eindeutigkeit geht verloren (Existenzsatz von
Peano).
12
Hinreichend fu r Existenz und Eindeutigkeit ist im u brigen bereits eine
abgeschwachte Form der Lipschitz-Bedingung (2.5), die lokale Lipschitz-
Stetigkeit. Sie ist gegeben, falls es zu jeder kompakten Teilmenge K D
eine Lipschitzkonstante LK gibt mit
f (x, y) f (x, z) LK y z (x, y), (x, z) K. (2.6)
Allerdings kann mit dieser schwacheren Voraussetzung nicht die Existenz
auf dem ganzen Intervall [a, b] garantiert werden. Die folgenden zwei Bei-
spiele beleuchten die Situation.
Beispiel 2.2: y = y 1/2 , y(0) = 1.
Die rechte Seite besitzt in y = 0 eine Singularitat, und die Voraussetzungen (2.5) wie auch (2.6)
werden in einer Umgebung der Singularitat verletzt. Die exakte Losung lautet
y(x) = (1 3x/2)2/3
und ist auf den ersten Blick nur auf x ] , 2/3[ deniert. F
ur x 2/3 lauft sie in die
Singularitat (sie kollabiert). Ist y auch rechts von der Singularitat deniert?
Als ersten Hinweis, dass sich ein Randwertproblem (RWP) ganzlich anders
verhalten kann als ein AWP (mehr in Kapitel 7), behandeln wir noch das
Beispiel 2.4: y + y = 0 mit RB
( )
a) y(0) = 0, y 2 = 1, b) y(0) = 0, y() = 0, c) y(0) = 0, y() = 1.
13
2.2 Einflu von St
orungen
und damit
x
y(x) z(x) y(x0 ) z(x0 ) + L y(t) z(t) dt
| {z } | {z }
=:m(x) x0 =m(t)
bzw.
x
m(x) y(x0 ) z(x0 ) + L m(t) dt. ()
x0
m(x) ist stetig und
x
Lx
q(x) := e m(t) dt
x0
stetig dierenzierbar. Also gilt
( )
m(x) = eLx q(x) = LeLx q(x) + eLx q (x). ()
Eingesetzt in ():
LeLx q(x) + eLx q (x) y0 z0 + LeLx q(x)
= q (x) y0 z0 eLx
x
( )
= q(x) = q (t) dt y0 z0 eLx0 eLx /L
x0
14
Abschatzungen eingesetzt in ():
m(x) y0 z0 (eL(xx0 ) 1 + 1)
z L(x-x 0)
y
0
0
} <e |y0 -z 0 |
x0 x
Im Beweis von Satz 2.2 ist eine vereinfachte Fassung des Lemmas von
Gronwall versteckt, das eine wichtige Rolle bei verschiedenen Abschatzungen
spielt. In der vollstandigen Fassung lautet es:
Mit dem Lemma von Gronwall kann man die Abschatzung u ber die Aus-
wirkung von Storungen in den Anfangswerten erweitern. Sei y eine exakte
Losung und z eine approximative Losung, die einen Defekt (x) aufweist,
15
Mit m(x) := y(x) z(x) und := y(x0 ) z(x0 ) sind die obigen Vor-
aussetzungen erfu
llt, denn
x x
y(x) z(x) = y(x0 ) z(x0 ) + (t) dt + (f (t, y) f (t, z)) dt.
x0 x0
Ableitung bezu
glich einem Parameter
16
Solution of van der Pol Equation, = 1 Solution of van der Pol Equation, = 3
2.5 2.5
2 2
1.5 1.5
1 1
0.5 0.5
solution y1
solution y1
0 0
0.5 0.5
1 1
1.5 1.5
2 2
2.5 2.5
0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20
time t time t
wobei p fur einen (oder mehrere) Parameter steht. Die partiellen Ableitun-
gen f /y sowie f /p seien stetig in einer Umgebung der Losung y. Als
Sensitivitat (bzw. Sensitivitatsmatrix ) bezeichnet man die Ableitung der
Losung nach dem Parameter,
y(x; p) y(x; p + p) y(x; p)
(x) := = lim .
p p0 p
Die Sensitivitat kann als Losung einer speziellen Dierentialgleichung cha-
rakterisiert werden. Dazu betrachten wir neben der Losung y zum Parame-
ter p eine weitere Losung z zum Parameter q = p + p,
17
falls wir A = max f /p setzen. Damit kann das Lemma von Gronwall
angewandt werden (beachte = z(x0 ) y(x0 ) = 0), und wir erhalten
(vergl. die Abschatzung (2.7))
A ( )
z(x) y(x) |q p| e L(xx0 )
1 .
L
Die Dierenz der beiden Losungen kann also u ber die Dierenz der Para-
meter p = q p beliebig klein gemacht werden. Im Grenzu bergang p 0
folgt so aus (2.9) die Dierentialgleichung oder Sensitivitatengleichung
f f
(x) = (x, y(x); p) (x) + (x, y(x); p). (2.10)
y p
Da (2.10) bei gegebenem y ein System linearer Dierentialgleichungen fu
r
die Sensitivitat darstellt, kann man die entsprechende Losungstheorie
einsetzen, um eine expliziten Ausdruck fu
r zu gewinnen.
Exkurs: Losung von y (x) = A(x)y(x) + b(x), A(x) Rn Rn , b(x) Rn stetig.
Ein Fundamentalsystem linear unabhangiger Losungen 1 , . . . , n der homogenen Gleichung
z = Az fassen wir zusammen als
und denieren
P (x, x0 ) := W (x)W 1 (x0 ) Propagationsmatrix.
Im homogenen Fall hat man dann f
ur die Losung zum Anfangswert z0 die Darstellung
(Folgt u
ber Variation der Konstanten.)
Anwendung von ( 2.12) auf die Gleichung (2.10) fu r die Sensitivitat ergibt
(beachte 0 = (x0 ) = 0)
x
f
(x) = P (x, s) (s, y(s); p) ds (2.13)
x0 p
18
mit der Propagationsmatrix P des linearen Systems z = f /y z. Der
Wie wirken sich nun Storungen im Parameter auf die Losung aus? Wegen
y
y(x; p + p) = y(x; p) + (x; p)p + o(|p|)
p
r y(x) = y(x; p + p) y(x; p) die Beziehung
gilt in linearer Theorie fu
.
y(x) = (x) p, (2.14)
und (x) ist die zugehorige (punktweise) Konditionszahl.
Abh
angigkeit der L
osung von den Anfangswerten
Die Anfangswerte kann man auch als Parameter interpretieren. Dazu be-
trachten wir die Dierentialgleichung
z (x) = F (x, z(x); y0 ) := f (x, z(x) + y0 ), z(x0 ) = 0.
Oensichtlich hangt z vom Parametervektor y0 ab, und y(x) := z(x) + y0
lost das AWP
y (x) = f (x, y(x)), y(x0 ) = y0 .
19
Damit kann die obige Argumentation wiederholt werden. Deniere die Sen-
sitivitatsmatrizen
z(x; y0 ) y(x; y0 )
(x) := , (x) := = (x) + Inn .
y0 y0
Die Dierentialgleichung fu
r kann von (2.10) u
bernommen werden und
ergibt sich zu
F F f
(x) = (x, z(x); y0 )(x)+ (x, z(x); y0 ) = (x, z(x)+y0 )((x)+I),
z y0 y
und wegen (x) = (x) sowie (x0 ) = 0nn folgt:
Die Aussagen von Satz 2.4 und Satz 2.5 komplettieren damit die in Satz 2.2
und (2.7) gegebenen quantitativen Abschatzungen.
20
2.3 Differentialgleichungen als Dynamische Systeme
Was kann man u ber Losungen von (2.16) qualitativ aussagen? Als Aus-
gangspunkt der Analyse betrachtet man kritische Punkte ys mit
f (ys ) = 0, d.h. ys (x) 0 .
In kritischen Punkten verschwindet die rechte Seite f , die Losung ist stati-
onar (Abb. 6). Im nachsten Schritt wird um ys linearisiert. Setze
y = ys + y = y = ys + y = y
und weiter
y = f (y) = f (ys + y) = f (ys ) +Df (ys )y + o(y)
|{z} | {z }
=y =0
y y
1
ys
ys
y
x 2
21
Also beschreibt die lineare Dierentialgleichung
.
y = Df (ys ) y
Nehmen wir an, die Eigenwerte 1,2 von Df (ys ) seien verschieden und un-
gleich Null. Die folgenden Falle/Phanomene konnen dann auftreten.
y
y(x)
y
s
y
x
y
1
y y
y
y
s
~
y
y~
x
22
c) EW reell, 0 < 1 < 2 : Dann ist ys instabiler Knoten, alle Losungen y
sind instabil.
y
y
y
y
s
y y
y
s
y y
y
s
23
f ) EW 1,2 = i, > 0: Dann ist ys instabiler Strudelpunkt.
y y
y
s
Wichtiger ist, dass neben den Fallen a) - f) noch weitere Phanomene existie-
ren, die nicht durch Linearisierung erfassbar sind. Ein Beispiel sind Grenz-
zyklen, das sind periodische Bewegungen, in die die Losungen hineinlau-
fen (stabiler Grenzzyklus, Abb. 13) oder die Losungen abstoen (instabiler
Grenzzyklus).
~
y
24
Existenz von Grenzzyklen
In der Ebene kann die Existenz von Grenzzyklen mithilfe des folgenden
Satzes (ohne Beweis) nachgewiesen werden:
Die Existenz von Grenzzyklen ist aber nicht auf die Ebene beschrankt. Im
Rn erfolgt der praktische Nachweis eines Grenzzyklus u
ber die Technik des
PoincareSchnitts.
(y0 )
P
y0
25
Beispiel 2.6: Anwendung der Aussagen auf den Brusselator (1.5)
Bedingung f
ur kritische Punkte
0 = A + y12 y2 (B + 1)y1 ,
0 = By1 y12 y2 .
( )
A
Einziger kritischer Punkt ist ys = B/A .
( )
B 1 A2
Die Eigenwertanalyse von DF (ys ) = ergibt:
B A2
ur B < A2 + 1
stabiler Knoten f
ys : Zentrum ur B = A2 + 1
f
instabil ur B > A2 + 1
f
Fall B A2 + 1 damit klar.
Fall B > A2 + 1: Man kann zeigen, dass alle Losungen beschrankt sind. Nach PoincareBendixson
mu daher ein Grenzzyklus existieren! (Fall (i) im Satz 2.6 kann nicht zutreen, Fall (ii) nur,
wenn der AW bereits auf dem Grenzzyklus liegt. Ein Zentrumsknoten mit periodischer Losung
kommt nicht in Frage.) Abbildung 15 zeigt qualitativ das Phasenportrait f
ur A = 1, B = 3, d.h.
B > A2 + 1.
y
2
3
y
~ s
y
^
y
1 y
1
Abbildung 15: Skizze des stabilen Grenzzyklus beim Brusselator. y wird von ys abgestossen,
lauft in Grenzzyklus hinein, dagegen lauft y von auen in Grenzzyklus.
Bemerkungen:
Der Ubergang
B A2 + 1 nach B > A2 + 1 markiert eine sogenannte
HopfBifurkation, eine Verzweigung, in der stabile Losungen in peri-
odische Losungen u
bergehen.
26
Ab Dimension n = 3 existieren noch komplexere Phanomene, die meist
nur durch numerische Studien zu entdecken sind.
Kritische Punkte und Grenzzyklen nennt man auch Attraktoren des dy-
namischen Systems. Es gibt daru
berhinaus strange attractors, die bei
chaotischen Systemen auftreten. Ein beru
hmtes Beispiel ist die Lorenz
Dierentialgleichung
y1 = by1 + y2 y3 ,
y2 = y2 + y3 ,
y3 = y1 y2 + ry2 y3 .
40
30
20
y(3)
10
10
20
10
50
0 40
30
10
20
y(2) 20 10
y(1)
27
Der Begriff der Evolution
x,x0 y0 := y(x),
Fu
r weitere Details siehe Deuhard/Bornemann, Kap. 2.
28
Kapitel 3
Einschritt- und Extrapolationsverfahren
Dieses Kapitel stellt eine erste und sehr wichtige Verfahrensklasse zur Losung
von Anfangswertproblemen vor, die sogenannten Einschrittverfahren. Allge-
mein basieren numerische Verfahren fu r AWPe auf folgendem Ansatz: Man
fu
hrt ein Gitter auf dem betrachteten Intervall [a, b] ein,
d.h. man diskretisiert das Kontinuum [a, b]. Die Schrittweiten hi = xi+1
xi mussen nicht aquidistant sein. Dann berechnet man, ausgehend von y0 ,
.
sukzessive Nahrungen yi = y(xi ) an den diskreten Punkten,
y0 y1 y2 . . . yn .
29
y
. yi+1
.y(xi+1)
yi .
h x
x x
i i+1
3.1 Beispiele fu
r Einschrittverfahren
y1 = y0 + hf (x0 , y0 ),
y2 = y1 + hf (x1 , y1 ),
.. ..
. .
yi+1 = yi + hf (xi , yi ). (3.1)
30
xi+1
(iii) Quadraturformel y(xi+1 ) = y(xi ) + f (t, y(t)) dt
x
|i {z }
.
=hf (xi ,yi )
bei der man deutlich die Verwandschaft zur numerischen Quadratur sieht.
Als Notation fu
r Einschrittverfahren (ESV) verwendet man das Schema
31
Wie gut approximiert yi+1 die exakte Losung y(xi+1 )? Zur Analyse der Ver-
fahren fu
hrt man verschiedene Begrie ein. Der erste ist eine lokale Eigen-
schaft:
Nach (3.5) gibt (h) die Dierenz von exakter Losung und numerischer
Approximation, skaliert mit h, an. Eine zweite Interpretation ist
y(x0 + h) y0 y1 y0
(h) = .
| h
{z h }
} | {z
Sehnensteigung Steigung
exakt Approx.
32
Definition 3.2 Konsistenzordnung
Ein Verfahren heit konsistent, falls der lokale Diskretisierungsfehler fu
r
h 0 ebenfalls gegen 0 strebt:
(h) = O(hp ) .
lim en (X) = 0 .
n
Zum diskreten Wert yi sei ui (x) die exakte Losung der Dierentialgleichung
mit Anfangswert ui (xi ) = yi fu
r i = 0, . . . , n. Insbesondere ist u0 = y die
exakte Losung zum Anfangswert y0 . Der globale Fehler lautet dann
y(xn ) yn = u0 (xn ) yn
33
y
. u 2 (xn )
. u 1 (x n )
y2 .
y . . u 0 (x n )
y0 . 1
x
x0 x1 x2 xn =X
Damit folgt
n2
u0 (xn ) yn un1 (xn ) yn + ui (xn ) ui+1 (xn )
i=0
n2
un1 (xn ) yn + ui (xi+1 ) ui+1 (xi+1 )eL|xn xi+1 | .
i=0
Bei der letzten Abschatzung wurde Satz 2.2 angewandt, um die Die-
renz ui (xn ) ui+1 (xn ) auf die Dierenz der Anfangswerte ui (xi+1 )
ui+1 (xi+1 ) zuru
ckzuspielen. Wegen yi+1 = ui+1 (xi+1 ) hat man nun
n1
u0 (xn ) yn ui (xi+1 ) yi+1 eL|xn xi+1 | .
i=0
Die Beitrage auf der rechten Seite sind lokale Fehler nach einem Schritt.
Fu
r ein konsistentes Verfahren ist demnach
ui (xi+1 ) yi+1 c hip+1 c hpmax hi
34
mit hmax = max hi . Schlielich bekommen wir fu
r den globalen Fehler
i
n1
u0 (xn ) yn c hpmax hi eL|xn xi+1 |
i=0
xn
c hpmax eL(xn t) dt
x0
unter der Annahme xn > x0 . Mit Ausintegrieren lautet das Resultat:
35
3.3 RungeKuttaVerfahren
K2
y
1
K1
y0
x x
0 1
Diskretisierungsvorschrift
y1 = y0 + 12 h(K1 + K2 )
mit Korrekturen K1 = f (x0 , y0 ),
K2 = f (x1 , y0 + h K1 ),
bzw. mit Verfahrensfunktion :
y1 = y0 + h (x0 , y0 , h) ,
(x, y, h) = 21 [f (x, y) + f (x + h, y + hf (x, y))]
= 21 (K1 + K2 )
36
Die Verfahrenskoezienten bi , ci fu
r i = 1, . . . , s sowie aij fu
r i = 1, . . . , s
und j i 1 legen die Methode fest, zusammen mit der Stufenzahl s.
Beispiel 3.3: HeunVerfahren
Stufenzahl s = 2 und Koezienten
1 1
b1 = , b2 = , c1 = 0, c2 = 1, a21 = 1.
2 2
y1 = y0 + h K2
mit K1 = f (x0 , y0 ),
( )
K2 = f x0 + h2 , y0 + h
2 K1 .
c1 0
c2 a21 ... 0
.. .. ... ... bzw. c A
. .
cs as1 . . . ass1 0 bT
b1 b2 . . . bs
0 0
1 1
1 1 Heun 2 2
mod. Euler
1 1
2 2 0 1
0
1 1
2 2
1 klassisches RungeKuttaVerfahren
2 0 12
(Runge 1895, Kutta 1901)
1 0 0 1
1 2 2 1
6 6 6 6
37
Koeffizientenbestimmung
Wie kommt man auf die Koezienten? Die Koezienten sind zunachst
freie Parameter des Verfahrens. Sie werden so bestimmt, dass das Verfahren
moglichst gut ist, d.h. eine moglichst hohe Ordnung erreicht. Wesentliches
Hilfsmittel bei der Bestimmung sind Taylorreihenentwicklungen der exakten
Losung und der numerischen Naherung.
Als Beispiel fu
r eine Koezientenbestimmung (den sogenannten Abgleich)
betrachte man den Fall s = 2:
y1 = y0 + h(b1 K1 + b2 K2 ) ,
K1 = f (x0 + c1 h, y0 ) , K2 = f (x0 + c2 h, y0 + ha21 K1 ).
Taylorentwicklung liefert
38
Oensichtlich wird bereits bei s = 3 Stufen der Abgleich sehr aufwandig.
Das Vorgehen vereinfacht sich etwas, falls man nur autonome Dierenti-
algleichungen betrachtet. Man wird von einem fu r autonome Dierential-
gleichungen hergeleiteten Verfahren natu
rlich verlangen, dass es auch im
nichtautonomen Fall wohldeniert ist (d.h. invariant unter Autonomisie-
rung ist).
i1
denn wegen der letzten Zeile in F ergibt sich x0 + ci h = x0 + h aij 1 als
j=1
Forderung.
Damit genu gt es, den Abgleich nur im autonomen Fall durchzufu hren und
nachtraglich mit den Koezienten aij die Zwischenstutzstellen oder Knoten
ci mittels (3.9) zu denieren. Aber diese Erleichterung hilft nicht allzuviel,
um Verfahren mit s 3 Stufen und entsprechend hoherer Ordnung zu
gewinnen. Auch wenn die einzelnen Schritte nur aus wiederholter Taylo-
rentwicklung bestehen, ist der Abgleich auf diese Art und Weise schnell
nicht mehr handhabbar. Der folgende Satz gibt die Bedingungsgleichungen
bis Ordnung 4 an.
39
den Bedingungsgleichungen erfu
llt sind:
s
p=1: bi = 1 (1)
i=1
s
p=2: bi ci = 1/2 (2) sowie (1)
i=1
s
s
p=3: bi c2i = 1/3, bi aij cj = 1/6 (3) sowie (1) (2)
i=1 i,j=1
s s
p=4: bi c3i = 1/4, bi ci aij cj = 1/8,
i=1 i,j=1
s s
bi aij c2j = 1/12, bi aij ajk ck = 1/24 sowie (1) (3)
i,j=1 i,j,k=1
Auf diese Art und Weise gelingt es, die fortgesetzte Taylorentwicklung syste-
matisch und elegant darzustellen, und zwar sowohl fu r die exakte Losung als
auch fu
r die numerische Approximation. (Naheres in Deuhard/Bornemann
und Hairer/Wanner)
40
Verfahrenskonstruktion
Damit ist a43 = 1 und a42 = 1. Aus der Denition der ci ergeben sich
schlielich noch a21 = 1/3, a31 = 1/3, a41 = 1.
41
Resultat ist die Kuttasche 3/8-Regel
0
1 1
3 3
2
3 13 1
1 1 1 1
1 3 3 1
8 8 8 8
(Zu uberpru
fen ist noch die u
briggebliebene Bedingung bi aij c2j = 1/12,
die von diesen Koezienten tatsachlich erfu
llt wird.)
Wie oben erwahnt kann man mit Hilfe des Abgleichs u ber die Butcher-
Reihen die Bedingungsgleichungen fu r Verfahren hoher Ordnung ( 5) zwar
elegant aufstellen, die Zahl der Bedingungsgleichungen nimmt jedoch stark
zu:
Ordnung p 1 2 3 4 5 6 7
Mindeststufen s 1 2 3 4 6 7 9
# Bed.gleich. 1 2 4 8 17 37 85
# Zahl Koe. 1 3 6 10 21 28 45
Die Koezientenbestimmung ist dann nur noch mit weiteren vereinfachen-
den Annahmen (vergl. den Zugang u ber die Quadratur) und mit Computer
AlgebraUnterstu tzung (Maple, Mathematica) moglich. In der folgenden
Tabelle ist aufgefu
hrt, wieviele Stufen notwendig sind, um eine bestimmte
Ordnung zu erreichen. Oensichtlich gilt die Ungleichung s p.
Stufenzahl s 1 2 3 4 5 6 7 8
max. Ordnung p 1 2 3 4 4 5 6 6
42
3.4 Ru
ckw
artsanalyse des globalen Diskretisierungsfehlers
Nachdem wir nun mit den Runge-Kutta-Verfahren eine groe Klasse von
Methoden zur Verfugung haben, soll die in Abschnitt 3.2 gemachte Konver-
genzaussage um eine relativ neue Sichtweise erganzt werden. Ziel ist eine
Interpretation des globalen Fehlers im Sinne der Ru ckwartsanalyse. Wir
stellen uns dazu die folgende Frage: Kann man die numerische Losung yn
nach n Schritten als exakte Losung zu einer gestorten Dierentialgleichung
auassen? Und wie sieht diese Dierentialgleichung aus?
Wir betrachten dazu den autonomen Fall mit konstanter Schrittweite. Die
Verfahrensvorschrift eines expliziten ESV lautet dann
Die Funktionen di sind aus dem Abgleich bekannt und setzen sich aus ele-
mentaren Dierentialen zusammen.
Beispiel 3.5: Beim expliziten Euler ist d2 = d3 = 0 da (y, h) = f (y). Beim Verfahren von
Heun dagegen ist
1 1
(y, h) = (f (y) + f (y + hf (y))) = (f (y) + f (y) + hf f (y) + O(h2 ))
2 2
y = f ( y ) + h2 f3 (
y ) + hf2 ( y) + . . .
43
u
ber Taylorentwicklung darstellen als
y(xi + h) = yi + h(f (yi ) + hf2 (yi ) + h2 f3 (yi ) + . . .)
h2
+ (f (yi ) + hf2 (yi ) + . . .)(f (yi ) + hf2 (yi ) + . . .) + . . (3.10)
.
2
Auf der anderen Seite ist yi+1 u
ber
yi+1 = yi + hf (yi ) + h2 d2 (yi ) + h3 d3 (yi ) + . . . (3.11)
gegeben, so dass aus der Forderung yi+1 = y(xi + h) durch Vergleich der
Reihenentwicklungen die unbekannten Funktionen fi folgen. So ist z.B.
1
f2 (y) = d2 (y) f f (y). (3.12)
2
Bei einem Verfahren der Ordnung p ist nun
yi+1 = y(xi + h) + hp+1 p+1 (y(xi )) + O(hp+2 ),
so dass die fu
hrenden Terme in der Dierenz beider Entwicklungen bis zum
Index p herausfallen: fi = 0 fu
r i = 2, . . . , p. Zusammenfassend folgt
1) Die in Satz 3.1 mithilfe einer Vorwartsanalyse und dem Lemma von
Gronwall gegebene Abschatzung fu r den globalen Diskretisierungsfeh-
L|Xx0 |
ler lautete en (X) chmax (e
p
1)/L. Fur groe Lipschitzkon-
stanten und lange Integrationsintervalle ist diese Abschatzung nicht be-
sonders aussagekraftig und oft viel zu pessimistisch. Die Ru
ckwartsana-
lyse erlaubt dagegen eine qualitative Interpretation.
44
2) Dieses relativ neue Resultat entfaltet seine ganze Wirkung erst im Zu-
sammenhang mit den sogenannten geometrischen Integratoren, siehe
Kap. 6. Literatur dazu: Hairer/Lubich/Wanner, Geometric Numerical
Integration, Springer 2002
3) Die Reihendarstellung des gestorten AWP mu nicht konvergieren,
vergl. die Problematik der asymptotischen Entwicklung. Sie kann aber
als Ausgangspunkt fu r den Nachweis einer solchen Entwicklung des
globalen Fehlers herangezogen werden.
3.5 Fehlersch
atzung und Schrittweitensteuerung
45
Gesamtfehler
p
e (x) = O(h )
n
Rundungsfehler
Schauen wir uns nun eine spezielle Losungskurve an, so wird man oensicht-
lich fordern, dass starke Anderungen in der Losung ein lokal feines Gitter
verlangen, wahrend annahernd konstante Abschnitte mit groen Schritten
durchlaufen werden konnen, Abb. 21.
x
groe Schrittweiten
kleine Schrittweiten
Da das Losungsverhalten aber a priori nicht bekannt ist (im gewissen Ge-
gensatz zur Quadratur, bei der man den Integranden kennt), kann die
Gitterstruktur nicht zu Beginn der numerischen Integration festgelegt wer-
den. Stattdessen wird man die Gitterpunkte wahrend des Losungsprozesses
schrittweise anpassen, d.h. ausgehend von (xi , yi ) wird eine neue Schrittwei-
te hi bestimmt und damit xi+1 festgelegt. Man spricht deswegen von einer
Schrittweitensteuerung, kurz SWS.
46
Das Ziel der SWS kann nun so formuliert werden: Wahle hi so, dass die
nachste Naherung yi+1 eine vorgegebene Fehlertoleranz erfu
llt.
Allerdings steuern alle gangigen Verfahren nur den lokalen Fehler und nicht
den globalen! D.h., die Schrittweite hi wird so gewahlt, dass fu
r den aktu-
ellen Schritt gilt
u(xi+1 )|u(xi )=yi yi+1 (3.13)
mit vorgegebener Toleranz . Die exakte Losung u(x) zum AW u(xi ) = yi ist
im Kriterium (3.13) natu rlich nicht bekannt und muss durch eine Approxi-
mation yi+1 ersetzt werden. Man spricht deswegen von der Fehlerschatzung,
vergleiche das Vorgehen beim adaptiven Simpson.
Formal im ButcherTableau
c A
s
bT
y1 = y0 + h bi Ki
i=1
bT s
bi Ki
y1 = y0 + h
i=1
Die Idee der Einbettung geht auf Fehlberg zuruck, und die von ihm ent-
wickelten adaptiven Verfahren nennt man deswegen RungeKuttaFehlberg
Verfahren. Details dazu und eine moderne Implementierung werden weiter
unten besprochen. Im folgenden konzentrieren wir uns auf die Grundidee
der Schrittweitensteuerung und nehmen an, eine zweite Approximation yi+1
der Ordnung p + 1 stehe zur Verfugung.
47
Mit den Naherungen yi+1 und yi+1 setzt man nun folgendermaen an, um
eine optimale Schrittweite hopt. zu bestimmen: Es gilt
.
yi+1 yi+1 = c hp+1 (yi+1 Konsistenzord. p).
c hp+1
opt.
=
Fur den Einsatz in der Praxis sind zahlreiche Heuristiken notwendig, um die
Leistungsfahigkeit der SWS zu optimieren und Sonderfalle abzufangen. We-
sentlich sind z.B. die Begrenzung der neuen Schrittweite nach oben/unten
und die Einfuhrung eines Sicherheitsfaktors.
48
Algorithmus 3.1 Einschrittverfahren adaptiv
Start: x0 , y0 , h0 , , xend gegeben;
i = 0;
while xi < xend
x = xi + hi ;
y = yi + hi (xi , yi , hi );
i , yi , hi );
y = yi + hi (x
e = y y;
hopt. = h p+1 ; % Sicherheitsf.
e
hopt. = min( hi , max(hi , hopt. )); % Schranken ,
if e
xi+1 = x;
yi+1 = y;
hi+1 = min(hopt. , xend xi+1 );
i = i + 1;
else
hi = hopt. ; % Wiederhole Schritt
end
end
Bemerkungen:
49
In Algorithmus 3.1 wird die weniger genaue Losung als neue Approximation
herangezogen, denn fu r sie liegt mit der Dierenz y y ein Fehlerschatzer
vor. Diese Argumentation gilt heute als u berholt, und man rechnet in den
gangigen Verfahren mit der besseren Losung weiter, ein Vorgehen, das man
auch als lokale Extrapolation bezeichnet. Begru ndung: Der Fehlerschatzer
bezieht sich auf den lokalen Fehler und liefert kaum Information u ber den
eigentlich wichtigeren globalen Fehler. Stattdessen wird u ber die SWS das
Gitter an den lokalen Losungsverlauf angepasst, und diese Anpassung funk-
tioniert, wie die Erfahrung zeigt, mit der lokalen Extrapolation genauso gut,
mit dem zusatzlichen Plus der hoheren Ordnung.
Damit fuhren wir die folgende Verbesserung im Algorithmus 3.1 ein: ist ein
Verfahren der Ordnung p + 1 und eines der Ordnung p. Dann arbeitet der
Algorithmus automatisch im Modus der lokalen Extrapolation; alle anderen
Schritte bleiben wie gehabt.
Eingebettete Verfahren
Wie oben bereits erwahnt, stellt die Einbettung eines zweiten RK-Verfahrens,
das moglichst wenig zusatzlichen Aufwand erfordert, die gangige Technik
zur SWS dar. Am Beispiel des klassischen RK-Verfahrens mit p = 4 wollen
wir die prinzipielle Vorgehensweise studieren. Im Tableau haben wir also
eine zusatzliche Zeile
0
1 1
2 2
1 1
2 0 2
1 0 0 1
1 2 2 1
6 6 6 6
b1 b2 b3 . . . bs
50
Erster Ansatz: Wir nehmen s = 4 an und versuchen, die vier Koezienten
b1 , . . . , b4 so zu bestimmen, dass
y1 = y0 + h(b1 K1 + . . . + b4 K4 )
wobei die Koezienten aij und ci bereits eingesetzt wurden. Einzige Losung
des linearen Gleichungssystems ist b = b, d.h. mit s = 4 ist es nicht moglich,
solch ein eingebettetes Verfahren zu konstruieren!
51
liefern. Ein Beispiel zeigt das folgende Butcher-Tableau:
0
1 1
2 2
1
20 12
1 0 0 1
1 16 26 26 1
6
1 2 2 1
6 6 6 6 0
1 2 2 1
6 6 6 0 6
Die Auswertung der Dierenz y y reduziert sich hier auf die einfache
Formel y y = (K1 (n.S.) K4 (a.S.))/6, d.h. die Groe y muss gar nicht
explizit berechnet werden.
Bemerkung:
Eine in der Herleitung einfache, aber in der Auswertung teure Alternative
zu den eingebetteten Verfahren ist die Richardson-Extrapolation, bei der
ein Schritt y1 = y0 + 2h(x0 , y0 , 2h) mit doppelter Schrittweite und zwei
Schritte y1/2 = y0 +h(x0 , y0 , h), y1 = y1/2 +h(x0 +h, y1/2 , h) mit einfacher
Schrittweite miteinander kombiniert werden, um mittels Extrapolation eine
Approximation hoherer Ordnung zu erhalten. Siehe dazu den Abschnitt
u
ber Extrapolation unten.
Das Verfahren hat die Ordnung 5 mit einem eingebetteten Verfahren der
Ordnung 4, daher stammt die Namensgebung DOPRI5(4). Erreicht wird
dies mit s = 6 und s = 7 Stufen, wobei durch den FSAL-Ansatz die
zusatzliche Stufe nicht ins Gewicht fallt. Man hat also eektiv nur 6 Stufen
52
pro Schritt auszuwerten. Der Koezientensatz lautet
0
1 1
5 5
3 3 9
10 40 40
4 44 56 32
5 45 15 9
8 19372 25360 64448 212
9 6561 2187 6561 729
9017 355 46732 49 5103
1 3168 33 5247 176 18656
35 500 125 2187 11
1 384 0 1113 192 6784 84
35 500 125 2187 11
384 0 1113 192 6784 84 0
5179 7571 393 92097 187 1
57600 0 16695 640 339200 2100 40
(i) Der fuhrende Fehlerterm 6 (y) des Verfahrens 5. Ordnung soll moglichst
hrende Fehlerterm 5 (y) des Verfahrens 4. Ordnung dage-
klein sein, der fu
gen soll moglichst dominant sein im Vergleich zu den weiteren Termen der
zugehorigen Fehlerentwicklung.
(ii) Der Term 5 (y), der sich ja aus elementaren Dierentialen von f zusam-
mensetzt, soll im Sonderfall der Quadratur f = f (x) nicht verschwinden.
(iii) Die Koezienten ci sollen in [0, 1] liegen und verschieden sein.
Einen guten Integrationscode zeichnen aber nicht nur die Koezienten aus.
Zwei weitere Aspekte sind von besonderer Bedeutung:
53
exakte Losung moglichst gut approximiert. Grundsatzlich geht man dabei
intervallweise vor und konstruiert nach jedem Schritt den Interpolanten auf
[xi , xi+1 ].
Die einfachste Methode bei ESV verwendet die Daten yi , fi sowie yi+1 , fi+1
und konstruiert das HermiteInterpolationspolynom vom Grad 3
u() = (1 )yi + yi+1 + ( 1) ((1 2)(yi+1 yi ) + ( 1)hfi + hfi+1 )
mit 0 1. Probe: u(0) = yi , u (0) = yi = fi , u(1) = yi+1 , u (1) = fi+1
Die Formel (3.14) ist in diesem Zusammenhang ein I-Regler , d.h. sie beruck-
sichtigt die Summe der bisherigen Abweichungen vom Sollwert (das Inte-
gral). Ein PI-Regler ist zusatzlich Proportional zur letzten Abweichung vom
Sollwert und aus regelungstechnischer Sicht besser. Man kommt so auf die
Schrittweitenformel
( ) ( )
yi yi
hi+1 = hi (3.16)
yi+1 yi+1
mit Konstanten = 0.7/p, = 0.4/p. Speziell fu
r DOPRI5(4) in der Im-
54
plementierung von Hairer kann wahlweise die PI-Steuerung (3.16) oder die
klassische Variante (3.14) eingesetzt werden, mit zusatzlichen Sicherheits-
faktoren.
Numerische Software
Beispiel 3.6: In Abbildung 22 sieht man die numerische Losung der Van-der-Pol-Gleichung in
der Fassung (2.2) mit = 2, AW y0 = (3, 0)T und Intervall [0, 14]. Der Integrator ode45 aus MAT-
LAB lost die Dierentialgleichung mit der Standardeinstellung ATOL=1.e6 und RTOL=1.e3.
Mit markiert sind die diskreten Approximationen, so dass man sich einen Eindruck von den
Schrittweitenanderungen und der Gitteranpassung machen kann. Integrationsstatistik: NST = 58
(Gesamtzahl Schritte), NEF = 8 (verworfene Schritte).
55
4
y2 1
4
3 2 1 0 1 2 3
y1
3.6 Extrapolationsverfahren
Das Prinzip der Extrapolation hatten wir schon in Numerik II bei der Rom-
bergquadratur kennengelernt. Grundlage jeder Extrapolationstechnik ist die
Existenz einer asymptotischen Entwicklung. Fu r die Klasse der ESV reicht
dafu
r bereits eine Darstellung des lokalen Fehlers als Reihe
y(x) + h(x, y(x), h) = y(x + h) + hp+1 p+1 (y(x)) + hp+2 p+2 (y(x))
+ . . . + hN +1 N +1 (y(x)) + O(hN +2 )
(3.17)
aus. Es gilt der folgende Satz:
56
r 0 h h0 beschrankt ist.
wobei das Restglied E(X, h) fu
Beweis: Wir schreiben yh (X) = yn f
ur die numerische Approximation und suchen im ersten
Schritt eine Funktion ep mit
Die Summe yh (x) := yh (x) + hp ep (x) interpretieren wir als neues Verfahren mit der diskreten
i , yi , h). Es folgt f
Vorschrift yi+1 = yi + h(x die Beziehung
ur die Inkrementfunktion
(x, y, h) = (x, y hp ep (x), h) + hp1 (ep (x + h) ep (x)),
ep (x) = f (x, y(x))ep (x) p+1 (y(x))
y
Die nachste Koezientenfunktionen ep+1 bestimmt man analog, ausgehend vom neuen gedach-
ten Verfahren .
57
1) Wahle Mikroschrittweiten hi = H/ni mit der Makroschrittweite H und
z.B. ni {1, 2, 3, 4, 5, . . .} (Harmonische Folge) oder
ni {1, 2, 3, 4, 6, 8, 12 . . .} (Bulirsch-Folge).
3) Das Polynom P (h) vom Grad k 1 wird so gewahlt, dass es die diskreten
Approximationen interpoliert, P (hi ) = Ti,1 fu
r i = 1, . . . , k. Die Extrapola-
tion zum Grenzwert h = 0 eliminiert alle Terme der Expansion und liefert
(Begru
ndung s.u.). Damit liegt ein Einschrittverfahren der Ordnung k vor.
4) Praktische Durchfu
hrung: Mit dem Aitken-Neville-Schema
Tj,l1 Tj1,l1
Tj,l = Tj,l1 +
nj /njl+1 1
wird das Extrapolationstableau
T1,1
T2,1 T2,2
T3,1 T3,2 T3,3
.. .. .. ..
. . . .
Tk,1 Tk,2 Tk,3 . . . Tk,k
aufgebaut. Die Spalte l enthalt dann Approximationen der Ordnung l, und
Tkk = y = P (0).
58
Begru
ndung fu
r (3.18): Nach Satz 3.4 ist
Wir gehen nun analog zur Rombergquadratur vor (dort allerdings war eine h2 -Entwicklung die
Basis). Das Polynom P hat die Lagrange-Darstellung
k
k
k
P (h) = Ti,1 li (h) = y(x0 + H) h0i li (h) e1 hi li (h) . . .
i=1 i=1 i=1
k
k
ek1 hk1
i li (h) E(x0 + H, hi )hki li (h).
i=1 i=1
Da hl sich selbst interpoliert, gilt ur l = 0, . . . , k 1, und daher
hli li (h) = hl f
k
P (0) = y(x0 + H) E(x0 + H, hi )hki li (0)
i=1
k
k hi
= y(x0 + H) + (1)k hi E(x0 + H, hi )
hi hl
i=1 i=1 l=i
Weil der letzte Term durch eine Konstante abgeschatzt werden kann, hi /(hi hl ) ck , und
E(x0 + H, hi ) = O(H), bekommt man schlielich die Aussage (3.18).
Der GBS-Algorithmus
Der Aufwand der Extrapolation lohnt im Vergleich zu den Runge-Kutta-
Verfahren nur, wenn eine h2 -Entwicklung (und nicht eine h-Entwicklung wie
oben) vorliegt. Das Standardverfahren in diesem Fall ist die GBSMethode
nach Gragg/Bulirsch/Stoer (1966).
59
erhalt so yni (hi ) als Naherung fu
r y(X). Fu
r genu
gend glattes f und ni gera-
de besitzt yn (h) eine asymptotische Entwicklung analog der Trapezsumme
in der Rombergquadratur (siehe Stoer-Bulirsch, Kap. 7.2)
y(X) yn (h) = h2 1 (X) + . . . + h2N N (X) + h2N +2 (X, h). (3.19)
Wiederum ist das Restglied (x, h) beschrankt. Aufgrund der h2 -Entwicklung
kann man mit wenigen Extrapolationsschritten eine sehr hohe Ordnung rea-
lisieren! Deswegen ist das auf dieser Entwicklung basierende Extrapolati-
onsverfahren, der GBS-Algorithmus, eine etablierte Methode zur Losung
gewohnlicher Dierentialgleichungen, vor allem bei hohen Genauigkeitsan-
forderungen. Details siehe Deuhard/Bornemann und Hairer/Wanner.
Bemerkungen:
1) Die explizite Mittelpunktregel ist ein 2SchrittVerfahren, siehe nachstes
Kapitel. Das aus ihr abgeleitete ExtrapolationsVerfahren kann dage-
gen wie oben als RungeKuttaVerfahren interpretiert werden!
2) Oft wird ein sogenannter Glattungsschritt nach Gragg durchgefu
hrt.
Statt yni nimmt man
1
yni := [yni 1 + 2yni + yni +1 ]
4
Die Bedeutung in der Praxis ist umstritten.
3) Da die Entwicklung (3.19) nur fur ni gerade gilt, muss die Schrittwei-
tenfolge ebenfalls auf geraden Unterteilungen von H basieren:
RombergFolge { 1, 2, 4, 8, . . .}
BulirschFolge { 1, 2, 3, 4, 6, 8, . . .}
60
Kapitel 4
Mehrschrittverfahren
Mehrschrittverfahren (MSV) bilden neben den ESV die zweite groe Verfah-
rensklasse. Bei ihnen ieen Daten / Losungspunkte aus der Vergangenheit
in die Berechnungvorschrift mit ein. Genauer:
.
In einem kSchrittverfahren berechnet man die Approximation yi+1 = y(xi+1 )
aus den Daten
Die Daten stammen aus den vorherigen Schritten oder einer Anlaufrech-
nung. Neben der Konsistenz ist die Eigenschaft der Stabilitat bei MSV
notwendig, um die Konvergenz zu zeigen. Insofern zeigt sich hier ein we-
sentlicher Unterschied zu den ESV.
4.1 Verfahrensklassen
61
Der Integrand in (4.1) hangt von der (unbekannten) Losung y ab und wird
nun durch ein Interpolationspolynom p ersetzt. Zwei Falle werden unter-
schieden:
fik+1
1111
0000
0000
1111
p(x)
fi 1111
0000
f i1
0000
1111
0000
1111 Integration y
0000
1111 i+1
0000
1111
x ik+1 x i1 xi0000
1111x x
i+1
1
k
k+s
= h ds.
=1
j
0
| {z }
=j
=:j1
62
Die Verfahrensvorschrift schreibt sich dann folgendermaen:
yi+1 = yi + h(0 fik+1 + 1 fik+2 + . . . + k1 fi ) (4.2)
As spezielle Verfahren hat man
k=1 yi+1 = yi + hfi ,
( )
k=2 yi+1 = yi + h 12 fi1 + 23 fi ,
(5 )
k=3 yi+1 = yi + h 12 fi2 12
16
fi1 + 23
12 fi .
Die so gewonnenen expliziten MSV nennt man AdamsBashforthVerfahren.
b) Ein alternativer Zugang fordert, dass das Polynom p die Werte fik+j
fu
r j = 1, . . . , k, k + 1 interpoliert. Da p damit von yi+1 abhangt, liegt ein
implizites Verfahren vor.
p(x)
111
000
f i+1
000
111
fi
000
111
000
111
000
111
000
111
x ik+1 xi x i+1 x
63
und als Verfahren
yi+1 = yi + h(0 fik+1 + . . . + k1
fi + k fi+1 ). (4.3)
Spezielle Vertreter sind
( )
k=1 yi+1 = yi + h 12 fi + 12 fi+1 Trapezregel,
( 1 )
k=2 yi+1 = yi + h 12 fi1 +, 12 fi + 12
8 5
fi+1 .
Diese Klasse von impliziten MSV bezeichnet man als AdamsMoulton
Verfahren. Der Fall k = 0 macht auch Sinn, daraus ergibt sich
yi+1 = yi + hfi+1 Impliziter Euler.
AdamsVerfahren sind in der Praxis von groer Bedeutung und werden
meist als PradiktorKorrektorSchema implementiert:
(0)
(i) yi+1 = yi + h(0 fik+1 + . . . + k1 fi )
(Startschritt mit explizitem AdamsBashforthVerfahren)
(l+1) (l)
(ii) yi+1 = (yi+1 ), l = 0, 1, . . .
Fixpunktiteration Korrektorschritt; mit implizitem AdamsMoulton
Verfahren, wobei
(y) := hk f (xi+1 , y) + yi + h(0 fik+1 + . . . + k1
fi ).
Verfahren nach Nystr om und Milne. Analog zur Herleitung der Adams
Verfahren kann man weitere MSV konstruieren, indem man y = f (x, y) von
xi1 bis xi+1 integriert, d.h. u
ber 2 Intervalle. Der Fall a) mit explizit gege-
benem Interpolanten p fu hrt dann auf die NystromVerfahren. Ein Beispiel
ist der Fall
k=2 yi+1 = yi1 + 2hfi ,
64
der die beim GBS-Algorithmus in Abschnitt 3.6 eingesetzte Mittelpunktre-
gel liefert.
Der Fall b) mit implizitem Ansatz deniert die Verfahren nach Milne. Als
Beispiel betrachte
=:j1
65
BDFVerfahren sind implizit und nden vor allem bei steifen Dierential-
gleichungen Anwendung (s. Kapitel 5). Beispiele sind
bzw. kurz
k
k
l yi+l = h l fi+l (4.5)
l=0 l=0
mit zunachst beliebigen reellen Verfahrenskoezienten l und l .
Analog zu den ESV untersuchen wir nun den Fehler eines MSV in einem
Schritt. Dazu gibt es verschiedene Ansatze, die im Wesentlichen aquivalent
sind. Wir verwenden hier den Defekt als lokale Fehlerdenition, vgl. Deni-
tion 3.1.
Definition 4.1 Lokaler Diskretisierungsfehler MSV
Sei y(x) exakte Losung des AWP y = f (x, y), y(x0 ) = y0 . Der lokale Dis-
kretisierungsfehler des MSV (4.5) mit konstanter Schrittweite h ist definiert
als Defekt
( k )
1
k
(h) := l y(xi+l ) h l f (xi+l , y(xi+l )) .
h
l=0 l=0
66
Das Fehlerma (h) ergibt sich also durch Einsetzen der exakten Losung in
die Dierenzengleichung (4.5). Im Sonderfall exakter Daten
hat man fu
r die numerische Losung unter der Annahme k = 0 (explizites
Verfahren)
k1
k1
k yi+k + l yi+l = h l fi+l
l=0 l=0
und fu
r die exakte Losung
k1
k1
k y(xi+k ) + l y(xi+l ) = h l fi+l + h (h)
l=0 l=0
k
= (h) = (y(xi+k ) yi+k ) . (4.6)
h
Die Darstellung (4.6) ist analog zur Denition 3.1 des lokalen Fehlers bei
ESV. Oft normiert man das MSV noch durch k := 1, so dass (4.6) un-
abhangig von Verfahrenskoezienten geschrieben werden kann.
(h) = O(hp ) .
Zur Bestimmung der Konsistenzordnung eines MSV fu hrt man ahnlich den
RKVerfahren einen Abgleich durch und ermittelt Bedingungsgleichungen
fu
r die Koezienten l , l . Man setzt an
( k )
1 k
l y(x + lh) h l y (x + lh) = (h)
h
l=0 l=0
67
1
mit y(x + lh)) = y(x) + y (x)lh + y (x)(lh)2 + . . .
2
1
y (x + lh)) = y (x) + y (x)lh + y (x)(lh)2 + . . .
2
Sortieren nach hPotenzen ergibt
1 k k
(h) = y(x) l + y (x) (l l l )
h
l=0 l=0
h k
+ y (x) (l l2 2l l)
2
l=0
..
. (4.7)
hp1 (p) k
+ y (x) (l lp p l lp1 )
p!
l=0
p
+ O(h )
Fur die Ordnung p 1 mussen auch die hoheren Terme in der Entwicklung
(4.7) herausfallen, d.h.
k
k
q
l l = q l lq1 fu
r q = 1, . . . , p. (4.9)
l=0 l=0
68
Also liegt ein 2SchrittVerfahren der Konsistenzordnung p = 3 vor.
() : = k k + . . . + 1 + 0 , (4.10)
() : = k k + . . . + 1 + 0 , (4.11)
1
k
+ h2 (l l2 2l l)
2
l=0
..
. (4.13)
hp
k
+ (l lp pl lp1 ) + O(hp+1 ).
p!
l=0
69
Satz 4.1 Konsistenzordnung MSV
Das MSV l yi+l = h l fi+l hat die Konsistenzordnung p, falls eine der
drei aquivalenten Bedingungen erfullt ist:
k
k
q
k
(i) l = 0 und l l = q l lq1 fu
r q = 1, . . . , p
l=0 l=0 l=0
Bemerkungen
1) In der Literatur schreibt man oft (x, y, h) statt (h), um die Abhangig-
keit vom Gitter und der Losung zu betonen. Bedingung (ii) oben be-
deutet dann (x, exp, h) = O(hp ) . Ein anderer Ansatz verwendet statt
(x, y, h) den Dierenzenoperator
k
k
L(x, y, h) := l y(x + lh) h l y (x + lh) = h (x, y, h).
l=0 l=0
2) Einordnung Verfahrensklassen:
70
Definition 4.3 Globaler Diskretisierungsfehler
Der globale Diskretisierungsfehler eines MSV ist die Differenz
yi = hf (yi ). (4.14)
4.3 Stabilit
at von MSV
Wie hangen lokaler Fehler (d.h. die Konsistenz) und globaler Fehler (d.h.
Konvergenz) zusammen? Anders als bei ESV braucht man bei MSV eine
zusatzliche Stabilitatsbedingung, um Konvergenz zu zeigen.
71
Beispiel 4.2 Als Beispiel f
ur die Problematik betrachte das MSV mit k = 2 und p = 3 aus
Beispiel 4.1,
yi+2 + 4yi+1 5yi = h(4fi+1 + 2fi ).
Erzeugende Polynome sind () = 2 + 4 5 sowie () = 4 + 2. Das MSV wird angewendet
auf die triviale (!) Dierentialgleichung
y = 0 , AW y0 = 1 , sowie y1 = 1 + h. (4.15)
Die im Beispiel diskutierte Problematik tritt bei vielen MSV auf. Verfahren
mit scheinbar guten Eigenschaften wie hoher Konsistenzordnung und kleiner
Fehlerkonstante liefern schon bei einfachsten Testbeispielen extrem schlecht
Resultate.
Zur Analyse des Stabilitatsproblems wird wie im Beispiel 4.2 die triviale
Testgleichung y = 0, y(x0 ) = 1 betrachtet. Ein darauf angewandtes MSV
liefert die homogene Rekursion oder Differenzengleichung
0 yi + 1 yi+1 + . . . + k yi+k = 0 (4.18)
72
mit Startwerten y0 , . . . , yk1 .
Satz 4.2
Sei eine m-fache Nullstelle des charakteristischen Polynoms () aus
(4.10), d.h. () = () = . . . = (m1) () = 0. Dann gilt:
(ii) Die allgemeine Losung der homogenen Rekursion (4.18) ist eine
Linearkombination der insgesamt k speziellen Losungen nach (i).
(j)
Beweis: Betrachte spezielle Losung yi mit 1 j m,
(j) (j) (j)
0 yi + 1 yi+1 + . . . + k yi+k = Dj1 (0 i + 1 i+1 + . . . + k i+k ) = Dj1 (() i ).
folgt
( )
j1 i j1 i j1
D (() ) = () D + ... + Ds () Dj1s i + . . . + Dj1 () i = 0.
|{z} s | {z } | {z }
=0 =0 =0
Damit ist (i) gezeigt. Oensichtlich gibt es insgesamt k spezielle Losungen entsprechend dem Grad
von , und jede Linearkombination ist ebenfalls Losung der homogenen Rekursion. Zu zeigen
bleibt, dass solch eine Linearkombination durch die Startwerte y0 , . . . , yk1 eindeutig festgelegt
ist. Beweis dazu ist elementar, aber m
uhselig (Stoer/Bulirsch, S. 148, vergl. Bsp. 4.2).
73
Falls eine NST von ist mit || > 1: {yi } = {i } wachst exponentiell!
Falls eine NST von ist mit || < 1: {yi } = {i } fallt exponentiell!
Falls () = 0 = || 1.
}
() = 0
Falls = ist einfache NST. (4.19)
und || = 1
k = 1 , k1 = 1 , k2 = . . . = 0 = 0
74
Im
x einfache NST
x evt. mehrfache NST
x
x
Re
x
und damit
() = k k2 = k2 ( 2 1) = k2 ( 1)( + 1)
r k 6, ansonsten instabil
Die BDF-Verfahren sind stark stabil nur fu
(s. Hairer/Norsett/Wanner, Kap. III.3).
75
bzw. 2k1, s.o. Diese Verfahren sind aber instabil, und es zeigt sich, dass die
in Abschnitt 4.1 eingefu
hrten Verfahrensklassen (Adams, Milne/Simpson)
bereits optimal sind.
k + 2, falls k gerade
k + 1, falls k ungerade
k, falls k /k 0 (insbesondere falls MSV explizit).
Die Ordnung k+2 kann nur erzielt werden, wenn alle Nullstellen des charak-
teristischen Polynoms auf dem Rand des Einheitskreises liegen (schwach
stabiles Verfahren).
Der Beweis der DahlquistBarriere ist sehr aufwendig und verlangt einige
funktionentheoretische Hilfsmittel, s. Hairer/Norsett/Wanner, Kap. III.3.
Wie sieht nun der allgemeine Fall, d.h. die Konvergenz fur y = f (x, y), aus?
Es zeigt sich, dass die anhand von y = 0 ermittelte Stabilitatsbedingung
(4.19) hinreichend ist, um die Konvergenz fu r ein konsistentes Verfahren zu
zeigen.
76
4.4 Konvergenz von MSV
Wie bei den Einschrittverfahren setzen wir im folgenden voraus, dass die
rechte Seite f der Dierentialgleichung global Lipschitzstetig ist und zudem
genugend glatt, damit der Begri der Konsistenzordnung Sinn macht. Im
Unterschied zu den ESV muss fu r die Konvergenz auch der Fehler in den
Startdaten y0 , . . . , yk1 beru
cksichtigt werden:
y(x0 + jh) yj 0 r h 0,
fu j = 0, . . . , k 1.
y(X) yn 0 r h 0,
fu X = x0 + nh fest.
y(X) yn chp
r Schrittweiten h h0 .
fu
(iii) Betrachte das AWP y = 1, y(0) = 0 mit exakter Losung y(x) = x. Das MSV mit Startdaten
77
ur j = 0 : k 1 lautet
yj = jh f
k1
k
k yk + h jj = h j k (yk hk) + h (1) = h(1)
j=0 j=0
und aus der Konvergenz yk kh ergibt sich (1) = (1). Zusammen mit (1) = 0 impliziert
dies die Konsistenz als notwendige Voraussetzung.
Stabilitat und Konsistenz sind also notwendig, aber sind sie auch hinrei-
chend fur die Konvergenz? Es heisst nun tief Luft holen, um den formalen
Rahmen fu r den Konvergenzbeweis bereit zu stellen. Wir schreiben das MSV
als
k1
k1
yi+k = j yi+j +h j f (xi + jh, yi+j ) + hk f (xi + kh, yi+k )
j=0 j=0
Yi := (yi+k1 , yi+k2 , . . . , yi )T ,
der alle Daten, die aus der Vergangenheit in die Berechnung von yi+k ein-
gehen, enthalt (die sog. Schleppe des MSV).
Die Grundidee fur die Konvergenzanalyse besteht nun darin, das MSV auf
ein Einschrittverfahren bez. Y zuru
ckzufu
hren. Ideal ware eine Formulie-
78
rung des MSV in der Gestalt Yi+1 = Yi + h(xi , Yi , h) mit einer Verfahrens-
funktion , doch dies ist nicht zu erreichen. Es gilt aber
k1
yi+k j=0 j yi+j + h(xi , Yi , h)
y
i+k1 yi+k1
Yi+1 = .. = .. .
. .
yi+1 yi+1
Im Fall m = 1 kann man dies neu schreiben als
mit
k1 k2 ... 1 0 (xi , Yi , h)
1 0 ... . 0 0
A :=
1 . 0 , (x, Yi , h) :=
0 .
... .. .. ..
. . .
1 0 0
Der Fall m > 1: Dann liefert das Kronecker-Produkt A B zweier Matrizen A, B die passende
Schreibweise. Es ist A B die aus Blocken ajl B zusammengesetzte Matrix,
a11 B . . . a1k B
..
A B = ... .
ak1 B . . . akk B
A und B m ussen dabei nicht unbedingt quadratische Matrizen sein. Die Verfahrensvorschrift
(4.21) lautet f
ur m > 1
Im folgenden wollen wir uns der Einfachheit halber auf den Fall m = 1 zur
uckziehen und uns mit
dem Hinweis begn
ugen, dass f
ur m > 1 die Argumentation mithilfe des Kroneckerproduktes ganz
analog ablauft.
79
praktisch fertig. Aber: Die Matrix A enthalt mit den j die Koezien-
ten, die die Stabilitat des Verfahrens festlegen. Es gilt nun: Falls die Stabi-
litatsbedingung (4.19) vom MSV erfu llt wird, dann gibt es eine Vektornorm
auf dem Rk , so dass
Ax
A = sup 1. (4.23)
x=0 x
Man rechnet schnell aus, dass die Eigenwerte von A genau die Wurzeln des erzeugenden Polynoms
sind. Aufgrund der Stabilitatsbedingung (4.19) erf ullen alle EW |j | 1, und die Jordan-
Normalform von A ist
l+1 l+1
T 1 AT = diag 1 , . . . , l , ..
.
k1 .
k
Die mehrfachen EW seien dabei l+1 , . . . , k mit j = 0 oder j = 1 je nach Vielfachheit. Es
gibt nun ein > 0 mit 1 |j | > f
ur j = l + 1 : k (mehrfache EW liegen im Innern des
Einheitskreises). Setze ( )
D := diag 1, , 2 , . . . , k1 .
Mit einer weiteren Transformation ergibt sich
l+1 l+1
D1 JD = diag 1 , . . . , l , ..
. k1 ,
k
und in der Maximumnorm B = maxj l |Bjl | ist
D1 JD 1.
Nun denieren wir als Vektornorm x := D1 T 1 x , womit folgt
Ax = D1 T 1 AT D D1 T 1 x = D1 JD D1 T 1 x x.
80
Beweis: Wir betrachten das MSV als ESV in der Vorschrift(4.21). F
ur zwei diskrete Approxi-
mationen Yi und Zi ist
und somit
Yi+1 Zi+1 AYi Zi + h(xi , Yi , h) (xi , Zi , h).
Falls f Lipschitzstetig ist, so auch aus (4.20) und folglich (bezuglich Y ), und zwar mit
Lipschitzkonstante L . Wir nehmen die Norm mit A 1, falls die Stabilitatsbedingung (4.19)
erf
ullt ist, und folgern
Yi+1 Zi+1 (1 + hL )Yi Zi . (4.24)
Abb. 26 zeigt die weitere Vorgehensweise, die ahnlich der beim Satz 3.1 uber die Konvergenz
von ESV ist. Allerdings wird hier nicht eine Folge exakter Losungen mit Gronwall abgesch
atzt,
sondern umgekehrt eine Folge numerischer Approximationen mit (4.24).
Y(x0 ) Y(xn )
Y0
Y1
Y2 Yn
Y3 Y4
x1 x2 x3 x4 xn
(0)
Die Fortpanzung der Dierenz Y (x0 ) Y0 = Y0 Y0 in den Startdaten ergibt sich mit (4.24)
zu
Yn(0) Yn (1 + hL )n Y0 Y0 .
(0)
81
Jetzt muss nur noch aufsummiert werden:
Y (xn ) Yn = Y (xn ) Yn(n1) + Yn(n1) Yn(n2) + . . . + Yn(1) Yn(0) + Yn(0) Yn
n1
(1 + hL )n Y0 Y0 + chp+1 (1 + hL )i
(0)
i=0
n1
exp(nhL )Y0 exp(ihL )
(0)
Y0 + chp+1
i=0n
exp(nhL )Y0 exp(shL )ds
(0)
Y0 + chp+1
s=0
hp
Y0 + c (exp(nhL ) 1) .
(0)
= exp(nhL )Y0
L
(0)
ur Startdaten Y0
F Y0 = Y (x0 ) Y0 = O(hp ) folgt Konvergenz der Ordnung p.
Bemerkungen
1. Die Konvergenztheorie zu MSV kurz und pragnant zusammengefasst:
Stabilitat
Konvergenz
+ Konsistenz
2. Die gesamte hier diskutierte Theorie ist auf den Fall konstanter Schritt-
weiten beschrankt. Die Erweiterung auf variable Schrittweiten ist moglich,
verlangt aber erheblichen zusatzlichen Aufwand. Alle drei Punkte Kon-
sistenz, Stabilitat und Konvergenz mussen dazu neu aufgerollt werden!
Zum Abschluss des Kapitels u ber MSV noch ein paar Worte zur Implemen-
tierung. Im Gegensatz zu den Runge-Kutta-Verfahren ist diese sehr komplex
und lasst viel Spielraum fu
r Heuristiken. Die heute verfu
gbaren Codes glei-
chen in gewisser Hinsicht nervosen Rennpferden.
82
Ein Integrationsschritt von xi nach xi+1 verlauft nach dem Schema:
For l = 0 : m 1
Correct: Fixpunktiteration
% kSchrittAdamsMoulton,
% Koezientenber. verwendet Information des AdamsBashforth-Verf.
(l+1) (l+1)
Evaluate: Auswertung f (xi+1 , yi+1 ) =: fi+1
83
Vergleich der MATLAB-Integratoren ode45 und ode113
Was ist nun vorzuziehen, ein Runge-Kutta-Integrator wie ode45 oder ein
Adams-Verfahren wie ode113? Wir testen die MATLAB-Implementierungen
am Beispiel des Arenstorf-Problems, einer Dierentialgleichung zweiter Ord-
nung aus der Himmelsmechanik (3-Korper-Problem). Das folgende m-le
wertet die Dierentialgleichung aus.
1
10
ode23 ode45
10
10
0
RelTol
10
CPU
ode113
ode45 ode113
5
10
1
10
0 5 10 0 5 10
10 10 10 10 10 10
Fehler Fehler
84
Zum Zeitpunkt xend ist die Losung in den Komponenten y1 , y2 periodisch.
Wir bestimmen dort den absoluten Fehler der Integratoren (2-Norm), und
zwar fur immer genauer werdende Toleranzvorgaben. Zudem wird die Re-
chenzeit gemessen (auf einem Siemens-Fujitsu S6010 Notebook mit Pentium
III 1 GHz Prozessor).
Abb. 27 zeigt die so gewonnenen Ergebnisse. Links ist der tatsachliche Feh-
ler gegen die Vorgabe RelTol aufgetragen, in doppelt logarithmischer Ska-
lierung. Oensichtlich liegt der globale Fehler in der Groenordnung von
RelTol, die Integratoren arbeiten zuverlassig. Rechts ist ein Aufwandsdia-
gramm oder Work-Precision-Diagram zu sehen, mit dem Fehler gegen die
Rechenzeit. Bei niedriger Genauigkeit ist ode45 am ezientesten, bei hoher
dagegen ode113. Zum Vergleich ist auch der Runge-Kutta-Integrator ode23
mit Ordnung 2(3) eingetragen; er ist bei diesem glatten Problem chancen-
los. Bei der Analyse von ode113 fallt auf, dass dieses Verfahren sehr wenige
Funktionsauswertungen benotigt, aber durch die Koezientenberechnung
zusatzlicher Aufwand entsteht, so dass ode45 bei niedriger Genauigkeit im
Vorteil ist. Das Zittern der Werte von ode113 ist auf die komplexere Heuri-
stik mit simultaner Ordnungs- und Schrittweitensteuerung zuru ckzufu
hren.
Dadurch verhalt sich das MSV nicht so vorhersehbar - eine hohere Tole-
ranzforderung resultiert nicht unbedingt in einer besseren numerischen Ap-
proximation.
85
Kapitel 5
Numerik steifer Differentialgleichungen
Nachdem wir mit den expliziten Ein- und Mehrschrittverfahren zwei groe
Klassen von Integratoren fur gewohnliche Dierentialgleichungen kennenge-
lernt haben, wenden wir uns nun den sogenannten steifen Differentialglei-
chungen zu. Steife Dierentialgleichungen treten in vielen Anwendungen auf
und sind im Prinzip mit den bereits eingefu hrten Methoden losbar. Aller-
dings mussen explizite Verfahren die Schrittweiten bei dieser Problemklasse
sehr klein wahlen und sind nicht mehr ezient. Warum das so ist, soll
zunachst geklart werden, bevor Alternativen in Form von impliziten bzw.
semi-impliziten Ein- und Mehrschrittverfahren eingefu hrt werden.
5.1 Das Ph
anomen der Steifheit (Stiffness)
86
Van der Pol, = 1.e3, ode45 Ausschnitt
1.99
2
1 1.995
0
y
y
2
1
2.005
2
288 verworfenen Schritten. Bei genauem Hinsehen erkennt man, dass der
Integrator extrem kleine Schrittweiten wahlt, obwohl sich die Losung kaum
andert. Im Ausschnitt rechts sieht man ein Zittern der Losung, was auf ein
Stabilitatsproblem hindeutet. Andere Toleranzvorgaben beeinussen die-
ses Verhalten nur wenig, und auch die Anwendung des Adams-Verfahrens
ode113 fu hrt auf ahnliche Eekte (bei 2471 bzw. 339 erfolgreichen / ver-
worfenen Schritten).
87
diesem Beispiel fu
r den impliziten Euler beliebig gro gewahlt werden, beim
expliziten Euler dagegen nicht.
Expliziter Euler
0 h = 0.1
1 0 1 2 3 4 5
0 h = 0.2
1 0 1 2 3 4 5
0 h=1
1 0 1 2 3 4 5
x
Impliziter Euler
0 h=0.1
1 0 1 2 3 4 5
0 h=0.2
1 0 1 2 3 4 5
0 h=1
1 0 1 2 3 4 5
x
88
Sowohl bei der Van-der-Pol-Gleichung (2.3) als auch bei der Testgleichung
(5.1) spricht man von einer steifen Differentialgleichung (wobei der Grad der
Steifheit von den Parametern bzw. abhangt). Bevor der Begri steif
naher erlautert wird, eine Analyse der Diskretisierung von (5.1):
fu
r alle h > 0 erfu
llt, d.h. die numerische Losung fallt stets! Das implizite
Verfahren ist hier also u
berlegen.
89
1) Implizit besser als Explizit (s.o.)
Ein System gewohnlicher Dierentialgleichungen ist steif, wenn expli-
zite Verfahren aus Stabilitatsgru
nden extrem kleine Schrittweiten ver-
wenden, obwohl sich die Losung kaum andert; implizite Verfahren da-
gegen groe Schrittweiten einsetzen konnen.
2) Eigenwertansatz
Ein lineares System y = Ay ist steif, wenn A einige Eigenwerte j
mit stark negativem Realteil besitzt und einige mit schwach negativem
Realteil, d.h. das Verhaltnis
ist sehr gro. Im nichtlinearen Fall y = f (y) betrachtet man die Jaco-
bimatrix Df (y).
3) Singul
ar gest orte Systeme
Dierentialgleichungen der Form
y = f (y, z)
z = g(y, z) (vgl. Van-der-Pol)
Grenzschicht
asymptotische Phasen
90
Zur Auosung einer Grenzschicht benotigt jedes numerische Verfahren
kleine Schrittweiten. In der asymptotischen Phase jedoch sind groe
Schrittweiten von Vorteil, und hier sind implizite Verfahren im steifen
Fall u
berlegen.
y = (y ) + (5.2)
mit 0 und exakter Losung y(x) = (y0 (x0 ))e(xx0 ) + (x). Die
gegebene Funktion sei langsam variierend.
y * y
*
y
0 * y
0
*
*
*
*
x x1 x 2 * x x x x2 x
0 0 1
91
Obwohl das Modellproblem (5.2) (und allgemein steife Dierentialgleichun-
gen) extrem gut konditioniert ist, zeigt ein Verfahren wie der explizite Euler
ein katastrophales Verhalten.
5.2 Einschrittverfahren
Zunachst betrachten wir nochmals den expliziten Euler. Anwendung auf die
Testgleichung (5.3) liefert
y1 = y0 + h y0 = R(z) y0 mit R(z) := 1 + z, z = h .
Die Stabilitatsfunktion R des Euler-Verfahrens charakterisiert das Verhalten.
Sei C fest mit Re 0. Dann gibt es Schrittweiten h > 0 mit
|R(z)| = |R(h )| = |1 + h| > 1
92
und demnach y0 < y1 < y2 < y3 < . . ., also nicht A-stabil! Wu
rde stattdes-
sen
|R(z)| 1 z C := {z C : Re z 0}
gelten, so ware das Verfahren A-stabil.
Fu
r alle gangigen Einschrittverfahren y1 = y0 + h(x0 , y0 , h) fu
hrt die An-
wendung auf die Testgleichung (5.3) auf eine Vorschrift
y1 = R(z) y0
S := {z C : |R(z)| 1}.
Man hat die Aquivalenz
Astabil |R(z)| 1 z C C S
Das Stabilitatsgebiet des expliziten Euler ist eine Kreisscheibe mit Radius
1 und Mittelpunkt 1, siehe Abb. 32. Falls z = h S, arbeitet das
Verfahren stabil, ansonsten aber instabil!
Die Stabilitatsfunktion ist ein Polynom in z vom Grad s mit der Stufen-
r ein Polynom R kann |R(z)| 1 in C nicht erfu
zahl s. Fu llt werden
also kann das Verfahren nicht Astabil sein!
93
Im
S
i
2 Re
i
Im
S
2i
Re
94
Aus der Kenntnis des Stabilitatsgebietes kann man Aussagen u ber die ma-
ximal moglichen Schrittweiten des Verfahrens zu einer gegebenen Dierenti-
algleichung treen, vorausgesetzt man hat eine zumindest grobe Naherung
fu
r die Eigenwerte.
Beispiel 5.1: Wir betrachten das lineare System
( )
998 1998
y = Ay mit A = .
999 1999
A hat die Eigenwerte 1 = 1 und 2 = 1000. Zum Anfangswert y0 = (1, 0)T ndet man als
exakte Losung ( ) ( )
2 1
y(x) = ex + e1000x .
1 1
Die Losung besteht also aus einem langsamen und einem sehr schnell abklingenden Term. Der
ur die Steifheit verantwortlich. Diagonalisierung von A auf T 1 AT = diag (1 , 2 )
letztere ist f
und Transformation auf neue Variable u = T 1 y f uhrt auf die entkoppelten Gleichungen
ui = i ui , i = 1, 2.
Ein Runge-Kutta-Verfahren ist invariant unter einer solchen linearen Transformation. D.h., es
spielt keine Rolle, ob man das urspr
ungliche System direkt diskretisiert oder das transformierte
und dann das Ergebnis zur ucktransformiert, das Ergebnis ist dasselbe.
Im transformierten System erkennen wir die Testgleichung (5.3), und entsprechend konnen wir
nun die Schrittweiteneinschrankungen fur numerische Verfahren u
ber die Stabilitatsfunktion bzw.
direkt u
ber das Stabilitatsgebiet bestimmen. Der betragsgrote Eigenwert 2 ist entscheidend.
Im Fall expliziter Euler mu h so gewahlt werden, dass
|R(h2 )| 1 h 2/|2 | = 0.002,
denn 2 markiert den linken Rand des Stabilitatsgebietes S. Bei DOPRI5 liegt der linke Rand
bei ca. 3.3, daher folgt f
ur dieses Verfahren die Einschrankung h 3.3/|2 | = 0.0033.
95
mit Schrittweitensteuerung direkt u ber die Reaktion der Schrittweitensteue-
rung, denn diese erkennt die mogliche Instabilitat und arbeitet hart an der
Stabilitatsgrenze, es sei denn, die Genauigkeitsanforderungen verlangen per
se kleinere Schrittweiten. Haug kommt es dabei zu sehr vielen Schrittver-
werfungen und einem Zittern der Schrittweiten (vergl. Hausaufgabe H18).
Implizite Runge-Kutta-Verfahren
Der implizite Euler wurde bereits oben analysiert. Angewandt auf die Test-
gleichung (5.3) erhalten wir
1
y1 = y0 + h y1 = y1 = R(z) y0 mit R(z) := .
1z
r Re z das gleiche
Die Stabilitatsfunktion ist rational und hat fu
asymptotische Verhalten wie exp. Wegen |R(z)| 1 fu r alle z C ist
dieses Verfahren A-stabil.
ur z C hat man
(F
|R(z)|2 = 1/((1 z)(1 z)) = 1/((1 x iy)(1 x + iy)) = 1/(1 + x2 + y 2 2x) 1/(1 + x2 + y 2 ),
Gibt es weitere implizite Verfahren, die wie der implizite Euler A-stabil sind
und damit fur steife Dierentialgleichungen geeignet?
96
C
S i
1
i
Dieses Verfahren ist von der Ordnung p = 2 (Herleitung z.B. als Sonderfall
k = 1 der Adams-Moulton-Verfahren). Seine Stabilitatsfunktion ergibt sich
zu
h 1 + h/2 1 + z/2
y1 = y0 + (y0 + y1 ) y1 = y0 , , d.h. R(z) = .
2 1 h/2 1 z/2
Wiederum ist R eine rationale Funktion und das Verfahren A-stabil:
97
aber nicht beschrankt! Deswegen muss man die Kompaktifizierung C {} betrachten und deren
Rand
(C {}) = {z C : Re z = 0} {}
Die Beispiele der Trapez- und Mittelpunktregel legen nahe, die bereits ein-
gefuhrte Klasse der Runge-Kutta-Verfahren geeignet zu erweitern, um wei-
tere A-stabile Verfahren zu gewinnen. Butcher (1964) fuhrte dazu die im-
pliziten RungeKuttaVerfahren ein, kurz IRK-Verfahren. Das allgemeine
Diskretisierungsschema fu
r einen Schritt lautet
y1 = y0 + h(b1 K1 + b2 K2 + . . . + bs Ks ) (5.5)
mit Korrekturen / Zuwachsen
( )
s
Ki = f x0 + ci h, y0 + h aij Kj , i = 1, . . . , s. (5.6)
j=1
98
Stabilitatseigenschaften festzulegen. Falls die Koezienten auf der Diago-
nalen und im rechten oberen Dreieck verschwinden, aij = 0 fu r j i, liegt
der bekannte explizite Fall vor.
Aufgrund der impliziten Denition ist nicht klar, ob solche Verfahren wohl-
deniert sind, d.h., ob das nichtlineare Gleichungssystem fu r die Stufen
K1 , . . . , Ks u
berhaupt eine eindeutige Losung besitzt. Wir verschieben die-
se Frage aber auf spater und bleiben vorerst bei der Frage der Stabilitat.
99
Die Cramersche Regel liefert
( ) ( )
I hA y0 1l I hA 0
y1 = det / det
hbT y0 hbT 1
Der Nenner reduziert sich auf det(I zA), der Zahler auf
( ) ( )
I hA y0 1l I hA + h1lbT 0
det = det = y0 det(I hA + h1lbT ).
hbT y0 hbT y0
Aus der Determinantendarstellung folgt schlielich die behauptete Struktur als rationale Funktion
in z.
100
Mit anderen Worten, Rjk (z) = g(z) + O(z j+k+1 ) fur z 0. Man kann zei-
gen (Ubung H21), dass zu g(z) = exp(z) eindeutige Pade-Approximationen
existieren, die in der Pade-Tafel, Tab. 1, zusammengefasst sind. Die erste
Zeile in der Tafel (Nennergrad j = 0) enthalt die Stabilitatsfunktionen der
expliziten Runge-Kutta-Verfahren mit p = s = 1, 2, 3 und ist daher uninter-
essant. Ferner erkennt man als (1, 0)-Approximation die Stabilitatsfunktion
des impliziten Euler, als (1, 1)-Approximation in der Diagonalen die Stabi-
litatsfunktion der Trapezregel (bzw. Mittelpunktregel).
1 1+z 1+z+z 2 /2
1 1 1
1 1+z/2 1+2z/3+z 2 /6
1z 1z/2 1z/3
1 1+z/3 1+z/2+z 2 /12
1z+z 2 /2 12z/3+z 2 /6 1z/2+z 2 /12
Tabelle 1: Pade-Tafel f
ur j, k = 0, 1, 2
Fu
r den steifen Fall wichtig:
101
Bemerkung: Wahrend die Exponentialfunktion entlang der imaginaren
Achse fu r z ein anderes Verhalten aufweist als entlang der rellen Ach-
se, ist dies bei den Pade-Approximationen nicht der Fall. Man hat entweder
|R()| = 1 oder R() = 0, mit den jeweiligen Vor- und Nachteilen.
lim R(z) = 0.
z
w(x0 ) = y0
(5.7)
w (x0 + ci h) = f (x0 + ci h, w(x0 + ci h)) i = 1, . . . , s.
In den gegebenen Stu tzstellen 0 c0 < c1 < . . . < cs 1 erfu
llt w demnach
die Dierentialgleichung. Die Punkte x0 +ci h heien Kollokationspunkte. Mit
w(x0 + h) liefert das Verfahren eine Approximation der Losung y(x0 + h).
102
Satz 5.2 Das Kollokationsverfahren (5.7) ist aquivalent zu einem s-
stufigen Runge-Kutta-Verfahren mit Koeffizienten
ci 1
aij = lj (x) dx, bj = lj (x) dx, i, j = 1 : s,
0 0
mit den Lagrange-Polynomen lj (x) = l=j (x cl )/(cj cl ).
Beweis: Das Polynom w besitzt die Darstellung
s
w (x0 + th) = w (x0 + cj h)lj (t), t [0, 1].
j=1
Zugleich ist ci
1
w (x0 + th) dt = w(x0 + th)|c0i = (w(x0 + ci h) w(x0 ))/h
0 h
und damit
ci
s
s
w(x0 + ci h) = w(x0 ) + h w (x0 + cj h)lj (t) dt = y0 + h aij f (x0 + cj h, w(x0 + cj h)).
0 j=1 j=1
s
Wir setzen Ki := f (x0 + ci h, w(x0 + ci h)) = f (x0 + ci h, y0 + h j=1 aij Kj ) und erhalten so das
nichtlineare System (5.6) f
ur die Korrekturen. Weiter ist
1
s
s
w(x0 + h) = w(x0 ) + h w (x0 + cj h)lj (t) dt = y0 + h bj Kj = y1 .
0 j=1 j=1
103
Daraus folgt
1
s 1
s
1 k1
= t dt = ck1
j lj (t) dt = bj ck1
j . (5.8)
k 0 j=1 0 j=1
Speziell ist daher fu
r k = 1 die Konsistenzbedingung bj = 1 erfullt.
Analog zeigt man
cki s
= aij ck1
j , (5.9)
k j=1
was fu
r k = 1 die bekannte Bedingung aij = ci liefert: Jedes Kollokati-
onsverfahren ist konsistent und invariant unter Autonomisierung.
Die Beziehungen (5.8) und (5.9) charakterisieren die Koezienten des Kol-
lokationsverfahrens. Allgemein deniert man bei IRK-Verfahren die verein-
fachenden Annahmen
s
1
B(p) : bj ck1
j = , k = 1 : p; (5.10)
j=1
k
s
1 k
C(l) : aij ck1
j = ci , i = 1 : s, k = 1 : l; (5.11)
j=1
k
s
1
D(m) : bi ck1
i aij = bj (1 ckj ), j = 1 : s, k = 1 : m. (5.12)
i=1
k
Aus (5.8) folgt die Eigenschaft B(s), und aus (5.9) folgt C(s) fu
r die Kolloka-
tionsverfahren. Damit sind diese weitgehend charakterisiert, doch es bleibt
die Frage der Konsistenzordnung: Wie soll man die ci wahlen, damit das
Verfahren moglichst hohe Konsistenzordnung hat?
104
Satz 5.3 Ein durch Kollokation erzeugtes Runge-Kutta-Verfahren hat
r f C p (D) die Konsistenzordnung p genau dann, wenn die durch die
fu
Stutzstellen ci und Gewichte bi gegebene Quadraturformel die Ordnung p
besitzt.
Beweisskizze: Wir betrachten das Integral
x0 +h 1
s
g(x) dx = h g(x0 + th) dt = h bj g(x0 + cj h) + e(g)
x0 0 j=1
samt seiner Approximation durch die Quadraturformel mit Fehler e(g). Der Fehler e(g) ist von
der Ordnung p, falls
|e(g)| Chp+1 max |g (p) ()|.
(Vergleiche Gauquadratur: Dort war
1
s
a(t) dt = bj a(j ) + C D2s a()
1 j=1
1
1 1 s
1 1+ 1 + j 1 d2s 1+
g(x0 +th) dt = g(x0 +h ) d = bj g(x0 +h )+ C 2s g(x0 +h )|=
0 2 1 2 2 2 2 d 2
j=1
w (x) = f (x, w(x)) + (w (x) f (x, w(x))) = f (x, w(x)) + r(x, w(x)),
wobei r als Defekt oder Storung aufgefasst wird. Verglichen mit der urspr unglichen Dierential-
gleichung y = f (x, y(x)) folgt f
ur die Dierenz
x
y
w(x) y(x) = (x, s, w(s)) r(s, w(s)) ds
x0 y 0
(Sogenannte Formel von Alekseev/Grobner oder nichtlineare Variation der Konstanten). Ersetzt
man die Integration durch die Quadraturformel, so fallen alle Terme in der Summe heraus, da
r(x0 + ci h, w(x0 + ci h)) = w (x0 + ci h) f (x0 + ci h, w(x0 + ci h)) = 0 aufgrund der Kollokation.
mit
y
g(t) = (x, t, w(t)) (w (t) f (t, w(t))).
y0
Aufgepasst: Die Beschranktheit von g (p) () f
ur h 0 muss noch gezeigt werden!
105
Bemerkung: Satz 5.3 liefert fu r alle Kollokationsverfahren die zugehorige
Konsistenzordnung (Beispiele s.u.). Alternativ kann man zu einem gegebe-
nen IRK-Verfahren (Kollokationsverf. oder andere) die Konsistenzordnung
mit den in Satz 3.2 angegebenen Ordnungsbedingungen bestimmen. D.h.,
der Abgleich u
ber Taylorreihenentwicklungen lauft im impliziten Fall genau-
so ab wie im expliziten, nur die Summationen in den Bedingungsgleichungen
laufen immer von 1 bis s.
(i) Gau-Verfahren: Wahle die ci als Nullstellen des (auf [0, 1] verschobenen)
Legendre-Polynoms ds /dxs [xs (x 1)s ].
(ii) Radau-Verfahren: Wahle die ci als NST von ds1 /dxs1 [xs (x 1)s1 ]
(Radau IA) oder ds1 /dxs1 [xs1 (x 1)s ] (Radau IIA). Die letzte Klasse ist
durch 0 < c1 und cs = 1 charakterisiert und vorteilhafter:
106
Mit den Gau- und Radau IIA-Verfahren stehen damit zwei spezielle Ver-
fahrensklassen mit sehr guten Eigenschaften bereit.
Losung des nichtlinearen Systems. Nun kommen wir zur lange aufge-
schobenen Frage der Losung des nichtlinearen Gleichungssystems (5.6) in
den Korrekturen,
( )
s
Ki = f x0 + ci h, y0 + h aij Kj , i = 1 : s.
j=1
108
Auf den ersten Blick sieht diese Abschatzung wegen (I hJ)1 nicht
besser aus, aber es gilt nun:
1
ist EW von J :=
ist EW von (I hJ)1 . (5.16)
1 h
Reelle, stark negative Eigenwerte von J, die steifen Anteile, werden also
auf positive Eigenwerte 1 transformiert. Die Matrix (I hJ)1 kann
demnach als Filter interpretiert werden, der die steifen Anteile herausl-
tert. Es ist (I hJ)1 1 + c mit einer von der Steifheit unabhangigen
Konstanten c; falls J symmetrisch ist und nur reelle, negative Eigenwerte
vorliegen gilt c = 0. Also
1
D(K) < 1 falls h< , (5.17)
(1 + c)g/y
wobei g/y die von der Steifheit unabhangige Nichtlinearitat des Problems
beschreibt. Im Fall linearer Dierentialgleichungen konvergiert das verein-
fachte Newtonverfahren in einem Schritt.
109
man nach dem Banachschen Fixpunktsatz
. Z k
Z k+1
Z
Z k , = k := .
1 Z k1
Dabei ist k eine Approximation an die Konvergenzrate .
c1 a11 0 0 c1 a11 0 0
.. .. . . . .. .. . . .
. . 0 . . 0
.
cs as1 . . . ass cs as1 . . . a11
DIRK b1 . . . bs SDIRK b1 . . . bs
110
Numerische Software. Die Codes RADAU5 und RADAU von Hairer und
Wanner sind die zur Zeit wohl besten IRK-Verfahren (ftp-Adresse Kap. 3,
siehe auch www-m3.ma.tum.de/m3/software.html). RADAU5 basiert auf
dem L-stabilen Kollokationsverfahren Radau IIA zu s = 3 Stufen und hat
die Ordnung 2s 1 = 5 ( y1 Ordnung 3). RADAU enthalt zusatzlich die
Varianten mit s = 5 und s = 7 Stufen und ist als Verfahren mit variabler
Ordnung (5, 7, 13) und variabler Schrittweite implementiert. Bei diesen IRK-
Verfahren ist die Matrix Imsms hAJ nicht von der Blockgestalt (5.18); es
existiert aber eine Transformation, die den Aufwand fur eine LR-Zerlegung
z.B. im Fall s = 3 von (3m)3 /3 auf 5m3 /3 Operationen reduziert.
Linear-implizite Einschrittverfahren
Da die Implementierung der IRK-Verfahren sehr aufwandig ist und bis En-
de der 80er Jahre nicht klar war, ob solche Verfahren uberhaupt ezient
realisiert werden konnen, hat man bereits fru
h nach Alternativen gesucht.
Die Grundidee der linear-impliziten ESV wollen wir am Beispiel des linear-
impliziten Eulerverfahrens studieren:
Wir beschranken uns auf den autonomen Fall y = f (y). Der implizite Euler
lautet hier y1 = y0 + hK, K = f (y0 + hK). Wir fu hren nun eine einzige
r 0 = K f (y0 + hK) aus,
Newtoniteration mit Startwert K 0 = 0 fu
K 1 = K 0 (I hDf (y0 ))1 (K 0 f (y0 + hK 0 )) = (I hDf (y0 ))1 f (y0 ),
und setzen als neue Approximation
y1 = y0 + hK 1 = y0 + h(I hDf (y0 ))1 f (y0 ). (5.19)
Dieses Verfahren ist der linear-implizite Euler. Die Multiplikation mit der
Matrix (I hDf (y0 ))1 erfullt hier wie bei der Newton-Iteration der IRK-
Verfahren eine Filterfunktion. Oensichtlich hat der linear-implizite Euler
diesselbe Stabilitatsfunktion wie der implizite Euler und ist demnach L-
stabil. Daru
berhinaus hat dieses Verfahren die Konsistenzordnung 1.
111
ur lineare Probleme wie die Testfunktion y = y stimmt der semi-implizite Euler mit dem
F
berein. Daher ist R(z) = 1/(1 z). Die Ordnung folgt (i) f
impliziten u ur hDf (y0 ) < 1 aus der
Neumannschen Reihe
(I hDf (y0 ))1 = I + hDf (y0 ) + (hDf (y0 ))2 + (hDf (y0 ))3 + . . . ,
y1 (0) + O(h2 )
y1 (h) = y1 (0) + h
mit y1 (0) = y0 und y1 (0) = (Df (y0 )y0 + Df (y0 ) y1 (h) + f (y0 ))|h=0 = f (y0 ).
y1 (h) + hDf (y0 )
ur C folgt Re (h) 0
= 1 h EW von I hDf (y0 ). F
Sei EW von Df (y0 ), dann ist
1. Vergleiche auch die Aquivalenz
und somit Re (5.16).
Der linear-implizite Euler ist damit eine deutlich weniger aufwandige Alter-
native zum impliziten Euler. In jedem Integrationsschritt ist nur ein lineares
Gleichungssystem
fu
r die neue Approximation y1 (jetzt ohne ) zu losen.
112
Verlust an Stabilitat rechnen: Ab der Spalte 3 des Extrapolationstableaus
liegt keine A-Stabilitat mehr vor, da Teile von C in der Nahe der ima-
ginaren Achse nicht mehr zum Stabilitatsgebiet gehoren. Man spricht auch
.
von A()-Stabilitat mit = 89, 7o , vergl. Abb. 35.
Numerische Software fu r diesen Zugang: Codes LIMEX (Deuhard/Nowak)
sowie SEULEX (Hairer/Wanner).
S
C
Abbildung 35: A()-Stabilitat: Der Sektor {z : |arg (z)| < } ist in S enthalten.
113
(f (y) Jy) explizit. Wir denieren so das ROW-Verfahren
s
y1 = y0 + h bj Kj , (5.20)
j=1
( )
i1
i1
Ki = hJii Ki + f y0 + h ij Kj + hJ ij Kj . (5.21)
j=1 j=1
Die Vorschrift fu
r die Korrekturen Ki fu
hrt auf das lineare Gleichungssystem
( )
i1 i1
(I hii J)Ki = f y0 + h ij Kj + hJ ij Kj , (5.22)
j=1 j=1
das in jeder Stufe zu losen ist. Die Koezienten bi , ij , ij legen die Methode
fest. (Eine andere Herleitung geht von einem DIRK-Verfahren aus, auf das
1 Newtoniteration angewandt wird.)
114
i1 i
mit Koezienten i = j=1 ij sowie i = j=1 ij .
fu
r p = 2 zusatzlich
bj (jk + jk ) = 1/2, (5.25)
j,k
und fu
r p = 3 zusatzlich
bj jk jl = 1/3, bj (jk + jk )(kl + kl ) = 1/6. (5.26)
j,k,l j,k,l
115
Algorithmus 5.1 ROW-Verfahren
fu
r i = 1 : s
( )
i1
i1
Lose (I hJ)Ki = f y0 + h ij Kj + hJ ij Kj ;
j=1 j=1
s
y1 = y0 + h j=1 bj Kj ;
Mit der FSAL-Strategie (siehe Abschnitt 3.5) wird eine weitere Stufe aus-
gewertet,
(I hJ)K3 = f (x1 , y1 ) (6 + 2)(K2 f (x0 + 21 h, y1 + h21 K1 ))
2(K1 f (x0 , y0 )) + hDx f (x0 , y0 ).
116
Als Fehlerschatzer verwendet man
h
y1 y1 = (K1 2K2 + K3 )
6
Dabei ist y1 ein L-stabiles Verfahren der Ordnung 2 und y1 ein A-stabiles
Verfahren der Ordnung 3.
Bemerkungen
Beim Verfahren ode23s wird keine lokale Extrapolation durchgefu hrt
und mit der Approximation 2.ter Ordnung weitergerechnet. Die bessere
Stabilitat wird wichtiger als die hohere Ordnung eingeschatzt.
5.3 Mehrschrittverfahren fu
r steife Differentialgleichungen
117
Das MSV j yi+j = h j fi+j , angewandt auf y = y, ergibt die homo-
gen Dierenzengleichung
k
(j zj )yi+j = 0 mit z = h. (5.28)
j=0
Mit den Techniken aus Abschnitt 4.3 (Satz 4.2) kann die homogene Rekur-
sion (5.28) abhangig vom komplexen Parameter z analysiert werden. Das
zu (5.28) gehorige charakteristische Polynom ist
() z() , () = l
l , () = l l .
Eine Losung {yi } von (5.28) ist Nullfolge genau dann, wenn fur jede Null-
stelle von () z() gilt || < 1, und zwar unabhangig von z!
Bei der AStabilitat nehmen wir den Fall || = 1 mit und fordern nur, dass
Losungen der Testgleichung nicht wachsen. Es gilt damit
Satz 5.4 Das MSV j yi+j = h j fi+j ist A-stabil
Jede NST von () z() erfu r alle z C .
llt || 1 fu
Falls mehrfache NST ist, gilt || < 1.
Die Abhangigkeit von z wird wie bei ESV durch den Begri des Stabi-
litatsgebietes genauer beschrieben:
Beispiele fu
r Stabilitatsgebiete sind in Abb. 36 dargestellt.
118
S
C 3i C
111111
000000
C
5
11111
00000
S
000000
111111
000000
111111
00000
11111
i
00000
11111 000000
111111
00000
11111 000000
111111
000000
111111
00000
11111 000000
111111
5
00000 i
11111
1
000000
111111
000000
111111 S
0000003i
111111
Adams-Bashforth k = 2 Adams-Moulton k = 2 BDF k = 2
Abbildung 36: Stabilitatsgebiete von MSV
Was kann man mit MSV bezu glich Stabilitat und Ordnung maximal errei-
chen? Die Antwort (ohne Beweis):
Trotz dieser groen Einschrankung gibt es eine Verfahrensklasse, die sehr er-
folgreich bei steifen Problemen eingesetzt wird, namlich die BDFVerfahren
(bis k = 6). Die negative reelle Halbachse ist bei den BDF-Verfahren immer
in S enhalten, daher hat man dort gute Stabilitatseigenschaften. Ab k = 3
sind die Verfahren aber nur noch A()-stabil und daher nicht geeignet fu r
oszillatorische Probleme.
119
k=6
k=5
k=3 k=4
Beginnend mit Gear (1971) waren die BDF-Verfahren die ersten erfolgrei-
chen Methoden fu r steife Dierentialgleichungen. Auch heute noch gelten
sie als Standard fu r diese Problemklasse. Da die Verfahrenskoezienten fest
sind, unterscheiden sich die eingesetzten Verfahren nur in Bezug auf die, wie
bei MSV u blich, sehr komplexe Implementierung mit variabler Ordnung und
variabler Schrittweite. Basis ist der Ansatz (vergl. Abschnitt 4.1)
120
Dazu noch ein paar Details: Die Verfahren werden als Pradiktor-Korrektor-Schema implemen-
tiert. Das Korrektorpolynom ist obiges Polynom q, das wir im folgenden mit qC bezeichnen. Das
Pradiktorpolynom qP interpoliert (xnk , ynk ) bis (xn , yn ) und liefert mit Extrapolation auf xn+1
den Pradiktor- oder Startwert
0
yn+1 := qP (xn+1 ).
Die Dierenz der beiden Polynome interpoliert (xnk+1 , 0) bis (xn , 0) sowie (xn+1 , yn+1 yn+1
0 ).
mit
1 1 1
n+1 := + + ... + .
xn+1 xn xn+1 xn1 xn+1 xnk+1
Der Faktor n+1 (nicht mit den j des allgemeinen MSV zu verwechseln!) wird f uhrender Ko-
effizient (leading coefficient) genannt. Nun kann die Gleichung (5.29) neu formuliert werden.
0
Mithilfe des Pradiktors yn+1 0
und seiner Ableitung yn+1 := qP (xn+1 ) ist (5.29) aquivalent mit
dem nichtlinearen System
0
0 = yn+1 + n+1 (yn+1 yn+1
0
) f (xn+1 , yn+1 ) =: F (yn+1 ).
Man beachte, dass hier ganz analog zu den IRK- und ROW-Verfahren die Matrix n+1 I Dy f =
1
n+1 (I n+1 Dy f ) auftaucht und das Konvergenzverhalten im Vergleich zur Fixpunktiteration
im steifen Fall entscheidend verbessert. Als Schatzer f
ur den Fehler
e = y(xn+1 ) yn+1
verwendet man die Dierenz zwischen Pradiktor und Korrektor in den Formen
. 1 1
(i) e = (y 0 yn+1 )
n+1 xn+1 xnk n+1
. 1
(ii) e = (n+1 I Dy f (xn+1 , yn+1
0
))1 (y 0 yn+1 ).
xn+1 xnk n+1
Die zweite Schatzung (ii) erfordert eine weitere Vorwarts-R uckwarts-Substitution, ist aber im
steifen Fall besser und liefert groere Schrittweiten. In den meisten Codes wird aber die billigere
Schatzung (i) eingesetzt.
121
BDF-2-Verfahren mit SWS
Gegeben: yn2 , yn1 , yn zu Zeitpunkten xn2 , xn1 , xn
(zur
uckliegende Schritte oder Anlaufrechnung, impl. Euler).
Pr
adiktor:
0
yn+1 = yn + (xn+1 xn )[yn , yn1 ] + (xn+1 xn )(xn+1 xn1 )[yn , yn1 , yn2 ];
0
yn+1 = [yn , yn1 ] + (2xn+1 xn xn1 )[yn , yn1 , yn2 ];
Korrektor - l
ose nach yn+1 :
0
0 = yn+1 n+1 yn+1
0 + n+1 yn+1 f (yn+1 , xn+1 );
Fehlersch
atzer und neue Schrittweite:
. 1 1
(0)
ERR =
yn+1 yn+1
;
n+1 xn+1 xn2
hnew = hn min(u , max(l , s ERR1/3 ));
1 1
Dabei sind n+1 = + der leading coecient,
xn+1 xn1 xn+1 xn
( )2
1 m yi
y2 := m die Fehlernorm mit Gewichtung w,
i=1
ATOL+RTOL |wi |
u 1 + 2 die obere Schranke f
ur die Anderung h,
l = 0.2 die untere Schranke f
ur die Anderung h und s = 0.9 der Sicherheitsfaktor.
122
Van der Pol, = 1.E6
2.5 Van der Pol, = 1.E6
2
10
2 y(2)
1.5
0.5
ode23s
CPU
1
0 10
y(1)
0.5
ode15s
1
1.5
0
2.5 10
1 2 3 4 5 6 7 8
0 2 4 6 8 10 12 10 10 10 10 10 10 10 10
x Fehler
Wir testen die in MATLAB verfu gbaren steifen Integratoren anhand des
Van-der-Pol-Oszillators in den Koordinaten von Lienhart (2.3)
1
y1 = y2 , y2 = (y1 y23 /3 + y2 )
mit Anfangswert y1 (0) = 2, y2 (0) = 0. Die Anwendung des ROW-Verfahrens
ode23s aus MATLAB liefert dann zum Parameterwert = 106 die in
Abb. 38 dargestellte numerische Losung im Intervall [0, 11]. Ein explizites
Verfahren ist hier chancenlos, da Schrittweiten in der Groenordnung von
notwendig sind.
Auf der rechten Seite von Abb. 38 sieht man ein Aufwandsdiagramm, in
dem die erreichte Genauigkeit gegen die Rechenzeit fu r verschiedene Tole-
ranzvorgaben aufgetragen ist (vergl. Abschnitt 4.5, Abb. 27). Der Fehler am
Endpunkt xend = 11 wurde dabei durch Vergleich mit einer mit hoher Ge-
nauigkeit berechneten Referenzlosung ermittelt. Wahrend ode23s bei nied-
rigen Toleranzvorgaben ezienter und auch etwas genauer ist, schneidet das
BDF-Verfahren ode15s bei hoheren Genauigkeiten klar besser ab.
123
Exkurs: Lininenmethode Als zweites Beispiel betrachten wir das N -dimensionale System
y = Ay + d(y) (5.30)
wobei
2 1
.
1 2 2 . .
A=N
.. ..
. . 1
1 2
und
D(1 + y1 ) exp(/y1 ) + N 2 y2
D(1 + y2 ) exp(/y2 )
d(y) =
..
. .
D(1 + yN 1 ) exp(/yN 1 )
D(1 + yN ) exp(/yN ) + N 2 1
Dieses Dierentialgleichungssystem entsteht durch Anwendung der Linienmethode auf die Reaktions-
Diusionsgleichung
Die Variable T (x, t) beschreibt dabei den Temperaturverlauf einer chemischen Reaktion im ein-
dimensionalen Kontinuum 0 x 1. ist die Warmefreisetzung, die Aktivierungsenergie und
D = R exp()/() mit Reaktionsrate R.
Bei der Linienmethode wird zunachst aquidistant in x diskretisiert und dann auf dem Gitter
0 = x1 < x2 < . . . xN < 1 die Ableitung Txx durch einen Dierenzenoperator approximiert,
. 1
Txx (xj , t) = (T (xj1 , t) 2T (xj , t) + T (xj+1 , t)), j = 1 : N, x = 1/N.
x2
Die Variable yj approximiert nun T (xj , t) und gen
ugt dem System (5.30), wobei die Randwerte
u
ber
. T (x, t) T (x, t)
T (xN +1 , t) = T (1, t) = 1 und 0 = Tx (x0 , t) = T (x, t) = y2
2x
eingearbeitet werden.
Abb. 39 zeigt die numerische Losung der im Ort diskretisierten Gleichung (5.30) zu N = 100
Gitterpunkten. Zahlwerte der Konstanten sind R = 5, = 20, = 1, der Anfangswert ist yj,0 = 1
f
ur alle Komponenten. Man sieht, wie sich nach einer gewissen Anlaufzeit von links nach rechts
eine Front ausbreitet.
124
2
Temperature T
1.5
1
0.4
0.3 1
0.8
0.2
0.6
0.1 0.4
0.2
0 0
Time t
x
Das Dierentialgleichungssystem (5.30) ist steif, da die Eigenwerte von A negativ sind und mit N 2
wachsen. Im Fall N = 100 liegt der betragsgrote Eigenwert bereits bei 3.9 104 , der kleinste bei
9.6. Die beiden Verfahren ode23s und ode15s arbeiten hier wiederum zuverlassig und berechnen
die Losung bei Standardeinstellung der Toleranzen (RelTol = 103 ) in 124 bzw. 129 Schritten
und 7.53 bzw. 2.62 Sekunden Rechenzeit (Pentium III 1 GHz). Falls N noch groer gewahlt wird,
kommt der linearen Algebra (Losung der linearen Gleichungssysteme im ROW-Verfahren bzw. in
der Newton-Iteration bei BDF) entscheidende Bedeutung zu. Ohne eine eziente lineare Algebra
kann auch der beste steife Integrator dann nichts mehr ausrichten (Option JPattern in odeset).
Der Begri der A-Stabilitat war bisher das zentrale Kriterium zur Analyse
und Einteilung von numerischen Verfahren fu r steife Dierentialgleichun-
gen. Da man bei der A-Stabilitat nur eine lineare autonome Testgleichung
125
betrachtet, stellt sich die Frage, wie dieses Konzept auf den nichtlinearen
Fall verallgemeinert werden kann. Wir beschranken uns in diesem Abschnitt
auf den autonomen Fall y = f (y).
Bei der linearen Theorie war unsere Modellvorstellung, dass alle Eigenwerte
negative Realteile haben und daher die Losungen zu einer stabilen Trajek-
torie hin streben. Eine nichtlineare Verallgemeinerung ist
mit 0 erfu
llt ist.
126
Eine naturliche Forderung an numerische Verfahren ist nun, dass sie die
Dissipativitat an die diskrete Approximation vererben. Auf Butcher (1975)
geht die folgende Denition zuru ck.
Aus der B-Stabilitat eines Verfahrens folgt die A-Stabilitat: Wir betrachten
y = y mit R, 0. Die Ungleichung (5.32) ist dann erfu llt, die
Gleichung ist dissipativ. Falls das Verfahren B-stabil ist, gilt
|y1 z1 | = |R(h)(y0 z0 )| |y0 z0 |
und damit |R(h)| 1 fur alle Schrittweiten h. Im komplexen Fall =
+ i mit 0 verwendet man die reellwertige Dierentialgleichung
( )
u = Au, A := ,
die wegen uT Au = (u21 + u22 ) 0 ebenfalls dissipativ ist.
Die Umkehrung dieser Aussage gilt nicht: So sind z.B. alle linear-impliziten
Einschrittverfahren nicht B-stabil. Nur in der Klasse der impliziten Runge-
Kutta-Verfahren ndet man B-stabile Vertreter.
y z, y z. Dann ist
m (x) = 2w w
, w w
= 2f (w(x0 + ci h)) f (w(x
0 + ci h)), w w
0.
Integration liefert
x0 +h
s
m(x0 + h) = m(x0 ) + m (x) dx = m(x0 ) + h bi m (x0 + ci h) m(x0 ).
x0 i=1
127
Dabei geht ein, dass m ein Polynom vom Grad 2s1 ist und durch die Quadratur exakt integriert
wird. Dar
uberhinaus sind die Gewichte bi positiv. Also folgt
y1 z1 2 = m(x0 + h) m(x0 ) = y0 z0 2
(i) bi 0, i = 1 : s;
(ii) M := (bi aij + bj aji bi bj )i,j=1:s ist positiv semidefinit
erfu
llen, ist das Verfahren B-stabil. Insbesondere sind RadauIIA-Verfahren
B-stabil.
128
legen. Trotzdem konnen wir die in Abschnitt 5.1 gegebene Charakterisierung
noch etwas prazisieren.
Betrachten wir zum Test den Fall y = Ay mit A symmetrisch und Eigenwerten m m1
. . . 1 . Im Fall der 2-Norm ist dann (A) = 1 und A = (max (AT A))1/2 = maxi |i |. Die
Bedingung (5.35) bedeutet dann, dass
1 (xend x0 ) max |i |,
i
y0 y(xend )
auf, die eine Kondition (x0 , xend ) besitzt.Im steifen Fall ist dieses Problem
129
sehr gut konditioniert. Die numerische Approximation ist eine weitere Ab-
bildung
y0 yn , y(xend ) yn Chp ,
die ebenfalls eine Kondition h (x0 , xend ) besitzt, die von der Schrittweite h
und dem Verfahren abhangt. Damit ist
wenn der Anfangswert mit y0 gestort wird. Man fordert nun sinnvollerweise
h (x0 , xend ) (x0 , xend ), was im steifen Fall bei expliziten Verfahren aber
nur mit sehr kleinen Schrittweiten h zu erreichen ist. Mit anderen Worten:
Das Anfangswertproblem ist steif, wenn die diskrete Kondition h (x0 , xend )
bei expliziten Verfahren nur fu r sehr kleine Schrittweiten h mit (x0 , xend )
u
bereinstimmt.
130
Kapitel 6
Stochastische Differentialgleichungen
6
Stochastische Dierentialgleichungen spielen eine groe Rolle in der Finanz-
mathematik, speziell bei der Modellierung von Aktienkursen. In diesem Ka-
pitel wollen wir die Grundlagen dazu bereitstellen und einfache numerische
Integrationsverfahren kennenlernen.
Weiterfu
hrende Literatur: Seydel, R.: Einfu
hrung in die numerische
Berechnung von Finanzderivaten. Springer, 2000
Gu
nther, M., Ju
ngel, A.: Finanzderivate mit MATLAB. Vieweg 2003
131
Dabei ist X die gesuchte (skalarwertige) Variable, die von der Zeit t abhangt,
die Funktion a ist der sogenannte Driftterm und b(X, t) t der Diffusions-
term mit einem verallgemeinerten stochastischen Prozess t (der Zufall).
Man hat also eine Uberlagerung von deterministischen und stochastischen
Eekten.
Wir betrachten ein unendlich langes Rohr, das mit Wasser oder Gas gef ullt ist. In der Mitte bei
x = 0 bringen wir einen Tropfen Farbe mit (normierter) Masse 1 ein und verfolgen, wie sich die
Farbmolekule nach links und rechts ausbreiten. Dieser Diusionsprozess kann als Random Walk
approximiert werden Wir fuhren dazu ein Gitter
xi = ix ur i = 0, 1, 2, . . . ;
f tj = jt f
ur j = 0, 1, 2, . . .
132
ein. Im Punkt (i, j) = (0, 0), also bei x = 0 und zum Zeitpunkt t = 0, startet nun ein Partikel
seine Bewegung nach folgendem Gesetz: Entweder es bewegt sich nach links zur Position i 1 mit
Wahrscheinlichkeit 1/2, oder es bewegt sich rechts zur Position i + 1 mit gleicher Wahrscheinlich-
keit 1/2.
Sei nun p(i, j) die Wahrscheinlichkeit, dass sich das Partikel zum Zeitpunkt tj = jt im Punkt
xi = ix bendet. Dann gelten
{
0 fur i = 0 1 1
p(i, 0) = und p(i, j + 1) = p(i 1, j) + p(i + 1, j) .
1 fur i = 0 2 2
In der Physik ergibt sich die Konstante aus D = RT /(NA F ) mit Gaskonstante R, Temperatur
T , Reibungskoezient F und Avogadro-Zahl NA .
In der Stochastik setzt man D = 1 und motiviert so den Wiener-Prozess. Die prazise Durchf
uhrung
des Grenzubergangs basiert im u
brigen auf dem Satz von Laplace - DeMoivre.
Fu
r eine numerische Realisierung des Wiener-Prozesses benotigen wir ein
diskretes Modell. Dazu sei h = t ein konstantes Zeitinkrement, mit dem
wir das Gitter tj = jh fu
r j = 0 : n festlegen. Dann ist
n
n
Wtn = (Wtj Wtj1 ) = Wj
j=1 j=1
133
Algorithmus 6.1 Diskreter Wiener-Prozess
t0 = 0; W0 = 0; h = T /n;
Fur j = 0 : n 1
tj+1 = tj + h;
Ziehe Zufallszahl Z N (0, 1);
Wj+1 = Wj + Z h;
In MATLAB stehen mit dem Befehl randn Zufallszahlen zur Verfu gung,
die N (0, 1)-normalverteilt sind. Diese werden mit einem speziellen Zufalls-
zahlengenerator erzeugt, und sukzessive Aufrufe von randn liefern eine Fol-
ge von entsprechenden Zufallszahlen. Eine Reproduktion dieser Folge ist
moglich, wenn man mit der Option state den Generator zuru cksetzt. Da-
mit lautet der Algorithmus 6.1 in MATLAB
randn(state,1);
T = 10; h = 0.01; n = T/h; W(1) = 0;
t = 0:h:T;
for i=1:n, W(i+1) = W(i)+randn*sqrt(h); end
134
h=0.01
3
1
Wt
6
0 1 2 3 4 5 6 7 8 9 10
t
(Fur Fragen zur Wohldeniertheit verweisen wir auf die entsprechende Li-
teratur). Mit diesem Integralbegri gilt
t
1 t
Ws dWs = Wt2 ,
0 2 2
im Unterschied zur klassischen Analysis, wo wir nur den ersten Term erwar-
ten wu
rden.
Bemerkung: Es existiert jedoch ein alternativer Integralbegri, das Stratonovich-Integral, bei
dem der obige Extraterm nicht auftritt. Es ist deniert als
t
n1
Xs dWs := lim X1/2(tk +tk+1 ) (Wtk+1 Wtk ). (6.3)
0 n
k=0
Mit dem Begri des Ito-Integrals konnen wir nun eine stochastische Diffe-
rentialgleichung nach Ito einfu
hren. Dazu schreiben wir (6.1) in Integralform
t t
Xt = X0 + a(Xs , s)ds + b(Xs , s)dWs (6.4)
0 0
135
mit dem gesuchten stochastischen Prozess Xt . Das erste Integral in (6.4) ist
ein Riemann- oder Lebesgue-Integral, das zweite jedoch ein Ito-Integral. Als
Kurzschreibweise verwendet man fur die stochastische Dierentialgleichung
(6.4) oft auch
dXt = a(Xt , t)dt + b(Xt , t)dWt (6.5)
mit Drift a und Diusion b dWt .
Fu
r verschwindenden Diusionsterm ist dies nichts anderes als der explizite
Euler. Die Implementierung des Verfahrens verknupft Euler-Zeitschritt mit
dem Wiener-Prozess:
Wir wollen dieses Verfahren nun anhand der speziellen stochastischen Dif-
ferentialgleichung (sog. geometrische Brownsche Bewegung)
136
Brownsche Bew., r=0.1, =0.4
14
12
10
8
S
0
0 1 2 3 4 5 6 7 8 9 10
t
studieren. Dabei ist r 0 z.B. ein konstanter risikofreier Zinssatz und > 0
die konstante Volatilitat. In MATLAB lautet die Simulation mit dem Euler-
Maruyama-Verfahren, mit gleichzeitiger Berechnung des Wiener-Prozesses,
In der Abbildung 41 sind zwei Pfade dargestellt, die mit diesem Algorithmus
berechnet wurden (r = 0.1, = 0.4, S0 = 1). Andere Realisierungen des
Wiener-Prozess fu
hren wiederum auf andere Pfade, und es ist oensichtlich,
dass es wenig Sinn macht, die Pfade punktweise zu interpretieren oder zu
bewerten. Stattdessen wendet man die Monte-Carlo-Simulation an:
137
1) Lose die stochastische Dierentialgleichung M mal, um Approximationen
ST,1 , . . . , ST,M zu erhalten.
1
M
.
E(ST ) = M = ST,i
M i=1
u
ber den Mittelwert M .
Im Fall der Brownschen Bewegung (6.7) ist die exakte Losung bekannt. Sie
lautet
St = S0 exp((r 2 /2)t + Wt )
und hangt naturlich von der Realisierung des Wiener-Prozess Wt ab. Zur
Bestimmung der Konvergenzordnung berechnet man nun M 1 Realisie-
rungen des Wiener-Prozess. Zu jeder Realisierung i berechnet man weiter N
h
numerische Approximationen ST,ij mit Schrittweiten hj = 21j h. Der Fehler
zur Schrittweite hj ergibt sich geschatzt zu
1
M
h
e(hj ) = |ST,i ST,ij |.
M i=1
138
Es ist nun
log e(hj ) log C + log hj ,
und zur Bestimmung der Parameter log C sowie ziehen wir die lineare
Ausgleichsrechnung heran. Das Ergebnis ist 0.5, was u
ber eine theore-
tische Analyse veriziert werden kann: Die starke Konvergenzordnung des
Euler-Maruyama-Verfahrens ist = 1/2.
Der Begri der schwachen Konvergenz basiert dagegen auf dem Fehler
Bezu
glich dieses Fehlermaes konvergiert das Euler-Maruyama-Verfahren
mit der Ordnung 1.
Im deterministischen Fall b = 0 ergibt sich diese Formel direkt aus der partiellen Integration
t t
f (X(s))|t0 = f (X(s))X (s) ds + 0 f (X(s)) ds
0 0
und X (s) = a(X(s)). Mit Dierentiation nach t folgt daraus die Kettenregel
d d d
f (X(t)) = f (X(t)) X(t).
dt dX dt
Die Ito-Formel ist also die Verallgemeinerung der Kettenregel auf den stochastischen Kalk
ul.
139
Mit der Ito-Formel als verallgemeinerter Kettenregel kann man nun wie bei
den Einschrittverfahren im deterministischen Fall die Taylorentwicklung der
exakten Losung aufstellen. Wir setzen f (X) = X und erhalten
t t
Xt = X0 + a(Xs )ds + b(Xs )dWs .
0 0
Im nachsten Schritt werden a und b mit der Ito-Formel substituiert,
t( s s )
1 2
Xt = X0 + a(X0 ) + (a a + a b )dz + a b dWz ds
2
t( 0
s 0
s 0
)
1
+ b(X0 ) + (b a + b b2 )dz + b b dWz dWs .
0 0 2 0
Beispiele:
1) Das Milstein-Verfahren basiert auf der Darstellung
t t s
Xt = X0 + a(X0 )t + b(X0 ) dWs + b (X0 )b(X0 ) dWz dWs + R
0 0 0
mit Termen hoherer Ordnung, die in R zusammegefasst sind. Wichtig: Das
Doppelintegral kann man explizit berechnen. In diskreter Form lautet das
Verfahren dann (jetzt nichtautonom angeschrieben)
1
X1 = X0 + a(X0 , t0 )t + b(X0 , t0 )W + b (X0 , t0 )b(X0 , t0 )((W )2 h).
2
Das Milstein-Verfahren ist von der starken Konvergenzordnung 1.
Bemerkungen:
Alles hier Dargestellte bezieht sich auf eine skalare stochastische Dif-
ferentialgleichung. Der mehrdimensionale Fall lauft ganz analog, falls
man es in allen Komponenten mit demselben Wiener-Prozess zu tun
hat. Ist dies nicht der Fall, treten zusatzliche gemischte Ableitungen
und Integrale auf, so dass sich die Verfahren hoherer Ordnung deutlich
komplizierter gestalten. Nur das Euler-Maruyama-Verfahren ist auch
in diesem Fall einfach zu handhaben.
141
Kapitel 7
Randwertprobleme
charakterisiert, mit unbekannter Periode T und Phasenxierung der Komponente zj . Wir trans-
formieren das System auf die neue Zeit x = Tt , so dass x [0, 1]. Das Ergebnis ist
} }
z = T g(z) z(0) = z(0 + 1)
DGL Randbed. (7.2)
T = 0 zj (0) =
Mit y(x) := (z T (x), T (x))T und a = 0, b = 1 liegt damit ein Randwertproblem (7.1) vor. Insge-
samt hat man m + 1 Dierentialgleichungen und m + 1 Randbedingungen.
142
7.1 EinfachSchieen
Die einfachste Technik zur Losung eines RWP ist das EinfachSchieen:
Betrachte dazu die skalare Dierentialgleichung zweiter Ordnung
u
u(x;s)
Steigung s
^
u(x;s)
x
a b
Abbildung 42: Einfachschieen
F (s) := u(b; s)
143
abhangig von s Rm und sucht den Anfangswert s mit r( s, y(b; s)) = 0. Es
ist also eine Losung des nichtlinearen Gleichungssystems
wobei
F (s1 , . . . , sk + k , . . . , sm ) F (s1 , . . . , sm )
ak = .
k
Notwendig dazu ist die Auswertung von F (s1 , . . . , sk + k , . . . , sm ) fu r k =
1, . . . , m, d.h. weitere m AWP der Form (7.3). Hinzu kommt schlielich der
Aufwand zur Losung des linearen Gleichungssystems fu r si . Der Haupt-
aufwand liegt damit in der Berechnung von DF (si ).
144
Dabei sind y1 , y2 die Komponenten des dynamischen Systems und y3 = T die gesuchte Periode.
Die erste Komponente wird am linken Rand xiert. Beim Einfach-Schieen ist nun das nichtlineare
Gleichungssystem
s1 y1 (1; s)
0 = F (s) = s2 y2 (1; s)
s1
zu losen, wobei die Variable s3 indirekt u
ber y1 (1; s) und y2 (1; s) eingeht. Als Jacobimatrix
erhalten wir
1 0 0 y1 /s1 y1 /s2 y1 /s3
DF (s) = 0 1 0 y2 /s1 y2 /s2 y2 /s3 .
1 0 0 0 0 0
Der zweite Term besteht aus Ableitungen der Losung nach den Anfangswerten und kann statt
u
ber numerische Dierentiation durch die Integration der Sensititivitatengleichung (Satz 2.5)
bestimmt werden:
f
(x) = (x, y(x; s)) (x), (x0 ) = I33 , (7.7)
y
mit
y1 /s1 y1 /s2 y1 /s3
(x) = y2 /s1 y2 /s2 y2 /s3 .
y3 /s1 y3 /s2 y3 /s3
Man lost nun simultan die Dierentialgleichung (7.6) und die Sensitivitatendierentialgleichung
(7.7) und erhalt so am rechten Rand x = 1 die Auswertung von F (s) sowie von DF (s).
Dieses Vorgehen erfordert seinerseits die Berechnung der Jacobimatrix f y der rechten Seite ent-
lang der Losung y(x; s). Die numerische Dierentiation hierf
ur ist aber langst nicht so empndlich.
Im Falle eines steifen Systems mit impliziter oder semi-impliziter Integration wird diese Informa-
tion sowieso benotigt.
4 4 4
y2
y2
y2
2 2 2
0 0 0
0 2 4 0 2 4 0 2 4
y y y
1 1 1
Die Abbildung 43 zeigt die Trajektorien zu den ersten 3 Iterationen beim Einfachschieen mit
dem Startwert s0 = (1, 2, 6) und = 1 sowie A = 1, B = 3. Das Verfahren konvergiert rasch
gegen die Periode T = 7.1569.
145
Im allgemeinen Fall F (s) = r(s, y(b; s)) hat man als Jacobimatrix
r(s, y(b; s)) r(s, y(b; s)) y(b; s)
DF (s) = + ,
s y s
so dass auch hier das in Beispiel 7.2 skizzierte Vorgehen u
ber die Sensiti-
vitaten zur Berechnung der Ableitungsinformation y(b; s)/s moglich ist.
Fu
r groe Lipschitzkonstanten L und lange Intervalle [a, b] kann eine kleine
Dierenz s1 s2 extrem verstarkt werden.
Abhilfe: Unterteile das Intervall [a, b] Mehrzielmethode, s.u.
146
7.2 Mehrzielmethode (MehrfachSchieen, Multiple Shooting)
Der Ansatz ergibt sich aus den Schwierigkeiten beim Einfach-Schieen: Wir
unterteilen das Intervall [a, b] durch ein Gitter
mit k Stu
tzpunkten oder Mehrzielknoten xj .
(x3 ,s3 )
a=x1 x2 x3 x x =b x
k1 k
Abbildung 44: Mehrzielmethode
In jedem Teilintervall [xj , xj+1 ] bezeichne y(x; xj , sj ) die Losung des An-
fangswertproblems y = f (x, y), y(xj ) = sj . Von den durch Einfach-Schieen
gewonnenen Losungsabschnitten fordern wir, dass die daraus stu ckweise zu-
sammengesetzte Funktion stetig ist.
147
r x [xj , xj+1 ] und j = 1 : k 1 stetig ist und die
Y (x) = y(x; xj , sj ) fu
Randbedingungen erfu llt. Also haben wir die Bedingungen
y(xj+1 ; xj , sj ) = sj+1 r j = 1, . . . , k 2,
fu
r(s1 , y(xk ; xk1 , sk1 ) = 0.
Dieses nichtlineare Gleichungssystem der Dimension m (k 1) schreiben
wir unter Hinzunahme des rechten Randwertes sk um als
y(x2 ; x1 , s1 ) s2
y(x3 ; x2 , s2 ) s3 s1
(S) :=
.
.. = 0,
S := ... . (7.8)
y(xk ; xk1 , sk1 ) sk sk
r(s1 , sk )
Zur Losung von (7.8) wird wiederum das Newtonverfahren eingesetzt. Wich-
tig ist dabei die spezielle Struktur der Jacobimatrix
G1 I
G2 I
.
D(S) =
.. ... ... ,
Gk1 I
A 0 ... 0 B
die man auch als Mehrzielmatrix bezeichnet. Die Blocke stellen ihrerseits
Funktionalmatrizen dar,
Gj = y(xj+1 ; xj , sj )/sj , A = r(s1 , sk )/s1 , B = r(s1 , sk )/sk .
Aufgrund dieser speziellen Struktur kann die Losung linearer Gleichungssy-
steme D(S i )S = (S i ) im Newtonverfahren mittels Blockelimination
erfolgen:
G1 s1 s2 = 1 ,
G2 s2 s3 = 2 ,
..
.
Gk1 sk1 sk = k1 ,
As1 + Bsk = k .
148
Durch sukzessives Eliminieren von s2 bis sk folgt das lineare Gleichungs-
system
(A + BGk1 Gk2 . . . G1 )s1 = r.S.
fu
r den Anfangswert s1 . Nach dessen Losung berechnet man s2 bis sk
u
ber die blockweise Ru
cksubstitution, so dass der Aufwand in der linearen
Algebra trotz Einfu
hrung der Mehrzielknoten nicht all zu viel hoher ist als
beim EinfachSchieen. Allerdings neigt das lineare Gleichungssystem fu r
s1 wie beim EinfachSchieen zu schlechter Kondition, weswegen man bei
empndlichen Problemen lieber die Gesamtmatrix D(S) mit einer QR-
Zerlegung behandelt (ist aufgrund der schlechten Kondition einer billigeren
LR-Zerlegung vorzuziehen).
Die lineare Algebra entscheidet jedoch weniger u ber die Ezienz als die Be-
rechnung der Mehrzielmatrix, die wie beim Einfach-Schieen u ber numeri-
sche Dierentiation oder die Sensitivitatengleichung bestimmt wird. Da die
Integrationsintervalle durch die intervallweise Vorgehensweise ku rzer sind,
ist der Aufwand kaum hoher als beim EinfachSchieen: Statt u ber ein lan-
ber k 1 Intervalle der Lange (b a)/(k 1)
ges Intervalls [a, b] wird u
integriert.
Numerische Software
Der Fortran-Code BOUNDS geht auf Bulirsch zuru ck und wurde seit 1971
fortlaufend weiterentwickelt. Eine aktuelle und vollstandige Neu-Implemen-
tierung ist der Code JANUS (Callies).
Die Idee der Kollokation (vergl. IRK-Verfahren) kann man auch zur Losung
von Randwertproblemen einsetzen. Ein bekannter Code hierzu ist COLSYS
(Ascher/Bader/Mattheij/Russel, ab 1981).
149