Reading Time: 9 minutes

Die Wahl einer Zeitintegrationsmethode ist eine der wichtigsten Entscheidungen in einer wissenschaftlichen Simulation. Das Verfahren bestimmt, wie sich die numerische Lösung von einer Zeitstufe zur nächsten bewegt, wie klein der Zeitschritt sein muss, wie viel jeder Schritt kostet und ob schnelle physikalische Prozesse aufgelöst oder unterdrückt werden.

Die gemeinsame Unterscheidung zwischen expliziten und impliziten Methoden ist nützlich, bietet jedoch keine vollständige Auswahlregel. Eine implizite Methode ist nicht automatisch genauer, und eine explizite Methode ist nicht automatisch für ernsthafte Simulationen geeignet. Stabilität, Genauigkeit, Rechenkosten, Steifigkeit, numerische Dämpfung und die physikalischen Zeitskalen des Problems müssen zusammen betrachtet werden.

Das zentrale Prinzip ist einfach: Die numerische Stabilität sagt uns nur, ob Fehler kontrolliert bleiben. Es sagt uns nicht, ob die berechnete Lösung der wahren physikalischen Lösung nahe kommt.

Was ist Zeitintegration?

Nachdem eine PDE im Raum diskretisiert wurde, wird sie oft zu einem System gewöhnlicher Differentialgleichungen:

du/dt = F(u, t)

Ein Zeitintegrator nähert sich der Veränderung des Vektors u über einen endlichen Schritt:

tⁿ → tⁿ⁺¹ = tⁿ + Δt

Die genaue Entwicklung ist normalerweise nicht verfügbar, daher erstellt der Algorithmus eine Annäherung aus bekannten Werten, abgeleiteten Bewertungen oder einem System, das den unbekannten zukünftigen Zustand betrifft.

Eine explizite Methode berechnet den neuen Zustand direkt aus bereits verfügbaren Informationen. Eine implizite Methode definiert den neuen Zustand durch eine Gleichung, die gelöst werden muss.

Explizite Zeitintegration

Forward Euler ist die einfachste explizite Methode:

uⁿ⁺¹ = uⁿ + Δt F(uⁿ, tⁿ)

Alles auf der rechten Seite ist bekannt. Es ist kein lineares oder nichtlineares System erforderlich. Dies macht jeden Schritt kostengünstig und einfach zu parallelisieren.

Explizite Runge-Kutta-Methoden höherer Ordnung berechnen mehrere Zwischenstufen. Die klassische Methode vierter Ordnung verwendet vier Derivateauswertungen:

k₁ = F(uⁿ, tⁿ)

k₂ = F(
    uⁿ + 0.5 Δt k₁,
    tⁿ + 0.5 Δt
)

k₃ = F(
    uⁿ + 0.5 Δt k₂,
    tⁿ + 0.5 Δt
)

k₄ = F(
    uⁿ + Δt k₃,
    tⁿ + Δt
)

uⁿ⁺¹ = uⁿ
      + Δt(k₁ + 2k₂ + 2k₃ + k₄) / 6

Explizite Methoden sind attraktiv, wenn jede Derivatbewertung erschwinglich ist und die Stabilitätsgrenze keine übermäßige Anzahl von Schritten erzwingt.

Implizite Zeitintegration

Rückwärts Euler bewertet das Derivat im unbekannten zukünftigen Zustand:

uⁿ⁺¹ = uⁿ + Δt F(uⁿ⁺¹, tⁿ⁺¹)

Der neue Wert erscheint auf beiden Seiten. Ein lineares Problem kann eine Matrixlösung erfordern, während ein nichtlineares Problem Newton-Iterationen oder einen anderen nichtlinearen Algorithmus erfordert.

Crank-Nicocolson mittelt die Ableitung zwischen dem aktuellen und dem zukünftigen Zustand:

uⁿ⁺¹ = uⁿ
      + 0.5 Δt [
          F(uⁿ, tⁿ)
          + F(uⁿ⁺¹, tⁿ⁺¹)
        ]

Implizite Methoden kosten pro Schritt mehr, aber geeignete Schemata können für Zeitschritte stabil bleiben, die eine explizite Methode divergieren würden. Dies ist besonders wertvoll für steife Systeme und feine räumliche Gitter.

Stabilität ist keine Genauigkeit

Eine stabile Berechnung liefert nicht unbedingt eine genaue physikalische Flugbahn. Die Unterscheidung kann mit der linearen Testgleichung untersucht werden:

dy/dt = λy

Die genaue Lösung nach einem Malschritt ist:

y(t + Δt) = exp(λΔt)y(t)

Eine numerische Methode erzeugt stattdessen:

yⁿ⁺¹ = R(z)yⁿ

z = λΔt

Die Funktion R(z) ist der Verstärkungsfaktor. Absolute Stabilität erfordert:

|R(z)| ≤ 1

Diese Bedingung verhindert ein unbegrenztes numerisches Wachstum für ein abklingendes Testproblem. Es garantiert nicht, dass R(z) exp(z) sehr nahe kommt.

Ein implizites Verfahren kann mit einem sehr großen Zeitschritt begrenzt bleiben, während die Zerfallsrate, Phase oder das vorübergehende Ansprechen schlecht wiedergegeben werden. Der NASA-Bericht explizite, implizite und hybride Methoden erörtert die Notwendigkeit, die Genauigkeit zu berücksichtigen, anstatt die Stabilität allein zu verwenden, um eine Methode zu rechtfertigen.

Stabilitätsbereiche verstehen

Der Stabilitätsbereich ist der Wertsatz von z = λΔt, für den der Amplifikationsfaktor begrenzt bleibt.

vorwärts Euler

Vorwärts hat Euler:

R(z) = 1 + z

Sein Stabilitätsbereich erfüllt:

|1 + z| ≤ 1

Dies bildet eine Scheibe mit dem Radius Eins bei −1 zentriert. Entlang der negativen realen Achse ist das stabile Intervall:

−2 ≤ z ≤ 0

Klassische RK4

Die klassische Runge-Kutta-Methode vierter Ordnung hat einen größeren, aber immer noch begrenzten Stabilitätsbereich. Entlang der negativen realen Achse bleibt es ungefähr stabil bis:

z ≈ −2.785

Dies ist erheblich größer als das Vorwärts-Euler-Intervall, aber keine explizite Runge-Kutta-Methode kann die gesamte linke Hälfte der komplexen Ebene umfassen.

Eine praktische Einführung in die absolute Stabilität, A-Stabilität und L-Stabilität finden Sie im Crash-Kurs über numerische Oden .

Rückwärts Euler

Rückwärts Euler hat:

R(z) = 1 / (1 - z)

Sein Stabilitätsbereich enthält die komplette linke Halbebene. Es ist daher a-stabil.

Wenn z zunehmend negativ wird, nähert sich der Verstärkungsfaktor Null. Stark abfallende Modi werden schnell unterdrückt. Dies macht rückwärts Euler L-stabil, obwohl es nur erste Ordnung genau ist.

Crank-Niconson

Crank-Nicocolson hat:

R(z) = (1 + z/2) / (1 - z/2)

Es ist auch a-stabil, da sein Stabilitätsbereich die linke Halbebene umfasst. Als z → −∞:

R(z) → −1

Hochsteife Modi werden nicht auf Null gefahren. Stattdessen können sie sich im Vorzeichen abwechseln, während sie nahezu konstante Größe beibehalten. Crank-Niconson ist daher nicht L-stabil und kann nichtphysische zeitliche Schwingungen erzeugen, wenn sehr große Schritte auf steife Systeme angewendet werden.

A-Stabilität und L-Stabilität

Eine a-stabile Methode ist für jeden Test-Gleichwert mit einem nicht-positiven Realteil stabil, unabhängig von der Zeitschrittgröße.

Eine L-stabile Methode ist a-stabil und erfüllt auch:

R(z) → 0 as z → −∞

Diese Unterscheidung ist wichtig für steife Systeme. A-Stabilität verhindert explosives Wachstum, während L-Stabilität sorgt dafür, dass ungelöste schnelle Abklingmodi stark gedämpft werden.

Nicht jede implizite Methode ist a-stabil, und nicht jede a-stabile Methode ist L-stabil. Die Eigenschaften gehören eher zum einzelnen Schema als zur gesamten impliziten Kategorie.

Reihenfolge der Genauigkeit

Die Reihenfolge eines Verfahrens bestimmt, wie schnell sein Fehler abnimmt, wenn der Zeitschritt kleiner wird.

Für eine Bestellmethode p:

Local truncation error = O(Δt^(p+1))
Global error           = O(Δt^p)
Verfahren Typ Bestellen Lokaler Fehler Globaler Fehler
vorwärts Euler Explizit 1 O(Δt²) O(Δt)
Rückwärts Euler Implizit 1 O(Δt²) O(Δt)
Crank-Niconson Implizit 2 O(Δt³) O(Δt²)
BDF2 Impliziter Multischritt 2 O(Δt³) O(Δt²)
Klassische RK4 Explizit 4 O(Δt⁵) O(Δt⁴)
Dormand–Prinz 5(4) Explizite eingebettete RK 5 mit Schätzer vierter Ordnung methodenabhängig ungefähr O(Δt⁵) für die Lösung fünfter Ordnung

Vorwärts und rückwärts haben Euler die gleiche formale Ordnung, obwohl sich ihre Stabilitätseigenschaften stark unterscheiden. RK4 kann bei gleicher Schrittweite viel genauer sein als rückwärts Euler, wenn die Stabilität den Einsatz ermöglicht.

Das explizite oder implizite Label beschreibt in erster Linie, wie ein Schritt berechnet wird. Es bestimmt nicht die formale Reihenfolge.

Was ist Steifheit?

Ein System ist steif, wenn es stark getrennte Zeitskalen enthält, und explizite Stabilitätsanforderungen erzwingen Schritte, die viel kleiner sind als die, die zur Auflösung des interessierenden Verhaltens erforderlich sind.

Bedenken Sie:

dy/dt = -1000(y - cos(t)) - sin(t)

Die gewünschte Lösung kann auf einer Zeitskala der Reihenfolge eins variieren, aber eine schnell abfallende Komponente hat eine Zeitskala nahe 0.001. Eine explizite Methode muss möglicherweise den schnellen Modus für Stabilität auflösen, auch wenn dieser Modus physisch unwichtig geworden ist.

Eine geeignete implizite Methode kann den schnellen Zerfall überschreiten und der langsameren Lösung folgen. Dies ist der Hauptgrund, warum implizite Methoden für steife Reaktionssysteme, Diffusionsgleichungen, elektrische Schaltkreise und eng gekoppelte Multiphysikmodelle verwendet werden.

Numerische Dämpfung

Implizitheit bedeutet nicht automatisch eine starke Dämpfung. Die Dämpfung wird durch den Verstärkungsfaktor der Methode gesteuert.

Rückwärts Euler unterdrückt die Modi stark, wenn |λΔt| groß ist. Dies kann erwünscht sein, wenn diese Moden eine ungelöste Steifheit darstellen. Es kann unerwünscht sein, wenn sie Wellen oder Transienten darstellen, die gemessen werden sollten.

Crank-Niconson führt viel weniger hochfrequente Dämpfung ein. Dies bewahrt ein gewisses Oszillationsverhalten, kann aber auch unerwünschte steife numerische Schwingungen erhalten.

Die Flow-3D-Diskussion von Implizite und explizite numerische Methoden veranschaulichen, wie große implizite Schritte das vorübergehende Verhalten verzerren können. Der Effekt sollte nicht als universeller fester Unterentspannungsfaktor interpretiert werden. Seine Größe hängt vom Integrationsschema, dem Zeitschritt, der Gleichung und dem iterativen Solver ab.

Unterentspannung ist ein separates Thema

Unterrelaxation wird häufig in nichtlinearen oder gekoppelten iterativen Lösern verwendet:

u(updated) =
    u(old)
    + α [
        u(computed)
        - u(old)
      ]

Der Parameter α liegt normalerweise zwischen Null und Eins. Kleinere Werte können eine iterative Lösung stabilisieren, verlangsamen jedoch ihre Konvergenz und verändern den scheinbaren Übergang, wenn die Iterationen vor der vollständigen Konvergenz gestoppt werden.

Unterentspannung ist nicht unvermeidbare Eigenschaft jedes impliziten Zeitintegrators. Es handelt sich um eine zusätzliche algorithmische Wahl, die innerhalb des nichtlinearen Lösungsprozesses auftreten kann.

CFL-Einschränkungen für explizite Methoden

Für eine Advektionsgleichung folgen explizite Methoden üblicherweise einer Courant-Bedingung:

Δt ≤ C Δx / |v|

Die Konstante C hängt von der räumlichen Methode und dem Zeitintegrator ab.

Bei einer explizit integrierten Diffusionsgleichung skaliert die Grenze typischerweise wie folgt:

Δt ≤ C Δx² / D

Diese quadratische Abhängigkeit kann auf feinen Maschen teuer werden. Die Halbierung der Zellengröße kann für ein diffusionskontrolliertes explizites Schema ungefähr viermal so viele Zeitschritte erfordern.

Diese Einschränkungen bedeuten nicht, dass explizite Methoden ungenau sind. Sie definieren einen Stabilitätsbereich. Bei hyperbolischen Simulationen erfordert die physikalische Notwendigkeit, die Wellenfahrt aufzulösen, bereits einen Schritt ähnlich dem CFL-Grenzwert.

Kosten pro Schritt

Explizite Methoden erfordern normalerweise Funktionsbewertungen, Flussberechnungen oder spärliche Matrix-Vektor-Produkte. Ihre Schritte sind relativ kostengünstig und skalieren auf paralleler Hardware oft gut.

Implizite Methoden können Folgendes erfordern:

  • Matrixbaugruppe
  • Jacobi-Konstruktion
  • Lineare Systemlösung
  • Preconditioner-Setup
  • Nichtlineare Newton-Iterationen
  • Konvergenzprüfungen

Ein implizites Verfahren ist nur dann effizient, wenn der größere, nutzbare Schritt die zusätzlichen Kosten jeder Lösung kompensiert.

Der Vergleich sollte daher die Gesamtkosten auf einem festen Fehlerniveau und nicht allein der Schrittzählung verwenden.

Adaptive explizite Integration

Eingebettete Runge-Kutta-Paare schätzen den Fehler, ohne zwei völlig unabhängige Integrationen abzuschließen.

Dormand-Prinz 5 (4), oft als RK45 bezeichnet, teilt eine Reihe von Zwischenstufen, um sowohl eine Approximation fünfter Ordnung als auch eine Fehlerschätzung niedrigerer Ordnung zu konstruieren.

Der normalisierte Fehler kann wie folgt bewertet werden:

error_ratio =
    estimated_error
    / (
        absolute_tolerance
        + relative_tolerance
          * solution_scale
      )

Wenn das Verhältnis unter eins liegt, kann der Schritt akzeptiert werden. Wenn es eins überschreitet, wird der Schritt abgelehnt und mit einem kleineren Δt wiederholt.

Ein typisches Update hat die Form:

Δt(new) =
    safety
    * Δt(old)
    * error_ratio^(-1/(p+1))

Praktische Implementierungen begrenzen auch, wie schnell der Schritt wachsen oder schrumpfen kann.

Adaptive implizite Integration

Implizite Löser können Fehler durch eingebettete Formeln, BDF-Methoden mit variabler Ordnung, Fehlerschätzungen oder Schrittverdoppelung schätzen.

Schrittverdoppelung vergleicht:

  • Ein Schritt der Länge Δt
  • Zwei Schritte der Länge Δt/2

Der Unterschied schätzt den zeitlichen Fehler. Dies kann mehrere implizite Lösungen erfordern, obwohl Matrixfaktorisierungen oder Vorkonditionierer manchmal wiederverwendet werden können, wenn der Operator ähnlich bleibt.

Große Produktionsmodelle können mehrere unabhängige Einschränkungen kombinieren. Die Pism-Zeitschritt-Dokumentation veranschaulicht, wie CFL, Diffusivität, Ausgabe und modellspezifische Grenzwerte interagieren in einem realen Simulationscode.

IMEX-Methoden

Implizite-explizite Methoden teilen die rechte Seite in steife und nicht steife Komponenten auf:

du/dt = Fexplicit(u) + Fimplicit(u)

Der kostengünstige nicht-steife Begriff wird explizit bewertet, während der steife Begriff implizit behandelt wird.

Für ein Konvektions-Diffusionsproblem:

∂u/∂t
+ v · ∇u
= D∇²u

Der Advektionsbegriff kann explizit und der Diffusionsbegriff implizit sein. Dies vermeidet eine globale nichtlineare Lösung für die vollständige Gleichung, während die strenge explizite Diffusionsbeschränkung entfernt wird.

IMEX-Schemata erfordern kompatible explizite und implizite Formeln. Ihre Reihenfolge und Stabilität hängen von der vollständigen gepaarten Methode ab, nicht nur von jeder Komponente isoliert.

Expliziter, impliziter und IMEX-Vergleich

Eigentum Explizit Implizit IMEX
Schrittberechnung direkt aus bekannten Staaten erfordert die Lösung für den zukünftigen Zustand Kombiniert direkte und implizite Stufen
Kosten pro Schritt Normalerweise niedrig normalerweise höher zwischen explizit und vollständig implizit
Stabilitätsbereich Für explizite RK-Methoden begrenzt kann sehr groß sein; methodenabhängig Abhängig von beiden Komponenten
Steife Systeme oft ineffizient normalerweise angemessen Geeignet, wenn die Steifheit getrennt werden kann
Wellenprobleme Oft effizient und verlustarm Erfordert sorgfältige Auswahl der Dämpfungs- und Phaseneigenschaften Nützlich für gemischte Wellen und steife Begriffe
Umsetzung Relativ einfach Benötigt lineare oder nichtlineare Löser Erfordert Operatortrennung und gepaarte Formeln

Auswahl einer Methode nach Physik

Problemtyp Gemeinsamer Ausgangspunkt Die Vernunft
nicht steife Ode Adaptive Explizite RK-Methode Niedrige Schrittkosten und zuverlässige eingebettete Fehlerkontrolle
Wellenausbreitung Explizite RK- oder Strukturerhaltungsmethode Die physikalische Auflösung setzt oft schon einen kleinen Schritt voraus
Explizite Diffusion auf einem feinen Netz implizite oder imex-Methode Vermeidet die schwere Einschränkung Δx²
stark steifes Reaktionssystem BDF, Radau oder ein anderer steifer Solver Explizite Stabilitätsanforderungen können unpraktisch sein
Advektions-Diffusionssystem IMEX oder Operator-Splitting Unterschiedliche Begriffe haben unterschiedliche numerische Eigenschaften
Leicht steifes System mit wichtigen Schwingungen Sorgfältig ausgewähltes implizites RK- oder Kurbel-Nicocolson-Schema erfordert Stabilität ohne übermäßige Dämpfung
Steady-State-Berechnung durch Pseudo-Zeit implizite oder beschleunigte Iteration vorübergehende Treue kann weniger wichtig sein als Konvergenz

Methodenreihenfolge ist immer noch wichtig

Ein implizites Verfahren erster Ordnung kann einen kleinen Zeitschritt für die Genauigkeit erfordern, selbst wenn die Stabilität einen großen ermöglicht. Eine explizite Methode vierter oder fünfter Ordnung kann für ein reibungsloses, nicht steifes Problem weitaus effizienter sein.

Umgekehrt kann eine explizite Methode höherer Ordnung die schwere Steifheit nicht überwinden, wenn ihr Stabilitätsbereich die relevanten Eigenwerte ausschließt.

Die Auswahl beinhaltet daher zwei getrennte Fragen:

  1. Ist der Stabilitätsbereich für das System und die vorgesehene Schrittweite geeignet?
  2. Ist die Bestellung hoch genug, um den erforderlichen Fehler zu erschwinglichen Kosten zu erfüllen?

Verifizierung durch Zeitschrittverfeinerung

Ein Zeitintegrator sollte getestet werden, indem die Simulation mit kleineren Schritten wiederholt wird. Vergleichen Sie physikalisch relevante Größen wie:

  • Spitzenamplitude
  • Wave Ankunftszeit
  • Schwingungsphase
  • Gesamtmasse oder Energie
  • Reaktionsausbeute
  • Schnittstellenposition
  • Steady-State-Wert

Ändert sich das Ergebnis nach der Halbierung des Zeitschritts signifikant, wurde der ursprüngliche Schritt nicht zeitlich konvergiert.

Zeitliche Tests sollten von der Netzkonvergenz getrennt werden. Die gleichzeitige Verfeinerung von Raum und Zeit macht es schwierig zu bestimmen, welche Fehlerquelle die Änderung verursacht hat.

Häufige Auswahlfehler

Ein häufiger Fehler ist, einen sehr großen impliziten Schritt zu tun, einfach weil die Methode stabil bleibt.

Andere häufige Fehler sind:

  • Verwirrung der lokalen und globalen Ordnung
  • Angenommen, jede implizite Methode ist a-stabil
  • Angenommen, jede A-stabile Methode dämpft steife Modi stark
  • Verwenden von expliziter Diffusion auf einem feinen Netz, ohne dessen Stabilitätsgrenze zu schätzen
  • Verwenden von Rückwärts-Euler, wenn Phasengenauigkeit wichtig ist
  • Verwenden von Crank-Nicocolson für starke Steifheit ohne Überprüfung der zeitlichen Schwingungen
  • Nichtlinear-Solver-Toleranzen in einer impliziten Methode ignorieren
  • Algorithmen mit unterschiedlichen Genauigkeitsstufen vergleichen
  • Stabilisierung ohne Durchführung von Zeitschritt-Konvergenztests
  • Anwenden einer generischen explizit-versus-impliziten Regel auf jede PDE

Ein praktischer Auswahlworkflow

  1. Identifizieren Sie die wichtigen physischen Zeitskalen.
  2. Bestimmen Sie, ob das halbdiskrete System steif ist.
  3. Schätzen Sie Advektion, Diffusion, Reaktion und wellenbasierte Einschränkungen ab.
  4. Entscheiden Sie, ob schnelle Modi aufgelöst oder gedämpft werden können.
  5. Wählen Sie ein Verfahren mit einem geeigneten Stabilitätsbereich.
  6. Wählen Sie eine Bestellung aus, die das Genauigkeitsziel erreichen kann.
  7. Schließen Sie die Kosten für Matrix und nichtlineare Lösung ein.
  8. Verwenden Sie gegebenenfalls adaptive Schritte.
  9. Wiederholen Sie die Simulation mit engeren Toleranzen oder kleineren Schritten.
  10. Vergleichen Sie die Gesamtlaufzeit mit demselben gemessenen Fehler.

Verwandte Anleitungen

weiterlesen

Schlussfolgerung

Explizite und implizite Zeitintegratoren lösen unterschiedliche numerische Probleme. Explizite Methoden bieten kostengünstige Schritte und sind für nicht steife Gleichungen, Wellen und Probleme wirksam, deren physikalische Auflösung bereits kleine Zeitschritte erfordert. Implizite Methoden können einschränkende Stabilitätsgrenzen vermeiden und sind häufig für steife Reaktionen, diffusionsdominierte Modelle und eng gekoppelte Systeme erforderlich.

Stabilität garantiert keine Genauigkeit. Eine implizite Berechnung kann begrenzt bleiben, während schnelle Transienten fehlen, Phasenfehler einführen oder einen Zeitschritt verwenden, der zu groß ist, um die physikalische Trajektorie zu reproduzieren.

Die individuelle Methode ist wichtiger als ihr breites Label. Rückwärts Euler ist stark dämpfend und erster Ordnung genau. Crank-Niconson ist zweiter Ordnung und A-stabil, unterdrückt jedoch keine extrem steifen Modi. RK4 bietet eine hohe Genauigkeit für nicht steife Systeme, hat jedoch einen begrenzten Stabilitätsbereich. IMEX-Methoden kombinieren explizite und implizite Behandlungen, wenn die Operatoren getrennt werden können.

Die richtige Methode ist diejenige, die den erforderlichen Fehler bei den niedrigsten glaubwürdigen Rechenkosten erfüllt. Diese Entscheidung sollte eher durch Stabilitätsanalyse, adaptive Fehlerkontrolle und Zeitschrittverfeinerung als durch die Wörter „explizit“ oder „implizit“ belegt werden.