病态矩阵的例子

从一个数值结果谈起 从这里至少可以看出三点，一理论上的满秩不能保证实际计算的满秩，有时甚至相差很大；二数学软件只能解决常态的问题，难以解决大量存在的有特殊性态的问题；三要研究非方阵、非满秩的问题。 * 以为基底，构造不超过n 次最佳平方逼近多项式时，其法方程的系数矩阵会出现所谓Hilbert矩阵。

n=5

>> A=hilb(5)

A =

1.0000 0.5000 0.3333 0.2500 0.2000

0.5000 0.3333 0.2500 0.2000 0.1667

0.3333 0.2500 0.2000 0.1667 0.1429

0.2500 0.2000 0.1667 0.1429 0.1250

0.2000 0.1667 0.1429 0.1250 0.1111

>> Eigenvalue=eig(A)

Eigenvalue =

0.0000

0.0003

0.0114

0.2085

1.5671

>> cond(A)

ans =

4.7661e+005

>> b=A*ones(5,1);

>> x=inv(A)*b;

>> x

solution =x

1.0000

1.0000

1.0000

1.0000

1.0000

>> norm(ones(5,1)-x)

ans =

3.0055e-011

>> norm(b-A*x)

ans =

6.9024e-012

n=8

>> A=hilb(8);

>> Eigenvalue=eig(A)'

Eigenvalue =

0.0000 0.0000 0.0000 0.0001 0.0015 0.0262 0.2981 1.6959

>> cond(A)

ans =

1.5258e+010

>> b=A*ones(8,1);

>> x=(inv(A)*b)'

x =

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

>> norm(x'-ones(8,1))

ans =

9.3866e-007

>> norm(A*x'-b)

ans =

7.0278e-008

就是说，对于系数矩阵为8阶的Hilbert矩阵，用Matlab中的”inv”语句，大致能达到6-7位有效数字，这时条件数为10的10次方量级。

n=10;

>> A=hilb(10);

>> cond(A)

ans =

1.6025e+013

> b=A*ones(10,1);

>> solution=inv(A)*b;

>> solution'

ans =

1.0000 1.0000 1.0000 1.0000 0.9999 1.0003 0.9991 1.0007 0.9996 1.0001

>> norm(ones(10,1)-solution)

ans =

0.0012

>> norm(A*solution-b)

ans =

7.6071e-005

n=12;

>> A=hilb(12);

>> cond(A)

ans =

1.7945e+016

>> x=inv(A)*(A*ones(12,1));

Warning: Matrix is close to singular or badly scaled.

Results may be inaccurate. RCOND = 2.632091e-017.

>> rank(A)

ans =

11

>> A=hilb(100);

>> rank(A)

ans =

18

>> A=hilb(1000);

>> rank(A)

ans =

24

>> rank(hilb(2000))

ans =

25

将系数矩阵作微调为Hilbert(2000)+7*10(-7)*

>> A=hilb(2000)+0eye(2000);

>> rank(A)

ans =

2000

广义逆矩阵

的理论与方法

南京大学计算机系

赵金熙

2009.10.15

各类数学软件的广泛使用，使得大量以前不能计算或者计算非常困难的问题成为可能或容易计算的问题。也就是说，现代数值方法丰硕成果，譬如上世纪70年代以来，样