Korteweg-De Vries方程n>1时数值不稳定的原因排查

我希望求解初始条件为$U(x,0) = \frac{n(n+1)}{\cosh^2(x)}$的Korteweg-De Vries（KdV）方程，但当整数$n>1$时，系统出现数值不稳定现象，请问可能的原因是什么？

复现代码 #

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

# 参数设置
L = 40    # 区间长度
N = 256   # 网格点数
dx = L / N  # 空间步长
x = np.linspace(-L/2, L/2, N)  # 离散空间域
dt = 0.1  # 时间步长
T = 20    # 终止时间
n = 1     # 初始条件参数n

# 初始条件 u(x,0) = n*(n+1)/cosh^2(x)
u0 = (n * (n + 1)) / np.cosh(x)**2

# 计算KdV方程的导数
def derivs(t, u):
    # 使用np.roll实现周期性边界条件
    uxx = np.roll(u, -1) - 2 * u + np.roll(u, 1)  # 二阶导数
    uxx /= dx**2
    uxxx = np.roll(uxx, -1) - np.roll(uxx, 1)  # 三阶导数
    uxxx /= 2 * dx
    # KdV方程：u的时间导数
    return -6 * u * (np.roll(u, -1) - np.roll(u, 1)) / (2 * dx) - uxxx

# 用四阶Runge-Kutta方法求解KdV方程
sol = solve_ivp(derivs, [0, T], u0, method='RK45', t_eval=np.arange(0, T, dt))

# 动画展示解
plt.figure()
for i in range(len(sol.t)):
    plt.clf()  # 清除上一帧图像
    plt.plot(x, sol.y[:, i])  # 绘制t[i]时刻的解
    plt.title(f"KdV方程解（t={sol.t[i]:.1f}）")
    plt.xlabel("x")  # x轴标签
    plt.ylabel("u(x,t)")  # y轴标签
    plt.ylim([-0.5, n*(n+1)+0.5])  # 设置y轴范围以优化显示
    plt.pause(0.1)  # 短暂停顿实现动画效果

plt.show()  # 显示最终图像窗口

可能的不稳定原因 #

• 空间离散格式精度不足：当前代码用二阶中心差分计算二阶导数，一阶中心差分（借助np.roll）计算一阶和三阶导数，这类低精度离散格式对高频数值噪声的抑制能力弱。当n>1时，初始解的峰值更高、梯度更陡，截断误差被快速放大，直接引发数值不稳定。KdV方程的三阶色散项对离散格式精度要求极高，低阶格式无法满足需求。
• 时间步长违反稳定性条件：KdV方程的显式数值解法存在严格的CFL类稳定性约束。n增大时，初始解的特征速度（与u的幅值正相关）显著提升，你固定的dt=0.1可能超过了当前空间分辨率（N=256）下的稳定时间步长上限，导致RK45方法出现数值发散。
• 边界条件假设不合理：代码用np.roll实现周期性边界条件，但初始条件$\frac{n(n+1)}{\cosh^2(x)}$在区间$[-20,20]$两端仅近似为0，并非严格周期连续。n越大，解的传播速度越快，会更早触及边界，这种非严格周期的边界会引入反射误差，进一步加剧数值振荡。
• 显式方法的固有局限：RK45属于显式时间积分方法，对于包含三阶导数的色散方程，其稳定区域有限。当n增大时，解的高阶导数幅值剧增，问题刚度显著提升，显式方法容易突破稳定边界，引发数值振荡甚至完全发散。

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