﻿ 130213013+李艳+矩阵计算+作业7 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 130213013+李艳+矩阵计算+作业7

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

error =

2.4825e-016
error =

1.1116e-015
error =

1.6120e-015
error =

3.4837e-016
error =

3.0860e-015
error =

6.1991e-015
error =

2.2847e-014
error =

3.9232e-014
error =

1.6496e-013
error =

4.3955e-014
error =

2.1440e-013
error =

4.3213e-013
error =

1.3134e-012
error =

1.9521e-012
error =

9.1956e-013
error =

1.2180e-011
error =

1.1582e-011
error =

3.3860e-011
error =

1.8725e-010
error =

6.1072e-011
error =

5.0874e-010
error =

7.8975e-010
error =

1.9822e-009
error =

2.5387e-009
error =

3.8170e-009

error =

7.9459e-009

2.gaussliexuan.m
function [L,U,q]=guassliexuan(A)
[m,n]=size(A);
p=[];
tol=1e-16;
if m==n
for k=1:n-1
[r,c]=max(abs(A(k:n,k)));
q(k)=c+k-1;
A([k,q(k)],:)=A([q(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
'A(k,k)=0!'
break;
end
end
else
'm~=n';
end
L=eye(n)+tril(A,-1);
U=triu(A);
Test2.m
clc
error1=[];
err1=[];
k=[5:30];
for n=k
[A,b,x]=matrix(n);
y=gauss(A,b);%计算解
norma=norm(inv(A),inf);
rnorm=norm(b-A*y,inf);
bnorm=norm(b,inf);
anorm=norm(A,inf);
error=(norma*anorm*rnorm)/bnorm;%计算解的精度
error1=[error1;error];
err=norm(abs(x-y),inf)/norm(x,inf);%真是相对误差
err1=[err1;err];
end
error1
err1
error=abs(error1-err1)

error1 =

1.0e-007 *

0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0003
0.0000
0.0002
0.0009
0.0009
0.0009
0.0040
0.0371
0.0856
0.0367
0.0730
0.2688

err1 =

1.0e-008 *

0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0002
0.0005
0.0001
0.0008
0.0043
0.0023
0.0020
0.0117
0.0511
0.3175
0.1402
0.1343
0.4631

error =

1.0e-007 *

0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0003
0.0000
0.0002
0.0005
0.0007
0.0007
0.0028
0.0320
0.0538
0.0226
0.0596
0.2225