8 Numerik gewöhnlicher Differentialgleichungen...• das in dieser Gleichung stehende Integral...

8 NUMERIK GEWÖHNLICHER DIFFERENTIALGLEICHUNGEN 203

8 Numerik gewöhnlicher Differentialgleichungen

8.1 Grundlagen

In der Numerik von gewöhnlichen Differentialgleichungen werden vorwiegend Aufgaben fol-gender Art behandelt:

Definition 8.1.1 (Aufgabenstellung)

• gegeben ist eine Anfangswertaufgabe:

– explizites Differentialgleichungssystem1–Ordnung:

~y ′ = ~f (x,~y) mit ~f : R × Rn → Rn

– Anfangswert:

~y(x0) = ~y0 =

(y01, . . . , y

0n

)T ∈ Rn

• gesucht ist der Funktionswert (–Vektor)~y(b) an der rechts von x0 liegenden Stelle b.

Es ist also nicht die Lösungsfunktion als solche gesucht, sondern der Funktionswert der Lösungan einer bestimmten Stelleb > x0.

Bemerkung 8.1.2

• Unter vernünftigen Voraussetzungen an die Funktion~f — etwa:~f ist auf[x0, b] × Rn

– stetig

– bzgl.~y gleichmäßig lipschitzstetig, d.h. es gibt eine KonstanteL > 0 mit

||~f (x,~y1) −~f (x,~y2)|| ≤ L · ||~y1 −~y2||

für alle x∈ [x0, b] und~y1,~y2 ∈ Rn

ist die Anfangswertaufgabe8.1.1eindeutig lösbar.

• In vielen Anwendungen gilt n= 1, d.h. es handelt sich um eine explizite Differentialglei-chung 1. Ordnung.

Bemerkung 8.1.3 Differentialgleichungen (oder auch Differentialgleichungssysteme) höhererOrdnung können immer auf ein äquivalentes Differentialgleichungssystem 1. Ordnung zurück-geführt werden:Die DGL

z(n) = f (x, z, z′, z′′, . . . , z(n−1)) (6)

n–ter Ordnung.

Fachhochschule Aachen 2. Januar 2008

204 8.1 Grundlagen

wird durch:y1 = zy2 = z′

...yn−1 = z(n−2)

yn = z(n−1)

überführt in das Differentialgleichungssystem 1. Ordnung:

y′1 = y2y′2 = y3

...y′n−1 = yny′n = f (x, y1, y2, . . . , yn)

(7)

Eine Lösung z der Differentialgleichung6 der Ordnung n führt zur Lösung

~y = (z, z′, . . . , z(n−1))T

des Differentialgleichungssystems7 und umgekehrt eine Lösung

~y = (y1, y2, . . . , yn)T

des Differentialgleichungssystems7 führt zu der Lösung

z = y1

der Differentialgleichung6 der Ordnung n.

Beispiel 8.1.4 (2–Körper–Problem)

Befindet sich im Punkt~P1(t) im Raum eine Masse M1 und bewegt sich mit der Geschwindigkeit~P1(t)′, sowie in Punkt~P2(t) eine Masse M2 mit Geschwindigkeit~P′2(t), so beeinflussen sich dieMassen gegenseitig aufgrund der Gravitation. Es gilt:

~P′′1 = γ · M2 ·(~P2 − ~P1)||~P2 − ~P1||3

~P′′2 = γ · M1 ·(~P1 − ~P2)||~P2 − ~P1||3

(Ableitungen nach der Zeit t) wobeiγ = 6.674 · 10−11 m3

Kg · s2 die Gravitationskonstante ist.Mit

~Pi(t) =

Pix(t)Piy(t)Piz(t)

, ~P′i(t) =

P′ix(t)P′iy(t)P′iz(t)

und ~P′′i (t) =

P′′ix(t)P′′iy(t)P′′iz(t)

für i = 1, 2 und

y1(t) = P1x(t), y2(t) = P1y(t), y3(t) = P1z(t), y4(t) = P′1x(t), y5(t) = P

′1y(t), y6(t) = P

′1z(t)

Mathematik III/Numerik Fachbereich 8


y7(t) = P2x(t), y8(t) = P2y(t), y9(t) = P2z(t), y10(t) = P′2x(t), y11(t) = P

′2y(t), y12(t) = P

′2z(t)

führt das auf das Differentialgleichungssystem 1. Ordnung:

y′1 = y4y′2 = y5y′3 = y6

y′4 = γ · M2 ·y7 − y1

(√

(y7 − y1)2 + (y8 − y2)2 + (y9 − y3)2)3

y′5 = γ · M2 ·y8 − y2

(√

(y7 − y1)2 + (y8 − y2)2 + (y9 − y3)2)3

y′6 = γ · M2 ·y9 − y3

(√

(y7 − y1)2 + (y8 − y2)2 + (y9 − y3)2)3

y′7 = y10y′8 = y11y′9 = y12

y′10 = γ · M1 ·y1 − y7

(√

(y7 − y1)2 + (y8 − y2)2 + (y9 − y3)2)3

y′11 = γ · M1 ·y2 − y8

(√

(y7 − y1)2 + (y8 − y2)2 + (y9 − y3)2)3

y′12 = γ · M1 ·y3 − y9

(√

(y7 − y1)2 + (y8 − y2)2 + (y9 − y3)2)3

Satz 8.1.5 (Rückführung der Anfangswertaufgabe auf eine Integralgleichung)Es sei eine Anfangswertaufgabe

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

gegeben.Integriert man die Gleichung~y ′ = ~f (x,~y), unter der Voraussetzung, dass~y und~y ′ eine Funk-tion von x ist, über dem Intervall[x0, x] (d.h. linker Randpunkt ist das feste x0 und der rechteRandpunkt das “variable“ x), so erhält man (wegen der oberenGrenze x wird hier t als Inte-grationsvariable eingeführt, es wird komponentenweise integriert):

x∫

x0

~y ′(t) dt =

x∫

x0

~f (t,~y(t)) dt

⇔ ~y(x) −~y(x0) =x∫

x0

~f (t,~y(t)) dt

⇔ ~y(x) = ~y(x0) +x∫

x0

~f (t,~y(t)) dt mit~y(x0) = ~y0

⇔ ~y(x) = ~y0 +x∫

x0

~f (t,~y(t)) dt


206 8.1 Grundlagen

Die Anfangswertaufgabe ist also äquivalent zur Integralgleichung

~y(x) = ~y0 +

x∫

x0

~f (t,~y(t)) dt (8)

Diese Integralgleichung ist mathematisch gesehen genausokompliziert wie die Anfangswert-aufgabe, denn die Anfangswertaufgabe besteht (im Wesentlichen) aus einer Gleichung, in dereine unbekannte Funktion~y und deren Ableitung~y ′ vorkommen — in der Integralgleichungtaucht dieselbe unbekannte Funktion~y und ein Integralüber diese Funktion auf.

Bemerkung 8.1.6 (Ansatz für numerische Verfahren)

• Intervallaufteilung:

x0 < x1 < . . . < xk < xk+1 < . . . < xn = b

• unter der Voraussetzung, dass man~yk = ~y(xk) schon kennt, die Integralgleichung8 lokalfür das Intervall von xk bis xk+1 anwenden:

~y(xk+1) = ~yk+1 = ~yk +

xk+1∫

xk

~f (t,~y(t)) dt

• das in dieser Gleichung stehende Integralxk+1∫

xk

~f (t,~y(t)) dt

näherungsweise mit einer (interpolatorischen) Quadraturformel berechnen.

xk−2 xk−1 xk xk+1

f (x, y(x))



Beispiel 8.1.7 (explizites Euler–Verfahren)

Nimmt man zur näherungsweisen Berechnung des Integrals die(linksseitige)Rechteckregel:

xk+1∫

xk

~f (t,~y(t)) dt ≈ ~f (xk,~yk) · (xk+1 − xk)

Skizze:


b

f (xk, yk)

bekannt

f (x, y(x))

so erhält man die Formel:

~yk+1 ≈ ~yk + hk ·~f (xk,~yk) mit hk = xk+1 − xk

Das auf dieser Formel beruhende Verfahren heißtexplizites Euler–Cauchy–Polygonzugverfahren.

Algorithmus 8.1.8 (explizites Euler–Verfahren)

Gegeben: Anfangswertaufgabe~y ′ = ~f (x,~y),~y(x0) = ~y0.Gesucht: “Funktionswert“ (eigentlich: Vektor)~y(b) mit b > x0.Vorgehen:Unterteile das Intervall[x0, b] in viele kleine Teilintervalle

x0 < x1 < x2 < . . . < xn = b

und setze hk = xk+1 − xk (Teilintervalllänge des k-ten Intervalls).Ausgehend vom Anfangswert(x0,~y0) berechne iterativ für k= 0, 1, 2, . . . , n− 1:

~yk+1 = ~yk + hk ·~f (xk,~yk)Dann gilt:

~yn ≈ ~y(xn) = ~y(b)Man kann zeigen: Sei h= max{ hk | 0 ≤ k < n } die größte der beteiligten Teilintervalllängen,so gilt für den Fehler:

||~yn −~y(b)|| ∈ O(h)d.h. das Euler–Cauchy–Verfahren hat die Ordnung1.


208 8.1 Grundlagen

Beispiel 8.1.9

gegeben: Anfangswertaufgabe y′ =√

x + y︸︷︷︸

f (x,y)

, y(0) = 1.0

gesucht: y(1.0)

Es soll das Euler–Cauchy–Polygonzugverfahren mit konstanter Schrittweite h= 0.1 durchge-führt werden.

• Anfangswert: x0 = 0.0, y0 = 1.0

• x1 = 0.1, y1 = y0 + h · f (x0, y0) = 1.1000000000

• x2 = 0.2, y2 = y1 + h · f (x1, y1) = 1.2095445115

• x3 = 0.3, y3 = y2 + h · f (x2, y2) = 1.3282687513

• x4 = 0.4, y4 = y3 + h · f (x3, y3) = 1.4558723857

• x5 = 0.5, y5 = y4 + h · f (x4, y4) = 1.5921027929

• x6 = 0.6, y6 = y5 + h · f (x5, y5) = 1.7367438242

• x7 = 0.7, y7 = y6 + h · f (x6, y6) = 1.8896079411

• x8 = 0.8, y8 = y7 + h · f (x7, y7) = 2.0505305294

• x9 = 0.9, y9 = y8 + h · f (x8, y8) = 2.2193656718

• x10 = 1.0, y10 = y9 + h · f (x9, y9) = 2.3959829323

Ein Näherungswert für den gesuchten Funktionswert y(1.0) ist y10 = 2.3959829323.

Bemerkung 8.1.10 (geometrische Beschreibung des Euler–Cauchy–Verfahrens)

Im eindimensionalen Fall kann das explizite Euler–Cauchy-Polygonzugverfahren geometrischanhand des Richtungsfeldes der Differentialgleichung y′ = f (x, y) hergeletitet werden:Zur Anfangswertaufgabe y′ = f (x, y), y(x0) = y0 sei der Funktionswert der Lösung y(x) an derStelle b> x0 gesucht.

Hierzu wird auf der x–Achse die Strecke von x0 nach b in kleine Teilintervalle unterteilt:

x0 < x1 < x2 < . . . < xn = b

Aufgrund des Anfangswertes kennen wir den Funktionswert y0 der Lösung y(x) an der Stelle x0— die Lösung y(x) soll ja durch den Punkt(x0, y0) gehen.Ausgehend von diesem Anfangswert(x0, y0) wird jetzt wie folgt ein Näherungswert y1 für denFunktionswert y(x1) der Lösung y(x) an der Stelle x1 ermittelt:

1. Man setzt den Punkt(x0, y0) in die rechte Seite der Differentialgleichung ein und erhältdie Steigung y′0 = f (x0, y0) der (ansonsten unbekannten) Lösungskurve y(x) im Punkt(x0, y0).

2. Hiermit ist die Tangente der Lösungskurve im Punkt(x0, y0) bekannt: diese geht durchden Punkt(x0, y0) und hat die Steigung y′0 = f (x0, y0).



3. In einer (kleinen) Umgebung von x0 wird die Lösung y(x) ganz gut durch ihre Tangenteangenähert:

Es gilt nach Taylor:

y(x) ≈ y(x0) + y′(x0) · (x− x0)︸︷︷︸

Tangente von y(x) in (x0,y0)

Durch Einsetzen von x1 für x folgt mit y(x0) = y0 und y′(x0) = y′0 = f (x0, y0):

y(x1) ≈ y1 := y0 + f (x0, y0) · (x1 − x0)Mit dieser Formel kann also anhand des Anfangswertes(x0, y0) ein Näherungswert y1 fürden Funktionswert der (unbekannten) Lösung y(x) an der Stelle x1 ermittelt werden!

Skizzen:

a.) Anfangswert:

-

6

x0 x1

y0 (x0, y0)q

b.) Linienelement mit Steigung f(x0, y0) ermitteln:

-

6

x0 x1

y0 q��(x0, y0)

c.) Linienelement zu Tangentevervollständigen:

-

6

x0 x1

y0 q��

��

��

��(x0, y0)

d.) Tangente an der Stelle x1 auswerten:

-

6

x0 x1

y0 q��

��

��

��(x0, y0)

qy1 (x1, y1)

Mit diesem Verfahren bzw. der Formel

y1 = y0 + (x1 − x0) · f (x0, y0)ist es also gelungen, ausgehend vom Anfangspunkt(x0, y0) einen weiter rechts liegenden Punkt(x1, y1) zu ermitteln mit y1 ≈ y(x1).Dies kann nun iterativ fortgesetzt werden:

y2 = y1 + (x2 − x1) · f (x1, y1) , y2 ≈ y(x2)y3 = y2 + (x3 − x2) · f (x2, y2) , y3 ≈ y(x3)

...yn = yn−1 + (xn − xn−1) · f (xn−1, yn−1) , yn ≈ y(xn) = y(b)


210 8.1 Grundlagen

Die Linienelemente in den Punkten(xk, yk) legen die Geradenstücke von(xk, yk) nach(xk+1, yk+1)fest und man gelangt mit diesen Geradenstücken von einem Punkt zum nächsten:

Skizze:

-

6

x0

y0 ry0�

��

x1

r

y1��

x2

r

y2��

x3

r

y3

x4

r

y4PPPP

x5

r

y5HHH

x6

r

y6PPPP

x7

r

y7

x8

r

y8��

x9 = b

r

y9 ≈ y(b)

Diese Skizze motiviert auch den Namen Polygonzugverfahren, der in der Literatur ebenfalls fürdieses Euler–Cauchy–Verfahren verwendet wird.

Beispiel 8.1.11 (implizites Euler–Verfahren)

Nimmt man zur näherungsweisen Berechnung des Integrals die(rechtsseitige)Rechteckregel:

xk+1∫

xk

~f (t,~y(t)) dt ≈ ~f (xk+1,~yk+1) · (xk+1 − xk)

Skizze:


b

f (xk+1, yk+1)

unbekannt

f (x, y(x))

so erhält man die Formel:

~yk+1 ≈ ~yk + hk ·~f (xk+1,~yk+1) mit hk = xk+1 − xk

In dieser Formel taucht der zu berechnende Wert (Vektor)~yk+1 sowohl auf der linken Seite, alsauch auf der rechten Seite als zweites Argument der Funktion~f auf — es handelt sich somit umeineimplizite Gleichung für~yk+1.Diese implizite Gleichung muss dann gelöst werden — im Allgemeinen iterativ.Das auf dieser Formel beruhende Verfahren heißtimplizites Euler–Verfahren.



Algorithmus 8.1.12 (implizites Euler–Verfahren)


x0 < x1 < x2 < . . . < xn = b


~yk+1 = ~yk + hk ·~f (xk+1,~yk+1)

Diese (für jeden Schritt von xk nach xk+1 anfallende)implizite Gleichung muss iterativ ausge-hend von einem geeigneten Startvektor gelöst werden.Dann gilt:

~yn ≈ ~y(xn) = ~y(b)Auch dieses Verfahren hat die Ordnung1, d.h. mit h= max{ hk | 0 ≤ k < n } gilt für denFehler:

||~yn −~y(b)|| ∈ O(h)

Bemerkung 8.1.13 (Durchführung der Fixpunktiteration in jedem Schritt)

Die in jedem Schritt des impliziten Euler–Verfahrens (iterativ) zu lösende Fixpunktgleichung:

~yk+1 = ~yk + hk ·~f (xk+1,~yk+1)

kann umgeschrieben werden zu

~yk+1 −~yk︸︷︷︸

=: ~z

= hk ·~f (xk+1,~yk +~yk+1 −~yk︸︷︷︸

=: ~z

)

also zu~z = hk ·~f (xk+1,~yk +~z)

wobei~z = ~0 ein vernünftiger Startvektor für die Fixpunktiteration ist.Hierbei kann

• die Gleichung~z = hk ·~f (xk+1,~yk +~z)

selbst als Fixpunktgleichung verwendet und mit dieser Gleichung die Fixpunktiterationdurchgeführt werden.

• die Gleichung umgeformt werden zu einer Nullstellengleichung:

~g(~z) = ~z− hk ·~f (xk+1,~yk +~z) = ~0

für welche dann das (ggf. mehrdimensionale) Newton–Verfahren zur Nullstellenbestim-mung (mit Startvektor~0) durchgeführt werden kann.


212 8.1 Grundlagen

Beispiel 8.1.14gegeben: Anfangswertaufgabe y′ =

√x + y

︸︷︷︸

f (x,y)

, y(0) = 1.0

gesucht: y(1.0)

Es soll das implizite Euler–Verfahren mit konstanter Schrittweite h= 0.1 durchgeführt werden.

Die in jedem Schritt zu lösende implizite Gleichung lautet:

z = h · f (xk+1, yk + z) = h ·√

xk+1 + yk + z

a.) Jeweilige Berechnung der impliziten Gleichung mit dem allgemeinen Iterationsverfahren

b.) Jeweilige Berechnung der impliziten Gleichung mit dem Newton–Verfahren

Startwert jeweils0 (die Iteration wird abgebrochen, sobald sich zwei iterierte Werte um wenigerals10−10 unterscheiden):

k xk yk Anzahl Iterationena.) b.)

0 0.000000 1.0000000000 − −1 0.100000 1.1100000000 8 42 0.200000 1.2295643924 8 43 0.300000 1.3583409811 8 34 0.400000 1.4960376646 8 35 0.500000 1.6424073100 8 36 0.600000 1.7972374575 8 37 0.700000 1.9603430380 8 38 0.800000 2.1315610586 8 39 0.900000 2.3107466223 8 310 1.000000 2.4977698797 7 3

Beispiel 8.1.15 (Trapezregel)

Die Berechnung des Integrals mit der Trapezregel:

xk+1∫

xk

~f (t,~y(t)) dt ≈ (xk+1 − xk)2

·(

~f (xk,~yk) +~f (xk+1,~yk+1))

Skizze:




b

bf (xk, yk)f (xk+1, yk+1)

bekannt unbekannt

f (x, y(x))

führt mit hk = xk+1 − xk auf die ebenfallsimplizite Gleichung:

~yi+k = ~yk +hk2·~f (xk,~yk) +

hk2·~f (xk+1,~yk+1)

für den zu berechnenden Wert~yk+1.Das auf dieser Formel beruhende (implizite) Verfahren für Differentialgleichungsn heißtTra-pezregel.

Algorithmus 8.1.16 (Trapezregel)


x0 < x1 < x2 < . . . < xn = b


~yk+1 = ~yk +hk2·~f (xk,~yk) +

hk2·~f (xk+1,~yk+1)

Diese (für jeden Schritt von xk nach xk+1 anfallende)implizite Gleichung muss iterativ ausge-hend von einem geeigneten Startvektor gelöst werden.Dann gilt:

~yn ≈ ~y(xn) = ~y(b)Dieses Verfahren hat die Ordnung2, d.h. mit h= max{ hk | 0 ≤ k < n } gilt für den Fehler:

||~yn −~y(b)|| ∈ O(h2)

Bemerkung 8.1.17Mit ~z = ~yk+1 −~yk kann diese implizite Gleichung umgeschrieben werden zu

~z =h2·~f (xk,~yk) +

h2·~f (xk+1,~yk +~z)

und diese Gleichung kann (ausgehend vom Startvektor~0)


214 8.1 Grundlagen

• entweder direkt zur Iteration verwendet werden,

• oder (nach Umformung in eine Nullstellengleichung) mit demNewton–Verfahren gelöstwerden.

Beispiel 8.1.18gegeben: Anfangswertaufgabe y′ =

√x + y

︸︷︷︸

f (x,y)

, y(0) = 1.0

gesucht: y(1.0)

Es soll die Trapezregel mit konstanter Schrittweite h= 0.1 durchgeführt werden.Die in jedem Schritt zu lösende implizite Gleichung lautet:

z =h2· f (xk, yk) +

h2· f (xk+1, yk + z) =

h2· √xk + yk +

h2· √xk+1 + yk + z

a.) Jeweilige Berechnung der impliziten Gleichung mit dem allgemeinen Iterationsverfahren

b.) Jeweilige Berechnung der impliziten Gleichung mit dem Newton–Verfahren

Startwert jeweils0 (die Iteration wird abgebrochen, sobald sich zwei iterierte Werte um wenigerals10−10 unterscheiden):

k xk yk Anzahl Iterationena.) b.)

0 0.000000 1.0000000000 − −1 0.100000 1.1048835949 7 32 0.200000 1.2193351156 7 33 0.300000 1.3429926794 7 34 0.400000 1.4755578207 7 35 0.500000 1.6167790988 7 36 0.600000 1.7664410793 7 37 0.700000 1.9243566154 7 38 0.800000 2.0903612584 7 39 0.900000 2.2643090958 6 310 1.000000 2.4460695821 6 3

Definition und Bemerkung 8.1.19

Zur Berechnung des neuen Wertes~yk+1 wird in allen Verfahren das Integral

xk+1∫

xk

~f (t,~y(t)) dt

über ein Interpolationspolynom~ϕ(t) für die unbekannte Funktion~f (t,~y(t)) näherungsweiseberechnet.Je nachdem, wo die Stützstellen dieses Interpolationspolynoms liegen, unterscheidet man dieVerfahrenstypen:



Einschrittverfahren: Bei den Einschrittverfahren liegenalle Stützstellen des Interpolations-polynoms~ϕ im aktuellen Integrationsintervall[xk, xk+1]:

xk−1 xk xk+1b b b

tk1 tk2 . . . tkm

bb

b

f (x, y)

ϕ(x)

Mehrschrittverfahren Bei den Mehrschrittverfahren sind die Stützstellen des Interpolations-polynoms ausschließlich unter den bereits durch die Intervallaufteilung gegebenen x–Werten:


b

b

b

bϕ(x)

f (xk−2, yk−2)f (xk−1, yk−1)

f (xk, yk)f (xk+1, yk+1)

f (x, y(x))

b b b b

bekannt unbekannt

Je nachdem, ob xk+1 zu den Stützstellen des Interpolations dazugehört oder nicht, sprichtman von einemexpliziten (gehört nicht dazu) oder einemimpliziten (gehört dazu) Mehr-schrittverfahren.

8.2 Einschrittverfahren

Definition und Bemerkung 8.2.1 (Runge–Kutta–Verfahren)

Zur Berechnung des nächsten Näherungswertes~yk+1 für~y(xk+1) ausgehend vomaktuellenWert~yk werden im Intervall[xk, xk+1] — die Länge des Interalls sei h= xk+1−xk — die m Stützstellenxkj = xk + aj · h für das Interpolationspolynom zugrunde gelegt.Die systematische Berechnung des Interpolationspolynoms~ϕ(x) (mit zum Teil noch unbekann-


216 8.2 Einschrittverfahren

ten Stützwerten) und die anschließende näherungsweise Berechnung des Integrals

xk+1∫

xk

~f (t,~y(t)) dt

über dieses Interpolatinspolynom~ϕ(x) führt auf Gleichungen folgender Form:

K1 = ~f(

xk + a1h,~yk + h ·m∑

j=1b1j · Kj

)

K2 = ~f(

xk + a2h,~yk + h ·m∑

j=1b2j · Kj

)

...

Km = ~f(

xk + amh,~yk + h ·m∑

j=1bmj · Kj

)

(9)

und

~yk+1 = ~yk + h ·m∑

j=1

cj · Kj

mit fest vorgegebenen, von den konkreten Stützstellen (undvom Verfahren) abhängenden Kon-stanten ai, ci und bij , 1 ≤ i, j ≤ m.Das auf diesen Formeln beruhende Verfahren heißtm–stufiges Runge–Kutta–Verfahren.Die Koeffizienten, welche dieses Verfahren vollständig beschreiben, werden häufig in einer Ta-belle, dem sogenanntenButcher–Diagramm, angegeben:

a1 b11 b12 . . . b1ma2 b21 b22 . . . b2m...

......

. . ....

am bm1 bm2 . . . bmmc1 c2 . . . cm

Bemerkung 8.2.2

Die bislang kennengelernten Verfahren sind von diesem Typ:

1. explizites Euler–Verfahren:~yk+1 = ~yk + h ·~f (xk,~yk)

︸︷︷︸

=: K1

Stufenzahl m ist1.K1 = ~f (xk + 0 · h,~yk + h · 0 · K1)

und~yk+1 = ~yk + h · 1 · K1

Butcher–Diagramm für das explizite Euler–Verfahren:

0 01



2. implizites Euler–Verfahren:

~yk+1 = ~yk + h ·~f (xk+1,~yk+1)︸︷︷︸

=: K1

Stufenzahl m ist1.K1 = ~f (xk + 1 · h,~yk + h · 1 · K1)

und~yk+1 = ~yk + h · 1 · K1

Butcher–Diagramm für das implizite Euler–Verfahren:

1 11

3. Trapezregel:

~yk+1 = ~yk + h ·(1

2·~f (xk,~yk)︸︷︷︸

=: K1

+1

2·~f (xk+1,~yk+1)︸︷︷︸

=: K2

)

Stufenzahl m ist2.

K1 = ~f(

xk + 0 · h,~yk + h · (0 · K1 + 0 · K2))

K2 = ~f(

xk + 1 · h,~yk + h · (12 ·K1 + 12 · K2))

und

~yk+1 = ~yk + h · (1

2· K1 +

1

2· K2)

Butcher–Diagramm für die Trapezregel:

0 0 0

1 12

12

12

12

Definition 8.2.3

1. Ein Runge–Kutta–Verfahren, bei dem die Matrix der Koeffizienten bij einestrikte untereDreiecksmatrix ist:

a1 0 0 0 . . . 0a2 b21 0 0 . . . 0a3 b31 b32 0 . . . 0...

......

. . . . . ....

am bm1 bm2 . . . bm(m−1) 0c1 c2 . . . cm−1 cm



heißt (m–stufiges)explizites Runge–Kutta–Verfahren.

Bei den meisten expliziten Runge–Kutta–Verfahren gilt darüber hinaus, dass a0 = 0 gilt.

In diesem Fall werden die Einträge gleich0 im Butcher–Diagramm fortgelassen, so dassein solches explizites Runge–Kutta–Verfahren dann durch das Schema

a2 b21a3 b31 b32...

......

. . .am bm1 bm2 . . . bm(m−1)

c1 c2 . . . cm−1 cm

beschrieben wird.

2. Ein Runge–Kutta–Verfahren, bei dem die Matrix der Koeffizienten bij eineuntere Drei-ecksmatrix ist:

a1 b11 0 . . . 0

a2 b21 b22. . .

......

......

. . . 0am bm1 bm2 . . . bmm

c1 c2 . . . cm

heißt (m–stufiges)diagonal–implizites Runge–Kutta–Verfahren.

Sind darüberhinaus alls Diagonalelemente bjj gleich, so spricht man von einemeinfachdiagonal–impliziten Runge–Kutta–Verfahren.

3. Ist die Matrix bij keineuntere Dreiecksmatrix, so spricht man von einem (m–stufigen) vollimpliziten Runge–Kutta–Verfahren.

Bemerkung 8.2.4

• Bei einemexpliziten Runge–Kutta–Verfahren ist das System9 zur Bestimmung derK’sein explizites Gleichungssystem, d.h.

1. zuerst kannK1 ausgerechnet werden,

2. anschließend kannK2 mit Hilfe vonK1 berechnet werden,

3. dann kannK3 mit Hilfe vonK1 undK2 berechnet werden

4. usw.

Nach der Berechnung allerK’s wird ~yk+1 mit Hilfe derK’s berechnet.

• Bei einemimpliziten Runge–Kutta–Verfahren ist das Gleichungssystem9 zur Bestim-mung derK’s ein implizites Gleichungssystem, welches iterativ zu lösen ist.

Hierbei kann



– entweder das System9 selbst zur Fixpunktiteration verwendet,

– oder (nach Umformung in eine Nullstellengleichung) ein Newton–Verfahren ange-wendet werden.

Als Startwert kann man jeweilsKj = ~f (xk,~yk) für alle j verwenden.

Es kann gezeigt werden, dass bei vernünftigen Aufgabenstellungen und hinreichend klei-nen Schrittweiten h das System9 für dieK’s eine eindeutige Lösung hat.

Definition 8.2.5 (Konsistenzordnung eines Einschrittverfahrens)

Sei~yk+1 der mittels des Einschrittverfahrens ausgehend von der Stelle xk und exaktem An-fangswert~yk = ~y(xk) berechnete Näherungswert für die exakte Lösung~y(xk+1) an der Stellexk+1 = xk + h und gilt:

~y(xk+1) −~yk+1h

∈ O(hp)

so heißt p dieKonsistenzordnungdes Verfahrens und das Verfahren heißtkonsistent, falls p≥ 1gilt.(Der Bruch in obiger Gleichung ist der durch h geteiltelokale Fehler, also der Fehler, der imPunkt xk+1 vorliegt, falls~yk im Punkt xk noch fehlerfrei war.)

Definition und Bemerkung 8.2.6 (Konvergenzordnung eines Einschrittverfahrens)

Sei~yn der mittels Intervallaufteilung

x0 < x1 < . . . < xn = b

ausgehend von dem Anfangswert~y(x0) = ~y0 bestimmte Näherungswert für die Lösung~y(xn) =~y(b), so heißt p dieKonvergenzordnungdes Verfahrens, wenn

(~y(xn) −~yn) ∈ O(hp)

gilt (globaler Fehler).Man kann zeigen, dass bei einem Einzelschrittverfahren dieKonvergenzordnung und die Kon-sitenzordnung gleich sind!

Algorithmus 8.2.7 (verbessertes Euler–Verfahren)

Butcher–Diagramm:

0 0 012

12

0

0 1

(es handelt sich somit um ein2–stufigesexplizitesVerfahren!)In Formeln:

K1 = ~f (xk,~yk)K2 = ~f (xk + h2 ,~yk +

h2K1)



und~yk+1 = ~yk + h · K2

oder kürzer in einerFormel:

~yk+1 = ~yk + h ·~f (xk +h2,~yk +

h2·~f (xk,~yk))

Die Konvergenzordnung des verbesserten Euler–Verfahrensist 2.

Algorithmus 8.2.8 (Verfahren von Heun)

Butcher–Diagramm:

0 0 0

1 1 012

12

(es handelt sich somit um einexplizitesVerfahren, die Stufenzahl m ist2!)In Formeln:

K1 = ~f (xk, yk)K2 = ~f (xk + h , ~yk + h · K1)

und

~yk+1 = ~yk + h ·1

2· (K1 + K2)

Die Ordnung des Verfahrens ist2.

Algorithmus 8.2.9 (Klassisches Runge–Kutta–Verfahren)

Butcher–Diagramm:

0 0 0 0 012

12

0 0 012

0 12

0 0

1 0 0 1 016

13

13

16

(Wiederum einexplizitesVerfahren, die Stufenzahl m ist4!)In Formeln:

K1 = ~f (xk,~yk)K2 = ~f (xk + h2 , ~yk +

h2· K1)

K3 = ~f (xk +h2

, ~yk +h2· K2)

K4 = ~f (xk + h , ~yk + h · K3)und

~yk+1 = ~yk + h · (1

6· K1 +

1

3· K2 +

1

3· K3 +

1

6· K4)

Das klassische Runge–Kutta–Verfahren hat die Ordnung4.



Beispiel 8.2.10

gegeben: Anfangswertaufgabe y′ =√

x + y︸︷︷︸

f (x,y)

, y(0) = 1.0

gesucht: y(1.0)

Es soll das klassische Runge–Kutta–Verfahren mit konstanter Schrittweite h= 0.2 durchgeführtwerden.

k xk yk K’s

0 0.0000000000 1.0000000000 K1 = 1.0000000000K2 = 1.0954451150K3 = 1.0997929403K4 = 1.1916201526

1 0.2000000000 1.2194032088 K1 = 1.1913870944K2 = 1.2800554356K3 = 1.2835142198K4 = 1.3697102076

2 0.4000000000 1.4756777625 K1 = 1.3695538553K2 = 1.4534899890K3 = 1.4563745265K4 = 1.5384903860

3 0.6000000000 1.7666035383 K1 = 1.5383769168K2 = 1.6187776963K3 = 1.6212591736K4 = 1.7002515617

4 0.8000000000 2.0905602789 K1 = 1.7001647799K2 = 1.7778011016K3 = 1.7799832553K4 = 1.8564904874

5 1.0000000000 2.4463010782

Der mit dem klassischen Runge–Kutta–Verfahren berechneteNäherungswert lautet:

y(1.0) ≈ 2.4463010782Algorithmus 8.2.11 (Gauß–Verfahren der Ordnung4)Butcher–Diagramm:

3−√

36

14

3−√

312

3+√

36

3+√

312

14

12

12

Es handelt sich somit um ein2–stufiges implitzites Verfahren.In Formeln:

K1 = ~f (xk +3 −

√3

6· h , ~yk + h · (

1

4· K1 +

3 − 2√

3

12· K2))

K2 = ~f (xk +3 +

√3

6· h , ~yk + h · (

3 + 2√

3

12· K1 +

1

4· K2))


222 8.3 Schrittweitensteuerung

und

~yk+1 = ~yk +h2· (K1 + K2)

Algorithmus 8.2.12 (Gauß–Verfahren der Ordnung6)Butcher–Diagramm:

12−

√15

10536

29−

√15

15536

−√

1530

12

536

+√

1524

29

536

−√

154

12

+√

1510

536

+√

1530

29

+√

1515

536

518

49

518

Es handelt sich somit um ein3–stufiges implitzites Verfahren.(Auf eine explizite Angabe der Formelm wird hier verzichtet!)

8.3 Schrittweitensteuerung

Ist eine Anfangswertaufgabe~y′ = ~f (x,~y) , ~y(x0) = ~y0

gegeben und ist der Funktionswert~y(b) der Lösung~y(x) an der Stelleb > x0 gesucht, so mussman, um mit einem numerischen Verfahren die Lösung einigermaßen genau zu bestimmen, dasIntervall in Teilintervalle aufteilen:

x0 < x1 < x2 < . . . < xn = b

um dann ausgehend vom Funktionswert~y0 an der Anfangsstellex0 schrittweise über~y1 ≈ ~y(x1),~y2 ≈ ~y(x2), . . . zum gewünschten Näherungswert~yn ≈ ~y(xn) = ~y(b) zu gelangen.Wählt man (zu) viele kleine Teilintervalle, so ist das Ergebnis zwar sehr genau, dafür ist derRechenaufwand aber sehr hoch.Wählt man (zu) wenige große Teilintervalle, ist der Rechenaufwand sehr gering, dafür dürftedas Ergebnis aber zu ungenau sein.Wie bei der numerischen Integration (Adaptive Quadratur) ist es anzustreben, bei einer vorge-gebenen Genauigkeitsschranke(b− x0) · ε > 0 die Schrittweiten der Intervallaufteilung jeweilsso groß zu wählen, dass die vorgegebene Genauigkeitsschranke gerade noch unterboten wird.Tunlichst sollte die Wahl der “passenden“ Schrittweiten dem numerischen Verfahren selbstüberlassen werden, so dass das Verfahren selbständig die lokalen Schrittweiten den lokalenBedürfnissen anpasst.Bei denadaptivenVerfahren ist es in jedem Berechnungsschritt (Berechnung von ~yk+1 aus-gehend vonxk und~yk) erforderlich, den hierbei gemachtenlokalen Fehlerzu schätzen, umentscheiden zu können,

• ob das Ergebnis genau genug ist und zum nächsten Teilintervall übergegangen werdenkann



• oder ob das Ergebnis zu ungenau ist und ausgehend vonxk mit einer kleineren Schrittweiteneu gerechnet werden muss.

Es gibt unterschiedliche Möglichkeiten, denlokalen Fehlerzu schätzen:

Bemerkung 8.3.1 (Schätzung des lokalen Fehlers)

1. Rechnung mit einemVerfahren und unterschiedlichen Schrittweiten:

Sei V ein Verfahren mit Konsistenzordnung (=Konvergenzordnung) q, weiter sei

• Ỹ der mittels dieses Verfahrens und Schrittweite h= xk+1 − xk in einemSchritt mitSchrittweite h ausgehend von(xk,~yk) berechnete Näherungswert für~y(xk+1),

• Y der mittels dieses Verfahrens mit Schrittweiteh2

in zweiSchritten ausgehend von(xk,~yk) berechnete Näherungswert für~y(xk+1)

so ist1

2q − 1(Y − Ỹ)

eine Schätzung für den lokalen Fehler der (besseren) NäherungY und

F = 12q − 1 ||(Y − Ỹ)||

ist ein Schätzwert für den Betrag eben dieses (lokalen) Fehlers.

2. Rechnung mit zwei unterschiedlichen Verfahren und gleicher Schrittweite:

Sei V ein Verfahren der Konsistenzordnung (=Konvergenzordnung) q und̃V ein Verfahrenmit Konsistenzordnung mindestens q+ 1, weiter sei

• Y der ausgehend von(xk,~yk) mit dem Verfahren V und Schrittweite h= xk+1 − xkberechnete Näherungswert für~y(xk+1), und

• Ỹ der ausgehend von(xk,~yk) mit dem VerfahreñV und derselben Schrittweite hberechnete Näherungswert für~y(xk+1),

so ist die DifferenzY − Ỹ

eine Schätzung für den lokalen Fehler der (schlechteren) NäherungY und

F = ||Y − Ỹ||

ist ein Schätzwert für den Betrag eben dieses (lokalen) Fehlers.

Ein möglicher Ansatz zur Schrittweitensteuerung ist dann (nach dem BuchNumerik–Algorithmenvon Frau Prof. Dr. G. Engeln–Müllges) der folgende:



Algorithmus 8.3.2 (Schrittweitensteuerung)

Es sei eine Fehlerschranke(b− x0) · ε > 0 vorgegeben.In folgendem Algorithmus wird sichergestellt, dass in jedem Integrationsschritt der Intervall-länge h der lokale Fehler ungefähr gleich h· ε ist.Ausgehend vom aktuellen Punkt(xk,~yk) und Schrittweite h= xk+1−xk werden zwei Näherungs-werteY bzw.Ỹ für den nächsten Wert~y(xk+1) berechnet,entweder

• mit ein und demselben Verfahren V der Ordnung q und den beidenSchrittweiten h (Nä-herungỸ und h

2(NäherungY)

• oder mit zwei Verfahren V der Ordnung q (NäherungY) undṼ der Ordnung mindestensq + 1 (NäherungỸ)

und der Betrag des hierbei gemachten FehlersF nach Bemerkung8.3.1geschätzt.Mit dieser FehlerschätzungF wird die folgende Steuergröße berechnet:

S= q√

h · εF

• Ist S ≥ 1, so wird der Fehler als klein genug angesehen, der WertY wird als Nähe-rungswert für y(xk+1) akzeptiert, man geht zum nächsten x–Wert xk+1 über und wählt alsneue Schrittweite:

h = min{ 2, S} · h

d.h. die Schrittweite wird nächsten Schritt geringfügig vergrößert,• Ist S< 1, so wird der Fehler als zu groß angesehen, der Schritt ausgehend von xk mussmit der kleineren Schrittweite:

h = max{ 12, S} · h

und neuem xk+1 = xk + h wiederholt werden.

8.3.1 Einbettungsformeln

erleichtern die Schrittweitensteuerung, wenn pro Schrittmit zwei unterschiedlichen VerfahrenNäherungen berechnet werden.Hierbei ist es hilfreich, wenn die beiden verwendeten VerfahrenV und Ṽ so in einem Zusam-menhang stehen, dass etwa die Zwischen– und Hilfswerte, welche zur Ausführung eines Schrit-tes mit dem VerfahrenV berechnet werden, bei der Durchführung des Schrittes mit VerfahrenṼ wiederverwendet werden können, so dass fürṼ nicht alles neu berechnet werden muss.Solche Paare zusammenhängender VerfahrenV undṼ heißenEinbettungsformeln.Einbettungsformeln sind vor allem bei expliziten Runge–Kutta–Verfahren weit verbreitet.



Beispiele hierzu sind das verbesserte Polygonzugverfahren (Ordnung2, übernimmt die Rollevon VerfahrenV) mit dem Koeffizientenschema:

12

12

0 1

ausgeschrieben: K1 = f (xi, yi)K2 = f (xi +

12· h, yi + h · 12 · K1)

undY = yi + h · (0 · K1 + 1 · K2)

und folgendes Runge–Kutta–Verfahren der Ornung3 (übernimmt die Rolle von VerfahreñV),Koeffizientenschema: 1

212

1 −1 216

23

16

ausgeschrieben:K1 = f (xi, yi)

K2 = f (xi + 12 · h, yi + h · 12 · K1)K3 = f (xi + 1 · h, yi + h · (−1 · K1 + 2 · K2))

undỸ = yi + h · (16 · K1 + 23 · K2 + 16 · K3)

Wie man sieht, werden in beiden VerfahrenV und Ṽ die GrößenK1 und K2 gleich berechnetund verwendet — bei VerfahreñV kommt zusätzlich einK3 hinzu.Bei den Formeln für die NäherungswerteY bzw. Ỹ treten dann unterschiedliche Koeffizientenauf.Man kann die Koeffizienten beider Verfahren in einSchema hineinschreiben:

12

12

1 −1 2V 0 1

Ṽ 16

23

16

Bemerkung 8.3.3 (Koeffizientenschema für Einbettungsformeln)

Das Koeffizientenschema einer allgemeinen (expliziten) Einbettungsformel lautet:

a2 b21a3 b31 b32...

.... . . . . .

am bm1 . . . . . . bm,m−1V c1 c2 . . . cm−1 cmṼ c̃1 c̃2 . . . c̃m−1 c̃m

(i. Allg. sind die letzten Koeffizienten cm, cm−1,. . . gleich0)



In jedem Schritt mit diesen Verfahren werden die Hilfsgrößen wie folgt berechnet:

K1 = f (xi, yi)K2 = f (xi + a2 · h, yi + h · b21 · K1)...

......

Km = f (xi + am · h, yi + h ·m−1∑

l=1

bmlKl)

und anschließend mit denselben Hilfsgrößen die beiden Näherungen:

Y = yi + h ·m∑

i=1

ci · Ki

und

Ỹ = yi + h ·m∑

i=1

c̃i · Ki

berechnet.

Bemerkung 8.3.4 (Einbettungsformel von England)

Eine bekanntes Formelpaar von diesem Typ sind die England–Formeln mit Ordnung q= 4 fürV und5 für Ṽ. Sie haben folgendes Koeffizientenschema:

12

12

12

14

14

1 0 −1 223

727

1027

0 127

15

28625

−125625

546625

54625

−378625

V 16

0 46

16

Ṽ 14336

0 0 35336

162336

125336

Beispiel 8.3.5 (Beispiel zur Schrittweitensteuerung)

Beute–Räuber–Differentialgleichungssystem (Lotka–Volterra):

y′1(t) = a · y1(t) · (1 − y2(t))y′2(t) = y2(t) · (y1(t) − 1)

y1(t) = “Beute“ und y2(t) = “Räuber“.Zeitliche Entwicklung der Populationen, berechnet mit derFormel von England(a = 10.0, y1(0) = 3.0, y2(0) = 1.0):



0

0.5

1

1.5

2

2.5

3

3.5

4

0 1 2 3 4 5

"beute""raeuber"

Entwicklung der Schrittweiten:

0

0.02

0.04

0.06

0.08

0.1

0 1 2 3 4 5

"schrittweite"


228 8.4 Stabilität, steife Differentialgleichungen

8.4 Stabilität, steife Differentialgleichungen

In der Praxis auftretende Differentialgleichungssystemehaben häufig die Eigenschaft, dass sichdie Lösungen aus “abklingenden“ e–Funktioneneλt zusammensetzen mitλ ∈ C und negativemRealteil Re(λ) < 0.Beachte: istλ = −a + b · i mit a, b ∈ R unda > 0, so ist:

eλt = e−at · (cos(bt) + i · sin(bt))

eine oszillierende, aber wegena > 0 betragsmäßig abklingende (komplexe) Funktion.Man möchte natürlich, dass auch eine numerisch bestimmte Näherungslösung für das Differen-tialgleichungssystem genauso “abklingende“ Bestandteile hat.Hierzu sind je nach Verfahren an die zu verwendende Schrittweiteh Bedingungen zu knüpfen!

Definition 8.4.1 (Testproblem)

Das Anfangswertproblemy′(x) = λ · y(x) , y(0) = 1

(λ ∈ C, Re(λ) < 0) heißtTestproblem.Die exakte Lösung des Testproblems lautet bekanntlich

y(x) = eλx

Da Re(λ) < 0 handelt es sich bei der Lösung um eine betragsmäßig monoton fallende Funktion.

Ist xk+1 = xk + h, so gilt für die exakte Lösung:

y(xk+1) = eλxk+1 = eλ(xk + h) = eλh · eλxk = eλh · y(xk)

d.h. die exakte Lösung y(xk+1) an der Stelle xk+1 erhält man aus der exakten Lösung y(xk) an

der Stelle xk durch Multiplikation mit dem Faktor eλh.Wegen Re(λ) < 0 ist der Betrag dieses Faktors eλh kleiner als1.

Bemerkung 8.4.2 (numerische Verfahren für die Testaufgabe)

Wendet man ein Runge–Kutta–Verfahren auf die Testdifferentialgleichung y′ = λy an, so gehtebenfalls der Näherungswert yk+1 für y(xk+1) durch Multiplikation mit einem Faktor F(λh) ausdem Näherungswert yk für y(xk) hervor, d.h. es gibt eine (komplexe) Funktion F(z), so dassdieser Faktor gleich F(λh) ist:

yk+1 = F(λh) · ykDiese (komplexe) Funktion F(z) heißtStabilitätsfunktion des numerischen Verfahrens.

Bemerkung 8.4.3 (Explizites Euler-Verfahren für das Testproblem)

Wendet man das explizite Euler–Verfahren mit Schrittweiteh auf das Testproblem an, so erhältman

y(xk+1) ≈ yk+1 = yk + h · f (xk, yk) = yk + hλ · yk = (1 + λh) · yk



d.h. die Stabilitätsfunktion des expliziten Euler–Verfahrens ist F(z) = 1 + z.

Damit das qualitative Verhalten der exakten Lösung (betragsmäßig abklingend) auch für dienumerisch bestimmte Näherungslösung gilt, ist zu fordern,dass auch der beim Euler–Verfahrenverwendete Faktor F(λh) = (1 + λh) betragsmäßig kleiner als1 ist, es ist also zu fordern:

|F(λh)| < 1 mit F(z) = 1 + z

Ist diese (vom Wertλ abhängende) Bedingung für die verwendete Schrittweite h verletzt, sokommt mit dem Euler–Verfahren ein unsinniges Ergebnis heraus!

Beispiel 8.4.4Es wird das Testproblem mitλ = −100 und das explizite Euler–Verfahren zu-grundegelegt.

Die exakte Lösung y(x) = e−100x ist eine extrem schnell abklingende e–Funktion.Die “Stailitätsbedingung“ für das Euler–Verfahren|1 + λh| < 1 führt mitλ = −100 auf

h < 0.02

• Wählt man die Schrittweite h= 0.005, so ist diese Bedingung erfüllt und die für dasIntervall [0, 1] mit dem Euler–Verfahren berechneten Lösung sieht wie folgtaus:

0

0.2

0.4

0.6

0.8

1

0 0.2 0.4 0.6 0.8 1

"h=0.005.txt"

Wie man sieht, stimmt die Lösung mit der exakten Lösung sehr gut überein.

• Wählt man als Schrittweite h= 0.025, so ist diese Bedingung verletzt und die für dasIntervall [0, 1] mit dem Euler–Verfahren berechnete Näherungslösung



-8e+06

-6e+06

-4e+06

-2e+06

0

2e+06

4e+06

6e+06

8e+06

1e+07

1.2e+07

0 0.2 0.4 0.6 0.8 1

"h=0.025.txt"

hat mit der exakten Lösung nicht viel zu tun.

Bemerkung 8.4.5 (Klassisches Runge–Kutta–Verfahren für das Testproblem)Wendet man das klassische Runge–Kutta–Verfahren auf die Testdifferentialgleichung y′ = λyan, so erhält man die Formel:

yk+1 =

(

1 + (λh) +(λh)2

2+

(λh)3

3!+

(λh)4

4!

)

· yk

D.h. die Stabilitätsfunktion für das klassiche Runge–Kutta–Verfahren lautet:

F(z) = 1 + z+z2

2+

z3

3!+

z4

4!

Um auch hierabklingendesVerhalten der Näherungslösung zu bekommen, muss für die zuverwendende Schrittweite h

|F(λh)| < 1

gefordert werden.

Beispiel 8.4.6Es wird das Testproblem mitλ = −100 und das explizite klassiche Runge–Kutta–Verfahren zugrundegelegt.Die exakte Lösung y(x) = e−100x ist eine extrem schnell abklingende e–Funktion.Die “Stailitätsbedingung“ für das klassiche Runge–Kutta–Verfahren führt mitλ = −100 auf

h < 0.0278

• Wählt man die Schrittweite h= 0.01, so ist diese Bedingung erfüllt und die für dasIntervall [0, 1] mit dem klassichen Runge–Kutta–Verfahren berechnete Lösung sieht wiefolgt aus:



0

0.2

0.4

0.6

0.8

1

0 0.2 0.4 0.6 0.8 1

"rk_0_01.txt"

Wie man sieht, stimmt die Lösung mit der exakten Lösung sehr gut überein.

• Wählt man als Schrittweite h= 0.03, so ist diese Bedingung verletzt und die für dasIntervall [0, 1] mit dem klassichen Runge–Kutta–Verfahren berechnete Näherungslösung

-10000

0

10000

20000

30000

40000

50000

60000

0 0.2 0.4 0.6 0.8 1

"rk_0_03.txt"

hat mit der exakten Lösung nicht viel zu tun.

Bemerkung 8.4.7 In den Anwendungen bei Differentialgleichungsen ist das Problem, dass sichdie Lösungen häufig aus abklingenden e–Funktionen mit unterschiedlichen “Abklingkonstan-ten“ λ1 undλ2 zusammensetzen, etwa:

y(t) = c1 · e−1 · t︸︷︷︸y1(t)

+c2 · e−100 · t︸︷︷︸y2(t)

+ . . .

hier hat manλ1 = −1 undλ2 = −1000.



Im Vergleich zu y1(t) ist der Bestandteil y2(t) “extrem schnell abklingend“ und schon für kleineWerte für t spielt y2(t) in der Lösung y(t) keinerlei Rolle mehr.Damit der vonλ2 = −1000 herrührende Bestandteil auch in der numerischen Lösung keineRolle spielt, muss dennoch für dass gewählte Verfahren mit Stabilitätsfunktion F(z):

|F(λ2h)| < 1

gelten, ansonsten sind die Ergebnisse unbrauchbar.

Definition 8.4.8 (steifes Differentialgleichungssystem)

Ein Differentialgleichungssystem heißtsteif, falls in der Lösung abklingende e–Funktionen mitstark unterschiedlichen Abklingkonstantenλk, k = 1, 2, . . . beteiligt sind.Quantitativ: die Zahl

St=max{|Re(λk)|}min{|Re(λk)|}

heißtSteifheitsmaß.Das System heißtsteif, falls das Steifheitsmaß St größer als1000 ist.

Beispiel 8.4.9 In der Lösung des Differentialgleichungssystems:(

y′1y′2

)

=

(−500.5 499.5

499.5 −500.5

)

·(

y1y2

)

+

(50 sin(100t)

3

)

mit Anfangswerten y1(0) = −1 und y2(0) = 2 tauchen die e–Funktionen e−1 · t (alsoλ1 = −1)und e−1000 · t (alsoλ2 = −1000) auf.Löst man dieses System mit dem Euler–Cauchy–Polygonzugverfahren, so muss

• für den “Hauptbestandteil“ e−1 · t der Lösung

|F(λ1 · h)| < 1 also |1 − 1 · h| < 1 also h< 2

gelten.

Damit also der Hauptbestandteil e−1 · t der Lösung qualitativ richtig erfasst wird, reichteine Schrittweite h< 2 aus.

• für den (eigentlich in der Lösung keine große Rolle spielenden) Bestandteil e−1000t muss

|F(λ2h)| < 1 also |1 − 1000h| < 1 also h< 0.002

gelten.

Obwohl dieser Bestandteil e−1000t in der exakten Lösung keine große Rolle spielt, mussim numerischen Verfahren dennoch eine winzig kleine Schrittweite zugrunde gelegt wer-den, damit dieser auch in der numerischen Lösung keine Rollespielt!

Je nachdem, ob und wie diese Bedingung eingehalten wird, sehen die berechneten Lösungen(etwa für y1(t)) aus:



1. Rechnung mit h= 0.0002, die Stabilitätsbedingung ist bei weitem erfüllt:

-3

-2

-1

0

1

2

3

4

5

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16

"h_0002.txt"

so in etwa sieht die exakte Lösung für y1(t) aus!

2. Rechnung mit h= 0.00192, die Stabilitätsbedingung ist so gerade noch erfüllt:

-3

-2

-1

0

1

2

3

4

5

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16

"h_00192.txt""h_0002.txt"

erst nach einiger Zeit klingt der vonλ2 = −1000 herrührende, stark oszilierende Be-standteil der numerischen Lösung ab — die berechnete Lösungnähert sich mehr undmehr der exakten Lösung.

3. Berechnung mit h= 0.00201, die Stabilitätsbedingung ist knapp nicht mehr erfüllt:



-3

-2

-1

0

1

2

3

4

5

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16

"h_00201.txt""h_0002.txt"

der vonλ2 = −1000 herrührende, stark oszilierende Bestandteil der numerischen Lösungschaukelt sich auf!

Die ermittelte Lösung stimmt immer weniger mit der exakten Lösung überein.

Definition 8.4.10 (Stabilitätsgebiet, Stabilitätsintervall)

Sei V ein numerisches Einschrittverfahren mit Stabilitätsfunktion F(z), d.h. bei Anwendung aufdie Testdifferentialgleichung y′ = λy erhält man die Iterationsvorschrift:

yk+1 = F(λh) · yk

• Die Menge komplexer Zahlen:

S := {z∈ C| |F(z)| < 1 }

heißt dasStabilitätsgebietdes Verfahrens V.

Ist bei der Lösung eines Differentialgleichungssystems eine (abklingende) e–Funktion mitAbklingkonstanteλ ∈ C (Re(λ) < 0) beteiligt, muss die Schrittweite so bemessen sein,dassλ · h im Stabilitätsgebiet liegt!

• Die Menge reeller ZahlenI := {x ∈ R| |F(x)| < 1 }

heißt dasStabilitätsintervall des Verfahrens V (es handelt sich meißt wirklich um einIntervall!).

Bemerkung 8.4.11 (Stabilitätsgebiet/Stabilitätsindervall einfacher Verfahren)

1. das Stabilitätsgebiet des Euler–Cauchy–Polygonzugverfahrens ist das Innere des Kreisesum−1 mit Radius1:



Skizze:

1−1−2−3

i

Das Stabilitätsintervall ist(−2, 0).

2. Skizze des Stabilitätsgebietes des klassischen Runge–Kutta–Verfahrens:

1−1−2−3

i

Das Stabilitätsintervall ist(−2.78, 0).

Man kann zeigen:

Bemerkung 8.4.12Für alle expliziten Runge–Kutta–Verfahren ist die Stabilitätsfunktion F(z)ein Polynom in z, das Stabilitätsgebiet ist ein (links des Nullpunktes) liegendes beschränktesGebiet und das Stabilitätsintervall ist ein endliches Intervall der Form(−a, 0) mit a > 0.

Bemerkung 8.4.13 (Stabilitätsfunktion/–Gebiet/–Interval des implizitenEuler–Verfahrens)

Die Anwendung des impliziten Euler–Verfahrens auf die Testdifferentialgleichung y′ = λy führtauf die explizite Gleichung

yk+1 =1

1 − λh · yk

Die Stabilitätsfunktion lautet somit:

F(z) =1

1 − z

Das Stabilitätsgebiet ist das Äußere des Kreises um1 mit Radius1:



1−1−2−3

i

Definition 8.4.14 (A–Stabilität)

Ein Runge–Kutta–Verfahren mit Stabilitätsgebiet S heißtA–stabil oderabsolut stabil, falls

C− = {z∈ C | Re(z) < 0 } ⊂ S

gilt.Bei einem absolut–stabilen Verfahren werden ohne Einschränkung für die Schrittweite alle ab-klingenden e–Funktionen qualitativ richtig berechnet.

Das implizite Euler–Verfahren ist absolut–stabil.

Beispiel 8.4.15Es wird das Testproblem mitλ = −100 und das implizite Euler–Verfahrenzugrundegelegt.Die exakte Lösung y(x) = e−100x ist eine extrem schnell abklingende e–Funktion.Bereits mit der recht großen Schrittweite h= 0.1 erhählt man qualitativ die richtige Lösung:

0

0.2

0.4

0.6

0.8

1

0 0.2 0.4 0.6 0.8 1

"ime_0_1.txt"

Es liegtλ · h = −100 · 0.1 = −10 im Stabilitätsgebiet.Für jedepositive Schrittweite liegtλ · h im Stabilitätsgebiet!

Bemerkung 8.4.16 (Stabilitätsfunktion/–Gebiet/–Intervall der Trapezregel)

Wird die (eindimensionale) Trapezregel



yk+1 = yk +h2

f (xk, yk) +h2

f (xk+1, yk+1)

auf die Testdifferentialgleichung y′ = λy angewendet, so erhält man

yk+1 = yk +λh2

yk +λh2

yk+1

also die (explizite) Gleichung:

yk+1 =1 + λh

2

1 − λh2

· yk

Die Stabilitätsfunktion ist somit

F(z) =1 + z

2

1 − z2

=2 + z2 − z

Die Bedingung|F(z)| < 1

ist für alle z mit Re(z) < 1 erfüllt,

Skizze des Stabilitätsgebietes der Trapezregel:

1−1−2−3

i

Somit folgt: auch die (ebenfalls implizite) Trapezregel ist absolut stabil.

Man kann zeigen:

Bemerkung 8.4.17 (Stabilität von Runge–Kutta–Verfahren)

• Alle impliziten Runge–Kutta–Verfahren sind absolut stabil.

• Ein explizitesRunge–Kutta–Verfahren kann nicht absolut stabil sein!Bemerkung 8.4.18

Für steife Differentialgleichungssysteme darf man kein explizites Runge–Kutta–Verfahrenverwenden!


Numerik gewöhnlicher DifferentialgleichungenGrundlagenEinschrittverfahrenSchrittweitensteuerungStabilität, steife Differentialgleichungen

8 Numerik gewöhnlicher Differentialgleichungen...• das in dieser Gleichung stehende Integral...

Documents

Transcript of 8 Numerik gewöhnlicher Differentialgleichungen...• das in dieser Gleichung stehende Integral...