[SciPy-user] Runge-Kutta ODE integrator in SciPy, odepack problems on SciPy Superpack for OSX?

Nathan Bell wnbell@gmail....
Mon Mar 17 01:43:28 CDT 2008


On Mon, Mar 17, 2008 at 12:00 AM, Anne Archibald
<peridot.faceted@gmail.com> wrote:
>  I have to disagree with this. Getting numerical code right, and
>  robust, and bug-free, requires a lot of knowledge, skill and care.

As does wrapping libraries.  Keep in mind that only a minority of
SciPy users/contributors will ever know SWIG or f2py, etc.  Keeping
things in Python greatly enhances accessibility of the implementation.

>  Even if it's only thirty lines. For example, suppose you implement
>  embedded adaptive RK 4/5 integration. How do you design a test case
>  that ensures you didn't get one of the coefficients wrong? After all,
>  the adaptiveness means that even quite serious miscalculations often
>  still converge to the right answer, they just take many more
>  iterations than they should.

How do you know that you didn't introduce a subtle error when wrapping
the library?

In either case you need comprehensive unit tests.  I don't see how
this is a win for existing codes.

> Or what about stiff ODEs?  Simple RK schemes will take *forever*.

Don't move the goalposts, the OP wanted an RK method.

In the case of stiff solvers, which can be considerably more
complicated, wrapping an existing code may prove to be the better
option.

>  Better, if you can, to use an ODE solver
>  somebody else has beaten the bugs out of over many years.  Especially
>  if all it takes is loading it in from scipy.

Wrapping a library and "loading it in from scipy" is not always a
trivial exercise.

What if the interface is a poor match for wrapping?
What if the library only supports double precision FP and you want
singles or extended precision?
What if the callback procedure is tedious/inflexible?
What if the library doesn't support RK23 and that's what you really want?
What if the library doesn't provide all the information that you might want?
What happens when the person who wraps the code abandons SciPy?
What if the original authors of the library abandon their code?
What if the library changes licenses?

All of these problems are easily solved by a Python implementation.

> Focus the effort that would have gone into debugging and testing your RK solver on solving
>  the quirks of your own problem.

You ignore all the potential problems that may arise when wrapping
someone else's code and the difficulty with providing a Pythonic
interface.

>  Once in a while you'll run into a problem when you need more
>  performance, or different abilities, or special techniques, that the
>  library doesn't cover. But until that happens, why waste your time
>  reinventing the wheel? Try the stock solutions first. If they don't
>  work, understand why, and whether reimplementing them will solve the
>  problem or just hide it.

In this case the wheel is 20 lines of dead-simple Python + NumPy code.
 I'd wager that the handling of parameters and callbacks is more
effort than the numerical component.

IMO the pure Python + NumPy approach wins hands down here.

-- 
Nathan Bell wnbell@gmail.com
http://graphics.cs.uiuc.edu/~wnbell/


More information about the SciPy-user mailing list