《线性代数计算方法.ppt》由会员分享，可在线阅读，更多相关《线性代数计算方法.ppt（128页珍藏版）》请在三一办公上搜索。
1、第3章 线性代数计算方法,1 高斯消去法,3 解实三对角线性方程组的追赶法,4 矩阵的三角分解,5 行列式和逆矩阵的计算,7 迭代法的收敛性,8 矩阵的特征值与特征向量的计算,2 高斯约当消去法,6 迭代法,在自然科学和工程技术中很多问题的解决常常归结为解线性代数方程组。例如:电学中的网络问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题等,直接法:经过有限次运算后可求得方程组精确解的方法(不计舍入误差!)迭代法：从解的某个近似值出发，通过构造一个无穷序列去逼近精确解的方法(一般有限步内得不到精确解),解线性方程组的两类方法:,1 高斯消去法,一、高斯消去法,且ai
2、i0,i=1,2,n,上三角阵,A,X,B,下三角阵?,化成上(下)三角!,怎么解?,高斯消去法,-2,-3,-2,例：,第2行：计算比例因子,（1）消去过程：,第一步：消，设,第2 行-l21 第1行，得到：,消去法的数值计算过程：,消去,-1,-1,第 i 行-li1 第1行，得到：,第 i 行：计算因子,第2行：,第 n 行：消去,第k步：消去,共进行 n 1步，得到,设，计算因子且计算,消元过程总体流程：,对于,对于,做,做,对于,做,（2）回代过程：,二、选主元消去法,为避免这种情况的发生，可通过交换方程的次序，选取绝对值大的元素作主元。基于这种思想导出了主元素法,在高斯消去法消去过
3、程中可能出现 的情 况，这时高斯消去法将无法进行；即使主因素 但很小，其作除数，也会导 致其它元素数量级的严重增长和舍误差的扩散。,例如：用高斯消去法求解下列方程组(用四位有效数字计算):,x2=0.5999959999,化简可得 x2=0.6000回代求得 x1=105(0.6-0.6000)=0而方程组的解应为 x1=0.4000 x2=0.6000,显然用上述方法求出的解x1与方程组的实际解相差很大。若改变两个方程的顺序,即 x1+x2=1 10-5 x1+x2=0.6-10-5得(1.000-1.00010-5)x2=0.6-1.00010-5 0.99999x2=0.59999 化简
4、得 x2=0.6000 回代求得 x1=(1-0.6000)=0.4000,x2=0.5999959999,高斯主元素消去法是顺序消去法的一种改进。它的基本思想是在逐次消元时总是选绝对值最大的元素(称之为主元)做除数，按顺序消去法的步骤消元。这里主要介绍求解线性方程组最常用的列主元素消去法和全主元素消去法。,所谓列主元素消去法就是在每一步消元过程中取系数子矩阵的第一列元素中绝对值最大者作主元。对线性方程组进行n-1次消元后，可得到上三角形方程组,列主元消去法,这种方法称为列主元Gauss消去法。,取四位有效数字计算。解 中-18为主元,交换和得,例1 用列主元素消去法解方程组,+12/18,+
5、1/18得,第二列消元时,主元为1.167,交换方程和得,+1/1167得,回代求得 x1=1.000,x2=2.000,x3=3.001方程组的实际解 x1=1,x2=2,x3=3,开始,F,T,F,T,找列最大值,当前列,换行,消元,回代、求解,找列最大值,换行,对于,对于,做,做,对于,做,/*本算法用高斯列主元素消去法求解矩阵方程AX=B。其中：A是NN矩阵；B是N1矩阵。输入：nA的行数；a二维矩阵A b矩阵B算法结束后，函数返回值为ERROR_CODE时，表示A是奇异的或病态的；否则，A代表行列式的值。a A消元后的上三角矩阵 b矩阵方程的解X*/,double Gaussian_
6、elimination(int n,double ann，double bn)int i,j,k,mk;double mm,f;for(k=0;kn-1;k+)mm=akk;mk=k;for(i=k+1;in;i+)if(fabs(mm)fabs(aik)mm=aik;mk=i;if(fabs(mm)=0)return(ERROR_CODE);,当前列,找列最大值,if(mk!=k)for(j=k;jn;j+)f=akj;akj=amkj;amkj=f;f=bk;bk=bmk;bmk=f;for(i=k+1;in;i+)mm=aik/akk;aik=0.0;for(j=k+1;jn;j+)ai
7、j=aij mm*aki;bi=bi mm*bk;,换行,消元,bn-1=bn-1/an-1 n-1;for(j=n-2;j=0;j-)for(k=j+1;kn;k+)bj=bj ajk*bk bj=bj/ajj;,回代、求解,所谓全主元素消去法，就是每步消元时选取系数子矩阵中绝对值最大的元素作主元。经过n-1次消元后,方程组可化为上三角形方程组,全主元消去法,例2 用全主元素消去法求解方程组,解这里主元为-18,交换方程与方程得,再全选主元，主元为2.333，交换x2和x3所在的两列，同时改变两未知数的排列号得：,-0.944/2.333得,已经化为三角方程组，回代求解：x1=1.000，x
8、2=3.000，x3=2.000 这里未知数x2与x3已对调，所以应恢复解的顺序，方程组的实际精确解为：x1=1.000，x2=2.000，x3=3.000,在全主元素的消元过程中，其系数矩阵可能进行了列对调，那么未知数也相应地作了对调。在得到计算结果之后，需要恢复未知量的原顺序。,值得注意的是：,图 3.2,找各子系数矩阵最大值,换行、换列,图 3.2,恢复解的顺序,2 高斯约当消去法,前面所述的消去法均要进行两个过程，即消元过程和回代过程。但对消元过程稍加改变可以把方程组化为对角形：,此时求解就不要回代了。这种无回代过程的主元素消去法称为 高斯约当(Jordan)消去法。特别是方程组还可化
9、为,显然等号右端即为方程组的解。对于n阶线性方程组，其增广矩阵为,首先把主元素(按列选主元或全选主元)调换到主对角线上，并化为1，再将主元素所在列的其它元素消为0，则第一次消元后增广矩阵化为,显见经n步消元后即得方程组的解。其具体计算步骤为:对 k=1,2,n 选主元(列选或全选主元)对 j=k,k+1,n+1计算,计算 还得指出,若用到全选主元,最后应注意恢复解的顺序。,对于,做,对于,做,4 矩阵的三角分解,复习,矩阵的初等变换,复习,矩阵的初等变换,4 矩阵的三角分解,i,j,复习,矩阵的初等变换,4 矩阵的三角分解,4 矩阵的三角分解,4 矩阵的三角分解,y,定理：若A的所有顺序主子式
10、 均不为0，则 A 的 LU 分解唯一（其中：L 或U为单位三角阵）。,下三角矩阵,上三角矩阵,（1）L 为单位下三角阵而 U 为一般上三角阵的分解称为Doolittle 分解（2）L 为一般下三角阵而 U 为单位上三角阵的分解称为Crout 分解。,Doolittle分解法：,例：用杜利特尔分解方法求解下列方程组的解:,1,2,3,4,5,Doolittle分解,L为下三角，U为单位上三角,比较第1行：,比较第1列：,计算L的各列与U的各行的次序如下图所示：,L为下三角，U为单位上三角,Courant分解法：,用克劳特分解方法求解下列方程组,解 令,利用矩阵乘法可得到,这样原方程组就化为依次
11、求下列两个三角形方程组,代入第二个方程组可求得原方程组的解为,这样,L、U中的元素都已求出。计算L的各列与U的各行的次序如图3.4所示。,图 3.4,用克劳特分解求解线性方程组的计算过程为:（1）LU分解过程:对于k=1,2,n依次计算,（2）求解过程:,求 l,求 u,求 y,图 3.5,求 x,从上述计算过程还可看出，求解yk的公式完全可以插在计算L、U元素的两个公式之间。,对方程组,做等价变换,与非线性方程的迭代方法一样，需要我们构造一个等价的方程，从而构造一个收敛序列，序列的极限值就是方程组的根,6 迭代法,2 赛德尔（Seidel）迭代,1 雅可比（Jacobi）迭代,1 雅可比Ja
12、cobi迭代,格式很简单：,例8 用简单迭代法解下列方程组,解将方程组写成等价形式,取初始值x(0)=0,按迭代公式,表 31,图 3.8,Jacobi迭代收敛条件,对于Jacobi迭代，我们有一些保证收敛的充分条件,定理：若M满足下列条件之一，则Jacobi迭代收敛。,M为行对角占优阵,M为列对角占优阵,M满足,2 赛德尔Seidel迭代,在Jacobi迭代中，使用最新计算出的分量值,例9 用赛德尔迭代法解方程组,解 将原方程组写成等价形式并按(375)构造赛德尔迭代公式,表 32,Jacobi迭代算法,1、输入系数矩阵M和向量b，和误差控制eps2、x1=0,0,.,0,x2=1,1,.,
13、1/赋初值3、while(|M*x2-b|eps)x1=x2;for(i=0;in;i+)x2i=0;for(j=0;ji;j+)x2i+=Mij*x1j for(j=i+1;jn;j+)x2i+=Mij*x1j x2i=(bi-x2i)/Mii 4、输出解x2,前次计算结果,当前计算结果,原系数矩阵,Siedel迭代算法,1、输入系数矩阵M和向量b，和误差控制eps2、x2=1,1,.,1/赋初值3、while(|M*x2-b|eps)for(i=0;in;i+)for(j=0;ji;j+)x2i+=Mij*x2j for(j=i+1;jn;j+)x2i+=Mij*x2j x2i=(bi-x
14、2i)/Mii 4、输出解x2,Siedel迭代收敛条件,定理：若A满足下列条件之一，则Seidel迭代收敛。,A为行或列对角占优阵,A对称正定阵,注：二种方法都存在收敛性问题。有例子表明：Seidel法收敛时，Jacobi法可能不收敛；而Jacobi法收敛时，Seidel法也可能不收敛。,7.1 向量范数 在三维空间中,常用三种方法来衡量一个向量r=(x,y,z)T的“长度”,即,7 迭代法的收敛性,定义1：设n维向量xR，记对应向量x的一个实数为x,若x满足下面三个性质:(1)非负性，即x0，当且仅当 x=0时，x=0(2)齐次性，x=x,为任意实数(3)三角不等式性，x+yx+y，yRn
15、。则称该实数x为向量x的范数。可验证，对于不小于1的正数p,n,具有向量范数的三条基本性质。,称上式为向量x的p范数,记为xp，即,在Rn中，常用的几种向量范数有:,定义2 设A为nn阶矩阵，定义,为矩阵A的范数。可以证明，矩阵范数满足下列的几个性质:(1)非负性，A0,当且仅当A=0时,A=0(2)齐次性，A=A,为任意实数(3)三角不等式性,A+BA+B,B与A为同阶矩阵(4)AxAx(5)ABAB,常见的矩阵范数,是 的最大特征值,谱半径：设nn阶矩阵A的特征值为i(i=1,2,n),则称,为矩阵A的谱半径。,定理3 设(A)是矩阵A的谱半径，则对于A的任何一种范数A,都有(A)A,定理
16、5 对于任意初始向量x(0)和常数项f,由迭代公式 x(k+1)=Ax(k)+f,k=0,1,2,收敛的充要条件是(A)1,通常(A)1很难证明，一般用 A 1来判断迭代是否收敛,8 矩阵的特征值与特征向量的计算,求矩阵的特征值与特征向量是自然科学许多分支中一个重要的问题。nn阶矩阵A的特征值是其特征方程 det(A-I)=0 的根,而属于特征值的特征向量是线性方程组(A-I)x=0 的非零解。,8.1 乘幂法 在实际问题中,往往不需要求矩阵A的全部特征值,而只要计算其按模最大或最小的特征值,也即只需求矩阵A的部分特征值。乘幂法便是一种行之有效的方法。设n阶实矩阵A有n个线性无关的特征向量,其
17、特征值按模排列如下:12n 其相应的特征向量分别记为 x1,x2,xn 现求按模最大的特征值1及其相应的特征向量x1。,取非零初始向量v0,作迭代构成一向量序列,因为特征向量x1,x2,xn线性无关,故v0可以表示成这n个线性无关的向量的线性组合,即 v0=1x1+2x2+nxn,于是 vk=Ak(1x1+2x2+nxn)=11kx1+22kx2+nnkxn 设1为实数,将上式改写成,若12,由上式有：,因此，当k充分大时，只要10,就有 vk1k1x1(1)(A-1I)vk11k(A-1I)x1=0 这说明vk是相应于特征值1的特征向量的近似向量。由式(1)，当k充分大时：v k+11k+1
18、1x1=1k1x11vk,例10 求矩阵,的按模最大的特征值与相应的特征向量。,(A6v0)1/(A5v0)13.4(A6v0)2/(A5v0)23.4(A6v0)3/(A5v0)33.4(A7v0)1/(A6v0)13.4(A7v0)2/(A6v0)23.4(A7v0)3/(A6v0)33.4,于是有 13.4相应的特征向量为 x1(792,-1120,792),若把x1的第一个分量化为1,则有 x1(1,-1.4,1)事实上,矩阵A的特征值与相应的特征向量为,习题：,要求:习题三 1.2，6请尽力试着编写程序,非线性方程组的迭代法,例 解下列非线性方程组,k,k,直接推广Newton迭代为
19、：,单个方程的牛顿迭代形式:,雅可比矩阵(Jacobi),Newton方法的迭代格式的推导,解非线性方程组的主要步骤：,(1)选定初值X0,(4)按公式计算X1,(2)计算F(X0),J(X0),(3)计算X0,(5)误差判断,例：,例 用牛顿法解下列非线性方程组,k,k,SolveNonlinear,程序框架,GetFun,agaus,GetDeltaX,GetJacobi,UpdateX,Main函数,。,。,。,。,。,。,#include stdlib.h#include math.h#include stdio.h#define N 3void SolveNonlinear(doub
20、le xN,double eps,int k);void GetFun(double xN,double bN);void GetJacobi(double xN,double JNN);double GetDeltaX(double JNN,double bN,double deltaX N);int agaus(double a NN,double bN,int n);void UpdateX(double xN,double deltaX N);,计算F(X0),计算F(X0),解非线性方程组主程序,计算X0主程序,更新计算结果（相当于计算X1）,非线性方程组求解程序（牛顿法）,/主函数
21、main()int i,k;double eps;double x N=1.0,1.0,1.0;double b N;eps=0.0000001;k=100;SolveNonlinear(x,eps,k);GetFun(x,b);printf(n);printf(n);for(i=0;iN;i+)printf(x(%d)=%13.7e,F(x)=%13.7en,i,xi,bi);printf(n);,/非线性方程组求解主程序void SolveNonlinear(double x0N，double eps，int k)int nIter=0;double bN,JNN,deltaX N,Nor
22、mDeltaX;while(1)nIter+;GetFun(x0,b);GetJacobi(x0,J);NormDeltaX=GetDeltaX(J,b,deltaX);UpdateX(x0,deltaX);if(nIterk)|(NormDeltaXeps)break;,计算F(X0),计算F(X0),计算X0,更新计算结果,循环截止条件,/计算F(X)，结果放在b里void GetFun(x,y)double x,y;y0=x0*x0+x1*x1+x2*x2-1.0;y1=2.0*x0*x0+x1*x1-4.0*x2;y2=3.0*x0*x0-4.0*x1+x2*x2;return;/计算
23、F(X)，结果放在J里void GetJacobi(x,J)double x N,JNN;J00=2*x0;J01=2*x1;J02=2*x2;J10=4*x0;J11=2*x1;J12=-4;J20=6*x0;J21=-4;J22=2*x2;return;,/计算X0主程序double GetDeltaX(double JNN,double bN,double deltaXN)int nResultFlag,i;double NormDeltaX;nResultFlag=agaus(J,b,N);if(nResultFlag)NormDeltaX=0;for(i=0;iN;i+)NormDe
24、ltaX=NormDeltaX+fabs(bi);deltaXi=bi;return NormDeltaX;,解线性方程组求X0，存于b中,/高斯全主元法解方程，结果放在b里int agaus(a,b,n)int n;double aNN,bN;int*js,l,k,i,j,is;double d,t;js=malloc(n*sizeof(int);/记载列交换标号l=1;for(k=0;kd)d=t;jsk=j;/记住列is=i;/记住行,求全主元,if(d+1.0=1.0)l=0;else if(jsk!=k)/列交换for(i=0;i=n-1;i+)t=aik;aik=aijsk;aij
25、sk=t;if(is!=k)for(j=k;j=n-1;j+)/行交换 t=akj;akj=aisj;aisj=t;t=bk;bk=bis;bis=t;/b交换,列交换,行交换,判断主元是否为零,/最大主元为零不能计算出正确结果 if(l=0)free(js);printf(failn);return(0);/以下为消元过程d=akk;akk=1;for(j=k+1;j=n-1;j+)/当前行除以akkakj=akj/d;bk=bk/d;for(i=k+1;i=n-1;i+)/消元for(j=k+1;j=n-1;j+)aij=aij-aik*akj;bi=bi-aik*bk;,主元为零终止计算
26、,消元,/以下过程为回代过程d=an-1n-1;if(fabs(d)+1.0=1.0)free(js);printf(failn);return(0);bn-1=bn-1/d;for(i=n-2;i=0;i-)t=0.0;for(j=i+1;j=n-1;j+)t=t+aij*bj;bi=bi-t;,回代、求解,/以下过程为恢复x正常顺序过程jsn-1=n-1;for(k=n-1;k=0;k-)if(jsk!=k)t=bk;bk=bjsk;bjsk=t;free(js);return(1);,恢复解的顺序,void UpdateX(double x0N,double deltaXN)int i;for(i=0;iN;i+)x0i=x0i-deltaXi;,/更新计算结果,（相当于计算X1）,计算雅可比矩阵的方法,拟牛顿法,数值法解非线性方程组的主要步骤：,(1)选定初值X0、h0,(5)按公式计算X1,(2)计算F(X0),F(X0),(4)计算X0,(6)误差判断,h,(3)计算z0,
