3、线性代数方程组的直接解法课件

第三章解线性方程组的直接法本章目标求解线性方程组：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