{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Lecture 11: Eigenvalue Problem Demo\n", "\n", "We solve the time-independent Schrödinger equation\n", "\n", "$$\n", "\\hat H\\phi_n(x) = -\\frac{1}{2}\\frac{d^2\\phi_n(x)}{dx^2} + V(x)\\,\\phi_n(x) = E_n\\,\\phi_n(x)\n", "$$\n", "\n", "on a grid, in units with $\\hbar = m = \\omega = 1$, for the harmonic oscillator $V(x) = \\frac{1}{2}x^2$.\n", "We know the exact answer, $E_n = n + \\frac{1}{2}$, so this is where we test the method before trusting it on a problem we cannot solve by hand." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import time\n", "from math import factorial\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import scipy.linalg\n", "import scipy.sparse\n", "import scipy.sparse.linalg\n", "import scipy.special\n", "\n", "COLORS = [\"#2a78d6\", \"#eb6834\", \"#1baf7a\", \"#eda100\"] # one color per state n = 0, 1, 2, 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. The Hamiltonian matrix\n", "\n", "On the grid $x_i = -L/2 + i\\,\\Delta x$, $i = 0, \\dots, N_D$, with $\\Delta x = L/N_D$,\n", "\n", "$$\n", "H_{ij} = -\\frac{1}{2\\Delta x^2}\\left(\\delta_{i+1,j} - 2\\delta_{ij} + \\delta_{i-1,j}\\right) + V(x_i)\\,\\delta_{ij}.\n", "$$\n", "\n", "Only three diagonals are nonzero, so we store $H$ as a sparse matrix: about $3(N_D + 1)$ numbers instead of $(N_D + 1)^2$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def V(x):\n", " return 0.5 * x**2\n", "\n", "\n", "def make_grid(L, ND):\n", " # ND + 1 points from -L/2 to L/2, and the spacing\n", " x = np.linspace(-L / 2, L / 2, ND + 1)\n", " return x, x[1] - x[0]\n", "\n", "\n", "def hamiltonian(x, dx, V):\n", " # sparse tridiagonal finite-difference Hamiltonian\n", " off = np.full(len(x) - 1, -0.5 / dx**2)\n", " diag = 1.0 / dx**2 + V(x)\n", " return scipy.sparse.diags([off, diag, off], [-1, 0, 1], format=\"csc\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check against the slides: $L = 8$ and $N_D = 4$ (so $\\Delta x = 2$) should give $H_{00} = 8.25$ and off-diagonal entries $-0.125$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x, dx = make_grid(8, 4)\n", "print(hamiltonian(x, dx, V).toarray())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Three ways to find the lowest eigenvalues\n", "\n", "**Predict first:** for $N_D = 3000$, which of these is fastest, and which is slowest?\n", "\n", "1. `scipy.linalg.eigh(H)`: dense matrix, computes **all** $N_D + 1$ eigenpairs.\n", "2. `scipy.sparse.linalg.eigsh(H, k=2, which='SM')`: ARPACK, asks directly for the 2 eigenvalues of smallest magnitude.\n", "3. `scipy.sparse.linalg.eigsh(H, k=2, sigma=0, which='LM')`: ARPACK in shift-invert mode, which finds the largest eigenvalues of $(H - \\sigma)^{-1}$, i.e. the $E_n$ closest to $\\sigma$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def best_time(f, repeat=3):\n", " # best wall-clock time over a few runs (like %timeit, but bounded)\n", " times = []\n", " for _ in range(repeat):\n", " t0 = time.perf_counter()\n", " f()\n", " times.append(time.perf_counter() - t0)\n", " return min(times)\n", "\n", "\n", "L = 8\n", "for ND in [600, 3000]:\n", " x, dx = make_grid(L, ND)\n", " H = hamiltonian(x, dx, V)\n", " H_dense = H.toarray()\n", " t_dense = best_time(lambda: scipy.linalg.eigh(H_dense))\n", " t_sm = best_time(lambda: scipy.sparse.linalg.eigsh(H, k=2, which=\"SM\"), repeat=1)\n", " t_si = best_time(lambda: scipy.sparse.linalg.eigsh(H, k=2, sigma=0, which=\"LM\"))\n", " print(f\"N_D = {ND}\")\n", " print(f\" 1. dense eigh (all eigenpairs) {t_dense:9.4f} s\")\n", " print(f\" 2. eigsh, which='SM' {t_sm:9.4f} s\")\n", " print(f\" 3. eigsh, shift-invert sigma=0 {t_si:9.4f} s\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Why is `which='SM'` so slow? ARPACK converges quickly to eigenvalues that stand out in magnitude from the rest of the spectrum.\n", "The largest eigenvalue of $H$ is about $2/\\Delta x^2$ (about $3 \\times 10^5$ at $N_D = 3000$), so on that scale $E_0 = 0.5$ and $E_1 = 1.5$ are bunched together near zero.\n", "Shift-invert maps $E_n \\to 1/(E_n - \\sigma)$. With $\\sigma = 0$ the ground state becomes the largest eigenvalue, 3 times larger than the next one, and ARPACK converges in a few iterations.\n", "Each iteration solves a tridiagonal linear system, which costs only $O(N_D)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. From eigenvectors to wave functions\n", "\n", "`eigsh` returns unit vectors, $\\sum_i |v_i|^2 = 1$. A wave function needs\n", "\n", "$$\n", "\\int |\\phi(x)|^2\\, dx \\approx \\sum_i |\\phi_i|^2\\, \\Delta x = 1\n", "\\qquad\\Rightarrow\\qquad \\phi_i = \\frac{v_i}{\\sqrt{\\Delta x}}.\n", "$$\n", "\n", "The overall sign of an eigenvector is arbitrary (ARPACK starts from a random vector), so a plot can flip from one run to the next.\n", "We fix it to match the exact Hermite functions, which are positive at large $x$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def lowest_states(H, dx, k, sigma=0):\n", " # k eigenvalues closest to sigma (ascending) and normalized wave functions (columns of psi)\n", " E, v = scipy.sparse.linalg.eigsh(H, k=k, sigma=sigma, which=\"LM\")\n", " order = np.argsort(E)\n", " E, psi = E[order], v[:, order] / np.sqrt(dx)\n", " for n in range(k):\n", " # sign convention: the outermost lobe on the right is positive\n", " lobe = np.flatnonzero(np.abs(psi[:, n]) > 0.01 * np.abs(psi[:, n]).max())[-1]\n", " psi[:, n] *= np.sign(psi[lobe, n])\n", " return E, psi\n", "\n", "\n", "def phi_exact(n, x):\n", " return scipy.special.eval_hermite(n, x) * np.exp(-x**2 / 2) / np.sqrt(2**n * factorial(n) * np.sqrt(np.pi))\n", "\n", "\n", "L, ND = 8, 600\n", "x, dx = make_grid(L, ND)\n", "E, psi = lowest_states(hamiltonian(x, dx, V), dx, k=6)\n", "\n", "# overlaps sum_i phi_m(x_i) phi_n(x_i) dx should be the identity matrix\n", "print(np.round(psi[:, :4].T @ psi[:, :4] * dx, 6) + 0.0) # + 0.0 turns -0.0 into 0.0" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(9, 5.5))\n", "ax.plot(x, V(x), color=\"0.25\", lw=2, label=\"$V(x)$\")\n", "for n in range(4):\n", " ax.axhline(E[n], color=COLORS[n], lw=1, ls=\"--\")\n", " ax.plot(x, phi_exact(n, x) + E[n], color=\"0.8\", lw=6) # exact, drawn underneath\n", " ax.plot(x, psi[:, n] + E[n], color=COLORS[n], lw=2, label=rf\"$\\phi_{n}(x) + E_{n}$\")\n", "ax.plot([], [], color=\"0.8\", lw=6, label=\"exact\")\n", "ax.set_xlabel(\"$x$\")\n", "ax.set_ylabel(r\"$V(x)$ or $\\phi_n(x) + E_n$\")\n", "ax.set_xlim(-L / 2, L / 2)\n", "ax.set_ylim(0, 4.5)\n", "ax.legend(loc=\"upper left\", bbox_to_anchor=(1.02, 1), frameon=False)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Validation: compare with the exact energies\n", "\n", "How close are we to $E_n = n + \\frac{1}{2}$?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\" n numerical exact error\")\n", "for n in range(6):\n", " print(f\"{n:2d} {E[n]:.8f} {n + 0.5:11.1f} {E[n] - (n + 0.5):+.1e}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The low states come out slightly **too low**, while the higher states come out **too high** by much more.\n", "Two different errors are at work. The next two cells separate them.\n", "\n", "### Error 1: the grid spacing $\\Delta x$\n", "\n", "Taylor expanding the finite difference,\n", "\n", "$$\n", "\\frac{\\phi_{i+1} - 2\\phi_i + \\phi_{i-1}}{\\Delta x^2} = \\phi''(x_i) + \\frac{\\Delta x^2}{12}\\,\\phi''''(x_i) + O(\\Delta x^4),\n", "$$\n", "\n", "so the discrete kinetic energy is off by $-\\frac{\\Delta x^2}{24}\\phi''''$.\n", "First-order perturbation theory then predicts $\\Delta E_n \\approx -\\frac{\\Delta x^2}{24}\\langle p^4\\rangle_n = -\\frac{\\Delta x^2}{32}\\left(2n^2 + 2n + 1\\right)$ for the harmonic oscillator.\n", "To see this error on its own, we use a large box ($L = 10$) and shrink $\\Delta x$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "L = 10\n", "NDs = np.array([50, 100, 200, 400, 800, 1600])\n", "dxs = L / NDs\n", "errors = []\n", "for ND in NDs:\n", " x, dx = make_grid(L, ND)\n", " E_ND, _ = lowest_states(hamiltonian(x, dx, V), dx, k=4)\n", " errors.append(E_ND - (np.arange(4) + 0.5))\n", "errors = np.array(errors)\n", "\n", "fig, ax = plt.subplots(figsize=(7, 5))\n", "for n in range(4):\n", " predicted = dxs**2 / 32 * (2 * n**2 + 2 * n + 1)\n", " ax.loglog(dxs, predicted, color=COLORS[n], lw=1, ls=\"--\")\n", " ax.loglog(dxs, np.abs(errors[:, n]), \"o\", color=COLORS[n], ms=8, label=f\"$n = {n}$\")\n", "ax.plot([], [], color=\"0.4\", lw=1, ls=\"--\", label=r\"predicted $\\Delta x^2 (2n^2 + 2n + 1)/32$\")\n", "ax.set_xlabel(r\"$\\Delta x$\")\n", "ax.set_ylabel(r\"$|E_n - (n + \\frac{1}{2})|$\")\n", "ax.legend(frameon=False)\n", "plt.show()\n", "\n", "print(\"ratio of errors when dx is halved:\")\n", "print(np.round(errors[:-1] / errors[1:], 2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Slope 2 on the log-log plot, and halving $\\Delta x$ cuts the error by a factor of 4: the method is second order, and the points sit on the prediction.\n", "\n", "### Error 2: the box size $L$\n", "\n", "The matrix is equivalent to hard walls just outside $x = \\pm L/2$, where $\\phi = 0$.\n", "The exact $\\phi_n$ leaks past the classical turning point $x_n = \\sqrt{2n + 1}$ and decays like $e^{-x^2/2}$.\n", "If a wall cuts into that tail, the state is squeezed and its energy goes **up**.\n", "Fix $\\Delta x = 0.02$ and vary $L$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"turning points x_n:\", np.round(np.sqrt(2 * np.arange(6) + 1), 2))\n", "print(\" L E_n - (n + 1/2) for n = 0, ..., 5\")\n", "for L in [6, 8, 10, 12]:\n", " x, dx = make_grid(L, round(L / 0.02))\n", " E_L, _ = lowest_states(hamiltonian(x, dx, V), dx, k=6)\n", " print(f\"{L:4d} \" + \" \".join(f\"{err:+.1e}\" for err in E_L - (np.arange(6) + 0.5)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The $L = 10$ and $L = 12$ rows agree to within $10^{-4}$: once the walls are far enough away, only the small, negative grid error from Error 1 is left.\n", "With $L = 8$ (walls at $\\pm 4$), $E_0$ and $E_1$ are fine, but the walls push $E_3$ up by about $10^{-3}$ and $E_5$ (turning point $x_5 \\approx 3.3$) by about $4 \\times 10^{-2}$.\n", "\n", "## 5. Your turn: Assignment 2\n", "\n", "Nothing above is specific to the harmonic oscillator except `V(x)` and the exact answers.\n", "For a potential without exact answers, check the results yourself before trusting them:\n", "\n", "- **Grid:** halve $\\Delta x$. The change in $E_n$ should shrink by about 4 each time.\n", "- **Box:** make $L$ larger. The energies you report should not move.\n", "- **Normalization:** divide by $\\sqrt{\\Delta x}$ and check $\\sum_i |\\phi_i|^2\\,\\Delta x = 1$ before plotting.\n", "- **Shift:** `sigma` should sit at or below the lowest energy you expect. Shift-invert finds the $E_n$ closest to $\\sigma$, so $\\sigma = 0$ works here only because every $E_n > 0$.\n", "\n", "**Looking ahead:** the eigenpairs give the propagator, $K(x, t; x', 0) = \\sum_n \\phi_n(x)\\,\\phi_n^*(x')\\,e^{-iE_n t/\\hbar}$, and the partition function, $Z(\\beta) = \\sum_n e^{-\\beta E_n}$." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.10" } }, "nbformat": 4, "nbformat_minor": 2 }