Risolvere equazioni polinomiali in Python

Question

Fino ad ora ho sempre Mathematica per risolvere equazioni analitiche. Ora però ho bisogno di risolvere alcune centinaia di equazioni di questo tipo (polinomio caratteristico)

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

in una volta, che i rendimenti tempi di calcolo terribilmente lunghi in Mathematica.

C'è un comando pronto all'uso in numpy o qualsiasi altro pacchetto per risolvere un'equazione di questo tipo? (fino ad ora ho usato Python solo per le simulazioni, quindi non so molto sugli strumenti analitici e non ho trovato nulla di utile nelle esercitazioni numpy).

Answer 1

Si usa numpy (apparentemente), ma non l'ho mai provato personalmente: http://docs.scipy.org/doc/numpy/reference/generated/numpy.roots.html#numpy.roots.

Numpy fornisce anche una classe polinomiale ... numpy.poly1d.

Questo trova le radici numericamente - se vuoi le radici analitiche, non penso che Numpy possa farlo per te.

Answer 2

Si consiglia di guardare SAGE che è una distribuzione python completa progettata per l'elaborazione matematica. Oltre a ciò, ho usato Sympy per questioni simili, come evidenziato da Marcin.

Answer 3

Ecco un esempio dalla documentazione 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)