[SciPy-dev] Test case for integrate/ode.py: banded systems

[Jesper Friis]

Andrew Straw wrote:
> And now, upon reading the comments in the code of that page, it seems
> you've uncovered what may be a bug in the VODE wrapper. Would you say it
> is? If so, we should report the bug, hopefully fix it either ourselves
> or have a developer do it, and fix up the page appropriately.
>
> Cheers!
> Andrew


I have to admit that I just started last week to use
python instead of Fortran, so I think that you are in a better
position to say if it should be considered a bug or not that the
nrowpd argument is missing in the python interface.

Searching through the function integrate/odepack/vode.f revealed that
the Jacobian function is only called in two places. For full systems
it is called with nrowpd=neq while for banded systems it is called
with nrowpd=ml*2+mu+1. So, as long as these values for nrowpd are
provided in the documentation there should be no problems using the
python interface.

Note however that the lower and upper bandwidths, ml and mu are
interchanged compared to what I wrote in the comment in the example
script. The reason for this turned out to be another little miss in
integrate/ode.py, which is corrected with the patch below.

It might also be interesting (at least for me) to add a way to obtain
information about the solution, like the number of steps that was
needed to reach the solution. What is the best way to do this?  Maybe
one could add a method called infodict() that returns a dictionary
containing all optional output from the solver?

In the end of the ode module there are two simple tests for non-banded
systems. Maybe they should be added to the Cookbook together with my
banded example. What do you think?


Regards
/Jesper



--------------------------------------------------------------

--- scipy-0.4.6/Lib/integrate/ode.py.org        2006-03-07
23:23:17.000000000 +0100
+++ scipy-0.4.6/Lib/integrate/ode.py    2006-03-07 23:23:33.000000000 +0100
@@ -276,8 +276,8 @@
          self.with_jacobian = with_jacobian
          self.rtol = rtol
          self.atol = atol
-        self.mu = lband
-        self.ml = uband
+        self.mu = uband
+        self.ml = lband

          self.order = order
          self.nsteps = nsteps