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

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

Beginnend mit der Advektionsgleichung ist konservative Form,

$$\frac{\partial u}{\partial t} = -\frac{\partial (\boldsymbol{v} u)}{\partial x} + s(x,t)$$

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

$$\frac{u_{j}^{n+1} - u_{j}^{n}}{\Delta t} = -\boldsymbol{v} \left[ \frac{1-\beta}{2\Delta x} \left( u_{j+1}^{n} - u_{j-1}^{n} \right) + \frac{\beta}{2\Delta x} \left( u_{j+1}^{n+1} - u_{j-1}^{n+1} \right) \right] + 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 = \frac{\boldsymbol{v}}{2}\frac{\Delta t}{\Delta x}$$

gibt,

$$-\beta r\phi_{j-1}^{n+1} + \phi_{j}^{n+1} + \beta r\phi_{j+1}^{n+1} = (1-\beta)r\phi_{j-1}^{n} + \phi_{j}^{n} - (1-\beta)r\phi_{j+1}^{n}$$

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

$$\begin{pmatrix} 1 & \beta r & & & 0 \\ -\beta r & 1 & \beta r & & \\ & \ddots & \ddots & \ddots & \\ & & -\beta r & 1 & \beta r \\ 0 & & & -\beta r & 1 \\ \end{pmatrix} \begin{pmatrix} u_1^{n+1} \\ u_2^{n+1} \\ \vdots \\ u_{J-1}^{n+1} \\ u_{J}^{n+1} \\ \end{pmatrix} = \begin{pmatrix} 1 & -(1 - \beta)r & & & 0 \\ (1 - \beta)r & 1 & -(1 - \beta)r & & \\ & \ddots & \ddots & \ddots & \\ & & (1 - \beta)r & 1 & -(1 - \beta)r \\ 0 & & &(1 - \beta)r & 1 \\ \end{pmatrix} \begin{pmatrix} u_1^{n} \\ u_2^{n} \\ \vdots \\ u_{J-1}^{n} \\ u_{J}^{n} \\ \end{pmatrix}$$ 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.$\beta=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.

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