共軛梯度法的預(yù)處理ppt課件

1、共軛梯度法的預(yù)處理,設(shè) A 為 n 階對(duì)稱正定矩陣, 考慮線性方程組 (1) 其中x 是未知向量,b是右端已知向量,一、知識(shí)回顧,1、共軛梯度法 如果我們對(duì)于固定的 ，在 空間中取 個(gè)線性無關(guān)的向量 使 滿足 當(dāng) 時(shí)， 一定是方程 (1)的精確解，這是共軛梯度法的出發(fā)點(diǎn),一、知識(shí)回顧,2、預(yù)處理共軛梯度法 （1）預(yù)處理共軛梯度法的引進(jìn) 共軛梯度法理論上是一種直接方法，即按照這個(gè)算法得到的 應(yīng)當(dāng)是方程 的精確解。但由于其數(shù)值的穩(wěn)定性不好，實(shí)際上由于舍入誤差的影響 不可能滿足原方程。如果把它看成一種迭代法，當(dāng)A條件數(shù)較大時(shí)，它收斂的很慢。因此引進(jìn)了改進(jìn)共軛梯度法的手段預(yù)處理共軛梯度法,一、知識(shí)回顧

2、,2）預(yù)處理共軛梯度法原理 將Ax=b變形為C-1Ax=C-1b。 令 則 要求 當(dāng) 是對(duì)稱正定矩陣時(shí)， 其中 是A的最大特征值， 是A的最小特征值，K(A)是A的條件數(shù),一、知識(shí)回顧,3）預(yù)處理共軛梯度方法 由于預(yù)處理矩陣的構(gòu)造方式很多，所以就形成了多種預(yù)處理方法，總的來說，大致可歸為如下幾種途徑： 1、取預(yù)處理矩陣（也稱預(yù)優(yōu)矩陣）為A的一個(gè)小帶寬部分（如三對(duì)角或?qū)蔷€部分）； 2、通過A的分裂，尤其是線性穩(wěn)定迭代法中的矩陣 A的分裂構(gòu)造預(yù)處理矩陣（如矩陣分裂法）； 3、通過A的各種近似分解得到預(yù)處理矩陣（如不完全分解）； 4、通過矩陣A的多項(xiàng)式構(gòu)造預(yù)處理矩陣； 5、子結(jié)構(gòu)、區(qū)域分裂、EBE

3、預(yù)處理途徑,一、知識(shí)回顧,對(duì)角預(yù)優(yōu)矩陣法 如果對(duì)稱正定線性方程組( 1) 的系數(shù)矩陣 A 的對(duì)角元素相差較大, 則可取預(yù)優(yōu)矩陣C-1為 這是最常用的預(yù)處理方陣,SSOR分裂 將A分裂為 ，其中 , 為嚴(yán)格下三角矩陣，則預(yù)處理矩陣C為： 其中 是參數(shù)，而且有,二、算法,1）共軛梯度法算法 令 對(duì) 直到收斂完成 結(jié)束,二、算法,2）預(yù)處理共軛梯度算法 令 對(duì) 直到收斂完成。 結(jié)束,三、代碼,1）共軛梯度法代碼 function x,n,tol=cg(A,b,ep) r=b;p=r;n=0;x=zeros(length(b),1); while n10 alpha=(r*r)/(p*A*p);r1=

4、r; x=x+alpha*p; tol=norm(r); if norm(r)ep break; end r=r-alpha*A*p; beta=(r*r)/(r1*r1); p=r+beta*p; n=n+1; end,x = 7.85971307544587 0.422926408294999 -0.0735922390236594 -0.540643016894632 0.0106261628540375 n = 5 tol = 6.745229690695623e-07,2）預(yù)處理共軛梯度法代碼 以對(duì)角預(yù)優(yōu)矩陣為例 clc;clear A=0.2 0.1 1 1 0;0.1 4 -1

5、1 -1;1 -1 60 0 -2;1 1 0 8 4;0 -1 -2 4 700; b=1 2 3 4 5;%SSOR此處程序改為 C=chuli1(A); A1=C*A*C; b1=C*b; ep=10(-2); x,n,tol=cg1(A1,b1,ep,C); function C=chuli1(A) D=diag(diag(A); C=D(-1/2,function C=chuli(A,w) D=diag(diag(A);L=tril(A,-1); C=(D-w*L)*D(-1/2,C=chuli(A,0.009); A1=inv(C)*A*inv(C); b1=inv(C)*b,fu

6、nction x,n,tol=cg1(A,b,ep,C) r=b;p=r;n=0;x=zeros(length(b),1); while n10 alpha=(r*r)/(p*A*p);r1=r; x=x+alpha*C*p; tol=norm(r); if norm(r)ep break; end r=r-alpha*A*p; beta=(r*r)/(r1*r1); p=r+beta*p; n=n+1; end,SSOR此處程序改為 x=x+alpha*inv(C)*p,SSOR分裂 x = 7.85971307544587 0.422926408295008 -0.073592239024

7、0462 -0.540643016894626 0.0106261628540363 m= 4 tol = 6.567520686375984e-05,對(duì)角預(yù)優(yōu)矩陣 x= 7.85971307544586 0.422926408295008 -0.0735922390240463 -0.540643016894626 0.0106261628540363 m= 4 tol= 4.731753657562764e-04,此矩陣并沒有顯示出預(yù)處理的優(yōu)勢(shì)，所以我們又取矩陣 A=diag(1 10 100 1000 10000)，輸出結(jié)果如下,五、參考文獻(xiàn),1】蔡大用，白峰杉.現(xiàn)代科學(xué)計(jì)算.北京：科學(xué)出版社，2000