Problem setting

Harmonic oscillator potential $ V (x) $ is

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

Given in.

Figure. Harmonic oscillator potential. The dotted line is E = 1.5Ry.

The corresponding steady-state Schrodinger equation is

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

Is.

Rydberg Atomic Unit Unit system called,

Electron mass $ m = 1/2 $ Dirac constant $ \ hbar = 1 $ Length to Bohr $ a_ {B} = (0.529177 Å) $ unit, Energy $ 1 Ry = 13.6058 eV = $

Using, the Schrodinger equation

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

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

Will be.

The exact solution ** of the ** th $ n (n = 0,1,2, ...) $ eigenvalue $ E_n $ for the harmonic oscillator potential is

E_n = (n+\frac{1}{2})(\hbar \omega), (n=0,1,2,...){\tag 5} Is. In this code, $ \ hbar = 1 $, $ \ omega = 1 $, so

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

Is.

** The exact solution of the normalized eigenfunction $ \ psi_n $ ** is as follows.

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

Where $ \ xi $ and $ A $ are given below.

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

Also, $ H_n $ is a Hermite polynomial.

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}

** The purpose of this paper is to obtain these by numerical calculation. ** **

** Evaluation method of solution **

In the code posted in this article, the evaluation method of the solution of the problem is set as follows. ** The wave function is a continuous function, and its first derivative is also continuous **. One requirement to satisfy this property of the wavefunction is ** logarithmic derivative of the wavefunction ** \frac{d (\ln\ \psi(x))}{dx} = \frac{\psi'(x)}{\psi(x)} {\tag 9}

You can impose the condition that ** is continuous. If the wave function is further standardized, the correct solution can be finally obtained. ** **

Therefore, As an evaluation function of the solution at the point $ x = x_m $ to check the solution integrated from the left side and the right side ($ \ psi_L (x) $ and $ \ psi_R (x) $, respectively)

\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}} think of. When $ \ Delta E (\ epsilon, x = x_m) = 0 $, $ \ epsilon $ is the correct energy eigenvalue. ** $ \ Delta E (\ epsilon, x = x_m) $ is calculated one after another from a certain energy range Emin to Emax as a function of $ \ epsilon $, and when $ \ Delta E $ becomes smaller than a certain value. Can be regarded as a solution. ** **

By the way, equation (10) may take a very large value, which may make numerical evaluation difficult. Therefore, in order to adjust the scale of $ \ Delta E (\ epsilon) $, in the code published in this article, the denominator of equation (10) multiplied by the scale adjustment factor (sum of logarithmic derivatives) is $ \ Delta E. It is adopted as $. Of course, the solution obtained in this way remains the same.

\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

I mentioned some points to note in the appendix about (1) how to select the value of x to check the consistency of the solution from both ends, and (2) how to select the symmetry of the solution and the initial condition, so I know the details of the code. If you want, please read it. ** **

"""
Shooting method+Numerov method for solving time-independent one-dimensional Schrodinger equation
Harmonic oscillator potential
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)))#Turnoff point
    return xa

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


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

i_match = int((xa-xL0)/delta_x)  #Index of the position to check the match between uL and uR. Select as a turning point.
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): #Harmonic oscillator potential
    v = m_elec*(omega**2)*(x**2)/2        
    return v

#Boundary condition / initial condition set
def set_condition_even(): #Even function
    uL[0]  = 0
    uR[0] = 0

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

def set_condition_odd(): #Odd function
    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):  #Calculation by Numerov method
    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(): #Evaluation function:formula(11)See
#    print("in E_eval")
#    print("delta_x=",delta_x)

    if uL[-1]*uR[-1] >0 : #If the signs are different, logderi may accidentally match, so add the condition of the same sign.(It doesn't matter if you just find the eigenvalues)

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






#Potential function plot
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-axis label
plt.ylabel('V (X) (Ry)') #y-axis label
plt.hlines([E], xL0,xR0, linestyles="dashed")  #Energy
plt.plot(XXX,POT,'-',color='blue')
plt.show()
#

# k^2(x)Plot
plt.xlabel('X (Bohr)') #x-axis label
plt.ylabel('k^2 (X) (Ry)') #y-axis label

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-axis label
    plt.ylabel('') #y-axis label
    plt.legend(loc='upper right')
    plt.show()


   



#Search for a solution

#Boundary condition 1(Even function)
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) #Turnoff point
    i_match = int((xa-xL0)/delta_x)  #Index of the position to check the match between uL and uR. Select at the turning point.
    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 :  #When a1 is True
        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) #Create a frame for the graph
plt.xlim(EEmin, EEmax) #The range of x to draw[xmin,xmax]To
plt.ylim(-10, 10) #The range of y to draw[ymin,ymax]To
plt.hlines([0], EEmin,EEmax, linestyles="dashed")  # y=Draw dashed lines on y1 and y2
plt.xlabel('Energy (Ry)') #x-axis label
plt.ylabel('Delta_E_function') #y-axis label
plt.show()

#Boundary condition 2(Odd function)
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) #Turnoff point
    i_match = int((xa-xL0)/delta_x)  #Index of the position to check the match between uL and uR. Select as a turning point.
    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 :  #When a1 is True
        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) #Create a frame for the graph
plt.xlim(EEmin, EEmax) #The range of x to draw[xmin,xmax]To
plt.ylim(-10, 10) #The range of y to draw[ymin,ymax]To
plt.hlines([0], EEmin,EEmax, linestyles="dashed")  # y=Draw dashed lines on y1 and y2
plt.xlabel('Energy (Ry)') #x-axis label
plt.ylabel('Delta_E_function') #y-axis label
plt.show()

result

(1) Evaluation function

The evaluation function (a plot of Eq. (9) is shown. The zero point is the solution of the energy eigenvalue. There are two types, the even number solution and the odd function solution.

Figure. For even functions

Figure. In the case of odd function

(2) Energy eigenvalues and wave functions

Some of the obtained eigenvalues and eigenfunctions are shown. ** Very well in line with the exact solution. ** **

The exact solution is E_n = (n+\frac{1}{2})(\hbar \omega), (n=0,1,2,...) Is. In this code, $ \ hbar = 1 $, $ \ omega = 1 $, so

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

Is.

Figure. Solution of n = 0 (ground state). Even function.

Figure. First excited state state (n = 1). Odd function.

Figure. Solution for n = 2. Even function.

Figure. Solution for n = 3. Odd function.

Addendum

** (1) Where to choose the point to check the solution integrated from both ends **

The point $ x_m $ for checking the consistency of the solution integrated from both ends was selected as the boundary (** called "turning point" **) where classical mechanical motion is prohibited. That is, for the given estimate of E

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

$ X $ that satisfies is defined as $ x_m $. That is, x_m= (\frac{2E}{m\omega^2})^\frac{1}{2} {\tag {13}}

Is

This is because if the integral exceeds $ x_m $, the solution of the ** differential equation will be mixed with a solution that diverges in the distance (special solution), and it will be numerically unstable. ** </ front>