基于奇异值分解的非线性系统流形解的数值计算

Numerical Computation of Manifold Solutions for Nonlinear Systems Based on Singular Value Decomposition

摘要: 本文利用边界矩阵技术，将非线性系统流形解的数值计算转化为增广非线性系统正则解的数值计算。利用雅可比矩阵的近似左零空间和右零空间构造增广矩阵，进而利用增广矩阵构造一族增广非线性系统，该族增广非线性系统中每个非线性系统均为给定非线性系统添加约束条件后的摄动系统，利用牛顿迭代更新增广非线性系统的解使得其逼近原非线性系统的流形解。

Abstract: In this paper, the numerical computation of manifold solutions of nonlinear systems is transformed into the numerical computation of regular solutions of augmented nonlinear systems using the boundary matrix technique. The approximate left-zero space and right-zero space of Jacobi matrix are used to construct the generalized matrix, and then the generalized matrix is used to construct a family of generalized nonlinear systems, and each nonlinear system in the family of generalized nonlinear systems is a regenerative system after the constraints are added to the given nonlinear system, and Newton iteration is utilized to update the solution of the generalized nonlinear system to make it approximate to the manifold solution of the original nonlinear system.

1. 引言

非线性系统的求解问题在科学计算领域占据基础性的核心地位。对于具有方阵结构的非线性系统，若在某一特定解处雅可比矩阵满秩且不存在奇异性，那么称该解为正则解。特别地，当非线性系统的解集构成一个流形时，这些解称为流形解。

流形解的求解方法已得到广泛研究。Allgower和Georg [1]提出的同伦延拓法通过构造连续的同伦路径，能够有效追踪非线性系统的多解问题，特别适用于解集构成流形的情况。Sommese和Wampler [2]利用数值代数几何方法，通过构造多项式系统的解空间，显著提高了流形解的求解精度。Tenenbaum等人[3]提出的流形学习方法通过提取非线性系统的低维流形结构，为流形解的降维求解提供了有效途径。Bates等人[4]研究了多项式系统的数值解流形，提出了基于代数几何的高效算法。Hauenstein和Sottile [5]利用数值同伦方法求解多项式系统的流形解，为流形解的全局求解提供了新思路。此外，Kuznetsov [6]基于分岔理论的方法通过分析系统参数变化时解的稳定性，揭示了流形解的分岔行为。O’Malley [7]引入的奇异摄动法通过分解系统的快慢尺度特性，能够求解具有多尺度行为的流形解。Raissi等人[8]将深度学习引入数值计算，通过训练神经网络逼近非线性系统的解，为高维流形解的求解提供了新的思路。

近年来，许多学者尝试使用牛顿法计算流形解。Griewank和Walther [9]提出自动微分技术，通过精确计算雅可比矩阵和高阶导数，提高了牛顿法在流形解附近的收敛性，但是计算量较大。Deuflhard [10]引入拟牛顿法，通过近似雅可比矩阵的逆矩阵，避免了直接处理奇异点带来的数值不稳定性，但存在收敛缓慢甚至不收敛的问题。Kelley [11]提出基于Krylov子空间的方法，通过迭代求解线性化系统的近似解，提高了大规模非线性系统中流形解的求解效率，但对初始向量的选择较为敏感。我们提出的算法解决了牛顿法在求解非线性系统流形解上的收敛性问题。

本文运用边界矩阵技术[12]，将一般系统的流形解问题转化为求解方形增广系统的正则解问题。利用雅可比矩阵的近似左零空间和右零空间构造增广矩阵并求逆，对传统牛顿迭代进行改进。经过理论推导证明了改进后的牛顿迭代具有局部二阶收敛性，并通过数值计算证明了该方法的有效性。

2. 符号和预备知识

列向量由小写黑体字母表示，如$b$ ，$x$ ，$y$ ，其中$0$ 代表零维向量。向量中的元素用普通小写字母表示，如$x,y$ 。我们分别用$ℝ$ ，$ℂ$ ，$\mathbb{Q}$ 和$\mathbb{N}$ 表示实数，复数，有理数和非负整数集合，矩阵由大写字母表示，如$J,A$ ，其中$O$ 和$I$ 分别表示零矩阵和单位矩阵。设$U\Sigma V^{\text{T}}$ 是$A$ 的奇异值分解，其中$U=\left[u_1,\cdots ,u_m\right]$ 是左奇异向量形成的酉矩阵，$V=\left[v_1,\cdots ,v_n\right]$ 是右奇异向量形成的酉矩阵。

如果$f:\Omega \subset {\mathbb{F}}^m\to {\mathbb{F}}^n$ 具有二阶连续导数，其中$\mathbb{F}=ℝ$ 或者$ℂ$ ，域$\Omega$ 是${\mathbb{F}}^m$ 的开子集，则称$f$ 是一个光滑映射。指定变量$x$ ，有$x\mapsto f\left(x\right)$ ，则$f$ 在任何特定$x_0$ 处的雅可比矩阵可表示为$J_f\left(x_0\right)$ 。

定义1：如果存在$\epsilon \in {ℝ}^{+}$ ，则$f\left(x\right)=0$ 的解$\widehat{x}$ 满足条件

$\left\{y\in {ℂ}^n:\Vert y-\widehat{x}\Vert <\epsilon \right\}\cap f^{-1}\left(0\right)=\left\{\widehat{x}\right\}$

则$\widehat{x}$ 是$f\left(x\right)=0$ 孤立解。

定义2：设$x_{*}$ 是$f\left(x\right)=0$ 的解，如果存在$x_{*}$ 的开邻域$\Delta$ ，使得$\Delta \cap f^{-1}\left(0\right)=\varphi \left(\Lambda \right)$ ，其中$z\mapsto \varphi \left(z\right)$ 是一个定义$\Delta$ 在内满足条件$\varphi \left(z_{*}\right)=x_{*}$ 的可微映射，且$rank\left({\varphi }_z\left(z_{*}\right)\right)=k$ ，则$x_{*}$ 作为$f\left(x\right)=0$ 零点其维数定义为

${\dim }_f\left(x_{*}\right)=\dim \left(range\left({\varphi }_z\left(z_{*}\right)\right)\right)\equiv rank\left({\varphi }_z\left(z_{*}\right)\right)=k$

3. 边界系统

对于给定的$x\in {ℝ}^n$ ，$J_f\left(x\right)$ 的奇异值分解表示为$J_f\left(x\right)=U\left(x\right)\Sigma \left(x\right){\left(V\left(x\right)\right)}^{\text{T}}$ ，其中 $U\left(x\right)=\left(U_1\left(x\right),U_2\left(x\right)\right)$ ，$U_2\left(x\right)\in {ℝ}^{m\times \left(m-k\right)}$ ，$V\left(x\right)=\left(V_1\left(x\right),V_2\left(x\right)\right)$ ，$V_2\left(x\right)\in {ℝ}^{n\times \left(n-k\right)}$ 。$U_2\left(x\right)$ 和$V_2\left(x\right)$ 分别是$J_f\left(x\right)$ 的左零和右零空间。如果近似解$\overline{x}$ 足够接近$x$ ，则有以下假设

假设1：假设有误差${\epsilon }_1\in {ℝ}^{+}$ ，满足

${\sigma }_1\left(\overline{x}\right)\ge \cdots \ge {\sigma }_k\left(\overline{x}\right)>{\epsilon }_1\ge {\sigma }_{k+1}\left(\overline{x}\right)\ge \cdots \ge {\sigma }_p\left(\overline{x}\right)$

那么有

$U_2\left(x\right)=\left(u_{k+1}\left(\overline{x}\right),\cdots ,u_m\left(\overline{x}\right)\right),V_2\left(x\right)=\left(v_{k+1}\left(\overline{x}\right),\cdots ,v_n\left(\overline{x}\right)\right)$

分别是$J_f\left(x\right)$ 的近似左零空间和右零空间。

定理1：设$J_f\left(x\right)$ 是光滑映射$f$ 在${ℝ}^m$ 的开域$\Omega$ 中任意$x$ 处的雅可比矩阵，假设在$x_{*}\in \Omega$ 处雅可比矩阵的秩是k，则矩阵

$J\left(x_{*}\right)={\left(\begin{array}{cc}{J}_f\left(x_{*}\right)& U_2\left(x_{*}\right)\\ V_2{\left(x_{*}\right)}^{\text{T}}& O\end{array}\right)}_{\left(m+n-k\right)\times \left(m+n-k\right)}$ (1)

是非奇异矩阵。

证明：假设矩阵(1)是奇异的，则存在非零向量$\left(\begin{array}{c}x\\ y\end{array}\right)$ ，使得

$\left(\begin{array}{cc}{J}_f\left(x_{*}\right)& U_2\left(x_{*}\right)\\ V_2{\left(x_{*}\right)}^{\text{T}}& O\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)$

根据奇异值分解的定义，有$V_2^{\text{T}}{V}_1=O_{n-k,k}$ ，因此存在向量$z$ 使得$x=V_{1}z$ ，故

$J_f\left(x_{*}\right)V_{1}z+U_2\left(x_{*}\right)y=0$

因为

$V_1^{\text{T}}{J}_f\left(x_{*}\right)=\left(\begin{array}{ccc}{\sigma }_1& & \\ & \ddots & \\ & & {\sigma }_k\end{array}\right)V_1^{\text{T}}$

所以

$V_1^{\text{T}}{J}_f\left(x_{*}\right)V_{1}z=\left(\begin{array}{ccc}{\sigma }_1& & \\ & \ddots & \\ & & {\sigma }_k\end{array}\right)V_1^{\text{T}}{V}_{1}z=\left(\begin{array}{ccc}{\sigma }_1& & \\ & \ddots & \\ & & {\sigma }_k\end{array}\right)z=\left(\begin{array}{c}0\\ ⋮\\ 0\end{array}\right)$

故$z=0$ ，与已知矛盾，所以矩阵(1)非奇异。 □

定理2：设$L$ 为Lipschitz常数，使得对于任何$x\in \Omega$ ，

$\Vert J_f\left(\overline{x}\right)-J_f\left(x\right)\Vert \le L\Vert \overline{x}-x\Vert$

如果

$\Vert \overline{x}-\widehat{x}\Vert <\frac{1}{L\Vert J{\left(\overline{x},\overline{x}\right)}^{-1}\Vert }$ (2)

那么存在一个开集$\overline{\Omega }\subset \Omega$ ，使得$\widehat{x},\overline{x}\in \overline{\Omega }$ ，并且对于任何$x\in \overline{\Omega }$ ，$J\left(x,\overline{x}\right)$ 是非奇异的。

证明：由于$\overline{x}$ 是足够接近$\widehat{x}$ 的近似解，所以存在一个开集$\overline{\Omega }\subset \Omega$ ，使得对于任何$x\in \overline{\Omega }$ ，有

$\Vert \overline{x}-x\Vert <\frac{1}{L\Vert J{\left(\overline{x},\overline{x}\right)}^{-1}\Vert }$

显然$\widehat{x},\overline{x}\in \overline{\Omega }$ 。因为

$\Vert J\left(\overline{x},\overline{x}\right)-J\left(x,\overline{x}\right)\Vert =\Vert J_f\left(\overline{x}\right)-J_f\left(x\right)\Vert \le L\Vert \overline{x}-x\Vert$

发现对于任何$x\in \overline{\Omega }$ ，有

$\Vert I_{m+n-k}-J{\left(\overline{x},\overline{x}\right)}^{-1}J\left(x,\overline{x}\right)\Vert \le \Vert J{\left(\overline{x},\overline{x}\right)}^{-1}\Vert \Vert J\left(\overline{x},\overline{x}\right)-J\left(x,\overline{x}\right)\Vert <1$

由定理2可知，对于任何$x\in \overline{\Omega }$ ，$J\left(x,\overline{x}\right)$ 是非奇异的。 □

4. 数值算法

对于非线性系统$f\left(x\right)=0$ ，设$J_f\left(x\right)=U\left(x\right)\Sigma \left(x\right)V{\left(x\right)}^{\text{T}}$ 为其对应雅可比矩阵的奇异值分解。对于$k\ge 1$ ，令$U\left(x\right)=\left(U_1\left(x\right),U_2\left(x\right)\right)$ ，其中$U_1\left(x\right)\in {ℝ}^{m\times k}$ ，$U_2\left(x\right)\in {ℝ}^{m\times \left(m-k\right)}$ 。令$V\left(x\right)=\left(V_1\left(x\right),V_2\left(x\right)\right)$ ，其中$V_1\left(x\right)\in {ℝ}^{n\times k}$ ，$V_2\left(x\right)\in {ℝ}^{n\times \left(n-k\right)}$ 。引入$\left(m-k\right)$ 维向量$y$ ，定义增广非线性系统

$g\left(x,y\right)=\{\begin{array}{l}f\left(x\right)+U_2\left(x\right)y=0\\ V_2{\left(x\right)}^{\text{T}}x=0\end{array}$ (3)

将求解系统$f\left(x\right)=0$ 的流形解问题转化为求解系统$g\left(x,y\right)=0$ 的正则解。

令$x_0$ 为系统$f\left(x\right)=0$ 的初始近似解，由于非线性系统(3)是高度非线性系统，本文利用交替迭代方法的思想，构造一系列非线性系统$g_k\left(x,y\right)=0$ ，其中

$g_1\left(x,y\right)=\{\begin{array}{l}f\left(x\right)+U_2\left(x_0\right)y=0\\ V_2{\left(x_0\right)}^{\text{T}}x=0\end{array}$

具体而言，设$x_k$ 是方程$g\left(x,y\right)=0$ 的第k步近似解，计算$f\left(x_k\right)$ 的雅可比矩阵$J_f\left(x_k\right)=U\left(x_k\right)\Sigma \left(x_k\right)V{\left(x_k\right)}^{\text{T}}$ ，得到$U_2\left(x_k\right)$ ，$V_2\left(x_k\right)$ ，则系统$g_k\left(x,y\right)$ 定义为

$\{\begin{array}{l}f\left(x\right)+U_2\left(x_k\right)y=0\\ V_2{\left(x_k\right)}^{\text{T}}x=0\end{array}$

利用牛顿迭代法求解系统$g_k\left(x,y\right)$ ，即

$\left(\begin{array}{c}{x}_{k+1}\\ y_{k+1}\end{array}\right)=\left(\begin{array}{c}{x}_k\\ y_k\end{array}\right)-{\left(\begin{array}{cc}{J}_f\left(x_k\right)& U_2\left(x_k\right)\\ V_2{\left(x_k\right)}^{\text{T}}& 0\end{array}\right)}^{-1}\left(\begin{array}{c}f\left(x_k\right)+U_2\left(x_k\right)y_k\\ V_2{\left(x_k\right)}^{\text{T}}{x}_k\end{array}\right)$ (4)

利用更新后的$x_{k+1}$ ，计算奇异值分解$J_f\left(x_{k+1}\right)=U\left(x_{k+1}\right)\Sigma \left(x_{k+1}\right)V{\left(x_{k+1}\right)}^{\text{T}}$ ，得到$U_2\left(x_{k+1}\right)$ ，$V_2\left(x_{k+1}\right)$ ，进而构造非线性系统$g_{k+1}\left(x,y\right)$ 。令$z_k=\left(x_k,y_k\right)$ ，增广后的雅可比矩阵表示为$J\left(z_k\right)$ ，则牛顿迭代(4)可以表示成

$z_{k+1}=z_k-J{\left(z_k\right)}^{-1}g\left(z_k\right)$ (5)

下面证明牛顿迭代(4)具有局部二阶收敛性。

定理3：设$g:{ℝ}^{m+n-k}\to {ℝ}^{m+n-k}$ 于开集$D\subseteq {ℝ}^{m+n-k}$ 上连续可微，$z_{*}$ 是$g_k\left(x,y\right)=0$ 的正则解，$J\left(z_{*}\right)$ 非奇异，且$J\left(z_k\right)$ 在$z_{*}$ 的某个邻域内满足Lipschitz条件，则存在$z_{*}$ 的一个邻域$\Omega \subseteq D$ ，使得对于任意初始点$z_0\in \Omega$ ，由牛顿迭代(4)产生的序列$\left\{z_k\right\}$ 收敛到$z_{*}$ ，且在$z_{*}$ 局部具有二阶收敛性，即当$\left\{z_k\right\}$ 收敛到$z_{*}$ 时，存在常数$C$ 和正整数$K$ ，使得对于$k\ge K$ ，满足不等式$\Vert z_{k+1}-z_{*}\Vert \le C{\Vert z_k-z_{*}\Vert }^2$ 。

证明：由式(5)可以得到

$z_{k+1}-z_{*}=z_k-z_{*}-J{\left(z_k\right)}^{-1}g\left(z_k\right)$ (6)

将$g\left(z_k\right)$ 进行泰勒展开得到

$g\left(z_k\right)=g\left(z_{*}\right)+J\left(z_{*}\right)\left(z_k-z_{*}\right)+O\left({\Vert z_k-z_{*}\Vert }^2\right)$ (7)

又$g\left(z_{*}\right)=0$ ，将式(7)代入式(6)有

$z_{k+1}-z_{*}=z_k-z_{*}-J{\left(z_k\right)}^{-1}\left(J\left(z_{*}\right)\left(z_k-z_{*}\right)+O\left({\Vert z_k-z_{*}\Vert }^2\right)\right)$ (8)

因为$J\left(z_{*}\right)$ 非奇异且满足Lipschitz条件，则有

$J{\left(z_k\right)}^{-1}J\left(z_{*}\right)=I+O\left(\Vert z_k-z_{*}\Vert \right)$

于是对式(8)化简可得

$z_{k+1}-z_{*}=O\left({\Vert z_k-z_{*}\Vert }^2\right)$

因此存在常数$C$ 和正整数$K$ ，使得对于$k\ge K$ ，满足

$\Vert z_{k+1}-z_{*}\Vert \le C{\Vert z_k-z_{*}\Vert }^2$ □

在特定的假设条件下，定理3证明了牛顿迭代(4)在正则解处的局部二阶收敛性。

5. 应用分析

非线性系统流形解的计算可以应用于数值问题中，比如亏损特征值的求解[13]。亏损特征值是指矩阵的几何重数小于代数重数的特征值。对于给定的矩阵$A$ 和近似特征值$\overline{\lambda }$ ，计算矩阵$A-\overline{\lambda }{I}_n$ 的奇异值分解，得到左零空间$U_2\left(\overline{\lambda }\right)$ 和右零空间$V_2\left(\overline{\lambda }\right)$ ，根据奇异值分布确定数值秩$q$ 。构造边界矩阵$B\left(\lambda ,\overline{\lambda }\right)$

$B\left(\lambda ,\overline{\lambda }\right)=\left(\begin{array}{cc}A-\lambda I_n& U_2\left(\overline{\lambda }\right)\\ V_2^{\text{T}}\left(\overline{\lambda }\right)& O_{q,q}\end{array}\right)$

使用边界矩阵的某些分量的解$F\left(\lambda \right)$ 来构造非线性系统$f\left(\lambda \right)=0$ ，具体来说

$f\left(\lambda \right)={\left(F_{1,1}\left(\lambda \right),\cdots ,F_{q,q}\left(\lambda \right)\right)}^{\text{T}}$

根据不同的$q$ 和$f^{\prime }\left(\widehat{\lambda }\right)$ 来构造增广系统$g_k\left(\lambda ,y\right)=0$ ，通过牛顿迭代

$\left({\lambda }_{k+1},y_{k+1}\right)=\left({\lambda }_k,y_k\right)-J{\left({\lambda }_k\right)}^{-1}{g}_k\left(\lambda ,y\right)$

求解出亏损特征值，其中$J\left({\lambda }_k\right)$ 是雅可比矩阵。通过奇异值分解和边界矩阵技术，可以将亏损特征值问题转化为增广非线性系统的流形解问题，利用牛顿迭代算法实现亏损特征值的求解。在优化理论与算法设计方面，许多复杂的优化问题同样也可转化为求解特定的非线性系统，利用文中算法进行求解。本文算法具有一定应用价值，但在大规模非线性系统的计算中，算法的计算复杂度较高，尤其是奇异值分解的计算成本较大，因此后续需要对算法进行研究完善，进一步拓展算法应用范围。

6. 数值计算

例1 考虑如下从${ℝ}^2\to {ℝ}^2$ 的映射$f$

$f\left(x_1,x_2\right)=\left(\begin{array}{c}{x}_1^3+x_1{x}_2^2-x_1+2x_1^2+2x_2^2-2\\ x_1^2{x}_2+x_2^3-x_2-3x_1^2-3x_2^2+3\end{array}\right)$

此代数系统的解是单位圆$x_1^2+x_2^2=1$ 。将非线性系统$f\left(x_1,x_2\right)=0$ 引入1维向量$y$ ，转化为高度非线性系统$g_k\left(x_1,x_2,y\right)$ ，利用牛顿迭代公式(4)进行求解。

从初始点$\left(x_1^{\left(0\right)},x_2^{\left(0\right)},y_0\right)=\left(1.8,0.6,0.1\right)$ 开始迭代，牛顿迭代最终在单位圆上收敛到点

$\left({\widehat{x}}_1,{\widehat{x}}_2,\widehat{y}\right)\approx \left(0.984871967607530,0.173283604016940,0.000000000000334\right)$

从初始点$\left(x_1^{\left(0\right)},x_2^{\left(0\right)},y_0\right)=\left(0.4,0.2,0.1\right)$ 开始，迭代收敛到点

$\left({\widehat{x}}_1,{\widehat{x}}_2,\widehat{y}\right)\approx \left(0.807454956739565,0.589929226972861,0.000000000000179\right)$

从初始点$\left(x_1^{\left(0\right)},x_2^{\left(0\right)},y_0\right)=\left(0.6,1.5,0.1\right)$ 开始，迭代收敛到点

$\left({\widehat{x}}_1,{\widehat{x}}_2,\widehat{y}\right)=\left(0.546169023159359,0.837674995532964,8.6875\text{e}-22\right)$

从初始点$\left(x_1^{\left(0\right)},x_2^{\left(0\right)},y_0\right)=\left(0.1,0.2,0.1\right)$ 开始，迭代收敛到点

$\left({\widehat{x}}_1,{\widehat{x}}_2,\widehat{y}\right)=\left(0.364515527645951,0.931197309975171,-2.4652\text{e}-32\right)$

可以发现，迭代结束后的引入向量$y$ 十分趋于0。经验算，迭代所求解是非线性系统$f\left(x_1,x_2\right)$ 的较好近似解。如图1所示，4个迭代序列都渐进地遵循解集的类法线方向朝向特定解。

例2

$f=\left\{x_1^2{x}_2+3x_1^2+x_2^3+3x_2^2-x_2-3,x_1^3-2x_1^2+x_1{x}_2^2-2x_2^2-x_1+2\right\}$

例3

$f=\left\{x_1^3\cos x_1+x_1^2{x}_2\cos x_2+x_1{x}_2^2\cos x_1+x_2^3\cos x_2-x_1\cos x_1-x_2\cos x_2,x_1^3+x_1^2{x}_2+x_1{x}_2^2+x_2^3-x_1-x_2\right\}$

例4

$f=\left\{x_1^3{x}_2+3x_1^2+x_1{x}_2^3+3x_2^2-x_1{x}_2-3,-x_1^4-x_1^2{x}_2^2-x_1^2+x_2^3-2x_2^2+x_1^2{x}_2-x_2+2\right\}$

易知例2~4代数系统的解包括单位圆$x_1^2+x_2^2=1$ 。将非线性系统引入一维向量$y$ ，转化为新的非线性系统$g\left(x_1,x_2,y\right)=0$ 。选取4个初始迭代点，迭代结果分别如表1~3所示，迭代序列效果分别如图2~4所示。从图中可以看出，迭代序列同样均遵循解集的类法线方向朝向特定解。

Figure 1. Iteration direction diagram of numerical Example 1

图1. 算例1迭代朝向图

Table 1. Example 2 iterative solution

表1. 例2迭代求解

$\left(x_1^{\left(0\right)},x_2^{\left(0\right)},y_0\right)$	$\left({\widehat{x}}_1,{\widehat{x}}_2\right)$	$\widehat{y}$
(1.8, 0.5, 0.1)	(0.878826438051687, 0.477141584628068)	3.5570e−33
(0.2, 0.4, 0.1)	(0.588654514186236, 0.808384724576217)	2.6402e−22
(0.6, 1.6, 0.1)	(0.269023744780287, 0.963133544605519)	7.0997e−30
(0.3, 0.2, 0.1)	(0.922612396045897, 0.385728358644330)	2.9969e−24

Figure 2. Iteration direction diagram of numerical Example 2

图2. 例2迭代朝向图

Table 2. Example 3 iterative solution

表2. 例3迭代求解

$\left(x_1^{\left(0\right)},x_2^{\left(0\right)},y_0\right)$	$\left({\widehat{x}}_1,{\widehat{x}}_2\right)$	$\widehat{y}$
(1.65, 0.6, 0.1)	(0.850885993243707, 0.525350384507016)	4.7980e−18
(0.4, 1.7, 0.1)	(0.457920680902070, 0.888993053967344)	3.8687e−18
(0.4, 0.5, 0.1)	(0.560400914857473, 0.828221476796459)	2.4652e−32
(0.6, 0.35, 0.1)	(0.923495769373285, 0.383608607762710)	0.00

Figure 3. Iteration direction diagram of numerical Example 3

图3. 例3迭代朝向图

Table 3. Example 4 iterative solution

表3. 例4迭代求解

$\left(x_1^{\left(0\right)},x_2^{\left(0\right)},y_0\right)$	$\left({\widehat{x}}_1,{\widehat{x}}_2\right)$	$\widehat{y}$
(1.8, 0.5, 0.1)	(0.983286424282257, 0.182065394356570)	−5.9714e−33
(0.2, 0.6, 0.1)	(0.241174587977960, 0.970481745378892)	−2.0251e−32
(0.55, 1.75, 0.1)	(0.556801891416083, 0.830645323658342)	−4.3141e−32
(0.5, 0.3, 0.1)	(0.829387581140966, 0.558673643775270)	0.00

例5 $f=\left\{\left(x_1-x_2\right){\left(x_1+x_2\right)}^2\left(x_1+x_2^2+x_2\right),\left(x_1-x_2\right){\left(x_1+x_2\right)}^2\left(x_1^2+x_2^2\right)\right\}$

例6 $f=\left\{x_1{x}_3+x_2^2-x_3,x_2+x_3,x_1{x}_2+x_3\right\}$

例7 $f=\left\{x_1^4\sin x_1+x_1^3{x}_2\sin x_1-x_1{x}_2\sin x_1-x_2^2\sin x_1,x_1^4\sin x_2+x_1^3{x}_2\sin x_2-x_1{x}_2\sin x_2-x_2^2\sin x_2\right\}$

例8 $f=\{x_1^5-x_1^2{x}_2+x_1^4{x}_2-x_1{x}_2^2,x_1^4{x}_2-x_1{x}_2^2+x_1^3{x}_2^2-x_2^3,x_1^4{x}_3-x_1{x}_2{x}_3+x_1^3{x}_2{x}_3-x_2^2{x}_3\}$

例9 $f=\left\{x_1^2+x_1{x}_2-x_1{x}_3-x_1-x_2+x_3,x_1{x}_2-x_2^2-x_2{x}_3-x_1+x_2+x_3,x_1{x}_3+x_2{x}_3-x_3^2-x_1-x_2+x_3\right\}$

Figure 4. Iteration direction diagram of numerical Example 4

图4. 例4迭代朝向图

使用MATLAB软件对例5~9进行迭代求解，结果如表4所示。从表中可以看出，引入向量$y$ 十分趋近于0，展示了改进牛顿迭代法对方形非线性系统正则解的良好收敛性。

Table 4. Iterative solution of nonlinear systems with regular solutions

表4. 非线性系统正则解迭代求解

system	$\left(x_0,y_0\right)$	$\widehat{x}$	$\widehat{y}$
Example 5	(0.01, 0.01, 0.01)	(−1.3346e−17, −6.6534e−18)	−3.7541e−84
Example 6	(1.01, 1.01, 1.01, 0.01)	(0.0398, 3.0949e−21, −1.3547e−21)	1.3547e−21
Example 7	(0.001, 0.001, 0.001)	(9.3354e−18, 2.2538e−17)	7.4967e−66
Example 8	(0.001, 0.001, 0.001, 0.001)	(1.2268e−17, 1.2268e−17, 1.2268e−17)	−8.0540e−66
Example 9	(1.01, 0.01, 1.01, 0.1)	(0.4680, −0.5466, 0.8948)	−0.5595

7. 结论

我们将非线性系统$f\left(x\right)=0$ 扩展至高维非线性系统$g_k\left(x,y\right)=0$ ，利用牛顿迭代公式

$\left(\begin{array}{c}{x}_{k+1}\\ y_{k+1}\end{array}\right)=\left(\begin{array}{c}{x}_k\\ y_k\end{array}\right)-{\left(\begin{array}{cc}{J}_f\left(x_k\right)& U_2\left(x_k\right)\\ V_2{\left(x_k\right)}^{\text{T}}& 0\end{array}\right)}^{-1}\left(\begin{array}{c}f\left(x_k\right)+U_2\left(x_k\right)y_k\\ V_2{\left(x_k\right)}^{\text{T}}{x}_k\end{array}\right)$

计算增广非线性系统的正则解，使其逼近原非线性系统的流形解。理论推导证明迭代在正则解附近局部二阶收敛，并通过数值计算证明了该方法的有效性。在未来的研究中，我们将探索更高效的矩阵分解技术和并行计算策略，以降低计算复杂度，提升算法在大规模问题上的适用性。

NOTES

^*通讯作者。