Fortran四阶Adams-Bashforth方法运行出现异常NaN输出问题排查

问题根因 #

1. 四阶Adams-Bashforth公式括号书写错误 #

四阶Adams-Bashforth的正确公式为：
y(r+1) = y(r) + (h/24) * (55*f(t_r,y_r) - 59*f(t_{r-1},y_{r-1}) + 37*f(t_{r-2},y_{r-2}) - 9*f(t_{r-3},y_{r-3}))
你代码里的括号只包裹了第一个55*f(t(r),y(r))项，相当于只有第一项乘了h/24，后面三个系数项没有被h/24缩放，直接参与计算，必然导致数值异常增大最终溢出为NAN。

2. t数组越界访问（未定义行为） #

你在adams_bashford_four_method子程序中定义的t是长度为4的局部数组，仅分配了t(1)~t(4)的内存空间：

real:: y(1000001), t(4), ab_calculation

你循环给t数组赋值的逻辑是把最近4个时间点存入t(1)~t(4)，但后续调用f时却用了t(r)、t(r-1)这类下标：当循环变量r >=5时，t(r)已经超出了数组的合法下标范围，属于未定义行为。
这也能解释你碰到的反常现象：加打印语句后程序运行正常，是因为打印语句改变了程序的栈内存布局，越界访问时刚好读取到了合法的数值，属于巧合，不具备可复现性。 #

修复方案 #

修改adams_bashford_four_method子程序的逻辑即可：

1. 修正公式括号位置，保证所有系数项都被h/24缩放
2. 调用f时使用t数组的合法下标1~4
   修正后的子程序代码如下：

subroutine adams_bashford_four_method(t0,y,h,n)
integer:: r, i
real:: y(1000001), t(4)
do r = 4,n
    ! t(4)对应r时刻，t(3)对应r-1时刻，t(2)对应r-2时刻，t(1)对应r-3时刻
    do i = 4,1,-1
        t(5-i) = t0 + (r-i)*h
    end do
    ! 修正括号位置，用t数组合法下标访问时间值
    y(r+1) = y(r) + (h/24) * (55*f(t(4),y(r)) - 59*f(t(3),y(r-1)) + 37*f(t(2),y(r-2)) - 9*f(t(1),y(r-3)))
end do
end subroutine adams_bashford_four_method

修复后使用测试参数t0 = 1、y0 = 0.5、t1 = 2、h = 0.02运行，即可得到正确结果，不会再出现NAN问题。 #

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