{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Quantum Numerov; an Improved Schrodinger Solver

\n", "\n", "The Numerov method is a highly efficient solver for ODEs involving only \n", "second derivatives, such as the Schrodinger equation:\n", "

\n", "$$ \n", "\\psi(x+h) \\simeq \\frac{2 [1-\n", "\\frac{5}{12}h^2k^2(x) ]\\psi(x) - [1+ \\frac{h^2} {12} k^2(x-h) ]\n", "\\psi(x-h)} { 1+h^2 k^2(x+h)/12} .\n", "$$\n", " \n", "Here it is used to match right and left wave functions via bisection algorihtm, and present results with Vpython." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/html": [ "

" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") { window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")};}else{ element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") { window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")};}else{ element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") { window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")};}else{ element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "\"\"\" From \"COMPUTATIONAL PHYSICS\" & \"COMPUTER PROBLEMS in PHYSICS\"\n", " by RH Landau, MJ Paez, and CC Bordeianu (deceased).\n", " Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, \n", " C Bordeianu, Univ Bucharest, 2021. \n", " Please respect copyright & acknowledge our work.\"\"\"\n", "\n", "# QuantumNumerovVP: Schodinger solver using Vpython \n", "\n", "import numpy as np\n", "\n", "psigr = canvas(width = 600, height = 300)\n", "psi = curve(color = color.yellow,canvas = psigr)\n", "psigr2 = canvas(y = 300,height = 200,width = 600)\n", "psi2 = curve(color = color.cyan,canvas = psigr2)\n", "energr = canvas(y = 500,width = 600,height = 200, title = 'Potential and Energy')\n", "poten = curve(color = color.cyan,canvas = energr) \n", "autoen = curve(x = list(range(0,1000)),canvas = energr)\n", "\n", "dl = 1e-6 # for bisection\n", "ul = [0.]*(1501)\n", "ur = [0.]*(1501)\n", "k2l = [0.]*(1501) # k**2 left wavefunc.\n", "k2r = [0.]*(1501) # k**2 Right wavefunc.\n", "n = 1501\n", "m = 5 # plot every 5 points\n", "imax = 100 \n", "xl0 = -1000 # leftmost x \n", "xr0 = 1000 # rightmost x \n", "h = (1.0*(xr0-xl0)/(n-1)) # 1.0 to float \n", "amin = -0.001 # root between amin and amax\n", "amax = -0.00085\n", "e = amin # initial guess for energy\n", "de = 0.01\n", "ul[0] = 0.0\n", "ul[1] = 0.00001\n", "ur[0] = 0.0\n", "ur[1] = 0.00001 \n", "im = 500 # match point \n", "nl = im+2 # for left psi\n", "nr = n-im+1 # for right psi\n", "istep = 0\n", "\n", "def V(x): # Square well Potential\n", " if (abs(x)< = 500): v = -0.001 \n", " else: v = 0\n", " return v\n", "\n", "def setk2(): # k2 = (sqrt(e-V))^2 \n", " for i in range(0,n): \n", " xl = xl0 + i*h\n", " xr = xr0 - i*h\n", " k2l[i] = e-V(xl)\n", " k2r[i] = e-V(xr)\n", "\n", "def numerov (n,h,k2,u): # Numerov algorithm \n", " b = (h**2)/12.0 \n", " for i in range( 1,n-1): \n", " u[i+1] = (2*u[i]*(1.0-5.*b*k2[i])-(1.+b*k2[i-1])*u[i-1])/(1.+b*k2[i+1])\n", " \n", "setk2()\n", "numerov (nl,h,k2l,ul) # left wave function\n", "numerov (nr,h,k2r,ur) # right wave function\n", "fact = ur[nr-2]/ul[im] # Rescale solution\n", "for i in range (0,nl): ul[i] = fact*ul[i]\n", "f0 = (ur[nr-1]+ul[nl-1]-ur[nr-3]-ul[nl-3])/(2*h*ur[nr-2]) # Log deriv\n", "dl = 1e-6 # for bisection\n", "ul = [0.]*(1501)\n", "ur = [0.]*(1501)\n", "k2l = [0.]*(1501) # k**2 L \n", "k2r = [0.]*(1501) # k**2 R\n", "n = 1501\n", "m = 5 # plot every 5 points\n", "imax = 100 \n", "xl0 = -1000 # leftmost x \n", "xr0 = 1000 # rightmost x \n", "h = (1.0*(xr0-xl0)/(n-1)) # 1.0 to float \n", "amin = -0.001 # root between amin and amax\n", "amax = -0.00085\n", "e = amin # initial energy guess\n", "de = 0.01\n", "ul[0] = 0.0\n", "ul[1] = 0.00001\n", "ur[0] = 0.0\n", "ur[1] = 0.00001 \n", "im = 500 \n", "nl = im+2 # match point L \n", "nr = n-im+1 # match point R \n", "istep = 0\n", "\n", "def normalize(): # normalize WF\n", " asum = 0\n", " for i in range( 0,n): \n", " if i > im :\n", " ul[i] = ur[n-i-1]\n", " asum = asum+ul[i]*ul[i]\n", " asum = sqrt(h*asum);\n", " elabel = label(pos = vector(700, 500,0), text = 'e = ', box = 0,display = psigr)\n", " elabel.text = (('e = %10.8f' )%e)\n", " ilabel = label(pos = vec(700, 400,0), tex15t = 'istep = ', box = 0,display = psigr)\n", " ilabel.text = (('istep = %4s') %istep)\n", " curve(vec(-1500,200,0),vec(-1000,200,0),vec(-1000,-200,0),vec(0,-200,0),vec(0,200,0),vec(1000,200,0))\n", " curve(vec(-1000,e*400000.0+200,0),vec(0,e*400000.0+200,0))\n", " label(pos = vec(-1150,-240,0),text = '0.001',box = 0,display = energr)\n", " label(pos = vec(-1000,300,0),text = '0',box = 0,display = energr)\n", " label(pos = vec(-900,180,0),text = '-500',box = 0,display = energr)\n", " label(pos = vec(-100,180,0),text = '500',box = 0,display = energr)\n", " label(pos = vec(-500,180,0),text = '0',box = 0,display = energr)\n", " label(pos = vec(900,120,0),text = 'r',box = 0,display = energr)\n", " j = 0\n", " psi2.clear()\n", " psi.clear()\n", " for i in range(0,n,m): \n", " xl = xl0+i*h\n", " ul[i] = ul[i]/asum #wave function normalized \n", " psi.append(pos = (vector(xl-500,10000*ul[i],0)))\n", " psi.visible = True\n", " rate(100)\n", " curve(vec(-830,-500,0),vec(-830,500,0),color = color.red,canvas = psigr)\n", " psi2.append(pos = (vec(xl-500, 1.0e5*ul[i]**2,0))) #plot psi\n", " psi2.vivible = True\n", " j += 1\n", " \n", "while abs(de) > dl and istep < imax : # bisection algorithm \n", " rate(100) # slowdown animation\n", " e1 = e # guessed root\n", " e = (amin+amax)/2 # half interval\n", " for i in range(0,n): \n", " k2l[i] = k2l[i]+ e-e1 \n", " k2r[i] = k2r[i]+ e-e1\n", " im = 500;\n", " nl = im+2\n", " nr = n-im+1;\n", " numerov (nl,h,k2l,ul) # Find wavefuntions for new k2l,k2r\n", " numerov (nr,h,k2r,ur) \n", " fact = ur[nr-2]/ul[im]\n", " for i in range(0,nl): ul[i] = fact*ul[i] \n", " f1 = (ur[nr-1]+ul[nl-1]-ur[nr-3]-ul[nl-3])/(2*h*ur[nr-2]) # Log deriv.\n", " if f0*f1 < 0: # bisection localize root\n", " amax = e\n", " de = amax-amin\n", " else:\n", " amin = e\n", " de = amax-amin\n", " f0 = f1\n", " normalize() \n", " istep = istep+1" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.13" } }, "nbformat": 4, "nbformat_minor": 2 }