v1.4 – Absolute Stabilität

 

1.4.1 Dahlquistsche Testgleichung

Die Dahlquistsche Testgleichung lautet:

\dot y = \lambda y,\quad \quad y\left( 0 \right) = 1

mit einem Parameter \lambda \in \mathbb{C}. Für die exakte Lösung y\left( t \right) = {e^{\lambda t}} gilt:

\mathop {\lim }\limits_{t \to \infty } \left| {y\left( t \right)} \right|\quad = \quad \left\{ {\begin{array}{*{20}{c}}  0 & {\operatorname{Re} \left\{ \lambda \right\} < 0} \\  \infty & {\operatorname{Re} \left\{ \lambda \right\} > 0} \\   \end{array} } \right.

Eine wünschenswerte Eigenschaft ist also, dass für \operatorname{Re} \left\{ \lambda \right\} < 0 gilt:

\left| {{u_n}} \right| \to 0 für {t_n} \to \infty

Ein numerisches Verfahren heißt für z = t\lambda absolut stabil, wenn diese Bedingung erfüllt ist.

1.4.2 Beispiel: Explizites Eulerverfahren

Für das explizite Eulerverfahren gilt:

{u_{k+1}} = {u_k}+{\tau _k}\Phi \left( {{t_k},{u_k},\tau ,f} \right)

\quad = {u_k}+{\tau _k}f\left( {{t_k},{u_k}} \right)

\quad = {u_k}+{\tau _k}\lambda {u_k}

\quad = {u_k}\left( {1+{\tau _k}\lambda } \right)

Bei einheitlicher Schrittweite ergibt sich daraus:

{u_{k+1}} = {u_k}\left( {1+\tau \lambda } \right)

\quad = {u_0}{\left( {1+\tau \lambda } \right)^k}

Für y\left( 0 \right) = 1\quad \Rightarrow \quad {u_0} = 1 folgt:

{u_{k+1}} = {u_0}{\left( {1+\tau \lambda } \right)^k} = {\left( {1+\tau \lambda } \right)^k}

Das Verfahren ist also genau dann stabil, wenn \left| {1+\tau \lambda } \right| < 1 ist. Daraus ergibt sich die Stabilitätsbedingung für das explizite Eulerverfahren:

0 < \tau < -\frac{{2\operatorname{Re} \left\{ \lambda \right\}}}{{{{\left| \lambda \right|}^2}}}

num-104-euler-polygonzug-absolut-stabil-bereich

Man sieht, dass diese Bedingung besonders ungünstig ist, wenn \operatorname{Im} \left\{ \lambda \right\} \gg \operatorname{Re} \left\{ \lambda \right\} oder wenn \left| \lambda \right| \gg 1.

1.4.3 Bereich der absoluten Stabilität, unbedingte absolute Stabilität

Der Bereich der absoluten Stabilität ist die Menge

A = \left\{ {z = \tau \lambda \in \mathbb{C}:\left| {{u_n}} \right|\xrightarrow{{{t_n} \to \infty }}0} \right\}

Das heißt A ist die Menge der Werte des Produktes \tau \lambda, für die die numerische Methode Lösungen liefert, die für {t_n} \to \infty gegen 0 gehen.

Ein numerisches Verfahren zur Approximation der Dahlquistschen Testgleichung heißt unbedingt absolut stabil, wenn \left| {{u_n}} \right|\xrightarrow{{{t_n} \to \infty }}0 für beliebige \operatorname{Re} \left\{ \lambda \right\} < 0 und beliebige Schrittweiten \tau > 0 gilt.
Eine Methode ist unbedingt absolut stabil, wenn {\mathbb{C}^-} = \left\{ {z \in \mathbb{C}:\operatorname{Re} \left\{ z \right\} < 0} \right\} \subset A gilt.

1.4.4 Beispiele für Stabilität und Instabilität

1.4.4.1 Numerische Tests

Sei \lambda = -5,\:\:\:T = 8. Dann lautet die Stabilitätsbedingung:

0 < \tau < -\frac{{2\operatorname{Re} \left\{ \lambda \right\}}}{{{{\left| \lambda \right|}^2}}} = -\frac{{2 \cdot \left( {-5} \right)}}{{{{\left| {-5} \right|}^2}}} = \frac{{10}}{{25}} = 0,4

Im numerischen Test beobachtet man folgendes Verhalten:

  • \tau = 0,41\quad \Rightarrow \quadLösung schaukelt sich auf
  • \tau = 0,40\quad \Rightarrow \quadLösung oszilliert unverändert
  • \tau = 0,39\quad \Rightarrow \quadLösung fällt
  • \tau = 0,15\quad \Rightarrow \quadsinnvolle Lösung

1.4.4.2 Anfangswertaufgabe

Die Anfangswertaufgabe

\dot y\left( t \right) = \lambda \left( {y-\cos \left( t \right)} \right)-\sin \left( t \right),\quad \quad y\left( 0 \right) = 1

hat für beliebige \lambda die Lösung

y\left( t \right) = \cos \left( t \right)

Trotzdem diktiert uns die Größe von \lambda = \frac{{\partial f}}{{\partial y}} die maximal verwendbare Schrittweite.
\lambda = \frac{{\partial f}}{{\partial y}} gilt, da \dot y = \lambda y als Linearisierung von \dot y = f\left( {t,y} \right) aufgefasst werden kann. Als Steigung der Geradengleichung muss \lambda daher die Ableitung von f nach y sein.

num-104-anfangswertaufgabe-euler-explizit-stabil-instabil

1.4.4.3 Verfahren von Heun

Wir wollen nun den Bereich der absoluten Stabilität für das Verfahren 2. Ordnung von Heun bestimmen. Erinnerung: Das Verfahren lautet:

\begin{array}{*{20}{c}}  1 &\vline & 1 & {} \\ \hline{} &\vline & {\frac{1}{2}} & {\frac{1}{2}} \\   \end{array}

\Phi = \frac{1}{2}{f_{1,k}}+\frac{1}{2}{f_{2,k}}

{f_{1,k}} = f\left( {{t_k},{u_k}} \right) = \lambda {u_k}

{f_{2,k}} = f\left( {{t_k}+1 \cdot \tau ,\:\:{u_k}+1 \cdot \tau {f_{1,k}}} \right) = \lambda {u_k}+\lambda \tau {f_{1,k}} = \lambda {u_k}+{\lambda ^2}\tau {u_k}

{u_{k+1}} = {u_k}+\tau \left[ {\frac{1}{2}\lambda {u_k}+\frac{1}{2}\left( {\lambda {u_k}+{\lambda ^2}\tau {u_k}} \right)} \right] = {u_k}\underbrace {\left( {1+\tau \lambda +\frac{1}{2}{{\left( {\tau \lambda } \right)}^2}} \right)}_{ = :R\left( {\tau \lambda } \right)}

1.4.4.4 Verfahren von Collatz

Dieses Verfahren ist ebenfalls 2. Ordnung. Die Vorschrift lautet:

\begin{array}{*{20}{c}}{\frac{1}{2}} &\vline & {\frac{1}{2}} & {} \\ \hline{} &\vline & 0 & 1 \\   \end{array}

\Phi = {f_{2,k}}

{f_{1,k}} = f\left( {{t_k},{u_k}} \right) = \lambda {u_k}

{f_{2,k}} = f\left( {{t_k}+\frac{\tau }{2},\:\:{u_k}+\frac{\tau }{2}{f_{1,k}}} \right) = \lambda \left( {{u_k}+\frac{\tau }{2}{f_{1,k}}} \right) = \lambda {u_k}+\frac{\tau }{2}{\lambda ^2}{u_k}

{u_{k+1}} = {u_k}+\tau {f_{2,k}} = {u_k}+\tau \left( {\lambda {u_k}+\frac{\tau }{2}{\lambda ^2}{u_k}} \right) = {u_k}\underbrace {\left( {1+\tau \lambda +\frac{1}{2}{{\left( {\tau \lambda } \right)}^2}} \right)}_{ = :R\left( {\tau \lambda } \right)}

Für beide Verfahren ist also

R\left( z \right) = 1+z+\frac{1}{2}{z^2}

Dies gilt allgemein für alle Verfahren 2. Ordnung. Für explizite Einschrittverfahren 3. Ordnung ergibt sich:

R\left( z \right) = 1+z+\frac{1}{2}{z^2}+\frac{1}{6}{z^3}

Die folgende Grafik stellt die Ränder der Stabilitätsbereiche für explizite Einschrittverfahren 1. (grün), 2. (türkis), 3. (blau) und 4. Ordnung (rot) dar:

num-104-stabilitatsbereiche-explizite-einschrittverfahren

Die bisher besprochenen Verfahren höherer Ordnung haben etwas größere, aber ähnliche Stabilitätsbereiche. Von unbedingter Stabilität ist man weit entfernt, hier helfen nur implizite Verfahren.

1.4.4.5 Implizites Eulerverfahren

Das implizite Eulerverfahren (Euler rückwärts) wird durch die Vorschrift

{u_0} = {y_0},\quad \quad {u_{k+1}} = {u_k}+\tau f\left( {{t_{k+1}},{u_{k+1}}} \right)

definiert. Es muss in jedem Zeitschritt ein Gleichungssystem gelöst werden. Wenn die Funktion f\left( {t,y} \right) nichtlinear von y abhängt, ist das Gleichungssystem nichtlinear. Implizite Verfahren sind also im Allgemeinen aufwändiger als explizite Verfahren.

Zur Analyse der absoluten Stabilität rechnen wir (mit {u_0} = 1):

{u_{k+1}} = {u_k}+\tau f\left( {{t_{k+1}},{u_{k+1}}} \right) = {u_k}+\tau \lambda {u_{k+1}}

\quad \Rightarrow \quad {u_{k+1}} = \frac{{{u_k}}}{{1-\tau \lambda }} = \frac{{{u_0}}}{{{{\left( {1-\tau \lambda } \right)}^k}}}

Es muss also gelten:

\left| {1-\tau \lambda } \right| > 1

num-104-euler-implizit-absolute-stabilitat-bereich

Das Verfahren ist unbedingt absolut stabil. Die A-Stabilität ist eine sinnvolle Eigenschaft für Aufgaben mit \operatorname{Re} \left\{ \lambda \right\} < 0, aber nicht für Aufgaben mit \operatorname{Re} \left\{ \lambda \right\} > 0.

1.4.5 Systeme von gewöhnlichen DGL

Bei Systemen von gewöhnlichen Differentialgleichungen

\dot{ \vec{ x }}= K\vec x,\quad \quad \vec x\left( {{t_0}} \right) = {\vec x_0}

mit einer diagonalisierbaren n \times n-Matrix K = X\Lambda {X^{-1}} gilt

{X^{-1}}\dot{ \vec{ x }}= \Lambda {X^{-1}}\vec x,\quad \quad \vec x\left( {{t_0}} \right) = {\vec x_0}

Für \vec y = {X^{-1}}\vec x gilt also:

\dot{ \vec{ y}} = \Lambda \vec y,\quad \quad \vec y\left( {{t_0}} \right) = {X^{-1}}\vec x\left( {{t_0}} \right) = {X^{-1}}{\vec x_0} = :{\vec y_0}

wobei \Lambda = diag\left( {{\lambda _1}, \ldots ,{\lambda _n}} \right) eine Diagonalmatrix ist. Das numerische Verfahren muss also gleichzeitig die Aufgaben {\dot y_i} = {\lambda _i}{y_i} gut approximieren können, d.h. für die Schrittweite muss gelten:

\tau {\lambda _i} \in A\quad \forall {\lambda _i}

1.4.6 Steifes Differentialgleichungssystem

Ein Differentialgleichungssystem

\dot{ \vec{ x }}= K\vec x,\quad \quad \vec x\left( {{t_0}} \right) = {\vec x_0}

bei dem die Eigenwerte \lambda von K, die negativen Realteil haben, weit auseinander liegen, bezeichnet man als steifes Differentialgleichungssystem. Derartige Systeme können nur mit Verfahren mit einem großen Bereich A absoluter Stabilität vernünftig gelöst werden. Explizite Verfahren sind wenig geeignet.

Beispiel: Feder-Dämpfer-Masse-Schwinger mit starker Dämpfung c = 100

\dot y = \frac{1}{m}\left( {\begin{array}{*{20}{c}}  0 & 1 \\{-1} & {-100} \\   \end{array} } \right)y,\quad \quad y\left( 0 \right) = \left( {\begin{array}{*{20}{c}}  0 \\  1 \\   \end{array} } \right)

Die Matrix hat die Eigenwerte {\lambda _1} = -99,99 und {\lambda _2} = -0,01, diese liegen also um etwa 4 Zehnerpotenzen auseinander. Während das implizite Eulerverfahren damit keine Schwierigkeiten hat, muss beim expliziten Eulerverfahren die Schrittweite kleiner als 2/99,99 sein.

num-104-euler-implizit-explizit-vergleich-stabilitat

Ähnliche Artikel

Kommentar verfassen