What's the best practice for using simplify?

[neozhaoliang]

Hi, I have a short code that computes the catacaustic of a circle, which should return a cardioid:

from sympy import *

t, X, Y = symbols("t X Y")
x = cos(t)
y = sin(t)
curve = Matrix([x, y])
light_source = Matrix([1, 0])
ray = curve - light_source
dx = diff(x, t)
dy = diff(y, t)
n = Matrix([dy, -dx])
reflected_ray = simplify(ray - 2 * ray.dot(n) * n / n.dot(n))
F = (Y - y) * reflected_ray[0] - (X - x) * reflected_ray[1]
dF = diff(F, t)
result = solve((F, dF), X, Y)
print(f"X(t)={simplify(result[X])}")
print(f"Y(t)={simplify(result[Y])}")

sympy returns the following result:

X(t)=(-3*cos(t) + cos(3*t) + 2)/(6*(cos(t) - 1))
Y(t)=2*sin(t)**3/(3 - 3*cos(t))

The result is correct but it's not in the most simplified form, which should be

$$x(t)=\frac{\cos(2t) + 2\cos(t)}{3}$$

$$ y(t)=\frac{\sin(2t) + 2\sin(t)}{3}.$$

I verified manually that the above is equivalent to the result by sympy.

My question is, is there any possible improvement of my usage of sympy, so that I can get a more simplified form of the result? For example, to eliminate the denominator?

[asmeurer]

simplify() tries various simplification functions, but not all of them. To get better trig simplification you should use trigsimp directly, and try the different options documented in the method argument https://docs.sympy.org/latest/modules/simplify/simplify.html#module-sympy.simplify.trigsimp. You can also use the exact simplification functions from the fu module if you know which ones you want. https://docs.sympy.org/latest/modules/simplify/fu.html

In this case, the fu function does a pretty bad job at simplifying this function:

>>> pprint(fu(result[X]))
                                         2
 2⋅sin(2⋅t) + sin(4⋅t)     (cos(2⋅t) + 1) ⋅sin(2⋅t)          5⋅sin(2⋅t)                  sin(6⋅t)                   4⋅sin(t)
──────────────────────── - ──────────────────────── - ──────────────────────── + ───────────────────────── + ───────────────────────
2⋅(-2⋅sin(t) + sin(2⋅t))   3⋅(-2⋅sin(t) + sin(2⋅t))   4⋅(-2⋅sin(t) + sin(2⋅t))   12⋅(-2⋅sin(t) + sin(2⋅t))   -12⋅sin(t) + 6⋅sin(2⋅t)

But the groebner method is able to give the result you are expecting.

>>> pprint(trigsimp(result[X], method='groebner'))
2⋅cos(t)   cos(2⋅t)
──────── + ────────
   3          3

No doubt fu could be improved here. I would say this is a pretty catastrophic failure in its default heuristics.

[neozhaoliang]

I love this solution!

[asmeurer]

I guess to be fair, result[X] is pretty complicated to begin with. fu() isn't making it worse, but it is completely missing the substitutions that actually make it simpler

>>> pprint(result[X])
                                     2                                                                                                                                    ↪
                                  sin (t)⋅cos(2⋅t)                                                                    2⋅sin(t)⋅sin(2⋅t)⋅cos(t)                            ↪
───────────────────────────────────────────────────────────────────────────────────── - ───────────────────────────────────────────────────────────────────────────────── ↪
   2                               2           2                               2           2                               2           2                               2  ↪
sin (t) - 3⋅sin(t)⋅sin(2⋅t) + 2⋅sin (2⋅t) + cos (t) - 3⋅cos(t)⋅cos(2⋅t) + 2⋅cos (2⋅t)   sin (t) - 3⋅sin(t)⋅sin(2⋅t) + 2⋅sin (2⋅t) + cos (t) - 3⋅cos(t)⋅cos(2⋅t) + 2⋅cos ( ↪

↪                                                                                                                                      2                                  ↪
↪                                      sin(t)⋅sin(2⋅t)⋅cos(2⋅t)                                                                   2⋅sin (2⋅t)⋅cos(t)                      ↪
↪ ──── - ───────────────────────────────────────────────────────────────────────────────────── + ──────────────────────────────────────────────────────────────────────── ↪
↪           2                               2           2                               2           2                               2           2                         ↪
↪ 2⋅t)   sin (t) - 3⋅sin(t)⋅sin(2⋅t) + 2⋅sin (2⋅t) + cos (t) - 3⋅cos(t)⋅cos(2⋅t) + 2⋅cos (2⋅t)   sin (t) - 3⋅sin(t)⋅sin(2⋅t) + 2⋅sin (2⋅t) + cos (t) - 3⋅cos(t)⋅cos(2⋅t)  ↪

↪                                                      2                                                                                              2                   ↪
↪                                                   cos (t)⋅cos(2⋅t)                                                                        cos(t)⋅cos (2⋅t)              ↪
↪ ───────────── - ───────────────────────────────────────────────────────────────────────────────────── + ─────────────────────────────────────────────────────────────── ↪
↪        2           2                               2           2                               2           2                               2           2                ↪
↪ + 2⋅cos (2⋅t)   sin (t) - 3⋅sin(t)⋅sin(2⋅t) + 2⋅sin (2⋅t) + cos (t) - 3⋅cos(t)⋅cos(2⋅t) + 2⋅cos (2⋅t)   sin (t) - 3⋅sin(t)⋅sin(2⋅t) + 2⋅sin (2⋅t) + cos (t) - 3⋅cos(t)⋅ ↪

↪
↪
↪ ──────────────────────
↪                 2
↪ cos(2⋅t) + 2⋅cos (2⋅t)

[oscarbenjamin]

Here is one possibility:

In [19]: result[X].simplify().expand(trig=True).factor().trigsimp().factor()
Out[19]: 
2⋅cos(t) + cos(2⋅t)
───────────────────
         3

[neozhaoliang]

This is awesome, but quite unintuitive!

[asmeurer]

The intuition behind this solution is that expand(trig=True) expands the double, triple, etc. trig identities (this is what it's documented to do), and trigsimp() by default applies the reverse identities (and to be fair the fact that that's one thing it does is somewhat arbitrary and might change in the future).