Ask Your Question
2

Problems and errors in solve an equation

asked 2019-06-06 17:38:46 -0600

pull_over93 gravatar image

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 flag offensive close merge delete

2 answers

Sort by ยป oldest newest most voted
1

answered 2019-06-07 04:52:42 -0600

Emmanuel Charpentier gravatar image

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,

edit flag offensive delete link more
1

answered 2019-06-07 03:24:14 -0600

rburing gravatar image

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.

edit flag offensive delete link more

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-06-06 17:38:46 -0600

Seen: 57 times

Last updated: Jun 07