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

方法一：Gauss列主元消去法

function [x]=gauss(A,b)

    n=size(A,1);x=zeros(n,1);

    for k=1:n-1

%looking for column max anSEO靠我d 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));

eSEO靠我nd

        end

        if Max==0

            disp(This matrix can not be solved by Gauss algorithm);

        end

Atemp=A(MaxIndex,:);btemp=b(SEO靠我MaxIndex);

        A(MaxIndex,:)=A(k,:);b(MaxIndex)=b(k);

        A(k,:)=Atemp;b(k)=btemp;     

        %diagonalization

        for i=k+1:n

MSEO靠我ik=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)/ASEO靠我(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,bSEO靠我,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 methoSEO靠我d can not converge;

    else

        x=B*x0+f;

        while norm(x-x0)>=e 

            x0=x;

            x=B*x0+f;

        end

    end

end

方法三：Gauss-Seidel迭代

functionSEO靠我 [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)))>=SEO靠我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

附：不使用SEO靠我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);fSEO靠我=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 kSEO靠我

            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,SEO靠我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)

                break;

            end

x0=xSEO靠我;k=k+1;

        end

        if k==N

            disp(迭代次数已到达上限!);

        end

        disp([迭代次数 k=,num2str(k)])

    end

end

function [x]=gs_nodlu(A,b,x0,e,SEO靠我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-SeidelSEO靠我 method can not converge;

    else

        n=size(A,1);x=zeros(n,1);k=0;

        while k

            for i=1:n

                if i==1

x(1)=(b(1)-A(1,2:n)SEO靠我*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-1SEO靠我)-A(i,i+1:n)*x0(i+1:n))/A(i,i);

                end  

            end

            if norm(x-x0)

                break;

            end

            x0=x;k=k+1;

        end

        if k==N

            disp(迭代次数已到达上限!);

endSEO靠我

        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];SEO靠我

disp(answer for Q1);

%guass method

x=gauss(A,b);

disp(Guass method:x=);disp(x);

%jacobi method & gs methSEO靠我od

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,SEO靠我N);

disp(Jacobi method:x=);disp(x);

x=gs(A,b,x0,e);

disp(Guass-Seidel method:x=);disp(x);

N=100;x=gs_nodSEO靠我lu(A,b,x0,e,N);

disp(Guass-Seidel method:x=);disp(x);

%set the problem

A=3*eye(100)-circshift(eye(100),SEO靠我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 QSEO靠我2);

%guass method

x=gauss(A,b);

disp(Guass method:x=);disp(x);

%jacobi method & gs method

x0=zeros(100,1)SEO靠我;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(JacoSEO靠我bi 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,SEO靠我N);

disp(Guass-Seidel method:x=);disp(x);