查找根/功能

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