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. )
