User-friendly Python interface to ODE solvers
Many Fortran and C libraries to solve ODEs are cumbersome to use due to numerous inputs and outputs. We developed a Python package that wraps existing ODE solvers in a consistent, easy-to-use class interface. Options are given sensible defaults, which may be overridden. By establishing a common interface, a user can easily try different solvers without having to learn the intricacies of each solver.
The ODE function can be written in Python, C, or Fortran. The solver directly calls the Fortran/C function, eliminating the Python overhead during integration. This yields the high performance required for large 2D and 3D problems.
Download
Source code is available under the GPL:
ode-0.1a4.tar.gz (242K)
This is a fairly early release, so there are undoubtedly bugs and rough edges. Feedback to <mrgates2 @ uiuc.edu> is welcome.
Presentations
- User-friendly Python interface to ODE solvers. Poster from my LLNL summer internship, 2005.
- Two-for-one deal: Explicit projective methods for stiff ODEs and User-friendly Python interface to ODE solvers. Presentation I gave at UIUC computer science seminar, 2005.
Current solvers
- LSODE — Livermore Solver for ODEs
- Adams and BDF methods with stiffness detection.
- RKC — Runge-Kutta-Chebyshev method
- Stabilized explicit method for use with mildly stiff problems.
- Forward Euler
- Backward Euler
- Trapezoid
- Fourth-order Runge-Kutta
- Backward Euler
- Suite of common integrators, implemented in Python.
Interface comparison
Here is a comparison of the Python and Fortran code to drive the LSODE integrator, with corresponding sections shaded in color. Note the simplicity and straight-forwardness of the Python implementation.
| Python | Fortran | |
|---|---|---|
| Interface | Standardized | Solver-specific |
| Options | Defaults, can override | Specify numerous options |
| Memory | Automatically managed | Size based on problem and method |
| Errors | Raises exceptions | Explicitly check |
| Sample code |
from ode.lsode import lsode
import myode
t0 = 0.
tf = 0.7
y0 = myode.initial_condition( t0 )
try:
solver = lsode( myode.f, t0, y0 )
solver.reltol = 1e-3
solver.set_jacobian_bandwidth( 19*19, 19*19 )
y = solver.integrate( tf )
print "Steps", solver.steps
print "F evals", solver.function_evals
except Exception, e:
print e
|
integer neq, lrw, liw parameter (neq=6859, liw=6879, lrw=7496909) integer istate, itask, itol, iopt, mf, iwork(liw) double precision atol, rtol, t0, tf, rwork(lrw), y(neq) itol = 1 rtol = 1.0d-3 atol = 1.0d-3 iopt = 0 mf = 25 itask = 1 istate = 1 iwork(1) = 19*19 iwork(2) = 19*19 t0 = 0.0 tf = 0.7 call initial_condition( neq, t0, y ) call dlsode( f, neq, y, t0, tf, itol, rtol, atol, itask, istate, iopt, rwork, lrw, iwork, liw, jac, mf ) if (istate .lt. 0) then write(6,*) 'Failed istate=', istate else write (6,*) 'Steps ', iwork(11) write (6,*) 'F evals', iwork(12) endif |
History
This work started as a summer internship at Lawrence Livermore National Lab in 2005. Steven Lee, my mentor, gave much of the vision for the project. Pat Miller was our Python guru who answered many questions. After leaving LLNL, I re-implemented the interface and added some integrators written in Python. The new interface is a bit slimmer and uses NumPy instead of Numeric.
