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 $y'(t)=\lambda y(t),\ y(0)=1$ mit $\lambda \in \mathbb {C}$ besitzt als Lösung die Exponentialfunktion $y(t)=e^{\lambda t}$. Bei den meisten Verfahren für gewöhnliche Differentialgleichungen kann man die berechnete Näherungslösung nach einem Zeitschritt mit einer Schrittweite $h$ ebenfalls als eine Funktion schreiben, die nur vom Produkt $z=h\lambda \in \mathbb {C}$ abhängt. Diese Funktion ist die Stabilitätsfunktion und wird oft mit $R(z)$ bezeichnet. Durch einen Vergleich mit der Exponentialfunktion $e^{z}=e^{h\lambda }$ bekommt man grundlegende Informationen über das numerische Verfahren. So beziehen sich einige Stabilitätsbegriffe auf die Eigenschaften von $R(z)$.

Stabilitätsgebiet und Stabilitätsbegriffe

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

$S=\{z\in \mathbb {C} :|R(z)|<1\}.$

Denn bei Einschrittverfahren gilt für die Näherungen $y_{n}$ zum Zeitpunkt $t_{n}=nh$ die Beziehung $y_{n}=R(z)y_{n-1}=\ldots =\left(R(z)\right)^{j}y_{n-j}=\ldots =\left(R(z)\right)^{n}y_{0}$ und daher gilt

$y_{n}{\xrightarrow {n\to \infty }}0\iff z\in S.$

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

Beispiel

Das explizite Euler-Verfahren $y_{n+1}=y_{n}+hf(t_{n},y_{n})$ ergibt für die Testgleichung mit $f(t,y)=\lambda y$ nach einem Schritt

$y_{1}=y_{0}+h\lambda y_{0}=(1+h\lambda )y_{0}$,

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

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

$y_{1}=y_{0}+h\lambda y_{1}\iff y_{1}={\frac {1}{1-h\lambda }}y_{0}$,

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

$|1-z|>1$

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

Die Stabilitätsfunktion von Runge-Kutta-Verfahren

Runge-Kutta-Verfahren sind vollständig durch die Koeffizienten $A,b,c$ aus ihrem Butcher-Tableau festgelegt. Bei der Testgleichung ist der Anfangswert $y_{0}=1$ und für die Stufen ergibt sich im ersten Zeitschritt

$k_{i}=\lambda \left(1+h\sum _{j=1}^{s}a_{ij}k_{j}\right),\quad i=1,\dotsc ,s.$

Dies ist ein quadratisches lineares Gleichungssystem für den Vektor $k=(k_{1},\dotsc ,k_{s})^{T}$ in der Form $(I-zA)k=\lambda e$ mit dem Vektor $e=(1,\dotsc ,1)^{T}.$ Mit dessen Lösung bekommt man dann die Runge-Kutta-Näherung $y_{1}\approx y(h)$ in der Form

$y_{1}=y_{0}+h\sum _{j=1}^{s}b_{j}k_{j}=1+hb^{T}k=1+zb^{T}(I-zA)^{-1}e=:R(z).$

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

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

$R(z)=1+zb^{T}(I-zA)^{-1}e=1+zb^{T}e+z^{2}b^{T}Ae+\dotsb +z^{s}b^{T}A^{s-1}e.$

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 $\sum _{j=0}^{m}\alpha _{j}y_{n-j}=h\sum _{j=0}^{m}\beta _{j}f(y_{n-j})$ auf die Testgleichung an, ergibt sich wieder mit $z=h\lambda$ die Gleichung

$\sum _{j=0}^{m}\alpha _{j}y_{n-j}-z\sum _{j=0}^{m}\beta _{j}y_{n-j}=\sum _{j=0}^{m}(\alpha _{j}-z\beta _{j})y_{n-j}=0.$

Dies ist eine lineare Differenzengleichung, die man einfach lösen kann. Denn die Folge $y_{n}=u^{n}$ ist eine nichttriviale Lösung dieser Differenzengleichung, wenn u eine Nullstelle des charakteristischen Polynoms

$0{\stackrel {!}{=}}\sum _{j=0}^{m}\alpha _{j}u^{m-j}-z\sum _{j=0}^{m}\beta _{j}u^{m-j}=\varrho (u)-z\sigma (u)$

ist, wobei man die Polynome

$\varrho (u)=\sum _{j=0}^{m}\alpha _{j}u^{m-j}$

$\sigma (u)=\sum _{j=0}^{m}\beta _{j}u^{m-j}$

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

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 $S$. Der Rand dieses Gebietes besteht aus denjenigen $z\in \mathbb {C}$, bei dem für die Nullstellen $|u|=1$ gilt, wo die Nullstellen also auf dem komplexen Einheitskreis liegen. Da $\varrho (u)-z\sigma (u)=0\Leftarrow z=\varrho (u)/\sigma (u)$ gilt, ist die Bestimmung des Stabilitätsgebiets bei Mehrschrittverfahren sogar besonders einfach, denn seinen Rand erhält man i. W. explizit durch

$\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 \}}.$

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 $Y$ der Stufenlösungen gilt

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

und der Zeitschritt wird daher zu

$y^{[n]}=zBY+Vy^{[n-1]}=(V+zB(I-zA)^{-1}U{\big )}y^{[n-1]}.$

In jedem Zeitschritt erfolgt also die Multiplikation mit derselben Matrix

$M(z)=V+zB(I-zA)^{-1}U.$

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

$y'(t)=Qy(t),\quad y(0)=y_{0},\quad Q\in \mathbb {R} ^{d\times d}.$

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

$y_{n}={\big (}R(hQ){\big )}^{n}y_{0}=P{\big (}R(hJ){\big )}^{n}P^{-1}y_{0}.$

Bei einer diagonalisierbaren Matrix $Q$ ist, ist $R(hJ)$ eine Diagonalmatrix mit den Diagonalelementen $R(h\lambda _{j})$. Wenn für alle Eigenwerte $\lambda _{j}$ von $Q$ gilt, dass $h\lambda _{j}\in S$ ist, dann konvergiert auch hier $y_{n}\to 0\,(n\to \infty )$. Bei dieser Differentialgleichung sieht man gleichzeitig, dass es sinnvoll ist, $S$ als offene Menge zu definieren. Denn im diagonalisierbaren Fall bleiben zwar Lösungen auf dem Rand mit $h\lambda _{j}\in \partial S$ 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.