﻿ 130213037+左芝+矩阵计算上机作业七 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 130213037+左芝+矩阵计算上机作业七

%列选主元

function y=gauss(A,b)
[m,n]=size(A);
p=[];
P=eye(n);
tol=1e-16;
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],:);
P([k,p(k)],:)=P([p(k),k],:);
if abs(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
L='A(k,k)=0!';
break;
end
end
else
L='m~=n';
end
L=eye(n)+tril(A,-1);
U=triu(A,0);
b=P*b;
for j=1:(n-1)
b(j)=b(j)/L(j,j);
b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);
end
b(n)=b(n)/L(n,n);
y=b;
for j=n:-1:2
y(j)=y(j)/U(j,j);
y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);
end
y(1)=y(1)/U(1,1);

%矩阵

function [A,b,x]=matrix(n)

A=ones(n,n);

A1=(-2)*tril(A,-1);

A2=(-1)*triu(A,1);

A=A+A1+A2;

A(:,n)=1;

x=rand(n,1);

b=A*x;

%验证

for n=5:30

[A,b,x]=matrix(n);

y=gauss(A,b);

error=norm(x-y,2)

end

%运行结果

error =

2.2204e-016

error =

1.5701e-016

error =

3.2311e-015

error =

5.3325e-015

error =

1.3521e-015

error =

4.4717e-015

error =

7.7664e-015

error =

2.6142e-014

error =

4.3588e-014

error =

2.4976e-013

error =

2.8337e-013

error =

8.0581e-013

error =

4.1528e-013

error =

5.7491e-012

error =

1.3495e-012

error =

2.0822e-011

error =

7.5540e-012

error =

1.0360e-010

error =

4.0775e-011

error =

1.4145e-010

error =

3.6417e-010

error =

1.0592e-009

error =

1.0323e-009