三項遞歸公式的數值穩定性分析

郭　浩，范秀香，付波，趙熙臨，趙遠陽，徐光輝

(湖北工業大學電氣與電子工程學院， 湖北 武漢 430068)

三項遞歸公式的數值穩定性分析

郭浩，范秀香，付波，趙熙臨，趙遠陽，徐光輝

(湖北工業大學電氣與電子工程學院， 湖北 武漢 430068)

[摘要]三項遞歸公式是計算正交多項式的高效方法，有些遞歸計算穩定，但有些遞歸計算具有數值不穩定性，可能影響正交多項式的實際應用效果。為判定三項遞歸計算的數值穩定性，將三項遞歸公式按階數轉化為離散狀態方程，基于離散控制理論，利用李雅普諾夫穩定性原理分析三項遞歸計算的數值穩定性，提出了一個判定三項遞歸計算數值穩定的充分條件，并通過實例驗證。

[關鍵詞]三項遞歸公式; 李雅普諾夫第二法; 離散時變線性系統; 穩定性

圖像正交矩面臨高階矩數值不穩定的問題。研究[1]表明，即使沒有基函數離散誤差，高階Chebichef矩仍會因三項遞歸計算的數值傳遞誤差過大而偏離真實值；究其原因，由于計算機存儲器長度限制，產生數值截斷誤差并且在迭代計算過程中逐步積累放大，導致Chebichef多項式在高階出現較大偏差，進而限制其對應的正交矩在分析大圖像時不能達預期效果。因此，R.Mukundan與Liang[2]、Zhang[3]等分別提出針對Tchebichef矩、Hahn矩以及Krawtchouk矩的對稱算法抑制數值誤差的傳遞，這些方法可以一定程度上降低數值誤差的積累速度，但未能從根本上找到決定正交多項式三遞歸計算數值穩定性的本質原因。

這種由于存儲器字長引起的截斷誤差不能消除，只能抑制，合理的數學模型與高效的算法流程可有效降低數值誤差的累積速度。G.R. Amayeh等[4]采用Arbitrary Precision Arithmetic(絕對精度)來提高Zernike多項式的計算精度。尤瑞霖[5]提出加權Krawtchouk多項式關于對稱線的對稱性能，減少了遞歸算法帶來的累積誤差。W.Gautschi[6-7]討論過離散正交多項式的遞歸收斂性，但該解釋的理論性稍弱。由于沒有一套嚴謹的數值誤差穩定性評價體系，目前的正交矩計算方法評價都是以少數試驗驗證，在什么條件下收斂不能確認，限制了正交矩在大型圖像模式識別與重構中的應用。

依筆者所掌握的文獻，目前還未發現有效評價三項遞歸算法數值穩定性的理論與方法。本文立足于尋找判定三項遞歸計算數值穩定性的內在規律，把三項遞歸公式轉變為一個離散時變線性控制系統，利用李雅普諾夫穩定性判別理論分析該系統的穩定性，給出了系統穩定的充分條件，從而確定三項遞歸計算的數值穩定性。

1三項遞歸公式

1.1三項遞歸公式

三項遞歸公式是一類實際使用較多的公式，尤其每一類經典正交多項式都存在一個三項遞歸公式，一般表示為：

Pk(x,N)=Ak(x,N)Pk-1(x,N)-

Bk(x,N)Pk-2(x,N)

(1)

其中P0(x,N)=f0(x,N), P1(x,N)=f1(x,N).公式(1)中的系數Ak(x,N),Bk(x,N)分別是關于k,x,N的多項式，Ak(x,N) = fA(k, x,N), Bk(x,N) = fB(k, x,N)。

1.2遞歸數值誤差

易得，遞歸誤差滿足：

ek(x)=Ak(x)ek-1(x)-Bkek-2(x)

(2)

可見，三項遞歸計算數值誤差分析的本質可以轉化為研究基于式(2)的二階誤差微分方程的數值穩定性。式(2)本質是一個離散時變線性系統，劉永清[11]等人給出了一個時變離散系統穩定性判斷的充分條件，但對于不具有極限定常系統的時變離散系統，此判據無能為力。本文運用李雅普諾夫第二法來研究誤差三項遞歸公式的數值穩定性。

2穩定性分析

2.1李雅普諾夫第二法

在現代控制理論中，李雅普諾夫第二法是研究穩定性的主要方法，該方法通過構造李雅普諾夫函數，從能量增減的角度直接判斷系統的穩定性。

1)V(x)對所有的x都具有連續的一階偏導數。

2)V(x)是正定的。即當x=0時，V(x)=0，當x≠0時，V(x)>0。

2.2離散時變系統穩定性判據

定理：三項遞歸公式(1)數值穩定的一個充分條件是：

存在有界實數序列：{P|P(k)>0},{Q|Q(k)>0},{M|M2(k)≤P(k)Q(k)},使得不等式 {S1+|S2|≤ 0; S3+|S2|≤ 0}成立。其中

(3)

證明：由上節可知，三項遞歸公式(1)的數值穩定性等價于判定二階誤差微分方程ek(x)=Ak(x)ek-1(x)-Bkek-2(x)的穩定性，做變量代換xk=ek-2(x)，yk=xk+1，將誤差微分方程(2)改寫為

(4)

構造李雅普諾夫函數

(5)

其中{P|P(k)>0}, {Q|Q(k)>0},{M| M2(k)≤P(k)Q(k)}，使得V(k,x,y)≥0。計算能量函數(5)的差分函數ΔV：

ΔV(k,x,y)=V(k+1,xk+1,

yk+1)-V(k,x,y)=

[P(k+1)+2M(k+1)A(k)+Q(k+1)A2(k)-

M(k)]2xkyk+[Q(k+1)B2(k)-

當S2≥ 0時，因-2xy≤x2+y2得：

ΔV(k,x,y)

當不等式集{S2≥ 0; S1+S2≤0; S3+S2≤0}滿足時，ΔV(k,x,y) ≤0,系統(4)收斂。同理，S2< 0時，因-2xy> -( x2+y2)得：

不等式集{S2< 0; S1+S2<0; S3+S2<0}滿足時，ΔV(k,x,y)<0,系統(4)收斂。

證畢

3實驗結果與分析

本文分別對經典的Legendre三項遞歸公式和不存在極限系統的三項遞歸公式做數值穩定性分析。

3.1Legendre三項遞歸公式

Legendre三項遞歸公式：

(6)

其中，L(0)=1,L(1)=0.5；x∈[-1,1]。

定義有界序列：令P(k)=e-k,M(k)=0,Q(k)=e-k，則有

把P(k),M(k),Q(k)帶入式(3)有：

S1=e-(k+1)+e-(k+1)A2(k)-e-k

S2=e-(k+1)A(k)B(k)

S3=e-(k+1)B2(k)-e-k

則由定理可求得系統(6)的一個收斂區間為x∈[-0.452, 0.452]。本實驗所求得的收斂區間與所取的P(k)，Q(k)，M(k)序列有關，若另定義一組序列，則可求得相應的收斂區間。

本文采用相對誤差E(k,x)表示系統的數值穩定性

其中，T(k,x)是采用Arbitrary Precision Arithmetic計算遞歸公式的精確值，F(k, x)表示由計算機16位雙精度模式計算的觀測值。若相對誤差發散，說明計算機16位雙精度模式執行的遞歸計算不可靠。在求得的收斂區域內隨機取一點x=0.4分析。圖1a表示在x=0.4時，L(k)精確值(Arbitrary Precision Arithmetic) 隨k的變化趨勢；圖1b表示x=0.4位置遞歸計算的相對誤差的變化趨勢。

從圖1中可知，在x=0.4時，遞歸計算值是收斂的，除偶爾有相對誤差達到10-12，其他絕大部分遞歸相對誤差保持在10-15左右，遞歸過程保持收斂。

(a)L(k)趨勢圖

(b)E(k)趨勢圖圖 1　x=0.4時數值分析

3.2無極限系統三項遞歸公式

T(k)=[0.01sin(kx)+0.35]T(k-1)+

[Vsin(k)+0.6]T(k-2)

(7)

易知：A(k,x)=0.05sin(kx)+0.39;B(k)=-[sin(k)+0.6]。

定義有界序列：P(k)=1,M(k)=1,Q(k)=2，則有

把P(k),M(k),Q(k)帶入式(3)有：

S1=2A2(k)+2A(k)-1

S2=B(k)+2A(k)B(k)+1

S3=2B2(k)-1

則由定理可求得系統(7)的收斂點為x=333和x=355。

在序列P(k)=1，Q(k)=1，M(k)=2條件下，取一收斂點做分析。圖2a表示在x=333時，遞歸公式的精確值(Arbitrary Precision Arithmetic)隨k的變化趨勢，圖2b表示遞歸公式在該點迭代過程中相對誤差的變化趨勢。

從圖2中可知，在x=333時，遞歸計算值是收斂的，除偶爾有相對誤差達到10-14，其他絕大部分遞歸相對誤差保持在10-16左右，遞歸過程保持收斂。

(a)T(k)趨勢圖

(b)E(k)趨勢圖圖 2　x=333時數值分析

4結論

本文對三項遞歸公式的數值穩定性做了分析和探究，構建解析三項遞歸公式的數值誤差等效二階離散時變微分方程，以李雅普諾夫第二法作為判斷依據，給出了數值誤差在多項式遞歸過程中收斂的充分條件。

[參考文獻]

[1]Mukundan R, Ramakrishnan K R. Fast computation of Legendre and Zernike moments[J]. Pattern Recognition,1995, 28(9): 1433-1442.

[2]Liang J, Shu H Z, Zhu H Q, et al. The image by discrete orthogonal moments[J]. Journal of Shanghai Jiaotong University, 2006, 40 (6): 796-800.

[3]Zhang, Guojun LuoZhu, FuBo, et al. A symmetry and bi-recurrence algorithm of accurate computing krawtchouk moments[J]. Pattern Recognition Letters, 2010, 31(7): 548-554.

[4]Amayeh G, Erol A, Bebis G,et al. Accurate and efficient computation of high order zernike moments[C]. ISVC 2005: 462-469.

[5]尤瑞霖. 基于Krawtchouk矩的圖像分析[D].南京:南京理工大學,2012.

[6]Gautschi W. Computational aspects of three-term recurrence relations[J]. SIAM Review, 1967, 9(1): 24-82.

[7]Gautschi W. Is the recurrence relation for orthogonal polynomials always stable[J]. BIT Numerical Mathematics, 1993, 33:277-284.

[8]段廣仁. 線性系統理論[M].黑龍江：哈爾濱工業大學出版社,1996.

[9]龔德國,李奇云.一類離散線性時變系統的穩定性[J].系統工程的理論與實踐,1993(1):15-18.

[10] Chen Xiaodan, Chen Shaobai. The feasible control lyapunov function of nonliner affine-control system[C]. Third International Symposisium on Intelligent Information Technology and Security Informatics, 2010:442-445.

[11] 劉永清, 鄔齊斌, 唐功友. 時變離散系統的穩定性(1)[J]. 山東化工學院學報, 1984(2):56-64.

[責任編校： 張巖芳]

Numerical Stability Analysis of the Three Recursion Relations

GUO Hao, FAN Xiuxiang, FU Bo, ZHAO Xilin, ZHAO Yuanyang, XU Guanghui

(SchoolofElectricalandElectronicEngin.,HubeiUniv.ofTech.,Wuhan430068,China)

Abstract：Three recursive relations are effective way to calculate the orthogonal polynomials. Some recursive calculations are stable, but some recursive calculations are not, which may affect the application effect of orthogonal polynomials. To determine the numerical stability of three recursive calculations, we transformed the three recursion formula into the discrete state equation of numerical errors according to the order number. Based on the discrete control theory, we analyzed numerical stability of three recursive calculations by using the lyapunov second method and put forward a set of sufficient conditions for the robust stability of the three recursion relations. An experiment was designed to verify the feasibility of our conditions.

Keywords：three recursion relations; lyapunov second method; discrete time varying linear system; stability

[收稿日期]2015-06-14

[基金項目]國家自然科學基金(61072130，51109088)；武漢市科技攻關計劃項目(2013012401010845)；湖北工業大學科研基金項目(BSQD12107)；廣東省工業攻關項目(2011B010100037)

[作者簡介]郭浩(1990-)， 男， 湖北武漢人，湖北工業大學碩士研究生，研究方向為圖像處理 [通訊作者] 范秀香(1974-)，女，江西臨川人，湖北工業大學實驗師，研究方向為信號處理

[文章編號]1003-4684(2016)02-0058-04

[中圖分類號]TP13

[文獻標識碼]：A