Contenu du code: potentiel de puits profond

Tout le code est [unité atomique de Rydberg](https://ja.wikipedia.org/wiki/%E5%8E%9F%E5%AD%90%E5%8D%98%E4%BD%8D%E7%B3% BB) est utilisé. c'est,

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

Est d'être.

Trouvez la fonction d'onde et la valeur propre de l'énergie dans l'état lié. L'énergie à rechercher est E = 0-20 (Ry).

Figure 1. Potentiel bien formé. La largeur du puits est de 2 Bohr et la hauteur du puits est de 1000 Ry.

Pour les puits infiniment profonds, la solution exacte de la valeur propre d'énergie est

E_n = n^2 \pi^2/4 {\tag 2} Et

E1=2.46740110027234 (Ry) E2 = 9.869604401089358 (Ry) E3 = 22.20660990245106 (Ry) ... [1].

Code: Potentiel de puits profond

En raison du calcul du test, les conditions de calcul sont relativement souples. Si vous souhaitez améliorer la précision Il est bon de réduire la largeur de grille delta_x.

"""
Méthode de prise de vue+Résolution d'équation de Schrödinger unidimensionnelle indépendante du temps par la méthode Numerov
Potentiel profond et bien formé
1 Sept. 2017
"""
import numpy as np
import matplotlib.pyplot as plt

delta_x=0.05

xa =1 #Eh bien limite. X=Il y a une fin de puits à ± xa.
eps_E=0.005 #Condition de convergence

nn=5 #Xa les coordonnées des deux extrémités lors de l'intégration des deux extrémités(Quand-xa0)Paramètre indiquant combien de fois
xL0, xR0  = -nn*xa, nn*xa

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

delta_x = (xR0-xL0)/(Nx)
i_match = int((xa-xL0)/delta_x)  #Un index de la position pour vérifier la correspondance entre uL et uR. Je l'ai choisi comme limite du puits.


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)

E=np.pi**2/4

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 de type puits
    if np.abs(x) >  xa :
        v = 1000.0
    else :
        v = 0
        
    return v

#Condition aux limites / jeu de conditions initiales
def set_condition():
    uL[0]  = 0
    uL[1] =1e-6

    uR[0] = 0
    uR[1] =1e-6
#
set_condition()


def setk2 (E): # for E<0
    for i in range(Nx+1):
        xxL = xL0 + i*delta_x
        xxR = xR0 -  i*delta_x
        k2L[i] = E-V(xxL) 
        k2R[i] = E-V(xxR) 
        
        
def Numerov (N,delta_x,k2,u):  #Développement 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():
    uLdash = (uL[-1]-uL[-2])/delta_x
    uRdash = (uR[-2]-uR[-1])/delta_x
    logderi_L=  uLdash/uL[-1]
    logderi_R=  uRdash/uR[-1]    
    return (logderi_L- logderi_R)/(logderi_L+logderi_R)





#Diagramme des fonctions potentielles
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
XXX= np.linspace(xL0,xR0, Nx+1)
plt.plot(XXX, k2L,'-')
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))

    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 = 20
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


    set_condition_even()
    setk2(EE)

    Numerov (nL,delta_x,k2L,uL)
    Numerov (nR,delta_x,k2R,uR)
 
    a1= E_eval()
  
    if a1 :  
        Elis.append(EE)
        check_Elis.append(a1)
        if np.abs(a1) <= eps_E :  #Tracer lors de la recherche d'une solution
            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 = 20
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
   
    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)

(1-A) Calcul de la fonction d'évaluation

Figure: Pour des fonctions paires. Le point zéro est la valeur propre.

Figure: Pour les fonctions impaires.

(1-B) Fonction d'onde obtenue et énergie intrinsèque

Figure.Fonction d'onde d'état de base\psi(x)Et densité de probabilité|\psi(x)|^2

Solution exacte Par rapport à $ E_1 = 2,4674011 $, cela correspond à une erreur d'environ 4%.

Figure. Premier état excité

C'est à moins de 6% de la solution exacte. E_2 = 9.8696044 (Ry)

Addendum: méthode de prise de vue (méthode de prise de vue)

(1.Tout d'abord

Équation de Schrödinger unidimensionnelle \frac{d^2}{d x^2}\psi(x) +k^2(x) \psi =0 {\tag 1}

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

Est ** dans un état lié **

\lim_{x \to \infty} \psi(x) = 0 {\tag 3} Rencontrer.

** C'est la condition aux limites. ** Résolvez l'équation ci-dessus (1,2) sous la condition aux limites (3). ** À ce moment-là, il n'y a pas toujours de solution pour tout $ E $. La solution n'existe que lorsqu'un certain $ E = E_n $. C'est ce qu'on appelle une valeur unique. Et la fonction à ce moment-là s'appelle une fonction propre. ** ** Autrement dit, l'équation (1) aboutit à un problème mathématique appelé ** équation propre **. Nous l'appelons ** "problème de valeur unique" **.

(2) Procédure de la méthode de prise de vue

La méthode de prise de vue est l'une des solutions au problème des valeurs propres des équations différentielles. Considérez une fonction qui satisfait la condition aux limites, intégrez l'équation différentielle normale intégrale des deux extrémités en tant que problème de valeur initiale et vérifiez si les deux solutions correspondent à un certain point. ** Si les valeurs propres d'énergie correctes sont données, les deux solutions concordent. ** **

La méthode de prise de vue utilise cette propriété [2]. Suivez les étapes ci-dessous pour résoudre le problème de valeur propre.

1.Résolvez l'équation différentielle des deux extrémités en donnant une estimation de la valeur propre de l'énergie
2.Vérifiez si les deux fonctions de solution correspondent à un point approprié

3-1.S'ils ne correspondent pas(Figure 1)Modifie la valeur propre estimée de l'énergie et revient à 1
3-2.S'ils correspondent(Figure 2)A obtenu des valeurs propres d'énergie et des fonctions propres.

Figure 1: Deux fonctions d'onde obtenues en résolvant des équations différentielles de gauche et de droite. Les solutions ne sont pas connectées en douceur (correspondant à 3-1).

Figure 2: Lorsque les solutions sont connectées en douceur. L'énergie donnée à ce moment et la fonction d'onde obtenue sont les valeurs propres et les fonctions propres correctes.

Addendum: Comment évaluer la solution

À propos de cela, c'était mélangé,

[[Calcul scientifique / technique par Python] Résolution de l'équation de Schrodinger unidimensionnelle en régime permanent par la méthode de tir (2), potentiel d'oscillateur harmonique, mécanique quantique](http://qiita.com/sci_Haru/items/edb475a6d2eb7e901905# Méthode d'évaluation)

Prière de se référer à.

Les références

[1] Article Qiita [Calcul scientifique / technique par Python] Résolution d'équations différentielles ordinaires du second ordre par la méthode Numerov, calcul numérique