第二章解线性方程组的直接方法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=RBend
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=RBRA=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=RBend
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]’；
升级会员 