Status
Description
This is a modern object-oriented Fortran implementation of the DDEABM Adams-Bashforth-Moulton ODE solver. The original Fortran 77 code was obtained from the SLATEC library. It has been extensively refactored.
DDEABM uses the Adams-Bashforth-Moulton predictor-corrector formulas of orders 1 through 12 to integrate a system of first order ordinary differential equations of the form dx/dt = f(t,x). Also included is an event-location capability, where the equations can be integrated until a specified function g(t,x) = 0. Dense output is also supported.
This project is hosted on GitHub.
Examples
The ddeabm_module provides a thread-safe and object-oriented interface to the DDEABM method. Some example use cases are presented below:
Basic integration
This example shows how to integrate a conic orbit (6 state equations) around the Earth from an initial time t0 to a final time tf:
program ddeabm_example
use ddeabm_module, wp => ddeabm_rk
implicit none
real(wp),parameter :: mu = 398600.436233_wp !! Earth gravitational parameter (km^3/s^2)
integer,parameter :: n = 6 !! number of state variables
type(ddeabm_class) :: s
real(wp),dimension(n) :: x0,x
real(wp) :: t0,tf,t
integer :: idid
call s%initialize(n,maxnum=10000,df=twobody,rtol=[1.0e-12_wp],atol=[1.0e-12_wp])
!initial conditions:
x0 = [10000.0_wp,10000.0_wp,10000.0_wp,& !initial state [r,v] (km,km/s)
1.0_wp,2.0_wp,3.0_wp]
t0 = 0.0_wp !initial time (sec)
tf = 1000.0_wp !final time (sec)
write(*,'(A/,*(F15.6/))') 'Initial time: ',t0
write(*,'(A/,*(F15.6/))') 'Initial state:',x0
t = t0
x = x0
call s%integrate(t,x,tf,idid=idid)
write(*,'(A/,*(F15.6/))') 'Final time: ',t
write(*,'(A/,*(F15.6/))') 'Final state:',x
contains
subroutine twobody(me,t,x,xdot)
!! derivative routine for two-body orbit propagation
implicit none
class(ddeabm_class),intent(inout) :: me
real(wp),intent(in) :: t
real(wp),dimension(:),intent(in) :: x
real(wp),dimension(:),intent(out) :: xdot
real(wp),dimension(3) :: r,v,a_grav
real(wp) :: rmag
r = x(1:3)
v = x(4:6)
rmag = norm2(r)
a_grav = -mu/rmag**3 * r ! acceleration due to gravity
xdot(1:3) = v
xdot(4:6) = a_grav
end subroutine twobody
end program ddeabm_example
It produces the following output:
Initial time:
0.000000
Initial state:
10000.000000
10000.000000
10000.000000
1.000000
2.000000
3.000000
Final time:
1000.000000
Initial time:
10667.963305
11658.055962
12648.148619
0.377639
1.350074
2.322509
The intermediate integration points can also be reported to a user-defined procedure. For the above example, the following subroutine could be defined:
subroutine twobody_report(me,t,x)
!! report function - write time,state to console
implicit none
class(ddeabm_class),intent(inout) :: me
real(wp),intent(in) :: t
real(wp),dimension(:),intent(in) :: x
write(*,'(*(F15.6,1X))') t,x
end subroutine twobody_report
Which can be added to the class on initialization:
call s%initialize(n,maxnum=10000,df=twobody,&
rtol=[1.0e-12_wp],atol=[1.0e-12_wp],&
report=twobody_report)
This function is then called at each time step
if the equations are integrated using the integration_mode=2 option like so:
call s%integrate(t,x,tf,idid=idid,integration_mode=2)
Event location
A user-defined event function g(t,x) can also be defined in order to stop the integration at a specified event (i.e., when g(t,x)=0). In the above example, say it is desired that the integration stop when z = x(3) = 12,000 km. The event function for this would be:
subroutine twobody_event(me,t,x,g)
!! event function for z = 12,000 km
implicit none
class(ddeabm_with_event_class),intent(inout) :: me
real(wp),intent(in) :: t
real(wp),dimension(:),intent(in) :: x
real(wp),intent(out) :: g
g = 12000.0_wp - x(3)
end subroutine twobody_event
For event finding, the ddeabm_with_event_class type is used (which is an extension of the main ddeabm_class). For example:
type(ddeabm_with_event_class) :: s
...
call s%initialize_event(n,maxnum=10000,df=twobody,&
rtol=[1.0e-12_wp],atol=[1.0e-12_wp],&
g=twobody_event,root_tol=1.0e-12_wp)
...
call s%integrate_to_event(t,x,tf,idid=idid,gval=gval)
In this case, root_tol is the tolerance for the event location, and gval is the value of the event function at the final time (note that the integration will stop when g(t,x)=0 or at t=tf, whichever occurs first).
A vector event function is also supported (in which case, the integration stops if any of the roots are found). This is done using the ddeabm_with_event_class_vec type.
Fixed time step
All of the integration methods have an optional argument (tstep) to enable a fixed time step, which can be used for dense output, or to specify a fixed step used for event finding (since the default step may be too large). For example, for performing a root-finding integration with the event function evaluated every 100 seconds:
call s%integrate_to_event(t,x,tf,idid=idid,gval=gval,tstep=100.0_wp)
The event function is also evaluated at the end of each integrator step. The points within a step are computed by interpolation, so tstep does not change the steps taken by the integrator.
Maximum step size
By default, the size of the integration steps is not limited. An optional max_step_size argument can be passed to initialize (or initialize_event) to limit the magnitude of every step. For example, to limit the step size to 60 seconds:
call s%initialize(n,maxnum=10000,df=twobody,&
rtol=[1.0e-12_wp],atol=[1.0e-12_wp],&
max_step_size=60.0_wp)
This can be useful when the solution has short-lived features that a large step could pass over (for example, two event roots closer together than the step size). Unlike tstep, it changes the steps taken by the integrator, so it can increase the number of derivative evaluations.
Maximum order
DDEABM is a variable-order method. By default, the maximum order is 12. A different maximum order (from 1 to ddeabm_max_order_limit, which is 20) can be specified using the optional max_order argument to initialize (or initialize_event):
call s%initialize(n,maxnum=10000,df=twobody,&
rtol=[1.0e-25_wp],atol=[1.0e-25_wp],&
max_order=16)
A lower maximum order can improve stability on some problems. A higher maximum order is mainly useful for very tight tolerances when the library is compiled with REAL128, where it can significantly reduce the number of steps. For the default REAL64, a maximum order above 12 makes little difference.
Continuing after an interrupted integration
The integration is interrupted (and can be continued) when the code returns one of these values of idid:
-1: the maximum number of steps (maxnum) has been taken.-2: the tolerances are too small, and have been increased by the code.-3: a solution component is zero and the corresponding absolute tolerance is zero.-4: the problem appears to be stiff.
To continue, call resume_integration and then call the integration routine again. For example:
call s%integrate(t,x,tf,idid=idid)
if (idid == -1) then
! allow another maxnum steps
call s%resume_integration()
call s%integrate(t,x,tf,idid=idid)
end if
For idid=-2, the increased tolerances are used for the rest of the problem. For idid=-3, the problem must first be corrected (by calling initialize again with a positive atol, which also restarts the integration).
Building DDEABM
DDEABM and the test programs will build with any modern Fortran compiler. A Fortran Package Manager (FPM) manifest file (fpm.toml) is included, so that the library and tests cases can be compiled with FPM. For example:
fpm build --profile release
fpm test --profile release
To generate the documentation using ford, run:
ford ford.md
To use ddeabm within your fpm project, add the following to your fpm.toml file:
[dependencies]
ddeabm = { git="https://github.com/jacobwilliams/ddeabm.git" }
A specific version can also be specified:
[dependencies]
ddeabm = { git="https://github.com/jacobwilliams/ddeabm.git", rev = "2.1.0" }
By default, the library is built with double precision (real64) real values. Explicitly specifying the real kind can be done using the following processor flags:
| Preprocessor flag |
Kind |
Number of bytes |
REAL32 |
real(kind=real32) |
4 |
REAL64 |
real(kind=real64) |
8 |
REAL128 |
real(kind=real128) |
16 |
For example, to build a single precision version of the library, use:
fpm build --profile release --flag "-DREAL32"
Dependencies
Building the library requires the roots-fortran module. Building the tests requires the pyplot-fortran module. FPM will automatically download the correct versions of both (see fpm.toml).
Documentation
The latest API documentation can be found here. This was generated from the source code using FORD.
License
The DDEABM source code and related files and documentation are distributed under a permissive free software license (BSD-style). The original DDEABM Fortran 77 code is public domain.
Keywords
Adams-Bashforth-Moulton Method, DEPAC, Initial Value Problems, ODE, Ordinary Differential Equations, Predictor-Corrector, SLATEC, Modern Fortran
References
- L. F. Shampine, M. K. Gordon, “Solving ordinary differential equations with ODE, STEP, and INTRP”, Report SLA-73-1060, Sandia Laboratories, 1973.
- L. F. Shampine, M. K. Gordon, “Computer solution of ordinary differential equations, the initial value problem”, W. H. Freeman and Company, 1975.
- L. F. Shampine, H. A. Watts, “DEPAC - Design of a user oriented package of ode solvers”, Report SAND79-2374, Sandia Laboratories, 1979.
- H. A. Watts, “A smoother interpolant for DE/STEP, INTRP and DEABM: II”, Report SAND84-0293, Sandia Laboratories, 1984.
- R. P. Brent, “An algorithm with guaranteed convergence for finding a zero of a function”, The Computer Journal, Vol 14, No. 4., 1971.
- R. P. Brent, “Algorithms for minimization without derivatives”, Prentice-Hall, Inc., 1973.
评论