Matplotlib - Error in plotting 2 values with same dimensions

Question

def Hamiltonian(alpha,h):

    Sx = np.array([[0,1],[1,0]])
    Sy = np.array([[0,-1j],[1j,0]])
    Sz = np.array([[1,0],[0,-1]])
    I  = np.array([[1,0],[0,1]])

    H = -1*((alpha*np.kron(np.kron(Sx,Sx),I))
       + (alpha*np.kron(np.kron(Sy,Sy),I))
       + (alpha*np.kron(np.kron(I,Sx),Sx))
       + (alpha*np.kron(np.kron(I,Sy),Sy))
       + (h*np.kron(np.kron(I,Sz),I)))

    return H
np.set_printoptions(linewidth=100)
Hamiltonian(1,0.5).real

Which returns the following matrix (just inputting this for clarity in what I'm trying to do below)

After defining Hamiltonian I want to look for entropy as a function of the h parameter. The physics behind the code does not matter for this type of question.

# von Neumann entropy as a function of h and beta - Complete
# Definition of a mixed state: [Thermal Density Matrix used]
h = np.arange(0,2.5,0.1)
beta = 2
for i in range(h.size):
    H = Hamiltonian(1.0, h[i] )
    rho_thermal = expm(-1.0 * beta * H)
    tr = np.trace(rho_thermal)
    rho_thermal = rho_thermal / tr
    np.set_printoptions(linewidth=100)
    eigvals_rho_thermal, eigvecs_rho_thermal = LA.eigh(rho_thermal)
    # Entropy
    s = 0.0
    for i in range(eigvals_rho_thermal.size):
        s += -1.0 * (eigvals_rho_thermal[i] * np.log(eigvals_rho_thermal[i]))
    print(s)
plt.plot(h,s)

My question is why do i get the following error:

My code returns 25 values for s and 25 for h why will it not plot them?

Answer 1

You need to store the value of s you obtain in each loopstep somewhere to be able to plot each s for each h later on. In general you may simplify your code a bit and use a function instead of a for-loop. Here is how I would do it.

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm

def Hamiltonian(alpha,h):

    Sx = np.array([[0,1],[1,0]])
    Sy = np.array([[0,-1j],[1j,0]])
    Sz = np.array([[1,0],[0,-1]])
    I  = np.array([[1,0],[0,1]])

    H = -1*((alpha*np.kron(np.kron(Sx,Sx),I))
       + (alpha*np.kron(np.kron(Sy,Sy),I))
       + (alpha*np.kron(np.kron(I,Sx),Sx))
       + (alpha*np.kron(np.kron(I,Sy),Sy))
       + (h*np.kron(np.kron(I,Sz),I)))

    return H

# von Neumann entropy as a function of h and beta - Complete
# Definition of a mixed state: [Thermal Density Matrix used]
def get_entropy(beta, h, alpha=1.0):
    H = Hamiltonian(alpha, h)
    rho_thermal = expm(-1.0 * beta * H)
    tr = np.trace(rho_thermal)
    rho_thermal = rho_thermal / tr
    eigvals_rho_thermal, eigvecs_rho_thermal = np.linalg.eigh(rho_thermal)
    # Entropy
    s = -np.sum(eigvals_rho_thermal*np.log(eigvals_rho_thermal))
    return s

h = np.arange(0,2.5,0.1)
beta = 2    
s = [get_entropy(beta, hi) for hi in h]
plt.plot(h,s)

plt.show()