Problemstellung

Oberschwingungsoszillatorpotential $ V (x) $

V(x) = \frac{m\omega^2x^2}{2} {\tag1}

Gegeben in.

Abbildung Harmonisches Oszillatorpotential. Die gepunktete Linie ist E = 1,5 Ry.

Die entsprechende stationäre Schrödinger-Gleichung lautet

\frac{d^2}{d x^2}\psi(x) +k^2(x) \psi =0 {\tag 2}

k^2(x) = \frac{2m}{\hbar^2}(E-\frac{m\omega^2x^2}{2}\) {\tag 3}

Ist.

Rydberg-Atomeinheit Einheitensystem genannt,

Elektronenmasse $ m = 1/2 $ Dirac-Konstante $ \ hbar = 1 $ Länge in Bohr $ a_ {B} = (0,529177 Å) $ Einheit, Energie $ 1 Ry = 13.6058 eV = $

Mit der Schrödinger-Gleichung

$\frac{d^2}{d x^2}\psi(x) +k^2(x) \psi =0 $

k^2(x) = (E-\frac{\omega^2x^2}{4}\) {\tag 4}

Wird sein.

Die genaue Lösung ** des ** th $ n (n = 0,1,2, ...) $ Eigenwerts $ E_n $ für das harmonische Oszillatorpotential ist

E_n = (n+\frac{1}{2})(\hbar \omega), (n=0,1,2,...){\tag 5} Ist. In diesem Code ist $ \ hbar = 1 $, $ \ omega = 1 $, also

E_0 = 0.5 E_1 = 1.5 E_2= 2.5 E_3= 3.5 ...

Ist.

** Die genaue Lösung der standardisierten Eigenfunktion $ \ psi_n $ ** lautet wie folgt.

\psi_n(x)=AH_n(\xi)\exp\left(-\frac{\xi^2}{2}\right) {\tag 6}

Wobei $ \ xi $ und $ A $ unten angegeben sind.

\xi=\sqrt{\frac{m\omega}{\hbar}}x,\ A=\sqrt{\frac{1}{n!2^n}\sqrt\frac{m\omega}{\pi\hbar}}, {\tag 7}

Außerdem ist $ H_n $ ein Hermeet-Polynom.

H_n(x)=(-1)^n\exp\left(x^2\right)\frac{\mathrm{d}^n}{\mathrm{d}x^n}\exp\left(-x^2\right) {\tag 8}

** Der Zweck dieses Papiers ist es, diese durch numerische Berechnung zu erhalten. ** **.

** Lösungsbewertungsmethode **

In dem in diesem Artikel veröffentlichten Code wird die Bewertungsmethode für die Lösung des Problems wie folgt festgelegt. ** Die Wellenfunktion ist eine kontinuierliche Funktion, und ihre Differenzierung erster Ordnung ist ebenfalls kontinuierlich **. Eine Anforderung, die diese Eigenschaft der Wellenfunktion erfüllt, ist ** Differential der logarithmischen Funktion der Wellenfunktion ** \frac{d (\ln\ \psi(x))}{dx} = \frac{\psi'(x)}{\psi(x)} {\tag 9}

Sie können die Bedingung auferlegen, dass ** kontinuierlich ist. Wenn die Wellenfunktion weiter standardisiert ist, kann schließlich die richtige Lösung erhalten werden. ** **.

Deshalb, Als Bewertungsfunktion der Lösung am Punkt $ x = x_m $ zur Überprüfung der von links und rechts integrierten Lösung ($ \ psi_L (x) $ bzw. $ \ psi_R (x) $)

\Delta E(\epsilon, x=x_m) =\frac{\psi_L'(x_m)}{\psi_L(x_m)}-\frac{\psi_R'(x_m)}{\psi_R(x_m)} {\tag {10}} Denk an. Wenn $ \ Delta E (\ epsilon, x = x_m) = 0 $ ist, ist $ \ epsilon $ der richtige Energieeigenwert. ** $ \ Delta E (\ epsilon, x = x_m) $ wird nacheinander aus einem bestimmten Energiebereich von Emin bis Emax als Funktion von $ \ epsilon $ berechnet, und wenn $ \ Delta E $ kleiner als ein bestimmter Wert wird Kann als Lösung angesehen werden. ** **.

Übrigens kann Gleichung (10) einen sehr großen Wert annehmen, was die numerische Auswertung schwierig machen kann. Um die Skala von $ \ Delta E (\ epsilon) $ anzupassen, ist in dem in diesem Artikel veröffentlichten Code der Nenner von Gleichung (10) multipliziert mit dem Skalenanpassungsfaktor (Summe der logarithmischen Differenzierung) $ \ Delta E. Es wird als $ angenommen. Natürlich bleibt die auf diese Weise erhaltene Lösung dieselbe.

\Delta E(\epsilon, x=x_m) =(\frac{\psi_L'(x_m)}{\psi_L(x_m)}-\frac{\psi_R'(x_m)}{\psi_R(x_m)})/(\frac{\psi_L'(x_m)}{\psi_L(x_m)}+\frac{\psi_R'(x_m)}{\psi_R(x_m)}) {\tag {11}}

Code

Ich habe im Anhang einige Punkte erwähnt, die zu beachten sind: (1) Auswahl des Werts von x zur Überprüfung der Konsistenz der Lösung an beiden Enden und (2) Auswahl der Symmetrie der Lösung und der Anfangsbedingung, damit ich die Details des Codes kenne. Wenn Sie möchten, lesen Sie es bitte. ** **.

"""
Aufnahmemethode+Zeitunabhängige eindimensionale Schrödinger-Gleichungslösung nach der Numerov-Methode
Harmonisches Oszillatorpotential
3 Sept. 2017
"""
import numpy as np
import matplotlib.pyplot as plt


iter_max = 100 
delta_x=0.01

xL0, xR0  = -10, 10

E=1.5

eps_E=0.01
#nn=4

hbar=1
omega=1
m_elec=0.5  

def calc_match_pos_osci(E):
    xa =np.sqrt(2*E/(m_elec*(omega**2)))#Wendepunkt
    return xa

xa = calc_match_pos_osci(E) #Wendepunkt
print("xa=",xa)
#xL0, xR0  = -nn*xa, nn*xa


Nx =  int((xR0-xL0)/delta_x)

i_match = int((xa-xL0)/delta_x)  #Ein Index der Position, um die Übereinstimmung zwischen uL und uR zu überprüfen. Als Wendepunkt auswählen.
nL = i_match
nR = Nx-nL

print(xL0,xR0, i_match, delta_x)
print(nL,nR)

uL = np.zeros([nL],float)
uR = np.zeros([nR],float)

print("E= ",E)
print("xL0,xR0, i_match, delta_x=",xL0,xR0, i_match, delta_x)
print("Nx, nL,nR=",Nx, nL,nR)




def V(x): #Harmonisches Oszillatorpotential
    v = m_elec*(omega**2)*(x**2)/2        
    return v

#Randbedingung / Anfangsbedingung eingestellt
def set_condition_even(): #Gleiche Funktion
    uL[0]  = 0
    uR[0] = 0

    uL[1] =  1e-12
    uR[1] =  1e-12

def set_condition_odd(): #Komische Funktion
    uL[0]  = 0
    uR[0] = 0

    uL[1] = -1e-12
    uR[1] =  1e-12

    
    

set_condition_even()


def setk2 (E): # for E<0
    for i in range(Nx+1):
        xxL = xL0 + i*delta_x
        xxR = xR0 -  i*delta_x
        k2L[i] = (2*m_elec/hbar**2)*(E-V(xxL))
        k2R[i] =(2*m_elec/hbar**2)*(E-V(xxR))
        
        
def Numerov (N,delta_x,k2,u):  #Berechnung nach der Numerov-Methode
    b = (delta_x**2)/12.0
    for i in range(1,N-1):
        u[i+1] = (2*u[i]*(1-5*b*k2[i])-(1+b*k2[i-1])*u[i-1])/(1+b*k2[i+1]) 
        
        
xL=np.zeros([Nx])
xR=np.zeros([Nx])

for i in range (Nx):
    xL[i] = xL0 + i*delta_x
    xR[i] = xR0 - i*delta_x
    
k2L=np.zeros([Nx+1])
k2R=np.zeros([Nx+1])

setk2(E)



def E_eval(): #Bewertungsfunktion:Formel(11)Sehen
#    print("in E_eval")
#    print("delta_x=",delta_x)

    if uL[-1]*uR[-1] >0 : #Wenn die Codes unterschiedlich sind, stimmen logderi möglicherweise versehentlich überein. Fügen Sie daher die Bedingung desselben Codes hinzu.(Es spielt keine Rolle, ob Sie nur den eindeutigen Wert finden)

        uLdash = (uL[-1]-uL[-2])/delta_x
        uRdash = (uR[-2]-uR[-1])/delta_x

        logderi_L=  uLdash/uL[-1]
        logderi_R=  uRdash/uR[-1]    
   #     print("logderi_L, R=",logderi_L,logderi_R)
    
        return (logderi_L- logderi_R)/(logderi_L+logderi_R) #Formel(9)
    else:
        return False






#Mögliches Funktionsdiagramm
XXX= np.linspace(xL0,xR0, Nx)
POT=np.zeros([Nx])
for i in range(Nx):
    POT[i] = V(xL0 + i*delta_x)
plt.xlabel('X (Bohr)') #x-Achsenbeschriftung
plt.ylabel('V (X) (Ry)') #y-Achsenbeschriftung
plt.hlines([E], xL0,xR0, linestyles="dashed")  #Energy
plt.plot(XXX,POT,'-',color='blue')
plt.show()
#

# k^2(x)Handlung
plt.xlabel('X (Bohr)') #x-Achsenbeschriftung
plt.ylabel('k^2 (X) (Ry)') #y-Achsenbeschriftung

XXX= np.linspace(xL0,xR0, len(k2L[:-2]))
plt.plot(XXX, k2L[:-2],'-')
plt.show()
#


def normarize_func(u):
    factor = ((xR0-xL0)/Nx)*(np.sum(u[1:-2]**2))
    return factor
def plot_eigenfunc(color_name):  
    uuu=np.concatenate([uL[0:nL-2],uR[::-1]],axis=0)
    XX=np.linspace(xL0,xR0, len(uuu))

    factor=np.sqrt(normarize_func(uuu))
 #   print("fcator = ",factor)
    plt.plot(XX,uuu/factor,'-',color=color_name,label='Psi')
    plt.plot(XX,(uuu/factor)**2,'-',color='red',label='| Psi |^2')
 
    plt.xlabel('X (Bohr)') #x-Achsenbeschriftung
    plt.ylabel('') #y-Achsenbeschriftung
    plt.legend(loc='upper right')
    plt.show()


   



#Eine Lösung finden

#Randbedingung 1(Gleiche Funktion)
EEmin=0.1
EEmax = 5
delta_EE=0.01

NE = int((EEmax-EEmin)/delta_EE)
Elis=[]
Solved_Eigenvalu=[]
check_Elis= []
for i in range(NE+1):
    EE=EEmin+i*(EEmax-EEmin)/NE
   

    xa = calc_match_pos_osci(EE) #Wendepunkt
    i_match = int((xa-xL0)/delta_x)  #Ein Index der Position, um die Übereinstimmung zwischen uL und uR zu überprüfen. Wählen Sie am Wendepunkt.
    nL = i_match
    nR = Nx-nL
    
    uL = np.zeros([nL],float)
    uR = np.zeros([nR],float)

    set_condition_even()
    setk2(EE)

    Numerov (nL,delta_x,k2L,uL)
    Numerov (nR,delta_x,k2R,uR)
 
    a1= E_eval()
    #print ("a1=",a1)
    if a1 :  #Wenn a1 wahr ist
        Elis.append(EE)
        check_Elis.append(a1)
        if np.abs(a1) <= eps_E :           
            print("Eigen_value = ", EE)
            Solved_Eigenvalu.append(EE)
            plot_eigenfunc("blue")
            
plt.plot(Elis, check_Elis, 'o',markersize=3, color='blue',linewidth=1)
plt.grid(True) #Erstellen Sie einen Diagrammrahmen
plt.xlim(EEmin, EEmax) #Der Bereich von x zum Zeichnen[xmin,xmax]Zu
plt.ylim(-10, 10) #Der Bereich von y zum Zeichnen[ymin,ymax]Zu
plt.hlines([0], EEmin,EEmax, linestyles="dashed")  # y=Zeichnen Sie eine gestrichelte Linie auf y1 und y2
plt.xlabel('Energy (Ry)') #x-Achsenbeschriftung
plt.ylabel('Delta_E_function') #y-Achsenbeschriftung
plt.show()

#Randbedingung 2(Komische Funktion)
EEmin=0.1
EEmax = 5
delta_EE=0.01

NE = int((EEmax-EEmin)/delta_EE)
Elis=[]
Solved_Eigenvalu=[]
check_Elis= []
for i in range(NE+1):
    EE=EEmin+i*(EEmax-EEmin)/NE
   

    xa = calc_match_pos_osci(EE) #Wendepunkt
    i_match = int((xa-xL0)/delta_x)  #Ein Index der Position, um die Übereinstimmung zwischen uL und uR zu überprüfen. Als Wendepunkt auswählen.
    nL = i_match
    nR = Nx-nL
    
    uL = np.zeros([nL],float)
    uR = np.zeros([nR],float)

    set_condition_odd()
    setk2(EE)

    Numerov (nL,delta_x,k2L,uL)
    Numerov (nR,delta_x,k2R,uR)
 
    a1= E_eval()
    #print ("a1=",a1)
    if a1 :  #Wenn a1 wahr ist
        Elis.append(EE)
        check_Elis.append(a1)
        if np.abs(a1) <= eps_E :           
            print("Eigen_value = ", EE)
            Solved_Eigenvalu.append(EE)
            plot_eigenfunc("blue")
    

plt.plot(Elis, check_Elis, 'o',markersize=3, color='red',linewidth=1)
plt.grid(True) #Erstellen Sie einen Diagrammrahmen
plt.xlim(EEmin, EEmax) #Der Bereich von x zum Zeichnen[xmin,xmax]Zu
plt.ylim(-10, 10) #Der Bereich von y zum Zeichnen[ymin,ymax]Zu
plt.hlines([0], EEmin,EEmax, linestyles="dashed")  # y=Zeichnen Sie eine gestrichelte Linie auf y1 und y2
plt.xlabel('Energy (Ry)') #x-Achsenbeschriftung
plt.ylabel('Delta_E_function') #y-Achsenbeschriftung
plt.show()

Ergebnis

(1) Bewertungsfunktion

Eine Auswertungsfunktion (ein Diagramm von Gleichung (9) ist gezeigt. Der Nullpunkt ist die Lösung des Energieeigenwerts. Es gibt zwei Typen, eine Lösung mit geraden Zahlen und eine Lösung mit ungeraden Funktionen.

Abbildung. Für gleichmäßige Funktionen

Abbildung: Bei ungerader Funktion

(2) Energieeigenwert und Wellenfunktion

Einige der erhaltenen Eigenwerte und Eigenfunktionen sind gezeigt. ** Sehr gut im Einklang mit der genauen Lösung. ** **.

Die genaue Lösung ist E_n = (n+\frac{1}{2})(\hbar \omega), (n=0,1,2,...) Ist. In diesem Code ist $ \ hbar = 1 $, $ \ omega = 1 $, also

E_0 = 0.5 E_1 = 1.5 E_2= 2.5 E_3= 3.5 ...

Ist.

Abbildung. Lösung von n = 0 (Grundzustand). Gleiche Funktion.

Abbildung Erster Zustand im angeregten Zustand (n = 1). Komische Funktion.

Abbildung Lösung für n = 2. Gleiche Funktion.

Abbildung Lösung für n = 3. Komische Funktion.

Nachtrag

** (1) Wo soll der Punkt ausgewählt werden, an dem die von beiden Seiten integrierte Lösung überprüft werden soll? **

Der Punkt $ x_m $, der die Konsistenz der von beiden Enden integrierten Lösung überprüft, wurde als Grenze (** als "Wendepunkt" ** bezeichnet) ausgewählt, an der klassische dynamische Bewegungen verboten sind. Das heißt, für die gegebene Schätzung von E.

E = m \omega^2 x^2 /2 {\tag {12}}

$ X $, das erfüllt, ist definiert als $ x_m $. Das ist, x_m= (\frac{2E}{m\omega^2})^\frac{1}{2} {\tag {13}}

Ist

Dies liegt daran, dass bei einer Integration von $ x_m $ die Lösung der Differentialgleichung ** mit einer Lösung (Sonderlösung) gemischt wird, die in der Entfernung divergiert, was zu einer numerischen Instabilität führt. ** </ front>