(a) Write a Python script that generates a complex random Hermitian matrix of size $N x N$ and calculates all its eigenvalues $\lambda$ The matrix’ diagonal elements should be independently identically distributed real random numbers with a Gaussian distribution with mean zero and variance one. The off-diagonal elements should have real and imaginary parts that are independently identically distributed random numbers with a Gaussian distribution with mean zero and variance one.
In [1]:
import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt
from math import factorial
from itertools import combinations_with_replacement
from scipy.integrate import quad
%matplotlib inline
In [20]:
def generate_random_hermitian(N):
A = np.random.normal(size=(N, N))
H = (np.tril(A) + np.triu(A.T.conj())) / np.sqrt(2)
return H
def check_if_hermitian(H):
return np.sum(np.matrix(H).getH() == np.matrix(H)) == H.size
(b) Plot a histogram of the scaled eigenvalues $\lambda_i / \sqrt{2N}$. Scale the histogram such that the histogram’s area is normalized to one. Set the number of bins to 64 to get a reasonableresolution.
In [23]:
N = 300
x = np.linspace(-1, 1, 1000)
plt.hist(LA.eigvalsh(generate_random_hermitian(N)) / np.sqrt(2 * N), bins=64, normed=True);
(c) Extend the script from the previous two parts to generate M matrices and plot the distribution of the N × M scaled eigen values. As one can show the distribution of the scaled eigenvalues converges in the limit N → ∞ to some limiting distribution that has a rather simple mathematical expression. Can you guess the mathematical law of this limiting distribution from the numerical data? Choose the parameters N and M such that you are close enough to the limiting case but the computing time remains limited such that the simulation can be carried out in a reasonable amount of time.
We chose N = 100 and M = 1000. The distribution is the Marčenko–Pastur distribution.
In [24]:
N = 100
M = 1000
eignval = []
for i in range(M):
eignval.extend(list(LA.eigvalsh(generate_random_hermitian(N)) / np.sqrt(2 * N)))
plt.hist(eignval, bins=64, normed=True);
plt.plot(x, .63 * np.sqrt(1 - x ** 2)); # .315 is an approximation
(d) The SciPy library has two functions to calculate the eigen values of a dense matrix, one that works for arbitrary quadratic matrices, one that is limited to Hermitian matrices. Give three reasons why one should prefer the specialized function over the more general in the case of Hermitian matrices.
My reasons are:
The normalized eigen functions of the quantum harmonic oscillator are given by
$$ \psi_n(x)\equiv \left\langle x \mid n \right\rangle = {1 \over \sqrt{2^n n! l}}~ \pi^{-1/4} e^{-x^2 / (2 l^2)} H_n(x/l), $$where $H_n(x)$ denotes the Hermite polynomials, where the first few are given by
$H_n(x) = 1$,
$H_n(x) = 2x$,
$H_n(x) = 4x^2 - 2$
(a) Plot the first five eigen functions $\Psi_n(x)$ into a single diagram for ℓ = 1. Annotate your plot appropriately.
In [5]:
def quantum_eignfunctions(x, n, l=1):
h = scipy.special.hermite(n)
coef = 1 / np.sqrt(2 ** n * factorial(n) * l) * np.pi ** -.25
return coef * np.exp(-x**2 / (2 * l * l)) * h(x / l)
In [6]:
x = np.linspace(-5, 5, 1000)
for n in range(0, 5):
plt.plot(x, quantum_eignfunctions(x, n), label=r'$n$ = ' + str(n))
plt.legend()
plt.xlabel(r'$x$')
plt.ylabel(r'$\Psi_n(x)$')
plt.grid(linestyle=':')
(b) Look into the SciPy documentation how to evaluate numerically one-dimensional integrals. Check numerically the orthogonality relation
$\int_{-\infty}^{\infty} \Psi_m(x)\Psi_n(x) = \delta_{m,n}$.
In [7]:
psi_norm = lambda x, n, m: quantum_eignfunctions(x, n) * quantum_eignfunctions(x, m)
for n_m in combinations_with_replacement(range(5), 2):
n, m = n_m
y, err = (quad(psi_norm, -np.inf, np.inf, args=(n, m)))
if y < err:
y = 0.
np.random.random((4, 4))np.random.random((4, 4)) print('For n =', n, 'and m =', m, ': integral =', round(y, 2))