解线性方程组的直接法

解线性方程组的直接法

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

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

则说矩阵是可逆的，并把矩阵称为的逆矩阵.使得第4页,共180页，2024年2月25日，星期天行列式性质(6)行列式按任一行展开其中Ai,j为ai,j的代数余子式，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用不等于零的数乘某行，加至另一行。第17页,共180页，2024年2月25日，星期天第18页,共180页，2024年2月25日，星期天第19页,共180页，2024年2月25日，星期天高斯消去法的主要思路：将系数矩阵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)放入一个二维数组.消去每一步算出的aij(k)放入aij的位置,bi(k)放入bi,mik放在aik上.消去完毕后,数组各位置上的数据如下示:回代过程的计算没有数据组织问题.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(新)说明:系数矩阵存放于数组A中,右端向量放于数组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.

即使高斯消元法可行，如果很小，运算中用它作分母会导致其它元素数量计的严重增长和舍入误差的扩散。第30页,共180页，2024年2月25日，星期天顺序高斯消去法的使用条件

使用条件之一

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

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

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

使用条件之二

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

定理

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

原因:求行乘数时用了”小主元”0.0001作除数.第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位有效数字。解(1)顺序消元：

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

前例方法2中所用的是在Gauss消去法中加入选主元过程,利用换行,避免”小主元”作除数,该方法称为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,列主元消去法能避免”小主元”作除数,误差分析与数值试验表明：（Gauss）列主元消去法”基本稳定”.第39页,共180页，2024年2月25日，星期天列主元消去法中方程对换问题引入向量P=(p1,p2,…,pn)t=(1,2,…,n)t即i

pi记录方程(5.1)的行号.当第k步消元若需第k个方程与第ik个方程对换,则只需交换P中当消元过程完后,向量P中的pi(i=1,…,n)指示初始方程(5.1)在作的i步消元时被选为主元的方程序号.这种序号交换避免了方程间的对换,n较大时,有其优越性.第40页,共180页，2024年2月25日，星期天例

用列主元高斯消去法求解方程组(用三位有效数字计算)解选主元选主元第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说明:系数矩阵存放于数组A中,右端向量放于数组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

的顺序，须记录交换次序，解完后再换回来。全主元高斯消去法具有很好的稳定性，但选全主元比较费时，故在实际计算中很少使用。第49页,共180页，2024年2月25日，星期天程序示例★列主元高斯消元法算法描述：

1.输入数据：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消去法是将系数矩阵化为上三角矩阵，再进行回代求解；Jordan法是将系数矩阵化为对角矩阵，从而无须再进行回代求解的一种直接解法。这种解法基本过程与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消去法可以执行）的条件是什么？第82页,共180页，2024年2月25日，星期天Gauss消去法可以执行定理

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

(k=1,2,…,n-1)，则第83页,共180页，2024年2月25日，星期天由上述对系数矩阵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为单位上三角矩阵。第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日，星期天对于线性方程组系数矩阵非奇异,经过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的计算过程发现:算法的数据组织:对于线性方程，

第103页,共180页，2024年2月25日，星期天Doolittle分解法的数据组织可以用以下紧凑格式表示:第104页,共180页，2024年2月25日，星期天UL第105页,共180页，2024年2月25日，星期天

Doolittle法解本节方程组(例1)的数据组织过程：解:第106页,共180页，2024年2月25日，星期天第107页,共180页，2024年2月25日，星期天所以x第108页,共180页，2024年2月25日，星期天例4用杜利特尔分解法求解方程组解先求增广系数矩阵的杜利特尔分解，即第109页,共180页，2024年2月25日，星期天再一例：第110页,共180页，2024年2月25日，星期天第111页,共180页，2024年2月25日，星期天开始输出错误信息停机1.Doolittle法解方程组流程图程序实现计算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.本实验是后续改进算法(平方根法、Cholesky分解法)的基础。对于大型方程组，一般采用迭代算法（下一章）.第114页,共180页，2024年2月25日，星期天问题:在Doolittle法(包括LU分解紧凑格式)中,会反复用到公式仍有可能为小主元做除数?为此,也应考虑在算法中加入选取列主元即紧凑格式的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日，星期天

实现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日，星期天平方根法的数值稳定性用平方根法求解对称正定方程组时不需选取主元由可知因此平方根法是数值稳定的事实上,对称正定方程组也可以用顺序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月