On Jul 29, 12:06 am, wflynny <wfly...@gmail.com> wrote:
> While finding eigenvalues for a matrix of mine, I ran into a
> polynomial decompose error. It seems that when a polynomial has both a
> cubic and quartic term and coefficients that consist of symbols,
> factortools.py runs into trouble. To simplify the error, here is some
> test code that throws the same error:
>
> --------------------------
> import sympy as sp
>
> x = sp.Symbol('x')
> z = sp.Symbol('z')
> g = cos(z)*x**4 + x**3 + 1
> res = sp.solve(g, x)
> print res
> --------------------------
>
I replied to Aaron but afterward tried your problem. There is a
problem with the current quartic_roots solver (Issue 1468) and the new
code is awaiting review. Also, the current distribution takes a long
time (>30 seconds for me) to simplify the results. It could be that
issue 1495 is affecting this. Running your problem through my patched
version of sympy, after "polifying" the original equation and running
the result through cse, gives the following:
>>> var('x z')
>>> r,p=polify(cos(z)*x**4 + x**3 + 1 )
>>> ans=solv(p,x)
>>> orig_ans=[tmp.subs(r) for tmp in ans];rr,ee=cse(orig_ans)
rr
>>> for ri in rr:
... print ri
...
(x0, cos(z))
(x1, x0**(-3))
(x2, x1/2)
(x3, x1**2)
(x4, x3/4)
(x5, 128*x2/27)
(x6, x4 - x5)
(x7, x6**(1/2))
(x8, x2 + x7)
(x9, x3**(1/6))
(x10, x8**(-1/3))
(x11, 1/x10)
(x12, x9**2)
(x13, 3**(1/2))
(x14, I*x13/2)
(x15, 1/2 + x14)
(x16, 1/2 - x14)
(x17, 1/x16)
(x18, x11**(-2))
(x19, 1/x15)
(x20, 16*x10*x9/3)
(x21, 4*x11)
(x22, x12 + x20 + x21)
(x23, x22**(1/2))
(x24, -x17*x20/4)
(x25, x19**2)
(x26, 16*x12*x18*x25/9)
(x27, 32*x12*x17*x18/9)
(x28, x17**2)
(x29, 16*x12*x18*x28/9)
(x30, 1/x18)
(x31, 4*x30)
(x32, 16*x12*x18/9)
(x33, 2*x19*x32)
(x34, 16*x17*x9/3)
(x35, -32*x9/3)
(x36, x19*x27)
(x37, -x19*x35/2)
(x38, x26 + x27 + x29 + x31 + x32 + x33 + x34 + x35 + x36 + x37)
(x39, x38**(1/2))
(x40, -x21/4)
(x41, x12/2)
(x42, -x19*x20/4)
(x43, x23/4)
(x44, -3*x35/128)
(x45, -x39)
(x46, -8*x43*x44)
(x47, x24 + x40 + x41 + x42 + x45 + x46)
(x48, x47**(1/2))
(x49, x48/2)
(x50, -x46)
(x51, x24 + x39 + x40 + x41 + x42 + x50)
(x52, x51**(1/2))
(x53, x52/2)
>>> for ei in ee:
... print ei
...
x43 + x49 - x44
x53 - x43 - x44
x43 - x44 - x49
-x43 - x44 - x53
>>>
The results can be verified at z=pi:
###
>>> var('x z');r,p=polify(cos(z)*x**4 + x**3 + 1 );assert x in
>>> p;ans=solv(p,x);anst=[tmp.subs(r) for tmp in
>>> ans];rr,ee=cse(anst);atpi=[tmp.subs(z,pi).evalf(n=4) for tmp in anst];print
>>> atpi;print [(-x**4+x**3+1).subs(x,tmp).evalf(n=4) for tmp in atpi]
(x, z)
[1.380, 2.194e-1 - 9.145e-1*I, -8.192e-1, 2.194e-1 + 9.145e-1*I]
[0, -3.205e-6 - 2.137e-6*I, 3.815e-6, -3.205e-6 + 2.137e-6*I]
>>>
###
So the roots satisfy the original equation. Using the current
distribution gives the following results (for z=pi) which are close,
but inconsistent with the original equation:
>>> var('x z');r,p=polify(cos(z)*x**4 + x**3 + 1 );assert x in
>>> p;ans=solve(p,x,simplified=False,simplify=False);anst=[tmp.subs(r) for tmp
>>> in ans];rr,ee=cse(anst);atpi=[tmp.subs(z,pi).evalf(n=4) for tmp in
>>> anst];print atpi;print [(-x**4+x**3+1).subs(x,tmp).evalf(n=4) for tmp in
>>> atpi]
(x, z)
[2.806e-1 + 8.976e-1*I, 2.806e-1 - 8.976e-1*I, -8.663e-1 - .0e-22*I,
1.305 + .0e-22*I]
[6.909e-2 + 2.211e-1*I, 6.909e-2 - 2.211e-1*I, -2.133e-1 + .0e-12*I,
3.214e-1 + .0e-9*I]
/c
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to sympy@googlegroups.com
To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sympy?hl=en
-~----------~----~----~----~------~----~------~--~---
