Post on 04-Sep-2019
transcript
Kapitel 2
Gewohnliche Differentialgleichungen
2.1 Vorbemerkungen
Betrachte skalare Gleichung
y′ =dy
dx= f(x, y(x)) (2.1)
• y(x) ist die gesuchte Lsg.-Funktion der unabhangigen Variablen x
• Mit der Anfangsbed. y(x0) = y0
• Anfangswertproblem: Lsg. heißt Integralkurve
Richtungsfeld: Trage an Pkt. (x, y) lokal die Steigungen auf. −→ Eindruck uber denVerlauf der Losung (Abb. 2.1). Durch die Wahl von x0 wird eine bestimmte Losungausgewahlt.
Abbildung 2.1: Richtungsfeld einer Differentialgleichung.
Existenz und Eindeutigkeit einer Losung: Satz von Picard-Lindelof.
13
14 KAPITEL 2. GEWOHNLICHE DIFFERENTIALGLEICHUNGEN
2.2 Einfache Einschrittverfahren
Idee: Schreibe dx und dy als endliche Intervalle ∆y und ∆x. “Schreite” vorwarts mit∆x = h
dy
dx=
∆y
∆x= f −→ yk+1 − yk = ∆xf
Allgemein konnen Einschrittverfahren geschrieben werden als
yk+1 = yk + hΦ(xk, yk; yk+1, h) (2.2)
• gehe von xk nach xk+1 und yk nach yk+1
• falls auf rhs.Φ = Φ(xk, yk;h)
explizites Verfahren, sonst implizit, falls auch von yk+1 abhangig.
• Φ ist die Inkrement-Funktion (Steigung)
2.2.1 Taylorreihen-Methoden
Einschrittverfahren konnen durch Taylor-Reihenentwicklung konstruiert werden. Die ge-suchte Losungsfunktion y(x) wird um den Startwert x0 entwickelt
y(x) = y(x0) + (x− x0)y′(x0) +(x− x0)2
2y′′(x0) + ...+
(x− x0)p
p!y(p)(x0) +Rp+1 (2.3)
Hier ist Rp+1 eine Restglied von O((x−x0)p+1). Werde Rp+1 vernachlassigt, und betrachteder Schritt von xk nach xk+1 mit h = xk+1 − xk und xk = x0.
yk+1 = yk + hy′k +h2
2y′′k + ...+
hp
p!y(p)k (2.4)
wobei y(p)k = y(p)(xk) die p-te Ableitung im Punkt xk ist.
BeispieleMit p = 1 folgt die Euler-Methode
yk+1 = yk + hf(xk, yk) (2.5)
als einfachstes Einschrittverfahren, manchmal als Euler-Cauchy bezeichnet (vgl. Fig. 2.2).Wie die Abb. 2.2 zeigt, ist dies Verfahren recht grob und nur fur kleine h einigermaßengeeignet. Die Ableitung wird hierbei nur an einem Punkt ausgewertet und extrapoliert.Das Verfahren ist von 1.Ordnung (s.u).Hohere Ordnungen lassen sich durch Verwendung hoherer Ableitungen gewinnen. Mit
y′ = f(x, y)
y′′ = fx + fyy′ = fx + ffy
folgt
yk+1 = yk + hf(xk, yk) +h2
2[fx(xk, yk) + f(xk, yk)fy(xk, yk)] (2.6)
und analog fur hohere Ordnungen.
c© W. Kley; Skript: Computational Astrophysics
2.2. EINFACHE EINSCHRITTVERFAHREN 15
Abbildung 2.2: Die Euler Methode (Polygonzugmethode) an einem Beispiel illustriert.
2.2.2 Fehler
a) Lokaler Diskretisierungfehleran der Stelle xk+1 mit exakter Losung y(xk) und y(xk+1) und numerischer yk.
dk+1 := y(xk+1)− yk+1
= y(xk+1)− y(xk)− Φ(xk, yk; yk+1, h) (2.7)
(2.8)
Der Lokaler Diskretisierungfehler ist also die Abweichung wahre - numerische Losungin einem Integrationsschritt.
b) Globaler Diskretisierungfehler
gk := y(xn)− yk (2.9)
nach k Schritten: x0 → x1 → x2, ....
Bei Euler
dk+1 := y(xk+1)− y(xk)− hy′(xk)
=h2
2y′′(xk + Θh) mit 0 < Θ < 1 (2.10)
d.h. das Euler-Verfahren ist von 1. Ordnung, das Restglied ist ∝ h2 oder O(h2), also
yk+1 = yk + hf(xk, yk) +O(h2)
Bemerkungen:
c© W. Kley; Skript: Computational Astrophysics
16 KAPITEL 2. GEWOHNLICHE DIFFERENTIALGLEICHUNGEN
1) AllgemeinesEin Einschrittverfahren hat die Fehlerordnung p, falls fur den lokalen Diskretisie-rungfehler gilt:
max1≤k≤n
|dk| ≤ D = const.× hp+1 = O(hp+1)
wobei D das Maximum der dk ist.
2) Einschrittverfahren:gehe von einem Schritt xk nach xk+1
3) Gobaler Fehler:ist (bei expliziten Methoden) gegeben durch
|gn| ≤D
hL
(enhL − 1
)≤ D
hLenhL (2.11)
wobei der maximale Betrag der lokalen Fehler dk ist. L Lipschitzkonstante.
Φ(x, y, h)− Φ(x, y∗, h) ≤ L|y − y∗|
mit x ∈ [a, b], y, y∗ ∈ [−∞,∞], h ∈ [0, h0].
4) Konsistenz:Einschrittverfahren sind konsistent (mit der DGL), falls die Fehlerordnung minde-stens 1 ist. (d.h. die numerische Losung konvergiert fur h → 0 gegen die ’wahre’Losung).
2.2.3 Methode von Heun
Das Richtungfeld andert sich von Punkt zu Punkt, daraus folgt eine einfache Verbesserung
a) Naherungswert (Schatzwert) mit Euler
yk+1 = yk + hf(xk, yk) (2.12)
b) Nachbesserung
yk+1 = yk +h
2[f(xk, yk) + f(xk+1, yk+1)] (2.13)
Die ist ein Beispiel fur eine einfaches Pradiktor-Korrektor Verfahren. Die Steigungen bei(xk, yk) und (xk+1, yk+1) werden arithmetisch gemittelt.Umgeschrieben lautet Heun:
K1 = f(xk, yk) (2.14)
K2 = f(xk + h, yk + hK1) (2.15)
yk+1 = yk +h
2(K1 +K2) (2.16)
(2.17)
Bei Heun ist
dk ≈h3
3!C +O(h2) (2.18)
d.h. die Ferhlerordnung ist 2.
c© W. Kley; Skript: Computational Astrophysics
2.3. RUNGE-KUTTA VERFAHREN 17
2.3 Runge-Kutta Verfahren
2.3.1 Methodik
Wandle die Differentialgleichung in eine Integralgleichung um.∫ xk+1
xk
dy
dxdx =
∫ xk+1
xk
fdx
liefert
y(xk+1) = y(xk) +
∫ xk+1
xk
f(x, y(x)) dx (2.19)
Approximiere das Integral durch Quadraturformel mit Gewichten Ai und Stutzstellenξi ∈ [xk, xk+1]. Sei h ≡ xk+1 − xk∫ xk+1
xk
f(x, y(x)) dx ≈ hm∑i=1
Ai f(ξi, y(ξi)) (2.20)
Hier sind Ai und ξi vorerst beliebig mit Ai ≥ 0. Aus dem Spezialfall
f = const. →m∑i=1
Ai = 1
Sei nun
ξi = xk + αi h 0 ≤ αi ≤ 1, α1 = 0
dann
y(xk + h) = y(xk) + hm∑i=1
Ai f(xk + αih, y(xk + αih)) (2.21)
Nun sind die y(xk + αih) unbekannt, falls αi 6= 0.Benutze wieder Integration
y(xk + αih) = y(xk) +
∫ αi
0
y′(xk + τh)︸ ︷︷ ︸f(x,y)
dτ i = 1, ...,m
wende erneut Quadraturformel an mit gleichen Knoten αi, d.h.∫ αi
0
g(w)dw ≈m∑j=1
βij g(αi) i = 1, ...,m
Aus
g = const. → αi =m∑j=1
βij
ansonsten sind die βij noch frei. Insgesamt folgt
yk+1 = y(xk + h) = yk + h
m∑i=1
AiKi (2.22)
c© W. Kley; Skript: Computational Astrophysics
18 KAPITEL 2. GEWOHNLICHE DIFFERENTIALGLEICHUNGEN
Ki = f
(xk + αih, yk + h
m∑j=1
βijKj
)(2.23)
Diese beiden Gleichungen bilden i.A. ein System von m nicht-linearen Gleichungen. DieGl. (2.22) und (2.23) definieren ein m-stufiges Runge-Kutta Verfahren.Zu den Stutzstellen: Man wird α1 = 0 wahlen, sodass xk die erste Stutzstelle ist.Bspl.
1) Euler-Methodem = 1, α1 = 0, β11 = 0, A1 = 1
Gl. (2.23): → K1 = f(xk, y(xk))Gl. (2.22): → yk+1 = yk + hK1
Dies ist das Euler-Verfahren mit einer Funktionsauswertung.
2) Heun-Methodem = 2, α1 = 0, β11 = β12 = 0,
α2 = 1, β21 = 1, β22 = 0, A1 = A2 = 1/2
Es folgt
K1 = f(xk, yk) (2.24)
K2 = f(xk + h, yk + hK1) (2.25)
yk+1 = yk +h
2
f(xk, yk)︸ ︷︷ ︸K1
+ f(xk + h, yk + hK1)︸ ︷︷ ︸K2
(2.26)
2 Funktionsauswertungen, (auch verbesserte Euler-Methode):
3) Mittelpunktsmethode
m = 2, α1 = 0, β11 = β12 = 0,
α2 = 1/2, β21 = 1/2, β22 = 0, A1 = 0, A2 = 1
Es folgt
K1 = f(xk, yk) (2.27)
K2 = f(xk +1
2h, yk +
h
2K1) (2.28)
yk+1 = yk + hK2 (2.29)
oder
yk+1 = yk + h f
(xk +
h
2, yk +
h
2f(xk, yk)
)(2.30)
Modifiziertes Euler-Verfahren, verbesserte Polygonzugmethode.
Bemerkungen
a) Verbessertes Euler-Verfahren, auch Heun-Verfahren (s.o.), Pradiktor-Korrektor Ver-fahren. Explizites Euler-Verfahren als Pradiktor, implizite Trapezmethode als Kor-rektor.
c© W. Kley; Skript: Computational Astrophysics
2.3. RUNGE-KUTTA VERFAHREN 19
b) Zur Trapezregel:
yk+1 = yk +h
2[f(xk, yk) + f(xk+1, yk+1)]
Dies ist eine implizite Integrationsmethode. Jeder Integrationsschritt verlangt dieLosung einer impliziten Gleichung fur f(xk+1, yk+1). (f ist i.A. nicht-linear). Losungdurch Fix-Punkt-Iteration
y(0)k+1 = yk + hf(xk, yk) (2.31)
y(n+1)k+1 = yk +
h
2
[f(xk, yk) + f(xk+1, y
(n)k+1)
]n = 0, 1, ... (2.32)
(2.33)
c) Es gibt zu jedem m ≥ 2 beliebig viele RK-Verfahren.
2.3.2 Anwendungen
1) Klassische Runge-Kutta-Verfahren durch Butcher-Schema
αi βij ...α1 β11 β12 ... β1mα2 β21 ... ... β2m
∑βij = αi
α3...
......
...αm βm1 ... ... βmm
A1 A2 ... Am∑Ai = 1
(2.34)
wobeim∑i=1
Ai = 1
undm∑j=1
βij = αi
Die Verfahren sind explizit, falls βij = 0, fur j ≥ i. Dann K1 = f(x, y), und K2
explizit aus K1 usw.
2) BeispieleEuler
α β0 0
1(2.35)
Heunαi βij0 0 01 1 0
1/2 1/2
(2.36)
c© W. Kley; Skript: Computational Astrophysics
20 KAPITEL 2. GEWOHNLICHE DIFFERENTIALGLEICHUNGEN
Modifiziertes Euler-Verfahren
αi βij1/2 0 01/2 1/2 0
0 1
(2.37)
3) RK4
αi βij ...0 0 0 0 01/2 1/2 0 0 01/2 0 1/2 0 01 0 0 1 0
1/6 1/3 1/3 1/6
(2.38)
Ausgeschrieben
K1 = f(xk, yk) (2.39)
K2 = f(xk +1
2h, yk +
h
2K1) (2.40)
K3 = f(xk +1
2h, yk +
h
2K2) (2.41)
K4 = f(xk+1, yk + hK3) (2.42)
Insgesamt
yk+1 = yk +h
6(K1 + 2K2 + 2K3 +K4) (2.43)
Abbildung 2.3: Skizzierte Lage der Schatzwerte bei RK4 (nach N.R.)
Brauche 4 Fkt. Auswertungen: 1 zu Beginn, zwei Trials in der Mitte, eine am En-de, siehe Abb. 2.3. Die RK-Methode ist in jedem Schritt gleich, jeder Pkt. konnteAnfangspkt sein, brauche keine vorherigen Pkt., d.h. “einfache Routine”.
4) Fehlerordnungen
Euler 1Heun 2Verbessertes Polygon 2RK4 4
c© W. Kley; Skript: Computational Astrophysics
2.4. SYSTEME VON DGL. ERSTER ORDNUNG 21
Max. Fehlerordnung p von RK-Verfahren
m 1 2 3 4 5 6 7p 1 2 3 4 4 5 6
Also ist RK4 besonders beliebt.
5) ZitateFor many scientific users, fourth order RK is not just the first word on ODE inte-grators, but the last word as well.
Altes Arbeitspferd, aber sehr solide (bei adaptivem Gebrauch). Bulirsch Stoer, Pradiktor-Korrektor, Rennpferde (empfindlich)
- RK nimmt man, wenn man es nicht besser weiß
- Probleme mit besseren Methoden
- Effizienz kein Problem
- RK4 funktioniert fast immer (mit adapt.).
6) Rechenbeispielsei
y′ = y2, y(0) = −4 0 ≤ x ≤ 0.3
Wahle h = 0.1, integriere in 3 Schritten. Die exakte Losung ist
y = − 1
x+ 1/4
Integration ergibt
xk exakt Euler Heun RK40.0 -4.000000 -4.000000 -4.000000 -4.000000
0.1 -2.857143 -2.400000 -2.912000 -2.857341
0.2 -2.222222 -1.824000 -2.275003 -2.222404
0.3 -1.818182 -1.491302 -1.861791 -1.818323
2.4 Systeme von Dgl. erster Ordnung
DGL hoherer Ordnung lassen sich auf Systeme erster Ordnung zuruckfuhrenBspl:
y′′ + q(x)y′ = r(x)
y′ = z(x)
Die zweite Gl. abgeleitet und die erste eingesetzt ergibt
z′ = r(x)− q(x)z(x)
D.h. habe Systemyi(x)
dx= fi(x, y1, ..., yN) i = 1, ..., N
oder als Vektorgleichungy′ = f(x,y(x)), y(x0) = y0
c© W. Kley; Skript: Computational Astrophysics
22 KAPITEL 2. GEWOHNLICHE DIFFERENTIALGLEICHUNGEN
Startvektor y0. Ublich ist hier
y =
y1y1......yN
f =
y2y3...yN−1f(x, y1, ..., yN)
y0 =
y0y(1)0
...
...
y(N−1)0
(2.44)
Bisherige (explizite) Methoden lassen sich direkt auf ein System umschreiben.Bspl.
1. Euler
K1 = f(xk,yk) (2.45)
yk+1 = yk + hK1 (2.46)
2. Heun
K1 = f(xk,yk) (2.47)
K2 = f(xk + h,yk + hK1) (2.48)
yk+1 = yk +h
2(K1 + 2K2) (2.49)
3. RK4
K1 = f(xk,yk) (2.50)
K2 = f(xk +1
2h,yk +
h
2K1) (2.51)
K3 = f(xk +1
2h,yk +
h
2K2) (2.52)
K4 = f(xk+1,yk + hK3) (2.53)
yk+1 = yk +h
6(K1 + 2K2 + 2K3 + K4) (2.54)
wobei yk,yk+1,K1,K2,K3,K4 und f Vektoren im Rn sind.
Bem: Der Fehler ist jetzt ein Vektor dk, brauche geeignete Norm (euklidisch).
2.5 Adaptive Schrittweiten-Kontrolle
Die Idee zur Verbesserung ist:
• Viele kleine Schritte in “gefahrlichem” Terrain
• Große in uninteressanter Landschaft.
• Die Genauigkeit kann z.B. mit Hilfe von Erhaltungsgroßen gemessen werden.
• Es ist eine Abschatzung des Diskretisierungsfehlers notig
c© W. Kley; Skript: Computational Astrophysics
2.5. ADAPTIVE SCHRITTWEITEN-KONTROLLE 23
Abbildung 2.4: Schrittgroße beim adaptive RK4-Verfahren (aus N.R.).
Wir wenden hier ein adaptives Verfahren auf RK4 an:Zunachst Schrittverdopplung:unabhangig, zwei Schritte mit h, einen mit 2h, brauche 11 Funktionsauswertungen (imVergleich zu 8), d.h. ein Overhead von 1.375.Also ist
y(x+ 2h) = y1 + (2h)5Φ +O(h6) (2.55)
y(x+ 2h) = y2 + 2(h)5Φ +O(h6) (2.56)
mit Φ ≈ y(5)(x)/5!. Hier ist y2 der Schatzwert fur die kleine, und y1 fur die große Schritt-weite. Ein guter Fehlerschatzer ist nun
∆ = y1 − y2 = 30h5 Φ +O(h6)
und
y(x+ 2h) = y2 +∆
15+ O(h6)
Dies ist eine Ordnung hoher als RK4, aber der Abschneidefehler ist nicht bekannt. Nehmetrotzdem Fehler ∆ und auch neue O(h6) Schatzung. (Aber merke: Hohere Ordnung nichtimmer hohere Genauigkeit.)
∆ ∝ h5
Habe Schritt h1 mit Fehler ∆1 durchgefuhrt. Es folgt, dass fur einen Schritt h0 mit Fehler∆0 gilt
h0 = h1
(∆0
∆1
)0.2
(2.57)
Hier sei ∆0 die gewunschte Genauigkeit. ∆ ist hier ein lokaler Fehler. Der schlimmste Falltritt ein, falls alle dk ≈ ∆ > 0, dann ergibt sich ein Aufsummieren der Fehler. Muss also beikleinerem h auch kleineres ∆0 wahlen (weil mehr Schritte benotigt werden). Ansonstenist die Wahl von ∆0 nicht bestimmt, richtet sich nach Bedarf (absoluter Fehler, nachSteigung ...). Eine gute Wahl (N.R.) ist
∆0 ≈ ε h
(dy
dx
)i
(2.58)
Die relative Genauigkeit fur Steigung von y, nicht y selbst (ist restriktiver).
c© W. Kley; Skript: Computational Astrophysics
24 KAPITEL 2. GEWOHNLICHE DIFFERENTIALGLEICHUNGEN
Es folgt
h0 = S h1
(∆0
∆1
)0.2
∆0 ≥ ∆1
bei einer Vergroßerung der Schrittweite, und
h0 = S h1
(∆0
∆1
)0.25
∆0 < ∆1
bei einer Verkleinerung der Schrittweite, wobei S ein Sicherheitsfaktor (z.B. 0.9) ist. DiePotenz der ∆ wird deswegen geandert, weil bei obiger Definition (2.58) ein Faktor hauftritt.- gut zu gebrauchen bei nicht glatten Funktionen (robust)- Exzellent durch steiniges, diskontinuierliches Terrain- Vorsicht bei inneren singularen PunktenEin adaptives Verfahren ist generell empfohlen (z.B. ODEINT in Numerical Recipes).
2.6 Stabilitat
Benutze Eigenschaften der DGL (rhs) zur Wahl des Losungsalgorithmus
2.6.1 Inharente Instabilitat
Betrachte Beispiel:
y′(x) = λ [y(x)− F (x)] + F ′(x) mit y(x0) = y0 (2.59)
Homogene Gleichungy′(x) = λy
Die allgemeine Losung der hom. Gl. lautet
yhom(x) = Ceλx
eine partikulare Lsg. der inhomogenen
ypart(x) = F (x)
es folgt die Gesamtlosung
y(x) = [y0 − F (x0)] eλ(x−x0) + F (x) (2.60)
Speziell fur y0 = F (x0) ist y(x) = F (x), aber fur die kleine Abweichung
y0 = F (x0) + ε
wird die Losungy(x) = εeλ(x−x0) + F (x)
Ist also λ > 0 so divergiert die Losung exponentiell von der ersten fur beliebig kleineAbweichungen ε.
c© W. Kley; Skript: Computational Astrophysics
2.6. STABILITAT 25
- starke Empfindlichkeit von den Anfangsbedingungen.- schlecht konditioniert
Ist Beispiel fur inharente Instab.
Methode mit hoher Genauigkeit (hohe Fehlerordnung, hohe Rechengenauigkeit).
→ Rundungs-, Diskretisierungsfehler klein.
Bspl.:
y′(x) = 10
y(x)− x2
1 + x2︸ ︷︷ ︸F (x)
+2x
(1 + x2)2︸ ︷︷ ︸F ′(x)
mit y(0) = y0 = 0 hat die Losung
y(x) =x2
1 + x2
ist also vom obigen allgemeinen Typ (2.59). Die Numerik (RK4 mit h = 0.01) zeigt großereAbweichungen aufgrund der abklingenden exponentiellen Lsg., vgl. Abb. 2.5.
Numerische
Loesung
analytisch
Abbildung 2.5: Beispiel fur eine inharente Instabilitat (aus Schwarz.)
2.6.2 Absolute Stabilitat
Der Begriff der absoluten Stabilitat wird anhand der folgenden Testgleichung definiert.
y′(x) = λy(x) (2.61)
Die Losung lautet eλx
c© W. Kley; Skript: Computational Astrophysics
26 KAPITEL 2. GEWOHNLICHE DIFFERENTIALGLEICHUNGEN
Betrachte RK4
K1 = λ yk (2.62)
K2 = λ (yk +1
2hK1) = (λ+
1
2hλ2)yk (2.63)
K3 = λ (yk +1
2hK2) = (λ+
1
2h(λ+
1
2hλ2))yk (2.64)
K4 = λ (yk + hK3) = ... (2.65)
yk+1 = yk +h
6(K1 + 2K2 + 2K3 +K4) (2.66)
= (1 + hλ+1
2h2λ2 +
1
6h3λ3 +
1
24h4λ4)yk
= F (hλ) yk (2.67)
alsoyk+1 = F (hλ) yk (2.68)
hier ist F (hλ) der Beginn der Taylorreihe fur ehλ. Fur die exakte Losung gilt y(xk+1) =ehλy(xk), also ist der Abschneidefehler dk = O(h5). Fur ein RK Verfahren der Ordnung pentspricht F (λh) fur Testbeispiel (2.61) der Summe der p+ 1 ersten Taylor-Terme.Wir unterscheiden zwei Falle:
a) reelle Werte von λλ > 0 :
→ hλ > 0→ F (hλ) > 1
d.h. yk werden (qualitativ) richtig berechnet.
λ < 0 : Die Losung klingt also ab (z.B. bei Systemen von DGL. eine Komponente),aber die numerische Losung klingt nur ab, falls
|F (hλ)| < 1
aber weil F (µ) ein Polynom 4. Grades ist, gilt
limhλ→−∞
F (hλ) = +∞ → nicht konvergent fur alle h, −→ problems
b) Komplexe Werte von λSysteme von DGL haben oft oszillierende Komponenten mit λ ∈ C. Brauche Bedin-gung |F (hλ)| < 1 (notwendig und hinreichend)
Def.:Fur ein Einschrittverfahren, das fur die Testaufgabe (2.61) auf
yk+1 = F (hλ) yk
fuhrt, heißt die MengeB := {µ ∈ C | |F (hλ)| < 1}
das Gebiet der absoluten Stabilitat, hier ist µ = hλ.
Muss also h so wahlen, dass |F (hλ)| < 1 gilt fur alle Komponenten im System.Bspl.:
c© W. Kley; Skript: Computational Astrophysics
2.6. STABILITAT 27
Im(µ)
µRe( )
Abbildung 2.6: Bereiche der absoluten Stabilitat fur verschiedene ODE-Integratoren. Do-pri5: Nach Dormand & Price (1980): Eingebettete RK5 mit Schrittweitenkontrolle.
1) Explizite RK-VerfahrenDie Bereiche der Stabilitat fur verschiedene Verfahren ist in Abb. 2.6 dargestellt.
D.h bei großen negativen λ braucht man eine kleine Schrittweiten.
Bei hoherer Fehler-Ordnung p −→ großere stabile Bereiche
2) Implizite Trapez-Methode
yk+1 = yk +h
2[f(xk, yk) + f(xk+1, yk+1)]
Fur das Testbeispiel gilt
yk+1 = yk +h
2[λyk + λyk+1]
Fur den Verstarkungsfaktor folgt mit µ = hλ
|F (µ)| =∣∣∣∣2 + µ
2− µ
∣∣∣∣ < 1 ∀µ mit Re(µ) < 0.
also ist Trapezmethode absolut stabil, wie das implizite Euler-Verfahren in der Ab-bildung 2.6.
c© W. Kley; Skript: Computational Astrophysics