多極展開廣義極小殘值算法

1、第4章 基于多極展開法的廣義極小殘值算法線性方程組的求解的諸多方法中，傳統(tǒng)的Gauss消去法或是由Gauss消去法派生的其他方法解方程所需的計(jì)算次數(shù)與方程秩N的3次方成比例，而在Arnoldi算法67的基礎(chǔ)上提出的GMRES迭代算法用于求解稠密系數(shù)矩陣線性方程組時(shí)，其計(jì)算次數(shù)與N成比例。GMRES方法因其存儲(chǔ)量、計(jì)算量較少，每步迭代具有最優(yōu)性并且迭代不會(huì)中斷等優(yōu)點(diǎn)，成為較為可靠的Krylov子空間投影法。但是當(dāng)問題規(guī)模較大時(shí)，使用GMRES法求解依然很困難。故將其與多極展開法(FMM)相結(jié)合更新計(jì)算結(jié)構(gòu)，提高求解效率。其要點(diǎn)是求解時(shí)可只用行矩陣求解，而無須生成系數(shù)總矩陣，故可減少運(yùn)算量及計(jì)算所

2、需內(nèi)存，提高計(jì)算速度。4.1 廣義極小殘值(GMRES)算法 Krylov子空間方法設(shè)待求線性代數(shù)方程組為 (4-1)式中，為階非對(duì)稱非奇異稠密實(shí)矩陣，為實(shí)數(shù)組成的向量。設(shè)是維線性空間中的一組線性無關(guān)向量，由其張成的線性子空間記為，。Krylov子空間投影方法的基本思想是在子空間中構(gòu)造滿足(4-2)的，即式(4-1)近似解。 (4-2)式中，是另一子空間中的任一向量。令，其中為一維向量，則式(4-2)可寫為 (4-3)故若令，則近似解可表示為 (4-4)其中，若選取，這樣選取的Krylov子空間方法稱之為Arnoldi算法。近似解的構(gòu)造方法如下：首先，取為任一向量，令，則可將(4-1) 式化為

3、 (4-5)之后，特別取，假設(shè)其中各向量線性無關(guān)，再由其構(gòu)造一組標(biāo)準(zhǔn)正交基，正交化過程的計(jì)算流程為：(1)任取，計(jì)算和(2)迭代 記迭代中形成矩陣為，且有，上述正交化過程也叫Arnoldi過程。下一步，構(gòu)造近似解，其推導(dǎo)過程如下：令，因?yàn)橛忠驗(yàn)橛勺涌臻g的一組標(biāo)準(zhǔn)正交基組成，故有，由第2步迭代可推出3個(gè)重要關(guān)系式 (4-6) (4-7) (4-8)(4-6)式和(4-7)式在算法中已直接給出，(4-8)式可由(4-7)式導(dǎo)出，在 (4-7) 式兩邊乘得 (4-9) (4-10)即 (4-11)兩邊取模 (4-12)即 (4-13)(4-6)式可寫成在(4-7)式中令則，所以 (4-14)推導(dǎo)完畢

4、。最后，求解最小二乘問題，判斷并使求得的近似解滿足精度要求。 GMRES算法但是上述Arnoldi算法在理論上難以分析其收斂性，并且存在中斷問題，所以人們轉(zhuǎn)而探索其它的Krylov子空間投影方法，如果選取，這樣的方法稱之為GMRES算法。此時(shí)，在中極小化就等價(jià)于在中極小化。原問題歸結(jié)為最小二乘問題。GMRES算法的計(jì)算步驟可以歸結(jié)為：(1)初始化 任取，計(jì)算；(2)迭代 For （直到滿足條件）doStep1 正交化Step2 標(biāo)準(zhǔn)化 Step3 更新與 為上Hessenberg矩陣，當(dāng)時(shí)第一列省略，并且有；(3)求解最小二乘問題，得到；(4)構(gòu)造近似解。 GMRES算法的實(shí)用化處理GMRES

5、算法應(yīng)用于求解邊界元方程組目前來說是一種比較理想的算法，但是經(jīng)典的算法存在著兩大不足:第一，就是隨著迭代次數(shù)的增加，由Arnoldi過程生成的正交基也不斷增加，從而造成存儲(chǔ)量和計(jì)算量的不斷增加；第二，就是隨著迭代次數(shù)的增加，由于誤差積累造成正交基正交性的喪失，從而使迭代的收斂速度放緩甚至發(fā)散。因此GMRES實(shí)用化的關(guān)鍵就是解決這兩個(gè)不足。對(duì)于第一個(gè)不足，傳統(tǒng)方法采用重啟技術(shù)，其基本思想就是計(jì)算過程中最多存儲(chǔ)m個(gè)基向量，每隔m次迭代就重啟一次，即用最后一次的余量作為新一輪迭代的初始余量，最后一次的迭代解作為新一輪迭代的初始解，即GMRES(m)算法。但是這種方法不僅要根據(jù)經(jīng)驗(yàn)確定m的值，而且也沒

6、有收斂性保證，經(jīng)驗(yàn)表明該法經(jīng)常出現(xiàn)停滯甚至發(fā)散的現(xiàn)象，為此一些學(xué)者做了大量的工作68。但實(shí)踐表明這些方法應(yīng)用于三維彈性靜力問題的邊界元法，效果依然不夠理想，反而非重啟型GMRES取得了比較好的效果69。而為了能使非重啟型GMRES實(shí)用化，必須盡量減少迭代的次數(shù)，以減少存儲(chǔ)量和計(jì)算量。本文采用的方法就是在經(jīng)典算法的基礎(chǔ)上引入預(yù)條件技術(shù)。對(duì)于非對(duì)稱矩陣，預(yù)條件的目的就是使預(yù)條件后的矩陣的特征值的分布盡量鄰近。目前最常用的預(yù)條件技術(shù)有兩類，一類是左預(yù)條件，另一類是右預(yù)條件，它們的變換式為：左預(yù)條件 右預(yù)條件 上兩式中矩陣稱為預(yù)條件陣，它的選取應(yīng)符合兩條原則，第一要使預(yù)條件處理后的矩陣的特征值分布盡量

7、靠近，第二就是便于求解。在邊界元方程組的求解中，經(jīng)常取為系數(shù)矩陣的對(duì)角元素組成的對(duì)角矩陣，稱為對(duì)角預(yù)條件陣，或者取由對(duì)角線元素和次對(duì)角線元素組成的帶狀矩陣，稱為帶狀預(yù)條件陣。變換后的方程組記為。對(duì)于左預(yù)條件，；對(duì)于右預(yù)條件，。具體算法的實(shí)施就是在計(jì)算過程中用預(yù)條件處理后的方程組代替原方程組，再進(jìn)行一定的整理。對(duì)于第二個(gè)不足，可采用雙精度存儲(chǔ)和運(yùn)算或修正的Gramm-Achmidt正交化過程。為提高所得向量組的正交性，本文采用前者使正交程度大大提高，同時(shí)加快收斂速度。為避免存儲(chǔ)量的增加，可以只引入少量雙精度計(jì)算，對(duì)及相應(yīng)計(jì)算采用雙精度存儲(chǔ)，精度混合運(yùn)算時(shí)先將單精度轉(zhuǎn)化為雙精度，再用雙精度計(jì)算。4

8、.2 快速多極展開法(FMM)4.2.1 FMM基本思想 FMM源于靜電場計(jì)算，用于計(jì)算大量粒子間兩兩相互作用的電場解析?；军c(diǎn)在于將粒子按空間位置分為不同的集合，分屬于距離較近的集合及集合內(nèi)的粒子間的相互作用兩兩直接計(jì)算；當(dāng)集合相距足夠遠(yuǎn)時(shí)，集合之間多個(gè)粒子的相互作用根據(jù)多極展開公式近似計(jì)算?？臻g中多個(gè)粒子對(duì)某一粒子作用疊加的直接計(jì)算式為 (4-15)式中，是一組影響系數(shù)。(4-15)式在數(shù)學(xué)、物理、化學(xué)及生物等不同學(xué)科中表達(dá)形式相近，普遍用于化學(xué)中分子動(dòng)力學(xué)和量子力學(xué)模擬、天文學(xué)中大規(guī)模引力系統(tǒng)的計(jì)算、電子工程中的電容和電感的計(jì)算和不可壓縮液體的動(dòng)力學(xué)分析70-73。可見，在使用上述公式求

9、解電荷間的相互作用問題時(shí)，當(dāng)需要同時(shí)計(jì)算的電荷數(shù)增加時(shí)，計(jì)算量將以階上升。由Greengard等提出的快速多極展開法，能夠基于球諧函數(shù)在空間中的多極展開，采用遞歸算法結(jié)構(gòu)，將多個(gè)粒子對(duì)某處的作用層層近似轉(zhuǎn)化為單一粒子的作用，將計(jì)算量降為階。直接計(jì)算如圖4-1所示，當(dāng)距離足夠遠(yuǎn)時(shí)，使用FMM方法相當(dāng)于把Q區(qū)域內(nèi)的所有粒子對(duì)位于P點(diǎn)的粒子的作用近似歸結(jié)為域內(nèi)某一點(diǎn)對(duì)P點(diǎn)粒子的作用，如圖4-2所示。圖4-1 直接計(jì)算Fig. 4-1 Direct Calculation圖4-2 采用多極展開法計(jì)算Fig. 4-2 Calculation with FMM在多層的遞歸算法中，對(duì)給定的區(qū)域，首先定義包含

10、所有源點(diǎn)的最小立方體(或平面)作為計(jì)算區(qū)域。先將計(jì)算區(qū)域按某種方式層層細(xì)化為越來越小的多個(gè)子區(qū)域。在0層，有一個(gè)立方體包含所有的計(jì)算主域。第層是把第層的每個(gè)立方體分成若干份形成的。一個(gè)立方體是否需要再細(xì)分，是根據(jù)其中包含源點(diǎn)的個(gè)數(shù)是否少于某一設(shè)定值，這個(gè)值是根據(jù)解題規(guī)模的大小而確定。之后再從最后一層開始，將相互作用力的計(jì)算分成兩個(gè)部分，第一部分是與該立方體相鄰的立方體，即有公共的邊界點(diǎn)，它們之間的作用力需要直接計(jì)算。第二部分是與該立方體不相鄰的立方體，即沒有公共的邊界點(diǎn)，它們之間的作用力用FMM計(jì)算，將其上的每一單元格中所有粒子的作用歸結(jié)為某一粒子的作用，之后再進(jìn)行其上一層的FMM近似計(jì)算，如

11、此遞歸下去直到最上層。其遞歸計(jì)算結(jié)構(gòu)如圖4-3所示。圖4-3 多層遞歸FMMFig. 4-3 Multi-level FMM4.2.2 FMM相關(guān)定理 FMM的所有公式，建立在球坐標(biāo)系上。為適用FMM，定義階次的球面調(diào)和函數(shù)：，式中，是連帶勒讓德函數(shù)，可以用羅巨格(Rodrigues)公式定義式中，是次勒讓德多項(xiàng)式。下面給出的一系列FMM相關(guān)定理適于域積分項(xiàng)的計(jì)算。其中，多極展開和局部展開定理，描述了電荷集在遠(yuǎn)處的電場只與電荷集等效的總電荷(球心處)大小有關(guān)，而與其分布形式無關(guān)，這與力學(xué)的圣維南原理相當(dāng)。FMM主要依賴多極展開和局部展開的平移。相關(guān)定理敘述如下：(1)多極展開定理假設(shè)有個(gè)強(qiáng)度為

12、的電荷分別位于球坐標(biāo)分別為的點(diǎn)且這些點(diǎn)均位于球心在原點(diǎn)，半徑為的球內(nèi), 如圖4-4所示。則在任意點(diǎn)，當(dāng)時(shí)，由電荷產(chǎn)生的位勢為 (4-16)式中。此外，對(duì)任意的有圖4-4 多極展開Fig. 4-4 Multipole Expansion(2)局部展開定理 假設(shè)有個(gè)強(qiáng)度為的電荷分別位于點(diǎn)上，其坐標(biāo)分別為，如圖4-5所示。且點(diǎn)位于球心在原點(diǎn)，半徑為的球外。則在任意點(diǎn)，當(dāng)時(shí)，由電荷產(chǎn)生的位勢為：，式中，。此外，對(duì)任意的有 (4-17)圖4-5 局部展開Fig. 4-5 Local Expansion(3)多極展開的平移假設(shè)有個(gè)強(qiáng)度為的電荷位于球心在，半徑為的球內(nèi)。則由這些電荷在任意點(diǎn)產(chǎn)生的位勢為 (4

13、-18)式中，是向量的球坐標(biāo)。由此，對(duì)任意位于球心在原點(diǎn)，半徑為的球的外部的點(diǎn)有式中式中 (4-19)對(duì)任意的有，如圖4-6所示。圖4-6 多極展開平移Fig. 4-6 Transiton of the Multipole Expansion(4)局部展開的平移假設(shè)是中的一對(duì)點(diǎn)，是向量的球坐標(biāo)，為自然數(shù)，為階局部展開的中心。如圖4-7所示。在點(diǎn)處的位勢表示為 那么，在中的任意點(diǎn)有 式中。圖4-7 局部展開的平移Fig. 4-7 Transiton of the Local Expansion(5)多極展開到局部展開的變換假設(shè)個(gè)強(qiáng)度為的電荷位于球心在，半徑為的球內(nèi)，并且有。則有對(duì)應(yīng)的多極展開式(

14、4-18)在球心在原點(diǎn)，半徑為的球內(nèi)收斂。因此，對(duì)在任意的點(diǎn)，由電荷產(chǎn)生的位勢用局部展開式寫成 (4-20)式中，由(4-19)式定義，此外，對(duì)任意的有 (4-21)如圖4-8所示。圖4-8 多極展開到局部展開的變換Fig. 4-8 Transformation from the Multipole Expansion to the Local Expansion4.2.3 FMM的使用條件用FMM級(jí)數(shù)展開式計(jì)算粒子集合之間相互作用時(shí)，必須保證參與運(yùn)算的集合相距足夠遠(yuǎn)，否則將嚴(yán)重影響其計(jì)算精度。如圖4-9所示，假設(shè)有兩組點(diǎn)電荷分別位于點(diǎn)集和，在下面條件下，定義點(diǎn)集的相互獨(dú)立性。.x1r.x3.

15、xm.x2.x0.y1.y2.y3.yny0rr圖 4-9 點(diǎn)集示意圖Fig. 4-9 Diagrammatic Sketch of the Set of Points點(diǎn)集獨(dú)立定義 若存在和實(shí)數(shù)，且 則稱，點(diǎn)集與是相互獨(dú)立的，或者稱點(diǎn)集與相距足夠遠(yuǎn)。用表示中心在，半徑為的球，表示中心在，半徑為的球，根據(jù)三角不等式，若，則對(duì)任意有點(diǎn)處由處電荷產(chǎn)生的為實(shí)可表示為，取截?cái)嚯A數(shù)，截?cái)嗾`差。是表征源粒子和場粒子之間間距的量。設(shè)兩集合中心的距離為，則截?cái)嗾`差為由上式得4.3 基于FMM的GMRES算法隨著對(duì)邊界元等問題研究的進(jìn)一步深入，由于研究對(duì)象結(jié)構(gòu)的復(fù)雜和尺寸的巨大，為獲得準(zhǔn)確的分析結(jié)果，需要用較多的

16、單元數(shù)去離散該結(jié)構(gòu)，所獲得的影響系數(shù)矩陣將為一個(gè)階數(shù)巨大系數(shù)稠密的矩陣，若用Gauss消去法一類的直接解法求解，所需的內(nèi)存容量需以幾十GB計(jì)算，運(yùn)行時(shí)間需以天計(jì)算，這是不現(xiàn)實(shí)的。同Gauss消去法相比，GMRES 方法所需計(jì)算時(shí)間由O() 減少至O()，但所需內(nèi)存容量相同，皆為O()。建立在FMM基礎(chǔ)上的GMRES 方法將計(jì)算所需內(nèi)存減少到了O()，并提高了計(jì)算速度，為這一難題的解決提供了可能。經(jīng)分析可以發(fā)現(xiàn)GMRES算法的運(yùn)算量主要來自兩個(gè)方面：(1)在開始迭代前，首先要初始化形成個(gè)矩陣元素(必需)；(2)在形成正交基的每次迭代過程中，計(jì)算矩陣與向量積需要次運(yùn)算?？梢?，提高GMRES 算法計(jì)

17、算效率的關(guān)鍵是矩陣與向量積運(yùn)算的加速。將FMM引入GMRES的基本思想，即隱式構(gòu)造一個(gè)與原離散積分方程相聯(lián)系的稠密矩陣的稀疏表示，再將這一稀疏表示用于GMRES方法所需的矩陣向量積的快速計(jì)算中，以減少運(yùn)算量及計(jì)算所需內(nèi)存，提高計(jì)算速度。在這里，我們將矩陣與向量的乘積理解為計(jì)算給定強(qiáng)度的奇點(diǎn)分布在空間任意一點(diǎn)所產(chǎn)生的勢。根據(jù)前面對(duì)FMM使用條件介紹，給定奇點(diǎn)對(duì)空間任意一點(diǎn)的影響應(yīng)首先分成近場影響和遠(yuǎn)場影響兩部分，之后遠(yuǎn)場部分影響可采用有效的近似計(jì)算勢的的方法近似計(jì)算。根據(jù)多極展開定理，在足夠遠(yuǎn)場的近似計(jì)算中，該影響只與給定奇點(diǎn)的強(qiáng)度有關(guān)，而與其具體分布無關(guān)。采用用于處理含有奇點(diǎn)的勢流邊值問題的F

18、MM方法，將問題域逐層次剖分，再對(duì)矩陣元素進(jìn)行形式處理，使其轉(zhuǎn)化為適用于多極展開公式的形式后，遞歸運(yùn)用勢函數(shù)在球坐標(biāo)系下的多極展開和局域展開表示，可將矩陣中一部分元素與運(yùn)算集中起來，轉(zhuǎn)化為向量與相對(duì)固定的多極系數(shù)的計(jì)算，而無需再分別生成并存儲(chǔ)這些元素。這一稀疏表示被用于快速近似勢的遠(yuǎn)場影響部分，即核函數(shù)與遠(yuǎn)場邊界變量乘積的邊界積分；近場的影響采用直接計(jì)算的方法得到，兩者相加就得到了考慮所有源點(diǎn)影響的場點(diǎn)的勢。由此，便隱式構(gòu)造一個(gè)與離散邊界積分方程相聯(lián)系的稠密系數(shù)矩陣的稀疏表示。這樣，由于計(jì)算過程中便不再需顯式生成，并存儲(chǔ)矩陣中的全部的影響系數(shù)，多極展開GMRES算法可使占用內(nèi)存的乘積最小化，進(jìn)而可使大型稠密矩陣線性方程組的求解具有較低的存儲(chǔ)需求。在現(xiàn)代計(jì)算機(jī)上，內(nèi)存容量的增長速度低于CPU和硬盤的增長速度。內(nèi)存是計(jì)算機(jī)上最“緊俏”的