[SciPy-Dev] Updates to VODE and ZVODE solvers in single step mode

Thu Mar 5 14:22:33 EST 2015

Paul, All,

in light of ODE solvers, on the sundials list an interesting article
comparing solvers came up. This compares also CVODE with DVODE,
http://www.sciencedirect.com/science/article/pii/S0010465515000715

In all cases CVODE outperformed DVODE 2 to 5 times for this case.
Comparison also with lsode, ddaspk. Another datapoint that VODE used in
scipy is outdated.

In case you want to give it a shot, I maintain
https://github.com/bmcage/odes, one of the interfaces to sundials out there
which has an API like ode of scipy. Would it not be possible to change z in
x+iy transparently in a wrapper layer?

Greetings,
Benny

2015-03-02 13:30 GMT+01:00 Paul Nation <nonhermitian at gmail.com>:

> Benny,
>
> Many thanks for your help.  It seems like perhaps my usage is a bit
> specialized, and maybe does not warrant a pull to SciPy.
>
> In our particular case we also do some runtime Cython code generation for
> doing the RHS sparse matvec (if time-dependent) that I think would be a
> tricky thing to cast into real and imag parts without some major
> rewriting.  As we do quantum mechanics stuff, the vector space is naturally
> complex.
>
> Best,
>
> Paul
>
> On Mar 2, 2015, at 6:20 PM, Benny Malengier <benny.malengier at gmail.com>
> wrote:
>
>
>
> 2015-03-02 9:59 GMT+01:00 Paul Nation <nonhermitian at gmail.com>:
>
>> The zvode docs say that mode 5 returns the exact time answer and not an
>> interpolated result.
>>
>> Paul
>>
>
> Yes, I did not say mode 5 was not like that. I noted that the other mode
> also give you a good value at the wanted output time, only interpolated. In
> practise this would not be very different.
>
> Changing RHS as you have is a valid use of stop time. Just being able to
> stop on stop time is quite limiting however, more complicated solver have
> rootfinding, and can stop on every root found (stoptime is just root of
> t-stoptime).
>
> As to ZVODE, as far as I'm aware, nobody developed this further? One
> should take real and imaginary and convert to a CVODE problem I believe,
> but I never did complex problems. I do see
> http://netlib.sandia.gov/ode/zvode.f mentions such an approach:
>
> or a complex stiff ODE system in which f is not analytic,
> ZVODE is likely to have convergence
> failures, and for this problem one should instead use DVODE on the
> equivalent real system (in the real and imaginary parts of y).
>
> I don't mind somebody adding extra capabilities to VODE/ZVODE in scipy, if
> you do a nice PR it probably would be accepted if there is no influence on
> old code.
> A somewhat better approach in my view would be to deprecate it and convert
> to code which still sees releases. The API would be more complex though ...
> That has been discussed before though with no movement, people seem hapy to
> use ode.integrate for simple things, and then use other bindings when the
> problem outgrows VODE.
>
> Note that next release of the Sundials suite should be in the coming
> months. Over the years, one can assume bugs in VODE/ZVODE have been fixed
> in CVODE only.
>
> Benny
>
>
>
>>
>>
>> On Mar 2, 2015, at 17:43, Benny Malengier <benny.malengier at gmail.com>
>> wrote:
>>
>>
>>
>> 2015-03-02 6:49 GMT+01:00 Paul Nation <nonhermitian at gmail.com>:
>>
>>> When using the single-step mode in either the VODE or ZVODE ode solver,
>>> the default mode (2) called in:
>>>
>>>     def step(self, *args):
>>>         itask = self.call_args[2]
>>>         self.call_args[2] = 2  # Step mode is here
>>>         r = self.run(*args)
>>>         self.call_args[2] = itask
>>>         return r
>>>
>>> results in taking a single step that (typically) goes beyond the output
>>> time requested in the solver.  When doing, for example, monte carlo
>>> algorithms, this leads to a big performance hit because one must take a
>>> step back, reset the solver  and then use the normal mode to go to the
>>> requested stop time.  Instead, these solvers support a mode (5) that will
>>> never step beyond the end time.  The modified step function is in that case:
>>>
>>
>> You do obtain the output at the requested time though, it is only an
>> interpolated value of the actually computed solutions.
>> So, the only reason you should do above is if something is happening at a
>> certain time and you want to change data or so. You mention monte carlo,
>> but I don't see how that is related to such a usecase in general. I suppose
>> in your application probably yes, but in general MC does not need this
>>
>> The docs say only to use endtime if you have changing RHS or Jacobian,
>> and to otherwise not try to outsmart the solver, as the solver needs extra
>> work in case you set the endtime.
>>
>> Note that (Z)VODE was replaced by CVODE by the authors of VODE, which has
>> many improvements and several python bindings, all of which expose setting
>> a stop time. In my view, VODE is only present in scipy as a first attempt
>> solver, to be replaced by more modern solvers for heavy lifting.
>>
>> Benny
>>
>>
>>>
>>> def step(self, *args):
>>>         itask = self.call_args[2]
>>>         self.rwork[0] = args[4]    #Set to stop time
>>>         self.call_args[2] = 5       #Set single step mode to stop at
>>> requested time.
>>>         r = self.run(*args)
>>>         self.call_args[2] = itask
>>>         return r
>>>
>>> Currently in order to implement this, one needs to create their own ODE
>>> integrator subclass of VODE or ZVODE, overload the step function, then
>>> create an ode instance and then finally add the custom integrator using
>>> ode._integrator.  I think supporting both options natively would be a nice
>>> thing to have in SciPy.
>>>
>>> In addition, often it is not necessary to do a full reset of the ode
>>> solver using ode.reset().  Often times one just needs to change the RHS
>>> vector (and possibly the time) and set the flag for the solver to start
>>> anew (ode._integrator.call_args[3] = 1).  This to results in a large
>>> performance benefit for things like monte carlo solving. Right now I need
>>> to call
>>>
>>> ode._y = new_vec
>>> ode._integrator.call_args[3] = 1
>>>
>>> when I want to accomplish this.  Adding support for a “fast reset” might
>>> also be a good thing to have in SciPy.
>>>
>>> All of the code to accomplish such things are already being used in the
>>> QuTiP monte carlo solver(
>>> https://github.com/qutip/qutip/blob/master/qutip/mcsolve.py) and would
>>> therefore be fairly painless to add to SciPy.
>>>
>>> Best regards,
>>>
>>> Paul
>>>
>>>
>>>
>>> _______________________________________________
>>> SciPy-Dev mailing list
>>> SciPy-Dev at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>>>
>>>
>> _______________________________________________
>> SciPy-Dev mailing list
>> SciPy-Dev at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>>
>>
>> _______________________________________________
>> SciPy-Dev mailing list
>> SciPy-Dev at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>>
>>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
>
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20150305/c8fa341c/attachment.html>