python牛顿法求函数极值 python牛顿迭代法求方程解

目录

• 简单迭代法

• 简单迭代法的Aitken加速算法
• 基于Pyhton实现的Aitken加速算法

• 牛顿迭代法

• 基于Pyhton实现的牛顿迭代法

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

---

简单迭代法

对于待求解方程，先把方程写成 $f(x)=0$ 的形式，然后改成如下同解形式：
$x = \varphi(x)$选一个初始值 $x_0$，然后做迭代：
$x_{k+1}=\varphi(x_{k}), \;\;\;\;k=1,2,3,\dots$如果迭代序列 $\{x_k\}$

简单迭代法的收敛条件

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

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

---

简单迭代法的Aitken加速算法

假设 $\varphi(x)$ 在 $\alpha$ 处可导，有
$x_{k+1} - \alpha=\varphi$假设$\varphi$,有
$\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}$记
$\widehat 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加速算法的迭代格式，序列 $\{\widehat x_k\}$ 要比序列 $\{x_k\}$ 更快地收敛于 $\alpha$

为方便计算可构造如下Aitken加速算法
$\left\{\begin{array}{l}y_k=\varphi(x_k)\\z_k=\varphi(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}\right.$

---

基于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\pi2$

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$用其线性主部近似 $f(x)$
$f(x) \approx f(x_0) + f$将非线性方程 $f(x)=0$ 近似化为线性方程
$f(x_0) + f$若 $f$ 不为0，则有
$x_1 = x_0-\frac{f(x_0)}{f$一般地，有
$x_{k+1} = x_k-\frac{f(x_k)}{f$以上迭代格式称为牛顿迭代格式。

基于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

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

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

实验结果

最终结果为 0.5671432904097838