現(xiàn)在科學(xué)計(jì)算方法

1、實(shí)驗(yàn)報(bào)告6 / 6實(shí)驗(yàn)時(shí)間：2012.10.23地點(diǎn)：理學(xué)院機(jī)房報(bào)告人：王振宇題目共軛梯度法解線性方程組實(shí)驗(yàn)?zāi)繕?biāo)運(yùn)用matlab實(shí)現(xiàn)共軛梯度法解線性方程組數(shù)學(xué)原理共軛梯度法又稱(chēng)共軛斜量法，是解線性代數(shù)方程組和非線性方程組的一種數(shù)值方法，例如對(duì)線性代數(shù)方程組Ax二？, (1)式中A為n階矩陣，x和？為n維列向量,當(dāng)A對(duì)稱(chēng)正定時(shí),可以證明求(1)的解X*和求二次泛函代)丄(他陽(yáng))-的極小值問(wèn)題是等價(jià)的。此處(x, y )表示向量x和y的內(nèi)積。由此，給定了初始向量x(0)，按某一方向去求 式取極小值的點(diǎn)x(1),就得到下一個(gè)迭代值x(2),再由x(2)出發(fā)，求x(3)等等，這樣來(lái)逼近 x*。若取求極

2、小值的方向?yàn)镕在x(k=1,2，)處的負(fù)梯度方向就是所謂最速下降法，然而理論和實(shí)際計(jì)算表明這個(gè)方法的收斂速度較慢， 共軛梯度法則是在 x(k-1)處的梯度方向r(k-1)和這一步的修正方向p(k-1)所構(gòu)成的二維平面內(nèi)，尋找使F減小最快的方向作為下一步的修正方向p(k),即求極小值的方向p (其第一步仍取負(fù)梯度方向)。計(jì)算公式為p 00 =產(chǎn)= / - Ax s公式再逐次計(jì)算富£防N黑f+p血*子門(mén)一血叫-W初>4射鞏*4叭尸維-。丄尹(2" - jJ fe>十些戸Z ,公式(k=1,2,)。可以證明當(dāng)i乞時(shí)，,朋)叫(月嚴(yán))叫公式從而平p(1),p(2)形成一

3、共軛向量組;r(0),r(1),形成一正交向量組。后者說(shuō)明若沒(méi)有舍入誤差的話，至多n次迭代就可得到(1)的精確解。然而在實(shí)際計(jì)算中，一般都有舍入誤差，所以r(0),r(1),并不真正互相正交，而r(0) , r(1),等也只是逐步逼近(1)的真解，故一般將共軛梯度法作為迭代法來(lái)使用。近來(lái)在解方程組(1)時(shí),常將共軛梯度法同其他一些迭代法結(jié)合作用。特別是對(duì)病態(tài)方程組這種方法往往能收到比較顯著的效果。其方法是選取一對(duì)稱(chēng)正定矩陣B并進(jìn)行三角分解，得 B=LLT。將方程組(1)化為h y =b, (3)此處 y=ITx,b=l-1?,h=l -1AI-T,而lt = (r1)1 =公式再對(duì)(3)用共軛

4、梯度法，計(jì)算公式為= jf - 蘇 Q i Pq "=尸 Q I如"耳"(環(huán)如),祝+1 心十阻現(xiàn).厲切滬畑皿(畑訕公式k=0,1,2,)適當(dāng)選取B,當(dāng)B很接近A時(shí)，h的條件數(shù)較之 A大大減小,從而可使共軛梯度法的收斂速度大為加快，由一些迭代法的矩陣分裂A=M-N，可選取M為這里的B，例如對(duì)稱(chēng)超松弛迭代(SSOR)，強(qiáng)隱式迭代(SIP)等,這類(lèi)方法常稱(chēng)為廣義共軛梯度法或 預(yù)條件共軛梯度法，它也可用于解代數(shù)特征值問(wèn)題。四.程序設(shè)計(jì)(請(qǐng)將Matlab代碼打印頁(yè)附在該報(bào)告扉頁(yè))原方法fun cti on x = cgm(A,b)n=len gth(b);x=zeros

5、(n ,1);r=b-(A*x);p=r;while (norm(r)=0)s=r;a=(r'*r)/(p'*A*p);x=x+a*p;r=r-a*A*p;c=(r'*r)/(s'*s);p=r+c*p;endend>> A=0 1 -1 2; 1 0 -1 1; -2 1 0 1;b=1 0 1'>> cgm(A,b)ans =0011>> A=23.6 -6.8 4.5 6.3 7.2 8.6;4.2 7.62 0.8 -3.1 -4.5 5.3;9.8 -6.8 45.23.07 -4.37 7.82;8.4 -

6、9.05 24.4 2.65 6.53 -3.64;-9.5 4.54 15.8 4.89 4.8 -7.84;26.32 3.08 4.6 -3.5 -7.8 25.6;b=-9.5 8.12 42.3 -21.3 4.7 3.47'>> cgm(A,b)進(jìn)入死循環(huán)修改后function x = rcgm(A,b,tol)x = b;r = b - A*x;p = r;for k = 1:numel(b);z = A*p;a= (r'*r)/(p'*z);x = x + a*p;s = r'*r;r = r - a*z;if ( norm(r) &

7、lt; tol ) return ;endB = (r'*r)/s;p = r + B*p;end>>A=23.6 -6.8 4.5 6.3 7.2 8.6;4.2 7.62 0.8 -3.1 -4.5 5.3;9.8 -6.8 45.23.07 -4.37 7.82;8.4 -9.05 24.4 2.65 6.53 -3.64;-9.5 4.54 15.8 4.89 4.8 -7.84;26.32 3.08 4.6 -3.5 -7.8 25.6;>>b=-9.5 8.12 42.3 -21.3 4.7 3.47'>>tol=1e-6 >> rcgm(A,b,tol)ans =1.3280783458807260.351356154412808-0.1573211537983757.614700330721686-8.551039777438135-2.808210189216937五結(jié)果分析（結(jié)論）依照理論直接寫(xiě)出的算法非常簡(jiǎn)單，但是循環(huán)過(guò)程時(shí)間復(fù)雜度非常高，求解 6 階方程