{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Approximate Integration: #\n",
    "\n",
    "<p>Implementations of the following numerical integration techniques are given below:\n",
    "<a href=\"#lh\">Left-hand Riemann sum</a>,\n",
    "<a href=\"#rh\">Right-hand Riemann sum</a>,\n",
    "<a href=\"#mr\">Midpoint Rule</a>,\n",
    "<a href=\"#tr\">Trapezoid Rule</a>, and\n",
    "<a href=\"#sr\">Simpson's Rule</a>.\n",
    "Modify and evaluate the python code as you wish.</p>\n",
    "\n",
    "<p>Each function takes as input a\n",
    "function $f$, an interval $[a,b]$, and an integer $n$. Recall $\\Delta x = \\frac{b-a}{n}$ and $x_i = a+i\\Delta x$ for each $0\\le i \\le n$.</p>\n",
    "\n",
    "<br>\n",
    " <aside class=\"callout panel panel-info\">\n",
    " <div class=\"panel-heading\">\n",
    " <h3 id=\"sage-syntax\">Syntax for entering mathematical expressions in Python (`math`)</h3>\n",
    " </div>\n",
    " <div class=\"panel-body\">\n",
    " <p><ul>\n",
    "\t <li><p> Use `*` for multiplication, for example 10x would be input as `10*x`. </p></li>\n",
    "\t <li><p> Powers of x are input using `**`.  For example the square function $x^2$ is input as `x**2`.</p></li>\n",
    "\t <li><p> The exponential function $e^x$ is input as `exp(x)`. </p></li>\n",
    "\t <li><p> The natural logarithm $\\ln(x)$ (i.e. base e) is input as `log(x)`.</p></li>\n",
    " </ul>\n",
    " </div>\n",
    " </aside>\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a name=\"lh\"> Left-Hand Riemann Sum: </a>\n",
    "\n",
    "The function `lefthand_rs` outputs the left-hand Riemann sum approximation of $\\int_a^b f(x) dx$ using n partitions of the interval:\n",
    "$$\\int_a^bf(x)dx \\approx \\sum_{i=1}^n f(x_{i-1})\\Delta x = \\Delta x(f(x_0)+f(x_1)+\\cdots +f(x_{n-1})).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.388510504059924\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "def lefthand_rs(fcn,a,b,n):\n",
    "    # output: left-hand riemann sum approx of int_a^n fcn(x)dx using n steps\n",
    "    Deltax = (b-a)*1.0/n\n",
    "    return Deltax*sum([fcn(a+Deltax*i) for i in range(n)])\n",
    "\n",
    "#Example\n",
    "n=6\n",
    "a=0\n",
    "b=1\n",
    "def f(x):\n",
    "    return(math.sin(x))\n",
    "print(lefthand_rs(f,a,b,n));"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a name=\"rh\"> Right-Hand Riemann Sum: </a>\n",
    "\n",
    "The function `righthand_rs` outputs the right-hand Riemann sum approximation of $\\int_a^b f(x) dx$\n",
    "using n partitions of the interval:\n",
    "\n",
    "$$\\int_a^bf(x)dx \\approx \\sum_{i=1}^n f(x_{i})\\Delta x = \\Delta x(f(x_1)+f(x_1)+\\cdots +f(x_{n})).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4806386944084447\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "def righthand_rs(fcn,a,b,n):\n",
    "# output: right-hand riemann sum approx of int_a^b fcn(x)dx using n steps\n",
    "    Deltax = (b-a)*1.0/n\n",
    "    return Deltax*sum([fcn(a+Deltax*(i+1)) for i in range(n)])\n",
    "\n",
    "#Example\n",
    "n=20\n",
    "a=0\n",
    "b=1\n",
    "def f(x):\n",
    "    return(math.sin(x))\n",
    "print(righthand_rs(f,a,b,n));"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a name=\"mr\"> Midpoint Rule: </a>\n",
    "\n",
    "The function `midpoint_rule` outputs the midpoint rule approximation of $\\int_a^b f(x) dx$ using n partitions of the interval:\n",
    "\n",
    "$$\\int_a^bf(x)dx \\approx \\sum_{i=1}^n f(\\overline{x}_{i-1})\\Delta x,$$\n",
    "where $\\overline{x}_i$ is the midpoint of inteval $[x_{i-1},x_i]$, which is $\\overline{x}_i=\\frac{x_{i-1}+x_i}{2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.45988929071851814\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "def midpoint_rule(fcn,a,b,n):\n",
    "# output: midpoint rule approx of int_a^b fcn(x)dx using n steps\n",
    "    Deltax = (b-a)*1.0/n\n",
    "    xs=[a+Deltax*i for i in range(n+1)]\n",
    "    ysmid=[fcn((xs[i]+xs[i+1])/2) for i in range(n)]\n",
    "    return Deltax*sum(ysmid)\n",
    "\n",
    "#Example\n",
    "n=10;\n",
    "a=0\n",
    "b=1\n",
    "def f(x):\n",
    "    return(math.sin(x))\n",
    "print(midpoint_rule(f,a,b,n));"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a name=\"tr\"> Trapezoid Rule: </a>\n",
    "The function `trapezoid_rule` outputs the trapezoid rule approximation of $\\int_a^b f(x) dx$ using n partitions of the interval:\n",
    "\n",
    "$$\\int_a^bf(x)dx \\approx \\frac{\\Delta x}{2}(f(x_0)+2f(x_1)+2f(x_2)+ \\cdots + 2f(x_{n-1})+f(x_n)).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4593145488579764\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "def trapezoid_rule(fcn,a,b,n):\n",
    "# output: trapezoid rule approx of int_a^b fcn(x)dx using n steps\n",
    "    Deltax = (b-a)*1.0/n\n",
    "    coeffs = [2]*(n-1)\n",
    "    coeffs = [1]+coeffs+[1]\n",
    "    valsf = [fcn(a+Deltax*i) for i in range(n+1)]\n",
    "    return (Deltax/2)*sum([coeffs[i]*valsf[i] for i in range(n+1)])\n",
    "\n",
    "#Example\n",
    "n=10;\n",
    "a=0\n",
    "b=1\n",
    "def f(x):\n",
    "    return(math.sin(x))\n",
    "print(trapezoid_rule(f,a,b,n));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a name=\"sr\"> Simpsons Rule:</a> ##\n",
    "\n",
    "The function `simpsons_rule` outputs the Simpson's Rule approximation of $\\int_a^b f(x) dx$ using n partitions\n",
    "of the interval:\n",
    "\t\t\t\t\n",
    "$$\\int_a^bf(x)dx \\approx \\frac{\\Delta x}{3}(f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+ \\cdots + 2f(x_{n-2})+4f(x_{n-1})+f(x_n)).$$\n",
    "\t\t\t\t\n",
    "**Note:** n must be even."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4596976941318664522384551673894748091698\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "def simpsons_rule(fcn,a,b,n):\n",
    "# output: simpsons rule approx of int_a^b fcn(x)dx using n steps (n must be an even integer)\n",
    "    Deltax = (b-a)*1.0/n\n",
    "    n2=int(n/2)\n",
    "    coeffs = [4,2]*n2\n",
    "    coeffs = [1] +coeffs[:n-1]+[1]\n",
    "    valsf = [fcn(a+Deltax*i) for i in range(n+1)]\n",
    "    return (Deltax/3)*sum([coeffs[i]*valsf[i] for i in range(n+1)])\n",
    "\n",
    "#Example\n",
    "n=800;\n",
    "a=0\n",
    "b=1\n",
    "def f(x):\n",
    "    return(math.sin(x))\n",
    "print('{:.40f}'.format(simpsons_rule(f,a,b,n)))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Piecewise Defined Functions: ###\n",
    "\n",
    "Sometimes you may find the need to define $f$ as a piecewise defined function.  For example, suppose you are trying to approximate the integral $\\int_{0}^{10} f(x)~dx$ where\n",
    "$$f(x)=\\begin{cases}\n",
    "\\displaystyle x,&x\\le 1\\\\\n",
    "x^2,&x > 1\\\\\n",
    "\\end{cases}\n",
    "$$ (not that you would not need to approximate in this case since you can easily find antiderivates, but we'll just use it as a simple example).\n",
    "\n",
    "\n",
    "We can replace the line `f(x) = ...` in the code above with a python function defining $f$:\n",
    "\n",
    "```python\n",
    "def f(x):\n",
    "    if x<=1:\n",
    "        return x;\n",
    "    elif x>1:   # elif means \"else if\"\n",
    "        return x**2;\n",
    "```\n",
    "\n",
    "Try to cut-and-paste this code to replace the function $f$ in Simpson's Rule above.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.8333333333333330372738600999582558870316\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    if x<=1:\n",
    "        return x;\n",
    "    elif x>1:   # elif means \"else if\"\n",
    "        return x**2;\n",
    "\n",
    "n=800;\n",
    "a=0\n",
    "b=2\n",
    "print('{:.40f}'.format(simpsons_rule(f,a,b,n)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.11.1"
  },
  "vscode": {
   "interpreter": {
    "hash": "1a1af0ee75eeea9e2e1ee996c87e7a2b11a0bebd85af04bb136d915cefc0abce"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}