Steifigkeit bei linearen Problemen.

Bild 1 zeigt die Lösungskurven zu diesem Problem, diese sind einmal mit linear, das andere mal mit logarithmisch skalierter x-Achse wiedergegeben.

Beide Eigenwerte sind negativ und liefern einen Steifigkeitsquotienten von 200. Die Eigenwerte und der Steifigkeitsquotient sind konstant, da das System linear ist.

Numerische Integration

Angegeben ist der Aufwand für die beiden expliziten Integratoren ode23 (rot, durchgezogen) und ode45 (blau, durchgezogen), sowie für die impliziten Integratoren ode23s (rot, gestrichelt), ode15s (blau, gestrichelt) und ode23t (schwarz, gestrichelt).

Das unterschiedliche Verhalten der beiden Gruppen von Integratoren ist offensichtlich. Bei den expliziten Integratoren wächst der Aufwand linear mit der Intervalllänge. Diese Integratoren müssen für ein bestimmtes Integrationsintervall von t nach t+delta unabhängig von t etwa denselben Aufwand treiben. Bei den impliziten Integratoren ist der Aufwand zur Integration von t nach t+delta dagegen für kleine t deutlich größer als für große t.

Bild 3 zeigt die Schrittweiten, mit denen die einzelnen Integratoren bei der Integration von 0 bis 30 arbeiten, auch hier sieht man deutlich das unterschiedliche Verhalten der beiden Gruppen. Am Anfang des Integrationsintervalles benutzen die Integratoren beider Gruppen in etwa dieselben Schrittweiten. Während die Schrittweiten der impliziten Integratoren aber im Verlaufe der Integration auf etwa 2.5 anwachsen, bleiben die der expliziten Integratoren im Bereich von 2.0e-2.

Die Schrittweite der expliziten Integratoren liegt im gesamten Integrationsintervall in der Größenordnung von 1/|µ1|, wobei µ1 durch den betragsgrößten Eigenwert µ1=-200 gegeben ist.

Die Schrittweite der impliziten Integratoren liegt nur zu Beginn der Integration in diesem Bereich und richtet sich dann nach 1/|µ2|, wobei µ2=-1 der zweite Eigenwert ist.

Bild 3a zeigt zum Vergleich die Schrittweiten bei einem linearen Differentialgleichungssystem y' = M*y + v mit der Matrix:

M =	0	+1
	-20	-21

Vektor v, Anfangswerte und Genauigkeitsanforderungen sind dieselben wie bisher. Eigenwerte dieser Matrix sind -1 und -20. Die nicht-steifen Integratoren erreichen hier, wie zu erwarten, etwa 10 mal längere Schrittweiten wie beim ursprünglichen System, die steifen Integratoren erreichen dieselben Schrittweiten. Da sich der betragskleinere Eigenwert nicht geändert hat, ist auch dies zu erwarten.

Schaut man sich die Lösungskurven des Problems genauer an, so sieht man, daß die steife Komponente exp(-200*t) nur am Anfang des Integrationsintervalles, etwa bis t=1.0e-2, einen wesentlichen Beitrag zu den Lösungswerten liefert. Bild 4 zeigt diesen Anteil für y₁ und y₂ in graphischer Form. In dieser Transientenphase kann die steife Komponente nicht vernachlässigt werden, jeder Integrator muß seine Schrittweite am Eigenwert µ1=-200 orientieren, also Schrittweiten in der Größenordnung von 1/|µ1| verwenden. Steife (implizite) Interatoren müssen hier dieselben kleinen Schritte machen wie die nicht-steifen (expliziten) Integratoren, da letztere weniger Aufwand pro Schritt erfordern sind sie vorzuziehen.

Nach der Transientenphase ist der Beitrag der steifen Komponente zu den Lösungswerten dagegen minimal. Die steifen Integratoren können dies ausnutzen und ihre Schrittweiten jenen Komponenten anpassen, die noch aktiv sind, also noch einen nennenswerten Beitrag zu den Lösungswerten leisten, die nicht-steifen Integratoren können dies nicht, ihre Schrittweite wird weiterhin durch µ1 bestimmt.

Eine Konsequenz daraus ist, daß die Wahl zwischen einem steifen und einem nicht-steifen Integrator nicht nur vom Steifigkeitsquotienten, sondern auch von der Länge des Integrationsintervalles und von den Anfangswerten abhängt.
Die folgende Tabelle listet die Laufzeiten zur Integration des Testproblems über dem Intervall von 0 bis 1.0e-2 auf, dieses Intervall liegt in der Transientenphase des Problems. Merkliche Unterschiede ergeben sich erst bei höheren Genauigkeitsanforderungen. Die expliziten Integratoren höherer Ordnung sind hier, wie zu erwarten, schneller als die impliziten Integratoren höherer Ordnung. Wirklich auffallend ist aber der Unterschied zwischen den Verfahren niedriger Ordnung (ode23, ode23s, ode23t), die für hohe Genauigkeitsanforderungen ungeeignet sind, und den Verfahren höherer Ordnung.

Tabelle 1: Laufzeiten, Integration bis 1.0e-2

Tabelle 1a: Schrittweiten, Integration bis 1.oe-2

• Die Verfahren niederer Ordnung sind die langsamsten, interessanterweise sind die steifen Integratoren (ode23s und ode23t) niederer Ordnung bei einer relativen Toleranz von 1.oe-11 sogar langsamer als die nicht-steifen Integratoren niederer Ordnung.
• Die expliziten Verfahren höherer Ordnung (ode113 nur mit etwas gutem Willen) liegen in der Mittelgruppe.
• Die impliziten Integratoren höherer Ordnung sind die schnellsten.

Die Frage, ob ein Problem besser mit einem impliziten oder expliziten Verfahren zu behandeln ist, hängt nicht nur vom Steifigkeitskoeffizienten und dem Integrationsintervall (Transientenphase, Anfangswerte) ab, sondern auch von den Genauigkeitsanforderungen. Werden diese so hoch gesetzt, daß sie die Schrittweitenwahl bestimmen, so verlieren die impliziten Integratoren ihre Vorteile gegenüber den expliziten.

Tabelle 2: Laufzeiten, Integration bis 100

Tabelle 2a: Schrittweiten, Integration bis 100

Ob man ein lineares Problem nun mit einem steifen oder nicht-steifen Integrator behandeln soll, hängt also nicht nur vom:

• Steifigkeitsquotienten
• Integrationsintervall (Anteil der Transientenphase)

• Genauigkeitsanforderungen

Außerhalb der Transientenphase haben die beiden Lösungskomponenten unseres Testproblems praktisch die Werte 1 und 0. Die Integratoren machen also im ganzen Intervall danach nichts anderes als aus vorliegenden 0 bzw. 1 Werten neue 0 bzw. 1 Werte zu berechnen. Man kann die Probleme, die die nicht-steifen Integratoren gegenüber den steifen Integratoren bei der Integration eines steifen Problems haben, verdeutlichen, wenn man die numerisch berechneten Lösungen in t-Richtung stark gestaucht darstellt. Bild 5 und Bild 5a zeigen dies.

Sie sehen dort jeweils den Verlauf der numerisch berechneten Lösungskurve y₂(t) über dem Intervall von t=29.5 bis 30, wobei die y-Achse immer einen Ausschnitt von 1-2*RelTol bis 1+2*RelTol zeigt. RelTol hat bei dieser Rechnung den Wert 1.0e-5. Während die steifen Integratoren trotz der Stauchung noch glatte Kurvenzüge liefern, liefern die nicht-steifen Integratoren Zickzack-Kurven. Die Schrittrichtung der steifen Integratoren sind also wesentlich besser als die der nicht-steifen Integratoren.

WARUM IST DAS SO ?

AUTOMATISCHES ERKENNEN WANN TRANSIENTENPHASE BEENDET

..................................................................

Zu diesen Anfangswerten [2 0] stellt:

• Bild 1 die exakte Lösung dar, sowie den Anteil, der von der langsamen Komponente alleine geliefert wird. Es hängt natürlich von der geforderten Lösungsgenauigkeit ab, ab wann man sagen kann, daß der Anteil der schnellen bzw. steifen Komponente vernachlässigt werden kann. Wenn man vom optischen Eindruck dieses Bildes ausgeht, wird man etwa 0.05 sagen.
• Bild 2 den Aufwand (flops) dar, den verschiedene Integratoren bei unterschiedlich langen Integrationsintervallen treiben. Man sieht, daß ab 1.6 alle steifen Integratoren weniger Gleitkommaoperationen als alle nicht-steifen Integratoren erfordern und daß der Aufwand der nicht-steifen Integratoren für noch längere Integrationsintervalle erheblich schneller wächst als der für die steifen.
• Bild 3 den Aufwand (flops) dar, den ode23 und ode23s mit wachsendem Integrationsintervall treiben müssen. Außerdem ist der Aufwand dargestellt, der entsteht, wenn man mit ode23 bis 0.05 bzw. 0.1 integriert und dann mit ode23s weiterintegriert. Obwohl der Start eines Integrators immer mit erheblichem Mehraufwand verbunden ist, sieht man, daß auf diese Weise für die meisten Integrationsintervalle bei den hier geforderten Rechengenauigkeiten (AbsTol=1e-4, RelTol=1e-3) der Aufwand für den hybriden Ansatz konkurrenzfähig ist.

• Bild 4 ist da Äquivalent zu Bild 2, jedoch mit dem Wert 1.0-6 für die Eigenschaft 'AbsTol' anstelle des Default-Wertes 1.0e-3. Die Grenze, ab der die steifen Integratoren weniger Aufwand erfordern als die nicht-steifen, ist hier gegenüber Bild 2 deutlich nach hinten verschoben.
• Bild 5 ist da Äquivalent zu Bild 3, jedoch mit den Werten AbsTol=1e-6 und RelTol=1e-6. Bei diesen Genauigkeite ist es offensichtlich von Vorteil, den Anfang des Intervalles mit ode23 und den Rest, ab etwa 0.05 mit ode23s zu integrieren.

• das Differentialgleichungssystem liefert Lösungen mit sehr unterschiedlich schnellen Lösungskomponenten und bringt dadurch die Voraussetzungen für steife Probleme mit. Manche Autoren nennen solche Systeme steif, wobei ein steifes System nicht unbedingt zu einem steifen Integrationsproblem führt. Im Verlaufe der Zeit weisen die Lösungen eines solchen Systems in der Regel Phasen auf, in denen die schnellste oder die schnellsten Komponenten nur einen unwesentlichen Beitrag zur Gesamtlösung beitragen. Die Integration über einen solchen Bereich ist ein steifes Problem. In der Regel gibt es aber auch Phasen, in denen der Beitrag der schnellsten Komponente zur Gesamtlösung nicht unwesentlich ist, die Integration über einen solchen Bereich ist nicht-steif, denn alle Integratoren müssen hier ihre Schrittweite nach der schnellten Komponente richten, also hinreichend kleine Schrittweiten wählen. Wir wollen diese Bereiche steife bzw. nicht-steife Bereiche nennen. Manchmal bezeichnet man steife Bereiche auch als transiente Phasen eines Problems.
  Bei einem linearen Anfangswertproblem ist die Problematik sehr viel einfacher als im nicht-linearen Fall. Die Lösungskomponenten lassen sich durch die Eigenwerte des Systems bestimmen. Wenn alle Eigenwerte negativ und darüberhinaus von unterschiedlicher Größenordnung sind, in unserem Beispiel sind sie -1 und -100, sind die Voraussetzungen für ein steifes Problem gegeben. Darüber hinaus kann die schnellste Komponente nur am Anfang des Integrationsintervalles einen wesentlichen Beitrag zur Gesamtlösung leisten. Ob dies überhaupt der Fall ist, hängt im linearen Fall von den Anfangswerten ab.
  In unserem Beispiel mit den Anfangswerten [2 0] haben wir am Anfang einen nicht-steifen Bereich, der dann allmählich in einen steifen Bereich übergeht und 'immer steifer' wird.
• das Integrationsintervall bestimmt, ob die Integration in einem steifen oder nicht-steifen Bereich oder aber über eine Folge von Bereichen durchzuführen ist, in denen sich steife und nicht-steife Bereiche abwechseln. Wobei zu anzumerken ist, daß der Übergang zwischen solchen Bereichen verschieden schnell ausfallen kann.
• bei unterschiedlichen Genauigkeitsanforderungen verschieben sich die Grenzen zwischen steifen und nicht-steifen Bereichen, weil unterschiedliche Integratoren (im Zweifelsfalle ode15s) unterschiedlich auf Genauigkeitsanforderungen reagieren.

Beispiel 2.3

   y1(t) = exp(-t)
   y2(t) = -exp(-t)

• bei geringen Genauigkeitsanforderungen sind die steifen Integratoren effizienter als die expliziten Integratoren, bei höheren Anforderungen verschwindet dieser Unterschied und der Unterschied zwischen Verfahren höherer und niedrigerer Ordnung ist wesentlich bedeutsamer. Dies ist nicht erstaunlich, da wir schon wissen, daß die Steifigkeit des Problems erst durch Rechenungenauigkeiten eingeschleppt wird.
• es gibt keinen einfachen Zusammenhang zwischen den benötigten Gleitkommaoperationen und der Durchlaufszeit, der über mehrere Verfahren hinweg Bestand hat.

2.3 Absolut stabile Integration

   y' = -k y,   y(0) = 1

Was passiert nun, wenn wir die Schrittweite erhöhen? Wenn es darum geht, allein die Testgleichung zu integrieren, werden wir das nicht tun. Wenn die Testgleichung aber nur eine Lösungskomponente beschreibt, die gegenüber einer anderen Lösungskomponente kaum in Erscheinung tritt, wie dies bei steifen Problemen der Fall ist, ist dies ja gerade das Verhalten, das wir von einem steifen Integrator erwarten.

Das explizite Eulerverfahren, liefert für unser Testproblem:

   y(istep+1) = ( 1 - h*k ) y(istep);

Ein numerisches Verfahren integriert unser Testproblem absolut stabil, wenn die Folge der gelieferten Werte | y_n | monoton gegen 0 fällt.

Solange h*k < 2 bleibt, wird unser Testproblem von allen betrachteten Integratoren absolut stabil integriert. Neben dem expliziten Euler-Verfahren also auch vom:

• impliziten Eulerverfahren

     yn+1 = yn/(1+kh)

• der (impliziten) Trapezregel

     yn+1 = yn*(2-kh)/(2+kh)

In den Bildern zu h*k=5 und h*k=20 werden nur noch das implizite Eulerverfahren und das Trapezverfahren berücksichtigt. Beide Verfahren integrieren das Testproblem absolut stabil, unabhängig davon, wie groß der Wert von h*k gewählt wird.

Zusammenfassend können wir zum Unterschied zwischen expliziten und steifen Integratoren sagen:

• das explizite Eulerverfahren integriert das Testproblem nur für relativ kleine Werte h*k absolut stabil, insbesondere bleibt dieser Bereich wegen der Bedingung h*k < 2 an k gebunden. Für ein steifes Problem bedeutet dies, daß die Schrittweite des expliziten Eulerverfahrens (und dies gilt allgemein für alle expliziten Verfahren) sich immer nach der schnellsten Lösungskomponente richten muß, wenn das System als ganzes absolut stabil integriert werden soll. Das bedeutet aber entsprechend kleine Schrittweiten.
• das implizite Eulerverfahren und die (implizite) Trapezregel integrieren das Testproblem für alle Werte h*k > 0 absolut stabil. Bei einem steifen Problem kann die Schrittweite unabhängig von k, also unabhängig von der schnellsten Lösungskomponente gewählt werden, ohne daß die stabile Integration gefährdet wird.

Wir haben bisher vorausgesetzt, daß der Wert k im Testproblem reell ist, es ist jedoch kein Problem, obige Überlegungen auf komplexe Werte k auszudehnen. Im Bild ist der Verlauf der Oberfläche von exp(-k*t) für k=-1+i dargestellt. Das Verhalten verschiedener Integratoren für dieses Testproblem bei 'größeren' Schrittweiten h wird in einem anderen Bild dargestellt. Sie sehen hier jeweils links den Realteil der gelieferten Lösung, rechts den Imaginärteil.

2.4 A-Stabilität

   y' = k y, y(0)=y0, k komplex

Dahlquist (1963) nennt einen Integrator A-Stabil, wenn sein absoluter Stabilitätsbereich die gesamte negative Halbebene der komplexen Zahlen umfaßt.

Wir haben gesehen, daß das explizite Eulerverfahren nicht A-stabil sein kann, der Bereich der absoluten Stabilität ist bei diesem Verfahren gegeben durch das Innere des Kreises um -1 mit Radius 1.

Das implizite Eulerverfahren und die (implizite) Trapezregel sind dagegen A-stabil, dasselbe gilt auch für das BFD-2 Verfahren:

yn+2 - 4/3 yn+1 + 1/3 yn = 2/3 h yn+2

A-stabile Verfahren sind offensichtlich gut geeignet zur Lösung steifer Probleme. Die Klasse der A-stabilen Verfahren ist allerdings klein. Dahlquist (1963) hat gezeigt:

• ein explizites lineares Mehrschrittverfahren kann nicht A-stabil sein,
• ein implizites, A-stabiles lineares Mehrschritt-Verfahren hat maximal die Ordnung 2,
• unter allen impliziten, A-stabilen linearen Mehrschritt-Verfahren hat die Trapezregel die kleinste Fehlerkonstante.

• ein A(0)-stabiler Integrator hat also einen Bereich der absoluten Stabilität, der die negative, reelle Halbachse der komplexen Zahlenebene umfaßt,
• ein A(45)-stabiler Integrator integriert alle Testgleichungen mit -imag(k) <= real(k) <= imag(k), real(k) < 0 absolut stabil,
• ein A(90)-stabiler Integrator ist A-stabil.

2.5 Dämpfung und L-Stabilität

2.6 Fallbeispiele

1. Ruby Laser Oszillator