{
 "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": 11,
   "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"
   ]
  },
  {
   "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": 1,
   "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": 13,
   "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\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))"
   ]
  },
  {
   "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": null,
   "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": [],
   "source": [
    "def test_trapz():\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "    \n",
    "def test_simps():\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\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": null,
   "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": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "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
}