3、線性代數(shù)方程組的直接解法課件

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

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

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

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

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

非奇異，即A1

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

這時(shí)，消元過(guò)程無(wú)法進(jìn)行

選主元的必要性2即使但如果其絕對(duì)值很小，也將會(huì)嚴(yán)重影響計(jì)算結(jié)果的精度。選主元消去法1在求解線性方程組時(shí)，其系數(shù)矩陣絕大多數(shù)9例：?jiǎn)尉冉夥匠探M（設(shè)機(jī)器字長(zhǎng)4位）精確解為例：?jiǎn)尉冉夥匠探M（設(shè)機(jī)器字長(zhǎng)4位）精確解為10用GaussianElimination計(jì)算：小主元可能導(dǎo)致計(jì)算失敗。用GaussianElimination計(jì)算：小主元可能11一、列主元消去法省去換列的步驟，每次僅選一列中最大的元。一、列主元消去法省去換列的步驟，每次僅選一列中最大的元。12二、全主元消去法每一步選絕對(duì)值最大的元素為主元素，保證。Stepk:①選?、贗fik

k

then交換第k行與第ik

行;Ifjk

k

then交換第k列與第jk

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

的順序，須記錄交換次序，解完后再換回來(lái)。二、全主元消去法每一步選絕對(duì)值最大的元素為主元素，保證13注：列主元法沒(méi)有全主元法穩(wěn)定。注意：這兩個(gè)方程組在數(shù)學(xué)上嚴(yán)格等價(jià)。例：列主元法例：列主元法例：全主元法注：列主元法沒(méi)有全主元法穩(wěn)定。注意：這兩個(gè)方程組在數(shù)學(xué)上嚴(yán)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三、若當(dāng)消去法若當(dāng)消去法與高斯消去法的主要區(qū)別：把a(bǔ)kk(k)

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

1:其中

Lk=Stepn1:其中21對(duì)于Lk

有如下性質(zhì)：1、記為L(zhǎng)3、2、對(duì)于Lk有如下性質(zhì)：1、記為L(zhǎng)3、2、22記為L(zhǎng)單位下三角陣記

U=A

的

LU

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

的

Doolittle--LU

分解唯一（其中L

為單位下三角陣，U為上三角陣）。注1：?jiǎn)挝幌氯?上三角稱為Doolittle分解

唯一

下三角*單位上三角稱為Crout分解。唯一下三角*上三角不唯一注2：矩陣的分解可以通過(guò)消元法來(lái)實(shí)現(xiàn)，L是所有乘積系數(shù)構(gòu)成的矩陣，U是消元后的上三角定理若A的所有順序主子式均不24問(wèn)題：前面講的是可以將A分解成為L(zhǎng)U的形式。那么，為什么要分解呢？

注3：如果將A化成為L(zhǎng)U的形式，則Ax=b求解很簡(jiǎn)單：LUx=bLy=bUx=y問(wèn)題又來(lái)了：用消去法解方程組就可以將A分解成為L(zhǎng)U的形式。而分解的目的也是為了解方程組。方程組不是已經(jīng)解出來(lái)了嗎？這不是循環(huán)邏輯嗎？如何解釋?！問(wèn)題：前面講的是可以將A分解成為L(zhǎng)U的形式。注3：如果將A化25二、直接分解法

——LU

分解的緊湊格式（杜立特分解法）思路：通過(guò)比較法直接導(dǎo)出L和

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

r=1,2,…,nstep1對(duì)計(jì)算計(jì)算出了U的r行的每個(gè)元素Doolittle分解算法：step2對(duì)計(jì)算計(jì)算出了L的r列的每個(gè)元素

在編程序時(shí)，為了節(jié)省存儲(chǔ)單元，將U的上三角元素存放在上三角部分，L的嚴(yán)格下三角元素存放在A的下三角部分，L的對(duì)角元素為1，不存儲(chǔ)。對(duì)r=1,2,…,nstep1對(duì)29Step3:方程Ly=b

的求解表達(dá)Step4:方程Ux=y

的求解表達(dá)下面求解，請(qǐng)寫(xiě)出求解表達(dá)式Step3:方程Ly=b的求解表達(dá)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例：利用系數(shù)矩陣的LU分解,求解方程組解：LU分解的緊湊格式為上一頁(yè)下一頁(yè)返回

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

稱為對(duì)稱陣，如果aij=aji

。定義一個(gè)矩陣A

稱為正定陣，如果對(duì)任意非零向量都成立。對(duì)稱正定陣的幾個(gè)重要性質(zhì)

A1

亦對(duì)稱正定，且aii>0

A

的順序主子陣Ak

亦對(duì)稱正定

A

的特征值i

>0

A

的全部順序主子式

det(Ak

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

證明將對(duì)稱

正定陣

A

做LU

分解=u11uij/uii111u22unnU=uij記為

A對(duì)稱即因?yàn)锳的順序主子式都大于0，所以u(píng)ii>0定理設(shè)矩陣A對(duì)稱正定，則存在唯一的對(duì)角元為正的下三角陣L，使得。定理設(shè)n階對(duì)稱正定矩陣A，則存在唯一的單位下三角陣L及對(duì)角陣D

使得。證明將對(duì)稱正定陣A做LU分解=u11uij38將展開(kāi)，可得因?yàn)橥?，?duì)j=i+1,i+2,…,n當(dāng)r=i時(shí)有：

i=1,2,…,n當(dāng)r=i=j時(shí)有：將展開(kāi)，可得因?yàn)橥恚瑢?duì)39算法

step1:的對(duì)A作Choleski分解

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

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

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

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

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

非奇異，即A1

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

這時(shí)，消元過(guò)程無(wú)法進(jìn)行

選主元的必要性2即使但如果其絕對(duì)值很小，也將會(huì)嚴(yán)重影響計(jì)算結(jié)果的精度。選主元消去法1在求解線性方程組時(shí)，其系數(shù)矩陣絕大多數(shù)50例：?jiǎn)尉冉夥匠探M（設(shè)機(jī)器字長(zhǎng)4位）精確解為例：?jiǎn)尉冉夥匠探M（設(shè)機(jī)器字長(zhǎng)4位）精確解為51用GaussianElimination計(jì)算：小主元可能導(dǎo)致計(jì)算失敗。用GaussianElimination計(jì)算：小主元可能52一、列主元消去法省去換列的步驟，每次僅選一列中最大的元。一、列主元消去法省去換列的步驟，每次僅選一列中最大的元。53二、全主元消去法每一步選絕對(duì)值最大的元素為主元素，保證。Stepk:①選?、贗fik

k

then交換第k行與第ik

行;Ifjk

k

then交換第k列與第jk

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

的順序，須記錄交換次序，解完后再換回來(lái)。二、全主元消去法每一步選絕對(duì)值最大的元素為主元素，保證54注：列主元法沒(méi)有全主元法穩(wěn)定。注意：這兩個(gè)方程組在數(shù)學(xué)上嚴(yán)格等價(jià)。例：列主元法例：列主元法例：全主元法注：列主元法沒(méi)有全主元法穩(wěn)定。注意：這兩個(gè)方程組在數(shù)學(xué)上嚴(yán)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三、若當(dāng)消去法若當(dāng)消去法與高斯消去法的主要區(qū)別：把a(bǔ)kk(k)

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

1:其中

Lk=Stepn1:其中62對(duì)于Lk

有如下性質(zhì)：1、記為L(zhǎng)3、2、對(duì)于Lk有如下性質(zhì)：1、記為L(zhǎng)3、2、63記為L(zhǎng)單位下三角陣記

U=A

的

LU

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

的

Doolittle--LU

分解唯一（其中L

為單位下三角陣，U為上三角陣）。注1：?jiǎn)挝幌氯?上三角稱為Doolittle分解

唯一

下三角*單位上三角稱為Crout分解。唯一下三角*上三角不唯一注2：矩陣的分解可以通過(guò)消元法來(lái)實(shí)現(xiàn)，L是所有乘積系數(shù)構(gòu)成的矩陣，U是消元后的上三角定理若A的所有順序主子式均不65問(wèn)題：前面講的是可以將A分解成為L(zhǎng)U的形式。那么，為什么要分解呢？

注3：如果將A化成為L(zhǎng)U的形式，則Ax=b求解很簡(jiǎn)單：LUx=bLy=bUx=y問(wèn)題又來(lái)了：用消去法解方程組就可以將A分解成為L(zhǎng)U的形式。而分解的目的也是為了解方程組。方程組不是已經(jīng)解出來(lái)了嗎？這不是循環(huán)邏輯嗎？如何解釋?！問(wèn)題：前面講的是可以將A分解成為L(zhǎng)U的形式。注3：如果將A化66二、直接分解法

——LU

分解的緊湊格式（杜立特分解法）思路：通過(guò)比較法直接導(dǎo)出L和

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

r=1,2,…,nstep1對(duì)計(jì)算計(jì)算出了U的r行的每個(gè)元素Doolittle分解算法：step2對(duì)計(jì)算計(jì)算出了L的r列的每個(gè)元素

在編程序時(shí)，為了節(jié)省存儲(chǔ)單元，將U的上三角元素存放在上三角部分，L的嚴(yán)格下三角元素存放在A的下三角部分，L的對(duì)角元素為1，不存儲(chǔ)。對(duì)r=1,2,…,nstep1對(duì)70Step3:方程Ly=b

的求解表達(dá)Step4:方程Ux=y

的求解表達(dá)下面求解，請(qǐng)寫(xiě)出求解表達(dá)式Step3:方程Ly=b的求解表達(dá)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