本文讨论如下弱非线性互补问题, 即寻求向量$ u\in {{\Bbb R}} ^n $使得

(1.1) $ \begin{eqnarray} \left\{ \begin{array}{l} u\geq0, \\ w:=Au+\Phi(u)\geq0, \\ w^Tu=0, \end{array} \right. \end{eqnarray} $

其中$ A\in {{\Bbb R}} ^{n\times n} $是给定的大规模稀疏矩阵, $ \Phi(u):{{\Bbb R}} ^n\rightarrow {{\Bbb R}} ^n $为非线性Lipschitz连续函数.

弱非线性互补问题在科学与工程领域有着广泛的应用. 例如, 运筹学, 均衡理论, 工程设计等. 当$ (1.1) $式中的$ \Phi(u)=q $, 弱非线性互补问题退化为线性互补问题

(1.2) $ \begin{eqnarray} \left\{ \begin{array}{l} u\geq0, \\ w:=Au+q\geq0, \\ w^Tu=0, \end{array} \right. \end{eqnarray} $

它常出现于工程或经济等各类实际问题, 如双矩阵博弈的纳什均衡, 径向轴承的自由边值问题等^{[11, 19, 21]}. 在文献[19] 中, Lemke通过研究双矩阵均衡点等问题引出线性互补问题. 关于该问题的算法构造与求解也一直是优化或数值代数领域的研究热点, 特别当矩阵$ A $是大规模且稀疏的情形, 研究成果也颇为丰硕. 例如, 求解线性互补问题的投影逐次超松弛迭代算法^[12]以及一般的不动点迭代方法^[24]. 近些年, 白等人利用矩阵分裂思想并通过隐不动点方程建立了一系列基于模的矩阵分裂迭代格式来求解线性互补问题^{[1, 5, 7, 9]}. 进而, 他们又提出了有效的多分裂迭代格式来求解线性互补问题^{[1-4, 6, 8]}. 后来, 基于这些分裂迭代格式又发展了一些更有效的加速技巧^{[27, 28]}.

对于弱非线性互补问题, 孙等人给出了一个单调的半光滑牛顿法^[23]. 熟知, 半光滑(光滑)牛顿法是非常有效的一类求解非光滑函数的方法, 它常用于互补问题或变分不等式等^{[10, 18, 22]}. 最近, 黄等人^[16]也提出了一个基于模的分裂算法来求解弱非线性互补问题$ (1.1) $, 该算法被证明是可行且高效的.

鉴于以上一些研究成果, 通过引入一个光滑函数及巧妙的等价变换, 我们提出了一个具有高阶收敛性的加速牛顿迭代法来求解一类大规模弱非线性互补问题, 其思想源于牛顿法至少具有二阶收敛性的优势. 数值实验部分将该算法与当前一些有效算法做了详细比较, 说明了所提出的算法是可行性且有效的.

为分析方便, 通篇文章给出如下一些记号: 设$ {\Bbb N}_k=\{1, 2, \cdots, k\} $表示含有$ k $个正整数的集合. 对任意$ x\in {{\Bbb R}} ^{n} $, $ \|x\| $表示欧氏范数. 给定两个$ n\times m $实矩阵$ A=(a_{ij}) $和$ B=(b_{ij}) $, 我们记$ A\geq B $ (或$ A>B $), 若$ a_{ij}\geq b_{ij} $ (或$ a_{ij}>b_{ij} $) 对于所有的$ i\in {\Bbb N}_{n} $及$ j\in {\Bbb N}_{m} $都成立. 符号$ |A| $及$ \rho(A) $分别表示$ A\in {{\Bbb R}} ^{n\times n} $的绝对值和谱半径. 对于可微函数$ F(x) $, $ F'(x) $表示其Jacobian矩阵. $ A^{-1} $表示方阵$ A $的逆矩阵. 此外, $ {\rm Diag}\{a_1, a_2, \cdots, a_n\} $代表对角矩阵, 其中$ a_i (i\in {\Bbb N}_n) $表示其主对角线上的元素.

本文剩余章节组织如下. 在第2节, 提出了一个求解弱非线性互补问题的带高阶收敛性的加速牛顿迭代法. 第3节详细分析了算法的收敛速率. 第4节给出了包括源于自由边值问题离散的弱非线性互补问题的一些数值结果, 充分验证了所构造的算法的有效性和可行性. 最后是结论部分.

这一节, 主要提出一个具有高阶收敛性的加速牛顿迭代法. 首先, 给出几个重要的引理和结论, 将弱非线性互补问题进行合理巧妙地转化, 从而能有助于利用所构造的算法进行有效求解, 并为下一节的收敛性分析做一些必要准备.

引理2.1  问题$ (1.1) $等价于如下非线性方程组

(2.1) $ \begin{equation} (I+A)x+(I-A)|x|-\Phi(|x|-x)=0, \end{equation} $

其中$ A\in {{\Bbb R}} ^{n\times n} $是给定的实矩阵, $ \Phi(u):{{\Bbb R}} ^n\rightarrow {{\Bbb R}} ^n $为非线性Lipschitz连续函数. $ I $为相应维数的单位矩阵. $ x \in {{\Bbb R}} ^n $为未知向量.

证  首先假设$ x $为非线性方程组$ (2.1) $的解, 那么

(2.2) $ \begin{equation} |x|+x=A(|x|-x)+\Phi(|x|-x). \end{equation} $

令

(2.3) $ \begin{equation} w:=|x|+x, u:=|x|-x. \end{equation} $

由$ (2.2) $及$ (2.3) $式可得

$ w=Au+\Phi(u), w\geq 0, u\geq 0, w^Tu=0, $

这意味着$ u $是问题$ (1.1) $的解.

另一方面, 假设$ u $是问题$ (1.1) $的解. 显然

$ w\geq 0, u\geq 0, w^Tu=0, $

其中$ w=Au+\Phi(u) $.

类似可得

$ w:=|x|+x, u:=|x|-x. $

从而

$ x=\frac{1}{2}(w-u) $

即为非线性方程组$ (2.1) $的解.

接下来, 引入一个光滑函数$ F: {{\Bbb R}} ^{n+1}\rightarrow {{\Bbb R}} ^{n} $

(2.4) $ \begin{eqnarray} F(z):=(I-A)\sqrt{x^2+\varepsilon^2 {\bf e}}+(I+A)x-\Phi(\sqrt{x^2+\varepsilon^2 {\bf e}}-x), \end{eqnarray} $

其中$ z=(\varepsilon, x^T)^T \in {{\Bbb R}} ^{n+1} $, $ \varepsilon >0 $, $ {\bf e}=(1, 1, \cdots, 1)^T\in {{\Bbb R}} ^n $,

(2.5) $ \begin{eqnarray} \sqrt{x^2+\varepsilon^2 {\bf e}}:=\Big[\sqrt{x_{1}^{2}+\varepsilon^2}, \sqrt{x_{2}^{2}+\varepsilon^2}, \cdots, \sqrt{x_{n}^{2}+\varepsilon^2}\Big]^T\in {{\Bbb R}} ^n. \end{eqnarray} $

如此一来, 可定义如下非线性可微函数

(2.6) $ \begin{eqnarray} \Psi(z):= \left( \begin{array}{cccccc} \varepsilon\\ F(z)\\ \end{array} \right)\in {{\Bbb R}} ^{n+1}. \end{eqnarray} $

引理2.2  问题$ (1.1) $等价于如下非线性方程组

(2.7) $ \begin{eqnarray} \Psi(z)=0, \end{eqnarray} $

其中$ z=(\varepsilon, x^T)^T \in {{\Bbb R}} ^{n+1} $, $ \Psi(z) $由$ (2.6) $式所定义.

证  结论容易由公式(2.4) 及引理$ 2.1 $推得. 此处不再赘述.

引理2.3  非线性函数$ (2.6) $的Jocabian矩阵具有如下形式

(2.8) $ \begin{eqnarray} \Psi'(z)= \left( \begin{array}{cccccc} 1{\quad} &{\bf 0}\\ (I-A-\Phi'_{\widetilde{z}}(\widetilde{z}))g_\varepsilon{\quad} &(I-A-\Phi'_ {\widetilde{z}}(\widetilde{z}))D_x+\Phi'_{\widetilde{z}}(\widetilde{z})+A+I \end{array} \right)\in{{\Bbb R}} ^{(n+1)\times(n+1)}, \end{eqnarray} $

其中, $ g_\varepsilon=(g_1, g_2, \cdots, g_n)^T\in{{\Bbb R}} ^{n}, $$ g_i=\frac{\varepsilon} {\sqrt{x_{i}^{2}+\varepsilon^2}}, $$ D_x={\rm Diag}(d_1, d_2, \cdots, d_n)\in{{\Bbb R}} ^{n\times n}, $$ d_i=\frac{x_{i}}{\sqrt{x_{i}^{2}+\varepsilon^2}} $$ ( i\in {\Bbb N}_n) $且$ \widetilde{z}:=\sqrt{x^2+\varepsilon^2 {\bf e}}-x $, $ {\bf 0}\in {{\Bbb R}} ^{n} $为零向量, $ \Phi(u) $为变量$ u $的某种非线性函数.

证  式(2.6) 两边对变量$ z $求导可得

(2.9) $ \begin{eqnarray} \Psi'(z)= \left( \begin{array}{cccccc} 1{\quad} &{\bf 0}\\ F'_\varepsilon{\quad} &F'_x \end{array} \right). \end{eqnarray} $

由(2.4) 式以及(2.5) 式可得

(2.10) $ \begin{eqnarray} F'_\varepsilon&=&(I-A)(\frac{\varepsilon}{\sqrt{x_{1}^{2}+\varepsilon^2}}, \frac{\varepsilon}{\sqrt{x_{2}^{2}+\varepsilon^2}}, \cdots, \frac{\varepsilon}{\sqrt{x_{n}^{2}+\varepsilon^2}})^T{}\\ &&-\Phi'_{\widetilde{z}}(\widetilde{z})(\frac{\varepsilon}{\sqrt{x_{1}^{2}+\varepsilon^2}}, \frac{\varepsilon}{\sqrt{x_{2}^{2}+\varepsilon^2}}, \cdots, \frac{\varepsilon}{\sqrt{x_{n}^{2}+\varepsilon^2}})^T{}\\ &=&(I-A-\Phi'_{\widetilde{z}}(\widetilde{z}))g_\varepsilon. \end{eqnarray} $

由于$ \widetilde{z}'_x=D_x-I, $则

(2.11) $ \begin{eqnarray} F'_x&=&(I-A)D_x+A+I-\Phi'_{\widetilde{z}}(\widetilde{z})\cdot\widetilde{z}'_x{}\\ &=&(I-A-\Phi'_{\widetilde{z}}(\widetilde{z}))D_x+\Phi'_{\widetilde{z}}(\widetilde{z})+I+A. \end{eqnarray} $

下面, 给出求解非线性光滑函数$ (2.7) $的一个具有高阶收敛性的加速牛顿型迭代算法. 具体算法描述如下.

算法2.1 (高阶牛顿迭代法(HONM))

步1   给定初值点$ z^0=(\varepsilon_0, (x^{0})^{T})^T, $适当维数的矩阵$ A $, 向量$ q $以及参数$ \sigma_1, \sigma_2\in(0, 1) $, $ m\geq2 $. 令$ k:=0 $.

步2   计算$ (2.7) $的$ \Psi(z^k) $, Jacobian矩阵$ \Psi'(z^k) $及矩阵$ A_k:=(\Psi'(z^k))^{-1}. $

步3   令$ z^{k, 1}=z^{k} $, 置$ j:=1. $

步4   计算$ \Psi(z^{k, j}) $, 更新向量序列

$ \begin{eqnarray*} z^{k, j+1}=z^{k, j}-A_k b, \end{eqnarray*} $

再计算$ \Psi(z^{k, j+1}) $, 其中$ b=\Psi(z^{k, j}). $

步5   令$ j:=j+1, z^{k, j}:=z^{k, j+1}, \Psi(z^{k, j}):=\Psi(z^{k, j+1}), y:=A_kb. $当$ j=m $, 返回步$ 6 $, 否则返回步$ 4 $.

步6   当$ \|y\|<\sigma_1 $或$ \|\Psi(z^{k, m})\|<\sigma_2, $令$ z^*=z^{k, m} $, 否则$ k:=k+1 $, 返回步3.

事实上, 高阶牛顿法(HONM) 亦可写成下列具体迭代格式

(2.12) $ \begin{eqnarray} \left\{ \begin{array}{l} z^{k, 0}=z^k, \\ z^{k, j}=z^{k, j-1}-(\Psi'(z^k))^{-1}\Psi(z^{k, j-1}), j=1, 2\cdots, m, \\ z^{k+1}=z^{k, m}, k=0, 1, 2, \cdots. \end{array} \right. \end{eqnarray} $

注2.1  根据引理$ 2.2 $, 我们易知由算法$ 2.1 $产生的解$ z^* $亦是问题$ (1.1) $的解$ u^* $.

注2.2  参数$ \varepsilon_k $可设置为$ \varepsilon_k:=\varepsilon_{k-1}^m. $由于算法$ 2.1 $中的$ m $为大于等于$ 2 $的常数, 因此序列$ \{\varepsilon_k\}_{0}^{\infty} $单调下降并趋向于$ 0 $.

注2.3  若取$ m=1 $, 则算法HONM将退化为一般的牛顿迭代法.

本节我们着重讨论高阶牛顿迭代法(HONM) 的收敛性能. 首先给出一个定义及一些重要引理.

定义3.1^[17]  令$ F: D\subset {{\Bbb R}} ^{n}\rightarrow {{\Bbb R}} ^{n}, $$ x^*\in D $为方程组$ F(x)=0 $的解. $ S\subset D $为解集的子集. 对任意的初始点$ x^0\in S, $若迭代序列$ \{x^k, k=0, 1, \cdots\} $是适定的并收敛于$ x^* $, 则称$ x^* $为迭代序列的吸引点.

鉴于经典的牛顿迭代法具有二阶收敛速度, 我们有下面的结论, 具体可参看文献[17] 以及其中涉及到的相关参考文献.

引理3.1^[17]  假设$ F: D\subset {{\Bbb R}} ^{n}\rightarrow {{\Bbb R}} ^{n} $在开区间$ S_0\in D $为Fréchet可微函数, 且$ F'(x^*) $非奇异, 其中$ x^* $为方程组$ F(x)=0 $的解. 那么, 映射$ G(x)=x-F'(x)^{-1}F(x) $在闭区域$ S=\bar{S}(x^*, \delta)\subset S_0 $上是适定的. 进一步, 若不等式

(3.1) $ \begin{eqnarray} \|F'(x)-F'(x^*)\|\leq\beta\|x-x^*\| \end{eqnarray} $

成立, 则牛顿迭代法至少二阶收敛, 其中常数$ \beta>0 $, $ x\in S $.

引理3.2^[17]  假设$ F: D\subset {{\Bbb R}} ^{n}\rightarrow {{\Bbb R}} ^{n} $在不动点$ x^*\in int(D) $上Fréchet可微且其谱半径

(3.2) $ \begin{eqnarray} \rho(F'(x^*))=\sigma<1. \end{eqnarray} $

则对任意初始点$ x^0\in S $, $ x^* $是迭代序列$ x^{k+1}=F(x^k) (k=0, 1, \cdots) $在开区域$ S=S(x^*, \delta)\subset D $上的吸引点.

引理3.3^[17]  假设$ F: D\subset {{\Bbb R}} ^{n}\rightarrow {{\Bbb R}} ^{m} $在凸集$ D_0\subset D $连续可微且

(3.3) $ \begin{eqnarray} \|F'(u)-F'(v)\|\leq\beta\|u-v\|^p, u, v\in D_0. \end{eqnarray} $

则

(3.4) $ \begin{eqnarray} \|F(y)-F(x)-F'(x)(y-x)\|\leq\frac{\beta}{1+p}\|y-x\|^{1+p}, x, y\in D_0 \end{eqnarray} $

成立, 其中常数$ \beta\geq0, p\geq0. $

下面将给出高阶加速牛顿迭代法的收敛性, 详细有如下定理.

定理3.1  假设$ \Psi:D\subset {{\Bbb R}} ^{n+1}\rightarrow {{\Bbb R}} ^{n+1} $在$ S(z^*, \delta)\in D $为Fréchet可微, $ \Psi'(z^*) $非奇异, 且存在$ \beta>0 $, 对任意的$ z\in S $有

(3.5) $ \begin{eqnarray} \|\Psi'(z)-\Psi'(z^*)\|\leq\beta\|z-z^*\| \end{eqnarray} $

成立, 则$ z^* $为算法$ 2.1 $产生的迭代序列$ \{z^k\}_{0}^{\infty} $的吸引点, 且

(3.6) $ \begin{eqnarray} \|z^{k+1}-z^*\|\leq L_2\|z^{k}-z^*\|^{m+1}, \end{eqnarray} $

其中$ z^* $为$ \Psi(z)=0 $的解, $ L_2 $为不依赖于迭代次数$ k $的常数.

证  首先, 考虑当$ m=2 $时的迭代格式(2.12). 即

(3.7) $ \begin{eqnarray} z^{k+1}=z^{k}-(\Psi'(z^k))^{-1}\Big[\Psi(z^k)-\Psi\big(z^k-(\Psi'(z^k))^{-1}\Psi(z^k)\big)\Big], k=0, 1, \cdots. \end{eqnarray} $

根据引理$ 3.1 $, $ P(z):=z-(\Psi'(z))^{-1}\Psi(z) $在$ S_1:=\bar{S}(z^*, \delta_1)\subset S $上是适定的且

(3.8) $ \begin{eqnarray} \|P(z)-z^*\|\leq \xi\|z-z^*\|^{2}, z\in S_1, \end{eqnarray} $

其中$ \xi>0 $. 因此, 映射$ M(z)=P(z)-(\Psi'(z))^{-1}\Psi(P(z)) $在闭区域$ S_2=\bar{S}(z^*, \delta_2)\subset S_1 $为适定的, 其中$ \delta_2\leq\frac{\delta_1}{\xi} $. 注意到$ \Psi(z^*)=0 $, 因此$ P(z^*)=z^* $. 经过适当计算可得

$ \begin{eqnarray*} P'(z) &=&I+(\Psi'(z))^{-2}\Psi''(z)\Psi(z)-(\Psi'(z))^{-1}\Psi'(z)\\ &=&(\Psi'(z))^{-2}\Psi''(z)\Psi(z), \end{eqnarray*} $

因此$ P'(z^*)=0. $进一步有,

$ \begin{eqnarray*} M'(z^*) &=&P'(z^*)+(\Psi'(z^*))^{-2}\Psi''(z^*)\Psi(P(z^*)) -(\Psi'(z^*))^{-1}\Psi'(P(z^*))P'(z^*)\\ &=&0, \end{eqnarray*} $

其中最后一个等式源于$ \Psi(P(z^*))=\Psi(z^*)=0 $的事实. 这意味着$ \rho(M'(z^*))=0<1. $那么, 由引理$ 3.2 $即可得$ z^* $为$ (3.7) $式的吸引点.

另外, 也注意到$ \Psi'(z^*) $的非奇异性, 因此对任意的$ z\in S_2 $, $ \|(\Psi'(z))^{-1}\|\leq\zeta $. 再由引理$ 3.3 $以及假设得

(3.9) $ \begin{eqnarray} \|M(z)-z^*\|&\leq&\big\|(\Psi'(z))^{-1}\big\|\big\|\Psi'(z)[P(z)-z^*]-\Psi\big(P(z)\big)\big\| &\leq&\zeta\Big[\big\|\Psi\big(P(z)\big)-\Psi(z^*)-\Psi'(z^*)(P(z)-z^*)\big\| \\ &&+\big\|\big(\Psi'(z^*)-\Psi'(z)\big)(P(z)-z^*)\big\|\Big] \\ &\leq&\zeta\Big[\frac{1}{2}\beta\big\|P(z)-z^*\big\|^2+\beta\|z-z^*\|\|P(z)-z^*\|\Big]\\ &\leq&\beta\zeta\Big[\frac{1}{2}\xi^2\|z-z^*\|+\xi\Big]\|z-z^*\|^3\\\nonumber &\leq&\beta\zeta\xi\Big[\frac{1}{2}\xi\delta_2+1\Big]\|z-z^*\|^3. \nonumber\end{eqnarray} $

令$ L_0:=\beta\zeta\xi\Big[\frac{1}{2}\xi\delta_2+1\Big] $, 显然它不依赖于迭代数$ k $.

若在不等式$ (3.9) $的左边取$ z^{k+1}=M(z^k) $, 则有

(3.10) $ \begin{eqnarray} \|z^{k+1}-z^*\|\leq L_0\|z^{k}-z^*\|^3, \end{eqnarray} $

这意味着对于$ m=2 $, 不等式$ (3.6) $是成立.

下面将用归纳法来证明迭代格式$ (2.12) $至少有$ m+1 $阶收敛速率.

由$ (3.9) $式可得

(3.11) $ \begin{eqnarray} \|z^{k, 2}-z^*\|\leq L_0\|z^{k}-z^*\|^3=O(\|z^{k}-z^*\|^3). \end{eqnarray} $

假设

(3.12) $ \begin{eqnarray} \|z^{k, m-1}-z^*\|\leq L_1\|z^{k}-z^*\|^m=O(\|z^{k}-z^*\|^m) \end{eqnarray} $

成立. 接下来证明

(3.13) $ \begin{eqnarray} \|z^{k, m}-z^*\|\leq L_2\|z^{k}-z^*\|^{m+1}=O(\|z^{k}-z^*\|^{m+1}), \end{eqnarray} $

其中$ L_1, L_2 $为不依赖于$ k $的常数.

由$ (2.12) $式并结合引理$ 3.3 $及$ (3.12) $式, 可得

(3.14) $ \begin{eqnarray} \|z^{k, m}-z^*\|&\leq&\big\|(\Psi'(z^k))^{-1}\big\|\Big[\big\|\Psi(z^{k, m-1})-\Psi(z^*)-\Psi'(z^*)(z^{k, m-1}-z^*)\big\|\\ &&+\big\|\big(\Psi'(z^*)-\Psi'(z^k)\big)(z^{k, m-1}-z^*)\big\|\Big] \\ &\leq&\zeta\Big[\frac{1}{2}\beta\big\|z^{k, m-1}-z^*\big\|^2+\beta\|z^k-z^*\|\|z^{k, m-1}-z^*\|\Big]\\ &\leq&\zeta\beta\Big[\frac{1}{2}L_{1}^{2}\|z^{k}-z^*\big\|^{m-1}+L_1\Big]\|z^{k}-z^*\|^{m+1} \\ &\leq&\zeta\beta\Big[\frac{1}{2}L_{1}^{2}\delta_{2}^{m-1}+L_1\Big]\|z^{k}-z^*\|^{m+1} \\ &=&L_2\|z^{k}-z^*\|^{m+1}=O(\|z^{k}-z^*\|^{m+1}), \end{eqnarray} $

其中$ L_2:=\zeta\beta\Big[\frac{1}{2}L_{1}^{2}\delta_{2}^{m-1}+L_1\Big] $为不依赖于$ k $的常数. 因此定理结论成立.

这一节, 给出几个求解弱非线性互补问题的重要例子来验证所提出的算法HONM的有效性. 我们从几个性能, 包括迭代步(记为'IT'), 所消耗的CPU时间(记为'CPU'), 残量(记为'RES') 方面, 比较了基于Fischer的半光滑牛顿法(记为'FBSN')^{[13, 20]}, 基于'cosh' 的光滑牛顿法(记为'CBSN ')^[26], 基于模的矩阵分裂迭代法(记为'MMSA')^[16]. 终止准则为

$ \begin{eqnarray*} {\rm RES}:=\|\Psi(z^{k})\|<10^{-10} \end{eqnarray*} $

或超过最大迭代次数$ k_{\max} $. 数值实验的机器环境是: Intel(R) Core(TM) i7-2670QM, CPU 2.20GHZ, 内存8 GB的Windows 10操作系统. 软件环境为: MATLAB R2017a.

由算法1可知HONM是一般牛顿法的加速版本且具有高阶收敛速率. 但同时, 也注意到当$ m $很大时, 由于内迭代次数增大将导致计算量也增加. 通常取$ m=2 $即可保证算法$ 2.1 $的高阶收敛效率又不增加过多的内迭代数. 具体看下面例子的结果. 在这些例子, 算法HONM的参数$ m=2 $, FBSN的参数$ \rho=0.15 $, $ \beta=0.25 $, $ p=3.0 $ (详细参见文献[20]). 所有测试中, 若$ {\rm RES}<10^{-17} $, 则置其为$ 0 $. 初始点的选取为$ z_{(1)}^{0}=[\varepsilon_0, q^T]^T $, $ z_{(2)}^{0}=(0, 0, \cdots, 0)^T $以及$ z_{(3)}^{0}=(1, 1, \cdots, 1)^T $, 其中$ \varepsilon_0=0.01^{10} $.

例4.1^[23]  考虑问题(1.2), 即弱非线性互补问题的特殊形式, 其中

$ \begin{eqnarray*} A={\rm tridiag}(-1, 4, -1), q=(-1, \cdots, -1)^T. \end{eqnarray*} $

例4.2  该例问题来源于自由边值问题离散. 令$ \Omega=(0, 1)\times(0, 1) $, 函数$ g $满足$ g(0, y)=y(1-y) $, $ g(x, y)=0, $当$ y=1 $, $ y=0 $或$ x=1 $. 考虑下列问题: 即寻找$ t $使得

(4.1) $ \begin{eqnarray} \left\{ \begin{array}{ll} t\geq 0, & {\rm in}\ \Omega, \\ -\triangle t+f(x, y, t)-8(y-0.5)\geq 0, {\quad} & {\rm in}\ \Omega, \\ t(-\triangle t+f(x, y, t)-8(y-0.5))=0, {\quad} & {\rm in}\ \Omega, \\ t=g, & {\rm on}\ \partial\Omega, \end{array} \right. \end{eqnarray} $

其中$ f(x, y, t) $连续可微且在$ \overline{\Omega}\times\{t:t\geq 0\} $, $ \frac{\partial f}{\partial t}\geq 0 $. 离散后可得到非线性互补问题(1.1) 的形式, 其中$ A={\rm spdiags}([-e -e 4*e -e -e], [-(n-1) -1 0 1 n-1], (n-1)^2, (n-1)^2) $为$ (n-1)^2 $阶大小的块对角矩阵, 其中$ e={\rm ones}((n-1)^2, 1) $, $ \Phi(u)=f(x, y, u) $.

在例$ 4.1 $, 针对不同阶数的$ n $比较了算法HONM, FBSN以及CBSN的不用性能. 从表 1–2, 我们发现HONM, FBSN及CBSN都是有效算法. 然而, 从迭代数和CPU时间上看, 算法FBSN以及CBSN的收敛性能都弱于算法HONM. 同时注意到随着规模的增加这种优势保持得较为明显. 例$ 4.2 $的测试结果显示在表 3–6以及图 1–2. 从这些数值性能以及残量变化趋势图可见, 无论在迭代次数还是CPU耗时, 算法MMSA比算法FBSN好. 然而, 本文提出的算法HONM又明显优于算法MMSA, 特别是当维数$ n=225 $增加到$ n=3969 $时, HONM在求解非线性问题时迭代次数不但没有大幅增加, 有些问题求解过程的迭代次数反而还有所减少, 仅仅在CPU时间上有所增加, 并且残量的下降速度也明显快于算法FBSN和MMSA. 这说明了在计算大规模弱非线性互补问题时算法HONM是相当稳定的. 总而言之, 所构造的HONM是具有高阶收敛性能的加速牛顿型算法.

本文提出了一个求解弱非线性互补问题的高效迭代算法HONM, 该算法可看做是对一般牛顿法的加速迭代版本. 在适当的假设下算法HONM被证明至少有$ m+1 $阶收敛速度. 所提出的这个算法非常适用于求解大规模弱非线性互补问题, 如自由边值离散问题. 数值实验部分详细比较了当前几个流行且有效的算法, 从而进一步验证了本文所设计算法的有效性和实用性. 事实上该方法也可用于其他非线性方程组及变分不等式等问题的求解, 特别是当前比较流行的张量互补问题或张量方程组等^{[15, 25]}, 这也是我们未来即将开展的相关研究工作.