Forward Difference Method (Explicit Method)
Let us consider how to numerically calculate the time evolution of the wave function for the Schrödinger equation:
\[ \small i\hbar\frac{\partial \psi(x,t)}{\partial t} = -\frac{\hbar^2}{2m}\frac{\partial^2 \psi(x,t)}{\partial x^2} +V(x)\psi(x,t). \]
This is because it is a hassle to solve the equations one by one, and there are many cases where an analytical solution cannot be found. For the function \(\small \psi(x,t)\), divide time and space into equal intervals and denote them as:
\[ \small \begin{align*} &x_j = 0,\Delta x, 2\Delta x,\cdots, (M-1)\Delta x \\ &t_k = 0,\Delta t, 2\Delta t,\cdots, (N-1)\Delta t. \end{align*} \]
Let us denote it as \(\small \psi_{j,k} = \psi(j\Delta x, k\Delta t)\).
We consider a method of calculating the time evolution by approximating partial differentials by difference, similar to ordinary differential equations. When actually calculating the differential approximation, we can put
\[ \small \begin{align*} &\frac{\partial \psi_{j,k}}{\partial t} = \frac{\psi_{j,k+1}-\psi_{j,k}}{\Delta t} \\ &\frac{\partial^2 \psi_{j,k}}{\partial x^2} = \frac{\psi_{j+1,k}+\psi_{j-1,k}-2\psi_{j,k}}{\Delta x^2}. \end{align*} \]
Let us use this to rewrite the partial differential equation into a difference equation. Since
\[ \small i\hbar\frac{\psi_{j,k+1}-\psi_{j,k}}{\Delta t} = -\frac{\hbar^2}{2m}\frac{\psi_{j+1,k}+\psi_{j-1,k}-2\psi_{j,k}}{\Delta x^2} +V(x_{j})\psi_{j,k} \]
holds, by rearranging the equation, it can be expressed as:
\[ \small \psi_{j,k+1} = \psi_{j,k}+i\Delta t\left(\frac{\hbar}{2m}\frac{\psi_{j+1,k}+\psi_{j-1,k}-2\psi_{j,k}}{\Delta x^2}-\frac{V(x_{j})}{\hbar}\psi_{j,k}\right). \]
By performing calculations in order from the previous point in time, the solution to the partial differential equation can be found. This can be thought of as a method equivalent to Euler’s method for ordinary differential equations. It may also be possible to perform calculations in multiple stages, such as the Runge-Kutta method. In the numerical solution of partial differential equations, this type of calculation method is called an explicit method. This method is not often used in practice because it has poorer convergence than the Crank-Nicolson method, which will be explained in the next section, but it can easily calculate even very complex partial differential equations and is a useful method for quickly getting an idea of the solution.
Crank-Nicholson Method
In the previous section, we solved the partial differential equation using forward differences with respect to time, but it is easy to guess that this will result in a bias with respect to time. Of course, such biases will decrease as the grid gets finer, but if it is possible to formulate the equation in a form that is free of bias, it would be better to use that.
Let us consider two approximations, a forward difference approximation and a backward difference approximation, and take the average. Let us assume that the time change range is \(\small \Delta t/2\), and consider two equations: one approximated by forward differences from the previous grid, and one approximated by backward differences from the next grid. In particular, we calculate
\[ \small \begin{align*} &i\hbar\frac{\psi_{j,k+1/2}-\psi_{j,k}}{\Delta t/2} = -\frac{\hbar^2}{2m}\left(\frac{\psi_{j+1,k}+\psi_{j-1,k}-2\psi_{j,k}}{\Delta x^2} \right) +V(x_{j})\psi_{j,k} \\ &i\hbar\frac{\psi_{j,k+1}-\psi_{j,k+1/2}}{\Delta t/2} = -\frac{\hbar^2}{2m}\left(\frac{\psi_{j+1,k+1}+\psi_{j-1,k+1}-2\psi_{j,k+1}}{\Delta x^2} \right) +V(x_{j})\psi_{j,k+1}. \end{align*} \]
Adding the two equations together gives us
\[ \small \begin{align*} &i\hbar\frac{\psi_{j,k+1}-\psi_{j,k}}{\Delta t} = -\frac{\hbar^2}{4m}\Biggl(\frac{\psi_{j+1,k+1}+\psi_{j-1,k+1}-2\psi_{j,k+1}}{\Delta x^2}+\frac{\psi_{j+1,k}+\psi_{j-1,k}-2\psi_{j,k}}{\Delta x^2} \Biggl) \\ &\qquad\qquad\qquad+V(x_{j})\frac{\psi_{j,k}+\psi_{j,k+1}}{2}. \end{align*} \]
This method of calculating derivatives with respect to coordinates is called the Crank-Nicholson method. Methods that use backward differences are called implicit methods in contrast to explicit methods, and the Crank-Nicholson method is a type of implicit method. By rearranging the equation,
\[ \small \begin{align*} &a_j\psi_{j-1,k+1}+b_j\psi_{j,k+1}+c_j\psi_{j+1,k+1} = -a_j\psi_{j-1,k}-d_j\psi_{j,k}-c_j\psi_{j+1,k}\\ &a_j = -\frac{\hbar\Delta t}{4m\Delta x^2}i\\ &b_j = 1+\frac{\hbar\Delta t}{2m\Delta x^2}i+i\frac{V(x_j)}{2\hbar}\Delta t \\ &c_j = -\frac{\hbar\Delta t}{4m\Delta x^2}i\\ &d_j = -1+\frac{\hbar\Delta t}{2m\Delta x^2}i+i\frac{V(x_j)}{2\hbar}\Delta t \end{align*} \]
can be obtained. Expressed in matrix form, this is:
\[ \small \left[ \begin{array}{cccccc} b_{1} & c_{1} & 0 & 0 & \cdots & 0 & 0 \\ a_{2} & b_{2} & c_{2} & 0 & \cdots & 0 & 0 \\ 0 & a_{3} & b_{3} & c_{3} & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \ddots & b_{M-3} & c_{M-3} \\ 0 & 0 & 0 & 0 & \cdots & a_{M-2} & b_{M-2} \end{array} \right] \left[ \begin{array}{c} \psi_{1,k+1} \\ \psi_{2,k+1} \\ \psi_{3,k+1} \\ \vdots \\ \psi_{M-3,k+1} \\ \psi_{M-2,k+1} \end{array} \right]=\qquad\qquad\qquad\qquad \\
\small \left[ \begin{array}{l} -a_{1} \psi_{0,k} – d_{1} \psi_{1,k} – c_{1} \psi_{2,k} – a_{1} \psi_{0, k+1} \\ -a_{2} \psi_{1,k} – d_{2} \psi_{2,k} – c_{2} \psi_{3,k} \\ -a_{3} \psi_{2,k} – d_{3} \psi_{3,k} – c_{3} \psi_{4,k} \\ \qquad\qquad\qquad\vdots \\ -a_{M-3} \psi_{M-4,k} – d_{M-3} \psi_{M-3,k} – c_{M-3} \psi_{M-2,k} \\ -a_{M-2} \psi_{M-3,k} – d_{M-2} \psi_{M-2,k} – c_{M-2} \psi_{M-1,k} – c_{M-2} \psi_{M-1, k+1} \end{array} \right]. \]
If the boundary conditions \(\small \psi_{j,0},\psi_{0,k},\psi_{M-1,k},\forall j,k\) are given, the time evolution of the wave function can be calculated by repeatedly solving this simultaneous equation starting from the front of the time grid. The above is the numerical solution method for partial differential equations known as the finite difference method. When solving the Schrödinger equation numerically, it is important to note that the boundary conditions and the values of the wave function on the grid must be calculated as complex numbers.
Free Particle
As an actual example, we show how the free-state Schrödinger equation is solved using the finite difference method. This is an exact solution even in analytical formulas, so it seems that solving it using finite difference method would also yield the same result as imagined. I tried making it into a video format, so I think it will be easier to get an idea of it.
One-dimensional Coulomb Potential
In the previous article:
I also tried to solve the Schrödinger equation for one-dimensional Coulomb potential, for which I presented a dubious approximate solution, using the finite difference method. Because the wave function mixes between the regions of \(\small x>0\) and \(\small x<0\), interference fringes arise, resulting in chaotic motion. I would like to say that my prediction was correct to some extent in the sense that the movement goes back and forth between the two extremes, but even so, the result was unexpected. I made a video of the time evolution, but it ended up looking kind of silly…
Comments