第三章解线性方程组的直接方法.docx
- 文档编号:8497782
- 上传时间:2023-01-31
- 格式:DOCX
- 页数:27
- 大小:180.04KB
第三章解线性方程组的直接方法.docx
《第三章解线性方程组的直接方法.docx》由会员分享,可在线阅读,更多相关《第三章解线性方程组的直接方法.docx(27页珍藏版)》请在冰豆网上搜索。
第三章解线性方程组的直接方法
在这章中我们要学习线性方程组的直接法,特别是适合用数学软件在计算机上求解的方法.
3.1方程组的逆矩阵解法及其MATLAB程序
3.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 例3.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 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 3.2三角形方程组的解法及其MATLAB程序 3.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 end end 例3.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]’ 3.3高斯(Gauss)消元法和列主元消元法及其MATLAB程序 3.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 例3.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 3.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 例3.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]’ 3.4LU分解法及其MATLAB程序 3.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 例3.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 3.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 例3.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 3.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 例3.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因此, 是对称矩阵. 3.5求解线性方程组的LU方法及其MATLAB程序 3.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 3.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 3.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]’; 3.6误差分析及其两种MATLAB程序 3.6.1
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第三章 解线性方程组的直接方法 第三 线性方程组 直接 方法