第三章解線性方程組的直接法本章目標求解線性方程組：1、高斯消去法

高斯消去法選主元消去法約當消去法2、矩陣三角分解法直接分解法平方根法追趕法第三章解線性方程組的直接法本章目標求解線性方程組：1、1§1高斯消元法

/*GaussianElimination*/高斯消元法：首先將A化為上三角陣，再回代求解=§1高斯消元法/*GaussianEliminat2例用高斯消元法解方程組例用高斯消元法解方程組3消元記Step1：設，計算因子將增廣矩陣/*augmentedmatrix*/第i行mi1

第1行，得其中消元記Step1：設，計算因子將增廣矩陣/*4Stepk：設，計算因子且計算共進行？步n

1Stepk：設，計算因子共進行？步n5回代問題1、方程組有解的條件；問題2、什么情況下消去法能求解；問題3、求解的誤差估計。回代問題1、方程組有解的條件；6定理若方程組系數矩陣A的所有順序主子式均不為0，則高斯消元法能求得方程組的唯一解。注：事實上，只要A

非奇異，即A1

存在，則可通過逐次消元及行交換，將方程組化為三角形方程組，求出唯一解。定理若方程組系數矩陣A的所有順7定義設矩陣每一行對角元素的絕對值都大于同行其它元素絕對值之和，則稱為嚴格對角占優陣。例設矩陣為嚴格對角占優，經過一步Gauss消去得到其中則也是嚴格對角占優定理設方程組，如果為嚴格對角占優矩陣。則用Gauss消去法求解時，全不為零。定義設矩陣每一行對8選主元消去法1在求解線性方程組時，其系數矩陣絕大多數是非奇異的，但可能出現主元素

這時，消元過程無法進行

選主元的必要性2即使但如果其絕對值很小，也將會嚴重影響計算結果的精度。選主元消去法1在求解線性方程組時，其系數矩陣絕大多數9例：單精度解方程組（設機器字長4位）精確解為例：單精度解方程組（設機器字長4位）精確解為10用GaussianElimination計算：小主元可能導致計算失敗。用GaussianElimination計算：小主元可能11一、列主元消去法省去換列的步驟，每次僅選一列中最大的元。一、列主元消去法省去換列的步驟，每次僅選一列中最大的元。12二、全主元消去法每一步選絕對值最大的元素為主元素，保證。Stepk:①選取②Ifik

k

then交換第k行與第ik

行;Ifjk

k

then交換第k列與第jk

列;③消元注：列交換改變了xi

的順序，須記錄交換次序，解完后再換回來。二、全主元消去法每一步選絕對值最大的元素為主元素，保證13注：列主元法沒有全主元法穩定。注意：這兩個方程組在數學上嚴格等價。例：列主元法例：列主元法例：全主元法注：列主元法沒有全主元法穩定。注意：這兩個方程組在數學上嚴14/*Gausseliminationstep*/for(k=0;k<DIM;k++)for(i=k+1;i<DIM;i++){A[i][k]=A[i][k]/A[k][k];for(j=k+1;j<DIM;j++)A[i][j]=A[i][j]-A[i][k]*A[k][j];xx[i]=xx[i]-A[i][k]*xx[k];}/*Gausseliminationstep*/15for(k=0;k<DIM;k++){

pelement=fabs(A[k][k]);i0=k;for(i=k;i<DIM;i++)if(fabs(A[i][k])>pelement){

pelement=fabs(A[i][k]);i0=i;}if(i0!=k)/*若i0不等于k，交換i0,k兩行*/

{for(j=0;j<DIM;j++){pelement=A[k][j];A[k][j]=A[i0][j];A[i0][j]=pelement;}pelement=xx[k];xx[k]=xx[i0];xx[i0]=pelement;

}for(i=k+1;i<DIM;i++)

{A[i][k]=A[i][k]/A[k][k];for(j=k+1;j<DIM;j++)A[i][j]=A[i][j]-A[i][k]*A[k][j]; xx[i]=xx[i]-A[i][k]*xx[k];}}/*列主元消去法*/for(k=0;k<DIM;k++)/*列主元消去法*/16for(k=0;k<DIM;k++){

pelement=fabs(A[k][k]);i0=k;for(i=k;i<DIM;i++)if(fabs(A[i][k])>pelement)

{

pelement=fabs(A[i][k]);i0=i;}if(i0!=k)/*若i0不等于k，交換i0,k兩行*/

{for(j=0;j<DIM;j++){pelement=A[k][j];A[k][j]=A[i0][j];A[i0][j]=pelement;}pelement=xx[k];xx[k]=xx[i0];xx[i0]=pelement;}…}/*列主元消去法*/for(k=0;k<DIM;k++)/*列主元消去法*/17for(k=0;k<DIM;k++){…

{A[i][k]=A[i][k]/A[k][k];for(j=k+1;j<DIM;j++)A[i][j]=A[i][j]-A[i][k]*A[k][j]; xx[i]=xx[i]-A[i][k]*xx[k];

}}/*列主元消去法*/for(k=0;k<DIM;k++)/*列主元消去法*/18三、若當消去法若當消去法與高斯消去法的主要區別：把akk(k)

所在列的上、下元素全消為0;相當于求逆矩陣。三、若當消去法若當消去法與高斯消去法的主要區別：把19§2矩陣三角分解法一、高斯消元法的矩陣形式：Step1:記L1=，則即Gauss消元過程，消去第一列相當于左乘矩陣L1§2矩陣三角分解法一、高斯消元法的矩陣形式：Step120Stepn

1:其中

Lk=Stepn1:其中21對于Lk

有如下性質：1、記為L3、2、對于Lk有如下性質：1、記為L3、2、22記為L單位下三角陣記

U=A

的

LU

分解記為L單位下三角陣記U=A的LU分解23定理若A的所有順序主子式均不為0，則A

的

Doolittle--LU

分解唯一（其中L

為單位下三角陣，U為上三角陣）。注1：單位下三角*上三角稱為Doolittle分解

唯一

下三角*單位上三角稱為Crout分解。唯一下三角*上三角不唯一注2：矩陣的分解可以通過消元法來實現，L是所有乘積系數構成的矩陣，U是消元后的上三角定理若A的所有順序主子式均不24問題：前面講的是可以將A分解成為LU的形式。那么，為什么要分解呢？

注3：如果將A化成為LU的形式，則Ax=b求解很簡單：LUx=bLy=bUx=y問題又來了：用消去法解方程組就可以將A分解成為LU的形式。而分解的目的也是為了解方程組。方程組不是已經解出來了嗎？這不是循環邏輯嗎？如何解釋?！問題：前面講的是可以將A分解成為LU的形式。注3：如果將A化25二、直接分解法

——LU

分解的緊湊格式（杜立特分解法）思路：通過比較法直接導出L和

U的計算公式。L的第1行乘以U的每一列：L的每一行乘以U的第1列：二、直接分解法思路：通過比較法直接導出L和U的計算公式26所以，有：L的第2行乘以U的第i列：L的第i行乘以U的第2列：因此有：因此有：所以，有：L的第2行乘以U的第i列：L的第i行乘以U的第2列27一般設已得出U的前r-1行和L的前r-1列，則對于一般設已得出U的前r-1行和L的前r-1列，則對于28對

r=1,2,…,nstep1對計算計算出了U的r行的每個元素Doolittle分解算法：step2對計算計算出了L的r列的每個元素

在編程序時，為了節省存儲單元，將U的上三角元素存放在上三角部分，L的嚴格下三角元素存放在A的下三角部分，L的對角元素為1，不存儲。對r=1,2,…,nstep1對29Step3:方程Ly=b

的求解表達Step4:方程Ux=y

的求解表達下面求解，請寫出求解表達式Step3:方程Ly=b的求解表達Step30主要程序段for(r=0;r<DIM;r++)

{for(i=r;i<DIM;i++)for(jj=0;jj<r;jj++)A[r][i]=A[r][i]-A[r][jj]*A[jj][i];for(i=r+1;i<DIM;i++)

{for(jj=0;jj<r;jj++)A[i][r]=A[i][r]-A[i][jj]*A[jj][r];A[i][r]=A[i][r]/A[r][r];}

}/*DoolittleFactorization*/主要程序段for(r=0;r<DIM;r++)/*Dooli31for(i=DIM-1;i>=0;i--){for(j=i+1;j<DIM;j++)xx[i]=xx[i]-A[i][j]*xx[j];xx[i]=xx[i]/A[i][i];}/*solvetheequationLy=b*/for(i=1;i<DIM;i++)for(j=0;j<i;j++)xx[i]=xx[i]-A[i][j]*xx[j];/*solvetheequationUx=y*/for(i=DIM-1;i>=0;i--)/*solvet32例：利用系數矩陣的LU分解,求解方程組解：LU分解的緊湊格式為上一頁下一頁返回

例：利用系數矩陣的LU分解,求解方程組解：LU分解的緊湊格33則求解原方程組等價于求解下面兩個方程組Ly=b及Ux=y則求解原方程組等價于求解下面兩個方程組Ly=b及Ux=y34三、三對角方程組追趕法若A滿足Gauss消去法可行的條件，則可用LU分解法求解其中：三、三對角方程組追趕法若A滿足Gauss消去法可行的條件，則35解方程組Ax=d分為兩步，即求解Ly=d和Ux=y，計算公式如下：上述方法為求解三對角方程組的追趕法，也稱Thomas算法.追趕法公式簡單，計算量和存儲量都小。解方程組Ax=d分為兩步，即求解Ly=d和Ux=y，計算公式36§3平方根法—對稱正定矩陣的分解法定義一個矩陣A=(aij)nn

稱為對稱陣，如果aij=aji

。定義一個矩陣A

稱為正定陣，如果對任意非零向量都成立。對稱正定陣的幾個重要性質

A1

亦對稱正定，且aii>0

A

的順序主子陣Ak

亦對稱正定

A

的特征值i

>0

A

的全部順序主子式

det(Ak

)>0§3平方根法—對稱正定矩陣的分解法定義一個矩陣A=37

證明將對稱

正定陣

A

做LU

分解=u11uij/uii111u22unnU=uij記為

A對稱即因為A的順序主子式都大于0，所以uii>0定理設矩陣A對稱正定，則存在唯一的對角元為正的下三角陣L，使得。定理設n階對稱正定矩陣A，則存在唯一的單位下三角陣L及對角陣D

使得。證明將對稱正定陣A做LU分解=u11uij38將展開，可得因為同理，對j=i+1,i+2,…,n當r=i時有：

i=1,2,…,n當r=i=j時有：將展開，可得因為同理，對39算法

step1:的對A作Choleski分解

Forj=i+1,i+2,…,nFori=1,2,…,n算法step1:的對A作Choleski分解40step3:解方程組step2:解方程組Ly=bstep3:解方程組step2:解方程組Ly41第三章解線性方程組的直接法本章目標求解線性方程組：1、高斯消去法

高斯消去法選主元消去法約當消去法2、矩陣三角分解法直接分解法平方根法追趕法第三章解線性方程組的直接法本章目標求解線性方程組：1、42§1高斯消元法

/*GaussianElimination*/高斯消元法：首先將A化為上三角陣，再回代求解=§1高斯消元法/*GaussianEliminat43例用高斯消元法解方程組例用高斯消元法解方程組44消元記Step1：設，計算因子將增廣矩陣/*augmentedmatrix*/第i行mi1

第1行，得其中消元記Step1：設，計算因子將增廣矩陣/*45Stepk：設，計算因子且計算共進行？步n

1Stepk：設，計算因子共進行？步n46回代問題1、方程組有解的條件；問題2、什么情況下消去法能求解；問題3、求解的誤差估計。回代問題1、方程組有解的條件；47定理若方程組系數矩陣A的所有順序主子式均不為0，則高斯消元法能求得方程組的唯一解。注：事實上，只要A

非奇異，即A1

存在，則可通過逐次消元及行交換，將方程組化為三角形方程組，求出唯一解。定理若方程組系數矩陣A的所有順48定義設矩陣每一行對角元素的絕對值都大于同行其它元素絕對值之和，則稱為嚴格對角占優陣。例設矩陣為嚴格對角占優，經過一步Gauss消去得到其中則也是嚴格對角占優定理設方程組，如果為嚴格對角占優矩陣。則用Gauss消去法求解時，全不為零。定義設矩陣每一行對49選主元消去法1在求解線性方程組時，其系數矩陣絕大多數是非奇異的，但可能出現主元素

這時，消元過程無法進行

選主元的必要性2即使但如果其絕對值很小，也將會嚴重影響計算結果的精度。選主元消去法1在求解線性方程組時，其系數矩陣絕大多數50例：單精度解方程組（設機器字長4位）精確解為例：單精度解方程組（設機器字長4位）精確解為51用GaussianElimination計算：小主元可能導致計算失敗。用GaussianElimination計算：小主元可能52一、列主元消去法省去換列的步驟，每次僅選一列中最大的元。一、列主元消去法省去換列的步驟，每次僅選一列中最大的元。53二、全主元消去法每一步選絕對值最大的元素為主元素，保證。Stepk:①選取②Ifik

k

then交換第k行與第ik

行;Ifjk

k

then交換第k列與第jk

列;③消元注：列交換改變了xi

的順序，須記錄交換次序，解完后再換回來。二、全主元消去法每一步選絕對值最大的元素為主元素，保證54注：列主元法沒有全主元法穩定。注意：這兩個方程組在數學上嚴格等價。例：列主元法例：列主元法例：全主元法注：列主元法沒有全主元法穩定。注意：這兩個方程組在數學上嚴55/*Gausseliminationstep*/for(k=0;k<DIM;k++)for(i=k+1;i<DIM;i++){A[i][k]=A[i][k]/A[k][k];for(j=k+1;j<DIM;j++)A[i][j]=A[i][j]-A[i][k]*A[k][j];xx[i]=xx[i]-A[i][k]*xx[k];}/*Gausseliminationstep*/56for(k=0;k<DIM;k++){

pelement=fabs(A[k][k]);i0=k;for(i=k;i<DIM;i++)if(fabs(A[i][k])>pelement){

pelement=fabs(A[i][k]);i0=i;}if(i0!=k)/*若i0不等于k，交換i0,k兩行*/

{for(j=0;j<DIM;j++){pelement=A[k][j];A[k][j]=A[i0][j];A[i0][j]=pelement;}pelement=xx[k];xx[k]=xx[i0];xx[i0]=pelement;

}for(i=k+1;i<DIM;i++)

{A[i][k]=A[i][k]/A[k][k];for(j=k+1;j<DIM;j++)A[i][j]=A[i][j]-A[i][k]*A[k][j]; xx[i]=xx[i]-A[i][k]*xx[k];}}/*列主元消去法*/for(k=0;k<DIM;k++)/*列主元消去法*/57for(k=0;k<DIM;k++){

pelement=fabs(A[k][k]);i0=k;for(i=k;i<DIM;i++)if(fabs(A[i][k])>pelement)

{

pelement=fabs(A[i][k]);i0=i;}if(i0!=k)/*若i0不等于k，交換i0,k兩行*/

{for(j=0;j<DIM;j++){pelement=A[k][j];A[k][j]=A[i0][j];A[i0][j]=pelement;}pelement=xx[k];xx[k]=xx[i0];xx[i0]=pelement;}…}/*列主元消去法*/for(k=0;k<DIM;k++)/*列主元消去法*/58for(k=0;k<DIM;k++){…

{A[i][k]=A[i][k]/A[k][k];for(j=k+1;j<DIM;j++)A[i][j]=A[i][j]-A[i][k]*A[k][j]; xx[i]=xx[i]-A[i][k]*xx[k];

}}/*列主元消去法*/for(k=0;k<DIM;k++)/*列主元消去法*/59三、若當消去法若當消去法與高斯消去法的主要區別：把akk(k)

所在列的上、下元素全消為0;相當于求逆矩陣。三、若當消去法若當消去法與高斯消去法的主要區別：把60§2矩陣三角分解法一、高斯消元法的矩陣形式：Step1:記L1=，則即Gauss消元過程，消去第一列相當于左乘矩陣L1§2矩陣三角分解法一、高斯消元法的矩陣形式：Step161Stepn

1:其中

Lk=Stepn1:其中62對于Lk

有如下性質：1、記為L3、2、對于Lk有如下性質：1、記為L3、2、63記為L單位下三角陣記

U=A

的

LU

分解記為L單位下三角陣記U=A的LU分解64定理若A的所有順序主子式均不為0，則A

的

Doolittle--LU

分解唯一（其中L

為單位下三角陣，U為上三角陣）。注1：單位下三角*上三角稱為Doolittle分解

唯一

下三角*單位上三角稱為Crout分解。唯一下三角*上三角不唯一注2：矩陣的分解可以通過消元法來實現，L是所有乘積系數構成的矩陣，U是消元后的上三角定理若A的所有順序主子式均不65問題：前面講的是可以將A分解成為LU的形式。那么，為什么要分解呢？

注3：如果將A化成為LU的形式，則Ax=b求解很簡單：LUx=bLy=bUx=y問題又來了：用消去法解方程組就可以將A分解成為LU的形式。而分解的目的也是為了解方程組。方程組不是已經解出來了嗎？這不是循環邏輯嗎？如何解釋?！問題：前面講的是可以將A分解成為LU的形式。注3：如果將A化66二、直接分解法

——LU

分解的緊湊格式（杜立特分解法）思路：通過比較法直接導出L和

U的計算公式。L的第1行乘以U的每一列：L的每一行乘以U的第1列：二、直接分解法思路：通過比較法直接導出L和U的計算公式67所以，有：L的第2行乘以U的第i列：L的第i行乘以U的第2列：因此有：因此有：所以，有：L的第2行乘以U的第i列：L的第i行乘以U的第2列68一般設已得出U的前r-1行和L的前r-1列，則對于一般設已得出U的前r-1行和L的前r-1列，則對于69對

r=1,2,…,nstep1對計算計算出了U的r行的每個元素Doolittle分解算法：step2對計算計算出了L的r列的每個元素

在編程序時，為了節省存儲單元，將U的上三角元素存放在上三角部分，L的嚴格下三角元素存放在A的下三角部分，L的對角元素為1，不存儲。對r=1,2,…,nstep1對70Step3:方程Ly=b

的求解表達Step4:方程Ux=y

的求解表達下面求解，請寫出求解表達式Step3:方程Ly=b的求解表達Step71主要程序段for(r=0;r<DIM;r++)

{for(i=r;i<DIM;i++)for(jj=0;jj<r;jj++)A[r][i]=A[r][i]-A[r][jj]*A[jj][i];for(i=r+1;i<DIM;i++)

{for(jj=0;jj<r;jj++)A[i][r]=A[i][r]-A[i][jj]*A[jj][r];A[i][r]=A[i][r]/A[r][r];}

}/*DoolittleFactorization*/主要程序段for(r=0;r<DIM;r++)/*Dooli72for(i=DIM-1;i>=0;i--){for(j=i+1;j<DIM;j++)xx[i]=xx[i]-A[i][j]*xx[j];xx[i]=xx[i]/A[i][i];}/*solvetheequationLy=b*/for(i=1;i<DIM;i++)for(j=0;j<i;j++)xx[i]=xx[i]-A[i][j]*xx