Technische Universität München WiSe 25/26

School of Computation, Information and Technology

Prof. Dr. Hans-Joachim Bungartz
Manuel Geiger
Marc Marot-Lassauzaie

Numerisches Programmieren, Übungen

Musterlösung 10. Übungsblatt: Gewöhnliche Differentialgleichungen (ODE) II

1) Begriffe
Was versteht man unter den folgenden Begriffen:
a) lokaler Diskretisierungsfehler, d) Konsistenz,

b) globaler Diskretisierungsfehler, e) Stabilität,

c) Konvergenz, f) Steifheit?
Machen Sie sich insbesondere den Unterschied zwischen lokalem und globalem Diskretisierungs-
fehler klar. Nutzen Sie hierfür eine Skizze!

Lösung:
a) lokaler Diskretisierungsfehler:
Größter Fehler, der durch die Verwendung des Differenzenquotentienten anstatt der
Ableitung entsteht.
Lokaler Diskretisierungsfehler der Eulermethode:

y(t + δt) − y(t)

( )
l(δt) := max − f (t, y(t))
a≤t≤b−δt δt

Oft wird anstatt dieser eher unhandlichen Definition des lokalen Diskretisierungsfehlers
(der Eulermethode) eine andere (nicht äquivalente) Darstellung verwendet, die in der Praxis
leichter handhabbar ist. Sie verwendet als lokalen Diskretisierungsfehler die Differenz
zwischen analytischer und numerischer Lösung am Ende eines einzelnen Zeitschritts:

|yk+1 − y(tk+1 )|,

wobei yk = y(tk ) gilt.

b) globaler Diskretisierungsfehler:
Größter Fehler zwischen den numerischen Approximationen yk und den korrespondierenden
analystischen Werten y(tk ) zu konkreten Zeitpunkten tk :

e(δt) := max {|yk − y(tk )|}

k=0,...,N

Eine in der Praxis eher gebräuchliche (aber nicht äquivalente) Form des globalen Diskreti-
sierungsfehlers betrachtet lediglich die Differenz zwischen numerischer und analytischer
 Lösung der Differentialgleichung am Ende des betrachteten Zeitintervalls (nach n Zeit-
schritten zum Endzeitpunkt t = b):
|yn − y(tn = b)|
wobei y0 = y(t0 ) gilt.
c) Konvergenz:
Konvergiert die durch ein Lösungsverfahren erzeugte Folge mit immer kleiner werdender
Schrittweite (h → 0) gegen die exakte Lösung, so heißt dieses Lösungsverfahren konvergent.
Es gilt also:
lim yk = y(tk )
h→0,n→∞

Der globale Fehler konvergiert folglich gegen 0.

d) Konsistenz:
Der Begriff der Konsistenz ist eng mit dem lokalen Diskretisierungsfehler verbunden.
Falls für δt → 0 der lokale Fehler l(δt) gegen 0 konvergiert (also l(δt) → 0), so heißt ein
Verfahren konsistent.
e) Stabilität:
• Stabilität ist die Eigenschaft eines Verfahrens/Algorithmus.
• Ein Verfahren ist stabil, wenn es gegenüber kleinen Störungen unempfindlich ist.
• Summieren sich kleine lokale Fehler nur zu kleinen globalen Fehlern auf, so ist ein
Verfahren stabil.
f) Steifheit:
Die Begriffe Konsistenz und Konvergenz sind eher analytischer Natur: Für ausreichend/be-
liebig kleine δt soll etwas passieren. Es gibt jedoch Differentialgleichungen, die konsistent
bzw. konvergent sind, jedoch nur für sehr kleine δt. Ein zu kleines δt kann die Diskretisie-
rung jedoch unpraktikabel machen.
Differentialgleichungen mit dieser Eigenschaft bezeichnet man als steif.
Somit ist Steifheit eine Problemeigenschaft, die bei der numerischen Lösung der Diffe-
rentialgleichung mit expliziten Einschrittverfahren die Verwendung einer sehr kleinen
Schrittweite δt erzwingt, die für eine vorgegebene Genauigkeit in der Approximation
eigentlich nicht nötig wäre.
Für die Beziehung zwischen Stabilität, Konsitenz und Konvergenz gilt bei den expliziten
Verfahren:
• ESV: Konvergenz ⇒ Konsistenz
• MSV: Konsistenz + Stabilität ⇒ Konvergenz

2) Implizites Euler-Verfahren (Rückwärts-Euler)

Gegeben sei das folgende AWP:
ẏ(t) = −12(y(t))2 ∀ t ≥ 1,
y(1) = 1.
Dabei ist bekannt, dass y(t) > 0 für alle t ≥ 1 gilt. Gesucht ist eine Näherungslösung y1 ≈ y(t1 )
zum Zeitpunkt t1 = 1.5.
 a) Lösen Sie das AWP analytisch mit Hilfe der Separation der Variablen und berechnen Sie
die exakte Lösung y(1.5)!
b) Berechnen Sie y1 mit Hilfe des expliziten Eulers.
c) Wenden Sie nun das implizite Euler-Verfahren
yk+1 = yk + hf (tk+1 , yk+1 )
zur Berechnung von y1 an, und vergleichen Sie das Ergebnis mit b).

Lösung:
a) Separation der Variablen:
dy
= −12y 2
dt
1
− dy = 12 dt
y2
Z y
1 Z t
− dη = 12 dτ
y(t0 ) η 2 t0
Z y
1 Z t
− dη = 12 dτ
1 η2 1
1
−1 = 12t − 12
y
1
y(t) =
12t − 11
Damit gilt y(1.5) = 1
7

b) Explizites Euler-Verfahren
y1 = y0 + hf (t0 , y0 )
= 1 + 0.5(−12 · 12 ) = −5

c) Implizites Euler-Verfahren
yk+1 = yk + h · f (yk+1 , tk+1 )
= yk − 12h · yk+1
2

Bemerkung: Die explizite Lösung ist gegeben durch

0 = 12hyk+1
2
+ yk+1 − yk
√
−1 ± 1 + 48hyk
yk+1 =
√24h
y(t)>0 −1 + 1 + 48hyk
=⇒ yk+1 =
24h
Ein Zeitschritt mit dem impliziten Euler-Verfahren:
y0 = 1
y1 = y0 − 12h · y12
√
−1 + 1 + 24 · 1
y1 = = 1
12 3
3) Quadratur und AWP-Lösung
In dieser Aufgabe wollen wir uns noch einmal die Analogie von Quadratur und numerischer
Lösung von Anfangswertproblemen (AWP) verdeutlichen. Wie in der Vorlesung definieren wir
ein AWP durch eine gewöhnliche Differentialgleichung (1) und einen Anfangswert (2)

ẏ(t) = f (t, y(t)) (1)

y(a) = y0 . (2)

Ziel ist die Approximation yi der zugehörigen Lösungsfunktion y(ti ) zu bestimmten Zeitpunkten
ti . Wenn wir Einschrittverfahren verwenden, interessiert uns immer ein Teilintervall [tk ; tk+1 ],
um aus einem alten Wert yk den neuen yk+1 zu berechnen. Dabei können wir den Hauptsatz
der Differential- und Integralrechnung benutzen, der uns folgenden Zusammenhang liefert:
Z tk+1 Z tk+1
(1)
y(tk+1 ) − y(tk ) = ẏ(t) dt = f (t, y(t)) dt . (3)
tk tk

Um nun approximierte Werte yk und yk+1 auf der linken Seite von (3) zu erhalten, integriert
man die rechte Seite numerisch.
a) Benutzen Sie die im Folgenden definierte »dumme Rechtecksregel«, um aus (3) das explizite
Euler-Verfahren herzuleiten!
Die »dumme Rechtecksregel« arbeitet wie die normale Rechtecksregel, nur wird die
Funktion am linken Rand des Intervalls ausgewertet und nicht in der Mitte:
Z b
f (t) dt ≈ (b − a) · f (a) .
a

b) Benutzen Sie die Trapezregel sowie einen zusätzlichen Approximationsschritt, um aus (3)
das Verfahren von Heun herzuleiten!

Lösung:
a) Die »dumme Rechtecksregel«
Z b
f (t) dt ≈ (b − a) · f (a) .
a

liefert direkt

yk+1 − yk = (tk+1 − tk ) ·f (tk , yk )

| {z }
=δt
⇒ yk+1 = yk + δt · f (tk , yk ) .

b) Die Trapezregel lautet

Z b
1
f (t) dt ≈ (b − a) · [f (a) + f (b)] .
a 2
Damit ergibt sich zunächst
1
yk+1 − yk = (tk+1 − tk ) · [f (tk , yk ) + f (tk+1 , yk+1 )]
| {z } 2
=δt
1
⇒ yk+1 = yk + δt · [f (tk , yk ) + f (tk+1 , yk+1 )] . (4)
2
 Allerdings taucht in Glg. (4) nun auch auf der rechten Seite das noch unbekannte yk+1
auf. Um die Lösung einer (i.A.) nichtlinearen Gleichung zu vermeiden, approximiert man
diesen Wert mit Hilfe eines einfach zu berechnenden anderen Verfahrens, wie z. B. des
expliziten Euler-Verfahrens. So erhalten wir:
 
1
yk+1 = yk + δt · f (tk , yk ) + f (tk+1 , yk + δt · f (tk , yk )) . (5)

2 | {z }
≈yk+1

4) Instabilität der Mittelpunktsregel

Neben den bisher betrachteten Einschrittverfahren gibt es eine weitere Klasse der sogenannten
Mehrschrittverfahren. Deren Grundidee ist, die in früheren als dem letzten Schritt bereits
berechneten Approximationen nicht wegzuwerfen sondern mitzunutzen.
Ein einfacher Vertreter der Mehrschrittverfahren ist die Mittelpunktsregel (MPR). Die MPR ist
ein Zweischrittverfahren der folgenden Form:

yk+1 = yk−1 + 2δt · f (tk , yk ), k = 1, . . . , N − 1

Damit ein Mehrschrittverfahren konvergiert, ist zusätzlich zur Konsistenz nun auch eine Stabili-
tätsbedingung notwendig. Wir wollen in dieser Aufgabe an der MPR beobachten, was passiert,
wenn diese Stabilität verletzt wird.
Dazu betrachten wir das AWP

ẏ(t) = −y(t)
y(0) = 1. (6)

Zu Beginn des Verfahrens liegt nur der Anfangswert y0 vor. Somit benötigen wir noch einen
weiteren Wert y1 , bevor die Mittelpunktsregel gestartet werden kann. Typischerweise benutzt
man zur Erzeugung dieses Wertes einfache Einschrittverfahren. Im Folgenden wollen wir daher
y1 mit Hilfe eines Schrittes des expliziten Euler-Verfahrens berechnen.
a) Berechnen Sie die analytische Lösung des AWPs (6)!
b) Wenden Sie die Mittelpunktsregel für t ∈ [0; 6], N = 3 und δt = 2 auf das AWP (6) an!
Was stellen Sie im Vergleich mit der analytischen Lösung fest?
c) Berechnen Sie nun die Ergebnisse mit halber Schrittweite (N = 6 und δt = 1) und
vergleichen Sie sie wieder mit der analytischen Lösung!

Lösung:
a) Die analytische Lösung des AWPs lautet y(t) = e−t . Das erkennt man entweder direkt
aus dem AWP oder rechnet es beispielsweise mit der Separation der Variablen nach
 (ẏ(t) = −y(t) = F (t) · G(y) mit F (t) = 1 und G(y) = −y, sowie y(t0 = 0) = 1):

dy
= F (t) · G(y)
dt
dy
= F (t)dt
G(y)
Z y
1 Z t
dη = F (τ )dτ
y0 G(η) t0
Z y
1 Z t
− dη = 1 dτ
1 η 0
−ln(y) + ln(1) = t
y(t) = e−t , t ∈ [0; b].

b) Für N = 3 und δt = 2 ergibt sich y1 mit einem Startschritt des expliziten Euler-Verfahrens
zu y1 = y0 + δt · f (t0 , y0 ) = 1 + 2(−1) = −1. Damit kann die Mittelpunktsregel als
Zweischrittverfahren starten und ergibt:

y2 = y0 + 2δt · f (t1 , y1 ) = 1 − 2 · 2(−1) = 5,

y3 = y1 + 2δt · f (t2 , y2 ) = −1 − 2 · 2(5) = −21.

Man erkennt, dass die Werte betragsmäßig größer werden, aber das Vorzeichen wech-
selt. Dagegen beschreibt die analytische Lösung ein gleichmäßiges Abklingen gegen Null.
Offensichtlich ist die Stabilitätsbedingung bei der Mittelpunktsregel hier verletzt.
c) Ein zu b) analoges Ergebnis erzielt man auch mit halb so großer Schrittweite δt = 1. Der
Euler-Startschritt liefert y1 = y0 + 1 · (−y0 ) = 0. Damit berechnet man über die Vorschrift
der Mittelpunktsregel:

y2 = y0 + 2δt · f (t1 , y1 ) = 1−2·0=1 ,

y3 = y1 + 2δt · f (t2 , y2 ) = 0 − 2 · 1 = −2 ,
y4 = y2 + 2δt · f (t3 , y3 ) = 1 − 2 · (−2) = 5 ,
y5 = y3 + 2δt · f (t4 , y4 ) = −2 − 2 · 5 = −12 ,
y6 = y4 + 2δt · f (t5 , y5 ) = 5 − 2 · (−12) = 29 .

Die Instabilität lässt sich also nicht durch Verringerung der der Schrittweite beheben,
sondern verschiebt lediglich den Zeitpunkt des Anschwellens der Werte nach hinten. Die
Betrachtung weiterer Verfeinerungslevel mit Hilfe des matlab-Programms vis_midpoint.m
macht dies nochmals deutlich (vgl. Abb. 1).
 visualisation of midpoint rule instability for dy/dt=−y: 5 timestep levels
15

10

0
y

−5

level 0 : dt = 2
−10 level 1 : dt = 1
level 2 : dt = 0.5
level 3 : dt = 0.25
level 4 : dt = 0.125
analytic solution y=exp(−t)
−15
0 1 2 3 4 5 6
t

Abbildung 1: Visualisierung der Instabilität der Mittelpunktsregel: Fünf verschiedene Schritt-

weitenlevel und analytische Lösung.