!**************************************************************************
!* Collection of Fortran 90 subroutines to Integrate a System of Ordinary *
!* Differential Equations By the Runge-Kutta-Fehlberg method (simple or *
!* double precision). *
!* ---------------------------------------------------------------------- *
!* REFERENCE: H A Watts and L F Shampine, *
!* Sandia Laboratories, *
!* Albuquerque, New Mexico. *
!* *
!* See demonstration program trkf45.f90. *
!**************************************************************************
subroutine fehl ( f, neqn, y, t, h, yp, f1, f2, f3, f4, f5, s )
!
!**************************************************************************
!
! FEHL takes one Fehlberg fourth-fifth order step (single precision).
!
! Discussion:
!
! FEHL integrates a system of NEQN first order ordinary differential
! equations of the form
! dY(i)/dT = F(T,Y(1),---,Y(NEQN))
! where the initial values Y and the initial derivatives
! YP are specified at the starting point T.
!
! FEHL advances the solution over the fixed step H and returns
! the fifth order (sixth order accurate locally) solution
! approximation at T+H in array S.
!
! The formulas have been grouped to control loss of significance.
! FEHL should be called with an H not smaller than 13 units of
! roundoff in T so that the various independent arguments can be
! distinguished.
!
! Author:
!
! H A Watts and L F Shampine,
! Sandia Laboratories,
! Albuquerque, New Mexico.
!
! RKF45 is primarily designed to solve non-stiff and mildly stiff
! differential equations when derivative evaluations are inexpensive.
! RKF45 should generally not be used when the user is demanding
! high accuracy.
!
! Parameters:
!
! Input, external F, a subroutine of the form
! subroutine f(t,y,yp)
! to evaluate the derivatives.
! YP(I) = dY(I) / dT
!
! Input, integer NEQN, the number of equations to be integrated.
!
! Input, real Y(NEQN), the current value of the dependent variable.
!
! Input, real T, the current value of the independent variable.
!
! Input, real H, the step size to take.
!
! Input, real YP(NEQN), the current value of the derivative of the
! dependent variable.
!
! Output, real F1(NEQN), real F2(NEQN), real F3(NEQN), real F4(NEQN),
! real F5(NEQN) are arrays of dimension NEQN which are needed
! for internal storage.
!
! Output, real S(NEQN), the computed estimate of the solution at T+H.
!
implicit none
!
integer neqn
!
real ch
external f
real f1(neqn)
real f2(neqn)
real f3(neqn)
real f4(neqn)
real f5(neqn)
real h
real s(neqn)
real t
real y(neqn)
real yp(neqn)
!
ch = h / 4.0E+00
f5(1:neqn) = y(1:neqn) + ch * yp(1:neqn)
call f ( t + ch, f5, f1 )
ch = 3.0E+00 * h / 32.0E+00
f5(1:neqn) = y(1:neqn) + ch * ( yp(1:neqn) + 3.0E+00 * f1(1:neqn) )
call f ( t + 3.0E+00 * h / 8.0E+00, f5, f2 )
ch = h / 2197.0E+00
f5(1:neqn) = y(1:neqn) + ch * ( 1932.0E+00 * yp(1:neqn) &
+ ( 7296.0E+00 * f2(1:neqn) - 7200.0E+00 * f1(1:neqn) ) )
call f ( t + 12.0E+00 * h / 13.0E+00, f5, f3 )
ch = h / 4104.0E+00
f5(1:neqn) = y(1:neqn) + ch * ( ( 8341.0E+00 * yp(1:neqn) &
- 845.0E+00 * f3(1:neqn) ) + ( 29440.0E+00 * f2(1:neqn) &
- 32832.0E+00 * f1(1:neqn) ) )
call f ( t + h, f5, f4 )
ch = h / 20520.0E+00
f1(1:neqn) = y(1:neqn) + ch * ( ( -6080.0E+00 * yp(1:neqn) &
+ ( 9295.0E+00 * f3(1:neqn) - 5643.0E+00 * f4(1:neqn) ) ) &
+ ( 41040.0E+00 * f1(1:neqn) - 28352.0E+00 * f2(1:neqn) ) )
call f ( t + h / 2.0E+00, f1, f5 )
!
! Ready to compute the approximate solution at T+H.
!
ch = h / 7618050.0E+00
s(1:neqn) = y(1:neqn) + ch * ( ( 902880.0E+00 * yp(1:neqn) &
+ ( 3855735.0E+00 * f3(1:neqn) - 1371249.0E+00 * f4(1:neqn) ) ) &
+ ( 3953664.0E+00 * f2(1:neqn) + 277020.0E+00 * f5(1:neqn) ) )
return
end
subroutine fehl_d ( f, neqn, y, t, h, yp, f1, f2, f3, f4, f5, s )
!
!*******************************************************************************
!
!! FEHL_D takes one Fehlberg fourth-fifth order step (double precision).
!
!
! Discussion:
!
! FEHL integrates a system of NEQN first order ordinary differential
! equations of the form
! dY(i)/dT = F(T,Y(1),---,Y(NEQN))
! where the initial values Y and the initial derivatives
! YP are specified at the starting point T.
!
! FEHL advances the solution over the fixed step H and returns
! the fifth order (sixth order accurate locally) solution
! approximation at T+H in array S.
!
! The formulas have been grouped to control loss of significance.
! FEHL should be called with an H not smaller than 13 units of
! roundoff in T so that the various independent arguments can be
! distinguished.
!
! Author:
!
! H A Watts and L F Shampine,
! Sandia Laboratories,
! Albuquerque, New Mexico.
!
! RKF45 is primarily designed to solve non-stiff and mildly stiff
! differential equations when derivative evaluations are inexpensive.
! RKF45 should generally not be used when the user is demanding
! high accuracy.
!
! Parameters:
!
! Input, external F, a subroutine of the form
! subroutine f(t,y,yp)
! to evaluate the derivatives.
! YP(I) = dY(I) / dT
!
! Input, integer NEQN, the number of equations to be integrated.
!
! Input, double precision Y(NEQN), the current value of the dependent
! variable.
!
! Input, double precision T, the current value of the independent variable.
!
! Input, double precision H, the step size to take.
!
! Input, double precision YP(NEQN), the current value of the derivative
! of the dependent variable.
!
! Output, double precision F1(NEQN), F2(NEQN), F3(NEQN), F4(NEQN),
! F5(NEQN) are arrays of dimension NEQN which are needed
! for internal storage.
!
! Output, double precision S(NEQN), the computed estimate of the solution
! at T+H.
!
implicit none
!
integer neqn
!
double precision ch
external f
double precision f1(neqn)
double precision f2(neqn)
double precision f3(neqn)
double precision f4(neqn)
double precision f5(neqn)
double precision h
double precision s(neqn)
double precision t
double precision y(neqn)
double precision yp(neqn)
!
ch = h / 4.0D+00
f5(1:neqn) = y(1:neqn) + ch * yp(1:neqn)
call f ( t + ch, f5, f1 )
ch = 3.0D+00 * h / 32.0D+00
f5(1:neqn) = y(1:neqn) + ch * ( yp(1:neqn) + 3.0D+00 * f1(1:neqn) )
call f ( t + 3.0D+00 * h / 8.0D+00, f5, f2 )
ch = h / 2197.0D+00
f5(1:neqn) = y(1:neqn) + ch * ( 1932.0D+00 * yp(1:neqn) &
+ ( 7296.0D+00 * f2(1:neqn) - 7200.0D+00 * f1(1:neqn) ) )
call f ( t + 12.0D+00 * h / 13.0D+00, f5, f3 )
ch = h / 4104.0D+00
f5(1:neqn) = y(1:neqn) + ch * ( ( 8341.0D+00 * yp(1:neqn) &
- 845.0D+00 * f3(1:neqn) ) + ( 29440.0D+00 * f2(1:neqn) &
- 32832.0D+00 * f1(1:neqn) ) )
call f ( t + h, f5, f4 )
ch = h / 20520.0D+00
f1(1:neqn) = y(1:neqn) + ch * ( ( -6080.0D+00 * yp(1:neqn) &
+ ( 9295.0D+00 * f3(1:neqn) - 5643.0D+00 * f4(1:neqn) ) ) &
+ ( 41040.0D+00 * f1(1:neqn) - 28352.0D+00 * f2(1:neqn) ) )
call f ( t + h / 2.0D+00, f1, f5 )
!
! Ready to compute the approximate solution at T+H.
!
ch = h / 7618050.0D+00
s(1:neqn) = y(1:neqn) + ch * ( ( 902880.0D+00 * yp(1:neqn) &
+ ( 3855735.0D+00 * f3(1:neqn) - 1371249.0D+00 * f4(1:neqn) ) ) &
+ ( 3953664.0D+00 * f2(1:neqn) + 277020.0D+00 * f5(1:neqn) ) )
return
end
subroutine rkf45 ( f, neqn, y, t, tout, relerr, abserr, iflag, work, iwork )
!
!*******************************************************************************
!
!! RKF45 carries out the Runge-Kutta-Fehlberg method (single precision).
!
!
! Author:
!
! H A Watts and L F Shampine,
! Sandia Laboratories,
! Albuquerque, New Mexico.
!
! RKF45 is primarily designed to solve non-stiff and mildly stiff
! differential equations when derivative evaluations are inexpensive.
! RKF45 should generally not be used when the user is demanding
! high accuracy.
!
! Abstract:
!
! RKF45 integrates a system of NEQN first order ordinary differential
! equations of the form:
!
! dY(i)/dT = F(T,Y(1),Y(2),...,Y(NEQN))
!
! where the Y(1:NEQN) are given at T.
!
! Typically the subroutine is used to integrate from T to TOUT but it
! can be used as a one-step integrator to advance the solution a
! single step in the direction of TOUT. On return, the parameters in
! the call list are set for continuing the integration. The user has
! only to call RKF45 again (and perhaps define a new value for TOUT).
! Actually, RKF45 is an interfacing routine which calls subroutine
! RKFS for the solution. RKFS in turn calls subroutine FEHL which
! computes an approximate solution over one step.
!
! Reference:
!
! E. Fehlberg,
! Low-order Classical Runge-Kutta Formulas with Stepsize Control,
! NASA Technical Report R-315.
!
! L F Shampine, H A Watts, S Davenport,
! Solving Non-stiff Ordinary Differential Equations - The State of the Art,
! Sandia Laboratories Report SAND75-0182,
! To appear in SIAM Review.
!
! Parameters:
!
! Input, external F, a subroutine of the form
! subroutine f ( t, y, yp )
! real t
! real y(*)
! real yp(*)
! to evaluate the derivatives YP(I) = dY(I) / dT
!
! Input, integer NEQN, the number of equations to be integrated.
!
! Input/output, real Y(NEQN), the solution vector at T.
!
! Input/output, real T, the independent variable.
!
! Input, real TOUT, the output point at which solution is desired.
!
! Input, real RELERR, ABSERR, the relative and absolute error tolerances
! for the local error test. At each step the code requires:
! abs(local error) <= relerr*abs(y) + abserr
! for each component of the local error and solution vectors
!
! Output, integer IFLAG, indicator for status of integration.
!
! Workspace, real WORK(3+6*NEQN), an array to hold information internal
! to RKF45 which is necessary for subsequent calls.
!
! Workspace, integer IWORK(5), an array used to hold information internal
! to RKF45 which is necessary for subsequent calls.
!
!
! first call
!
! The user must provide storage in his calling program for the arrays
! in the call list - y(neqn) , work(3+6*neqn) , iwork(5) ,
! declare f in an external statement, supply subroutine f(t,y,yp) and
! initialize the following parameters-
!
! neqn -- number of equations to be integrated. (neqn >= 1)
!
! y(*) -- vector of initial conditions
!
! t -- starting point of integration , must be a variable
!
! tout -- output point at which solution is desired.
! t=tout is allowed on the first call only, in which case
! rkf45 returns with iflag=2 if continuation is possible.
!
! relerr,abserr -- relative and absolute local error tolerances
! which must be non-negative. relerr must be a variable while
! abserr may be a constant. the code should normally not be
! used with relative error control smaller than about 1.e-8 .
! to avoid limiting precision difficulties the code requires
! relerr to be larger than an internally computed relative
! error parameter which is machine dependent. in particular,
! pure absolute error is not permitted. if a smaller than
! allowable value of relerr is attempted, rkf45 increases
! relerr appropriately and returns control to the user before
! continuing the integration.
!
! iflag -- +1,-1 indicator to initialize the code for each new
! problem. normal input is +1. the user should set iflag=-1
! only when one-step integrator control is essential. in this
! case, rkf45 attempts to advance the solution a single step
! in the direction of tout each time it is called. since this
! mode of operation results in extra computing overhead, it
! should be avoided unless needed.
!
!
! output from rkf45
!
! y(*) -- solution at t
! t -- last point reached in integration.
! iflag = 2 -- integration reached tout. indicates successful retur
! and is the normal mode for continuing integration.
! =-2 -- a single successful step in the direction of tout
! has been taken. normal mode for continuing
! integration one step at a time.
! = 3 -- integration was not completed because relative error
! tolerance was too small. relerr has been increased
! appropriately for continuing.
! = 4 -- integration was not completed because more than
! 3000 derivative evaluations were needed. this
! is approximately 500 steps.
! = 5 -- integration was not completed because solution
! vanished making a pure relative error test
! impossible. must use non-zero abserr to continue.
! using the one-step integration mode for one step
! is a good way to proceed.
! = 6 -- integration was not completed because requested
! accuracy could not be achieved using smallest
! allowable stepsize. user must increase the error
! tolerance before continued integration can be
! attempted.
! = 7 -- it is likely that rkf45 is inefficient for solving
! this problem. too much output is restricting the
! natural stepsize choice. use the one-step integrator
! mode.
! = 8 -- invalid input parameters
! this indicator occurs if any of the following is
! satisfied - neqn <= 0
! t=tout and iflag /= +1 or -1
! relerr or abserr < 0.
! iflag == 0 or < -2 or > 8
! work(*),iwork(*) -- information which is usually of no interest
! to the user but necessary for subsequent calls.
! work(1),...,work(neqn) contain the first derivatives
! of the solution vector y at t. work(neqn+1) contains
! the stepsize h to be attempted on the next step.
! iwork(1) contains the derivative evaluation counter.
!
!
! subsequent calls
!
! RKF45 returns with all information needed to continue
! the integration. if the integration reached tout, the user need onl
! define a new tout and call RKF45 again. In the one-step integrator
! mode (iflag=-2) the user must keep in mind that each step taken is
! in the direction of the current tout. Upon reaching tout (indicated
! by changing iflag to 2),the user must then define a new tout and
! reset iflag to -2 to continue in the one-step integrator mode.
!
! If the integration was not completed but the user still wants to
! continue (iflag=3,4 cases), he just calls rkf45 again. with iflag=3
! the relerr parameter has been adjusted appropriately for continuing
! the integration. in the case of iflag=4 the function counter will
! be reset to 0 and another 3000 function evaluations are allowed.
!
! However,in the case iflag=5, the user must first alter the error
! criterion to use a positive value of abserr before integration can
! proceed. if he does not,execution is terminated.
!
! Also,in the case iflag=6, it is necessary for the user to reset
! iflag to 2 (or -2 when the one-step integration mode is being used)
! as well as increasing either abserr,relerr or both before the
! integration can be continued. if this is not done, execution will
! be terminated. the occurrence of iflag=6 indicates a trouble spot
! (solution is changing rapidly,singularity may be present) and it
! often is inadvisable to continue.
!
! If iflag=7 is encountered, the user should use the one-step
! integration mode with the stepsize determined by the code.
! If the user insists upon continuing the integration with RKF45,
! he must reset iflag to 2 before calling RKF45 again. otherwise,
! execution will be terminated.
!
! If iflag=8 is obtained, integration can not be continued unless
! the invalid input parameters are corrected.
!
! The arrays work and iwork contain information
! required for subsequent integration, and should not be altered.
!
implicit none
!
integer neqn
!
real abserr
external f
integer iflag
integer iwork(5)
integer k1
integer k1m
integer k2
integer k3
integer k4
integer k5
integer k6
real relerr
real t
real tout
real work(6*neqn+3)
real y(neqn)
!
! Compute indices for the splitting of the work array
!
k1m = neqn + 1
k1 = k1m + 1
k2 = k1 + neqn
k3 = k2 + neqn
k4 = k3 + neqn
k5 = k4 + neqn
k6 = k5 + neqn
!
! This interfacing routine merely relieves the user of a long
! calling list via the splitting apart of two working storage arrays.
!
call rkfs ( f, neqn, y, t, tout, relerr, abserr, iflag, work(1), work(k1m), &
work(k1), work(k2), work(k3), work(k4), work(k5), work(k6), work(k6+1), &
iwork(1), iwork(2), iwork(3), iwork(4), iwork(5) )
return
end
subroutine rkf45_d ( f, neqn, y, t, tout, relerr, abserr, iflag, work, iwork )
!
!*******************************************************************************
!
!! RKF45_D carries out the Runge-Kutta-Fehlberg method (double precision).
!
!
! Author:
!
! H A Watts and L F Shampine,
! Sandia Laboratories,
! Albuquerque, New Mexico.
!
! RKF45 is primarily designed to solve non-stiff and mildly stiff
! differential equations when derivative evaluations are inexpensive.
! RKF45 should generally not be used when the user is demanding
! high accuracy.
!
! Abstract:
!
! RKF45 integrates a system of NEQN first order ordinary differential
! equations of the form:
!
! dY(i)/dT = F(T,Y(1),Y(2),...,Y(NEQN))
!
! where the Y(1:NEQN) are given at T.
!
! Typically the subroutine is used to integrate from T to TOUT but it
! can be used as a one-step integrator to advance the solution a
! single step in the direction of TOUT. On return, the parameters in
! the call list are set for continuing the integration. The user has
! only to call RKF45 again (and perhaps define a new value for TOUT).
! Actually, RKF45 is an interfacing routine which calls subroutine
! RKFS for the solution. RKFS in turn calls subroutine FEHL which
! computes an approximate solution over one step.
!
! Reference:
!
! E. Fehlberg,
! Low-order Classical Runge-Kutta Formulas with Stepsize Control,
! NASA Technical Report R-315.
!
! L F Shampine, H A Watts, S Davenport,
! Solving Non-stiff Ordinary Differential Equations - The State of the Art,
! Sandia Laboratories Report SAND75-0182,
! To appear in SIAM Review.
!
! Parameters:
!
! Input, external F, a subroutine of the form
! subroutine f ( t, y, yp )
! double precision t
! double precision y(*)
! double precision yp(*)
! to evaluate the derivatives YP(I) = dY(I) / dT
!
! Input, integer NEQN, the number of equations to be integrated.
!
! Input/output, double precision Y(NEQN), the solution vector at T.
!
! Input/output, double precision T, the independent variable.
!
! Input, double precision TOUT, the output point at which solution is
! desired.
!
! Input, double precision RELERR, ABSERR, the relative and absolute error
! tolerances for the local error test. At each step the code requires:
! abs(local error) <= relerr*abs(y) + abserr
! for each component of the local error and solution vectors
!
! Output, integer IFLAG, indicator for status of integration.
!
! Workspace, double precision WORK(3+6*NEQN), an array to hold information
! internal to RKF45 which is necessary for subsequent calls.
!
! Workspace, integer IWORK(5), an array used to hold information internal
! to RKF45 which is necessary for subsequent calls.
!
!
! first call
!
! The user must provide storage in his calling program for the arrays
! in the call list - y(neqn) , work(3+6*neqn) , iwork(5) ,
! declare f in an external statement, supply subroutine f(t,y,yp) and
! initialize the following parameters-
!
! neqn -- number of equations to be integrated. (neqn >= 1)
!
! y(*) -- vector of initial conditions
!
! t -- starting point of integration , must be a variable
!
! tout -- output point at which solution is desired.
! t=tout is allowed on the first call only, in which case
! rkf45 returns with iflag=2 if continuation is possible.
!
! relerr,abserr -- relative and absolute local error tolerances
! which must be non-negative. relerr must be a variable while
! abserr may be a constant. the code should normally not be
! used with relative error control smaller than about 1.e-8 .
! to avoid limiting precision difficulties the code requires
! relerr to be larger than an internally computed relative
! error parameter which is machine dependent. in particular,
! pure absolute error is not permitted. if a smaller than
! allowable value of relerr is attempted, rkf45 increases
! relerr appropriately and returns control to the user before
! continuing the integration.
!
! iflag -- +1,-1 indicator to initialize the code for each new
! problem. normal input is +1. the user should set iflag=-1
! only when one-step integrator control is essential. in this
! case, rkf45 attempts to advance the solution a single step
! in the direction of tout each time it is called. since this
! mode of operation results in extra computing overhead, it
! should be avoided unless needed.
!
!
! output from rkf45
!
! y(*) -- solution at t
! t -- last point reached in integration.
! iflag = 2 -- integration reached tout. indicates successful retur
! and is the normal mode for continuing integration.
! =-2 -- a single successful step in the direction of tout
! has been taken. normal mode for continuing
! integration one step at a time.
! = 3 -- integration was not completed because relative error
! tolerance was too small. relerr has been increased
! appropriately for continuing.
! = 4 -- integration was not completed because more than
! 3000 derivative evaluations were needed. this
! is approximately 500 steps.
! = 5 -- integration was not completed because solution
! vanished making a pure relative error test
! impossible. must use non-zero abserr to continue.
! using the one-step integration mode for one step
! is a good way to proceed.
! = 6 -- integration was not completed because requested
! accuracy could not be achieved using smallest
! allowable stepsize. user must increase the error
! tolerance before continued integration can be
! attempted.
! = 7 -- it is likely that rkf45 is inefficient for solving
! this problem. too much output is restricting the
! natural stepsize choice. use the one-step integrator
! mode.
! = 8 -- invalid input parameters
! this indicator occurs if any of the following is
! satisfied - neqn <= 0
! t=tout and iflag /= +1 or -1
! relerr or abserr < 0.
! iflag == 0 or < -2 or > 8
! work(*),iwork(*) -- information which is usually of no interest
! to the user but necessary for subsequent calls.
! work(1),...,work(neqn) contain the first derivatives
! of the solution vector y at t. work(neqn+1) contains
! the stepsize h to be attempted on the next step.
! iwork(1) contains the derivative evaluation counter.
!
!
! subsequent calls
!
! RKF45 returns with all information needed to continue
! the integration. if the integration reached tout, the user need onl
! define a new tout and call RKF45 again. In the one-step integrator
! mode (iflag=-2) the user must keep in mind that each step taken is
! in the direction of the current tout. Upon reaching tout (indicated
! by changing iflag to 2),the user must then define a new tout and
! reset iflag to -2 to continue in the one-step integrator mode.
!
! If the integration was not completed but the user still wants to
! continue (iflag=3,4 cases), he just calls rkf45 again. with iflag=3
! the relerr parameter has been adjusted appropriately for continuing
! the integration. in the case of iflag=4 the function counter will
! be reset to 0 and another 3000 function evaluations are allowed.
!
! However,in the case iflag=5, the user must first alter the error
! criterion to use a positive value of abserr before integration can
! proceed. if he does not,execution is terminated.
!
! Also,in the case iflag=6, it is necessary for the user to reset
! iflag to 2 (or -2 when the one-step integration mode is being used)
! as well as increasing either abserr,relerr or both before the
! integration can be continued. if this is not done, execution will
! be terminated. the occurrence of iflag=6 indicates a trouble spot
! (solution is changing rapidly,singularity may be present) and it
! often is inadvisable to continue.
!
! If iflag=7 is encountered, the user should use the one-step
! integration mode with the stepsize determined by the code.
! If the user insists upon continuing the integration with RKF45,
! he must reset iflag to 2 before calling RKF45 again. otherwise,
! execution will be terminated.
!
! If iflag=8 is obtained, integration can not be continued unless
! the invalid input parameters are corrected.
!
! The arrays work and iwork contain information
! required for subsequent integration, and should not be altered.
!
implicit none
!
integer neqn
!
double precision abserr
external f
integer iflag
integer iwork(5)
integer k1
integer k1m
integer k2
integer k3
integer k4
integer k5
integer k6
double precision relerr
double precision t
double precision tout
double precision work(6*neqn+3)
double precision y(neqn)
!
! Compute indices for the splitting of the work array
!
k1m = neqn + 1
k1 = k1m + 1
k2 = k1 + neqn
k3 = k2 + neqn
k4 = k3 + neqn
k5 = k4 + neqn
k6 = k5 + neqn
!
! This interfacing routine merely relieves the user of a long
! calling list via the splitting apart of two working storage arrays.
!
call rkfs_d ( f, neqn, y, t, tout, relerr, abserr, iflag, work(1), &
work(k1m), work(k1), work(k2), work(k3), work(k4), work(k5), work(k6), &
work(k6+1), iwork(1), iwork(2), iwork(3), iwork(4), iwork(5) )
return
end
subroutine rkfs ( f, neqn, y, t, tout, relerr, abserr, iflag, yp, h, &
f1, f2, f3, f4, f5, savre, savae, &
nfe, kop, init, jflag, kflag )
!
!*******************************************************************************
!
!! RKFS implements the Runge-Kutta-Fehlberg method (single precision).
!
!
! Discussion:
!
! RKFS integrates a system of first order ordinary differential
! equations as described in the comments for RKF45.
!
! The arrays yp, f1, f2, f3, f4, and f5 (of dimension at least neqn) and
! the variables h, savre, savae, nfe, kop, init, jflag and kflag are used
! internally by the code and appear in the call list to eliminate
! local retention of variables between calls. Accordingly, they
! should not be altered. Items of possible interest are
!
! YP - the derivative of the solution vector at T;
! H - an appropriate stepsize to be used for the next step;
! NFE - the number of derivative function evaluations.
!
! The expense is controlled by restricting the number
! of function evaluations to be approximately MAXNFE.
! As set, this corresponds to about 500 steps.
!
! REMIN is the minimum acceptable value of RELERR. Attempts
! to obtain higher accuracy with this subroutine are usually
! very expensive and often unsuccessful.
!
implicit none
!
integer neqn
!
real a
real abserr
real ae
real dt
real ee
real eeoet
real eps
real esttol
real et
external f
real f1(neqn)
real f2(neqn)
real f3(neqn)
real f4(neqn)
real f5(neqn)
real h
logical hfaild
real hmin
integer iflag
integer init
integer jflag
integer k
integer kflag
integer kop
integer, parameter :: maxnfe = 3000
integer mflag
integer nfe
logical output
real relerr
real, parameter :: remin = 1.0E-12
real rer
real s
real savae
real savre
real scale
real t
real tol
real toln
real tout
real y(neqn)
real yp(neqn)
real ypk
!
! Check the input parameters.
!
eps = epsilon ( eps )
if ( neqn < 1 ) then
iflag = 8
return
end if
if ( relerr < 0.0E+00 ) then
iflag = 8
return
end if
if ( abserr < 0.0E+00 ) then
iflag = 8
return
end if
mflag = abs ( iflag )
if ( abs ( iflag ) < 1 .or. abs ( iflag ) > 8 ) then
iflag = 8
return
end if
!
! Is this the first call?
!
if ( mflag == 1 ) then
go to 50
end if
!
! Check continuation possibilities
!
if ( t == tout .and. kflag /= 3 ) then
iflag = 8
return
end if
if ( mflag /= 2 ) then
go to 25
end if
!
! iflag = +2 or -2
!
if ( kflag == 3 ) go to 45
if ( init == 0 ) go to 45
if ( kflag == 4 ) go to 40
if ( kflag == 5 .and. abserr == 0.0E+00 ) then
stop
end if
if ( kflag == 6 .and. relerr <= savre .and. abserr <= savae ) then
stop
end if
go to 50
!
! iflag = 3,4,5,6,7 or 8
!
25 continue
if ( iflag == 3 ) go to 45
if ( iflag == 4 ) go to 40
if ( iflag == 5 .and. abserr > 0.0E+00 ) go to 45
!
! Integration cannot be continued since user did not respond to
! the instructions pertaining to iflag=5,6,7 or 8
!
stop
!
! Reset function evaluation counter
!
40 continue
nfe = 0
if ( mflag == 2 ) then
go to 50
end if
!
! Reset flag value from previous call
!
45 continue
iflag = jflag
if ( kflag == 3 ) then
mflag = abs ( iflag )
end if
!
! Save input iflag and set continuation flag for subsequent input checking.
!
50 continue
jflag = iflag
kflag = 0
!
! Save relerr and abserr for checking input on subsequent calls
!
savre = relerr
savae = abserr
!
! Restrict relative error tolerance to be at least as large as
! 2*eps+remin to avoid limiting precision difficulties arising
! from impossible accuracy requests
!
rer = 2.0E+00 * epsilon ( rer ) + remin
!
! The relative error tolerance is too small.
!
if ( relerr < rer ) then
relerr = rer
iflag = 3
kflag = 3
return
end if
dt = tout - t
if ( mflag == 1 ) go to 60
if ( init == 0 ) go to 65
go to 80
!
! Initialization:
! set initialization completion indicator, init
! set indicator for too many output points, kop
! evaluate initial derivatives
! set counter for function evaluations, nfe
! evaluate initial derivatives
! set counter for function evaluations, nfe
! estimate starting stepsize
!
60 continue
init = 0
kop = 0
a = t
call f ( a, y, yp )
nfe = 1
if ( t == tout ) then
iflag = 2
return
end if
65 continue
init = 1
h = abs ( dt )
toln = 0.0E+00
do k = 1, neqn
tol = relerr * abs ( y(k) ) + abserr
if ( tol > 0.0E+00 ) then
toln = tol
ypk = abs ( yp(k) )
if ( ypk * h**5 > tol) then
h = ( tol / ypk )**0.2E+00
end if
end if
end do
if ( toln <= 0.0E+00 ) then
h = 0.0E+00
end if
h = max ( h, 26.0E+00 * eps * max ( abs ( t ), abs ( dt ) ) )
jflag = sign ( 2, iflag )
!
! Set stepsize for integration in the direction from T to TOUT.
!
80 continue
h = sign ( h, dt )
!
! Test to see if RKF45 is being severely impacted by too many output points.
!
if ( abs ( h ) >= 2.0E+00 * abs ( dt ) ) then
kop = kop + 1
end if
!
! Unnecessary frequency of output.
!
if ( kop == 100 ) then
kop = 0
iflag = 7
return
end if
!
! If too close to output point, extrapolate and return.
!
if ( abs ( dt ) <= 26.0E+00 * eps * abs ( t ) ) then
y(1:neqn) = y(1:neqn) + dt * yp(1:neqn)
a = tout
call f ( a, y, yp )
nfe = nfe + 1
t = tout
iflag = 2
return
end if
!
! Initialize output point indicator.
!
output = .false.
!
! To avoid premature underflow in the error tolerance function,
! scale the error tolerances
!
scale = 2.0E+00 / relerr
ae = scale * abserr
!
! Step by step integration.
!
100 continue
hfaild = .false.
!
! Set smallest allowable stepsize.
!
hmin = 26.0E+00 * eps * abs ( t )
!
! Adjust stepsize if necessary to hit the output point.
! Look ahead two steps to avoid drastic changes in the stepsize and
! thus lessen the impact of output points on the code.
!
dt = tout - t
if ( abs ( dt ) >= 2.0E+00 * abs ( h ) ) go to 200
!
! The next successful step will complete the integration to the output point.
!
if ( abs ( dt ) <= abs ( h ) ) then
output = .true.
h = dt
go to 200
end if
h = 0.5E+00 * dt
!
! Core integrator for taking a single step
!
! The tolerances have been scaled to avoid premature underflow in
! computing the error tolerance function ET.
! To avoid problems with zero crossings, relative error is measured
! using the average of the magnitudes of the solution at the
! beginning and end of a step.
! The error estimate formula has been grouped to control loss of
! significance.
!
! To distinguish the various arguments, H is not permitted
! to become smaller than 26 units of roundoff in T.
! Practical limits on the change in the stepsize are enforced to
! smooth the stepsize selection process and to avoid excessive
! chattering on problems having discontinuities.
! To prevent unnecessary failures, the code uses 9/10 the stepsize
! it estimates will succeed.
!
! After a step failure, the stepsize is not allowed to increase for
! the next attempted step. This makes the code more efficient on
! problems having discontinuities and more effective in general
! since local extrapolation is being used and extra caution seems
! warranted.
!
! Test number of derivative function evaluations.
! If okay, try to advance the integration from T to T+H.
!
200 continue
!
! Too much work.
!
if ( nfe > maxnfe ) then
iflag = 4
kflag = 4
return
end if
!
! Advance an approximate solution over one step of length H.
!
call fehl ( f, neqn, y, t, h, yp, f1, f2, f3, f4, f5, f1 )
nfe = nfe + 5
!
! Compute and test allowable tolerances versus local error estimates
! and remove scaling of tolerances. Note that relative error is
! measured with respect to the average of the magnitudes of the
! solution at the beginning and end of the step.
!
eeoet = 0.0E+00
do k = 1, neqn
et = abs ( y(k) ) + abs ( f1(k) ) + ae
if ( et <= 0.0E+00 ) then
iflag = 5
return
end if
ee = abs ( ( -2090.0E+00 * yp(k) + ( 21970.0E+00 * f3(k) &
- 15048.0E+00 * f4(k) ) ) + ( 22528.0E+00 * f2(k) - 27360.0E+00 * f5(k) ) )
eeoet = max ( eeoet, ee / et )
end do
esttol = abs ( h ) * eeoet * scale / 752400.0E+00
if ( esttol <= 1.0E+00 ) then
go to 260
end if
!
! Unsuccessful step. Reduce the stepsize, try again.
! The decrease is limited to a factor of 1/10.
!
hfaild = .true.
output = .false.
if ( esttol < 59049.0E+00 ) then
s = 0.9E+00 / esttol**0.2E+00
else
s = 0.1E+00
end if
h = s * h
if ( abs ( h ) < hmin ) then
iflag = 6
kflag = 6
return
else
go to 200
end if
!
! Successful step. Store solution at T+H and evaluate derivatives there.
!
260 continue
t = t + h
y(1:neqn) = f1(1:neqn)
a = t
call f ( a, y, yp )
nfe = nfe + 1
!
! Choose next stepsize. The increase is limited to a factor of 5.
! If step failure has just occurred, next stepsize is not allowed to increase
!
if ( esttol > 0.0001889568E+00 ) then
s = 0.9E+00 / esttol**0.2E+00
else
s = 5.0E+00
end if
if ( hfaild ) then
s = min ( s, 1.0E+00 )
end if
h = sign ( max ( s * abs ( h ), hmin ), h )
!
! End of core integrator
!
! Should we take another step?
!
if ( output ) then
t = tout
iflag = 2
end if
if ( iflag > 0 ) go to 100
!
! Integration successfully completed
!
! one-step mode
!
iflag = - 2
return
end
subroutine rkfs_d ( f, neqn, y, t, tout, relerr, abserr, iflag, yp, h, &
f1, f2, f3, f4, f5, savre, savae, nfe, kop, init, jflag, kflag )
!
!*******************************************************************************
!
!! RKFS_D implements the Runge-Kutta-Fehlberg method (double precision).
!
!
! Discussion:
!
! RKFS integrates a system of first order ordinary differential
! equations as described in the comments for RKF45.
!
! The arrays yp, f1, f2, f3, f4, and f5 (of dimension at least neqn) and
! the variables h, savre, savae, nfe, kop, init, jflag and kflag are used
! internally by the code and appear in the call list to eliminate
! local retention of variables between calls. Accordingly, they
! should not be altered. Items of possible interest are
!
! YP - the derivative of the solution vector at T;
! H - an appropriate stepsize to be used for the next step;
! NFE - the number of derivative function evaluations.
!
! The expense is controlled by restricting the number
! of function evaluations to be approximately MAXNFE.
! As set, this corresponds to about 500 steps.
!
! REMIN is the minimum acceptable value of RELERR. Attempts
! to obtain higher accuracy with this subroutine are usually
! very expensive and often unsuccessful.
!
implicit none
!
integer neqn
!
double precision a
double precision abserr
double precision ae
double precision dt
double precision ee
double precision eeoet
double precision eps
double precision esttol
double precision et
external f
double precision f1(neqn)
double precision f2(neqn)
double precision f3(neqn)
double precision f4(neqn)
double precision f5(neqn)
double precision h
logical hfaild
double precision hmin
integer iflag
integer init
integer jflag
integer k
integer kflag
integer kop
integer, parameter :: maxnfe = 3000
integer mflag
integer nfe
logical output
double precision relerr
double precision, parameter :: remin = 1.0D-12
double precision rer
double precision s
double precision savae
double precision savre
double precision scale
double precision t
double precision tol
double precision toln
double precision tout
double precision y(neqn)
double precision yp(neqn)
double precision ypk
!
! Check the input parameters.
!
eps = epsilon ( eps )
if ( neqn < 1 ) then
iflag = 8
return
end if
if ( relerr < 0.0D+00 ) then
iflag = 8
return
end if
if ( abserr < 0.0D+00 ) then
iflag = 8
return
end if
mflag = abs ( iflag )
if ( abs ( iflag ) < 1 .or. abs ( iflag ) > 8 ) then
iflag = 8
return
end if
!
! Is this the first call?
!
if ( mflag == 1 ) then
go to 50
end if
!
! Check continuation possibilities
!
if ( t == tout .and. kflag /= 3 ) then
iflag = 8
return
end if
if ( mflag /= 2 ) then
go to 25
end if
!
! iflag = +2 or -2
!
if ( kflag == 3 ) go to 45
if ( init == 0 ) go to 45
if ( kflag == 4 ) go to 40
if ( kflag == 5 .and. abserr == 0.0D+00 ) then
stop
end if
if ( kflag == 6 .and. relerr <= savre .and. abserr <= savae ) then
stop
end if
go to 50
!
! iflag = 3,4,5,6,7 or 8
!
25 continue
if ( iflag == 3 ) go to 45
if ( iflag == 4 ) go to 40
if ( iflag == 5 .and. abserr > 0.0D+00 ) go to 45
!
! Integration cannot be continued since user did not respond to
! the instructions pertaining to iflag=5,6,7 or 8
!
stop
!
! Reset function evaluation counter
!
40 continue
nfe = 0
if ( mflag == 2 ) then
go to 50
end if
!
! Reset flag value from previous call
!
45 continue
iflag = jflag
if ( kflag == 3 ) then
mflag = abs ( iflag )
end if
!
! Save input iflag and set continuation flag for subsequent input checking.
!
50 continue
jflag = iflag
kflag = 0
!
! Save relerr and abserr for checking input on subsequent calls
!
savre = relerr
savae = abserr
!
! Restrict relative error tolerance to be at least as large as
! 2*eps+remin to avoid limiting precision difficulties arising
! from impossible accuracy requests
!
rer = 2.0D+00 * epsilon ( rer ) + remin
!
! The relative error tolerance is too small.
!
if ( relerr < rer ) then
relerr = rer
iflag = 3
kflag = 3
return
end if
dt = tout - t
if ( mflag == 1 ) go to 60
if ( init == 0 ) go to 65
go to 80
!
! Initialization:
! set initialization completion indicator, init
! set indicator for too many output points, kop
! evaluate initial derivatives
! set counter for function evaluations, nfe
! evaluate initial derivatives
! set counter for function evaluations, nfe
! estimate starting stepsize
!
60 continue
init = 0
kop = 0
a = t
call f ( a, y, yp )
nfe = 1
if ( t == tout ) then
iflag = 2
return
end if
65 continue
init = 1
h = abs ( dt )
toln = 0.0D+00
do k = 1, neqn
tol = relerr * abs ( y(k) ) + abserr
if ( tol > 0.0D+00 ) then
toln = tol
ypk = abs ( yp(k) )
if ( ypk * h**5 > tol) then
h = ( tol / ypk )**0.2D+00
end if
end if
end do
if ( toln <= 0.0D+00 ) then
h = 0.0D+00
end if
h = max ( h, 26.0D+00 * eps * max ( abs ( t ), abs ( dt ) ) )
jflag = sign ( 2, iflag )
!
! Set stepsize for integration in the direction from T to TOUT.
!
80 continue
h = sign ( h, dt )
!
! Test to see if RKF45 is being severely impacted by too many output points.
!
if ( abs ( h ) >= 2.0D+00 * abs ( dt ) ) then
kop = kop + 1
end if
!
! Unnecessary frequency of output.
!
if ( kop == 100 ) then
kop = 0
iflag = 7
return
end if
!
! If too close to output point, extrapolate and return.
!
if ( abs ( dt ) <= 26.0D+00 * eps * abs ( t ) ) then
y(1:neqn) = y(1:neqn) + dt * yp(1:neqn)
a = tout
call f ( a, y, yp )
nfe = nfe + 1
t = tout
iflag = 2
return
end if
!
! Initialize output point indicator.
!
output = .false.
!
! To avoid premature underflow in the error tolerance function,
! scale the error tolerances
!
scale = 2.0D+00 / relerr
ae = scale * abserr
!
! Step by step integration.
!
100 continue
hfaild = .false.
!
! Set smallest allowable stepsize.
!
hmin = 26.0D+00 * eps * abs ( t )
!
! Adjust stepsize if necessary to hit the output point.
! Look ahead two steps to avoid drastic changes in the stepsize and
! thus lessen the impact of output points on the code.
!
dt = tout - t
if ( abs ( dt ) >= 2.0D+00 * abs ( h ) ) go to 200
!
! The next successful step will complete the integration to the output point.
!
if ( abs ( dt ) <= abs ( h ) ) then
output = .true.
h = dt
go to 200
end if
h = 0.5D+00 * dt
!
! Core integrator for taking a single step
!
! The tolerances have been scaled to avoid premature underflow in
! computing the error tolerance function ET.
! To avoid problems with zero crossings, relative error is measured
! using the average of the magnitudes of the solution at the
! beginning and end of a step.
! The error estimate formula has been grouped to control loss of
! significance.
!
! To distinguish the various arguments, H is not permitted
! to become smaller than 26 units of roundoff in T.
! Practical limits on the change in the stepsize are enforced to
! smooth the stepsize selection process and to avoid excessive
! chattering on problems having discontinuities.
! To prevent unnecessary failures, the code uses 9/10 the stepsize
! it estimates will succeed.
!
! After a step failure, the stepsize is not allowed to increase for
! the next attempted step. This makes the code more efficient on
! problems having discontinuities and more effective in general
! since local extrapolation is being used and extra caution seems
! warranted.
!
! Test number of derivative function evaluations.
! If okay, try to advance the integration from T to T+H.
!
200 continue
!
! Too much work.
!
if ( nfe > maxnfe ) then
iflag = 4
kflag = 4
return
end if
!
! Advance an approximate solution over one step of length H.
!
call fehl_d ( f, neqn, y, t, h, yp, f1, f2, f3, f4, f5, f1 )
nfe = nfe + 5
!
! Compute and test allowable tolerances versus local error estimates
! and remove scaling of tolerances. Note that relative error is
! measured with respect to the average of the magnitudes of the
! solution at the beginning and end of the step.
!
eeoet = 0.0D+00
do k = 1, neqn
et = abs ( y(k) ) + abs ( f1(k) ) + ae
if ( et <= 0.0D+00 ) then
iflag = 5
return
end if
ee = abs ( ( -2090.0D+00 * yp(k) + ( 21970.0D+00 * f3(k) &
- 15048.0D+00 * f4(k) ) ) &
+ ( 22528.0D+00 * f2(k) - 27360.0D+00 * f5(k) ) )
eeoet = max ( eeoet, ee / et )
end do
esttol = abs ( h ) * eeoet * scale / 752400.0D+00
if ( esttol <= 1.0D+00 ) then
go to 260
end if
!
! Unsuccessful step. Reduce the stepsize, try again.
! The decrease is limited to a factor of 1/10.
!
hfaild = .true.
output = .false.
if ( esttol < 59049.0D+00 ) then
s = 0.9D+00 / esttol**0.2D+00
else
s = 0.1D+00
end if
h = s * h
if ( abs ( h ) < hmin ) then
iflag = 6
kflag = 6
return
else
go to 200
end if
!
! Successful step. Store solution at T+H and evaluate derivatives there.
!
260 continue
t = t + h
y(1:neqn) = f1(1:neqn)
a = t
call f ( a, y, yp )
nfe = nfe + 1
!
! Choose next stepsize. The increase is limited to a factor of 5.
! If step failure has just occurred, next stepsize is not allowed to increase
!
if ( esttol > 0.0001889568D+00 ) then
s = 0.9D+00 / esttol**0.2D+00
else
s = 5.0D+00
end if
if ( hfaild ) then
s = min ( s, 1.0D+00 )
end if
h = sign ( max ( s * abs ( h ), hmin ), h )
!
! End of core integrator
!
! Should we take another step?
!
if ( output ) then
t = tout
iflag = 2
end if
if ( iflag > 0 ) go to 100
!
! Integration successfully completed
!
! one-step mode
!
iflag = - 2
return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
! Example:
!
! May 31 2001 9:45:54.872 AM
!
! Modified:
!
! 31 May 2001
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! None
!
implicit none
!
character ( len = 8 ) ampm
integer d
character ( len = 8 ) date
integer h
integer m
integer mm
character ( len = 9 ), parameter, dimension(12) :: month = (/ &
'January ', 'February ', 'March ', 'April ', &
'May ', 'June ', 'July ', 'August ', &
'September', 'October ', 'November ', 'December ' /)
integer n
integer s
character ( len = 10 ) time
integer values(8)
integer y
character ( len = 5 ) zone
!
call date_and_time ( date, time, zone, values )
y = values(1)
m = values(2)
d = values(3)
h = values(5)
n = values(6)
s = values(7)
mm = values(8)
if ( h < 12 ) then
ampm = 'AM'
else if ( h == 12 ) then
if ( n == 0 .and. s == 0 ) then
ampm = 'Noon'
else
ampm = 'PM'
end if
else
h = h - 12
if ( h < 12 ) then
ampm = 'PM'
else if ( h == 12 ) then
if ( n == 0 .and. s == 0 ) then
ampm = 'Midnight'
else
ampm = 'AM'
end if
end if
end if
write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )
return
end
!end of file rkf45.f90