Resolviendo ecuaciones polinomiales en Python

Question

Hasta ahora siempre tengo Mathematica para resolver ecuaciones analíticas. Ahora, sin embargo tengo que resolver unos pocos cientos de ecuaciones de este tipo (polinomios característicos)

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

a la vez, que los rendimientos de los tiempos de cálculo muy largos en Mathematica.

¿Hay algún comando listo para usar en numpy o en cualquier otro paquete para resolver una ecuación de este tipo? (Hasta ahora he usado Python solo para simulaciones, así que no sé mucho sobre herramientas analíticas y no pude encontrar nada útil en los tutoriales numpy).

Answer 1

Usas numpy (aparentemente), pero nunca lo intenté yo mismo: http://docs.scipy.org/doc/numpy/reference/generated/numpy.roots.html#numpy.roots.

Numpy también proporciona una clase de polinomio ... numpy.poly1d.

Esto encuentra las raíces numéricamente, si quieres las raíces analíticas, no creo que numpy pueda hacer eso por ti.

Answer 2

Es posible que desee mirar SAGE que es una distribución completa de python diseñada para el procesamiento matemático. Más allá de eso, he usado Sympy para asuntos algo similares, como resaltó Marcin.

Answer 3

Aquí se muestra un ejemplo de docs 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)

Answer 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)