Python数值计算（6）

从本篇开始及接下来的几篇，将会介绍一下非线性方程的求根中最常见的几种算法及其实现。

1. 二分法

对于一个在区间[a,b]上连续的函数f(x)，如果f(a)*f(b)<0,则f(x)在[a,b]上至少有一个根，这个是根的保证，其次，我们可以每次缩小求根区间，记c=(a+b)/2，如果f(c)与f(a)符号相同，则根会落在[c,b]之间，如果f(c)与f(b)符号相同，则根会落在[a,c]之间，如果f(c)等于0呢？那更好办，c不就是根么……每次将求根区间缩小一半，这就是二分法的精髓，当然证明也很简单（反证法），这里就没必要去废话了。

直接丢出实现的代码：

def bisect(f,a,b,tol):
    c=(a+b)/2
    k=f(c)
    while math.fabs(k)>tol:
        if k*f(a)>0:
            a=c
        else:
            b=c
        c=(a+b)/2
        k=f(c)
    return c

如果我们要求:

的根，显然在[1,2]上存在根，计算根的代码为：

def func(x):
    return x**3-x-1
x=bisect(func,1,2,1e-6)
print(x,func(x)) # 1.3247179985046387 1.759583065918946e-07

我们可以为求根增加绘图，更加详细阐述迭代的过程（但是，个人不推荐这种方式，不建议将计算过程和呈现过程混合在一起）：

def bisect(f,a,b,tol,draw=False):
    if draw:
        x=np.linspace(a,b,100)
        fh=np.vectorize(f)
        y=fh(x)
        plt.plot(x,y,'r')
        plt.grid()
    
    ret=0
    c=(a+b)/2
    k=f(c)
    while True:
        c=(a+b)/2
        k=f(c)
        plt.plot([c,c],[0,k],'b--')
        if math.fabs(k)<tol:
            ret=c
            break
        if k*f(a)>0:
            a=c
        else:
            b=c
    if draw:
        plt.show()
    return ret

调用该函数，图形如下：

x=bisect(func,0,2,1e-6,True)
print(x) # 1.3247179985046387

2. 不动点迭代

不动点迭代就是用新的计算结果，作为下一轮计算的初始值，如果满足迭代收敛的话，则最后就可以得到根（在允许误差范围内），关于不动点求根的原理和迭代收敛的理论这里不打算展开说明，可以自行查阅数值计算的资料。

实现的代码非常简单：

def fixIterator(f,x,tol):
    while True:
        k=f(x)
        if math.fabs(k-x)<tol:
            break
        x=k
    return k

但是如果你的函数不收敛的话，程序就会跑飞的，例如还是上面那个函数的求根，如果写成:

就会发散，但是如果写成：

则迭代是收敛的。测试代码如下：

def func2(x):
    return math.pow(x+1,1/3)

x=fixIterator(func2,1.5,1e-6)
print(x,x-func2(x)) # 1.324718011988197 4.43451095843983e-08