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.
