Final: draft integration bullshite

Final/Final - Tight-binding propagation method.ipynb

@@ -289,7 +289,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 4,
    "metadata": {
     "deletable": false,
     "nbgrader": {
@@ -306,12 +306,43 @@
    },
    "outputs": [],
    "source": [
-    "def integrate(yk, x):\n",
-    "    # YOUR CODE HERE\n",
-    "    raise NotImplementedError()\n",
-    "#https://en.wikipedia.org/wiki/Numerical_integration#Quadrature_rules_based_on_interpolating_functions\n",
-    "#https://www.researchgate.net/post/Is_there_any_method_better_than_Gauss_Quadrature_for_numerical_integration\n",
-    "# \"Make sure to add all the docstrings and comments as to not lose points.\""
+    "def integrate(yk, x, b=None, c=None):\n",
+    "    \"\"\"\n",
+    "    Numerically integrates function yk over [x[0], x[-1]] using\n",
+    "    Simpson's rule over the grid provided by x.\n",
+    "    \n",
+    "    Args:\n",
+    "        yk: function of one numerical argument that returns a numeric\n",
+    "            or an array of function values such that x[i] corresponds to yk[i]\n",
+    "        x:  array of numerics as argument to yk\n",
+    "\n",
+    "    Returns:\n",
+    "        A numeric value for the quadrature of yk over x with error\n",
+    "        of order \n",
+    "    \"\"\"\n",
+    "    \n",
+    "    # If yk is callable, we use it to determine the function values\n",
+    "    # over array x.\n",
+    "    if callable(yk):\n",
+    "        yk = yk(x)\n",
+    "    \n",
+    "    a, b = yk[0], yk[-1]\n",
+    "    \n",
+    "    h = x[1:] - x[:-1]\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 += 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 = 3/8*( (x[1] - x[0])*yk[0] + np.sum(yk[]) + (x[-1] - x[-2])*yk[-1])\n",
+    "    return integral"
    ]
   },
   {