Python 中使用有限差分法的基本一维平流方程

Question

嗨，我正在尝试使用有限差分逆风法在 python 中编写一个简单的平流方程。我在网上查找了一个简单的示例，但我发现的代码比我预期的要复杂一些。我有一个用 python 编写的简单代码，其中使用了高斯初始条件。我正在使用的方程和我编写的代码如下：

方程

import numpy as np
import matplotlib.pyplot as plt

#Mesh
Nx = 211 #Number of x grid points
Nt = 1000 #Number of t grid points

dx = 1/(Nx - 1) #Dividing dx into equal dx intervales
dt = 1e-7 #Time step

#Boundary conditions
x= np.linspace(0.1, 1.0, num=Nx)  
t= np.linspace(0.0, Nt*dt, num=Nt)

#Initializing array to zero
u=np.zeros((Nx, Nt)) #Initializing array to zero

#Velocity
c = 1 

#Gaussian initial condition
mu = 0.5
sigma = 0.1

for i in range(1,Nx-1):
    u[i,0] = np.exp(-(x[i]-mu)**2/(2*sigma**2)) / np.sqrt(2*np.pi*sigma**2)


#Looping over all grid points excluding those at the boundary conditions
for j in range(0,Nt-1):   
    for i in range(1,Nx-1):
        
        u[i+1,j] = u[i,j] + c*(((1/x[i]) + (u[i+1,j] - (u[i,j])/dx)))*dt

plt.plot(x, u[:, 0])
plt.plot(x, u[:, 200])
plt.plot(x, u[:, 300])
plt.plot(x, u[:, 400])

由于速度 c>0，我预计绘图会向左移动，并且希望查看从初始高斯起点开始的 4 个不同时间的绘图。但我得到的没有任何意义。我还想指出，当我注释掉 for 循环并执行 plt.plot(x,u[:,0] 时，当我在取消注释 for 循环后再次键入它时，我会看到初始高斯曲线，如预期的那样，我得到了绘图如下所示，我将非常感谢任何人能给我的任何帮助。

输出图

Answer 1

你可以试试这个。特别注意：

我通过将源项 (c/r) 乘以 0 暂时关闭了它。然后，准确的解（增加 Nx）将保持不变地平流。

我已经更正了您的 du/dx 术语和 LHS 上的时间级别。

使用前向欧拉方案，您只需要两个时间级别，边进行边绘制。

我增加了 x 节点的数量以提高准确性。我也已经和dx保持一致了

你仍然需要考虑你的边界条件。

import numpy as np
import matplotlib.pyplot as plt

# Time dependence
Nt = 1000          # number of timesteps
dt = 1e-3

# Space dependence
Nx = 901           # number of nodes
xmin, xmax = 0.1, 1.0
dx = ( xmax - xmin ) / (Nx - 1)
x= np.linspace(xmin, xmax, Nx)

# Equation variable (speed)
c = 1 

# Dependent variable with gaussian initial condition
mu = 0.5
sigma = 0.1
u = np.exp(-(x-mu)**2/(2*sigma**2)) / np.sqrt(2*np.pi*sigma**2)
plt.plot(x, u)

# Loop over all grid points excluding those at the boundary conditions (this IS a problem)
for j in range( 1, Nt + 1 ):
    uold = u.copy()
    for i in range(1, Nx - 1):
        u[i] = uold[i] + c * dt * ( uold[i+1] - uold[i] ) / dx    + 0 * c / x[i]
    if j % 100 == 0: plt.plot( x, u )

plt.show()

输出：