With your example, it matters what order you substitute the variables. If you substitute vx=0 first, you get 0/sqrt(vy**2 + vz**2) == 0. If you substitute vy=vz=0, you get vx/sqrt(vx**2), which should simplify to sqrt(vx) for vx >= 0 (which it is, since vx=0). Can you assume nonnegative for these variables. Of course, if you substitute all three at once, you get 0/0.
There's no general solution to this problem, but you aren't working in a completely general situation. I think the solution you should adopt will depend on just how much structure you know about your expressions. Sometimes you can rewrite expressions so that they don't have 0/0 (like you learn in calculus). I don't see how to do it for this expression, though. Aaron Meurer On Tue, Sep 10, 2013 at 2:43 PM, Jason Moore <[email protected]> wrote: > Also, why does .evalf give zero in my previous example? Shouldn't that also > raise a DivideByZero error? > > > Jason > moorepants.info > +01 530-601-9791 > > > On Tue, Sep 10, 2013 at 4:41 PM, Jason Moore <[email protected]> wrote: >> >> I'll don't have a concrete example to show at the moment, but I feel like >> we get these singularity issues when we do use simplification and don't when >> we don't or when we expand everything out. I'll try to put together an >> example where this is the case. >> >> >> Jason >> moorepants.info >> +01 530-601-9791 >> >> >> On Tue, Sep 10, 2013 at 4:14 PM, Stefan Krastanov >> <[email protected]> wrote: >>> >>> > In this case, real and positive assumptions don't help for the lambdify >>> > call. I guess making lambdify assumption aware could be an option. >>> >>> I was actually expecting that one would always use `simplify` before >>> `lambdify`. This would make a big difference, because with the >>> appropriate assumptions `simplify` would remove all such removable >>> singularities[1], which is the technical term for these point. >>> Unfortunately a quick google search for "numerical removable >>> singularity" does not show anything promising. >>> >>> [1]: https://en.wikipedia.org/wiki/Removable_singularity >>> >>> On 10 September 2013 14:19, Jason Moore <[email protected]> wrote: >>> > In this case, real and positive assumptions don't help for the lambdify >>> > call. I guess making lambdify assumption aware could be an option. >>> > >>> > I think we should define dynamicsymbols in mechanics to have a default >>> > real >>> > assumption, this could help in simplification. >>> > >>> > The different basis is an interesting idea. Do you have an example of >>> > that? >>> > >>> > >>> > Jason >>> > moorepants.info >>> > +01 530-601-9791 >>> > >>> > >>> > On Tue, Sep 10, 2013 at 2:01 PM, Stefan Krastanov >>> > <[email protected]> wrote: >>> >> >>> >> The poorman's solution that I used when I had this problem was to do >>> >> the calculations in a different basis. I do not know how >>> >> general/automatic this solution would be. >>> >> >>> >> However, a common reason for bad/no simplification is lack of >>> >> appropriate assumptions. Define all you symbols as real/positive and >>> >> half of the issues will go away. >>> >> >>> >> On 10 September 2013 12:17, Jason Moore <[email protected]> wrote: >>> >> > 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. >>> >> >>> >> -- >>> >> 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. >>> > >>> > >>> > -- >>> > 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. >>> >>> -- >>> 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. >> >> > > -- > 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. -- 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.
