[Scientific / technical calculation by Python] Solving the Schrodinger equation in the steady state in the 3D isotropic harmonic oscillator potential by the matrix method, boundary value problem, quantum mechanics

Introduction

** 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

[Science and technology calculation by Python] Solving the boundary value problem of ordinary differential equations in matrix format, numerical calculation [1]

I would appreciate it if you could refer to.

Contents

3D isotropic harmonic oscillator potential V(r) = \frac{m_e \omega^2r^2}{2} {\tag {1}} The steady-state Schrodinger equation for the radial wavefunction $ R (r) $ with respect to

u(r) = r R(r) {\tag{2}}

As

-u''(r)+\frac{2 m_e}{\hbar^2}(V(r)+\frac{L(L+1)\hbar^2}{2m_e r^2}) = E \ u(r) {\tag {3}} Can be given.

** 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 **

E_n = \frac{\hbar \omega}{2} (n+\frac{3}{2}),\ ( \ n = 4N+2L+3 {\tag{4}})

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 E = \hbar \omega [1.5, 2.5, 3.5, ...] {\tag{5}})

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

code

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

result

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.

Addendum

** If the functional form of the potential changes, simply change the formula for v [i]. ** **

References

[1] [Science and technology calculation by Python] Solving the boundary value problem of ordinary differential equations in matrix format, numerical calculation

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