高斯消元法,列主元素消元法及LU分解法的matlab程序.doc

§2.2.1高斯消元法的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
运行命令及结果
a=[2.51 1.48 4.53;1.48 0.93 -1.30 ;2.68 3.04 -1.48];
b=[0.05;1.03;-0.53];[RA,RB,n,X] =gaus (A,b)
请注意：因为RA=RB=n，所以此方程组有唯一解.

RA =

3

RB =

3

n =

3

X =

1.4531
-1.5892
-0.2749

§2.2.2 列主元素消元法的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
运行命令及结果
a=[2.51 1.48 4.53;1.48 0.93 -1.30 ;2.68 3.04 -1.48];
b=[0.05;1.03;-0.53];[RA,RB,n,X]=liezhu(A,b)
请注意：因为RA=RB=n，所以此方程组有唯一解.

RA =

3

RB =

3

n =

3

X =

1.4531
-1.5892
-0.2749

§2.2.3 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
运行命令及结果
a=[2.51 1.48 4.53;1.48 0.93 -1.30 ;2.68 3.04 -1.48];hl=zhjLU(A)
请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：

RA =

3

U =

2.5100 1.4800 4.5300
0 0.9300 -3.9711
0 0 -0.0837

L =

1.0000 0 0
0.5896 1.0000 0
1.0677 1.5696 1.0000

hl =

2.5100 0.1439 13.6410

>> U=[2.5100 1.4800 4.5300
0 0.9300 -3.9711
0 0 -0.0837];
>>

L= [1.0000 0 0
0.5896 1.0000 0
1.0677 1.5696 1.0000];
>> b=[0.05;1.03;-0.53];U1=inv(U); L1=inv(L); X=U1*L1*b,x=A\b

X =

-111.8440
110.9531
25.7324

x =

1.4531
-1.5892
-0.2749

例2.1： 用高斯消元法求解下面的非齐次线性方程组。

解
1. 在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，所以此方程组有唯一解.
RA = 4，RB = 4，n = 4，X =[0 -0.5 0.5 0]’

例2.2 用列主元素消元法解线性方程组。
.
解
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.3 将矩阵直接进行LU分解，然后解矩阵方程。
，
解
1. 在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
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