計(jì)算方法線(xiàn)性方程組的直接解法

第六章線(xiàn)性方程組的直接法2/109§2.0引言§2.1Gauss消去法§2.2矩陣分解在解線(xiàn)性方程組中的作用補(bǔ)充向量的范數(shù)與矩陣的范數(shù)§2.3直接法的誤差分析§2.4線(xiàn)性方程組迭代解法

3/109§2.0引言在自然科學(xué)和工程技術(shù)中很多問(wèn)題的解決常常歸結(jié)為解線(xiàn)性代數(shù)方程組。例如電學(xué)中的網(wǎng)絡(luò)問(wèn)題，用最小二乘法求實(shí)驗(yàn)數(shù)據(jù)的曲線(xiàn)擬合問(wèn)題，解非線(xiàn)性方程組問(wèn)題，用差分法或者有限元方法解常微分方程、偏微分方程邊值問(wèn)題等都導(dǎo)致求解線(xiàn)性代數(shù)方程組。而這些方程組的系數(shù)矩陣大致分為兩種，一種是低階稠密矩陣(例如，階數(shù)大約為≤150)，另一種是大型稀疏矩陣(即矩陣階數(shù)高且零元素較多)。

4/109§2.0引言設(shè)有線(xiàn)性方程組Ax=b，其中為非奇異陣，，關(guān)于線(xiàn)性方程組的數(shù)值解法一般有兩類(lèi)：直接法與迭代法。5/1091.直接法就是經(jīng)過(guò)有限步算術(shù)運(yùn)算，可求得方程組精確解的方法(若計(jì)算過(guò)程中沒(méi)有舍入誤差)。但實(shí)際計(jì)算中由于舍入誤差的存在和影響，這種方法也只能求得線(xiàn)性方程組的近似解。本章將闡述這類(lèi)算法中最基本的高斯消去法及其某些變形。這類(lèi)方法是解低階稠密矩陣方程組的有效方法，近十幾年來(lái)直接法在求解具有較大型稀疏矩陣方程組方面取得了較大進(jìn)展。6/1092.迭代法迭代法是用某種極限過(guò)程去逐步逼近線(xiàn)性方程組精確解的方法。迭代法具有需要計(jì)算機(jī)的存貯單元較少、程序設(shè)計(jì)簡(jiǎn)單、原始系數(shù)矩陣在計(jì)算過(guò)程中始終不變等優(yōu)點(diǎn)，但存在收斂性及收斂速度問(wèn)題。迭代法是解大型稀疏矩陣方程組(尤其是由微分方程離散后得到的大型方程組)的重要方法。7/109§2.1Gauss消去法高斯（Gauss）消去法是解線(xiàn)性方程組最常用的方法之一。它的基本思想是通過(guò)逐步消元（行的初等變換），把方程組化為系數(shù)矩陣為三角形矩陣的同解方程組，然后用回代法解此三角形方程組（簡(jiǎn)單形式）得原方程組的解。

8/109§2.1Gauss消去法例：下面我們來(lái)討論一般的解n階方程組的高斯消去法。9/1091.消元將原方程組記為A(1)x=b(1)，其中A(1)=(aij(1))=(aij)，b(1)=b(1)第一次消元。其中注：若a11(1)=0，則在第1列中至少有一個(gè)元素不為0，可交換行后再消元10/109(2)第k次消元。其中注：為減少計(jì)算量，令，則11/109(3)當(dāng)k=n

–1時(shí)得完成第n

–1次消元后得到與原方程組等價(jià)的三角形方程組A(n)x=b(n)注：當(dāng)det(A)≠0時(shí)，顯然有aii(i)≠0，(i=1，…，n)，稱(chēng)為主元素。12/1092.回代求解三角形方程組A(n)x=b(n)，得到求解公式注：求解過(guò)程稱(chēng)為回代過(guò)程。13/1093.Gauss消去法的計(jì)算量以乘除法的次數(shù)為主(1)

消元過(guò)程：第k步時(shí)(n–k)+(n–k)(n–k+1)=(n–k)(n–k+2)共有14/1093.Gauss消去法的計(jì)算量(2)

回代過(guò)程：求xi中，乘n–i次，除1次，共n–i+1次（i=1，…，n–1）共有15/1093.Gauss消去法的計(jì)算量(3)總次數(shù)為注：當(dāng)n=20時(shí)約為2670次，比克萊姆法則9.71020次大大減少。16/1094.說(shuō)明(1)若消元過(guò)程中出現(xiàn)akk(k)=0，則無(wú)法繼續(xù)(2)若akk(k)≠0，但較小，則小主元做除數(shù)將產(chǎn)生大誤差(3)每次消元都選擇絕對(duì)值最大者作主元，稱(chēng)為高斯主元消去法(4)通常第k步選擇——第k列主對(duì)角元以下元素絕對(duì)值最大者作主元（該行與第k行對(duì)調(diào)），稱(chēng)為列主元消去法。17/109例1用舍入三位有效數(shù)字求解線(xiàn)性方程組（準(zhǔn)確解是x1=10.0，x2=1.00）解：1)

不選主元的Gauss消去法計(jì)算結(jié)果：

x1=-10.0，x2=1.01，此解無(wú)效；

2)

按列選主元的Gauss消去法計(jì)算結(jié)果：

x1=10.0，x2=1.00.18/109(5)行標(biāo)度化當(dāng)線(xiàn)性方程組的系數(shù)矩陣的元素大小相差很大時(shí)，則可能引進(jìn)因丟失有效數(shù)字而產(chǎn)生的誤差，并且舍入誤差的影響嚴(yán)重，即使用Gauss主元消去法得到的解也不可靠．為避免這一問(wèn)題，可將系數(shù)矩陣的行元素按比例增減以使元素的變化減?。缬妹啃性氐淖畲竽３撔懈髟?，使它們的模都不大于1，這叫行標(biāo)度化，其目的是要找到真正的主元．19/109(5)行標(biāo)度化如用每行元素的最大模除該行各元素，使它們的模都不大于1，這叫行標(biāo)度化，其目的是要找到真正的主元．消元過(guò)程仍是對(duì)原方程組進(jìn)行的，只不過(guò)在Gauss列主元消去法的算法中，按列選主元ck時(shí)，應(yīng)修改為其中20/1095.算法設(shè)計(jì)輸入A，n，epsFor(k=1ton-1)選主元A(I0，k)——確定行號(hào)p=A(I0，k)若|p|<=epsText=1，退出循環(huán)若I0kT換行（I0k）消元若|A(n，n)|<=epsText=1F回代若ext=1T無(wú)解F輸出解I0=kfor(i=k+1ton)若|A(i，k)|>|A(I0，k)|TI0=ifor(i=k+1ton)A(i，k)=A(i，k)/A(k，k)for(j=k+1ton+1)A(i,j)=A(i,j)-A(i,k)*A(k,j)A(n，n+1)=A(n，n+1)/A(n，n)for(k=n–1to1step-1)w=0for(j=k+1ton)w=w+A(k，j)*A(j，n+1)A(k，n+1)=A(k，n+1)–wA(k，n+1)=A(k，n+1)/A(k，k)for(j=kton)n=A(k，j),A(k，j)=A(I0，j),A(I0，j)=n21/1095.算法設(shè)計(jì)(1)

高斯消去法a=[5-1-1-1-4;-110-1-112;-1-15-18;-1-1-11034];x=[0;0;0;0]n=4fork=1:n-1fori=k+1:na(i,:)=a(i,:)-a(k,:)*a(i,k)/a(k,k)endendforj=n:-1:1

x(j)=(a(j,n+1)-a(j,j+1:n)*x(j+1:n))/a(j,j)end22/109(2)列主元素消去法a=[-0.040.040.123;0.56-1.560.321;-0.241.24-0.280];x=[0;0;0];n=3;fork=1:n-1[c,i]=max(abs(a(k:n,k)));

q=i+k-1ifq~=km=a(q,:);a(q,:)=a(k,:);a(k,:)=mendfori=k+1:na(i,:)=a(i,:)-a(k,:)*a(i,k)/a(k,k)endendforj=n:-1:1

x(j)=(a(j,n+1)-a(j,j+1:n)*x(j+1:n))/a(j,j)end23/109§2.2矩陣分解在解線(xiàn)性方程組中的作用1.原理若A=LR，則Ax=b

LRx=b

若其中L、R為三角形矩陣，則方程組易解定義：(1)稱(chēng)為單位下三角陣(2)設(shè)L為單位下三角陣，R為上三角陣，若A=LR，則稱(chēng)A可三角（LR）分解，并稱(chēng)A=LR為A的三角分解（杜利特爾分解）24/109定理：(1)A=(aij)nn有唯一LR分解

A的順序主子式k≠0，k=1，2，...，n

–1(2)若A=(aij)nn可逆，則存在置換陣P，使PA的n個(gè)順序主子式全不等于0注：由Ax=b

PAx=Pb知，經(jīng)適當(dāng)行交換后，A總可三角分解25/1092.直接三角分解法設(shè)A的各階順序主子式不為0，且A=LR，即(1)主對(duì)角線(xiàn)（含）上邊：當(dāng)列標(biāo)m

>

行標(biāo)k時(shí)，lkm=0

j=k，k+1，...，n；k=1，2，...，n(2)主對(duì)角線(xiàn)下邊：當(dāng)行標(biāo)m

>

列標(biāo)k時(shí)，rkm=0

i=k+1，k+2，...，n；k=1，2，...，n–126/109設(shè)A的各階順序主子式不為0，且A=LR，即(1)主對(duì)角線(xiàn)（含）上邊：當(dāng)列標(biāo)m

>

行標(biāo)k時(shí)，lkm=0

j=k,k+1,...,n；k=1,2,...,n(2)主對(duì)角線(xiàn)下邊：當(dāng)行標(biāo)m>列標(biāo)k時(shí)，rkm=0

i=k+1,k+2,...,n；k=1,2,...,n–1欲求lik和rkj27/109(1)主對(duì)角線(xiàn)（含）上邊：當(dāng)列標(biāo)m

>

行標(biāo)k時(shí)，lkm=0

j=k,k+1,...,n；k=1,2,...,n(2)主對(duì)角線(xiàn)下邊：當(dāng)行標(biāo)m>列標(biāo)k時(shí)，rkm=0

i=k+1,k+2,...,n；k=1,2,...,n–1欲求lik和rkj設(shè)k=1，a1j=l11r1j

r1j=a1j

j=1，2，...，nai1=li1r11

li1=ai1/r11

i=2，3，...，n28/109一般地，最后：r11r12...r1n第1步l21r22...r2n第2步l31l32......ln1ln2rnn第n步29/1093.求解方程組下三角Ly=b：上三角Rx=y：緊湊格式：30/109a=[2100-3;-3-4-1213;123-4;4149-13];b=[10;5;-2;7];y=[0000]’;x=[0000]'n=4;forr=1:na(r,r:n)=a(r,r:n)-a(r,1:r-1)*a(1:r-1,r:n)a(r+1:n,r)=(a(r+1:n,r)-a(r+1:n,1:r-1)*a(1:r-1,r))/a(r,r)endR=triu(a,0)L=eye(n)+tril(a,-1)forr=1:ny(r)=b(r)-L(r,1:r-1)*y(1:r-1)endforj=n:-1:1x(j)=(y(j)-R(j,j+1:n)*x(j+1:n))/R(j,j)end31/109緊湊格式a=[2100-310;-3-4-12135;123-4-2;4149-137];x=[0000]'n=4;forr=1:na(r,r:n+1)=a(r,r:n+1)-a(r,1:r-1)*a(1:r-1,r:n+1)a(r+1:n,r)=(a(r+1:n,r)-a(r+1:n,1:r-1)*a(1:r-1,r))/a(r,r)endR=triu(a,0)forj=n:-1:1x(j)=(a(j,n+1)-R(j,j+1:n)*x(j+1:n))/R(j,j)end32/1094.選主元的三角分解法5.平方根法(1)定理：若A正定，則A可唯一分解為A=LDLT。其中L為單位下三角，D為元素全為正數(shù)的對(duì)角陣。(2)設(shè)A=LDLT則有：33/109(3)為了避免重復(fù)計(jì)算，引進(jìn)中間量tik=likdk，將上式改寫(xiě)為：34/109(4)Ax=b

LDLTx=b

Ly=b，Dz=y，LTx=z求解公式：上述方法稱(chēng)為“改進(jìn)的平方根法”注：不能選主元作行交換——破壞對(duì)稱(chēng)性，必是數(shù)值穩(wěn)定的35/109a=[5-41;-46-4;1-46];b=[2;-1;-1]d=[000]';t=[000;000;000];L=[000;000;000];n=3fork=1:n

d(k)=a(k,k)-t(k,1:k-1)*L(k,1:k-1)'

t(k+1:n,k)=a(k+1:n,k)-t(k+1:n,1:k-1)*L(k,1:k-1)'L(k+1:n,k)=t(k+1:n,k)/d(k)endx=[000]';y=[000]';z=[000]';fori=1:ny(i)=b(i)-L(i,1:i-1)*y(1:i-1)endfori=1:nz(i)=y(i)/d(i)endfori=n:-1:1x(i)=z(i)-L(i+1:n,i)'*x(i+1:n)end36/1096.追趕法在實(shí)際問(wèn)題中，經(jīng)常遇到以下形式的方程組這種方程組的系數(shù)矩陣稱(chēng)為三對(duì)角矩陣。以下針對(duì)該方程組的特點(diǎn)提供一種簡(jiǎn)便有效的算法——

追趕法。37/109設(shè)直接三角分解為：A=LR其中容易看出pi=ci（1，2，…，n

–1），計(jì)算L和R的其余元素的公式為38/109容易看出pi=ci（1，2，…，n

–1），計(jì)算L和R的其余元素的公式為計(jì)算過(guò)程：r1→l2→r2→...→ln→rn39/109有了A的LR分解后，線(xiàn)性方程組Ax=d等價(jià)于兩個(gè)簡(jiǎn)單的方程組：Ly=d，Rx=y于是，計(jì)算y的公式是y1=d1，yi=di

–

liyi–1

（i=2，3，…，n）計(jì)算過(guò)程：r1→y1→l2→r2→y2→...→ln→rn→yn

（追）計(jì)算x的公式是xn=yn/rn，xi=(yi–cixi+1)/ri

（i=n–1，n–2，…，1）(趕）40/109定理：若A為對(duì)角占優(yōu)（diagonallydominant）的三對(duì)角陣，且滿(mǎn)足|b1|>|c1|>0，|bi|≥|ai|+|ci|，aici≠0，i=2，3，…，n-1|bn|>|an|>0.則追趕法可解以A為系數(shù)矩陣的方程組注：1.如果A是嚴(yán)格對(duì)角占優(yōu)陣，則不要求三對(duì)角線(xiàn)上的所有元素非零。2.根據(jù)不等式可知：分解過(guò)程中，矩陣元素不會(huì)過(guò)分增大，算法保證穩(wěn)定。3.運(yùn)算量為O(6n)。41/1097.QR算法42/109補(bǔ)充向量的范數(shù)與矩陣的范數(shù)在線(xiàn)性方程組的數(shù)值解法中，經(jīng)常需要分析解向量的誤差，需要比較誤差向量的“大小”或“長(zhǎng)度”。那么怎樣定義向量的長(zhǎng)度呢？我們?cè)诔醯葦?shù)學(xué)里知道，定義向量的長(zhǎng)度，實(shí)際上就是對(duì)每一個(gè)向量按一定的法則規(guī)定一個(gè)非負(fù)實(shí)數(shù)與之對(duì)應(yīng)，這一思想推廣到維線(xiàn)性空間里，就是向量的范數(shù)或模。用Rn表示n維實(shí)向量空間，用Cn表示n維復(fù)向量空間，首先將向量長(zhǎng)度概念推廣到Rn(或Cn)中。43/1091.向量的范數(shù)向量的范數(shù)可以看作是描述向量“大小”的一種度量.范數(shù)的最簡(jiǎn)單的例子，是絕對(duì)值函數(shù)：有三個(gè)熟知的性質(zhì)：(1)x

0

x

>0

x

=0當(dāng)且僅當(dāng)x=0(2)ax=

a

x

a為常數(shù)(3)

x+y

≤

x

+

y

44/1091.向量的范數(shù)范數(shù)的另一個(gè)簡(jiǎn)單例子是三維歐氏空間的長(zhǎng)度設(shè)x=(x1,x2,x3)，則x的歐氏范數(shù)定義為：歐氏范數(shù)也滿(mǎn)足三個(gè)條件：

x，y

R3，a為常數(shù)(1)

x

≥0，且x=0

x=0(2)

ax

=

a

x

(3)

x+y

≤

x

+

y

前兩個(gè)條件顯然，第三個(gè)條件在幾何上解釋為三角形一邊的長(zhǎng)度不大于其它兩邊長(zhǎng)度之和。因此，稱(chēng)為三角不等式。45/109向量范數(shù)的一般概念：定義1：設(shè)V是數(shù)域F上的向量空間，對(duì)V中任一向量α，都有唯一實(shí)數(shù)α

與之對(duì)應(yīng)，滿(mǎn)足如下三個(gè)條件：1)正定性：α≥0，且α=0

α=02)齊次性：kα=|k|α

，這里k

F3)三角不等式：α+

α

+

則稱(chēng)α為α的范數(shù)。定義了范數(shù)的向量空間稱(chēng)為賦范向量空間.簡(jiǎn)單性質(zhì)：(1)x

0——單位向量(2)||x||=||–x||(3)|||x||–||y|||||x–y||——當(dāng)x

y時(shí)，||x||||y||46/109Cn上的常見(jiàn)范數(shù)有：1)1-范數(shù)

2)2-范數(shù)稱(chēng)為歐氏范數(shù)3)-范數(shù)不難驗(yàn)證，上述三種范數(shù)都滿(mǎn)足定義的條件。注：上述形式的統(tǒng)一：

1

p

47/109定理5：定義在Cn上的向量范數(shù)||x||是變量x分量的一致連續(xù)函數(shù)。（f(x)=||x||）定理6：在Cn上定義的任何兩個(gè)范數(shù)都是等價(jià)的。即存在正數(shù)k1與k2(k1≥k2>0)，對(duì)一切xCn，不等式k1||x||b

||x||a

k2||x||b成立。對(duì)常用范數(shù)，容易驗(yàn)證下列不等式：48/109例1設(shè)A為正定矩陣，xCn，令，則||x||A為向量范數(shù)，稱(chēng)為加權(quán)范數(shù)證明：(1)顯然xA≥0，且xA=0

x=0(2)(3)因?yàn)锳正定，所以可逆P，使A=PTP

所以||x+y||A=||P(x+y)||2=||Px+Py||2||Px||+||Py||2=||x||A+||y||A綜上，||x+y||A為范數(shù)。注：當(dāng)0<p<1時(shí)，不是范數(shù)49/109有了范數(shù)的概念，就可以討論向量序列的收斂性問(wèn)題。定義2：設(shè)給定Cn中的向量序列{xk}，即x0，x1，…，xk，…其中若對(duì)任何i(i=1,2,…,n)都有則向量稱(chēng)為向量序列{xk}的極限，或者說(shuō)向量序列{xk}依坐標(biāo)收斂于向量x*，記為50/109定義2：設(shè)給定Cn中的向量序列{xk}，即x0，x1，…，xk，…其中若對(duì)任何i(i=1,2,…,n)都有則向量稱(chēng)為向量序列{xk}的極限，或者說(shuō)向量序列{xk}依坐標(biāo)收斂于向量x*，記為定理7：向量序列{xk}依坐標(biāo)收斂于x*的充要條件是如果一個(gè)向量序列{xk}與向量x*滿(mǎn)足上式，就說(shuō)向量序列{xk}依范數(shù)收斂于x*，于是便得：向量序列依范數(shù)收斂與依坐標(biāo)收斂是等價(jià)的。51/1092.矩陣的范數(shù)下面我們給出矩陣范數(shù)的一般定義。為簡(jiǎn)單起見(jiàn)，僅討論實(shí)數(shù)域的情形。定義3

若矩陣ARnn的某個(gè)非負(fù)的實(shí)值函數(shù)N(A)=||A||，滿(mǎn)足條件1)正定性：A≥0，且A=0

A=02)齊次性：kA=|k|||A||

kR3)三角不等式：||A+B||

||A||+||B||4)相容性：||AB||||A||||B||則稱(chēng)N(A)是Rnn上的一個(gè)矩陣范數(shù)(或模)。52/1092.矩陣的范數(shù)定義3

若矩陣ARnn的某個(gè)非負(fù)的實(shí)值函數(shù)N(A)=||A||，滿(mǎn)足條件1)正定性：A≥0，且A=0

A=02)齊次性：kA=|k|||A||

kR3)三角不等式：||A+B||

||A||+||B||4)相容性：||AB||||A||||B||則稱(chēng)N(A)是Rnn上的一個(gè)矩陣范數(shù)(或模)。如由Rnn上2-范數(shù)可以得到Rnn中矩陣的一種范數(shù)稱(chēng)為A的Frobenius范數(shù)。||A||F顯然滿(mǎn)足正定性、齊次性及三角不等式。53/109由于在大多數(shù)與估計(jì)有關(guān)的問(wèn)題中，矩陣和向量會(huì)同時(shí)參與討論，所以希望引進(jìn)一種矩陣的范數(shù)，它是和向量范數(shù)相聯(lián)系而且和向量范數(shù)相容的，即5)||Ax||||A||||x||對(duì)任何向量xRn及矩陣ARnn都成立。注：滿(mǎn)足相容性條件的范數(shù)很多，其中最重要的是算子范數(shù)54/109再引進(jìn)一種矩陣的范數(shù)。定義4

設(shè)xRn及矩陣ARnn，||x||v為一種向量范數(shù)(如v=1，2或∞)，相應(yīng)地定義一個(gè)矩陣的非負(fù)函數(shù)可驗(yàn)證||A||v滿(mǎn)足定義3，所以||A||v是Rnn上矩陣的一個(gè)范數(shù)，||A||v還滿(mǎn)足相容性條件5)，稱(chēng)為A的算子范數(shù)，或誘導(dǎo)范數(shù)。注：(1)可以證明：(2)||Ax||是x的連續(xù)函數(shù)，從而存在x*使||x*||=1且||Ax*||=||A||(3)若||.||為算子范數(shù)，I為單位矩陣，則||I||=155/109定理8：設(shè)n階方陣A=(aij)nn，則與||x||1相容的矩陣范數(shù)是與||x||2相容的矩陣范數(shù)是其中1為矩陣ATA的最大特征值。與||x||相容的矩陣范數(shù)是上述范數(shù)分別稱(chēng)為1-范數(shù)，2-范數(shù)和-范數(shù)又稱(chēng)為列模、譜模和行模56/109注：A的范數(shù)||A||與A的特征值間的關(guān)系設(shè)為矩陣A的任一特征值，e為相應(yīng)的特征向量，則e=Ae，因?yàn)閨|||e||=||Ae||||A||||e||，所以||||A||從而得定理9：矩陣A的任一特征值的絕對(duì)值不超過(guò)A的范數(shù)||A||。定義6：矩陣A的諸特征值的最大絕對(duì)值稱(chēng)為A的譜半徑，記為：由定理9可知：57/109例(p138)：設(shè)A=(aij)nn，||A||為其算子范數(shù)，則||A||<1

I–A可逆，且58/109§2.3直接法的誤差分析1.誤差向量與殘差向量設(shè)x為線(xiàn)性方程組Ax=b的準(zhǔn)確解，用數(shù)值算法得到的近似解為x*，則稱(chēng)向量e=x*-x為誤差向量。顯然，||e||可以表示近似解x*的準(zhǔn)確程度：||e||越小，近似解x*越準(zhǔn)確。但由于x不知道，所以無(wú)法計(jì)算e。變通的辦法是考慮殘差向量：r=b–Ax*的范數(shù)的大小。但是，在某些情況下，盡管||r||很小，||e||也可能很大。59/109例線(xiàn)性方程組的準(zhǔn)確解是(1，1)T，若用某種方法得到計(jì)算解(2.000，0.500)，則殘差向量其-范數(shù)||r||=0.0015，而誤差向量其-范數(shù)||e||=1，它是殘差向量的667倍。60/109能否建立殘差向量與誤差向量之間的關(guān)系？由于r=b–Ax*=Ax–Ax*=A(x–x*)=Ae，故有e=A–1r，||e||||A–1||||r||另一方面，由b=Ax得||b||||A||||x||，因此其中，是近似解x*的相對(duì)誤差，而表示相對(duì)殘差上式表明近似解x*的相對(duì)誤差大約是相對(duì)殘差的||A||||A–1||倍。61/1092.條件數(shù)與病態(tài)方程組定義5：設(shè)A為n階非奇矩陣，稱(chēng)數(shù)||A||||A–1||為矩陣A的條件數(shù)，記為cond(A)。顯然，當(dāng)A的條件數(shù)cond(A)較小時(shí)，只要?dú)埐钕蛄亢苄?，近似解x*的相對(duì)誤差就很小，即近似解的精度較高.

但若cond(A)較大時(shí)，則即使相對(duì)殘差較小，近似解x*的相對(duì)誤差也可能很大，近似解的精度就很低。62/1092.條件數(shù)與病態(tài)方程組上例中線(xiàn)性方程組系數(shù)矩陣A的條件數(shù)：||A||=3，||A-1||

1000cond(A)3000

很大！這就解釋了為什么殘差向量小未必保證近似解的誤差小。定義6：當(dāng)cond(A)相對(duì)很大時(shí)，稱(chēng)方程組Ax=b為病態(tài)的，否則稱(chēng)為良態(tài)的。63/1093.病態(tài)方程組的判斷理論上雖然可以用矩陣的條件數(shù)來(lái)度量線(xiàn)性方程組求解問(wèn)題的病態(tài)程度.但是，條件數(shù)大到怎樣才算大，并沒(méi)有客觀(guān)的適用的標(biāo)準(zhǔn)，且當(dāng)方陣的階n較大時(shí)，A-1的計(jì)算也非易事，所以一般地說(shuō)，遇到下列幾種情況之一就應(yīng)考慮線(xiàn)性方程組可能是病態(tài)的：

(1)系數(shù)矩陣的某兩行(列)幾乎近似相關(guān)；

(2)系數(shù)矩陣的行列式的值很小；

(3)用主元消去法解線(xiàn)性方程組時(shí)出現(xiàn)小主元；

(4)近似解x*已使殘差向量r=b-Ax*的范數(shù)很小，但該近似解仍不符合問(wèn)題的要求．64/1094.解的精度估計(jì)與改善我們知道，影響計(jì)算精度的一個(gè)重要因素是舍入誤差的傳播．如前所述，按計(jì)算過(guò)程逐步地分析舍入誤差的傳播是極為困難又不實(shí)用的．一種常用的辦法是向后誤差分析方法，它將計(jì)算過(guò)程中舍入誤差的影響歸結(jié)為原始數(shù)據(jù)的小擾動(dòng)，也就是說(shuō)，將計(jì)算所得的近似解x*看成是某個(gè)帶有小擾動(dòng)的方程組(A+A)x=b+b的準(zhǔn)確解．65/109定理：設(shè)x*是線(xiàn)性方程組(A+A)x=b+b的解，x是線(xiàn)性方程組Ax=b的解，且A的擾動(dòng)A滿(mǎn)足||A-1||||A||<1，則有上式表明了近似解x*的相對(duì)誤差與A及b的相對(duì)誤差之間的關(guān)系，其中A的條件數(shù)起著重要的作用．66/109估計(jì)式只是理論性的，實(shí)際應(yīng)用有很大的困難．因此，—般采用下述方法估計(jì)近似解的精確程度．對(duì)于線(xiàn)性方程組Ax=b，任取一個(gè)非零向量y，計(jì)算出d=Ay；然后采用與解Ax=b同一算法解線(xiàn)性方程組Ay=d，得到它的近似解y*，y*的相對(duì)誤差為由于求解Ay=d與求解Ax=b的算法相同，且系數(shù)矩陣同為A，所以可以認(rèn)為Ax=b的近似解x*的相對(duì)誤差大約也是當(dāng)然，對(duì)病態(tài)方程組來(lái)說(shuō)，這種估計(jì)方法是不大可靠的．67/109

病態(tài)方程組的求解問(wèn)題是計(jì)算方法研究的一個(gè)重要課題．這里介紹一下改善近似解的精度的方法：

1)采用雙倍字長(zhǎng)進(jìn)行運(yùn)算，特別是在做形如∑aibi的計(jì)算時(shí)，每一乘積項(xiàng)均取雙倍字長(zhǎng)的數(shù)，并與其他項(xiàng)累加后再舍入成單字長(zhǎng)，這樣能使結(jié)果的精度提高，同時(shí)也比全都用雙倍字長(zhǎng)進(jìn)行運(yùn)算節(jié)省內(nèi)存和機(jī)時(shí)．

68/1092)采用迭代改善的辦法．設(shè)已得到Ax＝b的近似解x(0)，采用雙倍字長(zhǎng)計(jì)算Ax(0)后再舍入成單字長(zhǎng)，并計(jì)算殘差向量r(0)=b–Ax(0)，記x=x–x(0)，則x應(yīng)滿(mǎn)足方程組Ax＝r(0)．（Ax＝A(

x–x(0))＝Ax–Ax(0)=b–Ax(0)=r(0)）采用原單字長(zhǎng)運(yùn)算的算法求出x的近似解x(0)(在A(yíng)=LR或A=QR分解已完成下，只需作順代和回代求解)

從而得到Ax＝b的改善了的近似解x(1)=x(0)+x(0)，若||x(0)||<，則x(1)就是滿(mǎn)足精度（）要求的解。否則再重復(fù)上述做法。69/109§2.4線(xiàn)性方程組迭代解法1.基本思想若方程組Ax=b可寫(xiě)成等價(jià)形式x=Bx+g，(2.4-1)則給定一個(gè)初始向量x(0)，可以得到迭代公式x(k+1)=Bx(k)+g，k=0，1，…(2.4-2)若上式確定的向量序列{x(k)}收斂于x，則x顯然是(2.4-1)的解，從而為方程組Ax=b的解。形如(2.4-2)的逐次逼近的方法稱(chēng)為簡(jiǎn)單迭代法，B稱(chēng)為該迭代法的迭代矩陣。70/1092.迭代法的收斂性

定理1：對(duì)任意初始向量x(0)及常向量g，迭代法(2.4-2)收斂的充分必要條件是迭代矩陣B的譜半徑(B)<1。這一結(jié)論在理論上是頗為重要的，但實(shí)際用起來(lái)不甚方便。當(dāng)n較大時(shí)，譜半徑的計(jì)算并非易事，但由于(B)||B||，所以只要B的某種范數(shù)小于1，迭代法(2.4-2)必收斂：

定理2：若迭代矩陣Ｂ的某種范數(shù)||B||<1，則(2.4-2)確定的迭代法對(duì)任意初值x(0)均收斂于方程組x=Bx+g的唯一解x。71/109

定理3：若||B||=q<1，則(2.4-2)確定的向量序列滿(mǎn)足：(2.4-3)(2.4-4)其中x是方程組(2.4-1)的準(zhǔn)確解。證明：(1)因?yàn)閤(k)–x=B(x(k–1)–x)=B(x(k–1)–x(k)+x(k)–x)，所以||x(k)–x||=||B[(x(k–1)–x(k))+(x(k)–x)]||

||B||||(x(k–1)–x(k))+(x(k)–x)||

q(||x(k–1)–x(k)||+||x(k)–x||)，從而得(2.4-3)。72/109(2)因?yàn)閤(k)–x(k–1)=B(x(k–1)+g)–B(x(k–2)+g)=B(x(k–1)–x(k–2))=…=B

k–1(x(1)–x(0))，所以有||x(k)–x(k–1)||=||B

k–1(x(1)–x(0))||

||B

k–1||||x(1)–x(0)||

||B||

k–1||x(1)–x(0)||代入(2.4-3)即得(2.4-4)(2.4-3)與(2.4-4)給出||x(k)–x||的估計(jì)分別稱(chēng)為后驗(yàn)估計(jì)與先驗(yàn)估計(jì)。73/109先驗(yàn)估計(jì)可以在求出x(1)后就可以估計(jì)出在精度要求為下大約需要的迭代次數(shù)k。從證明過(guò)程看出，先驗(yàn)估計(jì)是由后驗(yàn)估計(jì)放大得到，所以較粗。從(2.4-4)可以知道，q=||B||越小，{x(k)}收斂越快。x(1)–x(0)=Bx(0)+g–x(0)=g–(I–B)x(0)的范數(shù)越小，或x(0)滿(mǎn)足方程組Ax=b的程度越好，則x(k)與x就越接近。74/109后驗(yàn)估計(jì)用計(jì)算過(guò)程中前后相繼兩個(gè)迭代向量之差的范數(shù)來(lái)估計(jì)后一個(gè)近似值的誤差，較為可靠，并給出迭代終止準(zhǔn)則，例如在q

0.5的條件下，有(1)絕對(duì)誤差準(zhǔn)則當(dāng)||x(k)–x(k–1)||<時(shí)終止(2)相對(duì)誤差準(zhǔn)則當(dāng)||x(k)–x(k–1)||/||x(k)||<時(shí)終止注意當(dāng)q=||B||接近于1時(shí)，此準(zhǔn)則可能失效。75/1093.簡(jiǎn)單迭代法、Jacobi迭代與Seidel迭代若將方陣A分裂為A=M

–

N其中M可逆，則線(xiàn)性方程組Ax=b可寫(xiě)成等價(jià)形式Mx=Nx+b，或x=Bx+g其中B=M–1N，g=M–1b稱(chēng)x(k+1)=Bx(k)+g，k=0，1，…

為簡(jiǎn)單迭代法輸入矩陣A，向量b，初值x0，eps計(jì)算B，g計(jì)算x1=B*x0+gwhilenorm(x1-x0,inf)>epsx0=x1k=k+1x1=B*x0+g輸出k，x176/109若有迭代公式（分量式）(2.4-2')由于在計(jì)算新分量xi(k+1)時(shí)，x1(k+1)，x2(k+1)，…，xi-1(k+1)都已經(jīng)算出，故可將上式改為(2.4-7)稱(chēng)(2.4-7)為Seidel迭代法（迭代公式）。77/109將式中x1(k+1)，x2(k+1)，…，xn-1(k+1)項(xiàng)移到左邊：即有78/109Seidel迭代法的矩陣形式為：x(k+1)=BSx(k)+h其中，BS稱(chēng)為Seidel迭代的迭代矩陣。79/109例1設(shè)線(xiàn)性方程組問(wèn)p，q滿(mǎn)足怎樣條件時(shí)，對(duì)任意初始向量，簡(jiǎn)單迭代和Seidel迭代收斂？解：因?yàn)楹?jiǎn)單迭代法的迭代矩陣，而B(niǎo)的特征多項(xiàng)式：(p–)2=q2，特征值為

=p

iq，所以當(dāng)|

|=p2+q2<1時(shí)，簡(jiǎn)單迭代法對(duì)任意初始向量x(0)收斂。80/109Seidel迭代公式是即其迭代矩陣是BS的特征多項(xiàng)式是從而根據(jù)E.I.Jury準(zhǔn)則，若p2<1，q2<(1+p)2，則有(BS)<1，即有Seidel迭代對(duì)任意初始向量x(0)收斂。81/109例2用簡(jiǎn)單迭代法解線(xiàn)性方程組問(wèn)用后驗(yàn)估計(jì)與先驗(yàn)估計(jì)迭代次數(shù)k要多大才使近似解x(k)的誤差小于10-3？取x(0)=0.輸入矩陣A，向量b，初值x0，eps計(jì)算B，g計(jì)算x1=B*x0+gwhilenorm(x1-x0,inf)>epsx0=x1k=k+1x1=B*x0+g輸出k，x182/109解：在Matlab中B=[0.00000.1000-0.20000.0000;0.09090.00000.0909-0.2727;-0.20000.10000.00000.1000;0.0000-0.37500.12500.0000];g=[0.6000;2.2727;-1.1000;1.8750];norm(B,inf);x0=[0000]';x1=B*x0+g;k=1;whilenorm(x1-x0,inf)/norm(x1,inf)>0.001x0=x1;k=k+1x1=B*x0+gend對(duì)于先驗(yàn)估計(jì)，經(jīng)計(jì)算需k=1383/109在簡(jiǎn)單迭代法中若將系數(shù)矩陣分裂為A=D+L+R，其中D為A的對(duì)角線(xiàn)部分，L為A的嚴(yán)格下三角部分，R為A的嚴(yán)格上三角部分。則可將方程組化為等價(jià)方程組：(D+L+R)x=bDx=(–L–R)x+b

x=D-1(–L–R)x+D-1b，從而得到迭代公式x(k+1)=D-1(–L–R)x(k)+D-1b，稱(chēng)為Jacobi迭代法，其迭代矩陣為：BJ=D-1(–L–R)。其分量式為84/109Jacobi迭代算法流程圖：其中：BJ=D-1(–L–R)，g=D-1b輸入矩陣A，向量b，初值x0，eps計(jì)算D，L，R計(jì)算B，g計(jì)算x1=B*x0+gwhilenorm(x1-x0,inf)>epsx0=x1k=k+1x1=B*x0+g輸出k，x185/109例3-1用Jacobi迭代法解線(xiàn)性方程組取x(0)=0解：這里，，86/109Jacobi迭代法：87/109在Matlab中：A=[5-1-1-1;

-110-1-1;-1-15-1;-1-1-110];b=[-4;12;8;34];D0=diag(A);D=diag(D0);R=triu(A,1);L=tril(A,-1)BJ=inv(D)*(-L-R);g=inv(D)*bnorm(BJ,inf);x0=[0000]';x1=BJ*x0+g;k=1;whilenorm(x1-x0,inf)>0.001x0=x1;k=k+1x1=BJ*x0+gend88/109Jacobi迭代：x(k+1)=D-1(–L–R)x(k)+D-1b

x(k+1)=–D-1L

x(k)–D-1R

x(k)+D-1bSeidel迭代：x(k+1)=–D-1L

x(k+1)–D-1R

x(k)+D-1bDx(k+1)=–L

x(k+1)–R

x(k)+bDx(k+1)+L

x(k+1)=

–R

x(k)+bx(k+1)=

–(D+L)-1R

x(k)+(D+L)-1b即Seidel迭代的迭代矩陣：BS=-(D+L)-1R其分量式為89/109Seidel迭代的迭代矩陣：BS=-(D+L)-1R，g=(D+L)-1b算法：輸入矩陣A，向量b，初值x0，eps計(jì)算D，L，R計(jì)算B，g計(jì)算x1=B*x0+gwhilenorm(x1-x0,inf)>epsx0=x1k=k+1x1=B*x0+g輸出k，x190/109例3-2用Seidel迭代法解線(xiàn)性方程組取x(0)=0解：這里，，91/109Seidel迭代法：BS=-(D+L)-1R92/109在Matlab中：A=[5-1-1-1;

-110-1-1;-1-15-1;-1-1-110];b=[-4;12;8;34];D0=diag(A);D=diag(D0);R=triu(A,1);L=tril(A,-1);BS=-inv(D+L)*R;g=inv(D+L)*b;norm(BS,inf)x0=[0000]';x1=BS*x0+g;k=1;whilenorm(x1-x0,inf)>0.001x0=x1;k=k+1x1=BS*x0+gend93/1094.對(duì)角占優(yōu)方程組定義：設(shè)A=(aij)nn，若，i=1，2，…，n，則稱(chēng)A按行嚴(yán)格對(duì)角占優(yōu)

注：若A按行嚴(yán)格對(duì)角占優(yōu)，則A可逆94/109定理：設(shè)Ax=b，若A按行嚴(yán)格對(duì)角占優(yōu)，則對(duì)于任意的x(0)，J-迭代、S迭代均收斂。證明：(1)

，所以J-迭代收斂

95/109定理：設(shè)Ax=b，若A按行嚴(yán)格對(duì)角占優(yōu)，則對(duì)于任意的x(0)，J-迭代、S迭代均收斂。證明：(2)BS=-(D+L)-1R——欲證(BS)<1因?yàn)椋篒–BS=I+(D+L)-1R=(D+L)-1[(D+L)+R]所以：0=|I–BS|=|(D+L)-1||

(D+L)+R||(D+L)+R|=0令：有|G|=0.若|

|1，則

，i=1，2，...，n

G嚴(yán)格對(duì)角占優(yōu)

G可逆矛盾96/109注：若將方程組經(jīng)初等變換化為同解方程組，使其系數(shù)矩陣嚴(yán)格對(duì)角占優(yōu)，則可使迭代法收斂——p216例定理：若A正定，則S-迭代收斂若2D–A也正定，則J-迭代收斂，否則J-迭代不收斂97/109§2.5逐次超松弛迭代法和塊迭代法1.逐次超松弛迭代法使用迭代法的困難是計(jì)算量難以估計(jì)，有些方程組的迭代格式雖然收斂，但收斂速度慢而使計(jì)算量變得很大。

逐次超松弛迭代法(SuccessiveOverRelaxationMethod,簡(jiǎn)寫(xiě)為SOR)可以看作帶參數(shù)ω的高斯-塞德?tīng)柕?，是G-S方法的一種修正或加速.是求解大型稀疏矩陣方程組的有效方法之一.

這種方法將前一步的結(jié)果xi