基于Python实现Aitken迭代法和牛顿迭代法

对于非线性方程，我们可以使用迭代的方式求出近似解。下面介绍两种比较经典的算法：简单迭代法、牛顿法

---

简单迭代法

对于待求解方程，先把方程写成 $f(x)=0$ 的形式，然后改成如下同解形式：
$$x=\phi (x)$$选一个初始值 $x_0$，然后做迭代：
$$x_{k+1}=\phi (x_k),k=1,2,3,\dots$$如果迭代序列 { x k } \{x_k\} {xk​} 收敛于真解，则称简单迭代法使收敛的。

简单迭代法的收敛条件

根据压缩映射原理，如果 $\phi (x)$ 为定义域上的压缩映射，则 $\phi (x)$ 在定义域上有唯一的不动点，且对于任何初值，由迭代格式 $x_{k+1}=\phi (x_k)$ 产生的迭代序列 { x k } \{x_k\} {xk​} 收敛于该不动点。

另外我们还常常讨论简单迭代法的局部收敛性，感兴趣的同学查阅相关书籍，此处不再赘述。

---

简单迭代法的Aitken加速算法

假设 $\phi (x)$ 在 $\alpha$ 处可导，有
$$\begin{gathered} x_{k+1}-\alpha ={\phi }^{\prime }({\xi }_1)(x_k-\alpha ) \\ {x}_{k+2}-\alpha ={\phi }^{\prime }({\xi }_2)(x_{k+1}-\alpha ) \end{gathered}$$假设${\phi }^{\prime }({\xi }_1)\approx {\phi }^{\prime }({\xi }_2)$,有
$$\frac{x_{k+1}-\alpha }{x_{k+2}-\alpha }=\frac{x_k-\alpha }{x_{k+1}-\alpha }$$得到
$$\alpha =x_k-\frac{(x_{k+1}-x_k{)}^2}{x_{k+2}-2x_{k+1}+x_k}$$记
$${\hat{x}}_k=x_k-\frac{(x_{k+1}-x_k{)}^2}{x_{k+2}-2x_{k+1}+x_k}k=1,2,3,\dots$$以上即为Aitken加速算法的迭代格式，序列 { x ^ k } \{\widehat x_k\} {x k​} 要比序列 { x k } \{x_k\} {xk​} 更快地收敛于 $\alpha$ 。

为方便计算可构造如下Aitken加速算法
$$\{\begin{array}{l}{y}_k=\phi (x_k)\\ z_k=\phi (x_k)\\ x_{k+1}=x_k-\frac{{(y_k-x_k)}^2}{z_k-2y_k+x_k},k=0,1,2,\cdots \end{array}$$

---

基于Pyhton实现的Aitken加速算法

算法代码

import numpy as np
import math

def Aitken(x0,epsilon,iternum,phi):#初值，精度要求，最大迭代次数,迭代函数
    xk_1 = x0
    for i in range(iternum):
        y = phi(xk_1)
        z = phi(y)
        if (z - 2*y +xk_1)!= 0:
            xk = xk_1 - (y - xk_1)**2 / (z - 2*y +xk_1)
            print("第",i+1,"次迭代 ","xk=",xk,"  xk-1=",xk_1,"  |xk - xk-1|=",abs(xk-xk_1));
            if abs(xk-xk_1)<epsilon:
                return xk
            else:
                xk_1 = xk
        else:
            return x
    print("方法失败")
    return 0

实验

计算 $x_{k+1}=1.6+0.99cos(x)$ 在初值取 $x_0=\frac{\pi }{2}$ 时的根。

def phi(x): #φ
    return 1.6+0.99*math.cos(x)
    
Aitken(1.57080,1/(10)**7,10,phi)#初值1.57080，精度10^-7,最大迭代次数为10

实验结果

最终结果为 1.585471801521943

---

牛顿迭代法

设函数 $f(x)=0$ 在有根区间 $[a,b]$ 上二阶连续可微，$x_0$ 是 $\alpha$ 的近似值，将 $f(x)$ 在 $x_0$ 处作Taylor展开，有
$$f(x)=f(x_0)+f^{\prime }(x_0)(x-x_0)+f^{\prime \prime }({\xi }_0)\frac{(x-x_0{)}^2}{2}$$用其线性主部近似 $f(x)$
$$f(x)\approx f(x_0)+f^{\prime }(x_0)(x-x_0)$$将非线性方程 $f(x)=0$ 近似化为线性方程
$$f(x_0)+f^{\prime }(x_0)(x-x_0)=0$$若 $f^{\prime }(x_0)$ 不为0，则有
$$x_1=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$$一般地，有
$$x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime }(x_k)}$$以上迭代格式称为牛顿迭代格式。

基于Pyhton实现的牛顿迭代法

算法代码

import numpy as np
import math
from scipy.misc import derivative

def Newton(x0,epsilon,iternum,f):#初值，精度要求，最大迭代次数,迭代函数
    xk_1 = x0
    for i in range(iternum):
        fx = f(xk_1)
        fdx = derivative(f,xk_1,dx=1e-6)
        if fdx!= 0:
            xk = xk_1 - fx / fdx
            print("第",i+1,"次迭代 ","xk=",xk,"  xk-1=",xk_1,"  |xk - xk-1|=",abs(xk-xk_1));
            if abs(xk-xk_1)<epsilon:
                return xk
            else:
                xk_1 = xk
        else:
            break
    print("方法失败")
    return 0

实验
求方程 $x{e}^x-1=0$ 在 $0.5$ 处的根，要求精度为 $\epsilon =10^{-5}$ 。

def func(x): #f
    return x*math.exp( x )-1
    
Newton(0.5,1/(10**5),10,func)

实验结果

最终结果为 0.5671432904097838