如何在Python中实现二分法

Question

我希望能够编写一个Python程序，运用二分法来确定以下方程的根：

f(x) = -26 + 85x - 91x2 +44x3 -8x4 + x5

二分法是一种用于估算多项式f(x)的根的数值方法。

是否有可用的伪代码、算法或库可以告诉我答案？

- Scrubatpython

Answer 1

基本技术

以下是展示基本技术的代码：

>>> def samesign(a, b):
        return a * b > 0

>>> def bisect(func, low, high):
    'Find root of continuous function where f(low) and f(high) have opposite signs'

    assert not samesign(func(low), func(high))

    for i in range(54):
        midpoint = (low + high) / 2.0
        if samesign(func(low), func(midpoint)):
            low = midpoint
        else:
            high = midpoint

    return midpoint

>>> def f(x):
        return -26 + 85*x - 91*x**2 +44*x**3 -8*x**4 + x**5

>>> x = bisect(f, 0, 1)
>>> print(x, f(x))
0.557025516287 3.74700270811e-16

容差

为了在达到给定的容差时提前退出，可以在循环末尾添加一个测试：

def bisect(func, low, high, tolerance=None):
    assert not samesign(func(low), func(high))   
    for i in range(54):
        midpoint = (low + high) / 2.0
        if samesign(func(low), func(midpoint)):
            low = midpoint
        else:
            high = midpoint
        if tolerance is not None and abs(high - low) < tolerance:
            break   
    return midpoint

- Raymond Hettinger

---

我该如何添加容差，而不是迭代循环n次？ - ultimatetechie

Answer 2

您可以在早期的 Stack Overflow 问题这里看到解决方案，该方案使用了scipy.optimize.bisect。或者，如果您的目的是学习，那么维基百科上关于二分法的伪代码入口也是一个很好的指南，可以帮助您自行在 Python 中实现，正如早期问题的评论者所建议的那样。

- Simon

Answer 3

我的实现比其他解决方案更通用，同时更简单：（并且是公共领域）

def solve(func, x = 0.0, step = 1e3, prec = 1e-10):
    """Find a root of func(x) using the bisection method.

    The function may be rising or falling, or a boolean expression, as long as
    the end points have differing signs or boolean values.

    Examples:
        solve(lambda x: x**3 > 1000) to calculate the cubic root of 1000.
        solve(math.sin, x=6, step=1) to solve sin(x)=0 with x=[6,7).
    """
    test = lambda x: func(x) > 0  # Convert into bool function
    begin, end = test(x), test(x + step)
    assert begin is not end  # func(x) and func(x+step) must be on opposite sides
    while abs(step) > prec:
        step *= 0.5
        if test(x + step) is not end: x += step
    return x

- Tronic

Answer 4

# Defining Function
def f(x):
    return x**3-5*x-9

# Implementing Bisection Method
def bisection(x0,x1,e):
    step = 1
    print('\n\n*** BISECTION METHOD IMPLEMENTATION ***')
    condition = True
    while condition:
        x2 = (x0 + x1)/2
        print('Iteration-%d, x2 = %0.6f and f(x2) = %0.6f' % (step, x2, f(x2)))

        if f(x0) * f(x2) < 0:
            x1 = x2
        else:
            x0 = x2

        step = step + 1
        condition = abs(f(x2)) > e

    print('\nRequired Root is : %0.8f' % x2)


# Input Section
x0 = input('First Guess: ')
x1 = input('Second Guess: ')
e = input('Tolerable Error: ')

# Converting input to float
x0 = float(x0)
x1 = float(x1)
e = float(e)

#Note: You can combine above two section like this
# x0 = float(input('First Guess: '))
# x1 = float(input('Second Guess: '))
# e = float(input('Tolerable Error: '))


# Checking Correctness of initial guess values and bisecting
if f(x0) * f(x1) > 0.0:
    print('Given guess values do not bracket the root.')
    print('Try Again with different guess values.')
else:
    bisection(x0,x1,e)

代码和输出在这里

此外，codesansar.com/numerical-methods/拥有大量的算法、伪代码和使用不同编程语言进行数值分析的程序。

- saki

Answer 5

  def f(x):
    return -26 + 58 * x - 91 * x**2 + 44 * x**3 - 8 * x**4 + x**5


def relative_error(p, p_old):
    return abs((p - p_old) / p)


def bisection(a, b, no_iterations, accuracy):
    i = 1
    f_a = f(a)
    p_old = 0
    while i <= no_iterations:
        p = (a + b) / 2
        f_p = f(p)
        error = relative_error(p, p_old)
        print('%2d %4f %4f %6f %6f %4.6f' % (i, a, b, p, f_p, error))
        if error < accuracy or f_p == 0:
            break
        p_old = p
        i += 1
        if f_a * f_p > 0:
            a = p
            f_a = f_p
        else:
            b = p


bisection(1, 2, 20, 0.00001)

- Nirmal Sankalana

Answer 6

以容差为准：

# there is only one root
def fn(x):
 return x**3 + 5*x - 9
 # define bisection method
def bisection( eq, segment, app = 0.3 ):
 a, b = segment['a'], segment['b']
 Fa, Fb = eq(a), eq(b)
 if Fa * Fb > 0:
  raise Exception('No change of sign - bisection not possible')   
 while( b - a > app ): 
  x = ( a + b ) / 2.0
  f = eq(x)
  if f * Fa > 0: a = x
  else: b = x  
 return x
 #test it
print bisection(fn,{'a':0,'b':5}, 0.00003) # => 1.32974624634

直播链接：http://repl.it/k6q

- Roman Zhak

Answer 7

可以通过将容差作为停止条件来修改上述二分法：

def samesign(a, b):
        return a*b > 0

def bisect(func, low, high):
    assert not samesign(func(low), func(high))
    n = 20
    e = 0.01 # the epsilon or error we justify for tolerance
    for i in range(n):
        if abs(func(low)-func(high)) > e:
            midpoint = (low + high) / 2
            print(i, midpoint, f(midpoint))
            if samesign(func(low), func(midpoint)):
                low = midpoint
            else:
                high = midpoint
        else:
            return round(midpoint, 2)
    return round(midpoint, 2)

- KokoEfraim

Answer 8

我正在尝试进行一些修改：

1. While循环：算法执行的容差和迭代次数
2. 除了根方法外，保存由算法生成的点向量xn（所有c点）和所有图像f(c)的向量
3. 假设Xs是根的给定近似值，保存绝对误差np.linalg.norm(Xs-xn)
4. 保存近似误差：np.linalg.norm(xn+1 - xn)。

这是我的代码：

import numpy as np
def bisection3(f,x0,x1,tol,max_iter):
    c = (x0+x1)/2.0
    x0 = c
    Xs = 0.3604217029603
    err_Abs = np.linalg.norm(x0-Xs)
    itr = 0
    f_x0 = f(x0)
    f_x1 = f(x1)
    xp = [] # sucesion que converge a la raiz
    errores_Abs = []
    errores_Aprox = []
    fs = [] # esta sucecion debe converger a cero 
    while(((x1-x0)/2.0 > tol) and (itr< max_iter)):
        if f(c) == 0:
            return c
        elif f(x0)*f(c) < 0:
            x1 = c
        else :
            x0 = c
        c = (x0+x1)/2.0    
        itr +=1
        err_Abs = np.linalg.norm(x0-Xs)
        err_Aprox = np.linalg.norm(x1 - x0)
        fs.append(f(c))
        xp.append(c)
        errores_Abs.append(err_Abs)
        errores_Aprox.append(err_Aprox)
    return x0,errores_Abs, errores_Aprox,fs,xp

我有一个执行示例：

f  = lambda x : 3*x + np.sin(x) - np.exp(x)
X0_r1 ,  err_Abs_r1,err_Aprox_r1, fs_r1 , xp_r1 =   bisection3(f,0,0.5,1e-5,100)

- robintux