高等數(shù)值分析課程設(shè)計--求解線性方程組ax=b的極小化方法比較

1、<p><b>　　課程考核論文</b></p><p>　　課程名稱： 高等數(shù)值分析 </p><p>　　論文題目： 求解線性方程組Ax=b的極 </p><p>　　小化方法比較 </p><p>　　姓 名： <

2、/p><p>　　學(xué) 號： </p><p>　　成 績： </p><p><b>　　摘要</b></p><p>　　本文對求解線性方程組Ax=b的極小化方法進行了理論上的闡述，并且選取三種具有代表性的方法：最速下降法、共軛梯度法（CG）

3、、預(yù)處理共軛梯度法（PCG），使用Matlab編程并分別求解相同的線性方程組，在準(zhǔn)確性和收斂速度方面進行比較。結(jié)果表明，如果預(yù)處理矩陣選擇得當(dāng)，使用預(yù)處理共軛梯度法（PCG）效果最好。</p><p><b>　　極小化方法</b></p><p>　　極小化方法的基本原理是變分原理。</p><p>　　設(shè)A對稱正定，求解的方程組為</

4、p><p><b>　　(1.1)</b></p><p>　　其中，，?？紤]2次函數(shù)，定義為</p><p><b>　　(1.2)</b></p><p><b>　　有如下性質(zhì)</b></p><p><b>　?、?對一切，</b&g

5、t;</p><p><b>　　(1.3)</b></p><p><b>　?、?對一切，</b></p><p><b>　　(1.4)</b></p><p>　?、?設(shè)為(1.1)的解，則</p><p><b>　　而且對一切，有&

6、lt;/b></p><p><b>　　(1.5)</b></p><p>　　則有定理：設(shè)A對稱正定，則為（1.1）解的充分必要條件是滿足</p><p>　　求，使取最小，這就是求解等價于方程組（1.1）的變分問題。求解的方法一般是構(gòu)造一個向量序列，使。</p><p><b>　　最速下降法<

7、;/b></p><p>　　可通過求泛函的極小點來獲得式（1.1）的解。為此，可以從任一出發(fā)，沿著泛函在處下降最快的方向搜索下一個近似點，使得在該方向上達到極小值。</p><p>　　由多元微積分知識，在處下降最快的方向是在該點的負(fù)梯度方向，經(jīng)過簡單計算，得：</p><p><b>　　,</b></p><p&

8、gt;　　此處，也稱為近似點對應(yīng)的殘量。取，確定的值使取得極小值。由式（1.4），令</p><p><b>　　,</b></p><p><b>　　得</b></p><p><b>　　從上式求出并記為：</b></p><p><b>　　(1.6)<

9、;/b></p><p>　　綜合上述推導(dǎo)過程，可得最速下降法的算法描述：</p><p>　?、?給定初始點，容許誤差，置。</p><p>　　⑵ 計算，若，停算，輸出作為近似解。</p><p>　?、?按式（1.6）計算步長因子，置，，轉(zhuǎn)步驟⑵。</p><p>　　最速下降法在理論上是可以收斂的，但是收

10、斂可能會比較緩慢，而且由于舍入誤差存在，實際算出的與理論上的最速下降方向可能有很大差別，使計算時出現(xiàn)數(shù)值不穩(wěn)定性。最速下降法在現(xiàn)代科學(xué)與工程中不再具有實用價值。</p><p><b>　　共軛梯度法（CG）</b></p><p>　　當(dāng)線性方程組Ax=b的系數(shù)矩陣A是對稱正定矩陣，可以采用共軛梯度法對該方程組進行求解，可以證明，式(1.7)所示的n元二次函數(shù)<

11、;/p><p><b>　　(1.7)</b></p><p>　　取得極小值點x*是方程Ax=b的解。共軛梯度法是把求解線性方程組的問題轉(zhuǎn)化為求解一個與之等價的二次函數(shù)極小化的問題。從任意給定的初始點出發(fā)，沿一組關(guān)于矩陣A的共軛方向進行線性搜索，在無舍入誤差的假定下，最多迭代n次（其中n為矩陣A的階數(shù)），就可求得二次函數(shù)的極小點，也就求得線性方程組Ax = b的解。其迭

12、代格式為公式(1.8)。</p><p><b>　　(1.8)</b></p><p>　　共軛梯度法中關(guān)鍵的兩點是迭代格式(21)中最佳步長和搜索方向的確定。其中可以通過一元函數(shù)極小化來求得，其表達式為公式(22)；取，則，要求滿足，可得到的表達式(23).</p><p><b>　　(1.9)</b></p&

13、gt;<p><b>　　(1.10)</b></p><p>　　經(jīng)過一系列的證明和簡化，最終可得共軛梯度法的計算過程如下</p><p>　?、?給定初始點，容許誤差，計算，，置。</p><p><b>　?、?計算步長因子</b></p><p><b>　　置，。&

14、lt;/b></p><p>　?、?若，停算，輸出作為近似解。</p><p><b>　?、?計算</b></p><p><b>　　置，，轉(zhuǎn)步驟⑵。</b></p><p>　　1.3 預(yù)處理共軛梯度法（PCG）</p><p>　　由于計算機存在舍入誤差，故

15、前述的共軛梯度法得到的向量會逐漸失去正交性，因而不大可能在n步之內(nèi)得到原方程組的精確解。此外，遇到求解大規(guī)模的線性方程組，即使能夠在n步收斂，收斂速度也不盡如人意?？梢圆捎妙A(yù)處理技術(shù)改善收斂性質(zhì)，加快收斂速度。</p><p>　　預(yù)處理就是對原方程組進行顯式或隱式的修正，使得修正后的方程組能夠更有效地求解。即構(gòu)造一個預(yù)處理矩陣M，然后用迭代法求解如下的同解方程組</p><p><

16、b>　　(1.11)</b></p><p><b>　　或</b></p><p><b>　　(1.12)</b></p><p>　　相應(yīng)得到的算法分別稱為左預(yù)處理迭代法或右預(yù)處理迭代法。將預(yù)處理技術(shù)和共軛梯度法結(jié)合，可以得到預(yù)處理共軛梯度法，算法如下：</p><p>　　

17、⑴ 給定初始點，容許誤差，計算，，置。</p><p><b>　　⑵ 計算步長因子</b></p><p><b>　　置，。</b></p><p>　?、?若，停算，輸出作為近似解。</p><p><b>　?、?計算</b></p><p>&

18、lt;b>　　置，，轉(zhuǎn)步驟⑵。</b></p><p>　　預(yù)處理共軛梯度法成功與否，關(guān)鍵在于預(yù)處理矩陣M 是否選得合適。一個</p><p>　　好的預(yù)處理矩陣應(yīng)該具有如下特征：</p><p><b>　　⑴ M 對稱正定；</b></p><p>　　⑵ M 與A 的稀疏性相近；</p>

19、;<p>　?、?方程組MSI = rI</p><p><b>　　⑷ 易于求解；</b></p><p>　?、?的特征值分布比較“ 集中”，使較小。</p><p>　　預(yù)處理矩陣的構(gòu)造方法有很多，這里介紹本次實驗使用的兩種方法。</p><p><b>　　⑴ 對角預(yù)處理法</b&g

20、t;</p><p>　　若A 是嚴(yán)格對角占優(yōu)的矩陣，則可以選擇為預(yù)處理矩陣，若A是嚴(yán)格分塊對角占優(yōu)矩陣，則可以選擇為預(yù)處理矩陣；但是這樣簡單的選擇，往往不能帶來很好的加速收斂效果。</p><p><b>　　⑵ 矩陣分裂方法</b></p><p>　　此方法的基本思想是將A 分裂為：；通過矩陣和來構(gòu)造預(yù)處理矩陣M，一般的做法是取線性穩(wěn)定迭

21、代法中相應(yīng)的A 的分裂。比如：Jacobi 分裂（），Gauss - Seidei 分裂（，L 是A 的嚴(yán)格下三角部分），SSOR 分裂等，下面介紹SSOR分裂。將A分裂為，其中，為嚴(yán)格下三角矩陣，它的元素是由A 相應(yīng)部分元素取負(fù)號以后構(gòu)成的。取：</p><p><b>　　(1.13)</b></p><p>　　其中（）是參數(shù)，一般取為1，本次實驗中也取為1。&

22、lt;/p><p><b>　　則預(yù)處理矩陣為：</b></p><p><b>　　(1.14)</b></p><p>　　從理論上講，此預(yù)處理矩陣可以使等于的平方根。</p><p><b>　　2 實驗與分析</b></p><p>　　對于線性方程

23、組，首先使用matlab編寫最速下降法程序，并且隨機生成的對稱正定系數(shù)矩陣A和的矩陣b。多次運行之后得到可以求出收斂解的矩陣并保存，則分別使用最速下降法、共軛梯度法、預(yù)處理共軛梯度法解決同一個線性方程組。以下是矩陣部分截圖，由于矩陣太大，不方便放在論文中，但是生成矩陣的程序在第三章中。</p><p><b>　　圖1 系數(shù)矩陣A</b></p><p><b&

24、gt;　　圖2 矩陣b</b></p><p>　　最速下降法的實驗結(jié)果</p><p>　　由于結(jié)果矩陣x為的矩陣，在本文無法展示出截取計算出的最終收斂的結(jié)果，故截取前5個和后5個數(shù)值如下表：</p><p>　　表1 最速下降法的部分結(jié)果</p><p>　　圖3 Matlab顯示的工作區(qū)</p><p&g

25、t;　　可以看出，使用最速下降法，迭代1821次之后得到收斂的結(jié)果。</p><p>　　共軛梯度法的實驗結(jié)果</p><p>　　用共軛梯度法解同樣的線性方程組，截取部分結(jié)果如下：</p><p>　　表2 共軛梯度法的部分結(jié)果</p><p>　　圖4 Matlab顯示的工作區(qū)</p><p>　　可以看出，使用共

26、軛梯度法，迭代112次之后得到收斂的結(jié)果。</p><p>　　預(yù)處理共軛梯度法的實驗結(jié)果</p><p>　　首先將預(yù)處理矩陣M設(shè)為矩陣A的對角陣，得到數(shù)據(jù)如下：</p><p>　　表3 M為A的對角陣時預(yù)處理共軛梯度法的部分結(jié)果</p><p>　　圖5 Matlab顯示的工作區(qū)</p><p>　　再用第二種方

27、法：用SSOR分裂法得到預(yù)處理矩陣M，部分結(jié)果如下：</p><p>　　表4 用SSOR分裂法得到M的預(yù)處理共軛梯度法的部分結(jié)果</p><p>　　圖6 Matlab顯示的工作區(qū)</p><p>　　可以看出，使用預(yù)處理共軛梯度法，若預(yù)處理矩陣M取A的對角陣，需迭代123次之后得到收斂的結(jié)果；若預(yù)處理矩陣M用SSOR分裂法取得，需迭代90次之后得到收斂的結(jié)果。&

28、lt;/p><p>　　比較上述解線性方程組的三種方法：</p><p>　　①最速下降法法是以殘向量作為解的修正方向，收斂速度可能很慢。若對應(yīng)系數(shù)矩陣A的最大、最小特征值，當(dāng)時，收斂很慢；而且當(dāng)很小時，因舍入誤差影響，計算會不穩(wěn)定。由實驗也可以看出，最速下降法需要的迭代次數(shù)最多。</p><p>　?、诠曹椞荻确▽儆谧兎址椒ǖ姆秶笙禂?shù)矩陣A對稱正定。使用CG法求

29、解n階方程組，理論上最多n步便可得到精確解。實驗也證明，只需要112步便可得到精確解，比最速下降法快得多。</p><p>　　③ PCG法的使用效果，在于預(yù)處理矩陣的選擇；若預(yù)處理矩陣為系數(shù)矩陣的對角陣，用此方法效果可能不盡如人意，比如實驗中迭代次數(shù)為123次，比CG法還慢些；而使用SSOR法對矩陣進行預(yù)處理，再用CG法求解方程組會比較快，實驗也可看出，90步迭代后，便可得到達到精確度要求的解。</p&g

30、t;<p>　　3 實驗用到的matlab程序</p><p><b>　　生成隨機矩陣的程序</b></p><p><b>　　N=1000;</b></p><p>　　D = diag(rand(N,1));</p><p>　　U = orth(rand(N,N));<

31、/p><p>　　B = U' * D * U ;</p><p>　　save('date_A.mat','B');</p><p>　　b = rand(N,1);</p><p><b>　　最速下降法程序</b></p><p>　　A 最速下降法主程序

32、</p><p>　　load date_A.mat</p><p>　　load date_b.mat</p><p>　　[x,iter] = mgrad(B,b);</p><p><b>　　B 最速下降法模塊</b></p><p><b>　　%最速下降法</b>

33、</p><p>　　function [x,iter]=mgrad(A,b,x,ep,N)</p><p>　　%用途：用最速下降法解線性方程組Ax=b</p><p>　　%格式：[x,iter]=mgrad(A,b,x,ep,N) 其中A為系數(shù)矩陣，b為右端向量，x為初始向量</p><p>　　%（默認(rèn)零向量），ep為精度（默認(rèn) 1e

34、-6），N為最大迭代次數(shù)（默認(rèn)1000次），返回參數(shù)</p><p>　　%x，iter分別為近似解向量和迭代次數(shù)</p><p>　　if nargin<5, N=1e12;end</p><p>　　if nargin<4,ep=1e-6;end</p><p>　　if nargin<3,x=zeros(size(b)

35、);end</p><p>　　for iter=1:N</p><p><b>　　r=b-A*x;</b></p><p>　　if norm(r)<ep,break;end</p><p>　　a=r'*r/(r'*A*r);</p><p><b>　　x=

36、x+a*r;</b></p><p><b>　　iter</b></p><p><b>　　x</b></p><p><b>　　end</b></p><p><b>　　共軛梯度法程序</b></p><p>

37、　　A 共軛梯度法主程序</p><p>　　load date_A.mat</p><p>　　load date_b.mat</p><p>　　[x,iter]=mcg(B,b)</p><p><b>　　B 共軛梯度法模塊</b></p><p>　　function [x,iter]=

38、mcg(A,b,x,ep,N)</p><p>　　%用途：用共軛梯度法解線性方程組Ax=b</p><p>　　%格式：[x,iter]=mcg(A,b,x,ep,N)其中A為系數(shù)矩陣，b為右端向量，x為初始向量（默認(rèn)</p><p>　　%零向量）ep為精度（默認(rèn)1e-6），N為最大迭代次數(shù)（默認(rèn)10000次），返回參數(shù)x，iter</p><

39、;p>　　%分別為近似解向量和迭代次數(shù)</p><p>　　if nargin <5,</p><p><b>　　N=10000;</b></p><p><b>　　end</b></p><p>　　if nargin<4,</p><p><b

40、>　　ep=1e-6;</b></p><p><b>　　end</b></p><p>　　if nargin<3,x=zeros(size(b));</p><p><b>　　end</b></p><p><b>　　r=b-A*x;</b>&

41、lt;/p><p>　　for iter=1:N</p><p><b>　　rr=r'*r;</b></p><p>　　if iter==1</p><p><b>　　p=r;</b></p><p><b>　　else</b></p&

42、gt;<p>　　beta=rr/rr1;</p><p>　　p=r+beta*p;</p><p><b>　　end</b></p><p><b>　　q=A*p;</b></p><p>　　alpha=rr/(p'*q);</p><p>　

43、　x=x+alpha*p;</p><p>　　r=r-alpha*q;</p><p><b>　　rr1=rr;</b></p><p>　　if (norm(r)<ep),break;</p><p><b>　　end</b></p><p><b>

44、　　end</b></p><p>　　3.4 預(yù)處理共軛梯度法程序</p><p>　　A 得到預(yù)處理矩陣M的程序</p><p>　?、?M為A的對角矩陣</p><p>　　load date_A.mat</p><p>　　M=diag(diag(B));</p><p>

45、　　save('date_juzhenM1.mat','M');</p><p>　　② M由SSOR分裂法得到</p><p>　　load date_A.mat</p><p>　　D=diag(diag(B));</p><p>　　CL=tril(-B);</p><p>　　L

46、=((D-CL)*D)^1/2;</p><p><b>　　M=L*L';</b></p><p>　　save('date_juzhenM.mat','M');</p><p>　　B 預(yù)處理共軛梯度法主程序</p><p>　　load date_A.mat</p>

47、;<p>　　load date_b.mat</p><p>　　[x,iter]=mpcg(B,b)</p><p>　　C 預(yù)處理共軛梯度法模塊</p><p>　　function [x,iter]=mpcg(A,b,x,ep,N)</p><p>　　%用途：用預(yù)處理共軛梯度法解線性方程組Ax=b</p>

48、<p>　　%格式：[x,iter]=mpcg(A,b,x,ep,N)其中A為系數(shù)矩陣，b為右端向量，x為初始向量（默認(rèn)</p><p>　　%零向量）ep為精度（默認(rèn)1e-6），N為最大迭代次數(shù)（默認(rèn)10000次），返回參數(shù)x，iter</p><p>　　%分別為近似解向量和迭代次數(shù)</p><p>　　load date_juzhenM.mat<

49、;/p><p><b>　　L=inv(M);</b></p><p>　　if nargin <5,</p><p><b>　　N=10000;</b></p><p><b>　　end</b></p><p>　　if nargin<4,

50、</p><p><b>　　ep=1e-6;</b></p><p><b>　　end</b></p><p>　　if nargin<3,x=zeros(size(b));</p><p><b>　　end</b></p><p>　　r=

51、L*(b-A*x);</p><p>　　for iter=1:N</p><p><b>　　rr=r'*r;</b></p><p>　　if iter==1</p><p><b>　　p=r;</b></p><p><b>　　else</b

52、></p><p>　　beta=rr/rr1;</p><p>　　p=r+beta*p;</p><p><b>　　end</b></p><p><b>　　q=L*A*p;</b></p><p>　　alpha=rr/(p'*q);</p>

53、;<p>　　x=x+alpha*p;</p><p>　　r=r-alpha*q;</p><p><b>　　rr1=rr;</b></p><p>　　if (norm(r)<ep),break;</p><p><b>　　end</b></p><p