方法一：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

 

 

使用范例：

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);
x=gs(A,b,x0,e);
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);
x=gs(A,b,x0,e);
disp('Guass-Seidel method:x=');disp(x);