查找根/功能

Question

的零，我嘗試使用二分法指出，找到一個函數的根：

if f(a)*f(b) < 0 then a root exists, 
then you repeat with f(a)*f(c)<0 where c = (a+b)/2  

但林不知道如何修復代碼，以便它工作正常。 這是我的代碼，但它無法正常

from scipy import * 
from numpy import * 


def rootmethod(f, a, b, tol): 


    x = a 
    fa = sign(eval(f)) 

    x = b 
    fb = sign(eval(f)) 

    c = a + b 
    iterations = 0 

    if fa == 0: 
     return a 
    if fb == 0: 
     return b 

    calls = 0   
    fx = 1 

    while fx != 0: 
     iterations = iterations + 1 
     c *= 0.5 
     x = a + c 
     fc = sign(eval(f)) 
     calls = calls + 1 

     if fc*fa >= 0: 
      x = a 
      fx = sign(eval(f)) 
     if fc == 0 or abs(sign(fc)) < eps: 
      fx = sign(eval(f)) 
      return x, iterations, calls 





print rootmethod("(x-1)**3 - 1", 1, 3, 10*e-15) 

新的編輯工作..但仍然不起作用

if fa*fb < 0: 

     while fx != 0: 
      iterations = iterations + 1 
      c = (a + b)/2.0 
      x = c 
      fc = sign(eval(f)) 
      calls = calls + 1 

      if fc*fa >= 0: 
       x = c 
       fx = sign(eval(f)) 
      if fc == 0 or abs(sign(fc)) < tol: 
       fx = sign(eval(f)) 
       return x, iterations, calls 

編輯：改變C =（A + B）* 2到c =（A + b）/ 2在該方法的描述中。

Answer 1

坦白說你的代碼有點亂。這是一些有效的。閱讀循環中的註釋。 （BTW解決您給出的函數是2，而不是3.75）

from scipy import * 
from numpy import * 


def rootmethod(f, a, b, tol): 


    x = a 
    fa = sign(eval(f)) 

    x = b 
    fb = sign(eval(f)) 

    c = a + b 
    iterations = 0 

    if fa == 0: 
    return a 
    if fb == 0: 
    return b 

    calls = 0   
    fx = 1 

    while 1: 
    x = (a + b)/2 
    fx = eval(f) 

    if abs(fx) < tol: 
     return x 

    # Switch to new points. 
    # We have to replace either a or b, whichever one will 
    # provide us with a negative 
    old = b # backup variable 
    b = (a + b)/2.0 

    x = a 
    fa = eval(f) 

    x = b 
    fb = eval(f) 

    # If we replace a when we should have replaced b, replace a instead 
    if fa*fb > 0: 
     b = old 
     a = (a + b)/2.0 




print rootmethod("(x-1)**3 - 1", 1, 3, 0.01) 

Answer 2

我覺得 您 一個問題是：

x = a + c 

由於c = (a + b)*.5，你不需要在這裏a ...添加

更新

你不似乎沒有檢查是否fa * fb < 0開始，我也看不到你在哪裏縮小範圍：你應該重新分配a或b到c，然後重新計算c。

代碼這已經有一段時間，因爲我跟蟒蛇上次播放的，所以把它與鹽^ _^

x = a 
fa = sign(eval(f)) 

x = b 
fb = sign(eval(f)) 

iterations = 0 

if fa == 0: 
    return a 
if fb == 0: 
    return b 

calls = 0   
fx = 1 

while fa != fb: 
    iterations += 1 
    c = (a + b)/2.0 
    x = c 
    fc = eval(f) 
    calls += 1 

    if fc == 0 or abs(fc) < tol: 
     #fx = fc not needed since we return and don't use fx 
     return x, iterations, calls 
    fc = sign(fc) 
    if fc != fa: 
     b = c 
     fb = fc 
    else 
     a = c 
     fa = fc 
#error because no zero is expected to be found 

Answer 3

的糧食，我相信你的循環應該是這樣的（僞代碼，並留出一些檢查）：

before loop: 
a is lower bound 
b is upper bound 
Establish that f(a) * f(b) is < 0 

while True: 
    c = (a+b)/2 
    if f(c) is close enough to 0: 
     return c 
    if f(a) * f(c) > 0: 
     a = c 
    else 
     b = c 

換句話說，如果中點是沒有答案的，然後使它根據其簽署新的終端之一。

Answer 4

請注意，該代碼具有由四捨五入誤差簡單的缺乏

a=0.015707963267948963 
b=0.015707963267948967 
c=(a+b)*.5 

c再次成爲b（一探究竟！）。 如果像1e-16那樣有很小的公差，你可以在無限循環 中結束。

def FindRoot(fun, a, b, tol = 1e-16): 
    a = float(a) 
    b = float(b) 
    assert(sign(fun(a)) != sign(fun(b))) 
    c = (a+b)/2 
    while math.fabs(fun(c)) > tol: 
    if a == c or b == c: 
     break 
    if sign(fun(c)) == sign(fun(b)): 
     b = c 
    else: 
     a = c 
    c = (a+b)/2 
    return c 

現在，一遍又一遍地調用eval並不是很有效率。 這裏是你可以做什麼，而不是

expr = "(x-1.0)**3.0 - 1.0" 
fn = eval("lambda x: " + expr) 
print FindRoot(fn, 1, 3) 

或者你可以把FindRoot內EVAL和lambda定義。

它對您有幫助嗎？

Reson