Problem in solving non linear equation equation

[pull_over93]

Hi everybody, I'm tring to solve this equation without succes: omega_nf_eq = 0

import sympy as sym    

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 = sym.symbols('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', real  = True) 

t, omega_nf = sym.symbols('t, omega_nf', real = True)

omega_nf_eq = sym.Eq(omega_nf, -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)

solution = sym.solve(omega_nf_eq.rhs, omega, dict = True , force=True, manual=True, set=True)

Even after half an hour, sympy is unable to give a solution.

I also tried sage math, but I had the same failure.

Suggestions?

[Oscar Benjamin]

SymPy will find the solution eventually I think.

This is a polynomial of order 8 having only even powers so with a
substitution omega**2 -> x it's a quartic in x with complicated
symbolic coefficients. The general formula for a quartic is horrendous
and in this case your coefficients are already large expressions.

What this means is that the solution will take a while to compute and
will probably be too complicated to be useful when done anyway. Really
complicated analytic expressions can be difficult to use without
substituting numeric values in and it's much easier to find the roots
if you substitute the numbers in first.

I'm not sure why it is so slow though. You can get to the solutions
for omega^2 like this:

In [39]: coeffs = Poly(omega_nf_eq.rhs.subs(omega, sqrt(x)),
x).coeffs()

In [40]: a, b, c, d, e = symbols('a, b, c, d, e')

In [41]: p = Poly(a*x**4+b*x**3+c*x**2+d*x+e, x)

In [42]: rs = solve(p, x)

Print the first root:

In [43]: rs[0].subs(zip((a,b,c,d,e), coeffs))

(output omitted as it's long and complicated)


My suggestion is to consider why you need the solutions to this
equation. There may be a better way to reach your actual end goal then
getting the solutions in explicit analytic form.

--
Oscar

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/a4ebc67b-3e40-448b-b6a3-5500f6b9dab3%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Chris Smith]

As Oscar said, the equation is essentially a quartic with symbolic coefficients for which it is not possible to write a single generally valid solution (since the solutions depend on the relationship between the coefficients);

eq = Eq(omega_nf, C0 + C1*omega**2 + C2*omega**4 + C3*omega**6 + C4*omega**8)
d = Dict({C0: -k_phiphi*k_thetatheta*k_yy*k_zz + k_phiphi*k_yy*k_ztheta**2 + k_thetatheta*k_yphi**2*k_zz - k_yphi**2*k_ztheta**2,
C1: J_d*k_phiphi*k_yy*k_zz + J_d*k_thetatheta*k_yy*k_zz - J_d*k_yphi**2*k_zz - J_d*k_yy*k_ztheta**2 + J_p**2*Omega**2*k_yy*k_zz + 0.382*k_phiphi*k_thetatheta*k_yy + 0.382*k_phiphi*k_thetatheta*k_zz - 0.382*k_phiphi*k_ztheta**2 - 0.382*k_thetatheta*k_yphi**2,
C2: -J_d**2*k_yy*k_zz - 0.382*J_d*k_phiphi*k_yy - 0.382*J_d*k_phiphi*k_zz - 0.382*J_d*k_thetatheta*k_yy - 0.382*J_d*k_thetatheta*k_zz + 0.382*J_d*k_yphi**2 + 0.382*J_d*k_ztheta**2 - 0.382*J_p**2*Omega**2*k_yy - 0.382*J_p**2*Omega**2*k_zz - 0.145924*k_phiphi*k_thetatheta,
C3: 0.382*J_d**2*k_yy + 0.382*J_d**2*k_zz + 0.145924*J_d*k_phiphi + 0.145924*J_d*k_thetatheta + 0.145924*J_p**2*Omega**2,
C4: -0.145924*J_d**2, x: J_g1 + omega})

real_roots(eq.xreplace(d.subs({...}))) should very quickly give you the real solution when you put in for {...} the known values, e.g. `{J_d: 1, k_zz: 2, etc...}`. If the solution is written in terms of RootOf you can evaluate those to whateve precision is needed, e.g. `RootOf(x+x**3-4*x**4,1).n(3) -> 0.725`

BTW, I have been experimenting with a constant-extracting routine:

def cify(eq, v):
 from sympy.core.rules import Transform
 sym=numbered_symbols('C')
 def do(e, x):
 if e.is_Add:
 e = collect(e, x)
 c, fx = e.as_independent(x)
 s = next(sym)
 did[s] = c
 return e.func(s,fx)
 did={}
 e = None
 while e != eq:
 e = eq
 eq = eq.xreplace(Transform(lambda x: do(x, v), lambda x:(x.is_Add or x.is_Mul) and x.as_independent(v)[0].args))
 r, d = cse(did.values())
 did = dict(zip(did.keys(), d))
 r = dict(r)
 return r, did, eq
eq = omega_nf_eq.rhs
print(cify(eq,omega))
({
x0: k_ztheta**2,
x1: k_phiphi*x0,
x2: k_yphi**2,
x3: k_thetatheta*x2,
x4: k_phiphi*k_thetatheta,
x5: k_yy*k_zz,
x6: 0.382*k_yy,
x7: 0.382*k_zz,
x8: J_d*k_phiphi,
x9: J_d*k_thetatheta,
x10: J_d*x2,
x11: J_d*x0,
x12: J_p**2*Omega**2,
x13: 0.145924*k_phiphi,
x14: J_d**2,
x15: k_yy*x14
},
{
C0: k_yy*x1 + k_zz*x3 - x0*x2 - x4*x5,
C1: -k_yy*x11 - k_zz*x10 - 0.382*x1 + x12*x5 - 0.382*x3 + x4*x6 + x4*x7 + x5*x8 + x5*x9,
C2: -k_thetatheta*x13 - k_zz*x15 + 0.382*x10 + 0.382*x11 - x12*x6 - x12*x7 - x6*x8 - x6*x9 - x7*x8 - x7*x9,
C3: J_d*x13 + 0.145924*x12 + x14*x7 + 0.382*x15 + 0.145924*x9,
C4: -0.145924*x14
},
C0 + C1*omega**2 + C2*omega**4 + C3*omega**6 + C4*omega**8)

[pull_over93]

Thank you so much, for the answer. I have to get the solutions to plot them. The solutions represent the critical speed of a shaft and there's no other way.  Other students are using other software for their project like Maple or Mathematica, which they cause less problem. I was using python and I've written around 50 pages of code in jupyter notebook and i get stuck in this last calculation. Using Sagemath i've solved the equation and then i've copied into my python code.

> To unsubscribe from this group and stop receiving emails from it, send an email to sy...@googlegroups.com.

[pull_over93]

Thank you so much, for the answer. I have to try it . Really Thank you!