Python数值计算（7）

继续求根的数值算法。这次介绍牛顿法和割线法。

3. 牛顿迭代法

该算法是数值分析中最著名也是最强健的数值方法之一，其迭代公式为：

例如，还是求方程：

在[1,2]之间的根，其算法及测试代码为：

def f(x):
    return x**3-x-1

def df(x):
    return 3*x**2-1

def newton(f,df,p0,tol):
    cnt=0
    ret=0
    while True:
        p=p0-f(p0)/df(p0)
        if math.fabs(p-p0) < tol:
            ret=p
            break
        cnt+=1
        p0=p
    return ret,cnt
    
x,cnt=newton(f,df,1,1e-6)
print(cnt,x) # 4 1.3247179572447898

可以看到，只需要4次就可以计算出误差在1e-6的水平。

4. 割线法

牛顿法速度快效率高，但是有个很大的问题是必须要事先知道函数的导数，有时候计算导数是一件费时费力的事情，例如上例中是预选定义好了导函数。为了克服这个问题，可以逆向思维，导数的定义之一通过计算微商，然后取极限的方式，可以将处的导数用微商代替，这就是割线法的原理，实现代码如下：

def secant(f,a,b,tol):
    ret=0
    cnt=0
    p0,p1=a,b
    q0,q1=f(p0),f(p1)
    while True:        
        p=p1-q1*(p1-p0)/(q1-q0)
        if math.fabs(p-p1)<tol:
            ret=p
            break
        p0=p1
        q0=q1
        p1=p
        q1=f(p)
        cnt+=1
    return ret,cnt

x,cnt=secant(f,1,2,1e-6)
print(cnt,x) # 6 1.3247179572446703

测试显示通过6次迭代也可以实现误差在1e-6的水平。