Problems and errors in solve an equation

edit

This is a quartic polynomial in omega^2, so make the substitution and solve:

var('omega_squared')
solve(omega_nf_eq.subs({omega : sqrt(omega_squared)}), omega_squared)

Take plusminus the square roots of these solutions to get the solutions for omega.

Okay. Time to think a bit.

Let's try to solve this in SR, the symbolic ring. Let's get rid of the goddamn numeric constants :

sage: var("k_1","k_2")
(k_1, k_2)
sage: E1="-J_d^2*k_yy*k_zz*omega^4 + 0.382*J_d^2*k_yy*omega^6 + 0.382*J_d^2*k_zz
....: *omega^6 - 0.145924*J_d^2*omega^8 + J_d*k_phiphi*k_yy*k_zz*omega^2 - 0.382
....: *J_d*k_phiphi*k_yy*omega^4 - 0.382*J_d*k_phiphi*k_zz*omega^4 + 0.145924*J_
....: d*k_phiphi*omega^6 + J_d*k_thetatheta*k_yy*k_zz*omega^2 - 0.382*J_d*k_thet
....: atheta*k_yy*omega^4 - 0.382*J_d*k_thetatheta*k_zz*omega^4 + 0.145924*J_d*k
....: _thetatheta*omega^6 - J_d*k_yphi^2*k_zz*omega^2 + 0.382*J_d*k_yphi^2*omega
....: ^4 - J_d*k_yy*k_ztheta^2*omega^2 + 0.382*J_d*k_ztheta^2*omega^4 + J_p^2*Om
....: ega^2*k_yy*k_zz*omega^2 - 0.382*J_p^2*Omega^2*k_yy*omega^4 - 0.382*J_p^2*O
....: mega^2*k_zz*omega^4 + 0.145924*J_p^2*Omega^2*omega^6 - k_phiphi*k_thetathe
....: ta*k_yy*k_zz + 0.382*k_phiphi*k_thetatheta*k_yy*omega^2 + 0.382*k_phiphi*k
....: _thetatheta*k_zz*omega^2 - 0.145924*k_phiphi*k_thetatheta*omega^4 + k_phip
....: hi*k_yy*k_ztheta^2 - 0.382*k_phiphi*k_ztheta^2*omega^2 + k_thetatheta*k_yp
....: hi^2*k_zz - 0.382*k_thetatheta*k_yphi^2*omega^2 - k_yphi^2*k_ztheta^2"
sage: E2=SR(E1.replace("0.382","k_1").replace("0.145924","k_2"))

What's this ?

sage: E2.expand()
J_p^2*Omega^2*k_2*omega^6 - J_d^2*k_2*omega^8 - J_p^2*Omega^2*k_1*k_yy*omega^4 - J_p^2*Omega^2*k_1*k_zz*omega^4 + J_d^2*k_1*k_yy*omega^6 + J_d^2*k_1*k_zz*omega^6 + J_d*k_2*k_phiphi*omega^6 + J_d*k_2*k_thetatheta*omega^6 + J_p^2*Omega^2*k_yy*k_zz*omega^2 + J_d*k_1*k_yphi^2*omega^4 - J_d*k_1*k_phiphi*k_yy*omega^4 - J_d*k_1*k_thetatheta*k_yy*omega^4 + J_d*k_1*k_ztheta^2*omega^4 - J_d*k_1*k_phiphi*k_zz*omega^4 - J_d*k_1*k_thetatheta*k_zz*omega^4 - J_d^2*k_yy*k_zz*omega^4 - k_2*k_phiphi*k_thetatheta*omega^4 - k_1*k_thetatheta*k_yphi^2*omega^2 + k_1*k_phiphi*k_thetatheta*k_yy*omega^2 - k_1*k_phiphi*k_ztheta^2*omega^2 - J_d*k_yy*k_ztheta^2*omega^2 + k_1*k_phiphi*k_thetatheta*k_zz*omega^2 - J_d*k_yphi^2*k_zz*omega^2 + J_d*k_phiphi*k_yy*k_zz*omega^2 + J_d*k_thetatheta*k_yy*k_zz*omega^2 - k_yphi^2*k_ztheta^2 + k_phiphi*k_yy*k_ztheta^2 + k_thetatheta*k_yphi^2*k_zz - k_phiphi*k_thetatheta*k_yy*k_zz

Looks like a polynomial in omega. Let's check this : sage: LK=E2.expand().coefficients(omega); LK

[[-k_yphi^2*k_ztheta^2 + k_phiphi*k_yy*k_ztheta^2 + k_thetatheta*k_yphi^2*k_zz - k_phiphi*k_thetatheta*k_yy*k_zz,
  0],
 [J_p^2*Omega^2*k_yy*k_zz - k_1*k_thetatheta*k_yphi^2 + k_1*k_phiphi*k_thetatheta*k_yy - k_1*k_phiphi*k_ztheta^2 - J_d*k_yy*k_ztheta^2 + k_1*k_phiphi*k_thetatheta*k_zz - J_d*k_yphi^2*k_zz + J_d*k_phiphi*k_yy*k_zz + J_d*k_thetatheta*k_yy*k_zz,
  2],
 [-J_p^2*Omega^2*k_1*k_yy - J_p^2*Omega^2*k_1*k_zz + J_d*k_1*k_yphi^2 - J_d*k_1*k_phiphi*k_yy - J_d*k_1*k_thetatheta*k_yy + J_d*k_1*k_ztheta^2 - J_d*k_1*k_phiphi*k_zz - J_d*k_1*k_thetatheta*k_zz - J_d^2*k_yy*k_zz - k_2*k_phiphi*k_thetatheta,
  4],
 [J_p^2*Omega^2*k_2 + J_d^2*k_1*k_yy + J_d^2*k_1*k_zz + J_d*k_2*k_phiphi + J_d*k_2*k_thetatheta,
  6],
 [-J_d^2*k_2, 8]]

Hey ! it's indeed a polynomial of degree 4 in omega^2... Let's give its coefficients a more sy(n|mpa)thetic name (i. e. create a simpler-to-write polynomial having the same value given the right substitution) and check it :

sage: SK=[var("K_{}".format(u[1]))==u[0] for u in LK]
sage: E3=sum([SK[u].lhs()*omega^(2*u) for u in range(len(SK))])
sage: bool(E3.subs(SK)==E2)
True

Transform this in a polynomial in omega^2 :

sage: E4=E3.subs([omega^(2*u)==OmSq^u for u in range(len(SK))]);E4
K_8*OmSq^4 + K_6*OmSq^3 + K_4*OmSq^2 + K_2*OmSq + K_0

Now, this fourth-degree polynomial has four (possibly equal, possibly complex) roots:

sage: Sol4=[u.subs(SK) for u in solve(E4, OmSq)]
sage: len(Sol4)
4

Each of these solutions s for OmSq gives two possible solutions for omega (namely $-\sqrt{s}$ and $\sqrt{s}$) :

sage: Sol3=flatten([[omega==-sqrt(s.rhs()),omega==sqrt(s.rhs())] for s in Sol4])
sage: len(Sol3)
8

If one insists to reincorporate the goddamn numeric constants in the final solution, this is easy :

sage: Sol=[u.subs([k_1==0.382, k_2==0.145924]) for u in Sol3]

Printing the solutions is left as an exercise for the (masochistic) reader (fitted with an A0 printer and lots of time). Similarly, checking them (by substitution in E2) is quite slow...

An alternative (=more general (=better)) solution is to create a multivariate polynomial ring with the right unknowns (please choose easier names !) and compute the solutions by searching the "right" ideal. This would also allow you to create the set of roots without having to compute an explicit solution for them (which do not, in general, exists for degrees >=5).

HTH,