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)
On Thursday, June 6, 2019 at 6:13:23 PM UTC-5, pull_over93 wrote:
>
> 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?
>
--
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 [email protected].
To post to this group, send email to [email protected].
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/18a9fc5a-7569-4496-9f04-47c415e64f05%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
