[SciPy-Dev] ode and odeint callable parameters

[Benny Malengier]

I would like to come back to this issue.

I added a dopri wrapper to the odes scikit,
https://github.com/bmcage/odes/blob/master/scikits/odes/dopri5.py which
wraps the scipy.integrate calls to dopri5 and dop853.

When doing test with the Van der Pol example, I see this solver fails. The
seemingly good solution, is actually bad. In the test script you send, you
integrate via:

  sol1 = [r1.integrate(T) for T in tt[1:]]

What this does is actually raise an error after time 166.49 sec everytime,
but only one is printed to output, the others seemingly supressed, and you
keep restarting the solver via this call, even if not successful.
WIth the wrapper, after time 166.49, solutions never reach the output, they
go to
  sol1.errors
assuming you did
  sol1 = r1b.solve(tt, y0)
and as a user you must decide what to do.

This somewhat proves that returning output also on wrong output (not
satisfying atol/rtol) like ode.integrate does is not a good practice. In
the above case, with scipy.ode, you should test for r1.successful() after
every integrate, which you failed to do.

With the odes scikit approach, wrong output just stops the solve routine
and goes to sol.errors. The dopri warning output is to increase nsteps. If
you do that, and restart at the last computed solution, the error of dop853
becomes:
  UserWarning: dop853: problem is probably stiff (interrupted)
So an indication to use another solver.

Your van der pol example:
https://github.com/bmcage/odes/blob/master/docs/src/examples/ode/van_der_pol.py
which shows this behavior via the scikit API.

Benny


2016-01-26 15:08 GMT+01:00 Marcel Oliver <m.oliver at jacobs-university.de>:

> Benny Malengier writes:
>  >     I attach a small code snippet.  In fact, my memory failed me in that
>  >     dop853 indeed complains about insufficient nmax (nsteps in the
> Python
>  >     API).  What breaks is the accuracy of the solution, but it's still
>  >     remarkably good and can be controlled by nsteps.  It's obvious that
>  >     one should not solve this kind of problem with dop853, but I am
> still
>  >     amazed that it does not fall flat on its face.
>  >
>  >     In fact, vode/Adams dies completely (as it should) and vode/BDF
> gives
>  >     lots of warnings and computes a completely wrong solution (BDF
> schemes
>  >     are supposed to work here).  lsoda is fine, though.
>  >
>  > You forgot
>  > , with_jacobian=True
>  > when doing bdf metho via ode. So change into:
>  >
>  > r2 = ode(van_der_pol).set_integrator('vode', method='bdf',
> with_jacobian=True)
>  >
>  > Normally, as vode is for stiff methods, you set ml or mr to indicate how
>  > banded the jacobian is. Here you want full jacobian, so you should
> indicate
>  > this, see doc.
>
> Oops, I completely overlooked that vode is defaulting to simple
> fixed-point iteration in the BDF solver.  Thanks for pointing this
> out!
>
> I wonder if there is any situation where one would want to use BDF and
> still not want to use a full Newton or quasi-Newton iteration since
> the convergence condition for a simple fixed point iteration is
> basically as bad as using an explicit solver in the first place.  But
> there is certainly educational value to see this behavior...
>
>  > Can I use your code in an example under the scipy license ? I'll
>  > see what odes does.
>
> Absolutely.  It's pretty trivial of course, so feel free to use and/or
> adapt.
>
> Best,
> Marcel
>
> PS.: I'll try to get odes to compile later.
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> https://mail.scipy.org/mailman/listinfo/scipy-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20160129/9615f28c/attachment.html>