Wie kann die Advektionsgleichung mit der Crank-Nicolson-Methode diskretisiert werden?

8

Die Advektionsgleichung muss diskretisiert werden, um für die Crank-Nicolson-Methode verwendet zu werden. Kann mir jemand zeigen, wie das geht?

Pandoragami
quelle
1
Willkommen bei SciComp! Der Umfang Ihrer Frage passt sehr gut zu dieser Seite. Um gute Antworten zu erhalten, sollten Sie jedoch genauer sein. Bitte geben Sie an, was Sie insbesondere nicht verstehen. Ihr Code sieht gut strukturiert und dokumentiert aus, aber zur Beantwortung Ihrer Frage reicht möglicherweise ein kleineres Code-Snippet aus.
Jan
Dies kann das Debuggen vereinfachen, wenn Sie einen Testfall verwenden, in dem Ihr Eingabevektor beispielsweise 5 Elemente enthält, und den Code mit einem Debugger wie gdb und ddd durchlaufen. Dies kann dazu beitragen, die Fehlerquelle einzugrenzen. Ich denke, dass die meisten Fragen zum Debuggen von Code hier nicht sehr gut funktionieren, da der Großteil der Arbeit darin besteht, herauszufinden, wo der Fehler überhaupt liegt. Sobald Sie es gefunden haben, ist die Erklärung häufig (aber nicht immer) einfach. Können Sie einen Komponententest durchführen, um herauszufinden, ob es in Tridiagonal Fehler gibt, die dieses Verhalten verursachen könnten?
Geoff Oxberry
Schauen Sie sich dieses Beispiel im Wikipedia-Eintrag zur Crank-Nicolson-Methode an. Wenn Sie und auf Null setzen, wird dies zu einem einfachen Advektionsproblem. Es bleibt die Randbedingungen einzuarbeiten ...D xkDx
Jan
Sorry Leute, aber die Crank-Nicolson-Methode ist für ein Advektionsproblem völlig ungeeignet. Stabilität und Genauigkeit der lokalen Differentialnäherung garantieren leider keine Konsistenz. Sie erhalten die Lösung eines anderen Problems. In einigen trivialen Fällen haben Sie vielleicht Glück, aber die Advektionsgleichung ist im Allgemeinen unversöhnlich. Überprüfen Sie die Lehrbücher zur numerischen Analyse. Das Crank-Nicolson-Schema wird niemals verwendet. Die Leute müssen es wegen der Stabilität immer noch bei einigen technischen Problemen anwenden, aber dort müssen sie die Lösung kontrollieren, während sie sich entwickelt, und heuristisch sein, um erro

Antworten:

19

Beginnend mit der Advektionsgleichung ist konservative Form,

ut=(vu)x+s(x,t)

Die Crank-Nicolson-Methode besteht aus einer zeitlich gemittelten zentrierten Differenz.

ujn+1ujnΔt=v[1β2Δx(uj+1nuj1n)+β2Δx(uj+1n+1uj1n+1)]+s(x,t)

In Bezug auf die Notation beziehen sich Indizes auf Punkte im Raum und hochgestellte Indizes auf Zeitpunkte.

Die Punkte bei liegen in der Zukunft: Sie sind unbekannt. Wir müssen nun die obige Gleichung neu anordnen, so dass alle bekannten auf der rechten und die unbekannten auf der rechten Seite stehen.n+1

Die Substitution vornehmen,

r=v2ΔtΔx

gibt,

βrϕj1n+1+ϕjn+1+βrϕj+1n+1=(1β)rϕj1n+ϕjn(1β)rϕj+1n

Dies ist die Advektionsgleichung, die nach der Crank-Nicolson-Methode diskretisiert wurde. Sie können es als Matrixgleichung schreiben,

(1βr0βr1βrβr1βr0βr1)(u1n+1u2n+1uJ1n+1uJn+1)=(1(1β)r0(1β)r1(1β)r(1β)r1(1β)r0(1β)r1)(u1nu2nuJ1nuJn)
Wenn Sie 1/2 einstellen, erhalten Sie eine trapezförmige Integration in der Zeit. Für Crank-Nicolson ist dies genau das, was Sie wollen.β=1/2

Ein paar warnende Worte. Dies ist die grundlegende Lösung, die Sie wollten, aber Sie müssen eine Art Randbedingung für ein gut gestelltes Problem angeben. Außerdem ist Crank-Nicolson nicht unbedingt die beste Methode für die Advektionsgleichung. Es ist genau zweiter Ordnung und bedingungslos stabil , was fantastisch ist. Es wird jedoch (wie bei allen Schablonen mit zentrierter Differenz) eine Störschwingung erzeugen, wenn Sie sehr scharfe Spitzenlösungen oder Anfangsbedingungen haben.

Ich habe den folgenden Code für Sie in Python geschrieben, damit Sie loslegen können. Der Code löst die Advektionsgleichung für eine anfängliche Gaußsche Kurve, die sich mit konstanter Geschwindigkeit nach rechts bewegt.

Gaußsche Kurve mit konstanter Geschwindigkeit nach rechts

from __future__ import division
from scipy.sparse import spdiags
from scipy.sparse.linalg import spsolve
import numpy as np
import pylab

def make_advection_matrices(z, r):
    """Return matrices A and M for advection equations"""
    ones = np.ones(len(z))
    A = spdiags( [-beta*r, ones, beta*r], (-1,0,1), len(z), len(z) )
    M = spdiags( [(1-beta) * r, ones, -(1-beta) * r], (-1,0,1), len(z), len(z) )
    return A.tocsr(), M.tocsr()

def plot_iteration(z, u, iteration):
    """Plot the solver progress"""
    pylab.plot(z, u, label="Iteration %d" % iteration)

# Set up basic constants
beta = 0.5
J = 200 # total number of mesh points
z = np.linspace(-10,10,J) # vertices
dz = abs(z[1]-z[0]) # space step
dt = 0.2    # time step
v = 2 * np.ones(len(z)) # velocity field (constant)
r = v / 2 * dt / dz

# Initial conditions (peak function)
gaussian = lambda z, height, position, hwhm: height * np.exp(-np.log(2) * ((z - position)/hwhm)**2)
u_init = gaussian(z, 1, -3, 2)

A, M = make_advection_matrices(z, r)
u = u_init
for i in range(10):
    u = spsolve(A, M * u)
    plot_iteration(z, u, i)

pylab.legend()
pylab.show()
Boyfarrell
quelle
3
Eigentlich habe ich auch Ihre frühere Frage gesehen, aber sie war so allgemein, dass ich sie nicht beantworten konnte (als Sie eine Code-Seite gepostet haben). Nach meiner Erfahrung sind die Leute sehr hilfreich, wenn Sie auf dieser Seite gute Fragen stellen. Gute Reise!
Boyfarrell
Ich habe nur gescherzt.
Pandoragami
@boyfarrel Gibt es eine Chance, dass Sie eine C ++ / C-Version davon haben? Es ist in Ordnung, wenn nein. Ich benutze Matlab nicht viel und ich habe keine Lust, es zu lernen. Sogar Fortan wäre besser.
Pandoragami
2
<unangemessene Kommentare gelöscht>
Paul