冉育紅，李存吉，殷俊鋒

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

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

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

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

1 問題的來源

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

1.1 非局部擴散模型及其快速配置法

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

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

1.2 剛度矩陣的元素估計

1.3 剛度矩陣的結構分析

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

2 預處理GMRES算法

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

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

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

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

3 數(shù)值算例

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

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

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

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

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

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

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

4 總結與展望

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

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