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

11: Solve some TODOs, do some TeX

Browse Source
This commit is contained in:
Koen Vendrig
2022-03-18 09:39:48 +01:00
committed by Kees van Kempen
parent f7ac2c3ff7
commit 1d7dec786c
1 changed files with 55 additions and 16 deletions
Download Patch File Download Diff File Expand all files Collapse all files

71
Week 6/11 Parabolic PDEs: 1D Schrödinger Equation with Potential.ipynb Normal file → Executable file
View File

@ -160,12 +160,13 @@
"\n", "\n",
"with $\\delta_{ij} = 1 \\iff i = j$, and using the approximations $\\hbar = 1$, $m = 0.5$.\n", "with $\\delta_{ij} = 1 \\iff i = j$, and using the approximations $\\hbar = 1$, $m = 0.5$.\n",
"\n", "\n",
"\n", "$$\n",
"finite differences for second derivative to get a tridiagonal matrix (one lower and one higher than diagonal)\n", "%finite differences for second derivative to get a tridiagonal matrix (one lower and one higher than diagonal)\n",
"then plot eigenvalues\n", "%then plot eigenvalues\n",
"no time component -> not complex?\n", "%no time component -> not complex?\n",
"energy should be real\n", "%energy should be real\n",
"eigenvectors squared should be normalized" "%eigenvectors squared should be normalized\n",
"$$"
] ]
}, },
{ {
@ -201,7 +202,6 @@
" A tuple of the eigenvalues E and eigenfunctions \\Psi of the hamiltonian H.\n", " A tuple of the eigenvalues E and eigenfunctions \\Psi of the hamiltonian H.\n",
" \"\"\"\n", " \"\"\"\n",
" \n", " \n",
" # TODO: Is n the dimensionality of Psi?\n",
" n = len(x)\n", " n = len(x)\n",
" h = x[1] - x[0]\n", " h = x[1] - x[0]\n",
" \n", " \n",
@ -286,7 +286,6 @@
"fig = plt.figure(figsize=(12,12))\n", "fig = plt.figure(figsize=(12,12))\n",
"\n", "\n",
"for n in range(the_first_few):\n", "for n in range(the_first_few):\n",
" # TODO: Should SEQStat return normalized eigenfunctions?\n",
" # Normalize the found eigenfunctions using Riemann midpoint sums\n", " # Normalize the found eigenfunctions using Riemann midpoint sums\n",
" prob_mass = np.sum((x[1:] - x[:-1])*(Psi[1:, n]**2 + Psi[:-1, n]**2)/2)\n", " prob_mass = np.sum((x[1:] - x[:-1])*(Psi[1:, n]**2 + Psi[:-1, n]**2)/2)\n",
" Psi[:, n] /= np.sqrt(prob_mass)\n", " Psi[:, n] /= np.sqrt(prob_mass)\n",
@ -383,7 +382,50 @@
} }
}, },
"source": [ "source": [
"YOUR ANSWER HERE" "To do the derivation of the Crank-Nicolson approach following $\\textit{Numerical Analysis}$ by Burden and Faires, we start from averaging our 'Forward-Difference method' at the jth step in t:\n",
"\n",
"$$\n",
"\\frac{w_{i,j+1}-w_{i,j}}{k}-\\alpha^2\\frac{w_{i+1,j}-2w_{i,j}+w_{i-1,j}}{h^2} = 0.\n",
"$$\n",
"\n",
"The local truncation error is:\n",
"\n",
"$$\n",
"\\tau_F = \\frac{k}{2}\\frac{\\partial^2 u}{\\partial t^2}(x_i,\\mu_j)+O(h^2).\n",
"$$\n",
"\n",
"We do the same for the 'Backward-Difference method' to get:\n",
"\n",
"$$\n",
"\\frac{w_{i,j+1}-w_{i,j}}{k}-\\alpha^2\\frac{w_{i+1,j+1}-2w_{i,j+1}+w_{i-1,j+1}}{h^2} = 0,\n",
"$$\n",
"\n",
"again with a local truncation error:\n",
"\n",
"$$\n",
"\\tau_B = -\\frac{k}{2}\\frac{\\partial^2 u}{\\partial t^2}(x_i,\\hat{u} _j)+O(h^2).\n",
"$$\n",
"\n",
"If we make the following assumption:\n",
"\n",
"$$\n",
"\\frac{\\partial^2 u}{\\partial t^2}(x_i,\\hat{\\mu}_j) \\approx \\frac{\\partial^2 u}{\\partial t^2}(x_i,\\mu_j),\n",
"$$\n",
"\n",
"Then we can get the 'averaged-difference method' to write:\n",
"\n",
"$$\n",
"\\frac{w_{i,j+1}-w_{i,j}}{k}-\\frac{\\alpha^2}{2}\\left[\\frac{w_{i+1,j}-2w_{i,j}+w_{i-1,j}}{h^2}+\\frac{w_{i+1,j+1}-2w_{i,j+1}+w_{i-1,j+1}}{h^2}\\right] = 0,\n",
"$$\n",
"\n",
"which now has an truncation error of the order $O(k^2+h^2)$. This is known as the Crank-Nicolson method and is represented in the form described above.\n",
"\n",
"For $\\lambda_1$ and $\\lambda_2$ we have the following expressions:\n",
"\n",
"$$\n",
"\\lambda_1 = \\alpha^2\\frac{k}{h^2}\\\\\n",
"\\lambda_2 = 1.\n",
"$$"
] ]
}, },
{ {
@ -432,8 +474,7 @@
"outputs": [], "outputs": [],
"source": [ "source": [
"def SEQDyn(x, t, V, f):\n", "def SEQDyn(x, t, V, f):\n",
" # YOUR CODE HERE\n", " # YOUR CODE HERE"
" raise NotImplementedError()"
] ]
}, },
{ {
@ -457,8 +498,7 @@
"source": [ "source": [
"# Animate your solution here ...\n", "# Animate your solution here ...\n",
"\n", "\n",
"# YOUR CODE HERE\n", "# YOUR CODE HERE"
"raise NotImplementedError()"
] ]
}, },
{ {
@ -502,8 +542,7 @@
}, },
"outputs": [], "outputs": [],
"source": [ "source": [
"# YOUR CODE HERE\n", "# YOUR CODE HERE"
"raise NotImplementedError()"
] ]
}, },
{ {
@ -567,7 +606,7 @@
], ],
"metadata": { "metadata": {
"kernelspec": { "kernelspec": {
"display_name": "Python 3", "display_name": "Python 3 (ipykernel)",
"language": "python", "language": "python",
"name": "python3" "name": "python3"
}, },