Polynomial function cannot be solved by Python sympy

Question

I have problems by solving a polynomial function with sympy. The following example shows a case which gives an error message that I cannot manage. If the polynomial gets simpler the solver works properly. Please copy and paste the code to check the error on your system as well.

 import sympy
 from sympy import I
 omega = sympy.symbols('omega')

 def function(omega):
    return - 0.34*omega**4      \
           + 7.44*omega**3      \
           + 4.51*I*omega**3    \
           + 87705.64*omega**2  \
           - 53.08*I*omega**2   \
           - 144140.03*omega    \
           - 22959.95*I*omega   \
           + 42357.18 + 50317.77*I
 sympy.solve(function(omega), omega)

Do you know how I can achieve a result? Thank you for your help.

EDIT:

This is the error message:

---

TypeError                                 Traceback (most recent call last)
<ipython-input-7-512446a62fa9> in <module>()
      1 def function(omega):
      2     return - 0.34*omega**4                 + 7.44*omega**3                 + 4.51*I*omega**3               + 87705.64*omega**2             - 53.08*I*omega**2              - 144140.03*omega               - 22959.95*I*omega              + 42357.18 + 50317.77*I
----> 3 sympy.solve(function(omega), omega)

C:\Anaconda\lib\site-packages\sympy\solvers\solvers.pyc in solve(f, *symbols, **flags)
   1123     # restore floats
   1124     if floats and solution and flags.get('rational', None) is None:
-> 1125         solution = nfloat(solution, exponent=False)
   1126 
   1127     if check and solution:  # assumption checking

C:\Anaconda\lib\site-packages\sympy\core\function.pyc in nfloat(expr, n, exponent)
   2463             return type(expr)([(k, nfloat(v, n, exponent)) for k, v in
   2464                                list(expr.items())])
-> 2465         return type(expr)([nfloat(a, n, exponent) for a in expr])
   2466     rv = sympify(expr)
   2467 

C:\Anaconda\lib\site-packages\sympy\core\function.pyc in nfloat(expr, n, exponent)
   2497     return rv.xreplace(Transform(
   2498         lambda x: x.func(*nfloat(x.args, n, exponent)),
-> 2499         lambda x: isinstance(x, Function)))
   2500 
   2501 

C:\Anaconda\lib\site-packages\sympy\core\basic.pyc in xreplace(self, rule)
   1085 
   1086         """
-> 1087         value, _ = self._xreplace(rule)
   1088         return value
   1089 

C:\Anaconda\lib\site-packages\sympy\core\basic.pyc in _xreplace(self, rule)
   1093         """
   1094         if self in rule:
-> 1095             return rule[self], True
   1096         elif rule:
   1097             args = []

C:\Anaconda\lib\site-packages\sympy\core\rules.pyc in __getitem__(self, key)
     57     def __getitem__(self, key):
     58         if self._filter(key):
---> 59             return self._transform(key)
     60         else:
     61             raise KeyError(key)

C:\Anaconda\lib\site-packages\sympy\core\function.pyc in <lambda>(x)
   2496 
   2497     return rv.xreplace(Transform(
-> 2498         lambda x: x.func(*nfloat(x.args, n, exponent)),
   2499         lambda x: isinstance(x, Function)))
   2500 

C:\Anaconda\lib\site-packages\sympy\core\function.pyc in nfloat(expr, n, exponent)
   2463             return type(expr)([(k, nfloat(v, n, exponent)) for k, v in
   2464                                list(expr.items())])
-> 2465         return type(expr)([nfloat(a, n, exponent) for a in expr])
   2466     rv = sympify(expr)
   2467 

C:\Anaconda\lib\site-packages\sympy\core\function.pyc in nfloat(expr, n, exponent)
   2463             return type(expr)([(k, nfloat(v, n, exponent)) for k, v in
   2464                                list(expr.items())])
-> 2465         return type(expr)([nfloat(a, n, exponent) for a in expr])
   2466     rv = sympify(expr)
   2467 

TypeError: __new__() takes exactly 3 arguments (2 given)

Answer 1

As I mentioned in the comments I'm not familiar with sympy, but here's how to find the roots of your equation using the arbitrary-precision mpmath module.

In order to avoid precision issues with Python floats, it's usual to pass floats to mpmath in string form, unless it's convenient to construct them from integers. I guess it's not really an issue with your equation, since your coefficients have rather low precision, but anyway...

Here's your equation translated directly into mpmath syntax:

from mpmath import mp

I = mp.mpc(0, 1)

def func(x):
    return (-mp.mpf('0.34') * x ** 4
        + mp.mpf('7.44') * x ** 3
        + mp.mpf('4.51') * I * x ** 3
        + mp.mpf('87705.64') * x ** 2
        - mp.mpf('53.08') * I * x ** 2
        - mp.mpf('144140.03') * x
        - mp.mpf('22959.95') * I * x
        + mp.mpf('42357.18') + mp.mpf('50317.77') * I)

mpf is the arbitrary-precision float constructor, mpc is the complex number constructor. Please see the mpmath docs for info on how to call these constructors.

However, we don't need to mess about with I like that: we can just define the coefficients directly as complex numbers.

from __future__ import print_function
from mpmath import mp

# set precision to ~30 decimal places
mp.dps = 30

def func(x):
    return (mp.mpf('-0.34') * x ** 4
        + mp.mpc('7.44', '4.51') * x ** 3
        + mp.mpc('87705.64', '-53.08') * x ** 2
        + mp.mpc('-144140.03', '-22959.95') * x
        + mp.mpc('42357.18', '50317.77'))

x = mp.findroot(func, 1)
print(x)
print('test:', func(x))

output

(1.35717989161653180236360985534 - 0.202974596285109153971324419197j)
test: (-3.2311742677852643549664402034e-26 + 6.4623485355705287099328804068e-27j)

But how can we find the other roots? Simple!

Let u be a root of f(x). Then let f(x) = g(x)(x - u) and any root of g(x) is also a root of f(x). We can conveniently do this multiple times by using a for loop that saves each found root to a list and then builds a new function from the previous function, storing this new function in another list.

In this version, I use the "muller" solver, as that's recommended when looking for complex roots, but it actually gives the same answers as using the default "secant" solver.

from __future__ import print_function
from mpmath import mp

# set precision to ~30 decimal places
mp.dps = 30

def func(x):
    return (mp.mpf('-0.34') * x ** 4
        + mp.mpc('7.44', '4.51') * x ** 3
        + mp.mpc('87705.64', '-53.08') * x ** 2
        + mp.mpc('-144140.03', '-22959.95') * x
        + mp.mpc('42357.18', '50317.77'))

x = mp.findroot(func, 1)
print(x)
print('test:', func(x))

funcs = [func]
roots = []

#Find all 4 roots
for i in range(4):
    x = mp.findroot(funcs[i], 1, solver="muller")
    print('%d: %s' % (i, x))
    print('test: %s\n%s\n' % (funcs[i](x), funcs[0](x)))
    roots.append(x)
    funcs.append(lambda x,i=i: funcs[i](x) / (x - roots[i]))    

output

(1.35717989161653180236360985534 - 0.202974596285109153971324419197j)
test: (-3.2311742677852643549664402034e-26 + 6.4623485355705287099328804068e-27j)
0: (1.35717989161653180236360985534 - 0.202974596285109153971324419197j)
test: (-3.2311742677852643549664402034e-26 + 6.4623485355705287099328804068e-27j)
(-3.2311742677852643549664402034e-26 + 6.4623485355705287099328804068e-27j)

1: (0.2859569222439674364374376897 + 0.465618100581394597267702975072j)
test: (2.70967991111831205485691272044e-27 - 4.34146435347156317045282996313e-27j)
(0.0 + 6.4623485355705287099328804068e-27j)

2: (-497.86487129703641182172086688 + 6.49836193448774263077718499855j)
test: (1.25428695883356196194609388771e-26 - 3.46609896266051486795576778151e-28j)
(3.11655159180984988723836070362e-21 - 1.65830325771275337225587644119e-22j)

3: (518.104087424352383171155113452 + 6.50370044356891310239702468087j)
test: (-4.82755073209873082528016484574e-30 + 7.38353321095804877623117526215e-31j)
(-1.31713649147437845007238988587e-21 + 1.68350641700147843422461467477e-22j)

In

lambda x,i=i: funcs[i](x) / (x - roots[i])

we specify i as a default keyword argument, so that the value that i had when the function was defined is used. Otherwise, the current value of i is used when the function is called, which is not what we want.

---

This technique for finding multiple roots can be used for arbitrary functions. However, when we want to solve polynomials, mpmath has a better way that can find all roots simultaneously: the polyroots function. We just need to give it a list (or tuple) of the polynomial's coefficients.

from __future__ import print_function
from mpmath import mp

# set precision to ~30 decimal places
mp.dps = 30

coeff = (
    mp.mpf('-0.34'),
    mp.mpc('7.44', '4.51'),
    mp.mpc('87705.64', '-53.08'),
    mp.mpc('-144140.03', '-22959.95'),
    mp.mpc('42357.18', '50317.77'),
)

roots = mp.polyroots(coeff)
for i, x in enumerate(roots):
    print('%d: %s' % (i, x))
    print('test: %s\n' % mp.polyval(coeff, x))

output

0: (1.35717989161653180236360985534 - 0.202974596285109153971324419197j)
test: (6.4623485355705287099328804068e-27 - 6.4623485355705287099328804068e-27j)

1: (0.2859569222439674364374376897 + 0.465618100581394597267702975072j)
test: (0.0 + 0.0j)

2: (-497.86487129703641182172086688 + 6.49836193448774263077718499855j)
test: (-2.27689218423463552142807161949e-21 + 7.09242751778865525915133624646e-23j)

3: (518.104087424352383171155113452 + 6.50370044356891310239702468087j)
test: (-7.83663157514495734538720675411e-22 - 1.08373584941517766465574404422e-23j)

As you can see, the results are very close to those obtained by the previous technique. Using polyroots is not only more accurate, it has the big advantage that we don't need to specify a starting approximation for the root, mpmath is smart enough to create one for itself.