[SciPy-User] lsoda vs. Coulomb friction

[Warren Weckesser]

Hi Ryan,

I have two suggestions, but the second option requires using a 
capability that the scipy ODE solvers currently don't provide.

Since the underlying problem is that you are using discontinous 
equations with a solver that expects smooth equations, you can replace 
your discontinuous function with a smooth approximation.  The attached 
script, coulomb_approx.py, includes an approximation for your friction 
terms.  It generates the attached PNG file.

Also attached is lsoda_vs_coulomb_friction2.py, which is a modified 
version of your script.  It computes the solution twice, once as in the 
original, and once with the approximation.  With the smooth 
approximation, lsoda does not complain, and--at least at the resolution 
of the pylab plot--you can not distinguish the two solutions.   Whether 
an approximation like this is good enough depends on the ultimate goals 
of your analysis.  There maybe behavior very close to v=0 that you want 
to observe but that are not demonstrated with the smooth approximation.

The second option is to use an ODE solver with "root finding" (aka 
"event detection").  You would set the solver to stop when it attempted 
to cross v=0.  Depending on the direction of the crossing, it would 
change a parameter appropriately and then resume from the point where it 
stopped.  Unfortunately, the current ODE solvers in scipy do not include 
this capability (and "rolling your own" can be hard to get right).

Warren


Ryan Krauss wrote:
> I am trying to use odeint (i.e. lsoda) on a mechanical system that
> involves a mass and friction modeled as a viscous term plus Coulomb
> friction.  The discontinuity near 0 velocity makes lsoda mad.  It
> complains that it has to do excess work (lazy algorithm :).  In spite
> of its complaining, the result leads to fairly good agreement between
> model and experiment.  But is there a better way to handle this?
>
> This is the func I am passing to odeint (see attached example script
> for the details):
>
> def ydot_ol(x, t, u, C):
>     m = C[0]
>     b1 = C[1]
>     b2 = C[2]
>     y1 = x[0]
>     y2 = x[1]
>     ydot = zeros(2,)
>     ydot[0] = y2
>     ydot[1] = (u - b1*y2 - b2*sign(y2))/m
>     return ydot
>
> u is an external input.
>
> lsoda doesn't complain until the system is essentially stopping.
>
> Thanks,
>
> Ryan
>   
> ------------------------------------------------------------------------
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: coulomb_approx.py
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100203/71418db3/attachment.ksh>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: coulomb_approx.png
Type: image/png
Size: 20128 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100203/71418db3/attachment.png>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: lsoda_vs_coulomb_friction2.py
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100203/71418db3/attachment-0001.ksh>