** The steady-state Schrodinger equation can be reduced to the (general) eigenvalue problem of a matrix by the finite difference method. The purpose of this paper is to determine the eigenenergy and eigenfunctions of electrons in the three-dimensional isotropic harmonic oscillator potential using this method. ** **
About the outline of the matrix method
I would appreciate it if you could refer to.
3D isotropic harmonic oscillator potential
As
** This equation can be transformed into a confluent supergeometric differential equation, and the exact solution of the energy eigenvalue $ E $ of a regular solution at the origin is **
It is known to be ** [2]. ** Where L is the orbital quantum number and N is the natural number $ N = 0,1,2, ... $.
Therefore the energy spectrum
Will be. Most of them are degenerate.
In this paper, we solve equation (3) as a matrix eigenvalue problem and calculate the energy eigenvalue $ E_n $ and the eigenfunction $ u_n (r) $.
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 $ 1Ry = 13.6058 eV
Is to be.
"""
Boundary value problem by matrix method
3D isotropic harmonic oscillator potential
"""
import numpy as np
import scipy.linalg
import matplotlib.pyplot as plt
delta_x = 0.05
x0, x1 = 0.001, 10
N=int((x1-x0)/delta_x)
print("N=",N)
L=0 #Orbital quantum number. Change as appropriate.
hbar=1
m_elec=1/2
omega=1
y = np.zeros([N-1,N+1])
y[:,0] = 0
y[:,-1] = 0
A=np.zeros([N-1,N-1])
B=np.identity(N-1)
v = np.zeros([N-1])
vcent = np.zeros([N-1])
veff = np.zeros([N-1])
for i in range(N-1): #Harmonic oscillator ponshall
x = x0 + i*delta_x
vcent[i] = L*(L+1)/x**2 #Centrifugal force potential
v[i] = m_elec*omega**2*(x**2)/2
veff[i] = v[i] +vcent[i]
for i in range(N-1): #Be careful of the index position because it is a triple diagonal matrix.
if i == 0:
A[i,i] = 2/(delta_x**2) + veff[i]
A[i,i+1] = -1/(delta_x**2)
elif i == N-2:
A[i,i-1] = -1/(delta_x**2)
A[i,i] = 2/(delta_x**2) + veff[i]
else:
A[i,i-1] = -1/(delta_x**2)
A[i,i] = 2/(delta_x**2) + veff[i]
A[i,i+1] = -1/(delta_x**2)
eigen, vec= scipy.linalg.eigh(A,B)
print("eigen values_3points=",eigen[0:4])
for j in range(N-1):
for i in range(1,N):
y[j, i] = vec[i-1,j]
#
# for plot
X= np.linspace(x0,x1, N+1)
plt.plot(X, y[0,:],'-',markersize=5,label='y1')
plt.plot(X, y[1,:],'-',markersize=5,label='y2')
plt.plot(X, y[2,:],'-',markersize=5,label='y3')
plt.legend(loc='upper right')
plt.xlabel('X') #x-axis label
plt.ylabel('Y') #y-axis label
plt.show()
L=0 eigen values_3points= [ 1.46118173 3.44087896 5.42483146]
** Exact values of 1.5, 3.5 and 5.5 are in agreement within a few percent of error. ** **
The function form of $ u_n (r) $ n = 0, 2, 4 is shown in the figure below.
L=1 eigen values_3points= [ 2.4997526 4.49899205 6.49760609]
** Exact values of 2.5, 4.5 and 6.5 are in agreement within a few percent of error. ** **
The function form of $ u_n (r) $ n = 1, 3, 5 is shown in the figure below.
** If the functional form of the potential changes, simply change the formula for v [i]. ** **
[2] Kenichi Goto et al. ["Detailed Theory Applied Quantum Mechanics Exercise"](https://www.amazon.co.jp/%E8%A9%B3%E8%A7%A3%E7%90%86%E8 % AB% 96% E5% BF% 9C% E7% 94% A8% E9% 87% 8F% E5% AD% 90% E5% 8A% 9B% E5% AD% A6% E6% BC% 94% E7% BF % 92-% E5% BE% 8C% E8% 97% A4-% E6% 86% B2% E4% B8% 80 / dp / 4320031717), Kyoritsu Publishing, 1982.
Recommended Posts