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

Quantum Eigenvalue

\n", "

\n", " A nucleon is bound in a square well potential. Using solutions\n", " of Schrödinger equations from the left and from the right of the well, the bisection algorithm is used to search for an energy at which the log derivatives match, i.e. $\\Delta$ vanishes:\n", " \n", " $$\n", "\\Delta (E,x) = \\left. \\frac{\\psi'_L(x)/\\psi_L(x) -\n", "\\psi'_R(x)/\\psi_R(x)}{\\psi'_L(x)/\\psi_L(x)\n", " + \\psi'_R(x)/\\psi_R(x) }\\right|_{x=x_{match}} \n", " $$ " ] }, { "cell_type": "code", "execution_count": 1, "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": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/glow.min\");}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/glowcomm\");}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/jquery-ui.custom.min\");}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/glow.min\"], function(){console.log(\"GLOW LOADED\");});}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/glowcomm\"], function(){console.log(\"GLOWCOMM LOADED\");});}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/jquery-ui.custom.min\"], function(){console.log(\"JQUERY LOADED\");});}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, 2020. \n", " Please respect copyright & acknowledge our work.\"\"\"\n", "\n", "# QuantumEigen.py: Finds E and psi via rk4 + bisection\n", "\n", "%matplotlib notebook\n", "# mass/((hbar*c)**2)= 940MeV/(197.33MeV-fm)**2 =0.4829\n", "# well width=20.0 fm, depth 10 MeV, Wave function not normalized. \n", "\n", "from vpython import *\n", "import numpy as np\n", "\n", "psigr = display(x=0,y=0,width=600,height=300, title='R & L Wavefunc')\n", "Lwf = curve(x=list(range(502)),color=color.red)\n", "Rwf = curve(x=list(range(997)),color=color.yellow)\n", "eps = 1E-3 # Precision\n", "n_steps = 501 \n", "E = -17.0 # E guess\n", "h = 0.04\n", "count_max = 100\n", "Emax = 1.1*E # E limits\n", "Emin = E/1.1 \n", "F = np.zeros(2)\n", "y = np.zeros(2)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Final eigenvalue E = -15.952090037952768\n", "iterations, max = 15\n" ] } ], "source": [ "def f(x, y, F,E):\n", " F[0] = y[1]\n", " F[1] = -(0.4829)*(E-V(x))*y[0]\n", "\n", "def V(x):\n", " if (abs(x) < 10.): return (-16.0) # Well depth\n", " else: return (0.)\n", "\n", "def rk4(t, y,h,Neqs,E):\n", " F =np.zeros((Neqs),float)\n", " ydumb = np.zeros((Neqs),float)\n", " k1 = np.zeros((Neqs),float)\n", " k2 = np.zeros((Neqs),float)\n", " k3 = np.zeros((Neqs),float)\n", " k4 = np.zeros((Neqs),float)\n", " f(t, y, F,E)\n", " for i in range(0,Neqs):\n", " k1[i] = h*F[i]\n", " ydumb[i] = y[i] + k1[i]/2.\n", " f(t + h/2., ydumb, F,E)\n", " for i in range(0,Neqs):\n", " k2[i] = h*F[i]\n", " ydumb[i] = y[i] + k2[i]/2.\n", " f(t + h/2., ydumb, F,E)\n", " for i in range(0,Neqs):\n", " k3[i]= h*F[i]\n", " ydumb[i] = y[i] + k3[i]\n", " f(t + h, ydumb, F,E);\n", " for i in range(0,Neqs):\n", " k4[i]=h*F[i]\n", " y[i]=y[i]+(k1[i]+2*(k2[i]+k3[i])+k4[i])/6.0\n", "\n", "def diff(E, h):\n", " y = np.zeros((2),float)\n", " i_match = n_steps//3 # Matching radius\n", " nL = i_match + 1 \n", " y[0] = 1.E-15; # Initial left wf\n", " y[1] = y[0]*sqrt(-E*0.4829) \n", " for ix in range(0,nL + 1):\n", " x = h * (ix -n_steps/2)\n", " rk4(x, y, h, 2, E)\n", " left = y[1]/y[0] # Log derivative\n", " y[0] = 1.E-15; # For even; reverse for odd\n", " y[1] = -y[0]*sqrt(-E*0.4829) # Initialize R wf\n", " for ix in range( n_steps,nL+1,-1):\n", " x = h*(ix+1-n_steps/2)\n", " rk4(x,y, -h, 2, E)\n", " right = y[1]/y[0] # Log derivative\n", " return( (left - right)/(left + right) )\n", "\t\t\n", "def plot(E, h): # Repeat integrations for plot\n", " x = 0. \n", " n_steps = 1501 # # integration steps\n", " y = np.zeros((2),float)\n", " yL = np.zeros((2,505),float) \n", " i_match = 500 # Matching point\n", " nL = i_match + 1; \n", " y[0] = 1.E-40 # Initial left wf\n", " y[1] = -sqrt(-E*0.4829) *y[0]\n", " for ix in range(0,nL+1): \n", " yL[0][ix] = y[0]\n", " yL[1][ix] = y[1]\n", " x = h * (ix -n_steps/2)\n", " rk4(x, y, h, 2, E)\n", " y[0] = -1.E-15 # For even; reverse for odd\n", " y[1] = -sqrt(-E*0.4829)*y[0]\n", " j=0\n", " Rwf.clear() \n", " for ix in range(n_steps -1,nL + 2,-1): # right WF\n", " x = h * (ix + 1 -n_steps/2) # Integrate in\n", " rk4(x, y, -h, 2, E)\n", " Rwf.append(pos =vec( 2.*(ix + 1 -n_steps/2)-500.0, \n", " y[0]*35e-9 +200,0)) #plot new\n", " Rwf.visible=True #make visible new\n", " j +=1\n", " x = x-h \n", " normL = y[0]/yL[0][nL]\n", " j=0\n", " Lwf.clear() #erase previous plot left\n", " # Renormalize L wf & derivative\n", " for ix in range(0,nL+1):\n", " x = h * (ix-n_steps/2 + 1) \n", " y[0] = yL[0][ix]*normL \n", " y[1] = yL[1][ix]*normL\n", " Lwf.append(pos =vec( 2.*(ix -n_steps/2+1)-500.0,\n", " y[0]*35e-9+200,0)) # Factor for scale\n", " Lwf.visible=True\n", " j +=1\n", " \n", "for count in range(0,count_max+1):\n", " rate(1) # Slow rate to show changes\n", " # Iteration loop\n", " E = (Emax + Emin)/2. # Divide E range\n", " Diff = diff(E, h) \n", " if (diff(Emax, h)*Diff > 0): Emax = E # Bisection algor\n", " else: Emin = E \n", " if ( abs(Diff) < eps ): break\n", " if count >3: # First iterates too irregular\n", " rate(4)\n", " plot(E, h)\n", " elabel = label(pos=vec(700, 400,0), text='E=', box=0)\n", " elabel.text = 'E=%13.10f' %E\n", " ilabel = label(pos=vec(700, 600,0), text='istep=', box=0)\n", " ilabel.text = 'istep=%4s' %count\n", "elabel = label(pos=vec(700, 400,0), text='E=', box=0) # Last iteration\n", "elabel.text = 'E=%13.10f' %E\n", "ilabel = label(pos=vec(700, 600,0), text='istep=', box=0)\n", "ilabel.text = 'istep=%4s' %count\n", " \n", "print(\"Final eigenvalue E = \",E)\n", "print(\"iterations, max = \",count)" ] } ], "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 }