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) **
** 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>
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
E1=2.46740110027234 (Ry) E2 = 9.869604401089358 (Ry) E3 = 22.20660990245106 (Ry) ... [1].
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()
Figure. For even functions. The zero point is the eigenvalue.
Figure. For odd functions.
Figure.Ground state wavefunction
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.
One-dimensional Schrodinger equation
Is ** in bound state **
** 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" **.
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.
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.
[1] Qiita article [Scientific and technical calculation by Python] Solving second-order ordinary differential equations by Numerov method, numerical calculation
Recommended Posts