Kent-Andre and Dag Sverre playing with simple linear system

gistfile1.txt

from numpy import *
from scipy.linalg import eig, eigh

N = 100 
M = 40 

noise = 0.01

i = arange(M, dtype=float)
cinv = 1 + i**2
Cinv = matrix(diag(cinv))

p = noise * ones(N)
pinv = 1. / p
pinv[1*N//5:2*N//5] = 0
Pinv = matrix(diag(pinv))

Q = matrix(zeros([M, N]))
for i in range(M): 
  for j in range(N): 
    Q[i,j] = sin(i*pi*j/N)

A = Cinv + Q*Pinv*Q.transpose()

A_jacobi = matrix(diag(diagonal(A)))

e, v = eigh(A, A_jacobi)

e = sorted(e.real)

import pylab 
pylab.plot(e)
pylab.show(e)