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

Final: Implement better tests for integrate in task 2.2

Browse Source
This commit is contained in:
Kees van Kempen
2022-03-18 16:30:37 +01:00
parent d284f0c9a8
commit 11bcebe10d
1 changed files with 19 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files

46
Final/Final - Tight-binding propagation method.ipynb
View File

@ -308,8 +308,8 @@
"source": [ "source": [
"def integrate(yk, x):\n", "def integrate(yk, x):\n",
" \"\"\"\n", " \"\"\"\n",
" Numerically integrates function yk over [x[0], x[-1]] using\n", " Numerically integrates function yk over [x[0], x[-1]] using Simpson's 3/8\n",
" Simpson's rule over the grid provided by x.\n", " rule over the grid provided by x.\n",
" \n", " \n",
" Args:\n", " Args:\n",
" yk: function of one numerical argument that returns a numeric\n", " yk: function of one numerical argument that returns a numeric\n",
@ -326,23 +326,17 @@
" if callable(yk):\n", " if callable(yk):\n",
" yk = yk(x)\n", " yk = yk(x)\n",
" \n", " \n",
" #a, b = yk[0], yk[-1]\n", " # The distance h_i = x[i + 1] - x[i] is not necessarily constant. The choice of\n",
" \n", " # partitioning of the interval is subject to mathematical considerations I will\n",
" # not go into.\n",
" h = x[1:] - x[:-1]\n", " h = x[1:] - x[:-1]\n",
" \n", " \n",
" # Instead of summing over something with n/2, we use step size 2,\n",
" # thus avoiding any issues with 2 ∤ n.\n",
" #integral = h/3*(a + 2*np.sum(yk[2:-1:2]) + 4*np.sum(yk[1:-1:2]) + b)\n",
" \n",
" #integral = 3*h/8*(a + 3*np.sum(yk[2:-1:2]) + 3*np.sum(yk[1:-1:2]) + b)\n",
" #integral = (b - a)/8*(x[0] + x[-1]) + 3*(x[1:] - x[:-1])/8*( )\n",
" integral = 0\n", " integral = 0\n",
" integral += 3/8*(x[1] - x[0])*yk[0]\n", " integral += 3/8*(x[1] - x[0])*yk[0]\n",
" integral += 9/8*yk[1:-1:3]@h[1::3]\n", " integral += 9/8*h[1::3]@yk[1:-1:3]\n",
" integral += 9/8*yk[2:-1:3]@h[2::3]\n", " integral += 9/8*h[2::3]@yk[2:-1:3]\n",
" integral += 6/8*yk[:-1:3]@h[::3]\n", " integral += 6/8*h[ ::3]@yk[ :-1:3]\n",
" integral += 3/8*(x[-1] - x[-2])*yk[-1]\n", " integral += 3/8*(x[-1] - x[-2])*yk[-1]\n",
" #integral = 3/8*( (x[1] - x[0])*yk[0] + np.sum(yk[]) + (x[-1] - x[-2])*yk[-1])\n",
" return integral" " return integral"
] ]
}, },
@ -363,22 +357,20 @@
"task": false "task": false
} }
}, },
"outputs": [ "outputs": [],
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exact: 333333333.3333333\n",
"Integrated: 333333448.1291072\n"
]
}
],
"source": [ "source": [
"def test_integrate():\n", "def test_integrate():\n",
" # Test integral 1 of f with F its primitive with integration constant 0\n",
" f = lambda x: x**2\n", " f = lambda x: x**2\n",
" x = np.logspace(0, 3, 10000000)\n", " F = lambda x: x**3/3\n",
" print(\"Exact:\", 1000**3/3)\n", " x = np.logspace(0, 3, 1000000)\n",
" print(\"Integrated:\", integrate(f, x))\n", " assert np.isclose(integrate(f, x), F(x[-1]) - F(x[0]))\n",
" \n",
" # Test integral 2 of f with F its primitive with integration constant 0\n",
" f = lambda x: np.sin(2*x)/(2 + np.cos(2*x))\n",
" F = lambda x: -.5*np.log(np.cos(2*x) + 2)\n",
" x = np.linspace(0, 10, 1000)\n",
" assert np.isclose(integrate(f, x), F(x[-1]) - F(x[0]))\n",
" \n", " \n",
"test_integrate()" "test_integrate()"
] ]