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. )
(1 votes)
Loading...
Ich mag mir Deinen Code nicht ansehen und habe das einfach schnell mit Gnuplot gemacht:
Verschönerungen können immer noch gemacht werden.
P.S.: Die erste Version hatte den Faktor ½ im Potential verschustert. Beachte, daß das immer noch nur für eine einzige Wahl des Potentialparameters k und der reduzierten Masse m gilt, nämlich k=m=1 in atomaren Einheiten. Wenn Dir das nicht gefällt, dann mußt Du k und m in die Koordinaten einbauen und beachten, daß e0=sqrt(k/m).
Omg, du bist ein Engel! Danke danke..
Ich habe es mir tatsächlich schon gedacht, das auch etwas falsch ist an den Komponenten aber ich konnte nicht erklären warum lul