2.3 Algorithmen und Strategien in DAESOL
2.3.3 Schrittweiten- und Ordnungssteuerung
Zur Berechnung einer neuen Schrittweite und Ordnung f¨ur den n¨achsten Schritt ¨ ubertra-gen wir die Fehlerformeln aus Schritt n+ 1 auf Schritt n+ 2 und leiten eine Sch¨atzung f¨ur den Fehler im n¨achsten Schritt her, der nur auf Werten beruht, die im Schritt n+ 1 bekannt sind.
Mit Formel (2.19) als Sch¨atzung f¨ur den lokalen Diskretisierungsfehler im Schritt n + 1 erhalten wir f¨ur den akkumulierten Fehler im Schritt n+ 2
τn+2 = ψ1(n+ 2)·. . .·ψk(n+ 2)· ∇k+1[x(tn+2), x(tn+1), . . . , x(tn+1−k)]
+ψ1(n+ 2)·. . .·ψk+1(n+ 2)· ∇k+2[x(tn+2), x(tn+2), x(tn+1), . . . , x(tn+1−k)]
= ψ1(n+ 2)·. . .·ψk(n+ 2)·
∇k+1[x(tn+2), x(tn+1), . . . , x(tn+1−k)]
+ψk+1(n+ 2)· ∇k+2[x(tn+2), x(tn+2), x(tn+1), . . . , x(tn+1−k)]
= ψ1(n+ 2)·. . .·ψk(n+ 2)·
∇k+1[x(tn+1), x(tn), . . . , x(tn−k)] +ψk+2(n+ 2)· ∇k+2[x(tn+2), . . . , x(tn−k)]
+ψk+1(n+ 2)· ∇k+2[x(tn+2), . . . , x(tn−k)]
+ψk+1(n+ 2)·ψk+2(n+ 2)· ∇k+3[x(tn+2), x(tn+2), x(tn+1), . . . , x(tn−k)]
= ψ1(n+ 2)·. . .·ψk(n+ 2)·
∇k+1[x(tn+1), x(tn), . . . , x(tn−k)]
+(ψk+1(n+ 2) +ψk+2(n+ 2))· ∇k+2[x(tn+1), . . . , x(tn−k−1)]
+(ψk+1(n+ 2) +ψk+2(n+ 2))·ψk+3(n+ 2)· ∇k+3[x(tn+2), . . . , x(tn−k−1)]
+ψk+1(n+ 2)·ψk+2(n+ 2)· ∇k+3[x(tn+2), x(tn+2), x(tn+1), . . . , x(tn−k)]
=. ψ1(n+ 2)·. . .·ψk(n+ 2)·
∇k+1[x(tn+1), x(tn), . . . , x(tn−k)]
+(ψk+1(n+ 2) +ψk+2(n+ 2))· ∇k+2[x(tn+1), . . . , x(tn−k−1)]
. F¨ur den lokalen Fehler w¨ahlen wir eine Sch¨atzung analog zu (2.21):
Ek(n+ 2) = hn+2· kτn+2k
˙
= hn+2·ψ1(n+ 2)·. . .·ψk(n+ 2)· k∇k+1xn+1k (2.22)
=: h2n+2·q(hn+2)· k∇k+1xn+1k mit
q(hn+2) = (hn+2+ψ1(n+ 1))·. . .·(hn+2+ψk−1(n+ 1)).
Die Schrittweite sollte so gew¨ahlt werden, daß im n¨achsten Schritt Ek(n+ 2)≤T OL
gilt. Da aus Fehlerformel (2.22) nicht einfach eine neue Schrittweite zu berechnen ist, leiten wir zun¨achst eine vereinfachte Fehlerformel her:
Die analoge Formel zu (2.21) auf ¨aquidistantem Gitter ist Eˆk(n+ 2) :=k!hk+1k∇k+1xn+1k.
Mit ˆEk(n+ 2)≤T OL[ erh¨alt man eine erste Sch¨atzung f¨ur die neue Schrittweite hˆ = k+1
s T OL[
k!k∇k+1xn+1k, (2.23)
die sogenannte maximale gleichm¨aßige Schrittweite, siehe Bleser [Ble86].
Um eine neue Schrittweite und Ordnung f¨ur den n¨achsten Schritt zu bestimmen, wird diese Formel f¨ur verschiedene Ordnungen k0 =k −1, k, k+ 1 mit T OL < T OL, z. B.[ T OL[ =
1
2·T OLausgewertet. Wir verringern oder erh¨ohen die Ordnung um 1, falls die berechnete Schrittweite f¨ur die zugeh¨orige Ordnung sich signifikant gegen¨uber der Ordnung im letzten Schritt ge¨andert hat. Ansonsten wird die aktuelle Ordnung beibehalten.
Fehlerformel (2.23) scheint allerdings ungeeignet f¨ur die Berechnung im Fall variabler Schrittweiten. Zum einen basiert die Fehlerformel auf variablem Gitter (Faktor∇k+1xn+1), zum anderen wird die Schrittweite als konstant angenommen ( (k+ 1)-te Wurzel).
Wir betrachten ˆhals eine erste N¨aherung f¨ur die neue Schrittweite und testen, ob sie auch die genauere Fehlerformel (2.22) auf nicht¨aquidistantem Gitter erf¨ullt:
Ek(n+ 2) =h2q(h)k∇k+1xn+1k ≤T OL. (2.24) Falls obige Ungleichung erf¨ullt ist, so wird die Schrittweite akzeptiert, ansonsten wird sie mit Hilfe von Fehlerformel (2.24) reduziert:
h2 = 1
q(h∗) · T OL
k∇k+1xn+1k. (2.25)
Dabei bezeichne h∗ die zuvor gew¨ahlte Schrittweite.
Bemerkung 2.3.2 (Schrittweiten- und Ordnungssteuerung)
Shampine und Bogacki [SB89] haben gezeigt, daß eine Schrittweitensteuerung basierend auf Approximationen auf ¨aquidistantem Gitter bei der Berechnung einer Schrittweite f¨ur den n¨achsten Schritt zu pessimistisch ist, wenn eine Vergr¨oßerung der Schrittweite ange-strebt wird, andererseits bei Schrittzur¨uckweisungen die Verkleinerung oft zu optimistisch gew¨ahlt wird. Die in DAESOL verwendeten Schrittweitensch¨atzungen gehen zun¨achst von Formeln auf ¨aquidistantem Gitter aus, der sogenanntenmaximalen gleichm¨aßigen Schritt-weite, siehe Bleser [Ble86]. Sie werden aber im Anschluß mit Hilfe der Formeln auf va-riablem Gitter korrigiert. Gleichzeitig wird zugelassen, daß sich die Schrittweite in jedem Schritt ¨andern kann.
Auch die Ordnung wird relativ schnell hochgeschaltet. Dabei wird das von Gear und Watanabe [GW74] aus Stabilit¨atsgr¨unden geforderte Festhalten der Ordnung f¨ur ein bis drei Schritte (abh¨angig von der aktuellen Ordnung) bei einer Erh¨ohung der Ordnung um eins in der Praxis immer eingehalten.
Bleser [Ble86] hat f¨ur die Beispiele aus STIFF DETEST, einer Beispielsammlung von steifen gew¨ohnlichen Differentialgleichungssystemen von Enright et al. [EHL75], gezeigt, daß eine Schrittweiten- und Ordnungssteuerung, die auf den Formeln auf nicht¨aquidi-stantem Gitter beruht, zu einem schnelleren Hochschalten der Ordnung und zu weniger zur¨uckgewiesenen Schritten f¨uhrt als bei konventionellen Sch¨atzungen.
Bemerkung 2.3.3 (Schrittweitensteuerung und Monitor-Strategie)
In vielen Codes wird die Schrittweite ¨uber mehrere Schritte hinweg konstant gehalten, teil-weise auch, um eine neue Auswertung und Zerlegung der IterationsmatrixJ zu vermeiden.
Falls die Schrittweite dann ver¨andert wird, so meist relativ stark, was allerdings leicht zu Stabilit¨atsproblemen f¨uhren kann. Zum anderen kann es leicht vorkommen, daß die Itera-tionsmatrix J˜im Newton-Verfahren, die auch explizit von der Schrittweite abh¨angt, nach großen Schrittweiten¨anderungen keine gute N¨aherung mehr f¨urJ ist. Die Konvergenzraten im Newton-Verfahren werden schlecht und die Iterationsmatrix muß mit aktueller Schritt-weite neu zerlegt werden. Da in den meisten Codes der 2. Schritt aus der Monitor-Strategie fehlt, erfordert dies gleichzeitig eine Auswertung der Ableitungen der Modellfunktionen.
Die oben beschriebene Schrittweitensteuerung f¨uhrt auf eine gleichm¨aßige ¨Anderung der Schrittweiten. Wenn sich die Schrittweite nach einigen Schritten zu sehr ge¨andert hat, so wird zun¨achst nur die Iterationsmatrix J mit aktueller Schrittweite und BDF-Koeffizienten neu zerlegt.
Schrittweitensteuerung nach Schrittzur¨uckweisungen Schrittzur¨uckweisungen k¨onnen aus zwei Gr¨unden auftreten:
• zum einen, weil der gesch¨atzte Fehler zu groß ist
• und zum anderen, weil das Newton-Verfahren nicht innerhalb von drei Schritten konvergiert hat.
Hat das Newton-Verfahren konvergiert, aber der gesch¨atzte Fehler ist zu groß, so berech-nen wir eine neue Schrittweite aus Fehlerformel (2.21)
h(neu)n+1 = h(alt)n+1 ·T OL0
Ek(n+ 1) , (2.26)
mit T OL0 = cred·T OL, cred ≤ 1. Diese Formel beruht wiederum auf Approximationen basierend auf dem variablen Gitter.
Falls das Newton-Verfahren nicht konvergiert hat, so ist k∆x0k eine Approximation f¨ur kxCn+1−xPn+1k=ψ1(n+ 1)·. . .·ψk(n+ 1)· k∇k+1xn+1k
und wir erhalten eine Sch¨atzung f¨ur den lokalen Fehler mit kEkk= h
ψk+1(n+ 1) · k∆x(0)k. (2.27) Die Schrittweite sollte so weit verkleinert werden, daß das Newton-Verfahren innerhalb von zwei Schritten konvergiert (zur Sicherheit wird, falls n¨otig, noch eine weitere Iteration erlaubt).
Das heißt f¨ur die Konvergenzrate
δ0(neu)≤! 1 4.
Da die oben dargestellten Schritte 1. - 4. der Monitor-Strategie nacheinander ausgef¨uhrt wurden, wurde die Zerlegung von ˜J mit aktuellen Matrizen ∂f∂x, ∂A∂x und ∂g∂x, aktuellen BDF-Koeffizienten und aktueller Schrittweite durchgef¨uhrt. Die Gr¨oßeκ in Fromel (2.15) ist dann gleich Null und wir erhalten
δ0 = ω0
2 k∆x0k. Mit
δ0(neu)= ω0
2 k∆x0(neu)k≤! 1 4 haben wir somit
k∆x0(neu)k ≤ k∆x0k 4·δ0
.
Die Fehlerformel aus (2.21) f¨ur die neue Schrittweite sollte dann Ek(n+ 1) ≤ T OL[ an Stelle von Ek(n+ 1) ≤T OL erf¨ullen mit versch¨arfter Toleranz
T OL[ = h
ψk+1(n+ 1) · k∆x0k 4·δ0
.
Wir erhalten eine neue Schrittweite h(neu)n+1 aus Formel (2.26) mit T OL[ an Stelle vonT OL0.