Ask Your Question
0

How to find all roots with solve?

asked 2019-11-18 18:27:26 -0500

Tofi gravatar image

updated 2019-11-18 18:30:04 -0500

I have the following quantity:

c_p = 1 - (4*sin(th)^2+2*5*sin(th)/pi+(5/(2*pi))^2)

and I am trying to find its roots (the solutions for c_p==0).

When plotting c_p as a function of th between $-\pi$ and $\pi$, I can see the curve crosses the x-axis at four positions. However, solve(c_p==0, th) is only giving two roots: $$\left[theta = -\arcsin\left(\frac{5}{4 \ \pi} + \frac{1}{2}\right),\quad theta = -\arcsin\left(\frac{5}{4 \ \pi} - \frac{1}{2}\right)\right].$$ It appears that solve can only find the roots that are in the domain of the $\arcsin$ function, i.e in the interval $[-\dfrac{\pi}{2},\dfrac{\pi}{2}]$. Is there a way to get the other two roots?

edit retag flag offensive close merge delete

1 answer

Sort by ยป oldest newest most voted
2

answered 2019-11-19 13:40:05 -0500

Emmanuel Charpentier gravatar image

updated 2019-11-20 08:20:42 -0500

Sympy offers us a solution, but sage is not (yet) ready to parse it totally. But with a little Sympy's doc reading and a little manual unscrewing, one can get it:

sage: [u.args[0].args[1]._sage_() for u in c_p.solve(theta, algorithm="sympy")]
[2*pi*n + arctan((2*pi - 5)/(pi*sqrt(20/pi - 25/pi^2 + 12))),
 pi + 2*pi*n - arctan((2*pi - 5)/(pi*sqrt(20/pi - 25/pi^2 + 12))),
 -pi + 2*pi*n + arctan(1/2*(2*pi + 5)/(pi*sqrt(-5/pi - 25/4/pi^2 + 3))),
 2*pi + 2*pi*n - arctan(1/2*(2*pi + 5)/(pi*sqrt(-5/pi - 25/4/pi^2 + 3)))]

A bit of fooling around with latex and string splicing allows us to get the $\LaTeX$ form:

$$2 \pi n + \arctan\left(\frac{2 \pi - 5}{\pi \sqrt{\frac{20}{\pi} - \frac{25}{\pi^{2}} + 12}}\right)$$ $$\pi + 2 \pi n - \arctan\left(\frac{2 \pi - 5}{\pi \sqrt{\frac{20}{\pi} - \frac{25}{\pi^{2}} + 12}}\right)$$ $$-\pi + 2 \pi n + \arctan\left(\frac{2 \pi + 5}{2 \pi \sqrt{-\frac{5}{\pi} - \frac{25}{4 \pi^{2}} + 3}}\right)$$ $$2 \pi + 2 \pi n - \arctan\left(\frac{2 \pi + 5}{2 \pi \sqrt{-\frac{5}{\pi} - \frac{25}{4 \pi^{2}} + 3}}\right)$$

Note: Maxima could also get it (via to_poly_solve), but gives pompous pedantic nonsense involving logs of complex expressions, probably due to the domain:complex default...

Note 2: Substituting these values in c_p gets rid of the trig functions easily (via trig_reduce() after declaring n integer...). However, the resulting expressions involve sums ad products of scalar and square roots values ; proving them to be 0 is another possibly nontrivial problem... It can be "checked" numerically, and formally proved for the last two:

sage: L2=[u.args[0].args[1]._sage_() for u in c_p.solve(theta, algorithm="sympy")]
sage: [c_p.subs(theta==u).trig_reduce().n() for u in L2]
[2.22044604925031e-16,
 2.22044604925031e-16,
 2.22044604925031e-16,
 2.22044604925031e-16]

sage: [bool(c_p.subs(theta==u).trig_reduce()==0) for u in L2]
[False, False, True, True]

EDIT : A few notes for the terminally curious:

fricas and giac give the same results as maxima (i. e. Sage).

I have been unable to obtain no-nonsensical results from to_poly_solve.

One can also screw around Mathematica's answers to get another form of the results (closer to what Sage/Maxima gives) :

sage: var("k", domain="integer")
k
sage: LM=[u[1][2][1].sage(locals={"C":lambda x:k}) for u in mathematica.Solve(c_p
....: ==0, theta)]
sage: LM
[pi + 2*pi*k + arcsin(1/4*(2*pi + 5)/pi),
 2*pi*k - arcsin(1/4*(2*pi + 5)/pi),
 pi + 2*pi*k - arcsin(1/4*(2*pi - 5)/pi),
 2*pi*k + arcsin(1/4*(2*pi - 5)/pi)]

which is:

$$\pi + 2 \pi k + \arcsin\left(\frac{2 \pi + 5}{4 \pi}\right)$$ $$2 \pi k - \arcsin\left(\frac{2 \pi + 5}{4 \pi}\right)$$ $$\pi + 2 \pi k - \arcsin\left(\frac{2 \pi - 5}{4 \pi}\right)$$ $$2 \pi k + \arcsin\left(\frac{2 \pi - 5}{4 \pi}\right)$$

which turns out to be somewhat easier to check/prove:

sage: [c_p.subs(theta==u).trig_reduce().expand() for u in LM]
[0, 0, 0, 0]
edit flag offensive delete link more

Comments

Ticket ?

tmonteil gravatar imagetmonteil ( 2019-11-20 08:45:53 -0500 )edit

On what subject ? Getting all Sympy results in Sage u=is a large work, already undertaken notably by someone whose trac login is rws and whose name escapes me presently (getting old isn't fun...). The present occurrence is but a special case...

Do you suggest a "sympy usability" metaticket ?

The Mathematica interface has other problems, most notably the lack of correct Mathematica<-->Sage conversions. E. g., try to use arctan in a Mathematica call...

sage: mathematica.Integrate(arctan(x),x)
Integrate[Arctan[x], x]

This is true in the reverse direction also 'Sage complains of the lack of a function arcTan...)

Since Mathematica isn't free, I doubt that polishing its interface can be considered high-priority. Sympy's can, however...

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2019-11-20 11:27:07 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2019-11-18 18:27:26 -0500

Seen: 132 times

Last updated: Nov 20 '19