diff --git a/Week 2/6 Composite Numerical Integration: Trapezoid and Simpson Rules.ipynb b/Week 2/6 Composite Numerical Integration: Trapezoid and Simpson Rules.ipynb index a83c9cd..3dea41d 100644 --- a/Week 2/6 Composite Numerical Integration: Trapezoid and Simpson Rules.ipynb +++ b/Week 2/6 Composite Numerical Integration: Trapezoid and Simpson Rules.ipynb @@ -39,7 +39,7 @@ "cell_type": "raw", "metadata": {}, "source": [ - "team_members = \"\"" + "team_members = \"Koen Vendrig, Kees van Kempen\"" ] }, { @@ -75,7 +75,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "metadata": { "deletable": false, "nbgrader": { @@ -92,10 +92,8 @@ }, "outputs": [], "source": [ - "# Import packages here ...\n", - "\n", - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "import numpy as np\n", + "import scipy.integrate" ] }, { @@ -126,7 +124,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": { "deletable": false, "nbgrader": { @@ -144,17 +142,23 @@ "outputs": [], "source": [ "def trapz(yk, dx):\n", - " # YOUR CODE HERE\n", - " raise NotImplementedError()\n", + " a, b = yk[0], yk[-1]\n", + " h = dx\n", + " integral = h/2*(a + 2*np.sum(yk[1:-1]) + b)\n", + " return integral\n", " \n", "def simps(yk, dx):\n", - " # YOUR CODE HERE\n", - " raise NotImplementedError()" + " a, b = yk[0], yk[-1]\n", + " h = dx\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", + " return integral" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "metadata": { "deletable": false, "nbgrader": { @@ -171,10 +175,41 @@ }, "outputs": [], "source": [ - "# Compare your results here ...\n", + "# We need a function to integrate, so here we go.\n", + "f = lambda x: x**2\n", "\n", - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "n = 100001\n", + "a, b = 0, 1\n", + "h = (b - a)/n\n", + "xk = np.linspace(a, b, n + 1)\n", + "yk = f(xk)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "For function f(x) = x^2\n", + "for boundaries a = 0 , b = 1 and steps n = 100001 the algorithms say:\n", + "trapezoid:\t\t 0.33333333334999976\n", + "Simpson:\t\t 0.3333300000666658\n", + "scipy.integrate.trapz:\t 0.33333333334999965\n", + "scipy.integrate.simps:\t 0.3333333333333335\n" + ] + } + ], + "source": [ + "print(\"For function f(x) = x^2\")\n", + "print(\"for boundaries a =\", a, \", b =\", b, \"and steps n =\", n, \"the algorithms say:\")\n", + "print(\"trapezoid:\\t\\t\", trapz(yk, h))\n", + "print(\"Simpson:\\t\\t\", simps(yk, h))\n", + "print(\"scipy.integrate.trapz:\\t\", scipy.integrate.trapz(yk, xk))\n", + "print(\"scipy.integrate.simps:\\t\", scipy.integrate.simps(yk, xk))" ] }, { @@ -291,7 +326,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3", "language": "python", "name": "python3" },