《MATLAB线性方程组.ppt》由会员分享，可在线阅读，更多相关《MATLAB线性方程组.ppt（108页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第三章第三章 线性代数方程组及矩阵特征值线性代数方程组及矩阵特征值3.0 3.0 预备知识预备知识3.1 3.1 直接法直接法3.2 3.2 迭代法迭代法3.3 3.3 不可解问题不可解问题3.4 3.4 病态问题病态问题 3.0预备知识预备知识一、对角阵与三角阵一、对角阵与三角阵1、对角阵：、对角阵：udiag(A)提取提取mn的矩阵的矩阵A 的主对角线上元素，生的主对角线上元素，生 成一个具有成一个具有min(m,n)个元素的列向量个元素的列向量 diag(A,k)提取第提取第k条对角线的元素条对角线的元素udiag(V)设设V为具有为具有m个元素的向量，将产生一个以个元素的向量，将产生一
2、个以向量向量V的元素为主对角线元素的的元素为主对角线元素的mm对角矩阵对角矩阵diag(V,k)产生一个以向量产生一个以向量V的元素为第的元素为第k k条对角线的条对角线的元素的元素的nn(n=m+k)对角阵对角阵 2 2、矩阵的三角阵矩阵的三角阵u 下三角矩阵下三角矩阵 tril(A)提取矩阵提取矩阵A的下三角阵的下三角阵 tril(A,k)提取矩阵提取矩阵A的第的第k条对角线以下的元素条对角线以下的元素u上三角矩阵上三角矩阵 triu(A)提取矩阵提取矩阵A的上三角矩阵的上三角矩阵 triu(A,k)提取矩阵提取矩阵A的第的第k条对角线以下的元素条对角线以下的元素二、二、用于专门学科中的矩
3、阵用于专门学科中的矩阵(1)(1)魔方矩阵魔方矩阵 魔方矩阵是每行、每列及两条对角线上的魔方矩阵是每行、每列及两条对角线上的元素和都相等的矩阵。对于元素和都相等的矩阵。对于n阶魔方阵，其阶魔方阵，其元素由元素由1,2,3,n2共共n2 个整数组成个整数组成.magic(n)生成一个生成一个n阶魔方阵阶魔方阵,各行各列各行各列及两条及两条 对角线和为对角线和为(1+2+3+.+n2)/n 例例 magic(5)ans=17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9例：将例：将101125等等25个数填入一个个数填
4、入一个5行行5列的表格中，使列的表格中，使其每行每列及对角线的和均为其每行每列及对角线的和均为565。命令如下：。命令如下：M=100+magic(5)(2)(2)范得蒙矩阵范得蒙矩阵 范得蒙范得蒙(Vandermonde)矩阵最后一列全为矩阵最后一列全为1，倒，倒数第二列为一个指定的向量，其他各列是其后列与数第二列为一个指定的向量，其他各列是其后列与倒数第二列的点乘积。倒数第二列的点乘积。vander(V)生成以向量生成以向量V为基础向量的范得蒙矩阵。为基础向量的范得蒙矩阵。例例 A=vander(1;2;3;5)A=1 1 1 1 8 4 2 1 27 9 3 1 125 25 5 1 (
5、3)希尔伯特矩阵（希尔伯特矩阵（Hilbert）Hilbert矩阵的每个元素矩阵的每个元素hij=1/(i+j-1)hilb(n)生成生成n阶希尔伯特矩阵阶希尔伯特矩阵 invhilb(n)专求专求n阶的希尔伯特矩阵的逆矩阵阶的希尔伯特矩阵的逆矩阵 注意注意1：高阶：高阶Hilbert矩阵一般为病态矩阵，所以直矩阵一般为病态矩阵，所以直接求逆可能出现错误结论，故应该采用接求逆可能出现错误结论，故应该采用invhilb(n)注意注意2：由于：由于Hilbert矩阵本身接近奇异，所以建矩阵本身接近奇异，所以建议处理该矩阵时建议尽量采用符号运算工具箱，若议处理该矩阵时建议尽量采用符号运算工具箱，若采
6、用数值解时应该验证结果的正确性采用数值解时应该验证结果的正确性 (4)托普利兹矩阵托普利兹矩阵(toeplitz)toeplitz矩阵除第一行第一列外，其他每个元素都与矩阵除第一行第一列外，其他每个元素都与左上角的元素相同。左上角的元素相同。toeplitz(x,y)生成一个以生成一个以x为第一列，为第一列，y为第一行的为第一行的托普利兹矩阵。这里托普利兹矩阵。这里x,y均为向量，两者不必等长。均为向量，两者不必等长。toeplitz(x)用向量用向量x生成一个对称的托普利兹矩阵。生成一个对称的托普利兹矩阵。例例 T=toeplitz(1:5,1:7)T=1 2 3 4 5 6 7 2 1 2
7、 3 4 5 6 3 2 1 2 3 4 5 4 3 2 1 2 3 4 5 4 3 2 1 2 3 (5)帕斯卡矩阵帕斯卡矩阵 二次项二次项(x+y)n展开后的系数随展开后的系数随n的增大组成一的增大组成一个三角形表，称为杨辉三角形。由杨辉三角形表组个三角形表，称为杨辉三角形。由杨辉三角形表组成的矩阵称为帕斯卡成的矩阵称为帕斯卡(Pascal)矩阵。矩阵。pascal(n)生成一个生成一个n阶帕斯卡矩阵。阶帕斯卡矩阵。pascal(6)ans=1 1 1 1 1 1 1 2 3 4 5 6 1 3 6 10 15 21 1 4 10 20 35 5656 1 5 15 35 70 126 1
8、 6 21 56 126 252例例 求求(x+y)5的展开式。的展开式。pascal(6)次对角线上的元素次对角线上的元素1，5，10，10，5，1即为即为展开式的系数。展开式的系数。三、三、向量和矩阵的范数向量和矩阵的范数u norm(V)或或norm(V,2)求向量求向量V（或矩阵（或矩阵V)的的2范数范数u norm(V,1)求向量求向量V（或矩阵（或矩阵V）的）的1范数范数u norm(V,inf)求向量求向量V（或矩阵（或矩阵V)的的范数范数 例例 已知已知V，求，求V的的3种范数。种范数。命令如下：命令如下：V=-1,1/2,1;v1=norm(V,1)%求求V的的1范数范数 v
9、2=norm(V)%求求V的的2范数范数 vinf=norm(V,inf)%求求V的的范数范数常用的向量范数常用的向量范数:范数意义下的单位向量范数意义下的单位向量:X=x1,x2T1-11|X|1=111-1-1|X|2=1-111-1-1|X|=1常用的矩阵范数常用的矩阵范数:例例 求矩阵求矩阵A的三种范数的三种范数命令如下：命令如下：A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,21,3;11,18,25,2,19;a1=norm(A,1)%求求A的的1范数范数 a2=norm(A)%求求A的的2范数范数 ainf=norm(A,inf)
10、%求求A的的范数范数 四、四、矩阵的逆与伪逆矩阵的逆与伪逆1、矩阵的逆（后面研究完矩阵的逆（后面研究完Gauss消去法后还将给出求消去法后还将给出求逆的方法）逆的方法）求方阵求方阵A的逆可用的逆可用 inv(A)注意：该函数也适用于符号变量构成的矩阵的求逆注意：该函数也适用于符号变量构成的矩阵的求逆例例 用求逆矩阵的方法解线性方程组。用求逆矩阵的方法解线性方程组。命令如下：命令如下：A=1,2,3;1,4,9;1,8,27;b=5,2,6;x=inv(A)*b 一般情况下，用左除一般情况下，用左除x=Ab比求矩阵的逆的方法更比求矩阵的逆的方法更有效有效(因为因为A奇异或接近奇异时，用奇异或接近
11、奇异时，用inv(A)可能出错可能出错)例：例：Hilbert矩阵（非常有名的病态矩阵）：矩阵（非常有名的病态矩阵）：验证从验证从55到到1414的的Hilbert矩阵病态特征矩阵病态特征clearformat long;for n=5:14 H=hilb(n);norm1=norm(H*inv(H)-eye(size(H);H1=invhilb(n);norm2=norm(H*H1-eye(size(H);fprintf(n=%3.0f norm1=%e norm2=%en,n,norm1,norm2)end n=5 norm1=1.409442e-011 norm2=3.637979e-0
12、12 n=6 norm1=2.534893e-009 norm2=1.455203e-011 n=7 norm1=9.810538e-008 norm2=5.208793e-010 n=8 norm1=3.123671e-006 norm2=4.822187e-008 n=9 norm1=8.116595e-005 norm2=2.479206e-007 n=10 norm1=2.645008e-003 norm2=1.612897e-005 n=11 norm1=7.200720e-002 norm2=1.305122e-003Warning:Matrix is close to singu
13、lar or badly scaled.Results may be inaccurate.RCOND=2.632766e-017.n=12 norm1=1.176913e+001 norm2=5.576679e-002Warning:Matrix is close to singular or badly scaled.Results may be inaccurate.RCOND=2.339949e-018.n=13 norm1=5.323696e+001 norm2=1.137063e+001Warning:Matrix is close to singular or badly sca
14、led.Results may be inaccurate.RCOND=1.708191e-019.n=14 norm1=1.224232e+004 norm2=2.836298e+002说明说明1:对于对于Hilbert求逆时，不建议用求逆时，不建议用inv,可用可用 invhilb直接产生逆矩阵直接产生逆矩阵说明说明2：符号工具箱中也对符号矩阵定义了符号工具箱中也对符号矩阵定义了inv()函数函数,即使对更高阶的奇异矩阵也可以精确的求解出逆矩阵即使对更高阶的奇异矩阵也可以精确的求解出逆矩阵例：例：H=sym(hilb(30);norm(double(H*inv(H)-eye(size(H)
15、结果为结果为ans=0说明说明3：对于奇异对于奇异矩阵用数值方法的矩阵用数值方法的inv()函数，会给出函数，会给出错误的结果，但符号工具箱中的错误的结果，但符号工具箱中的inv()会给出正确的结论会给出正确的结论例例 A=16,2,3,13;5,11,10,8;9,7,6,12;4,14,15,1;det(A)B=inv(A)结果：结果：ans=0Warning:Matrix is close to singular or badly scaled.Results may be inaccurate.RCOND=1.306145e-017.B=1.0e+014*但用符号工具箱的但用符号工具箱
16、的inv:A=16,2,3,13;5,11,10,8;9,7,6,12;4,14,15,1;A=sym(A)inv(A)结果：结果：?Error using=sym.invError,(in inverse)singular matrix 2、矩阵的伪逆矩阵的伪逆 pinv(A)若若A不是一个方阵，或不是一个方阵，或A是一个非满秩的方阵时，矩阵是一个非满秩的方阵时，矩阵A没有逆矩阵，但可以找到一个与没有逆矩阵，但可以找到一个与A的转置矩阵的转置矩阵A同同型的矩阵型的矩阵B，使得：，使得：ABA=A,BAB=B此时称矩阵此时称矩阵B为矩阵为矩阵A的伪逆。的伪逆。例例 求求A的伪逆，并将结果送的伪
17、逆，并将结果送BA=3,1,1,1;1,3,1,1;1,1,3,1;B=pinv(A)例例 求矩阵求矩阵A的伪逆的伪逆 A=0,0,0;0,1,0;0,0,1;pinv(A)五、五、求方阵求方阵A的行列式的行列式:det(A)例例 用克莱姆用克莱姆(Cramer)方法求解线性方程组方法求解线性方程组(不建议使用）不建议使用）程序如下：程序如下：D=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2;%定义系数矩阵定义系数矩阵b=4;6;12;6;%定义常数项向量定义常数项向量D1=b,D(:,2:4);%用方程组的右端向量置换用方程组的右端向量置换D的第的第1列列 D2=D
18、(:,1),b,D(:,3:4);%用方程组的右端向量置换用方程组的右端向量置换D的第的第2列列D3=D(:,1:2),b,D(:,4:4);%用方程组的右端向量置换用方程组的右端向量置换D的的第第3列列D4=D(:,1:3),b;%用方程组的右端向量置换用方程组的右端向量置换D的第的第4列列DD=det(D);x1=det(D1)/DD;x2=det(D2)/DD;x3=det(D3)/DD;x4=det(D4)/DD;x1,x2,x3,x4 六、六、求矩阵求矩阵A的秩的秩 rank(A)七、七、求矩阵求矩阵A的迹的迹 trace(A)例：例：D=2,2,-1,1;4,3,-1,2;8,5,
19、-3,4;3,3,-2,2;r=rank(D)t=trace(A)结果：结果：r=4 t=6九、求矩阵特征多项式、特征值、特征向量九、求矩阵特征多项式、特征值、特征向量设设A是一个是一个 nn 矩阵，矩阵，为为A的特征多项式。的特征多项式。特征多项式的根称为矩阵特征多项式的根称为矩阵A的特征值。的特征值。c=poly(A)求矩阵求矩阵A的特征多项式的系数的特征多项式的系数roots(c)求多项式求多项式c的根的根八八、求矩阵求矩阵A的共轭矩阵的共轭矩阵 conj(A)eig(A)求矩阵求矩阵A的特征值的特征值 常用的调用格式有常用的调用格式有：uE=eig(A)求矩阵求矩阵A的全部特征值，构成
20、向量的全部特征值，构成向量E。uV,D=eig(A)求矩阵求矩阵A的全部特征值，构成对角阵的全部特征值，构成对角阵D，并求，并求A的特征向量构成的特征向量构成V的列向量。的列向量。uV,D=eig(A,nobalance)与第与第2种格式类似，但种格式类似，但第第2种格式中先对种格式中先对A作相似变换后求矩阵作相似变换后求矩阵A的特征值和的特征值和特征向量，而格式特征向量，而格式3直接求矩阵直接求矩阵A的特征值和特征向量，的特征值和特征向量，即即A中某项非常小，这样求出的特征值及特征向量更中某项非常小，这样求出的特征值及特征向量更精确。精确。uV,D=eig(A,B)计算广义特征值和特征向量，
21、使计算广义特征值和特征向量，使 AV=BVD。例：设矩阵例：设矩阵A=3 4-2;3-1 1;2 0 5;E=eig(A)V,D=eig(A)V,D=eig(A,nobalance)现求解线性方程组现求解线性方程组 情形1：m=n（恰定方程组）（恰定方程组）在MATLAB中的求解命令有：情形情形2 2：mn（超定方程），多用于曲线拟合（超定方程），多用于曲线拟合。解线性方程组的一般函数文件如下：解线性方程组的一般函数文件如下：function x,y=line_solution(A,b)m,n=size(A);y=;if norm(b)0%b非零，非零，非齐次方程组非齐次方程组 if rank
22、(A)=rank(A,b)%方程组相容方程组相容（有解有解）if rank(A)=n%有有唯一解唯一解 x=Ab;else%方程组有方程组有无穷多个解无穷多个解,基础解系基础解系 disp(原方程组有有无穷个解，其齐次方程组的基础原方程组有有无穷个解，其齐次方程组的基础 解系为解系为y，特解为，特解为x);y=null(A,r);x=Ab;end else%方程组不相容方程组不相容（无解）（无解），给出最小二乘解，给出最小二乘解 disp(方程组的最小二乘法解是：方程组的最小二乘法解是：);x=Ab;endelse%齐次方程组齐次方程组 if rank(A)=n%列满秩列满秩 x=zero(n
23、,1)%0解解 else%非非0解，无穷多个解解，无穷多个解 disp(方程组有无穷个解，基础解系为方程组有无穷个解，基础解系为x);x=null(A,r);end endreturn如在如在MATLAB命令窗口，输入命令命令窗口，输入命令 A=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2;b=4,6,12,6;x,y=line_solution(A,b)及：及：A=2,7,3,1;3,5,2,2;9,4,1,7;b=6,4,2;x,y=line_solution(A,b)分别显示其求解结果。分别显示其求解结果。求解线性方程组（m=n)用克莱姆法则，理论上可行，但实际
24、运算时工作量大，耗时。下面研究、一、一、GaussGauss简单消元法（简单消元法（m=n)m=n)设设 有有用线性代数中的克莱姆法则求解时，工作量非常大，用线性代数中的克莱姆法则求解时，工作量非常大，工作量小的方法是工作量小的方法是 Gauss Gauss消元法。消元法。3.1 3.1 解线性方程组的直接法解线性方程组的直接法消消 元：元：以此类推，最后方程组化为：回回 代：代：失效，故应选主元失效，故应选主元 二、列主元素消去法二、列主元素消去法-计算结果可靠计算结果可靠到此原方程组化为到此原方程组化为到此原方程组化为到此原方程组化为(上三角方程组上三角方程组)（3.2)3.2)(n-1)
25、原方程组化为原方程组化为以上为消元过程。以上为消元过程。(n)回代求解公式回代求解公式 （3.3)3.3)是回代过程。是回代过程。（3.3)3.3)例例1：在：在MATLAB上，用上，用Gauss列主元消去法求解方程组：列主元消去法求解方程组：cleara=-0.04 0.04 0.12 3;0.56-1.56 0.32 1;-0.24 1.24-0.28 0%先定义增广矩阵先定义增广矩阵x=0,0,0;%先将解设为零向量，后面重新赋值先将解设为零向量，后面重新赋值tempo=a(2,:);a(2,:)=a(1,:);a(1,:)=tempo;a%第一第一次选主元（第一行与第二行交换）次选主元
26、（第一行与第二行交换）a(2,:)=a(2,:)-a(1,:)*a(2,1)/a(1,1);%将第一个对角元将第一个对角元下面的数字消为下面的数字消为0a(3,:)=a(3,:)-a(1,:)*a(3,1)/a(1,1);aa=0.5600 -1.5600 0.3200 1.0000 0 -0.0714 0.1429 3.0714 0.0000 0.5714 -0.1429 0.4286tempo=a(3,:);a(3,:)=a(2,:);a(2,:)=tempo;a%第二第二次选主元次选主元a(3,:)=a(3,:)-a(2,:)*a(3,2)/a(2,2);a%第二次将第第二次将第二个对角
27、元下的数字消为二个对角元下的数字消为0 x(3)=a(3,4)/a(3,3);%消去完成，现在开始回代消去完成，现在开始回代x(2)=(a(2,4)-a(2,3)*x(3)/a(2,2);x(1)=(a(1,4)-a(1,2:3)*x(2:3)/a(1,1);x运行得方程组的解为：运行得方程组的解为：a=0.5600 -1.5600 0.3200 1.0000 0.0000 0.5714 -0.1429 0.4286 0.0000 0 0.1250 3.1250 x=7.0000 7.0000 25.0000也可直接用也可直接用 x=Ab说明：说明：(1)也也可可采采用用无无回回代代的的列列主
28、主元元消消去去法法(叫叫Gauss-Jordan消消去去法法)，该该法法同同时时消消去去对对角角元元上上下下的的元元素素，且且仍仍旧旧需需要要选选主主元元，但但比比有有回回代代的的列列主主元元消消去去法法的的乘乘除除运运算算次次数数多多。GaussJordan消消去去法法的的优优点点之一是用它来算逆矩阵的算法非常容易解释。之一是用它来算逆矩阵的算法非常容易解释。(2)有有回回代代的的列列主主元元消消去去法法所所进进行行的的乘乘除除运运算算次次数数为为 ，计算量很小。，计算量很小。例例2：用：用GaussJordan消去法求解上例中的矩阵消去法求解上例中的矩阵 的逆矩阵。的逆矩阵。clearA=
29、-0.04 0.04 0.12;0.56-1.56 0.32;-0.24 1.24-0.28a=A,eye(3);atempo=a(2,:);a(2,:)=a(1,:);a(1,:)=tempo;a%第一第一次选主元，并交换第一行和第二行次选主元，并交换第一行和第二行a(1,:)=a(1,:)/a(1,1)%将主元标准化将主元标准化for i=2:3;a(i,:)=a(i,:)-a(i,1)*a(1,:);end;%将主元下的元素消为零将主元下的元素消为零atempo=a(3,:);a(3,:)=a(2,:);a(2,:)=tempo;a%第二第二次选主元，交换第二行和第三行次选主元，交换第二
30、行和第三行a(2,:)=a(2,:)/a(2,2);a%第二个主元标准化第二个主元标准化a=0.5600 -1.5600 0.3200 1.0000 0 -0.0714 0.1429 3.0714 0.0000 0.5714 -0.1429 0.4286for i=1:3 if i=2,a(i,:)=a(i,:)-a(i,2)*a(2,:);endend%a(2,2)的上下元素消为的上下元素消为0aa=1.0000 0 -0.1250 0 3.8750 4.8750 0 1.0000 -0.2500 0 0.7500 1.7500 0 0 0.1250 1.0000 0.1250 0.1250
31、a(3,:)=a(3,:)/a(3,3)for i=1:3;if i=3,a(i,:)=a(i,:)-a(i,3)*a(3,:);end;end;aA_inv=a(:,4:6)A*A_inv结果为：结果为：A_inv=1.0000 4.0000 5.0000 2.0000 1.0000 2.0000 8.0000 1.0000 1.0000也可用也可用rref命令来求命令来求三、三、Gauss Gauss 全主元消去法：全主元消去法：优点优点-计算结果更可靠；计算结果更可靠；缺点缺点-挑主元花时间更多，挑主元花时间更多，次序有变动，程序复杂。次序有变动，程序复杂。四、四、Gauss消元法的应用
32、消元法的应用 1、求行列式：、求行列式：det(A)2、求逆矩阵：、求逆矩阵：inv（A）或）或 A(-1)rref(A,eye(size(A)将将(A,E)化化为为行行最最简简形形，其其实就是实就是Gauss-Jordon消去法消去法（3）求解方程组）求解方程组Ax=b u求逆矩阵的思想求逆矩阵的思想:x=inv(A)*b或或x=A(-1)*bu左左除除法法（原原理理上上是是运运用用高高斯斯消消元元法法求求解解,但但MatlabMatlab在在实际执行过程中是通过实际执行过程中是通过LULU分解法进行的）分解法进行的）：x=Abu符符号号矩矩阵阵法法（此此法法最最接接近近精精确确值值，但但运
33、运算算速速度度慢慢）x=sym(A)sym(b)u化为行最简形：化为行最简形：C=A,b rref(C)例例1：在：在MATLAB上，用上，用Gauss列主元消去法求解方程组：列主元消去法求解方程组：A=-0.04 0.04 0.12;0.56-1.56 0.32;-0.24 1.24-0.28;b=3;1;0;x11=inv(A)*b%法法1：求逆思想：求逆思想x12=A(-1)*b%法法1：求逆思想：求逆思想x3=Ab%法法2：左除法：左除法x4=sym(A)sym(b)%法法3：符号矩阵法：符号矩阵法C=A,b;rref(C)%法法4：化为行最简形：化为行最简形定义定义3.13.1 叫叫
34、的三角（因子）分解，其中的三角（因子）分解，其中 是是是上三角。是上三角。下三角下三角,若若L为单位下三角阵（对角元全为为单位下三角阵（对角元全为1 1），），U为上为上三角阵，则称三角阵，则称 A=LU 为为Doolittle分解分解;若若 L 是是下三角，下三角，U是单位上三角，则称是单位上三角，则称A=LU为为Crout分解分解。定理定理3.1 3.1 n 阶阵阶阵 有唯一有唯一Doolittle分解分解(Crout(Crout）的所有的所有顺序主子式不为顺序主子式不为0.三角分解不唯一三角分解不唯一,为此引入为此引入定义定义3.2 3.2 五利用矩阵三角分解法求解方程组五利用矩阵三角分
35、解法求解方程组 为什么要讨论三角分解？若在消元法进行前能实为什么要讨论三角分解？若在消元法进行前能实现三角分解现三角分解 ，则则 容易回代求解容易回代求解 在在Gauss消去法中，选主元改变了行的次序，消去法中，选主元改变了行的次序，尽管对于尽管对于Gauss消去法来说，这种次序的变换无法消去法来说，这种次序的变换无法事先知道，但是这个变化的影响却可以用一个算事先知道，但是这个变化的影响却可以用一个算子子P表示，其中表示，其中P是一个置换矩阵。用是一个置换矩阵。用P左乘原始左乘原始矩阵矩阵A得到：得到：PAx=Py 或或 对对 做做Gauss消去法不需要选主元，所以对消去法不需要选主元，所以对
36、 做做LU分解，同样不需要选主元。分解，同样不需要选主元。实际上，将选主元实际上，将选主元Gauss消去法里的行交换同消去法里的行交换同样作用于单位矩阵，所得矩阵即为样作用于单位矩阵，所得矩阵即为P。在在MATLAB中，中，LU分解的命令是分解的命令是lu，有两种格式：，有两种格式：u l,u,p=lu(A)其中其中A是待分解矩阵；是待分解矩阵；l,u,p分别代表分别代表L(主对角线元素主对角线元素为为1的下三角矩阵）、的下三角矩阵）、U（上三角矩阵）和由单位阵变上三角矩阵）和由单位阵变换出的置换矩阵换出的置换矩阵P 满足：满足：PA=LU，即，即 A=P-1LUu l,u=lu(A)其中其中
37、l=P-1L，所以，所以A=lU=P-1LU 用此格式求出来的用此格式求出来的 l 并不一定是真正的下三角矩阵，并不一定是真正的下三角矩阵，需要换行后才能是真正的下三角矩阵需要换行后才能是真正的下三角矩阵u实现实现LU分解后，分解后，若采用若采用l,u=lu(A)，则，则Ax=b的解为的解为 x=u(l b)若采用若采用l，u，p=lu(A)，则，则Ax=b的解为的解为 x=u(l(p*b)例例 利用利用LULU分解法求解方程组分解法求解方程组A=2 4 2 6;4 9 6 15;2 6 9 18;6 15 18 40%先录入系数矩阵先录入系数矩阵b=9 23 22 47L,U=lu(A)%对
38、对A矩阵做矩阵做LU分解分解y=Lb%求解方程组求解方程组Ly=bx=Uy%求解方程组求解方程组Ux=y得到方程组的最终解得到方程组的最终解故方程组的最终解为故方程组的最终解为:x=(0.5000,2.0000,3.0000,-1.0000)或或A=2 4 2 6;4 9 6 15;2 6 9 18;6 15 18 40%先录入系数矩阵先录入系数矩阵b=9 23 22 47L,U,P=lu(A)%对对A矩阵做矩阵做LU分解分解y=LP*b%求解方程组求解方程组Ly=Pbx=Uy%求解方程组求解方程组Ux=y得到方程组的最终解得到方程组的最终解lX=QR：X为方阵，为方阵，Q为正交矩阵，为正交矩
39、阵，R为上三角矩阵为上三角矩阵 六、利用矩阵六、利用矩阵QR分解法求解方程组分解法求解方程组lQ，R=qr(X)：产生一个一个正交矩阵：产生一个一个正交矩阵Q和一个上和一个上三角矩阵三角矩阵R，满足，满足X=QRlQ，R，E=qr(X)：产生一个一个正交矩阵：产生一个一个正交矩阵Q、一个、一个上三角矩阵上三角矩阵R、置换矩阵、置换矩阵E，满足，满足XE=QRl实现实现QR分解后，分解后，若采用若采用Q，R=qr(X)，则，则Ax=b的解为的解为 x=R(Qb)若采用若采用Q，R，E=qr(X)，则，则Ax=b的解为的解为 x=E(R(Qb)例：例：用用QR分解求解线性方程组。分解求解线性方程组
40、。A=3，1，-4，1;1，-3，0，2;0，2，1，-1;1，6，-1，-3;b=12，-6，4，0;Q，R=qr(A);x=R(Qb)或采用或采用QR分解的第分解的第2种格式，命令如下：种格式，命令如下：Q，R，E=qr(A);x=E*(R(Qb)两种得到结果都是两种得到结果都是 x=-16.4444 20.6667 -1.1111 36.2222七、利用矩阵七、利用矩阵Cholesky分解法求解方程组分解法求解方程组l若矩阵若矩阵X是对称正定的，是对称正定的，X=RR,R为上三角矩阵为上三角矩阵lR=chol(X)：产生一个上三角阵：产生一个上三角阵R，使，使RR=X 若若X为非对称正定
41、，则输出一个出错信息为非对称正定，则输出一个出错信息lR，p=chol(X)：这个命令格式将不输出出错信息。：这个命令格式将不输出出错信息。当当X为对称正定的，则为对称正定的，则p=0，R与上述格式得到的结与上述格式得到的结果相同；否则果相同；否则p为一个正整数。为一个正整数。如果如果X为满秩矩阵，则为满秩矩阵，则R为一个阶数为为一个阶数为q=p-1的上三的上三角阵，且满足角阵，且满足 RR=X(1:q，1:q)l实现实现Cholesky分解后，线性方程组分解后，线性方程组Ax=b变成变成RRx=b，所以，所以x=R(Rb)。【例】用【例】用Cholesky分解求解线性方程组。分解求解线性方程
42、组。A=3，1，-4，1;1，-3，0，2;0，2，1，-1;1，6，-1，-3;b=12，-6，4，0;R=chol(A)l结果：结果：?Error using=chol Matrix must be positive definitel命令执行时，出现错误信息，说明命令执行时，出现错误信息，说明A为非正为非正定矩阵。定矩阵。八、解三对角方程组八、解三对角方程组追赶法追赶法给定方程组给定方程组按行严格对角占优按行严格对角占优（三对角方程组）（三对角方程组）其求解算法是其求解算法是Gauss消去法的一种变形，称为三对角消去法的一种变形，称为三对角法（追赶法）。解法如下：法（追赶法）。解法如下：
43、1、由第一个方程、由第一个方程令令2、将其代入第二个方程、将其代入第二个方程 ，得：，得：再令再令3、将其代入第三个方程、将其代入第三个方程得得再令再令以此类推，以此类推，令令将将代入最后一个方程代入最后一个方程令令以上过程称为以上过程称为“追追”；总结总结“追追”的算法：的算法：Step1：（初始化变量）：（初始化变量）Step2：下面开始下面开始“赶赶”：Step3：先求：先求Step4：再求其他的：再求其他的三对角方程组的求解程序三对角方程组的求解程序 tri_diag.m 如下：如下：%tri_diag(a,b,c,d,n)solves a tridiagonal equation.f
44、unction x=tri_diag(a,b,c,d,n)for i=2:n r=a(i)/b(i-1);d(i)=d(i)-r*d(i-1);b(i)=b(i)-r*c(i-1);endd(n)=d(n)/b(n);for i=n-1:-1:1 d(i)=(d(i)-c(i)*d(i+1)/b(i);endx=d;迭代法适用于解大型稀疏方程组迭代法适用于解大型稀疏方程组(万阶以上的方万阶以上的方程组程组,系数矩阵中零元素占很大比例系数矩阵中零元素占很大比例,而非零元按某种而非零元按某种模式分布模式分布)背景背景:电路分析、边值问题的数值解和数学物理方程电路分析、边值问题的数值解和数学物理方程
45、问题问题:(1)如何构造迭代格式？如何构造迭代格式？(2)迭代格式是否收敛？迭代格式是否收敛？(3)收敛速度如何？收敛速度如何？(4)如何进行误差估计？如何进行误差估计？3.2 3.2 解线性方程组的迭代法解线性方程组的迭代法一、稀疏矩阵存储方式的产生与转化一、稀疏矩阵存储方式的产生与转化 1、由由完全存储完全存储方式方式 转为转为 稀疏存储稀疏存储方式命令：方式命令：B=sparse(A)将矩阵将矩阵A A转化为稀疏存储方式的矩阵转化为稀疏存储方式的矩阵B B sparse(m,n)生成一个生成一个mnmn的所有元素都是的所有元素都是0 0的稀疏的稀疏 矩阵矩阵sparse(u,v,S)u,
46、v,S是是3 3个等长的向量。个等长的向量。S是要建立的稀疏矩阵的非是要建立的稀疏矩阵的非0 0元素；元素；u(i)、v(i)分别是分别是S(i)的行标和列标；的行标和列标；该函数生成一个该函数生成一个max(u)行、行、max(v)列并以列并以S为稀疏元素为稀疏元素的稀疏矩阵的稀疏矩阵 u=1,1,4;v=1,2,1;s=1,3,7;sparse(u,v,s)结果结果：ans=(1,1)1 (4,1)7 (1,2)3U,V,S=find(A)返回矩阵返回矩阵A中非中非0 0元素的下标和元素，元素的下标和元素，这里产生的这里产生的U、V、S可作可作sparse(u,v,s)的参数的参数full
47、(A)返回稀疏存储矩阵返回稀疏存储矩阵A A对应的对应的完全存储方式完全存储方式矩阵矩阵相关操作的函数：相关操作的函数：2 2、产生一个稀疏矩阵产生一个稀疏矩阵 把要建立的稀疏矩阵的非把要建立的稀疏矩阵的非0元素及其所在行和列的元素及其所在行和列的位置表示出来后由位置表示出来后由MATLAB自己产生其稀疏存储方式，自己产生其稀疏存储方式，这需要使用这需要使用spconvert函数。调用格式为：函数。调用格式为：B=spconvert(A)其中其中A为一个为一个m3或或m4的矩阵，每行表示一个非的矩阵，每行表示一个非0元元素，素，m是非是非0元素的个数，元素的个数，A每个元素的意义是：每个元素的
48、意义是：(i,1)第第i个非个非0元素所在的行。元素所在的行。(i,2)第第i个非个非0元素所在的列。元素所在的列。(i,3)第第i个非个非0元素值的实部。元素值的实部。(i,4)第第i个非个非0元素值的虚部，若矩阵的全部元素都是元素值的虚部，若矩阵的全部元素都是实数，则无须第四列。实数，则无须第四列。A=1,1,4;1,2,1;1,3,7;spconvert(A)结果：结果：ans=(1,1)4 (1,2)1 (1,3)7 3 3、单位稀疏矩阵的产生、单位稀疏矩阵的产生 eye 产生一个完全存储方式的单位矩阵产生一个完全存储方式的单位矩阵 speye(m,n)产生产生一个一个mn的的稀疏存储
49、单位稀疏存储单位矩阵矩阵 speye(2,3)ans=(1,1)1 (2,2)1二、简单迭代法二、简单迭代法 1.1.迭代法建立迭代法建立.考虑考虑 Ax=b (矩阵矩阵B不唯一不唯一)对应写出对应写出 产生向量序列产生向量序列若收敛若收敛,记记则于则于(1)(1)两端取极限有：两端取极限有：上式说明：上式说明：是解向量是解向量 ，从而当，从而当k充分大时充分大时注意注意:迭代阵迭代阵B不唯一，不唯一，B的选取的选取影响收敛性。影响收敛性。解向量解向量(1)(1)叫简单迭代法叫简单迭代法,B叫迭代矩阵。叫迭代矩阵。2.2.收敛性收敛性.定义定义 称称 为矩阵为矩阵B B的谱半径。的谱半径。定理
50、 定理定理 对于简单迭代法，若迭代矩阵对于简单迭代法，若迭代矩阵设有方程组设有方程组(其中其中 )Ax=b，即，即 (2)(2)作等价变形作等价变形(3)(3)-Jacobi-Jacobi迭代法迭代法(k=0,1,2,)(4)(4)法法1 1：法法2：Ax=bA=D+(L+U)Dx=-(L+U)x+b=x=-D-1(L+U)x+D-1 bLU=x=BJ x+f=迭代公式：迭代公式：x(k+1)=BJ x(k)+f,(k=0,1,2)BJ=-D-1（L+U）,f=D-1bD编写实现编写实现Jacobo迭代法的函数迭代法的函数jacobi.m如下如下:function s=jacobi(a,b,x