SIR：一种改进牛顿法的非线性方程组全局收敛求解器

||=−|| =; −+=+SoftwareX 7（2018）59原始软件出版物SIR-一种有效的方程组求解器Jan Scheffel，Kristoffer Lindvall*瑞典斯德哥尔摩KTH皇家理工学院电气工程与计算机科学学院聚变等离子体物理系ar t i cl e i nf o文章历史记录：2017年5月17日收到2018年1月15日收到修订版，2018年关键词：牛顿法求根方程求解器MATLABa b st ra ct半隐式根求解器（SIR）是一种求解非线性方程组全局收敛的迭代方法在这里，我们提出了MATLAB和MAPLE代码SIR，可以很容易地实现在任何应用中，线性或非线性方程组需要有效地解决。该代码采用最近开发的高效稀疏矩阵算法和改进的数值微分。SIR收敛是准单调的，并在实根附近接近二阶。全局收敛性通常优于牛顿此外，该算法不能降落在局部最小值上，如具有线搜索的牛顿法的情况版权所有©2018作者.由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前软件版本v1.0用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX-D-17-00039法律代码许可BSD 3条款，请参阅http://opensource.org/licenses/BSD-3-Clause使用MATLA B R的软件代码语言、工具和服务关于MAPLETM问题支持电子邮件kfli@kth.se1. 动机和意义代数方程组通常通过计算求解，无论是直接法还是迭代法。本文提出的半隐式根求解器（SIR [1]）是为了改善牛顿法的全局收敛性而提出的。线搜索[2]通常与牛顿法结合SIR的开发最初受到半隐式偏微分方程（PDE）算法[3-5 ]的启发最近，SIR算法已经通过使用当前的软件（MATLAB和MAPLE）版本，以及优化的编程，以适应这两个软件平台变得更加有效。这包括正确的矩阵处理和使用数字而不是符号运算。* 通讯作者。电子邮件地址：kfli@kth.se（K. Lindvall）。https://doi.org/10.1016/j.softx.2018.01.003首先，在第2节中简要概述了SIR。对于SIR算法的全面解释，建议读者参考[1]。第3节给出了一个应用示例。在第4节中可以找到详细描述算法的伪代码。2. 软件描述单方程f（x）的根0，其中f（x）X在迭代形式的重新公式化之后，由SIR发现，xi+1+βxi+1=βxi+β（xi），（1）其中β是实参数。当量（1）将具有与原始方程相同的根我们把它塑造成xi+1=α（xi−（xi））+（xi），（ 2）其中α β/（1β）是用于优化全局的参数，局部收敛引入Φ（xα）α（x（x））可以证明，收敛性要求对于根的邻域中的迭代xi成立<$Φ/<$x<1。牛顿因此，在根附近潜在地实现最大化的二阶收敛2352-7110/©2018作者。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softx|| ≡=−==−[客户端]1=260J。Scheffel，K.Lindvall/SoftwareX 7（2018）59(a) 牛顿法（b）牛顿法/线搜索。(c)使用（d）结构系统报告。图1.一、（a）牛顿法，（b）带线搜索的牛顿法，（c）MATLAB内置求解器fsolve和（d）带子迭代的SIR的收敛图。一组51乘51的均匀分布的x0 ∈ [-5，5]根据根求解器迭代的数量进行颜色映射。解出的两个方程是f1= x1−cos（x2）和f2=x2−3 cos（x1），其中f=（f 2+ f 2）/2。 解为x=[−0]。6843，2. 3245]。然而，牛顿SIR通过强制单调收敛来解决这个问题。技巧是在每次迭代i时选择合适的值。因此SIR迭代方程xi+1=αi（xi−（xi））+（xi），（3）使用我Ri′（xi）α1− α′（xi）.（四）其中，X′（x）d/ dx。典型地， 初始值为间隔0。五，零。99，然后SIR自动地将其减小到零，以在根附近实现二阶收敛由于保证了单调收敛，SIR将连续找到方程的所有实根x推广到方程组，SIR现在求解xi+1=Ai（xi−（xi））+（xi），（5）哪里A=I+（R-I）J-1，（6）（R）mn= δmnRm。（七）这里雅可比矩阵J有分量Jmn<$（xm<$m（x））/<$xn，δmn是克罗内克增量，I是单位矩阵。通常通过求解线性矩阵关系来获得A是最经济的Φ′=（A-I）J+I，（8）采用对角矩阵法和稀疏矩阵法，当求解多维方程时，不能保证严格的单调收敛，因为没有已知的有效程序（如1D情况下的根括号）来安全地保持从起点到根的方向。因此，SIR [ 1 ]中采用了“准单调性”的过程SIR还防止了某些陷阱。在Eq中可以看到一个明显的问题。（4）;存在A矩阵可能变得奇异的风险。解决这个问题的方法是子迭代，Rm朝向统一被选择。参见[1]进一步讨论增强收敛的措施。总而言之，在每次迭代中，SIR算法将Rm值减小到零以接近二阶收敛。如果非单调收敛在任何维度上变得明显，则使用通过增加Rm值朝向统一的局部子迭代3. 说明性示例SIR已被比较牛顿在某些情况下，=-∈=：→∈x=1 toImaxdo= −+==+−||==+x0=xn=1nnnnnnnnNKNKnn nnn=面=[客户端]J. Scheffel，K. Lindvall / SoftwareX 7（2018）59-6261具有优越的收敛特性，特别是当奇异性出现在A矩阵的计算众所周知，NL方法有时可能落在局部极小值上而不是落在方程的根在[1]中提供了这方面的一个由于迭代求解器对起始点x0非常敏感，因此该图使用颜色标记显示使用51 x 51均匀分布x0[-5，5]的集合的收敛质量。解出的两个方程是x1 cos（x2）和x2 3 cos（x1），见图。1.一、收敛图清楚地表明，即使对于这个看似简单的问题，SIR和NL格式的全局收敛性质也是复杂的可以看出，在子迭代模式（SIR-s）中，几乎所有起始点的收敛都是快速的，而在这里，标准SIR收敛大约三分之一的起始点。点NL方法不能匹配SIR-s的收敛算法1半隐式根求解器SIRSIR程序输入：一个向量函数RNRN 初步估计x0的RN.输出：向量x，是矩阵方程x（x）的根。参数：N-待解方程数、解算精度、Imax最大迭代次数、允许/省略子迭代的子标志、IS-最大子迭代次数、K -单调性检查点数、α c -A矩阵元素的最大允许幅值、MC-单调性检查参数、d Φ dx 0 -初始值、控制每次迭代时的Δ Φ i/Δ xi求解R=（A-I）J+I得到A比计算J-1要便宜。conv=|xi −xi−1|−|xi−1−xi−2|R：= dΦI由于过度降落在非零局部极小值，并执行J={d×0，：=单位矩阵就像SIR一样。MATLAB的内置根求解器fsolve的直接实现 同样，我们看到了在非零局部最小值上着陆的常见陷阱。 通过选择更合适的算法，可以使结果略有改善。使用默认设置尝试了所有算法、信任区域、信任区域狗腿和Levenberg-Marquardt方法。Levenberg-Marquardt方法采用高斯-牛顿方向和最速下降方向的线搜索技术，证明是具有最佳全局收敛性的算法。用Levenberg- Marquardt方法得到的fsolve结果如图所示. 1（c）.必须指出的是，这些结果绝不是最终的和决定性的，因为不是所有的选项设置都已经尝试过了。然而，他们强调了流行的根求解器的全局收敛性差的紧迫问题，以及如何SIR与子迭代可以解决这个问题。最近，我们采取了一系列措施，对于i<$（xm−<$m（x））/<$xn}x0的A I（R I）J−1x1A（x）如果sub和 max（conv）>0，则对于j=1到1，（见正文）S1=单调性检验S2=序列|A NJ|对于所有尺寸nfalse如果max（S2）ac和min（S1）≥Mc则返回Mc保证拟单调性索引k在区间xi，xi+1中遍历K个等距点。大的K值会导致代价高昂的函数求值;然而，我们发现K = 1的值适用于大多数问题，包括[1]中表1的问题。值MC=0将对应于严格单调的5. 影响SIR通常用作GWRM应用中非线性代数方程组的鲁棒求解器[1]。GWRM是一个62J. Scheffel，K.Lindvall/SoftwareX 7（2018）59时间谱偏微分方程求解器，已应用于线性和非线性偏微分方程。后者包括与数值天气预报有关的Burger该算法具有扩展的全局收敛性相比，牛顿因此，SIR具有在大量的计算物理领域中使用的潜力6. 结论SIR是近年来发展起来的一种求解一般非线性方程组的方法.全局收敛性往往优于牛顿SIR的代码也很简单;紧凑的MATLAB和MAPLE代码可以在GitHub存储库中找到（链接在元数据中）。引用[1] Scheffel J， Håkansson C.解非线性方程组，半隐式方法. 59.第五十九章http://dx.doi.org/10.1016/j.apnum。2009年5月2日。[2] 按WH，Teukolsky SA，Vetterling WT，Flannery BP。 数字配方;科学计算的艺术。剑桥大学出版社;2007年。[3] Gottlieb D ， Orszag A. 谱 方 法 的 数 值 分 析 ： 理 论 与 应 用 。 Philadelphia ：SIAM;1977.[4] HarnedD，Kerner W. 三维可压缩磁流体力学模拟的半隐式方法。J ComputPhys1985;60：62-75.[5] Harned D，Schnack DD.三维长时间尺度磁流体动力学计算的半隐式方法。JComput Phys 1986;65：57-70.[6] Scheffel J.初值问题的时间谱解。偏微分方程：理论、分析与应用。新星科学出版公司; 2011年。p. 1-49号。