[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