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@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-0002.pl
-------------- 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-0001.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-0003.pl