在 Python 中求解多项式方程答案

【问题标题】：Solving polynomial equations in Python在 Python 中求解多项式方程
【发布时间】：2012-05-04 23:11:45
【问题描述】：

到目前为止，我一直使用 Mathematica 来求解解析方程。但是现在我需要解决数百个这种类型的方程（特征多项式）

a_20*x^20+a_19*x^19+...+a_1*x+a_0=0 (constant floats a_0,...a_20)

这会在 Mathematica 中产生非常长的计算时间。

在 numpy 或任何其他包中是否有一个现成的命令来解决这种类型的方程？（到目前为止，我只将 Python 用于模拟，因此我对分析工具了解不多，并且在 numpy 教程中找不到任何有用的东西。

【问题讨论】：

为什么你认为python会比mathematica更快？

标签： python wolfram-mathematica

【解决方案1】：

您使用 numpy（显然），但我自己从未尝试过：http://docs.scipy.org/doc/numpy/reference/generated/numpy.roots.html#numpy.roots。

Numpy 还提供多项式类...numpy.poly1d。

这会以数字方式找到根 -- 如果您想要分析根，我认为 numpy 无法为您做到这一点。

【讨论】：

【解决方案2】：

这是来自 simpy 文档的示例：

>>> from sympy import *
>>> x = symbols('x')
>>> from sympy import roots, solve_poly_system

>>> solve(x**3 + 2*x + 3, x)
           ____          ____
     1   \/ 11 *I  1   \/ 11 *I
[-1, - - --------, - + --------]
     2      2      2      2

>>> p = Symbol('p')
>>> q = Symbol('q')

>>> sorted(solve(x**2 + p*x + q, x))
          __________           __________
         /  2                 /  2
   p   \/  p  - 4*q     p   \/  p  - 4*q
[- - + -------------, - - - -------------]
   2         2          2         2

>>> solve_poly_system([y - x, x - 5], x, y)
[(5, 5)]

>>> solve_poly_system([y**2 - x**3 + 1, y*x], x, y)
                                   ___                 ___
                             1   \/ 3 *I         1   \/ 3 *I
[(0, I), (0, -I), (1, 0), (- - + -------, 0), (- - - -------, 0)]
                             2      2            2      2

(a link to the docs with this example)

【讨论】：

【解决方案3】：

您可能想查看SAGE，它是专为数学处理而设计的完整 python 发行版。除此之外，正如 Marcin 强调的那样，我使用 Sympy 处理一些类似的事情。

【讨论】：

是的，SAGE 非常好（尽管它很可能实际上使用 Numpy 来完成这类任务）。

【解决方案4】：

import decimal as dd
degree = int(input('What is the highest co-efficient of x? '))
coeffs = [0]* (degree + 1)
coeffs1 = {}
dd.getcontext().prec = 10
for ii in range(degree,-1,-1):
    if ii != 0:
        res=dd.Decimal(input('what is the coefficient of x^ %s ? '%ii))
        coeffs[ii] = res
        coeffs1.setdefault('x^ %s ' % ii, res)
    else:
        res=dd.Decimal(input('what is the constant term ? '))
        coeffs[ii] = res
        coeffs1.setdefault('CT', res)
coeffs = coeffs[::-1]
def contextmg(start,stop,step):
    r = start
    while r < stop:
        yield r
        r += step
def ell(a,b,c):
    vals=contextmg(a,b,c)
    context = ['%.10f' % it for it in vals]
    return context
labels = [0]*degree
for ll in range(degree):
    labels[ll] = 'x%s'%(ll+1)
roots = {}
context = ell(-20,20,0.0001)
for x in context:
    for xx in range(degree):
        if xx == 0:
            calculatoR = (coeffs[xx]* dd.Decimal(x)) + coeffs[xx+1]
        else:
            calculatoR = calculatoR * dd.Decimal(x) + coeffs[xx+1]

    func =round(float(calculatoR),2)
    xp = round(float(x),3)
    if func==0 and roots=={} :
        roots[labels[0]] = xp
        labels = labels[1:]
        p = xp
    elif func == 0 and xp >(0.25 + p):
        roots[labels[0]] = xp
        labels = labels[1:]
        p = xp

print(roots)

【讨论】：

这几行代码只是在 python 3 中使用了简单的逻辑，并且只导入了 1 个模块。这段代码可用于求解任意长度的多项式