matlab求方程組的解函式：

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

2.方法二：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);