[SciPy-Dev] GSoC' 15 Ideas: define general integrators; parallelize integration; higher-order ODEs

[Geordie McBain]

2015-03-17 3:37 GMT+11:00 Max Mertens <max.mail at dameweb.de>:
> Hello everyone,
>
> I am a bit late to write to this mailing list, therefore sorry for that.
>
> My name is Max Mertens; I am a Electrical Engineering student at the
> University of Ulm, Germany.
> I have done a few projects in C(++) on physical simulations of multibody
> dynamics with collision detection, and molecular dynamics with
> interatomic potentials.
>
> I got interested in SciPy as it facilitates physical simulations in an
> easy way, but have not worked with SciPy so far. I would like to
> contribute to SciPy as part of the GSoC, and am willing to contribute by
> providing an "easy-fix" to one of the open issues.
> Is there already a possibility to perform multibody/molecular dynamics
> in SciPy or are there contributions needed in this scope?
>
> The scipy.integrate module with ode and odeint provides several methods
> to integrate differential equations. For this I have the following ideas:
>
>   * add a interface/class/method to define general integration methods:
>     For now you can specify various integrators like "vode" and "dopri5"
>     as a string. The new code would allow to enter a Butcher tableau to
>     define implicit/explicit (embedded) Runge-Kutta methods, which would
>     cover "dopri5" and "dopri853" (see [0]); and possibly other general
>     integrator descriptions as needed.

Excellent idea.  I would like this very much.  I had to write
something like this with Butcher tableaux for general implicit
Runge-Kutta methods (for differential-algebraic equations rather than
ordinary differential equations), but it's a bit in-house and
application-specific for publishing.  Having this as part of SciPy
would be great.

>   * add distributed integration: Linear multistep integrators like
>     Runge-Kutta with multiple differential equations can be parallelized
>     to speed up calculations. This would distribute integrations to
>     multiple threads; and/or, if needed for complicated equations like
>     in multibody dynamics, distribute to multiple physical machines (I
>     have developed a few ideas on how to accomplish this).

Another excellent idea and something I've had on my own TODO list for
a while.  I'm not sure what the SciPy policy on parallelization is
though; SciPy doesn't seem to deal with it much.
   I would be very interested to read your ideas on this.  Ultimately
I suppose I would like my integrator to use a sparse parallel linear
solver.

>   * provide methods to integrate higher-order differential equations: is
>     this needed, or are users of the library expected to express these
>     as multiple first-order differential equations? Could this step be
>     automated?

I think this could be worthwhile too.  At the very least for
second-order equations which arise so often in physical applications
and for which there are popular special methods like
https://en.wikipedia.org/wiki/Newmark-beta_method and
https://en.wikipedia.org/wiki/Verlet_integration.

Another reason I had to write my own implicit Runge-Kutta methods was
because the differential equations that I was dealing with weren't
'ordinary' in the sense that they could be solved for the first
derivative so were 'implicit' or 'differential-algebraic'.  This
wasn't a particularly difficult generalization but is not currently
covered by anything in scipy.optimize; e.g. scipy.optimize.ode insists
on a form y' = f (t, y) but my equation was more like f (t, y, y') =
0, as in

Hairer, E., C. Lubich, & M. Roche (1989). The numerical solution of
differential-algebraic systems by Runge-Kutta methods, Volume 1409 of
Lecture Notes in Mathematics. Berlin: Springer

Brenan, K. E., S. L. Campbell, & L. R. Petzold (1996). Numerical
solution of initial-value problems in differential-algebraic
equations, Volume 14 of Classics in Applied Mathematics. Philadelphia:
Society for Industrial and Applied Mathematics

DAEs can be much trickier than ODEs though (as described in those two
books and Petzold's other papers), so it is harder to write robust
general-purpose programs for them for inclusion in something as
high-level as SciPy; e.g. this is a large part of why I describe my
solutions as too application-specific.

Good luck.  There is much worthwhile work to be done here.