牛顿法及Python实现

目录

1 原理

2 牛顿法求解步骤

3 牛顿法的几何解释

4 案例&Python实现

---

1 原理

牛顿法是基于泰勒公式来实现的。泰勒公式的意义：如果函数满足一定的条件，泰勒公式可以用函数在某一点的各阶导数值做系数构建一个多项式来近似表达这个函数。

设f(x)" role="presentation">f(x)$f(x)$在某邻域内n+1阶可导，则f(x)" role="presentation">f(x)$f(x)$的泰勒展开式为：

f(x)=f(x0)+f′(x0)(x−x0)+f″(x−x0)22!+......+f(n)(x−x0)nn!+Rn(x)" role="presentation">f(x)=f(x0)+f′(x0)(x−x0)+f′′(x−x0)22!+......+f(n)(x−x0)nn!+Rn(x)$f(x)=f(x_0)+f^{\prime }(x_0)(x-x_0)+\frac{f^{{}^{″}}(x-x_0{)}^2}{2!}+......+\frac{f^{(n)}(x-x_0{)}^n}{n!}+R_n(x)$

取其线性部分(即泰勒展开式的前两项)，并令其等于0，即：

f(x0)+f′(x0)(x−x0)=0" role="presentation">f(x0)+f′(x0)(x−x0)=0$f(x_0)+f^{\prime }(x_0)(x-x_0)=0$

以此可以得到非线性方程f(x)=0" role="presentation">f(x)=0$f(x)=0$的近似解。若f′(x0)≠0" role="presentation">f′(x0)≠0$f^{\prime }(x_0)\ne 0$，则可得到方程的一个解：

这样就可以得到牛顿法的迭代公式：

​​​​​​​

已经证明，如果是连续的，并且待求的零点是孤立的，那么在零点周围存在一个区域，只要初始值位于这个邻近区域内，那么牛顿法必定收敛。 并且，如果不为0, 那么牛顿法将具有平方收敛的性能. 粗略的说，这意味着每迭代一次，牛顿法结果的有效数字将增加一倍。

2 牛顿法求解步骤

牛顿法求解步骤如下：

① 求解迭代函数

② 构造迭代公式；

③ 确定迭代初值；

③ 根据精度要求确定迭代停止条件。

3 牛顿法的几何解释

 牛顿法的迭代过程如下：

• 迭代起点P(x0,f(x0))" role="presentation">(x0,f(x0))$(x_0,f(x_0))$；
• 过点(x0,f(x0))" role="presentation">(x0,f(x0))$(x_0,f(x_0))$作f(x)" role="presentation">f(x)$f(x)$的切线，切线与x轴的交点的横坐标为x1" role="presentation">x1$x_1$；
• 更新迭代起点为(x1,f(x1))" role="presentation">(x1,f(x1))$(x_1,f(x_1))$；
• 过点(x1,f(x1))" role="presentation">(x1,f(x1))$(x_1,f(x_1))$作f(x)" role="presentation">f(x)$f(x)$的切线，切线与x轴的交点的横轴标为；
• ......
• 如此循环，切线与x轴交点的横坐标xn" role="presentation">xn$x_n$会无限接近于真实解x∗" role="presentation">x∗$x\ast$。

4 案例&Python实现

用牛顿法求解方程的根，精度要求ε=10−9" role="presentation">ε=10−9$\epsilon ={10}^{-9}$。

输出近似解以及迭代次数。

程序流程如下：

 Python代码如下：

#----牛顿法求根-----#
import numpy as np
def f(x):
    y=x**3-x-1#求根方程的表达式
    return y
def g(x):
    y=3*x**2-1#求根方程的导函数
    return y
def main():
    x0=1.5 #取初值
    e=10**(-9) #误差要求
    L=0 #初始化迭代次数
    while abs((f(x0)-0))>e: #采用残差来判断
        x1=x0-f(x0)/g(x0) #迭代公式,x(n+1)=x(n)-f(x(n))/f'(x(n))
        x0=x1
        L=L+1 #统计迭代次数
    print(f"x1={x1}") #输出数值解
    print(f(x0)-0)  # 验证解的正确性
    print(f"L={L}") #输出迭代次数
if __name__ == '__main__':
   main()

运行结果如下：

x1=1.3247179572447898
1.865174681370263e-13
L=4