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

Bound States of the Momentum Space Schrödinger Equation

\n", "

\n", " The momentum space Schrödinger equation is solved for bound states by converting the integral equation into a set of linear algebraic equations:\n", " $$\n", "\\frac{k^{2}}{2\\mu} \\psi(k) + \\frac{2}{\\pi} \\int_{0}^{\\infty} dp\n", "p^{2} V(k,p) \\psi(p) = E \\psi(k), \\quad\n", "V(k,p) = \\frac{1}{kp} \\int_{0}^{\\infty} dr\\sin(kr) V(r) \\sin(pr),\n", "\\\\\n", "\\int_{0}^{\\infty} dp p^{2} V(k,p) \\psi(p) \\simeq \\sum_{j=1}^{N}\n", "w_jk_{j}^{2} V(k,k_j) \\psi(k_j) \\quad \\Rightarrow \\quad\n", "\\frac{k^{2}}{2\\mu} \\psi(k) + \\frac{2}{\\pi} \\sum_{j=1}^{N}w_j\n", "k_{j}^{2} V(k,k_j) \\psi(k_j) =E. \\label{BS.pse3}\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" } ], "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 (deceases), Univ Bucharest, 2020. \n", " Please respect copyright & acknowledge our work.\"\"\"\n", "\n", "# Bound.py: Bound state solutn of Lippmann-Schwinger equation in p space\n", "\n", "# from vpython import *\n", "import numpy as np\n", "from numpy.linalg import*\n", "\n", "min1 =0.; max1 =200.; u =0.5; b =10." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " M (size), lmbda, ReE = 16 -1024 -49908.19098862982\n", " M (size), lmbda, ReE = 16 -512 -22380.033477865014\n", " M (size), lmbda, ReE = 16 -256 -9284.341736932587\n", " M (size), lmbda, ReE = 16 -128 -3411.003181753296\n", " M (size), lmbda, ReE = 16 -64 -1075.959919874758\n", " M (size), lmbda, ReE = 16 -32 -319.43524496556023\n", " M (size), lmbda, ReE = 16 -16 -99.29251195626384\n", " M (size), lmbda, ReE = 16 -8 -30.86088382857418\n", " M (size), lmbda, ReE = 16 -4 -9.294986176821684\n", " M (size), lmbda, ReE = 16 -2 -2.7985328823847087\n", " M (size), lmbda, ReE = 24 -1024 -51290.69884814227\n", " M (size), lmbda, ReE = 24 -512 -21071.009025313448\n", " M (size), lmbda, ReE = 24 -256 -7609.78443706071\n", " M (size), lmbda, ReE = 24 -128 -2425.155449403763\n", " M (size), lmbda, ReE = 24 -64 -717.997555659577\n", " M (size), lmbda, ReE = 24 -32 -202.00811243397357\n", " M (size), lmbda, ReE = 24 -16 -53.92762488480716\n", " M (size), lmbda, ReE = 24 -8 -14.041875630397037\n", " M (size), lmbda, ReE = 24 -4 -4.413758333618474\n", " M (size), lmbda, ReE = 24 -2 -1.6542753342186143\n" ] } ], "source": [ "def gauss(npts,a,b,x,w):\n", " pp = 0.; m = (npts + 1)//2; eps = 3.E-10 # Accuracy: ADJUST!\n", " for i in range(1,m+1):\n", " t = cos(np.pi*(float(i)-0.25)/(float(npts) + 0.5))\n", " t1 = 1 \n", " while((abs(t-t1)) >= eps):\n", " p1 = 1. ; p2 = 0.; \n", " for j in range(1,npts+1):\n", " p3 = p2\n", " p2 = p1 \n", " p1=((2*j-1)*t*p2-(j-1)*p3)/j\n", " pp = npts*(t*p1-p2)/(t*t-1.) \n", " t1 = t; t = t1 - p1/pp \n", " x[i-1] = -t\n", " x[npts-i] = t \n", " w[i-1] = 2./((1.-t*t)*pp*pp) \n", " w[npts-i] = w[i-1] \n", " for i in range(0,npts):\n", " x[i] = x[i]*(b-a)/2. + (b + a)/2. \n", " w[i] = w[i]*(b-a)/2.\n", "\n", "for M in range(16, 32, 8):\n", " z=[-1024, -512, -256, -128, -64, -32, -16, -8, -4, -2]\n", " for lmbda in z:\n", " A = np.zeros((M,M), float) # Hamiltonian\n", " WR = np.zeros((M), float) # Eigenvalues, potential\n", " k = np.zeros((M), float); w = np.zeros((M),float); # Pts & wts\n", " gauss(M, min1, max1, k, w) # Call gauss points\n", " for i in range(0,M): # Set up Hamiltonian\n", " for j in range(0,M):\n", " VR = lmbda/2/u*sin(k[i]*b)/k[i]*sin(k[j]*b)/k[j]\n", " A[i,j] = 2./np.pi*VR*k[j]*k[j]*w[j] \n", " if (i == j): A[i,j] += k[i]*k[i]/2/u \n", " Es, evectors = eig(A) \n", " realev = Es.real # Real eigenvalues\n", " for j in range(0,M):\n", " if (realev[j]<0):\n", " print(\" M (size), lmbda, ReE = \",M,\" \",lmbda,\" \",realev[j])\n", " break" ] } ], "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 }