Stabilitätsfunktion

Die Stabilitätsfunktion ist in der Numerik ein Hilfsmittel, um Lösungsverfahren für gewöhnliche Differentialgleichungen zu analysieren. Die einfache Testgleichung von Germund Dahlquist <math>y'(t)=\lambda y(t),\ y(0)=1</math> mit <math>\lambda\in\Complex</math> besitzt als Lösung die Exponentialfunktion <math>y(t)=e^{\lambda t}</math>. Bei den meisten Verfahren für gewöhnliche Differentialgleichungen kann man die berechnete Näherungslösung nach einem Zeitschritt mit einer Schrittweite <math>h</math> ebenfalls als eine Funktion schreiben, die nur vom Produkt <math>z=h\lambda\in\Complex</math> abhängt. Diese Funktion ist die Stabilitätsfunktion und wird oft mit <math>R(z)</math> bezeichnet. Durch einen Vergleich mit der Exponentialfunktion <math>e^z=e^{h\lambda}</math> bekommt man grundlegende Informationen über das numerische Verfahren. So beziehen sich einige Stabilitätsbegriffe auf die Eigenschaften von <math>R(z)</math>.

Stabilitätsgebiet und Stabilitätsbegriffe

Mit Hilfe der Stabilitätsfunktion <math>R(z)</math> lässt sich das Stabilitätsgebiet <math>S</math> beschreiben und berechnen in der Form

<math>S=\{z\in\Complex: |R(z)|<1\}.</math>

Denn bei Einschrittverfahren gilt für die Näherungen <math>y_n</math> zum Zeitpunkt <math>t_n=nh</math> die Beziehung <math>y_n=R(z)y_{n-1}=\ldots=\left(R(z)\right)^jy_{n-j}=\ldots=\left(R(z)\right)^ny_0</math> und daher gilt

<math>y_n\xrightarrow{n\to\infty} 0 \iff z\in S.</math>

Wenn <math>S</math> die ganze linke komplexe Halbebene umfasst, heißt das Verfahren A-stabil. Dann ist der Betrag von <math>R</math> in der ganzen offenen linken Halbebene kleiner als 1. Besonders günstig für ein Verfahren ist es, wenn <math>R(z)</math> außerdem noch den Grenzwert 0 hat, wenn <math>z</math> auf der reellen Achse gegen <math>-\infty</math> strebt, sodass sich also der Betrag von <math>R(z)</math> dort asymptotisch wie die Exponentialfunktion verhält. Dann heißt das Verfahren L-stabil.

Beispiel

Das explizite Euler-Verfahren <math>y_{n+1} = y_n + hf(t_n, y_n)</math> ergibt für die Testgleichung mit <math>f(t,y) = \lambda y</math> nach einem Schritt

<math>y_1 = y_0 + h \lambda y_0 = (1 + h \lambda) y_0</math>,

also gilt für seine Stabilitätsfunktion <math>R(z) = 1+z</math>. Sein Stabilitätsgebiet besteht daher aus allen komplexen Zahlen <math>z</math> mit <math>|1+z| < 1</math>, was dem Inneren des Kreises mit Mittelpunkt <math>-1</math> und Radius <math>1</math> in der komplexen Zahlenebene entspricht.

Für das implizite Euler-Verfahren <math>y_{n+1} = y_n + hf(t_{n+1}, y_{n+1})</math> folgt dagegen mit <math>f(t,y) = \lambda y</math>

<math>y_1 = y_0 + h\lambda y_1 \iff y_1 = \frac{1}{1-h \lambda} y_0</math>,

also <math>R(z) = \frac{1}{1-z}</math>. Das Stabilitätsgebiet ist nun durch die Bedingung <math>|\tfrac{1}{1-z}| < 1</math> gegeben, die mit

gleichwertig ist, was dem Äußeren des Kreises mit Mittelpunkt <math>1</math> und Radius <math>1</math> entspricht. Es enthält daher die ganze offene linke Halbebene und somit ist das implizite Euler-Verfahren A-stabil. Wegen <math>\lim_{z \to -\infty} \frac{1}{1-z} = 0</math> ist es sogar L-stabil.

Die Stabilitätsfunktion von Runge-Kutta-Verfahren

Runge-Kutta-Verfahren sind vollständig durch die Koeffizienten <math>A, b, c</math> aus ihrem Butcher-Tableau festgelegt. Bei der Testgleichung ist der Anfangswert <math>y_0=1</math> und für die Stufen ergibt sich im ersten Zeitschritt

<math>k_i=\lambda\left(1+h\sum_{j=1}^sa_{ij}k_j\right),\quad i=1, \dotsc, s.</math>

Dies ist ein quadratisches lineares Gleichungssystem für den Vektor <math>k=(k_1, \dotsc, k_s)^T</math> in der Form <math>(I-zA)k=\lambda e</math> mit dem Vektor <math>e=(1, \dotsc, 1)^T.</math> Mit dessen Lösung bekommt man dann die Runge-Kutta-Näherung <math>y_1\approx y(h)</math> in der Form

Dies ist bei Runge-Kutta-Verfahren eine rationale Funktion, daher wird sie gerne mit <math>R(z)</math> bezeichnet.

Bei expliziten Runge-Kutta-Verfahren ist die Koeffizientenmatrix <math>A</math> eine strikt untere Dreiecksmatrix, daher bricht die Neumann-Reihe von <math>(I-zA)^{-1}</math> nach s<math></math> Summanden ab und man bekommt

<math>R(z)=1+zb^T(I-zA)^{-1}e=1+zb^T e+z^2b^TAe+\dotsb+z^sb^TA^{s-1}e.</math>

Daher ist die Stabilitätsfunktion eines expliziten Runge-Kutta-Verfahrens ein Polynom, solche Verfahren können nicht A-stabil sein.

Bei impliziten Runge-Kutta-Verfahren sind aber z. B. die Gauß-Legendre-Verfahren A-stabil. Die Stabilitätsfunktionen dieser speziellen Verfahren sind sogar sehr gute Approximationen an die Exponentialfunktion, nämlich die sogenannten Padé-Approximationen.

Die Stabilitätsfunktion von Mehrschrittverfahren

Wendet man ein lineares Mehrschrittverfahren <math>\sum_{j=0}^m\alpha_jy_{n-j}=h\sum_{j=0}^m\beta_jf(y_{n-j})</math> auf die Testgleichung an, ergibt sich wieder mit <math>z=h\lambda</math> die Gleichung

<math>\sum_{j=0}^m\alpha_jy_{n-j}-z\sum_{j=0}^m\beta_jy_{n-j}=\sum_{j=0}^m(\alpha_j-z\beta_j)y_{n-j}=0.</math>

Dies ist eine lineare Differenzengleichung, die man einfach lösen kann. Denn die Folge <math>y_n=u^n</math> ist eine nichttriviale Lösung dieser Differenzengleichung, wenn <math>u</math> eine Nullstelle des charakteristischen Polynoms

<math>0\stackrel!=\sum_{j=0}^m\alpha_ju^{m-j}-z\sum_{j=0}^m\beta_ju^{m-j}=\varrho(u)-z\sigma(u)</math>

ist, wobei man die Polynome

<math>\varrho(u)=\sum_{j=0}^m\alpha_ju^{m-j}</math>

<math>\sigma(u)=\sum_{j=0}^m\beta_ju^{m-j}</math>

eingeführt hat. Also bekommt man mit den von <math>z</math> abhängenden Nullstellen <math>u</math> des Polynoms <math>\varrho(u)-z\sigma(u)</math> die Lösungen <math>u^n</math> zur Testgleichung und daher liegt <math>z</math> im Stabilitätsgebiet des Verfahrens, wenn alle diese Lösungen gegen 0 gehen für <math>n\to\infty</math>. Daher kann man die betragsmaximale Nullstelle <math>|u(z)|</math> als Stabilitätsfunktion des Verfahrens ansehen.

Datei:Stability region for BDF6.svg

Stabilitätsgebiet für das 6-stufige BDF-Verfahren

Diese Interpretation erscheint sehr unhandlich. Allerdings interessiert man sich oft weniger für die Stabilitätsfunktion, sondern für das Stabilitätsgebiet <math>S</math>. Der Rand dieses Gebietes besteht aus denjenigen <math>z\in\Complex</math>, bei dem für die Nullstellen <math>|u|=1</math> gilt, wo die Nullstellen also auf dem komplexen Einheitskreis liegen. Da <math>\varrho(u)-z\sigma(u)=0\Leftarrow z=\varrho(u)/\sigma(u)</math> gilt, ist die Bestimmung des Stabilitätsgebiets bei Mehrschrittverfahren sogar besonders einfach, denn seinen Rand erhält man i. W. explizit durch

<math>\partial S=\Big\{\frac{\varrho(u)}{\sigma(u)}:\, |u|=1\Big\}=\Big\{\frac{\varrho(e^{i\varphi})}{\sigma(e^{i\varphi})}:\,\varphi\in[0,2\pi)\Big\}.</math>

Als Beispiel wird das Stabilitätsgebiet für das 6-stufige BDF-Verfahren gezeigt.

Die Stabilitätsfunktion von allgemeinen linearen Verfahren

Obwohl auch Mehrschrittverfahren in der Gestalt von allgemeinen linearen Verfahren geschrieben werden können, ist die Struktur ähnlich derjenigen der Runge-Kutta-Verfahren weiter oben. Daher bekommt man ein ähnliches Ergebnis. Für den Vektor <math>Y</math> der Stufenlösungen gilt

<math>Y=zAY+Uy^{[n-1]}\quad \Rightarrow Y=(I-zA)^{-1}Uy^{[n-1]}</math>

und der Zeitschritt wird daher zu

In jedem Zeitschritt erfolgt also die Multiplikation mit derselben Matrix

Es gilt daher <math>y^{[n]}=M(z)^ny^{[0]}\to0\,(n\to\infty)</math>, wenn die Potenzen von <math>M(z)</math> gegen 0 gehen, also alle Eigenwerte von <math>M(z)</math> innerhalb des komplexen Einheitskreises liegen. Daher kann man hier den Spektralradius von <math>M(z)</math> als Stabilitätsfunktion <math>R(z)</math> in der Definition des Stabilitätsgebiets <math>S</math> ansehen.

Weitergehende Bedeutung für lineare Systeme

Die obige Testgleichung von Dahlquist ist sehr einfach, hat aber eine weitergehende Bedeutung für Systeme von linearen, autonomen und homogenen Differentialgleichungen

<math>y'(t)=Q y(t),\quad y(0)=y_0,\quad Q\in\R^{d\times d}.</math>

Die exakte Lösung ist <math>y(t)=e^{tQ}y_0</math> mit dem Matrixexponential <math>e^{tQ}</math>. Die numerische Lösung <math>y_n</math> kann man jetzt mit der Matrix-Stabilitätsfunktion <math>R(tQ)</math> darstellen. Wenn dabei <math>J=P^{-1}QP</math> die Jordan-Normalform von <math>Q \ (=PJP^{-1})</math> ist, gilt

Bei einer diagonalisierbaren Matrix <math>Q</math> ist, ist <math>R(hJ)</math> eine Diagonalmatrix mit den Diagonalelementen <math>R(h\lambda_j)</math>. Wenn für alle Eigenwerte <math>\lambda_j</math> von <math>Q</math> gilt, dass <math>h\lambda_j\in S</math> ist, dann konvergiert auch hier <math>y_n\to 0\,(n\to\infty)</math>. Bei dieser Differentialgleichung sieht man gleichzeitig, dass es sinnvoll ist, <math>S</math> als offene Menge zu definieren. Denn im diagonalisierbaren Fall bleiben zwar Lösungen auf dem Rand mit <math>h\lambda_j\in\partial S</math> noch beschränkt, aber im Allgemeinen nicht mehr, wenn mehrfache Eigenwerte mit Jordanblöcken auftreten.

Literatur

E. Hairer, G. Wanner: Solving Ordinary Differential Equations II, Stiff problems. Springer Verlag.
K. Strehmel, R. Weiner, H. Podhaisky: Numerik gewöhnlicher Differentialgleichungen – Nichtsteife, steife und differential-algebraische Gleichungen. Springer Spektrum, 2012.