用matlab求方程组解的三种方法

方法一：Gauss列主元消去法

function [x]=gauss(A,b)
    n=size(A,1);x=zeros(n,1);
    for k=1:n-1
        %looking for column max and exchange rows
        Max=abs(A(k,k));MaxIndex=k;
        for u=k+1:n
            if(abs(A(u,k))>Max)
                MaxIndex=u;Max=abs(A(u,k));
            end
        end
        if Max==0
            disp('This matrix can not be solved by Gauss algorithm');
        end
        Atemp=A(MaxIndex,:);btemp=b(MaxIndex);
        A(MaxIndex,:)=A(k,:);b(MaxIndex)=b(k);
        A(k,:)=Atemp;b(k)=btemp;     
        %diagonalization
        for i=k+1:n
            Mik=A(i,k)/A(k,k);b(i)=b(i)-Mik*b(k);
            for j=k+1:n
                A(i,j)=A(i,j)-Mik*A(k,j);
            end   
        end
    end
    %solving
    x(n)=b(n)/A(n,n);
    for i=n-1:-1:1
        x(i)=(b(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);
    end
end

 

 

方法二：Jacobi迭代

function [x]=jacobi(A,b,x0,e)
    D=diag(diag(A));L=tril(A,-1);U=triu(A,1);
    B=-D\(L+U);f=D\b;
    if max(abs(eig(B)))>=1
        x='Jacobi method can not converge';
    else
        x=B*x0+f;
        while norm(x-x0)>=e 
            x0=x;
            x=B*x0+f;
        end
    end
end

 

 

方法三：Gauss-Seidel迭代

function [x]=gs(A,b,x0,e)
    D=diag(diag(A));L=tril(A,-1);U=triu(A,1);
    B=-(D+L)\U;f=(D+L)\b;
    if max(abs(eig(B)))>=1
        x='Guass-Seidel method can not converge';
    else
        x=B*x0+f;
        while norm(x-x0)>=e 
            x0=x;
            x=B*x0+f;
        end
    end
end

 

 

附：不使用DLU分解的迭代算法

function [x]=jacobi_nodlu(A,b,x0,e,N)
    D=diag(diag(A));L=tril(A,-1);U=triu(A,1);
    B=-D\(L+U);f=D\b;
    if max(abs(eig(B)))>=1
        x='Jacobi method can not converge';
    else
        n=size(A,1);x=zeros(n,1);k=0;
        while k<N
            for i=1:n
                if i==1
                    x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);
                elseif i==n
                    x(n)=(b(n)-A(n,1:n-1)*x0(1:n-1))/A(n,n);
                else
                    x(i)=(b(i)-A(i,1:i-1)*x0(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);       
                end  
            end
            if norm(x-x0)<e
                break;
            end
            x0=x;k=k+1;
        end
        if k==N
            disp('迭代次数已到达上限!');
        end
        disp(['迭代次数 k=',num2str(k)])
    end
end

 

 

function [x]=gs_nodlu(A,b,x0,e,N)
        D=diag(diag(A));L=tril(A,-1);U=triu(A,1);
    B=-(D+L)\U;f=(D+L)\b;
    if max(abs(eig(B)))>=1
        x='Guass-Seidel method can not converge';
    else
        n=size(A,1);x=zeros(n,1);k=0;
        while k<N
            for i=1:n
                if i==1
                    x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);
                elseif i==n
                    x(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);
                else
                    x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);       
                end  
            end
            if norm(x-x0)<e
                break;
            end
            x0=x;k=k+1;
        end
        if k==N
            disp('迭代次数已到达上限!');
        end
        disp(['迭代次数 k=',num2str(k)])
    end
end

 

 

使用范例：

clear
clc
%set the problem
A=[2,-2,-1;
   4,1,-2;
   -2,1,-1];
b=[-2;1;-3];
disp('answer for Q1');
%guass method
x=gauss(A,b);
disp('Guass method:x=');disp(x);
%jacobi method & gs method
x0=[0;0;0];e=1e-3;
x=jacobi(A,b,x0,e);
disp('Jacobi method:x=');disp(x);
N=100;x=jacobi_nodlu(A,b,x0,e,N);
disp('Jacobi method:x=');disp(x);
x=gs(A,b,x0,e);
disp('Guass-Seidel method:x=');disp(x);
N=100;x=gs_nodlu(A,b,x0,e,N);
disp('Guass-Seidel method:x=');disp(x);

%set the problem
A=3*eye(100)-circshift(eye(100),1)-circshift(eye(100),-1);A(1,100)=0;A(100,1)=0;
b([1:100],1)=1;b(1,1)=2;b(100,1)=2;
disp('answer for Q2');
%guass method
x=gauss(A,b);
disp('Guass method:x=');disp(x);
%jacobi method & gs method
x0=zeros(100,1);e=1e-4;
x=jacobi(A,b,x0,e);
disp('Jacobi method:x=');disp(x);
N=100;x=jacobi_nodlu(A,b,x0,e,N);
disp('Jacobi method:x=');disp(x);
x=gs(A,b,x0,e);
disp('Guass-Seidel method:x=');disp(x);
N=100;x=gs_nodlu(A,b,x0,e,N);
disp('Guass-Seidel method:x=');disp(x);