解线性方程组的直接法

解线性方程组的直接法第一页，共一百八十页，编辑于2023年，星期二

自然科学和工程计算中的很多问题的解决常常归结为求解线性方程组。如三次样条插值函数问题、用最小二乘远离确定拟合曲线、求解微分方程的数值解等，最终都要转化为求解线性方程组。求解线性方程组可采用：5.1引言与预备知识直接法——假定计算过程没有舍入误差的情况下，经过有限步算术运算后能求得线性方程组精确解的方法。经过有限步运算就能求得精确解的方法，但实际计算中由于舍入误差的影响，这类方法也只能求得近似解；例如：高斯消去法、三角分解法等。迭代法——构造适当的向量序列，用某种极限过程去逐步逼近精确解。例如:雅可比迭代法、高斯-赛德尔迭代法等。第二页，共一百八十页，编辑于2023年，星期二一般地：对低阶方程组用直接法，对高阶方程组和大型稀疏方程组用迭代法求解。预备知识:1.向量与矩阵：第三页，共一百八十页，编辑于2023年，星期二矩阵的基本运算(4)转置矩阵(5)定义:

对于阶矩阵，如果有一个阶矩阵,

则说矩阵是可逆的，并把矩阵称为的逆矩阵.使得第四页，共一百八十页，编辑于2023年，星期二行列式性质(6)行列式按任一行展开其中Ai,j为ai,j的代数余子式，Ai,j=(-1)i+jMi,j,Mi,j为ai,j的余子式.第五页，共一百八十页，编辑于2023年，星期二2.特殊矩阵(1)对角矩阵：如果ij时，ai,j=0.(2)三对角矩阵：如果|i-j|>1时，ai,j=0.形如(3)上三角矩阵：如果i>j时，ai,j=0.形如（类似有下三角矩阵）第六页，共一百八十页，编辑于2023年，星期二(4)上海森伯格(Hessenberg)矩阵：如果i>j+1时，ai,j=0.形如(5)对称矩阵：如果AT=A.(6)埃尔米特矩阵：设ACn×n,如果AH=A.(7)对称正定矩阵：如果AT=A且对任意非零向量第七页，共一百八十页，编辑于2023年，星期二(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.第八页，共一百八十页，编辑于2023年，星期二定理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为对称正定阵。第九页，共一百八十页，编辑于2023年，星期二定理4(Jordan标准型）：设A为n阶矩阵,则存在一个非奇异矩阵P使得其中(1)当A的若当标准型中所有若当块Ji均为一阶时，此标准型若当(Jordan)块变成对角阵。(2)若A的特征值各不相同，则其若当标准型必为对角阵diag(1,2,…,n).第十页，共一百八十页，编辑于2023年，星期二设线性方程组为或写成矩阵形式或简单地记为:5.2高斯消元法第十一页，共一百八十页，编辑于2023年，星期二第十二页，共一百八十页，编辑于2023年，星期二先看几类特殊三角形方程组的解

★对角形方程组若aii≠0,则

xi=bi/aii,i=1,2,…,n。运算量O(n)。第十三页，共一百八十页，编辑于2023年，星期二★下三角方程组（返回LU）(5.19)第十四页，共一百八十页，编辑于2023年，星期二运算量：因为计算xi需要i

次乘法和除法，所以★上三角方程组（返回Gauss）返回LU(5.20)第十五页，共一百八十页，编辑于2023年，星期二运算量：因为计算xi

需要n-i

次乘法和1次除法，即

n–i+1次乘除，所以第十六页，共一百八十页，编辑于2023年，星期二

消元法的基本思想就是通过对方程组作初等变换，把一般形式的线性方程组化为等价的易于求解的三角方程组。★高斯消元法（GaussianElimination）高斯消元法是一种古老的求解线性方程组的方法，它就是通过一系列的初等变换（消元），把线性方程组（5.1）化为等价的上三角方程组（5.5），然后通过回代方法求出原方程组的解。顺序高斯消去法消元过程:依从左到右、自上而下的次序将主对角元下方的元素化为零。1不作行交换。2用不等于零的数乘某行，加至另一行。第十七页，共一百八十页，编辑于2023年，星期二第十八页，共一百八十页，编辑于2023年，星期二第十九页，共一百八十页，编辑于2023年，星期二高斯消去法的主要思路：将系数矩阵A化为上三角矩阵，然后回代求解。考虑n阶线性方程组：矩阵形式=下面我们讨论一般线性方程组的高斯消去法第二十页，共一百八十页，编辑于2023年，星期二对于第二十一页，共一百八十页，编辑于2023年，星期二重复上述过程，可以得到与（5.1）等价的上三角方程组：第二十二页，共一百八十页，编辑于2023年，星期二或写成矩阵形式A(n-1)x=b(n-1),

其中，第二十三页，共一百八十页，编辑于2023年，星期二约定由（5.8）逐个回代，可得（5.1）的解第二十四页，共一百八十页，编辑于2023年，星期二此即Gauss消去法.

综述:若A非奇,总可通过带有行交换或不带有行交换的消元过程,将A化成非奇三角矩阵A(n-1).因此,回代求解过程也可进行到底.但实际中,A是否非奇异难以判定.算法应考虑:消元过程某一步找不到非零元,于是计算中断.消元可进行到底,但ann(n-1)=

0,回代求解过程无法进行.故算法设计要考虑以上情况,给出计算中断的信息.第二十五页，共一百八十页，编辑于2023年，星期二

算法组织将增广矩阵(A,b)放入一个二维数组.消去每一步算出的aij(k)放入aij的位置,bi(k)放入bi,mik放在aik上.消去完毕后,数组各位置上的数据如下示:回代过程的计算没有数据组织问题.Ab第二十六页，共一百八十页，编辑于2023年，星期二算法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/akkaik

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

aik

akj

aij(新)

1.2.3bi-aikbkbi(新)说明:系数矩阵存放于数组A中,右端向量放于数组b中*常用|akk|≤第二十七页，共一百八十页，编辑于2023年，星期二步骤2.bn/ann

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

3.2.1s-akjxjs3.3s/akkxk步骤4.End第二十八页，共一百八十页，编辑于2023年，星期二

★高斯消元法运算量：消元过程：乘：

除：回代过程：运算量：第二十九页，共一百八十页，编辑于2023年，星期二★高斯消元法的局限性：在高斯消元过程中，我们假定了对角元素由于每次消元时是按未知量的自然顺序进行的，而顺序消元不改变A的主子式的值，因此高斯消元法可行的充分必要条件为A的各阶主子式不为零。实际上只要detA≠0，方程组就有解。2.

即使高斯消元法可行，如果很小，运算中用它作分母会导致其它元素数量计的严重增长和舍入误差的扩散。第三十页，共一百八十页，编辑于2023年，星期二顺序高斯消去法的使用条件

使用条件之一

定理线性方程组系数矩阵A的顺序主子矩阵Ak

(k=1,2,…,n)非奇异，则顺序高斯消去法能实现方程组的求解。即方程组能用顺序高斯消去法求解的充要条件是系数行列式的顺序主子式非零。在原方程组的系数矩阵中如何反映出这个条件呢?我们知道高斯消去法能按顺序进行到底的充要条件应是：第三十一页，共一百八十页，编辑于2023年，星期二A的k阶顺序主子矩阵Ak的行列式推论：如果A的顺序主子矩阵的行列式detAk≠0

(k=1,2,…,n-1)，则第三十二页，共一百八十页，编辑于2023年，星期二

使用条件之二

n阶矩阵A为严格对角占优矩阵是指其每个主对角元的绝对值大于同一行其他元素绝对值之和，即一阶严格对角占优矩阵指一个非零数。

定理

方程组系数矩阵A为严格对角占优矩阵则可实现用顺序高斯消去法求解。第三十三页，共一百八十页，编辑于2023年，星期二Gauss主元素消元法例用Gauss消去法求解线性方程组解:[方法1]方程组的准确解为x*=(9998/9999,10000/9999)T.小主元a11回代求解x2=1.0000,x1=0.0000.计算结果相当糟糕.

原因:求行乘数时用了”小主元”0.0001作除数.第三十四页，共一百八十页，编辑于2023年，星期二解:[方法2]若上题在求解时将1,2行交换,即回代后得到

x=(1.0000,1.0000)T与准确解x*=(9998/9999,10000/9999)T相差不大.主元a11第三十五页，共一百八十页，编辑于2023年，星期二再如：解方程组

精确解为：x1=1/3，x2=2/3。计算取5位有效数字。解(1)顺序消元：

(2)交换方程的顺序：经高斯消元：法（2）较好第三十六页，共一百八十页，编辑于2023年，星期二可见，由于计算机和方程组本身的限制，在Gauss消元过程中会出现零主元和小主元，而使Gauss消去法无法进行或出现较大舍入误差，为避免此类现象，可以通过行交换来选取绝对值较大的主元素.

前例方法2中所用的是在Gauss消去法中加入选主元过程,利用换行,避免”小主元”作除数,该方法称为Gauss列主元消去法。第三十七页，共一百八十页，编辑于2023年，星期二★Gauss列主元消元法高斯消元过程中的称为主元素。如果在一列中选取按模最大的元素作为主元素，将其调换到主干方程位置再做消元，则称为列主元消元法。第一步,

先在第1列选出绝对值最大的元素,

交换第1行和第m行的所有元素，再做消元操作。第二步，对每个k=1，2，…，n做消元前，选出中绝对值最大的元素，对k行和m行交换后，再作消元。第三十八页，共一百八十页，编辑于2023年，星期二列以下(含主对角元)元素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,列主元消去法能避免”小主元”作除数,误差分析与数值试验表明：（Gauss）列主元消去法”基本稳定”.第三十九页，共一百八十页，编辑于2023年，星期二列主元消去法中方程对换问题引入向量P=(p1,p2,…,pn)t=(1,2,…,n)t即ipi记录方程(5.1)的行号.当第k步消元若需第k个方程与第ik个方程对换,则只需交换P中当消元过程完后,向量P中的pi(i=1,…,n)指示初始方程(5.1)在作的i步消元时被选为主元的方程序号.这种序号交换避免了方程间的对换,n较大时,有其优越性.第四十页，共一百八十页，编辑于2023年，星期二例

用列主元高斯消去法求解方程组(用三位有效数字计算)解选主元选主元第四十一页，共一百八十页，编辑于2023年，星期二消元过程完成，得到上三角形方程组再作回代可求得第四十二页，共一百八十页，编辑于2023年，星期二开始输出无解信息消元换行停机回代求解Gauss列主元消去法的算法设计流程图第四十三页，共一百八十页，编辑于2023年，星期二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说明:系数矩阵存放于数组A中,右端向量放于数组b中.算法存储,I/O同Guass消去法

第四十四页，共一百八十页，编辑于2023年，星期二步骤3.(回代)步骤4.Fori=n-1,…,13.13.2Forj=k+1,…,n

3.2.13.3步骤5.End第四十五页，共一百八十页，编辑于2023年，星期二

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

Then输出”方程奇异”stopIf

k≠rThen

bk

br;akj

arj(j=k,…,n)说明:列主元算法在实际中应用较广“全面选主元”算法第四十六页，共一百八十页，编辑于2023年，星期二算法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/akkaik

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

aik

akj

aij(新)

1.2.3bi-aikbkbi(新)|akk|≤第四十七页，共一百八十页，编辑于2023年，星期二步骤2.bn/ann

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

3.2.1s-akjxjs3.3s/akkxk步骤4.End第四十八页，共一百八十页，编辑于2023年，星期二全主元高斯消去法：第k步消元时选A(k)中绝对值最大的元素为主元，即①先选取全主元：ifik

k

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

k

then交换第k列和第jk列；③消元列交换改变了

xi

的顺序，须记录交换次序，解完后再换回来。全主元高斯消去法具有很好的稳定性，但选全主元比较费时，故在实际计算中很少使用。第四十九页，共一百八十页，编辑于2023年，星期二程序示例★列主元高斯消元法算法描述：

1.输入数据：n，A，b2.分解过程：（1）对k列选主元，k=1，2，…，n-1

（2）交换第k行和第m行方程；（3）生成高斯消元过程的新矩阵和右端项第五十页，共一百八十页，编辑于2023年，星期二

3.回代过程：对于i=n,n-1,…,2,1，解上三角方程组

4.输出

第五十一页，共一百八十页，编辑于2023年，星期二程序源码主要部分：输入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]);}第五十二页，共一百八十页，编辑于2023年，星期二找主元素并进行交换两行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的元素

} }第五十三页，共一百八十页，编辑于2023年，星期二高斯消元后生成新的矩阵元素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;}第五十四页，共一百八十页，编辑于2023年，星期二回代解方程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];}第五十五页，共一百八十页，编辑于2023年，星期二2.Gauss列主元消去法消元过程的矩阵描述由于Gauss列主元消去法每一步都要选取列主元,因此不可避免要进行行交换即表示不换行初等矩阵第五十六页，共一百八十页，编辑于2023年，星期二因此,Gauss列主元消去法的消元过程为:显然上三角阵仍然为单位下三角矩阵第五十七页，共一百八十页，编辑于2023年，星期二初等矩阵的乘积,称为排列阵则推广到一般情形令仍然为单位下三角矩阵则单位下三角阵与上三角阵的乘积第五十八页，共一百八十页，编辑于2023年，星期二综合以上讨论,有定理.(列主元素的三角分解定理):

第五十九页，共一百八十页，编辑于2023年，星期二Gauss消去法是将系数矩阵化为上三角矩阵，再进行回代求解；Jordan法是将系数矩阵化为对角矩阵，从而无须再进行回代求解的一种直接解法。这种解法基本过程与Gauss法相同，只是在消去过程中，不仅将主对角线以下的元素消成0，而且同时将主对角线以上的元素也消成0。高斯-若当(Gauss-Jordan)消元法第六十页，共一百八十页，编辑于2023年，星期二第4行分别乘分别加到第3，2，1行上以n=4为例说明高斯-若当消元过程。首先高斯消元：第六十一页，共一百八十页，编辑于2023年，星期二第六十二页，共一百八十页，编辑于2023年，星期二第六十三页，共一百八十页，编辑于2023年，星期二第六十四页，共一百八十页，编辑于2023年，星期二第六十五页，共一百八十页，编辑于2023年，星期二归一方法的例子第六十六页，共一百八十页，编辑于2023年，星期二列选主高斯-约当(Jordan)消去法。第六十七页，共一百八十页，编辑于2023年，星期二第六十八页，共一百八十页，编辑于2023年，星期二高斯-约当(Jordan)消去法求逆。第六十九页，共一百八十页，编辑于2023年，星期二第七十页，共一百八十页，编辑于2023年，星期二5.3

矩阵三角分解法5.3.1

直接三角分解法

5.3.2平方根法5.3.3

追赶法第七十一页，共一百八十页，编辑于2023年，星期二

定义将矩阵A分解成一个下三角阵L和一个上三角阵U的乘积，即A=LU，称为A的三角分解或LU分解。第七十二页，共一百八十页，编辑于2023年，星期二例如A=这里有A的两种不同的三角分解，类似可举出很多，一般，若A=LU是一个三角分解，任取与A同阶的非奇异对角矩阵D，则A=(LD)(D-1U)=L1U1也是A的三角分解。第七十三页，共一百八十页，编辑于2023年，星期二。杜利特尔分解(Doolittle)常用的两种三角分解克洛特Crout分解有的叫库朗(Courant)分解第七十四页，共一百八十页，编辑于2023年，星期二5.3.1

直接三角分解法

用矩阵相乘来解释高斯消元过程，取n=4。记(5.1)化为方程组(5.6)的过程，即等价于L1A=A(1),L1b=b(1),即第七十五页，共一百八十页，编辑于2023年，星期二记(5.6)化为方程组(5.7)的过程等价于L2A(1)

=A(2),L2b(1)=b(2),即第七十六页，共一百八十页，编辑于2023年，星期二记第七十七页，共一百八十页，编辑于2023年，星期二因此其中，=LU第七十八页，共一百八十页，编辑于2023年，星期二第七十九页，共一百八十页，编辑于2023年，星期二一般地，高斯消元完成后,有：故从而U=A(n-1)第八十页，共一百八十页，编辑于2023年，星期二即且顺序主元第八十一页，共一百八十页，编辑于2023年，星期二解:由上述分析不难得到

问题：矩阵A存在LU分解（即Gauss消去法可以执行）的条件是什么？第八十二页，共一百八十页，编辑于2023年，星期二Gauss消去法可以执行定理

证明（略）.推论(P172)：如果A的顺序主子矩阵的行列式detAk≠0

(k=1,2,…,n-1)，则第八十三页，共一百八十页，编辑于2023年，星期二由上述对系数矩阵A左乘一系列的三角初等矩阵知，A可以分解为一个单位下三角矩阵L和一个上三角矩阵U。这个过程派生出解线性方程组的直接分解法。一旦实现A=LU，则Ax=b可以化为LUx=b。令Ux=y，则Ly=b。由Ly=b（即(5.4)）解出y；再由Ux=y（即(5.5)）解出x。如果A的各阶主子式不为零，则可以实现：杜利特尔(Doolittle)分解：如果L为单位下三角矩阵，U为上三角矩阵；克洛特Crout分解：如果L为下三角矩阵,U为单位上三角矩阵。第八十四页，共一百八十页，编辑于2023年，星期二杜利特尔分解法设A

的各阶主子式不为零，A分解为A=LU，即先计算U的元素，后计算L的元素：U的第1行：第八十五页，共一百八十页，编辑于2023年，星期二L的第1列：第八十六页，共一百八十页，编辑于2023年，星期二再计算U的第2行元素；计算L的第2列元素；……计算U的第k行元素：第八十七页，共一百八十页，编辑于2023年，星期二计算L的第k列元素：第八十八页，共一百八十页，编辑于2023年，星期二综合以上分析,可以推导出：U的第一行L的第一列------(1)------(2)U的第r行------(3)(逐行算出)L的第r列

------(4)(逐列算出)上述(1)~(4)式所表示的分解过程称为A的Doolittle分解LU分解的运算量：第八十九页，共一百八十页，编辑于2023年，星期二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)求解总运算量：第九十页，共一百八十页，编辑于2023年，星期二举例第九十一页，共一百八十页，编辑于2023年，星期二用紧凑格式解：第九十二页，共一百八十页，编辑于2023年，星期二再如：第九十三页，共一百八十页，编辑于2023年，星期二紧凑格式分解A=LU总结：1.先确定U中第一行元素（即等于A中第一行元素）.2.再确定L中第一列元素：3.确定U中第r行元素:4.再确定L中第r列元素:第九十四页，共一百八十页，编辑于2023年，星期二再如：作分解A=LU第九十五页，共一百八十页，编辑于2023年，星期二再如：第九十六页，共一百八十页，编辑于2023年，星期二对于线性方程组系数矩阵非奇异,经过Doolittle分解后Ax=L(Ux)=b可化为下面两个三角形方程组消去回代L=第九十七页，共一百八十页，编辑于2023年，星期二上述解线性方程组的方法称为三角(LU)分解的Doolittle法.第九十八页，共一百八十页，编辑于2023年，星期二Doolittle法的特点:紧凑(无中间过程);内积计算(精度高)例1：

用Doolittle法解方程组解:由Doolittle分解逐行算出U的元素逐列算出L的元素第九十九页，共一百八十页，编辑于2023年，星期二逐行算出U的元素逐列算出L的元素第一百页，共一百八十页，编辑于2023年，星期二例2用多利特尔分解求解方程组：解

设A=LU，即

第一百零一页，共一百八十页，编辑于2023年，星期二解下三角方程组Ly=b，即解上三角方程组Ux=y，即第一百零二页，共一百八十页，编辑于2023年，星期二A,b,x,L,U,y需要单独存储,而从lij,uij的计算过程发现:算法的数据组织:对于线性方程，第一百零三页，共一百八十页，编辑于2023年，星期二Doolittle分解法的数据组织可以用以下紧凑格式表示:第一百零四页，共一百八十页，编辑于2023年，星期二UL第一百零五页，共一百八十页，编辑于2023年，星期二Doolittle法解本节方程组(例1)的数据组织过程：解:第一百零六页，共一百八十页，编辑于2023年，星期二第一百零七页，共一百八十页，编辑于2023年，星期二所以x第一百零八页，共一百八十页，编辑于2023年，星期二例4用杜利特尔分解法求解方程组解先求增广系数矩阵的杜利特尔分解，即第一百零九页，共一百八十页，编辑于2023年，星期二再一例：第一百一十页，共一百八十页，编辑于2023年，星期二第一百一十一页，共一百八十页，编辑于2023年，星期二开始输出错误信息停机1.Doolittle法解方程组流程图程序实现计算U的第k行计算L的第k列输出x由Ly=b计算yk输出Uk,Lk,yk由Ux=y计算xk第一百一十二页，共一百八十页，编辑于2023年，星期二2.Doolittle分解主要算法程序：输入：n，A，b计算U的第k行元素和L的第k列元素计算U的第k行计算L的第k列3.程序源代码：（Doolim.C，Dooli.C）回代求解第一百一十三页，共一百八十页，编辑于2023年，星期二小结

Doolittle求解法的特点:紧凑，比高斯消去法精度高.2.实际计算中，由于舍入误差的影响，不可避免地往往只能求得近似解.3.为保证LU分解能够进行，要求A的各阶顺序主子式不为零，为取消此限制，可采用类似列主元消元法处理，程序稍复杂，或采用其它算法（如迭代法）。4.本实验是后续改进算法(平方根法、Cholesky分解法)的基础。对于大型方程组，一般采用迭代算法（下一章）.第一百一十四页，共一百八十页，编辑于2023年，星期二问题:在Doolittle法(包括LU分解紧凑格式)中,会反复用到公式仍有可能为小主元做除数?为此,也应考虑在算法中加入选取列主元即紧凑格式的Doolittle列主元法.第一百一十五页，共一百八十页，编辑于2023年，星期二列主元Doolittle法步骤:第一步:第一百一十六页，共一百八十页，编辑于2023年，星期二例.用列主元Doolittle法解线性方程组解:第一百一十七页，共一百八十页，编辑于2023年，星期二所以原方程组的解为第一百一十八页，共一百八十页，编辑于2023年，星期二克洛特Crout分解*(如果L为下三角矩阵，U为单位上三角矩阵)

先计算L的第一列，再计算U的第一行，其余类此。类似于多利特尔分解法，得到：对k=1，2，…，n返回(5.26)第一百一十九页，共一百八十页，编辑于2023年，星期二

实现A=LU，则Ax=b可以化为LUx=b。令Ux=y，则Ly=b。由Ly=b解出yi:

再由Ux=y解出xi:

第一百二十页，共一百八十页，编辑于2023年，星期二注*：为保证LU分解能够进行，要求，既要求A的

所有顺序主子式不为零，而线性方程组有解，只须|A|≠0

即可。为取消此限制，采用类似列主元消元法处理。克洛特Crout分解列主元直接分解算法*输入：n，A，b计算L的第k列元素和U的第k行元素第一百二十一页，共一百八十页，编辑于2023年，星期二第一百二十二页，共一百八十页，编辑于2023年，星期二例4*

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

设A=LU，即

第一百二十三页，共一百八十页，编辑于2023年，星期二解下三角方程组Ly=b，即解(单位)上三角方程组Ux=y，即第一百二十四页，共一百八十页，编辑于2023年，星期二★克洛特Crout分解算法描述第一百二十五页，共一百八十页，编辑于2023年，星期二由Ly=b解出yi:

再由Ux=y解出xi:第一百二十六页，共一百八十页，编辑于2023年，星期二程序源码主要部分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]); }}第一百二十七页，共一百八十页，编辑于2023年，星期二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];}第一百二十八页，共一百八十页，编辑于2023年，星期二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];}第一百二十九页，共一百八十页，编辑于2023年，星期二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]);}第一百三十页，共一百八十页，编辑于2023年，星期二5.3.2平方根法第一百三十一页，共一百八十页，编辑于2023年，星期二5．A对称正定，则A的对角元素aii>0，i=1,2,…,n。第一百三十二页，共一百八十页，编辑于2023年，星期二第一百三十三页，共一百八十页，编辑于2023年，星期二定理1若A为对称矩阵，且A的各阶顺序主子式不为零，则A可唯一分解成第一百三十四页，共一百八十页，编辑于2023年，星期二证明第一百三十五页，共一百八十页，编辑于2023年，星期二引理

如果定理1中的A正定，则对角矩阵D中的元素非负。定理2(Cholesky分解)

若A为对称正定矩阵，则存在下三角（乔列斯基）矩阵U，使A=UUT。（平方根法）证明：因为A对称，由定理1知，其中L为单位下三角阵，D为对角阵。如果A正定，则由引理知，D中的元素非负。记

第一百三十六页，共一百八十页，编辑于2023年，星期二下面我们用直接分解法来确定计算L元素的递推公式：几点说明：(1)L

为对角线元素全为正的下三角矩阵；(2)只有对称正定矩阵(SPD)才存在Cholesky分解；(3)SPD矩阵的Cholesky分解存在且唯一将展开第一百三十七页，共一百八十页，编辑于2023年，星期二-------------(1)-------------(2)-------------(3)第一百三十八页，共一百八十页，编辑于2023年，星期二第一百三十九页，共一百八十页，编辑于2023年，星期二紧凑格式乔列斯基(Cholesky)分解A=LLT总结：1.先确定L中第一列元素：2.确定L中对角线元素:4.确定L中第i行元素:(r=2,3,…,i-1)第一百四十页，共一百八十页，编辑于2023年，星期二第一百四十一页，共一百八十页，编辑于2023年，星期二第一百四十二页，共一百八十页，编辑于2023年，星期二紧凑格式解法：第一百四十三页，共一百八十页，编辑于2023年，星期二第一百四十四页，共一百八十页，编辑于2023年，星期二对称正定线性方程组的解法线性方程组-------------(5)-------------(6)则线性方程组(10)可化为两个三角形方程组-------------(7)-------------(8)第一百四十五页，共一百八十页，编辑于2023年，星期二------(9)------(10)对称正定方程组的平方根法第一百四十六页，共一百八十页，编辑于2023年，星期二例1.用平方根法解对称正定方程组解:第一百四十七页，共一百八十页，编辑于2023年，星期二第一百四十八页，共一百八十页，编辑于2023年，星期二即所以原方程组的解为第一百四十九页，共一百八十页，编辑于2023年，星期二第一百五十页，共一百八十页，编辑于2023年，星期二第一百五十一页，共一百八十页，编辑于2023年，星期二平方根法的数值稳定性用平方根法求解对称正定方程组时不需选取主元由可知因此平方根法是数值稳定的事实上,对称正定方程组也可以用顺序Gauss消去法求解而不必加入选主元步骤第一百五十二页，共一百八十页，编辑于2023年，星期二改进的平方根法5.325.33第一百五十三页，共一百八十页，编辑于2023年，星期二例如：

用LDLT解方程组解第一百五十四页，共一百八十页，编辑于2023年，星期二于是，Ax=bLDLTx=b.解方程组有：第一百五十五页，共一百八十页，编辑于2023年，星期二第一百五十六页，共一百八十页，编辑于2023年，星期二第一百五十七页，共一百