寻找函数的根/零点

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

我觉得你遇到的问题之一是：

    x = a + c

因为 c = (a + b)*.5，所以这里不需要再加 a 了...

更新

你似乎没有先检查 fa * fb < 0，另外我也没看到你是怎么缩小范围的：你应该在循环中把 a 或 b 重新赋值为 c，然后再重新计算 c。

代码 我已经有一段时间没玩 Python 了，所以请你多包涵 ^_^

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 2

我觉得你的循环应该像这样（用伪代码表示，省略了一些检查）：

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 3

老实说，你的代码有点乱。这里有一些可以正常工作的代码。请仔细阅读循环中的注释。
(顺便说一下，你给出的函数的正确答案是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)