[SciPy-Dev] Loss of precision using lsoda f2py interface or ode class

[Juan Luis Cano]

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"]