数值分析 非线性方程与方程组的数值解法(二)

牛顿法及其变形方法

牛顿法

牛顿法的构造

  如果方程 $f(x)=0$ 是线性方程，则它的根是容易求解的. 牛顿法是一种线性化的近似方法，其基本思想是将非线性方程 $f(x)=0$转化为线性方程来迭代求解.
  设 $f(x)$在 $x^{\ast }$ 附近二次连续可微， 是附近的一个近似值（假定$f^{\prime }(x_0)\ne 0$ ),将函数在点$x_0$处做一阶 Taylor 展开，有$$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$$其根记为 $x_1=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$，即为 $x^{\ast }$的新近似值. 重复上述过程，将函数$f(x)$在点 $x_1$展开得
$$f(x)\approx f(x_1)+f^{\prime }(x_1)(x-x_1)$$则方程$f(x)=0$ 近似为线性方程$$f(x_1)+f^{\prime }(x_1)(x-x_1)=0$$
得到新的近似解 $x_2=x_1-\frac{f(x_1)}{f^{\prime }(x_1)}$
如此继续下去，得到迭代序列$$x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime }(x_k)},k=0,1,2,\cdots$$
这种方法称为牛顿法，上式称为牛顿迭代公式.

几何解释

牛顿法在几何上是运用曲线$y=f(x)$的切线与$x$轴的交点坐标来近似曲线与$x$轴交点坐标的，因此牛顿法也称为切线法。

计算步骤

函数代码

import random
import sympy

def Newton_itera(a, b, f, R):
    x = sympy.symbols('x')
    X_k = [random.uniform(a, b), "x2"] #构造XK序列
    X_k[-1] = float(X_k[-2] - f.subs(x, X_k[-2]) / (sympy.diff(f, x).subs(x, X_k[-2])))
    while (abs(X_k[-1] - X_k[-2]) >= R):
        x_k = float(X_k[-1] - f.subs(x, X_k[-1]) / (sympy.diff(f, x).subs(x, X_k[-1])))
        X_k.append(x_k)
    return X_k

例题

我们可以通过牛顿法求解分数和根式

1. $\sqrt{C}$
2. $\frac{1}{C}$

当$C=115时，求解\sqrt{115}的近似值$
  解 由于$\sqrt{C}$是方程$x^2-C=0$的正根，因此取$f(x)=x^2-C,f^{\prime }(x)=2x$构造迭代公式$$x_{k+1}=\frac{1}{2}(x_k+\frac{C}{x_k}),k=0,1,2,\cdots$$
计算结果如下

k	0	1	2	3	4
$x_k$	10	10.750000	10.723837	10.723805	10.723805

def Newton_itera(a, b, f, R):
    x = sympy.symbols('x')
    X_k = [10, "x2"]
    X_k[-1] = float(X_k[-2] - f.subs(x, X_k[-2]) / (sympy.diff(f, x).subs(x, X_k[-2])))
    while (abs(X_k[-1] - X_k[-2]) >= R):
        x_k = float(X_k[-1] - f.subs(x, X_k[-1]) / (sympy.diff(f, x).subs(x, X_k[-1])))
        X_k.append(x_k)
    return X_k


def main():
    x = sympy.symbols('x')
    f = x**2-115
    print(Newton_itera(10,11,f,10**-7))



if __name__ == "__main__":
    main()

结果

[10, 10.75, 10.723837209302326, 10.723805294811097, 10.723805294763608]

牛顿法的变形

简化牛顿法

  为了避免牛顿迭代法中的导数计算，取某常数$M\ne 0$代替牛顿迭代公式中的$f^{\prime }(x_k)$,得到简化牛顿迭代公式$$x_{k+1}=f(x_k)-\frac{x_k}{M},k=0,1,2.\cdots$$显然，$M最好取与f^{\prime }(x^{\ast })$较为接近的常数，有时$M取f^{\prime }(x_0)$,此时简化牛顿法为线性收敛。
  简化牛顿法的几何解释：过曲线$y=f(x)$上点$P_k(x_k,f(x_k))$作斜率$M$的平行弦，与$x$轴的交点横坐标$x_{k+1}$作$x^{\ast }$的近似值，因此简化牛顿法也称为平行弦法

牛顿下山法

  由牛顿法局部收敛性可知，初始值$x_0$应选取根在$x^{\ast }$附近，但对于有些问题，很难检验满足条件的初始值$x_0$,此时可利用下山法扩大初始值的选取范围，将牛顿迭代公式修改为$$x_{k+1}=x_k-\lambda \frac{f(x_k)}{f^{\prime }(x_k)},k=0,1,2,\cdots$$式中$\lambda$是一个参数，其选取应使得$$|f(x_{k+1})|<|f(x_k)|$$成立，这种方法称为牛顿下山法，$\lambda$称为下山因子，满足$0<{\epsilon }_r\le \lambda \le 1$,${\epsilon }_r$为下山因子下界，为方便，一般开始时可简单的取$\lambda =1$,然后逐步分半减小，即选取$\lambda =1,\frac{1}{2}.\frac{1}{4},\cdots ,\lambda \ge {\epsilon }_r$,直至$|f(x_{k+1})|<|f(x_k)|$成立。

割线法

  为避免计算f"(x)，由导数定义，有$f^{\prime }(x)\approx \frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}$，代入牛顿迭代公式得$$x_{k+1}=x_k-\frac{f(x_k)}{f(x_k)-f(x_{k-1})}(x_k-x_{k-1})$$使用式(2.17）求解非线性方程的方法称为割线法（或弦截法）.
  割线法的几何解释:过曲线$y=f(x)$上点$P_{k-1}(x_{k-1},f(x_{k-1}))$及$P_k(x_k,f(x_k))$作割线，记割线与x轴交点的横坐标为$x_{k+1}$，用割线与$x$轴的交点近似切线与x轴的交点.
  利用割线法计算$x_{k+1}$,时要用到前两步的结果$x_k$,和$x_{k-1}$,这类算法称为两步迭代法．而牛顿法在计算$x_{k+1}$,时只用到前一步的值$x_k$，这类算法称为单步方法.因此，在使用割线法公式时必须先给出两个初始值$x_0$和$x_1$.

多项式方程求根法

  很多实际应用问题需要求多项式的全部零点，它等价于求多项式方程（ n>1 )
$$P(x)=a_0{x}^n+a_1{x}^{n-1}+\cdots +a_{n-1}x+a_n=0\text{\hspace{1em}(2.20)}$$
的全部根，式中系数$a_0,a_1,\cdots ,a_n$是实数.
  前面介绍的任何非线性方程求根方法均适用于多项式方程求根，但由于多项式方程的特殊性，因此可以针对其特点提供更有效的算法，通常使用牛顿法最好.本节介绍求解多项式方程( n>1）的牛顿法.
  定理：对于任意多项式$c(x)$和$d(x)$，其中$d(x)\ne 0$，存在唯一的多项式$Q(x)$和$r(x)$，满足$c(x)=Q(x)d(x)+r(x)$且$r(x)$的次数小于$d(x)$的次数. $Q(x)$称为商多项式，$r(x)$称为余多项式.
  由定理可知，当$c(x)$为$n$次，$d(x)$为$m$次，且$m<n$，$Q(x)$为$n-m$次，$r(x)$次数小于$m$次.
  由于牛顿迭代公式中涉及函数和函数的导函数运算，因此，对于代数多项式方程，其计算量较大.多项式的牛顿法基本思想是:将牛顿法应用于求解多项式方程,再利用秦九韶方法求出牛顿公式中的函数值$f(x_k)$和$f^{\prime }(x_k)$ .
  将牛顿法应用于多项式方程，有迭代公式
$$x_{k+1}=x_k-\frac{P_n(x_k)}{P_n^{\prime }(x_k)},k==0,1,2,\cdots \text{\hspace{1em}(2.21)}$$
  借助秦九韶方法，不直接使用牛顿迭代公式以及其导数表达式，分别建立$P(x_k)$，$P^{\prime }(x_k)$与多项式系数的关系，得到$P_n(x_k)$，$P_n^{\prime }(x_k)$，从而简化计算.
  用$(x-x_k)$除多项式$P_n(x)$，设商$$Q(x)=b_0{x}^{n-1}+b_1{x}^{n-2}+\cdots +b_{n-2}x+b_{n-1}\text{\hspace{1em}(2.22)}$$
余数为$b_n$,，则$$P_n(x)=Q(x)(x-x_k)+b_n\text{\hspace{1em}(2.23)}$$
将多项式方程的基本形式和$Q(x)$带入上式得到如下递推关系$$\{\begin{array}{ll}{b}_0=a_0& \\ b_1=a_1+b_{i-1}{x}_k,& i=1,2,\cdots ,n\end{array}\text{\hspace{1em}(2.24)}$$

同理，再用(x一x,)除Q(x)，得到商$$H(x)=c_0{x}^{n-2}+c_1{x}^{n-3}+\cdots ＋c_{n-3}x+c_{n-2}\text{\hspace{1em}(2.25)}$$
余数为$c_{n-1}$，则$$Q(x)=H(x)(x-x_k)+c_{n-1}\text{\hspace{1em}(2.26)}$$
将式（2.22）和式(2.25）代入上式，比较系数，得如下递推关系式$$\{\begin{array}{ll}{c}_0=b_0& \\ c_i=b_i+c_{i-1}{x}_k,& i=1,2,\dots ,n一1\end{array}\text{\hspace{1em}(2.27)}$$
由式可得$Q(x_k)=c_{n-1}$,对基本形式求导，可得到$P^{\prime }(x_k)=Q(x_k)=c_{n-1}$因此，通过递推式(2.27）计算$c_{n-1}$，可得到$P_n^{\prime }(x_k)$.

计算步骤

综上所述，牛顿法求多项式方程根的计算步骤归纳如下.

1. 取初始值$x_0$,通常可以取$x_0=0$.
2. 对$k=0,1,2,\dots$计算
   $$b_0=a_0,b_i=a_i+b_{i-1}{x}_k,(i=1,2,\cdots ,n)$$$$c_0=b_0,c_j=b_j+c_{j-1}{x}_k,(j=2,3,\cdots ,n-1)$$$$x_{k+1}=x_k-\frac{b_n}{c_{n-1}}$$
3. 当误差$|b_n|<\epsilon$或$|x_{k+1}一x_k|=|\frac{b_n}{c_{n-1}}|<\epsilon$时，迭代终止;否则，转步骤2继续迭代.

代码实现

def multinomial_root(coef,R,x0):
    '''多项式求根公式的迭代实现'''
    #读入多项式阶数和迭代初值x0
    n= len(coef)
    X_k = [x0]
    
    while True:
        B, C = [coef[0]], [coef[0]]                  #初始化b0 = c0 =a0
        for i in range(1, n):
            B.append(coef[i] + B[i - 1] * X_k[-1])   #计算bi
            C.append(B[i] + C[i - 1] * X_k[-1])      #计算cj
        X_k.append(X_k[-1] - B[n - 1] / C[n - 2])    #计算x{k+1}
        if (abs(X_k[-1]-X_k[-2]))< R:                # 判断 |x{k+1}-xk| 是否达到精度
            return X_k