[SciPy-User] lsoda vs. Coulomb friction

[Warren Weckesser]

Looks like Anne's email arrived while I was finishing mine.  I could 
have save myself some time by simply saying "Yeah, what Anne said."  :)

Warren


Anne Archibald wrote:
> 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
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>