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 9. Übungsblatt: Gewöhnliche Differentialgleichungen (ODE)

1) Kondition von Anfangswertproblemen

a) Berechnen Sie die analytische Lösung y(t) der beiden folgenden gewöhnlichen Differential-
gleichungen (ODE) mithilfe der Separation der Variablen:
i) ẏ(t) = 2y(t) ii) ẏ(t) = −2y(t)
b) Gegeben seien folgende Anfangswertbedingungen für die jeweiligen ODEs aus a). Berechnen
Sie nun die analytische Lösung der daraus resultierenden Anfangswertprobleme (AWP):
i) ẏ(t) = 2y(t), y(0) = 3, t ≥ 0 ii) ẏ(t) = −2y(t), y(0) = y0 , t ≥ 0.
c) Diskutieren Sie jeweils die Kondition der beiden AWP aus b).

Lösung:
a) i) Die Lösung des ODEs bestimmen wir mittels Separation der Variablen:
dy
= 2y
dt
dy
= dt
2y
Z y
1 Z t
dη = dτ
y0 2η t0
1
(ln |y| − ln |y0 |) = t − t0
2
y(t)
= e2(t−t0 )
y0
y(t) = ±y0 · e2(t−t0 )
y(t0)=y0
y(t) = y0 · e2(t−t0 )

ii) Analog zum ODE (a) ergibt sich hier die Lösung
y(t) = y0 · e−2(t−t0 ) .

b) i) Da für t = 0 = t0 als Anfangswert y(t = 0) = 3 gegeben ist, muss für die AWP
Lösung
y(t) = 3 · e2t
gelten.
 ii) Analog zu i) ergibt sich mit t0 = und y(0) = y0 :
y(t) = y0 · e−2t .

c) Wenden wir uns nun der Betrachtung der (absoluten) Kondition zu:
i) Im Fall von gestörten Anfangswerten
yε (0) = y0 + ε
erhalten wir als Lösung des gestörten AWPs
yε (t) = (y0 + ε) · e2t .
Der Fehler nimmt folglich mit fortschreitender Zeit t zu:
|yε (t) − y(t)| = |ε| · e2t → ∞ für t → ∞
⇒ schlecht-konditioniertes Problem!
ii) y(t) = y0 · e−2t .
Wir erhalten als Lösung des gestörten AWPs
yε (t) = (y0 + ε) · e−2t .
Der Fehler wird gedämpft:
|yε (t) − y(t)| = |ε| · e−2t → 0 für t → ∞
⇒ gut-konditioniertes Problem!

2) Einschrittverfahren
Gegeben sei das folgende Anfangswertproblem:
ẏ(t) = t · y(t), y(0) = 1, t ≥ 0.
a) Berechnen Sie die analytische Lösung y(t) des AWPs mit Hilfe der Separation der Variablen!
b) Berechnen Sie im Intervall [0; 4] numerische Lösungswerte yk , welche die Lösung y(t) in den
Stellen tk approximieren, d.h. yk ≈ y(tk ). Rechnen Sie mit Schrittweite tk+1 − tk = δt = 1
und verwenden Sie die folgenden Verfahren:
i) Explizites Euler-Verfahren:
Bei diesem Verfahren wird lediglich die Steigung im aktuellen Punkt tk zur Berechnung
der numerischen Lösung betrachtet:
tk = t0 + k · δt;
yk+1 = yk + δt · f (tk , yk );

ii) Verfahren von Heun:

Analog zur Trapezregel bei der Quadratur wird bei diesem Verfahren im Vergleich
zum Euler-Verfahren ein zweiter f -Wert zur numerischen Berechnung der Steigung
hinzugezogen:
tk = t0 + k · δt;
δt
yk+1 = yk + · (f (tk , yk ) + f (tk+1 , yk + δt · f (tk , yk ))) ;
2
 iii) Klassisches Runge-Kutta-Verfahren (Zusatz):
Analog zur Fassregel werden hier Zwischenwerte Ti für die Näherung der Steigung
berechnet und mit 1/6 gewichtet:

tk = t0 + k · δt;
T1 = f (tk , yk );
δt δt
T2 = f (tk + , yk + T1 );
2 2
δt δt
T3 = f (tk + , yk + T2 );
2 2
T4 = f (tk+1 , yk + δtT3 );
δt
yk+1 = yk + · (T1 + 2T2 + 2T3 + T4 ) ;
6

Vergleichen Sie die Ergebnisse der einzelnen Verfahren mit der analytischen Lösung. Was
stellen Sie fest und wie lässt sich dies deuten?

Lösung:
a) Wir berechnen die Separation (Trennung) der Variablen für ẏ(t) = t · y(t) = F (t) · G(y)
mit F (t) = t und G(y) = y:

dy
= F (t) · G(y)
dt
dy
= F (t)dt
G(y)
Z y
1 Z t
dη = F (τ )dτ
y(0) G(η) 0
Z y
1 Z t
dη = τ dτ
y(0) η 0
2
t
ln|y| − ln|y(0)| =
| {z } 2
=0
t2
|y(t)| = e 2 , t ∈ [0; b]
t2
y(t) = ±e , 2 t ∈ [0; b]
t2
Wegen y(0) = 1 > 0 gilt y(t) = +e 2 , t > 0.

b) Wir berechnen im Intervall [a; b] = [0; 4] zu den äquidistanten Zeitpunkten tk = t0 + k · δt,

k = 1, 2, 3, 4 eine numerische Approximation.
i) Explizites Euler-Verfahren:

y1 = y0 + δt · f (t0 , y0 ) = 1 + 1 · 0 · 1 = 1,
y2 = y1 + δt · f (t1 , y1 ) = 1 + 1 · 1 · 1 = 2,
y3 = y2 + δt · f (t2 , y2 ) = 2 + 1 · 2 · 2 = 6,
y4 = y3 + δt · f (t3 , y3 ) = 6 + 1 · 3 · 6 = 24.
 t2
Die Werte der analytischen Lösung y(t) = e 2 betragen dagegen
1
y(t1 ) = e 2 = 1.6487 . . . ,
4
y(t2 ) = e 2 = 7.3890 . . . ,
9
y(t3 ) = e 2 = 90.017 . . . ,
16
y(t4 ) = e 2 = 2980.9 . . . .

ii) Verfahren von Heun:

δt 1 3
y1 = y0 + · (t0 · y0 + f (t1 , y0 + δt · t0 · y0 )) = 1 + · (1 · 1) = ,
2 2 2
δt 3 1 3 3 3
y2 = y1 + · (t1 · y1 + f (t2 , y1 + δt · t1 · y1 )) = + · ( + 2 · ( + ))
2 2 2 2 2 2
21
= = 5.25,
4
δt 21 1 21 21 42
y3 = y2 + · (t2 · y2 + f (t3 , y2 + δt · t2 · y2 )) = + · ( + 3 · ( + ))
2 4 2 2 4 4
21 1 189 42 273
= + ·( + ) = = 34.125,
4 2 4 4 8
δt
y4 = y3 + · (t3 · y3 + f (t4 , y3 + δt · t3 · y3 ))
2
273 1 273 273 273
= + · (3 · +4·( +3· ))
8 2 8 8 8
273 1 273 21
= + · (3 · + 2 · 273) = · 273 = 358, 3 . . . .
8 2 8 16

Im Vergleich zu den Werten der analytischen Lösung erkennt man, dass das Ver-
fahren von Heun für diese (relativ großen) Zeitschritte durchaus noch einen nicht
t2
unerheblichen Fehler im Vergleich zur analytischen Lösung y(t) = e 2 macht, aber
schon besser ist als das explizite Euler-Verfahren.
iii) Klassisches Runge-Kutta-Verfahren:
Um einen globalen Runge-Kutta-Schritt zu machen, sind jeweils vier Subschritte
T1 , . . . , T4 auszuführen:
k = 0:

T1 = f (t0 , y0 ) = t0 · y0 = 0,
δt δt 1 1
T2 = f (t0 + , y0 + · T1 ) = ·1= ,
2 2 2 2
δt δt 1 5 5
T3 = f (t0 + , y0 + · T2 ) = · = ,
2 2 2 4 8
13 13
T4 = f (t1 , y0 + δt · T3 ) = 1 · = .
8 8
Damit ergibt sich für den ersten Gesamtschritt:
δt 1 5 13
y1 = y0 + (T1 + 2T2 + 2T3 + T4 ) = 1 + · (0 + 1 + + )
6 6 4 8
1 31 79
= 1+ ·( )= = 1, 645 . . . .
6 8 48
k = 1:

79
T1 = f (t1 , y1 ) = t1 · y1 = ,
48
δt δt 3 79 79 9 79 711 237
T2 = f (t1 + , y1 + · T1 ) = ·( + )= · = = ,
2 2 2 48 2 · 48 4 48 192 64
δt δt 3 79 4 1 711
T3 = f (t1 + , y1 + · T2 ) = ·( · + · )
2 2 2 48 4 2 192
3 2 · 4 · 79 + 711 3 · 1343
= ·( )= ,
2 2 · 192 4 · 192
79 3 · 1343 16 · 79 + 3 · 1343 5293
T4 = f (t2 , y1 + δt · T3 ) = 2 · ( + ) = 2·( )= .
48 4 · 192 4 · 192 2 · 192
Damit ergibt sich für den zweiten Gesamtschritt:
δt 79 1 79 711 3 1343 5293
y2 = y1 + (T1 + 2T2 + 2T3 + T4 ) = + ·( + + · + )
6 48 6 48 2 · 48 2 192 2 · 192
79 1 8 · 79 + 4 · 711 + 3 · 1343 + 5293
= + ·( )
48 6 2 · 192
79 12798 16590
= + = = 7, 200 . . . .
48 2304 48 · 48

Weitere Gesamtschritte ersparen wir uns in Anbetracht der relativ unschönen Brüche
und überlassen diesem dem Computer (z.B. kurzes Matlab-Programm schreiben).
Wir erhalten dann als Ergebnis

y3 = 77.70 . . . ,

y4 = 1856.8 . . .

Das Klassische Runge-Kutta-Verfahren approximiert die analytischen Werte also

noch etwas besser als das explizite Euler-Verfahren und das Verfahren von Heun.

Je größer die Ordnung des Verfahrens (Euler: 1, Heun: 2, kl. Runge-Kutta: 4), desto
besser die numerischen Ergebnisse.
3) Zusatzaufgabe: Euler-Verfahren und Zinsberechnung
Die Verzinsung eines Guthabens, das am Anfang den Wert y0 habe und pro Jahr p% Zinsen
erhält, kann man als Anwendung des expliziten Euler-Verfahrens auf ein Anfangswertproblem
ẏ = f (y) interpretieren.
a) Wiederholung: Wie sieht die Berechnungsvorschrift des expliziten Euler-Verfahrens mit
Schrittweite δt für das AWP ẏ = f (y) aus?
b) Die Diga-Bank schüttet Zinsen immer zum Jahresende aus. Geben Sie die rechte Seite
f (y) an, so dass das Eulerverfahren mit δt=1 (Jahr) gerade diese Verzinsung beschreibt.
c) Berechnen Sie die analytische Lösung des AWPs!
d) Nun macht die Spaßkasse das Angebot, Zinsausschüttungen vierteljährlich statt nur am En-
de des Jahres durchzuführen. Geben Sie die Euler-Formel für diese neue Schrittweite an und
berechnen Sie daraus eine explizite Vorschrift, um direkt den Wert des Gesamtguthabens
nach einem vollen Jahr (yend ) zu ermitteln.
e) Vergleichen Sie für den konkreten Fall eines Startguthabens von y0 =100000 Euro und
einer Verzinsung von p = 2 Prozent p.a. die verschiedenen Guthaben yend nach einem
Jahr bei der Diga-Bank, der Spaßkasse und einer virtuellen Bank mit Verzinsung nach der
analytischen Lösung des AWPs! Was stellen Sie fest?

Lösung:
a) Explizites Euler-Verfahren für das angegebene AWP:

yk+1 = yk + δt · f (yk ) .

b) Ausgehend vom Startkapital y0 soll mit δt = 1 das Startkapital plus Zinsen bei y1
herauskommen. Die Zinsen für den Zeitraum eines Jahres sind δt · p/100 · y0 . Damit ergibt
sich insgesamt:
y1 = y0 + δt · p/100 · y0 ,

und die Funktion f der rechten Seite ist f (y) = p/100 · y.

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

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
p
dη = dτ
y0 η t0 100
p
ln(y) − ln(y0 ) = ·t
100
pt
y(t) = y0 · e 100 , t ∈ [0; b].
d) Wir verwenden nun 4 Schritte des Euler-Verfahrens zur kleineren Schrittweite δt = 1/4
(Jahr). Die Berechnungsvorschrift für yk+1 lautet dann:
!
p δt p
yk+1 = yk + δt · · yk = 1+ · yk .
100 100

Wenn wir diese Definition rekursiv einsetzen, so erhalten wir die explizite Formel zur
direkten Berechnung von yk+1 aus y0 :
! !2 !k+1
δt p δt p δt p
yk+1 = 1+ · yk = 1+ · yk−1 = . . . = 1+ · y0 .
100 100 100

e) Für die konkreten Werte y0 =100000 Euro und p = 2 erhalten wir folgende Gesamtgutha-
benswerte yend nach einem Jahr (b = 1):

Diga-Bank : yend = 105 + 1 · 2/100 · 105 = 102000, 00 ,

1
!4
·2
Spaßkasse : yend = y4 = 1+ 4
· 105 = 102015, 05 ,
100
2
Anal. Lsg. : yend = 105 · e 100 = 102020, 13 .

Man erkennt, dass mit ca. 15 Euro (0,015%) ein nicht unerheblicher Unterschied im Endka-
pital zwischen Diga und Spaßkasse vorliegt. Die Differenz der Spaßkasse zur analytischen
Lösung beträgt nur etwa 5 Euro.