Explizite gewohnliche Differentialgleichung 1. Ordnung mit

Anfangsbedingung

Gesucht ist eine Funktion y(x), welche erfullt

y′ = f(x, y)

y(x0) = y0

Differentialgleichung

Anfangsbedingung

Wenn f in x stetig ist und einer Lipschitzbedingung genugt, dann

existiert eine eindeutige Losung in der Umgebung des Anfangspunktes

x0.

1

Differentialgleichung gibt Richtungsfeld vor

y′ = xy/4− 1

Drei Losungen zu verschie-

denen Anfangsbedingungen

sind eingetragen 0 1 2 3 4

0

1

2

3

4

2

Eulersches Polygonzugverfahren

Fur die Differentialgleichung und

Anfangsbedingung

y′ = xy/4− 1

y(0) = 3

sind die exakte Losung sowie

drei Naherungen mit Schrittwei-

ten h = 1; 12;14 eingetragen 0 1 2 3 4

0

1

2

3

4

3

Einschrittverfahren: Ablaufschema

Wahle Schrittweite h und maximale Schrittzahl N ;

setze x0 und y0 laut Anfangsbedingung;

fur i = 0,1, . . . , N

xi+1 = xi + h ;

yi+1 = yi + hF(xi, yi, h) .

Die Einschrittverfahren unterscheiden sich in der Wahl der Fortschrei-

tungsrichtung F

Euler-Verfahren: F(x, y, h) = f(x, y),

Modifiziertes Euler-Verfahren: F(x, y, h) = f(

x+ h2, y + h

2f(x, y))

4

Fortschreitungsrichtung F(x, y, h)

5

Modifiziertes Euler-Verfahren

xm

ym

xm+1

ym+1

ym+ ’ym_h2

xm+h2_

F (x, y, h) =

f

(

x+h

2, y +

h

2f(x, y)

)

Verfahren von Heun

xm

ym

xm+1

ym+1

ym+hym’

F (x, y, h) =1

2(k1 + k2) mit

k1 = f(x, y)

k2 = f(x+ h, y + hf(x, y)),

Fortschreitungsrichtung F(x, y, h) beim klassischen Ver-

fahren von Runge-Kutta

F(x, y, h) =1

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

mit

k1 = f(x, y)

k2 = f

(

x+h

2, y +

h

2k1

)

k3 = f

(

x+h

2, y +

h

2k2

)

k4 = f(x+ h, y + hk3).

6

Lokaler Diskretisierungsfehler d(xm, ym, h)

Unterschied zwischen Fort-

schreitungsrichtung F eines

Einschrittverfahrens und ex-

akte Richtung D.

d(xm, ym, h) =

F(xm, ym, h)−D(xm, ym, h) xm

ym

xm+1

ym+1

ym+1

h

D

F

7

Globaler Diskretisierungsfehler

Ist Y die exakte Losung der Anfangswertaufgabe

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

und ym die Naherungslosung an der Stelle xm, so nennt man die

Differenz

e(xm, h) = ym − Y (xm)

den globalen Diskretisierungsfehler.

8

Ordnung eines Einschrittverfahrens

Der lokale Diskretisierungsfehler wird fur h → 0 immer kleiner.

Wie rasch?

Die großte naturliche Zahl p mit

d(xm, ym, h) = O(hp)

heißt Ordnung des Verfahrens.

Interpretation

Ordnung 2 bedeutet, der Fehler nimmt quadratisch in h ab: halbe

Schrittweite viertelt den Fehler

9

Konvergenz des Einschrittverfahrens

Ist der lokale Diskretisierungsfehler von der Ordnung p ≥ 1 und genugt

F einer Lipschitzbedingung, so geht auch der globale Diskretisierungs-

fehler mit dieser Ordnung nach Null: Das Einschrittverfahren ist kon-

vergent von der Ordnung p.

10

Numerische Losung: Wichtige Einschrittverfahren sind

• Eulersche Polygonzugverfahren (weil es das einfachste ist: Ver-

fahren 1. Ordnung)

• Verfahren von Heun, modifiziertes Euler-Verfahren (weil sie ge-

nauer sind: Verfahren 2. Ordnung)

• Klassische Runge-Kutta-Verfahren (weil man damit in der Praxis

oft rechnet; Verfahren 4. Ordnung).

• RK-Verfahren mit der Dormand-Prince-Formel (weil Matlabs ode45

damit rechnet, Ordnung 5 mit Kontrollrechnung 4. Ordnung).

11

Wann soll man welche Methode verwenden ?

12

Gegeben ist die Differentialgleichung

y(t) = −5 y(t)

y(0) = 1.

Das explizite Eulerverfahren und das Heun Ver-

fahren mit Schrittweite 0.41 divergiert. Erst mit

Schrittweite kleiner als 0.4 konvergieren beide.

Das implizite Eulerverfahren konvergiert fur alle

oben angegebenen Schrittweiten.0 1 2 3 4 5 6 7 8

−3

−2

−1

0

1

2

3

t

y(t)

Heun: Schrittweiten 0.41, 2/5,0.39,0.3

0 1 2 3 4 5 6 7 8−3

−2

−1

0

1

2

3

t

y(t)

Implizit: Schrittweiten 0.41, 2/5,0.39,0.3

0 1 2 3 4 5 6 7 8−3

−2

−1

0

1

2

3

t

y(t)

Explizit: Schrittweiten 0.41, 2/5,0.39,0.3

13

Wann soll man welche Methode verwenden ?

Stabilitat:

Es gibt verschiedene Definitionen fur Stabilitat eines numerischen Verfahrens. All-gemein kann man sich merken, dass bei einen impliziten Verfahren die Schrittweiteunabhangig von den Problem gewahlt werden kann, bei expliziten Verfahren dieSchrittweite den Problem angepasst werden muss.

• Explizites Euler (1− 5× 0.41)k y0 divergiert (und oszilliert) fur k → ∞;

• Heun Euler(

1− 5× 0.41+ 25× 0.412)k

y0 divergiert fur k → ∞;

• Explizites Euler(

[1 + 5× 0.41]−1)k

y0 konvergiert fur k → ∞.

14

• Robertson ODE

Gegeben ist die Differentialgleichung y(0) = (1,0,0)und

y1(t) = −0.04 y1(t)

+104y2(t) y3(t),

y2(t) = 0.04 y1(t)− 104y2(t) y3(t)

−3× 107y22(t),

y3(t) = 3× 107y22(t).

Wie man an den Plot sieht, funktioniert der ode45

schlechter als der ode15s, obwohl der ode45 hohereOrdnung hat !

0 1 2 30

0.5

1

1.5

2

2.5

3

3.5

4x 10

−5

t

y 2(t)

ode45

0 1 2 30

0.5

1

1.5

2

2.5

3

3.5

4x 10

−5

t

y 2(t)

ode15s

15

Wann soll man welche Methode verwenden ?

• Steifheit einer Differenzialgleichung:

Sei y(t) = f(y(t)), t ≥ 0, y(0) = y0 gegeben.

Sei Jf(z) = ∂f∂y

∣

∣

∣

∣

y=z.

Dann ist eine mogliche Definition der Steifigkeit

S :=maxj |λj|

minj |λj|

wobei λ1, . . ., λn die Eigenwerte der Matrix Jf(z) ist.

16

Wann soll man welche Methode verwenden ?

• Robertson ODE

Jf(z)

∣

∣

∣

∣

z=y=

y1 104 y3 104y20.04 −104 y3 − 6× 107y2 −104y2 − 6× 107y20 6× 107y2 0

an der Stelle (1,0.1,0.1) erhalt man als Eigenwerte −6, 0.001, 0.00001.

17

Moderne Runge-Kutta-Verfahren

Klassisches RK-Verfahren wertet f(x, y) viermal pro Schritt aus:

f(x, y), f

(

x+h

2, y +

h

2k1

)

, f

(

x+h

2, y +

h

2k2

)

, f(x+ h, y + hk3).

Neuere Verfahren werten f an speziell gunstigen Zwischenstellen aus

und liefern gleichzeitig zwei Werte mit unterschiedlicher Fehlerord-

nung (Differenz ≈ Fehler).

Das Verfahren RK5(4) von Dorman und Prince (MATLAB: ode45)

wertet f sechsmal aus und liefert Ergebnis mit Fehlerordnung 5, ver-

wendet Ergebnis mit Fehlerordnung 4 zur Differenzbildung und Feh-

lerabschatzung

18

Ein- und Mehrschrittverfahren

Runge-Kutta-Verfahren sind Einschritt-Verfahren: um y(x+h) zu berechnen, brau-chen sie die Losung nur am unmittelbar vorhergehenden Punkt y(x).

Mehrschritt-Verfahren verwenden zur Berechnung von y(x+h) die Werte von meh-reren zuruckliegenden Punkten y(x), y(x− h), y(x− 2h) . . ..

Beispiel: Adams-Bashforth-Moulton-Verfahren. Eine Variante davon ist als ode113

in MATLAB verfugbar.

Vorteil: hohe Genauigkeit im Verhaltnis zum Rechenaufwand, besonders bei”teurer“

Auswertung von f .

Nachteil: Braucht Anlaufphase. Nicht einfach bei variabler Schrittweite.

19

Fehlerkontrolle, Schrittweitensteuerung

Eine Schatzung des tatsachlichen Fehlers ǫ1 (absolut oder relativ) bei Schrittweiteh1 sei bekannt. Schrittweite und Fehler stehen bei Fehlerordnung p im Verhaltnis

ǫ2

ǫ1=

(

h2

h1

)p

Um eine gewunschte Fehlerschranke ǫ2 zu erreichen: Andere Schrittweite gemaß

h2 = h1

(

ǫ2

ǫ1

)1

p

Steuerung in Matlab: Schranken fur relativen und absoluten Fehleroptions=odeset(’RelTol’,1.e-7,’AbsTol’,1.e-10)

20

System gewohnlicher Differentialgleichungen 1. Ordnung,

Anfangswertproblem

Vektorielle Schreibweise (im Skriptum komponentenweise ausgefuhrt) zeigt die Ana-logie zum Anfangswertproblem fur eine explizite DG 1. Ordnung

Gesucht ist eine vektorwertige Funktion y(x), welche erfullt

y′ = f(x,y)

y(x0) = y0

Einschrittverfahren (Eulerscher Polygonzug, Heun, etc. ) lassen sich direkt verall-gemeinern.Numer. Losung im MATLAB mit ode45 analog zu einer DG.

21

Gewohnliche Differentialgleichung hoherer Ordnung

Eine DG hoherer Ordnung (oder ein System solcher DG) lasst sich

durch Einfuhren von Hilfsfunktionen in ein aquivalentes System von

DGs 1. Ordnung transformieren.

Beispiel: Mathematisches Pendel, y = − sin(y)

setzey = y1

y = y2−→ neues Glsyst.

y1 = y2y2 = − sin(y1)

22

Beispiel: die Blasius-Gleichung y′′′ = −yy′′/2

beschreibt Stromung in laminarer Grenzschicht

setze

y = y1

y′ = y2y′′ = y3

−→ neues Glsyst.

y′1 = y2

y′2 = y3y′3 = −y1y3/2

Weitere Beispiele in den Ubungsunterlagen!

23

Allgemein: Umformen von y(d) = f(x, y, y′, . . . , y(d−1))

Man setzt y1 = y, y2 = y′, . . . , yd = y(d−1) und schreibt

y′1 = y2

y′2 = y3...

y′d−1 = yd

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

24