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

Thu Aug 22 21:03:06 EDT 2013

On 8/22/13, Juan Luis Cano <juanlu001 at gmail.com> wrote:
> 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
>

Luis, thanks for working on this.

In this line:
    https://github.com/Juanlu001/scipy/blob/integrate2/scipy/integrate/tests/test_lsoda2.py#L6
the coefficient of `y[1]*y[2]` should be `1.0e4` (not `1.0`).

Warren

> 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"]
>
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>