In the mechanics package we sometimes end up with symbolic expressions that
can give a divide by zero error if we use the expression in a numerical
function. Below is a simple example of an expression that could be
generated by the mechanics package. It should evaluate to zero (i.e the
limit as vx -> zero. evalf says it's zero but if I use lambdify (or
probably any code generator) to transform the expression to something that
is purely numeric we get divide by zeros. I feel like this must be a common
issue when transforming symbolic expressions to numerical functions. Does
anyone know how or if this is dealt with? It be nice if we could generate
equations of motion in the mechanics package and generate code from them
that wouldn't have divide by zero errors, requiring some fix on the
numerical side.
In [1]: from sympy import symbols, sqrt, lambdify
In [2]: k, vx, vy, vz = symbols('k vx vy vz')
In [3]: f = k * sqrt(vx ** 2 + vy ** 2 + vz ** 2)
In [4]: dfdvx = f.diff(vx)
In [5]: dfdvx
Out[5]: k*vx/sqrt(vx**2 + vy**2 + vz**2)
In [6]: dfdvx.evalf(subs={k: 1.0, vx: 0.0, vy: 0.0, vz: 0.0})
Out[6]: 0
In [7]: numeric_dfdvx = lambdify((k, vx, vy, vz), dfdvx)
In [8]: numeric_dfdvx(1.0, 0.0, 0.0, 0.0)
---------------------------------------------------------------------------
ZeroDivisionError Traceback (most recent call last)
<ipython-input-8-55e44f8ab1bf> in <module>()
----> 1 numeric_dfdvx(1.0, 0.0, 0.0, 0.0)
/usr/local/lib/python2.7/dist-packages/numpy/__init__.pyc in
<lambda>(_Dummy_26, _Dummy_27, _Dummy_28, _Dummy_29)
ZeroDivisionError: float division by zero
In [9]: numeric_dfdvx = lambdify((k, vx, vy, vz), dfdvx, modules='numpy')
In [10]: numeric_dfdvx(1.0, 0.0, 0.0, 0.0)
/usr/local/lib/python2.7/dist-packages/numpy/__init__.py:1: RuntimeWarning:
invalid value encountered in double_scalars
"""
Out[10]: nan
In [11]: numeric_dfdvx = lambdify((k, vx, vy, vz), dfdvx, modules='mpmath')
In [12]: numeric_dfdvx(1.0, 0.0, 0.0, 0.0)
---------------------------------------------------------------------------
ZeroDivisionError Traceback (most recent call last)
<ipython-input-12-55e44f8ab1bf> in <module>()
----> 1 numeric_dfdvx(1.0, 0.0, 0.0, 0.0)
<string> in <lambda>(k, vx, vy, vz)
/usr/local/lib/python2.7/dist-packages/sympy/mpmath/ctx_mp_python.pyc in
__rdiv__(s, t)
206 if t is NotImplemented:
207 return t
--> 208 return t / s
209
210 def __rpow__(s, t):
/usr/local/lib/python2.7/dist-packages/sympy/mpmath/ctx_mp_python.pyc in
__div__(self, other)
/usr/local/lib/python2.7/dist-packages/sympy/mpmath/libmp/libmpf.pyc in
mpf_div(s, t, prec, rnd)
928 if not sman or not tman:
929 if s == fzero:
--> 930 if t == fzero: raise ZeroDivisionError
931 if t == fnan: return fnan
932 return fzero
ZeroDivisionError:
In [19]: dfdvx.limit(vx, 0.0)
Out[19]: 0
Jason
moorepants.info
+01 530-601-9791
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/sympy.
For more options, visit https://groups.google.com/groups/opt_out.
