用matlab实现高斯列主元消去法解线性方程及LU分解

 

function x=gaussLinearEquation(A,b)

 

%高斯法解线性方程Ax=b

 

disp('原方程为AX=b：')

 

A

 

b

 

disp('------------------------')

 

n=length(b);

 

eps=10^-2;

 

for k=1:n-1

 

    %找列主元

 

    [mainElement,index]=max(abs(A(k:n,k)));

 

    index=index+k-1;%index在A(k:n,k)中的行号转换为在A中的行号

 

    if abs(mainElement)<eps

 

        disp('列元素太小！！');

 

        break;

 

    elseif index>k

 

        %列主元所在行不是当前行，将当前行与列主元所在行交换

 

        temp=A(k,:);

 

        A(k,:)=A(index,:);

 

        A(index,:)=temp;

 

    end

 

    %消元

 

    for i=k+1:n

 

        m(i,k)=A(i,k)/A(k,k);%A(k,k)将A(i,k)消为0所乘系数

 

        A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n);%第i行消元处理

 

        b(i)=b(i)-m(i,k)*b(k);%还有b也需要处理！！       

 

    end

 

end

 

disp('消元后所得到的上三角阵是')

 

A

 

%回代

 

b(n)=b(n)/A(n,n);

 

for i=n-1:-1:1

 

    %sum(A(i,i+1:n).*b(i+1:n)')表示已知

 

    b(i)=(b(i)-sum(A(i,i+1:n).*b(i+1:n)'))/A(i,i);

 

end

 

clear x;

 

x=b;

 

disp('AX=b的解x是')

 

x

 

用法：

 

在控制台输入：

 

A=[1.003 0.333 1.504 -0.333;

 

   -2.011 1.455 0.506 2.956;

 

   4.329 -1.952 0.006 2.087;

 

   5.113 -4.004 3.332 -1.112];

 

b=[ 3.005,5.407,0.136,3.772 ]';

 

执行gaussLinearEquation(A,b);即可得到结果。

使用matlab实现矩阵的LU分解

function [L,U]=myLU(A)

%实现对矩阵A的LU分解，L为下三角矩阵

A

[n,n]=size(A);

L=zeros(n,n);

U=zeros(n,n);

for i=1:n

    L(i,i)=1;

end

for k=1:n

    for j=k:n

        U(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j)');

       

    end

    for i=k+1:n

        L(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)'))/U(k,k);       

    end

   

end

用法，在控制台输入

A=[1 2 3 -4;-3 -4 -12 13;2 10 0 -3;4 14 9 -13];

然后执行[l,u]=myLU(A);

这样得到l和u，可以通过l*u与A比较来验证LU分解的正确性。

 

注意：

[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就行了。