# Problems and errors in solve an equation

Hi everybody, I want to solve this non linear equation: omega_nf_eq = 0.

m,J_d,J_p,y,Y,omega,Omega,phi,Phi,z,Z,theta,Theta,k_yy,k_zz,k_phiphi,k_yphi,k_ztheta,k_thetatheta,plane_xy1,plane_xy2,plane_xz1,plane_xz2 = var('m J_d J_p y Y omega Omega phi Phi z Z theta Theta k_yy k_zz k_phiphi k_yphi k_ztheta k_thetatheta plane_xy1 plane_xy2 plane_xz1 plane_xz2')
t = var('t')

omega_nf_eq = -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_thetatheta*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*Omega^2*k_yy*k_zz*omega^2 - 0.382*J_p^2*Omega^2*k_yy*omega^4 - 0.382*J_p^2*Omega^2*k_zz*omega^4 + 0.145924*J_p^2*Omega^2*omega^6 - k_phiphi*k_thetatheta*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_phiphi*k_yy*k_ztheta^2 - 0.382*k_phiphi*k_ztheta^2*omega^2 + k_thetatheta*k_yphi^2*k_zz - 0.382*k_thetatheta*k_yphi^2*omega^2 - k_yphi^2*k_ztheta^2 == 0
solve(omega_nf_eq, omega)


But the Sage is unable to find the solution, damn it. I get this error message:

TypeError: ECL says: Memory limit reached. Please jump to an outer pointer, quit program and enlarge the memory limits before executing the program again.


I know the equation is big but i didn't expected such many problems. I've already tried sympy but nothing. Is there another way?

edit retag close merge delete

Sort by ยป oldest newest most voted

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,

more

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.

more