Code content: Deep well-shaped potential

All code is [Rydberg atomic unit](https://ja.wikipedia.org/wiki/%E5%8E%9F%E5%AD%90%E5%8D%98%E4%BD%8D%E7%B3% BB) is used. this is,

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

Is to be.

Find the wavefunction and energy eigenvalues of the bound state. The energy to be searched is E = 0-20 (Ry).

Figure 1. Well-shaped potential. The width of the well is 2 Bohr and the height of the well is 1000 Ry.

For infinitely deep wells, the exact solution to the energy eigenvalues is

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

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

Code: Deep well potential

Since it is a test calculation, the calculation conditions are relatively loose. If you want to improve the accuracy It is good to reduce the grid width delta_x.

"""
Shooting method+Numerov method for solving time-independent one-dimensional Schrodinger equation
Deep well-shaped potential
1 Sept. 2017
"""
import numpy as np
import matplotlib.pyplot as plt

delta_x=0.05

xa =1 #Well boundary. x=There is a well end at ± xa.
eps_E=0.005 #Convergence condition

nn=5 #Xa the coordinates of both ends when integrating from both ends(When-xa0)Parameter indicating how many times
xL0, xR0  = -nn*xa, nn*xa

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

delta_x = (xR0-xL0)/(Nx)
i_match = int((xa-xL0)/delta_x)  #Index of the position to check the match between uL and uR. I chose it as the boundary of the well.


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): #Setting of well type potential
    if np.abs(x) >  xa :
        v = 1000.0
    else :
        v = 0
        
    return v

#Boundary condition / initial condition set
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):  #Development by the 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():
    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)





#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
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)') #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 = 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 :  #Plot when finding a 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) #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 = 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 :  #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)

(1-A) Calculation of evaluation function

Figure. For even functions. The zero point is the eigenvalue.

Figure. For odd functions.

(1-B) Obtained wavefunction and intrinsic energy

Figure.Ground state wavefunction\psi(x)And probability density|\psi(x)|^2

Exact solution It matches with an error of about 4% compared to $ E_1 = 2.4674011 $.

Figure. First excited state

It is within 6% of the exact solution. E_2 = 9.8696044 (Ry)

Addendum: Shooting method (shooting method)

(1.First of all

One-dimensional Schrodinger equation \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}

Is ** in bound state **

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

** This is a boundary condition. ** Solve the above equation (1,2) under the boundary condition (3). ** At that time, there is not always a solution for any $ E $. The solution exists only when a certain $ E = E_n $. It is called an eigenvalue. And the function at that time is called an eigenfunction. ** ** That is, equation (1) results in a mathematical problem called ** eigen equation **. We call it ** "eigenvalue problem" **.

(2) Shooting method procedure

The shooting method is one of the solutions to the eigenvalue problem of differential equations. Consider a function that satisfies the boundary conditions, integrate the ordinary differential equations from both ends as an initial value problem, and investigate whether the solutions of the two match at a certain point. ** If the correct energy eigenvalues are given, both solutions will agree. ** **

The shooting method utilizes this property [2]. Follow the steps below to solve the eigenvalue problem.

1.Solve the differential equation from both ends by giving an estimate of the energy eigenvalue
2.Find out if both solution functions match at an appropriate point

3-1.If they do not match(Figure 1)Changes the estimated energy eigenvalues and returns to 1
3-2.If they match(Figure 2)Has obtained energy eigenvalues and eigenfunctions.

Figure 1: Two wavefunctions obtained by solving differential equations from the left and right. The solutions are not connected smoothly (corresponding to 3-1).

Figure 2: When the solutions are connected smoothly. The energy given at this time and the obtained wave function are the correct eigenvalues and eigenfunctions.

Addendum: How to evaluate the solution

About this, it was mixed up,

[[Science and technology calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (2), Harmonic oscillator potential, Quantum mechanics](http://qiita.com/sci_Haru/items/edb475a6d2eb7e901905#Solution Evaluation method)

Please refer to.

References

[1] Qiita article [Scientific and technical calculation by Python] Solving second-order ordinary differential equations by Numerov method, numerical calculation