Problème de réglage

Potentiel d'oscillateur harmonique $ V (x) $

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

Donné dans.

Figure: Potentiel de l'oscillateur harmonique. La ligne pointillée est E = 1,5 Ry.

L'équation de Schrödinger en régime permanent correspondante est

\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}

Est.

Unité atomique de Rydberg Système d'unité appelé,

Masse électronique $ m = 1/2 $ Constante de Dirac $ \ hbar = 1 $ Longueur en Bohr $ a_ {B} = (0,529177 Å) $ unité, Énergie 1 $ Ry = 13,6058 eV = $

En utilisant, l'équation de Schrödinger

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

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

Sera.

La solution exacte ** de la ** th $ n (n = 0,1,2, ...) $ valeur propre $ E_n $ pour le potentiel d'oscillateur harmonique est

E_n = (n+\frac{1}{2})(\hbar \omega), (n=0,1,2,...){\tag 5} Est. Dans ce code, $ \ hbar = 1 $, $ \ omega = 1 $, donc

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

Est.

** La solution exacte de la fonction propre normalisée $ \ psi_n $ ** est la suivante.

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

Où $ \ xi $ et $ A $ sont donnés ci-dessous.

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

De plus, $ H_n $ est un polynôme Hermeet.

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}

** Le but de cet article est de les obtenir par calcul numérique. ** **

** Méthode d'évaluation de la solution **

Dans le code publié dans cet article, la méthode d'évaluation de la solution du problème est définie comme suit. ** La fonction d'onde est une fonction continue, et sa différenciation du premier ordre est également continue **. Une condition qui satisfait cette propriété de la fonction d'onde est ** différentielle de la fonction logarithmique de la fonction d'onde ** \frac{d (\ln\ \psi(x))}{dx} = \frac{\psi'(x)}{\psi(x)} {\tag 9}

Vous pouvez imposer la condition que ** est continue. Si la fonction d'onde est davantage normalisée, la solution correcte peut finalement être obtenue. ** **

Donc, En tant que fonction d'évaluation de la solution au point $ x = x_m $ pour vérifier la solution intégrée du côté gauche et du côté droit ($ \ psi_L (x) $ et $ \ psi_R (x) $, respectivement)

\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}} penser à. Lorsque $ \ Delta E (\ epsilon, x = x_m) = 0 $, $ \ epsilon $ est la valeur propre d'énergie correcte. ** $ \ Delta E (\ epsilon, x = x_m) $ est calculé l'un après l'autre à partir d'une certaine plage d'énergie Emin à Emax en fonction de $ \ epsilon $, et lorsque $ \ Delta E $ devient plus petit qu'une certaine valeur Peut être considéré comme une solution. ** **

À propos, l'équation (10) peut prendre une valeur très grande, ce qui peut rendre l'évaluation numérique difficile. Par conséquent, pour ajuster l'échelle de $ \ Delta E (\ epsilon) $, dans le code publié dans cet article, le dénominateur de l'équation (10) multiplié par le facteur d'ajustement d'échelle (somme de la différenciation logarithmique) est $ \ Delta E. Il est adopté comme $. Bien entendu, la solution ainsi obtenue reste la même.

\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

J'ai mentionné quelques points à noter dans l'annexe sur (1) comment sélectionner la valeur de x pour vérifier la cohérence de la solution des deux côtés, et (2) comment sélectionner la symétrie de la solution et la condition initiale, donc je connais les détails du code. Si vous le souhaitez, veuillez le lire. ** **

"""
Méthode de prise de vue+Résolution d'équation de Schrödinger unidimensionnelle indépendante du temps par la méthode Numerov
Potentiel d'oscillateur harmonique
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)))#Tournant
    return xa

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


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

i_match = int((xa-xL0)/delta_x)  #Un index de la position pour vérifier la correspondance entre uL et uR. Sélectionnez comme point de virage.
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): #Potentiel d'oscillateur harmonique
    v = m_elec*(omega**2)*(x**2)/2        
    return v

#Condition aux limites / jeu de conditions initiales
def set_condition_even(): #Même fonction
    uL[0]  = 0
    uR[0] = 0

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

def set_condition_odd(): #Fonction impaire
    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):  #Calcul par la méthode Numerov
    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(): #Fonction d'évaluation:formule(11)Voir
#    print("in E_eval")
#    print("delta_x=",delta_x)

    if uL[-1]*uR[-1] >0 : #Si les codes sont différents, logderi peut accidentellement correspondre, alors ajoutez la condition du même code.(Peu importe si vous trouvez simplement la valeur unique)

        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) #formule(9)
    else:
        return False






#Diagramme de fonction potentielle
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)') #étiquette de l'axe des x
plt.ylabel('V (X) (Ry)') #étiquette de l'axe y
plt.hlines([E], xL0,xR0, linestyles="dashed")  #Energy
plt.plot(XXX,POT,'-',color='blue')
plt.show()
#

# k^2(x)Terrain
plt.xlabel('X (Bohr)') #étiquette de l'axe des x
plt.ylabel('k^2 (X) (Ry)') #étiquette de l'axe y

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)') #étiquette de l'axe des x
    plt.ylabel('') #étiquette de l'axe y
    plt.legend(loc='upper right')
    plt.show()


   



#Trouver une solution

#Condition aux limites 1(Même fonction)
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) #Tournant
    i_match = int((xa-xL0)/delta_x)  #Un index de la position pour vérifier la correspondance entre uL et uR. Sélectionnez au tournant.
    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 :  #Quand a1 est vrai
        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) #Créer un cadre graphique
plt.xlim(EEmin, EEmax) #La plage de x à dessiner[xmin,xmax]À
plt.ylim(-10, 10) #La plage de y à dessiner[ymin,ymax]À
plt.hlines([0], EEmin,EEmax, linestyles="dashed")  # y=Tracez une ligne brisée sur y1 et y2
plt.xlabel('Energy (Ry)') #étiquette de l'axe des x
plt.ylabel('Delta_E_function') #étiquette de l'axe y
plt.show()

#Condition aux limites 2(Fonction impaire)
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) #Tournant
    i_match = int((xa-xL0)/delta_x)  #Un index de la position pour vérifier la correspondance entre uL et uR. Sélectionnez comme point de virage.
    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 :  #Quand a1 est vrai
        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) #Créer un cadre graphique
plt.xlim(EEmin, EEmax) #La plage de x à dessiner[xmin,xmax]À
plt.ylim(-10, 10) #La plage de y à dessiner[ymin,ymax]À
plt.hlines([0], EEmin,EEmax, linestyles="dashed")  # y=Tracez une ligne brisée sur y1 et y2
plt.xlabel('Energy (Ry)') #étiquette de l'axe des x
plt.ylabel('Delta_E_function') #étiquette de l'axe y
plt.show()

résultat

(1) Fonction d'évaluation

Une fonction d'évaluation (un tracé de l'équation (9) est représentée. Le point zéro est la solution de la valeur propre de l'énergie. Il existe deux types, une solution de nombre pair et une solution de fonction impaire.

Figure. Pour des fonctions paires

Figure. En cas de fonction impaire

(2) Valeur propre de l'énergie et fonction d'onde

Certaines des valeurs propres et des fonctions propres obtenues sont affichées. ** Très bien en ligne avec la solution exacte. ** **

La solution exacte est E_n = (n+\frac{1}{2})(\hbar \omega), (n=0,1,2,...) Est. Dans ce code, $ \ hbar = 1 $, $ \ omega = 1 $, donc

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

Est.

Figure: Solution de n = 0 (état basal). Même fonction.

Figure: Premier état d'état excité (n = 1). Fonction étrange.

Figure: Solution pour n = 2. Même fonction.

Figure: Solution pour n = 3. Fonction étrange.

Addenda

** (1) Où choisir le point pour vérifier la solution intégrée des deux côtés **

Le point $ x_m $, qui vérifie la cohérence de la solution intégrée des deux côtés, a été sélectionné comme frontière (** appelée "point de retournement" **) où le mouvement dynamique classique est interdit. Autrement dit, pour l'estimation donnée de E

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

$ X $ qui satisfait est défini comme $ x_m $. C'est, x_m= (\frac{2E}{m\omega^2})^\frac{1}{2} {\tag {13}}

Est

En effet, si l'intégration dépasse $ x_m $, la solution de l'équation différentielle ** sera mélangée à une solution (solution spéciale) qui diverge dans la distance, entraînant une instabilité numérique. ** </ avant>