In this case, my rewrite idea works. For instance, the result from
trigsimp no longer has the removable singularity
In [10]: trigsimp(a)
Out[10]:
___
╲╱ 2 ⋅sin(q₂(t))⋅cos(q₁(t))
───────────────────────────
⎛ π⎞
2⋅sin⎜q₂(t) + ─⎟
⎝ 4⎠
In [12]: trigsimp(a).subs(S("q1(t)"), 0).subs(S("q2(t)"), 0)
Out[12]: 0
The Fu trigsimp algorithm is pretty customizable if you dig down into
it, so if you can figure out a reliable trigonometric transformation
that fixes this issue, you can apply it without the full power of
trigsimp. Or rewrite your algorithms to avoid them in the first place,
if you can pinpoint the source.
Aaron Meurer
On Wed, Sep 11, 2013 at 12:44 PM, Gilbert Gede <[email protected]> wrote:
> I think a really common example that came up in mechanics was with trig
> expressions. Here's a simple example (simple relative to actual mechanics
> expressions):
> sin(q1(t))*sin(q2(t))*cos(q1(t)) / (sin(q1(t))*sin(q2(t)) +
> sin(q1(t))*cos(q2(t)))
> When I substitute in q1(t) = 0, q2(t) = 0, I get nan.
>
> Generally, we don't know enough about the structure to know which order to
> substitute things in, i.e. of the coordinates (q's), which ones will end up
> in the denominator or even if you'll have an expression where both variables
> end up being able to cause a nan:
> sin(q1(t))*sin(q2(t))*cos(q1(t)) / (sin(q1(t))*sin(q2(t)) +
> sin(q1(t))*cos(q2(t))) +
> sin(q2(t))*sin(q1(t))*cos(q2(t)) / (sin(q2(t))*sin(q1(t)) +
> sin(q2(t))*cos(q1(t)))
>
> Some inspection/simplification seems to be required, which isn't practical
> for the larger expressions we can generate.
>
> Another question we can ask, relevant to this problem in mechanics, is when
> these nan's that pop up, should the result always be 0? That knowledge might
> be useful.
>
> -Gilbert
>
>
> On Wed, Sep 11, 2013 at 11:07 AM, Jason Moore <[email protected]> wrote:
>>
>> Thanks Aaron. The main assumption that we can make for our package is that
>> these variables are real numbers, so that will still leave us with divide by
>> zero issues. Maybe once I have a list of expressions that are common to our
>> problems, we can figure out how to deal with that subset.
>>
>>
>> Jason
>> moorepants.info
>> +01 530-601-9791
>>
>>
>> On Wed, Sep 11, 2013 at 1:45 PM, Aaron Meurer <[email protected]> wrote:
>>>
>>> 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.
>>
>>
>> --
>> 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.
