Finite Difference Method #

The finite difference method approximates differential equations with difference equations for numerical solutions.

A web tool implementing this method is available at the link below:

Solving the 1D Time-Independent Schrödinger Equation Using the Finite Difference Method

1. Schrödinger Equation #

Consider the Schrödinger equation for a quantum system in atomic units ($m=1, \hbar=1$):

$$ \begin{align}\label{eq1} \left[-\frac{1}{2}\frac{d^2}{dx^2}+V(x)\right]\psi(x)=E\psi(x) \end{align} $$

ここで、$V(x)$はポテンシャル、$E$はエネルギー固有値を表します。また、$x$は区間$[x_a, x_b]$で定義され、境界条件として固定端条件を考えます。つまり、

where ( V(x) ) is the potential, ( E ) is the energy eigenvalue, and ( x ) is defined in the interval ([x_a, x_b]) with fixed boundary conditions:

$$ \begin{align} \psi(x_a)=\psi(x_b)=0 \end{align} $$

が波動関数に課されている場合を考えます。

1.1. Discretization with Finite Differences #

To solve the Eq. \eqref{eq1} numerically, the wavefunction is treated discretely. The spatial coordinate ( x ) is discretized as¹:

$$ \begin{gather} x_n = n\Delta x + x_a,~~~(n=0,1,\cdots,N) \\ \Delta x \equiv \frac{x_b-x_a}{N} \end{gather} $$

Boundary conditions impose $\psi(x_0) = \psi(x_{N+1}) = 0$.

Using discrete points, the second derivative can be approximated as:

$$ \begin{align} \frac{d^2f}{dx^2}\approx\frac{f(x+\Delta x)-2f(x)+f(x-\Delta x)}{\Delta^2} + O(\Delta^2) \end{align} $$

Substituting this into the Schrödinger equation \eqref{eq1} yields:

$$ \begin{align} -\frac{1}{2\Delta x^2}\psi(x_{n+1})+\left(\frac{1}{\Delta x^2}+V(x_n)\right)\psi(x_n)-\frac{1}{2\Delta x^2}\psi(x_{n-1})=E\psi(x_n) \end{align} $$

For each discrete point ( n ), this can be rewritten as a matrix equation with $\psi_n\equiv\psi(x_n)$:

$$ \begin{gather} \psi_0 = \psi_{N} = 0\\ -\frac{1}{2\Delta x^2}\psi_{n+1}+\left[\frac{1}{\Delta x^2}+V(x_n)\right]\psi_n-\frac{1}{2\Delta x^2}\psi_{n-1}=E\psi_n,~~~(n=1,2,\cdots,N-1) \label{eq2} \end{gather} $$

$$ \begin{align} \begin{pmatrix} d_1 & e_2 & 0 & \cdots & 0 \\ e_2 & d_2 & e_3 & \cdots & 0 \\ 0 & e_3 & d_3 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & e_{N-1} \\ 0 & 0 & 0 & e_{N-1} & d_{N-1} \end{pmatrix} \begin{pmatrix} \psi_1 \\ \psi_2 \\ \psi_3 \\ \vdots \\ \psi_{N-1} \end{pmatrix} =E \begin{pmatrix} \psi_1 \\ \psi_2 \\ \psi_3 \\ \vdots \\ \psi_{N-1} \end{pmatrix}\label{eq3} \end{align} $$

Here,

$$ \begin{align} d_i &\equiv \frac{1}{\Delta x^2}+V(x_i), &&(i=1,2,\cdots, N-1) \\ e_i &\equiv -\frac{1}{2\Delta x^2}, &&(i=2,3,\cdots, N-1). \end{align} $$

Diagonalizing this matrix gives the energy eigenvalues and eigenfunctions.

1.2. Interpretation via Basis Functions #

The method can also be interpreted using basis functions. Representing the wavefunction as a linear combination of basis functions with the coefficients $c_i$:

$$ \begin{align}\label{eq4} \psi(x)=\sum_i c_i \varphi_i(x) \end{align} $$

Additionally, the eigenfunctions $\varphi_i(x)$ are defined for $i=0, 1, \cdots, N$ as follows:

$$ \begin{eqnarray} \varphi_i(x)= \left\{ \begin{aligned} & \frac{1}{\sqrt{\Delta x}} &&(x_i-\Delta x/2 \lt x \lt x_i+\Delta x/2) \\ & \frac{1}{2\sqrt{\Delta x}} &&(x = x_i\pm \Delta x/2) \\ & 0 &&(\text{otherwise}) \end{aligned} \right. \end{eqnarray} $$

Furthermore, integrals involving $\varphi_i(x)$ can be calculated as follows:

$$ \begin{align} \int_{x_a}^{x_b} f(x)\varphi_i(x)dx = \frac{1}{\sqrt{\Delta x}}\int_{x_i-\Delta x/2}^{x_i+\Delta x/2} f(x)dx \end{align} $$

Here, it is assumed that $f(x)$ is an arbitrary function.

1.2.1. Inner Product #

To compute the matrix elements, we first calculate the inner product for an arbitrary function $f(x)$. As a result,

$$ \begin{align} \int_{x_a}^{x_b} \varphi_i(x)f(x)\varphi_j(x)dx &=\frac{1}{\sqrt{\Delta x}}\int_{x_i-\Delta x/2}^{x_i+\Delta x/2} f(x)\varphi_j(x)dx \nonumber\\ &=\frac{1}{\Delta x}\int_{x_i-\Delta x/2}^{x_i+\Delta x/2} f(x)dx \delta_{i,j}\nonumber\\ &\approx \frac{1}{\Delta x}\int_{x_i-\Delta x/2}^{x_i+\Delta x/2} f(x_i)dx \delta_{i,j}\nonumber\\ &= f_i\delta_{i,j} \label{eq5} \end{align} $$

In particular, considering the case where $f(x) = 1$,

$$ \begin{align} \int_{x_a}^{x_b} \varphi_i(x)\varphi_j(x)dx = \delta_{i,j} \end{align} $$

it follows that the basis functions $\varphi_i(x)$ form an orthonormal basis.

The inner product involving the differential operator is given as follows:

$$ \begin{align} \int_{-\infty}^{\infty} \varphi_i(x)\frac{d}{dx}\varphi_j(x)dx &=\frac{1}{\sqrt{\Delta x}}\int_{x_i-\frac{\Delta x}{2}}^{x_i+\frac{\Delta x}{2}} \frac{d}{dx}\varphi_j(x)dx \nonumber\\ &=\frac{1}{\sqrt{\Delta x}}\int_{\varphi(x_i-\frac{\Delta x}{2})}^{\varphi(x_i+\frac{\Delta x}{2})} d\varphi_j(x) \nonumber\\ &=\frac{1}{\sqrt{\Delta x}}\left[\varphi_j(x_i+\frac{\Delta x}{2})-\varphi_j(x_i-\frac{\Delta x}{2})\right] \nonumber\\ &=\frac{1}{\sqrt{\Delta x}}\left[\frac{1}{2\sqrt{\Delta x}}\left(\delta_{j,i+1}+\delta_{j,i}-\delta_{j,i-1}-\delta_{j,i}\right)\right] \nonumber\\ &=\frac{1}{2\Delta x}\left(\delta_{j,i+1}-\delta_{j,i-1}\right) \label{eq6} \end{align} $$

The inner product involving the second-order differential operator is given as follows:

$$ \begin{align} \int_{-\infty}^{\infty} \varphi_i(x)\frac{d^2}{dx^2}\varphi_j(x)dx &=\frac{1}{\sqrt{\Delta x}}\biggl[\frac{d\varphi_j}{dx}\biggr|_{x=x_i+\Delta x/2}-\frac{d\varphi_j}{dx}\biggr|_{x=x_i-\Delta x/2}\biggr] \nonumber\\ &=\frac{1}{\Delta x}\Bigl[\delta(x_i-x_j+\Delta x)-2\delta(x_i-x_j)+\delta(x_i-x_j-\Delta x)\Bigr] \nonumber\\ &\to \frac{1}{\Delta x^2}\left(\delta_{j,i+1}-2\delta_{j,i}+\delta_{j,i-1}\right)\label{eq7} \end{align} $$

Here, the following result was used:

$$ \begin{gather} \frac{d\varphi_j}{dx}=\frac{1}{\sqrt{\Delta x}}\Bigl[\delta\left(x-(x_j-\Delta x/2)\right)-\delta(x-(x_j+\Delta x/2))\Bigr] \\ \delta(x-x’)\xrightarrow{離散化}\frac{\delta_{i,j}}{\Delta x} \end{gather} $$

1.2.2. Matrix Elements of the Hamiltonian #

Using the results derived above, we consider the matrix elements of the Hamiltonian $\hat{H} = -\frac{1}{2}\frac{d^2}{dx^2} + V(x)$ for the coefficients $c_n$ in accordance with the variational principle from Equation \eqref{eq4}. By applying Equations \eqref{eq5} and \eqref{eq7}, we obtain:

$$ \begin{align} H_{ij}\equiv\langle\varphi_i|\hat{H}|\varphi_j\rangle &= -\frac{1}{2}\langle\varphi_i|\frac{d^2}{dx^2}|\varphi_j\rangle + \langle\varphi_i|V(x)|\varphi_j\rangle \nonumber \\ &= -\frac{1}{2\Delta x^2}\left(\delta_{j,i+1}-2\delta_{j,i}+\delta_{j,i-1}\right) + V(x_i)\delta_{i,j} \end{align} $$

Thus,

$$ \begin{align} \begin{pmatrix} d_1 & e_2 & 0 & \cdots & 0 \\ e_2 & d_2 & e_3 & \cdots & 0 \\ 0 & e_3 & d_3 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & e_{N-1} \\ 0 & 0 & 0 & e_{N-1} & d_{N-1} \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_{N-1} \end{pmatrix} =E \begin{pmatrix} c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_{N-1} \end{pmatrix}\label{eq8} \end{align} $$

This result is identical to Equation \eqref{eq3}.

2. Numerical Calculation #

The implementation for solving the problem using the finite difference method with JavaScript and HTML has been created on the following page. Users can specify the potential and execute the computation:

Solving the 1D Time-Independent Schrödinger Equation via the Finite Difference Method

The appearance of the page is as shown below:

Due to either JavaScript’s limited speed or inefficient implementation, you may notice slower computation times when the matrix size reaches around 600x600 elements.