下边代码段是关于N lowest eigenvalues of the tridiagonal matrix in python的代码,希望能对码农们有较大用。
''' r = lamRange(d,c,N).
Returns the sequence {r[0],r[1],...,r[N]} that
separates the N lowest eigenvalues of the tridiagonal
matrix [A] = [cdc]; that is, r[i] < lam[i] < r[i+1].
'''
from numpy import ones
def lamRange(d,c,N):
lamMin,lamMax = gerschgorin(d,c)
r = ones(N+1)
r[0] = lamMin
# Search for eigenvalues in descending order
for k in range(N,0,-1):
# First bisection of interval(lamMin,lamMax)
lam = (lamMax + lamMin)/2.0
h = (lamMax - lamMin)/2.0
for i in range(1000):
# Find number of eigenvalues less than lam
p = sturmSeq(d,c,lam)
numLam = numLambdas(p)
# Bisect again & find the half containing lam
h = h/2.0
if numLam < k: lam = lam + h
elif numLam > k: lam = lam - h
else: break
# If eigenvalue located, change the upper limit
# of search and record it in [r]
lamMax = lam
r[k] = lam
return r