1D Schrödinger Equation: Finite Difference Method

Recommended resources

The freely available resources this lesson will direct you to.

Hello! Welcome to the first lesson of our new module on numerical methods in quantum mechanics.

In our previous module, we explored fascinating collective phenomena like the Fractional Quantum Hall Effect, where the breakdown of analytical methods for strongly-correlated systems necessitated brilliant intuitive leaps like the Laughlin wavefunction. This highlights a crucial theme in modern physics: many, if not most, interesting problems cannot be solved with pen and paper alone. This is where computational physics becomes an indispensable tool.

Your background in writing simulations makes this module a natural fit. We will begin building a practical toolkit for solving quantum problems numerically. Our goal for today is to:

Implement a finite difference method to find the bound-state energies of the 1D time-independent Schrödinger equation.

We will transform the familiar differential equation into a matrix problem that a computer can solve, allowing us to find the quantum states for virtually any one-dimensional potential you can imagine.

1. The Schrödinger Equation as a Matrix Problem

The time-independent Schrödinger equation (TISE) is the cornerstone for finding stationary states:

H^ψ (x) = E ψ (x)

where the Hamiltonian for a particle in one dimension is:

H^= - \frac{ℏ ^{2}}{2 m} \frac{d ^{2}}{d x ^{2}} + V (x)

For potentials like the infinite square well or the harmonic oscillator, we can solve this differential equation analytically. But what about a more complex potential, like an electron in a disordered lattice or the double-well potential we'll see later? The analytical path is often blocked.

The core idea of the finite difference method is to discretize the problem. Instead of trying to find the continuous function $ψ (x)$ , we will find its value at a series of discrete points $x_{i}$ separated by a small step $Δ x$ . We represent our wavefunction as a vector:

ψ = ψ (x_{1}) ψ (x_{2}) ⋮ ψ (x_{N})

To discretize the TISE, we need to approximate the second derivative. We can do this using the central difference formula:

Finite Difference Approximation for Quantum Mechanics — This diagram illustrates how the second derivative of the wavefunction at a point $i$, $\psi_i''$, can be approximated using the values at its neighboring points, $\psi_{i-1}$ and $\psi_{i+1}$. The distance between points is the grid spacing, $\Delta x$.

The formula for the second derivative at point $x_{i}$ is:

\frac{d ^{2} ψ}{d x ^{2}}_{x_{i}} \approx \frac{ψ ( x _{i + 1} ) - 2 ψ ( x _{i} ) + ψ ( x _{i - 1} )}{( Δ x ) ^{2}} = \frac{ψ _{i + 1} - 2 ψ _{i} + ψ _{i - 1}}{( Δ x ) ^{2}}

By substituting this into the TISE for each point $x_{i}$ , we convert the single differential equation into a system of $N$ coupled algebraic equations:

- \frac{ℏ ^{2}}{2 m} (\frac{ψ _{i + 1} - 2 ψ _{i} + ψ _{i - 1}}{( Δ x ) ^{2}}) + V (x_{i}) ψ_{i} = E ψ_{i}

Rearranging this gives:

(- \frac{ℏ ^{2}}{2 m ( Δ x ) ^{2}}) ψ_{i - 1} + (\frac{ℏ ^{2}}{m ( Δ x ) ^{2}} + V (x_{i})) ψ_{i} + (- \frac{ℏ ^{2}}{2 m ( Δ x ) ^{2}}) ψ_{i + 1} = E ψ_{i}

This equation has a very specific structure: the value $ψ_{i}$ is related only to its nearest neighbors, $ψ_{i - 1}$ and $ψ_{i + 1}$ . When we write the entire system of equations for all $i = 1, \dots, N$ , it takes the form of a matrix eigenvalue equation:

H ψ = E ψ

Here, $ψ$ is the vector of wavefunction values, and $H$ is the Hamiltonian, now represented as a matrix.

2. Constructing the Hamiltonian Matrix

The equation above defines the elements of our Hamiltonian matrix $H$ .

Diagonal elements ( $H_{ii}$ ): These correspond to the terms multiplying $ψ_{i}$ .
$H_{ii} = \frac{ℏ ^{2}}{m ( Δ x ) ^{2}} + V (x_{i})$
Off-diagonal elements ( $H_{i, i \pm 1}$ ): These correspond to the terms multiplying the neighbors $ψ_{i - 1}$ and $ψ_{i + 1}$ .
$H_{i, i - 1} = H_{i, i + 1} = - \frac{ℏ ^{2}}{2 m ( Δ x ) ^{2}}$
All other elements are zero.

This structure results in a tridiagonal matrix, which is a sparse matrix with non-zero elements only on the main diagonal, the super-diagonal (above the main one), and the sub-diagonal (below the main one).

For bound-state problems, we solve the equation on a finite grid $[x_{min}, x_{ma x}]$ and assume the wavefunction vanishes at the boundaries: $ψ (x_{min}) = ψ (x_{ma x}) = 0$ . This "Dirichlet boundary condition" is naturally incorporated into our matrix equation if we construct the matrix for the interior points only.

The problem is now perfectly set up for a numerical solution: construct the matrix $H$ and find its eigenvalues and eigenvectors. The eigenvalues will give us the quantized energy levels $E$ , and the corresponding eigenvectors will give us the shape of the wavefunctions $ψ (x)$ .

3. Python Implementation

Let's see this in action. The following video provides an excellent walkthrough of the entire process, from the initial equation to a Python implementation that can handle any potential.

Eigenstates of ANY 1D Potential in PYTHON

This video by Mr. P Solver demonstrates how to translate the theory we just discussed into a working Python script. It covers the discretization, the matrix formulation, and the use of SciPy's efficient solvers.

Please watch the following segments: Discretization & Dimensionless Form (03:20 - 06:23): This covers the finite difference approximation, similar to our discussion. The video also introduces the useful practice of making the equation dimensionless, which simplifies the constants involved. The Matrix Eigenvalue Problem (06:23 - 09:45): This part visualizes how the system of equations becomes a tridiagonal matrix eigenvalue problem. Python Code: Setup and Matrix Construction (09:45 - 014:59): Follow how the spatial grid, potential, and the diagonals of the Hamiltonian matrix are defined in NumPy. Python Code: Solving and Plotting (14:59 - 17:28): See how scipy.linalg.eigh_tridiagonal is used to find the energies and wavefunctions, and how to plot the results.

The video brilliantly demonstrates the power of this method. Once the framework is built, swapping out the harmonic oscillator for a more complex potential—like the one with a Gaussian bump—is trivial and yields fascinating physical results that would be very difficult to obtain analytically.

The key steps of the implementation are:

Define Constants and Grid: Set up physical constants ( $m, ℏ$ ) and the spatial grid using np.linspace. Calculate the grid spacing $Δ x$ .
Define Potential: Create a Python function V(x) that returns the potential energy at any point x. This makes the code general.
Construct Matrix Diagonals: Create two NumPy arrays: one for the main diagonal (d) and one for the off-diagonals (e).
- d_i = \hbar^2 / (m * dx^2) + V(x_i)
- e_i = -\hbar^2 / (2 * m * dx^2)
Solve the Eigenvalue Problem: Use a numerical library to find the eigenvalues and eigenvectors. scipy.linalg.eigh_tridiagonal(d, e) is highly optimized for this specific type of matrix.
Analyze Results: The returned eigenvalues are the bound-state energies. The eigenvectors are the discretized wavefunctions, which can be plotted. Remember to normalize them!

4. Alternative Approaches and Advanced Considerations

While the method in the video is excellent, it's useful to see other ways to construct and solve the matrix problem, especially as you might encounter them in different contexts.

Quantum Mechanics with Python

This blog post, 'Quantum Mechanics with Python', provides a nice textual summary of the method and introduces an important concept for large-scale computations: sparse matrices.

Please read the sections "Discretizing the 1D Schrödinger Equation" and "Solving the Eigenvalue Problem Numerically". Focus on how the Hamiltonian is constructed and the motivation for using sparse matrix libraries.

As the article mentions, when your grid size N becomes very large (e.g., in 2D or 3D problems, or for very high precision), storing the full $N \times N$ Hamiltonian matrix in memory becomes inefficient or impossible. Since our matrix is mostly zeros, we can use sparse matrix formats that only store the non-zero elements. Libraries like scipy.sparse are designed for this. The function scipy.sparse.linalg.eigs is a sparse matrix eigenvalue solver, which is essential for more advanced computational work, including the exact diagonalization you'll encounter later.

For another perspective on building the matrix, the "Dot Physics" video (LINK) shows how to construct the full matrix using np.diag, which can be more explicit if you want to visualize the matrix structure itself before optimizing with a tridiagonal or sparse solver.

Test your understanding! Apply your knowledge.

Now it's your turn to get your hands dirty. Your task is to implement the finite difference solver yourself. You can use the code from the "Mr. P Solver" video as a starting point.

Problem: Solve the TISE for a double-well potential, given by:
$V (x) = V_{0} [(\frac{x}{a})^{4} - 2 (\frac{x}{a})^{2} + 1] = V_{0} ((\frac{x}{a})^{2} - 1)^{2}$

This potential has two minima at $x = \pm a$ , separated by a barrier of height $V_{0}$ .

Instructions:

Write a Python script to find and plot the first 4-5 energy eigenvalues and their corresponding wavefunctions ( $∣ ψ (x) ∣^{2}$ ).
Choose sensible parameters. For example, you could work in atomic units ( $ℏ = 1, m = 1$ ) and try $a = 2$ and $V_{0} = 10$ . You'll need to choose a grid that is large enough to contain the wavefunctions (e.g., from $x = - 5$ to $x = 5$ ).
Observe the first two eigenstates (the ground state and the first excited state). What do you notice about their shapes and their energies?
Experiment by changing the barrier height $V_{0}$ . How does the energy difference between the first two states change as you increase or decrease $V_{0}$ ? What does this imply about quantum tunneling?

Show hints on what to expect

You should find that the ground state $ψ_{0}$ is an even function (symmetric) with probability density peaked in both wells. The first excited state $ψ_{1}$ should be an odd function (antisymmetric) with a node at $x = 0$ .

These two states will be very close in energy, forming a "doublet". This is because they correspond to the symmetric and antisymmetric combinations of a particle being in the left well or the right well. As you increase the barrier height $V_{0}$ , the tunneling between the wells becomes less likely. This should cause the energy splitting $Δ E = E_{1} - E_{0}$ to decrease dramatically, as the two states become nearly degenerate.

Conclusion

In this lesson, we have built a powerful and general tool. By discretizing space and approximating the kinetic energy operator, we transformed the time-independent Schrödinger equation from a challenging differential equation into a standard matrix eigenvalue problem.

Key Takeaways:

The finite difference method approximates derivatives on a discrete grid, turning differential equations into algebraic ones.
The 1D TISE becomes a matrix eigenvalue problem $H ψ = E ψ$ .
For a 1D problem with nearest-neighbor interactions, the Hamiltonian matrix $H$ is tridiagonal.
Numerical libraries like NumPy and SciPy provide efficient functions (eigh_tridiagonal, eigs) to solve this problem for the energy eigenvalues and wavefunction eigenvectors.
This method is versatile and can be applied to almost any arbitrary potential $V (x)$ .

Preview of the Next Lesson:

We've successfully tackled the time-independent Schrödinger equation. The natural next step is to consider dynamics. In our next lesson, you will implement the Crank-Nicolson method to simulate wavepacket evolution for the 1D time-dependent Schrödinger equation. This will allow us to watch quantum systems evolve in time.