本文详解如何通过在数值求解器中直接施加齐次狄利克雷边界条件（u=0），替代后处理反射模拟，从而彻底消除伪反射、相位错误和延迟回波，获得物理准确的波传播与反射动画。
本文详解如何通过在数值求解器中直接施加齐次狄利克雷边界条件（u=0），替代后处理反射模拟，从而彻底消除伪反射、相位错误和延迟回波，获得物理准确的波传播与反射动画。
在使用 scipy.integrate.solve_ivp 求解一维波动方程初值问题（IVP）时，若未显式处理空间边界，数值解会因截断域外信息而产生非物理的伪反射（unwanted reflection）——表现为波抵达边界后“延迟反弹”、相位不翻转（应为 π 相位跃变）、甚至出现多重重叠干扰。这并非 solve_ivp 的缺陷，而是源于边界条件缺失：默认情况下，有限差分近似将边界点视为内部节点，导致其导数计算依赖于不存在的外部网格点，等效于施加了隐式（且错误的）周期性或自由边界条件。
正确解法是：在 ODE 右端函数 wave_eq 中，对状态变量的时间导数显式强加物理边界约束。对于固定端（如弦两端固定），即齐次狄利克雷条件 $ u(0,t) = u(L,t) = 0 $，其直接推论是边界速度为零：
$$
\frac{\partial u}{\partial t}(0,t) = 0, \quad \frac{\partial u}{\partial t}(L,t) = 0.
$$
因此，在构造 dudt = y[N:] 后，只需两行代码即可强制满足该条件：
dudt[0] = 0 dudt[-1] = 0 # 等价于 dudt[N-1] = 0
此举确保了数值解在每一步积分中严格满足 $ u(0,t) \equiv 0 $ 和 $ u(L,t) \equiv 0 $，反射行为由波动方程本身自然导出，具备正确的 $ \pi $ 相位反转与瞬时响应。
✅ 关键优势：
- 消除所有延迟伪反射（如原代码中 T=20s 时出现的“二次回波”）；
- 自动保证反射波振幅守恒、相位正确（固定端反射必反相）；
- 无需手动滚动、镜像或叠加反射波形，代码简洁鲁棒；
- 可无缝扩展至其他边界条件：例如自由端（$ \partial u/\partial x = 0 $）只需修改 uxx 计算时采用单侧差分或设置 uxx[0] = uxx[-1] = 0。
以下为精简后的核心求解函数，已集成正确边界处理：
def wave_eq(t, y, c, L, dx):
N = len(y) // 2
u = y[:N]
# 使用二阶中心差分近似 u_xx（推荐用 np.gradient 两次，或自定义二阶差分）
uxx = np.gradient(np.gradient(u, dx), dx)
dudt = y[N:] # ∂u/∂t 当前值
dudt[0] = 0 # 固定左端：∂u/∂t(0,t) = 0
dudt[-1] = 0 # 固定右端：∂u/∂t(L,t) = 0
du_tdt = c**2 * uxx # ∂²u/∂t² = c² ∂²u/∂x²
return np.concatenate([dudt, du_tdt])
注意事项：
- 空间步长 dx 需满足 CFL 条件：$ c \cdot dt / dx \leq 1 $（当前 c=4, dt=0.1, dx=0.01 得 CFL=40，严重超限！建议将 dt 降至 0.002 或更小，否则高频误差会激化数值振荡）；
- 初始条件 u_t0 必须与边界条件兼容（如 u_t0[0]=u_t0[-1]=0），否则初始瞬间违反物理约束；
- 若需模拟行波持续输入（如左端驱动正弦波），则应改用非齐次边界条件，并在 dudt[0] 处赋值为驱动信号的导数（如 dudt[0] = omega * A * cos(omega * t)）。
通过此方法，您可稳定生成从单高斯脉冲到连续正弦驱动、再到驻波形成的全过程动画——所有反射均由偏微分方程与边界条件自洽决定，真正还原波动物理本质。
