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^{\prime }(t)=\lambda y(t),y(0)=1$ mit $\lambda \in ℂ$ 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 ℂ$ 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 ℂ:|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}=\dots ={(R(z))}^j{y}_{n-j}=\dots ={(R(z))}^n{y}_0$ und daher gilt

$y_n\stackrel{n\to \mathrm{\infty }}{\to }0⟺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 $-\mathrm{\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⟺y_1=\frac{1}{1-h\lambda }{y}_0$,

also $R(z)=\frac{1}{1-z}$. Das Stabilitätsgebiet ist nun durch die Bedingung $|\frac{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 $\underset{z\to -\mathrm{\infty }}{\lim }\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 (1+h\sum _{j=1}^s{a}_{ij}{k}_j),i=1,\dots ,s.$

Dies ist ein quadratisches lineares Gleichungssystem für den Vektor $k=(k_1,\dots ,k_s{)}^T$ in der Form $(I-zA)k=\lambda e$ mit dem Vektor $e=(1,\dots ,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+h{b}^{T}k=1+z{b}^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+z{b}^T(I-zA{)}^{-1}e=1+z{b}^{T}e+z^2{b}^{T}Ae+\cdots +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 \mathrm{\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 ℂ$, 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=\{\frac{\varrho (u)}{\sigma (u)}:|u|=1\}=\{\frac{\varrho (e^{i\phi })}{\sigma (e^{i\phi })}:\phi \in [0,2\pi )\}.$

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+U{y}^{[n-1]}\Rightarrow Y=(I-zA{)}^{-1}U{y}^{[n-1]}$

und der Zeitschritt wird daher zu

$y^{[n]}=zBY+V{y}^{[n-1]}=(V+zB(I-zA{)}^{-1}U)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 \mathrm{\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^{\prime }(t)=Qy(t),y(0)=y_0,Q\in {ℝ}^{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(=PJ{P}^{-1})$ ist, gilt

$y_n=(R(hQ){)}^n{y}_0=P(R(hJ){)}^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 \mathrm{\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.