<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 23-10-2017
<center>  <h1> Computational Quantum Dynamics  </h1> </center>
Exercise 1
3.  Eigen values of random Hermitian matrices.
(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.
Solution
Step1:  (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.
Step2:  (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.
Solution
We chose N = 100 and M = 1000. The distribution is the Marčenko–Pastur distribution.
Step3:  (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.
Solution
My reasons are
Step4:  (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}$.
<ASSISTANT_TASK:>
Python Code:
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  
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
Explanation: 23-10-2017
<center>  <h1> Computational Quantum Dynamics  </h1> </center>
Exercise 1
3.  Eigen values of random Hermitian matrices.
(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.
Solution
End of explanation
N = 300
x = np.linspace(-1, 1, 1000)
plt.hist(LA.eigvalsh(generate_random_hermitian(N)) / np.sqrt(2 * N), bins=64, normed=True);
Explanation: (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.
End of explanation
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
Explanation: (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.
Solution
We chose N = 100 and M = 1000. The distribution is the Marčenko–Pastur distribution.
End of explanation
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)
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=':')
Explanation: (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.
Solution
My reasons are:
* A function specific for hermitian matrices is probably optimized to work faster than in the general case.
* It returns the eignvalues as real numbers and not complex with the immaginary part close to 0.
    * For this it returns the values sorted
* It helps to keep the code clean.
3.  Eigen functions of the quantum harmonic oscillator.
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.
End of explanation
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))
Explanation: (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}$.
End of explanation
<END_TASK>