Final: Fix integrate(yk, x) and implement test

Final/Final - Tight-binding propagation method.ipynb

@@ -306,7 +306,7 @@
    },
    "outputs": [],
    "source": [
-    "def integrate(yk, x, b=None, c=None):\n",
+    "def integrate(yk, x):\n",
     "    \"\"\"\n",
     "    Numerically integrates function yk over [x[0], x[-1]] using\n",
     "    Simpson's rule over the grid provided by x.\n",
@@ -326,7 +326,7 @@
     "    if callable(yk):\n",
     "        yk = yk(x)\n",
     "    \n",
-    "    a, b = yk[0], yk[-1]\n",
+    "    #a, b = yk[0], yk[-1]\n",
     "    \n",
     "    h = x[1:] - x[:-1]\n",
     "    \n",
@@ -334,20 +334,21 @@
     "    # 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 = 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 += 3/8*(x[1] - x[0])*yk[0]\n",
-    "    integral += 3/8*np.sum(yk[])\n",
-    "    \n",
-    "    integral += 3/8*(x[-1] - x[-2])*yk[-1])\n",
+    "    integral += 9/8*yk[1:-1:3]@h[1::3]\n",
+    "    integral += 9/8*yk[2:-1:3]@h[2::3]\n",
+    "    integral += 6/8*yk[:-1:3]@h[::3]\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"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 5,
    "metadata": {
     "deletable": false,
     "nbgrader": {
@@ -362,11 +363,22 @@
      "task": false
     }
    },
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Exact: 333333333.3333333\n",
+      "Integrated: 333333448.1291072\n"
+     ]
+    }
+   ],
    "source": [
     "def test_integrate():\n",
-    "    # YOUR CODE HERE\n",
-    "    raise NotImplementedError()\n",
+    "    f = lambda x: x**2\n",
+    "    x = np.logspace(0, 3, 10000000)\n",
+    "    print(\"Exact:\", 1000**3/3)\n",
+    "    print(\"Integrated:\", integrate(f, x))\n",
     "    \n",
     "test_integrate()"
    ]