﻿ 130213003戴斌作业7 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 130213003戴斌作业7

1.估计5到20阶Hlibert矩阵的无穷条件数

k=zeros(16,1);
for n=5:1:20
for i=1:n
for j=1:n
A(i,j)=1.0/(i+j-1);
end
end
k((n-4),1)=cond(A,inf)
end

k =

1.0e+018 *

0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0012
0.0371
0.3650
4.0846
0.9756
1.0286
1.0381
9.2880
6.2804
8.6548

2.

function x=gausslu_pp(A,b)
n=length(b);
max=A(1,1);
p=1;
for k=1:n-1
for i=k:n
if A(i,k)>max
max=A(i,k);
p=i;
end
end
A([p,k],:)=A([k,p],:);

z=b(p);
b(p)=b(k);
b(k)=z;
end

tol=10^(-16);
for k=1:n-1
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);
end
end
L=eye(n)+tril(A,-1);
U=triu(A,0);

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=zeros(n,1);
x(n)=b(n)/U(n,n);
for j=n-1:-1:1
x(j)=(b(j)-U(j,j+1:n)*x(j+1:n))/U(j,j);
end

function [A,b,X]=improve(n)
A=zeros(n,n);
A=tril(-ones(n,n),-1);
for i=1:n
A(i,i)=1;
end
A(:,n)=1;
X=rand(n,1);
b=A*X;

for n=5:1:20;
[A,b,X]=improve(n);
x=gausslu_pp(A,b);
while ((norm((X-x),inf)/norm(X,inf))>eps
r=b-A*x;
z=gausslu_pp(A,r);
X=x+z;
x=X;
end
x
end

x =

0.1835
0.3685
0.6256
0.7802
0.0811
0.9294
0.7757
0.4868
0.4359
0.4468
0.3063
0.5085
0.5108
0.8176
0.7948
0.6443
0.3786
0.8116
0.5328
0.3507