7.2: Discretizing Partial Differential Equations
- Page ID
- 34843
With discretized derivatives, differential equations can be formulated as discrete systems of equations. We will discuss this using a specific example: the discretization of the time-independent Schrödinger wave equation in 1D.
7.2.1 Deriving a Finite-Difference Equation
The 1D time-independent Schrödinger wave equation is the second-order ordinary differential equation
\[-\frac{\hbar^2}{2m} \, \frac{d^2\psi}{dx^2} + V(x) \psi(x) = E \psi(x),\]
where \(\hbar\) is Planck's constant divided by \(2\pi\), \(m\) is the mass of the particle, \(V(x)\) is the potential, \(\psi(x)\) is the quantum wavefunction of an energy eigenstate of the particle, and \(E\) is the corresponding energy. The differential equation is usually treated as an eigenproblem, in the sense that we are given \(V(x)\) and seek to find the possible values of the eigenfunction \(\psi(x)\) and the energy eigenvalue \(E\). For convenience, we will adopt units where \(\hbar=m = 1\):
\[-\frac{1}{2}\, \frac{d^2\psi}{dx^2} + V(x) \psi(x) = E \psi(x).\]
To discretize this differential equation, we simply evaluate it at \(x = x_n\):
\[-\frac{1}{2}\, \psi''(x_n) + V_n \psi_n = E \psi_n,\]
where, for conciseness, we denote
\[V_n \equiv V(x_n).\]
We then replace the second derivative \(\psi''(x_n)\) with a discrete approximation, specifically the three-point rule:
\[-\frac{1}{2h^2}\, \Big[\psi_{n+1} - 2\psi_n + \psi_{n-1} \Big] + V_n \psi_n = E \psi_n.\]
This result is called a finite-difference equation, and it would be valid for all \(n\) if the number of discretization points is infinite. However, if there is a finite number of discretization points, \(\{x_0, x_1, \dots, x_{N-1}\}\), then the finite-difference formula fails at the boundary points, \(n=0\) and \(n=N-1\), where it involves the value of the function at the "non-existent" points \(x_{-1}\) and \(x_{N}\). We'll see how to handle this problem in the next section.
Boundaries aside, the finite-difference equation describes a matrix equation:
\[\left\{-\frac{1}{2h^2}\begin{bmatrix} \ddots & \ddots \\ \ddots & -2 & 1\\ & 1 & -2 & 1\\ & & 1 & -2 & \ddots \\ & & & \ddots & \ddots & \end{bmatrix} + \begin{bmatrix} & \ddots & \\ & & V_{n-1} & \\ & & & V_n & \\ & & & & V_{n+1} & \\ & & & & & \ddots \end{bmatrix} \right\}\begin{bmatrix}\vdots \\ \psi_{n-1} \\ \psi_n \\ \psi_{n+1} \\ \vdots\end{bmatrix} = E \begin{bmatrix}\vdots \\ \psi_{n-1} \\ \psi_n \\ \psi_{n+1} \\ \vdots\end{bmatrix}.\]
The second-derivative operator is represented by a tridiagonal matrix with \(-2\) in each diagonal element, and \(1\) in the elements directly above and below the diagonal. The potential operator is represented by a diagonal matrix, where the elements along the diagonal are the values of the potential at each discretization point. In this way, the Schrödinger wave equation is reduced to a discrete eigenvalue problem.
7.2.2 Boundary Conditions
We now have to figure out how to handle the boundaries. Let us suppose \(\psi(x)\) is defined over a finite interval, \(a \le x \le b\). As we recall from the theory of differential equations, the solution to a differential equation is not wholly determined by the differential equation itself, but also by the boundary conditions that are imposed. Thus, we have to specify how \(\psi(x)\) behaves at the end-points of the interval. We will show how this is done for a couple of the most common boundary conditions; other choices of boundary conditions can be handled using the same kind of reasoning.
Dirichlet Boundary Conditions
Under Dirichlet boundary conditions, the wavefunction vanishes at the boundaries:
\[\psi(a) = \psi(b) = 0.\]
Physically, these boundary conditions apply if we let the potential blow up in the external regions, \(x>b\) and \(x<a\), thus forcing the wavefunction to be strictly confined to the interval \(a \le x \le b\).
We have not yet stated how the discretization points \(\{x_0, \dots, x_{N-1}\}\) are distributed within the interval; we will make this decision in tandem with the implementation of the boundary conditions. Consider the first discretization point, \(x_{0}\), wherever it is. The finite-difference equation at this point is
\[-\frac{1}{2h^2}\, \Big[\psi_{-1} - 2\psi_0 + \psi_{1} \Big] + V_0 \psi_0 = E \psi_0.\]
This involves the wavefunction at \(x_{-1}\), which lies just outside our set of discretization points. But if we choose the discretization points so that \(x_{-1}=a\), then \(\psi_{-1}=0\) under Dirichlet boundary conditions, so the above finite-difference formula reduces to
\[-\frac{1}{2h^2}\, \Big[- 2\psi_0 + \psi_{1} \Big] + V_0 \psi_0 = E \psi_0.\]
As for the other boundary, the finite-difference equation at \(x_{N-1}\) involves \(\psi_{N}\). If we choose the discretization points so that \(x_{N}=b\), then the finite-difference formula becomes
\[-\frac{1}{2h^2}\, \Big[ \psi_{N-2} - 2\psi_{N-1} \Big] + V_{N-1} \psi_{N-1} = E \psi_{N-1}.\]
From this, we conclude that the discretization points ought to be equally spaced, with \(x_{0}\) at a distance \(h\) to the right of the left boundary \(a\) and \(x_{N-1}\) a distance \(h\) to the left of the right boundary \(b\). This is shown in the following figure:
Since there are \(N\) discretization points, the interval should contain \((N+1)\) multiples of \(h\). Hence,
\[h = \frac{b - a}{N + 1} \;\; \Rightarrow \;\; x_n \,=\, a + h (n+1) \,=\, \frac{a(N-n)+b(n+1)}{N+1}.\]
Having made the above choices, the matrix equation becomes
\[\left\{-\frac{1}{2h^2}\begin{bmatrix} -2 & 1 \\ 1 & -2 & \ddots \\ & \ddots & \ddots & 1 \\ & & 1 & -2\end{bmatrix} + \begin{bmatrix} V_0 \\ & V_1 \\ & & \ddots\\ & & & V_{N-1} \end{bmatrix} \right\}\begin{bmatrix}\psi_0 \\ \psi_1 \\ \vdots \\ \psi_{N-1}\end{bmatrix} = E \begin{bmatrix}\psi_0 \\ \psi_1 \\ \vdots \\ \psi_{N-1}\end{bmatrix}.\]
You can check for yourself that the first and last rows of this equation are the correct finite-difference equations at the boundary points, corresponding to Dirichlet boundary conditions.
Neuman boundary conditions
Neumann boundary conditions are another common choice of boundary conditions. They state that the first derivatives vanish at the boundaries:
\[\psi'(a) = \psi'(b) = 0.\]
An example of such a boundary condition is encountered in electrostatics, where the first derivative of the electric potential goes to zero at the surface of a charged metallic surface.
We follow the same strategy as before, figuring out the discretization points in tandem with the boundary conditions. Consider again the finite-difference equation at the first discretization point:
\[-\frac{1}{2h^2}\, \Big[\psi_{-1} - 2\psi_0 + \psi_{1} \Big] + V_0 \psi_0 = E \psi_0.\]
To implement the condition that first derivative vanishes at the boundary, we invoke the mid-point rule. Suppose the boundary point \(x=a\) falls in between the points \(x_{-1}\) and \(x_{0}\). Then, according to the mid-point rule,
\[\frac{\psi_0 - \psi_{-1}}{h} \approx \psi'(a) = 0.\]
With this choice, therefore, we can make the replacement \(\psi_{-1} = \psi_0\) in the finite-difference equation, which then becomes
\[-\frac{1}{2h^2}\, \Big[- \psi_0 + \psi_{1} \Big] + V_0 \psi_0 = E \psi_0.\]
Similarly, to apply the Neumann boundary condition at \(x=b\), we let the boundary fall between \(x_{N-1}\) and \(x_{N}\), so that the finite-difference equation becomes
\[-\frac{1}{2h^2}\, \Big[ \psi_{N-2} - \psi_{N-1} \Big] + V_{N-1} \psi_{N-1} = E \psi_{N-1}.\]
The resulting distribution of discretization points is shown in the following figure:
Unlike the Dirichlet case, the interval contains \(N\) multiples of \(h\). Hence, we get a different formula for the positions of the discretization points
\[h = \frac{b - a}{N} \;\; \Rightarrow \;\; x_n \,=\, a + h \left(n+\frac{1}{2}\right) \,=\, \frac{a(N-n-\tfrac{1}{2})+b(n+\tfrac{1}{2})}{N}.\]
The matrix equation is:
\[\left\{-\frac{1}{2h^2}\begin{bmatrix} -1 & 1 \\ 1 & -2 & \ddots \\ & \ddots & \ddots & \ddots \\ && \ddots & -2 & 1 \\ & & & 1 & -1\end{bmatrix} + \begin{bmatrix} V_0 \\ & V_1 \\ & & \ddots\\ & & & V_{N-1} \end{bmatrix} \right\}\begin{bmatrix}\psi_0 \\ \psi_1 \\ \vdots \\ \psi_{N-1}\end{bmatrix} = E \begin{bmatrix}\psi_0 \\ \psi_1 \\ \vdots \\ \psi_{N-1}\end{bmatrix}. \]
Due to the Neumann boundary conditions and the mid-point rule, the tridiagonal matrix has \(-1\) instead of \(-2\) on its corner entries. Again, you can verify that the first and last rows of this matrix equation correspond to the correct finite-difference equations for the boundary points.