Theorie und Numerik von Erhaltungsgleichungen

Aus Wikibooks

Wechseln zu: Navigation, Suche

Dieses Buch steht im Regal Mathematik

[Bearbeiten] Zusammenfassung des Projekts

Zielgruppe: Studierende der Mathematik, Physik oder des Maschinenbaus ab dem 5. Semester mit grundlegenden Kenntnissen in Numerik.

Lernziele: Grundlegende theoretische Aussagen zu Erhaltungungsgleichungen und der zugehörigen Numerik mit Schwerpunkt nichtlinearer Gleichungen der Strömungsmechanik. Darüberhinaus grundlegendes Wissen zu Finite-Volumen-Verfahren und nichtlinearen Lösungsverfahren.

Buchpatenschaft / Ansprechperson: Benutzer:P. Birken

Sind Co-Autoren gegenwärtig erwünscht? Nein, aber die Korrektur von Tippfehlern oder Anmerkungen zu Fehlern auf der Diskussionsseite ist erwünscht.

Projektumfang und Abgrenzung zu anderen Wikibooks: Nein.

Aufbau des Buches:

Was sind Erhaltungsgleichungen?
Grundlegende Analysis nichtlinearer Erhaltungsgleichungen
Erste Aussagen zur Numerik
Finite-Volumen-Verfahren
Speziellere Aussagen zu Finite-Volumen-Verfahren
Implizite Zeitintegration
Lösung nichtlinearer Gleichungssysteme

Inhaltsverzeichnis

1 Zusammenfassung des Projekts
2 Literatur
3 Erhaltungssätze
- 3.1 Erhaltungsgleichungen
- 3.2 Beispiele für Erhaltungsgleichungen
4 Analysis von Erhaltungsgleichungen
- 4.1 Entropielösungen
- 4.2 Die Lösung der verschwindenden Viskosität
5 Erste Aussagen zur Numerik
6 Zeitintegration
7 Nichtlineare Gleichungssystemlöser
- 7.1 Newton-Verfahren
- 7.2 Mehrgitter-Verfahren
8 Fluid-Struktur-Interaktion
- 8.1 Modellierung
- 8.2 Kopplungsverfahren

[Bearbeiten] Literatur

Als Literatur ist zu empfehlen:

C. Hirsch: Numerical computation of internal and external flows, 2 Bände, 1988, Wiley & Sons, New York. Dieses Buch beschreibt umfassend den Stand der Wissenschaft zum damaligen Zeitpunkt.
Randall J. LeVeque: Numerical Methods for Conservation Laws: Lectures in Mathematics, 1992, Birkhäuser Verlag. Dies ist ein recht kurzes Buch, dass aus einer Vorlesung LeVeques an der ETH Zürich entstanden ist und stellt die vermutlich beste Einführung in das Thema dar.
Randall J. LeVeque: Finite Volume Methods for Hyperbolic Problems, 2002, Cambridge Texts in Applied Mathematics, ISBN 0-521-00924-3. Zehn Jahre später das ganze nochmal, aber in ausführlich.
P. Wesseling: Principles of Computational Fluid Dynamics, 2000, Springer Verlag. Wesseling geht sehr gewissenhaft vor und behandelt Strömungsmechanik und Verfahren für inkompressible und kompressible Strömungen.

[Bearbeiten] Erhaltungssätze

Erhaltungssätze sind physikalische Gesetze, die aussagen, dass eine bestimmte Größe, die Erhaltungsgröße, in einem abgeschlossenen System konstant ist, also erhalten bleibt. Sie spielen eine zentrale Rolle in vielen Bereichen der Physik. Neben starker experimenteller Absicherung lassen sie sich nach dem Noether-Theorem mit bestimmten grundlegenden Symmetrien in Verbindung bringen und sind somit auch theoretisch stark fundiert. Etwa gehört zur Energieerhaltung die Translationsinvarianz physikalischer Aussagen in der Zeit. Beispiele sind

Energieerhaltung
Impulserhaltung
Drehimpulserhaltung
Erhaltung von elektrischer Ladung

Kein echter Erhaltungssatz ist die Massenerhaltung. Diese gilt nur bei nichtrelativistischen Geschwindigkeiten; bei hohen Geschwindigkeiten muss die spezielle Relativitätstheorie in Betracht gezogen werden. Die berühmte Formel

E = m c 2

wird in der Physik als Teil der Energieerhaltung betrachtet: Masse ist dann eine spezielle Form der Erhaltungsgröße „totale Energie“.

Keine Erhaltungsgrößen sind beispielsweise Druck, Temperatur oder elektrische Spannung.

[Bearbeiten] Erhaltungsgleichungen

Die mathematische Formulierung eines Erhaltungssatzes wird als Erhaltungsgleichung bezeichnet. Hier werden wir uns mit Strömungsmechanik beschäftigen. Eine der großen wissenschaftlichen Leistungen des 19. Jahrhunderts war die Erkenntnis, dass sich Strömungen durch Kopplung der Erhaltungsgleichungen für Masse, Energie und Impuls, zusammen mit einer Zustandsgleichung, komplett beschreiben lassen. Ausgangspunkt ist dabei die integrale Form. Wir betrachten als Beispiel die Massenerhaltung. Gegeben sei ein Volumen $Ω$ . Die darin enthaltene Masse $M$ ist dann gegeben durch das Integral der Dichte $ρ$ über $Ω$ :

M =	∫	ρdΩ.
	Ω

Aussage des Erhaltungssatzes der Masse ist, dass Masse nicht zerstört oder erschaffen wird. Änderungen der Masse M im Volumen $Ω$ mit der Zeit sind also Resultat von Massentransport über den Rand von $Ω$ . Dies lässt sich formulieren als

$M_t = \frac{d}{dt} \int_{\Omega} \rho d\Omega = -\int_{\partial \Omega} f \cdot n dS,$

wobei $f$ die Massenflussdichte beschreibt. Diese lässt sich präzisieren als $f = ρ v$ , wobei $v$ die Geschwindigkeit eines kleinen Massenvolumens beschreibt. Damit ergibt sich

$M_t = \frac{d}{dt} \int_{\Omega} \rho d\Omega = -\int_{\partial \Omega} \rho v \cdot n dS.$

Zu beachten ist hierbei, dass dabei außer der Integrierbarkeit von $ρ$ und $ρ v$ keine Voraussetzungen an die Funktionen $ρ$ und $v$ oder an das Volumen $Ω$ gestellt wurden. Eine zweite wichtige Form einer Erhaltungsgleichung ist die differentielle Form. Dazu wird die letzte Gleichung mittels des Integralsatzes von Gauß umgeformt:

$\frac{d}{dt} \int_{\Omega} \rho d\Omega + \int_{\Omega} \nabla \cdot \rho v d\Omega=0.$

Ist das Volumen $Ω$ nicht zeitabhängig, lässt sich die Zeitdifferentiation in das Integral reinziehen. Da die Aussage für beliebige Volumina gilt, ergibt sich die differentielle Form der Massenerhaltung

$\rho_t + \nabla \cdot \rho v=0.$

Allgemeiner haben Erhaltungsgleichungen die Form

$u_t + \nabla \cdot f(u, \nabla u)=0$

oder, wenn zusätzlich Volumenkräfte modelliert werden müssen,

$u_t + \nabla \cdot f(u, \nabla u)=g(u).$

Letztere Form wird im Englischen auch als Balance law bezeichnet (deutsch: Bilanzgleichung).

Eine dritte wichtige Form der Erhaltungsgleichung ist die schwache Form. Dazu multiplizieren wir die differentielle Form (hier für den eindimensionalen Fall) mit glatten Testfunktionen $\phi \in C_0^1(\mathbb{R}\times \mathbb{R})$ , also aus dem Raum der stetig differenzierbaren Funktionen mit kompaktem Träger. Integration liefert

$\int_0^{\infty}\int_{-\infty}^{\infty}[\phi u_t + \phi f(u)_x]dxdt=0.$

Partielle Integration ergibt aufgrund des kompakten Trägers von $φ$

$\int_0^{\infty}\int_{-\infty}^{\infty}[\phi_t u + \phi_xf(u)]dxdt=-\int_{-\infty}^{\infty}\phi(x,0)u(x,0).$

Eine Funktion $u (x, t)$ wird dann schwache Lösung der Erhaltungsgleichung genannt, wenn die obige Gleichung für alle $\phi \in C_0^1(\mathbb{R}\times \mathbb{R}^+)$ erfüllt ist.

[Bearbeiten] Beispiele für Erhaltungsgleichungen

Das wichtigste Beispiel der Strömungsmechanik ist das System der Navier-Stokes-Gleichungen

$\partial_t \rho + \nabla \cdot {\textbf m} = 0,$

$\partial_t m_i + \sum_{j=1}^3 \partial_{x_j} (m_i v_j + p \delta_{ij}) = \sum_{j=1}^3 \partial_{x_j} S_{ij}, \qquad i=1,2,3$

$\partial_t \rho E + \nabla \cdot (H {\textbf m}) = \sum_{j=1}^3 \partial_{x_j} \left (\sum_{i=1}^3 S_{ij} v_i - W_j \right) + q.$

Als Vereinfachung ergeben sich bei Vernachlässigung der Terme 2. Ordnung die Euler-Gleichungen, bei denen die Flüsse $f$ nicht von $\nabla u$ abhängen, also ein System erster Ordnung:

$\partial_t \rho + \nabla \cdot {\textbf m} = 0,$

$\partial_t m_i + \sum_{j=1}^3 \partial_{x_j} (m_i v_j + p \delta_{ij}) = 0, \qquad i=1,2,3$

$\partial_t \rho E + \nabla \cdot (H {\textbf m}) = 0.$

[Bearbeiten] Analysis von Erhaltungsgleichungen

Erhaltungsgleichungen in differentieller Form sind partielle Differentialgleichungen (PDEs) und lassen sich somit mit deren Maschinerie analysieren. Dazu ist es notwendig, verschiedene Typen von PDEs auseinanderzuhalten. Zunächst bezeichnet man als Ordnung einer PDE den höchsten Grad der auftretenden Ableitungen. Man spricht von einem System, wenn $u$ ein Vektor ist, ansonsten von einer skalaren Gleichung. Sind $f$ und $g$ linear heißt die Gleichung linear, ansonsten nichtlinear. Die Navier-Stokes-Gleichungen sind also ein nichtlineares System zweiter Ordnung, die Euler-Gleichung ein nichtlineares System erster Ordnung. Eine Gleichung heißt instationär, wenn sie eine Zeitabhängigkeit enthält und stationär, wenn das nicht der Fall ist. Sie heißt mehrdimensional, wenn es mehrere Ortsvariablen gibt.

Schließlich teilt man die Gleichungen noch in hyperbolisch, elliptisch, parabolisch und gemischte Formen ein. Grob kann man sagen, dass diese sich durch die Arten der erlaubten Randbedingungen unterscheiden.

Elliptische Gleichungen erlauben nur Randbedingungen, keine Anfangsbedingungen. Damit gibt es für elliptische Gleichungen keine Variable, die man einer Zeit zuordnen könnte; Sie entsprechen stationären Zuständen und Störungen der Nebenbedingungen machen sich sofort in der gesamten Lösung bemerkbar.

Parabolische Gleichungen erlauben Anfangs- und Randbedingungen. Sie haben eine Zeitabhängigkeit, Störungen der Nebenbedingungen können sich sofort in der kompletten Lösung bemerkbar machen. Ferner haben parabolische Gleichungen eine Glättungseigenschaft: Sind die Anfangsdaten zum Zeitpunkt Null endlich oft stetig differenzierbar, ist es möglich, dass die Lösung für Zeiten größer Null unendlich oft stetig differenzierbar ist.

Hyperbolische Gleichungen erlauben ebenfalls Anfangs- und Randbedingungen. Sie haben eine Zeitabhängigkeit, Störungen der Nebenbedingungen breiten sich mit endlicher Geschwindigkeit in ausgezeichnete Richtungen aus. Im Gegensatz zu parabolischen Gleichungen, die Anfangsdaten glätten, können hier auch bei beliebig glatten Anfangsdaten unstetige Lösungen auftreten.

Wir werden hier vor allem hyperbolische Erhaltungsgleichungen betrachten und nebenbei noch parabolische. Aus Gründen der Einfachheit untersuchen wir zunächst ausschließlich eindimensionale Probleme. Dieses Vorgehensweise ist deswegen sinnvoll, weil ein numerisches Verfahren, was schon im eindimensionalen nicht funktioniert, in mehreren Dimensionen erst recht nicht funktioniert. Zunächst betrachten wir die skalare nichtlineare hyperbolische Erhaltungsgleichung mit $f(u) \in C^1(\mathbb{R})$ :

u t + f (u) x = 0

auf dem Gebiet $\mathbb{R}\times\mathbb{R}^+$ und

u (x,0) = u 0 (x)

auf $\mathbb{R}.$

[Bearbeiten] Entropielösungen

Nichtlineare hyperbolische Erhaltungsgleichungen lassen in der Regel mehrere Lösungen zu. Da Erhaltungsgleichungen physikalische Gesetze beschreiben sollen, ist es sinnvoll, aus diesen mit einer Zusatzbedingung jene Lösungen herauszuziehen, die physikalisch sinnvoll sind. Hierbei hat sich die Entropiebedingung als nützlich erwiesen. Nach dem zweiten Hauptsatz der Thermodynamik ist die Gesamtentropie eines Systems nichtsinkend. Ferner steigt die Entropie eines Gases wenn man sich über einen physikalischen Schock bewegt. Dies führt zur Entropiebedingung

$\int_{x_1}^{x_2} \eta(u(x,t_2))dx \leq \int_{x_1}^{x_2} \eta(u(x,t_1))dx + \int_{t_1}^{t_2} \psi(u(x_1,t))dt - \int_{t_1}^{t_2} \psi(u(x_2,t))dt.$

$η(u)$ ist die so genannte Entropiefunktion und $ψ(u)$ der Entropiefluss. An dieses Entropie-Entropiefluss-Paare stellt man die Bedingung, dass $η$ glatt und konvex ist, $ψ$ glatt und ferner die Beziehung

ψ'(u) = η'(u) f'(u)

gilt. Dann nennt man eine Lösung, die die obige Ungleichung erfüllt Entropielösung der Erhaltungsgleichung. Analog kann man eine schwache Form herleiten:

$\int_0^{\infty}\int_{-\infty}^{\infty}[\phi_t\eta(u) + \phi_x\psi(u)]dxdt + \int_{-\infty}^{\infty}\phi(x,0)\eta(u(x,0))dx \geq 0.$

Eine schwache Lösung die die obige Gleichung für alle zulässigen Entropiepaare erfüllt, heißt schwache Entropielösung der Erhaltungsgleichung. Diese Lösung ist dann sogar eindeutig und hängt in der $L 1$ -Norm stetig von den Anfangsdaten ab!

[Bearbeiten] Die Lösung der verschwindenden Viskosität

Eine alternative Herangehensweise ist die Betrachtung der Gleichung

u t + f (u) x = ε u x x

für $ε > 0$ . Diese Gleichung ist parabolisch und hat eine eindeutige Lösung $u ε$ . Im Grenzwert $\epsilon \rightarrow 0$ ergibt sich die Lösung der verschwindenen Viskosität (vanishing viscosity solution). Diese erfüllt die Entropiebedingung.

[Bearbeiten] Erste Aussagen zur Numerik

Dazu sei ein Gitter aus Gitterzellen $Ω i$ gegeben mit $Ω i = (x i - 1, x i)$ und diskrete Zeitpunkte $t 0, t 1,...$ . Die Schrittweiten beschreiben wir mit $Δ x i = x i - x i - 1$ und $Δ t i = t i + 1 - t i$ . Mit $u_i^n$ bezeichnen wir die numerische Approximation an die exakte Lösung $u (x, t n)$ auf $Ω i$ .

Die zentrale physikalische Eigenschaft die Erhaltungsgleichungen modellieren, ist offensichtlich die Erhaltung. Das numerische Verfahren muss dies berücksichtigen, wenn Hoffnung bestehen soll, eine physikalisch sinnvolle Approximation zu berechnen. Wir werden dies später im Satz von Lax-Wendroff präzisieren, aber zunächst die Definition:

[Bearbeiten] Definition: Konservatives Verfahren

Ein numerisches Verfahren zur Lösung der Erhaltungsgleichung $u t + f (u) x = 0$ heißt konservativ, wenn es die Form

$u_i^{n+1} = u_i^n - \frac{\Delta t_n}{\Delta x_i}[F_{i+1/2} - F_{i-1/2}]$

hat. Die Funktion $F$ wird auch als numerische Flussfunktion bezeichnet. Hintergrund ist, dass die numerische Näherung an die Erhaltungsgröße gegeben ist durch

$\sum_{i=I}^J \Delta x_i u_i^{n+1} = \sum_{i=I}^J \Delta x_i u_i^n - \Delta t_n [F_I - F_J].$

Sind die Flüsse an den Randpunkten also physikalisch sinnvoll, erhält das numerische Verfahren die Erhaltungsgröße. Entscheidend für ein sinnvolles numerisches Verfahren ist die Wahl der Flüsse. Die Definition ist offensichtlich auch erfüllt, wenn wir die Nullfunktion als Fluss wählen, eine gute numerische Approximation ist davon aber nicht zu erwarten. Um dem abzuhelfen, braucht es eine Eigenschaft, die verlangt, dass unsere numerische Flussfunktion $F$ den physikalischen Fluss $f$ approximiert.

[Bearbeiten] Definition: Konsistentes Verfahren

Eine Flussfunktion $F$ mit den Argumenten $u 1$ , $u 2$ bis $u j$ heißt konsistent, wenn die beiden folgenden Eigenschaften erfüllt sind:

$F (u, u,..., u) = f (u)$
$F$ ist Lipschitz-stetig in allen Komponenten.

Ausgestattet mit diesen beiden Eigenschaften lässt sich bereits ein wichtiges Resultat beweisen.

[Bearbeiten] Satz von Lax-Wendroff

Da nichtlineare Erhaltungsgleichungen in der Regel mehrere schwache Lösungen zulassen, ist die Frage der Konvergenz nicht einfach zu beantworten. Es kann sogar sein, dass ein numerisches Verfahren konvergiert, aber nicht gegen eine schwache Lösung der Erhaltungsgleichung. Mit diesem Problem beschäftigt sich der Satz von Lax-Wendroff.

Satz von Lax-Wendroff (Lax und Wendroff, 1960)

Gegeben sei eine Folge von Gittern mit festen Maschenweiten $Δ x l$ und $Δ t l$ mit $\Delta x_l, \Delta t_l \rightarrow 0$ für $l \rightarrow \infty$ . Gegeben sei ein konsistentes und konservatives numerisches Verfahren, die entsprechende numerische Lösung auf Gitter $l$ sei $U l (x, t)$ . Sei ferner die totale Variation von $U l (x, t)$ beschränkt in dem Sinne, dass zu jedem $T$ ein $R > 0$ existiert, so dass

$TV(U_l(\cdot,t)):=\sup \sum_{j=1}^N|U_l(\xi_j,t)-U_l(\xi_{j-1},t)|<R$ für alle $0 \leq t <T, l=1,2,\ldots$

und das Supremum über alle Unterteilungen $(ξ) n$ der reellen Achse genommen wird. Schliesslich konvergiere $U l (x, t)$ gegen eine Funktion $u (x, t)$ in folgendem Sinne: Auf jeder beschränkten Menge $\Omega = [a,b] \times [0,T]$ gelte

$\int_0^T\int_a^b|U_l(x,t) - u(x,t)|dxdt \rightarrow 0$ , für $l \rightarrow \infty .$

Dann ist $u (x, t)$ eine schwache Lösung der Erhaltungsgleichung.

Beweis

Zu zeigen ist, dass $u (x, t)$ die schwache Form der Erhaltungsgleichung erfüllt:

$\int_0^{\infty} \int_{-\infty}^{\infty}[\phi_tu+\phi_xf(u)]dxdt = -\int_{\infty}^{\infty}\phi(x,0)u(x,0)dx$ für alle $\phi \in C_0^1.$

Wir betrachten zunächst das feste Gitter $l$ und vernachlässigen zur Vereinfachung der Notation die zugehörigen Indizes. Die numerische Lösung $U l (x, t)$ ist für alle $Δ x$ , $Δ t$ gegeben durch die konservative Vorschrift

$u_i^{n+1} = u_i^n + \frac{\Delta t}{\Delta x}(F_{i+1/2}^n-F_{i-1/2}^n).$

Dies multiplizieren wir mit einer beliebigen Testfunktion $φ$ und erhalten mit der Bezeichnung $\phi_i^n$ für die Einschränkung von $φ$ auf die Zelle $i$ und das Zeitintervall $[t n, t n + 1]$ :

$\phi_i^n u_i^{n+1} = \phi_i^n u_i^n + \frac{\Delta t}{\Delta x}\phi_i^n (F_{i+1/2}^n-F_{i-1/2}^n).$

Summieren wir die obige Formel über alle $i$ und $n \geq 0$ auf, ergibt sich

$\sum_{n=0}^{\infty} \sum_{i=-\infty}^{\infty} \phi_i^n(u_i^{n+1} - u_i^n) = -\frac{\Delta t}{\Delta x} \sum_{n=0}^{\infty} \sum_{i=-\infty}^{\infty} \phi_i^n (F_{i+1/2}^n-F_{i-1/2}^n).$

Diese Summen schreiben wir nun so um, dass die Differenzen nun von Testfunktionen genommen werden. Hierbei ist zu beachten, dass die Summen nur formal unendlich sind, da die Testfunktionen einen kompakten Träger haben und damit für $|i|\geq i_0$ und für $n \geq n_0$ Null sind. Was wir also tun ist, Summen der Form $\sum_{i=1}^m a_i(b_i-b_{i-1})$ umzuschreiben als $a_mb_m - a_1b_0 - \sum_{i=1}^{m-1} (a_{i+1}-a_i)b_i$ . Damit ergibt sich

$-\sum_{i=-\infty}^{\infty} \phi_i^0u_i^0 - \sum_{n=1}^{\infty} \sum_{i=-\infty}^{\infty} (\phi_i^n - \phi_i^{n-1})u_i^n = \frac{\Delta t}{\Delta x} \sum_{n=0}^{\infty} \sum_{i=-\infty}^{\infty} (\phi_{i+1}^n - \phi_i^n)F_{i-1/2}^n$

und nach weiteren Umformungen

$\Delta t \Delta x\left[\sum_{n=1}^{\infty} \sum_{i=-\infty}^{\infty} \left(\frac{\phi_i^n - \phi_i^{n-1}}{\Delta t}\right)u_i^n + \sum_{n=0}^{\infty} \sum_{i=-\infty}^{\infty} \left(\frac{\phi_{i+1}^n - \phi_i^n}{\Delta x}\right)F_{i-1/2}^n\right] = -\Delta x \sum_{i=-\infty}^{\infty} \phi_i^0u_i^0.$

Diese Gleichung gilt auf jedem Gitter $l$ . Wir betrachten nun den Grenzwert $l \rightarrow \infty$ . Definieren wir Anfangsdaten beispielsweise derart, dass $u_i^0 = 1/\Delta x \int_{x_i-1/2}^{x_i+1/2}u(x,0) dx$ , so konvergiert die rechte Seite aufgrund der Glattheit von $φ$ und der Konvergenz von $U l$ gegen $u$ gegen $-\int_{-\infty}^{\infty} \phi(x,0)u(x,0)dx$ . Der erste Term in der Summe auf der linken Seite konvergiert ebenfalls wegen der Glattheit von $φ$ und der Konvergenz von $U l$ gegen $\int_0^{\infty} \int_{-\infty}^{\infty} \phi_t(x,t)u(x,t)dxdt$ .

Der zweite Term in der Klammer enthält auch die Flussfunktion und ist deswegen etwas weniger glatt und schwieriger. Die Lipschitz-Stetigkeit von $F$ , gemeinsam mit der beschränkten Variation von $U l$ ist allerdings ausreichend, um die Konvergenz von $F$ gegen $f$ fast überall zu garantieren. Damit ergibt sich dann auch für diesen Term Konvergenz gegen das gewünschte Integral $\int_0^{\infty} \int_{-\infty}^{\infty} \phi_x(x,t)f(u(x,t))dxdt$ und damit der Beweis des Satzes.

[Bearbeiten] Entropiebedingung

Das Problem des Satzes von Lax-Wendroff ist, dass die schwache Lösung der Erhaltungsgleichung nicht eindeutig sein muss und er macht keine Aussagen darüber, zu welcher dieser Lösungen man konvergiert. Noch genauer kann es sein, dass unterschiedliche Folgen von Gittern zu unterschiedlichen schwachen Lösungen konvergieren. Ebensowenig sagt er etwas über den Abschneidefehler aus, mit dem wir in Anbetracht von endlich großen Maschenweiten konfrontiert sind oder liefert überhaupt eine Konvergenzaussage. Für das erste Problem ist die Lösung die so genannte Entropiebedingung. Im kontinuierlichen Fall stellt diese eine eindeutige Lösung bereit und auch im diskreten Fall lässt sich eine Entropiebedingung angeben. Gegeben sei also ein Entropiepaar $η$ , $ψ$ und dazu ein diskretes Analogon $Φ$ , das zu $p s i$ konsistent ist (im bekannten Sinne). Dann ist die diskrete Entropiebedingung:

$\eta(u_i^{n+1})\leq \eta(u_i^n) - \frac{\Delta t}{\Delta x}(\Phi^n_{i+1/2}-\Phi^n_{i-1/2}).$

Erfüllt die numerische Flussfunktion diese Bedingung lässt sich der Beweis für den Satz von Lax-Wendroff in analoger Weise durchführen und zeigen, dass die Grenzfunktion dann die schwache Entropielösung ist.

[Bearbeiten] Stabilität

Wir haben nun Bedingungen an ein numerisches Verfahren an der Hand, die im Falle der Konvergenz die Konvergenz gegen eine physikalisch sinnvolle Lösung garantieren. Was fehlt, ist ein Konvergenzbeweis und dafür ist Stabilität des Verfahrens notwendig. Genauer liefert Konsistenz, dass der lokale Abschneidefehler $τ n$ bei verschwindender Zeitschrittweite gegen Null geht. Was zur Konvergenz noch fehlt ist dann eine Eigenschaft die besagt, dass sich die bereits gemachten Fehler nicht aufschaukeln. Dies bedeutet, dass der durch das numerische Verfahren produzierte Fehler bei gegebenen Schrittweiten $Δ x$ und $Δ t$ für $t \to \infty$ in der entsprechenden Norm beschränkt bleibt.

Zur numerischen Lösung $u n$ zum Zeitpunkt $t n$ gehöre der Fehler $E n$ . Das numerische Verfahren $N (u)$ liefert dann eine Näherung $u n + 1 = N (u n)$ zum Zeitpunkt $t n + 1 = t n + Δ t$ mit Fehler $E n + 1$ . Also haben wir

$\begin{align} E^{n+1} &= u^{n+1}-u(x,t^{n+1}) = N(u(x,t^n) + E^n) - u(x,t^{n+1})\\ &= N(u(x,t^n) + E^n) -N(u(x,t^n)) + N(u(x,t^n)) - u(x,t^{n+1})\\ &= [N(u(x,t^n) + E^n) -N(u(x,t^n))] - \Delta t \tau^n.\end{align}$

Der zweite Term ist der lokale Abschneidefehler und Stabilität ist die Eigenschaft, die den ersten Term kontrolliert. Wir betrachten zunächst Fall einer linearen Gleichung und dafür ein lineares numerisches Verfahren. Dann ist

N (u (x, t n) + E n) - N (u (x, t n)) = N (e n)

und es ist ausreichend die Wirkung des numerischen Verfahrens auf den Fehler $e n$ zu betrachen. Bleibt der Fehler beschränkt, lässt sich Konvergenz zeigen (Äuivalenzsatz von Lax).

[Bearbeiten] Von-Neumann-Stabilitätsanalyse

Wir betrachten zunächst Stabilität in der $L 2$ -Norm. Die theoretische Rechtfertigung dieser Form der Stabilitätsanalyse basiert auf linearen räumlich periodischen Problemen, für die sich dann sehr starke Aussagen treffen lassen. Allerdings lässt sich die Technik auch sehr gut auf andere Probleme anwenden und liefert sehr gute praktische Ergebnisse, womit sie das Standardanalyseverfahren für Stabilität von numerischen Verfahren für zeitabhängige PDEs darstellt. Gegeben sei auf dem Intervall $[0, L]$ die Gleichung

u t + a u x = 0

mit $a \neq 0$ .

mit periodischen Randbedingungen und die auf ganz $\mathbb{R}$ periodisch fortgesetzte Lösung $u$ . Das numerische Verfahren definiert dann eine Evolution des Fehlers in der Zeit und lässt sich über eine Amplifikationsmatrix darstellen, deren Dimension allerdings mit kleiner werdender Maschenweite ansteigt. Die Idee ist nun, den periodischen Fehler zum Zeitpunkt $t n$ im Diskretisierungspunkt $x i = i Δ x$ in eine w:Fourier-Reihe

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

zu entwickelt. Dies diagonalisiert die Amplifikationsmatrix und wir können dann die Entwicklung der einzelnen Fourierkoeffizienten in der Zeit analysieren. Nach der parsevalschen Ungleichung gilt $\|e_i\|_{L_2} = \|c_j\|_{L_2}$ . Die $L 2$ -Stabilitätsbedingung reduziert sich dann darauf, dass das numerischer 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 + a u x = 0,

mit $a\in \mathbb{R}$

In der Zeit nehmen wir das explizite Euler-Verfahren $u n + 1 = u n + Δ t (a u x) n$ 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

$\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 i j Δ x I$ 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 2$ -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}{\Delta x} (u^n_i - u^n_{i-1})$

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

$\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. Unter Verwendung eines leicht anderen nichtäquivalenten Stabilitätsbegriffes wird manchmal eine andere Bedingung an STabilität gestellt:

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

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

[Bearbeiten] TV-Stabilität

Im nichtlinearen Fall braucht man in der Regel nichtlineare numerische Verfahren und dann funktioniert die obige Fehlerdarstellung nicht. Ist das numerische Verfahren eine Kontraktion, also gilt

$\|N(P)-N(Q)\|\leq N \|P-Q\|,$

so ergibt sich immerhin

$\|E^{n+1}\|\leq \|N(u(x,t^n) + E^n) -N(u(x,t^n))\| + \Delta t \|\tau^n\| \leq \|E^n\| + \Delta t \|\tau^n\|,$

womit sich Stabilität zeigen lässt. Diese Eigenschaft gilt insbesondere für die später erwähnten monotonen Verfahren, die aber maximal erster Ordnung sind. Entsprechend brauchen wir einen besseren Stabilitätsbegriff, nämlich die TV-Stabilität, bei der es um Konvergenz in der $L 1$ -Norm geht. Hintergrund ist, dass Mengen der Form

$\{v \in L_1: \; TV(v) \leq R$ und $supp(v(\cdot,t)) \subset [-M,M] \forall t\in [0,T]\}$

in $L 1$ kompakt sind.

Definition TV-stabil

Wir nennen ein numerisches Verfahren nun TV-stabil, falls die numerischen Lösungen für alle $Δ t < Δ t 0$ in einer Menge der obigen Form sind mit M und R fix.

Hierzu gilt folgender Satz:

Satz

Gegeben sei ein konservatives numerisches Verfahren mit einem Lipschitz-stetigen numerischen Fluss und für die Anfangsdaten $u 0$ existiere ein $Δ t 0$ und ein $R > 0$ , so dass für alle n und alle $Δ t$ mit $Δ t < Δ t 0$ und $n\Delta t\leq T$ gelte:

$TV(u^n) \leq R.$

Dann ist das numerische Verfahren TV-stabil.

Beweis

Siehe LeVeque, Seite 250.

Hiermit können wir nun ein Konvergenzresultat beweisen:

Satz

Die numerische Lösung U sei von einem TV-stabilen, konservativen und konsistenten numerischen Verfahren mit fester Zeitschrittweite $Δ t$ generiert. Dann konvergiert die Lösung für $\Delta t \rightarrow 0$ in der $L 1$ -Norm gegen eine schwache Lösung der Erhaltungsgleichung.

Beweis

Wir nehmen an, das die Aussage falsch wäre, es gäbe also eine Nullfolge $(Δ t) j$ , so dass die zugehörige Folge von numerischen Lösungen nicht gegen eine schwache Lösung konvergiert. Da das Verfahren TV-stabil ist, liegt die Folge $u (Δ t j)$ der numerischen Lösungen in einer in $L 1$ kompakten Menge und entsprechend existiert eine $L 1$ -konvergente Teilfolge, die die Voraussetzungen des Satzes von Lax-Wendroff erfüllt. Diese konvergiert also gegen eine schwache Lösung der Erhaltungsgleichung.

[Bearbeiten] Monotone Verfahren

Eine Eigenschaft der hier betrachteten Erhaltungsgleichungen ist, dass der Loesungsoperator invers monoton ist: gilt fuer zwei Anfangsdaten $u^0 \geq v^0$ , dann gilt fuer die entsprechenden Loesungen $u(t,\cdot)\geq v(t,\cdot)$ . Ein konservatives Verfahren, das die entsprechende diskrete Eigenschaft hat, nennt man monoton. Genauer, falls fuer die numerische Loesung aus

$u_j^n\geq v_j^n$ fuer alle

j

folgt, dass

$u_j^{n+1}\geq v_j^{n+1}$ fuer alle

j

,

so heisst das Verfahren monoton. Dies ist aequivalent mit

$\frac{\partial H_j}{\partial u_{j+l}}\geq 0$ fuer alle

l

,

wobei $H$ mittels $u_j^{n+1} = H_j(u^n)$ definiert ist.

Monotone Verfahren haben den grossen Vorteil, dass sie umfassende Konvergenzaussagen erlauben. Harten, Hyman und Lax bewiesen 1976, dass monotone Verfahren im Falle der Konvergenz gegen die schwache Entropieloesung konvergieren. 1980 bewiesen dann Crandall und Majda folgendes Resultat:

Satz

Die durch ein monotones Verfahren generierte Loesung konvergiert in $L 1$ gegen die schwache Entropieloesung der Erhaltungsgleichung fuer $Δ t$ gegen Null mit $Δ t / D e l t a x$ fix. Dieses Resultat gilt auch im mehrdimensionalen.

Kroener und Ohlberger bewiesen 2000 eine grundlegende Fehlerabschaetzung fuer monotone Verfahren. Um 2000 wurde von verschiedenen Autoren folgendes mehrdimensionale Resultat erarbeitet, was es auf der einen Seite tatsaechlich erlaubt, a priori Fehlerschaetzer zu erarbeiten, aber auf der anderen Seite das grosse Problem monotoner Verfahren vor Augen fuehrt.

Satz

Sei $u 0 (x)$ in $BV(\mathbb{R}^d)$ und $u (x, t)$ eine schwache Entropieloesung der Erhaltungsgleichung. Sei $U$ eine durch ein monotones numerisches Verfahren mit dem expliziten Eulerverfahren als Zeitintegration generierte numerische Loesung, mit einem in bestimmter Weise beschraenkten Zeitschritt. Sei $K$ eine bestimmte Teilmenge von $\mathbb{R}^d \times \mathbb{R}^+$ . Dann existiert eine Konstante $C > 0$ , so dass

$\|u-U\|_{L_1(K)} \leq Ch^{1/4}.$

Im eindimensionalen Fall gilt die Aussage mit $h 1 / 2$ , diese Konvergenzrate ist optimal.

Numerisch beobachtet man haeufig in O(h)-Verhalten, nichtsdestotrotz ist diese Konvergenzrate zu niedrig. Der Grund fuer die niedrige Rate liegt in der Eigenschaft der Monotonie. Diese ist sehr weitgehend, was die starken Konvergenzrate erst erlaubt, allerdings auch eine extrem starke Einschraenkung an die Verfahren. Ein monotones Verfahren hat nebenbei noch die folgenden Eigenschaften:

Die $l 1$ -Norm der numerischen Loesung ist mit der Zeit nichtsteigend ( $l 1$ -kontraktiv).
Damit gilt auch: Die Totale Variation ist nicht steigend (TVD).
Und daraus folgt: Monotonie der Loesung bleibt erhalten (Monotonieerhaltend).

Diese Eigenschaften bilden eine Hierarchie von Verfahrensklassen, die echte Teilmengen bilden. Es liegt also auf der Hand, die Anforderungen an das Verfahren abzuschwaechen, um hoehere Konvergenzordnung zu erzielen.

[Bearbeiten] Monotonieerhaltende Verfahren

Die Lösung $u (x, t)$ einer skalaren hyperbolischen Erhaltungsgleichung hat die Eigenschaft, dass aus Monotonie der Anfangsdaten Monotonie der Funktionen $u(x,\cdot)$ für alle Zeiten folgt. Das diskrete Analogon davon sind die monotonieerhaltenden Verfahren, die bei monotonen Anfangsdaten $(u_j^0)_j$ monotone Näherungslösungen $(u_j^n)_j$ für alle $n$ liefern.

Godunow bewies in den 1950er Jahren, dass lineare Einschrittverfahren zweiter Ordnung für die lineare Advektionsgleichung nur in uninteressanten Spezialfällen monotonieerhaltend sein können. Später wurde die Aussage von anderen auch für lineare Mehrschrittverfahren bewiesen. Wir zeigen den Beweis von Godunow.

Jedes lineare Einschrittverfahren lässt sich in der Form

$u_j^{n+1} = \sum_m \gamma_m u_{j+m}^n$

schreiben (das Verfahren ist eine lineare Abbildung von $(u_j^n)_j$ auf $(u_j^{n+1})_j$ . Zunächst ein Hilfssatz

Satz

Das obige Verfahren ist genau dann monotonieerhaltend, wenn

$\gamma_m \geq 0$ fuer alle

m

.

Beweis

ObdA sei $n = 0$ und $(u_j^n)_j$ nichtfallend.

Hinrichtung: Sei das Verfahren monotonieerhaltend, aber es existiere ein p mit $γ p < 0.$ Wir wählen nun Anfangsdaten:

$u_j^0=0$ , $j \leq 0$ und $u_j^0=1$ ,

j > 0.

Dann ist $u_{1-p}^1-u_{-p}^1 = \gamma_p<0$ und damit $(u_j^1)_j$ nicht nichtfallend, was einen Widerspruch ergibt.

Rückrichtung: Seien alle $\gamma_m \geq 0$ . Dann ist für alle j:

$u_j^1-u_{j-1}^1 = \sum_m\gamma_m (u_{j+m}^0-u_{j+m-1}^0)\geq 0.$

also ist das Verfahren monotonieerhaltend.

[Bearbeiten] Satz von Godunow

Wir können nun den Satz von Godunow beweisen:

Lineare Einschrittverfahren zweiter Ordnung für die lineare Advektionsgleichung können nicht monotonieerhaltend sein, es sei denn $|c|\Delta t/\Delta x \in \mathbb{N}.$

Beweis

Sei $u(0,x) = (\frac{x}{\Delta x}-\frac{1}{2})^2 - \frac{1}{4}$ und die entsprechenden numerischen Anfangsdaten $u^0_j = (j-\frac{1}{2})^2 - \frac{1}{4}$ . Die exakte Lösung ist ein Polynom zweiten Grades und ergibt sich durch Advektion der Anfangsdaten:

$u(x,t) = (\frac{x-ct}{\Delta x} - \frac{1}{2})^2 - \frac{1}{4}.$

In Godunows Sinne bedeutet, dass ein Verfahren zweiter Ordnung diese Lösung exakt reproduziert, es muss also gelten:

$u_j^1 = \sum_m \gamma_m [(j+m-\frac{1}{2})^2-\frac{1}{4}]= (j-c\Delta t/Delta x - \frac{1}{2})^2 - \frac{1}{4}.$

Nehmen wir an, das Verfahren sei Monotonieerhaltend, also $\gamma_m \geq 0.$ Wegen $(j+m-\frac{1}{2})^2-\frac{1}{4} \geq 0$ folgt

$(j-c\Delta t/\Delta x - \frac{1}{2})^2 - \frac{1}{4} \geq 0$ fuer alle j.

Fuer $c\Delta t/\Delta x \in \mathbb{N}$ und groesser Null waehlen wir nun j mit $j > c Δ t / Δ x > j - 1$ . Dann gilt

$(j-c\Delta t/\Delta x - \frac{1}{2})^2 - \frac{1}{4} = (j-c\Delta t/\Delta x)(j-c\Delta t/\Delta x -1)<0$

und damit ein Widerspruch.

Ein etwas stärkerer Beweis stammt von Roe aus dem Jahr 1986, siehe Wesseling, Seite 346 für eine Skizze. Der Satz gilt ebenfalls für lineare Mehrschrittverfahren. Die Aussage ist also: Lineare Verfahren können nicht zweite Ordnung erreichen und um hohe Ordnung zu erreichen sind nichtlineare Verfahren unabdingbar. Hier haben sich insbesondere die TVD-Verfahren als nützlich erwiesen.

[Bearbeiten] TVD-Verfahren

Ein numerisches Verfahren hat die TVD-Eigenschaft (von Total Variation Diminishing), wenn für alle $n$

$TV(u^n)\leq TV(u^{n+1})$

gilt. Diese Eigenschaft gilt analog für die kontinuierliche Lösung und so ist es sinnvoll, dies auch vom numerischen Verfahren zu verlangen. Ist die totale Variation zum Anfangszeitpunkt beschränkt, so ist sie das bei einem TVD-Verfahren auch für alle späteren Zeitpunkte. Damit gilt nach den bereits gebrachten Aussagen zur TV-Stabilität: Konservative Verfahren mit Lipschitz-stetigem Fluss die die TVD-Eigenschaft haben sind TV-stabil und damit konvergent gegen eine schwache Lösung der Erhaltungsgleichung.

Die TVD-Eigenschaft stellt eine Bedingung an die kombinierte Raum-Zeit-Diskretisierung dar. Hierzu gibt es nützliche Sätze von Harten:

Satz

Gegeben sei das explizite Drei-Punkt-Verfahren

$u_j^{n+1} = u_j^n + a_j (u_{j+1}^n - u_j^n) + b_j(u_{j-1}^n-u_j^n),$

wobei $a j$ und $b j$ von $u n$ abhängen. Wenn die Koeffizienten für alle $j$ die Bedingung $a_j \geq 0$ , $b_j \geq 0$ und $a_j+b_{j+1} \leq 1$ erfüllen, ist das Verfahren TVD.

Beweis

Wir betrachten

$u_{j+1}^{n+1} - u_j^{n+1} = b_j (u_j^n - u_{j-1}^n) + (1-b_{j+1}-a_j)(u_{j+1}^n - u_j^n) + a_{j+1}(u_{j+2}^n- u_{j+1}^n).$

Nach Voraussetzung ist dann

$|u_{j+1}^{n+1} - u_j^{n+1}| \leq b_j |u_j^n - u_{j-1}^n| + (1-b_{j+1}-a_j) |u_{j+1}^n - u_j^n| + a_{j+1}|u_{j+2}^n- u_{j+1}^n|.$

Summation über $j$ ergibt genau die totale Variation zum Zeitpunkt $t n + 1$ :

$TV(u^{n+1}) = \sum_j|u_{j+1}^{n+1} - u_j^{n+1}| \leq \sum_j (b_{j+1}+1-b_{j+1}-a_j+a_j)|u_{j+1}^n - u_j^n| = TV(u^n)$

und damit die Aussage des Satzes.

Hierbei ist anzumerken, dass die Bedingung $a_j+b_{j+1} \leq 1$ eine Art CFL-Bedingung ist, damit gibt es also keinen Wiederspruch zur Konvergenz des Verfahrens. Für implizite Verfahren ergibt sich die Aussage

Satz

Gegeben sei das implizite Drei-Punkt-Verfahren

$u_j^{n+1} = u_j^n + a_j(u_{j+1}^n - u_j^n)+ b_j(u_{j-1}^n-u_j^n)+ c_j(u_{j-1}^{n+1}-u_j^{n+1})+ d_j(u_{j-1}^{n+1}-u_j^{n+1}).$

wobei $a j$ und $b j$ von $u n$ abhängen, während $c j$ und $d j$ von $u n + 1$ abhängen. Wenn die Koeffizienten die Bedingung $a_j\geq 0, b_j \geq 0, a_j+b_j \leq 1, d_j \leq 0$ und $-\infty <C\leq c_j$ erfüllen, ist das Verfahren TVD.

Der Beweis teilt das Verfahren in seinen expliziten und den impliziten Anteil auf und nutzt fuer den expliziten Anteil die vorherige Aussage. Die zusaetzlichen Bedingungen sind also dazu da, die totale Variation des impliziten Anteils ausreichend zu beschraenken.

Beide Saetze lassen sich auf breitere Diskretisierungssterne erweitern. Ferner kann der explizite Satz genutzt werden, um Bedingungen an TVD-Runge-Kutta-Verfahren zu ermitteln und solche zu konstruieren. Dazu spaeter mehr.

Wir wissen aus der vorherigen Vorlesung, dass lineare monotonieerhaltende Verfahren nicht von zweiter Ordnung sein können, und zwar global. Entsprechend sucht man den Ausweg über nichtlineare Verfahren. Leider stellt sich heraus, dass auch diese Grenzen haben: Osher bewies 1984, dass TVD-Verfahren an glatten Extrema zu erster Ordnung degenerieren und Goodman und LeVeque bewiesen 1985, dass ein TVD-Verfahren in mehreren Dimensionen maximal erste Ordnung haben kann. Nichtsdestotrotz sind TVD-Verfahren, die hoehere Ordnung ueber Limiter oder ueber ENO-Konstruktionen erzielen, sehr nuetzliche Hilfsmittel.

[Bearbeiten] Jameson-Schmidt-Turkel-Verfahren (JST-scheme), Teil 1

Neben den erwaehnten Limitern und ENO-Verfahren besteht eine andere Moeglichkeit darin, die instabile zentrale Diskretisierung als Ausgangspunkt zu nehmen und über nichtlineare Zusatzterme zu stabilisieren. Jameson, Schmidt und Turkel entwickelten hierzu das nach ihnen benannte JST-Verfahren, zum Einsatz in nichtlinearen Mehrgitterverfahren, da Limiter in diesem Kontext die Konvergenz verlangsamen können. Betrachten wir das semidiskrete Verfahren

$\partial_t u_j + \frac{1}{\Delta x}[F(u) -d]_{j-1/2}^{j+1/2}=0$

mit $F(u_{j+1/2})=\frac{1}{2}(f(u_j)+f(u_{j+1}))$ . Der Term künstlicher Dissipation wird definiert mittels

$d_{j+1/2} = \epsilon^{(2)}_{j+1/2}(u_{j+1}-u_j) - \epsilon_{j+1/2}^{(4)}(u_{j+2}-3u_{j+1}+3u_j-u_{j-1}).$

Der Koeffizient $\epsilon^{(2)}_{j+1/2}$ wird dabei so konstruiert, dass er in Regionen mit glatter Loesung von der Ordnung $Δ x 2$ ist, waehrend $\epsilon_{j+1/2}^{(4)}$ dort O(1) ist, so dass beide Terme dort von der Ordnung $Δ x 3$ sind, womit sich in glatten Regionen eine niedrige Diffusion ergibt. In der Naehe eines Schocks sollte $\epsilon^{(2)}_{j+1/2} \approx 1$ sein und damit das Verfahren erster Ordnung.

Das originale JST-Verfahren war nicht TVD, funktionierte aber gut in der Praxis, seitdem wurde es erheblich verbessert. Zunaechst ein Einschub.

[Bearbeiten] Positive Verfahren/LED-Verfahren

Wie bereits erwähnt sind TVD-Verfahren in mehreren Dimensionen maximal erster Ordnung, es ist also zu überlegen, ob man die TVD-Eigenschaft noch einmal relaxiert oder eine zumindest im eindimensionalen äquivalente Eigenschaft sucht. Dies führt auf positive Verfahren (Spekreijse 1987), beziehungsweise Local extremum diminishing schemes (LED, Jameson). Sei also das semidiskrete Verfahren

$\partial_t u_j + \sum_{k \neq j} c_{jk}(u_k-u_j)=0$

gegeben. Befindet sich im $u j$ ein lokales Maximum, so ist damit $u_k-u_j\leq 0$ . Dann gilt $\partial_t u_j\leq 0$ und damit, dass ein lokales Maximum nicht weiter wächst, wenn alle $c_{jk}\geq 0$ . Dies gilt offensichtlich auch im mehrdimensionalen, für unstrukturierte Gitter und für lokale Minima. Ein Verfahren mit der ersten Eigenschaft heisst LED, ein Verfahren mit der zweiten positiv. Manchmal redet man statt LED auch von einem diskreten Maximumsprinzip. In 1-D ist die LED-Eigenschaft äquivalent mit TVD (Tadmor 1988), nicht aber im mehrdimensionalen. Spekrijse nutzte positive Verfahren, um mehrdimensionale MUSCL-Verfahren zweiter Ordnung zu konstruieren, ein Konvergenzbeweis im mehrdimensionalen steht aber noch aus.

Für das JST-Verfahren gilt folgender Satz:

Satz

Es gelte immer dann, wenn $u j + 1$ oder $u j$ ein lokales Extremum ist, für die Koeffizienten des JST-Verfahrens:

$\epsilon_{j+1/2}^{(2)} \geq \frac{1}{4}|a_{j+1/2}|$ und $\epsilon_{j+1/2}^{(4)}=0,$

wobei der Koeffizient $a j + 1 / 2$ gegeben sei durch

$\frac{f(u_{j+1})-f(u_j)}{u_{j+1}-u_j}$ fuer $u_{j+1} \neq u_j$ und

$\partial_u f|u_j$ fuer

u j + 1 = u j

.

Dann ist das Verfahren LED.

Beweis

Sei $v j$ ein Maximum und damit die Koeffizienten $\epsilon_{j-1/2}^{(4)}$ und $\epsilon_{j+1/2}^{(4)}$ Null. Dann ergibt sich das semidiskrete Verfahren zu

$\partial_t u_j - (\epsilon_{j+1/2}^{(2)} - \frac{1}{2}a_{j+1/2})(u_{j+1}-u_j) + (\epsilon_{j-1/2}^{(2)}+ \frac{1}{2}a_{j-1/2})(u_j-u_{j-1})=0$

und damit, dass alle Koeffizienten positiv sind.

Damit wissen wir, wie die Koeffizienten $ε (2)$ und $ε (4)$ zu konstruieren sind.

[Bearbeiten] JST-Verfahren, Teil 2

Das JST-Verfahren wurde im Laufe der Jahrzehnte immer wieder angepasst, die aktuelle Version ist folgende, bei der die Koeffizienten definiert werden als:

$\epsilon_{j+1/2}^{(2)} = \alpha_{j+1/2}Q_{j+1/2}$ und $\epsilon_{j+1/2}^{(4)} = \frac{1}{2}\alpha_{j+1/2} (1-Q_{j+1/2})$

An einem Extremum sollte $Q$ also Eins sein. Dazu sei

Q j = R (u j + 1 - u j, u j - u j - 1)

und

Q j + 1 / 2 = max(Q j, Q j + 1)

mit

$R(x,y) = \left\{\begin{matrix} |\frac{x-y}{|x|+|y|}|^q & \mbox{falls } x\neq 0 \mbox{ oder } y\neq 0\\ 0 & \mbox{falls } x=y=0 \end{matrix}\right.$

$q$ ist dabei eine positive ganze Zahl. Wenn $x$ und $y$ entgegengesetzte Vorzeichen haben (entspricht einem Extremum), ist $R (x, y) = 1$ , ansonsten kleiner.

[Bearbeiten] Systeme von Gleichungen

Bisher haben wir ausschliesslich skalare Gleichungen betrachtet, manchmal in mehreren Dimensionen. Nun sind die meisten praxisrelevanten Gleichungen Systeme, die Frage ist also, welche Aussagen noch uebrig bleiben. Zunaechst gilt der Satz von Lax-Wendroff weiterhin. Schon schwieriger wird es bei der Entropiebedingung, diese muss man bei Systemen auf so genannte genuinely nonlinear systems einschraenken, was in etwa einer Konvexitaetsbedingung an die Flussfunktion entspricht. Damit ist noch keine Konvergenz des numerischen Verfahrens gezeigt, was der Situation im kontinuierlichen entspricht. Ausser fuer Spezialfaelle (starke Einschraenkungen Anfangsdaten oder das System) gibt es keine Aussagen zu Loesungen fuer nichtlineare Systeme.

Fuer Konvergenz des numerischen Verfahren brauchen wir Stabilitaet und wir wollen die Probleme an der totalen Variation erlaeutern. Fuer lineare Systeme ist es moeglich die totale Variation so zu definieren, dass das diagonalisierte System betrachtet wird und dann die totale Variation die Summe der totalen Variationen der einzelnen Loesungskomponenten. In jeder einzelnen Komponente ist die skalare Theorie anwendbar und damit auch die gesamte totale Variation nichtsteigend.

Fuer nichtlineare Systeme funktioniert das nicht mehr, da sich die Diagonalisierungsmatrizen mit der Loesung aendern. In der Tat kann die totale Variation schon fuer die kontinuerliche Loesung ansteigend sein. Ein Beispiel sind zwei kollidierende Schockwellen bei den Euler-Gleichungen, die nach der Kollision zwei auseinanderlaufende Schockwelle unter Erhoehung der totalen Variation erzeugen. Gerade bei den Euler-Gleichungen in mehreren Dimensionen tritt noch das Problem auf, dass ein Eigenwert mehrfach vorkommt, bei der Erweiterung der MHD-Gleichung ist dann sogar die schwache Entropieloesung des Riemann-Problems nicht mehr eindeutig. Einer der wenigen Konvergenzbeweise fuer nichtlineare Systeme ist Glimm's Random Choice Method.

[Bearbeiten] Zeitintegration

Wir betrachten ab jetzt die semidiskrete Gleichung

U t + F (U) = 0,

wobei $U$ ein Gittervektor ist und $F (u)$ aus der räumlichen Diskretisierung stammt. Es handelt sich also um ein nichtlineares System gewöhnlicher Differentialgleichungen, die mit einer der vielen Verfahren für Anfangswertprobleme gelöst werden können.

[Bearbeiten] TVD-Runge-Kutta-Verfahren

Eine Frage ist, unter welchen Bedingungen die Zeitdiskretisierung die TVD-Eigenschaft der räumlichen Diskretisierung erhält bzw. wann das Gesamtverfahren TVD ist. Wir hatten bereits Bedingungen an Verfahren bei Nutzung des expliziten Euler-Verfahrens hergeleitet. Diese Technik kann genutzt werden, um zu zeigen, dass das Heun-Verfahren die TVD-Eigenschaft erhält. Dieses ist ein Runge-Kutta-Verfahren 2. Ordnung der Form:

$U * = U n + Δ t F (U n),$
$U^{n+1} = U^n + \frac{\Delta t}{2}[F(U^n) + F(u^*)].$

Dies lässt sich als Folge zweier Schritte des expliziten Euler-Verfahrens schreiben:

$U * = U n + Δ t F (U n),$
$U * * = U * + Δ t F (U *),$
$U^{n+1} = \frac{1}{2}[U^n + U^{**}].$

Wenn nun das Verfahren, dass durch $F (U)$ TVD ist (etwa die Bedingungen des Satzes von Harten für das explizite Euler-Verfahren erfüllt), sind die einzelnen Schritte TVD und es gilt

$TV(U^{**}) \leq TV(U^*) \leq TV(U^n)$

und damit

$TV(U^{n+1}) \leq \frac{1}{2}TV(U^n) +\frac{1}{2}TV(U^{**})\leq TV(U^n),$

das Gesamtverfahren ist dann also TVD und hat dasselbe Stabilitätsgebiet. Dieses ist sogar bei dritter Ordnung möglich (Shu und Osher 1988):

$U * = U n + Δ t F (U n),$
$U^{**} = \frac{3}{4}U^n + \frac{1}{4}U^* + \frac{\Delta t}{4} F(U^*),$
$U^{n+1} = \frac{1}{3}U^n + \frac{2}{3}U^{**} + \frac{2\Delta t}{3}F(U^{**}).$

[Bearbeiten] Strong stability preserving methods

Allgemeiner heißt ein Verfahren strong stability preserving (SSP), wenn es genau wie das explizite Euler-Verfahren in der $L_{\infty}$ -Norm stabil ist, also

$\|U^{n+1}\|_{\infty} \leq \|U^n\|_{\infty}$

unter einer CFL-artigen Stabilitätsbedindung. Gottlieb, Shu und Tadmor bewiesen 2001, dass in der Klasse der SSP-Runge-Kutta-Verfahren $m$ -ter Ordnung mit $m$ Stufen keine im Vergleich zum expliziten Euler relaxierte Stabilitätsbedingung möglich ist. Das Verfahren von Heun und das genannte Verfahren 3. Ordnung sind in dieser Hinsicht also optimal.

[Bearbeiten] Stationäre Rechnungen, explizite Zeitintegration

Bei stationären Rechnungen ist es nicht von Belang, wie die Zwischenergebnisse aussehen sondern nur, ob die numerische Näherung an die stationäre Lösung sinnvoll ist, entsprechend kann dann auf die TVD-Eigenschaft der Zeitintegration verzichtet werden. Eine zweite Idee ist es, die konvektiven und die diffusiven Terme mit unterschiedlichen Verfahren zu behandeln. Zunächst betrachten wir explizite Zeitintegration. Die Idee ist es, die Koeffizienten der Runge-Kutta-Verfahren so zu wählen, dass das Stabilitätsgebiet maximal wird. Jameson schlägt folgendes Verfahren vor. Die rechte Seite wird aufgespalten in den konvektiven und den diffusiven Anteil

F (U) = Q (U) + D (U)

.

Dann wird das Runge-Kutta-Verfahren definiert mittels

u (0) = u n .

u (k) = u n - α k Δ t (Q (k - 1) + D (k - 1))

, für

k = 1, m

u n + 1 = u (m) .

Hierbei ist $Q (0) = Q (u n)$ , $D (0) = D (u n)$ und ferner $Q (k) = Q (u (k))$ sowie $D (k) = β k D (u (k)) + (1 - β k) D (u (k - 1)).$

Die Koeffizienten $α k$ für den konvektiven Anteil werden nun so gewählt, dass das Stabilitätsgebiet entlang der imaginären Achse vergrößert wird. Für ein 5-stufiges Schema schlägt Jameson folgendes vor:

$\alpha_1 = \frac{1}{4}$ , $\alpha_2 = \frac{1}{6}$ , $\alpha_3 = \frac{3}{8}$ , $\alpha_4= \frac{1}{2}$ ,

α 5 = 1

Die Koeffizienten $β k$ für den diffusiven Anteil werden so gewählt, dass das Stabilitätsgebiet entlang der negativen reellen Achse vergrößert wird. Für ein 5-stufiges Schema schlägt Jameson folgendes vor:

β 1 = 1

,

β 2 = 0

,

β 3 = 0.56

,

β 4 = 0

,

β 5 = 0.44

.

Im Vergleich zum klassischen Runge-Kutta-Verfahren ist das Stabilitätsgebiet auf der negativen reellen Achse etwa dreimal so groß.

[Bearbeiten] Implizite Verfahren

Als zweite Klasse wollen wir hier noch implizite Verfahren erwähnen, die dann zum Einsatz kommen können, wenn explizite Verfahren aufgrund ihres Stabilitätsgebiets einen zu kleinen Zeitschritt haben. In der numerischen Strömungsmechanik sind die betrachteten Gleichungen in der Regel sehr steif, interessant sind also solche Verfahren, ein möglichst großes Stabilitätsgebiet haben. Das ist für A-stabile Verfahren erfüllt und so sind BDF-2, die implizite Mittelpunktsregel und das Crank-Nicholson-Verfahren sehr beliebt. Diese Verfahren sind alle zweiter Ordnung, was nach der Zweite Dahlquist-Barriere auch das Optimum für A-stabile lineare Mehrschritt-Verfahren ist. Bei Verfahren höherer Ordnung gibt es im wesentlichen zwei Ansätze. Entweder werden BDF-Verfahren höherer Ordnung genommen, in der Hoffnung, dass die Exlusion der imaginären Achse nicht so schlimm ist. Sowohl BFD-3 als auch BDF-4 wurden erfolgreich eingesetzt.

Die zweite Möglichkeit ist, implizite Runge-Kutta-Verfahren zu nutzen. Diese können grundsätzlich A-stabil sein, allerdings ist es dann erforderlich, ein nichtlineares Gleichungssystem in jeder Stufe zu lösen. Bijl et al (2002) haben so genannte ESDIRK-Verfahren genutzt und gezeigt, dass diese Verfahren kompetitiv sind und besser als die vergleichbaren BDF-Verfahren sein können.

Gemein ist allen Verfahren, dass nichtlineare Gleichungssysteme gelöst werden muessen.

[Bearbeiten] Nichtlineare Gleichungssystemlöser

Zur Lösung von Systemen nichtlinearer Gleichungssysteme wurden zahlreiche Gleichungssystemlöser entwickelt. Wir werden hier die Löser betrachten, die zur numerischen Lösung der Euler- und Navier-Stokes-Gleichungen eingesetzt werden.

[Bearbeiten] Newton-Verfahren

Das Newton-Verfahren (auch Newton-Raphson-Verfahren, aber auch Simpson trug wesentlich zur Entwicklung bei) löst Nullstellengleichungen $F (x) = 0$ und beruht auf einer Linearisierung der Funktion $F (x)$ . Die Iterationsvorschrift ist gegeben durch

x k + 1 = x k - J (x k) - 1 F (x k),

wobei $J (x k)$ die Jacobi-Matrix von $F$ an der Stelle $x k$ ist. Da die Inverse in der Regel nicht zur Verfügung steht, wird diese Iteration in zwei Schritten durchgeführt, nämlich der Lösung eines linearen Gleichungssystems und dem Lösungsupdate:

Löse $J (x k)Δ x k = - F (x k)$
Berechne $x k + 1 = x k + Δ x k$

Das Newton-Verfahren konvergiert lokal quadratisch. Die konkrete Form des Jacobi-Matrix hängt insbesondere vom Zeitintegrationsverfahren ab und wie dieses formuliert ist. Betrachten wir das $Θ$ -Verfahren ( $\Theta \in [0,1]$ ):

u n + 1 = u n + Δ t ((1 - Θ) f (u n + 1) + Θ f (u n)).

Nach einer Finite-Volumen-Diskretisierung erhalten wir für das komplette System auf dem gesamten Rechengitter

Ω U n + 1 = Ω U n + Δ t ((1 - Θ) F (U n + 1) + Θ F (U n)),

wobei $Ω$ die Diagonalmatrix bezeichnet, bei der auf der Diagonalen die entsprechenden Volumina zu finden sind. Dies lässt sich auf verschiedene Weisen in eine Nullstellengleichung für den Vektor der Unbekannten $U = U n + 1$ umformen:

Ω U n + 1 + Ω U n - Δ t ((1 - Θ) F (U n + 1) + Θ F (U n)) = 0.

Die Struktur der Jacobi-Matrix ist (auch für andere Zeitintegrationsverfahren) immer die, dass sie aus der Summe der Ableitung des Vektors $U$ und der Ableitung der numerischen Flussfunktion besteht, wobei die Faktoren $Δ t$ und $Ω$ nur bei einem der Terme auftauchen. Es ergibt sich dann für das lineare Gleichungssystem in diesem Fall

$A \Delta u = [\Omega \partial U - \Delta t (1-\Theta )\partial F(U)]\Delta u.$

Hierbei hat man noch die Auswahl, in welchem Variablensatz man das Newton-Verfahren durchführen will (konservative oder primitive Variablen etwa). Bei konservativen Variablen ist $\partial U$ eine Diagonalmatrix, dafür sind die Ableitungen der numerischen Flussfunktion in primitiven Variablen häufig etwas einfacher.

Die obige Matrix ist aufgrund der Eigenschaften von $F (U)$ unsymmetrisch, indefinit, keine M-Matrix, für hinreichend große $Δ t$ schlecht konditioniert. Sie ist eine Blockmatrix, wobei die Größe der Blöcke der Anzahl der Variablen des Systems der Strömungsgleichungen entspricht (je nach Raumdimension 3-5).

Die Aussage der lokal quadratischen Konvergenz gilt für den Fall, wo die Ableitung exakt ausgerechnet wurde und das lineare Gleichungssystem in jedem Schritt exakt berechnet wird. Wird auf eins oder beides verzichtet, so heißt das Verfahren inexakt und verliert man die quadratische Konvergenz, dafür kann die benötigte Rechenzeit drastisch reduziert werden. Erstes Beispiel ist das vereinfachte Newton-Verfahren, bei dem die Ableitung nur im ersten Schritt exakte berechnet wird. Dieses Verfahren kann immerhin superlineare konvergieren. Bei inexakten Newton-Verfahren ist eine quantitative Konvergenzaussage schwierig, grob gilt, je schlechter approximiert wird, desto schlechter ist das Konvergenzverhalten. Wichtigster Faktor ist die Geschwindigkeit, mit der die linearen Gleichungssysteme auf befriedigende Genauigkeit gelöst werden können.

[Bearbeiten] Newton-Krylow-Verfahren

Die Variante, bei der das lineare Gleichungssystem mittels eine Krylow-Unterraum-Verfahrens gelöst werden, nennt man Newton-Krylow-Verfahren. Aufgrund der Eigenschaften der Systemmatrix kommen nicht alle Verfahren in Frage, sondern nur solche für allgemeine Systeme, insbesondere GMRES und BiCGSTAB. Nach praktischer Erfahrung ist BiCGSTAB hier die bessere Wahl (Meister, Vömel 2001).

[Bearbeiten] Matrixfreie Verfahren

Krylow-Unterraum-Verfahren haben die Eigenschaft, dass sie nie die Matrix $A$ explizit benötigen, sondern immer nur Matrix-Vektor-Produkte. Ist die Matrix, wie in diesem Fall eine Jacobi-Matrix, repräsentiert das Matrix-Vektor-Produkt eine Richtungsableitung und kann damit durch Funktionsauswertungen approximiert werden:

$\frac{\partial R(u^*)}{\partial u} x = \lim_{\epsilon\rightarrow 0} \frac{R(u^*+\epsilon x)-R(u^*)}{\epsilon} \approx \frac{R(u^*+\epsilon x)-R(u^*)}{\epsilon}$

für $ε > 0$ . Das sich ergebende Verfahren braucht somit deutlich weniger Speicher und kann schneller sein als ein Verfahren bei dem die Matrix-Vektor-Produkte ausgerechnet werden. Als grundlegendes Krylow-Unterraum-Verfahren ist hier GMRES die beste Wahl. Dies liegt daran, dass hier letztlich eine Störung der Matrix vorliegt und damit der aufgespannte Krylow-Unterraum ${z, A z, A 2 z,..., A n - 1 z}$ ebenfalls gestört wird. GMRES ist dabei aufgrund seiner Minimierungseigenschaft robuster als etwa BiCGSTAB gegen solche Störungen: In jedem Schritt wird die Norm $\|Ax-b\|_2$ über den aufgespannten Krylow-Unterraum minimiert, GMRES liefert also immer in jedem Schritt etwas sinnvolles.

In unserem konkreten Fall ergibt sich:

R (U) = Ω U + Ω U n - Δ t ((1 - Θ) F (U) + Θ F (U n))

und damit

$\frac{\partial R(u^*)}{\partial u} x \approx \Omega \epsilon x - \frac{\Delta t (1-\Theta )}{\epsilon}F(u^*)-F(u^*+\epsilon v).$

Pro Matrix-Vektor-Multiplikation muss also eine weitere Auswertung der rechten Seite erfolgen, da $F (u *)$ während eines Newton-Schritts konstant ist.

[Bearbeiten] Vorkonditionierung

Entscheidender als die Wahl des Krylow-Unterraum-Verfahrens ist die Wahl der Vorkonditionierung. Unterschieden wird zwischen Links- und Rechtsvorkonditionierung (und beidseitiger):

P - 1 A x = P - 1 b

entspricht Linksvorkonditionierung und

A P - 1 y = b

,

x = P - 1 y

Rechtsvorkonditionierung. Es geht darum, einen Vorkonditionierer zu finden, der

In der Anwendung wenig kostet.
Die inverse von A gut approximiert.

Jedes lineare iterative Verfahren ist als Vorkonditionierer geeignet, da eine Anwendung eines iterativen Verfahrens einer Approximation von $A - 1$ entspricht. Für Strömungsprobleme haben sich hier die unvollständige LU-Zerlegung (ILU), bei der die Matrix A selbst approximiert wird und das symmetrische Gauß-Seidel-Verfahren (SGS, $A^{-1}\approx (D+L)D^{-1}(D+R)$ ) etabliert. SGS ist im supersonischen Bereich extrem gut, da dann die Matrix fast tridiagonal wird und SGS damit ein exakter Löser wäre. Im Bereich kleiner Machzahlen ist dies nicht der Fall und die Performance sinkt. Ferner hat sich herausgestellt, dass Rechtsvorkonditionierung effektiver ist als Linksvorkonditionierung.

Rechtsvorkonditioniertes GMRES

GMRES basiert darauf, mit Hilfe der Arnoldi-Prozedur eine orthogonale Basis des Krylow-Unterraums zu berechnen und dann über das Minimierungsproblem die Lösung auszurechnen. GMRES rechnet im Krylow-Unterraum ${r 0, A P - 1 r 0,...,(A P - 1) m - 1 r 0}$ und wendet bei der Berechnung der Lösung einmalig den Vorkonditionier nochmal an:

Gegeben $x_0 \in \mathbb{R}^n$ , initialisiere $H \in \mathbb{R}^{(m+1)\times m}$ mit Nullen.

Fuehre dann den Arnoldi-Prozess durch:

Berechne $r 0 = b - A x 0$ . If $r 0 = 0$ , then END. $v_1 = \frac{r_0}{\|r_0\|_2}$ . $\beta=\|r_0\|_2$ .
For $j = 1,..., m$

Berechne $z j = P - 1 v j$ .

Berechne $w j = A z j$

For $i = 1,..., j$ do $h_{ij} = v_i^T w_j$

$w_j = w_j - \sum_{i=1}^j h_{ij} v_i$

Berechne $h_{j+1,j}=\|w_j\|_2$ und $v j + 1 = w / h j + 1, j .$

Berechne die Norm der Naeherungsloesung ueber $y_j = argmin_y \|\beta e_1 - H_jy\|$ . Ist diese kleine genug, STOP.
Definiere $V m = [v 1,.., v m]$ .

Berechne anschliessend die Naeherungsloesung:

x m = x 0 + P - 1 V m y m

mit

y m

wie oben. Wenn vorgesehen: RESTART, sonst ENDE.

Im letzten Schritt wird die Loesung als Linearkombination der vorkonditionierten Arnoldi-Basisvektoren $z i = M - 1 v i$ berechnet.

FGMRES

Eine weitere Variante sind Vorkonditionierer, die sich mit jedem Schritt ändern, dazu gehören insbesondere nichtlineare Mehrgitterverfahren. Dazu ist es nötig, das Krylow-Unterraum-Verfahren etwas anzupassen: Oben haben wir die Vektoren $z i$ nach ihrer Berechnung nicht weiter abspeichern muessen. Es reichte, die $v i$ und den Vorkonditionierer abzuspeichern. Dies aendert sich bei varieblen Vorkonditionierern, bei denen $z_i=P^{-1}_iv_i$ gilt. Eine Möglichkeit so vorzugehen bietet FGMRES (flexible GMRES) (Saad 1993).

Gegeben $x_0 \in \mathbb{R}^n$ , initialisiere $H \in \mathbb{R}^{(m+1)\times m}$ mit Nullen.

Fuehre dann den Arnoldi-Prozess durch:

Berechne $r 0 = b - A x 0$ . If $r 0 = 0$ , then END. $v_1 = \frac{r_0}{\|r_0\|_2}$ . $\beta=\|r_0\|_2$ .
For $j = 1,..., m$

Berechne $z_j = P^{-1}_jv_j$ .

Berechne $w j = A z j$

For $i = 1,..., j$ do $h_{ij} = v_i^T w_j$

$w_j = w_j - \sum_{i=1}^j h_{ij} v_i$

Berechne $h_{j+1,j}=\|w_j\|_2$ und $v j + 1 = w / h j + 1, j .$

Berechne die Norm der Näherungslösung ueber $y_j = argmin_y \|\beta e_1 - H_jy\|$ . Ist diese kleine genug, STOP.
Definiere $Z m = [z 1,.., z m]$ .

Berechne anschliessend die Näherungslösung:

x m = x 0 + Z m y m

mit

y m

wie oben. Wenn vorgesehen: RESTART, sonst ENDE.

Der Zusatzaufwand ist hier das Abspeichern der Vektoren $z i$ (die $v i$ ) müssen eh gespeichert werden). Dafür ist es nun möglich, flexible Vorkondonditionierer einzusetzen. Eine andere Möglichkeit ist GMRES-R von van der Vorst und Vuik, 1994.

[Bearbeiten] Mehrgitter-Verfahren

Die wesentliche Alternative zu Newton-Verfahren sind Mehrgitter-Verfahren. Die schnellsten Verfahren zur Berechnung der stationären Euler-Gleichungen stammen aus dieser Klasse. Die Grundidee ist, eine Hierarchie von gröber werdenden Gittern aufzubauen und diese zu nutzen, um eine Approximation des unbekannten Fehlers einer Anfangsnäherung der Lösung zu berechnen.

[Bearbeiten] Lineare Mehrgitter-Verfahren

Zunächst sei auf dem feinsten Gitter ein lineares Gleichungssystem $A l x = b l$ mit regulärer Matrix $A_l \in \mathbb{R}^{n \times n}$ gegeben. Eine Iteration des Mehrgitter-Verfahrens sieht dann folgendermaßen aus:

M G (x l, b l, l)

if

(l = 0)

, $x_l = A_l^{-1}b_l$ (Löse exakt auf gröbstem Gitter)

else

$x_l = S_l^{\nu_1}(x_l,b_l)$ (Vorglättung)

r l - 1 = R l - 1, l (b l - A l x l)

(Restriktion)

v l - 1 = 0

Für $(j=0; \; j<\gamma; \; j++)$

M G (v l - 1, r l - 1, l - 1)

(Berechnung der Grobgitterkorrektur)

x l = x l + P l, l - 1 v l - 1

(Prolongation der Grobgitterkorrektur)

$x_l = S_l^{\nu_2}(x_l,b_l)$ (Nachglättung)

end if

End

Was ein sinnvoller Glätter $S l$ ist, hängt von der zu lösenden partiellen Differentialgleichung ab. Für viele sind Gauß-Seidel- oder Jacobi-Relaxation geeignet. Da niederfrequente Fehleranteile auf feinen Gittern hochfrequenten Fehleranteilen auf gröberen Gittern entsprechen, ergibt sich mit der Residuumsgleichung

A l - 1 e = A l - 1 (x l - 1 - x *) = A l - 1 x l - 1 - b l - 1 = r l - 1

ein Problem mit ähnlicher Struktur wie das Ursprungsproblem, allerdings mit kleinerer Dimension. Dies wird rekursiv bis zum gröbsten Gitter wiederholt ( $γ = 1$ entspricht einem V-Zyklus, $γ = 2$ einem W-Zyklus), wo das Gleichungssystem direkt gelöst werden kann.

V-Zyklus

W-Zyklus

Der berechnete Fehler wird dann sukzessive mittels Prolongation P auf die feineren Gitter rückprolongiert und die ursprüngliche Näherung der Lösung korrigiert. Dies funktioniert bei einem linearen Problem $A x = b$ mit Lösung $x *$ , da dann der Fehler $e = x k - x *$ der Näherungslösung $x k$ über die Residuumsgleichung $A e = r = A x k - b$ berechnet werden kann. Es folgt noch Nachglättung.

Für viele elliptische und parabolische Gleichungen gibt es umfassende Konvergenztheorie zu linearen Mehrgitterverfahren. Dies ist bei den Euler- und Navier-Stokes-Gleichungen nicht der Fall, sie finden aber manchmal als Vorkonditionierer oder als Löser in Newton-Verfahren Verwendung.

[Bearbeiten] Full Approximation Scheme

Auf ein nichtlineares Problem $L (u) = f$ lässt sich das obige Vorgehen nicht direkt übertragen, da die Residuumsgleichung $L (e) = r$ gar nicht lösbar sein muss. Deswegen löst man da auf dem groben Gitter statt dessen $L (u + e) = L (u) + r$ , was im linearen Fall äquivalent wäre. Es ergibt sich dann

$MG(\tilde{u}_l,u_l, f_l, l)$

if

(l = 0)

, $u_l = L_l^{-1}f_l$

else

$u_l = S_l^{\nu_1}(u_l,f_l)$

r l = f l - L l u l

Wähle $\tilde{u}_{l-1}$ und

s l - 1

$f_{l-1} = L_{l-1}(\tilde{u}_{l-1}) + s_{l-1}R_{l-1,l}r_l$

Für $(j=0; \; j<\gamma; \; j++)$ $MG(\tilde{u}_{l-1},u_{l-1},f_{l-1},l-1)$

$u_l = u_l + (1/s_{l-1})P_{l-1,l}(u_{l-1}-\tilde{u}_{l-1})$

$u_l = S_l^{\nu_2}(u_l,f_l)$

end if

End

$\tilde{u}_l$ beschreibt dabei eine Approximation an die Lösung und $s l$ wird so klein gewählt, dass das entsprechende nichtlineare Gleichungssystem lösbar ist. $s = 1$ entspricht dem Full Approximation Scheme (FAS) von Achi Brandt (1977). In Hackbusch (81, 82) und seinem Buch Multigrid Methods and Applications (85) findet sich eine ausführlichere Beschreibung von Auswahlmöglichkeiten.