第三章解线性方程组的直接方法.doc
蜗牛文库
-高品质阅读，高比例分成-
下载文档
温馨提示：
1. 部分包含数学公式或PPT动画的文件，查看预览时可能会显示错乱或异常，文件下载后无此问题，请放心下载。
2. 本文档由用户上传，版权归属用户，汇文网负责整理代发布。如果您对本文档版权有争议请及时联系客服。
3. 下载前请仔细阅读文档内容，确认文档内容符合您的需求后进行下载，若出现内容与标题不符可向本站投诉处理。
4. 下载文档时可能由于网络波动等原因无法下载或下载错误，付费完成后未能成功下载的用户请联系客服处理。
网站客服：3074922707
1. 部分包含数学公式或PPT动画的文件，查看预览时可能会显示错乱或异常，文件下载后无此问题，请放心下载。
2. 本文档由用户上传，版权归属用户，汇文网负责整理代发布。如果您对本文档版权有争议请及时联系客服。
3. 下载前请仔细阅读文档内容，确认文档内容符合您的需求后进行下载，若出现内容与标题不符可向本站投诉处理。
4. 下载文档时可能由于网络波动等原因无法下载或下载错误，付费完成后未能成功下载的用户请联系客服处理。
网站客服：3074922707
第三章
解线性方程组的直接方法
第三
线性方程组
直接
方法
5.1 方程组的逆矩阵解法及其MATLAB程序
5.1.1 线性方程组有解的判定条件及其MATLAB程序
判定线性方程组是否有解的MATLAB程序
function [RA,RB,n]=jiepb(A,b)
B=[A b];n=length(b); RA=rank(A);
RB=rank(B);zhica=RB-RA;
if zhica>0,
disp('请注意：因为RA~=RB，所以此方程组无解.')
return
end
if RA==RB
if RA==n
disp('请注意：因为RA=RB=n，所以此方程组有唯一解.')
else
disp('请注意：因为RA=RB<n，所以此方程组有无穷多解.')
end
end
例5.1.1 判断下列线性方程组解的情况.如果有唯一解，则求解.
(1) (2)
(3) (4)
解 在MATLAB工作窗口输入程序
>> 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，所以此方程组有唯一解.
RA = 4,RB =4,n =4
在MATLAB工作窗口输入
>>X=A\b,
运行后输出结果为 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<n，所以此方程组有无穷多解.
RA =2,RB =2,n =4
(3) 在MATLAB工作窗口输入程序
>> 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 1;4 2 -2 1;2 1 -1 -1];
b=[1; 2; 1]; [RA,RB,n]=jiepb(A,b)
运行后输出结果
请注意：因为RA=RB<n，所以此方程组有无穷多解.
RA =2,RB =2,n =3
5.2 三角形方程组的解法及其MATLAB程序
5.2.1 解三角形方程组的MATLAB程序
解上三角形线性方程组的MATLAB程序
function [RA,RB,n,X]=shangsan(A,b)
B=[A b]; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA;
if zhica>0,
disp('请注意：因为RA~=RB，所以此方程组无解.')
return
end
if RA==RB
if RA==n
disp('请注意：因为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=RB<n，所以此方程组有无穷多解.')
end
end
例5.2.1 用解上三角形线性方程组的MATLAB程序解方程组
.
解 在MATLAB工作窗口输入程序
>>A=[5 -1 2 3;0 -2 7 -4;0 0 6 5;0 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.0]’
5.3 高斯（Gauss）消元法和列主元消元法及其MATLAB程序
5.3.1 高斯消元法及其MATLAB程序
用高斯消元法解线性方程组的MATLAB程序
function [RA,RB,n,X]=gaus(A,b)
B=[A b]; n=length(b); RA=rank(A);
RB=rank(B);zhica=RB-RA;
if zhica>0,
disp('请注意：因为RA~=RB，所以此方程组无解.')
return
end
if RA==RB
if RA==n
disp('请注意：因为RA=RB=n，所以此方程组有唯一解.')
X=zeros(n,1); C=zeros(1,n+1);
for p= 1:n-1
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);
end
end
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);
end
else
disp('请注意：因为RA=RB<n，所以此方程组有无穷多解.')
end
end
例5.3.1 用高斯消元法和MATLAB程序求解下面的非齐次线性方程组，并且用逆矩阵解方程组的方法验证.
解 在MATLAB工作窗口输入程序
>> 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
0
RA =
4
RB =
4
n =
4
5.3.2 列主元消元法及其MATLAB程序
用列主元消元法解线性方程组的MATLAB程序
function [RA,RB,n,X]=liezhu(A,b)
B=[A b]; n=length(b); RA=rank(A);
RB=rank(B);zhica=RB-RA;
if zhica>0,
disp('请注意：因为RA~=RB，所以此方程组无解.')
return
end
if RA==RB
if RA==n
disp('请注意：因为RA=RB=n，所以此方程组有唯一解.')
X=zeros(n,1); C=zeros(1,n+1);
for p= 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;
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);
end
end
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);
end
else
disp('请注意：因为RA=RB<n，所以此方程组有无穷多解.')
end
end
例5.3.2 用列主元消元法解线性方程组的MATLAB程序解方程组
.
解 在MATLAB工作窗口输入程序
>> 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 0]’
5.4 LU分解法及其MATLAB程序
5.4.1判断矩阵LU分解的充要条件及其MATLAB程序
判断矩阵能否进行LU分解的MATLAB程序
function hl=pdLUfj(A)
[n n] =size(A); RA=rank(A);
if RA~=n
disp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下：'), RA,hl=det(A); return
end
if RA==n
for p=1:n,h(p)=det(A(1:p, 1:p));, end
hl=h(1:n);
for i=1:n
if h(1,i)==0
disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：'),hl;RA,return
end
end
if h(1,i)~=0
disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：')
hl;RA
end
end
例5.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)
运行后输出结果为
请注意：因为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 = 0
5.4.2 直接LU分解法及其MATLAB程序
将矩阵进行直接LU分解的MATLAB程序
function hl=zhjLU(A)
[n n] =size(A); RA=rank(A);
if RA~=n
disp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下:'), RA,hl=det(A);
return
end
if RA==n
for p=1:n
h(p)=det(A(1:p, 1:p));
end
hl=h(1:n);
for i=1:n
if h(1,i)==0
disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：'), hl;RA
return
end
end
if h(1,i)~=0
disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：')
for j=1:n
U(1,j)=A(1,j);
end
for k=2:n
for i=2:n
for j=2:n
L(1,1)=1;L(i,i)=1;
if i>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
例5.4.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
5.4.4 判断正定对称矩阵的方法及其MATLAB程序
判断矩阵是否是正定对称矩阵的MATLAB程序
function hl=zddc(A)
[n n] =size(A);
for p=1:n
h(p)=det(A(1:p, 1:p));
end
hl=h(1:n);zA=A';
for i=1:n
if h(1,i)<=0
disp('请注意：因为A的各阶顺序主子式hl不全大于零，所以A不是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：'), hl;zA,return
end
end
if h(1,i)>0
disp('请注意：因为A的各阶顺序主子式hl都大于零，所以A是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：')
hl;zA
end
例5.4.4 判断下列矩阵是否是正定对称矩阵：
（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不是对称矩阵
请注意：因为A的各阶顺序主子式hl不全大于零，所以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是对称矩阵
请注意：因为A的各阶顺序主子式hl都大于零，所以A是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：
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
请注意： A是对称矩阵
请注意：因为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是对称矩阵
请注意：因为A的各阶顺序主子式hl不全大于零，所以A不是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：
zA = -2 1 1 hl = -2 11 -38
1 -6 0
1 0 -4
可见不是正定矩阵，是负定矩阵；因为 = T 因此，是对称矩阵.
5.5 求解线性方程组的LU方法及其MATLAB程序
5.5.1 解线性方程组的直接LU分解法及其MATLAB程序
例5.5.1 首先将矩阵直接进行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依次如下：
L = 1 0 0 0
0 1 0 0
1 2 1 0
0 1 0 1
hl = 1 1 2 4
RA = 4
U = 1 0 2 0
0 1 0 1
0 0 2 1
0 0 0 2
分解为一个单位下三角形矩阵和一个上三角形矩阵的积 .
(2)在工作窗口输入
>> 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=A\b
运行后输出方程组的解
X = 8.50000000000000
0.50000000000000
-3.75000000000000
1.50000000000000
5.5.2 解线性方程组的选主元的LU方法及其MATLAB程序
例5.5.2 先将矩阵进行LU分解，然后解矩阵方程 其中
，.
解 方法1 编写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*b
P = 0 0 1 0
0 1 0 0
1 0 0 0
0 0 0 1
X =[-1.2013 3.3677 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.0000
U =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
方法2 编写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=inv(U); F1=inv(F); X=U1*F1*b
U=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.0000
X =[-1.2013 3.3677 0.0536 -1.4440]’
用LU分解法解线性方程组的MATLAB程序
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:n
h(p)=det(A(1:p, 1:p));
end
hl=h(1:n);
for i=1:n
if h(1,i)==0
disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：')
hl;RA
return
end
end
if h(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;
for p=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;
for k=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));
for k=2:n
Y(k)= B(r(k))- A(k,1:k-1)* Y(1:k-1);
end
X(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);
end
end
[RA,RB,n,X,Y]’;
33．
