2015-08-31 37 views
-1

我是新來的python,我很難通過使用平分法找到多項式的根。到目前爲止,我有2種方法。一個用於評估值x使用平分多項式的根

def eval(x, poly): 
""" 
Evaluate the polynomial at the value x. 
poly is a list of coefficients from lowest to highest. 

:param x:  Argument at which to evaluate 
:param poly: The polynomial coefficients, lowest order to highest 
:return:  The result of evaluating the polynomial at x 
""" 

result = poly[0] 
for i in range(1, len(poly)): 
    result = result + poly[i] * x**i 

return result 

下一個方法多項式應該用二分法查找多項式的根給出

def bisection(a, b, poly, tolerance): 
poly(a) <= 0 
poly(b) >= 0 

try: 
    if 









""" 
Assume that poly(a) <= 0 and poly(b) >= 0. 

:param a: poly(a) <= 0 Raises an exception if not true 
:param b: poly(b) >= 0 Raises an exception if not true 
:param poly: polynomial coefficients, low order first 
:param tolerance: greater than 0 
:return: a value between a and b that is within tolerance of a root of the polynomial 
""" 

我怎麼會發現使用二分法的根源?我已經提供了一個測試腳本來測試這些。

編輯:我跟着僞代碼,並結束了與此:

def bisection(a, b, poly, tolerance): 
#poly(a) <= 0 
#poly(b) >= 0 
difference = abs(a-b) 
xmid = (a-b)/2 
n = 1 
nmax = 60 




while n <= nmax: 
mid = (a-b)/2 
if poly(mid) == 0 or (b - a)/2 < tolerance: 
     print(mid) 

n = n + 1 
if sign(poly(mid)) == sign(poly(a)): 
    a = mid 
else: 
    b = mid 


return xmid 

這是正確的?我沒有能夠測試它,因爲返回xmid語句的縮進錯誤。

+0

只需按照僞代碼在這裏:https://en.wikipedia.org/wiki/Bisection_method? –

+0

你爲什麼要重新發明輪子? http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.bisect.html – Moritz

回答

0

你的代碼似乎很好,除了與xmidmid混亂。 mid = (a + b)/2而不是mid = (a - b)/2,您不需要difference變量。

清理了一下:

def sign(x): 
    return -1 if x < 0 else (1 if x > 0 else 0) 

def bisection(a, b, poly, tolerance): 
    mid = a # will be overwritten 
    for i in range(60): 
    mid = (a+b)/2 
    if poly(mid) == 0 or (b - a)/2 < tolerance: 
     return mid 
    if sign(poly(mid)) == sign(poly(a)): 
     a = mid 
    else: 
     b = mid 
    return mid 

print(bisection(-10**10, 10**10, lambda x: x**5 - x**4 - x**3 - x**2 - x + 9271, 0.00001))