Ist Crank-Nicolson ein stabiles Diskretisierungsschema für die Reaktions-Diffusions-Advektions-Gleichung (Konvektionsgleichung)?

26

Ich bin mit den üblichen Diskretisierungsverfahren für PDEs nicht sehr vertraut. Ich weiß, dass Crank-Nicolson ein beliebtes Verfahren zur Diskretisierung der Diffusionsgleichung ist. Ist das auch eine gute Wahl für den Advektionssemester?

Ich interessiere mich für die Lösung der Reaktions-Diffusions-Advektions- Gleichung,

ut+(vuDu)=f

Dabei ist der Diffusionskoeffizient der Substanz und die Geschwindigkeit.Duv

Für meine spezifische Anwendung kann die Gleichung in der Form geschrieben werden,

ut=D2ux2Diffusion+vuxAdvection (convection)+f(x,t)Reaction

Hier ist das Crank-Nicolson-Schema, das ich angewendet habe,

ujn+1ujnΔt=D[1β(Δx)2(uj1n2ujn+uj+1n)+β(Δx)2(uj1n+12ujn+1+uj+1n+1)]+v[1α2Δx(uj+1nuj1n)+α2Δx(uj+1n+1uj1n+1)]+f(x,t)

Beachten Sie die Begriffe und . Dies ermöglicht es dem Schema, zwischen folgenden Optionen zu wechseln:αβ

  • β=α=1/2 1/2 Crank-Niscolson,
  • β=α=1 ist vollständig implizit
  • β=α=0 ist vollständig explizit

Die Werte können unterschiedlich sein, so dass der Diffusionsbegriff Crank-Nicolson und der Advektionsbegriff etwas anderes sein kann. Was ist der stabilste Ansatz, was würden Sie empfehlen?

Boyfarrell
quelle

Antworten:

15

Dies ist eine gut umschriebene Frage, die sehr nützlich ist, um sie zu verstehen. Korrok verweist Sie zu Recht auf die von Neumann-Analyse und das Buch von LeVeque. Ich kann dem noch ein bisschen mehr hinzufügen. Ich würde gerne eine ausführliche Antwort schreiben, aber im Moment habe ich nur Zeit für eine kurze:

Mit 1/2 erhalten Sie eine Methode, die für beliebig große Schrittgrößen absolut stabil und genau zweiter Ordnung ist. Das Verfahren ist jedoch nicht L- stabil, so dass sehr hohe Frequenzen nicht gedämpft werden, was unphysikalisch ist.α=β=1/2

Mit Sie eine Methode, die ebenfalls bedingungslos stabil, aber nur genau 1. Ordnung ist. Diese Methode ist sehr dissipativ. Es ist L- stabil.α=β=1

Wenn Sie , kann Ihre Methode so verstanden werden, dass Sie eine additive Runge-Kutta-Methode auf die Mittendifferenz-Semidiskretisierung anwenden . Die Stabilitäts- und Genauigkeitsanalyse für solche Methoden ist erheblich komplizierter. Eine sehr schöne Abhandlung über solche Methoden finden Sie hier .αβ

Welcher Ansatz zu empfehlen ist, hängt stark von der Größe von , der Art der Ausgangsdaten und der Genauigkeit ab, die Sie suchen. Wenn eine sehr geringe Genauigkeit akzeptabel ist, ist ein sehr robuster Ansatz. Wenn mittel oder groß ist, ist das Problem diffusionsdominiert und sehr steif. In der Regel liefert 1/2 gute Ergebnisse. Wenn sehr klein ist, kann es vorteilhaft sein, eine explizite Methode und ein Upwinding höherer Ordnung für die konvektiven Terme zu verwenden.Dα=β=1Dα=β=1/2D

David Ketcheson
quelle
Eine sehr aufschlussreiche Antwort, danke! Gibt es eine Möglichkeit, die verschiedenen Regime zu definieren, in denen Diffusion und Advektion dominiert werden? Anders als die Größe der Begriffe zu vergleichen? Zum Beispiel nur durch den Vergleich von Koeffizienten? Was bedeutet der Fachbegriff L-Stabilität? Jeder empfiehlt dieses Buch, ich muss es kaufen!
Boyfarrell
Das Kriterium, das ich Ihnen gegeben habe, betrifft nur die Koeffizienten. Kurz gesagt bedeutet L-Stabilität, dass hohe Frequenzen stark gedämpft werden.
David Ketcheson
Wenn also eine glatte Funktion ist (wie in dem Sinne, dass es keine hochfrequenten Fourier-Komponenten hat), ist Crank-Nicolson eine gute Wahl. Wenn jedoch scharfe Kanten hat, ist eine gute Wahl. u(x)u(x)β=1
Boyfarrell
Das ist eine vernünftige, wenn auch sehr grobe Verallgemeinerung. Diese Auswahl funktioniert zumindest, wenn Sie nicht sehr genau sein müssen.
David Ketcheson
10

Im Allgemeinen sollten Sie eine implizite Methode für parabolische Gleichungen (den Diffusionsteil) verwenden - explizite Schemata für parabolische PDE benötigen einen sehr kurzen Zeitschritt, um stabil zu sein. Umgekehrt sollten Sie für den hyperbolischen Teil (Advektion) eine explizite Methode verwenden, da diese billiger ist und die Symmetrie des zu lösenden linearen Systems nicht durch die Verwendung eines impliziten Diffusionsschemas stört. In diesem Fall möchten Sie zentrierte Differenzen wie vermeiden und zu einseitigen Differenzen wechseln. aus Stabilitätsgründen.(uj+1uj1)/2Δt(ujuj1)/Δt

Ich würde vorschlagen, dass Sie sich Randy Leveques Buch oder Dale Durrans Buch für "von Neumann-Stabilitätsanalyse" ansehen . Dies ist ein allgemeiner Ansatz zur Feststellung der Stabilität Ihres Diskretisierungsschemas, vorausgesetzt, Sie haben periodische Randbedingungen. (Es gibt auch einen guten Wiki - Artikel hier .)

Die Grundidee ist anzunehmen, dass Ihre diskrete Approximation eine Summe von ebenen Wellen , wobei die Wellenzahl und die Frequenz ist. Sie stopfen eine ebene Welle in Ihre Annäherung an die PDE und beten, dass sie nicht explodiert. Wir können die ebene Welle als und wollen sicherstellen, dass .ei(kjΔxωnΔt)kωξneikjΔx|ξ|1

Betrachten Sie zur Veranschaulichung die gewöhnliche Diffusionsgleichung mit vollständig impliziter Differenzierung:

ujn+1ujnΔt=Duj1n+12ujn+1+uj+1n+1Δx2

Wenn wir in einer ebenen Welle substituieren, dann dividieren durch und , erhalten wir die GleichungξneikjΔx

ξ1Δt=DeikΔx2+eikΔxΔx2ξ

Räumen Sie jetzt ein bisschen auf und wir bekommen:

ξ=11+2DΔtΔx2(1coskΔx) .

Dies ist immer weniger als eins, so dass Sie im klaren sind. Wenden Sie dies für das explizite, zentrierte Schema für die Advektionsgleichung an:

ujn+1ujnΔt=vuj1nuj+1n2Δx

und sehen , was Sie erhalten. (Diesmal wird es einen Imaginärteil geben.) Sie werden feststellen, dass , was traurige Zeiten sind. Daher meine Ermahnung, dass du es nicht benutzt. Wenn Sie dies tun können, sollten Sie keine großen Schwierigkeiten haben, ein stabiles Schema für die vollständige Advektions-Diffusions-Gleichung zu finden.ξ|ξ|2>1

Das heißt, ich würde ein vollständig implizites Schema für den Diffusionsteil verwenden. Ändern Sie die Differenzierung im Advektionsteil in wenn und wenn und wählen Sie einen Zeitschritt, sodass . (Dies ist die Courant-Friedrichs-Lewy-Bedingung .) Sie ist nur in erster Ordnung korrekt, daher möchten Sie möglicherweise Diskretisierungsschemata höherer Ordnung nachschlagen, wenn dies Sie betrifft.ujuj1v>0ujuj+1v<0VΔt/Δx1

Daniel Shapero
quelle
Das ist eine wirklich detaillierte Antwort, danke.
Boyfarrell
Diese Antwort berücksichtigt nur Diskretisierungen, die auf den zeitlichen Vorwärts- und Rückwärts-Euler-Methoden basieren. Die Frage ist über Crank-Nicholson.
David Ketcheson