第三章解线性方程组的直接方法的matlab程序.doc

1、高等教育出版社 教育电子音像出版社 作者：任玉杰 第三章 解线性方程组的直接方法的MATLAB程序第三章 解线性方程组的直接方法在科技、工程、医学和经济等各个领域中，经常遇到求解包含个未知数、由个方程构成的线性方程组 （3.1）的问题.线性方程组（3.1）的解法一般有直接法和迭代法.迭代法是求解大型稀疏矩阵方程组，尤其是由微分方程离散后得到的大型方程组的重要方法.在这章中我们要学习线性方程组的直接法，特别是适合用数学软件在计算机上求解的方法.3.1 方程组的逆矩阵解法及其MATLAB程序3.1.1 逆矩阵、行列式及其MATLAB命令当矩阵为方阵时，的行列式通常表示为或.当时，可逆，且可以用表

2、3-1中的MATLAB命令求的逆矩阵和行列式. 表 3-1命 令功 能NA=inv(A)输入矩阵，运行后输出的逆矩阵；HA=det(A)输入矩阵，运行后输出的行列式的值.3.1.2 方程组的逆矩阵解法及其MATLAB程序在MATLAB中引进了矩阵除法的概念，它包括左除和右除.设两矩阵为阶方阵，且，皆可逆，则矩阵方程，,的解都可以用表 3-2中的命令求出. 表 3-2命 令相同功能的命令功 能X=ACX= inv(A)*C输入矩阵和C，运行后输出矩阵方程的解Y=D/BY=D*inv(B)输入矩阵和，运行后输出矩阵方程的解Z=A F/BZ= inv(A)*F* inv(B)输入矩阵，和，运行后输出

3、矩阵方程的解3.1.3 线性方程组有解的判定条件及其MATLAB程序解非齐次线性方程组，将其增广矩阵化成行阶梯形矩阵，便可判断其是否有解若有解，化成行最简形矩阵，便可写出其通解.判定线性方程组是否有解的MATLAB程序输入的量：系数矩阵和常数向量；输出的量：RA和RB分别是系数矩阵和增广矩阵的秩, n是方程组中未知量的个数，另外还输出有关方程组解的信息.说明：(1)求矩阵的秩用MATLAB命令：RA=rank(A).(2)如果有唯一解，则用表 3-2方法求解; 如果有无穷多解，则将其增广矩阵化成行最简形矩阵，便可写出其通解.根据线性方程组有解的判定条件，现提供判定线性方程组是否有解的MATLA

4、B程序function RA,RB,n=jiepb(A,b)B=A b;n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA;if zhica0,disp(请注意：因为RA=RB，所以此方程组无解.)returnendif RA=RB if RA=ndisp(请注意：因为RA=RB=n，所以此方程组有唯一解.) else disp(请注意：因为RA=RB A=2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7; b= 0; 0; 0; 0; RA,RB,n=jiepb(A,b)运行后输出结果为请注意：因为RA=RB=n，所以此方程

5、组有唯一解.RA = 4,RB =4,n =4在MATLAB工作窗口输入X=Ab, 运行后输出结果为 X =(0 0 0 0).(2) 在MATLAB工作窗口输入程序 A=3 4 -5 7;2 -3 3 -2;4 11 -13 16;7 -2 1 3;b= 0; 0; 0; 0;RA,RB,n=jiepb(A,b)运行后输出结果请注意：因为RA=RB A=4 2 -1;3 -1 2;11 3 0; b=2;10;8; RA,RB,n=jiepb(A,B)运行后输出结果请注意：因为RA=RB，所以此方程组无解.RA =2,RB =3,n =3(4)在MATLAB工作窗口输入程序 A=2 1 -1

6、 1;4 2 -2 1;2 1 -1 -1; b=1; 2; 1; RA,RB,n=jiepb(A,b)运行后输出结果请注意：因为RA=RB0,disp(请注意：因为RA=RB，所以此方程组无解.)returnendif RA=RB if RA=ndisp(请注意：因为RA=RB=n，所以此方程组有唯一解.) X=zeros(n,1); X(n)=b(n)/A(n,n);for k=n-1:-1:1 X(k)=(b(k)-sum(A(k,k+1:n)*X(k+1:n)/A(k,k); end else disp(请注意：因为RA=RBA=5 -1 2 3;0 -2 7 -4;0 0 6 5;0

7、 0 0 3;b=20; -7; 4;6; RA,RB,n,X=shangsan(A,b)运行后输出结果请注意：因为RA=RB=n，所以此方程组有唯一解.RA = RB =4, 4,n = 4,X =2.4 -4.0 -1.0 2.03.3 高斯（Gauss）消元法和列主元消元法及其MATLAB程序3.3.1 高斯消元法及其MATLAB程序现提供的MATLAB程序是利用高斯消元法解线性方程组的MATLAB程序.用高斯消元法解线性方程组的MATLAB程序输入的量：系数矩阵和常系数向量；输出的量：系数矩阵和增广矩阵的秩RA,RB, 方程组中未知量的个数n和有关方程组解及其解的信息.function

8、 RA,RB,n,X=gaus(A,b)B=A b; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA;if zhica0,disp(请注意：因为RA=RB，所以此方程组无解.)returnendif RA=RB if RA=ndisp(请注意：因为RA=RB=n，所以此方程组有唯一解.) X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1for k=p+1:n m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endend b=B(1:n,n+1);A

9、=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); endelse disp(请注意：因为RA=RB A=1 -1 1 -3; 0 -1 -1 1;2 -2 -4 6;1 -2 -4 1; b=1;0; -1;-1; RA,RB,n,X =gaus (A,b)运行后输出结果请注意：因为RA=RB=n，所以此方程组有唯一解.X = 0 -0.5000 0.5000 0RA = 4RB = 4n = 43.3.2 列主元消元法及其MATLAB程序现提供的MATLAB程序是利

10、用回代法解上三角形线性方程组.关于解下三角形线性方程组的程序与之类似.用列主元消元法解线性方程组的MATLAB程序输入的量：系数矩阵和常系数向量；输出的量：系数矩阵和增广矩阵的秩RA,RB, 方程组中未知量的个数n和有关方程组解及其解的信息.function RA,RB,n,X=liezhu(A,b)B=A b; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA;if zhica0,disp(请注意：因为RA=RB，所以此方程组无解.)returnendif RA=RB if RA=ndisp(请注意：因为RA=RB=n，所以此方程组有唯一解.)

11、X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1Y,j=max(abs(B(p:n,p); C=B(p,:);B(p,:)= B(j+p-1,:); B(j+p-1,:)=C;for k=p+1:n m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endend b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); endelse disp(请注意

12、：因为RA=RB A=0 -1 -1 1;1 -1 1 -3;2 -2 -4 6;1 -2 -4 1; b=0;1;-1;-1; RA,RB,n,X=liezhu(A,b)运行后输出结果请注意：因为RA=RB=n，所以此方程组有唯一解.RA = 4，RB = 4，n = 4，X =0 -0.5 0.5 03.4 LU分解法及其MATLAB程序3.4.1判断矩阵LU分解的充要条件及其MATLAB程序根据定理3.5和（3.16）式，现提供MATLAB程序如下：判断矩阵能否进行LU分解的MATLAB程序输入的量：系数矩阵；运行后输出的量：矩阵的秩R和各阶顺序主子式） 的值hl及其相关信息，能否进行L

13、U分解的信息.function hl=pdLUfj(A)n n =size(A); RA=rank(A); if RA=ndisp(请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下：), RA,hl=det(A); returnendif RA=n for p=1:n,h(p)=det(A(1:p, 1:p);, endhl=h(1:n);for i=1:nif h(1,i)=0disp(请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：),hl;RA,returnendend if h(1,i)=0 disp(请注

14、意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：)hl;RAendend例3.4.1 判断下列矩阵能否进行LU分解，并求矩阵的秩.（1）；（2）；（3）.解 （1）在MATLAB工作窗口输入程序 A=1 2 3;1 12 7;4 5 6;hl=pdLUfj(A)运行后输出结果为请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：RA = 3， hl = 1 10 -48（2）在MATLAB工作窗口输入程序 A=1 2 3;1 2 7;4 5 6;hl=pdLUfj(A)运行后输出结果为请注意：因

15、为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：RA = 3， hl =1 0 12（3）在MATLAB工作窗口输入程序 A=1 2 3;1 2 3;4 5 6;hl=pdLUfj(A)运行后输出结果为请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下RA = 2， hl = 03.4.2 直接LU分解法及其MATLAB程序根据定理3.5和（3.16）式，现提供MATLAB程序如下：将矩阵进行直接LU分解的MATLAB程序输入的量：系数矩阵；输出的量：矩阵的秩RA和各阶顺序主子式） 的值hl及其相关信息，能否进行LU分解的信息

16、.如果能进行LU分解，则输出和.function hl=zhjLU(A)n n =size(A); RA=rank(A); if RA=ndisp(请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下), RA,hl=det(A);returnendif RA=n for p=1:nh(p)=det(A(1:p, 1:p);endhl=h(1:n);for i=1:nif h(1,i)=0disp(请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：), hl;RAreturnendend if h(1,i)=0 disp

17、(请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：)for j=1:nU(1,j)=A(1,j);endfor k=2:nfor i=2:n for j=2:n L(1,1)=1;L(i,i)=1; if ijL(1,1)=1;L(2,1)=A(2,1)/U(1,1); L(i,1)=A(i,1)/U(1,1);L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k)/U(k,k);elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);endendendendhl;RA,U,Lendend例3.4

18、.3 用矩阵进行直接LU分解的MATLAB程序分解矩阵.解 在MATLAB工作窗口输入程序 A=1 0 2 0;0 1 0 1;1 2 4 3;0 1 0 3; hl=zhjLU(A)运行后输出结果L = 1 0 0 0 0 1 0 0 1 2 1 0 0 1 0 1 hl = 1 1 2 4 请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：RA = 4 U = 1 0 2 0 0 1 0 1 0 0 2 1 0 0 0 2 由此可见，例3.4.3和例3.4.2中的LU分解式相同.3.4.3 含交换矩阵的LU分解法及其MATLAB程序含交换

19、矩阵的LU分解法的MATLAB程序是使用列主元消元法进行的，其命令和功能如下：MATLAB命令功 能X U= lu (A) 输出X为一交换矩阵与单位下三角形矩阵L之积，U为上三角形矩阵，使得, 其中.L U P= lu (A) 输出L为单位下三角形矩阵，U为上三角形矩阵，P为一交换矩阵，使 .3.4.4 判断正定对称矩阵的方法及其MATLAB程序根据定理3.5和（3.16）式，现提供MATLAB程序如下：判断矩阵是否是正定对称矩阵的MATLAB程序输入的量：系数矩阵；输出的量： 矩阵的转置矩阵zA和的各阶顺序主子式 的值hl及其相关信息，能否进行LU分解的信息.function hl=zddc

20、(A)n n =size(A);for p=1:nh(p)=det(A(1:p, 1:p);endhl=h(1:n);zA=A;for i=1:n if h(1,i)0disp(请注意：因为A的各阶顺序主子式hl都大于零，所以A是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：) hl;zAend例3.4.5 判断下列矩阵是否是正定对称矩阵：（1）；(2) ； (3) ；（4）.解 （1）在MATLAB工作窗口输入程序 A=0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9;hl=zddc(A)运行后输出结果请注意：因为A的各阶顺序主子式hl不全大于零，所以

21、A不是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：zA = 1/10 -1 11 5 2 2 21 7 3 -3 13 8 4 4 41 9 hl = 1/10 11/5 -1601/10 3696/5 因此,即不是正定矩阵，也不是对称矩阵.（2）在MATLAB工作窗口输入程序 A=1 -1 2 1;-1 3 0 -3;2 0 9 -6;1 -3 -6 19,hl=zddc(A)运行后输出结果A = 1 -1 2 1 -1 3 0 -3 2 0 9 -6 1 -3 -6 19 请注意：因为A的各阶顺序主子式hl都大于零，所以A是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下

22、：zA = 1 -1 2 1 -1 3 0 -3 2 0 9 -6 1 -3 -6 19 hl = 1 2 6 24 （3）在MATLAB工作窗口输入程序 A=1/sqrt(2) -1/sqrt(2) 0 0; -1/sqrt(2) 1/sqrt(2) 0 0; 0 0 1/sqrt(2) -1/sqrt(2); 0 0 -1/sqrt(2) 1/sqrt(2), hl=zddc(A)运行后输出结果A= 985/1393 -985/1393 0 0 -985/1393 985/1393 0 0 0 0 985/1393 -985/1393 0 0 -985/1393 985/1393 请注意：

23、因为A的各阶顺序主子式hl不全大于零，所以A不是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：zA = 985/1393 -985/1393 0 0 -985/1393 985/1393 0 0 0 0 985/1393 -985/1393 0 0 -985/1393 985/1393 hl = 985/1393 0 0 0 可见，不是正定矩阵，是半正定矩阵；因为= T 因此，是对称矩阵.（4）在MATLAB工作窗口输入程序 A=-2 1 1;1 -6 0;1 0 -4;hl=zddc(A)运行后输出结果 A = -2 1 1 1 -6 0 1 0 -4请注意：因为A的各阶顺序主子式

24、hl不全大于零，所以A不是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：zA = -2 1 1 hl = -2 11 -38 1 -6 0 1 0 -4不是正定矩阵，是负定矩阵；因为 = T 因此，是对称矩阵.3.4.5 正定对称矩阵的楚列斯基（Cholesky）分解及其MATLAB程序阶正定或半正定的对称矩阵的楚列斯基分解的MATLAB程序如下：MATLAB命令功 能R=chol(A)输出R为上三角形矩阵U，使R,p =chol(A)输出R为上三角形矩阵U，使 = 如果读者想了解关于chol更多的使用信息，可以在MATLAB工作窗口输入查询命令help chol.3.5 求解线性方

25、程组的LU方法及其MATLAB程序3.5.1 解线性方程组的楚列斯基（Cholesky）分解法及其MATLAB程序例3.5.1 先将矩阵进行楚列斯基分解,然后解矩阵方程，并用其他方法验证.解 根据（3.26）式编写MATLAB程序，然后在工作窗口输入A=1 -1 2 1;-1 3 0 -3; 2 0 9 -6;1 -3 -6 19;b1=1:2:7; b=b1; R=chol(A);C=A-R*R,R1=inv(R);R2=R1; x=R1*R2*b,Rx=Ab-x运行后输出方程组的解和验证结果x = Rx = 1.0e-014 * C = 1.0e-015 * -8.0000 -0.7105

26、 0 0 0 0 0.3333 -0.0833 0 -0.4441 0 0 3.6667 0.2220 0 0 0 0 2.0000 0.1332 0 0 0 0 3.5.2 解线性方程组的直接LU分解法及其MATLAB程序例3.5.2 首先将矩阵直接进行LU分解，然后解矩阵方程，.解 (1) 首先将矩阵直接进行LU分解.在MATLAB工作窗口输入程序 A=1 0 2 0;0 1 0 1;1 2 4 3;0 1 0 3;b=1;2;-1;5; hl=zhjLU(A),A-L*U运行后输出LU分解请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下

27、：L = 1 0 0 0 0 1 0 0 1 2 1 00 1 0 1hl = 1 1 2 4RA = 4U = 1 0 2 0 0 1 0 1 0 0 2 1 0 0 0 2分解为一个单位下三角形矩阵和一个上三角形矩阵的积 .(2) 然后根据（3.27）式编写MATLAB程序，然后在工作窗口输入 U=1 0 2 0;0 1 0 1;0 0 2 1;0 0 0 2; L=1 0 0 0;0 1 0 0;1 2 1 0;0 1 0 1;b=1;2;-1;5;U1=inv(U); L1=inv(L); X=U1*L1*b,x=Ab运行后输出方程组的解X = 8.50000000000000 0.5

28、0000000000000 -3.75000000000000 1.500000000000003.5.3 解线性方程组的选主元的LU方法及其MATLAB程序例3.5.3 先将矩阵进行LU分解，然后解矩阵方程 其中，.解 方法1 根据（3.28）式编写MATLAB程序，然后在工作窗口输入 A=0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9; b=1;2;-1;5; L U P=LU(A), U1=inv(U); L1=inv(L); X=U1* L1*P*bP = 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 1X =-1.2013 3.3677

29、0.0536 -1.4440运行后输出结果L = 1.0000 0 0 0 -0.0909 1.0000 0 0 0.0091 0.4628 1.0000 0 0.4545 -0.6512 0.2436 1.0000U =11.0000 21.0000 13.0000 41.0000 0 3.9091 -1.8182 7.7273 0 0 3.7233 0.05120 0 0 -4.6171方法2 根据（3.29）式编写MATLAB程序，然后在工作窗口输入 A=0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9;b=1;2;-1;5; F U=LU(A), U1=i

30、nv(U); F1=inv(F); X=U1*F1*bU=11.0000 21.0000 13.0000 41.0000 0 3.9091 -1.8182 7.7273 0 0 3.7233 0.0512 0 0 0 -4.6171运行后输出结果F=0.0091 0.4628 1.0000 0 -0.0909 1.0000 0 0 1.0000 0 0 0 0.4545 -0.6512 0.2436 1.0000X =-1.2013 3.3677 0.0536 -1.4440解线性方程组的选主元的LU方法，现提供MATLAB程序如下：用LU分解法解线性方程组的MATLAB程序输入的量：系数矩阵

31、A和常数向量b；输出的量：系数矩阵A和增广矩阵B的秩RA,RB, 方程组中未知量的个数n和方程组解及其有关的信息.function RA,RB,n,X,Y=LUjfcz(A,b)n n =size(A);B=A b; RA=rank(A); RB=rank(B); for p=1:nh(p)=det(A(1:p, 1:p);endhl=h(1:n);for i=1:nif h(1,i)=0disp(请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：) hl;RAreturnendendif h(1,i)=0 disp(请注意：因为A的各阶主子式

32、都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：)X=zeros(n,1); Y=zeros(n,1); C=zeros(1,n);r=1:1;for p=1:n-1max1,j=max(abs(A(p:n,p); C=A(p,:); A(p,:)= A(j+P1,:); C= A(j+P1,:);g=r(p); r(p)= r(j+P1); r(j+P1)=g;for k=p+1:nH= A(k,p)/A(p,p); A(k,p) = H; A(k,p+1:n)=A(k,p+1:n)- H* A(p,p+1:n);endendY(1)=B(r(1);for k=2

33、:nY(k)= B(r(k)- A(k,1:k-1)* Y(1:k-1); endX(n)= Y(n)/ A(n,n);for i=n-1:-1:1 X(i)= (Y(i)- A(i, i+1:n) * X (i+1:n)/ A(i,i);endendRA,RB,n,X,Y;3.6 误差分析及其两种MATLAB程序对于实际问题导出的线性方程组 的系数矩阵和右端向量往往带有误差（扰动），下面讨论当或的微小变化对解的影响.3.6.1 用MATLAB软件作误差分析例3.6.2 解下列矩阵方程，并比较方程（1）和（2）有何区别，它们的解有何变化.其中解 (1) 矩阵方程的系数矩阵为7阶希尔伯特（Hil

34、bert）矩阵，我们可以用下列命令计算阶希尔伯特矩阵 h=hilb(n) % 输出h为n阶Hilbert矩阵在MATLAB工作窗口输入程序 A=hilb(7);b=1;2;2;2;2;2;2;X=Ab运行后输出的解为 X =（-35 504 -1260 -4200 20790 -27720 12012. （2）在MATLAB工作窗口输入程序 B =0.001,zeros(1,6);zeros(6,1),zeros(6,6);A=(B+hilb(7); b=1;2;2;2;2;2;2;X=Ab 运行后输出方程的解为 X =（-33 465 -966 -5181 22409 -29015 1241

35、3. 在MATLAB工作窗口输入程序 X =-33, 465,-966,-5181,22409,-29015,12413; X1 =-35,504,-1260,-4200,20790,-27720,12012; wu=X1- X运行后输出方程（1）和（2）的解的误差为 .方程（1）和（2）的系数矩阵的差为，常数向量相同，则的解的差为.的微小变化，引起的很大变化，即对的扰动是敏感的.3.6.2 求P 条件数和讨论解的性态的MATLAB程序根据定理3.10可知，判断线性方程组（3.38）是否病态需要计算条件数，在第一章中我们已经给出了计算矩阵的2范数条件数、1范数条件数、范数条件数、弗罗贝尼乌斯（Frobenius）范数条件数、条件数倒数等MATLAB命令.用这些命令计算各种条件数非常方便.现提供求P条件数和讨论解的性态的MATLAB程序.求P条件数和讨论解的性态的MATLAB程序输入的量：