关于连续性Sylvester方程的广义Richardson迭代

On the Generalized Richardson Iteration of the Continuous Sylvester Equation

摘要: 本文研究当系数矩阵A和B是正半定矩阵，且它们至少有一个是正定时求解连续Sylvester方程的迭代解AX+XB=C的广义Richardson迭代。我们首先分析了求解这类Sylvester方程的广义理查森迭代的收敛性，然后推导了它的最小谱半径的上界以及参数ω的最佳值，通过于HSS方法的比较，强调了所提方法的有效性。

Abstract: In this paper, we study the generalized Richardson iteration for solving the continuous Sylvester equationAX+XB=C, where the coefficient matrices A and B are assumed to be positive semidefinite and at least one of them is positive definite. We first analyze the convergence of the generalized Richardson iteration for solving such a class of Sylvester equations, then derive the upper bound of the minimum spectral radius and the best value of the parameter ω, and emphasize the effectiveness of the proposed method by comparing it with the HSS method.

1. 引言

矩阵方程的迭代算法研究经过经典迭代算法阶段到现代迭代算法阶段再到高级迭代算法阶段。由最初的Jacobi迭代法、Gauss-Seidel迭代法、Richardson迭代法到共轭梯度法、GMRES、BiCGStab、预处理子方法再到并行迭代算法、非线性迭代算法。从经典方法到现代高效算法的发展，反映了数值计算领域不断追求更高效、更稳定的数值方法的进程。这些方法不仅在科学计算中有重要应用，也在工程、金融和数据科学等领域发挥着关键作用。

求解Sylvester方程的方法有直接法和迭代法。常见的直接法是Bartels-Stewart方法[1]和Hessenberg-Schur方法[2]，它们通常通过正交相似变换将矩阵A和B转换为上三角或上Hessenberg形式，然后直接求解。这些方法主要适用于小规模问题。但是对于大型问题，直接求解方法计算量过大且数值不稳定，而迭代方法则具有数值稳定性且便于并行计算。现在对于大型稀疏Sylvester方程，一般采用迭代法求解。常见的方法为古典迭代法(如Smith迭代法[3]、ADI迭代等)和投影迭代法(如Krylov子空间法等)。另外，还有一类方法是利用Kronecker积将矩阵方程转化为大型线性方程组再求解，但这类方法计算量较大，一般是避免使用。

1845年，普鲁士数学家Jacobi提出了将矩阵方程Ax = b中的系数矩阵A分解为D + L + U的方法，其中D为对角阵，L和U分别为严格下三角矩阵和严格上三角矩阵。这一分解形式衍生出了Jacobi迭代法。

$D{x}_{k+1}=-\left(L+U\right)x_k+b$ .

Gauss-Seidel迭代法同样对系数矩阵A进行了类似的分解，其迭代格式为：

$\left(D+L\right)x_{k+1}=-U{x}_k+b$ .

而Richardson迭代法并未对系数矩阵A做分解，其迭代格式为：

$x_{k+1}=\left(I-A\right)x_k+b=x_k+r_k$ .

随后，Richardson迭代法的收敛性能，人们主要通过引入松弛因子和采用预处理方法两种途径来实现。例如：文献[4] [5]引入了松弛因子$\omega$ ，通过采用加权平均构建了逐次超松弛(SOR)迭代算法，显著提高了算法的收敛速度。文献[6]介绍了在系数矩阵为非奇异的M-矩阵或奇异的三对角M-矩阵的情况下，通过适当处理系数矩阵，再利用Gauss-Seidel迭代算法或Jacobi迭代算法，进一步加快了收敛速度。

在2011年，Bai在文献[7]中提出了求解连续Sylvester方程的Hermitian和反Hermitian分裂迭代法(HSS)。在2015年，基于Hermitian和反Hermitian分裂(HSS)迭代法，Li和Wu在文献[8]中提出了求解非Hermitian正定线性系统(Ax = b)的单步HSS迭代法(SHSS)。在2016年，Zeng等人在文献[9]中提出了一种求解复对称线性方程组(Ax = b)的参数化SHSS迭代方法。在近几年也有提出松弛版本的ADI迭代法(RADI)并分析其收敛性；构建非精确的ADI迭代法[10] (2020)年(IADI)并分析其收敛性；构建非精确的RADI迭代法(IRADI)并分析其收敛性等。

本文考虑连续Sylvester方程：

$AX+XB=C$ ，(1.1)

其中，$A\in C^{m\times m},B\in C,C\in C^{m\times n}$ 是给定的，要找到一个符合这个方程的矩阵$X\in C^{m\times n}$ 。并且系数矩阵满足以下条件：

1) 矩阵A，B，C都是大型稀疏矩阵；

2) 矩阵A，B中至少有一个为非Hermitian矩阵；

3) 矩阵A，B都是半正定，且其中至少有一个是正定矩阵；

众所周知，当A和−B不具有共同特征值的时方程(1.1)对于X具有唯一解。

在理论上，方程(1.1)等价于线性方程组

$Ax=c$ ，(1.2)

这里的$A=I\otimes A+B^{\text{T}}\otimes I$ ，其中$B^{\text{T}}$ 表示矩阵B的转置，x和c分别表示矩阵X和C的列向量，$\otimes$ 表示Kronecker积。

受到前面迭代方法的启发，本章将提出广义Richardson迭代来求解矩阵方程(1.1)。

2. Richardson迭代

Richardson迭代法是一种数值计算方法，用于计算连续函数的导数或求解微分方程，该方法的基本思想是通过迭代，以逐步逼近目标函数的导数或微分方程的解。

这里考虑线性方程组$Ax=b$ ，这里$A\in R^{n\times n}$ ，且A为对称正定矩阵，设矩阵A的特征值为$0<{\lambda }_{\min }={\lambda }_1\le {\lambda }_2\le {\lambda }_3\le \cdots \le {\lambda }_{n-1}\le {\lambda }_n\le {\lambda }_{\max }$ ，对任意的$\alpha >0$ ，方程$Ax=b$ 等价于$x=\left(I-\alpha A\right)x+\alpha b$ ，构造迭代：

$x_{k+1}=\left(I-\alpha A\right)x_k+\alpha b$ ，$k=0,1,2,\cdots$ (2.1)

上述迭代为线性方程组$Ax=b$ 的Richardson迭代。考虑迭代格式的收敛性，等价于考虑迭代过程中误差的收敛性。

定理1

设方程组有精确解x，则：

$x=\left(I-\alpha A\right)x+\alpha b$ (2.2)

由(2.1)~(2.2)得：

$e_{i+1}=\left(I-\alpha A\right)e_i, i=0,1,\cdots$

则有：

$e_i=\left(I-\alpha A\right)e_{i-1}={\left(I-\alpha A\right)}^i{e}_0$ ，

假设A为实对称矩阵，则$B=I-\alpha A$ 为实对称矩阵，其特征值${\lambda }_i$ 为实数，特征向量$u_i$ 两两正交，将误差在B特征向量构成的基上展开，有：

$e_i={\displaystyle \sum }_{j=1}^n{c}_j^i{u}_i={\left(I-\alpha A\right)}^i{e}_0={\left(I-\alpha A\right)}^i{\displaystyle \sum }_{j=i}^n{c}_j^0{u}_i={\displaystyle \sum }_{j=1}^n{c}_j^0{\lambda }_j^i{u}_j,$

随着迭代过程的进行，误差逐渐收敛，应有：

$-1<{\lambda }_j<1,i=1,2,\cdots ,n,$

即，

$\rho \left(I-\alpha A\right)<1\Rightarrow \{\begin{array}{l}1-\alpha {\lambda }_{\min }<1\Rightarrow {\lambda }_{\min }>0\\ 1-\alpha {\lambda }_{\max }>-1\Rightarrow 0<\alpha <\frac{2}{{\lambda }_{\max }}\end{array}$ ，

上式给出Richardson迭代的收敛条件：$\rho \left(I-\alpha A\right)<1$ 。

3. 广义Richardson迭代

关于$AX+XB=C$ 的广义Richardson迭代，先将迭代格式化为：

$ver\left(c\right)=\left(I\otimes A+B^{\text{T}}\otimes I\right)ver\left(x\right)$ ，(3.1)

其中的$\omega >0$ ，这个迭代格式可写为：

$X^{\left(k+1\right)}=X^{\left(k\right)}+\omega \left(C-A{X}^{\left(k\right)}-X^{\left(k\right)}B\right)$ ，(3.2)

这就是广义Richardson迭代，其中$\omega$ 是松弛参数，用$\omega$ 取值控制Richardson迭代算法的收敛性和收敛速度。

4. 广义Richardson迭代的收敛性分析

设迭代矩阵为$G_{\omega }=I-\omega U$ ，其中U为矩阵A和B特征值和，其收敛半径为$\rho \left(I-\omega U\right)$ ，假设所有特征值为$u_i,i=1,\cdots ,n$ ，且都是实数，则$u_{\min }\le u_i\le u_{\max }$ 。

然而，所有关于$G_{\omega }$ 的特征值${\lambda }_i$ 满足以下条件：

$1-\omega u_{\max }\le {\lambda }_i\le 1-\omega u_{\min }$ ，

特别地，如果$u_{\min }<0$ 且$u_{\max }>0$ ，则至少有一个特征值大于1，因此对于任意$\omega$ ，$\rho \left(G_{\omega }\right)>1$ 。在这种情况下，方法总是会偏离一些初始猜测。

在给出广义Richardon迭代的收敛定理之前，我们回顾一下以下已知的结果：

引理1：$\rho \left(AB\right)=\rho \left(BA\right)$ ，$\rho \left(A\otimes B\right)=\rho \left(A\right)\cdot \rho \left(B\right)$ ，$\forall A,B\in C^{n\times n}$ 。

引理2：使$K\in C^{n\times n}$ ，${\lambda }_j\in \lambda \left(K\right)$ ，则$\rho \left[{\left(\alpha I+K\right)}^{-1}\left(\beta I-K\right)\right]=\underset{1\le j\le n}{\max }\left|\frac{\beta -{\lambda }_j}{\alpha +{\lambda }_j}\right|$ 。

现在，我们给出广义Richardson迭代的收敛定理。

推论1：

当A、B的特征值为正实数时，即$u_{\min }>0$ ，当$\omega =\frac{2}{u_{\max }+u_{\min }}$ 时，(1.2.1)式收敛，$G_{\omega }$ 的谱半径达到最小为$\rho =\frac{u_{\max }-u_{\min }}{u_{\max }+u_{\min }}$ 。

证明：当特征值都是正的，为了使方法收敛，必须满足以下条件：

$1-\omega u_{\min }<1$ ，

$1-\omega u_{\max }>-1$ ，

第一个条件要求$\omega >0$ ，第二个条件要求$\omega <2/u_{\max }$ 。换句话说，该方法收敛于任何满足的标量$\omega$ 。

$0<\omega <\frac{2}{u_{\max }}$ ，

接下来讨论参数$\omega$ 的最佳$\omega$ 选择值是多少，也就是$\rho \left(G_{\omega }\right)$ 最小的$\omega$ 值是多少，关于$G_{\omega }$ 的谱半径为$\rho \left(G_{\omega }\right)=\max \left\{\left|1-\omega u_{\min }\right|,\left|1-\omega u_{\max }\right|\right\}$ 。

其中$\omega$ 的函数如图1所示。由曲线可知，当斜率为正的曲线$\left|1-\omega u_{\max }\right|$ 与斜率为负的曲线$\left|1-\omega u_{\min }\right|$ 相交时，$\omega$ 值达到最佳。

$-\left(1-\omega u_{\min }\right)=1-\omega u_{\max },$

Figure 1. The curve of $\rho \left(G_{\omega }\right)$ about $\omega$

图1. $\rho \left(G_{\omega }\right)$ 关于$\omega$ 的曲线

为${\omega }_{opt}=\frac{2}{u_{\max }+u_{\min }}$ ，

将其替换为两条曲线中的一条，可以得到相应的最佳谱半径：

${\rho }_{opt}=\frac{u_{\max }-u_{\min }}{u_{\max }+u_{\min }}$ ，

这个表达式显示了小特征值和大特征值存在的困难。对于实际问题，收敛速度可能非常小。此外，为了获得良好的收敛性，需要对特征值进行估计，以获得最优或接近最优的$\omega$ ，这可能会造成困难。最后，由于$u_{\max }$ 可以很大，曲线$\rho \left(G_{\omega }\right)$ 在$\omega$ 的最优值附近可以非常敏感。这些观测结果对于依赖于加速度参数的许多迭代方法都是常见的。

定理2

当特征值都为复数时，假设$u_s={\alpha }_s+{\beta }_{s}i$ ，$s=1\left(1\right)n$ 。是$G_{\omega }$ 的实部特征值，其中${\alpha }_m=\underset{1\le j\le n}{\min }\left\{{\alpha }_j\right\}$ ，${\alpha }_M=\underset{1\le j\le n}{\max }\left\{{\alpha }_j\right\}$ ，${\beta }_M=\underset{1\le j\le n}{\max }\left\{{\beta }_j\right\}$ ，$\Lambda ={\min }_s\left\{2{\alpha }_s/\left({\alpha }_s^2+{\beta }_s^2\right)\right\}$ ，$A={\alpha }_m\left({\alpha }_M-{\alpha }_m\right)$ ，$B=2{\beta }_M^2$ ，${\omega }^{\text{*}}=\frac{2}{{\alpha }_M+{\alpha }_m}$ ，${\omega }_1=\frac{{\alpha }_m}{{\alpha }_m^2+{\beta }_M^2}$ 。

当$0<\omega <\text{Λ}$ 时，广义Richardson迭代的收敛。

当${\omega }_{opt}=\{\begin{array}{ll}{\omega }_1\hfill & A\le B\hfill \\ {\omega }^{\text{*}}\hfill & A\ge B\hfill \end{array}$ ，

(1.2.2)迭代格式的谱半径上界达到最小，如下式：

$\underset{\omega }{\min }\rho \left(G_{\omega }\right)\le \{\begin{array}{ll}\frac{{\beta }_M}{{\left({\beta }_M^2+{\alpha }_m^2\right)}^{1/2}}\hfill & A\le B\hfill \\ \frac{{\left[{\left({\alpha }_M-{\alpha }_m\right)}^2+4{\beta }_M^2\right]}^{1/2}}{{\alpha }_M+{\alpha }_m}\hfill & A\ge B\hfill \end{array}$

证明：

设${\alpha }_m\le {\alpha }_s\le {\alpha }_M<1,s=1\left(1\right)n$ ，则$G_{\omega }$ 的特征值为：

${\lambda }_s=\left(1-\omega {\alpha }_s\right)+\omega {\beta }_{s}i,s=1\left(1\right)n$ ，

其中$|{\lambda }_s|<1$ ，则${\left(1-\omega {\alpha }_s\right)}^2+\omega {\beta }_s^2<1$ ，其中谱半径为$\rho \left(G_{\omega }\right)=\underset{s}{\max }\left|{\lambda }_s\right|=\underset{s}{\text{max}}\left|\left(1-\omega {\alpha }_s\right)+\omega {\beta }_{s}i\right|$ ，

将${\lambda }_s$ 平方可得：${\left|{\lambda }_s\right|}^2={\left(1-\omega {\alpha }_s\right)}^2+\omega {\beta }_s^2={\left|1-\omega {\alpha }_s\right|}^2+{\omega }^2{\left|{\beta }_s\right|}^2$ 。

我们只对$\omega >0$ 的情况感兴趣。设

$\underset{s}{\max }\left|{\beta }_s\right|={\beta }_M\ge 0,$

我们可得：

${\rho }^2\left(G_{\omega }\right)=\underset{s}{\max }\left\{{\left(1-\omega {\alpha }_s\right)}^2+{\omega }^2{\beta }_s^2\right\}\le \underset{s}{\max }\left\{{\left(1-\omega {\alpha }_s\right)}^2+{\omega }^2{\beta }_M^2\right\}$

进一步可得到：

$\underset{s}{\max }\left|1-\omega {\alpha }_s\right|=\{\begin{array}{ll}1-\omega {\alpha }_m\hfill & 0\le \omega \le {\omega }^{\text{*}}\hfill \\ \omega {\alpha }_M-1\hfill & \omega \ge {\omega }^{\text{*}}\hfill \end{array}$ ，

其中${\omega }^{\text{*}}=\frac{2}{{\alpha }_M+{\alpha }_m}$ 。

因此可得如下：

${\rho }^2\left(G_{\omega }\right)\le \{\begin{array}{ll}{\left(1-\omega {\alpha }_m\right)}^2+{\omega }^2{\beta }_M^2=g_1\left(\omega \right)\hfill & 0\le \omega \le {\omega }^{\text{*}}\hfill \\ {\left(\omega {\alpha }_M-1\right)}^2+{\omega }^2{\beta }_M^2=g_2\left(\omega \right)\hfill & \omega \ge {\omega }^{\text{*}}\hfill \end{array}$ ，

因此可得：

$\underset{\omega }{\min }{\rho }^2\left(G_{\omega }\right)\le \{\begin{array}{ll}\underset{\omega }{\min }{g}_1\left(\omega \right)\hfill & 0\le \omega \le {\omega }^{\text{*}}\hfill \\ \underset{\omega }{\min }{g}_2\left(\omega \right)\hfill & \omega \ge {\omega }^{\text{*}}\hfill \end{array}$ ，

为了得到${\min }_{0\le \omega \le {\omega }^{*}}{g}_1\left(\omega \right)$ 和${\min }_{\omega \ge {\omega }^{*}}{g}_2\left(\omega \right)$ 我们可通过对$g_1\left(\omega \right)$ 和$g_2\left(\omega \right)$ 求导，当$\omega \ge {\omega }_1$ 时$g_1\left(\omega \right)\ge 0$ ；当$\omega \ge {\omega }_2$ 时$g_2\left(\omega \right)\ge 0$ ，其中

${\omega }_1=\frac{{\alpha }_m}{{\alpha }_m^2+{\beta }_M^2},$

${\omega }_2=\frac{{\alpha }_M}{{\alpha }_M^2+{\beta }_M^2}.$

1) 当${\omega }_1\le {\omega }^{\text{*}}$ 时，即${\alpha }_m\left({\alpha }_M-{\alpha }_m\right)\le 2{\beta }_M^2$ 时，则：

${\min }_{0\le \omega \le {\omega }^{*}}{g}_1\left(\omega \right)=g_1\left({\omega }_1\right)$ ，

当${\omega }_1\ge {\omega }^{\text{*}}$ 时，即${\alpha }_m\left({\alpha }_M-{\alpha }_m\right)\ge 2{\beta }_M^2$ 时，则：

${\min }_{0\le \omega \le {\omega }^{*}}{g}_1\left(\omega \right)=g_1\left({\omega }^{*}\right)$ .

2) 由于${\omega }_2\le {\omega }^{\text{*}}$ ，故$\underset{\omega \ge {\omega }^{\text{*}}}{\min }{g}_2\left(\omega \right)=g_2\left({\omega }^{\text{*}}\right)=g_1\left({\omega }^{\text{*}}\right)$ 。

可得：

$\underset{\omega }{\min }{\rho }^2\left(G_{\omega }\right)\le \{\begin{array}{ll}{g}_1\left({\omega }_1\right)\hfill & A\le B\hfill \\ g_1\left({\omega }^{\text{*}}\right)\hfill & A\ge B\hfill \end{array}$ ，

$A={\alpha }_m\left({\alpha }_M-{\alpha }_m\right)$ ，$B=2{\beta }_M^2$ .

进一步可得到：

$\underset{\omega }{\min }\rho \left(G_{\omega }\right)\le \{\begin{array}{ll}\frac{{\beta }_M}{{\left({\beta }_M^2+{\alpha }_m^2\right)}^{1/2}}\hfill & A\le B\hfill \\ \frac{{\left[{\left({\alpha }_M-{\alpha }_m\right)}^2+4{\beta }_M^2\right]}^{1/2}}{{\alpha }_M+{\alpha }_m}\hfill & A\ge B\hfill \end{array}$ ，

对于上述$\underset{\omega }{\text{min}}\rho \left(G_{\omega }\right)$ 的上界，我们可以做如下的说明：

1) 如果$A\le B$ ，那么上界总是小于1。而且，等式成立——就是说：

$\underset{\omega }{\text{min}}\rho \left(G_{\omega }\right)=\frac{{\beta }_M}{{\left({\beta }_M^2+{\alpha }_m^2\right)}^{1/2}}$ ，

在${\alpha }_M+i{\beta }_M$ 为G的特征值情况下，有${\omega }_1={\omega }_{opt}$ 。

2) 如果$A\ge B$ ，${\alpha }_m{\alpha }_M>2{\beta }_M^2\ge {\beta }_M^2$ ，可得${\left({\alpha }_M-{\alpha }_m\right)}^2+4{\beta }_M^2<{\left({\alpha }_M+{\alpha }_m\right)}^2$ 然后我们也可以证明上界总是小于1。对此可以证明：

$\frac{{\left[{\left({\alpha }_M-{\alpha }_m\right)}^2+4{\beta }_M^2\right]}^{1/2}}{{\alpha }_M+{\alpha }_m}<1$ .

故在${\alpha }_m+i{\beta }_M$ 和${\alpha }_M+i{\beta }_M$ 为G的特征值情况下，有${\omega }^{\text{*}}={\omega }_{opt}$ 。

5. 数值实验

在本节中，我们比较了HSS和广义Richardson方法的计算性能来说明单步迭代求解Sylvester方程(1)的有效性。

其中由广义Richardson迭代方法的格式可得以下的浮点运算：

$R=C-A\ast X-X\ast B$ $2m^{2}n+2m{n}^2$ ，

$X = X + \omega \ast R$ $2n^2$ ，

得到最后的浮点运算次数是$f_{Ri}\left(m,n\right)=2m^{2}n+2m{n}^2+2n^2$ 。

在例子中的HSS算法是运用的Hermitian和反Hermitian分裂来迭代，其中的相关浮点运算如下：

第一步：将系数矩阵A和B进行如下分裂：

$A=H\left(A\right)+S\left(A\right)m^2,$

$B=H\left(B\right)+S\left(B\right)n^2,$

第二步：

计算交替方向隐式：

$\left(\alpha I+H\left(A\right)\right)X^{k+\frac{1}{2}}+X^{k+\frac{1}{2}}\left(\beta I+H\left(B\right)\right)=\left(\alpha I-S\left(A\right)\right)X^k+X^k\left(\beta I-S\left(B\right)\right)+C,$

$2m^{2}n+4mn+2m{n}^2,$

$\left(\alpha I+S\left(A\right)\right)X^{k+1}+ X^{k+1}\left(\beta I+S\left(B\right)\right)= \left(\alpha I-H\left(A\right)\right)X^{k+\frac{1}{2}}+X^{k+\frac{1}{2}}\left(\beta I-H\left(B\right)\right)+C,$

$2m^{2}n+4mn+2m{n}^2,$

得到最后的浮点运算次数是$f_{HSS}\left(m,n\right)=4m^{2}n+4m{n}^2+m^2+n^2+4mn$ 。

则有如下：

$\frac{f_{HSS}\left(m,n\right)}{f_{Ri}\left(m,n\right)}=\frac{4m^{2}n+4m{n}^2+m^2+n^2+4mn}{2m^{2}n+2m{n}^2+2n^2}$ .

在实际计算中，所有迭代都是从零矩阵开始，在MATLAB中以机器精度10⁻¹⁶进行，当前残差矩阵的范数满足$\Vert \left(R^{\prime },fr{o}^{\prime }\right)\Vert /\Vert C\Vert \le {10}^{-6}$ 时停止。为了方便，在我们的数值实验中，我们分别用IT表示迭代步骤数，用CPU表示计算时间(以秒为单位)，用$\omega$ 表示实验发现的单步迭代迭代参数的局部最优值。

例1

我们考虑二维对流扩散方程$-\left(u_{xx}+u_{yy}\right)+\sigma \left(x,y\right)u_x+\tau \left(x,y\right)u_y=f\left(x,y\right)$ 在单位平方$\left(0,1\right)\times \left(0,1\right)$ 上，具有dirichlet型边界条件，系数$\sigma$ 和$\tau$ 分别表示沿x和y方向的速度分量，这里我们取$\sigma$ 和$\tau$ 为常数的情况。拉普拉斯算子$u_{xx}+u_{yy}$ 为一阶导数$u_x$ 和$u_y$ 的不同的线性系统，自然等效于标准的Sylvester方程，在[11]中，描述了如何使用二阶有限中心差分算子对系数$\sigma$ 和$\tau$ 的任意值获得一般Sylvester方程$AX+XB=C$ ，当系数$\sigma$ 和$\tau$ 为常数时，A和B是三对角Toeplitz矩阵，定义为：

$A=tridiag\left(-1+\frac{\tau h}{2},2,-1-\frac{\tau h}{2}\right)$ 和$B=tridiag\left(-1+\frac{\sigma h}{2},2,-1-\frac{\sigma h}{2}\right)$

其中步长$h=\frac{1}{n+1}$ ，矩阵A对应y方向离散化，矩阵B对应x方向离散化，如果$\tau =\sigma$ ，则A = B。

这里我们考虑Sylvester方程(1.1)，矩阵$A,B\in R^{n\times n}$ 。参数$\tau$ 和$\sigma$ 取决于具有齐次Dirichlet边界条件

的对流扩散问题的性质。函数f定义为$f\left(x,y\right)={\text{e}}^{x+y}$ ，并测试了$\tau$ 、$\sigma$ 和$h=\frac{1}{n+1}$ 的不同值。在此例子中我们的A和B阶数一样，则$\frac{f_{HSS}\left(m,n\right)}{f_{Ri}\left(m,n\right)}=2$ ，HSS的浮点运算大约是广义Richardson的两倍。随着阶数

n的增加，其运算时间和最优的$\omega$ 如下表1所示，用广义Richardson迭代和HSS方法求解了本例中的方程。表1给出了我们在迭代次数(iter)和以秒为单位的CPU时间t(CPU)方面的实验总结，这些实验针对不同的参数$\tau$ 和$\sigma$ 以及矩阵的n阶。

Table 1. Performance of G-Richardson and HSS methods for the Example 1

表1. 广义 Richardson迭代和HSS迭代的比较

G-Richardson					HSS
h = 1/(n + 1)		ω	IT	CPU	α	IT	CPU
τ = 10 σ = 100	0.04	0.138	56	0.0002	0.75	23	0.006
	0.02	0.31	26	0.0003	0.53	30	0.041
	0.01	0.251	109	0.0069	0.31	52	0.419
	0.005	0.249	332	0.0813	0.15	104	4.39
τ = 1 σ = 100	0.04	0.13	53	0.0005	0.625	31	0.007
	0.02	0.248	32	0.0003	0.5	40	0.042
τ = 1 σ = 100	0.01	0.252	106	0.004	0.26	76	0.631
	0.005	0.25	320	0.085	0.19	104	4.044
τ = 50 σ = 0.1	0.04	0.251	85	0.0003	0.45	35	0.008
	0.02	0.249	375	0.005	0.37	40	0.039
	0.01	0.25	1411	0.073	0.17	88	0.685
	0.005	0.249	5068	1.389	0.08	216	8.563

由表1可知，无论当我们的系数$\sigma$ 和$\tau$ 取什么时，我们的广义Richasrdson迭代算法的运行时间远快于HSS算法。

6. 结论

本文给出了广义Richasrdson迭代算法中松弛参数选取的证明与研究，当选取最优的参数，运行时间与此迭代次数都会相对于少很多。当系数矩阵不是对称的时候与HSS算法进行了对比，数值实验表明：无论是收敛速度还是计算复杂度，广义Richardson算法解矩阵方程(1.1)的效果都好于HSS算法。