宋 巖,王小利,鐘瑞華

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

本文考慮如下來(lái)源于中子運(yùn)輸理論的非對(duì)稱(chēng)代數(shù)Riccati方程：

其中X∈?n×n是未知矩陣,O是零矩陣,已知矩陣Α=Δ-eqT∈?n×n,B=eeT∈?n×n,C=qqT∈?n×n,D =Γ -qeT∈?n×n,且分別是定義在區(qū)間[0,1]上的Gauss-Legendre權(quán)重和節(jié)點(diǎn).

式(1)可能存在多個(gè)解,但只有最小正解(即元素均為正數(shù)的矩陣)才具有物理意義[1].本文關(guān)注于求解式(1)最小正解的數(shù)值算法.Lu[2]在2005年首次將該矩陣方程轉(zhuǎn)化為一種非線性向量方程,并指出只要通過(guò)計(jì)算向量方程的最小正解,即可得到矩陣方程的最小正解.具體來(lái)說(shuō),若X∈?n×n是式(1)的最小正解,則滿(mǎn)足:

于是通過(guò)計(jì)算式(3)的最小正解,可以得到式(1)的最小正解.

關(guān)于式(3)的數(shù)值算法,目前主要有基于不動(dòng)點(diǎn)迭代的研究,如文獻(xiàn)[2-4],及基于Newton 型迭代的研究,如文獻(xiàn)[5-11].由于Newton型迭代具有快速收斂性,故得到了廣泛的關(guān)注.例如,Lu[7]在2005年提出了基于Newton法的一種迭代算法.進(jìn)一步,Lin等[8]在2008年提出了基于一種三階收斂的兩步Newton法.

本文基于文獻(xiàn)[12]所提出的一種具有三階收斂的兩步Newton法,提出了新的計(jì)算式(3)的數(shù)值迭代算法, 并證明了由該算法迭代產(chǎn)生的向量序列的單調(diào)收斂性. 最后, 通過(guò)一些數(shù)值實(shí)驗(yàn)來(lái)比較該算法和Newton法以及文獻(xiàn)[8]中具有三階收斂的兩步Newton法的計(jì)算效能.

1 一種三階收斂的兩步Newton型方法

對(duì)于u,v∈?n,記x=[uT,vT]T,則式(3)等價(jià)于f(x)= 0.給定初始點(diǎn)x0=[uT0,vT0]T∈?2n,考慮如下形式的兩步Newton迭代：

其中xk=[uTk,vTk]T,yk=[uˉTk,vˉTk]T.式(4)的標(biāo)量形式由Frontini和Sormani[12]在2003年將Newton法和文獻(xiàn)[13]結(jié)合而得到.此外,由式(3)可得f的Jacobi矩陣:

基于迭代格式(4),得到如下用于計(jì)算式(3)的兩步Newton 型算法.

算法1 一種計(jì)算式(3)的兩步Newton 型算法選擇初始向量u0,v0 ∈?n.對(duì)于k = 1,2,…,重復(fù)以下迭代過(guò)程,直至收斂:Step 1.求解關(guān)于uˉk,vˉk ∈?n的線性方程組:(In - G2(uk)- ZkH1(uk))vˉk = Zkgk - G2(uk)vk - H2(vk)uk +1 2 (vk + vk °P?uk + e),(In - G1(vk))uˉk = H1(uk)vˉk + gk,其中Ζk = H2(vk)(In - G1(vk))-1,gk=1 2 (uk + uk °Pvk + e)- G1(vk)uk - H1(uk)vk.Step 2.求解關(guān)于uk + 1,vk + 1 ∈?n的線性方程組:(In - G2(uˉk)- Ζˉk(H1(uˉk))vk + 1 = Ζˉkgˉk - G2(uˉk)vk - H2(vˉk)uk + vk °P?uk + e,(In - G1(vˉk))uk + 1 = H1(uˉk)vk + 1 + gˉk,其中Zˉk = H2(vˉk)(In - G1(vˉk))-1,gˉk = uk °Pvk + e - G1(vˉk)uk - H1(uˉk)vk.

在初始向量u0,v0∈?n滿(mǎn)足一定條件時(shí),可得到由上述算法迭代產(chǎn)生的向量序列單調(diào)收斂于式(3)的最小正解.該收斂性分析將在下一節(jié)給出.

2 算法的單調(diào)收斂性分析

對(duì)于任意的矩陣A,B∈?m×n,A°B=(aijbij)m×n表示A和B的Hadamard 內(nèi)積. 若對(duì)i= 1,2,…,m,j=1,2,…,n,有aij≥bij(aij＞bij),則記A≥B(A＞B);若aij＞0(aij≥0),則稱(chēng)A是正矩陣(非負(fù)矩陣).如果一個(gè)向量的所有分量都是正(負(fù)),我們稱(chēng)它為正(負(fù))向量.如果矩陣A的所有非對(duì)角元素都是非正的,則稱(chēng)實(shí)方陣A是Z-矩陣.Z-矩陣A可以表示成A=sI-B,其中s∈?,B≥O,I是n階單位矩陣.若s＞ρ(B),則稱(chēng)Z-矩陣A是非奇異矩陣,其中ρ(B)是B的譜半徑.

引理1[8]A∈?n×n是一個(gè)Z-矩陣,則以下條件等價(jià):

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

2)A-1≥O;

3)若Av＞0 ∈?n,則v＞0.

引理2[8]對(duì)?x,h1,h2∈?2n,有f″(x)(h1,h2)=[hT1LT1h2,hT1LT2h2,…,hT1LTnh2]∈?2n,其中

pi和是矩陣P和?第i列向量.特別地f?(x)=O,故由Taylor公式有

引理3[7]對(duì)任意h,h1,h2∈?n{0},且h1＜h2,則

引理4[7]若ρ(G(x?))≤1,則f′(x?)=I2n-G(x?)是M-矩陣,其中G(x)由式(5)定義.此外對(duì)任意x滿(mǎn)足0 ≤x＜x?,f′(x)是非奇異的M-矩陣.

引理5對(duì)?v1,v2,u1,u2∈?n,記x=[u1T,v1T]T,y=[u2T,v2T]T,若v1≤v2,u1≤u2,則f′(y)≤f′(x).

證明

故

定理記x?是向量方程f(x)= 0的最小正解,若初始向量x0∈?2n,滿(mǎn)足0≤x0＜x?,且f(x0)＜0,則由算法1迭代產(chǎn)生的向量序列{xk}k≥0和{yk}k≥0有定義,且有以下結(jié)論:

a)f(xk)＜0,k≥0;

b)0≤xk＜yk＜xk+1＜x?,k≥0;

證明利用數(shù)學(xué)歸納法證明a)和b).當(dāng)k= 0 時(shí),利用式(4)及引理1-4 易得a)和b)成立.假設(shè)k=l時(shí),a)和b)成立,下面考慮k=l+ 1的情形.由式(4)和引理5有

即xl+1-yl＞yl-xl.進(jìn)而,由引理3得

故由式(6)得

故當(dāng)k=l+ 1時(shí),f(xl+1)＜0,即a)成立.

對(duì)于b),由式(4)和式(6)得

則

得yl+1＜x*,故由引理4知f′(yl+1)是非奇異的M-矩陣,從而由引理1得f′(yl+1)-1＞O.當(dāng)k=l+ 1時(shí),由式(4)得xl+1＜yl+1,xl+1＜xl+2.由式(6)得及引理5得

進(jìn)而由式(4)和式(6)得

注意到f′(yl+1)是非奇異的M-矩陣,故由引理1知xl+2＞yl+1.于是由式(4)和式(6)得

即x*-yl+1＞yl+1-xl+1,則由引理3知

由式(4)和式(6)得

則

故xl+2＜x*.由此可知,當(dāng)k=l+ 1時(shí),b)成立.因此,對(duì)任意k≥0時(shí),a)和b)都成立.

3 數(shù)值實(shí)驗(yàn)

本節(jié)通過(guò)一些數(shù)值實(shí)驗(yàn)來(lái)比較算法1(記為T(mén)SNM)和文獻(xiàn)[7]中的基于Newton 法(記為NM)以及文獻(xiàn)[8]中的三階收斂的兩步Newton 法(記為MNI)的計(jì)算效能.分別選取(α,c)=(10-7,1 - 5×10-7)和(α,c)=(0.3,0.7)進(jìn)行測(cè)試,在每個(gè)測(cè)試中,給出迭代次數(shù)(記為IT)和相對(duì)殘差(記為ERR).

維數(shù)n分別取512,1 024,2 048,4 096,節(jié)點(diǎn)ci和權(quán)重ωi,i= 1,2,…,n,從區(qū)間[0,1]作數(shù)值積分得出.具體步驟,將區(qū)間[0,1]分為n4個(gè)長(zhǎng)度相等的區(qū)間,并對(duì)每個(gè)區(qū)間應(yīng)用Gauss-Legendre 四點(diǎn)求積公式.

我們選取初始點(diǎn)u0=v0=0∈?n,算法終止的條件是

其中eps = 2-52≈2.2204×10-16,范數(shù)取無(wú)窮范數(shù), 所有實(shí)驗(yàn)均在CPU - 2.30GHz（Intel(R)Core(TM)i5-6200），RAM - 4GB環(huán)境下進(jìn)行,MATLAΒ版本為2018b.數(shù)值實(shí)驗(yàn)結(jié)果如下表1-表2所示:

表1 α = 0.3,c = 0.7 時(shí)的數(shù)值模擬Tab.1 The numerical result for α = 0.3,c = 0.7

續(xù)表1

表2 α = 10-7,c = 1 - 5×10-7 時(shí)的數(shù)值模擬Tab.2 The numerical result for α = 10-7,c = 1 - 5×10-7

從以上數(shù)值實(shí)驗(yàn)結(jié)果可以發(fā)現(xiàn)TSNM 法也適用于接近奇異的情況,并且無(wú)論n取512,1 024,2 048還是4 096,TSNM 法比NM 法和MNI法的迭代次數(shù)少,殘差更小.因此,可以認(rèn)為T(mén)SNM 法在求式(3)的最小正解的問(wèn)題上更優(yōu)于NM法和MNI法.

圖1-圖2 是固定n取4 096,(α,c)為(0.3,0.7)和(10-7,1 - 5×10-7)的TSNM 法,NM 法和MNI 法收斂過(guò)程.由圖可知,在非奇異與接近奇異的情形下,TSNM 法的收斂速度都比Newton法和MNI法更快,特別是接近奇異時(shí),優(yōu)勢(shì)更明顯.

圖1 α = 0.3,c = 0.7,n = 4 096Fig.1 α = 0.3,c = 0.7,n = 4 096

圖2 α = 10-7,c = 1 - 5×10-7,n= 4 096Fig.2 α = 10-7,c = 1 - 5×10-7,n = 4 096