二分法

方程的求解是科学计算中不可或缺的一个部分，这里主要对单一的方程进行求解，针对超越方程这种非多项式构成的方程，由于其复杂度较高，获得解析解基本不太可能。一般常用的数值方法有二分法、不动点原理迭代法和最为著名的牛顿-拉夫逊法。同时还会指出每一种方法的代价是什么，谈到迭代技术，则收敛速度必须要得到关注；此外，对应的代码会利用Python给出。

二分法是一种十分重要的计算或查找逻辑，也许二分法这一种方法不是最优秀的求解方法，但是二分法的这一种求解或搜寻的逻辑思路是一种应用广泛的思想，不仅仅在数值计算领域，在计算机编程算法中，二分思想是一个很有效的算法。在数值分析中，二分法的产生的数学基础是依赖于函数的中值定理以及函数的连续性。更通俗一点来看，二分法的思路就是将正确的解“括住”，寻找到一个足够小的区间，满足事先设定好的精度即可。

函数的中值定理：

若$f(x)$是闭区间$[a, b]$上的连续函数，如果说该函数$f(x)$在区间边界处满足这个关系式:$f(a) \cdot f(b) < 0$，那么该函数在该区间上必定存在一个根，即使得函数为0的自变量取值必定在区间$(a, b)$内。

思路与过程

这一定理的本质意义可以这样理解。首先，函数的连续性可以让函数曲线的变化是平滑、不间断的，就排除了函数值跳跃的可能性；再者，这一定理中给出了一个先决条件$f(a) \cdot f(b) < 0$，即左端点，右端点处的函数值是异号的，前面由于连续性的要求，对于这一情况的描述之可能存在唯一的一种形式，就是：如果从一个函数值小于0的端点处开始，函数的变化必须平滑、连续，那么函数值要向从负值，变化到正值，一定会穿过$f(x) = 0$这一水平线，那么此时的$x$值就是对应方程的根。

上述只是阐述了“二分法”的计算思路，还没有体现出“二分”的概念；让想得到足够精确的解，还需对这个区间的搜索做一点控制工作，此时“二分”的概念得以体现。在寻找“括住根”的区间过程中，很难一步得到一个很精确的区间，因此，迭代思想必定会在这个过程中得到应用，第一步区间可以大一点，但是必须满足中值定理，只要中值定理满足，则这个区间就是具有意义的，为了减小区间范围，对区间长度进行对半分，此时一个关键问题不可避免，这两个区间到底使用哪一个？这一定是一个选择性问题。当我们对二分法求解的数学思路有了认识后，应当可以感觉到，要想把根括住，必须要让函数在区间端点的值是异号，也是说如何选择两个区间中的哪一个，唯一的判断依据就是“继续保持中值定理的满足”，只要这个区间还满足中值定理，那么就选择它。在进行多次的循环后，由于每一次区间长度都减半，一定会得到一个较小的区间。每一次完整的计算后，解的不确定性都会减小一半，得到的最后一个区间，取这个区间的中点作为估计值即可。

精度与收敛速度

对于二分法的精度的确定，完全可以从整个求解过程中得到，每一次计算都会使得区间长度减半。如果初始的区间是$[a, b]$，则完成第一次完整的计算后，区间长度变为$(b - a) / 2$；第二次计算完成后区间长度变为$(b - a) / 4$，以此类推可以得到进行了n次迭代后的区间长度为$(b - a) / 2^{n}$；注意在最后得到的区间长度是$(b - a) / 2^{n}$，而方程的根是取这个区间的中点作为估计值，由此可以得到二分法的精度为:

$$ eps = \frac{b - a}{2 ^{n + 1}} $$

对于收敛速度或者效率来说，可以自然得到为$(b - a) / 2$。

二分法的实现代码

前述二分法过程中，包含了两个重要信息；第一需要多次的迭代计算，使得估计值越来越精确，这意味着需要使用循环结构来实现；此外，在任何的一次计算中，需要对当前区间的二分后，要取选择下一次的区间时，要利用中值定理的条件去做判断，所以分支结构也一定会被利用。

伪代码对求解的描述

初始区间为$[a, b]$，其已经满足端点处函数值为异号的要求。

最外层：利用循环结构控制:循环控制条件为(b - a) / 2 与预先设定的tolerance的关系
    获得初始区间的中点：c = (b - a) / 2（中点作为估计值）
    进行二分后区间的选择：利用分支语句
    if f(c) == 0:
        这个就是正确解，直接可以结束循环
    if f(a) * f(c) < 0:
        左侧端点和中点c满足中值定理的要求，则选择它
        b = c 直接利用c作为右侧端点即可
    else: 此时意味着左侧的a处函数值与c处的值不满足中值定理，选择另一半
        a = c c作为左侧端点使用即可

Input的确定:

• 初始区间（这个区间必须满足中值定理的描述），以端点值最为输入即可
• tolerance值得给定，用于迭代计算得终止
• 所要求解的方程对应的函数$f(x)$，以def定义的函数形式引入即可

Output：最终的估计值是必须要返回的。其他的可以根据自己的需求而设定

CODE

# 库的调用
import numpy as np  # 导入numpy库，且重命名为np
import matplotlib.pyplot as plt

# 定义了一个二分法求解方程的函数
def bisection_method(f, a, b, tol=0.0001):
    # f:待求的方程对应的函数
    # a:初始区间的左侧端点值
    # b:初始区间的右侧端点值
    # tol: 当区间长度小于这个
    # 值时，结果可以接受，可结
    # 束计算，默认值为0.0001
    
    # 确定端点对应的函数值
    fa = f(a)
    fb = f(b)
    # 利用while循环控制
    while (b - a) / 2 > tol:  # 入口条件满足就要执行
        mid = (b + a) / 2  # 这一步实际上就是进了区间二分
        # 中点处的函数值f(mid)
        fmid = f(mid)
        if fmid == 0:
            return mid  # 此时c就是真正的解
        # 利用if语句进行区间的选择
        # 注意，区间的左、右始终用变量a，b来表示
        if fa * fmid < 0:
            # 选择中点左侧的区间
            b = mid
            fb = fmid
        else:
            a = mid
            fa = fmid
    # 上述循环结束后，区间a，b的值均已经发生了
    # 变化，取此时的中点作为估计值即可。
    return (a + b) / 2

eg:

求方程$\sin(x) - 6x -5=0$的根。

首先确定出方程根大致的区间

# 先定义出函数
def fun(x):
    return np.sin(x) - 6 * x - 5

# 利用函数图像来确定初始区间范围
# 产生x的离散点，范围可自行尝试。
# 此为-10，10之间的，产生100个点
x = np.linspace(-10, 10, 100)
# 带入定义好的函数fun
y = fun(x)

fig, ax = plt.subplots()
ax.plot(x, y,linestyle='--', linewidth=2, color='tab:blue',label='fun')
ax.set_title('$\sin(x) - 6x -5$')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.grid('True')
ax.legend(loc='best')
plt.show()

由此，先可以大致确定出初始区间可以选择为$[-2.5, 0]$。

# 调用bisection函数进行求解
# 前述的代码已经完成了函数的定义以及方程对应函数的确定
# 此时，要对方程进行计算即可
sol = bisection_method(fun, -2.5, 0, 0.0001)
print(sol)

结果为$-0.9708404541015625$。

OOP方式的编程实现

为什么要采用面向对象的编程？上述采用面向过程的编程范式，通过函数定义的实现计算过程，虽然也很好地解决了问题，但是更进一步地想，我们真正面对的问题是“求方程的根”，至于用什么数值方法可能不是我们的目的，并且方法的种类又不唯一。如果按照之前定义函数的方式来实现，求解对象与数值方法函数之间时存在割裂的。而采用OOP时，可以通过类的定义，确定出方程这个对象，然后将方程的具体形式（也就是这个方程对应的函数）作为对象的属性，初始区间、精度都作为对象的属性，而方程的解（或者说是求解）定义为方程对象的方法，而方法的数量又不是唯一的，因此，可以定义多个方法，如二分法、不动点迭代法，以及牛顿-拉夫逊法，这样的封装形式更加方便。

import numpy as np

# 定义一个方程对象
class myequation(object):
    
    def __init__(self, fun, a, b, tol=0.0001):
        self.fun = fun
        self.a = a
        self.b = b
        self.tol = tol
        
    def bisection(self):
        a_ = self.a
        b_ = self.b
        fa = self.fun(a_)
        fb = self.fun(b_)
    
        while (b_ - a_) / 2 > 0:
            mid_ = (b_ + a_) / 2
            fmid = self.fun(mid_)
            if self.fun(mid_) == 0:
                return mid_
            if fa * fmid < 0:
                b_ = mid_
                fb = fmid
            else:
                a_ = mid_
                fa = fmid
        return (a_ + b_) / 2
            
        
# 应用
def fun(x):
    return np.sin(x) - 6 * x - 5

eq = myequation(fun, -2.5, 0)
sol1 = eq.bisection()
print(sol1)

结果为$-0.9708989235042559$。

对于二分法的理解核心是二分这种逼近、缩小范围的思路，对于方程的求解，二分的实现方式被限定在满足函数中值定理的条件下，通过多次的迭代计算获取一个能够接受的小区间，真实解是存在于这个区间上。

如果你觉得这些对你还有些收获！当然可以打赏我。

请点击右下方Tip进行打赏！

本站文章如非特别说明，均为原创。未经本人同意，请勿转载！转载请务必注明出处！