Infinite loop while implementing Runge-Kutta method with step correction

Question

I'm implementing the Runge-Kutta-4 method for ODE approximation with a step correction procedure. Here's the code:

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

Where h is the step size, v is a vector of cant_ec components that corresponds to the order of the ODE (that is: v(0) = x,v(1) = y,v(2) = y', etc), tol is the tolerance of error and xf is the end of the x interval (assuming it starts at 0). All these values are inputted by the user before the subroutine call. The initial values given for this particular function are y(0) = -1. Everything else is defined by the user when running the script. The differential equation is given by:

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

noting that on the main program this assignment occurs:

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

where x and y are the initial values of the function, given by the user.

My issue is, the script seems to get stuck on an infinite loop when I call RK4h1.

Any help, or clue, would be appreciated. Thank you.

We need a complete program; otherwise, we're just guessing. For example, my initial guess is that your requested tolerance is too small. — steve, May 01 '21 at 03:26
I didn't add the complete program because this is pretty much it (at least to the issue, the rest of the file are the declaration of other methods, like Euler and Runge-Kutta-Fuhlberg). Still, I added some details about the initial values to the question which I realize now are important to know. — Morganuz, May 01 '21 at 03:36
Expecting others to add the missing main program and guessing the exact set of values that are causing your issues is a fairly good way of getting help. Good luck. — steve, May 01 '21 at 06:36
Please see [mcve]. The point is to enable others to try your code easily without unnecessary work. Also, be aware that `real(8)` is ugly and non-portable https://stackoverflow.com/questions/838310/fortran-90-kind-parameter — Vladimir F Героям слава, May 01 '21 at 09:22
I understand. I'll edit the post accordingly. Thank you for the correction. — Morganuz, May 01 '21 at 14:46

Lutz Lehmann · Accepted Answer · 2021-05-01T11:45:02.767

2

v2 = v2 + RK4(v+v2,h/2,cant_ec) is wrong, it should be v2 = RK4(v2,h/2,cant_ec), as the result of RK4 is the next point, not the update to the next point. As the error calculation is thus wrong, the step size gets reduced indefinitely. After some 50 reductions, the RK4 step will no longer advance the state, the increment being too small.

It will lead to problems if you have a fixed step number with a variable step size.

The inner loop does not make any sense whatsoever. The overall effect is that after each step size reduction i gets incremented by one. So theoretically, if n<50 the program should terminate, but with a final state very close to the initial state.

The local error should be compared to tol*h, so that the global error becomes proportional to tol.

There should also be an instruction that increases h if the local error becomes too small.

See How to perform adaptive step size using Runge-Kutta fourth order (Matlab)? for another example of using RK4 with step-size control by double steps.

edited May 01 '21 at 11:45

answered May 01 '21 at 07:09

Lutz Lehmann

25,219
2
22
51

Could you expand on "So theoretically, if n<50 the program should terminate, but with a final state very close to the initial state."? Should I control the end of the iterations with the end value of x? – Morganuz May 01 '21 at 14:52
Yes, the time loop should be something like `do while(v(0) – Lutz Lehmann May 01 '21 at 15:26
Thank you. Another question, what would be considered a "too small" local error? Is there anyway to have a lower bound? – Morganuz May 01 '21 at 15:45
Because of order 4, the unit step error changes with factor 16 if the step size changes with factor 2. So a gap of 20 or 50 between upper and lower bound would be advisable. Or use a narrower gap for a correspondingly smaller step size factor. – Lutz Lehmann May 01 '21 at 20:06

Infinite loop while implementing Runge-Kutta method with step correction

1 Answers1