In [1]:

import numpy as np
import numpy.linalg as la

In [2]:

n = 4 # crank up--8 is a good example

nums = np.arange(n).astype(np.float64)
nums

Out[2]:

array([ 0.,  1.,  2.,  3.])

In [3]:

B = 1 + nums[:, np.newaxis] + nums
B

Out[3]:

array([[ 1.,  2.,  3.,  4.],
       [ 2.,  3.,  4.,  5.],
       [ 3.,  4.,  5.,  6.],
       [ 4.,  5.,  6.,  7.]])

In [4]:

A = 1/B

from scipy.linalg import hilbert
A2 = hilbert(n)

print "difference: %g" % la.norm(A-A2)
print "cond: %g" % np.linalg.cond(A)

difference: 0
cond: 15513.7

/usr/lib/pymodules/python2.7/numpy/oldnumeric/__init__.py:11: ModuleDeprecationWarning: The oldnumeric module will be dropped in Numpy 1.9
  warnings.warn(_msg, ModuleDeprecationWarning)

In [5]:

b = np.ones((n,))

b1 = b.copy()
b1[-1] = 2.0/3.0
x1 = la.solve(A, b1)
print "residual of x1: %g"%(la.norm(b - np.dot(A, x1)))

b2 = b.copy()
b2[-1] = 0.667
x2 = la.solve(A, b2)
print "residual of x2: %g"%(la.norm(b - np.dot(A, x2)))

print "rel err in RHS: %g"%(la.norm(b2-b1)/la.norm(b1))
print "rel err in solution: %g"%(la.norm(x2-x1)/la.norm(x1))

residual of x1: 0.333333
residual of x2: 0.333
rel err in RHS: 0.000179605
rel err in solution: 0.0011524

In [6]:

x = np.ones((n,))
b = np.dot(A, x)
x3 = la.solve(A,b)
print "residual of x3: %g"%(la.norm(b - np.dot(A, x3)))
print "rel err in solution: %g"%(la.norm(x3-x)/la.norm(x))

residual of x3: 1.11022e-16
rel err in solution: 3.86179e-13

In []: