请求协助修正Runge-Kutta-Fehlberg法三阶ODE的Fortran95求解代码

DevOps 小满

2026-6-20

关于Runge-Kutta-Fehlberg法求解三阶线性ODE的Fortran代码问题

我需要用Runge-Kutta-Fehlberg（RKF45）方法求解三阶常微分方程（ODE）：
$$y''' = -2y''+y'+2y$$
初始条件为 $y(0)=3$、$y'(0)=-2$、$y''(0)=6$，精确解为 $y=\exp(-2t)+\exp(-t)+\exp(t)$。

我绘制了计算结果与精确解的对比图，但发现作为线性ODE，计算结果和精确解的匹配度没有达到预期。以下是我编写的Fortran 95代码，请帮忙指出问题并给出修正建议。

program ft

    !testing 45 rk method by examining y''' = -2y''+y'+2y, y''(0)=6, y'(0)=-2, and y(0)=3.
    !the exact solution should be y = exp(-2t)+exp(-t)+exp(t)
    implicit none
    integer, parameter :: i = SELECTED_REAL_KIND (10,200)
    REAL(i) :: t, val(3), dt, newval(3), newdt
    logical :: redo
    external runge

    t = 0.0

    !val = (y,y',y'')
    !initial condition
    val = (/3.0,-2.0,6.0/)
    !step size
    dt = 1.0E-4

    open (10, file = 'RK_test.dat', status = 'new')
    do while (t<4.0)
        redo = .true.
        write(10,*) t, val(1), val(2), val(3)

        do while (redo)
            call runge(t, dt, val, newdt, redo, newval)
            dt = newdt
        end do
        val = newval
        t = t + dt
    end do
    close(10)

end program ft

subroutine runge(t, dt, val, newdt, recalc, valp1)
    implicit none
    integer, parameter :: i = SELECTED_REAL_KIND (10,200)
    REAL(i), parameter :: deltacalc(2,6) = reshape([16.0/135, 0.0, 6656.0/12825, 28561.0/56430, -9.0/50, 2.0/55, &amp;
                25.0/216, 0.0, 1408.0/2565,   2197.0/4104,  -1.0/5,  0.0], shape(deltacalc), order=[2,1])
    REAL(i), parameter :: coe(6,6) = reshape ([0.0,         0.0,          0.0,          0.0,         0.0,      0.0, &amp;
                                           1.0/4,       0.0,          0.0,          0.0,         0.0,      0.0, &amp;
                                           3.0/32,      9.0/32,       0.0,          0.0,         0.0,      0.0, &amp;
                                           1932.0/2197, -7200.0/2197, 7296.0/2197,  0.0,         0.0,      0.0, &amp;
                                           439.0/216,   -8.0,         3680.0/513,   -845.0/4104, 0.0,      0.0, &amp;
                                           -8.0/27,     2.0,          -3544.0/2565, 1859.0/4104, -11.0/40, 0.0], &amp;
                                           shape(coe), order=[2,1])
    REAL(i) :: valbar(3), k(6,3), delta(3), error(3), var(3), f(3)
    REAL(i), intent(in) :: t, dt, val(3)
    REAL(i), intent(out) :: newdt, valp1(3)
    logical, intent(out) :: recalc
    integer :: a, b, c, d
    External func

    k = 0
    valbar = val
    valp1 = val
    var = val

    do a = 1,6
        do b = 1,3
            do c = 1,6
                do d = 1,3
                    var(d) = var(d) + coe(a,c)*k(c,d)
                end do
            end do
            call func(var, f)
            k(a,b) = dt*f(b)
        end do
    end do
    recalc = .false.

    do b = 1, 3
        if (val(b)==0.0) then
            error(b) = 1.0E-3
        else
            error(b) = 10.0**(log10(abs(val(b)))-3.0)
        end if


        do a = 1, 6
            valbar(b) = valbar(b) + k(a,b)*deltacalc(1,a)
            valp1(b) = valp1(b) + k(a,b)*deltacalc(2,a)
        end do

        delta(b) = 0.84*(error(b)/(abs(valbar(b) - valp1(b))/dt))**(1.0/4) 
        
        if (abs(valbar(b) - valp1(b))/dt>error(b)) then
            recalc = .true.
            !print*, &quot;oooops&quot;
        end if
    end do


    newdt = minval(delta)*dt
    if (newdt > 0.01) then
        newdt = 0.01
    end if

end subroutine runge

subroutine func(var, out)
    implicit none
    integer, parameter :: i = SELECTED_REAL_KIND (10,200)
    REAL(i), intent(in) :: var(3)
    REAL(i), intent(out) :: out(3)

    out(1) = var(2) 
    out(2) = var(3)
    out(3) = -2*var(3) + var(2) + 2*var(1)
end subroutine func

问题分析与修正建议

1. RKF45系数矩阵维度与映射错误

你定义的coe(6,6)和deltacalc(2,6)在reshape时使用order=[2,1]，导致系数的行/列映射完全错位。RKF45的Butcher表中，coe应为6个阶段对应前5个k值的权重（第1阶段无前置k），正确的系数定义应去掉冗余的0列，并使用默认行优先的reshape顺序。

修正后的系数定义：

REAL(dp), parameter :: b4(6) = [16.0_dp/135, 0.0_dp, 6656.0_dp/12825, 28561.0_dp/56430, -9.0_dp/50, 2.0_dp/55]  ! 四阶解权重
REAL(dp), parameter :: b5(6) = [25.0_dp/216, 0.0_dp, 1408.0_dp/2565, 2197.0_dp/4104, -1.0_dp/5, 0.0_dp]        ! 五阶解权重
REAL(dp), parameter :: coe(6,5) = [ &
    0.0_dp,         0.0_dp,          0.0_dp,          0.0_dp,         0.0_dp, &
    1.0_dp/4,       0.0_dp,          0.0_dp,          0.0_dp,         0.0_dp, &
    3.0_dp/32,      9.0_dp/32,       0.0_dp,          0.0_dp,         0.0_dp, &
    1932.0_dp/2197, -7200.0_dp/2197, 7296.0_dp/2197,  0.0_dp,         0.0_dp, &
    439.0_dp/216,   -8.0_dp,         3680.0_dp/513,   -845.0_dp/4104, 0.0_dp, &
    -8.0_dp/27,     2.0_dp,          -3544.0_dp/2565, 1859.0_dp/4104, -11.0_dp/40 &
]

2. k值计算循环逻辑错误

当前代码中var变量未在每个阶段开始时重置为初始状态，而是持续累加所有k值，导致每个阶段的状态计算完全错误。正确逻辑应为：每个阶段开始时将var重置为初始状态，再累加当前阶段对应的前置k值。

修正后的k值计算循环：

do a = 1,6
    var = val  ! 每个阶段重置为初始状态
    do c = 1,a-1  ! 仅累加前a-1个k值
        var = var + coe(a,c)*k(c,:)
    end do
    call func(var, f)
    k(a,:) = dt*f(:)
end do

3. 误差估计与步长调整错误

误差计算中不需要除以dt，RKF45的误差是四阶解与五阶解的直接差值：err = abs(valbar(b) - valp1(b))
步长调整公式的指数应为1/5（五阶方法的局部截断误差为O(dt^5)），而非1/4
误差容限计算可简化为error(b) = 1e-3 * max(abs(val(b)), 1e-6)，避免对数运算的精度问题

修正后的误差与步长调整部分：

do a = 1,3
    error(a) = 1e-3_dp * max(abs(val(a)), 1e-6_dp)
    valbar = val + matmul(k, b4)
    valp1 = val + matmul(k, b5)
    err = abs(valbar(a) - valp1(a))
    
    if (err > error(a)) then
        recalc = .true.
    end if
    
    delta(a) = 0.8_dp * (error(a)/err)**(1.0_dp/5.0_dp)
end do

4. 输出逻辑问题

当前代码在循环开始时写入上一步的时间与状态，应调整为计算完新状态后，写入新的时间t+dt和新状态newval，确保输出数据的时间与状态对应。

修正后的完整代码

program ft
    ! Testing RKF45 method for y''' = -2y''+y'+2y, y(0)=3, y'(0)=-2, y''(0)=6
    ! Exact solution: y = exp(-2t)+exp(-t)+exp(t)
    implicit none
    integer, parameter :: dp = SELECTED_REAL_KIND(10,200)
    REAL(dp) :: t, val(3), dt, newval(3), newdt
    logical :: redo
    external runge

    t = 0.0_dp
    val = [3.0_dp, -2.0_dp, 6.0_dp]  ! (y, y', y'')
    dt = 1.0E-4_dp

    open(10, file='RK_test.dat', status='replace')
    write(10,*) t, val(1), val(2), val(3)  ! 写入初始条件

    do while (t < 4.0_dp)
        redo = .true.
        do while (redo)
            call runge(t, dt, val, newdt, redo, newval)
            dt = newdt
        end do
        t = t + dt
        val = newval
        write(10,*) t, val(1), val(2), val(3)  ! 写入新的时间和状态
    end do
    close(10)

end program ft

subroutine runge(t, dt, val, newdt, recalc, valp1)
    implicit none
    integer, parameter :: dp = SELECTED_REAL_KIND(10,200)
    ! RKF45 Butcher table coefficients
    REAL(dp), parameter :: b4(6) = [16.0_dp/135, 0.0_dp, 6656.0_dp/12825, 28561.0_dp/56430, -9.0_dp/50, 2.0_dp/55]
    REAL(dp), parameter :: b5(6) = [25.0_dp/216, 0.0_dp, 1408.0_dp/2565, 2197.0_dp/4104, -1.0_dp/5, 0.0_dp]
    REAL(dp), parameter :: coe(6,5) = [ &
        0.0_dp,         0.0_dp,          0.0_dp,          0.0_dp,         0.0_dp, &
        1.0_dp/4,       0.0_dp,          0.0_dp,          0.0_dp,         0.0_dp, &
        3.0_dp/32,      9.0_dp/32,       0.0_dp,          0.0_dp,         0.0_dp, &
        1932.0_dp/2197, -7200.0_dp/2197, 7296.0_dp/2197,  0.0_dp,         0.0_dp, &
        439.0_dp/216,   -8.0_dp,         3680.0_dp/513,   -845.0_dp/4104, 0.0_dp, &
        -8.0_dp/27,     2.0_dp,          -3544.0_dp/2565, 1859.0_dp/4104, -11.0_dp/40 &
    ]
    REAL(dp) :: k(6,3), err, error(3), delta(3), var(3), f(3)
    REAL(dp), intent(in) :: t, dt, val(3)
    REAL(dp), intent(out) :: newdt, valp1(3)
    logical, intent(out) :: recalc
    integer :: a, c
    external func

    recalc = .false.
    k = 0.0_dp

    ! Compute all k stages
    do a = 1,6
        var = val
        do c = 1,a-1
            var = var + coe(a,c)*k(c,:)
        end do
        call func(var, f)
        k(a,:) = dt*f(:)
    end do

    ! Compute 4th order and 5th order solutions
    valp1 = val + matmul(k, b5)
    do a = 1,3
        error(a) = 1e-3_dp * max(abs(val(a)), 1e-6_dp)
        err = abs( (val + matmul(k, b4))(a) - valp1(a) )
        
        if (err > error(a)) then
            recalc = .true.
        end if
        
        delta(a) = 0.8_dp * (error(a)/err)**(1.0_dp/5.0_dp)
    end do

    ! Adjust step size, clamp to max 0.01 and avoid overshooting t=4.0
    newdt = minval(delta)*dt
    newdt = min(newdt, 0.01_dp)
    if (t + newdt > 4.0_d

请求协助修正Runge-Kutta-Fehlberg法三阶ODE的Fortran95求解代码

问题分析与修正建议 #

1. RKF45系数矩阵维度与映射错误 #

2. k值计算循环逻辑错误 #

3. 误差估计与步长调整错误 #

4. 输出逻辑问题 #