import numpy as np
import matplotlib.pyplot as plt

# Define a function to create a random Hermitian matrix of size n
def create_random_hermitian_matrix(n):
    A = np.random.rand(n, n) + 1j * np.random.rand(n, n)
    return A + A.conj().T

# Define the size of the matrix
n = 10

# Create the unperturbed density matrix (ground state)
rho_0_large = create_random_hermitian_matrix(n)
rho_0_large = (rho_0_large + rho_0_large.conj().T) / 2  # Ensure it's Hermitian

# Create a small perturbation matrix
perturbation_large = create_random_hermitian_matrix(n) * 0.05  # Scale it to be a small perturbation

# Linear response: rho(t) = rho_0 + perturbation
rho_linear_large = rho_0_large + perturbation_large

# Multiplicative response: rho(t) = rho_0 * perturbation (element-wise multiplication)
rho_multiplicative_large = rho_0_large * perturbation_large

# Compute the eigenvalues
eigenvalues_rho_0_large = np.linalg.eigvals(rho_0_large)
eigenvalues_linear_large = np.linalg.eigvals(rho_linear_large)
eigenvalues_multiplicative_large = np.linalg.eigvals(rho_multiplicative_large)

# Plot the eigenvalues
fig, ax = plt.subplots(figsize=(12, 6))
x = np.arange(n)

ax.plot(x, eigenvalues_rho_0_large.real, 'o-', label='rho_0')
ax.plot(x, eigenvalues_linear_large.real, 's-', label='Linear (rho_0 + perturbation)')
ax.plot(x, eigenvalues_multiplicative_large.real, 'd-', label='Multiplicative (rho_0 * perturbation)')

ax.set_xticks(x)
ax.set_xticklabels([f'Eigenvalue {i+1}' for i in x])
ax.set_ylabel('Eigenvalues (Real Part)')
ax.set_title('Eigenvalues of the Density Matrix for a Large System')
ax.legend()

plt.show()


