$$\require{cancel}$$

# 7.2: Discretizing Partial Differential Equations


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.

7.2: Discretizing Partial Differential Equations is shared under a CC BY-SA 4.0 license and was authored, remixed, and/or curated by Y. D. Chong via source content that was edited to conform to the style and standards of the LibreTexts platform; a detailed edit history is available upon request.