Search code examples
pythonscipyeigenvaluearpack

Strange behaviour of ARPACK for hermitian matrix

I want to obtain numerically the energy of the ground state of some hermitian matrix (see the definition of this matrix in the following code) and plot it in terms of the matrix-parametter "phase". 

import scipy.sparse as sparse
import scipy
import numpy
import numpy as np
import math
from scipy.special import binom
import cmath
import sympy
import matplotlib.pyplot as plt
import pylab
from copy import *
from numpy import linalg as LA



M=5#DIMENSION OF THE MATRIX


def tunneling(phase):#HAMILTONIAN MATRIX
    Matrix_hop = [[0 for x in range(M)] for y in range(M)]
    for i in range(M):
            if i+1==M:
                Matrix_hop[i][0] = -1.0
                Matrix_hop[i][i-1] = -1.0
        else:
                Matrix_hop[i][i+1] = -1.0
                Matrix_hop[i][i-1] = -1.0
    Matrix_hop[0][M-1]=-1.0*cmath.exp(1j*phase)
    Matrix_hop[M-1][0]=-1.0*cmath.exp(-1j*phase)
    return Matrix_hop 

def eigen_system(H):
    values, vectors = sparse.linalg.eigs(H,2,which='SR') #ARPACK!!
    energy_ground = values[0]
    return  vectors[:,0], energy_ground


init = 0.0
points = 1000  
final_value = 2*math.pi
steep = (final_value-init)/points
list_values_phase = np.arange(init,final_value,steep)
f1 = open("ground_state_energy.dat", "w")
for i in list_values_phase:
    phase = i
    f1.write(str(phase)+" ")
    H = np.asarray(tunneling(i))
    f1.write(str(np.real(eigen_system(H)[1]))+" ")
    f1.write("\n")
f1.close()



datalist = pylab.loadtxt("ground_state_energy.dat")
pylab.plot( datalist[:,0], datalist[:,1],label="ground state" )
pylab.legend()
pylab.xlabel("phase")
pylab.ylabel("Energy")
pylab.show()

I have used the ARPACK in Python for hermitian matrices, which is done using sparse.linalg.eigs. The problem is that, as you can see in the following figure, the ground state energy is not correctly computed, there is a lot of peaks, which means that the ground state is not correctly found. Actually seems that for this peaks, ARPACK doens't find the ground state and it obtain the first excitated state. enter image description here It's a very strange problem, since this matrix that I am using (which comes from quantum mechanics) can be solved analitycally as well as using Mathematica, and with ARPACK in Python doesn't work. Someone has some idea why this happens and how can be solved?? Thank you

I am using the last versión of scipy 0.19.1

Solution

  • In this function

    def eigen_system(H):
    values, vectors = sparse.linalg.eigs(H,2,which='SR') #ARPACK!!
    energy_ground = values[0]
    return  vectors[:,0], energy_ground

    you find the first two eigenvalues, and then take the first. The function eigs does not guarantee that the eigenvalues that it finds are ordered, and sometimes the first one is not the smallest.

    Instead of finding the two smallest, why not just find the smallest?

        values, vectors = sparse.linalg.eigs(H, 1, which='SR')  # ARPACK!!

    When I make that change, I get this plot:

    plot