{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "4ec40081b048ce2f34f3f4fedbb0be10",
     "grade": false,
     "grade_id": "cell-98f724ece1aacb67",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "# CDS: Numerical Methods Assignments\n",
    "\n",
    "- See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.\n",
    "\n",
    "- Solutions must be submitted via the Jupyter Hub.\n",
    "\n",
    "- Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\".\n",
    "\n",
    "## Submission\n",
    "\n",
    "1. Name all team members in the the cell below\n",
    "2. make sure everything runs as expected\n",
    "3. **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart)\n",
    "4. **run all cells** (in the menubar, select Cell$\\rightarrow$Run All)\n",
    "5. Check all outputs (Out[\\*]) for errors and **resolve them if necessary**\n",
    "6. submit your solutions  **in time (before the deadline)**"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "team_members = \"Koen Vendrig, Kees van Kempen\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "177d489a4104e3c95a4de1a4c7768c01",
     "grade": false,
     "grade_id": "cell-1e89a94d71771bb6",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "## Composite Numerical Integration: Trapezoid and Simpson Rules\n",
    "\n",
    "In the following we will implement the composite trapezoid and Simpson rules to calculate definite integrals. These rules are defined by\n",
    "\n",
    "\\begin{align}\n",
    "\t\\int_a^b \\, f(x)\\, dx &\\approx \\frac{h}{2} \\left[ f(a) + 2 \\sum_{j=1}^{n-1} f(x_j) + f(b) \\right] \n",
    "                          &\\text{trapezoid} \\\\\n",
    "                          &\\approx \\frac{h}{3} \\left[ f(a) + 2 \\sum_{j=1}^{n/2-1} f(x_{2j}) + 4 \\sum_{j=1}^{n/2} f(x_{2j-1}) + f(b) \\right]\t \n",
    "                          &\\text{Simpson}\n",
    "\\end{align}\n",
    "    \n",
    "with $a = x_0 < x_1 < \\dots < x_{n-1} < x_n = b$ and $x_k = a + kh$. Here $k = 0, \\dots, n$ and $h = (b-a) / n$ is the step size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "a60a63e0450a3157dd421b394288f18a",
     "grade": true,
     "grade_id": "cell-44d29c12deac7ed7",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate\n",
    "\n",
    "# And for printing the lambdas:\n",
    "import inspect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "e11bb7c15d840e7a9397f209769ebb66",
     "grade": false,
     "grade_id": "cell-ce9a56067e726f36",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 1\n",
    "\n",
    "Implement both integration schemes as Python functions $\\text{trapz(yk, dx)}$ and $\\text{simps(yk, dx)}$. The argument $\\text{yk}$ is an array of length $n+1$ representing $y_k = f(x_k)$ and $\\text{dx}$ is the step size $h$. Compare your results with Scipy's functions $\\text{scipy.integrate.trapz(yk, xk)}$ and $\\text{scipy.integrate.simps(yk, xk)}$ for a $f(x_k)$ of your choice.\n",
    "\n",
    "Try both even and odd $n$. What do you see? Why?\n",
    "\n",
    "Hint: go to the Scipy documentation. How are even and odd $n$ handled there?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "2bc6bd3985c2b7ab4ab051ebe94496f9",
     "grade": true,
     "grade_id": "cell-59f0de06f77dce3e",
     "locked": false,
     "points": 6,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def trapz(yk, dx):\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",
    "    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": 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": {
     "cell_type": "code",
     "checksum": "9599217f233689affb19148157e62b41",
     "grade": true,
     "grade_id": "cell-ff04b1d785ea895f",
     "locked": false,
     "points": 1,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# We need a function to integrate, so here we go.\n",
    "f = lambda x: x**2\n",
    "\n",
    "n    = 100001\n",
    "a, b = 0, 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "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",
      "\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": [
    "compare_integration(f, a, b, n)\n",
    "compare_integration(f, a, b, n + 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "f3a2f1f2b9ba3ffeb8646c346797d95a",
     "grade": false,
     "grade_id": "cell-1a7e33464e3be83b",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 2\n",
    "\n",
    "Implement at least one test function for each of your integration functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "c3b7ee0d3ced97054590e89bba97e031",
     "grade": true,
     "grade_id": "cell-d8f2e0aa55860e08",
     "locked": false,
     "points": 6,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "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",
    "    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",
    "    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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "ead1d68798b82e5c9c5dba354a255abb",
     "grade": false,
     "grade_id": "cell-71d20f6b94c6ed05",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 3\n",
    "\n",
    "Study the accuracy of these integration routines by calculating the following integrals for a variety of step sizes $h$:\n",
    "\n",
    "- $\\int_0^1 \\, x\\, dx$\n",
    "- $\\int_0^1 \\, x^2\\, dx$\n",
    "- $\\int_0^1 \\, x^\\frac{1}{2}\\, dx$\n",
    "\n",
    "The integration error is defined as the difference (not the absolute difference) between your numerical results and the exact results. Plot the integration error as a function of $h$ for both integration routines and all listed functions. Comment on the comparison between both integration routines. Does the sign of the error match your expectations? Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "90eaf3d9beb2347589518aba4e8ad3c4",
     "grade": true,
     "grade_id": "cell-b0bb51b7eae7769b",
     "locked": false,
     "points": 4,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "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": [
    "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.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}