Python sympy 无法解决多项式函数 [英] Polynomial function cannot be solved by Python sympy
问题描述
我在用 sympy 求解多项式函数时遇到了问题.以下示例显示了一个给出我无法管理的错误消息的案例.如果多项式变得更简单，则求解器可以正常工作.请复制并粘贴代码以检查系统上的错误.
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.
这是错误信息:
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)
推荐答案
正如我在评论中提到的，我不熟悉 sympy，但这里是如何使用任意精度 mpmath 模块.
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.
为了避免 Python 浮点数的精度问题，通常以字符串形式将浮点数传递给 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...
这是您的等式直接转换为 mpmath 语法:
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 是任意精度浮点构造函数，mpc 是复数构造函数.有关如何调用这些构造函数的信息，请参阅 mpmath 文档.
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.
然而，我们不需要像这样处理 I:我们可以直接将系数定义为复数.
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))
输出
(1.35717989161653180236360985534 - 0.202974596285109153971324419197j)
test: (-3.2311742677852643549664402034e-26 + 6.4623485355705287099328804068e-27j)
但是我们怎样才能找到其他的根呢?简单的！
But how can we find the other roots? Simple!
让 u 成为 f(x) 的根.然后让 f(x) = g(x)(x - u) 并且 g(x) 的任何根也是 f(x) 的根.我们可以通过使用 for 循环方便地多次执行此操作，该循环将每个找到的根保存到一个列表中，然后从前一个函数构建一个新函数，将这个新函数存储在另一个列表中.
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.
在这个版本中，我使用muller"求解器，这是在寻找复杂根时推荐的，但它实际上给出了与使用默认割线"求解器相同的答案.
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]))
输出
(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)
在
lambda x,i=i: funcs[i](x) / (x - roots[i])
我们指定 i 作为默认关键字参数，以便使用定义函数时 i 的值.否则，在调用函数时会使用i的current值，这不是我们想要的.
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.
这种寻找多个根的技术可用于任意函数.但是，当我们想要求解多项式时，mpmath 有一个更好的方法可以同时找到所有根:polyroots 函数.我们只需要给它一个多项式系数的列表(或元组).
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))
输出
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)
如您所见，结果与通过先前技术获得的结果非常接近.使用 polyroots 不仅更准确，它还有一个很大的优势，我们不需要为根指定一个起始近似值，mpmath 足够聪明，可以为自己创建一个.
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.
这篇关于Python sympy 无法解决多项式函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
