多層快速多極子技術(shù)MLFMM

1、FEKO Suite 5.5基本算法介紹基本算法介紹快速多極子（快速多極子（FMM）多層快速多極子技術(shù)（多層快速多極子技術(shù)（MLFMM）FMA overview我們知道解線性方程組的方法可分為兩類：一類是直接法，如高斯消元法等；一類是迭代法，如：共軛梯度法等用矩量法（MoM）求解線性方程組，它的系數(shù)矩陣是滿秩的。如用直接法求解，則計算機內(nèi)存需要O（N2），運算量達O（N3）；如用迭代法求解，內(nèi)存一樣需要O（N2），而每次迭代的運算量達O(N2).如此之多的內(nèi)存需要量，如此之大的運算量，大大限制了矩量法的應用范圍，在90年代以前，矩量法僅僅適用于電小尺寸物體（物理尺寸/工作波長 10）。20世紀

2、90年代以后，情況發(fā)生了改變，目前矩量法已可以計算相當大的電大尺寸物體，這主要歸歸功于Rokhlin提出的快速多極子算法，這是一種減少內(nèi)存需減少內(nèi)存需求求，加快矩陣和矢量相乘的技術(shù)?？焖俣鄻O子技術(shù)的基本思路快速多極子技術(shù)的基本思路大家都知道，迭代法的運算量主要取決于矩陣與矢量相乘的運算量。從物理意義上看，矩陣與矢量相乘雖實際上是源點對場點的作用，然由于整個考慮的源點和場點是重合的，因此也可以形式地認為是等效電流之間的相互作用?？焖俣鄻O子技術(shù)的基本思路是首先將未知等效電流分成小組。分小組可以按如下方式進行：首先用一個適當大小的長方體將物體剛好包住，然后將此長方體分成小長方體（小長方體究竟多大合適

3、，下面再作具體分析），并將非空小長方體標出儲存。此處非空小長方體是指其內(nèi)有未知等效電流的小長方體，也就是被物體邊界相割的小長方體。對任何一個非空小長方體，其他的非空小長方體可以分成兩類：一類為近相互作用，一類是遠相互作用。通常，兩小長方體中心之間的距離小于半個波長的為近相互作用，否則為遠相互作用。近遠相互作用介紹近遠相互作用介紹下邊來分析兩小長方體A和B的遠相互作用。設(shè)A和B內(nèi)分別都有100個未知數(shù)，如圖1所示。如果用通常方式來執(zhí)行他們之間的相互作用，則需要100*100次計算機操作。而快速多極子技術(shù)是用一種新的方式來執(zhí)行A和B之間的遠相互作用。其基本思路是將整個相互作用過程分解成三步：聚集、

4、轉(zhuǎn)移、發(fā)散。聚集就是將分布在A內(nèi)的100個未知數(shù)所對應的等效電流聚集在A的中心。其目的是獲得一組具有下列轉(zhuǎn)移特性的新函數(shù)：A內(nèi)所有等效電流對遠處的作用可以由執(zhí)行這組函數(shù)的轉(zhuǎn)移完成；轉(zhuǎn)移就是將聚集過程中得到的一組函數(shù)由A的中心轉(zhuǎn)移到B的中心；發(fā)散就是將轉(zhuǎn)移到B中心的那組函數(shù)發(fā)散到B內(nèi)所有100個未知數(shù)所對應的等效電流上，從而完成A和B的遠相互作用。此種作用方式由圖2表示。下邊會闡述平面波函數(shù)具有上述轉(zhuǎn)移特性，而且在能夠保證高精度情況下，所需平面波個數(shù)少于原未知數(shù)個數(shù)。這就是說，完成新函數(shù)從A中心到B中心的轉(zhuǎn)移，只需要少于100次的計算機操作。這就是快速多極子技術(shù)能夠加快完成A和B遠相互作用的原因

5、。作用過程的分解來源于積分方程中格林函數(shù)的多極子展開，故此項技術(shù)稱為快速多極子技術(shù)。由于格林函數(shù)的多極子展開在近相互作用時很難達到滿意精度，則這種新作用方式只適用于遠相互作用。這也就是我們將相互作用分成近相互作用和遠相互作用的原因。遠相互作用示意圖遠相互作用示意圖圖1：遠相互作用常規(guī)實施辦法示意圖100 個未知數(shù)100 個未知數(shù)100 x 100次計算機操作遠相互作用遠相互作用FMM實施方法實施方法圖2：遠相互作用FMA實施辦法示意圖小于100個平面波小于100個平面波聚集過程發(fā)散過程轉(zhuǎn)移過程100次計算機操作100 x100次計算機操作100 x100次計算機操作OO小長方體多大合適小長方體

6、多大合適這里我們討論小長方體多大合適的問題？由上述分析可以知道：轉(zhuǎn)移步驟所需計算量很小。然而，聚集和發(fā)散步驟并非如此。實際上，將原來的100個未知數(shù)所對應的100個基函數(shù)聚集成大致100個平面波函數(shù)，需要100 x100次計算機操作；將100個平面波函數(shù)發(fā)散給100個未知數(shù)所對應的100個基函數(shù)，也需要100 x100次計算機操作。因此，如果只考慮兩個小長方體的遠相互作用，快速多極子技術(shù)所需計算量是超過通常方式的計算量。然而，當我們考慮100個小長方體的遠相互作用時，情況就不同了。此時有100 x100=10000個未知數(shù)，用通常方式完成它們的相互作用需10000 x10000次計算機操作。如

7、用快速多極子技術(shù)，每個小長方體中的聚集需100 x100次計算機操作，現(xiàn)有100個小長方體，因此整個聚集需100 x100 x100次計算機操作。同樣道理，整個發(fā)散需100 x100 x100次計算機操作。至于轉(zhuǎn)移步驟，因為完成一次轉(zhuǎn)移需100次計算機操作，現(xiàn)有100個小長方體，需100 x100次轉(zhuǎn)移，因此整個轉(zhuǎn)移也需100 x100 x100次計算機操作，故用快速多極子技術(shù)完成整個相互作用需3*100 x100 x100次計算機操作，大大少于通常方式的計算量。由此可見，小長方體尺寸不能過大，因為過大會導致聚集和發(fā)散兩步驟地計算量過大；小長方體尺寸不能過小，因為過小會導致轉(zhuǎn)移步驟的計算量過大

8、。嚴格來說，假如有N個未知量，分成M組，這樣每組大致有N/M個未知數(shù)。根據(jù)上面分析，聚集和發(fā)散步驟所需計算機操作都是O（N2/M），轉(zhuǎn)移步驟需O（NM）次計算機操作。因此完成整個相互作用需要O（N2/M+MN）次計算機操作。不難知道，在M=N1/2時，完成整個相互作用需要的計算機操作次數(shù)最少為O（N3/2）。不難分析，此時的內(nèi)存需要量也為O（N3/2）快速多極子技術(shù)的數(shù)學原理快速多極子技術(shù)的數(shù)學原理離散積分方程系數(shù)矩陣的元素可表示成 （2.71）這里P（G）表示作用在格林函數(shù)G上的算子。假設(shè)x和y分別代表相距較遠的兩個小長方體A、B中的的未知數(shù)。那么x對y的作用可表示成 y=Zx （2.72）

9、快速多極子技術(shù)將此矩陣和矢量相乘分解成聚集、轉(zhuǎn)移、發(fā)散三步驟進行。下邊具體介紹此分解過程。很簡單，主要靠下面兩個數(shù)學恒等式。第一個便是關(guān)于格林函數(shù)的加法定律 （2.73）這里jl是第一類球面Bessel函數(shù)，hl(2)是第二類球面Hankel函數(shù)，Pl是Legendre多項式，以及dr.值得注意的是，在lz時，函數(shù)jl(z)和hl(z)幅度大致保持常數(shù)；在lz時，jl(z)衰減非?？欤鴋l2(z)遞增非?？?。這樣當dr時，式（2.73）能在保證精度下截斷。這樣展開（2.73）便可以寫成 （2.74）通常取L=kd+2ln(kd+pi)就能保證較高精度了。第二個恒等式便是式（2.74）中jlP

10、l的平面波展開 （2.75）|(2)0( 1) (21) ()()()|jk r dlllllejklj kd hkr P d rrd ( )ijijZg P G g dS dS|(2)0( 1) (21) ()() ()|jk r dLlllllejklj kd hkr P d rrd 204 ()()()()ljk dllljj k d P d reP k r d k 這里的積分符號“ ”表示是積分在整個單位球面上進行。此積分可用高斯面積分方法進行。具體來說,就是在區(qū)間0，pi上取L點，使得cos( )在區(qū)間-1，1上滿足Gauss-Legendre L點積分公式。這L點的theta值及積

11、分權(quán)因子可以直接調(diào)用可以直接調(diào)用文獻7中的子程序“”得到。對于水平方向phi值得選取，可以在區(qū)間0，2pi上等間隔選取2L個值。于是式（2.75）右邊的積分便可寫成： （2.76） 這里K=2L2，wp是權(quán)因子， , , 表示單位球面上取樣點的球坐標。將式（2.76）代入式（2.74），并將求和次序交換可得 其中 （2.77） （2.78） （2.79） 注意式（2.77）右邊中的 與r無關(guān)，而 與d無關(guān)。這表明式(2.77)已將格林函數(shù)表示的直接相互作用分解成遠距離的轉(zhuǎn)移和近距離的聚集或發(fā)散。為了更簡明的闡述，不失一般性，以P(G)=G為例說明。如圖3所示，取 ，利用（2.77）便可得矩陣元

12、素（2.71）的FMM表達式快速多極子技術(shù)的數(shù)學原理（續(xù)）快速多極子技術(shù)的數(shù)學原理（續(xù)）2d k(2)0(,)() (21)()()LlppllplTkr k rjlhkr P kr()ppVkd00 ,mom orrrdrr(,)pT kr kr (sincos,sinsin,cos)ppppppkppkkk(,)pp|1()(,)|4j r dKppppppejkV kd T kr krrd 快速多極子技術(shù)的數(shù)學原理（續(xù)）快速多極子技術(shù)的數(shù)學原理（續(xù)） （2.80）其中 （2.81） (2.82)注意Dip、Apj是矢量，Tp是標量。于是式（2.72）的矩陣和矢量相乘便可表示成 (2.83

13、)矩陣D和A的元素已由式（2.81），（2.82）給出。矩陣T的元素可寫成Tpq= 。顯然，T是對角陣。表達式(2.83)便是快速多極子技術(shù)將直接相互作用分解成聚集、轉(zhuǎn)移、發(fā)散三步驟的數(shù)學表達式。在這表達式中，Ax表示聚集過程。它將A中基函數(shù)聚集得到的平面波函數(shù)，所得結(jié)果x1表示K個平面波；Tx1表示轉(zhuǎn)移過程。它將聚集得到的平面波從A中心轉(zhuǎn)移到B中心，所得結(jié)果y1表示在B中心的K個平面波；最后Dy1就是發(fā)散過程。它將K個平面波發(fā)散到B中的基函數(shù)，從而得到最終結(jié)果y。modpjkrpjjSAg eSmodpjkrippiSDg eS1KijipppjpZD T A yD TAxpp qT快速多極

14、子技術(shù)的數(shù)學原理（續(xù)）快速多極子技術(shù)的數(shù)學原理（續(xù)）OOmm圖3：快速多極子技術(shù)具體實施示意圖多層快速多極子技術(shù)的基本思路多層快速多極子技術(shù)的基本思路由前邊的討論可知，快速多極子技術(shù)中的組不能太大，因為那樣轉(zhuǎn)移過程雖然能非常有效地計算，但聚集和發(fā)散過程都不能有效進行；組也不能太小，因為那樣聚集和發(fā)散過程能夠雖能有效進行，然轉(zhuǎn)移過程又不能有效計算。為此，我們是通過組的恰當大小來獲得快速多極子技術(shù)的最佳效率。這里我們將介紹一種新的方式來更有效地實現(xiàn)快速多極子技術(shù)。其基本思路就是將未知數(shù)分成不同層次的組，低層組大，高層組小，讓聚集和發(fā)散過程先在最高層進行，后通過移置、插值完成底層中的聚集和發(fā)散，而轉(zhuǎn)

15、移過程只在每層的部分組之間進行。這種實現(xiàn)方式被稱為多層快速多極子技術(shù)。這里，我們以聚集過程為例來具體闡述這一實現(xiàn)方式。如圖4所示，假設(shè)大組中有4m個未知數(shù)，這樣實現(xiàn)聚集需16m2次計算機操作。如果聚集先在小組進行，需4m2次計算機操作。后將所得的四類以小組中心為起點的平面波移到以大組中心為起點，并相加得到m個以大組中心為起點的平面波，這又需4m次計算機操作。接著再將m個以大組中心為起點的平面波插值，得到4m個大組中心平面波，從而完成大組聚集過程。后邊我們會說明此插值過程需64m次計算機操作。因此這種實現(xiàn)方式總共需4m2+4m+64m次計算機操作。在m很大時，明顯少于原來的16m2次計算機操作。

16、多層快速多極子技術(shù)實施辦法示意圖多層快速多極子技術(shù)實施辦法示意圖mmmmmmmm聚集過程發(fā)散過程轉(zhuǎn)移過程圖4 多層快速多極子技術(shù)實施辦法示意圖多層快速多極子技術(shù)的基本思路多層快速多極子技術(shù)的基本思路為了給出一個完整的多層快速多極子技術(shù)實現(xiàn)步驟，不失一般性，下面以求解域是正方形為例來具體說明。如圖5所示，將正方形分成四個小正方形，這是第一層分組，其中“*”號表示小組中心。后再將每一小正方形一分為四，此過程反復進行，直到最小正方形邊長在半個波長為止。為了敘述方便，這里只給出三層分組。圖5(b)和圖5(c)分別給出第二層和第三層的分組情況，其中符號“。”和“.”分別表示第二層和第三層的小組中心。和快

17、速多極子技術(shù)一樣，多層快速多極子技術(shù)也是將矩陣與矢量相乘分解成聚集、轉(zhuǎn)移、發(fā)散三步驟。然聚集和發(fā)散過程是先在最高層進行，后通過移置中心、以及插值完成低層的聚集和發(fā)散。至于轉(zhuǎn)移過程，與快速多極子技術(shù)中的執(zhí)行方式完全一樣，只是多層快速多極子技術(shù)只在同一層的次相鄰中心（后面將解釋）進行。先考慮聚集過程。以第二層的聚集為例。這層的聚集是通過將其高一層即第三層平面波的中心先移置、后插值完成的。假設(shè)在第三層中的每一小正方體有N3個未知數(shù)，這樣便有M=N/N3個小正方體。因為聚集在小正方體中心“.”的平面波個數(shù)通常小于未知數(shù)個數(shù)N3，因此從中心“.”到中心“?！钡囊浦盟栌嬎懔恳∮贜3xM=N。同樣可以分

18、析得到其他層的移置過程所需計算量也要小于N。下邊會闡述每一層中插值過程所需計算量也大致為O O（N N）。這樣執(zhí)行每一層的聚集過程計算量也就在O O（N N）量級。發(fā)散過程相似，計算量也在O O（N N）量級。三層快速多極子技術(shù)示意圖三層快速多極子技術(shù)示意圖圖圖5 5：三層快速多極子技術(shù)示意圖：三層快速多極子技術(shù)示意圖圖圖5 5（a a）第一層）第一層圖圖5 5（b b）第二層）第二層A A三層快速多極子技術(shù)示意圖三層快速多極子技術(shù)示意圖圖圖5 5：三層快速多極子技術(shù)示意圖：三層快速多極子技術(shù)示意圖圖圖5 5（c c）第三層）第三層a aA A多層快速多極子技術(shù)的基本思路多層快速多極子技術(shù)的基

19、本思路再來考慮轉(zhuǎn)移過程。為了清楚說明多層快速多極子技術(shù)轉(zhuǎn)移過程的實施，閑來解釋上述提及的次相鄰中心概念。次相鄰中心次相鄰中心是指在本層不屬于相鄰中心，但在低一層中，它們隸屬的組是相鄰組的那些中心。譬如在圖5（c）所示的第三層中，淺灰色陰影所蓋正方形中心都是中心a的次相鄰中心，其他的都不是。不難想象，任一中心的次相鄰中心個數(shù)不會超過一個固定常數(shù)C，對于此例而言，應不會超過40。前面已經(jīng)說過，轉(zhuǎn)移過程只在次相鄰中心，因為更遠處未知數(shù)的作用已通過低一層的發(fā)散過程得到。譬如，只有淺灰陰影所蓋正方形中心聚集的平面波需轉(zhuǎn)移到中心a，而深灰陰影所蓋正方形對a的作用已通過A處平面波的發(fā)散得到。假設(shè)這層每一小正

20、方形的平面波個數(shù)為N3，那么就有M=N/N3個小正方形。因為聚集在每一小正方形的平面波個數(shù)小于未知數(shù)個數(shù)N3，又轉(zhuǎn)移只在次相鄰中心進行，故在此層中完成轉(zhuǎn)移過程需大致 CxN3xM=CxN。這樣多層快速多極子技術(shù)中的聚集、轉(zhuǎn)移、發(fā)散，在每一層中的計算量都為N。又對于N個未知數(shù)問題，通常可分 lglgN N 層，故整個計算量為O O（N N lglgN N）。不難分析，多層快速多極子技術(shù)所需內(nèi)存不會超過O O（N NlglgN N）。多層快速多極子技術(shù)的數(shù)學原理多層快速多極子技術(shù)的數(shù)學原理與快速多極子技術(shù)比較，多層快速多極子技術(shù)增加了四個數(shù)學操作：聚集過程中兩個：一個是將以高層組中心為起點的平面波

21、移置到以低層組中心為起點的平面波，一個是小數(shù)目平面波插值得到大數(shù)目平面波；發(fā)散過程中兩個：一個是將以低層組中心為起點的平面波移置到以高層組中心為起點的平面波，一個是大數(shù)目平面波反插值得到小數(shù)目平面波。很明顯，發(fā)散中的兩個操作是聚集中兩個的逆操作，因此下面著重講述聚集中兩操作的數(shù)學表達式。平面波中心移置操作比較簡單，只要先將平面波乘上對應于中心間距相位因子，后將相同平面波方向 的平面波相加，便可得到一列以低層組中心為起點的平面波。但此列平面波的個數(shù)少于低層組所需數(shù)目。由于電場或磁場積分方程中的核載遠處表現(xiàn)得相當光滑，故額外所需平面波可通過插值得到。數(shù)學上有很多插值方法，這里只介紹一種較為實用的雙

22、立方插值。假設(shè)在一粗網(wǎng)格上，yc為一列平面波幅值，此列平面波的方向 對應于一個粗網(wǎng)格，如圖6所示。下面要計算出一列平面波幅值yp，此列平面波的方向 對應于一個細網(wǎng)格。譬如，要計算出對應于方向 的平面波幅值ya，如圖6所示。其中ya可表示成(*)： (2.84) * 參考“W H Press, S A Teukolsky,W T Vetterling and B P Flannery. Numerical Recipes in Fortran. New York: CAMBRIDGE University Press,1992”),(jipk),(ji）（aa,411141jjiijiautcy

23、多層快速多極子技術(shù)中插值實施示意圖多層快速多極子技術(shù)中插值實施示意圖123456789101112131415163d3d2d1d1d2d),(aa圖圖6 6：多層快速多極子技術(shù)中插值實施示意圖：多層快速多極子技術(shù)中插值實施示意圖多層快速多極子技術(shù)的數(shù)學原理多層快速多極子技術(shù)的數(shù)學原理其中 （2.85） （2.86）這里 cij=cl(k)(k=4x(i-1)+j),可由下式計算 (2.87)這里z(i)(i=1,2,4)代表著方向?qū)?，7，10，11點的平面波幅值；z(i)(i=5,6,8), z(i)(i=9,10,12),z(i)(i=13,14,16)分別代表著它們對 的導數(shù), 的導數(shù)，以及它們的交叉導數(shù)，并且26/ )(dta26/ )dua161)()(ikiizwkcl多層快速多極子技術(shù)的數(shù)學原理多層快速多極子技術(shù)的數(shù)學原理1 0 -3 2 0 0 0 0 -3 0 9 -6 2 0 -6 4 0 0 0 0 0 0 0 0 3 0 -9 6 -2 0 6 -40 0 0 0 0 0 0 0 0 0 9 -6 0 0 -6 40 0 3 -2 0 0 0 0 0 0 -9 6 0 0 6 -40 0 0 0 1 0 -3