第二章 解线性方程组的直接方法的MATLAB程序
?x2?x3?x4?0,??x?x?x?3x?1,?1234. ??2x1?2x2?4x3?6x4??1,??x1?2x2?4x3?x4??1.解 在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]’
2.4 LU分解法及其MATLAB程序
2.4.1判断矩阵LU分解的充要条件及其MATLAB程序 判断矩阵A能否进行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
例2.4.1 判断下列矩阵能否进行LU分解，并求矩阵的秩.
?123??123??123???????（1）?1127?；（2）?127?；（3）?123?.
?456??456??456???????解 （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依次如下：
29．
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
2.4.2 直接LU分解法及其MATLAB程序 将矩阵A进行直接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
例2.4.3 用矩阵进行直接LU分解的MATLAB程序分解矩阵
第二章 解线性方程组的直接方法的MATLAB程序
?1020????0101?. A??1243????0103???解 在MATLAB工作窗口输入程序
>> A=[1 0 2 0;0 1 0 1;1 2 4 3;0 1 0 3]; hl=zhjLU(A)
运行后输出结果
请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA
和各阶顺序主子式值hl依次如下：
L = 1 0 0 0 RA = 4
0 1 0 0 U = 1 0 2 0
1 2 1 0 0 1 0 1
0 1 0 1 0 0 2 1
0 0 0 2 hl = 1 1 2 4
2.4.4 判断正定对称矩阵的方法及其MATLAB程序 判断矩阵A是否是正定对称矩阵的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
例2.4.5 判断下列矩阵是否是正定对称矩阵：
?0.1234??1?121?1?1???211??????00???2?；（1）??12?34?；(2) ??130?3?； (3) ?2（4）?1?60??. 1?1??11211341??209?6?00????10?4?22?????????5789??1?3?619?11??0????0?????0?0?2122?1???2?解 （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
因此,A即不是正定矩阵，也不是对称矩阵.
（2）在MATLAB工作窗口输入程序
>> A=[1 -1 2 1;-1 3 0 -3;2 0 9 -6;1 -3 -6 19],hl=zddc(A)
31．
运行后输出结果
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
可见，A不是正定矩阵，是半正定矩阵；因为A= A 因此，A是对称矩阵.
（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
可见A不是正定矩阵，是负定矩阵；因为A = A 因此，A是对称矩阵.
2.5 求解线性方程组的LU方法及其MATLAB程序
2.5.1 解线性方程组的楚列斯基（Cholesky）分解法及其MATLAB程序
例3.5.1 先将矩阵A进行楚列斯基分解,然后解矩阵方程AX?b，并用其他方法验证.
1??1?12?1??????130?3???3?. A??,b??5?209?6??????1?3?619??7?????第二章 解线性方程组的直接方法的MATLAB程序
解 在工作窗口输入
>>A=[1 -1 2 1;-1 3 0 -3; 2 0 9 -6;1 -3 -6 19];
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.7105 0 0 0 0 0.3333 -0.0833 0 -0.4441 0 0 3.6667 0.2220 0 0 0 0 2.0000 0.1332 0 0 0 0
2.5.2 解线性方程组的直接LU分解法及其MATLAB程序
例3.5.2 首先将矩阵A直接进行LU分解，然后解矩阵方程AX?b
?1??1020?????2?. ?0101?，b??A????1?1243??????5??0103?????解 (1) 首先将矩阵A直接进行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 RA = 4
0 1 0 0 U = 1 0 2 0
1 2 1 0 0 1 0 1
0 1 0 1 0 0 2 1
hl = 1 1 2 4 0 0 0 2
A分解为一个单位下三角形矩阵L和一个上三角形矩阵U的积 A?LU.
(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.000 0.000 -3.000 1.000
2.5.3 解线性方程组的选主元的LU方法及其MATLAB程序
例3.5.3 先将矩阵A进行LU分解，然后解矩阵方程AX?b 其中
34??0.12?1??????12?34??，?2?. A??b???1?11211341??????5??5?789????解 方法1 根据（3.28）式编写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
运行后输出结果
L = 1.0000 0 0 0 P = 0 0 1 0
0 1 0 0
1 0 0 0 33．
0 0 0 1
X =[-1.2013 3.3677 0.0536
-1.4440]’