Numerik gewöhnlicherDifferentialgleichungen

MA2304 - SS16

Übungsblatt 11 – Musterlösung

Aufgabe 50 (Stabilität von BDF-Verfahren)

a) Weisen Sie nach, dass das explizite Differenzenverfahren der Ordnung 3

yi+2 + 4yi+1 − 5yi = h(4fi+1 + 2fi)

nicht Null-stabil ist, indem Sie das entsprechende charakteristische Polynom unter-suchen.

b) Ein BDF-Verfahren habe die Formq∑

l=0

alyi+l = hfi+q, q ≤ 6.

Beweisen Sie, dass ein solches Verfahren mit aq 6= 0 in einem unbeschränkten Be-reich absolut stabil ist.Hinweis: Betrachten Sie das charakteristische Polynom Πz(r) = ρ(r)− zσ(r). Wäh-len Sie M so groß, dass für alle z = λh < −M gilt |Πz(r)| > 0.

Lösung 50 (Stabilität von BDF-Verfahren)

a) Das 1. charakteristische Polynom zu diesem Verfahren lautet

ρ(r) = r2 + 4r − 5

und besitzt die Nullstellen r0 = 1, r1 = −5. Da |r1| > 1 ist die Wurzelbedingungnicht erfüllt, also ist das Verfahren nicht Null-stabil.

b) Wir betrachten das Polynom

Πz(r) := ρ(r)− zσ(r) =

q∑

l=0

alrl − zrq = (aq − z)rq +

q−1∑

l=0

alrl.

Wählen wir M > 0 so groß, dass für z = λh < −M die Relation

aq − z >

q−1∑

l=0

|al|

erfüllt ist, dann gilt für jedes r ∈ C mit |r| ≥ 1 die Beziehung

|Πz(r)| ≥ (aq − z)|r|q −q−1∑

l=0

|al||r|q > 0.

Folglich liegen für alle z < −M die Nullstellen von Πz(r) im Einheitskreis, d.h.|ri(z)| < 1, und das Verfahren erfüllt die absolute Wurzelbedingung und ist damitin dem unbeschränkten Bereich z < −M absolut stabil.

1

Aufgabe 51 (Prediktor-Korrektor Verfahren)

Zeigen Sie, dass das Prediktor-Korrektor Verfahren, welches durch Kombination des Adams–Bashforth Verfahren der Ordnung p und des Adams–Moulton Verfahren der Ordnung p+1entsteht, die Ordnung p+1 aufweist. Berechnen Sie hierzu den Konsistenzfehler unter derAnnahme, dass yAB

i−j = yAMi−j für j = 0, 1, 2, . . . und der Lipschitz-Stetigkeit von f(t, y) im

zweiten Argument.

Lösung 51 (Prediktor-Korrektor Verfahren)

Wir haben für den Konsistenzfehler ‖yABMi+1 −y(ti+h)‖ zu betrachten. Sei weiter yAM

i+1 dasErgebnis des Adams–Moulton Verfahren der Ordnung p+ 1. Dann gilt

∥∥yABM

i+1 − y(ti + h)∥∥ ≤

∥∥yABM

i+1 − yAMi+1

∥∥+

∥∥yAM

i+1 − y(ti + h)∥∥

≤∥∥yABM

i+1 − yAMi+1

∥∥+O(hp+2),

da das Adams–Moulton Verfahren die Ordnung p+ 1 besitzt.Weiter betrachten wir nun

yABMi+1 = yi + h

(

b−1 f(ti+1, y

(P ))+

p−1∑

j=0

bjfi−j

)

mit dem Prediktor des AB–Verfahrens

y(P ) = yi + h

p−1∑

j=0

bjfi−j.

Für die Lösung yABMi+1 des AM–Verfahrens gilt

yAMi+1 = yi + h

(

b−1 f(ti+1, y

AMi+1

)+

p−1∑

j=0

bjfi−j

)

.

Dann folgt daraus∥∥yABM

i+1 − yAMi+1

∥∥ = h

∥∥b−1

(f(ti+1, y

(P ))− f

(ti+1, y

AMi+1

))∥∥

≤ h |b−1|L(∥∥y(P ) − y(ti + h)

∥∥

︸ ︷︷ ︸

O(hp+1)

+∥∥yAM

i+1 − y(ti + h)∥∥

︸ ︷︷ ︸

O(hp+2)

)

= O(hp+2).

Daraus folgt schließlich∥∥yABM

i+1 − y(ti + h)∥∥ = O(hp+2).

Aufgabe 52 (BDF-Verfahren im Vergleich)

Betrachten Sie das Anfangswertproblem

y = −αy, y(0) = 1, (1)

mit der exakten Lösung y(t) = e−αt.

2

a) Schreiben Sie ein ODE-File mit der Signatur function f = ODERHS(t,y,flag,var),welches die rechte Seite der ODE (1) berechnet. Der Parameter var enthalte denWert von α. flag muss im Code nicht verwendet werden.

b) Implementieren Sie einen Schritt des BDF-Verfahrens der Schrittzahl q fürq ∈ {1, . . . , 4}. Erstellen Sie dazu eine Routinefunction [tnew,ynew,iter] = BDFall_step_fixpoint(FUNC,told,yold,h,var,q),die ausgehend von den alten Iterierten yold = (yi−q+1, . . . , yi−1, yi) die neue Iterier-te ynew berechnet. Die rechte Seite der ODE soll als FUNC übergeben werden, dieSchrittweite mit h und der aktuelle Zeitschritt mit told.Die Lösung des impliziten Gleichungssystems soll über eine Fixpunktiteration miteinem geeigneten Abbruchkriterium (maxiter = 500, tol = 1e-12) berechnet wer-den. var stellt wie in a) zusätzliche Parameter dar, die an die Funktion FUNC überge-ben werden. Neben der neuen Iterierten (tnew,ynew) = (ti+h, yi+1) soll die Anzahliter der durchgeführten Fixpunktiterationen zurückgegeben werden.Hinweis: Achten Sie darauf, dass das Verfahren auf eine vektorwertige Differenzial-gleichung angewendet werden kann. Zur Implementierung kann der Befehl switchhilfreich sein.

c) Schreiben Sie eine Funktionfunction [tvec,yvec,itervec] = BDFall(FUNC,t0,tN,y0,N,var,q,ANLAUF),die mit Hilfe der in b) implementierten Routine die ODE mit der rechten Seite FUNCnäherungsweise löst. Hierbei gibt (t0,tN) das entsprechende Zeitintervall, y0 denAnfangsvektor und N die Anzahl der Zeitschritte an. Das Funktionshandle ANLAUF

soll für q > 1 zum Berechnen der Anfangswerte y0, . . . , yq−1 verwendet werden, bei-spielsweise durch den Aufruf[tvec(1:q),yvec(:,1:q)] = feval(ANLAUF,FUNC,t0,t0+(q-1)*h,y0,q-1,var).

d) Auf der Vorlesungswebseite finden Sie die Routinen EULER.m, HEUN.m, RK3.m, RK4.msowie das Skript-File BDFall_test.m. Testen Sie Ihre Funktionen aus b) und c),indem Sie BDFall_test aufrufen. Führen Sie die Anlaufrechnung auch mit denanderen drei Einschrittverfahren durch. Wie verhält sich jeweils die Fehlerabnahme?Erklären Sie die Ergebnisse.

Lösung 52 (BDF-Verfahren im Vergleich)

102 103

10-14

10-12

10-10

10-8

10-6

10-4

10-2Fehler gegen Anzahl der Schritte mit Anlauf EULER

BDF1BDF2BDF3BDF4

102 103

10-14

10-12

10-10

10-8

10-6

10-4

10-2Fehler gegen Anzahl der Schritte mit Anlauf HEUN

BDF1BDF2BDF3BDF4

3

102 103

10-14

10-12

10-10

10-8

10-6

10-4

10-2Fehler gegen Anzahl der Schritte mit Anlauf RK3

BDF1BDF2BDF3BDF4

102 103

10-14

10-12

10-10

10-8

10-6

10-4

10-2Fehler gegen Anzahl der Schritte mit Anlauf RK4

BDF1BDF2BDF3BDF4

siehe Matlab files.

4

Aufgaben zum Selbststudium

Aufgabe 53 (Maximale Konsistenzordnung)

Betrachten Sie das Mehrschrittverfahren

yi+1 − 2yi + yi−1 = h2(afi+1 + bfi + afi−1

)

zur Approximation der skalaren Differenzialgleichung y = f(t, y).Welche Konsistenzordnung können Sie für a = 0, b ∈ R bzw. für a, b ∈ R maximalerreichen? Wie lauten die entsprechenden Koeffizienten?

Lösung 53 (Maximale Konsistenzordnung)

Umschreiben des Verfahrens liefert

yi+1 = 2yi − yi−1 + h2(afi+1 + bfi + afi−1

).

Mit der exakten Lösung y der Differenzialgleichung erhalten wir für den lokalen Abbruch-fehler e

e = y(ti + h)−(

2y(ti)− y(ti − h) + h2(ay(ti + h) + by(ti) + ay(ti − h)

))

. (2)

Für y ∈ Cp+2([t0, T ]) gelten die Taylorentwicklungen

y(t− jh) =

p∑

l=0

y(l)(t)

l!(−jh)l +O(hp+1),

y(t− jh) =

p−2∑

l=0

y(l+2)(t)

l!(−jh)l +O(hp−1) =

p∑

k=0

k(k − 1)y(k)(t)

k!(−jh)k−2 +O(hp−1).

Somit haben wir für den lokalen Abbruchfehler e aus (2) mit y(l)i := y(l)(ti) die Darstellung

e = y(0)i + hy

(1)i +

1

2h2y

(2)i +

1

3!h3y

(3)i +

1

4!h4y

(4)i +

1

5!h5y

(5)i +

1

6!h6y

(6)i

− 2y(0)i

+ y(0)i − hy

(1)i +

1

2h2y

(2)i − 1

3!h3y

(3)i +

1

4!h4y

(4)i − 1

5!h5y

(5)i +

1

6!h6y

(6)i

− h2a

(

y(2)i + hy

(3)i +

1

2h2y

(4)i +

1

3!h3y

(5)i +

1

4!h4y

(6)i

)

− h2b(

y(2)i

)

− h2a

(

y(2)i − hy

(3)i +

1

2h2y

(4)i − 1

3!h3y

(5)i +

1

4!h4y

(6)i

)

+O(h7)

= h2(1− 2a− b)y(2)i + h4

(1

12− a

)

y(4)i + h6

(2

6!− 2

4!a

)

y(6)i +O(h7).

Für a = 0 kann durch die Wahl b = 1 die Abschätzung e = O(h4) erreicht werden, so dassdas Verfahren

yi+1 − 2yi + yi−1 = h2f(ti, yi)

5

die Konsistenzordnung 3 besitzt.Damit e = O(h6) gilt, müssen die Gleichungen

1− 2a− b = 0,

1

12− a = 0

erfüllt sein, was a = 112

und b = 56

liefert. Mit a = 112

ist jedoch(26!− 2

4!a)6= 0, und das

Verfahren

yi+1 − 2yi + yi−1 =1

12h2(

f(ti+1, yi+1) + 10f(ti, yi) + f(ti−1, yi−1))

besitzt die Ordnung p = 5.

Aufgabe 54 (Stabilitätsgebiete von MSV)

Berechnen Sie die Stabilitätsgebiete der beiden folgenden impliziten 2-Schrittverfahren

BDF: yi+1 =4

3yi −

1

3yi−1 +

2

3hfi+1,

Adams–Moulton: yi+1 = yi +1

12h (5fi+1 + 8fi − fi−1),

und stellen Sie sie jeweils graphisch dar.

Lösung 54 (Stabilitätsgebiet von MSV)

Für die angegebenen Verfahren erhalten wir die charakteristischen Polynome

ρBDF (r) = r2 − 4

3r +

1

3, σBDF (r) =

2

3r2,

ρAM(r) = r2 − r, σAM(r) =1

12(5r2 + 8r − 1).

Wir betrachten nun das Polynom Πz(r) = ρ(r)− zσ(r). Ein MSV ist absolut stabil, fallsfür alle Nullstellen r(z) von Πz(r) die Relation |r(z)| < 1 gilt. Da die r(z) stetig vonz abhängen, folgt für z ∈ ∂S (S bezeichnet das Stabilitätsgebiet), dass mindestens eineNullstelle r(z) mit |r(z)| = 1 existiert. Mit r(z) = eiφ folgt aus Πz(r) = 0 die Relation

z =ρ(eiφ)

σ(eiφ).

Also haben wir

∂S ⊂ C =

{ρ(eiφ)

σ(eiφ): φ ∈ [0, 2π]

}

.

Für die vorliegenden Verfahren lautet die Parametrisierung der Kurve C

ΓBDF (φ) =3e2iφ − 4eiφ + 1

2e2iφ=

3− 4e−iφ + e−2iφ

2=

(e−iφ − 1)(e−iφ − 3)

2, φ ∈ [0, 2π],

ΓAM(φ) = 12e2iφ − eiφ

5e2iφ + 8eiφ − 1= 12

eiφ(eiφ − 1)

5e2iφ + 8eiφ − 1, φ ∈ [0, 2π].

Nachfolgende Grafiken zeigen die Kurven ΓBDF und ΓAM .

6

−6 −4 −2 0 2 4−4

−3

−2

−1

0

1

2

3

4Stab.gebiet BDF2

−6 −4 −2 0 2 4−4

−3

−2

−1

0

1

2

3

4Stab.gebiet AM2

Wegen ∂S ⊂ Γ folgt, da beide Kurven keine Knoten besitzen, ∂S = Γ. Somit ist jeweilsentweder das Innere der Kurve oder das gesamte äußere Gebiet der Kurve das Stabilitäts-gebiet. Da das Stabilitätsgebiet eines BDF-Vefahrens stets unbeschränkt ist, stellt dasäußere Gebiet der Kurve ΓBDF das Stabilitätsgebiet des vorliegenden BDF-Verfahrensdar. Beim Adams–Moulton–Verfahren erhält man beispielsweise durch die Überprüfungdes Wertes z = −4 für die Nullstellen r1/2(−4) folgendes Ergebnis:

Π−4(r) = r2 − r +1

3(5r2 + 8r − 1) =

1

3

(8r2 + 5r − 1

),

⇒ r1/2(−4) =1

16

(

−5±√25 + 32

)

⇒ |r1/2(−4)| ≤ 13

16< 1.

Damit ist gezeigt, dass das Innere der Kurve ΓAM das Stabilitätsgebiet sein muss.

7