求解IVP波动方程：非预期反射排查及驻波动画实现

一、非预期延迟反射的成因排查 #

• 并非solve_ivp核心机制问题，而是数值边界处理不当导致的伪反射：
  • 波动方程离散化后，若边界处未设置严格的物理约束（如固定边界的u=0），solve_ivp在求解离散ODE系统时，边界节点的导数计算会依赖邻域节点，当波传播到边界时，数值上会出现“伪泄漏”，之后由于离散格式的特性，这些泄漏的波会被数值边界反射回来，延迟恰好等于波往返边界的时间。
  • T=10s时，波尚未完成往返，伪反射未显现；T=20s时往返完成，伪反射叠加到物理反射上，导致相位偏移异常。
• 验证方向：检查边界节点的微分方程是否被正确约束——固定边界需强制u[0] = 0和u[-1] = 0，而非让其参与默认的差分计算。

二、双边界反射驻波的实现代码 #

以下是基于solve_ivp的波动方程驻波动画实现，左侧注入正弦入射波，双边界设为固定边界（反射带π相位）：

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# 物理参数
L = 10.0  # 空间长度
c = 1.0   # 波速
nx = 200  # 空间节点数
x = np.linspace(0, L, nx)
dx = x[1] - x[0]

# 波动方程转化为一阶ODE系统：y = [u, du/dt]，dy/dt = [du/dt, d²u/dt²] = [v, c² d²u/dx²]
def wave_eq(t, y):
    u = y[:nx]
    v = y[nx:]
    
    # 计算二阶空间导数（中心差分，边界用特殊处理保证固定约束）
    d2u_dx2 = np.zeros_like(u)
    # 内部节点用中心差分
    d2u_dx2[1:-1] = (u[2:] - 2*u[1:-1] + u[:-2]) / dx**2
    # 左侧固定边界：强制u[0]为入射波，导数匹配入射波变化率
    ω = 2*np.pi*c/(2*L)  # 驻波基频
    u_in = np.sin(ω*t)
    u[0] = u_in
    d2u_dx2[0] = (u[1] - 2*u_in + u_in) / dx**2  # 避免边界伪反射
    # 右侧固定边界：强制u[-1] = 0
    d2u_dx2[-1] = (0 - 2*u[-1] + u[-2]) / dx**2
    
    # 左侧入射波的速度项修正
    v[0] = ω*np.cos(ω*t)
    
    return np.concatenate([v, c**2 * d2u_dx2])

# 初始条件：静止无位移
u0 = np.zeros(nx)
v0 = np.zeros(nx)
y0 = np.concatenate([u0, v0])

# 求解时间范围
T = 40.0
t_eval = np.linspace(0, T, 200)

# 用Radau求解器处理刚性波动方程
sol = solve_ivp(wave_eq, [0, T], y0, t_eval=t_eval, method='Radau')

# 动画绘制
fig, ax = plt.subplots()
line, = ax.plot(x, sol.y[:nx, 0])
ax.set_ylim(-2, 2)
ax.set_xlabel('空间位置')
ax.set_ylabel('波幅')
ax.set_title('双边界固定的正弦入射波驻波动画')

def update(frame):
    line.set_ydata(sol.y[:nx, frame])
    return line,

ani = FuncAnimation(fig, update, frames=len(t_eval), interval=50, blit=True)
plt.show()

关键说明 #

• 把波动方程拆解为一阶ODE系统，将位移u和速度v=du/dt合并为状态向量，适配solve_ivp的求解要求。
• 固定边界通过强制节点位移值、修正二阶导数计算逻辑，彻底避免数值伪反射。
• 选择Radau求解器，针对波动方程这类刚性系统的求解稳定性更强。
• 入射波频率设为驻波基频，确保反射波与入射波叠加后形成稳定驻波。

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