{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Shooting Method with Numerov Algorithm\n", "\n", "\n", "**Authors:** Dou Du, Taylor James Baird and Giovanni Pizzi \n", "\n", " Go back to index\n", "\n", "**Source code:** https://github.com/osscar-org/quantum-mechanics/blob/master/notebook/quantum-mechancis/shooting_method.ipynb\n", "\n", "This notebook demonstrates the shooting method with the Numerov algorithm to search the \n", "eigenfunctions (wavefunctions) and eigenvalues for a one-dimensional quantum well.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## **Goals**\n", "\n", "* Understand the mathematical method to solve the Schrödinger equation \n", "numerically and the boundary condition for the 1D quantum well.\n", "* Understand the Numerov algorithm and how to improve the accuracy with \n", "high order correction.\n", "* Know how to use the shooting method with the Numerov algorithm to get \n", "the eigenvalues and eigenfunctions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## **Background theory**\n", "\n", "[More on the background theory.](./theory/theory_shooting_method.ipynb)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## **Tasks and exercises**\n", "\n", "1. Move the sliders for the width and depth of the quantum well. Do you understand\n", "the concept of quantum confinement? Do you know any numerical method to solve \n", "the Schrödinger equation for 1D quantum well?\n", "\n", "
\n", " Solution\n", " Please check the previous notebooks for the \n", " 1D quantum well.\n", " In that notebook, the one-dimensional Shrödinger equation was solved\n", " by numerical matrix diagonalization.\n", "
\n", "\n", "2. With the default width (1.20) and depth (0.20), move the sliders \n", "(on the left side) to the targeted energies. Report the energy when the tail \n", "of the wavefunction on the right converge to zero (line color turns to green). \n", "Is the energy the same as the eigenvalue shown in the right plot? You can also \n", "use the \"auto search\" button to get the eigenvalues, which searches the next \n", "solution when increasing the energy (i.e. it searches always upwards).\n", "
\n", " Solution\n", " The 1st eigenvalue is about 0.0092. You may need to click the \"Flip\n", " eigenfunctions\" button to make the comparsion. Check the exact eigenvalue \n", " by clicking on the eigenfunction in the plot.\n", "
\n", " \n", "3. Follow the same step to get all the eigenvalues, and make a table to compare \n", "the results with the eigenvalues from the figure. Compare the results with \n", "and without using the 4th derivative correction (checkbox). Which values \n", "should be more accurate and why?\n", "
\n", " Solution\n", " Please check the background theory section for the Numerov algorithm.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "## Interactive visualization\n", "(be patient, it might take a few seconds to load)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib widget\n", "\n", "from numpy import linspace, sqrt, ones, arange, diag, argsort, zeros\n", "from scipy.linalg import eigh_tridiagonal\n", "import matplotlib.pyplot as plt\n", "from ipywidgets import FloatSlider, jslink, VBox, HBox, Button, Label, Layout, Checkbox, Output\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "colors = ['#66c2a5','#fc8d62','#8da0cb','#e78ac3','#a6d854','#ffd92f'];\n", "ixx = 0\n", "\n", "def singlewell_potential(x, width, depth):\n", " x1 = ones(len(x))*depth\n", " for i in range(len(x)):\n", " if x[i] > - width/2.0 and x[i] < width/2.0:\n", " x1[i] =0\n", " return x1\n", "\n", "\n", "def diagonalization(hquer, L, N, pot=singlewell_potential, width = 0.1, depth = 0.0):\n", " \"\"\"Calculated sorted eigenvalues and eigenfunctions. \n", "\n", " Input:\n", " hquer: Planck constant\n", " L: set viewed interval [-L,L] \n", " N: number of grid points i.e. size of the matrix \n", " pot: potential function of the form pot\n", " x0: center of the quantum well\n", " width: the width of the quantum well\n", " depth: the depth of the quantum well\n", " Ouput:\n", " ew: sorted eigenvalues (array of length N)\n", " ef: sorted eigenfunctions, ef[:,i] (size N*N)\n", " x: grid points (arry of length N)\n", " dx: grid space\n", " V: Potential at positions x (array of length N)\n", " \"\"\"\n", " x = linspace(-L, L, N+2)[1:N+1] # grid points\n", " dx = x[1] - x[0] # grid spacing\n", " V = pot(x, width, depth)\n", " z = hquer**2 /2.0/dx**2 # second diagonals\n", "\n", " ew, ef = eigh_tridiagonal(V+2.0*z, -z*ones(N-1))\n", " ew = ew.real # real part of the eigenvalues\n", " ind = argsort(ew) # Indizes f. sort. Array\n", " ew = ew[ind] # Sort the ew by ind\n", " ef = ef[:, ind] # Sort the columns \n", " ef = ef/sqrt(np.sum(ef[0]*ef[0]*dx)) # Correct standardization \n", " return ew, ef, x, dx, V\n", "\n", "\n", "def plot_eigenfunctions(ax, ew, ef, x, V, width=1, updateTarget=True):\n", " \"\"\"Plot of the lowest eigenfunctions 'ef' in the potential 'V (x)'\n", " at the level of the eigenvalues 'ew' in the plot area 'ax'.\n", " \"\"\"\n", " global lnum, lax1, lspan\n", " \n", " fak = sfak.value/(50.0);\n", " \n", " try:\n", " lspan.remove()\n", " except:\n", " pass\n", " \n", " lspan = ax[0].axhspan(max(V), max(V)+0.05, facecolor='lightgrey')\n", " \n", " ax[0].set_xlim([min(x), max(x)])\n", " ax[0].set_ylim([min(V)-0.05, max(V)+0.05])\n", " \n", " ax[0].set_xlabel(r'$x/a$', fontsize = 10)\n", " ax[0].set_ylabel(r'$V(x)/V_0\\ \\rm{, Eigenfunctions\\ with\\ Eigenvalues}$', fontsize = 10)\n", " \n", " ax[1].set_xlim([min(x), max(x)])\n", " ax[1].set_ylim([min(V)-0.05, max(V) + 0.05])\n", " \n", " if updateTarget:\n", " loop1.min = min(V)-0.03\n", " loop1.min = int(loop1.min*100)/100.0\n", " loop1.value = loop1.min\n", " loop2.value = loop2.min\n", " loop3.value = loop3.min\n", " loop4.value = loop4.min\n", " \n", " ax[1].yaxis.set_label_position(\"right\")\n", " ax[1].yaxis.tick_right()\n", " \n", " ax[1].get_xaxis().set_visible(False)\n", " #ax[1].set_ylabel(r'$\\rm{\\ Eigenvalues}$', fontsize = 10)\n", " \n", " indmax = sum(ew < max(V)) \n", " if not hasattr(width, \"__iter__\"): \n", " width = width*ones(indmax) \n", " for i in arange(indmax): \n", " ax[0].plot(x, fak*ef[:, i]+ew[i], linewidth=width[i]+.1, color=colors[i%len(colors)])\n", " ax[1].plot(x, x*0.0+ew[i], linewidth=width[i]+2.5, color=colors[i%len(colors)])\n", " \n", " ax[0].plot(x, V, c='k', linewidth=1.6)\n", " lnum, = ax[0].plot(x, x*0 + loop1.value,'r--', linewidth=1.0)\n", " lax1, = ax[1].plot(x, x*0 + loop1.value,'r--', linewidth=1.0)\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mu = 0.06 # Potential parameter\n", "L = 1.5 # x range [-L,L]\n", "N = 200 # Number of grid points\n", "hquer = 0.06 # Planck constant\n", "sigma_x = 0.1 # Width of the Gaussian function\n", "zeiten = linspace(0.0, 10.0, 400) # time\n", "Flip = False # Flip the eigenfunction\n", "\n", "\n", "swidth = FloatSlider(value = 1.2, min = 0.1, max = 2.0, description = 'Width: ')\n", "sdepth = FloatSlider(value = 0.2, min = 0.05, max = 1.0, step = 0.05, description = 'Depth: ')\n", "sfak = FloatSlider(value = 3.0, min = 1.0, max = 5.0, step = 0.5, description = r'Zoom factor: ')\n", "\n", "output = Output()\n", "\n", "update = Button(description=\"Show all\")\n", "flip = Button(description=\"Flip eigenfunction\")\n", "search = Button(description=\"Auto search\")\n", "\n", "order = Checkbox(value=True, description=\"incl. 4th derivative\", indent=False,\n", " layout=Layout(width='180px'))\n", "\n", "loop1 = FloatSlider(value = -0.03, min = -0.03, max = 0.2,\n", " layout=Layout(height='450px', width='30px'), step = 0.01, readout_format=\".2f\", orientation='vertical')\n", "loop2 = FloatSlider(value = 0, min = 0, max = 99, \n", " layout=Layout(height='450px', width='30px'), step =1.0, readout_format='02d', orientation='vertical')\n", "loop3 = FloatSlider(value = 0, min = 0, max = 99,\n", " layout=Layout(height='450px', width='30px'), step =1.0, readout_format='02d', orientation='vertical')\n", "loop4 = FloatSlider(value = 0, min = 0, max = 99,\n", " layout=Layout(height='450px', width='30px'), step =1.0, readout_format='02d', orientation='vertical')\n", "\n", "\n", "Leng = Label('')\n", "Evalue = loop1.value + loop2.value/10000.0 + loop3.value/1000000.0 + loop4.value/100000000.0;\n", "Leng.value = \"Current value: \" + \"{:.8f}\".format(Evalue)\n", "\n", "width = 1.2\n", "depth = 0.2\n", "fak = 5.0\n", "\n", "ew, ef, x, dx, V = diagonalization(hquer, L, N, width = width, depth = depth)\n", " \n", "with output:\n", " global fig\n", " fig, ax = plt.subplots(1, 2, figsize=(6,6), gridspec_kw={'width_ratios': [10, 1]})\n", " fig.canvas.header_visible = False\n", " fig.canvas.layout.width = \"750px\"\n", "\n", " fig.suptitle('Numerical Solution ($\\psi$) of the Schrödinger Equation \\n for 1D Quantum Well', fontsize = 12)\n", " plot_eigenfunctions(ax, ew, ef, x, V)\n", " plt.show()\n", "\n", "def Numerov(y, E, Vn, dxx):\n", " y = zeros(len(y));\n", " y[0] = 0.0;\n", " \n", " if Flip:\n", " y[1] = -0.00000001\n", " else:\n", " y[1] = 0.00000001\n", " \n", " k2 = 2.0/(hquer**2)*(E-Vn)*dxx*dxx;\n", " \n", " for i in arange(2, len(y)):\n", " if order.value:\n", " y[i] = (2*(12.0-5.0*k2[i-1])*y[i-1] - (12+k2[i-2])*y[i-2])/(12+k2[i]);\n", " else:\n", " y[i] = 2*y[i-1] - k2[i-1]*y[i-1] - y[i-2]\n", " return y/(sqrt(np.sum(abs(y)**2*dxx))*50.0)*sfak.value\n", "\n", "def plot_numerov(c):\n", " Nn = 1000\n", " xx = linspace(-L, L, Nn+2)[1:Nn+1]\n", " dxx = xx[1] - xx[0];\n", " Vn = singlewell_potential(xx, width = swidth.value, depth = sdepth.value)\n", " yy = zeros(len(xx));\n", " \n", " Evalue = loop1.value + loop2.value/10000.0 + loop3.value/1000000.0 + loop4.value/100000000.0;\n", " \n", " yy = Numerov(yy, Evalue, Vn, dxx);\n", " \n", " if abs(yy[-1]) < 0.001:\n", " lnum.set_color(\"green\")\n", " lax1.set_color(\"green\")\n", " else:\n", " lnum.set_color(\"red\")\n", " lax1.set_color(\"red\")\n", "\n", " Leng.value = \"Current value: \" + \"{:.8f}\".format(Evalue) \n", " lnum.set_data(xx, yy + Evalue)\n", " lax1.set_data(xx, xx*0 + Evalue)\n", " \n", "def on_auto_search(b):\n", " Nn = 1000\n", " xx = linspace(-L, L, Nn+2)[1:Nn+1]\n", " dxx = xx[1] - xx[0];\n", " Vn = singlewell_potential(xx, width = swidth.value, depth = sdepth.value)\n", " yy = zeros(len(xx));\n", " \n", " Evalue = loop1.value + loop2.value/10000.0 + loop3.value/1000000.0 + loop4.value/100000000.0;\n", " yy = Numerov(yy, Evalue, Vn, dxx);\n", " \n", " increment = 0.01\n", " while abs(yy[-1]) > 0.001:\n", " tail_old = yy[-1]\n", " Evalue += increment;\n", " yy = Numerov(yy, Evalue, Vn, dxx);\n", " tail_new = yy[-1]\n", " \n", " if tail_old*tail_new < 0:\n", " Evalue -= increment\n", " increment /= 100.0\n", " yy = Numerov(yy, Evalue, Vn, dxx);\n", " \n", " Leng.value = \"Current value: \" + \"{:.8f}\".format(Evalue) \n", " lnum.set_data(xx, yy + Evalue)\n", " lax1.set_data(xx, xx*0 + Evalue)\n", " \n", " loop1.value = int(Evalue*100)/100.0;\n", " loop2.value = int((Evalue-loop1.value)*10000); \n", " loop3.value = int((Evalue-loop1.value-loop2.value/10000)*1000000);\n", " loop4.value = int((Evalue-loop1.value-loop2.value/10000-loop3.value/1000000)*100000000)\n", " \n", "\n", "def on_update_click(b):\n", " for i in ax[0].lines:\n", " i.set_alpha(1.0)\n", " for i in ax[1].lines:\n", " i.set_alpha(1.0)\n", "\n", " try:\n", " ann.remove()\n", " ann1.remove()\n", " except:\n", " pass\n", "\n", "def on_width_change(change):\n", " global ew, ef, x, dx, V\n", " ax[0].lines = []\n", " ax[1].lines = []\n", " \n", " try:\n", " ann.remove()\n", " ann1.remove()\n", " except:\n", " pass\n", "\n", " ew, ef, x, dx, V = diagonalization(hquer, L, N, width = swidth.value, depth = sdepth.value)\n", " plot_eigenfunctions(ax, ew, ef, x, V)\n", "\n", "def on_depth_change(change):\n", " global ew, ef, x, dx, V\n", " ax[0].lines = []\n", " ax[1].lines = []\n", " \n", " try:\n", " ann.remove()\n", " ann1.remove()\n", " except:\n", " pass\n", "\n", " ew, ef, x, dx, V = diagonalization(hquer, L, N, width = swidth.value, depth = sdepth.value)\n", " plot_eigenfunctions(ax, ew, ef, x, V)\n", " loop1.max = max(V)\n", " \n", "def on_xfak_change(change):\n", " ax[0].lines = []\n", " ax[1].lines = []\n", " \n", " try:\n", " ann.remove()\n", " ann1.remove()\n", " except:\n", " pass\n", "\n", " plot_eigenfunctions(ax, ew, ef, x, V, updateTarget=False)\n", " plot_numerov('test')\n", "\n", "def on_press(event):\n", " global ann, ann1, ixx\n", " \n", " ixx = min(enumerate(ew), key = lambda x: abs(x[1]-event.ydata))[0]\n", " \n", " for i in range(len(ax[1].lines)-1):\n", " ax[0].lines[i].set_alpha(0.1)\n", " ax[1].lines[i].set_alpha(0.1)\n", " ax[0].lines[i].set_linewidth(1.1)\n", " \n", " ax[0].lines[ixx].set_alpha(0.5)\n", " ax[1].lines[ixx].set_alpha(0.5)\n", " ax[0].lines[ixx].set_linewidth(2.0)\n", " \n", " try:\n", " ann.remove()\n", " ann1.remove()\n", " except:\n", " pass\n", " \n", " ann = ax[0].annotate(s = 'n = ' + str(ixx+1), xy = (0, ew[ixx]), xytext = (-0.15, ew[ixx]), xycoords = 'data', color='k', size=15)\n", " ann1 = ax[1].annotate(s = str(\"{:.3f}\".format(ew[ixx])), xy = (0, ew[ixx]), xytext = (-1.2, ew[ixx]+0.005), xycoords = 'data', color='k', size=9)\n", "\n", "def on_flip_eigenfunctions(b):\n", " global Flip\n", " x = lnum.get_xdata();\n", " y = lnum.get_ydata();\n", " lnum.set_data(x, -y+2.0*y[0])\n", " Flip = not Flip\n", "\n", " \n", "cid = fig.canvas.mpl_connect('button_press_event', on_press)\n", "\n", "swidth.observe(on_width_change, names = 'value')\n", "sdepth.observe(on_depth_change, names = 'value')\n", "sfak.observe(on_xfak_change, names = 'value')\n", "\n", "update.on_click(on_update_click)\n", "flip.on_click(on_flip_eigenfunctions)\n", "search.on_click(on_auto_search)\n", "\n", "loop1.observe(plot_numerov, names = 'value')\n", "loop2.observe(plot_numerov, names = 'value')\n", "loop3.observe(plot_numerov, names = 'value')\n", "loop4.observe(plot_numerov, names = 'value')\n", "\n", "\n", "label1 = Label(value=\"Targeted eigenvalue\")\n", "label2 = Label(value=\"Click to flip the eigenfunction\")\n", "label3 = Label(value=\"(click on a state to select it)\")\n", "label4 = Label(value=\"(tune to zoom the eigenfunctions)\")\n", "\n", "display(HBox([VBox([label1, HBox([loop1, loop2, loop3, loop4]), Leng, search, order, label2, flip]), output]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set the **width** and **depth** of the quantum well:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "display(HBox([swidth, sdepth]), VBox([HBox([sfak, label4]), HBox([update, label3])]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "# Legend\n", "\n", "(How to use the interactive visualization)\n", "\n", "## Interactive figures\n", "\n", "In the interative figure, the soild lines show the wavefunctions and their\n", "corresponding eigenvalues, which are solved by matrix diagonalization. \n", "There is a red dash line at the bottom of the figure, which shows the \n", "eigenfunction solved by Numerov algorithm. \n", "\n", "## Controls\n", "\n", "There are four vertical sliders to control the targeted eigenvalue E. The first \n", "slider controls the precision for tenths ($10^{-1}$) and hundredths ($10^{-2}$). \n", "The second slider controls thousandths ($10^{-3}$) and ten thousandths decimal ($10^{-4}$). The third slider controls hundred thousandths ($10^{-5}$) and \n", "millionths ($10^{-6}$). The last slider controls ten millionths ($10^{-7}$) \n", "and hundred millionths ($10^{-8}$). The current value is also displayed under \n", "the sliders.\n", "\n", "You need slowly move the 1st slider and observe the tail of the dashed line on \n", "the right edge. Once you see the tail change directions (up or down), the true \n", "value should be between these two values. You need to go back to a smaller value \n", "and start to tune the 2nd slider. Then the same procedure is for the 3rd and 4th\n", "slider. When the absolute value at the right edge is smaller than 0.001, the \n", "dashed red line will turn green. It reaches the desired accuracy for the \n", "wavefunction. Then, you can read out the current targeted value, which is the\n", "corresponding eigenvalue.\n", "\n", "You can also use the `Auto search` button, which finds the closest eigenvalue \n", "and eigenfunction (search in the upward direction). In order to make a comparison, \n", "you may also need to click the `Flip eigenfunctions` button." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.9.7" }, "voila": { "authors": "Dou Du and Giovanni Pizzi" } }, "nbformat": 4, "nbformat_minor": 4 }