[SciPy-User] lsoda vs. Coulomb friction

[Anne Archibald]

2010/2/3 David Goldsmith <d.l.goldsmith at gmail.com>:

> "   :
>     :
> Returns
>  -------
>     :
>     :
>  infodict : dict, only returned if full_output == True
>      Dictionary containing additional output information
>         :
>         :
>      'mused'  a vector of method indicators for each successful time step:
>               1: adams (nonstiff), 2: bdf (stiff)"
>
> Note that it's a _return_ value: apparently with odeint, you can't specify
> which method to use, you can only "hear" which method the algorithm decided
> to use for you.

In particular, odeint always (?) starts with a non-stiff method and
switches to a stiff method if the problem appears to warrant it. It
does not seem that odeint ever switches back to a non-stiff method
even if the problem enters a non-stiff region.

In any case, I'm not sure there's any reason to think that a stiff
solver will help with the OP's problem, which I think is coming from
the discontinuity. I suspect what's happening is that initially the
losses are dominated by the viscosity term, so that the discontinuity
doesn't bother the solver much. But when the velocities drop and the
discontinuous term begins to dominate, the solver tries to match it
better and better. Since the solver tries to approximate the solution
with a 4th or 5th order polynomial but the actual solution is probably
only differentiable once (or not at all), the solver goes crazy trying
to fit a kink with a curve, and you get zillions of basically useless
steps.

I can see three ways to solve the problem:

* Declare that when the discontinuity becomes important, the Coulomb
friction model becomes too crude an approximation, and stop.
(Obviously if the results agree with experiment this is unnecessary.)

* Come up with a model that describes in more detail what happens at
that transition; perhaps a smooth roll-off of forces from positive to
negative, perhaps an additional "stuck area" variable that increases
near zero velocity and that provides resistance to movement, or
perhaps some more sophisticated friction model from the literature.
(Unfortunately the tribology prof I used to work with went back to
Germany, or I'd ask if he had some suggestions.)

* Model it piecewise: run the solver with a RHS that assumes the
object is moving to the right, and stop when the velocity hits zero.
Then figure out whether and which way the object starts moving again,
and start up a new solver. This requires you to have a solver that can
stop when a condition is reached, which is a pain with the scipy
solvers. (I think that odeint can be used in a single-step mode that
allows you to "back up" using interpolation within the last step, but
really the right solution is to use a wrapper for lsodar, which
supports stopping conditions directly.)

Anne