Tutorial 19. Interactive -- Solve 1D Schrödinger equation numerically as a finite difference boundary value problem

Learning objectives

Learn the numerical algorithm for solving Schrödinger equation as a 1D boundary value problem
Apply the algorithm to solve 1D Harmonic Oscillator

Background

Earlier we learned that Schrödinger equation can be solved numerically using the numerical techniques for partial differential equations (PDEs) or ordinary differential equations (ODEs). We have learned about the Numerov method in our Interactive tutorial , where we solved the 1D Harmonic Oscillator problem as a finite difference initial value problem: we need to start from one grid point, propagate point by point to obtain the wave function. Obtaining the eigenvalues from the Numerov method is also not very convenient.

In this tutorial, we are going to show that 1D Schrödinger equation can also be solved numerically as a finite difference boundary value problem. As you will see, it is more convenient than the Numerov method. We are going to use the algorithm given in the following paper:

Donald Truhlar, “Finite difference boundary value method for solving one-dimensional eigenvalue equations”, Journal of Computational Physics 10, (1972) 123-132

For a 1D Schrödinger equation in atomic unit ($\hbar=1$)

\[-\frac{1}{2\mu}\frac{d^2}{d R^2}\phi_k(R) + V(R)\phi_k(R) = \epsilon_k \phi_k(R)\]

where $\phi_k(R)$ and $\epsilon_k$ are the $k$-th eigenfunction and eigenvalues of the Hamiltonian. The overall idea of this method is still discretizing the 1D space into $N$ grid points and replacing 2nd order derivative (kinetic energy part) with finite difference. The difference from the Numerov method is that, now we find a set of linear equations the wave function values at the grid points need to satisfy: $\sum_j F_{ij} \phi_{j,k} = \lambda_k^n \phi_{i,k}$

where

\[F_{ij} = \frac{1}{\mu}\left[\delta_{i-1,j}(1-\delta_{i,1})\right] + \delta_{ij}\left[-\frac{2}{\mu}+U^h_{ii}\right] + \frac{1}{\mu}\left[\delta_{i+1,j}(1-\delta_{i,N})\right]\]

$U^h_{ii} = -2h^2 V(R_i)$ $\lambda_k^h = -2h^2 \epsilon_k$

Here, $h$ is the discretization size, $i, j$ denote the indices of the grid points.

Then the problem is transformed into a linear algebra problem. We can just build the F matrix and diagonalize it to obtain the eigenvalues and values of the eigenfunctions of the Hamiltonian simultaneously at all grid points. No point-wise propagation anymore.

Try and learn

We will demonstrate how we can solve the 1D harmonic oscillator problem using the finite difference boundary value techniques. For simplicity, we use atomic units. Then $\hbar=1$. Mass unit is electron mass $m_e$=1. The length unit is bohr. The mass of the harmonic oscillator $\mu=0.5$, and angular frequency $\omega=2$.

Run the following code to explore how the grid size can influence the result

Instructions

First, click “Activate” to activate the code block.
Once you see the buttons to “run” at the bottom left corner of the code block, click “run” to run the code.
Please be patient. Starting the kernel can be slow sometimes.
Select the quantum number $v$ to want to visualize the solution to the wave function of that state. The energy of that state will also be printed out in the title of the plot.
Drag the slider bar to adjust the grid size (how many discrete points to use in the range [-10,10])


%matplotlib widget
import ipywidgets as widgets
import numpy as np
import matplotlib.pyplot as plt
import time

hbar = 1
class PDE_BV:
    def __init__(self,m,omega, L, npoints=1000):
        self.m = m
        self.omega = omega
        self.L = L
        self.xlower = -L*0.5
        self.xupper = L*0.5
        self.npoints = npoints
        self.x = np.linspace(self.xlower,self.xupper,self.npoints)
        self.h = self.x[1]-self.x[0]
        self.psi = [None]*self.npoints
        self.E = np.zeros(self.npoints)

    def getV(self, x):
        return 0.5*self.m*self.omega**2*x**2

    def getUii(self, i):
        return  -2*(self.h**2)*self.getV(self.x[i])

    def getF(self):
        F = np.zeros([self.npoints,self.npoints])
        for i in range(0,self.npoints):
            F[i,i] = self.getUii(i) - 2/self.m
            if i > 0:
                F[i,i-1] = 1/self.m
            if i < self.npoints-1:
                F[i,i+1] = 1/self.m

        return F

    def diagonalize(self):
        F = self.getF()
        w,v = np.linalg.eig(F)
        return w,v


    def calc(self):
        w,v = self.diagonalize()
        self.E = - w/2.0/self.h**2
        #print(w)
        #print(self.E[0])
        for k in range(0,len(w)):
            self.psi[k] = v[:,k]
            integral = self.h*np.dot(self.psi[k],self.psi[k])
            self.psi[k] = self.psi[k]/integral**0.5



    def plotWFN(self, v):
        plt.plot(self.x,self.psi[v],label=r'$\psi_v(x)$, k = ' + str(v))
        plt.title(r'$v=$'+ str(v) + r', $E_v$=' + '{:10.4f}'.format(self.E[v]))
        plt.legend()
        plt.xlabel(r'$x$(bohr)')
        plt.ylabel(r'$\psi(x)$')

### Now define the interactive widgets
# Number of discrete points to use in the range
n_options = widgets.IntSlider(
    value=100,
    min=50,
    max=501,
    step=10,
    description=r'npoints:',
    disabled=False,
)

v_options = widgets.IntSlider(
    value=0,
    min=0,
    max=10,
    step=1,
    description=r'state $v$:',
    disabled=False,
)


m1 = 0.5
omega1 = 2
L1 = 20

@widgets.interact(n=n_options, v=v_options)
def update(n,v):
    plt.clf()
    solver = PDE_BV(m1,omega1, L1, npoints=n)
    solver.calc()
    solver.plotWFN(v)

Liu Lab Computational Chemistry Tutorials

Tutorial 19. Interactive -- Solve 1D Schrödinger equation numerically as a finite difference boundary value problem

Learning objectives

Background

Try and learn

Run the following code to explore how the grid size can influence the result

Further reading