[SciPy-user] integrate ODE to steady-state?

[Anne Archibald]

2008/6/24 Zachary Pincus <zachary.pincus at yale.edu>:

> My problem is similar, actually, but in 2D. I have an outline of a
> bacteria as a 2D polygon and want to lay down a coordinate system on
> it, with a midline running from pole to pole ("x-axis"), and well-
> defined lines normal to the midline ("y-axis"). The classic solution
> would be to use the medial axis transform, which I had been (using
> Martin Held's excellent if quirky VRONI code). However, the medial
> axis is only a partial midline, it may branch, and there's no
> guarantee that it is possible to find reasonable straight-line normals
> to the axis (you get a clearance radius at each point instead). So I
> formulated the problem as trying to find pairs of positions along the
> outline of the shape where the distance between the pairs is
> minimized, but adjacent points don't "clump" too much, which naturally
> led to this sort of spring-system picture. (This is a very brief
> explanation -- sorry if unclear.) It doesn't seem particularly elegant
> to solve it out like this, but I wanted to see if it worked at all...
>
> Also, thanks for the references -- interesting indeed.
>
> Rob and Lou:
>> [try solving for dy/dt=0]
>
> Good idea too! I'll try that one as well, and see which seems to work
> better...

I'd set it up as a numerical minimization procedure - minimizing the
"potential energy" of all those springs. If you give it a good initial
guess and are lucky, the minimizer will follow a path down to the
optimum not too different from the one your ODEs would have followed,
but since the optimizer knows you don't care how it gets there it can
be much more efficient. The danger with both approaches is that there
may be some rather bad local minimum they can get trapped in. The best
approach is to choose a decent initial guess - maybe using the medial
axis transform, since you already have the code working - so that the
nearest, most obvious solution is the real optimum. Failing that, you
can look at some of the global optimizers in OpenOpt.

Better, of course, is to find an algorithm specifically adapted to
solving your particular problem, but I gather you've already done some
considerable searching.

Anne