非线性方程组重根的可信验证方法

桑海风，万保成

（1.吉林大学 数学学院，长春130012；2.北华大学 数学与统计学院，吉林 吉林132013；3.吉林农业大学 信息技术学院，长春130118）

在火箭喷口受力、核磁共振机设计和数码机床的控制等高风险应用领域，计算结果的可靠性至关重要.Rump［1］给出了标准的可信验证方法，该方法将浮点运算用于严格证明，解决了非奇异问题的验证.如判断非线性方程组单根的存在性与唯一性问题，利用可信验证方法可以得到该问题的一个近似解及其相应的误差界，使得在近似解的误差界范围内必存在一个精确解.但验证非线性方程组多重根的存在性是一个奇异问题，因为对非线性方程组的系数做微小扰动，一个孤立奇异解重数就可能发生改变，而且在验证非线性方程组是否具有重根的过程中存在舍入误差，故验证非常困难.因此，对于非线性方程组多重根的验证，首先要将该方程组正则化，得到一个新方程组，使原来方程组的多重根成为新方程组单根的一部分，然后再利用标准可信验证方法验证新方程组的单根.Rump等［2］给出了Jacobi矩阵秩亏为1的非线性方程组二重根验证方法，该方法通过将原方程组增加一个光滑变量，并增加一个含n－1个变元、n个方程的方程组，得到一个含2n个变元、2n个方程的新方程组，使得新方程组的Jacobi矩阵非奇异.将验证原方程组Jacobi矩阵秩亏为1的奇异问题转化为验证新方程组Jacobi矩阵满秩的非奇异问题.Li等［3］给出了Jacobi矩阵秩亏为1的多重根的验证方法，该方法通过将原方程组增加μ－1个光滑变量，并增加μ－1个方程，得到一个含n＋μ－1个变元及n＋μ－1个方程的新系统，将验证原方程组Jacobi矩阵秩亏为1的奇异问题转化为验证新方程组Jacobi矩阵满秩的非奇异问题.

本文利用Shen等［4］的数值算法及区间算法，研究Jacobi矩阵秩亏为q的非线性方程组重根的可信性验证方法.该方法通过将原方程组增加q个光滑变量，并增加q个方程，得到一个Jacobi矩阵非奇异的新方程组，将验证原方程组多重根的奇异问题转化为验证新方程组单根的非奇异问题.基于此方法，提出了一种可信验证算法，该算法输出一个近似解及其相应的误差界，使得在近似解的误差界范围内必存在一个精确解.本文推广了文献［2］中Jacobi矩阵秩亏为1的二重根验证方法与文献［3］中Jacobi矩阵秩亏为1的多重根验证方法.

1 预备知识

记实数区间集合为I（ℝ），矩阵A 的第i 行为Ai，：＝（Ai，1，Ai，2，…，Ai，n），q 阶单位阵为Iq.令f：D⊂ℝn→ℝn为非线性系统，x＊∈ℝn为f（x）＝0的解，Jf（x＊）为f（x）在x＊处的Jacobi矩阵.如果其秩为r＝n－q，则称Jf（x＊）秩亏为q（1≤q≤n）.本文假设所讨论的非线性系统f满足下列条件：

（H1）f在x＊的邻域内满足C2－Lipschitz条件；

（H2）Jf（x＊）秩亏为q；

（H3）存在非零向量μ＊∈Null（（Jf（x＊））T）和Null（Jf（x＊））的一组基η＊＝｛η＊1，…，η＊q｝，使得q×q阶矩阵

非奇异.其中q×n阶矩阵（η＊）T（（μ＊）TJJf（x＊））的第k列为

对足够接近x＊的x，令η（x）和h（x）为

的唯一解，μ（x）和g（x）为

的唯一解.其中η（x）和h（x）分别为n×q和q×q阶矩阵，μ（x）与g（x）分别为n维和q维向量，α为随机选取的q维向量.

定义q×q阶矩阵

基于上述符号，定义边界系统

其中λ为q维变元.由F（x，λ）的Jacobi矩阵为

可得：

由上述讨论及引理3可得：

证明：由式（1）及条件（H3）（η＊为 Null（Jf（x＊））的一组基），有

进而

由式（5）知g（x＊）＝hT（x＊）·α＝0.于是有

故（x＊，0）为边界系统F（x，λ）＝0的根.再由引理3知（x＊，0）为边界系统F（x，λ）＝0的单根.

2 主要结果

问题1 设非线性系统f：D⊂ℝn→ℝn满足条件（H1）～（H3），给定其近似解∈ℝn，如何确定区间向量X，使得f（x）的精确解x＊∈＋X.

Rump［5］给出了系数阵为一般稠密矩阵的线性方程组求解的可信验证方法，其中系数阵可以是浮点矩阵，也可以是区间矩阵.

定理2［1］给定矩阵A，T∈ℝn×n，向量b∈ℝn×n，区间向量X⊂I（ℝn×n）.如果

成立，则矩阵A，T均非奇异，且A－1b∈Tb＋（I－TA）X.

实现线性方程组求解的可信验证算法函数是INTLAB函数包中的Verifylss函数［1］.对于系数阵为区间矩阵的线性方程组，Verifylss函数输出区间向量，该区间向量包含此区间线性方程组所有可能的解.

定理3［1］给定区间矩阵∈I（ℝn×n）和区间向量∈I（ℝn×n），如果函数Verifylss运行成功，则该函数计算得到的区间向量X⊂I（ℝn×n）满足

区间Newton法利用区间运算，通过在点Newton法的基础上引进区间变量，构成了Newton法的区间变形，使所得的迭代程序在每次迭代过程中都产生解的界限，从而不仅得到解的近似，同时还可得到相应的误差.区间Newton法的这个特点，显然是一般点迭代不具备的.Krawczyk［6］针对区间Newton程序在计算量方面的缺陷，提出一种改进的区间Newton法，建立了不需要计算区间矩阵逆的Krawczyk算子.

定义1 设f：ℝn→ℝ，若存在区间值映象F：I（ℝn）→I（ℝ），使得

成立，则称F为函数f的区间扩展.

定义2 设F：I（ℝn）→I（ℝ），X，Y∈I（ℝn）且满足X⊆Y，若有F（X）⊆F（Y），则称区间值函数F具有包含单调性.

定义3 设函数f：D⊂ℝn→ℝn连续可微，X∈I（D），JF是Jf的具包含单调性的区间扩展，Y为任意非奇异矩阵，称

为Krawczyk算子.

Krawczyk算子的性质：如果X∩K（y，X）＝Ø，则X中不包含非线性方程组f（x）＝0的解.这个性质给出了解存在与否的检验条件.在区间Newton迭代过程中，如果遇到这种情况，则迭代停止.

在Krawczyk算子K（y，X）中，除y可在X中任取外，非奇异矩阵T也是任意的.且实矩阵T的选择恰当与否直接关系到区间迭代法的收敛速度.取y＝m（X），T＝［m（JF（X））］－1，则相应的K（y，X）为

式（6）称为Krawczyk算子的Moore形式.由此可构造Krawczyk区间Newton迭代法，取X（0）＝X，

其中k＝0，1，2，….

定理4［7］设函数f：D⊂ℝn→ℝn连续可微，JF是Jf的具包含单调性的区间扩展，若初始区间向量X（0）⊆I（D）是一个n维立方体，K（X（0））满足K（X（0））⊂int（X（0）），则非线性方程组f（x）＝0在X（0）中有唯一解，且序列｛X（k）｝至少线性地收敛于该唯一解.

应用可信验证方法的前提是已知条件能在计算机上经过验证.在Krawczyk［6］给出了验证非线性系统解存在性的区间Newton法基础上，Rump［8］做了进一步的研究，改进区间Newton法使其能更便于实际应用.

定理5［8］设函数f：D⊂ℝn→ℝn，其中f＝（f1，…，fn）∈C1.给定向量∈ℝn，区间向量X∈I（ℝn），且0∈X，矩阵T∈ℝn×n，且给定的区间矩阵M∈I（ℝn×n）满足条件：

如果

成立，则存在唯一的向量x＊∈＋X，使得f（x＊）＝0，并且每个矩阵∈M都是非奇异的，其中int（X）表示区间X的内部，▽表示梯度.

本文将区间Newton算法应用于文献［4］中非线性方程组的数值解法，以达到数值计算结果的可信性验证.令X∈I（ℝn）且0∈X，（）∈ℝn＋q为F（x，λ）＝0的近似解，T为JF（）的近似逆矩阵.首先要找到区间矩阵N∈I（ℝ（n＋q）×（n＋q））满足

由式（4）知

利用区间Newton算子的性质，如果

基于上述理论，设计算法如下.

算法1

2）初始化：利用区间转换函数intval，得到初始区间向量

返回（X，Λ），算法终止；

由定理5，可得下述命题.

命题1 如果算法1成功返回的包含区间＋X和的包含区间＋Λ，则必存在唯一的向量使得（x＊，0）为F（x，λ）＝0的精确解.进而存在唯一的向量使得x＊为f（x）＝0的精确解.

证明：由数值Newton算法可计算出边界系统F（x，λ）＝0的数值迭代增量‖Δ（x（k），λ（k））‖＜ε2的近似解如果算法1成功返回的包含区间和的包含区间，则由定理5可知，必存在唯一的向量使得F（x＊，0）＝0.再由边界系统的定义知

3 数值算例

本文数值实验基于 Windows7操作系统，软件分别是 MAPLE 15（Digits∶＝14）和 MATLAB R2011a（INTLAB V6）.在MAPLE和MATLAB中执行可信验证算法1，可计算出非线性系统的近似解及相应的误差界，使得在近似解的误差界范围内必存在一个精确解.

输出：

例1的解x＊＝（0，0）为孤立根.

例2 考虑非线性系统

输出：

例2的解x＊＝（1，1，0，0，0）为孤立根.

［1］Rump S M.Verification Methods：Rigorous Results Using Floating－Point Arithmetic［J］.Acta Numerica，2010，19：287－449.

［2］Rump S M，Graillat S.Verified Error Bounds for Multiple Roots of Systems of Nonlinear Equations ［J］.Numerical Algorithms，2010，54（3）：359－377.

［3］LI Nan，ZHI Lihong.Verified Error Bounds for Isolated Singular Solutions of Polynomial Systems：Case of Breadth One［J］.Theoretical Computer Science，2013，479（1）：163－173.

［4］SHEN Yunqiu，Ypma T J.Newton’s Method for Singular Nonlinear Equations Using Approximate Left and Right Nullspaces of the Jacobian［J］.Applied Numerical Mathematics，2005，54（2）：256－265.

［5］Rump S M.Kleine Fehlerschranken bei Matrixproblemen［D］.Karlsruhe：Universität Karlsruhe，1980.

［6］Krawczyk R.Newton－Algorithmen zur Bestimmung von Nullstellen mit Fehlerschranken［J］.Computing，1969，4（3）：187－201.

［7］Moore R E.A Computational Test for Convergence of Iterative Methods for Nonlinear Systems ［J］.SIAM Journal on Numerical Analysis，1978，15（6）：1194－1196.

［8］Rump S M.Solving Algebraic Problems with High Accuracy ［C］／／Proceedings of the Symposium on a New Approach to Scientific Computation.San Diego：Academic Press，1983：51－120.