Polynomial-Funktion kann nicht von Python gelöst werden sympy

Question

Ich habe Probleme, die durch eine Polynomfunktion mit sympy zu lösen. Das folgende Beispiel zeigt einen Fall, der eine Fehlermeldung gibt, die ich nicht verwalten kann. Wenn das Polynom einfacher wird, funktioniert der Solver ordnungsgemäß. Bitte kopieren und fügen Sie den Code ein, um den Fehler auf Ihrem System zu überprüfen.

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) 

Wissen Sie, wie ich ein Ergebnis erzielen kann? Vielen Dank für Ihre Hilfe.

EDIT:

Dies ist die Fehlermeldung:

---

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) 

Answer 1

Wie ich in den Kommentaren erwähnt ich mit sympy nicht vertraut bin, aber hier ist, wie die Wurzeln der Gleichung zu finden Verwendung der beliebige Genauigkeit mpmath Modul.

Um Präzision Probleme zu vermeiden, die mit Python schwimmt, ist es üblich, schwimmt auf mpmath in String-Form zu übergeben, es sei denn, es bequem ist sie aus ganzen Zahlen zu konstruieren. Ich denke, es ist nicht wirklich ein Problem mit Ihrer Gleichung, da Ihre Koeffizienten eher geringe Genauigkeit haben, aber trotzdem ...

Hier ist Ihre Gleichung direkt in mpmath Syntax übersetzt:

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 ist die beliebige Genauigkeit float-Konstruktor mpc ist die komplexe Zahl Konstruktor. Informationen zum Aufrufen dieser Konstruktoren finden Sie in den mpmath-Dokumenten.

Allerdings müssen wir nicht mit I so herumspielen: wir können nur die Koeffizienten direkt als komplexe Zahlen definieren.

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

Ausgang

(1.35717989161653180236360985534 - 0.202974596285109153971324419197j) 
test: (-3.2311742677852643549664402034e-26 + 6.4623485355705287099328804068e-27j) 

Aber wie können wir die anderen Wurzeln finden? Einfach!

Lassen Sie u eine Wurzel von f (x) sein. Dann lasse man f (x) = g (x) (x - u) und jede Wurzel von g (x) ist auch eine Wurzel von f (x). Wir können dies mehrere Male bequem tun, indem Sie eine for Schleife verwenden, die auf eine Liste jede gefunden Wurzel speichert und erstellt dann eine neue Funktion von der vorherige Funktion, diese neue Funktion in einer anderen Liste zu speichern.

In dieser Version verwende ich den "Muller" Solver, wie es für die Suche nach komplexen Wurzeln empfohlen wird, aber es gibt tatsächlich die gleichen Antworten wie mit dem Standard "Sekant" 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]))  

Ausgang

(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.95804877623117526215e-31j) 
(-1.31713649147437845007238988587e-21 + 1.68350641700147843422461467477e-22j) 

In

lambda x,i=i: funcs[i](x)/(x - roots[i]) 

wir i als Standard-Schlüsselwort Argument angeben, so dass der Wert, i hatte, als die Funktion verwendet wird, wurde festgelegt. Ansonsten wird der aktuelle Wert von i verwendet, wenn die Funktion aufgerufen wird, was nicht das ist, was wir wollen.

---

Diese Technik zum Finden mehrerer Wurzeln kann für beliebige Funktionen verwendet werden. Wenn wir jedoch Polynome lösen wollen, hat mpmath einen besseren Weg, alle Wurzeln gleichzeitig zu finden: die polyroots Funktion. Wir müssen ihm nur eine Liste (oder ein Tupel) der Koeffizienten des Polynoms geben.

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

Ausgang

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) 

Wie Sie sehen können, sind die Ergebnisse sehr nahe denen, durch die vorherige Technik erhalten. Die Verwendung von polyroots ist nicht nur genauer, es hat den großen Vorteil, dass wir keine Startapproximation für die Wurzel angeben müssen, mpmath ist schlau genug, um eine für sich selbst zu erstellen.

Answer 2

Dank der Hilfe von PM2Ring habe ich ein komplettes Codebeispiel erstellt, das die Koeffizienten eines beliebigen gegebenen Sympy-Polynoms vierten Grades extrahiert und löst.

import sympy 
from sympy import I 
import mpmath as mp 
import numpy as np 

omega = sympy.symbols('omega') 

# Define polynomial in sympy 
def function(omega): 
    return (- 0.34     *omega**4  \ 
      + (7.44 + 4.51*I)   *omega**3  \ 
      + (87705.64 - 53.08*I) *omega**2 \ 
      - (144140.03 + 22959.95*I)*omega \ 
      + 42357.18 + 50317.77*I) 

# Show defined polynomial 
print''+ str(function(omega).expand()) +'\n' 
result = sympy.Poly(function(omega).expand(), omega) 

# Obtain coefficients of the polynomial 
coeff = (
    mp.mpc(str(np.real(sympy.lambdify((I),result.coeffs()[0])(1j))) , str(np.imag(sympy.lambdify((I),result.coeffs()[0])(1j)))), 
    mp.mpc(str(np.real(sympy.lambdify((I),result.coeffs()[1])(1j))) , str(np.imag(sympy.lambdify((I),result.coeffs()[1])(1j)))), 
    mp.mpc(str(np.real(sympy.lambdify((I),result.coeffs()[2])(1j))) , str(np.imag(sympy.lambdify((I),result.coeffs()[2])(1j)))), 
    mp.mpc(str(np.real(sympy.lambdify((I),result.coeffs()[3])(1j))) , str(np.imag(sympy.lambdify((I),result.coeffs()[3])(1j)))), 
    mp.mpc(str(np.real(sympy.lambdify((I),result.coeffs()[4])(1j))) , str(np.imag(sympy.lambdify((I),result.coeffs()[4])(1j)))), 
) 

# Calculate roots of the polynomial 
roots = mp.polyroots(coeff) 
for no, frequency in enumerate(roots): 
    print frequency 

Ausgabe

-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 

(1.35717989161653 - 0.202974596285109j) 
(0.285956922243967 + 0.465618100581395j) 
(-497.864871297036 + 6.49836193448774j) 
(518.104087424352 + 6.50370044356891j)