Matlab线性方程组求解

线性方程组求解

1. 直接法

Gauss 消元法：

function x=DelGauss(a,b)

% Gauss 消去法

[n,m]=size(a);

nb=length(b);

det=1;% 存储行列式值

x=zeros(n,1);

for k=1:n-1

    for i=k+1:n

        if a(k,k)==0

            return

        end

     m=a(i,k)/a(k,k);

     for j=k+1:n

         a(i,j)=a(i,j)-m*a(k,j);

     end

     b(i)=b(i)-m*b(k);

  end

  det=det*a(k,k);

end

det=det*a(n,n);

   

 

for k=n:-1:1  % 回代

    for j=k+1:n

        b(k)=b(k)-a(k,j)*x(j);

    end

    x(k)=b(k)/a(k,k);

end

Example:

>> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898];

>> b=[1 0 1]';

>> x=DelGauss(A,b)

 

x =

 

    0.9739

   -0.0047

    1.0010

列主元 Gauss 消去法：

function x=detGauss(a,b)

% Gauss 列主元消去法

[n,m]=size(a);

nb=length(b);

det=1;% 存储行列式值

x=zeros(n,1);

for k=1:n-1

    amax=0;% 选主元

    for i=k:n

        if abs(a(i,k))>amax

            amax=abs(a(i,k));r=i;

        end

    end

    if amax<1e-10

       return;

    end

   

   

    if r>k  % 交换两行

        for j=k:n

            z=a(k,j);a(k,j)=a(r,j);a(r,j)=z;

        end

        z=b(k);b(k)=b(r);b(r)=z;det=-det;

     end

   

    for i=k+1:n   % 进行消元

        m=a(i,k)/a(k,k);

        for j=k+1:n

            a(i,j)=a(i,j)-m*a(k,j);

        end

        b(i)=b(i)-m*b(k);

    end

    det=det*a(k,k);

end

det=det*a(n,n);

 

for k=n:-1:1  % 回代

    for j=k+1:n

        b(k)=b(k)-a(k,j)*x(j);

    end

    x(k)=b(k)/a(k,k);

end

Example:

>> x=detGauss(A,b)

 

x =

 

    0.9739

   -0.0047

1.0010

 

Gauss-Jordan 消去法 ：

function x=GaussJacobi(a,b)

% Gauss-Jacobi 消去法

[n,m]=size(a);

nb=length(b);

x=zeros(n,1);

for k=1:n

    amax=0;% 选主元

    for i=k:n

        if abs(a(i,k))>amax

            amax=abs(a(i,k));r=i;

        end

    end

    if amax<1e-10

       return;

    end

   

   

    if r>k  % 交换两行

        for j=k:n

            z=a(k,j);a(k,j)=a(r,j);a(r,j)=z;

        end

        z=b(k);b(k)=b(r);b(r)=z;

    end

     % 进行消元

     b(k)=b(k)/a(k,k);

        for j=k+1:n

            a(k,j)=a(k,j)/a(k,k);

        end

        for i=1:n

            if i~=k

                for j=k+1:n

                    a(i,j)=a(i,j)-a(i,k)*a(k,j);

                end

                 b(i)=b(i)-a(i,k)*b(k);

            end

        end

    end

   

    for i=1:n

        x(i)=b(i);

    end

Example:

>> x=GaussJacobi(A,b)

 

x =

 

    0.9739

   -0.0047

    1.0010

LU 分解法：

function [l,u]=lu(a)

%LU 分解

n=length(a);

l=eye(n);

u=zeros(n);

for i=1:n

    u(1,i)=a(1,i);

end

for i=2:n

    l(i,1)=a(i,1)/u(1,1);

end

 

for r=2:n

    %%%%

    for i=r:n

        uu=0;

        for k=1:r-1

            uu=uu+l(r,k)*u(k,i);

        end

        u(r,i)=a(r,i)-uu;

    end

    %%%%

    for i=r+1:n

        ll=0;

        for k=1:r-1

                ll=ll+l(i,k)*u(k,r);

            end

            l(i,r)=(a(i,r)-ll)/u(r,r);

        end

     %%%%

End

function x=lusolv(a,b)

%LU 分解求解线性方程组 aX=b

if length(a)~=length(b)

    error('Error in inputing!')

    return;

end

n=length(a);

[l,u]=lu(a);

<