Python und Physik: Harmonischer Oszillator graphisch bis n= 4?
Ich würde gerne die Vorlesungen nacharbeiten von der Uni. Wir hatten das Thema Harmonischer Oszillator und ich würde gerne so ein Plot erstellen, was ungefähr so aussieht Allerdings bis n=4. Also einmal die stationären (psi_n(x)) und einmal die Aufenthaltswahrscheinlichkeit |psi_n(x)|^2 Wobei die Parabel das Potential darstellen soll. Das problem ist, dass ich das irgendwie nicht richtig gesplottet kriege. Meine Graphen sehen immer wieder komisch aus und weiß einfach nicht warum. Kann jemand mir helfen das zu Plotten? Also bis jetzt sieht mein Code so aus: import numpy as np import matplotlib.pyplot as plt from scipy.special import hermite # Konstanten e = 1.602176634e-19 # Elementarladung [C] m = 9.10938356e-31 # Masse des Elektrons hq = 1.0545718e-34 # Wirkungsquantum # Parameter x0 = 1e-10 # Klassischer Umkehrpunkt in m V0 = 7 * e # Potential bei einer Auslenkung um x0 NN = 4 # Anzahl von Wellenfunktionen N = np.arange(NN) c = V0 / x0**2 # Kraftkonstante Oszillators omega = np.sqrt(c / m) # Frequenz des Oszillators # Maximales benötigtes x xmax = np.sqrt(2 * hq / (m * omega) * (NN + 1/2)) x = np.linspace(-xmax, xmax, 200) y = np.sqrt(m * omega / hq) * x # Potential V = 0.5 * c * x**2 plt.figure(figsize=(7, 10)) plt.plot(x * 1e10, V / e, linewidth=1.5) plt.xlabel('x [Angstrom]', fontsize=20) plt.ylabel('V [eV]', fontsize=20) plt.grid(True) # Energieeigenwerte E = np.zeros((NN, len(x))) for n in N: E[n, :] = (n + 0.5) * hq * omega / e plt.plot(x * 1e10, E[N, :].T, 'g', linewidth=1.5) plt.ylim(0, E[-1, 0] * 1.2) #geraden Potenzreihen a = np.zeros((len(N) + 2, len(N) + 1)) for n in range(0, len(N), 2): a[0, n] = 1 for j in range(0, len(N), 2): a[j + 2, n] = 1 # ungerade Potenzreihen for n in range(1, len(N), 2): a[1, n] = 1 for j in range(1, len(N), 2): a[j + 2, n] = 1 # Wellenfunktionen psi(n,y) psi = np.zeros((len(N) + 1, len(x))) for n in range(len(N)): psi[n, :] = 0 for j in range(len(N)): psi[n, :] += a[j, n] * hermite(j)(y) psi[n, :] *= np.exp(-y**2 / 2) # Normierung der Wellenfunktionen dx = x[1]-x[0] psi[n, :] /= np.sqrt(dx * np.sum(psi[n, :]**2)) # Skalierung der Wellenfunktionen dE = E[1, 0] - E[0, 0] psi_max = np.max(psi) fact = dE / psi_max / 2 plt.plot(x * 1e10, fact * psi[N, :].T + E[N, :].T, 'r', linewidth=2) # Aufenthaltswahrscheinlichkeit plt.figure(figsize=(7, 10)) plt.plot(x * 1e10, V / e, linewidth=1.5) plt.xlabel('x [Angstrom]', fontsize=20) plt.ylabel('V [eV]', fontsize=20) plt.grid(True) plt.ylim(0, E[-1, 0] * 1.2) plt.plot(x * 1e10, E[N, :].T, 'g', linewidth=1.5) # Skalierung der Aufenthaltswahrscheinlichkeiten fact = dE / psi_max**2 / 1.2 plt.plot(x * 1e10, fact * (psi[N, :].T)**2 + E[N, :].T, 'k', linewidth=2) plt.show() (Bitte ignoriert erstmal die Zeilenabstände. )
