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

Answer 1

[indiachinacook]

Ich mag mir Deinen Code nicht ansehen und habe das einfach schnell mit Gnuplot gemacht:

Verschönerungen können immer noch gemacht werden.

piq=1/(sqrt(sqrt(pi)))
cutoff(x) = abs(x)<0.01 ? NaN : x
psi0(x)=cutoff(piq*exp(-0.5*x**2))
psi1(x)=cutoff(sqrt(2)*piq*x * exp(-0.5*x**2))
psi2(x)=cutoff(1/sqrt(2)*piq * (2*x**2-1) * exp(-0.5*x**2))
psi3(x)=cutoff(1/sqrt(3)*piq * (2*x**3-3*x) * exp(-0.5*x**2))
psi4(x)=cutoff(0.5/sqrt(6)*piq * (4*x**4-12*x**2+3) * exp(-0.5*x**2))
psi5(x)=cutoff(0.5/sqrt(15)*piq * (4*x**5-20*x**3+15*x) * exp(-0.5*x**2))

cutoff2(x) = x**2<0.01 ? NaN : 2*x**2
psi_0(x)=cutoff2(piq*exp(-0.5*x**2))
psi_1(x)=cutoff2(sqrt(2)*piq*x * exp(-0.5*x**2))
psi_2(x)=cutoff2(1/sqrt(2)*piq * (2*x**2-1) * exp(-0.5*x**2))
psi_3(x)=cutoff2(1/sqrt(3)*piq * (2*x**3-3*x) * exp(-0.5*x**2))
psi_4(x)=cutoff2(0.5/sqrt(6)*piq * (4*x**4-12*x**2+3) * exp(-0.5*x**2))
psi_5(x)=cutoff2(0.5/sqrt(15)*piq * (4*x**5-20*x**3+15*x) * exp(-0.5*x**2))

e0=1.0
niveau(x,e)=abs(x)<sqrt(2*e) ? e : NaN
bound=3.5

plot [-bound:bound] [0:bound**2/2] x**2/2 lw 3 title "", \
            e0/2*psi0(x)+0.5*e0 ls 2 title "", \
            e0/2*psi1(x)+1.5*e0 ls 3 title "", \
            e0/2*psi2(x)+2.5*e0 ls 4 title "", \
            e0/2*psi3(x)+3.5*e0 ls 5 title "", \
            e0/2*psi4(x)+4.5*e0 ls 6 title "", \
            e0/2*psi5(x)+5.5*e0 ls 7 title "", \
            niveau (x,0.5*e0) ls 2 lw 2 title "", \
            niveau (x,1.5*e0) ls 3 lw 2 title "", \
            niveau (x,2.5*e0) ls 4 lw 2 title "", \
            niveau (x,3.5*e0) ls 5 lw 2 title "", \
            niveau (x,4.5*e0) ls 6 lw 2 title "", \
            niveau (x,5.5*e0) ls 7 lw 2 title ""

pause -1

plot [-bound:bound]  [0:bound**2/2] x**2/2 lw 3 title "", \
            e0/2*psi_0(x)+0.5*e0 ls 2 title "", \
            e0/2*psi_1(x)+1.5*e0 ls 3 title "", \
            e0/2*psi_2(x)+2.5*e0 ls 4 title "", \
            e0/2*psi_3(x)+3.5*e0 ls 5 title "", \
            e0/2*psi_4(x)+4.5*e0 ls 6 title "", \
            e0/2*psi_5(x)+5.5*e0 ls 7 title "", \
            niveau (x,0.5*e0) ls 2 lw 2 title "", \
            niveau (x,1.5*e0) ls 3 lw 2 title "", \
            niveau (x,2.5*e0) ls 4 lw 2 title "", \
            niveau (x,3.5*e0) ls 5 lw 2 title "", \
            niveau (x,4.5*e0) ls 6 lw 2 title "", \
            niveau (x,5.5*e0) ls 7 lw 2 title ""

pause -1

  

  

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

Comments

[UnknownUser1192]

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