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)

Nunhood

"gives an error message that I cannot manage" and what is that error message? – Ballyhoo

Welcome to SO! You should add the error message to your question so people with the same problem can google it and find your question which hopefully gets solved. – Casillas

I added the final line of the error message. Hopefully that helps. – Nunhood

FWIW, a root of your function, found using mpmath.findroot is ~= 1.35717989161653180236360985534 - 0.202974596285109153971324419197i – Grondin

@PM2Ring: Which method and starting point did you use for mpmath.findroot? – Nunhood

I just used the default solver=Secant, with 1 as the starting approximation. But I converted all of the float constants to strings before passing them to mp.mpf(). I can post my code as an answer if you like. – Grondin

@PM2Ring: Thank you for your answer. Unfortunately, I get only one root with your workaround, but there should be four roots. If I change the starting point I can find the other roots as well. However, this procedure is not convenient for an automation... – Nunhood

@Aloex: You should've mentioned earlier that you want all 4 roots. :) There's a fairly simple way to do that; I'll post an answer that illustrates how. – Grondin

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.

Grondin

Thanks a lot! ... sry, that I missed to mention the other three roots. :) Finally, I putted your code into my script and got a nice workaround. I will post a full example if you come along with a sympy function at the beginning. – Nunhood

@Aloex: I'm glad my answer has helped you, even though it's not using sympy. Please consider accepting it. – Grondin

You can also use sympy.nsolve, which is a wrapper around mpmath.findroot. – Chemotaxis

Answer 2

Thanks to the help of PM2Ring, I created a complete code example that extracts the coefficients of an arbitrary given fourth order sympy polynomial and solves it.

import sympy
from sympy import I
import mpmath as mp
import numpy as np

omega = sympy.symbols('omega')

# Define polynomial in sympy
def function(omega):
    return ( -  0.34                   *omega**4      \
             + (7.44 + 4.51*I)         *omega**3      \
             + (87705.64 - 53.08*I)    *omega**2  \
             - (144140.03 + 22959.95*I)*omega \
             + 42357.18 + 50317.77*I)

# Show defined polynomial
print''+ str(function(omega).expand()) +'\n'
result = sympy.Poly(function(omega).expand(), omega)

# Obtain coefficients of the polynomial
coeff = (
    mp.mpc(str(np.real(sympy.lambdify( (I),result.coeffs()[0])(1j))) , str(np.imag(sympy.lambdify( (I),result.coeffs()[0])(1j)))),
    mp.mpc(str(np.real(sympy.lambdify( (I),result.coeffs()[1])(1j))) , str(np.imag(sympy.lambdify( (I),result.coeffs()[1])(1j)))),
    mp.mpc(str(np.real(sympy.lambdify( (I),result.coeffs()[2])(1j))) , str(np.imag(sympy.lambdify( (I),result.coeffs()[2])(1j)))),
    mp.mpc(str(np.real(sympy.lambdify( (I),result.coeffs()[3])(1j))) , str(np.imag(sympy.lambdify( (I),result.coeffs()[3])(1j)))),
    mp.mpc(str(np.real(sympy.lambdify( (I),result.coeffs()[4])(1j))) , str(np.imag(sympy.lambdify( (I),result.coeffs()[4])(1j))) ),
  )

# Calculate roots of the polynomial
roots = mp.polyroots(coeff)
for no, frequency in enumerate(roots):
    print frequency

Output

-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

(1.35717989161653 - 0.202974596285109j)
(0.285956922243967 + 0.465618100581395j)
(-497.864871297036 + 6.49836193448774j)
(518.104087424352 + 6.50370044356891j)

Nunhood

Answer 3

import sympy as *

x = symbols('x')

F = - 0.34*x**4+ 7.44*x**3+ 4.51*I*x**3+ 87705.64*x**2- 53.08*I*x**2- 144140.03*x- 22959.95*I*x+ 42357.18 + 50317.77*I

solve(F,x,dict = True)

[{x: -497.864871297036 + 6.49836193448774*I}, {x: 0.285956922243967 + 0.465618100581395*I}, {x: 1.35717989161653 - 0.202974596285109*I}, {x: 518.104087424352 + 6.50370044356891*I}]

Neu

Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center. – Riotous