How to programmatically calculate real eigenvalues and (optionally) complex eigenvectors for big matrix

Question

I am trying to solve 1D timeindependent schrodinger equation $$ -\frac{\hbar^{2}}{2 m} \frac{\partial^{2} \psi}{\partial x^{2}}+V(x) \psi=E \psi $$ for periodic potential to simulate crystal lattice using computational mathematics methods: $$\dfrac{\partial^2\psi}{\partial x^2}=\dfrac{\psi_{i+1}-2\psi_i+\psi_{i-1}}{\Delta x^2}$$ As I understand, I have to find eigenvalues for this matrix $$ \left[\begin{array}{cccc} \frac{1}{\Delta y^{2}}+L^{2} m V_{1} & -\frac{1}{2 \Delta y^{2}} & 0 & 0 \ldots \\ -\frac{1}{2 \Delta y^{2}} & \frac{1}{\Delta y^{2}}+L^{2} m V_{2} & -\frac{1}{2 \Delta y^{2}} & 0 \ldots \\ \ldots & \ldots & \ldots & -\frac{1}{2 \Delta y^{2}} \\ \ldots 0 & 0 & -\frac{1}{2 \Delta y^{2}} & \frac{1}{\Delta y^{2}}+L^{2} m V_{N-1} \end{array}\right]\left[\begin{array}{cc} \psi_1\\ \psi_2\\ \ldots\\ \psi_{N-1}\end{array}\right]=L^2mE\left[\begin{array}{cc} \psi_1\\ \psi_2\\ \ldots\\ \psi_{N-1}\end{array}\right];\quad L^2=2/\hbar^2 $$ Perturbation theory for weak connection gives continuous spectrum, so probably I have to take a really big matrix to get characteristic polynomial of huge order. I am looking for methods or at least python library to solve such a big equation. If it is important, I have nearly 5 GB of RAM and AMD® Athlon gold 3150u CPU, but probably I'll be able to use my university server

Answer 1

Ok, I don't need to calculate complex eigenvectors because all eigenstates are independent and are real in some point in time. So with numpy for potential $$U=\dfrac{U_0}{10(cos(i-5)\pi/5)^2+1};\quad U_0\approx 1\;eV$$ where $U_0$ I took nearly $1,6$ eV because experimentally silicium has 1,21 eV conduction band and where $i-$number if node (10 nodes per atom and 100 atoms). I used this code to calculate eigenvectors

import numpy as np
import math
import scipy.linalg as la
from numpy.linalg import eigvalsh
def potential(i):
    return -1/(10*(math.cos((i-5)*math.pi/5)-1)*(math.cos((i-5)*math.pi/5)-1)+1)
def element(i,j):
    if i==j:
        return 0.1524+potential(i)
    if abs(i-j)==1:
        return -0.1524/2
    return 0
I = 1j
a=[]
for i in range(0,5000):
    b=[]
    for j in range(0,5000):
        b.append(element(i,j))
    a.append(b)
arr = np.array(a)
print('matrix created')
eigvals, eigvecs = la.eig(arr)
print('eigenvalues and eigenvectors calculated')
eigvals = eigvals.real
f = open('Energies.txt', 'w')
for element in eigvals:
    f.write(str(element)+'\n')
print("eigenvalues stored")
f = open('Psi_functions.txt', 'w')
for element in eigvecs:
    f.write(str(element)+'\n') 

If anybody need Psi functions for selfpurposes you can get them here.

It's much simplier to calculate eigenvalues so I made 200 nodes per atom: $$U=\dfrac{U_0}{10(cos(i/20-5)\pi/5)^2+1};\quad U_0\approx 1\;eV$$ and used eigvalsh(arr) instead off la.eig(arr). So I got these energies: