Rozwiązywanie równań wielomianowych w Pythonie

Question

Do tej pory zawsze miałem Mathematica do rozwiązywania równań analitycznych. Teraz jednak muszę rozwiązać kilkaset tego typu równań (wielomian charakterystyczny)

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

naraz co daje strasznie długie czasy obliczeń w Mathematica.

Czy jest jak gotowa do użycia komenda w numpy lub jakiejkolwiek innej paczce do rozwiązania równania tego typu? (do tej pory używałem Pythona tylko do symulacji, więc nie wiem zbyt wiele o narzędziach analitycznych i nie mogłem znaleźć niczego przydatnego w numpy tutoriali).

Answer 1

Używasz numpy (podobno), ale ja sam nigdy tego nie próbowałem: http://docs.scipy.org/doc/numpy/reference/generated/numpy.roots.html#numpy.roots.

Numpy udostępnia także klasę wielomianową ... numpy.poly1d.

To odnajdzie korzenie liczbowo - jeśli chcesz mieć korzenie analityczne, nie sądzę, żeby numpy mógł to dla ciebie zrobić.

Answer 2

Możesz zajrzeć na SAGE, która jest kompletną dystrybucją pythona przeznaczoną do przetwarzania matematycznego. Poza tym użyłem Sympy do nieco podobnych spraw, jak podkreślił Marcin.

Answer 3

Oto przykład z Simpy dokumentów:

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