执行带阶跃修正的 Runge-Kutta 方法时的无限循环

Question

我正在为 ODE 逼近实现 Runge-Kutta-4 方法，并使用阶跃校正程序。这是代码：

function RK4 (v,h,cant_ec) !Runge-Kutta 4to Orden
    real(8),dimension(0:cant_ec)::RK4,v
    real::h
    integer::cant_ec

    real(8),dimension(0:cant_ec)::k1,k2,k3,k4

    k1 = h*vprima(v)
    k2 = h*vprima(v+k1/2.0)
    k3 = h*vprima(v+k2/2.0)
    k4 = h*vprima(v+k3)
    v = v + (k1+2.0*k2+2.0*k3+k4)/6.0 !la aproximación actual con paso h

    RK4 = v 

    end function RK4

    
    subroutine RK4h1(v,h,xf,tol,cant_ec) !Runge-Kutta con corrección de paso por método 1
    real(8),dimension(0:cant_ec)::v
    real::h,tol,xf
    integer::cant_ec,i,n

    real(8),dimension(0:cant_ec)::v1,v2
    real(8)::error

    n = int((xf-v(0))/h +0.5)
    open(2,file = "derivada.txt",status="replace")
    error = 2*tol

    do i = 1,n, 1
        do while(error > tol)
            v1 = RK4(v,h,cant_ec)
            v2 = RK4(v,h/2,cant_ec)
            v2 = v2 + RK4(v+v2,h/2,cant_ec)
            error = MAXVAL(ABS(v1-v2))
            
            if (error > tol) then 
                h = h/2
            end if
        
        end do 
    end do
    
    write(*,*)v1 
    write(2,*)v1
    close(2,status="keep")

    call system("gnuplot -persist 'derivada.p'")
    end subroutine Rk4h1

其中 h 是步长，v 是 cant_ec 分量的向量，对应于 ODE 的阶数（即：v(0) = x,v(1) = y,v(2) = y', etc), tol 是误差容限，xf 是 x 区间的终点（假设它从 0 开始）。所有这些值都是在子程序调用之前由用户输入的。此特定函数的初始值为 y(0) = -1。其他一切都由用户在运行脚本时定义。微分方程由下式给出：

function vprima(v,x,y) !definición de la función derivada
    real(8),dimension(0:cant_ec)::v,vprima

    
    vprima(0) = 1.0
    vprima(1) = (-v(1)/(v(0)**2+1))
    end function vprima

注意到在主程序上发生了这个分配：

v(0) = x
v(1) = y

其中 x 和 y 是函数的初始值，由用户提供。

我的问题是，当我调用 RK4h1 时，脚本似乎陷入无限循环。

如有任何帮助或线索，我们将不胜感激。谢谢。

Answer 1

v2 = v2 + RK4(v+v2,h/2,cant_ec)是错误的，应该是v2 = RK4(v2,h/2,cant_ec)，因为RK4的结果是下一个点，而不是更新到下一个点。由于错误计算因此是错误的，因此步长会无限减少。约50次减少后，RK4步不再推进状态，增量太小

固定步长可变步长会出问题

内部循环没有任何意义。整体效果是每次步长减小后 i 都会增加 1。所以理论上，如果 n<50 程序应该终止，但最终状态非常接近初始状态。

局部误差应该与tol*h进行比较，这样全局误差就与tol成正比。

如果局部误差变得太小，还应该有增加h的指令。

请参阅 How to perform adaptive step size using Runge-Kutta fourth order (Matlab)? 以了解使用 RK4 和双步步长控制的另一个示例。

执行带阶跃修正的 Runge-Kutta 方法时的无限循环

Infinite loop while implementing Runge-Kutta method with step correction

fortran

gfortran

numerical-methods