kees-study/cds-numerical-methods
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

11: Draft things. I dunno whether this is good at all.

Browse Source
This commit is contained in:
Kees van Kempen
2022-03-15 16:24:24 +01:00
parent f56c5c26df
commit 5be960e2ef
1 changed files with 111 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

125
Week 6/11 Parabolic PDEs: 1D Schrödinger Equation with Potential.ipynb
View File

@ -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"
},