+
Using the Finite-Difference Approximation and Hamiltonians to solve 1D Quantum Mechanics Problems

Developed by Marie Lopez del Puerto - Published May 16, 2017

This set of exercises are an introduction to the finite-difference approximation and its use in solving differential equations. The method is implemented in order to write the Schrödinger equation as a matrix (Hamiltonian). The numerical solution to the one-dimensional infinite well is calculated in order to explore issues of convergence by comparing numeric to analytic solutions. Finally, other one-dimensional problems are included, some of which have exact solutions and some of which can only be solved numerically.
Subject Areas Modern Physics and Quantum Mechanics Beyond the First Year MATLAB and IPython/Jupyter Notebook At the end of this exercise set students should: * Be able to explain truncation errors when finding numerical solutions to a differential equation. (Exercises 1-2) * Be able to describe what it means to discretize a differential equation and why it is necessary to do so in order to solve the equation using a computer. (Exercises 1-3) * Be able to explain convergence as a different number of points are used, $N$, and as different solutions are calculated, $n$. (Exercises 4-7) * Have practiced plotting, array manipulation, and data analysis. (Exercises 4-7) * Be able to change the Hamiltonian in an appropriate manner when faced with a different problem. (Exercises 6-8) 210 min

These exercises are not tied to a specific programming language. Example implementations are provided under the Code tab, but the Exercises can be implemented in whatever platform you wish to use (e.g., Excel, Python, MATLAB, etc.).

#Exercise 1# Taylor’s theorem for expanding a function $f(x)$ around a point $x_0$ is given by: $f(x)= \sum_{k=0}^\infty\frac{f^k (x_0)}{k!} (x-x_0)^k$ a) Write down Taylor’s theorem for a function $f(x_0+d)$ around a point $x_0$, $f(x_0+d)= \sum_{k=0}^\infty\frac{f^k (x_0)}{k!} d^k$, up to the third derivative. b) Rearrange the equation to solve for the first derivative, $f ’(x_0)$. The error is given by the terms on the right hand side that involve derivatives since we cannot calculate these without having an analytic expression for the function. Since $d$ is small, the leading term (the largest one) is the one that involves the second derivative: $-\frac{1}{2} d f''(x)$. The error approaches zero at the same rate that $d$ does; we call this error of order $d$, or $O(d)$. This is a truncation error, present even if the calculations are performed with infinite precision. #Exercise 2# a) Write down Taylor’s theorem for $f(x_0-d)$ up to the third derivative. b) Using the expressions you found in 1(a) and 2(a), calculate the first derivative using the two-point finite difference approximation: $f'(x_0)\simeq \frac{f(x_0+d)-f(x_0-d)}{2d}$. c) What is the order of the truncation error for the two-point finite difference approximation? d) One of your friends just showed up (late!). Briefly explain to him/her the advantages of using the finite-difference approximation instead of the definition of derivative to calculate derivatives numerically. #Exercise 3# For the Schrödinger equation, we need the second derivative of the wave function. In this case the finite difference approximation is: $\psi''(x)\simeq \frac{\psi(x+d)-2\psi(x)+\psi(x-d)}{d^2}$ plus a truncation error of order $d^2$. Or using our previous notation: $\psi''(x_j)\simeq \frac{\psi(x_{j-1})-2\psi(x_j)+\psi(x_{j+1})}{d^2}$ Substituting the finite difference approximation for the second derivative of the wave function in the Schrödinger equation, we obtain a set of $N$ equations, one for each point on the discretized wave function: $-\frac{\hbar^2}{2m} \frac{ \psi(x_{j-1} )-2\psi(x_j )+\psi(x_{j+1} )}{d^2} +V_j\psi(x_j)=E \psi(x_j )$ In other words, the discretized wave function must satisfy the Schrödinger equation at each point. a) Consider an infinite square well of length $L$. Inside the well, $V=V_j=0$. The Schrödinger equation becomes $-\frac{\hbar^2}{2m} \frac{ \psi(x_{j-1} )-2\psi(x_j )+\psi(x_{j+1} )}{d^2} =E \psi(x_j )$. Write out the Schrödinger equation for the first three ($j = 1, 2, 3$) and the last three ($j = N-2, N-1, N$) points. b) It is convenient to write these $N$ equations in matrix form (one equation per row): ${\bf H}\psi(x)=E \psi(x)$ Matrix $\bf H$ is called a Hamiltonian. The equation is an eigenvalue equation, where $\psi(x)$ are the eigenvectors and $E$ are the eigenvalues. In this case it is a symmetric tri-diagonal matrix of size $N \times N$, and the boundary conditions are $\psi(x_0)= \psi(x_{N+1}) =0$. Write out the Hamiltonian $\bf H$ in matrix form so you see what it looks like. It is useful to think of how matrix multiplication works and the fact that each row must represent each of the equations in (a). ![](images/1DQM/Hmatrix.jpg "") #Exercise 4# a) Input matrix $\bf H$ into your computer program in such a way that it is easy to change the matrix size by changing the value of $N$. Use $m = 9.11 \times 10^{-31} kg$, $L = 10 nm$, and $N = 100$. Be careful with units! b) Find the five lowest eigenvalues (allowed energies for $n = 1$ through $n = 5$) and eigenvectors (the wave functions) of matrix $\bf{H}$. In MATLAB you can use the eigs() function. In Python you can use the eig() function in the numpy.linalg package. In Mathematica you can use Eigensystems[ ]. c) For the infinite well, how do the numerical and exact eigenvalues compare? Change the number of points $N$ to 10, 100, and 1000 points. How do the numeric and analytic eigenvalues (allowed energies) compare as you change $N$? Are the numeric and analytic allowed energies closer for low $n$ or high $n$? [Note: notice that $N$ and $n$ mean different things!] Explain. #Exercise 5# a) Make a plot with 10 subplots (two columns with five rows each. In the first column plots the first five allowed wave functions, $\psi (x)$ vs. $x$, for both the numeric solution and the analytic solution (for each value of $n$ plot the numeric and analytic function on the same plot so you can compare them). For the analytic solution, use a large enough number of points to get a smooth curve. Values on the x-axis should be in nanometers [nm]. In the second column plot the first five probability densities, $|\psi (x)|^2$, for both the numeric solution and the analytic solution (for each value of $n$ plot the numeric and analytic function on the same plot so you can compare them). For the analytic solution, use a large enough number of points to get a smooth curve. Note: The probability densities should be normalized to 1. You may need to multiply the eigenvectors by an appropriate normalization constant and multiply some by -1 for the phases of the numeric and exact solutions to match. b) Change the number of points $N$ to 10, 100, and 1000 points. Notice what happens to your plots. How do the numeric and analytic wave functions and probability densities compare as you change $N$? Are the numeric and analytic wave functions closer for low $n$ or high $n$? [Note: notice that $N$ and $n$ mean different things!] Explain. #Exercise 6# a) Plot potential energy vs. x for an infinite well with a step: $V \rightarrow \infty$ for $x < 0$, $V = 0$ for $0 < x < \frac{L}{2}$, $V = 5 \frac{\pi^2 \hbar^2 }{2mL^2}$ for $\frac{L}{2} < x < L$, $V \rightarrow \infty$ for $x > L$. b) Modify matrix $H$ for this potential energy and find the first allowed energies. Compare to the first five analytic eigenvalues which are: 0.0087 eV, 0.0271 eV, 0.0429 eV, 0.0705 eV, and 0.1032 eV. c) Plot the first five wave functions and probability densities. Do not worry about normalizing the wave functions and probability densities in this case. Comment on the shape of the plots. #Exercise 7# a) Plot potential energy vs. x for a parabolic well: $V = \frac{1}{2}m\omega^2 (x-\frac{L}{2})^2$ for $0 < x < L$. Take $m = 9.11 \times 10^{-31} kg$ and $\omega = 1 \times 10^{14} Hz.$ b) Modify matrix $H$ for this potential energy and find the first allowed energies. Compare to the first five analytic eigenvalues which are given by $E_n = (n+\frac{1}{2})\hbar\omega$, for $n=0,1,2,\dots$. c) Plot the first five wave functions and probability densities. Do not worry about normalizing the wave functions and probability densities in this case. Comment on the shape of the plots. #Exercise 8# a) Plot potential energy vs. x for half a parabolic well: $V \rightarrow \infty$ for $x < 0$, $V = \frac{1}{2}m\omega^2 x^2$ for $x > 0$. Take $m = 9.11 \times 10^{-31} kg$ and $\omega = 1 \times 10^{14} Hz.$ b) Modify matrix $H$ for this potential energy and find the first allowed energies. Compare to the first five analytic eigenvalues which are given by $E_n = (2n-\frac{1}{2})\hbar\omega$, for $n=1,2,3,\dots$. c) Plot the first five wave functions and probability densities. Do not worry about normalizing the wave functions and probability densities in this case. Comment on the shape of the plots.