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

[Anne Archibald]

On 16/03/2008, Nathan Bell <wnbell at gmail.com> wrote:
> On Thu, Mar 13, 2008 at 6:07 PM, Zane Selvans <zane at ideotrope.org> wrote:
>  >
>  >  The only reason I ask about Runge-Kutta specifically is I know two
>  >  people who have the solution to my problem coded up already, one in
>  >  Fortran, and one in C, and they both used a Runge-Kutta integrator.  I
>  >  want to open-source my model code, but it depends on their codes, and if
>  >  I can't get them to let me publicize their work, I'm going to have to
>  >  re-write it from scratch... unless someone else has already done it.
>  >
>
> >  I don't know why they would have both chosen to write their own
>  >  numerical solutions from scratch if something publicly available would
>  >  have worked... but I guess it's possible.  A lot of people don't seem to
>  >  like to build on the work of others.
>
>
> As Andrew mentioned, writing a generic RK method (that accepts used
>  defined derivatives) is pretty trivial.  If you stored the
>  coefficients in the so-called "Butcher arrays" as arrays, as opposed
>  to hard coding for a particular order (e.g. RK2 or RK4), then you can
>  unify all (explicit) RK methods easily too.  Furthermore, adaptive RK
>  methods (e.g. RK45) are simply a post-process after a standard RK step
>  that can be handled in a uniform fashion.
>
>  I suspect that you could write a generic RK scheme in 20-30 lines by
>  using numpy for the vector operations .  There's really no reason to
>  use someone else's C/Fortran code for this problem (and many reasons
>  not to).

I have to disagree with this. Getting numerical code right, and
robust, and bug-free, requires a lot of knowledge, skill and care.
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. Or what about stiff ODEs? Simple RK
schemes will take *forever*. 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. Focus the effort that
would have gone into debugging and testing your RK solver on solving
the quirks of your own problem.

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.

Anne

P.S. your problem is really making me curious - are you really
simulating self-gravitating viscoelastic material? Is this geophysics?
Somehow that sounds more like a PDE type of problem than ODE; PDE
solvers are much less standardized and I don't know that there are
stock tools to attack them.  -A