特殊鞍點(diǎn)問題的高效預(yù)處理方法
【技術(shù)領(lǐng)域】
[0001 ] 本發(fā)明涉及一種特殊鞍點(diǎn)問題的高效預(yù)處理方法。
【背景技術(shù)】
[0002] 大多數(shù)科學(xué)與工程計(jì)算中的模型問題經(jīng)常需要求解一個(gè)或一系列大規(guī)模線性系 統(tǒng)。大規(guī)模線性方程組的求解往往占據(jù)了整個(gè)數(shù)值計(jì)算過程的主要部分，并且已成為科學(xué) 計(jì)算中一個(gè)突出的重要問題，同時(shí)也對實(shí)時(shí)計(jì)算和高精度分析提出了挑戰(zhàn)。因此，如何使 用合理的計(jì)算量求解一個(gè)線性方程組就成了數(shù)值計(jì)算方法的一個(gè)重要課題。當(dāng)今在計(jì)算數(shù) 學(xué)中十分重要的研究課題之一就是高效求解如下大型稀疏線性鞍點(diǎn)系統(tǒng)。
[0003]
[0004] 其中J e R"'C e
，5 e (可能m < < η)。若A對稱正定且C=0,我們稱 上式為經(jīng)典的鞍點(diǎn)問題，否則稱為廣義鞍點(diǎn)問題。
[0005] 鞍點(diǎn)問題在工程和科學(xué)計(jì)算上有著及其廣泛的應(yīng)用，在工程領(lǐng)域，隨著有限元 方法的日益普及，流體力學(xué)和固體力學(xué)成為了鞍點(diǎn)問題的主要源泉之一；在約束最優(yōu)化領(lǐng) 域，內(nèi)點(diǎn)法的每一次迭代都涉及到鞍點(diǎn)問題的求解；線性彈力學(xué)、電磁學(xué)、電路與計(jì)算機(jī)網(wǎng) 絡(luò)等各方面的許多問題也都?xì)w結(jié)為大型稀疏鞍點(diǎn)問題的求解。極其廣泛的應(yīng)用背景使得有 效地求解鞍點(diǎn)問題一直成為工程和科學(xué)計(jì)算的熱點(diǎn)，受到各個(gè)領(lǐng)域眾多學(xué)者的廣泛重視。 因此，鞍點(diǎn)問題的數(shù)值求解是丞待解決的重點(diǎn)課題之一。這不僅因?yàn)樗旧砭哂袠O高的應(yīng) 用背景，而且因?yàn)樗怯?jì)算數(shù)學(xué)與現(xiàn)代科學(xué)計(jì)算領(lǐng)域中極富挑戰(zhàn)性的難題，無論是從理論 方面還是算法本身方面，都存在很多沒有解決的問題。所以，繼續(xù)從事這方面的研究非常 重要，并且非常必要。
[0006] 由于鞍點(diǎn)問題自身特殊結(jié)構(gòu)，使得現(xiàn)有的數(shù)值計(jì)算方法求解這類問題依然存 在計(jì)算量和存儲量大，迭代求解不收斂或收斂很慢。因此，一些經(jīng)典的迭代算法如 Gauss-Seidel, S0R等方法均失效。故需要對鞍點(diǎn)問題進(jìn)行預(yù)處理，即將鞍點(diǎn)問題Ιχ =厶化 為等價(jià)的具有優(yōu)良性質(zhì)的線性系統(tǒng)=/>7)。預(yù)處理Krylov子空間方法是求解這類具 有特殊結(jié)構(gòu)線性方程組的基本方法之一。因此，針對鞍點(diǎn)問題的具體結(jié)構(gòu)和特殊性質(zhì)，設(shè) 計(jì)可行且高效的預(yù)處理子具有重要的理論意義和很高的實(shí)用價(jià)值。
[0007] 國內(nèi)外現(xiàn)狀，鞍點(diǎn)問題廣泛的產(chǎn)生于計(jì)算流體動力學(xué)、限制和加權(quán)最小二乘問題、 限制優(yōu)化、經(jīng)濟(jì)學(xué)、計(jì)算電磁學(xué)、橢圓型偏微分方程的混合有限元離散等應(yīng)用領(lǐng)域中。目前 已存在很多種方法求解鞍點(diǎn)問題，其中包括：直接法、零空間方法、投影法、HSS-型算法和 Uzawa-型算法等。由于來源于實(shí)際問題的鞍點(diǎn)系統(tǒng)一般涉及超大型不定的離散線性方程 組的求解，受限于存儲空間和巨大的計(jì)算費(fèi)用，直接解法難以有效地解決問題，學(xué)者們著 重于探討鞍點(diǎn)問題的迭代求解。近些年來，用于求解鞍點(diǎn)問題的數(shù)值方法相繼取得很大 進(jìn)展，例如，應(yīng)用Schur縮減方法時(shí)，針對Schur補(bǔ)系統(tǒng)求解的巨大計(jì)算費(fèi)用問題，逼近 Schur補(bǔ)的預(yù)處理技術(shù)被研究；將零空間方法與約束預(yù)處理相結(jié)合；預(yù)處理和不精確Uzawa算法被研究，并由線性迭代發(fā)展到非線性迭代；各種穩(wěn)定迭代算法與Krylov子空間方法結(jié) 合的預(yù)處理Krylov迭代法在很大程度上提高了鞍點(diǎn)系統(tǒng)的收斂速度。但由于鞍點(diǎn)系統(tǒng)結(jié) 構(gòu)特殊，且一般為大型稀疏、病態(tài)的線性系統(tǒng)，這種特殊性使得現(xiàn)有的數(shù)值計(jì)算方法求解 這類問題依然存在計(jì)算量和存儲量大，迭代求解不收斂或收斂很慢，通常需要對鞍點(diǎn)問題 進(jìn)行預(yù)處理化為具有優(yōu)良性質(zhì)的等價(jià)線性方程組，但進(jìn)行有效的預(yù)處理是很困難的。
[0008]國內(nèi)外很多學(xué)者做了大量的工作，提出了各種形式的預(yù)條件矩：塊對角形式、塊 三角形式和限制預(yù)條件矩陣等。Golub和Greif等人在2003年提出了一種增廣型預(yù)處理 鞍點(diǎn)問題。接著，各種增廣型預(yù)條件子也相繼被提出。塊三角預(yù)條件子首先由Bramble和 Pasciak于1988年提出：（1，1)和（1，2)塊分別為A和Βτ, (2, 2)塊為C的Schur余。顯 然，預(yù)處理后矩陣的特征值全為1，因此用GMRES等迭代法很快就能收斂到近似解。同塊 對角預(yù)條件子一樣，求（1，1)塊的逆是非常耗時(shí)的，因此近似的塊三角預(yù)條件子相繼被提 出。Simoncini和Cao對塊三角預(yù)處理矩陣進(jìn)行了細(xì)致的譜分析。Keller在2000年提出了 約束預(yù)條件子，對約束預(yù)處理后矩陣的特征值分布，特征向量等進(jìn)行分析，描述了GMRES 等Krylov子空間迭代法的收斂行為。約束預(yù)條件子的結(jié)構(gòu)和原系數(shù)矩陣結(jié)構(gòu)一致，但是 (1，1)塊被A的近似所取代。Cao和Dollar進(jìn)一步研究了約束預(yù)條件子的性質(zhì)，并推廣 到廣義鞍點(diǎn)問題。另外，Bai等人對廣義鞍點(diǎn)系統(tǒng)提出了HSS和PSS預(yù)條件子，并用預(yù)處 理的Krylov子空間方法求解。Pan等人根據(jù)Bai等人的PSS方法提出了一種DPSS預(yù)條件 子，并指出當(dāng)?shù)蜃于呌讴枙r(shí)，預(yù)處理后的矩陣的特征值一部分趨于〇,其余的趨于2。 本發(fā)明主要基于兩類特殊線性鞍點(diǎn)問題特殊結(jié)構(gòu)，設(shè)計(jì)有效的預(yù)處理子，對已有文獻(xiàn)所做 的工作大致總結(jié)如下：
[0009] 1)離散化混合型時(shí)諧Maxwell方程離散產(chǎn)生的鞍點(diǎn)問題
[0010] 由有限元離散的混合型時(shí)諧Maxwell方程：求u和p使得
[0011]
[0012]
[0013] w
[0014] uXn = 0 in Ω，
[0015] p = 0 in Ω，
[0016] 其中，ΩεΚ3是具有連通邊界0O的單連通多面體區(qū)域，n是3Ω上的外單位法向； u和ρ分別是電磁域和Lagrange乘子；波數(shù)滿足k2= ω2ε μ，其中ω彡〇，ε和μ分別 是頻率，電容率和滲透參數(shù)。
[0017] 解決Maxwell方程⑵通常采用離散的Helmholtz分解方法，代數(shù)多重網(wǎng)格方法 和鞍點(diǎn)系統(tǒng)求解方法等?；诎包c(diǎn)系統(tǒng)求解方法，利用第一類最低階的N6d6leC有限元近 似電磁域，用標(biāo)準(zhǔn)的點(diǎn)有限元來近似乘子，離散Maxwell方程可得如下鞍點(diǎn)問題：
[0018]
(3)
[0019] 其中wef和eIT,且對應(yīng)于離散的curl-curl算子的矩陣Je 對稱半正定 具有m維零空間，對應(yīng)于離散的發(fā)散算子的忍是滿秩矩陣，是向量質(zhì)量矩 陣。
[0020] 常用的Krylov子空間方法有CG法，MINRES(minimumresidualmethod)以及 GMRES法等。若波數(shù)k2>0,顯然、的（1，1)塊是不定矩陣，當(dāng)A-k2M變得越不定時(shí)，問題變 得越難解。最近，利用類似文中的譜等價(jià)性質(zhì)，Greif和Sch0_tzau對鞍點(diǎn)線性系統(tǒng)（3) 提出了兩種免增廣和免Schur余的塊對角預(yù)處理子
[0021]
[0022]
[0023]
[0024] 其中是離散的拉普拉斯算子，k2〈l。顯然，?和荇是免增廣和免Schur 余的，作者證明了預(yù)處理矩陣特征值的分布非常聚集；數(shù)值試驗(yàn)也說明了預(yù)條件的MINRES 迭代法迭代次數(shù)的變化對網(wǎng)格的細(xì)化和較小的波數(shù)不敏感。但是提出的預(yù)處理矩陣必須 k2〈l，這是非?？量痰?，而且，對具有非奇異（1，1)塊的通常的鞍點(diǎn)問題，快三角預(yù)處理 技術(shù)[4-6]要比塊對角預(yù)處理技術(shù)更有效，基于此，本發(fā)明將根據(jù)鞍點(diǎn)問題（3)的結(jié)構(gòu)特 征，進(jìn)一步研究并設(shè)計(jì)結(jié)構(gòu)化的免增廣和免Schur余的塊三角預(yù)處理子。
[0025]2)具有奇異（1，1)塊的非對稱鞍點(diǎn)問題
[0026] 考慮如下大型稀疏的2 X 2塊線性方程組：
[0027]
(6)
[0028] 其中奇異且具有高維的零空間，m彡η。鞍點(diǎn)系統(tǒng)（6)的具 體應(yīng)用請參考文獻(xiàn)。這里總假設(shè)^是非奇異的。
[0029] 特別指出，矩陣為是非奇異的當(dāng)且僅當(dāng)rank(B) =rank(C) =m和K/2AVC2非奇 異，其中VB2和Ve2分別是矩陣B和C的零空間的一組基。而且，的非奇異意味著null(A)Πnull(C) = {0}和null(AT)Πnull(B) = {0}。
[0030] 最近，基于鞍點(diǎn)系統(tǒng)（6)是對稱的情況，Greif和Sch0t,zau給出了 一種免 Schur余的塊對角預(yù)處理子：
[0031]
(7)
[0032] 而且他們指出假如A有m維的零空間時(shí)，預(yù)處理子僅有兩種不同的特征值 1和-1。接著，Cao基于預(yù)處理子心將鞍點(diǎn)系統(tǒng)推廣到求解非對稱鞍點(diǎn)系
[0033] 統(tǒng)（6)，提出的預(yù)處理子為：
[0034]
(8)
[0035] 其中妒6:藤_非奇異，且滿足A+BTWt可逆。而且Cao指出假如A有m維的零空 間時(shí)，預(yù)處理子化1;僅有兩種不同的特征值1和-1。而且，Cao也分析了另外兩種廣義 的塊三角預(yù)處理子：
[0036]
[0037] 進(jìn)一步指出，假如A有m維的零空間時(shí)，預(yù)處理子Γ^+孓和分別僅有三種 不同的特征值
> 而且，預(yù)處理子TAug比Τ_+更有效。
[0038] 2012年，Cao基于文獻(xiàn)的思想，給出了一種廣義的塊三角預(yù)處理子：
[0039]
(10)
[0040]其中k和j(辛0)是兩個(gè)實(shí)參數(shù)。而且，Cao分析了預(yù)處理子的特征值分 布，相應(yīng)的特征向量和最小多項(xiàng)式。最后指出，只有當(dāng)k= 2,j= -1時(shí)，預(yù)處理子 才僅有一種特征值1，這樣為了達(dá)到最優(yōu)效果時(shí)就限制了參數(shù)的選取?；诖?，本發(fā)明將 根據(jù)鞍點(diǎn)問題（6)的結(jié)構(gòu)特征，繼續(xù)研究并設(shè)計(jì)一種更加廣義的結(jié)構(gòu)化預(yù)處理子，并且只 需參數(shù)滿足一定條件時(shí)，預(yù)處理子的特征值更加聚集。

【發(fā)明內(nèi)容】

[0041] 針對上述情況，為克服現(xiàn)有技術(shù)之缺陷，本發(fā)明提供一種特殊鞍點(diǎn)問題的高效預(yù) 處理方法。
[0042] 本發(fā)明的技術(shù)方案是，（1)對由離散化混合型時(shí)諧Maxwell方程離散產(chǎn)生的鞍點(diǎn) 問遒
[0043] 根據(jù)此鞍點(diǎn)問題的結(jié)構(gòu)特征，設(shè)計(jì)兩種結(jié)構(gòu)化免增廣和免Schur余的塊三角預(yù)處 理子：
[0044]
[0045]
[0046]
[0047] (2)對設(shè)計(jì)的預(yù)處理子（11)和（12)，分析其構(gòu)造及應(yīng)用代價(jià)和已有的免增廣和免 Schur余的塊對角