計算方法線性方程組的直接解法

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

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

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

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

–1時得完成第n

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

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

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

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

x1=-10.0，x2=1.01，此解無效；

2)

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

x1=10.0，x2=1.00.18/109(5)行標(biāo)度化當(dāng)線性方程組的系數(shù)矩陣的元素大小相差很大時，則可能引進(jìn)因丟失有效數(shù)字而產(chǎn)生的誤差，并且舍入誤差的影響嚴(yán)重，即使用Gauss主元消去法得到的解也不可靠．為避免這一問題，可將系數(shù)矩陣的行元素按比例增減以使元素的變化減小．如用每行元素的最大模除該行各元素，使它們的模都不大于1，這叫行標(biāo)度化，其目的是要找到真正的主元．19/109(5)行標(biāo)度化如用每行元素的最大模除該行各元素，使它們的模都不大于1，這叫行標(biāo)度化，其目的是要找到真正的主元．消元過程仍是對原方程組進(jìn)行的，只不過在Gauss列主元消去法的算法中，按列選主元ck時，應(yīng)修改為其中20/1095.算法設(shè)計輸入A，n，epsFor(k=1ton-1)選主元A(I0，k)——確定行號p=A(I0，k)若|p|<=epsText=1，退出循環(huán)若I0kT換行（I0k）消元若|A(n，n)|<=epsText=1F回代若ext=1T無解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è)計(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矩陣分解在解線性方程組中的作用1.原理若A=LR，則Ax=b

LRx=b

若其中L、R為三角形矩陣，則方程組易解定義：(1)稱為單位下三角陣(2)設(shè)L為單位下三角陣，R為上三角陣，若A=LR，則稱A可三角（LR）分解，并稱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個順序主子式全不等于0注：由Ax=b

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

>

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

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

>

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

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

>

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

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

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

>

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

j=k,k+1,...,n；k=1,2,...,n(2)主對角線下邊：當(dāng)行標(biāo)m>列標(biāo)k時，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ù)的對角陣。(2)設(shè)A=LDLT則有：33/109(3)為了避免重復(fù)計算，引進(jìn)中間量tik=likdk，將上式改寫為：34/109(4)Ax=b

LDLTx=b

Ly=b，Dz=y，LTx=z求解公式：上述方法稱為“改進(jìn)的平方根法”注：不能選主元作行交換——破壞對稱性，必是數(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.追趕法在實際問題中，經(jīng)常遇到以下形式的方程組這種方程組的系數(shù)矩陣稱為三對角矩陣。以下針對該方程組的特點提供一種簡便有效的算法——

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

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

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

–

liyi–1

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

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

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

x，y

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

x

≥0，且x=0

x=0(2)

ax

=

a

x

(3)

x+y

≤

x

+

y

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

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

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

，這里k

F3)三角不等式：α+

α

+

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

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

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

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

1

p

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

||x||a

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

x=0(2)(3)因為A正定，所以可逆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ù)49/109有了范數(shù)的概念，就可以討論向量序列的收斂性問題。定義2：設(shè)給定Cn中的向量序列{xk}，即x0，x1，…，xk，…其中若對任何i(i=1,2,…,n)都有則向量稱為向量序列{xk}的極限，或者說向量序列{xk}依坐標(biāo)收斂于向量x*，記為50/109定義2：設(shè)給定Cn中的向量序列{xk}，即x0，x1，…，xk，…其中若對任何i(i=1,2,…,n)都有則向量稱為向量序列{xk}的極限，或者說向量序列{xk}依坐標(biāo)收斂于向量x*，記為定理7：向量序列{xk}依坐標(biāo)收斂于x*的充要條件是如果一個向量序列{xk}與向量x*滿足上式，就說向量序列{xk}依范數(shù)收斂于x*，于是便得：向量序列依范數(shù)收斂與依坐標(biāo)收斂是等價的。51/1092.矩陣的范數(shù)下面我們給出矩陣范數(shù)的一般定義。為簡單起見，僅討論實數(shù)域的情形。定義3

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

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

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

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

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

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

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

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

設(shè)xRn及矩陣ARnn，||x||v為一種向量范數(shù)(如v=1，2或∞)，相應(yīng)地定義一個矩陣的非負(fù)函數(shù)可驗證||A||v滿足定義3，所以||A||v是Rnn上矩陣的一個范數(shù)，||A||v還滿足相容性條件5)，稱為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ù)分別稱為1-范數(shù)，2-范數(shù)和-范數(shù)又稱為列模、譜模和行模56/109注：A的范數(shù)||A||與A的特征值間的關(guān)系設(shè)為矩陣A的任一特征值，e為相應(yīng)的特征向量，則e=Ae，因為||||e||=||Ae||||A||||e||，所以||||A||從而得定理9：矩陣A的任一特征值的絕對值不超過A的范數(shù)||A||。定義6：矩陣A的諸特征值的最大絕對值稱為A的譜半徑，記為：由定理9可知：57/109例(p138)：設(shè)A=(aij)nn，||A||為其算子范數(shù)，則||A||<1

I–A可逆，且58/109§2.3直接法的誤差分析1.誤差向量與殘差向量設(shè)x為線性方程組Ax=b的準(zhǔn)確解，用數(shù)值算法得到的近似解為x*，則稱向量e=x*-x為誤差向量。顯然，||e||可以表示近似解x*的準(zhǔn)確程度：||e||越小，近似解x*越準(zhǔn)確。但由于x不知道，所以無法計算e。變通的辦法是考慮殘差向量：r=b–Ax*的范數(shù)的大小。但是，在某些情況下，盡管||r||很小，||e||也可能很大。59/109例線性方程組的準(zhǔn)確解是(1，1)T，若用某種方法得到計算解(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*的相對誤差，而表示相對殘差上式表明近似解x*的相對誤差大約是相對殘差的||A||||A–1||倍。61/1092.條件數(shù)與病態(tài)方程組定義5：設(shè)A為n階非奇矩陣，稱數(shù)||A||||A–1||為矩陣A的條件數(shù)，記為cond(A)。顯然，當(dāng)A的條件數(shù)cond(A)較小時，只要殘差向量很小，近似解x*的相對誤差就很小，即近似解的精度較高.

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

1000cond(A)3000

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

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

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

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

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

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

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

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

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

從而得到Ax＝b的改善了的近似解x(1)=x(0)+x(0)，若||x(0)||<，則x(1)就是滿足精度（）要求的解。否則再重復(fù)上述做法。69/109§2.4線性方程組迭代解法1.基本思想若方程組Ax=b可寫成等價形式x=Bx+g，(2.4-1)則給定一個初始向量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)的逐次逼近的方法稱為簡單迭代法，B稱為該迭代法的迭代矩陣。70/1092.迭代法的收斂性

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

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

定理3：若||B||=q<1，則(2.4-2)確定的向量序列滿足：(2.4-3)(2.4-4)其中x是方程組(2.4-1)的準(zhǔn)確解。證明：(1)因為x(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)因為x(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||的估計分別稱為后驗估計與先驗估計。73/109先驗估計可以在求出x(1)后就可以估計出在精度要求為下大約需要的迭代次數(shù)k。從證明過程看出，先驗估計是由后驗估計放大得到，所以較粗。從(2.4-4)可以知道，q=||B||越小，{x(k)}收斂越快。x(1)–x(0)=Bx(0)+g–x(0)=g–(I–B)x(0)的范數(shù)越小，或x(0)滿足方程組Ax=b的程度越好，則x(k)與x就越接近。74/109后驗估計用計算過程中前后相繼兩個迭代向量之差的范數(shù)來估計后一個近似值的誤差，較為可靠，并給出迭代終止準(zhǔn)則，例如在q

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

–

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

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

=p

iq，所以當(dāng)|

|=p2+q2<1時，簡單迭代法對任意初始向量x(0)收斂。80/109Seidel迭代公式是即其迭代矩陣是BS的特征多項式是從而根據(jù)E.I.Jury準(zhǔn)則，若p2<1，q2<(1+p)2，則有(BS)<1，即有Seidel迭代對任意初始向量x(0)收斂。81/109例2用簡單迭代法解線性方程組問用后驗估計與先驗估計迭代次數(shù)k要多大才使近似解x(k)的誤差小于10-3？取x(0)=0.輸入矩陣A，向量b，初值x0，eps計算B，g計算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對于先驗估計，經(jīng)計算需k=1383/109在簡單迭代法中若將系數(shù)矩陣分裂為A=D+L+R，其中D為A的對角線部分，L為A的嚴(yán)格下三角部分，R為A的嚴(yán)格上三角部分。則可將方程組化為等價方程組：(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，稱為Jacobi迭代法，其迭代矩陣為：BJ=D-1(–L–R)。其分量式為84/109Jacobi迭代算法流程圖：其中：BJ=D-1(–L–R)，g=D-1b輸入矩陣A，向量b，初值x0，eps計算D，L，R計算B，g計算x1=B*x0+gwhilenorm(x1-x0,inf)>epsx0=x1k=k+1x1=B*x0+g輸出k，x185/109例3-1用Jacobi迭代法解線性方程組取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計算D，L，R計算B，g計算x1=B*x0+gwhilenorm(x1-x0,inf)>epsx0=x1k=k+1x1=B*x0+g輸出k，x190/109例3-2用Seidel迭代法解線性方程組取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.對角占優(yōu)方程組定義：設(shè)A=(aij)nn，若，i=1，2，…，n，則稱A按行嚴(yán)格對角占優(yōu)

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

，所以J-迭代收斂

95/109定理：設(shè)Ax=b，若A按行嚴(yán)格對角占優(yōu)，則對于任意的x(0)，J-迭代、S迭代均收斂。證明：(2)BS=-(D+L)-1R——欲證(BS)<1因為：I–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)格對角占優(yōu)

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

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

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