基于正弦余弦算法的非线性方程组求根

雍龙泉

(陕西理工大学 数学与计算机科学学院, 陕西 汉中 723000)

考虑非线性方程组

(1)

记向量x=(x1,x2,…,xn)，向量函数F(x)=(f1(x),f2(x),…,fn(x))T，则式(1)等价于非线性方程

F(x)=0。

(2)

很多物理、化工、工业里复杂的问题都可以转化为非线性方程组，求解非线性方程组最基本的方法是牛顿法以及各类修正牛顿法等[1-8]。多步迭代的高阶牛顿迭代法(诸如5阶、7阶、8阶、9阶牛顿迭代法等[9-19])成为当前研究的热点。当然，收敛阶越高，计算代价也就越大。

实际计算过程中，若其雅克比矩阵奇异，即行列式det(F′(x)) =0，可以采用阻尼牛顿法或延拓法进行处理。目前多数牛顿迭代法都依赖于初始点的选取和F(x)的性态；此外针对存在多个根的非线性方程组，牛顿迭代法无法同时找到多个根。本文提出了一种求非线性方程组的新方法，具体做法是把非线性方程组转化为一个无约束优化，采用正弦余弦算法求解与之等价的无约束优化问题而获得原问题尽可能多的根。该方法的优点是无需计算F(x)的雅克比矩阵，且能找到多个根。

1 建立优化模型

为了求解非线性方程组，建立如下无约束优化模型

(3)

(3)式是一个标准的无约束优化问题，为了避免计算F(x)的雅克比矩阵，下面采用正弦余弦算法求解无约束优化(3)式。

2 正弦余弦算法

澳大利亚学者Mirjalili[20]在2016年提出的一种采用正弦函数和余弦函数对种群进行更新的智能优化算法，之后便形成了正弦余弦算法(Sine Cosine Algorithm，SCA)。

SCA 算法结构简单，容易实现，其最显著的特点是基于正弦函数和余弦函数值的变化来达到寻优目的。在SCA算法中，主要有r1、r2、r3、r44个参数。其中，最关键的参数是r1，控制算法从全局搜索到局部开发的转换。当r1的值较大时，算法倾向于全局探索；当r1的值较小时，算法偏向于局部开发。

下面给出SCA算法的主要步骤。

初始化算法参数: 种群规模N，空间维数D，控制参数a，最大迭代次数Tmax；

在可行域空间中随机初始化N个个体组成初始种群；t=1；

计算当前每个个体的适应值，并记录最优个体位置Pit；

while (t

fori=1 toNdo

forj=1 toDdo

根据式r1=a-at/Tmax计算r1的值；

随机产生r2∈U[0,2π],r3∈U[0,2],r4∈U[0,1]；

ifr4<0.5

else

end if

end for

end for

越界处理；

计算每个个体的适应值并更新整个种群的最优值；t=t+1；

end while

图1分别给出了a=2、Tmax=100与a=2、Tmax=500时r1sinr2与r1cosr2的图像。且只有当r1>1时，函数r1sinr2与r1cosr2值才有可能大于1或者小于-1；当r1≤1时，函数r1sinr2与r1cosr2值必在-1和1之间。

各级人影部门均按需设岗，市级人影业务人员分为：作业天气过程预报岗、作业条件潜力预报岗、作业条件监测岗、联合作业方案制定岗、作业跟踪指挥岗、登机作业指挥岗、登机技术保障岗、外场地面保障岗、作业信息管理岗、作业效益评估岗、云水资源评估岗、指挥平台运行保障岗、探测设备保障岗、网络运行管理岗、信息上报及效果评估岗；旗县级人影业务人员分为：综合管理岗、火箭作业岗、应急保障岗、高炮管理岗、信息上报及效果评估岗等；作业点人员分为现场指挥岗、作业实施岗、空域请示及信息上报岗。这样调整以后，做到了每个岗位均由专人负责，专职专岗、岗岗有人、岗岗有责。

根据SCA算法设计原理，算法先进行全局搜索再进行局部开发。如图2所示,当|r1sinr2|>1或者|r1cosr2|>1时，算法进行全局搜索；当|r1sinr2|<1或者|r1cosr2|<1时，算法进行局部开发。

刘勇等[21]对其参数r1展开研究，提出转换参数r1非线性递减的正弦余弦算法;曲良东等[22]提出了正弦算法。目前SCA算法已成功应用到工农业中，并获得了较好的结果。更多的SCA算法见文献[23-26]。下面应用SCA算法求解(3)式。

3 数值算例

给出几个非线性方程组，算例1—3的共同特点是其雅克比矩阵存在奇异点，故牛顿法难以进行。通过将其转化为无约束优化，采用SCA算法求解，程序用MATLABR2009a编写，参数设置N=30,a=2,Tmax=1000，算例1—3搜索空间为-100≤xi≤100,i=1,2；算例4搜索空间为-0.5≤xi≤1.5,i=1,2,…,n。为消除随机数对算法的影响，算法运行10次。

(a) a=2,Tmax=100 (b) a=2,Tmax=500图1 r1sinr2与r1cosr2的图像

图2 SCA算法搜索原理(r1=2,r3=1)

算例1 考虑非线性方程组

(a) 适应值收敛曲线 (b) 平均适应值收敛曲线图3 算例1运行10次的适应值和平均适应值收敛曲线

算例2 考虑非线性方程组

其根为x*=(3,0.5)T，x*=(81/32，-1/3)T。该方程的雅克比矩阵存在奇异点(说明：奇异点分布在x1=0,x2=1,x2=0,x2=-2 四条直线上)，因此牛顿法难以进行。图4给出了SCA算法运行10次的适应值收敛曲线和平均适应值收敛曲线图；表1中给出了SCA算法运行10次后得到的结果。

(a) 适应值收敛曲线 (b) 平均适应值收敛曲线图4 算例2运行10次的适应值和平均适应值收敛曲线

算例3 考虑非线性方程组

(a) 适应值收敛曲线 (b) 平均适应值收敛曲线图5 算例3运行10次的适应值和平均适应值收敛曲线

算例4 考虑非线性方程组

其根为x*=(1,1,…,1)T。图6分别给出了n=5时SCA算法运行10次的适应值收敛曲线和平均适应值收敛曲线图。

(a) 适应值收敛曲线 (b) 平均适应值收敛曲线图6 算例4运行10次的适应值和平均适应值收敛曲线

算例5 考虑非线性方程组

其解为x*=(1,1)T,x*=(0,0)T。该方程非光滑，因此牛顿法难以进行。图7分别给出了SCA算法运行10次的适应值收敛曲线和平均适应值收敛曲线图；表1中给出了SCA算法运行10次后得到的结果。

(a) 适应值收敛曲线 (b) 平均适应值收敛曲线图7 算例5运行10次的适应值和平均适应值收敛曲线

表1 算例1、2、3、5计算10次的结果

4 结束语

通过把非线性方程组转化为无约束优化，采用SCA算法求解，针对唯一根的非线性方程组，该方法能够收敛到其唯一根；针对具有多个根的非线性方程组，该方法能够找到尽可能多的根；且该算法无需计算雅克比矩阵，因此对雅克比矩阵是否存在奇异点无限制，适用于多数非线性方程组。计算结果表明, 该方法计算效率高, 对非线性方程组求解较为有效。如果要得到较高精度的数值解，可以通过增加迭代次数实现(如需源代码，请与作者联系)。