subroutine fehl ( f, neqn, y, t, h, yp, f1, f2, f3, f4, f5, s )
!*****************************************************************************80
!
!! FEHL takes one Fehlberg fourth-fifth order step (double precision).
!
! Discussion:
!
! This routine integrates a system of NEQN first order ordinary differential
! equations of the form
! dY(i)/dT = F(T,Y(1:NEQN))
! where the initial values Y and the initial derivatives
! YP are specified at the starting point T.
!
! The routine 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.
! The routine should be called with an H not smaller than 13 units of
! roundoff in T so that the various independent arguments can be
! distinguished.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 27 March 2004
!
! Author:
!
! Original FORTRAN77 version by Herman Watts, Lawrence Shampine.
! FORTRAN90 version by John Burkardt.
!
! Reference:
!
! Erwin Fehlberg,
! Low-order Classical Runge-Kutta Formulas with Stepsize Control,
! NASA Technical Report R-315, 1969.
!
! Lawrence Shampine, Herman Watts, S Davenport,
! Solving Non-stiff Ordinary Differential Equations - The State of the Art,
! SIAM Review,
! Volume 18, pages 376-411, 1976.
!
! Parameters:
!
! Input, external F, a user-supplied subroutine to evaluate the
! derivatives Y'(T), of the form:
!
! subroutine f ( t, y, yp )
! real ( kind = rk ) t
! real ( kind = rk ) y(*)
! real ( kind = rk ) yp(*)
!
! Input, integer NEQN, the number of equations to be integrated.
!
! Input, real ( kind = rk ) Y(NEQN), the current value of the
! dependent variable.
!
! Input, real ( kind = rk ) T, the current value of the independent
! variable.
!
! Input, real ( kind = rk ) H, the step size to take.
!
! Input, real ( kind = rk ) YP(NEQN), the current value of the
! derivative of the dependent variable.
!
! Output, real ( kind = rk ) F1(NEQN), F2(NEQN), F3(NEQN), F4(NEQN),
! F5(NEQN), derivative values needed for the computation.
!
! Output, real ( kind = rk ) S(NEQN), the estimate of the solution at T+H.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer neqn
real ( kind = rk ) ch
external f
real ( kind = rk ) f1(neqn)
real ( kind = rk ) f2(neqn)
real ( kind = rk ) f3(neqn)
real ( kind = rk ) f4(neqn)
real ( kind = rk ) f5(neqn)
real ( kind = rk ) h
real ( kind = rk ) s(neqn)
real ( kind = rk ) t
real ( kind = rk ) y(neqn)
real ( kind = rk ) 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, yp, t, tout, relerr, abserr, flag )
!*****************************************************************************80
!
!! RKF45 carries out the Runge-Kutta-Fehlberg method (double precision).
!
! Discussion:
!
! This routine is primarily designed to solve non-stiff and mildly stiff
! differential equations when derivative evaluations are inexpensive.
! It should generally not be used when the user is demanding
! high accuracy.
!
! This routine 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 again (and perhaps define a new value for TOUT).
!
! Before the first call, the user must
!
! * supply the subroutine F(T,Y,YP) to evaluate the right hand side;
! and declare F in an EXTERNAL statement;
!
! * initialize the parameters:
! NEQN, Y(1:NEQN), T, TOUT, RELERR, ABSERR, FLAG.
! In particular, T should initially be the starting point for integration,
! Y should be the value of the initial conditions, and FLAG should
! normally be +1.
!
! Normally, the user only sets the value of FLAG before the first call, and
! thereafter, the program manages the value. On the first call, FLAG should
! normally be +1 (or -1 for single step mode.) On normal return, FLAG will
! have been reset by the program to the value of 2 (or -2 in single
! step mode), and the user can continue to call the routine with that
! value of FLAG.
!
! (When the input magnitude of FLAG is 1, this indicates to the program
! that it is necessary to do some initialization work. An input magnitude
! of 2 lets the program know that that initialization can be skipped,
! and that useful information was computed earlier.)
!
! The routine returns with all the information needed to continue
! the integration. If the integration reached TOUT, the user need only
! define a new TOUT and call again. In the one-step integrator
! mode, returning with FLAG = -2, the user must keep in mind that
! each step taken is in the direction of the current TOUT. Upon
! reaching TOUT, indicated by the output value of FLAG switching to 2,
! the user must define a new TOUT and reset FLAG to -2 to continue
! in the one-step integrator mode.
!
! In some cases, an error or difficulty occurs during a call. In that case,
! the output value of FLAG is used to indicate that there is a problem
! that the user must address. These values include:
!
! * 3, integration was not completed because the input value of RELERR, the
! relative error tolerance, was too small. RELERR has been increased
! appropriately for continuing. If the user accepts the output value of
! RELERR, then simply reset FLAG to 2 and continue.
!
! * 4, integration was not completed because more than MAXNFE derivative
! evaluations were needed. This is approximately (MAXNFE/6) steps.
! The user may continue by simply calling again. The function counter
! will be reset to 0, and another MAXNFE function evaluations are allowed.
!
! * 5, integration was not completed because the solution vanished,
! making a pure relative error test impossible. The user must use
! a 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 the requested accuracy
! could not be achieved, even using the smallest allowable stepsize.
! The user must increase the error tolerances ABSERR or RELERR before
! continuing. It is also necessary to reset FLAG to 2 (or -2 when
! the one-step integration mode is being used). The occurrence of
! FLAG = 6 indicates a trouble spot. The solution is changing
! rapidly, or a singularity may be present. It often is inadvisable
! to continue.
!
! * 7, it is likely that this routine is inefficient for solving
! this problem. Too much output is restricting the natural stepsize
! choice. The user should use the one-step integration mode with
! the stepsize determined by the code. If the user insists upon
! continuing the integration, reset FLAG to 2 before calling
! again. Otherwise, execution will be terminated.
!
! * 8, invalid input parameters, indicates one of the following:
! NEQN <= 0;
! T = TOUT and |FLAG| /= 1;
! RELERR < 0 or ABSERR < 0;
! FLAG == 0 or FLAG < -2 or 8 < FLAG.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 April 2011
!
! Author:
!
! Original FORTRAN77 version by Herman Watts, Lawrence Shampine.
! FORTRAN90 version by John Burkardt.
!
! Reference:
!
! Erwin Fehlberg,
! Low-order Classical Runge-Kutta Formulas with Stepsize Control,
! NASA Technical Report R-315, 1969.
!
! Lawrence Shampine, Herman Watts, S Davenport,
! Solving Non-stiff Ordinary Differential Equations - The State of the Art,
! SIAM Review,
! Volume 18, pages 376-411, 1976.
!
! Parameters:
!
! Input, external F, a user-supplied subroutine to evaluate the
! derivatives Y'(T), of the form:
!
! subroutine f ( t, y, yp )
! real ( kind = rk ) t
! real ( kind = rk ) y(*)
! real ( kind = rk ) yp(*)
!
! Input, integer NEQN, the number of equations to be integrated.
!
! Input/output, real ( kind = rk ) Y(NEQN), the current solution vector at T.
!
! Input/output, real ( kind = rk ) YP(NEQN), the current value of the
! derivative of the dependent variable. The user should not set or alter
! this information!
!
! Input/output, real ( kind = rk ) T, the current value of the independent
! variable.
!
! Input, real ( kind = rk ) TOUT, the output point at which solution is
! desired. TOUT = T is allowed on the first call only, in which case
! the routine returns with FLAG = 2 if continuation is possible.
!
! Input, real ( kind = rk ) 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 the solution vector Y.
! RELERR cannot be "too small". If the routine believes RELERR has been
! set too small, it will reset RELERR to an acceptable value and return
! immediately for user action.
!
! Input/output, integer FLAG, indicator for status of
! integration. On the first call, set FLAG to +1 for normal use, or to -1
! for single step mode. On return, a value of 2 or -2 indicates normal
! progress, while any other value indicates a problem that should
! be addressed.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer neqn
real ( kind = rk ) abserr
real ( kind = rk ), save :: abserr_save = -1.0D+00
real ( kind = rk ) ae
real ( kind = rk ) dt
real ( kind = rk ) ee
real ( kind = rk ) eeoet
real ( kind = rk ) eps
real ( kind = rk ) esttol
real ( kind = rk ) et
external f
real ( kind = rk ) f1(neqn)
real ( kind = rk ) f2(neqn)
real ( kind = rk ) f3(neqn)
real ( kind = rk ) f4(neqn)
real ( kind = rk ) f5(neqn)
real ( kind = rk ), save :: h = -1.0D+00
logical hfaild
real ( kind = rk ) hmin
integer flag
integer, save :: flag_save = -1000
integer, save :: init = -1000
integer k
integer, save :: kflag = -1000
integer, save :: kop = -1
integer, parameter :: maxnfe = 3000
integer mflag
integer, save :: nfe = -1
logical output
real ( kind = rk ) relerr
real ( kind = rk ) relerr_min
real ( kind = rk ), save :: relerr_save = -1.0D+00
real ( kind = rk ), parameter :: remin = 1.0E-12
real ( kind = rk ) s
real ( kind = rk ) scale
real ( kind = rk ) t
real ( kind = rk ) tol
real ( kind = rk ) toln
real ( kind = rk ) tout
real ( kind = rk ) y(neqn)
real ( kind = rk ) yp(neqn)
real ( kind = rk ) ypk
!
! Check the input parameters.
!
eps = epsilon ( eps )
if ( neqn < 1 ) then
flag = 8
return
end if
if ( relerr < 0.0D+00 ) then
flag = 8
return
end if
if ( abserr < 0.0D+00 ) then
flag = 8
return
end if
if ( flag == 0 .or. 8 < flag .or. flag < -2 ) then
flag = 8
return
end if
mflag = abs ( flag )
!
! Is this a continuation call?
!
if ( mflag /= 1 ) then
if ( t == tout .and. kflag /= 3 ) then
flag = 8
return
end if
if ( mflag == 2 ) then
if ( kflag == 3 ) then
flag = flag_save
mflag = abs ( flag )
else if ( init == 0 ) then
flag = flag_save
else if ( kflag == 4 ) then
nfe = 0
else if ( kflag == 5 .and. abserr == 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'RKF45 - Fatal error!'
write ( *, '(a)' ) ' KFLAG = 5 and ABSERR = 0.0'
stop
else if ( &
kflag == 6 .and. &
relerr <= relerr_save .and. &
abserr <= abserr_save ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'RKF45 - Fatal error!'
write ( *, '(a)' ) ' KFLAG = 6 and'
write ( *, '(a)' ) ' RELERR <= RELERR_SAVE and'
write ( *, '(a)' ) ' ABSERR <= ABSERR_SAVE'
stop
end if
!
! FLAG = 3, 4, 5, 6, 7 or 8.
!
else
if ( flag == 3 ) then
flag = flag_save
if ( kflag == 3 ) then
mflag = abs ( flag )
end if
else if ( flag == 4 ) then
nfe = 0
flag = flag_save
if ( kflag == 3 ) then
mflag = abs ( flag )
end if
else if ( flag == 5 .and. 0.0D+00 < abserr ) then
flag = flag_save
if ( kflag == 3 ) then
mflag = abs ( flag )
end if
!
! Integration cannot be continued because the user did not respond to
! the instructions pertaining to FLAG = 5, 6, 7 or 8.
!
else
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'RKF45 - Fatal error!'
write ( *, '(a)' ) ' Integration cannot be continued.'
write ( *, '(a)' ) ' The user did not respond to the output'
write ( *, '(a)' ) ' value FLAG = 5, 6, 7, or 8.'
stop
end if
end if
end if
!
! Save the input value of FLAG.
! Set the continuation flag KFLAG for subsequent input checking.
!
flag_save = flag
kflag = 0
!
! Save RELERR and ABSERR for checking input on subsequent calls.
!
relerr_save = relerr
abserr_save = abserr
!
! Restrict the relative error tolerance to be at least
!
! 2 * EPS + REMIN
!
! to avoid limiting precision difficulties arising from impossible
! accuracy requests.
!
relerr_min = 2.0D+00 * epsilon ( relerr_min ) + remin
!
! Is the relative error tolerance too small?
!
if ( relerr < relerr_min ) then
relerr = relerr_min
flag = 3
kflag = 3
return
end if
dt = tout - t
!
! Initialization:
!
! Set the initialization completion indicator, INIT;
! set the indicator for too many output points, KOP;
! evaluate the initial derivatives;
! set the counter for function evaluations, NFE;
! estimate the starting stepsize.
!
if ( mflag == 1 ) then
init = 0
kop = 0
call f ( t, y, yp )
nfe = 1
if ( t == tout ) then
flag = 2
return
end if
end if
if ( init == 0 ) then
init = 1
h = abs ( dt )
toln = 0.0D+00
do k = 1, neqn
tol = relerr * abs ( y(k) ) + abserr
if ( 0.0D+00 < tol ) then
toln = tol
ypk = abs ( yp(k) )
if ( tol < ypk * h**5 ) 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 ) ) )
flag_save = sign ( 2, flag )
end if
!
! Set the stepsize for integration in the direction from T to TOUT.
!
h = sign ( h, dt )
!
! Test to see if too may output points are being requested.
!
if ( 2.0D+00 * abs ( dt ) <= abs ( h ) ) then
kop = kop + 1
end if
!
! Unnecessary frequency of output.
! Seems like an error I'm willing to tolerate!
!
! if ( kop == 100 ) then
if ( kop == 10000 ) then
kop = 0
flag = 7
return
end if
!
! If we are too close to the output point, then simply extrapolate and return.
!
if ( abs ( dt ) <= 26.0D+00 * eps * abs ( t ) ) then
t = tout
y(1:neqn) = y(1:neqn) + dt * yp(1:neqn)
call f ( t, y, yp )
nfe = nfe + 1
flag = 2
return
end if
!
! Initialize the 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.
!
do
hfaild = .false.
!
! Set the smallest allowable stepsize.
!
hmin = 26.0D+00 * eps * abs ( t )
!
! Adjust the 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 ( 2.0D+00 * abs ( h ) <= abs ( dt ) ) then
else
!
! Will the next successful step complete the integration to the output point?
!
if ( abs ( dt ) <= abs ( h ) ) then
output = .true.
h = dt
else
h = 0.5D+00 * dt
end if
end if
!
! Here begins the 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 the number of derivative function evaluations.
! If okay, try to advance the integration from T to T+H.
!
do
!
! Have we done too much work?
!
if ( maxnfe < nfe ) then
flag = 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. The 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
flag = 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
exit
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
flag = 6
kflag = 6
return
end if
end do
!
! We exited the loop because we took a successful step.
! Store the solution for T+H, and evaluate the derivative there.
!
t = t + h
y(1:neqn) = f1(1:neqn)
call f ( t, y, yp )
nfe = nfe + 1
!
! Choose the next stepsize. The increase is limited to a factor of 5.
! If the step failed, the next stepsize is not allowed to increase.
!
if ( 0.0001889568D+00 < esttol ) 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
flag = 2
return
end if
if ( flag <= 0 ) then
exit
end if
end do
!
! One step integration mode.
!
flag = -2
return
end
subroutine timestamp ( )
!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
! Example:
!
! 31 May 2001 9:45:54.872 AM
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 18 May 2013
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! None
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
character ( len = 8 ) ampm
integer d
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
integer values(8)
integer y
call date_and_time ( values = 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 ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )
return
end