張 茹，徐春偉，靳聰明

(浙江理工大學(xué)理學(xué)院，杭州 310018)

?

最大熵方法在計(jì)算二維不變測(cè)度中的應(yīng)用

張 茹，徐春偉，靳聰明

(浙江理工大學(xué)理學(xué)院，杭州 310018)

在最大熵方法中采用三角單元上的分片線性基函數(shù)作為矩函數(shù)，用于二維映射的不變測(cè)度計(jì)算。有限元單元上分片線性基函數(shù)有局部支集和單位分割性質(zhì)，因此最大熵方法得到的非線性方程組的Jacobi矩陣是正定的帶狀矩陣，保證了非線性方程組解的唯一性和穩(wěn)定性。文中從理論上給出該方法的收斂階，并且數(shù)值實(shí)驗(yàn)結(jié)果與理論一致。

有限元單元；不變測(cè)度；分片線性基函數(shù)；最大熵

0 引 言

很多科學(xué)和工程問(wèn)題，都可以歸納為對(duì)離散動(dòng)力系統(tǒng)漸近性質(zhì)的研究。例如在分子動(dòng)力學(xué)模擬[1]中可得到每個(gè)原子任意時(shí)刻的位置和速度，但原子的運(yùn)動(dòng)軌跡和速度是雜亂無(wú)章的，只有通過(guò)統(tǒng)計(jì)平均才能有效預(yù)測(cè)分子的溫度、壓力等熱力學(xué)性質(zhì)和結(jié)構(gòu)性質(zhì)。映射的不變測(cè)度[2]反映了動(dòng)力系統(tǒng)的統(tǒng)計(jì)性質(zhì),因此對(duì)映射不變測(cè)度的研究非常重要。

假設(shè)映射S具有絕對(duì)連續(xù)不變測(cè)度，它的密度函數(shù)是映射S對(duì)應(yīng)的Frobenius-Perron算子的不動(dòng)點(diǎn)。關(guān)于不變密度的計(jì)算，Ulam[3]首先提出了分片常數(shù)逼近法，Ding等[4]提出了分片線性Markov方法，Boyarsky等[5]利用Markov變化將不變密度問(wèn)題等價(jià)為一個(gè)矩陣特征值問(wèn)題。最大熵方法由Jaynes[6]首先提出，已廣泛應(yīng)用于物理、化學(xué)、生物等領(lǐng)域,如一維映射的不變測(cè)度計(jì)算[7-8]、物種分布分析[9]、Fredholm積分方程求解[10-11]、化學(xué)反應(yīng)能量依賴(lài)關(guān)系分析[12]、圖像處理[13-14]、光譜分析[15]等。在傳統(tǒng)的最大熵方法中，用標(biāo)準(zhǔn)單項(xiàng)基底1,x,x2,…,xn作為矩函數(shù)，會(huì)產(chǎn)生病態(tài)非線性方程組。近年來(lái)，對(duì)一維映射的不變測(cè)度計(jì)算，Ding等[8]提出了基于分片線性基函數(shù)的最大熵方法。這種方法克服了產(chǎn)生病態(tài)方程組的缺陷，理論和數(shù)值結(jié)果表明該方法收斂。

高維映射的不變測(cè)度可用于研究分子構(gòu)象變化[1]等問(wèn)題。本文基于有限元方法中的三角元上的分片線性基函數(shù)，用最大熵方法計(jì)算二維映射的不變測(cè)度。分片線性基函數(shù)滿(mǎn)足單位分割性質(zhì)和支集性質(zhì)，因此由最大熵方法得到的非線性方程組的Jacobi矩陣是正定的帶狀矩陣。該方法簡(jiǎn)化了非線性方程組的計(jì)算，保證了解的唯一性和穩(wěn)定性。

1 基于分片線性基函數(shù)的最大熵方法

1.1 不變測(cè)度及最大熵原理

定義1[16]設(shè)(X,∑,μ)是一個(gè)有限測(cè)度空間，可測(cè)映射S:X→X為非奇異變換，即由μ(S-1(B))=0，有μ(B)=0。若對(duì)于所有的B∈∑，有μ(S-1(B))=μ(B)，則μ是S的不變測(cè)度。

定義2[16]令(X,∑,μ)為一個(gè)測(cè)度空間，可測(cè)變換S:X→X是非奇異變換。由公式

定義的算子Ps:L1(X)→L1(X)稱(chēng)為對(duì)應(yīng)于S的Frobenius-Perron(F-P)算子。

定理1[16]令Ps是可測(cè)變換S的F-P算子，并且f∈L1。則對(duì)于絕對(duì)連續(xù)有限測(cè)度

在S下不變，當(dāng)且僅當(dāng)f*是Ps的一個(gè)不動(dòng)點(diǎn)。

(1)

其中，當(dāng)f(x)=0時(shí)，有flnf=0。

對(duì)于密度函數(shù)f(x)，若矩

(2)

已知gi(i=1,…,n)是矩函數(shù)，密度函數(shù)f(x)的解并不唯一。考慮如下約束優(yōu)化問(wèn)題：

i=1,…,n}

(3)

命題1[16]約束優(yōu)化問(wèn)題(3)的解為：

(4)

其中，λi(i=1,…,n)滿(mǎn)足

(5)

1.2 二維映射不變測(cè)度的計(jì)算

其中,

對(duì)任意(x,y)∈X，基函數(shù){φi,j(x,y)}滿(mǎn)足單位分割性質(zhì)，即：

(6)

圖1 X的三角剖分

在最大熵方法中，選擇φi,j(x,y),(i,j=0,1,…,n)為矩函數(shù)，則矩為

mij=?Xf*(x,y)φi,j(x,y)dxdy=?suppφi,j(x,y)f*(x,y)φi,j(x,y)dxdy

(7)

(8)

由式(8)可得命題2。

命題2 若f∈L1(X)是一個(gè)非負(fù)函數(shù)，且滿(mǎn)足

?Xf(x,y)φi,j(x,y)dxdy=mij,i,j=0,1,2,…,n

(9)

則f是概率密度，即?Xf(x,y)dxdy=1。

證明：由式(6)和式(8)知，

由命題2和最大熵方法知，若λi,j(i,j=0,1,…,n)滿(mǎn)足

=mij,i,j=0,1,…,n

(10)

則

(11)

是近似不變測(cè)度。這里原最大熵方法式(5)簡(jiǎn)化為式(10)，基函數(shù)的單位分割性質(zhì)使計(jì)算變得簡(jiǎn)單。

對(duì)于方程(10)的求解，考慮等式左邊，則有

(12)

△k的一個(gè)頂點(diǎn)為(xi,yj)。在本文的三角剖分中，可能存在兩種情況：

=?△kφs(x,y)exp(λi,jφi,j(x,y)+λi+1,jφi+1,j(x,y)+

λi,j+1φi,j+1(x,y))dxdy

(13)

其中s取下標(biāo)(i,j),(i+1,j),(i,j+1)；或者

=?△kφs(x,y)exp(λi,jφi,j(x,y)+

λi-1,jφi-1,j(x,y)+λi,j-1φi,j-1(x,y))dxdy

(14)

其中s取下標(biāo)(i,j),(i-1,j),(i,j-1)，(i,j)，是三角形△k的直角頂點(diǎn)。此積分可用數(shù)值積分方法計(jì)算，得到一個(gè)關(guān)于λi,j(i,j=0,1,2,…,n)的非線性方程組。

紅色節(jié)日其背后是許多影響歷史的關(guān)鍵節(jié)點(diǎn)和歷史事件的原型，它猶如紅色文化中的顆顆珍珠，串起我黨紅色歷史與紅色文化。通過(guò)在這些節(jié)日舉辦的各種紀(jì)念活動(dòng)，可以讓大學(xué)生重溫紅色歷史，加深對(duì)中國(guó)共產(chǎn)黨的理解與情感，因此可以以制度化的形式來(lái)推進(jìn)這些紅色節(jié)日，以促進(jìn)和弘揚(yáng)紅色文化的傳承。

若X∈R2被分割成任意三角形單元，假設(shè)任意三角形單元的頂點(diǎn)分別是(xi,yi),(xj,yj),(xk,yk)，如圖2。采用有限元方法中的坐標(biāo)變換

0≤ξ,η≤1

(15)

則任意三角形單元就可以轉(zhuǎn)換為標(biāo)準(zhǔn)三角形單元(如圖3)，其頂點(diǎn)關(guān)系為：

(xi,yi)?(0,0)，(xj,yj)?(1,0)，

(xk,yk)?(0,1).

圖2 任意三角形單元

圖3 標(biāo)準(zhǔn)三角形單元

考慮標(biāo)準(zhǔn)三角形單元上的積分，頂點(diǎn)(0,0),(0,1),(1,0)對(duì)應(yīng)的分片線性基函數(shù)分別表示為

(16)

那么三角形△k上的積分可以表示為：

(17)

=(xj-xi)(yk-yi)-(xk-xi)(yj-yi)

(18)

利用數(shù)值積分可以得到關(guān)于λi,j(i,j=0,1,2,…,n)的非線性方程組(10)，可以用牛頓法或擬牛頓法迭代求解。

令G:R(n+1)2→R(n+1)2，它的元素定義為：

Gij(λ0,0,λ0,1,…,λn,n)=?Xφi,j(x,y)

(19)

其中，i,j=0,1,2,…,n.則方程組(10)可以寫(xiě)為

Gij(λ0,0,λ0,1,…,λn,n)=0

(20)

G的Jacobi矩陣

φm,l(x,y)dxdy

(21)

其中，i,j,m,l均從0,1,…,n取值。與一維問(wèn)題類(lèi)似，很容易證明該Jacobi矩陣是正定的，由于φi,j(x,y)只在以(xi,yj)為頂點(diǎn)的三角形單元上非零，所以Jacobi矩陣是帶狀矩陣。這些性質(zhì)保證了非線性方程組可以用擬牛頓方法有效求解。

2 收斂性分析

本節(jié)采用矩問(wèn)題的收斂理論[17-18]分析了計(jì)算二維映射不變測(cè)度的最大熵方法的收斂性，并給出誤差分析。

設(shè)X是凸拓?fù)湎蛄靠臻g，令W:L1(X)→[-∞,∞)是X上具有緊水平集的泛函，{Dn}是L1(X)上的閉集序列，且對(duì)于所有的n都有Dn+1?Dn?？紤]最大值問(wèn)題

(22)

和極限問(wèn)題

(23)

設(shè)對(duì)于每個(gè)n，fn是式(22)的解。

引理1[18]如果f*是式(23)的唯一解且滿(mǎn)足W(f*)>-∞。那么在拓?fù)鋁下，函數(shù)序列fn收斂到f*，序列W(fn)收斂到W(f*)。

熵H(f)=-?Xf(x,y)lnf(x,y)dxdy滿(mǎn)足一般理論的所有條件。如果將X剖分成2*2n(n=1,2,…)，如圖1中的嵌套直角三角形即Dn+1?Dn，式(10)的可行集對(duì)于包含關(guān)系是單調(diào)遞減的。根據(jù)引理1有下列性質(zhì)：

性質(zhì)1 若H(f*)>-∞，那么對(duì)于Dn的最大熵的解fn有

a)當(dāng)n→∞時(shí)，fn弱收斂到f*。

b)當(dāng)n→∞時(shí)，H(fn)→H(f*)。

此方法的收斂率依賴(lài)于函數(shù)lnf*到由所有φi張成的子空間的最小距離[19]

(24)

由收斂理論和二維空間中連續(xù)函數(shù)的分片線性插值理論，可以得到下面的誤差估計(jì)。

定理3 假設(shè)X上的f*(x)≥c，c是一個(gè)正數(shù)，f*在[a,b]上是連續(xù)可微的。則

(25)

其中，δ是三角形的最大直徑。

3 數(shù)值實(shí)驗(yàn)結(jié)果

本節(jié)把基于三角單元的分片線性最大熵方法應(yīng)用到具體例子。數(shù)值結(jié)果表明收斂階與理論估計(jì)一致。

例1 考慮二維映射

S(x,y)=(S1(x),S2(y))，

其中：

表1 S(x,y)和T(x,y)不變測(cè)度f(wàn)n的L1-誤差

對(duì)于大多數(shù)映射，其不變測(cè)度是未知的，則假設(shè)映射是遍歷的，由Birkhoff的逐點(diǎn)遍歷定理[20]有：

選擇一個(gè)足夠大的正整數(shù)M，可得矩的近似值為：

例2 考慮二維映射

T(x,y)=(T1(x),T2(y))，

其中:

T2(y)=4y(1-y).

4 結(jié) 論

本文將基于有限元中分片線性三角元基函數(shù)的最大熵方法進(jìn)行推廣,用于計(jì)算二維映射的不變測(cè)度。選用三角單元上的分片線性基函數(shù)作為矩函數(shù)，其滿(mǎn)足單位分割性質(zhì)和支集性質(zhì)，使得非線性方程組的Jacobi矩陣是正定的帶狀矩陣，簡(jiǎn)化了方程組的求解，提高了解的精度，保證了非線性方程組解的穩(wěn)定性。數(shù)值實(shí)驗(yàn)結(jié)果表明，用基于分片線性基函數(shù)的最大熵方法計(jì)算二維映射的不變測(cè)度是有效并且收斂的。按照同樣的思路，可把此方法推廣到高維問(wèn)題，但計(jì)算量會(huì)隨之增大。

[1]FRENKELD,SMITB.Understandingmolecularsimulation:fromalgorithmstoapplications[J] . 2nded.PhysicsToday,1996,50(7):66.

[2]GORAP,BOYARSKYA,PROPPEH.Onthenumberofinvariantmeasuresforhigher-dimensionalchaotictransformations[J].JournalofStatisticalPhysics,1991,62(3):709-728.

[3]UlamSM.Bookreviews:Acollectionofmathematicalproblems[J].Science,1960,132(3428):665-666.

[4]DINGJ,ZHOUAH.ConstructiveapproximationsofMarkovoperators[J].JournalofStatisticalPhysics,2001,105(5):863-878.

[5]BOYARSKYA,GORAP.LawsofChaos[M].Boston:Birkh?userBoston,1997:45-50.

[6]JAYNESET.Informationtheoryandstatisticalmechanics[J].PhysicalReview,1957,106(4):620-630.

[7]PEYMANESLAMI,MICALMISIUREWICZ.Singularlimitsofabsolutelycontinuousinvariantmeasuresforfamiliesoftransitivemaps[J].JournalofDifferenceEquations&Applications,2012,18(4):739-750.

[8]DINGJ,NOAHHRHEE.Birkhoff'sergodictheoremandthepiecewise-constantmaximumentropymethodforFrobenius-Perronoperators[J].InternationalJournalofComputerMathematics,2012,89(8):1-9.

[9]BARNHARTPR,GILLAMEH.Theimpactofsamplingmethodonmaximumentropyspeciesdistributionmodelingforbats[J].ActaChiropterologica,2015,16(8):241-248.

[10]JINCM,DINGJ.SolvingFredholmintegralequationsviaapiecewiselinearmaximumentropymethod[J].JournalofComputationalandAppliedMathematics,2016,304:130-137.

[11]MEADLR.ApproximatesolutionofFredholmintegralequationsbythemaximumentropymethod[J].JournalofMathematicalPhysics,1986,27(12):2903-2907.

[12]GALLISMA,HARVEYJK.Maximumentropyanalysisofchemicalreactionenergydependence[J].JournalofThermophysics&HeatTransfer,2015,10(2):217-223.

[13]SIDDIQIMH,ALAMMGR,HONGCS,etal.AnovelmaximumentropyMarkovmodelforhumanfacialexpressionrecognition[J].PlosOne,2016,11(9):3-8.

[14]GULLSF,SKILLINGJ.Maximumentropymethodinimageprocessing[J].CommunicationsRadar&SignalProcessingIeeProceedingsF,1984,131(6):646-659.

[15]FERNANDEZJE,SCOTV,GIULIOED.SpectrumunfoldinginX-rayspectrometryusingthemaximumentropymethod[J].RadiationPhysics&Chemistry,2014,95(4):154-157.

[16] 文蘭.動(dòng)力系統(tǒng)簡(jiǎn)介(迎接ICM2002特約文章)[J].數(shù)學(xué)進(jìn)展,2002,31(4):293-294.

[17]TEBOULLEM,VAJDAI.Convergenceofbestφ-entropyestimates[J].IEEETransactionsonInformationTheory,1993,39(1):297-301.

[18]LEWISAS.Theconvergenceofentropy-basedapproximationsformomentproblems[J].Optimization,2010,28(3):383-395.

[19]MACKEY,MICHAELC.Chaos,Fractals,andNoise[M].Berlin:Springer-Verlag,1994:97-100.

[20]BISWASP,SHIMOYAMAH,MEADLR.Lyapunovexponentsandthenaturalinvariantdensitydeterminationofchaoticmaps:aniterativemaximumentropyansatz[J].JournalofPhysicsAMathematical&Theoretical,2009,43(12):835-842.

(責(zé)任編輯： 康 鋒)

Application of Maximum Entropy Method in Calculating Two-Dimensional Invariant Measure

ZHANGRu,XUChunwei,JINCongming

(School of Sciences, Zhejiang Sci-Tech University, Hangzhou 310018, China)

The piecewise linear basis function on the triangle unit is used as the moment function in maximum entropy method to calculate invariant measure of two-dimensional mapping. Since the piecewise linear basis function on the finite element unit has local support set and unit partition property, Jacobi matrix of nonlinear equations obtained from the maximum entropy method is positive band matrix, which guarantees the uniqueness and stability of the solution. In theory, this paper gives convergence degree of this method, and numerical experiment results are consistent with the theory.

finite element unit; invariant measure; piecewise linear basis function; maximum entropy

10.3969/j.issn.1673-3851.2017.07.017

2016-11-14 網(wǎng)絡(luò)出版日期： 2017-03-28

國(guó)家自然科學(xué)基金項(xiàng)目(11571314)；浙江省自然科學(xué)基金項(xiàng)目(LY13A010014)；浙江理工大學(xué)521人才培養(yǎng)計(jì)劃項(xiàng)目

張茹(1990-)，女，河北石家莊人，碩士研究生，主要從事動(dòng)力系統(tǒng)和分子模擬等方面的研究。

靳聰明，E-mail：jincm@lsec.cc.ac.cn

O242.21

A

1673- 3851 (2017) 04- 0569- 06