From 5be960e2ef5adc3ddf730b550e2ad7f756ac3b5b Mon Sep 17 00:00:00 2001 From: Kees van Kempen Date: Tue, 15 Mar 2022 16:24:24 +0100 Subject: [PATCH] 11: Draft things. I dunno whether this is good at all. --- ...Schrödinger Equation with Potential.ipynb | 125 ++++++++++++++++-- 1 file changed, 111 insertions(+), 14 deletions(-) diff --git a/Week 6/11 Parabolic PDEs: 1D Schrödinger Equation with Potential.ipynb b/Week 6/11 Parabolic PDEs: 1D Schrödinger Equation with Potential.ipynb index 9d61ec1..4268513 100644 --- a/Week 6/11 Parabolic PDEs: 1D Schrödinger Equation with Potential.ipynb +++ b/Week 6/11 Parabolic PDEs: 1D Schrödinger Equation with Potential.ipynb @@ -39,7 +39,7 @@ "cell_type": "raw", "metadata": {}, "source": [ - "team_members = \"\"" + "team_members = \"Koen Vendrig, Kees van Kempen\"" ] }, { @@ -72,7 +72,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": { "deletable": false, "nbgrader": { @@ -89,10 +89,7 @@ }, "outputs": [], "source": [ - "# Import packages here ...\n", - "\n", - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "import numpy as np" ] }, { @@ -144,7 +141,24 @@ } }, "source": [ - "YOUR ANSWER HERE" + "The potential $V(x)$ is a scalar quantity, and will therefore only occur on the diagonal of the matrix $H$.\n", + "The second derivative operator will however be matrix valued, as $x$ is a vector, so that\n", + "\n", + "$$\n", + "\\partial^2/\\partial{}x^2 = \\begin{bmatrix} \\partial^2/\\partial{}x_1^2 & \\ldots & \\partial^2/\\partial{}x_1x_n \\\\ \\vdots & \\ddots & \\vdots \\\\ \\partial^2/\\partial{}x_nx_1 & \\ldots & \\partial^2/\\partial{}x_n^2 \\end{bmatrix},\n", + "$$\n", + "\n", + "or more cleanly, that element $ij$ equals\n", + "$$\n", + "(\\partial^2/\\partial{}x^2)_{ij} = \\partial^2/\\partial{}x_ix_j.\n", + "$$\n", + "\n", + "This leads to a hamiltonian matrix with elements\n", + "$$\n", + "H_{ij} = \\frac{-\\hbar^2}{2m}\\frac{\\partial^2}{\\partial{}x_ix_j} + V(x) \\delta_{ij} = -\\frac{\\partial^2}{\\partial{}x_ix_j} + V(x) \\delta_{ij},\n", + "$$\n", + "\n", + "with $\\delta_{ij} = 1 \\iff i = j$, and using the approximations $\\hbar = 1$, $m = 0.5$." ] }, { @@ -167,8 +181,62 @@ "outputs": [], "source": [ "def SEQStat(x, V):\n", - " # YOUR CODE HERE\n", - " raise NotImplementedError()" + " \"\"\"\n", + " Numerically solves the stationary Schrödinger equation over\n", + " the grid of x values, in a potential given by V, using\n", + " the approximations ℏ = 1, 𝑚 = 0.5.\n", + " \n", + " Args:\n", + " x: array of evenly spaced space value vectors x\n", + " V: function that takes in a location x and returns a scalar potential value\n", + "\n", + " Returns:\n", + " A tuple of the eigenvalues E and eigenfunctions \\Psi of the hamiltonian H.\n", + " \"\"\"\n", + " \n", + " n = len(x)\n", + " h = x[1] - x[0]\n", + " \n", + " # In the stationary case, \n", + " H = -1/h**2*np.ones(n) + np.eye(n)*V(x)\n", + " \n", + " # TODO: I think H could be complex, so I used eigh instead of eig.\n", + " return np.linalg.eigh(H)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[-2., -1., 0.],\n", + " [-1., -2., -1.],\n", + " [ 0., -1., -2.]])" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def stationary_hamiltonian(x, V, h, m): #construct the hamiltonian in tridiagonal matrix form for the stationary case (time-independent schrodinger equation) using the amount of x-steps m and the small x-parameter h (since hbar = 1 and m = 0.5 --> 2m = 1 all factors due to hbar or m will disappear)\n", + " H = np.zeros((m,m))\n", + " off_diag = -1/(h**2) * np.ones(m-1)\n", + " diag = [V[i] - 2*h**2 for i in range(m)]\n", + " np.fill_diagonal(H[1:,:], off_diag)\n", + " np.fill_diagonal(H[:,:], diag)\n", + " np.fill_diagonal(H[:,1:], off_diag)\n", + " return H\n", + "\n", + "V = lambda x: -x**2*0\n", + "x = np.array([1,2,3])\n", + "h = 1\n", + "m = 3\n", + "stationary_hamiltonian(x, V(x), h, m)" ] }, { @@ -195,7 +263,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": { "deletable": false, "nbgrader": { @@ -210,10 +278,39 @@ "task": false } }, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "array([-5. , -4.8989899 , -4.7979798 , -4.6969697 , -4.5959596 ,\n", + " -4.49494949, -4.39393939, -4.29292929, -4.19191919, -4.09090909,\n", + " -3.98989899, -3.88888889, -3.78787879, -3.68686869, -3.58585859,\n", + " -3.48484848, -3.38383838, -3.28282828, -3.18181818, -3.08080808,\n", + " -2.97979798, -2.87878788, -2.77777778, -2.67676768, -2.57575758,\n", + " -2.47474747, -2.37373737, -2.27272727, -2.17171717, -2.07070707,\n", + " -1.96969697, -1.86868687, -1.76767677, -1.66666667, -1.56565657,\n", + " -1.46464646, -1.36363636, -1.26262626, -1.16161616, -1.06060606,\n", + " -0.95959596, -0.85858586, -0.75757576, -0.65656566, -0.55555556,\n", + " -0.45454545, -0.35353535, -0.25252525, -0.15151515, -0.05050505,\n", + " 0.05050505, 0.15151515, 0.25252525, 0.35353535, 0.45454545,\n", + " 0.55555556, 0.65656566, 0.75757576, 0.85858586, 0.95959596,\n", + " 1.06060606, 1.16161616, 1.26262626, 1.36363636, 1.46464646,\n", + " 1.56565657, 1.66666667, 1.76767677, 1.86868687, 1.96969697,\n", + " 2.07070707, 2.17171717, 2.27272727, 2.37373737, 2.47474747,\n", + " 2.57575758, 2.67676768, 2.77777778, 2.87878788, 2.97979798,\n", + " 3.08080808, 3.18181818, 3.28282828, 3.38383838, 3.48484848,\n", + " 3.58585859, 3.68686869, 3.78787879, 3.88888889, 3.98989899,\n", + " 4.09090909, 4.19191919, 4.29292929, 4.39393939, 4.49494949,\n", + " 4.5959596 , 4.6969697 , 4.7979798 , 4.8989899 , 5. ])" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "x = np.linspace(-5, 5, 100)\n" ] }, { @@ -467,7 +564,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3", "language": "python", "name": "python3" },