﻿ 130213025+王颖洁+作业3+Gauss消去法列选全选主元 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 130213025+王颖洁+作业3+Gauss消去法列选全选主元

function x=gauss(A,b)
[m,n]=size(A);
p=[];
q=[];
tol=1e-14;
x=zeros(n,1);
if m==n
for k=1:n-1
[r,c]=find(A==max(max(abs(A(k:n,k:n)))));
p(k)=r(1);
q(k)=c(1);
A([k,p(k)],:)=A([p(k),k],:);
A(:,[k,q(k)])=A(:,[q(k),k]);
if A(k,k)>tol
A(k+1:n,k)=A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);
else
break;
end
end
L=tril(A,-1)+eye(n)
U=triu(A,0)
for k=1:length(p)
b([k,p(k)])=b([p(k),k]);
end
b(1)=b(1)/L(1,1);
for j=2:n
b(j)=(b(j)-L(j,1:j-1)*b(1:j-1))/L(j,j);
end
x(n)=b(n)/U(n,n);
for j=n-1:-1:1
x(j)=(b(j)-x(j+1:n)*U(j,j+1:n))/U(j,j);
end
for k=length(q):-1:1
x([k,q(k)])=x([q(k),k]);
end
else
'A is not a square!'
end

A=[0.001 1.00;1.00 2.00];
b=[1.00;3.00];
x=gauss(A,b)
err=A*x-b

L =
1.0000         0
0.5000    1.0000
U =
2.0000    1.0000
0   -0.4990
x =
1.0020
0.9990
err =
0
0

function x=gauss(A,b)
[m,n]=size(A);
p=[];
tol=1e-14;
x=zeros(n,1);
if m==n
for k=1:n-1
[r,c]=find(A==max(max(abs(A(k:n,k:n)))));
p(k)=r(1);
A([k,p(k)],:)=A([p(k),k],:);
if A(k,k)>tol
A(k+1:n,k)=A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);
else
break;
end
end
L=tril(A,-1)+eye(n)
U=triu(A,0)
for k=1:length(p)
b([k,p(k)])=b([p(k),k]);
end
b(1)=b(1)/L(1,1);
for j=2:n
b(j)=(b(j)-L(j,1:j-1)*b(1:j-1))/L(j,j);
end
x(n)=b(n)/U(n,n);
for j=n-1:-1:1
x(j)=(b(j)-x(j+1:n)*U(j,j+1:n))/U(j,j);
end
else
'A is not a square!'
end

A=[0.001 1.00;1.00 2.00];
b=[1.00;3.00];
x=gauss(A,b)
err=A*x-b

L =
1.0000         0
0.0010    1.0000
U =
1.0000    2.0000
0    0.9980
x =
1.0020
0.9990
err =
0

0