解線性方程組的直接法

解線性方程組的直接法

自然科學和工程計算中的很多問題的解決常常歸結為求解線性方程組。如三次樣條插值函數(shù)問題、用最小二乘遠離確定擬合曲線、求解微分方程的數(shù)值解等，最終都要轉化為求解線性方程組。求解線性方程組可采用：5.1引言與預備知識直接法——假定計算過程沒有舍入誤差的情況下，經過有限步算術運算后能求得線性方程組精確解的方法。經過有限步運算就能求得精確解的方法，但實際計算中由于舍入誤差的影響，這類方法也只能求得近似解；例如：高斯消去法、三角分解法等。迭代法——構造適當?shù)南蛄啃蛄校媚撤N極限過程去逐步逼近精確解。例如:雅可比迭代法、高斯-賽德爾迭代法等。第2頁,共180頁，2024年2月25日，星期天一般地：對低階方程組用直接法，對高階方程組和大型稀疏方程組用迭代法求解。預備知識:1.向量與矩陣：第3頁,共180頁，2024年2月25日，星期天矩陣的基本運算(4)轉置矩陣(5)定義:

對于階矩陣，如果有一個階矩陣,

則說矩陣是可逆的，并把矩陣稱為的逆矩陣.使得第4頁,共180頁，2024年2月25日，星期天行列式性質(6)行列式按任一行展開其中Ai,j為ai,j的代數(shù)余子式，Ai,j=(-1)i+jMi,j,Mi,j為ai,j的余子式.第5頁,共180頁，2024年2月25日，星期天2.特殊矩陣(1)對角矩陣：如果ij時，ai,j=0.(2)三對角矩陣：如果|i-j|>1時，ai,j=0.形如(3)上三角矩陣：如果i>j時，ai,j=0.形如（類似有下三角矩陣）第6頁,共180頁，2024年2月25日，星期天(4)上海森伯格(Hessenberg)矩陣：如果i>j+1時，ai,j=0.形如(5)對稱矩陣：如果AT=A.(6)埃爾米特矩陣：設ACn×n,如果AH=A.(7)對稱正定矩陣：如果AT=A且對任意非零向量第7頁,共180頁，2024年2月25日，星期天(8)正交矩陣：如果A-1=AT.(9)酉矩陣：設ACn×n,如果A-1=AH.(10)初等置換陣：由單位矩陣I交換第i行與第j行（或交換第i列與第j列）得到的矩陣。(11)置換陣：由初等置換陣的乘積得到的矩陣。定理1：設ARn×n,則下述命題等價：1.對任何b

Rn,方程組AX=b有唯一解。2.齊次方程組AX=0只有唯一解X=0。3.det(A)≠0.4.A-1存在.5.A的秩rank(A)=n.第8頁,共180頁，2024年2月25日，星期天定理2：設ARn×n為對稱正定陣,則：1.A為非奇異矩陣，則A-1亦是對稱正定陣。2.記Ak為A的順序主子式，則Ak亦是對稱正定陣。3.A的特征值

i>0(i=1,2,…,n).4.A的順序主子式都大于零，即det(Ak)>0(k=1,2,…,n).定理3：設ARn×n為對稱陣,如果det(Ak)>0(k=1,2,…,n)，或A的特征值

i>0(i=1,2,…,n),則A為對稱正定陣。第9頁,共180頁，2024年2月25日，星期天定理4(Jordan標準型）：設A為n階矩陣,則存在一個非奇異矩陣P使得其中(1)當A的若當標準型中所有若當塊Ji均為一階時，此標準型若當(Jordan)塊變成對角陣。(2)若A的特征值各不相同，則其若當標準型必為對角陣diag(

1,2,…,n).第10頁,共180頁，2024年2月25日，星期天設線性方程組為或寫成矩陣形式或簡單地記為:5.2高斯消元法第11頁,共180頁，2024年2月25日，星期天第12頁,共180頁，2024年2月25日，星期天先看幾類特殊三角形方程組的解

★對角形方程組若aii≠0,則

xi=bi/aii,i=1,2,…,n。運算量O(n)。第13頁,共180頁，2024年2月25日，星期天★下三角方程組（返回LU）(5.19)第14頁,共180頁，2024年2月25日，星期天運算量：因為計算xi需要i

次乘法和除法，所以★上三角方程組（返回Gauss）返回LU(5.20)第15頁,共180頁，2024年2月25日，星期天運算量：因為計算xi

需要n-i

次乘法和1次除法，即

n–i+1次乘除，所以第16頁,共180頁，2024年2月25日，星期天

消元法的基本思想就是通過對方程組作初等變換，把一般形式的線性方程組化為等價的易于求解的三角方程組。★高斯消元法（GaussianElimination）高斯消元法是一種古老的求解線性方程組的方法，它就是通過一系列的初等變換（消元），把線性方程組（5.1）化為等價的上三角方程組（5.5），然后通過回代方法求出原方程組的解。順序高斯消去法消元過程:依從左到右、自上而下的次序將主對角元下方的元素化為零。1不作行交換。2用不等于零的數(shù)乘某行，加至另一行。第17頁,共180頁，2024年2月25日，星期天第18頁,共180頁，2024年2月25日，星期天第19頁,共180頁，2024年2月25日，星期天高斯消去法的主要思路：將系數(shù)矩陣A化為上三角矩陣，然后回代求解。考慮n階線性方程組：矩陣形式=下面我們討論一般線性方程組的高斯消去法第20頁,共180頁，2024年2月25日，星期天對于第21頁,共180頁，2024年2月25日，星期天重復上述過程，可以得到與（5.1）等價的上三角方程組：第22頁,共180頁，2024年2月25日，星期天或寫成矩陣形式A(n-1)x=b(n-1),

其中，第23頁,共180頁，2024年2月25日，星期天約定由（5.8）逐個回代，可得（5.1）的解第24頁,共180頁，2024年2月25日，星期天此即Gauss消去法.

綜述:若A非奇,總可通過帶有行交換或不帶有行交換的消元過程,將A化成非奇三角矩陣A(n-1).因此,回代求解過程也可進行到底.但實際中,A是否非奇異難以判定.算法應考慮:消元過程某一步找不到非零元,于是計算中斷.消元可進行到底,但ann(n-1)=

0,回代求解過程無法進行.故算法設計要考慮以上情況,給出計算中斷的信息.第25頁,共180頁，2024年2月25日，星期天

算法組織將增廣矩陣(A,b)放入一個二維數(shù)組.消去每一步算出的aij(k)放入aij的位置,bi(k)放入bi,mik放在aik上.消去完畢后,數(shù)組各位置上的數(shù)據如下示:回代過程的計算沒有數(shù)據組織問題.Ab第26頁,共180頁，2024年2月25日，星期天算法Guass消去法輸入n,aij,bi,(i,j=1,..,n)輸出解x1,x2

…,xn步驟1.Fork=1,2,…,n-1(消元)

1.1ifakk=0*Then輸出”不能消元”stop1.2fori=k+1,…,n1.2.1aik/akk

aik

1.2.2Forj=k+1,…,n1.2.2.1aij-

aik

akj

aij(新)

1.2.3bi-aik

bk

bi(新)說明:系數(shù)矩陣存放于數(shù)組A中,右端向量放于數(shù)組b中*常用|akk|≤

第27頁,共180頁，2024年2月25日，星期天步驟2.bn/ann

xn(回代)步驟3.Fork=n-1,…,13.1bks3.2Forj=k+1,…,n

3.2.1s-akj

xjs3.3s/akk

xk步驟4.End第28頁,共180頁，2024年2月25日，星期天

★高斯消元法運算量：消元過程：乘：

除：回代過程：運算量：第29頁,共180頁，2024年2月25日，星期天★高斯消元法的局限性：在高斯消元過程中，我們假定了對角元素由于每次消元時是按未知量的自然順序進行的，而順序消元不改變A的主子式的值，因此高斯消元法可行的充分必要條件為A的各階主子式不為零。實際上只要detA≠0，方程組就有解。2.

即使高斯消元法可行，如果很小，運算中用它作分母會導致其它元素數(shù)量計的嚴重增長和舍入誤差的擴散。第30頁,共180頁，2024年2月25日，星期天順序高斯消去法的使用條件

使用條件之一

定理線性方程組系數(shù)矩陣A的順序主子矩陣Ak

(k=1,2,…,n)非奇異，則順序高斯消去法能實現(xiàn)方程組的求解。即方程組能用順序高斯消去法求解的充要條件是系數(shù)行列式的順序主子式非零。在原方程組的系數(shù)矩陣中如何反映出這個條件呢?我們知道高斯消去法能按順序進行到底的充要條件應是：第31頁,共180頁，2024年2月25日，星期天A的k階順序主子矩陣Ak的行列式推論：如果A的順序主子矩陣的行列式detAk≠0

(k=1,2,…,n-1)，則第32頁,共180頁，2024年2月25日，星期天

使用條件之二

n階矩陣A為嚴格對角占優(yōu)矩陣是指其每個主對角元的絕對值大于同一行其他元素絕對值之和，即一階嚴格對角占優(yōu)矩陣指一個非零數(shù)。

定理

方程組系數(shù)矩陣A為嚴格對角占優(yōu)矩陣則可實現(xiàn)用順序高斯消去法求解。第33頁,共180頁，2024年2月25日，星期天Gauss主元素消元法例用Gauss消去法求解線性方程組解:[方法1]方程組的準確解為x*=(9998/9999,10000/9999)T.小主元a11回代求解x2=1.0000,x1=0.0000.計算結果相當糟糕.

原因:求行乘數(shù)時用了”小主元”0.0001作除數(shù).第34頁,共180頁，2024年2月25日，星期天解:[方法2]若上題在求解時將1,2行交換,即回代后得到

x=(1.0000,1.0000)T與準確解x*=(9998/9999,10000/9999)T相差不大.主元a11第35頁,共180頁，2024年2月25日，星期天再如：解方程組

精確解為：x1=1/3，x2=2/3。計算取5位有效數(shù)字。解(1)順序消元：

(2)交換方程的順序：經高斯消元：法（2）較好第36頁,共180頁，2024年2月25日，星期天可見，由于計算機和方程組本身的限制，在Gauss消元過程中會出現(xiàn)零主元和小主元，而使Gauss消去法無法進行或出現(xiàn)較大舍入誤差，為避免此類現(xiàn)象，可以通過行交換來選取絕對值較大的主元素.

前例方法2中所用的是在Gauss消去法中加入選主元過程,利用換行,避免”小主元”作除數(shù),該方法稱為Gauss列主元消去法。第37頁,共180頁，2024年2月25日，星期天★Gauss列主元消元法高斯消元過程中的稱為主元素。如果在一列中選取按模最大的元素作為主元素，將其調換到主干方程位置再做消元，則稱為列主元消元法。第一步,

先在第1列選出絕對值最大的元素,

交換第1行和第m行的所有元素，再做消元操作。第二步，對每個k=1，2，…，n做消元前，選出中絕對值最大的元素，對k行和m行交換后，再作消元。第38頁,共180頁，2024年2月25日，星期天列以下(含主對角元)元素ai,k(k-1)(k≤i≤n)中選絕對值最大者即:Max(|ai,k(k-1)|)(k≤i≤n)并通過行交換(A

(k-1),b

(k-1))具體地：在Gauss消去法的第k步(k=1,…,n-1)消元時,先在A(k-1)的第k的k行與第ik行對應元素,使位于主對角線上仍記為akk(k-1),稱之為第k步消元的列主元,然后再進行消元計算.因為|mi,k|≤1,列主元消去法能避免”小主元”作除數(shù),誤差分析與數(shù)值試驗表明：（Gauss）列主元消去法”基本穩(wěn)定”.第39頁,共180頁，2024年2月25日，星期天列主元消去法中方程對換問題引入向量P=(p1,p2,…,pn)t=(1,2,…,n)t即i

pi記錄方程(5.1)的行號.當?shù)趉步消元若需第k個方程與第ik個方程對換,則只需交換P中當消元過程完后,向量P中的pi(i=1,…,n)指示初始方程(5.1)在作的i步消元時被選為主元的方程序號.這種序號交換避免了方程間的對換,n較大時,有其優(yōu)越性.第40頁,共180頁，2024年2月25日，星期天例

用列主元高斯消去法求解方程組(用三位有效數(shù)字計算)解選主元選主元第41頁,共180頁，2024年2月25日，星期天消元過程完成，得到上三角形方程組再作回代可求得第42頁,共180頁，2024年2月25日，星期天開始輸出無解信息消元換行停機回代求解Gauss列主元消去法的算法設計流程圖第43頁,共180頁，2024年2月25日，星期天3.2Gauss列主元消去法步驟1.Fori=1,2,…,n(消元)1.1i

pi

(記錄方程(5.1)的行號)步驟2.Fork=1,2,…,n-12.1求最小的rk,使得

2.2ifThen輸出”方程奇異”stop2.3ifk≠rThenpk

pr(記錄主行所在方程的行號)2.4Fori=k+1,…,n2.4.12.4.2Forj=k+1,…,n2.4.2.12.4.3說明:系數(shù)矩陣存放于數(shù)組A中,右端向量放于數(shù)組b中.算法存儲,I/O同Guass消去法

第44頁,共180頁，2024年2月25日，星期天步驟3.(回代)步驟4.Fori=n-1,…,13.13.2Forj=k+1,…,n

3.2.13.3步驟5.End第45頁,共180頁，2024年2月25日，星期天

補充:若不引入向量P,只需將算法中的1.1步改為:求最小的rk使得if

Then輸出”方程奇異”stopIf

k≠rThen

bk

br;akj

arj(j=k,…,n)說明:列主元算法在實際中應用較廣“全面選主元”算法第46頁,共180頁，2024年2月25日，星期天算法3.2’（帶方程對換的）列主元Guass消去法輸入n,aij,bi,(i,j=1,..,n)輸出

解x1,x2

…,xn步驟1.

Fork=1,2,…,n-1(消元)

1.1a.

求最小的rk使得

b.

ifThen輸出”方程奇異”stop

c.

ifk≠rThen

bk

br;akj

arj(j=k,…,n)1.2fori=k+1,…,n1.2.1aik/akk

aik

1.2.2Forj=k+1,…,n1.2.2.1aij-

aik

akj

aij(新)

1.2.3bi-aik

bk

bi(新)|akk|≤

第47頁,共180頁，2024年2月25日，星期天步驟2.bn/ann

xn(回代)步驟3.Fork=n-1,…,13.1bks3.2Forj=k+1,…,n

3.2.1s-akj

xjs3.3s/akk

xk步驟4.End第48頁,共180頁，2024年2月25日，星期天全主元高斯消去法：第k步消元時選A(k)中絕對值最大的元素為主元，即①先選取全主元：ifik

k

then交換第k行和第ik行；ifjk

k

then交換第k列和第jk列；③消元列交換改變了

xi

的順序，須記錄交換次序，解完后再換回來。全主元高斯消去法具有很好的穩(wěn)定性，但選全主元比較費時，故在實際計算中很少使用。第49頁,共180頁，2024年2月25日，星期天程序示例★列主元高斯消元法算法描述：

1.輸入數(shù)據：n，A，b2.分解過程：（1）對k列選主元，k=1，2，…，n-1

（2）交換第k行和第m行方程；（3）生成高斯消元過程的新矩陣和右端項第50頁,共180頁，2024年2月25日，星期天

3.回代過程：對于i=n,n-1,…,2,1，解上三角方程組

4.輸出

第51頁,共180頁，2024年2月25日，星期天程序源碼主要部分：輸入Ax=b的A和b

printf("Nowinputthematrixa(i,j),i,j=0,1,...,%d:\n",n-1);for(i=0;i<n;i++){ for(j=0;j<n;j++){scanf("%lf",&a[i][j]); printf(“a[%d][%d]=%f”,i,j,a[i][j]);//輸出交換前A的元素

}}printf("Nowinputthematrixb(i),i=0,...,%d:\n",n-1);for(i=0;i<n;i++){scanf("%lf",&b[i]); printf("b[%d]=%f",i,b[i]);}第52頁,共180頁，2024年2月25日，星期天找主元素并進行交換兩行for(i=0;i<n-1;i++)//findmainelements{ for(j=i+1,mi=i,mx=fabs(a[i][i]);j<n;j++) if(fabs(a[j][i])>mx) {mi=j; mx=fabs(a[j][i]);} if(i<mi) {tem=b[i];b[i]=b[mi];b[mi]=tem; for(j=i;j<n;j++) {tem=a[i][j]; a[i][j]=a[mi][j]; a[mi][j]=tem;printf(“a[%d][%d]=%f”,mi,j,a[i][j]);//輸出交換后

A的元素

} }第53頁,共180頁，2024年2月25日，星期天高斯消元后生成新的矩陣元素for(j=i+1;j<n;j++){tem=-a[j][i]/a[i][i]; b[j]+=b[i]*tem; for(k=i;k<n;k++) a[j][k]+=a[i][k]*tem;}第54頁,共180頁，2024年2月25日，星期天回代解方程x[n-1]=b[n-1]/a[n-1][n-1];for(i=n-2;i>=0;i--){ x[i]=b[i]; for(j=i+1;j<n;j++) x[i]-=a[i][j]*x[j]; x[i]/=a[i][i];}第55頁,共180頁，2024年2月25日，星期天2.Gauss列主元消去法消元過程的矩陣描述由于Gauss列主元消去法每一步都要選取列主元,因此不可避免要進行行交換即表示不換行初等矩陣第56頁,共180頁，2024年2月25日，星期天因此,Gauss列主元消去法的消元過程為:顯然上三角陣仍然為單位下三角矩陣第57頁,共180頁，2024年2月25日，星期天初等矩陣的乘積,稱為排列陣則推廣到一般情形令仍然為單位下三角矩陣則單位下三角陣與上三角陣的乘積第58頁,共180頁，2024年2月25日，星期天綜合以上討論,有定理.(列主元素的三角分解定理):

第59頁,共180頁，2024年2月25日，星期天Gauss消去法是將系數(shù)矩陣化為上三角矩陣，再進行回代求解；Jordan法是將系數(shù)矩陣化為對角矩陣，從而無須再進行回代求解的一種直接解法。這種解法基本過程與Gauss法相同，只是在消去過程中，不僅將主對角線以下的元素消成0，而且同時將主對角線以上的元素也消成0。高斯-若當(Gauss-Jordan)消元法第60頁,共180頁，2024年2月25日，星期天第4行分別乘分別加到第3，2，1行上以n=4為例說明高斯-若當消元過程。首先高斯消元：第61頁,共180頁，2024年2月25日，星期天第62頁,共180頁，2024年2月25日，星期天第63頁,共180頁，2024年2月25日，星期天第64頁,共180頁，2024年2月25日，星期天第65頁,共180頁，2024年2月25日，星期天歸一方法的例子第66頁,共180頁，2024年2月25日，星期天列選主高斯-約當(Jordan)消去法。第67頁,共180頁，2024年2月25日，星期天第68頁,共180頁，2024年2月25日，星期天高斯-約當(Jordan)消去法求逆。第69頁,共180頁，2024年2月25日，星期天第70頁,共180頁，2024年2月25日，星期天5.3

矩陣三角分解法5.3.1

直接三角分解法

5.3.2

平方根法5.3.3

追趕法第71頁,共180頁，2024年2月25日，星期天

定義將矩陣A分解成一個下三角陣L和一個上三角陣U的乘積，即A=LU，稱為A的三角分解或LU分解。第72頁,共180頁，2024年2月25日，星期天例如A=這里有A的兩種不同的三角分解，類似可舉出很多，一般，若A=LU是一個三角分解，任取與A同階的非奇異對角矩陣D，則A=(LD)(D-1U)=L1U1也是A的三角分解。第73頁,共180頁，2024年2月25日，星期天。杜利特爾分解(Doolittle)常用的兩種三角分解克洛特Crout分解有的叫庫朗(Courant)分解第74頁,共180頁，2024年2月25日，星期天5.3.1

直接三角分解法

用矩陣相乘來解釋高斯消元過程，取n=4。記(5.1)化為方程組(5.6)的過程，即等價于L1A=A(1),L1b=b(1),即第75頁,共180頁，2024年2月25日，星期天記(5.6)化為方程組(5.7)的過程等價于L2A(1)

=A(2),L2b(1)=b(2),即第76頁,共180頁，2024年2月25日，星期天記第77頁,共180頁，2024年2月25日，星期天因此其中，=LU第78頁,共180頁，2024年2月25日，星期天第79頁,共180頁，2024年2月25日，星期天一般地，高斯消元完成后,有：故從而U=A(n-1)第80頁,共180頁，2024年2月25日，星期天即且順序主元第81頁,共180頁，2024年2月25日，星期天解:由上述分析不難得到

問題：矩陣A存在LU分解（即Gauss消去法可以執(zhí)行）的條件是什么？第82頁,共180頁，2024年2月25日，星期天Gauss消去法可以執(zhí)行定理

證明（略）.推論(P172)：如果A的順序主子矩陣的行列式detAk≠0

(k=1,2,…,n-1)，則第83頁,共180頁，2024年2月25日，星期天由上述對系數(shù)矩陣A左乘一系列的三角初等矩陣知，A可以分解為一個單位下三角矩陣L和一個上三角矩陣U。這個過程派生出解線性方程組的直接分解法。一旦實現(xiàn)A=LU，則Ax=b可以化為LUx=b。令Ux=y，則Ly=b。由Ly=b（即(5.4)）解出y；再由Ux=y（即(5.5)）解出x。如果A的各階主子式不為零，則可以實現(xiàn)：杜利特爾(Doolittle)分解：如果L為單位下三角矩陣，U為上三角矩陣；克洛特Crout分解：如果L為下三角矩陣,U為單位上三角矩陣。第84頁,共180頁，2024年2月25日，星期天杜利特爾分解法設A

的各階主子式不為零，A分解為A=LU，即先計算U的元素，后計算L的元素：U的第1行：第85頁,共180頁，2024年2月25日，星期天L的第1列：第86頁,共180頁，2024年2月25日，星期天再計算U的第2行元素；計算L的第2列元素；……計算U的第k行元素：第87頁,共180頁，2024年2月25日，星期天計算L的第k列元素：第88頁,共180頁，2024年2月25日，星期天綜合以上分析,可以推導出：U的第一行L的第一列------(1)------(2)U的第r行------(3)(逐行算出)L的第r列

------(4)(逐列算出)上述(1)~(4)式所表示的分解過程稱為A的Doolittle分解LU分解的運算量：第89頁,共180頁，2024年2月25日，星期天LU分解的總運算量：由（5.17）知，計算U的第2行n-1個元素：（n-1)×1（1項乘積）第3行n-2個元素：（n-2)×2（2項乘積）

………

第n行1個元素：1×（n-1）（n-1項乘積）(n-1)求解總運算量：第90頁,共180頁，2024年2月25日，星期天舉例第91頁,共180頁，2024年2月25日，星期天用緊湊格式解：第92頁,共180頁，2024年2月25日，星期天再如：第93頁,共180頁，2024年2月25日，星期天緊湊格式分解A=LU總結：1.先確定U中第一行元素（即等于A中第一行元素）.2.再確定L中第一列元素：3.確定U中第r行元素:4.再確定L中第r列元素:第94頁,共180頁，2024年2月25日，星期天再如：作分解A=LU第95頁,共180頁，2024年2月25日，星期天再如：第96頁,共180頁，2024年2月25日，星期天對于線性方程組系數(shù)矩陣非奇異,經過Doolittle分解后Ax=L(Ux)=b可化為下面兩個三角形方程組消去回代L=第97頁,共180頁，2024年2月25日，星期天上述解線性方程組的方法稱為三角(LU)分解的Doolittle法.第98頁,共180頁，2024年2月25日，星期天Doolittle法的特點:緊湊(無中間過程);內積計算(精度高)例1：

用Doolittle法解方程組解:由Doolittle分解逐行算出U的元素逐列算出L的元素第99頁,共180頁，2024年2月25日，星期天逐行算出U的元素逐列算出L的元素第100頁,共180頁，2024年2月25日，星期天例2用多利特爾分解求解方程組：解

設A=LU，即

第101頁,共180頁，2024年2月25日，星期天解下三角方程組Ly=b，即解上三角方程組Ux=y，即第102頁,共180頁，2024年2月25日，星期天A,b,x,L,U,y需要單獨存儲,而從lij,uij的計算過程發(fā)現(xiàn):算法的數(shù)據組織:對于線性方程，

第103頁,共180頁，2024年2月25日，星期天Doolittle分解法的數(shù)據組織可以用以下緊湊格式表示:第104頁,共180頁，2024年2月25日，星期天UL第105頁,共180頁，2024年2月25日，星期天

Doolittle法解本節(jié)方程組(例1)的數(shù)據組織過程：解:第106頁,共180頁，2024年2月25日，星期天第107頁,共180頁，2024年2月25日，星期天所以x第108頁,共180頁，2024年2月25日，星期天例4用杜利特爾分解法求解方程組解先求增廣系數(shù)矩陣的杜利特爾分解，即第109頁,共180頁，2024年2月25日，星期天再一例：第110頁,共180頁，2024年2月25日，星期天第111頁,共180頁，2024年2月25日，星期天開始輸出錯誤信息停機1.Doolittle法解方程組流程圖程序實現(xiàn)計算U的第k行計算L的第k列輸出x由Ly=b計算yk輸出Uk,Lk,yk由Ux=y計算xk第112頁,共180頁，2024年2月25日，星期天2.Doolittle分解主要算法程序：輸入：n，A，b計算U的第k行元素和L的第k列元素計算U的第k行計算L的第k列3.程序源代碼：（Doolim.C，Dooli.C）回代求解第113頁,共180頁，2024年2月25日，星期天小結

Doolittle求解法的特點:緊湊，比高斯消去法精度高.2.實際計算中，由于舍入誤差的影響，不可避免地往往只能求得近似解.3.為保證LU分解能夠進行，要求A的各階順序主子式不為零，為取消此限制，可采用類似列主元消元法處理，程序稍復雜，或采用其它算法（如迭代法）。4.本實驗是后續(xù)改進算法(平方根法、Cholesky分解法)的基礎。對于大型方程組，一般采用迭代算法（下一章）.第114頁,共180頁，2024年2月25日，星期天問題:在Doolittle法(包括LU分解緊湊格式)中,會反復用到公式仍有可能為小主元做除數(shù)?為此,也應考慮在算法中加入選取列主元即緊湊格式的Doolittle列主元法.第115頁,共180頁，2024年2月25日，星期天列主元Doolittle法步驟:第一步:第116頁,共180頁，2024年2月25日，星期天例.用列主元Doolittle法解線性方程組解:第117頁,共180頁，2024年2月25日，星期天所以原方程組的解為第118頁,共180頁，2024年2月25日，星期天克洛特Crout分解*(如果L為下三角矩陣，U為單位上三角矩陣)

先計算L的第一列，再計算U的第一行，其余類此。類似于多利特爾分解法，得到：對k=1，2，…，n返回(5.26)第119頁,共180頁，2024年2月25日，星期天

實現(xiàn)A=LU，則Ax=b可以化為LUx=b。令Ux=y，則Ly=b。由Ly=b解出yi:

再由Ux=y解出xi:

第120頁,共180頁，2024年2月25日，星期天注*：為保證LU分解能夠進行，要求，既要求A的

所有順序主子式不為零，而線性方程組有解，只須|A|≠0

即可。為取消此限制，采用類似列主元消元法處理。克洛特Crout分解列主元直接分解算法*輸入：n，A，b計算L的第k列元素和U的第k行元素第121頁,共180頁，2024年2月25日，星期天第122頁,共180頁，2024年2月25日，星期天例4*

用克洛特分解求解方程組：解

設A=LU，即

第123頁,共180頁，2024年2月25日，星期天解下三角方程組Ly=b，即解(單位)上三角方程組Ux=y，即第124頁,共180頁，2024年2月25日，星期天★克洛特Crout分解算法描述第125頁,共180頁，2024年2月25日，星期天由Ly=b解出yi:

再由Ux=y解出xi:第126頁,共180頁，2024年2月25日，星期天程序源碼主要部分for(i=0;i<n;i++)u[i][i]=1;//U矩陣對角元素為1for(k=0;k<n;k++){for(i=k;i<n;i++)//計算L矩陣的第k列元素

{ l[i][k]=a[i][k]; for(j=0;j<=k-1;j++) l[i][k]-=(l[i][j]*u[j][k]); }}第127頁,共180頁，2024年2月25日，星期天for(j=k+1;j<n;j++)//計算矩陣U的第k行元素{u[k][j]=a[k][j]; for(i=0;i<=k-1;i++){u[k][j]-=(l[k][i]*u[i][j]);}u[k][j]/=l[k][k];}第128頁,共180頁，2024年2月25日，星期天for(i=0;i<n;i++)//由Ly=b計算y{ y[i]=b[i]; for(j=0;j<=i-1;j++) y[i]-=(l[i][j]*y[j]);y[i]/=l[i][i];}第129頁,共180頁，2024年2月25日，星期天for(i=n-1;i>=0;i--)//由Ux=y計算x{ x[i]=y[i]; for(j=i+1;j<n;j++) x[i]-=(u[i][j]*x[j]);}第130頁,共180頁，2024年2月25日，星期天5.3.2平方根法第131頁,共180頁，2024年2月25日，星期天5．A對稱正定，則A的對角元素aii>0，i=1,2,…,n。第132頁,共180頁，2024年2月25日，星期天第133頁,共180頁，2024年2月25日，星期天定理1若A為對稱矩陣，且A的各階順序主子式不為零，則A可唯一分解成第134頁,共180頁，2024年2月25日，星期天證明第135頁,共180頁，2024年2月25日，星期天引理

如果定理1中的A正定，則對角矩陣D中的元素非負。定理2(Cholesky分解)

若A為對稱正定矩陣，則存在下三角（喬列斯基）矩陣U，使A=UUT。（平方根法）證明：因為A對稱，由定理1知，其中L為單位下三角陣，D為對角陣。如果A正定，則由引理知，D中的元素非負。記

第136頁,共180頁，2024年2月25日，星期天下面我們用直接分解法來確定計算L元素的遞推公式：幾點說明：(1)L

為對角線元素全為正的下三角矩陣；(2)只有對稱正定矩陣(SPD)才存在Cholesky分解；(3)SPD矩陣的Cholesky分解存在且唯一將展開第137頁,共180頁，2024年2月25日，星期天-------------(1)-------------(2)-------------(3)第138頁,共180頁，2024年2月25日，星期天第139頁,共180頁，2024年2月25日，星期天緊湊格式喬列斯基(Cholesky)分解A=LLT總結：1.先確定L中第一列元素：2.確定L中對角線元素:4.確定L中第i行元素:(r=2,3,…,i-1)第140頁,共180頁，2024年2月25日，星期天第141頁,共180頁，2024年2月25日，星期天第142頁,共180頁，2024年2月25日，星期天緊湊格式解法：第143頁,共180頁，2024年2月25日，星期天第144頁,共180頁，2024年2月25日，星期天對稱正定線性方程組的解法線性方程組-------------(5)-------------(6)則線性方程組(10)可化為兩個三角形方程組-------------(7)-------------(8)第145頁,共180頁，2024年2月25日，星期天------(9)------(10)對稱正定方程組的平方根法第146頁,共180頁，2024年2月25日，星期天例1.用平方根法解對稱正定方程組解:第147頁,共180頁，2024年2月25日，星期天第148頁,共180頁，2024年2月25日，星期天即所以原方程組的解為第149頁,共180頁，2024年2月25日，星期天第150頁,共180頁，2024年2月25日，星期天第151頁,共180頁，2024年2月25日，星期天平方根法的數(shù)值穩(wěn)定性用平方根法求解對稱正定方程組時不需選取主元由可知因此平方根法是數(shù)值穩(wěn)定的事實上,對稱正定方程組也可以用順序Gauss消去法求解而不必加入選主元步驟第152頁,共180頁，2024年2月25日，星期天改進的平方根法5.325.33第153頁,共180頁，2024年2月25日，星期天例如：

用LDLT解方程組解第154頁,共180頁，2024年2月25日，星期天于是，Ax=b

LDLTx=b.解方程組有：第155頁,共180頁，2024年2月25日，星期天第156頁,共180頁，2024年2月25日，星期天第157頁,共180頁，2024年2月25日，星期天第158頁,共180頁，2024年2月