基于有限差分迎风格式的一维平流方程Python实现问题求助

我尝试用有限差分迎风格式在Python中实现一维平流方程，初始条件设为高斯分布。在线找到的示例代码都比较复杂，所以自己写了一份。预期当平流速度c>0时，高斯曲线应该向左移动，但实际运行结果不符合预期——注释掉计算循环时，初始高斯曲线显示正常；取消注释后，输出结果异常。

常见错误点及修正方案 #

以下是这类问题最容易踩的坑，结合代码逐一说明：

1. 迎风格式离散方向错误

当c>0时，平流信息从左向右传播，迎风格式需要取左侧差分，正确的离散公式应为：
u[n+1][i] = u[n][i] - c*dt/dx*(u[n][i] - u[n][i-1])
如果差分方向搞反，波形移动方向会完全不符合预期。

2. 未满足CFL稳定性条件

迎风格式要求Courant数C = c*dt/dx ≤ 1，如果C>1会直接导致数值不稳定，出现振荡、发散等异常结果。

3. 数组覆盖问题

循环中直接修改原数组u，会导致后续计算提前使用更新后的值，破坏迭代逻辑，必须用临时数组存储每一步的新状态。

4. 边界条件未处理

左边界i=0无法直接使用i-1索引，必须补充边界条件（比如周期性边界、Dirichlet边界），否则会导致边界点数值异常，影响整体波形。

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt

# 参数设置
nx = 100
L = 2.0
dx = L / (nx - 1)
c = 1.0  # 平流速度>0
dt = 0.001
nt = 100

# 检查Courant数，确保满足稳定性条件
C = c * dt / dx
print(f"Courant数: {C:.3f}")
assert C <= 1, "Courant数大于1，数值不稳定，请调小dt或增大dx"

# 初始条件：高斯分布
x = np.linspace(0, L, nx)
u = np.exp(-((x - 0.5)**2)/0.02)
u_initial = u.copy()

# 初始化临时数组存储新状态
u_new = u.copy()

# 迎风格式迭代计算
for n in range(nt):
    # 内部点计算：左差分格式
    u_new[1:] = u[1:] - c * dt / dx * (u[1:] - u[:-1])
    # 周期性边界条件（模拟无限域，保证波形连续）
    u_new[0] = u[0] - c * dt / dx * (u[0] - u[-1])
    # 交换数组，准备下一轮迭代
    u, u_new = u_new, u

# 绘图对比
plt.figure(figsize=(8, 4))
plt.plot(x, u_initial, label='初始状态', linestyle='--')
plt.plot(x, u, label=f'迭代{nt}步后 (c>0)')
plt.xlabel('x')
plt.ylabel('u')
plt.legend()
plt.title('一维平流方程迎风格式模拟（修正后）')
plt.show()

修正说明 #

• 新增Courant数检查，从根源避免数值不稳定；
• 使用临时数组u_new存储每一步结果，防止原数组覆盖导致的计算错误；
• 添加周期性边界处理，解决左边界无索引的问题，保证波形移动的连续性；
• 严格遵循c>0时的左差分格式，匹配信息传播方向，确保高斯曲线向左移动。

内容的提问来源于stack exchange，提问作者user24293678