Estou desenvolvendo um código maior para realizar cálculos de autovalor de enormes matrizes esparsas, no contexto da física computacional. Testo minhas rotinas contra o oscilador harmônico simples em uma dimensão, uma vez que os autovalores são bem conhecidos analiticamente. Fazendo isso e comparando minhas próprias rotinas aos solucionadores embutidos do SciPy, deparei-me com a estranheza exibida no gráfico abaixo. Aqui você pode ver os 100 primeiros autovalores computados numericamente e autovalores analíticos λ a n a
Em torno do valor próprio 40, os resultados numéricos começam a divergir dos analíticos. Isso não me surpreende (não vou explicar o porquê aqui, a menos que apareça na discussão). No entanto, o que é surpreendente para mim é que eigsh () produz autovalores degenerados (em torno do autovalor número 80). Por que o eigsh () se comporta assim mesmo com um número tão pequeno de valores próprios?
import numpy as np
from scipy.sparse.linalg import eigsh
import myFunctions as myFunc
import matplotlib.pyplot as plt
#discretize x-axis
N = 100
xmin = -10.
xmax = 10.
accuracy = 1e-5
#stepsize
h = (xmax - xmin) / (N + 1.)
#exclude first and last points since we force wave function to be zero there
x = np.linspace(-10. + h,10. - h,N)
#create potential
V = x**2
def fivePoint(N,h,V):
C0 = (np.ones(N))*30. / (12. * h * h) + V
C1 = (np.ones(N)) * (-16.) / (12. * h * h)
C2 = (np.ones(N)) / (12. * h * h)
H = sp.spdiags([C2, C1, C0, C1, C2],[-2, -1, 0, 1, 2],N,N)
return H
H = myFunc.fivePoint(N,h,V)
eigval,eigvec = eigsh(H, k=N-1, which='SM', tol=accuracy)
#comparison analytical and numerical eigenvalues
xAxes = np.linspace(0,len(eigval)-1,len(eigval))
analyticalEigval = 2. * (xAxes + 0.5)
plt.figure()
plt.plot(xAxes,eigval, '+', label=r"$\lambda_{num}$")
plt.plot(xAxes,analyticalEigval, label=r"$\lambda_{ana}$")
plt.xlabel("Number of Eigenvalue")
plt.ylabel("Eigenvalue")
plt.legend(loc=4)
plt.title("eigsh()-method: Comparison of $\lambda_{num}$ and $\lambda_{ana}$")
plt.show()
Respostas:
A degeneração de alguns autovalores me parece a marca registrada da quebra do algoritmo de Lanczos . O algoritmo de Lanczos é um dos métodos mais comumente usados para aproximar os autovalores e autovetores das matrizes hermitianas; é o que scipy.eigsh () usa, através de uma chamada para a biblioteca ARPACK .
Na aritmética exata, o algoritmo de Lanczos produz um conjunto de vetores ortogonais, mas na aritmética de ponto flutuante eles podem falhar em ser ortogonais e até se tornar linearmente dependentes. O fato realmente irritante é que essa perda de ortogonalidade ocorre precisamente quando um dos valores próprios aproximados converge para um dos valores reais reais - o algoritmo sabota a si mesmo, por assim dizer. O resultado é que você obterá pares espúrios de valores próprios próximos. Existem várias correções para isso, por exemplo, usando Gram-Schmidt para forçar qualquer vetor próprio convergente a ser ortogonal a cada passo.
No entanto, nenhum método é perfeito, especialmente se você estiver tentando calcular todo o espectro da sua matriz . Portanto, se você estiver tentando obter os 50 autovalores menores, é melhor aproximar a função de onda por um vetor com 100 elementos e solicitar apenas
eigsh()os primeiros 50 níveis de energia, em vez de usar um vetor com 50 pontos e pedir tudo dos autovalores.Se você quiser ler mais, consulte os Métodos numéricos de Yousef Saad para grandes problemas de autovalor .
fonte