關(guān)于非局部擴(kuò)散模型的一種快速預(yù)處理算法

冉育紅，李存吉，殷俊鋒

（1. 西北大學(xué)數(shù)學(xué)學(xué)院，陜西西安710127；2. 西北大學(xué)數(shù)學(xué)學(xué)院，陜西西安710127；3. 同濟(jì)大學(xué)數(shù)學(xué)學(xué)院，上海200092）

對(duì)于自然界那些用整數(shù)階擴(kuò)散方程不太好描述的現(xiàn)象，如反常擴(kuò)散和多孔介質(zhì)流等[1?3]，除了可以用分?jǐn)?shù)階擴(kuò)散方程建模以外，非局部擴(kuò)散模型也是一種可供選擇的建模方法。非局部擴(kuò)散理論對(duì)間斷性問題及其他連續(xù)體變形問題給出了一個(gè)恰當(dāng)?shù)拿枋?，而裂縫、斷裂等奇異性問題在經(jīng)典理論方程下無法給出恰當(dāng)描述。非局部理論被引入到連續(xù)介質(zhì)力學(xué)當(dāng)中，2000 年Silling 提出了一種叫做peridynamic的非局部擴(kuò)散模型[4?7]，在非局部擴(kuò)散模型框架下，物體內(nèi)部不再有接觸力，并且可以保證建立的方程都保持同一積分形式，對(duì)間斷或其他奇異性連續(xù)體的變形問題提供了新的研究思路。

然而，利用一般的數(shù)值方法數(shù)值離散非局部擴(kuò)散模型，通常得到一個(gè)稠密的剛度矩陣，利用直接法求解的計(jì)算量和存儲(chǔ)量非常大，從而阻礙了其廣泛應(yīng)用。最近Wang C 等提出了一種快速的配置法數(shù)值離散變系數(shù)非局部擴(kuò)散模型[8]，分析了剛度矩陣的結(jié)構(gòu)，并利用直接法和共軛梯度平方法求解了變系數(shù)非局部擴(kuò)散模型經(jīng)過快速配置方法數(shù)值離散得到的線性方程組。因?yàn)橹苯臃ǖ挠?jì)算量太大，共軛梯度平方法不穩(wěn)定，所以這兩個(gè)方法不實(shí)用。

由于變系數(shù)非局部擴(kuò)散模型經(jīng)過快速配置方法數(shù)值離散得到的線性方程組系數(shù)矩陣為非對(duì)稱的，因此可以用間接法中的GMRES 方法[9]求解。在實(shí)際求解過程中，系數(shù)矩陣的壞條件數(shù)會(huì)導(dǎo)致GMRES方法收斂的很慢甚至不收斂，而提高收斂速度最常用的方法是預(yù)處理技術(shù)[10?13]。一個(gè)好的預(yù)處理子應(yīng)該至少保持某種特殊結(jié)構(gòu)，并使得預(yù)處理后的系數(shù)矩陣的條件數(shù)很小。在文獻(xiàn)［8］的研究基礎(chǔ)上，本文采用了Toeplitz 及循環(huán)預(yù)處理GMRES 方法求解了變系數(shù)非局部擴(kuò)散模型經(jīng)過快速配置方法數(shù)值離散得到的線性方程組。

1 問題的來源

本節(jié)給出了所需要求解的線性方程組，即變系數(shù)非局部擴(kuò)散模型經(jīng)過快速配置法數(shù)值離散后所得到的離散方程。

1.1 非局部擴(kuò)散模型及其快速配置法

變系數(shù)非局部擴(kuò)散模型的一類體積約束非局部Dirichlet邊值問題，可以表示如下[5]：

式中：δ ＞0表示非局部擴(kuò)散現(xiàn)象擴(kuò)散的水平方向的范圍；α(?)用來表示可變的擴(kuò)散系數(shù)，并假設(shè)其有正的上下界且光滑。核被定義為

1.2 剛度矩陣的元素估計(jì)

1.3 剛度矩陣的結(jié)構(gòu)分析

從參考文獻(xiàn)［8］可知，剛度矩陣A可以分解為

2 預(yù)處理GMRES算法

本文利用快速配置法數(shù)值求解變系數(shù)非局部擴(kuò)散模型最終歸結(jié)為線性方程組式（6）的求解，其中系數(shù)矩陣A由式（18）所定義，右端項(xiàng)f由式（7）所定義。由于系數(shù)矩陣A 為非對(duì)稱矩陣，所以可以利用GMRES 方法求解線性方程組式（6）。為了提高GMRES 方法的收斂率，本節(jié)構(gòu)造了系數(shù)矩陣A 的Toeplitz 及循環(huán)預(yù)處理子，提出了預(yù)處理GMRES 算法。首先，本文給出幾個(gè)基本定義：

由于P1為Toeplitz矩陣，在PGMRES算法中，求解以P1為系數(shù)矩陣的殘量方程P1z=r的計(jì)算量太大。為減少計(jì)算量，我們進(jìn)一步用P1的Strang 循環(huán)近似矩陣逼近P1，令C為Toeplitz 矩陣T的Strang 循環(huán)逼近矩陣，即T≈C，則得到剛度矩陣A的如下預(yù)處理子

由于預(yù)處理子P2為循環(huán)矩陣，則P2可以通過離散的傅里葉變換進(jìn)行對(duì)角化P2=F*ΛFn，所以只需要一次FFT 和一次逆FFT 變換可求出，PGMRES算法中殘量方程P2z=r的解z=F*Λ?1Fn r，計(jì)算量?jī)H為O(NlogN).

同時(shí)，在以上PGMRES 算法中還涉及到Toeplitz 矩陣與向量相乘，由于Toeplitz 矩陣與向量的乘積可以嵌入到2(N- 1) ×2(N- 1)階循環(huán)矩陣與向量相乘中，也可以通過FFT計(jì)算量得到，計(jì)算量仍是O(NlogN),所以PGMRES 算法每次迭代的計(jì)算量為O(NlogN)[17?18]。

3 數(shù)值算例

通過數(shù)值實(shí)驗(yàn)來驗(yàn)證本文構(gòu)造的剛度矩陣A的循環(huán)預(yù)處理子P2的有效性，分別利用GMRES方法、預(yù)處理GMRES 方法求解線性方程組式（6），其中系數(shù)矩陣和右端項(xiàng)分別由式（18）和式（7）定義。所有迭代法都以零向量作為初始向量，迭代終止條件為

表1 算例1的數(shù)值結(jié)果Tab.1 The numerical results for Example 1

從表1和表2可以得出以下結(jié)論：三種算法算出來的數(shù)值近似解是相同的；預(yù)處理GMRES算法需要的迭代步數(shù)和計(jì)算時(shí)間比GMRES 算法需要的迭代步數(shù)和計(jì)算時(shí)間少；P2做預(yù)處理子時(shí)需要的計(jì)算時(shí)間在這三種算法中是最短的。

圖1 和圖2 當(dāng)中的橫坐標(biāo)表示矩陣特征值的實(shí)部，縱坐標(biāo)表示特征值的虛部。從這兩個(gè)圖中，發(fā)現(xiàn)預(yù)處理矩陣的特征值分布比原始矩陣的特征值分布集中，說明預(yù)處理矩陣的條件數(shù)比原始矩陣的條件數(shù)小，故預(yù)處理GMRES方法的迭代步數(shù)比GMRES方法的迭代步數(shù)少。

表2 算例2的數(shù)值結(jié)果Tab.2 The numerical results for Example 2

圖1 算例1中N為300時(shí)，矩陣特征值分布圖Fig.1 When N is 300 in Example 1,the eigenvalue distribution of the matrix

圖2 算例2中N為300時(shí)，矩陣特征值分布圖Fig.2 When N is 300 in Example 2,the eigenvalue distribution of the matrix

4 總結(jié)與展望

利用快速配置法數(shù)值離散變系數(shù)非局部擴(kuò)散模型所得到的一類具有Toeplitz 結(jié)構(gòu)的線性方程組的求解，針對(duì)系數(shù)矩陣的特殊結(jié)構(gòu)，構(gòu)造了一個(gè)Toeplitz 預(yù)處理子和循環(huán)預(yù)處理子，并將其應(yīng)用于GMRES方法中，即得到預(yù)處理GMRES方法。最后利用數(shù)值實(shí)驗(yàn)驗(yàn)證了預(yù)處理GMRES方法的有效性。

非局部擴(kuò)散模型的數(shù)值方法有待進(jìn)一步研究，接下來可以考慮從數(shù)值離散方法入手，提出新的離散方法，保證線性系統(tǒng)更易求解。