[SciPy-Dev] Loss of precision using lsoda f2py interface or ode class
Juan Luis Cano
juanlu001 at gmail.com
Thu Aug 22 20:15:34 EDT 2013
I've been interested in solving some of the bugs of
scipy.integrate.odeint, and there was some discussion about how to do it
here
https://github.com/scipy/scipy/issues/2570
However, I've been trying to wrap my head around lsoda.f and its f2py
interface and things are being very difficult. In lsoda.f there is a
sample program with the expected output, and I tried to replicate it to
begin with.
I'm working on this branch:
https://github.com/Juanlu001/scipy/tree/integrate2/scipy/integrate
These tests were added as three files in tests/: a Fortran program, a
Python program using lsoda.pyf directly and another Python program using
the ode class. See below for the outputs.
For some reason the three programs show different outputs. The
difference becomes bigger and bigger with each step.
I'm asking for some advice on how to debug this effectively, because
lsoda.f is a hell of GO TO and, to be honest, I'm not smart enough to
stay sane while trying to follow the logic. I already tried f2py
--debug-capi but I see too much noise and nothing useful. I also tried
changing the dtypes of the arrays to float, or even dtype('f8') to make
sure I'm using double precision everywhere, and though the results
change there are still large errors.
There may be also errors in translating the original from Fortran to
Python, but I am not able to spot them.
I fear odeint/ode bugs are starting to pile up (gh-1567, gh-1801,
gh-1976, gh-2515, gh-2570), and as many have suggested in the past a
rewrite or redesign would be quite helpful.
Maybe someone wants to take this and work a little during the upcoming
sprint.
Thank you in advance
Juan Luis Cano
--
This is the output of the Fortran program (good):
$ ./test_lsoda
at t = 4.0000E-01 y = 9.851712E-01 3.386380E-05 1.479493E-02
at t = 4.0000E+00 y = 9.055333E-01 2.240655E-05 9.444430E-02
at t = 4.0000E+01 y = 7.158403E-01 9.186334E-06 2.841505E-01
at t = 4.0000E+02 y = 4.505250E-01 3.222964E-06 5.494717E-01
at t = 4.0000E+03 y = 1.831976E-01 8.941773E-07 8.168015E-01
at t = 4.0000E+04 y = 3.898729E-02 1.621940E-07 9.610125E-01
at t = 4.0000E+05 y = 4.936362E-03 1.984221E-08 9.950636E-01
at t = 4.0000E+06 y = 5.161833E-04 2.065787E-09 9.994838E-01
at t = 4.0000E+07 y = 5.179804E-05 2.072027E-10 9.999482E-01
at t = 4.0000E+08 y = 5.283679E-06 2.113483E-11 9.999947E-01
at t = 4.0000E+09 y = 4.658659E-07 1.863464E-12 9.999995E-01
at t = 4.0000E+10 y = 1.431333E-08 5.725337E-14 1.000000E+00
no. steps = 361 no. f-s = 693 no. j-s = 64
method last used = 2 last switch was at t = 6.0092E-03
This is the output of the Python program using the Fortran interface
directly:
$ python test_lsoda.py
At t = 4.0000e-01 y = [ 9.84128935e-01 3.62239836e-05 1.58348410e-02]
At t = 4.0000e+00 y = [ 8.52155561e-01 3.37055417e-05 1.47810734e-01]
At t = 4.0000e+01 y = [ 2.02104281e-01 1.64041139e-05 7.97879315e-01]
At t = 4.0000e+02 y = [ 1.06122917e-06 2.44976665e-08 9.99998914e-01]
At t = 4.0000e+03 y = [ 6.31790031e-09 2.50796401e-10 9.99999993e-01]
At t = 4.0000e+04 y = [ 6.31740542e-10 2.52488308e-11 9.99999999e-01]
At t = 4.0000e+05 y = [ 4.98669737e-10 1.99540660e-11 9.99999999e-01]
At t = 4.0000e+06 y = [ -7.52067881e-10 -3.00837491e-11 1.00000000e+00]
At t = 4.0000e+07 y = [ -5.25052749e-10 -2.10020683e-11 1.00000000e+00]
At t = 4.0000e+08 y = [ 5.77160490e-10 2.30864142e-11 9.99999999e-01]
At t = 4.0000e+09 y = [ -3.33188378e-09 -1.33275352e-10 1.00000000e+00]
At t = 4.0000e+10 y = [ -2.20370350e-09 -8.81481397e-11 1.00000000e+00]
No. steps = 251, no. f-s = 614, no. j-s = 72
method last used = 2, last switch at t = 6.0089e-03
(note: the last two lines are garbage with current SciPy master - see
https://github.com/Juanlu001/scipy/commit/e11e978232c14b65f800922e7390b08bd8298734)
And this is the output of the program using the ode class in
scipy.integrate:
$ python test_lsoda2.py
At t = 4.0000e-01 y = [ 9.84127434e-01 3.62239556e-05 1.58363421e-02]
At t = 4.0000e+00 y = [ 8.52153740e-01 3.37055100e-05 1.47812555e-01]
At t = 4.0000e+01 y = [ 2.02145824e-01 1.64042718e-05 7.97837771e-01]
At t = 4.0000e+02 y = [ 1.05528994e-06 2.44971521e-08 9.99998920e-01]
At t = 4.0000e+03 y = [ 6.39513550e-09 2.53944135e-10 9.99999993e-01]
At t = 4.0000e+04 y = [ 5.53276503e-10 2.21169099e-11 9.99999999e-01]
At t = 4.0000e+05 y = [ 5.47740551e-11 2.19082310e-12 1.00000000e+00]
At t = 4.0000e+06 y = [ 6.03428535e-12 2.41369408e-13 1.00000000e+00]
At t = 4.0000e+07 y = [ 7.15628192e-13 2.86251003e-14 1.00000000e+00]
["Excess work done" errors"]
More information about the SciPy-Dev
mailing list