模拟振动弦：数值发散问题能否避免？

问题背景 #

我正在模拟一维弦上传播的高斯波：初始波速为1.0，x=0位置右侧波速变为0.5。对应的Python代码如下：

import numpy as np
import matplotlib.pyplot as plt

# Parameters
N = 301
x_max = 10.0
x_min = -x_max

v_left = 1.0
v_right = 0.5

t0 = -6.0
dt = 0.02

# Initial conditions
x = np.linspace(x_min, x_max, N)
v = (x < 0) * v_left + (x >= 0) * v_right
y = np.exp(-(x - v_left * t0) ** 2)
dy_dt = y * 2 * (x - v_left * t0) * v_left

dx = x[1] - x[0]

# Prepare for iteration
fig = plt.figure()
d2y_dx2 = np.zeros_like(x)
t = t0
iteration = 0

# Iterate
while t <= 4.0:
    # Approximate wave equation: d^2/dt^2 y = v^2 * d^2/dx^2
    d2y_dx2[1:-1] = np.diff(y, 2) / dx**2
    d2y_dt2 = v**2 * d2y_dx2
    
    y += dt * dy_dt + 0.5 * dt**2 * d2y_dt2
    dy_dt += dt * d2y_dt2
    
    # Check if we should plot
    for test_t in [-4.0, -2.0, -0.5, 0.5, 2.0, 4.0]:
        if abs(t - test_t) < (0.1 * dt):
            plt.plot(x, y, label=f"t={t:.2f}")
    
    # Prepare for next iteration
    iteration += 1
    t = t0 + (iteration * dt)

# Adjust plot
plt.xlim(-5.0, 5.0)
plt.ylim(-1.1, 1.1)
plt.plot([0.0, 0.0], [-1.1, 1.1], color="black")
plt.legend()

plt.show()

问题描述 #

该算法初始阶段运行正常，但最终会出现数值发散的异常结果。尝试提高x轴采样密度或减小dt，却发现发散出现得更早。请问：

1. 这是浮点数精度限制导致的吗？
2. 如何避免该问题？
3. 若无法避免，如何在发散前最大化x/t的采样密度？

1. 发散原因：并非浮点数精度，而是显式积分的稳定性问题 #

你当前使用的是二阶泰勒展开式的显式积分求解波动方程，这类方法受Courant-Friedrichs-Lewy (CFL) 稳定性条件严格约束：

对于波动方程，显式方法稳定的必要条件是 v_max * dt / dx ≤ 1，其中v_max是计算域内的最大波速。

你调整参数后发散更早的核心原因：

• 当提高x采样密度（减小dx）时，若dt不变，v_max * dt / dx的比值会增大；即使减小dt，若dt的减小幅度跟不上dx的减小幅度，比值仍会突破1。
• 一旦该比值超过1，显式方法的数值误差会呈指数级增长，直接导致发散。

2. 避免发散的解决方法 #

方法一：严格遵守CFL条件，同步调整dt与dx #

先计算当前dx下的最大稳定dt：

v_max = max(v_left, v_right)
dt = 0.9 * dx / v_max  # 取0.9而非1是为了留安全余量，避免临界状态的不稳定

确保在调整dx（比如增大N提高采样密度）时，dt同步按比例减小，始终维持v_max * dt / dx < 1。

方法二：改用波动方程专用的稳定显式格式——蛙跳法（Leapfrog Method） #

蛙跳法是针对波动方程设计的显式格式，稳定性更好，且同样满足CFL条件。核心更新逻辑如下：

# 初始化前一步的y值（利用初始条件和dt推导）
d2y_dx2[1:-1] = np.diff(y, 2) / dx**2
y_prev = y - dt * dy_dt + 0.5 * (dt**2) * (v**2 * d2y_dx2)

while t <= 4.0:
    # 蛙跳法更新y值
    y_next = np.zeros_like(y)
    # 内部点更新
    y_next[1:-1] = 2*y[1:-1] - y_prev[1:-1] + (v[1:-1]*dt/dx)**2 * (y[2:] - 2*y[1:-1] + y[:-2])
    # 边界条件（此处采用镜像边界，可根据需求调整）
    y_next[0] = y_next[1]
    y_next[-1] = y_next[-2]
    
    # 更新变量
    y_prev, y = y, y_next
    dy_dt = (y - y_prev) / (2*dt)  # 可选，计算速度项
    
    # 绘图逻辑保持不变...

方法三：改用无条件稳定的隐式格式（如Crank-Nicolson） #

如果需要长时间模拟且不想受CFL条件限制，可采用隐式格式。这类格式需求解线性方程组，实现复杂度更高，但能保证无条件稳定，适合波速突变、长时间演化的场景。

3. 发散前最大化采样密度的方法 #

若暂时不修改数值格式，最大化采样密度需遵循以下规则：

• 固定v_max * dt / dx为一个略小于1的常数（比如0.9），保证稳定性
• 当提高x采样密度（增大N，减小dx）时，同步按比例减小dt，即dt = 0.9 * dx / v_max
• 例如：N从301增至601，dx减半，dt也需减半，这样稳定性条件始终满足，同时x和时间的采样密度均翻倍

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