第二章解线性方程组的直接方法2017秋

1高斯消去法主元素法直接三角分解法平方根法与改进平方根法误差分析第二章线性方程组的直接方法2讨论n元线性方程组的直接解法.3其中方程组的矩阵形式为若矩阵A非奇异,即det(A)≠0,则方程组有唯一解.4直接解法是指,若不考虑计算过程中的舍入误差，经过有限次算术运算就能求出线性方程组的精确解的方法.但由于实际计算中舍入误差的存在,用直接解法一般也只能求出方程组的近似解.Cramer法则是一种不实用的直接法，下面介绍几种实用的直接法.解线性方程组的方法:直接方法和迭代法.迭代法是从解的某个近似值出发,通过构造一个无穷序列去逼近精确解的方法.一般地,有限步内得不到精确解.5Gauss消元法是一种规则化的消元法，其基本思想是通过逐次消元计算，把一般线性方程组的求解转化为等价的上三角形方程组的求解。§1Gauss消去法6消去后两个方程中的x1得再消去最后一个方程的x2得消元结束.经过回代得解:例考虑线性方程组顺序Gauss消去法7消元过程:先逐次消去变量x1,x2,将方程组化为同解的上三角形方程组.上过方法可推广到一般情况回代过程:按方程相反的顺序求解上三角形方程组.8第一步.设依次用乘矩阵的第1行加到第i行上,得到矩阵:则线性方程组的增广矩阵为记9其中第二步.设依次用乘矩阵的第2行加到第i行,得到矩阵:10其中11这就完成了消元过程.如此继续消元下去,第n1步结束后得到矩阵:12对此方程组进行回代，就可求出方程组的解.对应的方程组变成：13能用顺序Gauss消去法求解的条件是在消元过程中得到的主元必须全不为0，即顺序Gauss消去法通常也简称为Gauss消去法.主元素都不为零矩阵A的各阶顺序主子式都不为零.

顺序Gauss消去法中的称为主元素.14顺序Gauss消去法求解n元线性方程组的乘除运算量

第1次消元乘除运算量:

消元过程乘除运算量求li1:

(n－1)求aij(2):(n－1)2求bi(2):(n－1)共(n2－1)次15

回代过程乘除运算量：求xn:

1求xn－1:2求x1:n….16n=30时,顺序Gauss消去法只需9890次乘除法运算.顺序Gauss消去法求解n元线性方程组的乘除运算量17高斯消去法优缺点：简单易行要求主元均不为零,因而适用范围小数值稳定性差18例：单精度解方程组

精确解8个8个用顺序Gauss消去法计算：8个小主元可能导致计算失败.§2主元素法19若将方程组改写成:用顺序Gauss消去法,消元得回代得解:x2=1,x1=1与准确解非常接近.可见,第一种算法是不稳定的,第二种算法是稳定的.20此例说明,在消元过程中,应避免选取绝对值较小的数作主元.如例中的第二种解法,通过交换方程次序,选取绝对值较大的元素作为主元.基于这种想法导出了主元法.为了提高计算的数值稳定性,在消元过程中采用选择主元的方法.常采用的是列主元消去法和全主元消去法.21给定线性方程组Ax=b,记A(1)=A,b(1)=b,列主元Gauss消去法的具体过程如下:

首先在增广矩阵B(1)=(A(1),b(1))的第一列元素中,取

然后进行第一步消元得增广矩阵B(2)=(A(2),b(2)).列主元消去法22

然后进行第二步消元得增广矩阵B(3)=(A(3),b(3)).按此方法继续进行下去,经过n1步选主元和消元运算,得到增广矩阵B(n)=(A(n),b(n)).则方程组A(n)x=b(n)是与原方程组等价的上三角形方程组,可进行回代求解.只要|A|0,列主元Gauss消去法就可顺利进行.

再在矩阵B(2)=(A(2),b(2))的第二列元素中,取23

全主元素法每一步选绝对值最大的元素为主元素，保证Stepk:①选取②Ifik

k

then交换第k行与第ik

行;Ifjk

k

then交换第k列与第jk

列;③消元注：列交换改变了xi

的顺序，须记录交换次序，解完后再换回来。24例用主元素法求解线性方程组计算过程保留三位小数,方程的精确解为x1*=1,x2*=2,x3*=3.25解1.按列主元素法，求解过程如下消元回代得26解2.按全主元素法，求解过程如下回代得27全主元素法的精度优于列主元素法,这是由于全主元素是在全体系数中选主元,故它对控制舍入误差十分有效.但全主元素法在计算过程中,需同时作行与列的互换,因而程序比较复杂,计算时间较长.

列主元素法的精度虽然稍低于全主元素法,但其计算简单,工作量大为减少,且计算经验与理论实践均表明,它与全主元素法同样具有良好的数值稳定性.列主元素法是求解中小型稠密线性方程组的最好方法之一.

例3的计算结果表明28Gauss消元法的矩阵表示§3直接三角分解法两者等价29其中两者等价n=3时Gauss消元法的矩阵表示30其中两者等价3132条件:矩阵A经Gauss消元法后得到的上三角矩阵.33例求矩阵的三角分解.解34上述n=3的情形可以推广到一般情形,例如n=43536条件:矩阵A经Gauss消元法后得到的上三角矩阵.37例38条件:矩阵A经Gauss消元法后得到的上三角矩阵.Gauss消元法的矩阵表示Doolittle分解39矩阵的三角分解定理

设A为n阶方阵,若A

的前n阶顺序主子式Ai

(i=1,2,…,n)均不为0,则矩阵A存在唯一的Doolittle分解.存在性：唯一性：反证法40

Doolittle分解中LU

元素的求解次序Doolittle分解中LU

的另一求法41U的第一行L的第一列

Doolittle分解中LU

元素的求解右端用矩阵乘法展开,比较两边的第1行和第1列得42假设已求得U的前(r

1)行和L的前(r

1)列,r>1.下面求U的第r行和L的第r列.右端用矩阵乘法,比较两边的第r行的后(n－r+1)个元素arj

=L的第r行行向量与U的第j列列向量的内积(j

r)U的第r行43U的第r行44右端用矩阵乘法,比较两边的第r列的后(n－r)个元素air=L的第i行行向量与U的第r列列向量的内积(i

>r)L的第r列45L的第r列46对r=2,3,…,n,矩阵三角分解的紧凑格式47例用紧凑格式求矩阵的三角分解.解48先求Ly=b,得y;再求Ux=y,得x.直接三角分解法或Doolittle分解法.Doolittle分解法49Ly=b乘除运算量为50Ux=y乘除运算量为51例用直接三角分解法解解52先求Ly=b得y再求Ux=y

得x53由于在求出uij,lij和yi后,aij和bi就无需保留了,故上机计算时,可把L,U和y存在A,b所占单元,回代时x取代y,整个计算过程中不需要增加新的存储单元.在求一系列系数矩阵相同而右端项不同的线性方程组Ax=b(k),(k=1,2,…,m)时(如求逆矩阵),用三角分解法更为简便.每解一个方程组Ax=b(k)

仅需要增加n2

次乘除法运算.解线性方程组Ax=b的Doolittle三角分解法的计算量约为n3/3,与Gauss消去法相同.54Crout分解定理

设A为n阶方阵,若A

的前n阶顺序主子式Ai

(i=1,2,…,n)均不为0,则矩阵A可以唯一分解为其中L为下三角阵,U为单位上三角阵.55特殊的稀疏矩阵解三对角方程组的追赶法解三对角方程组56追赶法是求三对角线性方程组的三角分解法.

追赶法本质上是Gauss消元法.二阶常微分方程边值问题的差分离散方程组,热传导方程以及船体数学中建立的三次样条函数的三转角或三弯矩方程组均为三对角方程组.57例

4阶三对角矩阵的三角分解58例

4阶三对角矩阵的三角分解单位下二对角阵上二对角阵59定理设三对角矩阵A满足下列条件则它可以分解成A=LU三对角矩阵的三角分解对角占优阵单位下二对角阵上二对角阵工程中得到的三对角阵多数满足此条件60三对角矩阵三角分解中LU的求解次序61对k=3,…,n

三对角矩阵三角分解中LU的求解右端用矩阵乘法展开,比较两边元素得乘除运算量2n－262解三对角方程组的追赶法先求Ly=d,得y;再求Ux=y,得x.追：消元过程赶：回代过程63Ly=d乘除运算量为n－164Ux=y乘除运算量为2n－1追赶法总的乘除运算量5n－465追赶法的实质就是Gauss消元法,只是由于系数中出现了大量的零,在计算过程中将它们撇开,从而使计算公式大大简化,也大大减少了计算量.为节省计算机存储单元,计算得到的lk,uk分别存放在ak,bk的存储单元内,而yk,xk

存放在dk

的存储单元内.当系数矩阵为满足定理条件的严格对角占优阵时,追赶法具有良好的数值稳定性.66§4平方根法与改进的平方根法求解对称正定方程组,即其中A为对称正定阵.利用对称性,平方根法的计算量是Gauss消去法的一半.67对称正定阵的Cholesky分解定理

若A

是对称正定阵,则存在唯一的非奇异下三角阵L,使得且L的对角元素皆为正数,即lii>0(i=1,2,…,n).证明

A可进行Doolittle分解.A为对称正定阵,它的n个顺序主子式均大于0,设其中为单位下三角阵,U为上三角矩阵.用反证法68对U进一步分解.

A对称

A正定对角矩阵D的对角元均为正.令则其中为非奇异下三角阵,且对角元均为正.

A=LLT.69对称

A为3阶对称正定阵,A=LLT,怎样求L?对称70例对下列矩阵进行Cholesky分解71例对下列矩阵进行Cholesky分解Matlab解法：A=[4-26;-2175;6522]R=chol(A)L=R’72

A为n阶对称正定阵,A=LLT,怎样求L?L中元素的求解次序依次求L的第一列,第二列,…,第n列.73………...对称

A为n阶对称正定阵,A=LLT,怎样求L?li为L的第i

个行向量74对称

A为n阶对称正定阵,A=LLT,怎样求L?75求L的第一列求L的第二列

A为n阶对称正定阵,A=LLT,怎样求L?76

A为n阶对称正定阵,A=LLT,怎样求L?设已经求得L的前k－1列,现求L的第k列(k=3,4,…,n)77例求正定阵的Cholesky分解.解78例求正定阵的Cholesky分解.解79先求Ly=b,得y,再求LTx=y,得x.解正定线性方程组的平方根法或Cholesky分解法.平方根法或Cholesky分解法设A为对称正定阵80定理

若A

是对称正定阵,则存在唯一的单位下三角阵L和对角阵D,使得且D的对角元素皆为正数.对称正定阵的LDLT分解证明

A对称

A正定对角矩阵D的对角元均为正.81对称正定阵的LDLT分解本质上是对A作Doolittle分解,即LU分解.LDLT分解中的D=LU分解中的U的对角部分LDLT分解中的L=LU分解中的L82对称正定矩阵A的LU分解,计算量可以节省一半求U的第1行求L的第1列对称正定阵的LDLT分解中L,D的计算先对对称正定阵A作LU分解83求U的第k行(k=2,3,…,n)求L的第k列(k=2,3,…,n)对称正定阵的LDLT分解中L,D的计算节省了计算量求D84例求矩阵的LDLT分解.解8586解正定线性方程组的改进平方根法或LDLT分解法.先求Ly=b,得y,再求LTx=D1y,得x.改进平方根法或LDLT分解法设A为对称正定阵87平方根法与改进的平方根法的优点计算无须选主元,由于正定性,计算过程是数值稳定的计算量是Gauss消元法的一半由于对称性,实际计算可存储一半是求解中小型稠密正定线性方程组的好算法88§5误差分析用直接法解线性方程组，初始数据会有误差，计算过程同样会产生误差，这就需要对这些误差作一些分析.主要数学工具：向量范数，矩阵范数，条件数89向量范数向量范数是用来度量向量长度的,它可以看成是二、三维解析几何中向量长度概念的推广.定义(向量范数)对任一向量xRn，按照一定规则确定一个实数与它对应，该实数记为||x||，若||x||满足下面三个性质：1)||x||0；||x||=0当且仅当x=0(零向量)

正定性2)对任意实数，||x||=||||x||齐次性3)对任意向量x,yRn，||x+y||||x||+||y||三角不等式则称该实数||x||为向量x的范数.90Rn中常用的三种范数其中x1,x2,…,xn分别是x的n个分量.1-范数2-范数-范数用向量范数的定义来验证,即验证它们满足向量范数定义中的三个性质.向量的模91定义(范数等价)

在Rn中有两个范数||||和||||,若存在实数M,m>0,使得对任意的n维向量x,都有则称这两个范数等价.Rn中范数的重要性质:范数等价定理92范数等价定理:

Rn中任意两个范数等价.Rn中范数的重要性质:范数等价定理例1-范数,2-范数和-范数是两两等价的.93当不需要指明使用哪一种向量范数时，就用记号||.||泛指任何一种向量范数.有了向量的范数就可以用它来衡量向量的大小和表示向量的误差.设x

为Ax=b

的精确解，x*为其近似解绝对误差相对误差94矩阵范数是用于定义矩阵“大小”的量，类似于向量范数，可以定义n阶方阵A的范数.定义(矩阵范数)

设A为n阶方阵，按照一定规则有一实数与之对应，记为||A||,若||A||满足：1)||A||0,||A||=0当且仅当A=0时；2)对任意实数，||A||=||||A||；3)对任意两个n阶方阵A,B,都有||A+B||||A||+||B||;4)||AB||||A||||B||(相容性条件)则称||A||为矩阵A的范数.矩阵范数95常用的三种矩阵范数1-范数或列范数2-范数或谱范数-范数或行范数用矩阵范数的定义来验证,即验证它们满足矩阵范数定义中的四个性质.96定理设A为n阶方阵，||||是Rn中的向量范数,则是一种矩阵范数,称其为由向量范数||||诱导出的矩阵范数.由向量范数诱导出的矩阵范数矩阵范数与向量范数的相容性性质：设矩阵范数是由向量范数诱导出的，则对任意n阶方阵A,以及任意的n维向量x,有979899常用的三种矩阵范数均是由向量范数诱导出的.

对于给定的向量范数1-范数,2-范数及-范数,可以证明由它们诱导出的矩阵范数分别为1-范数或列范数2-范数或谱范数-范数或行范数100101实际计算常用1-范数与-范数因其计算比较简单,理论证明常用2-范数因为它有一些好性质.由矩阵范数与向量范数的相容性性质知对任意n阶方阵A,以及任意的n维向量x,有误差分析中常用102矩阵的误差可用矩阵范数表示.

设A*是A的近似矩阵绝对误差相对误差矩阵范数的等价定理也成立.103Matlab函数：normNORM(X)isthelargestsingularvalueofX,max(svd(X)).NORM(X,2)isthesameasNORM(X).NORM(X,1)isthe1-normofX,thelargestcolumnsum,=max(sum(abs((X)))).NORM(X,inf)istheinfinitynormofX,thelargestrowsum,=max(sum(abs((X')))).NORM(X,'fro')istheFrobeniusnorm,sqrt(sum(diag(X'*X))).NORM(X,P)isavailableformatrixXonlyifPis1,2,infor'fro'.

104NormofAMatrixMatlabcommand:norm(A,1)norm(A,2)=norm(A)norm(A,inf)norm(A,’fro’)105方程组的状态与条件数例方程组I方程组II右端项有0.00001的差别,最大相对误差为0.5105,但解分量的相对误差为50%.平面上两条接近于平行的直线的交点，当其中一条直线稍有变化时,新的交点可与原交点相差甚远.106当一个方程组,由于系数矩阵或右端项有微小扰动,而引起解发生巨大变化时,则称该方程组是“病态”的.分以下两种情形加以讨论只有右端项有扰动只有系数矩阵有扰动怎样刻画线性方程组“病态”程度.用条件数来描绘.107只有右端项有扰动原方程组扰动后方程组右端项有扰动b,引起解的变化为x,其相对误差为它究竟有多大?108故又这表明:当右端项有扰动b时,解的相对误差不超过右端项的相对误差的||A||||A1||倍.只有右端项有扰动从而109只有系数矩阵有扰动原方程组扰动后方程组系数矩阵有扰动A,引起解的变化为x,其相对误差为它究竟有多大?110故这表明:当系数矩阵有扰动A,解的扰动不超过系数矩阵相对误差的||A|