Python同情不能解决多项式函数

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.

通过使用sympy求解多项式函数，我遇到了问题。以下示例显示了一个案例，该案例提供了我无法管理的错误消息。如果多项式变得更简单，则求解器正常工作。请复制并粘贴代码以检查系统上的错误。

 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)

2 个解决方案

#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.

正如我在评论中提到的，我不熟悉同情，但这里是如何使用任意精度mpmath模块找到等式的根。

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...

为了避免Python浮点数的精度问题，通常以字符串形式将浮点数传递给mpmath，除非从整数构造它们很方便。我想这对你的方程式来说并不是一个问题，因为你的系数精度相当低，但无论如何......

Here's your equation translated directly into mpmath syntax:

这是你的等式直接翻译成mpmath语法：

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.

mpf是任意精度的浮点构造函数，mpc是复数构造函数。有关如何调用这些构造函数的信息，请参阅mpmath文档。

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.

让你成为f（x）的根。然后令f（x）= g（x）（x-u），g（x）的任何根也是f（x）的根。通过使用for循环将每个找到的根保存到列表然后从前一个函数构建一个新函数，将这个新函数存储在另一个列表中，我们可以方便地多次执行此操作。

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.

在这个版本中，我使用“muller”解算器，因为在寻找复杂的根时建议使用它，但它实际上给出了与使用默认的“割线”求解器相同的答案。

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.

我们将i指定为默认关键字参数，以便使用我在定义函数时所具有的值。否则，在调用函数时使用i的当前值，这不是我们想要的。

---

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.

这种用于查找多个根的技术可用于任意函数。但是，当我们想要求解多项式时，mpmath有一种更好的方法可以同时找到所有的根：多根函数。我们只需要给它一个多项式系数的列表（或元组）。

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.

如您所见，结果与先前技术获得的结果非常接近。使用多根不仅更准确，它具有很大的优势，我们不需要为根指定起始近似，mpmath足够聪明，可以为自己创建一个。

#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.

感谢PM2Ring的帮助，我创建了一个完整的代码示例，它提取任意给定的四阶sympy多项式的系数并求解它。

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)