﻿ 陈思130213001作业七 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 陈思130213001作业七

Gauss: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);
matrix: 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;
test: for n=5:30
[A,b,x]=matrix(n);
y=gauss(A,b);
error=norm(x-y,2)
end

5.1179e-016

error =

3.0375e-016

error =

2.3149e-016

error =

2.9710e-015

error =

5.5174e-015

error =

1.7947e-015

error =

1.3352e-014

error =

1.4152e-014

error =

5.2444e-014

error =

1.9273e-013

error =

4.5993e-013

error =

7.4337e-013

error =

1.4850e-012

error =

2.0571e-012

error =

6.4771e-012

error =

1.0776e-012

error =

7.8049e-012

error =

4.6851e-011

error =

2.4160e-011

error =

2.5830e-011

error =

1.2246e-010

error =

7.5735e-010

error =

3.2478e-009

error =

1.5314e-009

error =

1.9735e-009

error =

6.1606e-009