I m trying to find the eigenstates of the Hamiltonian of the hydrogene atom in a 1D problem. I'm defining my space grid, then my Hamiltonian and finally I use the function eigh from scipy.linalg to find the eigenvalues and eigenvector (see after). But, when I'm doing so, I'm finding some weird stuff: the eigenvector are coherent with quantum theory while the half the eigenvalues are wrong (some should not exist and some should exist but are linked to the wrong eigenvector). On the photo, please note that the units are in atomic units. Therefore my question is: why do I get such weird stuff? You will find the code I've used below. Thanks a lot for your answer import numpy as np from scipy.sparse import diags from scipy.linalg import eigh import matplotlib.pyplot as plt # Parameters L = 50.0 # domain size N = 1000 # number of grid points x = np.linspace(-L/2, L/2, N) dx = x[1] - x[0] # 1D Hydrogen-like potential (softened to avoid singularity at x=0) Z = 1.0 epsilon = 1e-2 V = -Z / np.sqrt(x**2 + epsilon**2) # Kinetic energy operator (finite difference, central) main_diag = np.full(N, -2.0) off_diag = np.full(N-1, 1.0) laplacian = diags([off_diag, main_diag, off_diag], [-1, 0, 1]) / dx**2 # Hamiltonian H = -0.5 * laplacian.toarray() + np.diag(V) # Solve eigenvalue problem num_states = 10 eigvals, eigvecs = eigh(H) eigvals = eigvals[:num_states] eigvecs = eigvecs[:, :num_states] # Normalize eigenfunctions eigvecs /= np.sqrt(np.sum(np.abs(eigvecs)**2, axis=0) * dx) # Plot plt.figure(figsize=(8,6)) for i in range(num_states): plt.plot(x, eigvecs[:, i] + eigvals, label=f'n={i+1}, E={eigvals:.3f}') plt.xlabel('x') plt.ylabel('Wavefunction + Energy') plt.title('1D Hydrogen Atom Eigenstates') plt.legend() plt.grid() plt.figure() plt.plot(x, eigvecs[:, 0]) plt.plot(x, np.exp(-np.abs(x)), label='Analytical 1s state', linestyle='--', color='red') plt.show() Continue reading...
