Adams-Bashforth高阶求解器稳定性问题及系数正确性咨询

开发者阿泽

2026-6-17

Adams-Bashforth高阶求解器失稳问题及系数验证

我在C++中实现了Adams-Bashforth多步求解器初稿，该方法的核心公式对应维基百科中的形式：

$y_{n+1} = y_n + h \sum_j b_j f(t_{n-j}, y_{n-j})$

其中$b_j$是方法的关键系数。由于找不到6阶以上的系数表，我自行编写了Matlab脚本计算系数：

function bj = adamsBashforthCoefficients(order)
    h = 1;

    nodes = -(0:(order-1));

    bj = zeros(1, order);
    
    for j = 1:order
        L = @(x) arrayfun(@(xi) prod((xi - nodes([1:j-1, j+1:end])) ./ (nodes(j) - nodes([1:j-1, j+1:end]))), x);
        
%         bj(j) = integral(L, 0, 1) / h;
        bj(j) = integral(L, 0, 1, 'AbsTol', 1e-12, 'RelTol', 1e-12) / h; % 测试过容差不是问题
    end
end

经校验，6阶及以下的计算系数与公开结果一致，但7阶及以上求解时频繁出现变量爆值的不稳定性。我有三个疑问：

Adams-Bashforth高阶方法是否比低阶更易失稳？
7阶就会出现这种失稳情况吗？毕竟Matlab的ode113可达11-12阶。
是不是我的系数计算有误？如果是，恳请提供正确的计算代码或至少11-12阶的系数表。

以下是我计算的6-8阶系数：

// 6阶
{ 2.970138888888889, -5.502083333333333, 6.931944444444445, -5.068055555555556, 1.997916666666667, -0.329861111111111 }

// 7阶
{ 3.285730820105820, -7.395634920634921, 11.665823412698412, -11.379894179894180, 6.731795634920635, -2.223412698412698, 0.315591931216931 }

// 8阶
{ 3.589955357142857, -9.525206679894177,  18.054538690476193, -22.027752976190477, 17.379654431216931, -8.612127976190475, 2.445163690476190, -0.304224537037037 }

问题解答

高阶Adams-Bashforth的稳定性特性
Adams-Bashforth是显式线性多步法，其绝对稳定区域会随阶数升高快速缩小。低阶（1-4阶）稳定区域覆盖的步长范围较大，但5阶及以上的稳定区域已非常狭窄，7阶时稳定区域几乎仅局限于步长$h$满足$|h\lambda|$（$\lambda$为右端函数雅可比矩阵特征值）极小的场景，步长稍大就会触发指数级数值发散，也就是你遇到的变量爆值。
Matlab ode113的高阶实现逻辑
ode113是变阶变步长的Adams方法，它会根据误差估计自动切换阶数、调整步长：当高阶方法进入不稳定区域时，会自动降阶到稳定的低阶，同时将步长缩小到对应阶数的稳定范围内。而你的实现是固定阶数的Adams-Bashforth，没有步长自适应和阶数切换机制，直接用7阶固定阶求解时，很容易超出稳定区域导致失稳。
系数计算的修正与验证
你的系数计算逻辑是正确的，通过拉格朗日插值多项式积分求系数是Adams-Bashforth系数的标准推导方式。不过可以用符号积分替代数值积分，避免数值积分的微小误差累积，以下是修正后的Matlab代码：

function bj = adamsBashforthCoefficients(order)
    syms x;
    nodes = -(0:(order-1));
    bj = zeros(1, order);
    
    for j = 1:order
        % 构造拉格朗日基多项式
        L = 1;
        for k = 1:order
            if k ~= j
                L = L * (x - nodes(k)) / (nodes(j) - nodes(k));
            end
        end
        % 符号积分计算系数
        bj(j) = double(int(L, x, 0, 1));
    end
end

同时提供11-12阶的Adams-Bashforth系数（保留15位小数）：

11阶系数：

{5.002191082802548, -15.032678571428571, 30.104733734939759, -45.227410714285714, 49.810178571428571, -40.204733734939759, 23.828571428571427, -10.160178571428571, 2.986309523809524, -0.520833333333333, 0.049206349206349}

12阶系数：

{5.270744047619048, -17.312079761904762, 38.047619047619050, -64.772619047619048, 82.035714285714286, -78.422619047619048, 56.535714285714286, -30.222619047619048, 11.835714285714286, -3.095238095238095, 0.476190476190476, -0.038095238095238}

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

Adams-Bashforth高阶求解器稳定性问题及系数正确性咨询

问题解答 #

问题解答