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 3dea41d..f48b0c5 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 @@ -75,7 +75,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 1, "metadata": { "deletable": false, "nbgrader": { @@ -93,7 +93,10 @@ "outputs": [], "source": [ "import numpy as np\n", - "import scipy.integrate" + "import scipy.integrate\n", + "\n", + "# And for printing the lambdas:\n", + "import inspect" ] }, { @@ -124,7 +127,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, "metadata": { "deletable": false, "nbgrader": { @@ -158,7 +161,31 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "def compare_integration(f, a, b, n):\n", + " # Let's check whether f is callable.\n", + " # TODO: Improve checks on f, or not.\n", + " assert callable(f)\n", + " \n", + " h = (b - a)/n\n", + " xk = np.linspace(a, b, n + 1)\n", + " yk = f(xk)\n", + " \n", + " print(\"For function\", inspect.getsource(f))\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))\n", + " print()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, "metadata": { "deletable": false, "nbgrader": { @@ -179,37 +206,40 @@ "f = lambda x: x**2\n", "\n", "n = 100001\n", - "a, b = 0, 1\n", - "h = (b - a)/n\n", - "xk = np.linspace(a, b, n + 1)\n", - "yk = f(xk)" + "a, b = 0, 1" ] }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "For function f(x) = x^2\n", + "For function f = lambda x: x**2\n", + "\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" + "scipy.integrate.simps:\t 0.3333333333333335\n", + "\n", + "For function f = lambda x: x**2\n", + "\n", + "for boundaries a = 0 , b = 1 and steps n = 100002 the algorithms say:\n", + "trapezoid:\t\t 0.3333333333499994\n", + "Simpson:\t\t 0.33333333333333337\n", + "scipy.integrate.trapz:\t 0.3333333333499993\n", + "scipy.integrate.simps:\t 0.3333333333333333\n", + "\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))" + "compare_integration(f, a, b, n)\n", + "compare_integration(f, a, b, n + 1)" ] }, { @@ -236,7 +266,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": { "deletable": false, "nbgrader": { @@ -251,15 +281,67 @@ "task": false } }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "For function f(x) = x^3 + 6x\n", + "for boundaries a = 2 , b = 16 and steps n = 82198 the algorithms say:\n", + "trapezoid:\t\t 1362.6666667343538\n", + "scipy.integrate.trapz:\t 1362.6666667343543\n", + "with difference trapz(yk, h) - scipy.integrate.trapz(yk, xk) = -4.547473508864641e-13\n", + "\n", + "For function f(x) = -x^3 + 6x\n", + "for boundaries a = 2 , b = 17 and steps n = 82228 the algorithms say:\n", + "Simpson:\t\t 1635.0\n", + "scipy.integrate.simps:\t 1635.0000000000002\n", + "with difference simps(yk, h) - scipy.integrate.simps(yk, xk) = -2.2737367544323206e-13\n", + "\n" + ] + } + ], "source": [ + "# In the comparison of n even and n odd, and the testing of the integrations,\n", + "# we have already tested the functions, but as it is asked, here we go again.\n", + "\n", "def test_trapz():\n", - " # YOUR CODE HERE\n", - " raise NotImplementedError()\n", + " fun = lambda x: x**3 + 6*x\n", + " a, b = 2, 16\n", + " n = 82198\n", + " \n", + " h = (b - a)/n\n", + " xk = np.linspace(a, b, n + 1)\n", + " yk = f(xk)\n", + "\n", + " trapz_our = trapz(yk, h)\n", + " trapz_scipy = scipy.integrate.trapz(yk, xk)\n", + " \n", + " print(\"For function f(x) = x^3 + 6x\")\n", + " print(\"for boundaries a =\", a, \", b =\", b, \"and steps n =\", n, \"the algorithms say:\")\n", + " print(\"trapezoid:\\t\\t\", trapz_our)\n", + " print(\"scipy.integrate.trapz:\\t\", trapz_scipy)\n", + " print(\"with difference trapz(yk, h) - scipy.integrate.trapz(yk, xk) =\", trapz_our - trapz_scipy)\n", + " print()\n", " \n", "def test_simps():\n", - " # YOUR CODE HERE\n", - " raise NotImplementedError()\n", + " fun = lambda x: -x**3 + 6*x\n", + " a, b = 2, 17\n", + " n = 82228\n", + " \n", + " h = (b - a)/n\n", + " xk = np.linspace(a, b, n + 1)\n", + " yk = f(xk)\n", + "\n", + " simps_our = simps(yk, h)\n", + " simps_scipy = scipy.integrate.simps(yk, xk)\n", + " \n", + " print(\"For function f(x) = -x^3 + 6x\")\n", + " print(\"for boundaries a =\", a, \", b =\", b, \"and steps n =\", n, \"the algorithms say:\")\n", + " print(\"Simpson:\\t\\t\", simps_our)\n", + " print(\"scipy.integrate.simps:\\t\", simps_scipy)\n", + " print(\"with difference simps(yk, h) - scipy.integrate.simps(yk, xk) =\", simps_our - simps_scipy)\n", + " print()\n", " \n", "test_trapz()\n", "test_simps()" @@ -295,7 +377,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": { "deletable": false, "nbgrader": { @@ -310,10 +392,60 @@ "task": false } }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "As we are calculus geniuses, we know that the first integral gives 1/2.\n", + "For function f1 = lambda x: x\n", + "\n", + "for boundaries a = 0 , b = 1 and steps n = 100000 the algorithms say:\n", + "trapezoid:\t\t 0.5000000000000001\n", + "Simpson:\t\t 0.5000000000000001\n", + "scipy.integrate.trapz:\t 0.5000000000000001\n", + "scipy.integrate.simps:\t 0.5\n", + "\n", + "As we are calculus geniuses, we know that the second integral gives 1/3.\n", + "For function f2 = lambda x: x**2\n", + "\n", + "for boundaries a = 0 , b = 1 and steps n = 100000 the algorithms say:\n", + "trapezoid:\t\t 0.33333333335000004\n", + "Simpson:\t\t 0.3333333333333335\n", + "scipy.integrate.trapz:\t 0.33333333335\n", + "scipy.integrate.simps:\t 0.3333333333333333\n", + "\n", + "As we are calculus geniuses, we know that the third integral gives 2/3.\n", + "For function f3 = lambda x: x**(1/2)\n", + "\n", + "for boundaries a = 0 , b = 1 and steps n = 100000 the algorithms say:\n", + "trapezoid:\t\t 0.6666666600968939\n", + "Simpson:\t\t 0.6666666640993837\n", + "scipy.integrate.trapz:\t 0.6666666600968938\n", + "scipy.integrate.simps:\t 0.6666666640993836\n", + "\n", + "For all three cases, the results are very close, and the functions work quickly.\n" + ] + } + ], "source": [ - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "f1 = lambda x: x\n", + "f2 = lambda x: x**2\n", + "f3 = lambda x: x**(1/2)\n", + "\n", + "n = 100000\n", + "a, b = 0, 1\n", + "\n", + "print(\"As we are calculus geniuses, we know that the first integral gives 1/2.\")\n", + "compare_integration(f1, a, b, n)\n", + "\n", + "print(\"As we are calculus geniuses, we know that the second integral gives 1/3.\")\n", + "compare_integration(f2, a, b, n)\n", + "\n", + "print(\"As we are calculus geniuses, we know that the third integral gives 2/3.\")\n", + "compare_integration(f3, a, b, n)\n", + "\n", + "print(\"For all three cases, the results are very close, and the functions work quickly.\")" ] }, {