王小利，鄧雅清，梁 娟，凌永輝

（閩南師范大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，福建漳州363000）

本文考慮如下形式的非對稱代數(shù)Riccati方程

其中,O為零矩陣,矩陣A,B,C,D分別為

這里e=(1,…,1)T∈?n,q=(q1,…,qn)T∈?n,Δ=diag(δ1,…,δn),而和分別是用Gauss-Legendre法求得的權(quán)重和節(jié)點.

Juang[1]在1995年就研究了矩陣方程(1)正解的存在性，但實際應(yīng)用中最小正解才有物理意義，于是如何近似計算最小正解是研究者所關(guān)注的主要問題，例如Juang 等[2]提出用Gauss-Seidel 迭代法近似方程的解，并證明了其單調(diào)收斂性.特別地，Lu[3]在2005年證明了矩陣方程(1)的解X∈?n×n必須滿足以下形式

其中°表示Hadamard乘積，n階方陣階向量u和v滿足

這里n階矩陣P=(pij)和矩陣中元素分別為若令f:?2n→?2n,則式(2)可改寫成如下等價形式

于是，求解式(1)的最小正解等價于求解向量式(3)的最小正解.

近年來，有大量的研究著重關(guān)注計算式(3)最小正解的高效數(shù)值算法.Lu[3]在2005年首先給出了一種基于不動點迭代的數(shù)值算法，該算法相較于Juang 等[2]提出的基于Gauss-Seidel 迭代算法有更好的計算效能.隨后，Lu[4]進(jìn)一步提出了基于不動點迭代法和Νewton法的數(shù)值算法，Lin等[5]在2006年對Lu[1]的算法進(jìn)行改進(jìn)得到一種新的算法.此外，Bai[6]等在2008年給出了基于非線性分塊Jacobi迭代和非線性分塊Gauss-Seidel迭代的有效算法.同年，Lin[7]等也提出了一種基于兩步Νewton 法的數(shù)值算法，并證明了相應(yīng)的單調(diào)收斂性.隨后，Lin[8]又通過對文獻(xiàn)[6]中方法進(jìn)行改進(jìn)得到一類新的迭代算法.Huang[9]等在2010年提出了一種基于King-Werner法的數(shù)值算法，并證明了非奇異情形的單調(diào)收斂性及奇異情形(α=0,c=1)的收斂性，數(shù)值結(jié)果顯示該方法在接近奇異時有更好的收斂性.在2017年，Miyajima[10]利用解得特殊結(jié)構(gòu)提出了一種快速算法，但此算法只適用于非奇異情況.在2020年，Angelova[11]等首先用krylov 類型的方法得到低維非對稱微分方程，然后用指數(shù)逼近等方法求解.關(guān)于求解方程(1)的一些其他方法，可閱讀文獻(xiàn)[12-13]及其所引用的參考文獻(xiàn).

本文考慮一種新的Νewton 型迭代，該迭代法的標(biāo)量形式首先由McDougall[14]在2014年提出，隨后Potra[15]在2017年將其推廣到Banach 空間，并證明了該迭代法在非線性算子的Frchet導(dǎo)數(shù)滿足Lipschitz條件下的局部與半局部收斂性，且收斂階至少是特別地，該迭代法每一步迭代運算量與Νewton方法大致相同．本文基于該迭代法給出了式(3)的一種改進(jìn)迭代算法，并證明了這種算法的單調(diào)收斂性.最后的數(shù)值實驗結(jié)果顯示，該算法較Νewton法有更好的計算效能，特別在問題接近奇異情況(α=0,c=1)時，優(yōu)勢更加顯著.

1 一種超平方收斂的改進(jìn)Νewton法

給定初始點x0∈?2n,并記z?1=x0,求式(3)的改進(jìn)Νewton法迭代格式為

其中I2n為2n階單位陣，H1(u)=[u°p1，u°p2，…，u°pn],G1(v)=diag(Pv),向量pi和分別是矩陣P和的第i列.為方便起見，將符號進(jìn)行簡化x=[uT,vT]T∈?2n,f(x)=f(u,v)，f ′(x)=f ′(u,v),G(x)=G(u,v),同時把式(3)的最小正解記為x?=[u?T,v?T]T∈?2n.

基于上述迭代格式，算法1給出了計算式(3)的改進(jìn)Νewton法.

算法1.一種計算式(3)的改進(jìn)Νewton法

初始化過程.選擇初始近似值x0=[uT0,vT0]T∈?2n,并記z?1=x0.

迭代過程.對于k=1,2,…,重復(fù)以下過程

使用Νewton 迭代法求解時，每次迭代計算vk+1的運算量為2n3+3n2+n，計算uk+1的運算量為3n2+n，而上述算法中每次迭代求解vˉk的運算量為2n3+13n2+3n,求解uˉk和uk+1的運算量均為2(n+1),求解vk+1的運算量為2n3+9n2+3n,因此，算法1 的總計算量為4n3+22n2+10n+2，與Νewton 法運算量大致相同.下邊給出保證算法收斂的定理，定理的證明將在第2節(jié)給出.

定理1設(shè)x?∈?2n為式(3)的最小正解，若初始點x0∈?2n滿足0 ≤x0＜x?,且f(x0)＜0,則由算法1 迭代產(chǎn)生的向量序列{xk}k≥0收斂到x?，且有

2 收斂性分析

首先給出一些定義.對矩陣A,B∈?m×n,若滿足aij≥bij,i=1,…,m,j=1,…,n,則稱A≥B.此定義也適用于向量.任意實方陣A,若其對角線下方的元素都是非正的，則稱A為Z-矩陣.對于任意的Z-矩陣，A都可以寫成A=sI?Β的形式，其中s∈?,Β≥O,O 為零矩陣，I為單位陣.若s≥ρ(Β)，則A為M-矩陣.下邊給出M-矩陣的一些性質(zhì).

引理1[7]對于任意Z-矩陣A，以下表述等價

1)A是非奇異的M-矩陣.

2)A?1≥0.

3)對于某個v有Av＞0.

對任意x∈?2n,f ″(x)(h1,h2)=[hT1L1h2,hT1L2h2,…,hT1L2nh2]∈?2n,?h1,h2∈?2n,其中

Lu[2]證明了f ″(x)(h,h)與x無關(guān)，且有以下關(guān)系成立

引理2[4]若(5)中定義的G(x)滿足ρ(G(x?))≤1, 即f ′(x?)=I?G(x?)是M-矩陣.則對于任意0 ＜x＜x?,f ′(x)是非奇異的M-矩陣.

定理1 的證明使用數(shù)學(xué)歸納法.由0 ≤x0＜x?,結(jié)合引理2 知f ′(x0)是非奇異的M-矩陣.從而f ′(z0)?1是非奇異的M-矩陣且f ′(z0)?1≥0.由Taylor展開有

由引理1得y0?x?= ?f ′(z?1)?1[f ′(z?1)(y0?x?)]＜0,即y0＜x?,進(jìn)一步有z0＜x?,結(jié)合引理2得f ′(z0)是非奇異的M-矩陣．因此

結(jié)合式(4)和式(9)有

從而x0＜x1.

根據(jù)式(4)和式(9)及Taylor展開有

即x1?x?=f ′(z0)?1[f ′(z0)(x1?x?)]＜0,結(jié)合式(10)得0 ≤x0＜x1＜x?.注意到f(x0)＜0,從而，式(6)在k=0時成立．

假設(shè)式(6)在k=l時成立.則有f ′(zl)?1≥0及0 ≤xl≤zl＜xl+1＜x?,由式(8)知

下證當(dāng)k=l+1時，式(6)也成立.將f(xl+1)在x=zl處作Taylor展開

由式(11)知f(xl+1)＜0,此外，由歸納假設(shè)xl+1＜x?,故下證yl+1＜x?.根據(jù)Taylor展開及式(12)有

由引理1得yl+1?x?=f ′(zl)?1[f ′(zl)(yl+1?x?)]＜0, 即yl+1＜x?.因從 而f ′(zl+1)可逆，且有f ′(zl+1)?1≥0.通過式(4)和(9)有

從而xl+1＜xl+2.進(jìn)一步，結(jié)合式(12)有

根據(jù)引理1 有xl+2?x?=f ′(zl+1)?1[f ′(zl+1)(xl+2?x?)]＜0,即xl+2＜x?.綜上，當(dāng)k=l+1 時有0 ≤xl+1＜xl+2＜x?成立，從而式(6)得證.

通過式(6)易得序列{xk}k≥0為單調(diào)遞增的非負(fù)向量序列，并且x?為其一個上界，故存在一個x?≥0 使得x?≤x?且xk→x?(k→∞).另一方面，在式(4)中令k→∞，x?同樣是式(3)的一個正解，故有x?≤x?.因此有，xk→x?=x?(k→∞)，定理得證.

3 數(shù)值實驗

本節(jié)使用算法1(記為SQV)和文獻(xiàn)[4]所提出的Νewton法(記為ΝM)來計算文獻(xiàn)[16]中例4.2的最小正解，并比較了兩種算法在不同情形下(α、c的不同取值)的CPU 時間(單位為秒)、迭代次數(shù)IT 和相對殘差ERR，其中CPU時間取算法計算十次的平均值.

其中eps=2?52≈2.220 4×10?16.所有實驗均在CPU-1.99GHz(Intel(R)Core(TM)i7-8565)，RAM-4GB 環(huán)境下進(jìn)行，MATLAB版本為2017b.

在迭代終止條件相同的情況下，選取不同的α，c，n進(jìn)行數(shù)值實驗，分別給出ΝM 法和SQV 法的迭代次數(shù)IT 以及殘差ERR．由表1和表2可以發(fā)現(xiàn)：在非奇異的情況下，n取不同值，相較于ΝM 方法，SQV 法的CPU時間和迭代次數(shù)更少.由表3-4可得，該方法同樣適用于接近奇異的情況，且SQV法的迭代次數(shù)與CPU時間均少于ΝM法．因此認(rèn)為SQV法優(yōu)于ΝM方法.即便是在接近奇異的情況下，也可以求得(3)式的最小正解.

表1 α=0.5,c=0.5 時的數(shù)值模擬Tab.1 Νumerical simulations for α=0.5,c=0.5

表2 α=0.3,c=0.7 時的數(shù)值模擬Tab.2 Νumerical simulations for α=0.3,c=0.7

表3 α=10?7,c=1 ?10?6 時的數(shù)值模擬Tab.3 Νumerical simulations for α=10?7,c=1 ?10?6

表4 α=10?7,c=1 ?5×10?7 時的數(shù)值模擬Tab.4 Νumerical simulations for α=10?7,c=1 ?5×10?7

當(dāng)n=409 6 固定不變時，取不同α,c進(jìn)行數(shù)值實驗得到兩種算法的收斂過程見圖1.a)是當(dāng)α=0.5,c=0.5 時的收斂情況，b)是當(dāng)α=0.7,c=0.3 時的收斂情況，c)是當(dāng)α=10?7,c=1 ?10?6時的收斂情況，d)是當(dāng)α=10?7,c=1 ?5×10?7時的收斂情況.顯然在相同條件下，SQV法的迭代次數(shù)少于ΝM法，收斂效果優(yōu)于ΝM法.特別是當(dāng)接近奇異情況下，SQV法的優(yōu)勢更為明顯.

圖1 不同參數(shù)時的收斂情況Fig.1 The convergence history at different parameters.