1. Steife Probleme.

Steife Probleme tauchen bei Systemen von Differentialgleichungen auf, deren Jacobi-Matrix durchweg Eigenwerte mit negativem Realteilen aufweisen, die darüberhinaus noch von stark unterschiedlicher Größenordnung sind. Sei µ1 der betragsgrößte und µ2 der betragskleinste Realteil der Eigenwerte einer Jacobi-Matrix, dann wird der Steifigkeitsquotient q dieser Matrix definiert durch:

Durchweg negative Realteile bei den Eigenwerten und ein hoher Steifigkeitskoeffizient sind charakteristisch für steife Probleme und wir setzen dies im folgenden immer voraus. Diese Eigenschaften sind aber noch nicht hinreichend dafür, daß implizite Verfahren zur Lösung besser geeignet sind als explizite.

Bei einem linearen System ergeben sich die Lösungen als gewichtete Summen von Termen der Form exp(µ_it), wobei die µ_i die Eigenwerte der Jacobi-Matrix durchlaufen. Die Komponenten zu den Eigenwerten mit relativ betragsgroßen Realteilen nennen wir steife oder schnelle Komponenten.

Bei nichtlinearen Systemen bezieht sich der Begriff steife Komponente auf die Jacobi-Matrix der jeweils lokalen Linearisierung.

Da der Wert einer steifen Komponente viel schneller gegen 0 geht als der einer vergleichsweise nicht steifen Komponente, ergibt sich häufig der Fall, daß eine steife Komponente nur in Teilen des Integrationsintervalles einen wesentlichen Beitrag zur Gesamtsumme leistet. Ist dies der Fall, so sprechen wir davon, daß die steife Komponente aktiv ist, anderenfalls nennen wir sie passiv, latent oder schlafend. Bei linearen Problemen hat man typischerweise den Fall, daß die steifen Komponenten nur zu Beginn des Integrationsintervalles aktiv sind und dann latent werden, wir nennen diesen teil des Integrationsintervalles die Transientenphase.

Als Beispiel betrachten wir das lineare System 1. Ordnung:

y1'(t) y2'(t)

=

y2(t) - 1 -200 y1(t) - 201 y2(t) + 1

oder anders formuliert:

y'(t) =	0	-1	*	y	+	-1
	-200	-201				-1

mit den Anfangswerten:

y1(0) = 0.99 y2(0) = 1.00

Die Eigenwerte der Jacobi-Matrix sind -1 und -200. Der Steifigkeitsquotient ist also 200. Die Lösung des Problems, auf das wir noch genauer eingehen werden, ist:

y1(t) = 2*exp(-t) - 0.01*exp(-200*t) - 1 y2(t)=-2*exp(-t) + 2.00*exp(-200*t) + 1

Bild 1 zeigt die Lösungskurven zu diesem Problem, diese sind einmal mit linear, das andere mal mit logarithmisch skalierter x-Achse wiedergegeben. Der Eigenwert -200 liefert die steife Komponente. Für die Komponenten von y₁ haben wir:

Für die Komponenten von y₂ haben wir:

Hier ist der Beitrag der steifen Komponente für kleine t sicher nicht zu vernachlässigen, ab t=1.0e-2 fällt er aber auch hier nicht mehr ins Gewicht. Wir haben bei diesem Anfangswertproblem also eine Transientenphase, die von 0 bis höchstens 1.0e-1 reicht.

Die Schrittweiten eines Integrators richten sich im allgemeinen nach dem betragskleinsten Eigenwert der Jacobi-Matrix, alo nach den steifen Komponenten des Systems. Solange diese aktiv sind, in der Transientenphase, muß jeder Integrator entsprechend kleine Integrationsschritte durchführen. Ist µ der zugehörige Eigenwert, so liefert die Zeitkonstante 1/|µ| einen Anhaltspunkt für die Größenordnung des Integrationsschrittes. Wenn allerdings die steifen Komponenten passiv sind, wünscht man sich, daß die Schrittweiten nur von den aktiven Komponenten bestimmt wird. Integratoren, die ihre Schrittweiten so steuern können, nennen wir steife Integratoren. Wir werden sehen, daß explizite Integrationsverfahren nicht als steife Integratoren geeignet sind.

2. Matlab Integratoren für steife Probleme.

Name	Beschreibung
ode15s	der Standard-Integrator für steife Probleme im Matlab Angebot, ode15s ist insbesondere bei höheren Genauigkeitsanforderungen bei steifen Problemen die erste Wahl under den steifen Matlab Integratoren. ode15s verwendet Formeln variabler Ordnung aus der Familie der 'numerical differentiation formulas (NDF)'. Optional können stattdessen die BDF Formeln (Gear'sche Formeln) benutzt werden, die im allgemeinen aber weniger effizient sind. ode15s führt mit jeder Option ein Mehrschritt-Verfahren durch.
ode23s	Einschrittverfahren niedriger Ordnung für steife Systeme, basierend auf einer modifizierten Rosenbrock-Formel. Bei geringen Genauigkeitsanforderungen häufig effizienter als ode15s, für hohe Genauigkeitsanforderungen jedoch weniger geeignet.
ode23t	Integrator für moderat steife Probleme, basierend auf der Trapezregel. Verwendbar bei Problemen, bei denen numerische Dämpfung unerwünscht ist. Für hohe Genauigkeitsanforderungen ist ode23t wenig geeignet.
ode23tb	Implementation des TR-BDF2 Verfahrens, das auch als implizites Runge-Kutta Verfahren gedeutet werden kann. Bei geringen Genauigkeitsanforderungen häufig effizienter als ode15s, für hohe Genauigkeitsanforderungen weniger geeignet.

3. Beschreibung der Jacobi-Matrix

Die Integration unseres linearen Testbeispieles mit einem steifen Integrator unterscheidet sich im Aufruf daher nur durch den anderen Integratornamen von der Integration mit einem expliziten Namen.

Beispiel 3.1:

[t,y] = ode15s('LinODE',[0 10],[0.99 1]); plot(t,y);

mit dem ODE-File LinODE:

function dy=LinODE(t,y) % Lineares, steifes Problem dy = [y(2)-1; -200*y(1)-201*y(2)+1;]

Manchmal ist es jedoch vorteilhaft, die Jacobi-Matrix explizit im ODE-File zu beschreiben. Diese wird dann mit einem dritten Parameter aufgerufen. Dieser Parameter fehlt wenn nur die Ableitungen dy/dt in zu einem Zeitpunkt t abgerufen werden und hat den Wert 'jacobian', wenn die Jacobi-Matrix abgerufen wird. Der Integrator wird mit Hilfe der Funktion odeset darüber informiert, daß eine Beschreibung der Jacobi-Matrix verfügbar ist.

Beispiel 3.2:

opt=odeset('Jacobian','on'); [t,y] = ode15s('LinODEj',[0 10],[0.99 1],opt); plot(t,y);

mit dem ODE-File LinODEj:

function dy=LinODEj(t,y,flag) % Lineares, steifes Problem if flag eq '' dy = [y(2)-1; -200*y(1)-201*y(2)+1;] else if flag eq 'jacobian' dy = [0,1;-200,-201]; end

Es gibt weitere Optionen, mit denen man für den Integrator hilfreiche Eigenschaften der Jacobi-Matrix angeben kann. In Beispiel 3.2 ist die Jacobi-Matrix konstant, was durch die Eigenschaft 'JConstant' mit dem Wert 'on' angegeben werden kann. Falls die Jacobi-Matrix dünnbesetzt ist, kann man das Besetzungsmuster mit Hilfe der Eigenschaft 'JPattern' beschreiben.

Die Eigenschaft 'Vectorized' ist ebenfalls im Zusammenhang mit impliziten Integratoren von Bedeutung, bezieht sich aber auf die Beschreibung der Ableitungen dy/dt und nicht die der Jacobi-Matrix.
Die Integratoren berechnen zu jeder Lösungskomponente y_i(t) eine Folge von Werten y(i,j), wobei j dem Zeitpunkt t_j entspricht. Wenn das ODE-File so geschrieben ist, daß es in einem Aufruf die gesamte Folge y(i,1:n) zurückliefern kann, so nennt man das ODE-File vektorisiert. Beispiel 3.3:

opt=odeset('Vectorized','on'); [t,y] = ode15s('LinODEv',[0 10],[0.99 1],opt); plot(t,y);

mit dem ODE-File LinODEv:

function dy=LinODEv(t,y) % Lineares, steifes Problem, vektorisiert dy = [y(2,:)-1; -200*y(1,:)-201*y(2,:)+1;];

Die nachfolgende Tabelle zeigt, daß der implizite Integrator ode15s beim vorliegenden Problem deutlich schneller ist als der explizite Integrator ode15s, sofern deas Problem hinreichend lange steif ist. In der Transientenphase ist dagegen ode45 der schnellere Integrator. Die anderen Optionen wie 'Jacobian', 'JConstant' und 'Vectorized' haben bei derart einfachen Problemen keine verkürzende Wirkung auf die Durchlaufzeit.

Intervall	1.0e-2	1.0e-1	1.0	1.0e1	1,0e2
ode45	0.109	0.172	0.438	3.094	30.391
ode15s	0.28	0.52	0.61	0.83	0.95

• Steifigkeit bei linearen Problemen
• Ruby Laser Oszillator
• Van der Pool'sche Differentialgleichung