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 Crank und Nicholson publiziert. 1968 bewies Heinz-Otto Kreiss weitere zentrale Eigenschaften des Analyseverfahrens.

Der lineare, eindimensionale Fall

Gegeben sei auf einem Intervall [0,L] eine lineare, partielle Differentialgleichung mit konstanten Koeffizienten der Form

u_t + Au_x = 0,

Anfangsdaten u₀ = 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 Δx und Δt für $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 $E_i^n$ zum Zeitpunkt t_n im Diskretisierungspunkt x_i kann nun in eine Fourier-Reihe

$E_i^n = \sum_{-N}^{N}c_j^n e^{ij\Delta x I}$

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

Beispiel

Der einfachste Fall ist die lineare Advektionsgleichung

u_t + au_x = 0,

wobei a eine reelle Zahl ist. Eines der einfachsten vorstellbaren numerischen Verfahren zur Lösung solcher Gleichungen ist das explizite Euler-Verfahren u^{n + 1} = uⁿ − Δtau_x zur Zeitintegration gekoppelt mit zentralen Differenzen auf einem äquidistanten Gitter im Raum. Der zweite Term wird also mittels

$au_x \approx \frac{a}{2\Delta x} (u^n_{i+1} - u^n_{i-1})$

approximiert. Insgesamt ergibt sich das Verfahren

$\frac{u_i^{n+1} - u_i^n}{\Delta t} + \frac{a}{2\Delta x} (u^n_{i+1} - u^n_{i-1})=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 e^ijΔxI mit ϕ = jΔx:

$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

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

Das Verfahren ist L₂-stabil, falls $|G| \leq 1$ für alle ϕ, was hier nicht der Fall ist, da $|G|^2 = 1 + (\frac{a\Delta t}{\Delta x})^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

$au_x \approx \frac{a}{2\Delta x} (u^n_i - u^n_{i-1})$

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

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

also bedingte Stabilität.

Andere Gleichungen

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:

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

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

Weblinks

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