江蘇省南京林業(yè)大學(xué)理學(xué)院信息與計(jì)算科學(xué)系 錢耀飛

雅可比迭代法求解稀疏矩陣

江蘇省南京林業(yè)大學(xué)理學(xué)院信息與計(jì)算科學(xué)系 錢耀飛

本文在求解大型稀疏線性代數(shù)方程組時(shí)，通過(guò)Mathematica軟件，判斷雅可比迭代法的譜半徑和迭代精度達(dá)到的迭代次數(shù)。

雅可比迭代；稀疏矩陣；譜半徑；迭代次數(shù)

給定n階線性代數(shù)方程組Ax=b，其中：

當(dāng)n=5，10，20，50，100，200時(shí)，判斷雅可比迭代法的譜半徑和迭代精度達(dá)到時(shí)的迭代次數(shù)。

對(duì)于本題中給定的三對(duì)角線性方程組Ax=b，雅可比迭代法的迭代格式為：

由此可以計(jì)算編寫程序?yàn)椋?/p>

n=10；

b=ConstantArray[0，n]；

b[[1]]=1；b[[n]]=1；

x0=ConstantArray[0，n]；

x1=x0；err=1；k=0；

While[err＞10^(-8)&&k＜25000，

x1[[1]]=(b[[1]]+x0[[2]])/2；

Do[x1[[i]]=(b[[i]]+x0[[i-1]]+x0[[i+1]])/2，{i，2，n-1}]；

x1[[n]]=(b[[n]]+x0[[n-1]])/2；

err=Max[Abs[x1-x0]]；

x0=x1；

k++]；Print[“k=”，k，”err=”，N[err]]；N[x1]

k=375err=9.74235＊10^-9(＊迭代次數(shù)以此迭代的誤差＊)

{1.，1.，1.，1.，1.，1.，1.，1.，1.，1.}（＊求得的解＊）

我們可以看到，當(dāng)n=10時(shí)，雅可比迭代法在迭代了k=375次后，達(dá)到了題目給的精度要求。我們假設(shè)A=D-L-U，其中D，-L，-U分別為A的對(duì)角線元素，嚴(yán)格下三角部分和嚴(yán)格上三角部分，則可以得到雅可比迭代法的迭代矩陣為J=D-1(L+D)。

給出下面的程序，計(jì)算上述雅可比迭代的迭代矩陣的譜半徑。

n=10；

A=SparseArray[{Band[{2，1}]-＞-1，Band[{1，1}]-＞2，

Band[{1，2}]-＞-1}，n]；

d=Diagonal[A]；

（＊構(gòu)造矩陣A，對(duì)角線元素為2，下次和上次對(duì)角線元素均為-1＊）

Diag=DiagonalMatrix[d]；

（＊Diag表示構(gòu)造的對(duì)角矩陣＊）

L=-A＊SparseArray[{i_，j_}/；j＜i-＞1，{n，n}]；

U=-A＊SparseArray[{i_，j_}/；j＞i-＞1，{n，n}]；

A==Diag-L-U(＊檢驗(yàn)矩陣分解是否正確＊)

True(＊結(jié)果為True，說(shuō)明矩陣分解是正確的＊)

J=Inverse[Diag].(L+U)；（＊計(jì)算雅可比迭代矩陣＊）

Max[Abs[Eigenvalues[N[J]]]]（＊計(jì)算雅可比迭代矩陣的譜半徑＊）

0.959493

當(dāng)n等于其他值時(shí)，雅可比方法的迭代次數(shù)和迭代矩陣譜半徑如下表所示：

雅可比迭代法的迭代次數(shù)和迭代矩陣的譜半徑與方程組階數(shù)n的關(guān)系

?

由上表可知，當(dāng)雅可比迭代法的迭代次數(shù)隨著n的增大而急劇增加，誤差也越來(lái)越小，原因在于迭代矩陣半徑越來(lái)越小，接近于1。因此我們現(xiàn)在核心的問(wèn)題就是減少譜半徑，從而減少迭代的次數(shù)。

[1]李慶陽(yáng)，王能超，易大義.數(shù)值分析[M].北京：清華大學(xué)出版社，2008.

[2]王同科，張東麗，王彩華.Mathematica與數(shù)值分析實(shí)驗(yàn)[M].北京：清華大學(xué)出版社，2011.

book=38,ebook=40