- #1

- 4

- 0

## 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 :

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

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

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 :

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

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: