第二章解线性方程组的直接方法matlab用法.docx
- 文档编号:29372046
- 上传时间:2023-07-22
- 格式:DOCX
- 页数:35
- 大小:167.92KB
第二章解线性方程组的直接方法matlab用法.docx
《第二章解线性方程组的直接方法matlab用法.docx》由会员分享,可在线阅读,更多相关《第二章解线性方程组的直接方法matlab用法.docx(35页珍藏版)》请在冰豆网上搜索。
第二章解线性方程组的直接方法matlab用法
第二章解线性方程组的直接方法
在这章中我们要学习线性方程组的直接法,特别是适合用数学软件在计算机上求解的方法。
2.1方程组的逆矩阵解法及其MATLAB程序
2。
1。
3线性方程组有解的判定条件及其MATLAB程序
判定线性方程组
是否有解的MATLAB程序
function[RA,RB,n]=jiepb(A,b)
B=[Ab];n=length(b);RA=rank(A);
RB=rank(B);zhica=RB-RA;
ifzhica〉0,
disp('请注意:
因为RA~=RB,所以此方程组无解。
')
return
end
ifRA==RB
ifRA==n
disp(’请注意:
因为RA=RB=n,所以此方程组有唯一解.')
else
disp(’请注意:
因为RA=RB end end 例2。 1.4判断下列线性方程组解的情况.如果有唯一解,则用表3—2方法求解。 (1) (2) (3) (4) 解在MATLAB工作窗口输入程序 〉〉A=[23-15;312-7;41-36;1-24—7]; b=[0;0;0;0];[RA,RB,n]=jiepb(A,b) 运行后输出结果为 请注意: 因为RA=RB=n,所以此方程组有唯一解。 RA=4,RB=4,n=4 在MATLAB工作窗口输入 〉〉X=A\b, 运行后输出结果为X=(0000)’. (2)在MATLAB工作窗口输入程序 〉>A=[34—57;2-33-2;411—1316;7-213];b=[0;0;0;0]; [RA,RB,n]=jiepb(A,b) 运行后输出结果 请注意: 因为RA=RB〈n,所以此方程组有无穷多解。 RA=2,RB=2,n=4 (3) 在MATLAB工作窗口输入程序 〉〉A=[42-1;3—12;1130];b=[2;10;8];[RA,RB,n]=jiepb(A,B) 运行后输出结果 请注意: 因为RA~=RB,所以此方程组无解。 RA=2,RB=3,n=3 (4)在MATLAB工作窗口输入程序 >〉A=[21-11;42-21;21—1-1]; b=[1;2;1];[RA,RB,n]=jiepb(A,b) 运行后输出结果 请注意: 因为RA=RB RA=2,RB=2,n=3 2。 2三角形方程组的解法及其MATLAB程序 2。 2。 2解三角形方程组的MATLAB程序 解上三角形线性方程组 的MATLAB程序 function[RA,RB,n,X]=shangsan(A,b) B=[Ab];n=length(b);RA=rank(A);RB=rank(B);zhica=RB—RA; ifzhica〉0, disp(’请注意: 因为RA~=RB,所以此方程组无解。 ') return end ifRA==RB ifRA==n disp('请注意: 因为RA=RB=n,所以此方程组有唯一解。 ’) X=zeros(n,1);X(n)=b(n)/A(n,n); fork=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=RB〈n,所以此方程组有无穷多解.') end end 例2。 2。 2用解上三角形线性方程组的MATLAB程序解方程组 . 解在MATLAB工作窗口输入程序 >〉A=[5—123;0-27—4;0065;0003]; 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。 02。 0]’ 2。 3高斯(Gauss)消元法和列主元消元法及其MATLAB程序 2。 3。 1高斯消元法及其MATLAB程序 用高斯消元法解线性方程组 的MATLAB程序 function[RA,RB,n,X]=gaus(A,b) B=[Ab];n=length(b);RA=rank(A); RB=rank(B);zhica=RB—RA; ifzhica>0, disp(’请注意: 因为RA~=RB,所以此方程组无解。 ') return end ifRA==RB ifRA==n disp('请注意: 因为RA=RB=n,所以此方程组有唯一解。 ') X=zeros(n,1);C=zeros(1,n+1); forp=1: n-1 fork=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); end end b=B(1: n,n+1);A=B(1: n,1: n);X(n)=b(n)/A(n,n); forq=n-1: -1: 1 X(q)=(b(q)-sum(A(q,q+1: n)*X(q+1: n)))/A(q,q); end else disp(’请注意: 因为RA=RB end end 例2.3.2用高斯消元法和MATLAB程序求解下面的非齐次线性方程组,并且用逆矩阵解方程组的方法验证. 解在MATLAB工作窗口输入程序 〉〉A=[1—11-3;0-1-11;2—2—46;1—2-41]; b=[1;0;—1;—1];[RA,RB,n,X]=gaus(A,b) 运行后输出结果 请注意: 因为RA=RB=n,所以此方程组有唯一解. X= 0 -0。 5000 0。 5000 0 RA= 4 RB= 4 n= 4 2。 3.2列主元消元法及其MATLAB程序 用列主元消元法解线性方程组 的MATLAB程序 function[RA,RB,n,X]=liezhu(A,b) B=[Ab];n=length(b);RA=rank(A); RB=rank(B);zhica=RB—RA; ifzhica>0, disp('请注意: 因为RA~=RB,所以此方程组无解.’) return end ifRA==RB ifRA==n disp(’请注意: 因为RA=RB=n,所以此方程组有唯一解.’) X=zeros(n,1);C=zeros(1,n+1); forp=1: n-1 [Y,j]=max(abs(B(p: n,p)));C=B(p,: ); B(p,: )=B(j+p—1,: );B(j+p-1,: )=C; fork=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); end end b=B(1: n,n+1);A=B(1: n,1: n);X(n)=b(n)/A(n,n); forq=n-1: —1: 1 X(q)=(b(q)—sum(A(q,q+1: n)*X(q+1: n)))/A(q,q); end else disp(’请注意: 因为RA=RB ’) end end 例2.3.3用列主元消元法解线性方程组的MATLAB程序解方程组 。 解在MATLAB工作窗口输入程序 >〉A=[0—1—11;1—11-3;2-2-46;1-2-41]; b=[0;1;-1;—1];[RA,RB,n,X]=liezhu(A,b) 运行后输出结果 请注意: 因为RA=RB=n,所以此方程组有唯一解. RA=4,RB=4,n=4,X=[0—0。 50。 50]' 2.4LU分解法及其MATLAB程序 2.4。 1判断矩阵LU分解的充要条件及其MATLAB程序 判断矩阵 能否进行LU分解的MATLAB程序 functionhl=pdLUfj(A) [nn]=size(A);RA=rank(A); ifRA~=n disp(’请注意: 因为A的n阶行列式hl等于零,所以A不能进行LU分解.A的秩RA如下: ’),RA,hl=det(A);return end ifRA==n forp=1: n,h(p)=det(A(1: p,1: p));,end hl=h(1: n); fori=1: n ifh(1,i)==0 disp('请注意: 因为A的r阶主子式等于零,所以A不能进行LU分解。 A的秩RA和各阶顺序主子式值hl依次如下: ’),hl;RA,return end end ifh(1,i)~=0 disp(’请注意: 因为A的各阶主子式都不等于零,所以A能进行LU分解。 A的秩RA和各阶顺序主子式值hl依次如下: ’) hl;RA end end 例2。 4。 1判断下列矩阵能否进行LU分解,并求矩阵的秩. (1) ; (2) ;(3) . 解 (1)在MATLAB工作窗口输入程序 〉>A=[123;1127;456];hl=pdLUfj(A) 运行后输出结果为 请注意: 因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下: RA=3,hl=110—48 (2)在MATLAB工作窗口输入程序 〉>A=[123;127;456];hl=pdLUfj(A) 运行后输出结果为 请注意: 因为A的r阶主子式等于零,所以A不能进行LU分解。 A的秩RA和各阶顺序主子式值hl依次如下: RA=3,hl=1012 (3)在MATLAB工作窗口输入程序 >>A=[123;123;456];hl=pdLUfj(A) 运行后输出结果为 请注意: 因为A的n阶行列式hl等于零,所以A不能进行LU分解。 A的秩RA如下 RA=2,hl=0 2.4.2直接LU分解法及其MATLAB程序 将矩阵 进行直接LU分解的MATLAB程序 functionhl=zhjLU(A) [nn]=size(A);RA=rank(A); ifRA~=n disp('请注意: 因为A的n阶行列式hl等于零,所以A不能进行LU分解。 A的秩RA如下: '),RA,hl=det(A); return end ifRA==n forp=1: n h(p)=det(A(1: p,1: p)); end hl=h(1: n); fori=1: n ifh(1,i)==0 disp('请注意: 因为A的r阶主子式等于零,所以A不能进行LU分解。 A的秩RA和各阶顺序主子式值hl依次如下: ’),hl;RA return end end ifh(1,i)~=0 disp('请注意: 因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下: ') forj=1: n U(1,j)=A(1,j); end fork=2: n fori=2: n forj=2: n L(1,1)=1;L(i,i)=1; ifi>j L(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); else U(k,j)=A(k,j)—L(k,1: k—1)*U(1: k—1,j); end end end end hl;RA,U,L end end 例2。 4.3用矩阵进行直接LU分解的MATLAB程序分解矩阵 。 解在MATLAB工作窗口输入程序 >>A=[1020;0101;1243;0103];hl=zhjLU(A) 运行后输出结果 L=1000 0100 1210 0101 hl=1124 请注意: 因为A的各阶主子式都不等于零,所以A能进行LU分解。 A的秩RA和各阶顺序主子式值hl依次如下: RA=4 U=1020 0101 0021 0002 2.4.4判断正定对称矩阵的方法及其MATLAB程序 判断矩阵 是否是正定对称矩阵的MATLAB程序 functionhl=zddc(A) [nn]=size(A); forp=1: n h(p)=det(A(1: p,1: p)); end hl=h(1: n);zA=A'; fori=1: n ifh(1,i)〈=0 disp(’请注意: 因为A的各阶顺序主子式hl不全大于零,所以A不是正定的。 A的转置矩阵zA和各阶顺序主子式值hl依次如下: ’),hl;zA,return end end ifh(1,i)>0 disp('请注意: 因为A的各阶顺序主子式hl都大于零,所以A是正定的。 A的转置矩阵zA和各阶顺序主子式值hl依次如下: ') hl;zA end 例2。 4。 5判断下列矩阵是否是正定对称矩阵: (1) ; (2) ;(3) ;(4) . 解 (1)在MATLAB工作窗口输入程序 〉〉A=[0.1234;—12—34;11211341;5789];hl=zddc(A) 运行后输出结果 请注意: A不是对称矩阵 请注意: 因为A的各阶顺序主子式hl不全大于零,所以A不是正定的。 A的转置矩阵zA和各阶顺序主子式值hl依次如下: zA=1/10-1115 22217 3—3138 44419 hl=1/1011/5—1601/103696/5 因此, 即不是正定矩阵,也不是对称矩阵。 (2)在MATLAB工作窗口输入程序 〉>A=[1—121;-130-3;209-6;1-3—619],hl=zddc(A) 运行后输出结果 A=1—121 -130—3 209—6 1-3—619 请注意: A是对称矩阵 请注意: 因为A的各阶顺序主子式hl都大于零,所以A是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下: zA=1-121 —130—3 209-6 1—3-619 hl=12624 (3)在MATLAB工作窗口输入程序 >〉A=[1/sqrt (2)-1/sqrt (2)00;-1/sqrt (2)1/sqrt (2)00;001/sqrt (2)—1/sqrt (2);00—1/sqrt (2)1/sqrt (2)],hl=zddc(A) 运行后输出结果 A=985/1393-985/139300 —985/1393985/139300 00985/1393-985/1393 00-985/1393985/1393 请注意: A是对称矩阵 请注意: 因为A的各阶顺序主子式hl不全大于零,所以A不是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下: zA=985/1393-985/139300 -985/1393985/139300 00985/1393—985/1393 00—985/1393985/1393 hl=985/1393000 可见, 不是正定矩阵,是半正定矩阵;因为 = T因此, 是对称矩阵。 (4)在MATLAB工作窗口输入程序 〉>A=[—211;1-60;10-4];hl=zddc(A) 运行后输出结果 A=—211 1-60 10-4 请注意: A是对称矩阵 请注意: 因为A的各阶顺序主子式hl不全大于零,所以A不是正定的。 A的转置矩阵zA和各阶顺序主子式值hl依次如下: zA=—211hl=—211-38 1—60 10—4 可见 不是正定矩阵,是负定矩阵;因为 = T因此, 是对称矩阵. 2.5求解线性方程组的LU方法及其MATLAB程序 2.5。 1解线性方程组的楚列斯基(Cholesky)分解法及其MATLAB程序 例3。 5。 1先将矩阵 进行楚列斯基分解,然后解矩阵方程 ,并用其他方法验证。 。 解在工作窗口输入 >〉A=[1—121;-130—3;209—6;1-3-619]; b1=1: 2: 7;b=b1';R=chol(A);C=A-R’*R,R1=inv(R);R2=R1'; x=R1*R2*b,Rx=A\b—x 运行后输出方程组的解和验证结果 x=Rx=1.0e-014*C=1.0e-015* -8。 0000-0。 71050000 0。 3333—0。 08330—0。 444100 3.66670.22200000 2。 00000。 13320000 2.5.2解线性方程组的直接LU分解法及其MATLAB程序 例3.5。 2首先将矩阵 直接进行LU分解,然后解矩阵方程 , 。 解 (1)首先将矩阵 直接进行LU分解.在MATLAB工作窗口输入程序 〉〉A=[1020;0101;1243;0103];b=[1;2;—1;5];hl=zhjLU(A),A-L*U 运行后输出LU分解 请注意: 因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下: L=1000 0100 1210 0101 hl=1124 RA=4 U=1020 0101 0021 0002 分解为一个单位下三角形矩阵 和一个上三角形矩阵 的积 。 (2)在工作窗口输入 〉>U=[1020;0101;0021;0002];L=[1000;0100;1210;0101]; b=[1;2;-1;5];U1=inv(U);L1=inv(L);X=U1*L1*b,x=A\b 运行后输出方程组的解 X=8.50000000000000 0.50000000000000 —3。 75000000000000 1.50000000000000 2。 5。 3解线性方程组的选主元的LU方法及其MATLAB程序 例3。 5。 3先将矩阵 进行LU分解,然后解矩阵方程 其中 , 。 解方法1根据(3。 28)式编写MATLAB程序,然后在工作窗口输入 〉>A=[0.1234;-12—34;11211341;5789]; b=[1;2;-1;5];[LUP]=LU(A),U1=inv(U);L1=inv(L);X=U1*L1*P*b P=0010 0100 1000 0001 X=[—1。 20133。 36770。 0536—1。 4440]’ 运行后输出结果 L=1.0000000 —0。 09091.000000 0。 00910。 46281.00000 0.4545-0.65120.24361.0000 U=11.000021。 000013.000041.0000 03。 9091-1。 81827.7273 003.72330。 0512 000-4。 6171 方法2根据(3。 29)式编写MATLAB程序,然后在工作窗口输入 〉〉A=[0。 1234;—12-34;11211341;5789]; b=[1;2;-1;5];[FU]=LU(A),U1=inv(U);F1=inv(F);X=U1*F1*b U=11.000021。 000013。 000041。 0000 03.9091-1。 81827。 7273 003。 72330。 0512 000—4.6171 运行后输出结果 F=0.00910。 46281.00000 -0。 09091。 000000 1。 0000000 0.4545-0.65120。 24361.0000 X=[—1。 20133。 36770。 0536-1。 4440]' 用LU分解法解线性方程组 的MATLAB程序 function[RA,RB,n,X,Y]=LUjfcz(A,b) [nn]=size(A);B=[Ab];RA=rank(A);RB=rank(B); forp=1: n h(p)=det(A(1: p,1: p)); end hl=h(1: n); fori=1: n ifh(1,i)==0 disp(’请注意: 因为A的r阶主子式等于零,所以A不能进行LU分解。 A的秩RA和各阶顺序主子式值hl依次如下: ') hl;RA return end end ifh(1,i)~=0 disp('请注意: 因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下: ’) X=zeros(n,1);Y=zeros(n,1);C=zeros(1,n);r=1: 1; forp=1: n-1 [max1,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; fork=p+1: n H=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); end end Y (1)=B(r (1)); fork=2: n Y(k)=B(r(k))—A(k,1: k—1)*Y(1: k-1); end X(n)=Y(n)/A(n,n); fori=n—1: —1: 1 X(i)=(Y(i)-A(i,i+1: n)*X(i+1: n))/A(i,i); end end [RA,RB,n,X,Y]’;
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第二 线性方程组 直接 方法 matlab 用法