System of ODEs with RK4 & step doubling in Fortran : damping

[Tabatta]

Hello,

I'm recently trying to code a solver for a system of differential equations u'(t) = F(t,u), using a Runge Kutta 4 method with an adaptative stepsize. For this, I'm using the 'step doubling' method, which is the following : suppose that we now the solution u(i) at time t(i). Then, the procedure is the following :

1. Compute the solution u_h (i+1) at time t(i+1) = t(i)+h using the RK4 formula with a stepsize h,
2. Compute the solution u_h/2 (i+1) at time t(i+1) = t(i)+h using the RK4 formula with two steps of size h/2,
3. Compare error = abs(u_h/2 (i+1)-u_h (i+1)) with a chosen precision parameter p : if error>p, then start over with h=h/2, if error<p, then start over with h=2*h.

I added in my code some variable so that the code knows the value of error at the previous value of h, so that there is no loop. Here is my code written in fortran90, I added a lot of comments so I hope it is readable.

!__________________________________________________________________________________________________
!This subroutine solves the differential equation \dot{u} = F(t,u) between (ti,tf), with u(ti)=ui.
!It uses a Runge Kutta 4 method with a variable step size h with a precision pmin that can be specified by the user.
!Description of the arguments :
!1) number_equation, type : integer.
!Gives the number of equation of the system that needs to be solved.
!2) ti, type : double precision.
!Initial time.
!3) tf, type : double precision.
!Final time.
!4) ui, type : double precision, dimension(number_equation).
!Initial value of u.
!5) pmin, type : double precision.
!Precision of the calculation
!6) u, type : double precision, dimension(number_equation,0:10000000).
!Solution of the ODE.
!7) t, type : double precision, dimension(0:10000000).
!Time solution of the ODE
!Idea : computes u(t+1) from u(t) using RK4 with a step h, and with a step h/2. If the difference between the two solution
!is bigger than the precision required, the step h is reduced h=h/2. If it is smaller, then h is increased h=2*h.
!__________________________________________________________________________________________________
SUBROUTINE rk4_variable(number_equation,ti,tf,ui,pmin,u,t)

    IMPLICIT NONE
    INTEGER, INTENT(IN) :: number_equation
    REAL*8 :: t1h,t2h,t3h,t4h,t1hhalf,t2hhalf,t3hhalf,t4hhalf  !Intermediate times in the RK4 formulae
!xh denotes x computed with one step h while xhhalf denotes x computed with two steps h/2
    REAL*8, DIMENSION(number_equation) :: uh,uhhalf !Intermediate solutions using RK4 before choosing h
    REAL*8, DIMENSION(number_equation) :: u1h,u2h,u3h,u4h,u1hhalf,u2hhalf,u3hhalf,u4hhalf
    REAL*8 :: h, hhalf
    REAL*8, DIMENSION(number_equation) :: err,err1
    REAL*8, INTENT(IN) :: ti,tf
    REAL*8, DIMENSION(number_equation), INTENT(IN) :: ui
    INTEGER :: i,j,k
    REAL*8, DIMENSION(0:10000000), INTENT(OUT) :: t
    REAL*8, DIMENSION(number_equation,0:10000000), INTENT(OUT) :: u
    REAL*8, INTENT(IN):: pmin
    LOGICAL :: var_bool1,var_bool2,var_bool3
    t(0)=ti
    u(:,0)=ui
    h=(tf-ti)/999
    i=0    DO WHILE (t(i)<=tf)
        err1=0
        Looph : DO
            k=1
                t1h=t(i)
                t2h=t(i)+h*0.5
                t3h=t(i)+h*0.5
                t4h=t(i)+h

                u1h=u(:,i)
                u2h=u(:,i)+h*0.5*f(t1h,u1h,number_equation)
                u3h=u(:,i)+h*0.5*f(t2h,u2h,number_equation)
                u4h=u(:,i)+h*f(t3h,u3h,number_equation)

                uh=u(:,i)+(h/6)*f(t1h,u1h,number_equation)+(h/3)*f(t2h,u2h,number_equation)&
                +(h/3)*f(t3h,u3h,number_equation)+(h/6)*f(t4h,u4h,number_equation)
!RK4 formula for one step h

                 uhhalf=u(:,i)
                 hhalf=h/2

                 DO j=1,2

                                 t1hhalf=t(i)+hhalf*(j-1)
                                 t2hhalf=t(i)+hhalf*0.5+hhalf*(j-1)
                                 t3hhalf=t(i)+hhalf*0.5+hhalf*(j-1)
                                 t4hhalf=t(i)+hhalf+hhalf*(j-1)

                                 u1hhalf=uhhalf
                                 u2hhalf=uhhalf+hhalf*0.5*f(t1hhalf,u1hhalf,number_equation)
                                 u3hhalf=uhhalf+hhalf*0.5*f(t2hhalf,u2hhalf,number_equation)
                                 u4hhalf=uhhalf+hhalf*f(t3hhalf,u3hhalf,number_equation)

                                  uhhalf=uhhalf+(hhalf/6)*f(t1hhalf,u1hhalf,number_equation)+&
                                     (hhalf/3)*f(t2hhalf,u2hhalf,number_equation)&
                                    +(hhalf/3)*f(t3hhalf,u3hhalf,number_equation)&
                                    +(hhalf/6)*f(t4hhalf,u4hhalf,number_equation)
                 END DO
!RK4 formula for two steps h/2

                err=ABS(uh-uhhalf)

                k=1
                var_bool1=.FALSE. !Var_bool1 indicates if there is some value of k such that err(k)>pmin
                Loopk : DO WHILE (k<=number_equation)
                    IF (err(k)>pmin) THEN
                        var_bool1=.TRUE.

                        EXIT Loopk
                    ELSE
                        k=k+1
                    END IF
                END DO Loopk

                IF (var_bool1.eqv..TRUE.) THEN
                    var_bool2=.FALSE.
                    !Var_bool2 indicates if for k such that err(k)>pmin, this is either the first iteration
                    !of the h loop (ie err1(k)=0) or, in the previous iteration, the error was still
                    !larger that the precision pmin.
                    k=1
                    Loopk2 : DO WHILE (k<=number_equation)
                        IF (err(k)>pmin) THEN
                            IF ((err1(k)==0d0).OR.(err1(k)>pmin)) THEN
                                var_bool2=.TRUE.
                                EXIT Loopk2
                            ELSEIF (err1(k)<=pmin.AND.err1(k)/=0d0) THEN
                                k=k+1
                            END IF
                        ELSE
                            k=k+1
                        END IF
                    END DO Loopk2

                    IF (var_bool2.eqv..TRUE.) THEN
                    !If var_bool2 is true, then the stepsize is too big : reduce it and remember the
                    !value of the error.
                        err1=err
                        h=h/2
                    ELSE
                    !Otherwise, the stepsize is too big but h/2 was perfect : take the value of u
                    !computed with two steps h/2, and remember for the next time step that h=h/2
                    !is fine.
                        u(:,i+1)=uhhalf
                        t(i+1)=t(i)+h
                        h=h/2
                        EXIT Looph
                    END IF
                ELSE !There is no k such that err(k)>pmin
                    var_bool3=.FALSE.
                    !var_bool3 indicates if there is some k such that the error at the previous value
                    !of h (err1(k)) was larger that pmin
                    k=1
                    Loopk3 : DO WHILE (k<=number_equation)
                        IF (err1(k)>pmin) THEN
                            var_bool3=.TRUE.
                            EXIT Loopk3
                        ELSE
                            k=k+1
                        END IF
                    END DO Loopk3

                    IF (var_bool3.eqv..TRUE.) THEN
                    !If var_bool3 is true, then the stepsize is perfect : take the value of u computed
                    !with two steps h/2 since it has the best precision, and keep the value of h.
                        u(:,i+1)=uhhalf
                        t(i+1)=t(i)+h
                        EXIT Looph
                    ELSE
                    !Otherwise, it means that h can be enlarged
                        err1=err
                        h=2*h
                    END IF
                END IF

        END DO Looph
        i=i+1
        print*,'Pourcentage :',t(i)*100/tf,'%'

    END DO

END SUBROUTINE rk4_variable

My problem is that I am currently using this code to solve a system of 4 equations whose solution I know : one of the component should be a perfect sin^2. However, when I use the code such that it is written above, this solution is damped quite rapidly (less than 10 oscillations). The weird part is that when I replace all the "u(:,i+1)=uhhalf" when the value of h is chosen by "u(:,i+1)=uh", there is no damping anymore. I tried a lot of test but I really don't understand where this comes from ... Does anyone have an idea or an explanation ?
Thank you a lot for your time,

Amélie

EDIT : I modified my code in agreement with DrClaude's response, however, there is still a damping (but much slower). Any idea ?
EDIT2 : I modified the precision I was using, and the problem is solved. Thanks !

PS : I apologize if I made any english mistake, it is not my mothertongue !

 

[DrClaude]

Note: I've added "fortran" in the code tags to make it easier to read.

Your coefficients in the half-step version are incorrect:

uhhalf=uhhalf+(h/12)*f(t1hhalf,u1hhalf,number_equation)+(h/6)*f(t2hhalf,u2hhalf,number_equation)&
                    +(h/12)*f(t3hhalf,u3hhalf,number_equation)+(h/6)*f(t4hhalf,u4hhalf,number_equation)

Instead of modifying h many times in that loop, I would have a variable hhalf and use that in the loop. Then you won't make mistakes by incorrectly modifying h.

 

[Tabatta]

Thank you for the "fortran" tag, and I feel kind of silly for not having noticing the error in the coefficients.
I changed them in my code (I switched the 3rd (h/12) and the 4th (h/6)). The result is 'improved', the damping is much slower, but it still exists while it doesn't if I use uh everywhere. Could this be because the error on uhhalf is not controlled ? I will also change the hhalf.

EDIT I found the problem : the precision that I used wasn't good enough. Thank you for your time, you've made my day ahah.