布谷鸟牛顿迭代算法求解非线性方程组

Cuckoo Newton Iterative Algorithm for Solving Nonlinear Equations

摘要: 近年来，非线性方程组问题越来越多的出现在科学与工程计算领域中。例如机器学习、人工智能、金融计算、石油地质探测、卫星轨道预测等各个领域都涉及到了非线性方程组求解的问题，如何有效快速的求解各类非线性方程组问题受到人们的普遍关注。其中牛顿迭代是求解非线性方程(组)的一个重要方法。本文首先利用布谷鸟搜索算法得到初值，然后利用牛顿迭代法求解非线性方程组。理论分析了这些格式的收敛阶。最后给出一些数值算例对三种迭代格式进行了分析比较，验证了理论并分析了结论。

Abstract: In recent years, the problem of nonlinear equations appears more and more in the field of science and engineering calculation. Such as, machine learning, artificial intelligence, financial computing, petroleum geological exploration, satellite orbit prediction and other fields are involved in nonlinear equations. How to effectively and quickly solve various nonlinear equations has been widely concerned by many researchers. Newton’s method is an important method for solving nonlinear equations. The method in this paper is to use the variation of Newton’s method to improve the order of convergence, and combine with cuckoo algorithm to find a better initial point, so as to achieve the purpose of solving the nonlinear equations. Finally, some numerical examples are used to analyze and compare the different iterative schemes, and support the method.

1. 引言

非线性方程组问题 [1]，在机器学习、人工智能、金融计算、石油地质勘探、电力系统计算、卫星轨道预测等各个领域都有着非常广泛的应用。非线性方程组问题是非线性问题中经常出现的问题之一，其具体表现形式如下所示

$\{\begin{array}{l}{\phi }_1\left(x\right)={\phi }_1\left(x_1,x_2,\cdots ,x_m\right)=0\hfill \\ {\phi }_2\left(x\right)={\phi }_2\left(x_1,x_2,\cdots ,x_m\right)=0\hfill \\ ⋮\hfill \\ {\phi }_m\left(x\right)={\phi }_m\left(x_1,x_2,\cdots ,x_m\right)=0\hfill \end{array}$

用向量形式表示为

$\phi \left(x\right)=0,$

其中 $x={\left(x_1,x_2,\cdots ,x_m\right)}^{\text{T}}$， $\phi ={\left({\phi }_1,{\phi }_2,\cdots ,{\phi }_m\right)}^{\text{T}}$， ${\phi }_i:R^m\to R$， $i=1,2,\cdots ,m$。

非线性方程组在科学与工程计算中处于非常重要的地位。如何有效地利用计算机求解各类非线性问题一直以来都是困扰人们的一个难题。对于大多数非线性方程组而言，我们只能得到解的存在性，而解的形式推导过程相当困难甚至无法得到。因此，非线性方程组的数值求解具有很重要的实际意义。

牛顿迭代法是求解非线性方程(组)的经典方法，其最突出的优点是当初始点充分靠近问题的解时能得到很好的收敛效果。但该方法需要好的初始点。近年来，有无数的科学工作者在研究求解非线性方程及方程组，为该领域提供了很多的算法，杨新社等 [2] 提出一种布谷鸟搜索算法，也叫杜鹃搜索。赵鹏军 [3] 将求解非线性方程组的问题转化为优化问题进行求解。 [4] 提出了3种两步三阶收敛精度的牛顿迭代的变形形式。S. Weerakoon和T. G. I. Fernando [5] 同样提出了一种两步三阶收敛精度的牛顿迭代格式。

本文把杨新社等 [2] 提出一种布谷鸟搜索算法与 [4] 提出的两步三阶收敛的变形牛顿迭代法相结合，求解非线性方程组的数值解。

2. 相关算法描述

2.1. 布谷鸟算法 [2] 描述

基于布谷鸟的行为习惯，模拟布谷鸟寻找宿主鸟巢下蛋的行为，对非线性方程组的解在合理范围[c, d]内进行搜索，做出以下三个理想假设：

1) 每一个布谷鸟一次只下一个蛋，并将其下到随机的一个宿主鸟巢内；

2) 好的宿主鸟巢孵出的高质量的蛋会持续到下一代；

3) 宿主鸟巢的数量是固定的，并且宿主有一定概率 $p_a$ 发现布谷鸟的蛋，此时宿主会选择将蛋舍弃或舍弃鸟巢并在一个新的地方建立一个新鸟巢。

在以上3个假设条件下，当第 $i$ 代布谷鸟 $x_i^{(t+1)}$ 寻找新宿主鸟巢时由模拟莱维飞行产生新的宿主鸟巢的位置：

$x_i^{(t+1)}=x_i^{(t)}+\alpha \oplus L\left(\lambda \right)$ (2.1)

此处 $\alpha >0$，用来对新的宿主鸟巢的位置范围进行控制， $\oplus$ 表示点对点乘法，在大多数情况下， $\alpha$ 取1，其中莱维飞行 $L\left(\lambda \right)$ 表示为：

$L\left(\lambda \right)\sim u=t^{-\lambda },\left(1<\lambda \le 3\right)$ (2.2)

新的位置确定后，判断宿主是否发现布谷鸟的鸟蛋，将随机数 $r\in \left[0,1\right]$ 与 $p_a$ 进行比较，若 $r<p_a$ 则在范围 $\left[c,d\right]$ 内随机产生新的鸟巢位置，反之不变。经过多代延续，找到最具适应性的鸟巢位置，作为结果。

2.2. 布谷鸟算法求解非线性方程组 [3]

设有m个方程形成的非线性方程组有唯一解：

$\{\begin{array}{l}{\phi }_1\left(x\right)={\phi }_1\left(x_1,x_2,\cdots ,x_m\right)=0\hfill \\ {\phi }_2\left(x\right)={\phi }_2\left(x_1,x_2,\cdots ,x_m\right)=0\hfill \\ ⋮\hfill \\ {\phi }_m\left(x\right)={\phi }_n\left(x_1,x_2,\cdots ,x_m\right)=0\hfill \end{array}$ (2.3)

由布谷鸟算法找到的近似解 $x=\left(x_1,x_2,\cdots ,x_m\right)$ 带入

$Sum={\displaystyle \underset{i=1}{\overset{m}{\sum }}\left|{\phi }_i\left(x\right)\right|}$ (2.4)

作为鸟巢优劣的判断标准，即解优劣的标准，Sum越小，鸟巢的位置越合适，近似解 $x=\left(x_1,x_2,\cdots ,x_m\right)$ 就越合适。

2.3. 两步三阶精度牛顿迭代法 [4]

经典牛顿迭代法求解非线性方程组(2.3)的格式如下：给定初始值 $x_0$，

$x_{n+1}=x_n-F^{\prime }{\left(x_n\right)}^{-1}F\left(x_n\right),n=0,1,2,\cdots$ (2.5)

下面给出三种经典牛顿迭代法(2.5)的变形形式

两步算术平均牛顿迭代法(arithmetic mean Newton’ (AN) method)

给定初始值 $x_0$

$\{\begin{array}{l}{z}_{n+1}=x_n-F^{\prime }{\left(x_n\right)}^{-1}F\left(x_n\right),\hfill \\ x_{n+1}=x_n-{\left(\frac{F^{\prime }\left(x_n\right)+F^{\prime }\left(z_{n+1}\right)}{2}\right)}^{-1}F\left(x_n\right),n=0,1,2\cdots \hfill \end{array}$ (2.6)

该格式可以看成取 $F^{\prime }\left(x_n\right)$ 和 $F^{\prime }\left(z_{n+1}\right)$ 的算术平均 $\frac{F^{\prime }\left(x_n\right)+F^{\prime }\left(z_{n+1}\right)}{2}$，代替经典牛顿迭代法(2.5)中的 $F^{\prime }\left(x_n\right)$，因此称该格式为算术平均迭代格式(AN)。

两步调和平均牛顿迭代法(harmonic mean Newton’s method (HN))

给定初始值 $x_0$

$\{\begin{array}{l}{z}_{n+1}=x_n-F^{\prime }{\left(x_n\right)}^{-1}F\left(x_n\right),\hfill \\ x_{n+1}=x_n-{\left(2F^{\prime }\left(x_n\right)F^{\prime }\left(z_{n+1}\right)\right)}^{-1}\left(F^{\prime }\left(x_n\right)+F^{\prime }\left(z_{n+1}\right)\right)F\left(x_n\right),n=0,1,2\cdots \hfill \end{array}$ (2.7)

该格式可以看成取 $F^{\prime }\left(x_n\right)$ 和 $F^{\prime }\left(z_{n+1}\right)$ 的调和平均 $\frac{2F^{\prime }\left(x_n\right)F^{\prime }\left(z_{n+1}\right)}{F^{\prime }\left(x_n\right)+F^{\prime }\left(z_{n+1}\right)}$，代替经典牛顿迭代法(2.5)中的 $F^{\prime }\left(x_n\right)$，因此称该格式为调和平均牛顿迭代(HN)。

两步中点牛顿迭代法(midpoint Newton’s method) (MN)

给定初始值 $x_0$

$\{\begin{array}{l}{z}_{n+1}=x_n-F^{\prime }{\left(x_n\right)}^{-1}F\left(x_n\right),\hfill \\ x_{n+1}=x_n-{\left(F^{\prime }\left(\frac{1}{2}\left(x_n+z_{n+1}\right)\right)\right)}^{-1}F\left(x_n\right),n=0,1,2\cdots \hfill \end{array}$ (2.8)

该格式可以看做用 $x_n$ 和 $z_{n+1}$ 的平均 $\frac{x_n+z_{n+1}}{2}$ 的导数值 $F^{\prime }\left(\frac{x_n+z_{n+1}}{2}\right)$ 代替经典牛顿迭代法(2.5)中的 $F^{\prime }\left(x_n\right)$，因此称该格式为两步中点牛顿迭代法(MN)。

下面将布谷鸟搜索算法(2.1)~(2.2)分别与两步牛顿迭代法(2.6)，(2.7)，(2.8)结合起来。

2.4. 布谷鸟搜索两步牛顿迭代法

2.4.1. 首先估计解的大致范围

通过调整解可能存在的区间范围[c, d]，改变布谷鸟算法的搜索范围，使布谷鸟算法可以更快搜索出一个误差较小的解，记录下此时的区间范围[c, d]，以便多次实验。

2.4.2 应用布谷鸟搜索算法求出一组近似解

使用调整好的区间范围[c, d]，以及步长参数 $\alpha$ 和鸟巢数量，然后对非线性方程组的解进行搜索，根据方程组的复杂程度和实际需求对迭代次数或精度阈值进行调整，最终求出一个精度较高的近似解 $X_0=\left(x_1,x_2,x_3,\cdots ,x_m\right)$。

2.4.3. 应用两步牛顿迭代法进行计算

将上面利用布谷鸟搜索算法得到的近似解 $X_0=\left(x_1,x_2,x_3,\cdots ,x_m\right)$ 作为两步牛顿迭代法(2.6)，(2.7)，(2.8)的初值，然后分别用(2.6)，(2.7)，(2.8)进行计算。例如代入两步算术平均牛顿迭代法(2.6)

$\{\begin{array}{l}{Z}_1=X_0-F{\left(X_0\right)}^{-1}F\left(X_0\right)，\hfill \\ X_1=X_0-2{\left(F^{\prime }\left(X_0\right)+F^{\prime }\left(Z_1\right)\right)}^{-1}F\left(X_0\right)\hfill \end{array}$ (2.9)

得到 $X_1=\left(x_1,x_2,x_3,\cdots ,x_m\right)$，再将 $X_1$ 看做 $X_0$ 代入(2.9)式中，求出 $X_2=\left(x_1,x_2,x_3,\cdots ,x_m\right)$，如此反复循环计算，最后得到精度较高的数值解 $X_n=\left(x_1,x_2,x_3,\cdots ,x_m\right)$。其它两个牛顿迭代方法(2.7)，(2.8)与之类似。

下面给出两步牛顿法(2.6)，(2.7)，(2.8)的收敛性分析。

3. 收敛性分析

定义3.1. [4] 如果一个序列 $\left\{x_n|n>0\right\}$ 收敛于 $x_{\ast }$ 满足：

$\underset{x\to x_{\ast }}{\lim }\frac{x_{n+1}-x_{\ast }}{{\left(x_n-x_{\ast }\right)}^p}=C$ (3.1)

其中 $C\ne 0$，且 $p\ge 1$，则该序列的为 $p$ 阶收敛。设 $e_n=x_n-x_{\ast }$，则

$e_{n+1}=C{e}_n^p+O\left(e_n^{p+1}\right)$ (3.2)

称为该迭代法的误差计算公式， $p$ 为迭代法的收敛阶。

定义3.2. [6] 设 $x_{\ast }$ 为方程组 $F\left(x\right)$ 的一组解， $x_{n-2}$， $x_{n-1}$， $x_n$ 和 $x_{n+1}$ 为收敛于 $x_{\ast }$ 的四个相邻迭代近似值，则该方法的收敛阶 $p$ 可通过下式近似计算得到：

$p=\frac{\ln \left(\Vert x_{n+1}-x_n\Vert /\Vert x_n-x_{n-1}\Vert \right)}{\ln \left(\Vert x_n-x_{n-1}\Vert /\Vert x_{n-1}-x_{n-2}\Vert \right)}$ (3.3)

定理3.1. 设函数 $F\left(x\right):\mathbb{D}\subseteq {ℝ}^n\to {ℝ}^n$，非线性方程组 $F\left(x\right)$ 在邻域 $\mathbb{D}$ 上连续可微，且存在一组解 $x_{\ast }\in \mathbb{D}$ 使得 $F\left(x_{\ast }\right)=0$。 $F\left(x\right)$ 的雅可比矩阵 $F^{\prime }\left(x\right)$ 在 $D$ 上是连续且非奇异，且 $F\left(x\right)$ 在 $x_{\ast }$ 的某局部邻域内具有连续的三阶导数，则两步牛顿迭代法(2.6)是三阶收敛的，且满足误差方程：

$e_{k+1}=e_k-\left(C_2^2+\frac{1}{2}{C}_3\right)e_k^3+O\left(e_k^4\right)$ (3.4)

其中 $e_k=x_k-x_{\ast }$， $F\left(x_{\ast }\right)=0$， $C_k=k!F^{\prime }{\left(x_{\ast }\right)}^{-1}{F}^{(k)}\left(x_{\ast }\right)$， $k=2,\cdots ,n\left(C_1=F^{\prime }\left(x_{\ast }\right)\right)$。

证 分别对 $F\left(x_k\right)$， $F^{\prime }\left(x_k\right)$ 在 $x_{\ast }$ 处做Taylor级数展开：

$\begin{array}{c}F\left(x_k\right)=F\left(x_{\ast }\right)+F^{\prime }\left(x_{\ast }\right)e_k+\frac{1}{2!}{F}^{″}\left(x_{\ast }\right)e_k^2+\frac{1}{3!}{F}^{‴}\left(x_{\ast }\right)e_k^3+O\left(e_k^4\right)\\ =F^{\prime }\left(x_{\ast }\right)\left[e_k+C_2{e}_k^2+C_3{e}_k^3\right]+O\left(e_k^4\right),\end{array}$ (3.5)

$F^{\prime }\left(x_k\right)=F^{\prime }\left(x_{\ast }\right)\left[I+2C_2{e}_k+3C_3{e}_k^2\right]+O\left(e_k^3\right)$ (3.6)

由(3.5)和(3.6)可得

$F^{\prime }{\left(x_k\right)}^{-1}F\left(x_k\right)=e_k-C_2{e}_k^2+2\left(C_2^2-C_3\right)e_k^3+O\left(e_k^4\right)$ (3.7)

则由(2.6)和(3.7)可得

$z_{k+1}=x_{\ast }+C_2{e}_k^2+2\left(C_2^2-C_3\right)e_k^3+O\left(e_k^4\right)$ (3.8)

对 $F\left(z_{k+1}\right)$ 在点 $x_{*}$ 运用Taylor公式展开，整理可得

$F^{\prime }\left(z_{k+1}\right)=F^{\prime }\left(x_{\ast }\right)\left[I+2C_2^2{e}_k^2+4\left(C_2{C}_3-C_2^3\right)e_k^3\right]+O\left(e_k^4\right)$ (3.9)

由(3.5)和(3.9)可得

$F^{\prime }\left(x_k\right)+F^{\prime }\left(z_{k+1}\right)=2F^{\prime }\left(x_{\ast }\right)\left(I+C_2{e}_k+\left(C_2^2+\frac{3}{2}{C}_3\right)e_k^2+2\left(C_2{C}_3-C_2^3+C_4\right)e_k^3\right)+O\left(e_k^4\right)$ (3.10)

对于两步算术平均牛顿迭代法AN：由于

${\left(\frac{1}{2}\left(F^{\prime }\left(x_k\right)+F^{\prime }\left(z_{k+1}\right)\right)\right)}^{-1}F\left(x_k\right)=e_k-\left(C_2^2+\frac{1}{2}{C}_3\right)e_k^3+O\left(e_k^4\right)$

所以

$\begin{array}{c}{x}_{k+1}=x_k-\frac{2F\left(x_k\right)}{F^{\prime }\left(x_k\right)+F^{\prime }\left(z_{k+1}\right)}\\ =x_k-\left(e_k-\left(C_2^2+\frac{1}{2}{C}_3\right)e_k^3+O\left(e_k^4\right)\right)\end{array}$

则

$e_{k+1}=e_k-\left(C_2^2+\frac{1}{2}{C}_3\right)e_k^3+O\left(e_k^4\right)$

于是两步算术平均牛顿迭代法(2.6)具有三阶收敛精度。

类似的，可证明两步调和平均牛顿迭代法(2.7)和两步中点牛顿迭代法(2.8)具有三阶收敛精度。

4. 数值实验

本节通过3个数值算例来验证两步算术平均牛顿迭代法(2.6)，两步调和平均牛顿迭代法(2.7)和两步中点牛顿迭代法(2.8)的3阶收敛。设定实验初始解的存在区间均为[−10, 10]，给出实验迭代的终止条件为：

1) 迭代次数最大为100次，即迭代次数 $n<100$ ；

2) 相邻两次迭代 $x_{n+1}$ 与 $x_n$ 差的2范数小于 $1\times {10}^{-5}$，即 ${\Vert x_{n+1}-x_n\Vert }_2<1\times {10}^{-5}$。

若迭代达到最大迭代次数100，且最后一次迭代结果与近似解相差过多，则认为该算法不收敛。

例4.1. 求非线性方程组

$\{\begin{array}{l}3x_1^5+x_1^7+x_2^4+20=0\hfill \\ x_1^3+2x_2^2+2e^{x_2}-40=0\hfill \end{array}$ (4.1)

的近似值。其中近似值[−1.664274278843389, 2.706644847598123]^T看作准确值。

例4.2. 求非线性方程组

$\{\begin{array}{l}{x}_1^2-10x_1+x_2^2+8=0\hfill \\ x_1{x}_2^2+x_1-10x_2+8=0\hfill \end{array}$ (4.2)

已知其精确解是[1, 1]^T

例4.3. 求非线性方程组

$\{\begin{array}{l}{x}_2{x}_3+x_2{x}_4+x_3{x}_4=0\hfill \\ x_1{x}_3+x_1{x}_4+x_3{x}_4=0\hfill \\ x_1{x}_2+x_1{x}_4+x_2{x}_4=0\hfill \\ x_1{x}_2+x_1{x}_3+x_2{x}_3-1=0\hfill \end{array}$ (4.3)

[0.5773502691896, 0.57735026918962, 0.57735026918962, −0.2886751345948]^T认为是该方程组的准确值。

表1和表2是算例4.1，算例4.2，算例4.3的数值结果。表1的第2列 $X_0$ 是通过布谷鸟搜索算法(2.1)~(2.2)求出的初值。第3列AN，HN，MN分别表示两步算术平均牛顿迭代法(2.6)，两步调和平均牛顿迭代法(2.7)和两步中点牛顿迭代法(2.8)。第4列是达到给定误差精度所需要的迭代次数n，以及在给定的迭代终止条件下收敛阶的近似值p，其中NC表示不收敛。第6列 ${\Vert x_{n+1}-x_n\Vert }_2$ 表示迭代终止时最后两次迭代结果差的2范数。从表1可以看出，对于例4.1和例4.2，HN方法计算效果不好。

Table 1. Results of numerical experiment

表1. 数值实验结果

表2给出了三种迭代格式在相同初值，达到相同迭代终止条件，所需要的计算机运算时间。(注：表1和表2中初值 $X_0$ 不相同)。从表2可以看出，MN方法用时最少，而HN方法用时最多。

Table 2. Computing time

表2. 运算时间

通过表1和表2可以看出，由布谷鸟搜索算法提供的值是可以作为牛顿迭代法初始值使用。通过对比这三个方法在求解非线性方程组的数值结果，我们发现两步中点牛顿迭代法(2.8)无论在精度上，还是运算时间上都比另两个方法具有一定的优势。

5. 小结

本文首先用布谷鸟搜索算法得到非线性方程组的一个值，该值作为牛顿迭代的初值，然后分别用算术平均牛顿迭代法，几何平均牛顿迭代法，中点牛顿迭代法，进行计算得到非线性方程组三阶精度的近似值。该方法避免了牛顿迭代初值选取的影响。通过3个数值例子对这三个格式，进行了精度和运行效率(电脑计算时间)的比较。对这三个方法更详细的理论分析和数值比较是下一步的工作。

基金项目

河南省高等学校重点科研项目：21B110003；河南科技大学实验技术开发基金项目：SY2021033；河南科技大学大学生研究训练计划(SRTP)：2020175；河南省自然科学基金面上项目：202300410156。

参考文献

NOTES

^*通讯作者。