Finite Difference Solution of the Schrodinger Equation

Introduction

The time-independent Schrodinger equation is a linear ordinary differential equation that describes the wavefunction or state function of a quantum-mechanical system. The solution of the Schrodinger equation yields quantized energy levels as well as wavefunctions of a given quantum system. The Schrodinger equation can be used to model the behavior of elementary particles, and atoms. When combined with the superposition principle, the Schrodinger equation accounts for the formation of chemical bonds, and hence can be used to model molecular systems and periodic systems such as crystalline materials.

In this article, we focus on the simple one-dimensional Schrodinger equation. We show how the equation can be solved using the finite difference method. We shall illustrate our example using the quantum harmonic oscillator. By changing the potential energy of the system, the code can also be applied to other fundamental one-dimensional models such as square well systems, and the hydrogen atom.

General Formalism

Finite Difference Implementation in Python

import necessary libraries

import numpy as np
import matplotlib.pyplot as plt

define potential energy function

def Vpot(x):
    return x**2

enter computational parameters

a = float(input('enter lower limit of the domain: '))
b = float(input('enter upper limit of the domain: '))
N = int(input('enter number of grid points: '))

define x grid and step size

x = np.linspace(a,b,N)
h = x[1]-x[0]

create kinetic energy matrix

T = np.zeros((N-2)**2).reshape(N-2,N-2)
for i in range(N-2):
    for j in range(N-2):
        if i==j:
            T[i,j]= -2
        elif np.abs(i-j)==1:
            T[i,j]=1
        else:
            T[i,j]=0

create potential energy matrix

V = np.zeros((N-2)**2).reshape(N-2,N-2)
for i in range(N-2):
    for j in range(N-2):
        if i==j:
            V[i,j]= Vpot(x[i+1])
        else:
            V[i,j]=0

create hamiltonian matrix

H = -T/(2*h**2) + V

find eigenvalues and eigenvectors, then sort them in ascending order

val,vec=np.linalg.eig(H)
z = np.argsort(val)
z = z[0:4]
energies=(val[z]/val[z][0])
print(energies)

plot wavefunctions for first 4 lowest states

plt.figure(figsize=(12,10))
for i in range(len(z)):
    y = []
    y = np.append(y,vec[:,z[i]])
    y = np.append(y,0)
    y = np.insert(y,0,0)
    plt.plot(x,y,lw=3, label="{} ".format(i))
    plt.xlabel('x', size=14)
    plt.ylabel('$\psi$(x)',size=14)
plt.legend()
plt.title('normalized wavefunctions for a harmonic oscillator using finite difference method',size=14)
plt.show()

Sample Output for the Quantum Harmonic Oscillator

Using a = -6, b = 6, N = 1001, we obtain the following:

Numerical and Exact energies for first 4 lowest states.

In summary, we’ve shown that the finite difference scheme is a very useful method for solving an eigenvalue equation such as the Schrodinger equation. We illustrated our implementation using the harmonic oscillator system. By changing the potential energy of the system, the code can also be applied to other fundamental one-dimensional models such as square well systems, and the hydrogen atom.