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

## Main Question or Discussion Point

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 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.

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)

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 !

Last edited:

Related Differential Equations News on Phys.org
DrClaude
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:
Fortran:
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.

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: