我试图通过使用分步傅里叶方法求解非线性薛定谔方程来模拟孤子传播。以下是从教科书上复制的 Matlab 代码示例。

alpha = 0;
beta_2 = 1;
gamma = 1;

A = N*sech(T);   % initial pulse shape
L = max(size(A));
delta_omega = 1/L/delta_T*2*pi;
omega = (-L/2:1:L/2-1)*delta_omega;

solution = A;   % analytical solution (same as original it is periodic)

step_num_hist = round(10.^linspace(1,3.5,20));   % total number of steps
avg_err = zeros(3,max(size(step_num_hist)));

end_pt = 0.5*pi;   % the end of simulation, one soliton period

for step_num = step_num_hist
    h = end_pt/step_num;   % step size
    A_t = A;
    for n=1:step_num
        A_f = fftshift(fft(A_t));     
        A_f = A_f.*exp(-alpha*(h/2)-1i*beta_2/2*omega.^2*(h/2));
        A_t = ifft(ifftshift(A_f));

        A_t = A_t.*exp(1i*gamma*(abs(A_t).^2*h));

        A_f = fftshift(fft(A_t));
        A_f = A_f.*exp(-alpha*(h/2)-1i*beta_2/2*omega.^2*(h/2));
        A_t = ifft(ifftshift(A_f));
    end
    comp_result = abs(A_t);
    avg_err(1,m) = sum(abs(comp_result.^2 - abs(solution).^2))/max(size(solution))/max(abs(solution).^2);   % avg relative intensity error, found in a paper
    m = m+1
end

该代码似乎产生了一些有意义的东西，所以我认为它没有错。问题是，当步数增加（步长h减小）超过某个数字（在我的情况下约为 1e-3）时，计算解和解析解之间的误差不会进一步减小。这发生在 Matlab 和 Python 中。我还尝试实现了一些我在文献中找到的高阶方法。它们确实收敛得更快，但“平均相对强度误差”不会低于 1e-3 左右，而我发现比较各种算法的论文可以下降到 1e-7 甚至 1e-10（总步数也是 10e4) 的数量级。

所以我的问题是为什么会这样？为什么误差最初随着步长减小，但在步长达到某个点后停止减小？什么可能出错？