Theorie und Numerik von Erhaltungsgleichungen

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

1. Energieerhaltung
2. Impulserhaltung
3. Drehimpulserhaltung
4. 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“. Ebenso ergeben sich bei zerfallenden Atomkernen nach dem Zerfall kleinere Massensummen der Tochterteilchen.

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

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 $\mathrm{\Omega }$. Die darin enthaltene Masse $M$ ist dann gegeben durch das Integral der Dichte $\rho$ über $\mathrm{\Omega }$:

$M=\underset{\mathrm{\Omega }}{\int }\rho d\mathrm{\Omega }.$

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

$M_t=\frac{d}{dt}\underset{\mathrm{\Omega }}{\int }\rho d\mathrm{\Omega }=-\underset{\partial \mathrm{\Omega }}{\int }f\cdot ndS,$

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

$M_t=\frac{d}{dt}\underset{\mathrm{\Omega }}{\int }\rho d\mathrm{\Omega }=-\underset{\partial \mathrm{\Omega }}{\int }\rho v\cdot ndS.$

Zu beachten ist hierbei, dass dabei außer der Integrierbarkeit von $\rho$ und $\rho v$ keine Voraussetzungen an die Funktionen $\rho$ und $v$ oder an das Volumen $\mathrm{\Omega }$ 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}\underset{\mathrm{\Omega }}{\int }\rho d\mathrm{\Omega }+\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }\cdot \rho vd\mathrm{\Omega }=0.$

Ist das Volumen $\mathrm{\Omega }$ nicht zeitabhängig, lässt sich Zeitdifferentiation und Integration vertauschen. Da die Aussage für beliebige Volumina gilt, ergibt sich die differentielle Form der Massenerhaltung

${\rho }_t+\mathrm{\nabla }\cdot \rho v=0.$

Allgemeiner haben Erhaltungsgleichungen die Form

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

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

$u_t+\mathrm{\nabla }\cdot f(u,\mathrm{\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 $\varphi \in C_0^1(ℝ\times ℝ)$, also aus dem Raum der stetig differenzierbaren Funktionen mit kompaktem Träger. Integration liefert

$\underset{0}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}[\varphi u_t+\varphi f(u{)}_x]dxdt=0.$

Partielle Integration ergibt aufgrund des kompakten Trägers von $\varphi$

$\underset{0}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}[{\varphi }_{t}u+{\varphi }_{x}f(u)]dxdt=-\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\varphi (x,0)u(x,0)dx.$

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

Beispiele für Erhaltungsgleichungen

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

${\partial }_t\rho +\mathrm{\nabla }\cdot \mathbf{\text{<mi fromhbox="1">m</mi>}}=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},i=1,2,3$

${\partial }_t\rho E+\mathrm{\nabla }\cdot (H\mathbf{\text{<mi fromhbox="1">m</mi>}})=\sum _{j=1}^3{\partial }_{x_j}(\sum _{i=1}^3{S}_{ij}{v}_i-W_j)+q.$

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

${\partial }_t\rho +\mathrm{\nabla }\cdot \mathbf{\text{<mi fromhbox="1">m</mi>}}=0,$

${\partial }_t{m}_i+\sum _{j=1}^3{\partial }_{x_j}(m_i{v}_j+p{\delta }_{ij})=0,i=1,2,3$

${\partial }_t\rho E+\mathrm{\nabla }\cdot (H\mathbf{\text{<mi fromhbox="1">m</mi>}})=0.$

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, das 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(ℝ)$:

$u_t+f(u{)}_x=0$ auf dem Gebiet $ℝ\times {ℝ}^{+}$ und

$u(x,0)=u_0(x)$ auf $ℝ.$

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

$\underset{x_1}{\overset{x_2}{\int }}\eta (u(x,t_2))dx\le \underset{x_1}{\overset{x_2}{\int }}\eta (u(x,t_1))dx+\underset{t_1}{\overset{t_2}{\int }}\psi (u(x_1,t))dt-\underset{t_1}{\overset{t_2}{\int }}\psi (u(x_2,t))dt.$

$\eta (u)$ ist die sogenannte Entropiefunktion und $\psi (u)$ der Entropiefluss. An dieses Entropie-Entropiefluss-Paar stellt man die Bedingung, dass $\eta$ glatt und konvex ist, $\psi$ glatt und ferner die Beziehung

${\psi }^{\prime }(u)={\eta }^{\prime }(u)f^{\prime }(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:

$\underset{0}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}[{\varphi }_t\eta (u)+{\varphi }_x\psi (u)]dxdt+\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\varphi (x,0)\eta (u(x,0))dx\ge 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!

Die Lösung der verschwindenden Viskosität

Eine alternative Herangehensweise ist die Betrachtung der Gleichung

$u_t+f(u{)}_x=ϵu_{xx}$

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

Erste Aussagen zur Numerik

Dazu sei ein Gitter aus Gitterzellen ${\mathrm{\Omega }}_i$ gegeben mit ${\mathrm{\Omega }}_i=(x_{i-1},x_i)$ und diskrete Zeitpunkte $t_0,t_1,\dots$. Die Schrittweiten beschreiben wir mit $\mathrm{\Delta }{x}_i=x_i-x_{i-1}$ und $\mathrm{\Delta }{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 ${\mathrm{\Omega }}_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:

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{\mathrm{\Delta }{t}_n}{\mathrm{\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\mathrm{\Delta }{x}_i{u}_i^{n+1}=\sum _{i=I}^J\mathrm{\Delta }{x}_i{u}_i^n-\mathrm{\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.

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:

1. $F(u_1,\dots ,u_j)=f(u)$
2. $F$ ist Lipschitz-stetig in allen Komponenten.

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

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 $\mathrm{\Delta }{x}_l$ und $\mathrm{\Delta }{t}_l$ mit $\mathrm{\Delta }{x}_l,\mathrm{\Delta }{t}_l\to 0$ für $l\to \mathrm{\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\le t<T,l=1,2,\dots$

und das Supremum über alle Unterteilungen $(\xi {)}_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 $\mathrm{\Omega }=[a,b]\times [0,T]$ gelte

$\underset{0}{\overset{T}{\int }}\underset{a}{\overset{b}{\int }}|U_l(x,t)-u(x,t)|dxdt\to 0$, für $l\to \mathrm{\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:

$\underset{0}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}[{\varphi }_{t}u+{\varphi }_{x}f(u)]dxdt=-\underset{\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\varphi (x,0)u(x,0)dx$ für alle $\varphi \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 $\mathrm{\Delta }x$, $\mathrm{\Delta }t$ gegeben durch die konservative Vorschrift

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

Dies multiplizieren wir mit einer beliebigen Testfunktion $\varphi$ und erhalten mit der Bezeichnung ${\varphi }_i^n$ für die Einschränkung von $\varphi$ auf die Zelle $i$ und das Zeitintervall $[t_n,t_{n+1}]$:

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

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

$\sum _{n=0}^{\mathrm{\infty }}\sum _{i=-\mathrm{\infty }}^{\mathrm{\infty }}{\varphi }_i^n(u_i^{n+1}-u_i^n)=-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}\sum _{n=0}^{\mathrm{\infty }}\sum _{i=-\mathrm{\infty }}^{\mathrm{\infty }}{\varphi }_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|\ge i_0$ und für $n\ge 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_m{b}_m-a_1{b}_0-\sum _{i=1}^{m-1}(a_{i+1}-a_i)b_i$. Damit ergibt sich

$-\sum _{i=-\mathrm{\infty }}^{\mathrm{\infty }}{\varphi }_i^0{u}_i^0-\sum _{n=1}^{\mathrm{\infty }}\sum _{i=-\mathrm{\infty }}^{\mathrm{\infty }}({\varphi }_i^n-{\varphi }_i^{n-1})u_i^n=\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}\sum _{n=0}^{\mathrm{\infty }}\sum _{i=-\mathrm{\infty }}^{\mathrm{\infty }}({\varphi }_{i+1}^n-{\varphi }_i^n)F_{i-1/2}^n$

und nach weiteren Umformungen

$\mathrm{\Delta }t\mathrm{\Delta }x[\sum _{n=1}^{\mathrm{\infty }}\sum _{i=-\mathrm{\infty }}^{\mathrm{\infty }}\left(\frac{{\varphi }_i^n-{\varphi }_i^{n-1}}{\mathrm{\Delta }t}\right)u_i^n+\sum _{n=0}^{\mathrm{\infty }}\sum _{i=-\mathrm{\infty }}^{\mathrm{\infty }}\left(\frac{{\varphi }_{i+1}^n-{\varphi }_i^n}{\mathrm{\Delta }x}\right)F_{i-1/2}^n]=-\mathrm{\Delta }x\sum _{i=-\mathrm{\infty }}^{\mathrm{\infty }}{\varphi }_i^0{u}_i^0.$

Diese Gleichung gilt auf jedem Gitter $l$. Wir betrachten nun den Grenzwert $l\to \mathrm{\infty }$. Definieren wir Anfangsdaten beispielsweise derart, dass $u_i^0=1/\mathrm{\Delta }x\underset{x_i-1/2}{\overset{x_i+1/2}{\int }}u(x,0)dx$, so konvergiert die rechte Seite aufgrund der Glattheit von $\varphi$ und der Konvergenz von $U_l$ gegen $u$ gegen $-\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\varphi (x,0)u(x,0)dx$. Der erste Term in der Summe auf der linken Seite konvergiert ebenfalls wegen der Glattheit von $\varphi$ und der Konvergenz von $U_l$ gegen $\underset{0}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\varphi }_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 $\underset{0}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\varphi }_x(x,t)f(u(x,t))dxdt$ und damit der Beweis des Satzes.

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 $\eta$, $\psi$ und dazu ein diskretes Analogon $\mathrm{\Phi }$, das zu $\psi$ konsistent ist (im bekannten Sinne). Dann ist die diskrete Entropiebedingung:

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

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.

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 ${\tau }_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 $\mathrm{\Delta }x$ und $\mathrm{\Delta }t$ für $t\to \mathrm{\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+\mathrm{\Delta }t$ mit Fehler $E^{n+1}$. Also haben wir

$\begin{array}{cc}{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))]-\mathrm{\Delta }t{\tau }^n.\end{array}$

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 (Äquivalenzsatz von Lax).

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\ne 0$.

mit periodischen Randbedingungen und die auf ganz $ℝ$ 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\mathrm{\Delta }x$ in eine Fourier-Reihe

$E_i^n=\sum _{-N}^N{c}_j^n{e}^{ij\mathrm{\Delta }xI}$

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 $\Vert e_i{\Vert }_{L_2}=\Vert c_j{\Vert }_{L_2}$. Die $L_2$-Stabilitätsbedingung reduziert sich dann darauf, dass das numerischer Verfahren genau dann stabil ist, wenn der Spektralradius der Amplifikationsmatrix $\rho (G)$ betragsmäßig kleiner gleich eins ist.

Beispiel

Der einfachste Fall ist die lineare Advektionsgleichung

$u_t+a{u}_x=0,$ mit $a\in ℝ$

In der Zeit nehmen wir das explizite Euler-Verfahren $u^{n+1}=u^n+\mathrm{\Delta }t(a{u}_x{)}^n$ gekoppelt mit zentralen Differenzen auf einem äquidistanten Gitter im Raum. Der zweite Term wird also mittels

$a{u}_x\approx \frac{a}{2\mathrm{\Delta }x}(u_{i+1}^n-u_{i-1}^n)$

approximiert. Insgesamt ergibt sich

$\frac{u_i^{n+1}-u_i^n}{\mathrm{\Delta }t}+\frac{a}{2\mathrm{\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 $e^{ij\mathrm{\Delta }xI}$ mit $\varphi =j\mathrm{\Delta }x$:

$c_j^{n+1}-c_j^n+\frac{a\mathrm{\Delta }t}{2\mathrm{\Delta }x}{c}_j^n(e^{I\varphi }-e^{-I\varphi })=0.$

Die Amplifikationsmatrix ist nun gegeben durch

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

Das Verfahren ist $L_2$-stabil, falls $|G|\le 1$ für alle $\varphi$, was hier nicht der Fall ist, da $|G{|}^2=1+{\left(\frac{a\mathrm{\Delta }t}{\mathrm{\Delta }x}\right)}^2{\sin }^2\varphi .$ 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

$a{u}_x\approx \frac{a}{\mathrm{\Delta }x}(u_i^n-u_{i-1}^n)$

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

$\frac{a\mathrm{\Delta }t}{\mathrm{\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)\le 1+c\mathrm{\Delta }t$.

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

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

$\Vert N(P)-N(Q)\Vert \le N\Vert P-Q\Vert ,$

so ergibt sich immerhin

$\Vert E^{n+1}\Vert \le \Vert N(u(x,t^n)+E^n)-N(u(x,t^n))\Vert +\mathrm{\Delta }t\Vert {\tau }^n\Vert \le \Vert E^n\Vert +\mathrm{\Delta }t\Vert {\tau }^n\Vert ,$

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)\le R$ und $\mathrm{supp}(v(\cdot ,t))\subset [-M,M]\mathrm{\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 $\mathrm{\Delta }t<\mathrm{\Delta }{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 $\mathrm{\Delta }{t}_0$ und ein $R>0$, so dass für alle n und alle $\mathrm{\Delta }t$ mit $\mathrm{\Delta }t<\mathrm{\Delta }{t}_0$ und $n\mathrm{\Delta }t\le T$ gelte:

$TV(u^n)\le 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 $\mathrm{\Delta }t$ generiert. Dann konvergiert die Lösung für $\mathrm{\Delta }t\to 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 $(\mathrm{\Delta }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(\mathrm{\Delta }{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.

Monotone Verfahren

Eine Eigenschaft der hier betrachteten Erhaltungsgleichungen ist, dass der Lösungsoperator invers monoton ist: gilt für zwei Anfangsdaten $u^0\ge v^0$, dann gilt für die entsprechenden Lösungen $u(t,\cdot )\ge v(t,\cdot )$. Ein konservatives Verfahren, das die entsprechende diskrete Eigenschaft hat, nennt man monoton. Genauer heisst das Verfahren monoton, falls für die numerische Lösung aus

$u_j^n\ge v_j^n$ für alle $j$ folgt, dass

$u_j^{n+1}\ge v_j^{n+1}$ für alle $j$

gilt. Dies ist äquivalent mit

$\frac{\partial H_j}{\partial u_{j+l}}\ge 0$ für 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 Entropielösung konvergieren. 1980 bewiesen dann Crandall und Majda folgendes Resultat:

Satz

Die durch ein monotones Verfahren generierte Lösung konvergiert in $L_1$ gegen die schwache Entropielösung der Erhaltungsgleichung für $\mathrm{\Delta }t$ gegen Null mit $\mathrm{\Delta }t/\mathrm{\Delta }x$ fix. Dieses Resultat gilt auch im mehrdimensionalen.

Kröner und Ohlberger bewiesen 2000 eine grundlegende Fehlerabschätzung für monotone Verfahren. Um 2000 wurde von verschiedenen Autoren folgendes mehrdimensionale Resultat erarbeitet, was es auf der einen Seite tatsächlich erlaubt, a priori Fehlerschätzer zu erarbeiten, aber auf der anderen Seite das große Problem monotoner Verfahren vor Augen führt.

Satz

Sei $u_0(x)$ in $BV({ℝ}^d)$ und $u(x,t)$ eine schwache Entropielösung der Erhaltungsgleichung. Sei $U$ eine durch ein monotones numerisches Verfahren mit dem expliziten Eulerverfahren als Zeitintegration generierte numerische Lösung, mit einem in bestimmter Weise beschränkten Zeitschritt. Sei $K$ eine bestimmte Teilmenge von ${ℝ}^d\times {ℝ}^{+}$. Dann existiert eine Konstante $C>0$, so dass

$\Vert u-U{\Vert }_{L_1(K)}\le C{h}^{1/4}.$

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

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

1. Die $l_1$-Norm der numerischen Lösung ist mit der Zeit nichtsteigend ($l_1$-kontraktiv).
2. Damit gilt auch: Die Totale Variation ist nicht steigend (TVD).
3. Und daraus folgt: Monotonie der Lösung 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 abzuschwächen, um eine höhere Konvergenzordnung zu erzielen.

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\ge 0$ für alle $m$.

Beweis

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

Hinrichtung: Sei das Verfahren monotonieerhaltend, aber es existiere ein p mit ${\gamma }_p<0.$ Wir wählen nun Anfangsdaten:

$u_j^0=0$, $j\le 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\ge 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)\ge 0.$

also ist das Verfahren monotonieerhaltend.

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|\mathrm{\Delta }t/\mathrm{\Delta }x\in \mathbb{N}.$

Beweis

Sei $u(0,x)={(\frac{x}{\mathrm{\Delta }x}-\frac{1}{2})}^2-\frac{1}{4}$ und die entsprechenden numerischen Anfangsdaten $u_j^0={(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}{\mathrm{\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\mathrm{\Delta }t/\mathrm{\Delta }x-\frac{1}{2})}^2-\frac{1}{4}.$

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

${(j-c\mathrm{\Delta }t/\mathrm{\Delta }x-\frac{1}{2})}^2-\frac{1}{4}\ge 0$ für alle $j$.

Für $c\mathrm{\Delta }t/\mathrm{\Delta }x\in \mathbb{N}$ und größer Null wählen wir nun $j$ mit $j>c\mathrm{\Delta }t/\mathrm{\Delta }x>j-1$. Dann gilt

${(j-c\mathrm{\Delta }t/\mathrm{\Delta }x-\frac{1}{2})}^2-\frac{1}{4}=(j-c\mathrm{\Delta }t/\mathrm{\Delta }x)(j-c\mathrm{\Delta }t/\mathrm{\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.

TVD-Verfahren

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

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

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\ge 0$, $b_j\ge 0$ und $a_j+b_{j+1}\le 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}|\le 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}|\le \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}\le 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\ge 0,b_j\ge 0,a_j+b_j\le 1,d_j\le 0$ und $-\mathrm{\infty }<C\le c_j$ erfüllen, ist das Verfahren TVD.

Der Beweis teilt das Verfahren in seinen expliziten und den impliziten Anteil auf und nutzt für den expliziten Anteil die vorherige Aussage. Die zusätzlichen Bedingungen sind also dazu da, die totale Variation des impliziten Anteils ausreichend zu beschränken.

Beide Sätze 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 später mehr.

Wir wissen aufgrund des Satzes von Godunow, 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 höhere Ordnung über Limiter oder über ENO-Konstruktionen erzielen, sehr nützliche Hilfsmittel.

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

Neben den erwähnten Limitern und ENO-Verfahren besteht eine andere Möglichkeit 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}{\mathrm{\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}={ϵ}_{j+1/2}^{(2)}(u_{j+1}-u_j)-{ϵ}_{j+1/2}^{(4)}(u_{j+2}-3u_{j+1}+3u_j-u_{j-1}).$

Der Koeffizient ${ϵ}_{j+1/2}^{(2)}$ wird dabei so konstruiert, dass er in Regionen mit glatter Lösung von der Ordnung $\mathrm{\Delta }{x}^2$ ist, während ${ϵ}_{j+1/2}^{(4)}$ dort O(1) ist, so dass beide Terme dort von der Ordnung $\mathrm{\Delta }{x}^3$ sind, womit sich in glatten Regionen eine niedrige Diffusion ergibt. In der Nähe eines Schocks sollte ${ϵ}_{j+1/2}^{(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. Zunächst ein Einschub.

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\ne j}{c}_{jk}(u_k-u_j)=0$

gegeben. Befindet sich im $u_j$ ein lokales Maximum, so ist damit $u_k-u_j\le 0$. Dann gilt ${\partial }_t{u}_j\le 0$ und damit, dass ein lokales Maximum nicht weiter wächst, wenn alle $c_{jk}\ge 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 folgt aus der LED-Eigenschaft die TVD-Eigenschaft (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:

${ϵ}_{j+1/2}^{(2)}\ge \frac{1}{4}|a_{j+1/2}|$ und ${ϵ}_{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}$ für $u_{j+1}\ne u_j$ und

${\partial }_{u}f|u_j$ für $u_{j+1}=u_j$.

Dann ist das Verfahren LED.

Beweis

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

${\partial }_t{u}_j-({ϵ}_{j+1/2}^{(2)}-\frac{1}{2}{a}_{j+1/2})(u_{j+1}-u_j)+({ϵ}_{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.

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:

${ϵ}_{j+1/2}^{(2)}={\alpha }_{j+1/2}{Q}_{j+1/2}$ und ${ϵ}_{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)=\{\begin{array}{cc}|\frac{x-y}{|x|+|y|}{|}^q& \text{falls }x\ne 0\text{ oder }y\ne 0\\ 0& \text{falls }x=y=0\end{array}$

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

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 übrig bleiben. Zunächst 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 einschränken, was in etwa einer Konvexitätsbedingung an die Flussfunktion entspricht. Damit ist noch keine Konvergenz des numerischen Verfahrens gezeigt, was der Situation im kontinuierlichen entspricht. Ausser für Spezialfälle (starke Einschränkungen an die Anfangsdaten oder das System) gibt es keine Aussagen zu Lösungen für nichtlineare Systeme.

Für Konvergenz des numerischen Verfahren brauchen wir Stabilität und wir wollen die Probleme an der totalen Variation erläutern. Für lineare Systeme ist es möglich die totale Variation so zu definieren, dass das diagonalisierte System betrachtet wird und dann die totale Variation die Summe der totalen Variationen der einzelnen Lösungskomponenten. In jeder einzelnen Komponente ist die skalare Theorie anwendbar und damit auch die gesamte totale Variation nichtsteigend.

Für nichtlineare Systeme funktioniert das nicht mehr, da sich die Diagonalisierungsmatrizen mit der Lösung ändern. In der Tat kann die totale Variation schon für die kontinuierliche Lösung ansteigend sein. Ein Beispiel sind zwei kollidierende Schockwellen bei den Euler-Gleichungen, die nach der Kollision zwei auseinanderlaufende Schockwellen unter Erhöhung 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 Entropielösung des Riemann-Problems nicht mehr eindeutig. Einer der wenigen Konvergenzbeweise für nichtlineare Systeme ist Glimm's Random Choice Method.

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.

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:

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

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

1. $U^{*}=U^n+\mathrm{\Delta }tF(U^n),$
2. $U{**}^{=}{U}^{*}+\mathrm{\Delta }tF(U^{*}),$
3. $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^{**})\le TV(U^{*})\le TV(U^n)$

und damit

$TV(U^{n+1})\le \frac{1}{2}TV(U^n)+\frac{1}{2}TV(U^{**})\le 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):

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

Strong stability preserving methods

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

$\Vert U^{n+1}{\Vert }_{\mathrm{\infty }}\le \Vert U^n{\Vert }_{\mathrm{\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.

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-{\alpha }_k\mathrm{\Delta }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)}={\beta }_{k}D(u^{(k)})+(1-{\beta }_k)D(u^{(k-1)}).$

Die Koeffizienten ${\alpha }_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},{\alpha }_5=1$

Die Koeffizienten ${\beta }_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:

${\beta }_1=1,{\beta }_2=0,{\beta }_3=0,56,{\beta }_4=0,{\beta }_5=0,44$.

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

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 müssen.

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.

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:

1. Löse $J(x_k)\mathrm{\Delta }{x}_k=-F(x_k)$
2. Berechne $x_{k+1}=x_k+\mathrm{\Delta }{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 $\mathrm{\Theta }$-Verfahren ($\mathrm{\Theta }\in [0,1]$):

$u^{n+1}=u^n+\mathrm{\Delta }t((1-\mathrm{\Theta })f(u^{n+1})+\mathrm{\Theta }f(u^n)).$

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

$\mathrm{\Omega }{U}^{n+1}=\mathrm{\Omega }{U}^n+\mathrm{\Delta }t((1-\mathrm{\Theta })F(U^{n+1})+\mathrm{\Theta }F(U^n)),$

wobei $\mathrm{\Omega }$ 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:

$\mathrm{\Omega }{U}^{n+1}+\mathrm{\Omega }{U}^n-\mathrm{\Delta }t((1-\mathrm{\Theta })F(U^{n+1})+\mathrm{\Theta }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 $\mathrm{\Delta }t$ und $\mathrm{\Omega }$ nur bei einem der Terme auftauchen. Es ergibt sich dann für das lineare Gleichungssystem in diesem Fall

$A\mathrm{\Delta }u=[\mathrm{\Omega }\partial U-\mathrm{\Delta }t(1-\mathrm{\Theta })\partial F(U)]\mathrm{\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 $\mathrm{\Delta }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 exakt 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.

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).

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=\underset{ϵ\to 0}{\lim }\frac{R(u^{*}+ϵx)-R(u^{*})}{ϵ}\approx \frac{R(u^{*}+ϵx)-R(u^{*})}{ϵ}$

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,Az,A^{2}z,\dots ,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 $\Vert Ax-b{\Vert }_2$ über den aufgespannten Krylow-Unterraum minimiert, GMRES liefert also immer in jedem Schritt etwas sinnvolles.

In unserem konkreten Fall ergibt sich:

$R(U)=\mathrm{\Omega }U+\mathrm{\Omega }{U}^n-\mathrm{\Delta }t((1-\mathrm{\Theta })F(U)+\mathrm{\Theta }F(U^n))$

und damit

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

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

Vorkonditionierung

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

$P^{-1}Ax=P^{-1}b$

entspricht Linksvorkonditionierung und

$A{P}^{-1}y=b$, $x=P^{-1}y$

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

1. In der Anwendung wenig kostet.
2. 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,\dots ,(A{P}^{-1}{)}^{m-1}{r}_0\}$ und wendet bei der Berechnung der Lösung einmalig den Vorkonditionier nochmal an:

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

Führe dann den Arnoldi-Prozess durch:

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

   Berechne $z_j=P^{-1}{v}_j$.

   Berechne $w_j=A{z}_j$

   For $i=1,\dots ,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}=\Vert w_j{\Vert }_2$ und $v_{j+1}=w/h_{j+1,j}.$

   Berechne die Norm der Näherungslösung über $y_j={\mathrm{argmin}}_y\Vert \beta e_1-H_{j}y\Vert$. Ist diese kleine genug, STOP.

3. Definiere $V_m=[v_1,\dots ,v_m]$.

Berechne anschliessend die Näherungslösung:

$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 Lösung 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 müssen. Es reichte, die $v_i$ und den Vorkonditionierer abzuspeichern. Dies ändert sich bei varieblen Vorkonditionierern, bei denen $z_i=P_i^{-1}{v}_i$ gilt. Eine Möglichkeit so vorzugehen bietet FGMRES (flexible GMRES) (Saad 1993).

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

Führe dann den Arnoldi-Prozess durch:

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

   Berechne $z_j=P_j^{-1}{v}_j$.

   Berechne $w_j=A{z}_j$

   For $i=1,\dots ,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}=\Vert w_j{\Vert }_2$ und $v_{j+1}=w/h_{j+1,j}.$

   Berechne die Norm der Näherungslösung über $y_j={\mathrm{argmin}}_y\Vert \beta e_1-H_{j}y\Vert$. Ist diese kleine genug, STOP.

3. Definiere $Z_m=[z_1,\dots ,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.

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.

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 {ℝ}^{n\times n}$ gegeben. Eine Iteration des Mehrgitter-Verfahrens sieht dann folgendermaßen aus:

$MG(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++)$ $MG(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 ($\gamma =1$ entspricht einem V-Zyklus, $\gamma =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 $Ax=b$ mit Lösung $x^{*}$, da dann der Fehler $e=x^k-x^{*}$ der Näherungslösung $x^k$ über die Residuumsgleichung $Ae=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.

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.

Nichtlineares Mehrgitterverfahren von Jameson & al

Eine weitere Variante wurde von Jameson und Turkel im Laufe von Jahrzehnten für die stationären Euler-Gleichungen entwickelt. Dabei wird auf Nachglättung verzichtet (${\nu }_2=0$) und als Vorglätter wird ein Schritt mit dem oben erwähnten speziellen Runge-Kutta-Verfahrens gemacht, das nämlich neben einem grossen Stabilitätsgebiet noch Glättungseigenschaften hat. Da die Restriktion auf ein gröberes Gitter führt, ist es dort möglich, den Zeitschritt deutlich zu erhöhen und so sehr rasche Konvergenz zum stationären Zustand zu erzielen. Dazu wird das Runge-Kutta-Verfahren umformuliert als

$U_l^{(1)}=U_l^{(0)}-{\alpha }_1\mathrm{\Delta }{t}_l[F_l^{(0)}+C_l]$

$\dots$

$U_l^{(q+1)}=U_l^{(0)}-{\alpha }_{q+1}\mathrm{\Delta }{t}_l[F_l^{(q)}+C_l].$

Dabei ist $C_l=R_{l+1,l}{F}_{l+1}(U_{l+1})-F_l(U_l^{(0)})$.

Als Zyklus wird ein 4- oder 5-Gitter-W-Zyklus gewählt. Die Konvergenz wird dabei durch weitere Maßnahmen noch beschleunigt. Zunächst indem nicht in jeder Zelle mit demselben Zeitschritt gerechnet wird, sondern mit einer global gültigen CFL-Zahl, mittels derer dann der lokale Zeitschritt ausgerechnet wird (local time stepping). Allgemeiner haben Lee, Turkel, Roe und Van Leer nichtlineare Vorkonditionierer entwickelt, mit denen die Zeitableitung multipliziert wird.

Als zweites wird Residuenglättung genutzt (Jameson, Baker 1984), was den Glättungseffekt erhöht. Dazu werden in jeder Stufe des Runge-Kutta-Verfahrens die Flüsse $F(U)$ vor der Aufdatierung von $U^{(k)}$ durch $\overline{F}(U)$ ersetzt, wobei $\overline{F}(U)$ durch die folgende Gleichung mit $ϵ$ klein gegeben ist:

$-ϵ{\overline{F}}_{j-1}^{(k)}+(1+2ϵ){\overline{F}}_j^{(k)}-ϵ{\overline{F}}_{j+1}^{(k)}=F(U{)}_j^{(k)}.$

Diese Gleichung lässt sich leicht lösen und erhält den stationären Zustand, ist aber für die Berechnung instationärer Probleme komplett ungeeignet.

Es ist möglich, mit diesem Verfahren stationäre Lösungen der Euler-Gleichungen für Tragflächenumströmungen mit drei bis fünf W-Zyklen zu berechnen (Jameson, Caughey 2001). Die stationären Navier-Stokes-Gleichungen sind schwieriger zu lösen, je nach Problem sind 30 bis 100 Schritte notwendig.

Dual Timestepping

Das Problem bei dem oben beschriebenen Mehrgitter-Verfahren ist, dass ein Großteil der Beschleunigungstechniken nur für stationäre Gleichungen funktioniert. Um die Technik auf instationäre Probleme zu übertragen wird das Dual Time Stepping Verfahren genutzt (Jameson 91). Gegeben sei dazu die Semidiskrete instationäre Gleichung

$U_t+F(U)=0.$

Wird diese mittels eines impliziten Verfahrens integriert, ergibt sich ein nichtlineares Gleichungssystem

$F^{*}(U)=0,$

wobei $F^{*}$ von der konkreten Wahl der Zeitintegration abhängt. Wird dazu eine zusätzliche Zeitableitung addiert in Dualzeit, ergibt sich das modifizierte Problem

$\frac{\partial U}{\partial t^{*}}+F^{*}(U)=0.$

Die Lösung des eigentlichen Problems ist dann die Lösung des letzten Problems wenn die duale Zeitableitung verschwindet, also der stationäre Zustand in Dualzeit erreicht ist. Dafür wird das oben erwähnte Mehrgitterverfahren angewandt. Aufgrund der zusätzlichen Terme wird dabei das Spektrum des Operators $F^{*}$ etwas verschoben, was die Konvergenz negativ beeinflusst. Im Fall der RANS-Gleichungen sind 50 bis 250 Iterationen notwendig.

Konvergenztheorie im hyperbolischen

Eine Konvergenztheorie im hyperbolischen Fall ist für keines der Verfahren bekannt. Ein wesentliches Problem sind Schocks. Durch die Kopplung auf dem groben Gittern werden Informationen von beiden Seiten des Schocks theoretisch nach Prolongation über das komplette Gitter verteilt. Dies erschwert die theoretische Analyse bisher entscheidend.

Fluid-Struktur-Interaktion

Als letztes betrachten wir ein Thema, bei dem hyperbolische Erhaltungsgleichungen nur eine Nebenrolle spielen, nämlich gekoppelte Systeme von partiellen Differentialgleichungen, bei denen eine die Strömung eines Fluids beschreibt und eine zweite eine Struktur. Das meistuntersuchte Problem ist dabei die Kombination, bei der Fluid und Struktur über Kräfte interagieren, handelt es sich um Luft, spricht man von Aeroelastizität. Das Fluid verformt also durch mechanische Kräfte das Fluid, was wiederum Kräfte in der Struktur bewirkt. Diese können sich in nichtlinearer Weise aufschaukeln. Beispiele sind

• Flattern von Flugzeugtragflächen oder
• der Zusammensturz der Tacoma-Narrows-Brücke.

Eine andere Form ist die thermomechanische Kopplung, bei der an der Struktur aufgrund von Hitze plastische Verformungen zu befürchten sind oder hervorgerufen werden sollen. Beispiele sind

• der Wiedereintritt von Raumfahrzeugen in die Atmosphäre, bei denen die erhitzte Atmosphäre das Gefährt erhitzt,
• thermische Belastungen in Motoren und Triebwerken oder
• Abkühlungsprozesse mittels Gasen bei Stahlumformung.

Modellierung

Gegeben sind zwei partielle Differentialgleichungen auf verschiedenen Gebieten ${\mathrm{\Omega }}_1$ und ${\mathrm{\Omega }}_2$. Dazu seien Kopplungsbedingungen gegeben für die Menge ${\mathrm{\Omega }}_1\cap {\mathrm{\Omega }}_2$. Konkreter haben wir ein System für das Fluid

$\frac{d}{dt}{x}_1=f(x_1,x_2)$

und eines für die Struktur

$\frac{d}{dt}{x}_2=g(x_1,x_2).$

Die Variablensätze 1 und 2 sind dabei den verschiedenen Gebieten zugeordnet und die Kopplungsbedingung ist in den einzelnen Abbildungen enthalten. Zum Beispiel könnte f den Navier-Stokes- oder Euler-Gleichungen entsprechen und g der Wärmeleitungsgleichung in einer Struktur. Kopplungsbedingung wäre dann, dass auf dem gemeinsamen Rand die Temperatur T und der Wärmefluss q übereinstimmen. Dies ist ein erstes nützliches Modellproblem vor der echten thermomechanischen Kopplung.

Kopplungsverfahren

Grundsätzlich wird zwischen zwei Lösungsansätzen unterschieden: monolithischen, bei denen das gekoppelte System als ein Gleichungssystem gelöst wird, und partitionierten, bei denen bestehende Lösungsverfahren zur Lösung der Einzelsysteme eingesetzt werden. Wesentlicher Vorteil der letzteren Ansätze ist, dass bestehende Löser für die jeweiligen Probleme verwendet werden können, was den Entwicklungsaufwand drastisch reduziert. Für gekoppelte Probleme existiert dagegen in der Regel nichts, so dass die Implementierung eines monolithischen Verfahrens sehr hohen Aufwand bedeuten kann. Wir werden uns hier auf partitionierte Verfahren beschränken.

Die Einzelsysteme seien also mit dem jeweiligen Verfahren im Raum diskretisiert, wobei wir der Einfachheit halber die Notation für $x_1$ und $x_2$ beibehalten:

$\frac{d}{dt}{x}_1=F(x_1,x_2),$

$\frac{d}{dt}{x}_2=G(x_1,x_2).$

Lose Kopplung

Einfachster Fall ist nun die lose oder schwache Kopplung. Gegeben sei für die jeweiligen Systeme ein Zeitintegrationsverfahren ${\mathrm{\Phi }}_i$, was aus der Lösung $x_i^n$ zum Zeitpunkt $t_n$ die Lösung $x_i^{n+1}$ zum Zeitpunkt $t_{n+1}$ generiert. Die erste Möglichkeit ist dann folgender Algorithmus, der parallel ausgeführt werden kann:

$x_1^{n+1}={\mathrm{\Phi }}_1(x_1^n,x_2^n),$

$x_2^{n+1}={\mathrm{\Phi }}_2(x_1^n,x_2^n).$

Eine üblichere Variante, die nur sequentiell funktioniert, ist:

$x_1^{n+1}={\mathrm{\Phi }}_1(x_1^n,x_2^n),$

$x_2^{n+1}={\mathrm{\Phi }}_2(x_1^{n+1},x_2^n).$

Beide Kopplungsverfahren sind maximal erster Ordnung in der Zeit und explixit, so dass sich erhebliche Zeitschritteinschränkungen ergeben können. Beide Varianten kann man innerhalb eines Zeitschritts iterieren, um die Genauigkeit zu erhöhen.

Einschub Thermische Kopplung

Wir betrachten nun konkret die Kopplung der Navier-Stokes-Gleichungen mit der Wärmeleitungsgleichung. Im diskreten Fall ist es nicht möglich, in Struktur und Fluid sowohl Temperatur als auch Wärmefluss am Rand vorzuschreiben, da dies für keine der Gleichungen sinnvolle Randbedingungen ergibt. Es gibt also die Möglichkeit, dem Fluid und der Struktur Temperatur oder Wärmefluss vorzuschreiben.

Die eine Möglichkeit bedeutet: Die Wärmeleitungsgleichung hat als Randbedingung eine Temperaturvorgabe und gibt dem Fluid einen Wärmefluss, während die Strömung als Randbedingung einen Wärmefluss erhält und eine Randtemperatur an die Struktur zurückgibt. Laut Giles 1997 ist dieses Verfahren instabil.

Deswegen erhält die Strömung am Rand Temperaturdaten und die Struktur einen Wärmefluss.

Starke Kopplung

Schwache Kopplung hat verschiedene Probleme: Geringe Ordnung, möglicherweise starke Zeitschrittbeschränkungen und schließlich erfüllt das Verfahren keine zusätzlichen Bedingungen wie Energieerhaltung. Eine Lösung ist eine implizite Herangehensweise, wobei dann die Frage ist, wie ein partitioniertes Lösungsverfahren aussieht. Folgende Formulierung bietet sich an (siehe beispielsweise Mathies, Steindorf 2002):

$x_1^{n+1}={\mathrm{\Phi }}_1(x_1^{n+1},x_2^{n+1}),$

$x_2^{n+1}={\mathrm{\Phi }}_2(x_1^{n+1},x_2^{n+1}).$

Ein Newton-Schritt für dieses System wäre dann:

$\left(\begin{array}{cc}I-D_{x_1}{\mathrm{\Phi }}_1& -D_{x_2}{\mathrm{\Phi }}_1\\ -D_{x_1}{\mathrm{\Phi }}_2& I-D_{x_2}{\mathrm{\Phi }}_2\end{array}\right)\left(\begin{array}{c}\mathrm{\Delta }{x}_1\\ \mathrm{\Delta }{x}_2\end{array}\right)=-\left(\begin{array}{c}{x}_1-{\mathrm{\Phi }}_1(x_1,x_2)\\ x_2-{\mathrm{\Phi }}_2(x_1,x_2)\end{array}\right).$

Nun liegen die entsprechenden partiellen Ableitungen nicht vor. Diese können allerdings entweder durch eine matrixfreie Technik in einem Krylov-Unterraum-Verfahren wie oben bereitgestellt werden oder durch Newton-Schritte im lokalen Verfahren.

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, das 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.

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.