[Scientific / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (1), well-type potential, quantum mechanics

Introduction

Using Python, ** [Time-independent one-dimensional Schrodinger equation] for a given potential $ V (x) $ (https://ja.wikipedia.org/wiki/%E3%82%B7%E3 % 83% A5% E3% 83% AC% E3% 83% BC% E3% 83% 87% E3% 82% A3% E3% 83% B3% E3% 82% AC% E3% 83% BC% E6% 96 % B9% E7% A8% 8B% E5% BC% 8F # .E6.99.82.E9.96.93.E3.81.AB.E4.BE.9D.E5.AD.98.E3.81.97.E3.81. AA.E3.81.84.E3.82.B7.E3.83.A5.E3.83.AC.E3.83.BC.E3.83.87.E3.82.A3.E3.83.B3.E3.82. AC.E3.83.BC.E6.96.B9.E7.A8.8B.E5.BC.8F) **

\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 solved by a numerical calculation method called the shooting method [addendum] to obtain the wave function and energy eigenvalues. ** **

In this article, ** a typical problem in the field [well-shaped potential](https://ja.wikipedia.org/wiki/%E4%BA%95%E6%88%B8%E5%9E%8B% Consider E3% 83% 9D% E3% 83% 86% E3% 83% B3% E3% 82% B7% E3% 83% A3% E3% 83% AB) **. The ** Numerov method [1] ** is used to solve the differential equations.

The shooting method is a conventional method that solves the boundary value problem of ordinary differential equations as an initial value problem, and has been used for a long time.

** Python code using the shooting method is hard to find on the net (as of September 4, 2017). We hope that this article will contribute to the reader's understanding of quantum mechanics. ** </ font>


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).

t1.png 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

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

s2.png Figure. For odd functions.


(1-B) Obtained wavefunction and intrinsic energy

スクリーンショット 2017-09-04 0.25.47.png 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 $.

スクリーンショット 2017-09-04 0.25.56.png 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.

スクリーンショット 2017-09-04 1.28.38.png Figure 1: Two wavefunctions obtained by solving differential equations from the left and right. The solutions are not connected smoothly (corresponding to 3-1).

スクリーンショット 2017-09-04 1.28.56.png 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

Recommended Posts

[Scientific / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (1), well-type potential, quantum mechanics
[Scientific / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (2), harmonic oscillator potential, quantum mechanics
[Scientific / technical calculation by Python] Solving one-dimensional Newton equation by the 4th-order Runge-Kutta method
[Scientific / technical calculation by Python] Fitting by nonlinear function, equation of state, scipy
[Scientific / technical calculation by Python] Solving second-order ordinary differential equations by Numerov method, numerical calculation
Solving one-dimensional wave equation by finite difference method (python)
Peeping in Python, not scary quantum mechanics 2: Finite well-type potential
[Scientific / technical calculation by Python] Numerical solution of one-dimensional harmonic oscillator problem by velocity Verlet method
[Scientific / technical calculation by Python] Solving simultaneous linear equations, numerical calculation, numpy
[Scientific / technical calculation by Python] Sum calculation, numerical calculation
[Scientific / technical calculation by Python] Numerical solution of two-dimensional Laplace-Poisson equation by Jacobi method for electrostatic potential, elliptic partial differential equation, boundary value problem
[Scientific / technical calculation by Python] Numerical solution of one-dimensional unsteady heat conduction equation by Crank-Nicholson method (implicit method) and FTCS method (positive solution method), parabolic partial differential equation
[Scientific / technical calculation by Python] Solving the boundary value problem of ordinary differential equations in matrix format, numerical calculation
[Scientific / technical calculation by Python] Solving (generalized) eigenvalue problem using numpy / scipy, using library
[Scientific / technical calculation by Python] Analytical solution to find the solution of equation sympy
[Scientific / technical calculation by Python] Inverse matrix calculation, numpy
[Scientific / technical calculation by Python] Histogram, visualization, matplotlib
[Scientific / technical calculation by Python] Lagrange interpolation, numerical calculation
[Scientific / technical calculation by Python] Numerical solution of one-dimensional and two-dimensional wave equations by FTCS method (explicit method), hyperbolic partial differential equations
[Scientific / technical calculation by Python] Logarithmic graph, visualization, matplotlib
[Scientific / technical calculation by Python] Polar graph, visualization, matplotlib
[Scientific / technical calculation by Python] List of matrices that appear in Hinpan in numerical linear algebra
[Scientific / technical calculation by Python] List of usage of (special) functions used in physics by using scipy
[Scientific / technical calculation by Python] Numerical solution of eigenvalue problem of matrix by power method, numerical linear algebra
[Scientific / technical calculation by Python] 3rd order spline interpolation, scipy
Peeping in Python, not scary quantum mechanics 1: Infinite particle potential
[Numerical calculation method, python] Solving ordinary differential equations by Eular method
[Scientific / technical calculation by Python] Basic operation of arrays, numpy
[Scientific / technical calculation by Python] Monte Carlo simulation of thermodynamics of 2D Ising spin system by Metropolis method
[Scientific / technical calculation by Python] Numerical integration, trapezoidal / Simpson law, numerical calculation, scipy
[Scientific / technical calculation by Python] Vector field visualization example, electrostatic field, matplotlib
[Scientific / technical calculation by Python] 1D-3D discrete fast Fourier transform, scipy
[Scientific / technical calculation by Python] 2D random walk (drunk walk problem), numerical calculation