Saya sedang mengembangkan beberapa kode yang lebih besar untuk melakukan perhitungan nilai eigen dari matriks jarang besar, dalam konteks fisika komputasi. Saya menguji rutinitas saya terhadap osilator harmonik sederhana dalam satu dimensi, karena nilai eigen dikenal secara analitis. Melakukannya dan membandingkan rutinitas saya sendiri dengan pemecah bawaan SciPy, saya telah menemukan keanehan yang ditampilkan dalam plot di bawah ini. Di sini Anda dapat melihat 100 nilai eigen yang dihitung secara numerik pertama dan nilai eigen analitik λ a n a
Sekitar nilai eigen angka 40, hasil numersial mulai menyimpang dari yang analitis. Ini tidak mengejutkan saya (saya tidak akan membahas mengapa di sini, kecuali itu muncul dalam diskusi). Namun, yang mengejutkan bagi saya adalah bahwa eigsh () menghasilkan nilai eigen yang merosot (sekitar nilai eigen angka 80). Mengapa eigsh () berperilaku seperti itu bahkan untuk sejumlah kecil nilai eigen?
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()
Jawaban:
Bagi saya, degenerasi beberapa nilai eigen seperti ciri khas dari kegagalan algoritma Lanczos . Algoritma Lanczos adalah salah satu metode yang lebih umum digunakan untuk memperkirakan nilai eigen dan vektor eigen dari matriks Hermitian; itulah yang digunakan oleh scipy.eigsh (), melalui panggilan ke perpustakaan ARPACK .
Dalam aritmatika yang tepat, algoritma Lanczos menghasilkan satu set vektor ortogonal, tetapi dalam aritmatika titik mengambang ini bisa gagal menjadi ortogonal dan bahkan menjadi tergantung secara linear. Yang benar-benar menyebalkan adalah hilangnya ortogonalitas ini terjadi tepat ketika salah satu dari perkiraan nilai eigen telah konvergen ke salah satu nilai eigen nyata - algoritma itu menyabotase itu sendiri, sehingga bisa dikatakan. Hasilnya adalah Anda akan mendapatkan beberapa pasangan nilai eigen terdekat. Ada berbagai perbaikan untuk ini, misalnya menggunakan Gram-Schmidt untuk memaksa vektor eigen konvergen untuk menjadi ortogonal di setiap langkah.
Meskipun demikian, tidak ada metode yang sempurna, terutama jika Anda mencoba menghitung seluruh spektrum matriks Anda . Jadi jika Anda mencoba untuk mendapatkan 50 nilai eigen terkecil, Anda mungkin lebih baik mendekati fungsi gelombang dengan vektor dengan 100 elemen dan hanya meminta
eigsh()50 tingkat energi pertama, daripada menggunakan vektor dengan 50 poin dan meminta semua dari nilai eigen.Jika Anda ingin membaca lebih lanjut, lihat Metode Numerik Yousef Saad untuk Masalah Nilai Eigen Besar .
sumber