濾波對(duì)角化方法的簡(jiǎn)單探討及應(yīng)用

傅霞　劉楠

摘要：介紹了濾波對(duì)角化計(jì)算本征模式的方法，其中一個(gè)重要的過程是求解廣義逆矩陣。通過matlab編程實(shí)現(xiàn)了諧振腔電磁場(chǎng)本征值和本征模式的計(jì)算，并對(duì)比了采用奇異值分解和高斯消元法求解廣義逆矩陣對(duì)本征值計(jì)算結(jié)果的影響。

關(guān)鍵詞：濾波對(duì)角化；本征值；本征模式；奇異值分解；高斯消元法

中圖分類號(hào)：TB

文獻(xiàn)標(biāo)識(shí)碼：A

doi：10.19311/j.cnki.16723198.2017.09.098

0引言

濾波對(duì)角化算法（Filter Diagonalization Method，簡(jiǎn)稱FDM）是在1990年由Neuhauser提出的用于計(jì)算量子力學(xué)本征值的算法，對(duì)求解各種系統(tǒng)本征值得到較好的效果。隨著FDM算法在應(yīng)用中的長(zhǎng)時(shí)間探索與研究，其求解過程簡(jiǎn)化得到重大發(fā)展，并且在事件分析、收斂性、可靠性等方面都有很大改善，F(xiàn)DM算法在物理、化學(xué)、生物科學(xué)等領(lǐng)域得到迅速應(yīng)用。

諧振腔本征值和本征模式問題是電磁領(lǐng)域的一個(gè)重要問題，多數(shù)諧振腔的優(yōu)化都要以此為基礎(chǔ)，獲得本征值和本征模式后也利于進(jìn)行工作模式分析。將FDM方法應(yīng)用于電磁波本征模式的計(jì)算，其基本思想是：首先通過給定激勵(lì)源，采用時(shí)域有限差分（FDTD）求解電磁場(chǎng)，由于腔體本身的諧振而產(chǎn)生進(jìn)行濾波作用，然后利于矩陣對(duì)角化方法處理電磁場(chǎng)時(shí)域結(jié)果，求得特定結(jié)構(gòu)中電磁波的本征值及其對(duì)應(yīng)的本征模式。該方法是一種混合方法，在彌補(bǔ)了解析方法和數(shù)值方法的弱點(diǎn)的同時(shí)發(fā)揮它們各自的長(zhǎng)處。

1濾波

對(duì)于存在本征振蕩模式的系統(tǒng)，其參量具有exp（-iωmt）形式，其本征頻率為ωm，該系統(tǒng)可以用狀態(tài)矢量s描述，則s滿足的方程的一般形式為

ids（t）dt=Ls（t）（1）

其差分離散形式為

s（t+Δt）-s（t）Δt=-iLs（t）（2）

如果給定一個(gè)0時(shí)刻的初始值s（0），通過式（2），可以得到s（Δt），然后可以再由s（Δt）得到s（2Δt），一直計(jì)算下去，在沒有外加干擾的條件下，式（3）經(jīng)過若干時(shí)間步長(zhǎng)Δt的演變，系統(tǒng)中只會(huì)存在本征模式的線性組合，也就是說通過方程（2）實(shí)現(xiàn)了濾波過程，s將演變?yōu)楸菊髂Ｊ降木€性組合。從初始時(shí)刻開始，經(jīng)過一個(gè)特定時(shí)間段T后，s（T）、s（T+Δt）、s（T+2Δt）…都會(huì)是共振模式的線性組合。前面t在0～T之間的推進(jìn)過程即為濾波過程。

一般只關(guān)心某個(gè)頻段的本征模式，可以通過使用具有頻域特性的激勵(lì)來區(qū)分希望的和不希望的模式，即通過激勵(lì)，促使某些頻率范圍內(nèi)的共振模式形成。只要在式（2）右端增加一個(gè)激勵(lì)gf（t）即可

s（t+Δt）-s（t）Δt=-iHs（t）+gf（t）（3）

其中g(shù)=g（x，y，z），為空間位置的函數(shù)，gf（t）將激勵(lì)頻率在f（t）內(nèi)且具有與g（x，y，z）相同形式的本征模式。在初始時(shí)間段內(nèi)，經(jīng)過充分的激勵(lì)后，再采用方程（2）濾波，則經(jīng)過較短的時(shí)間后，將使s成為希望存在的頻率范圍內(nèi)的共振模式的線性組合。

2對(duì)角化

對(duì)角化部分是使用小規(guī)模的線性代數(shù)方法在一定的范圍內(nèi)區(qū)分相近的本征值。對(duì)于滿足等式Lvm=λvm的vm稱為H的本征模式，λ稱為本征值。找出L個(gè)向量s1，…，sL（L必須大于本征模式數(shù)M），sl是所有本征模式的線性組合。

采用FDTD法計(jì)算出的電場(chǎng)E或B構(gòu)造R和S矩陣，再根據(jù)對(duì)角化方法采用Matlab編程實(shí)現(xiàn)上述算法，其中關(guān)鍵的一個(gè)步驟就是求解S的廣義逆矩陣S，其計(jì)算方法有多種，包括采用奇異值分解法（SVD）、迭代法。本文分別SVD和高斯消元法兩種方法實(shí)現(xiàn)廣義逆的求解。在Matlab中廣義逆直接可以采用調(diào)用pinv函數(shù)實(shí)現(xiàn)，其本質(zhì)是SVD算法。

另外，也可以根據(jù)下式求解S：

S=ST（SST）-1（12）

但是由于先求SST的逆矩陣，數(shù)值誤差較大，在Matlab中通常轉(zhuǎn)換為采用以下公式求解

S=ST/（SST）（13）

其中“/”為Matlab的右除，其本質(zhì)是采用高斯消元法。

采用頻率在c（23.6m-1）～c（25.3m-1）之間的激勵(lì)，對(duì)1m×1.01m二維矩形腔體進(jìn)行FDTD求解，實(shí)現(xiàn)濾波過程。

采用的網(wǎng)格剖分為500×500。取P和L分別為100和20，采用SVD和高斯消元法兩種方法計(jì)算了本征值，計(jì)算結(jié)果及誤差如表1所示。可以看出采用高斯消元法，對(duì)角化過程的誤差更小，其計(jì)算廣義逆矩陣的時(shí)間為1.7ms，而采用SVD計(jì)算廣義逆矩陣的時(shí)間為0.9ms，相比濾波過程的計(jì)算時(shí)間（與網(wǎng)格數(shù)量有關(guān)，一般需要幾十分鐘），它們之間計(jì)算時(shí)間的差別可以忽略。同時(shí)可以看到兩種計(jì)算得到的本征值與理論值誤差基本一致，這是由于本例中采用FDTD進(jìn)行濾波過程的數(shù)值誤差更大。

4結(jié)論

本文介紹了采用濾波對(duì)角化計(jì)算本征模式的方法，通過matlab編程實(shí)現(xiàn)了諧振腔電磁場(chǎng)本征值和本征模式的計(jì)算，通過對(duì)比奇異值分解法和高斯消元法，發(fā)現(xiàn)高斯消元法用于對(duì)角化更加準(zhǔn)確。

參考文獻(xiàn)

[1]Neuhauser D.Time‐dependent reactive scattering in the presence of narrow resonances：Avoiding long propagation times[J].Journal of Chemical Physics，1991，95（7）：49274932.

[2]方劍委.基于濾波對(duì)角化方法提高傅立葉變換質(zhì)譜數(shù)據(jù)質(zhì)量[D].長(zhǎng)沙：國(guó)防科學(xué)技術(shù)大學(xué)，2013.

[3]曾亞文.微波器件諧振腔特征模式分析的算法研究[D].北京：電子科技大學(xué)，2008.

[4]Werner G R，Cary J R.Extracting degenerate modes and frequencies from time-domain simulations with filter-diagonalization[J].Journal of Computational Physics，2008，227（10）：52005214.

[5]郭文彬.奇異值分解及其在廣義逆理論中的應(yīng)用[D].上海：華東師范大學(xué)，2004.

[6]施吉林，史曉發(fā).廣義逆矩陣的兩個(gè)迭代算法[J].大連理工大學(xué)學(xué)報(bào)，1986，（S1）：3745.