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 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 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 die ganze linke komplexe Halbebene umfasst, heißt das Verfahren A-stabil. Dann ist der Betrag von 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 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 mit $|1+z|<1$ , was dem Inneren des Kreises mit Mittelpunkt und Radius 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 und Radius 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 eine strikt untere Dreiecksmatrix, daher bricht die Neumann-Reihe von $(I-zA)^{-1}$ nach s 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 abhängenden Nullstellen des Polynoms $\varrho (u)-z\sigma (u)$ die Lösungen $u^{n}$ zur Testgleichung und daher liegt 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 . 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 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 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 ist, ist $R(hJ)$ eine Diagonalmatrix mit den Diagonalelementen $R(h\lambda _{j})$ . Wenn für alle Eigenwerte $\lambda _{j}$ von 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, 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.

Basierend auf einem Artikel in:

Wikipedia.de