Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

  1. Nov 12, 2015 #1

    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 uh (i+1) at time t(i+1) = t(i)+h using the RK4 formula with a stepsize h,
    2. Compute the solution uh/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(uh/2 (i+1)-uh (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.

    Code (Fortran):

    !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)

        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

        DO WHILE (t(i)<=tf)
            Looph : DO


    !RK4 formula for one step h


                     DO j=1,2



                     END DO
    !RK4 formula for two steps h/2


                    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

                            EXIT Loopk
                        END IF
                    END DO Loopk

                    IF (var_bool1.eqv..TRUE.) THEN
                        !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.
                        Loopk2 : DO WHILE (k<=number_equation)
                            IF (err(k)>pmin) THEN
                                IF ((err1(k)==0d0).OR.(err1(k)>pmin)) THEN
                                    EXIT Loopk2
                                ELSEIF (err1(k)<=pmin.AND.err1(k)/=0d0) THEN
                                END IF
                            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.
                        !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.
                            EXIT Looph
                        END IF
                    ELSE !There is no k such that err(k)>pmin
                        !var_bool3 indicates if there is some k such that the error at the previous value
                        !of h (err1(k)) was larger that pmin
                        Loopk3 : DO WHILE (k<=number_equation)
                            IF (err1(k)>pmin) THEN
                                EXIT Loopk3
                            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.
                            EXIT Looph
                        !Otherwise, it means that h can be enlarged
                        END IF
                    END IF

            END DO Looph
            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,


    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 !
    Last edited: Nov 12, 2015
  2. jcsd
  3. Nov 12, 2015 #2


    User Avatar

    Staff: Mentor

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

    Your coefficients in the half-step version are incorrect:
    Code (Fortran):

    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.
  4. Nov 12, 2015 #3
    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.
    Last edited: Nov 12, 2015
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook