高斯消元法求解线性方程,包括把增广矩阵转换为三角矩阵形式的过程,消去阶段工作
步骤是把矩阵A分解成为下三角L和上三角U的乘积。这种计算L,U的过程称为LU分解法。
lu实现对矩阵的分解。
[L,U] = lu(A)
%%将矩阵A分解的上三角矩阵保存在U当中,将一个“心理学上的”下三角矩阵(例如一个下三角矩阵和置换矩阵的乘积)保存在L中,满足A=L*U,注意A不必须是方阵。
[L,U,P] = lu(A)
%%返回三个矩阵,下三角矩阵L、上三角矩阵U和一个置换矩阵P,满足P*A=L*U。
[L,U,p] = lu(A,'vector')
%%将使用向量而不是矩阵的方式返回置换信息,p是一个行向量,满足A(p,:) = L*U,如果使用[L,U,P] = lu(A,'matrix')则返回的P是向量,这也是默认的形式。
[L1,U1]=lu(x)
[L2,U2,P]=lu(x)
[L3,U3,P,Q] = lu(X)
MATLAB中[L1,U1]=lu(x)的结果:L是下三角的置换矩阵即L1=p*L2,U是上三角阵。Matlab中LU分解采用高斯消元法,结果是不唯一的,只要[L1,U1]=lu(x)满足L1*U1=x, [L2,U2,P]=lu(x)满足L2*U2=p*x,[L3,U3,P,Q] = lu(X)满足 L3*U3= P*X*Q就行了
例方程组:
x1+0.2*(x2)+0.5*(x3)=1
0.2*(x1)+0.4*(x2)+x3=2
0.5*(x1)+0.1*(x2)+0.6*(x3)=1.5
根据LU分解有LUx=b,记Y=L^(-1)b,则x=U^(-1)Y为方程组Ax=b解。
>> a=[1 0.2 0.5;0.2 0.4 1; 0.5 0.1 0.6];
>> b=[1 2 1.5];
>> [L,U]=lu(a)
L =
1.0000 0 0
0.2000 1.0000 0
0.5000 0 1.0000
U =
1.0000 0.2000 0.5000
0 0.3600 0.9000
0 0 0.3500
>> y=inv(L)*b'
y =
1.0000
1.8000
1.0000
或者
>> y=L^(-1)*b'
y =
1.0000
1.8000
1.0000
>> x=U^(-1)*y
x =
0
-2.1429
2.8571
或
>> x=inv(U)*y
x =
0
-2.1429
2.8571
下面就是lu函数的厂商代码,相信这段程序,对大家编写线性方程组直接解法的诸多算法都有启发作用。
function [L,U,x]=lux(A,b)
%LU 分解法解线性方程组(列主元LU分解)
[n,n]=size(A);
p=eye(n);%p记录了选择主元时候所进行的行变换
for k=:n-
[r,m]=max(abs(A(k:n,k))); %选列主元
m=m+k-;
if(A(m,k)~=)
if(m~=k)
A([k m],:)=A([m k],:);
p([k m])=p([m k]);
end
for i=k+:n
A(i,k)=A(i,k)/A(k,k);
j=k+:n; A(i,j)=A(i,j)-A(i,k)*A(k,j);
end
end
end
L=tril(A,-)+eye(n,n);
U=triu(A);
%解下三角矩阵 Ly=b
newb=p*b;
y=zeros(n,);
for k=:n
j=:k-;
y(k)=(newb(k)-L(k,j)*y(j))/L(k,k);
end
%解上三角方程组 Ux=y
x=zeros(n,);
for k=n:-:
j=k+:n;
x(k)=(y(k)-U(k,j)*x(j))/U(k,k);
end