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

Warren Weckesser warren.weckesser@gmail....
Thu Aug 22 20:10:41 CDT 2013


On 8/22/13, Warren Weckesser <warren.weckesser@gmail.com> wrote:
> On 8/22/13, Juan Luis Cano <juanlu001@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.
>

*Juan*, that is--sorry.

The other python file (test_lsoda.py) also needs the correction to the
coefficient.

Warren

> 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@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>>
>


More information about the SciPy-Dev mailing list