Von-Neumann-Stabilitätsanalyse

Die Von-Neumann-Stabilitätsanalyse (nach John von Neumann), manchmal auch L2-Stabilitätsanalyse, ist das Standardverfahren zur Untersuchung der Stabilität von numerischen Verfahren zur Lösung zeitabhängiger partieller Differentialgleichungen.

Das Verfahren wurde von John von Neumann in Los Alamos im Rahmen des Manhattan-Projekts entwickelt. Während des Krieges wurde die Methode unter Verschluss gehalten und erst 1947 von John Crank und Phyllis Nicolson publiziert. 1968 bewies Heinz-Otto Kreiss weitere zentrale Eigenschaften des Analyseverfahrens.

Gegeben sei auf einem Intervall ${\displaystyle [0,L]}$ eine lineare, partielle Differentialgleichung mit konstanten Koeffizienten der Form

${\displaystyle u_{t}+Au_{x}=0,}$

Anfangsdaten ${\displaystyle u_{0}=u(0)}$ sowie ein numerisches Verfahren zur Lösung. Die Bedingung, dass das Verfahren in der L2-Norm stabil ist, besagt dann, dass der durch das numerische Verfahren produzierte Fehler bei gegebenen Schrittweiten ${\displaystyle \Delta x}$ und ${\displaystyle \Delta t}$ für ${\displaystyle t\to \infty }$ beschränkt bleibt. Der erste Schritt der Von-Neumann-Stabilitätsanalyse besteht nun darin, die Lösung periodisch auf die kompletten reellen Zahlen fortzusetzen.

Der periodische Fehler ${\displaystyle E_{i}^{n}}$ zum Zeitpunkt ${\displaystyle t_{n}}$ im Diskretisierungspunkt ${\displaystyle x_{i}}$ kann nun in eine Fourier-Reihe

${\displaystyle E_{i}^{n}=\sum _{-N}^{N}c_{j}^{n}e^{ij\Delta xI}}$

entwickelt werden. Hierbei bezeichnet ${\displaystyle I}$ die imaginäre Einheit. Das numerische Verfahren definiert dann eine Evolution der Koeffizienten der Fourierreihe in der Zeit mittels einer so genannten Amplifikationsmatrix ${\displaystyle G}$. Die L2-Stabilitätsbedingung reduziert sich dann darauf, dass das numerische Verfahren genau dann stabil ist, wenn der Spektralradius der Amplifikationsmatrix ${\displaystyle \rho (G)}$ betragsmäßig kleiner gleich eins ist.

Der einfachste Fall ist die lineare Advektionsgleichung

${\displaystyle u_{t}+au_{x}=0,}$

wobei ${\displaystyle a}$ eine reelle Zahl ist. Eines der einfachsten vorstellbaren numerischen Verfahren zur Lösung solcher Gleichungen ist das explizite Euler-Verfahren ${\displaystyle u^{n+1}=u^{n}-(\Delta t)au_{x}}$ zur Zeitintegration gekoppelt mit zentralen Differenzen auf einem äquidistanten Gitter im Raum. Der zweite Term wird also mittels

${\displaystyle au_{x}\approx {\frac {a}{2\Delta x}}(u_{i+1}^{n}-u_{i-1}^{n})}$

approximiert. Insgesamt ergibt sich das Verfahren

${\displaystyle {\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}+{\frac {a}{2\Delta x}}(u_{i+1}^{n}-u_{i-1}^{n})=0,}$

was auch die Entwicklung der Fehler definiert und auch jedes einzelnen Terms der Fourier-Reihenentwicklung. Betrachten wir den j-ten Summanden, so ergibt Einsetzen in die obige Formel und Division durch ${\displaystyle e^{ij\Delta xI}}$ mit ${\displaystyle \phi =j\Delta x}$:

${\displaystyle c_{j}^{n+1}-c_{j}^{n}+{\frac {a\Delta t}{2\Delta x}}c_{j}^{n}(e^{I\phi }-e^{-I\phi })=0.}$

Die Amplifikationsmatrix ist nun gegeben durch

${\displaystyle G(\phi )={\frac {c_{j}^{n+1}}{c_{j}^{n}}}=1-I{\frac {a\Delta t}{\Delta x}}\sin \phi .}$

Das Verfahren ist ${\displaystyle L_{2}}$-stabil, falls ${\displaystyle |G|\leq 1}$ für alle ${\displaystyle \phi }$, was hier nicht der Fall ist, da ${\displaystyle |G|^{2}=1+\left({\frac {a\Delta t}{\Delta x}}\right)^{2}\sin ^{2}\phi .}$ Damit ist das Verfahren unabhängig von der Wahl der Schrittweiten instabil. Dieses Verhalten beobachteten die Mitarbeiter des Manhattan Projekts, was von Neumann zur Entwicklung der Stabilitätsanalyse führte. Wird im Raum das Upwind-Verfahren

${\displaystyle au_{x}\approx {\frac {a}{\Delta x}}(u_{i}^{n}-u_{i-1}^{n})}$

benutzt, so ergibt sich die Courant-Friedrichs-Lewy-Bedingung

${\displaystyle {\frac {a\Delta t}{\Delta x}}<1}$,

also bedingte Stabilität.

Im nichtlinearen Fall oder im Falle variabler Koeffizienten kann das Verfahren durch Linearisierung und Einfrieren der Koeffizienten angewandt werden, allerdings liefert die Analyse im allgemeinen Fall nur noch eine notwendige Bedingung für Stabilität, in Spezialfällen auch eine hinreichende. Ferner ist die Bedingung an den Spektralradius nun:

${\displaystyle \rho (G)\leq 1+c\Delta t}$.

Ein allgemeines Verfahren zur vollständigen Stabilitätsanalyse nichtlinearer Gleichungen ist nicht bekannt.

• Martin Neumann: Erklärung der Analyse in einem Skript Computational Physics I: Grundlagen der Uni Wien