{ "metadata": { "name": "", "signature": "sha256:dfc4bc93a3dfa6e710f000ae8ca8389f0091c899d5820fe703a32c54f16d8113" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt, mpld3\n", "import matplotlib.cm as cm\n", "from numpy import sqrt, exp, sinh\n", "plt.rcParams['figure.figsize'] = 12, 9\n", "mpld3.enable_notebook()" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 43 }, { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Numerov's Method" ] }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Single Step" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def numerov_step(k0, k1, k2, psi0, psi1, dx):\n", " \"\"\"compute psi2 via a single numerov step\"\"\"\n", " dx_sqr = dx**2\n", " c0 = (1 + 1/12. * dx_sqr * k0)\n", " c1 = 2 * (1 - 5/12. * dx_sqr * k1)\n", " c2 = (1 + 1/12. * dx_sqr * k2)\n", " return (c1 * psi1 - c0 * psi0) / c2" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 44 }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a basic test, let's check the \"symmetry property\" of a single Numerov step:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "assert np.isclose(numerov_step(2, 4, 7, numerov_step(7, 4, 2, 2, 4, .25), 4, .25), 2)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 45 }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Low-Level Numerov" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def numerov_iter(ks, psi0, psi1, dx):\n", " \"\"\"compute psis = [psi0, psi1, psi2, ...] for ks = [k0, k1, ...] via iterated numerov steps\"\"\"\n", " n = len(ks)\n", " psis = np.zeros(n, dtype=np.complex128)\n", " psis[0] = psi0\n", " psis[1] = psi1\n", " for i in range(2, n):\n", " psis[i] = numerov_step(ks[i-2], ks[i-1], ks[i], psis[i-2], psis[i-1], dx)\n", " return psis" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 46 }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "High-Level Numerov" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def numerov(k, psi0, psi1, x_min, x_max, n):\n", " \"\"\"compute psis = [psi0, psi1, ...] for k = k(x) according to given discretization\"\"\"\n", " xs, dx = np.linspace(x_min, x_max, num=n, retstep=True)\n", " ks = np.vectorize(k)(xs)\n", " return xs, numerov_iter(ks, psi0, psi1, dx)\n", "\n", "def numerov_right_to_left(k, psi_2ndlast, psi_last, x_min, x_max, n):\n", " \"\"\"compute psis = [..., psi_2ndlast, psi_last] for k = k(x) according to given discretization\"\"\"\n", " xs, dx = np.linspace(x_min, x_max, num=n, retstep=True)\n", " ks = np.vectorize(k)(xs)\n", " return xs, numerov_iter(ks[::-1], psi_last, psi_2ndlast, dx)[::-1]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 47 }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Scattering Coefficients" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def scatter(V, E, x_min, x_max, n, return_xs_psis=False):\n", " \"\"\"compute transmission & reflection coefficient for given potential and energy (and m = hbar = 1)\"\"\"\n", " # compute step size\n", " dx = (x_max - x_min) / float(n)\n", "\n", " # start with right-travling wave at energy E (w/ phase fixed to 1 at x_max)\n", " q_free = sqrt(2. * E)\n", " psi_2nd_last = 1\n", " psi_last = exp(1j * q_free * dx)\n", "\n", " # compute resulting solution using Numerov\n", " def k(x):\n", " return 2. * (E - V(x))\n", " xs, psis = numerov_right_to_left(k, psi_2nd_last, psi_last, x_min - dx, x_max + dx, n + 2)\n", " \n", " # fit psis[] to A * exp(iqx) + B * exp(-iqx)\n", " A, B = np.linalg.solve(\n", " [[exp(1j*q_free*(-dx)), exp(-1j*q_free*(-dx))],\n", " [1, 1]],\n", " [psis[0], psis[1]])\n", " \n", " # extract transmission and reflection coefficients\n", " T = 1 / abs(A)**2\n", " R = abs(B)**2 / abs(A)**2\n", " return (T, R, xs, psis) if return_xs_psis else (T, R)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 48 }, { "cell_type": "markdown", "metadata": {}, "source": [ "\"Test\" with a free particle:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "assert np.allclose(scatter(V=lambda _: 0, E=1, x_min=0, x_max=1, n=1e5), [1, 0])" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 49 }, { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Rectangular Potential" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def rect_potential(a):\n", " def V(x):\n", " return float(0 < x <= a)\n", " return V\n", "\n", "xs = np.linspace(-0.5, 1.5, num=100)\n", "plt.figure(figsize=(4, 3))\n", "plt.plot(xs, [rect_potential(a=1)(x) for x in xs])\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", "\n", "\n", "\n", "
\n", "" ], "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAQcAAADICAYAAADlYiRyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEuZJREFUeJzt3V9sk+XiB/Bvob0AQmAI7oy2yWQtpxNYmSlOYmbqnzmY\nh4p/LuoVgWXuEBZiYgyJXrgZJcybc2EvGAYR/LPMC+M4x1EVQyWCox5AMBnBbrrYNbo4YIEjF2Pd\n+7tY7I/xPm87+7xd+7x+PwlxdQ/vnvpu337f531cbZqmaSAiusO8Yk+AiEoTw4GIhBgORCTEcCAi\nIYYDEQkxHIhIKGc47NixA+Xl5Vi3bp3hmN27d8Pr9cLv9+P8+fOmTpCIiiNnOGzfvh3RaNTw8319\nfRgcHEQikcCBAwewc+dOUydIRMWRMxzq6+tRVlZm+PmjR49i27ZtAIC6ujqMj49jdHTUvBkSUVFI\nrzmkUim43e7MY5fLhZGREdnDElGR2c04yJ07sG02m26Mx+PB0NCQGV+OiP6EqqoqDA4O/um/J90c\nnE4nkslk5vHIyAicTqdu3NDQEDRNs+SfeFzD3/72atHnUcg/r75qzef3xBMawmFrPrc//uT7oiwd\nDqFQCEeOHAEA9Pf3Y+nSpSgvL5c9rFLSaWAebworyW4HNP6vh0I5Lyuee+45fPXVVxgbG4Pb7UZH\nRwdu3boFAGhtbUVTUxP6+vrg8XiwaNEiHDp0qOCTLjUMB3XNnw9MTRV7FqUpZzh0d3fnPEgkEjFl\nMqqanATKyoLFnkZBBYPBYk+hIOx24O9/DxZ7GiWJr3cmSKeBu+4KFnsaBWXVcJg/H/D5gsWeRkli\nOJhgcnL6FYjUY7dPnz/SYziYIJ2efgUi9cyfP33+SI/hYILJSYaDqubPZ3MwwnAwQTrNywpV2e1s\nDkYYDiZgc1AXm4MxhoMJ2BzUxeZgjOFgAi5IqovNwRjDwQS8lakuNgdjDAcTsDmoi7cyjTEcTMDm\noC5ugjLGcDABm4O62ByMMRxMwOagLi5IGmM4mIDNQV1ckDTGcDABN0Gpi83BGMPBBNwEpS42B2MM\nBxPwskJdbA7GGA4m4IKkutgcjDEcTMDmoC42B2MMBxOwOaiLzcEYw8EEbA7q4iYoYwwHE7A5qIuX\nFcYYDiZgc1AXLyuMMRxMwE1Q6mJzMMZwMAE3QamLzcFYznCIRqPw+Xzwer3o7OzUfX5sbAybNm3C\n+vXrsXbtWrz77ruFmGdJ42WFutgcjGUNh3Q6jba2NkSjUQwMDKC7uxuXLl2aMSYSiaC2thbfffcd\nYrEYXnzxRUz+xf5rc0FSXWwOxrKGQzweh8fjQWVlJRwOB8LhMHp7e2eMqaiowPXr1wEA169fx113\n3QX7X+wnhc1BXWwOxrL+FKdSKbjd7sxjl8uFM2fOzBjT0tKCRx55BCtXrsSNGzfw0UcfFWamJYzN\nQV1sDsayfkvbbLacB9i7dy/Wr1+PWCyGoaEhNDQ04MKFC1i8eLFubHt7e+bjYDBomTdnZXNQlxU3\nQcViMcRiMenjZA0Hp9OJZDKZeZxMJuFyuWaMOX36NF555RUAQFVVFe655x5cvnwZgUBAd7zbw8FK\n2BzUZcXLijtfeDs6OvI6TtY1h0AggEQigeHhYUxMTKCnpwehUGjGGJ/Ph+PHjwMARkdHcfnyZaxa\ntSqvyaiKzUFdvKwwlvX1zm63IxKJoLGxEel0Gs3NzaiurkZXVxcAoLW1FS+//DK2b98Ov9+Pqakp\nvPnmm1i2bNmcTL5UsDmoy4rNwSw2TdO0OflCNhvm6EvNuX/8A2htBbZsKfZM6M/673+Bf/5z+p9W\nle/PHndImoCXFepiczDGcDABLyvUxTUHYwwHE7A5qIvNwRjDwQRsDupiczDGcDABm4O62ByMMRxM\nwOagLivukDQLw8EEbA7q4mWFMYaDCdgc1MXLCmMMBxOwOaiLzcEYw8EE/B2S6mJzMMZwMAF/h6S6\n2ByMMRxMwMsKdbE5GGM4mIALkupiczDGcDABm4O62ByMMRxMwOagrvnzgakpwKK/TUAKw8EEbA7q\nstmAefOmA4JmYjiYgM1Bbby0EGM4mIDNQW1clBRjOJiAzUFtbA5iDAcTsDmojc1BjOFgAoaD2tgc\nxBgOkjSN4aA6NgcxhoOkqan/vx1GamJzEOO3tCQuRqqPvw1KLGc4RKNR+Hw+eL1edHZ2CsfEYjHU\n1tZi7dq1lnlz3NniJYX67HY2B5Gsr3npdBptbW04fvw4nE4nNmzYgFAohOrq6syY8fFx7Nq1C599\n9hlcLhfGxsYKPulSwuagPjYHsazNIR6Pw+PxoLKyEg6HA+FwGL29vTPGfPjhh3jmmWcy7769fPny\nws22BLE5qI8LkmJZwyGVSsHtdmceu1wupFKpGWMSiQSuXr2Khx9+GIFAAO+9915hZlqi2BzUxwVJ\nsazf1jabLecBbt26hXPnzuHLL7/EzZs3sXHjRjzwwAPwer26se3t7ZmPg8GgJdYn2BzUZ7XmEIvF\nEIvFpI+TNRycTieSyWTmcTKZzFw+/MHtdmP58uVYsGABFixYgIceeggXLlzIGQ5WwV8Rpz6rNYc7\nX3g7OjryOk7Wy4pAIIBEIoHh4WFMTEygp6cHoVBoxpgnn3wSX3/9NdLpNG7evIkzZ87g3nvvzWsy\nKuIvl1Wf1ZqDWbK+5tntdkQiETQ2NiKdTqO5uRnV1dXo6uoCALS2tsLn82HTpk2oqanBvHnz0NLS\n8pcKB15WqM9qzcEsNk2bm9+BY7PZMEdfak798APwxBNAIlHsmVC+6uuBN94AHnqo2DMpjHx/9rhD\nUhKbg/q4CUqM4SCJtzLVx01QYgwHSWwO6mNzEGM4SGJzUB+bgxjDQRKbg/p4K1OM4SCJm6DUx1uZ\nYgwHSdwEpT42BzGGgyReVqiPzUGM4SCJC5Lq44KkGMNBEpuD+ngrU4zhIInNQX1sDmIMB0lsDupj\ncxBjOEhic1Afm4MYw0ESm4P6eCtTjOEgiZug1MdbmWIMB0ncBKU+NgcxhoMkNgf1sTmIMRwksTmo\njwuSYgwHSVyQVB9vZYoxHCTxVqb62BzEGA6S2BzUx+YgxnCQxOagPjYHMYaDJDYH9bE5iDEcJPFW\npvrYHMQYDpJ4K1N93AQlljMcotEofD4fvF4vOjs7Dcd9++23sNvt+Pjjj02dYKljc1AfN0GJZQ2H\ndDqNtrY2RKNRDAwMoLu7G5cuXRKO27NnDzZt2mTJt7zLhs1BfbysEMsaDvF4HB6PB5WVlXA4HAiH\nw+jt7dWNe+utt/Dss89ixYoVBZtoqeKCpPq4ICmWNRxSqRTcbnfmscvlQiqV0o3p7e3Fzp07AUy/\naedfCW9lqo/NQSzrt/VsftBfeOEF7Nu3L/NOvtkuK9rb2zMfB4NBBIPBWU+0VLE5qM9qzSEWiyEW\ni0kfJ2s4OJ1OJJPJzONkMgmXyzVjzNmzZxEOhwEAY2NjOHbsGBwOB0KhkO54t4eDVbA5qM9qzeHO\nF96Ojo68jpP12zoQCCCRSGB4eBgrV65ET08Puru7Z4z58ccfMx9v374dW7ZsEQaDVbE5qM9qzcEs\nWcPBbrcjEomgsbER6XQazc3NqK6uRldXFwCgtbV1TiZZyngrU31Waw5msWlzdO/xjzUJq2lpATZs\nAJ5/vtgzoXz9+9/AgQPT/7SifH/2uENSEpuD+rgJSozhIImboNTHywoxhoMkNgf1cUFSjOEgic1B\nfWwOYgwHSbyVqT42BzGGgyRuglIfm4MYw0ESm4P62BzEGA6SuCCpPjYHMYaDJC5Iqo/NQYzhIInN\nQX1sDmIMB0lsDupjOIgxHCSxOaiPlxViDAdJbA7qY3MQYzhIYnNQH5uDGMNBEpuD+tgcxBgOkrgJ\nSn1sDmIMB0m8rFAfm4MYw0ESLyvUx+YgxnCQxOagPjYHMYaDJDYH9fHXxIkxHCSxOaiP77ItxnCQ\nxOagvj8uKyz4y9GlMBwksTmob948wGYDpqaKPZPSwnCQxOZgDVyU1JtVOESjUfh8Pni9XnR2duo+\n/8EHH8Dv96OmpgYPPvggLl68aPpESxU3QVkDb2fq5SzE6XQabW1tOH78OJxOJzZs2IBQKITq6urM\nmFWrVuHkyZNYsmQJotEonn/+efT39xd04qWClxXWwOagl7M5xONxeDweVFZWwuFwIBwOo7e3d8aY\njRs3YsmSJQCAuro6jIyMFGa2JYiXFdbA5qCXMxxSqRTcbnfmscvlQiqVMhx/8OBBNDU1mTM7BbA5\nWAObg17Ob2ubzTbrg504cQLvvPMOTp06Jfx8e3t75uNgMIhgMDjrY5cqNgdrsNJGqFgshlgsJn2c\nnOHgdDqRTCYzj5PJJFwul27cxYsX0dLSgmg0irKyMuGxbg8HK5iamr4FNo/3fJRnpY1Qd77wdnR0\n5HWcnN/WgUAAiUQCw8PDmJiYQE9PD0Kh0IwxP//8M55++mm8//778Hg8eU1ERWwN1sHLCr2czcFu\ntyMSiaCxsRHpdBrNzc2orq5GV1cXAKC1tRWvvfYarl27hp07dwIAHA4H4vF4YWdeArjeYB1ckNSz\nadrcbBq12WyYoy81Z27cACoqgP/9r9gzIVkeD3DsGOD1Fnsm5sv3Z49XyxLYHKyDzUGP4SCBuyOt\ng2sOegwHCVyQtA42Bz2GgwReVlgHm4Mew0ECm4N1WGkTlFkYDhLYHKzDSpugzMJwkMDmYB1sDnoM\nBwlsDtbB5qDHcJDA5mAdXJDUYzhIYHOwDt7K1GM4SOAmKOtgc9BjOEjgZYV1sDnoMRwk8LLCOtgc\n9BgOEtgcrIO3MvUYDhLYHKyDtzL1GA4S2Bysg81Bj+Eggc3BOtgc9BgOEtgcrIPNQY/hIIHNwTrY\nHPQYDhK4Cco6eCtTj+EgYXKSzcEquAlKj+Eggc3BOtgc9BgOErggaR1ckNRjOEjggqR1cEFSj+Eg\ngc3BOtgc9HKGQzQahc/ng9frRWdnp3DM7t274fV64ff7cf78edMnWapubw5mvKtxKbPy84vFYmwO\nAlnDIZ1Oo62tDdFoFAMDA+ju7salS5dmjOnr68Pg4CASiQQOHDiQeb/Mv4Lbm4OVf3gAaz+/WCzG\n5iCQNRzi8Tg8Hg8qKyvhcDgQDofR29s7Y8zRo0exbds2AEBdXR3Gx8cxOjpauBmXEK45WAebg17W\nb+1UKgW325157HK5cObMmZxjRkZGUF5erjveli2y0y0tP/4INDYWexZkBocD+Ogj4Pvviz2T0pE1\nHGw226wOcuc7+Ir+XlVVFf7zn9kdTyUDA8C//jX9cUdHR3EnU2DWfn7Tz21oqMjTKICqqqq8/l7W\ncHA6nUgmk5nHyWQSLpcr65iRkRE4nU7dsQYHB/OaIBEVR9Y1h0AggEQigeHhYUxMTKCnpwehUGjG\nmFAohCNHjgAA+vv7sXTpUuElBRGpJWtzsNvtiEQiaGxsRDqdRnNzM6qrq9HV1QUAaG1tRVNTE/r6\n+uDxeLBo0SIcOnRoTiZORIVl0+5cMCAiQgF3SF69ehUNDQ1YvXo1Hn/8cYyPjwvHVVZWoqamBrW1\ntbj//vsLNR1TWH1DWK7nF4vFsGTJEtTW1qK2thavv/56EWaZnx07dqC8vBzr1q0zHKPyucv1/PI6\nd1qBvPTSS1pnZ6emaZq2b98+bc+ePcJxlZWV2pUrVwo1DdNMTk5qVVVV2k8//aRNTExofr9fGxgY\nmDHm008/1TZv3qxpmqb19/drdXV1xZhqXmbz/E6cOKFt2bKlSDOUc/LkSe3cuXPa2rVrhZ9X+dxp\nWu7nl8+5K1hzuH1z1LZt2/DJJ59kC6hCTcM0Vt8QNpvnB6hxrkTq6+tRVlZm+HmVzx2Q+/kBf/7c\nFSwcRkdHM3ctysvLDf9D22w2PPbYYwgEAnj77bcLNR1pos1eqVQq55iRkZE5m6OM2Tw/m82G06dP\nw+/3o6mpCQMDA3M9zYJR+dzNRj7nTmrzb0NDA3799Vfdv3/jjTd0EzPaUHXq1ClUVFTgt99+Q0ND\nA3w+H+rr62WmVRBmbggrRbOZ53333YdkMomFCxfi2LFj2Lp1K3744Yc5mN3cUPXczUY+506qOXzx\nxRf4/vvvdX9CoRDKy8szwfHLL7/g7rvvFh6joqICALBixQo89dRTiMfjMlMqGDM3hJWi2Ty/xYsX\nY+HChQCAzZs349atW7h69eqczrNQVD53s5HPuSvYZUUoFMLhw4cBAIcPH8bWrVt1Y27evIkbN24A\nAH7//Xd8/vnnWVeTi8nqG8Jm8/xGR0czr67xeByapmHZsmXFmK7pVD53s5HXuZNYIM3qypUr2qOP\nPqp5vV6toaFBu3btmqZpmpZKpbSmpiZN0zRtaGhI8/v9mt/v19asWaPt3bu3UNMxRV9fn7Z69Wqt\nqqoqM9f9+/dr+/fvz4zZtWuXVlVVpdXU1Ghnz54t1lTzkuv5RSIRbc2aNZrf79c2btyoffPNN8Wc\n7p8SDoe1iooKzeFwaC6XSzt48KClzl2u55fPueMmKCIS4q+JIyIhhgMRCTEciEiI4UBEQgwHIhJi\nOBCREMOBiIT+Dx/NA7Yxm2gtAAAAAElFTkSuQmCC\n", "text": [ "