MATLAB解线性方程组的直接方法.docx
《MATLAB解线性方程组的直接方法.docx》由会员分享,可在线阅读,更多相关《MATLAB解线性方程组的直接方法.docx(31页珍藏版)》请在冰豆网上搜索。
MATLAB解线性方程组的直接方法
在这章中我们要学习线性方程组的直接法,特别是适合用数学软件在计算机上求解的方法.
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=RBend
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=RBRA=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=RBRA=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=RBend
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=RBend
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=RBend
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.