ch7Integrals.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "664a693f-3080-4696-8ca5-13e80f798568",
   "metadata": {},
   "source": [
    "(chapter:integrals)=\n",
    "\n",
    "# Numerical integration in R\n",
    "\n",
    "This chapter is a quick guide to numerical integration in R. Following the discussion you should understand how numerical integration works and be able to use R to compute numerically definite integrals of univariate functions.\n",
    "\n",
    ":::{tip}\n",
    "The discussion builds upon the material from chapters {ref}`chapter:plotting` and {ref}`chapter:functions` and assumes that you are able to work with user defined functions in R.\n",
    ":::\n",
    "\n",
    "The best way to work through the chapter is to read the discussion and then reproduce all code in your own R/RStudio installation. You can do this by copying code from a notebook cell to your R script and run the code line by line.\n",
    "\n",
    "The chapter includes some exercises that allow you to test your understanding. It is recommended that you try these exercises on your own before choosing the option to reveal the solutions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83a3764f-eaba-42bd-9bcc-bfee0ffedcc3",
   "metadata": {},
   "source": [
    "## Numerical integration - primer\n",
    "\n",
    "Consider a function of one variable, $f(x)$, and its definite integral from $a$ to $b$ denoted\n",
    "\n",
    ":::{math}\n",
    ":label: def-int\n",
    "\\int_{a}^{b} f(x) dx\n",
    ":::\n",
    "\n",
    "which corresponds geometrically to the signed area between the graph of of the function and the horizontal axis.\n",
    "\n",
    "### What is numerical integration?\n",
    "\n",
    "We have seen that the analytical approach to evaluating {eq}`def-int` is to find an anti-derivative of $f$, say $F$, and then evaluate the integral as\n",
    "\n",
    ":::{math}\n",
    "\\int_{a}^{b} f(x) dx = F(b) - F(a)\n",
    ":::\n",
    "\n",
    "By using the analytical approach we can find the exact value of  {eq}`def-int`. \n",
    "\n",
    "On the other hand numerical integration is the use of computational techniques for **approximately** evaluating {eq}`def-int` without needing to find anti-derivatives. Numerical integration is especially useful when we have to integrate functions who do not have closed-form anti-derivatives (or if we don't know the anti-derivative), but also in many cases where we are satisfied by using an approximate solution.\n",
    "\n",
    "### How does numerical integration work?\n",
    "\n",
    "We have seen in the lectures, that formally the definite integral \n",
    "\n",
    ":::{math}\n",
    "\\int_{a}^{b} f(x) dx\n",
    ":::\n",
    "\n",
    "can be defined as the limit of the sum of the areas of rectangles fitted *\"just under\"* (or *\"just over\"*) the graph of the function between $a$ and $b$ when the number of these rectangles converges to infinity (or equivalently when the *width* of a rectangle converges to zero).  \n",
    "\n",
    "While numerical integration can be done in many different ways, in general it involves evaluating {eq}`def-int` by calculating the are of rectangles fitted *\"just under\"* (or *\"just over\"*) the graph of the function between $a$ and $b$ when the number of these rectangles is specified as a \"sufficiently\" large (but not infinite - computers can't deal with infinity). The idea of numerical integration therefore draws on the formal definition of integration above, but since the number of used rectangles cannot be infinite, it leads to an approximate, rather than exact solution. However, by varying the number of rectangles fitted we can control the absolute value of the numerical error. \n",
    "\n",
    "(section:integrals:byhand)=\n",
    "\n",
    "## Numerical integration in R from first principles\n",
    "\n",
    "While R has in-built functions for calculating definite integrals numerically to high precision (discussed in section {ref}`section:integrals:integrate`), we will start by conducting numerical integration \"by hand\" in the sense of only by using basic R functions and numerical methods. The purpose of this section is to illustrate the idea of numerical integration and build intuition about how the result of numerical integration functions should be understood and interpreted. \n",
    "\n",
    "We will do this for a specific function, and specific limits, but the same procedure can be applied to most basic functions with known functional form and any integration limits.\n",
    "\n",
    ":::{note}\n",
    "The following discussion includes some graphs which are produced in R. The code for producing the graphs is hidden as it is not central to the discussion, however, you can unhide it to see how to produce these graphs yourself. You should be familiar with most of the plotting functions used from the discussion in chapter {ref}`chapter:plotting` but the codes below also use some functionality from `graphics` not discussed before (in particular, the functions `polygon()` and `rect()`).  \n",
    ":::\n",
    "\n",
    "Consider the function \n",
    "\n",
    "$$\n",
    "f(x) = x^2 + 1\n",
    "$$\n",
    "\n",
    "and suppose that we want to obtain the value of \n",
    "\n",
    "$$\n",
    "\\int_0^5 f(x) dx\n",
    "$$\n",
    "\n",
    "For context, note that this integral is easy to find analytically as follows. \n",
    "One anti-derivative of $f$ is $F(x)=\\frac{1}{3}x^3+x$. Then \n",
    "\n",
    "$$\n",
    "\\int_0^5 f(x) dx = F(5)-F(0) = \\left.\\left[\\frac{1}{3}x^3+x\\right]\\right|_0^5=140/3\\approx46.666\n",
    "$$\n",
    "\n",
    "We can represent the value of the integral graphically as the value of the shaded area in the graph below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "49e53596-eef0-4dc4-8f4c-55e7084c99f2",
   "metadata": {
    "tags": [
     "remove-cell"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=5, repr.plot.height=4, repr.plot.res=240)\n",
    "options(jupyter.plot_mimetypes = \"image/png\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "127abdfd-351d-4e6d-8e4e-873a53190ecb",
   "metadata": {
    "tags": [
     "hide-input"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABLAAAAPACAMAAAASenCyAAAANlBMVEUAAABNTU1oaGh8fHyMjIyampqnp6et2OaysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD/AAD///+9zLTcAAAACXBIWXMAACToAAAk6AGCYwUcAAAgAElEQVR4nO3d63ra2LKGUTtO+rBywvd/s6vBh0BsHIxnTan0jfFnr+5nPy0Kl96AEM7NPUATN0s/AIBLCRbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQhWoRvPLqOFL1X29MXCd4sK4UuVPX2x8N2iQvhSZU9fLHy3qBC+VNnTFwvfLSqEL1X29MXCd4sK4UuVPX2x8N2iQvhSZU9fLHy3qBC+VNnTFwvfLSqEL1X29MXCd4sK4UuVPX2x8N2iQvhSZU9fLHy3qBC+VNnTFwvfLSqEL1X29MXCd4sK4UuVPX2x8N2iQvhSZU9fLHy3qBC+VNnTFwvfLSqEL1X29MXCd4sKCy3Vbrdb5Li/c0oVEiyGW2apdruVFMspVUiwGG6RpdoJVgLBYjjBoopgMdwSS7WeXglWJcFiuAWWakW9EqxKgsVwgkUVwWK4+Uu1pl4JViXBYrjpS7WqXglWJcFiOMGiimAx3OylWlevBKuSYDHc5KVaWa8Eq5JgMZxgUUWwGG7uUq2tV4JVSbAYbupSra5XglVJsBhOsKgiWAw3c6nW1yvBqiRYDDdxqVbYK8GqJFgMN2+pdoIVRrAYbnawZh3tQk6pQoLFcNOWapW9EqxKgsVws5Zqnb0SrEqCxXCCRRXBYrhJS7XSXglWJcFiuDlLtdZeCVYlwWK4KUu1yjsaDpxShQSL4WYs1Xp7JViVBIvh5gWr/jhXcEqNcHPO0g+MzZmwVCvulWANIVjMUr9Ua+6VYFUSLIYrX6oVX8C6F6xSgsVw1Uu17l4JViXBYrg5wao9xgc4pQoJFsMVL9XKeyVYlQSL4WqXau29EqxKgsVwpUu18gtY94JVSrAYrnKp1t8rwaokWAxXuFQNeiVYlQSL4cqDVfafH8IpVUiwGK5uqTr0SrAqCRbDlS1Vi14JViXBYriqpepwAetesEoJFsMVLVWTXglWJcFiuNJglfynh3JKFRIshqtZqi69EqxKgsVwJUvVpleCVUmwGK5iqbpcwLoXrFKCxXAFS9WoV4JVSbAYrixYw/+zJZxShQSL4cYvVadeCVYlwWK44UvVqleCVUmwGG70UnW6gHUvWKUEi+EGL1WzXglWJcFiuJJgDf1PlnJKFRIshhu7VN16JViVBIvhhi5Vu14JViXBYriRS9XtAta9YJUSLIYbuFQNeyVYlQSL4cYtVcdeCVYlwWK40cEa9V+bxClVSLAYbthSteyVYFUSLIYbtVQ9eyVYlQSL4QYtVcsLWPeCVUqwGG7MUnXtlWBVEiyGGxmsEf+lyZxShQSL4YYsVdteCVYlwWK4EUvVt1eCVUmwGG7AUrW9gHUvWKUEi+GGBWvAY1mAU6qQYDHcx5eqc68Eq5JgMdyHl6p1rwSrkmAx3EeXqvMFrHvBKiVYDDcmWGMeywKcUoUEi+E+uFTNeyVYlQSL4T62VN17JViVBIvhBIsqgsVwH1qq9r0SrEqCxXAfWar+vRKsSoLFcB9YquZ3NBw4pQoJFsNdv1Rb6JVgVRIshrt6qTbRK8GqJFgMd+1SbaNXglVJsBjuY8Ea+1gW4JQqJFgMd+VSbaRXglVJsBjuuqXaSq8Eq5JgMdxVS7WRC1j3glVKsBjumqXaTq8Eq5JgMdwVS7WhXglWJcFiuPcv1ZZ6JViVBIvh3r1Um+qVYFUSLIa7Mlglj2UBTqlCgsVw712qbfVKsCoJFsO9c6k21ivBqiRYDPe+pdrWBax7wSolWAz3rqXaXK8Ea4ibc5Z+YGzOFcEqeywLcEqNIFjM8p6l2l6vBKuSYDHcO5Zqg70SrEqCxXCXL9X2LmDdC1YpwWK4i5dqk70SrEqCxXDvDFbpY1mAU6qQYDHcpUu1zV4JViXBYrgLl2qjvRKsSoLFcJct1TYvYN0LVinBYriLlmqzvRKsSoLFcO8IVvljWYBTqpBgMdwlS7XdXglWJcFiuAuWasO9EqxKgsVwf16q7V7AuhesUoLFcH9cqk33SrAqCRbD/Wmptt0rwaokWAz3h6XaeK8Eq5JgMdxFwZr0WBbglCokWAz39lJtvVeCVUmwGO7Npdp8rwSrkmAx3FtLtfULWPeCVUqwGO6NpQrolWBVEiyGO79UCb0SrEqCxXBnlyqiV4JVSbAY7txSZfRKsCoJFsOdWaqQXglWJcFiuNeXKqVXglVJsBjurWDNfiwLcEoVEiyGe3WpYnolWJUEi+FeW6qcXglWJcFiuFeWKuYC1r1glRIshnu5VEm9EqxKgsVwL5YqqleCVUmwGO73pcrqlWBVEiyGez1YyzyWBTilCgkWw/22VGG9EqxKgsVwp0uV1ivBqiRYDHeyVGEXsO4Fq5RgMdzxUuX1SrAqCRbDvQzWco9lAU6pQoLFcEdLFdgrwaokWAz3a6kSeyVYlQSL4Z6XKvAC1r1glRIshntaqsxeCVYlwWK4x6UK7ZVgVRIshntYqtReCVYlwWK4w1LF9kqwKgkWw+2XKrdXglVJsBjuv6UK7pVgVRIshru5Se6VYFUSLIZ7CtbSj2MhTqlCgsVwj8Fa+mEsxSlVSLAYbfcQrKUfxmKcUoUEi9EegrX0o1iOU6qQYDHY7hCspR/FgpxShQSLsXaHYC39KJbklCokWAx16JVgUUSwGOmQq/Clyp6+WPhuMdghWH9nL1X29MUEi4F2u7/3spcqe/pigsU4O8G6F6xSgsUwT70SLKoIFqM890qwqCJYDPKrV4JFFcFijP3ng0+9yl6q7OmLhe8WwxzuZxCse8EqFb5bjHLUK8Hi427OWfqBsQXHvRIsPk6wqHN0AUuwsqcvFr5bDHHaK8GiTPhuMcJOsI5lT18sfLcY4Nf9V4K1lz19sfDd4uN2gnUqe/pi4bvFh73slTvdKSNYfMhrvRIsqggWH/Hb5XbBuhesUoLFB+xObhgVrIPs6YsJFtc70yvBoopgcbVzvRIsqggW1zrbK8GiimBxpd/vbz/qVfZSZU9fLHy3uNru1RsaBEuwSoXvFtd6o1eCRZnw3eJab/RKsCgTvltc6Y1cCdbSD2DLwneL65y73C5Y94JVKny3uMrZzwcF616wSoXvFtf4Q68EizLhu8UVzt8wKlh72dMXC98t3m+fqzd75U53yggW7/PH11eCtfQD2DLB4l3euv9KsA6ypy8mWLzLJb0SLKoIFu9xUa8EiyqCxTtc1ivBoopgcbkLeyVYVBEsLnbRBfe/3Ye19APYsvDd4h0uuJ9BsO4Fq1T4bnG5C24YFay97OmLhe8Wl7rkflHBOsievlj4bnGh9/RKsCgTvltc5vxfOCFYL2RPXyx8t7jI+3olWJQJ3y0usROs98ievlj4bnGBS2+/EqwH2dMXC98t/uzdvXKnO2UEi7e99RcQCtZrsqcvJli86ZpeCRZVBIu3vPNyu2DdC1YpweItl38fR7CeZE9fTLB4w3W9EiyqCBbnvef7OIL1JHv6YoLFWVddv/rbfVhLP4AtC98t3nBtrwSLMuG7xXnX3M8gWPeCVSp8tzjr+l4JFmXCd4szPpArwVr6AWxZ+G7xuqvubxesB9nTFwvfLV513f3tgvUge/pi4bvFa3bv+oXIgvWb7OmLhe8Wr/horwSLMuG7xUu7K7+Pc1yspWdYVPb0xQSLUwN6JVhUESxOfPB6u2AJVinB4tjH7mcQrL3s6YsJFkfG9EqwqCJYHBnTK8GiimDxy6BeCRZVBItnAy63P/Qqe6mypy8Wvlsc+eDtooL1KHv6YuG7xS8j7mcQrHvBKhW+Wzwb1yvBokz4bvHko98fFKxn2dOPcnPO0g+MNRjZK8Hi4wSL84Z8IUewHmRPXyx8tzgYdPuVYB1kT18sfLfYG9wrwaJM+G5xP+7+9qNiLT3SorKnLyZY8cb3SrBOfP32c5HHsUmCla6gV4J14u7mRrFGEaxwI29nEKyD36d/+ijeR/IDeA6jjfh9yIL1G8Eq5DlMNvR2UcF69Eqwfjz+3wUezcZ4DoNV9UqwTny6ubnbF0uwBvAc5qrrVfZS/T79X2e/ZSJh7+YZizX26ziC9ez36b8L1jiesVR1vRKs3/wrWMN4xkJV3H4lWAcvp//2WbAG8YxlquyVYJ359+HPyxCew0ilvRKsM/8+/HkZwnOYqLZXgnXm34c/L0N4DgNV3c4gWAfZ0xcL361INV/HEawn2dMXC9+tRNWvr/52pztlBCtN3e1XgvUge/pighVmRq8EiyqClaXs+4OC9Sx7+mKClWRX+H0cwXqSPX0xwQpSfPuVYD3Inr6YYOWY1ivBoopgxZjYq+ylyp6+WPhuBZnXK8GiTPhu5Zjy6aBg7WVPXyx8t2KUfx1HsJ5lT18sfLdCzLn7SrAeZE9fLHy3Mky6+0qwHmRPXyx8tyLM7pVgUSZ8txLMur1dsB5lT18sfLcC7CbezyBYe9nTFwvfre1boFfudKeMYG3bErkSLMoI1qYt0yvBoopgbdj0q+2CtZc9fTHB2q7d5PtFBetB9vTFBGuzdlO/jiNYz7KnLyZYW7VgrwSLKoK1UYu9HfzbfVhLP4AtC9+tzVrsertgCVal8N3aqoVuZxCsvezpi4Xv1kYt2yvBokz4bm3Twr0SLMqE79YmLd0rwaJM+G5t0ZKX2wXrXrBKhe/WBi15P4Ng7WVPXyx8t7ZnBb0SLMqE79bWLHl7+3Gxln4eFpU9fTHB2pKV9EqwqCJYG7Lo7e3Hspcqe/pigrUdi9/O8Cx7qbKnLyZYW7HIL28/I3upsqcvJlgbsaZeCRZVBGsb1pQrwaKMYG3CWq62Pwhfquzpi4Xv1kas4W7RI+FLlT19sfDd2oTl/rKJM8KXKnv6YuG7tQVruVv0l/Clyp6+WPhubcD6eiVYlAnfrf5Wc3f7kfClyp6+WPhutbeu2xkehS9V9vTFwneru1X2SrD4uJtzln5gXG+dvRIsPk6wNmdV38Y5kb1U2dMXE6y2Vni1/Un2UmVPX0ywulrb3aLHspcqe/pigtXUmnslWFQRrJ5W3SvBoopgdbTCu9tPZC9V9vTFBKuh1X3b+XfZS5U9fTHB6me1dzM8CV+q7OmLhe9WR6vvlWBRJny3+lnv3aK/hC9V9vTFwnernQ69EizKhO9WNyu+u/1I+FJlT18sfLd6Wf2ng4/Clyp7+mLhu9VKl14JFmXCd6uTNr0SLMqE71Yfu7Xf3X4kfKmypy8Wvltt7Pq8vvrbne6UEawWWtzN8Ev2UmVPX0ywGmiWK8GijGCtX7teCRZVBGv1etwseiJ7qbKnLyZYa9fpYvuT7KXKnr6YYK1bqw8Hn2UvVfb0xQRr1Zr2KnupsqcvFr5bK9foZtFj4UuVPX2x8N1atV43ix4JX6rs6YuF79aa9bub4Un4UmVPXyx8t1asb68EizLhu7Vau8a9EizKhO/WWrXulWBRJny3Vqrpp4NPwpcqe/pi4bu1Tg+9EqymsqcvFr5bq9T0ZoYj2UuVPX0xwVqbrjdfHctequzpiwnWyjzcLbp0cT4oe6mypy8mWOvS+cPBX7KXKnv6YoK1KtvolWBRRbBWpPXNV8eylyp7+mKCtR6b6ZVgUUWw1qL5zaLHwpcqe/pi4bu1Ht1vFj0WvlTZ0xcL36212MLNV7+EL1X29MXCd2slttUrwaJM+G6twyZuFj0SvlTZ0xcL36012M6Hg0/Clyp7+mLhu7UC2+uVYFEmfLeWt8FeCRZlwndradu62P4kfKmypy8WvlsL29DNoieylyp7+mKCtaBtvr76W7AoI1iL2djNV8eylyp7+mKCtZQN90qwqCJYC9nazaInspcqe/pigrWIDd7LcCx7qbKnLyZYS9h4rwSLKoI13xbvFT0RvlTZ0xcL360lbPni1YPwpcqevlj4bs23pV/Ud074UmVPXyx8t6ZL6JVgUSZ8tybb8r1XR8KXKnv6YuG7Ndem7706Er5U2dMXC9+tmXab/3TwSfhSZU9fLHy3Jkqp1d+CtfQD2LLw3Zom5tXVXvhSZU9fLHy3Zgm5ePUke6mypy8mWBPstvqL+s7JXqrs6Ue5OWfpB7Z9EfdenchequzpRxGshcS9vPpbsCgjWMVS7r06kb1U2dMXE6xaUZ8OPstequzpiwlWpcha/S1YlBGsQrG9yl6q7OmLhe9WpcSLVw/Clyp7+mLhu1Un716GX8KXKnv6YuG7VSa5V4JFmfDdKhLye6/OCV+q7OmLhe9WiV3kvVdHwpcqe/pi4btVIPPOqxPhS5U9fbHw3RpPrwRr6QewZeG7NVr4e8FH4UuVPX2x8N0aLPxi+7PspcqevphgjZP4exlel71U2dMXE6xh9OpZ9lJlT19MsAbZeTv4S/ZSZU9fTLDG8OHgseylyp6+mGANoVcnspcqe/pigjWAWv0me6mypy8mWB/n3qvfhC9V9vTFwndrgORfy3BG+FJlT18sfLc+LPrXyJwTvlTZ0xcL362PcivDa8KXKnv6YuG79TH7i+169VL4UmVPXyx8tz7CrQznhC9V9vTFwnfrenJ1XvhSZU9fLHy3ruXXyLwlfKmypy8WvlvX8dHg28KXKnv6YuG7dY2dXv1J9lJlT19MsN7Jb2W4QPZSZU9fTLDe5fAX4riV4U+ylyp7+mKC9R47nw1eJHupsqcvJliXU6tLZS9V9vTFBOtS6X876ntkL1X29MUE6zI+GnyP7KXKnr6YYF1Er94jfKmypy8WvluXcSvD+4QvVfb0xcJ36yIutr9T+FJlT18sfLcuIFfvFr5U2dMXC9+tP/LJ4BXClyp7+mLhu/UHrl1dJXypsqcvFr5bb/LJ4JXClyp7+mLhu/UGdzJcLXypsqcvFr5bZ+3c2H698KXKnr5Y+G6d4UvOH5O9VNnTFxOsV6jVB2UvVfb0xQTrBZ8Mflj2UmVPX0ywfvN4rX3pU7637KXKnr6YYJ2QqyGylyp7+mKCdWSnV2NkL1X29MUE65lPBofJXqrs6YsJ1iN3MowTvlTZ0xcL360nPhkcKXypsqcvFr5bD1y6Git8qbKnLxa+W3tyNVr4UmVPXyx8t3wyWCF8qbKnLxa+W660VwhfquzpiyXvlg8GiyQv1b1glcrdLb9NtEzuUh1kT18sdbf8er5CqUv1KHv6Ypm7JVe1IpfqWfb0xRKD5YPBaoFLdSR7+mJ5wXKlvV7cUp3Inr5YWrB8MjhD2FL9Jnv6UW7OWfqBTeSDwUmSluql7OlHiQ/W04V2waoXs1Svyp6+WEqwdj4YnChkqc7Inr5YRrAe/5ZBtZojY6nOyp6+WMBu7Vxonyxgqd6SPX2xze+WWs23+aV6W/b0xTa+W+4QXcLGl+pPsqcvtuXd8rngQra8VBfInr7YZndr5ws4i9nsUl0me/piG92tnc8FF7TRpbpU9vTFNrlbLrQva5NLdbns6YttcLfkamkbXKr3yJ6+2OZ2y4WrFdjaUr1P9vTFthUsnwuuw6aW6t2ypy+2oWD5XHA1trNU18ievthWgvX4/Ru5WoWNLNWVsqcvto1guc6+LptYqqtlT19sA8Fyw9Xq9F+qj8ievlj3YPk9V2vUfKk+KHv6Yr2DJVer1HupPix7+mKNd+v5MrtcrUzjpRohe/piXXfLb+Vbsa5LNUj29MVa7pbL7OvWcqnGyZ6+WLfder47VLDWq9tSDZY9fbFWu7Xb7Vxmb6DVUo2XPX2xPrv1WKvd/nxY+ozkTX2WqkT29MWa7NbjK6ulz0Qu0mSpqmRPX6zDbrls1UyHpSqUPX2x1e/WzmX2fta+VLWypy+27mA9XbYSq15WvVTlsqcvtuJgPV9kp5v1LtUM2dMXW2mw3L3Q2jqXapbs6YutMFhHF60Eq6f1LdVM2dMXW1mwnm62EqvW1rVUs2VPX2xFwXqKlatW7a1nqZaQPX2xdQTr6HXV0ucaH7eOpVpM9vTFlt+tk3eBgrUFyy/VorKnL7bsbu1+vQ1Uq+0QLKost1vPL6xcs9oawQqwP3kXOOwyu+VN4KYJ1vY9vzeafNz5u/U8qV5tlWBt3+7YxONO3a2nTj3cuaBVWyVY27f73aTjTtstL6pyCFaAF8Wak60pu3X6JnDp04lqgpXh9WbVdqt4t54+CPQmMErMKfuqqOlnR6smWC6sZ4s6ZV+Im/58tMZna3SwTh+tXmWKO2VPRE7/VrRGhmtksI4jtfQpw5IiT9lnwdP/KVsfL9eYYHlRxbHgU/Y+OliPLujWteH6ULBOEvXwDRutQrC4LFpXlOs9wXr9eH43KL/LPmWzpz91cbYuDdglwXq7UFLFb9yHxW+uCNfrDvdI/fH/5/FN3+7hL4r3koo3CRbnjAjWcZB+D9TRPyoUlxEs/uQDwfo9SALFxwgWF3tvsJZebrZHsLiaYDGbYFFFsBhOsKgiWAwnWFQRLIYTLKr4LJDxsk/Z7OmLCRbjZZ+y2dMXEyzGyz5ls6cvJliMl33KZk9fTLAYL/uUzZ6+mGAxXvYpmz19McFivOxTNmX6mwU9b5p/9s8D/nkhS5/CD1byMIot/KOG/pY+iR+s5GEUW+jZzjqscR223koeRrGsH7IzeMvHzTrsCyt5GMWyfsjO4C0fN+uwL6zkYRTL+iE7g7d83KzDvrCSh1Es64fsDN7ycbMO+8JKHkaxrB+yM3jLx8067AsreRjFsn7IzuAtHzfrsC+s5GEUy/ohO4O3fNysw76wkodRLOuH7Aze8nGzDvvCSh5GsawfsjN4y8fNOuwLK3kYxbJ+yM7gLR8367AvrORhFMv6ITuDt3zcrMO+sJKHUSzrh+wM3vJxsw77wkoeRrGsH7IzeMvHzTrsCyt5GMWyfsjO4C0fN+uwL6zkYRTL+iE7g7d83KzDvrCSh1Es64fsDN7ycbMO+8JKHkaxrB+yM3jLx8067AsreRjFsn7IzuAtHzfrsC+s5GEUy/ohO4O3fNysw76wkodRLOuH7Aze8nGzDvvCSh5GsawfsjN4y8fNOuwLK3kYxbJ+yM7gLR8367AvrORhFMv6ITuDt3zcrMO+sJKHAfBnggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG1HB+nQzcdyvX25vbj59+TrviM+mzvng65e7m/24/5t72O9//Tfr7SJP8n++3Ux9on/eHJt44Psf/+x/und/fZ950NclBet/E3/MPz49LdanH7OO+WTmnA++3T6fR7f/Tjzu56ej3k1/kvdu5z7R3xYK1s/np/nm87yjnhEUrB8T9+vH7dFqTT6ZZs754K+TM+mvace9+3XQ2wWK9WXyC51/lwnWyS7fTjvsGTnBenjeJx3s7ni1Pk066IOpcx787+bUP5OO+2W5J3nv2+QXOr9eT0498M/bk8Mu/RorJliPf07MOdjXw4/2239LfSjXtzlHPZg654P9EW//2V/f+P7wWuvnlMN+3x/qy3+H/XEo1/TrWNOf6NtpT+2xz88/3a+fZi/zK1KC9XSVZc7R7n69M9r/vL/MOere3DkP/nf8juzwsmPOm8K/fj2znxf4o//L7Cf65yLvyL4fR+pu7jK/JiRY++26m7ZfP4+O9HPqVs+d88Hnkz92/5329uzTrxccP2a/qPx1BXzuEee/Ifty/AfQj0WaeSIjWPuXHXf30/br2/EfRF8mvoyePOeD08PNDfTrj2GG2+nB2v9RMPMj2AfLJ+pURrAergRP26+T1Zq5Z5PnPPcYFgnW3Kvu+z+GPs0d9fPEP/mefZv2Dv9CIcH6tL9oOG2/9qv1fI/d94mv5CfP+Zqf09txv8BpdXhD+H3uE3245v7zn/8yeTfvhda/S3yc8ZaMYD3cf739YE2e8zX/LvFH8t3sz89uD69lpz7Rh2vuXx/fid7OisjjKn/ff23j7p/5n1G+kBGsB9P26/TOzen9WDRYt8exnuNw78jUlwFfHl5GTn2i9y/qjm6JmjTv58OIT3eATevkeYJVfqCoYN1N/yjrcDJ9mtrIhzeEk5/o0/vcZ/2pcPiz96iTs24KPkuwyszPHMsAAAW6SURBVA+UFKz9vVGTvyVzuIvjbuqXrm8fz9ypT/QhzPv7ZO9//Hs77ULh4ck97uTC940KVv2BgoL118Q3K08+ffnrbu5rrC9PuZj6RN8ePbU/pqXjoVL7L23cf//n8K5wxlHfekALH38mwSp2t0CvDn7cTfyz/+kN4dJ/Mkz5aOPQq6cjfV/q53tEsMoPFBOsxXp1f3j9cTvpM6ynN4RLBmufjrsZB7o5uSZ59G2opWwzWEfvuX/7t1MO/+JTwrmvo5c5jx6+db3Un79fp10Qfn5DuOjFwlmHPn3v+X3atbOzBKvAYvdhPVjkPDrcIbTcx94/Z73kON6tV3ZsllmHvVvoFDpLsAoEButwuX2J36L3ZOr7/ZxgfRas5Ux7tve3zDx/zr7/7Stzv7O6wFY93Fsw+aAnsoI17QtQ+1U+ujYoWDNNe7YX+20ND+Zv1aFXs7+Qc3N8Lu1fx865UrhUsE7GPdmwSt+O/+xd6DfcnBCs6iOdLNr0o0+xxO1XD+9Wnl+6LvL51dQn+vPxc/x52hN+8gfB/lle+F53wapw9GtGP89/rzQ7WIdvuc3/G6AOv+n08Y+CZe4QmvpE/+8oHV+nvZ48+ZPg6/w/fF8QrAqHn+zd8+90n3wmzQ7W/vdCLfE31h2+oLJ/bn8cvmc3//P2uU/0b+NOui7682mV738cXkj7Syjmmbhfz38r4RJn0uRgfV3qKvRvB57/EeXcJ/r028/TXrWf/p1Is+7OPUuwShz/nZfTP+yfHKyTv21rZrBO/j7EJW4Bm/xEH4878SrDP0vu8guCVeP786/kmP83P08+j07/3rqZwbr/+nzsRf7q58lP9NFrrKmfL3xb9lk+JVhVvu5/SePtl+3/wf+yVxOP/r/9hxqf/lrml57MDtb9z8O4d3/NvmJ42OVPS+zyC0nBApoTLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECwW9OXm5ubT0z98++8fbr4v+XBYPcFiSbf/NeqfV/43vEqwWNL/9q+qfhz+58mrLXiVYLGoz/9l6m7/P7wh5AKCxaJ+7jv17703hFxEsFjWv/ti/fSGkIsIFgu7+y9Vn70h5CKCxcK+3zzyhpA/EiyW9s9Dr7wh5M8Ei8V9uvGGkMsIFov7uu/VX0s/CjoQLBb3xTtCLiRYLO2ba+5cSrBY2u3jp4QuYvFHgsXC9m8Ib5+/oQNvESyWdXhD+O3n7eM3dOAtgsWyHr9D+PgNHXiTYLGo5+8QHr6hs/SjYe0EiyX9+g7hj/3/+t/Sj4eVEyyWdPRLZfbf0Lld9tGweoLFgk5+qcz+Gzpflnw0rJ9gAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtDG/wGxyuTMD0vAjwAAAABJRU5ErkJggg=="
     },
     "metadata": {
      "image/png": {
       "height": 480,
       "width": 600
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f <- function(x){x^2 + 1}\n",
    "\n",
    "a <- 0\n",
    "b <- 5\n",
    "\n",
    "plot(f,from = -1, to = 6 , col=\"red\", lwd=2, axes=FALSE)\n",
    "axis(1, pos=0)\n",
    "axis(2, pos=0, labels=FALSE, las=2)\n",
    "abline(v=c(0,5))\n",
    "\n",
    "xpol <- c(a,seq(a,b,by=0.01),b)\n",
    "ypol <- c(0,f(seq(a,b,by=0.01)),0)\n",
    "polygon(xpol,ypol, col=\"lightblue\", lty=\"dashed\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42829459-f0e8-4464-848f-1588440625ae",
   "metadata": {},
   "source": [
    "Now let's start computing this definite integral numerically. \n",
    "\n",
    "Let's start by approximating the area under the graph of the function with the area of two rectangles of equal width just below the graph of the function, as in the graph below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "707be6d2-136e-469c-ab65-abfcf1e1b48b",
   "metadata": {
    "tags": [
     "hide-input"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABLAAAAPACAMAAAASenCyAAAANlBMVEUAAABNTU1oaGh8fHyMjIyampqnp6et2OaysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD/AAD///+9zLTcAAAACXBIWXMAACToAAAk6AGCYwUcAAAgAElEQVR4nO3d7WIT15atYRMgvfOJff8325FjiI0la5W0ZlXNGs/z52T39gme1UNvW4Dh4QmgiYetPwGAUYIFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCVejB02W28FFlX18sfFtUCB9V9vXFwrdFhfBRZV9fLHxbVAgfVfb1xcK3RYXwUWVfXyx8W1QIH1X29cXCt0WF8FFlX18sfFtUCB9V9vXFwrdFhfBRZV9fLHxbVAgfVfb1xcK3RYXwUWVfXyx8W1QIH1X29cXCt0WF8FFlX18sfFtUCB9V9vXFwrdFhfBRZV9fLHxbVAgfVfb1xcK3RYXwUWVfXyx8W1QIH1X29cXCt0WF8FFlX18sfFtU2GhUj4+b/LDveEkVEiymEyyqCBbTbTOqvfRKsCoJFtNtMqpHwUogWEy3UbA2+EHP8pIqJFhMt8Wo9tMrwaokWEy3waj284ZQsEoJFtNtEqzVf8iLvKQKCRbTrT+qPfVKsCoJFtOtPqo9vSEUrFKCxXQbBGvlH/BDXlKFBIvp1h7VvnolWJUEi+lWHtW+3hAKVinBYrrVg7XqD3eVl1QhwWK6dUe1t14JViXBYrpVR7W3N4SCVUqwmG7NUe2vV4JVSbCYbt1grfdjDfKSKiRYTLfiqHbYK8GqJFhMt96odviGULBKCRbTrTaqXfZKsCoJFtOtGKyVfqBFvKQKCRbTrTWqffZKsCoJFtOtNKp9viEUrFKCxXTrjGqvvRKsSoLFdGsFa40f5QZeUjM8XLL1J8bhrDKq3fZKsKYQLNayxqh2+4ZQsEoJFtOtMKod90qwKgkW09WPas+9EqxKgsV05aPada8Eq5JgMd0KwSr+Ae7iJVVIsJiuelT77pVgVRIspise1b7fEApWKcFiutpR7b1XglVJsJiudFS775VgVRIspisOVuG/fAovqUKCxXSVo9p/rwSrkmAxXeGo9v+GULBKCRbT1Y2qQ68Eq5JgMV3ZqFr0SrAqCRbTFQar6F88lZdUIcFiuqpR9eiVYFUSLKYrGlWPN4SCVUqwmK5mVF16JViVBIvpSkbVpleCVUmwmK4oWAX/0hJeUoUEi+kqRtWnV4JVSbCYrmBUfd4QClYpwWK6+aPq1CvBqiRYTDd9VK16JViVBIvpZo+qV68Eq5JgMd38YM399xXzkiokWEw3eVTNeiVYlQSL6eaOqtkbQsEqJVhMN3VU7XolWJUEi+lmjqpfrwSrkmAx3dxgzft3rcRLqpBgMd3EUTXslWBVEiymmzeqhm8IBauUYDHdtFG17JVgVRIspps1qp69EqxKgsV084I159+zMi+pQoLFdJNG1bRXglVJsJhuzqiaviEUrFKCxXRTRtW2V4JVSbCYblKwJvxLNuElVUiwmG7GqPr2SrAqCRbTTRhV3zeEglVKsJju/lF17pVgVRIsppsRrBmfx0a8pAoJFtPdParWvRKsSoLFdPeOqvUbQsEqJVhMd+eomvdKsCoJFtPdHaxJn8dGvKQKCRbT3Teq7r0SrEqCxXR3jar7G0LBKiVYTHfPqPr3SrAqCRbT3TGqA/RKsCoJFtPdPqoj9EqwKgkW0908qkP0SrAqCRbT3RGsqZ/HRrykCgkW0906qmP0SrAqCRbT3TiqY7whFKxSgsV0t43qKL0SrEqCxXQ3jeowvRKsSoLFdLeM6ji9EqxKgsV0N4zqQL0SrEqCxXQ3Bavg89iIl1QhwWK65aM6Uq8Eq5JgMd3iUR3pDaFglRIspls6qmP1SrAqCRbTLRzVwXolWFM8XLL1J8bhLA5W0eexES+pGQSLtSwb1dF6JViVBIvpFo3qaG8IBauUYDHdklEdr1eCVUmwmG7BqA7YK8GqJFhMtyhYhZ/HRrykCgkW042P6oi9EqxKgsV0w6M64htCwSolWEw3Oqpj9kqwKgkW0w2O6qC9EqxKgsV0w8Eq/jw24iVVSLCYbmxUR+2VYFUSLKYbGtVR3xAKVinBYrqRUR23V4JVSbCYbmBUB+6VYFUSLKa7Pqoj90qwKgkW010d1aF7JViVBIvpro3q2L0SrEqCxXTXg7XO57ERL6lCgsV0V0Z18F4JViXBYrqPR3XwN4SCVUqwmO7DUR2+V4JVSbCY7qNRHb9XglVJsJjug1EF9EqwKgkW010eVUKvBKuSYDHdxVFF9EqwKgkW010aVUavBKuSYDHd5WCt+3lsxEuqkGAx3YVRhfRKsCoJFtOdH1XIG0LBKiVYTHd2VDG9EqxKgsV050aV0yvBqiRYTHdmVEG9EqxKgsV070eV1CvBqiRYTHcuWFt8HhvxkiokWEz3blRRvRKsSoLFdD+PKuoNoWCVEiym+2lUYb0SrEqCxXRvR5XWK8GqJFhM93Owtvo8NuIlVUiwmO7NqOJ6JViVBIvpXo8q7g2hYJUSLKZ7NarAXglWJcFiuv9GldgrwaokWEz3Y1SRvRKsSoLFdN9HldkrwaokWEz3MqrQXglWJcFiun9HldorwaokWEz3PKrYXglWJcFiutOocnslWJUEi+n+GVVwrwSrkmAx3cNDcq8Eq5JgMd0pWFt/DhvykiokWEz3T7C2/hS25CVVSLCYLnxU2dcXC98WBR7DR5V9fbHwbTHfo2BRJnxbTPcoWFt/AkcWvi2me3xMH1X29cXCt8Vsp9/QED6q7OuLhW+LyR4FS7AqhW+Luf79He7ho8q+vlj4tpjq5TtywkeVfX2x8G0x0/fvIAwfVfb1xcK3xUQ/vuM5fFTZ1xcL3xYT/fiO5/BRZV9fLHxbzPPfn9AQPqrs62d5uGTrT4xjeBSsF9nXzyJYVHr9R/aFjyr7+mLh22KSN3/EaPiosq8vFr4t5nj7RyKHjyr7+mLh22KKn/4I9/BRZV9fLHxbzPDzXzkRPqrs64uFb4sJ3v0VOeGjyr6+WPi2uN/7v9IrfFTZ1xcL3xZ3O/NXEIaPKvv6YuHb4l7n/srU8FFlX18sfFvc6exf8Rw+quzri4Vvi/uc/yvpw0eVfX2x8G1xl/O9Sh9V9vXFwrfFPS70Kn1U2dcXC98Wd7jUq/RRZV9fLHxb3O5ir9JHlX19sfBtcbuLvUofVfb1xcK3xc0u9yp9VNnXFwvfFre6/IYwflTZ1xcL3xY3+qhX6aPKvr5Y+La4zYe9Sh9V9vXFwrfFTT7uVfqosq8vFr4tbnGlV+mjyr6+WPi2uMG1XqWPKvv6YuHbYrmrvUofVfb1xcK3xXJXe5U+quzri4Vvi8Wu9yp9VNnXFwvfFktdf0MYP6rs64uFb4uFRnqVPqrs64uFb4tlhnqVPqrs64uFb4tFxnqVPqrs64uFb4slBnuVPqrs64uFb4sFRnuVPqrs64uFb4thj8O9Sh9V9vXFwrfFqAW9Sh9V9vXFwrfFoAW5ih9V9vXFwrfFmEW9Sh9V9vXFwrfFkGW9Sh9V9vXFwrfFiIW9Sh9V9vXFwrfFgKW9Sh9V9vXFwrfFdYt7lT6q7OuLhW+Lq5b3Kn1U2dcXC98WVy3vVfqosq8vFr4trrmhV+mjyr6+WPi2uOKGN4Txo8q+vlj4tvjYTb1KH1X29cXCt8WHbutV+qiyry8Wvi0+cmOv0keVfX2x8G3xgVt7lT6q7OuLhW+Ly27uVfqosq8vFr4tLlny51/9LHxU2dcXC98WF9zTq/RRZV9fLHxbnHdPruJHlX19sfBtcdZ9vUofVfb1xcK3xTl39ip9VNnXFwvfFmfc26v0UWVfXyx8W7x3d6/SR5V9fbHwbfHO/b1KH1X29cXCt8XPJvQqfVTZ1xcL3xY/m9Cr9FFlX18sfFv8ZEav0keVfX2x8G3x1ow3hPGjyr6+WPi2eGNOr9JHlX19sfBt8dqkXqWPKvv6YuHb4pVZvUofVfb1xcK3xX+m9Sp9VNnXz/JwydafGPswr1eCxf0Ei49M7JVgUSZ8W7yY2av0UWVfXyx8W/xraq/SR5V9fbHwbfFsbq/SR5V9fbHwbXEyuVfpo/r5+j/+/LbJ53FI4dviaX6v0kf18/WfHx4Ua5bwbVHQq/RR/Xz991+K90vyE3iG6eb3Kn1UglXIM8x2198/eEn4qM4E6++X/3eDz+ZgPMNoJb1KH9XP1//y8PD5VCzBmsAzTFaSq/hR/Xz9rxe/y0TCFvPEghX1Kn1UP1//l2DN44nlqupV+qjeXf+bYE3jicUq61X6qN5f/+cXwZrEE0tV16v0UV26Xp4m8AxDFfYqfVSCVcgzzFTZq/RRCVYhzzBSaa/SR5V9fbHwbYWq7VX6qLKvLxa+rUzFvUofVfb1xcK3Fam6V+mjyr6+WPi2EpX3Kn1U2dcXC99WoPpepY8q+/pi4dvKs0Kv0keVfX2x8G2lqfnjZH4WPqrs64uFbyvMOr1KH1X29cXCt5VlnVzFjyr7+mLh24qyVq/SR5V9fbHwbSVZrVfpo8q+vlj4toKs16v0UWVfXyx8WzlW7FX6qLKvLxa+rRQr/fLgi/BRZV9fLHxbIdbtVfqosq8vFr6tDOvmKn5U2dcXC99WhLV7lT6q7OuLhW8rweq9Sh9V9vXFwrcVYP1epY8q+/pi4ds6vg16lT6q7OuLhW/r8LboVfqosq8vFr6tg1v5tzN8Fz6q7OuLhW/r2DbqVfqosq8vFr6tQ9soV/Gjyr6+WPi2jmyzXqWPKvv6YuHbOrDtepU+quzri4Vv67g27FX6qLKvLxa+rcPaslfpo8q+vlj4to5q016ljyr7+mLh2zqobXuVPqrs64uFb+uYNu5V+qiyry8Wvq1D2rpX6aPKvr5Y+LaOaPNepY8q+/pi4ds6oO17lT6q7OuLhW/rcLb69sE3wkeVfX2x8G0dzS56lT6q7OuLhW/rYHaRq/hRZV9fLHxbx7KTXqWPKvv6YuHbOpJ9vB08CR9V9vXFwrd1IPvpVfqosq8vFr6t49hPruJHlX19sfBtHcaeepU+quzri4Vv6yh21av0UWVfXyx8W8ewo5++ehY+quzri4Vv6xD21qv0UWVfXyx8W0ewt1zFjyr7+mLh2zqA/fUqfVTZ1xcL31Z/O+xV+qiyr5/l4ZKtPzHuscdeCRb3E6wj2mWvBIsy4dvqbXe/PPgifFTZ1xcL31Zre+1V+qiyry8Wvq3O9pqr+FFlX18sfFuN7bdX6aPKvr5Y+Lb62nGv0keVfX2x8G11tdufvnoWPqrs64uFb6upffcqfVTZ1xcL31ZP+85V/Kiyry8Wvq2W9t6r9FFlX18sfFsN7fzt4En4qLKvLxa+rX4a9Cp9VNnXFwvfVjsNchU/quzri4Vvq5kOX149xY8q+/pi4dvqpUmv0keVfX2x8G210iRX8aPKvr5Y+LYa6fLl1VP8qLKvLxa+rT4a9Sp9VNnXFwvfVhuNchU/quzri4Vvq4lOX149xY8q+/pi4dvqoVmv0keVfX2x8G210CxX8aPKvr5Y+LY6aNer9FFlX18sfFv71+3t4En4qLKvLxa+rd3r2Kv0UWVfXyx8W3vXMVfxo8q+vlj4tvat5ZdXT/Gjyr6+WPi2dq1rr9JHlX19sfBt7VnXXMWPKvv6YuHb2q+2X149xY8q+/pi4dvarc69Sh9V9vXFwre1V51zFT+q7OuLhW9rp3r3Kn1U2dcXC9/WPjXvVfqosq8vFr6tPWr901fPwkeVfX2x8G3tUP9epY8q+/pi4dvan/65ih9V9vXFwre1O0foVfqosq8vFr6tnTnA28GT8FFlX18sfFv7cpBepY8q+/pi4dvak6PkKn5U2dcXC9/WjhynV+mjyr6+WPi2duNAuYofVfb1xcK3tReH6lX6qLKvLxa+rZ04VK7iR5V9fbHwbe3Csb68eoofVfb1xcK3tQeH61X6qLKvLxa+rR04XK7iR5V9fbHwbW3ueF9ePcWPKvv6YuHb2tohe5U+quzri4Vva2OHzFX8qLKvLxa+rU0d88urp/hRZV9fLHxbWzpsr9JHlX19sfBtbeiwuYofVfb1xcK3tZnjfnn1FD+q7OuLhW9rK4fuVfqosq8vFr6tbRw7V/Gjyr6+WPi2NnH0XqWPKvv6YuHb2sDhcxU/quzri4Vva30BvUofVfb1xcK3tbaEXMWPKvv6YuHbWllGr9JHlX19sfBtrSokV/Gjyr6+WPi21hTTq/RRZV9fLHxb68nJVfyosq8vFr6t1ST1Kn1U2dfP8nDJ1p9YgqhcCdbWn8AhCNZ2wnolWJQJ39Ya0nIVP6rs64uFb2sFeb1KH1X29cXCt1UvL1fxo8q+vlj4tqoFfnn1FD+q7OuLhW+rWGav0keVfX2x8G2VCs1V/Kiyry8Wvq1Cj7G9Sh9V9vXFwrdVJzdX8aPKvr5Y+LaqJOcqflTZ1xcL31aN4HeDz8JHlX19sfBtVUjPVfyosq8vFr6tAvG5ih9V9vXFwrc1nVw9xY8q+/pi4duaTa9OwkeVfX2x8G3NJVf/Ch9V9vXFwrc1lV69CB9V9vXFwrc1kVz9ED6q7OuLhW9rHr36T/iosq8vFr6taeTqlfBRZV9fLHxbk/jy6o3wUWVfXyx8W3Po1Vvho8q+vlj4tmaQq5+Fjyr7+mLh27qf7xx8L3xU2dcXC9/W3eTqjPBRZV9fLHxbd5Krs8JHlX19sfBt3cW7wQvCR5V9fbHwbd1Bri4KH1X29cXCt3UzufpA+Kiyry8Wvq0bydWHwkeVfX2x8G3dRK6uCB9V9vXFwrd1A7m6KnxU2dcXC9/WYnI1IHxU2dcXC9/WUnI1InxU2dcXC9/WMnI1JnxU2dcXC9/WEo8PjNr6f1ebyr6+WPi2xsnVElv/b2tT2dcXC9/WsMd/gvV/jAkfVfb1xcK3Nej0k1eCNSx8VNnXFwvf1pjnn2wXrGHho8q+vlj4tka8/NqgYA0LH1X29cXCt3Xdj98pKljDwkeVfX2x8G1d8+o3tgvWsPBRZV9fLHxbH3vzfTiCNSx8VNnXFwvf1kd++rZBwRoWPqrs64uFb+uyd9/lLFjDwkeVfX2x8G1dcuYPZRCsYeGjyr6+WPi2Ljj3Xc6CNSx8VNnXFwvf1lnn/1AGwRoWPqrs64uFb+uMS39En2ANCx9V9vXFwrf1zuU/UVSwhoWPKvv6YuHb+slHfwCyYA0LH1X29cXCt/XGx39eu2ANCx9V9vXFwrf1yrW/XkKwhoWPKvv6YuHb+uH634YjWMPCR5V9fbHwbb0Y+cu7BGtY+Kiyry8Wvq1nY3/XoGANCx9V9vXFwrf1NP5XowrWsPBRZV9fLHxbC/4mZ8EaFj6q7OuLRW/rcclfPC9Yw6JHJVilgre1KFeCtUDwqE6yry8Wu61ltXoSrAViR/Wv7OuLhW5rca4Ea4HQUX2XfX2xyG3dkCvBWiByVP/Jvr5Y4LZuypVgLRA4qteyr5/l4ZKtP7GV3ZgrwVogblRvZV8/i2A9Lf6FwTcEa1jWqN7Jvr5Y0rbuyZVgLZA0qjOyry8Ws637avUkWAvEjOq87OuLhWzr7lwJ1gIho7ok+/piCduaUKsnwVogYVQfyL6+2PG3NSdXgrXA8Uf1oezrix19W5Nq9SRYCxx9VFdkX1/s0Nua9cXVM8EaduhRXZd9fbEDb2tqrgRrgQOPakT29cWOuq3JtXoSrAWOOqpB2dcXO+a25udKsBY45qiGZV9f7IjbKqjVk2AtcMRRLZB9fbHjbasmV4K1wPFGtUj29cWOtq2qXAnWAkcb1ULZ1xc71LYqfurqB8EadqhRLZd9fbHjbOuxNFeCtcBxRnWT7OuLHWRb1bV6EqwFDjKqW2VfX+wI21qhVk+CtcARRnWH7OuL9d/WKrV6EqwF+o/qLtnXF2u+rbVq9SRYCzQf1b2yry/WeVvrvBX8TrCGdR7VBNnXF+u7rVVr9SRYC/Qd1RTZ1xdruq21a/UkWAs0HdUs2dcX67etx8d13wp+J1jD+o1qquzri/Xa1uNWtXoSrAV6jWq67OuLNdrWdq16JljDGo2qQvb1xdpsa9NYnQjWsDajqpF9fbEe29q8Vk+CtUCPUZXJvr5Yg21t+07wB8Ea1mBUlbKvL7b3be2kVk+CtcDeR1Us+/piu97Wfmr1JFgL7HpU9bKvL7bbbW38a4LvCdaw3Y5qHdnXF9vltjb83VaXCdawXY5qPdnXF9vdtnYZqxPBGra7Ua0r+/piu9rWbmN1IljDdjWq9WVfX2wv29rym27GCNawvYxqI9nXF9vDtvYfqxPBGraHUW0o+/piW2+rRaueCdawrUe1sezri225rT6xOhGsYYJFla221StWJ4I1TLASbPPy3WJb/WJ1IljDBCvBNi/jlbfV4+fXzxKsYYKVYv2X84rbahyrE8EaJlhJ1n1Zr7St3q16JljDBCvN4xuVP9IK2zpArE4Ea5hgJVonWoXbWi266xCsYYIVrfRFP39bj+/M/hG2IVjDBCteWQDmbuuIofpOsIYJFk9FX7rM29ZBO/WDYA0TLF6ZGq77t3XgL6reEKxhgsU5M8q1eFvvf9CAWJ0I1jDB4iP3NGR0W5d/jKOH6jvBGiZYDPmgKpc8PCz9/7H1kVsRrGGCxU2GgsWorTvQhmBRZesIdLJ1B9oQLKp4GY7ypIYJFlW8DEd5UsMEiypehqM8qWGCRRUvw1Ge1DDBooqX4ShPaphgUcXLcJQnNUywqOJlOMqTGiZYVPEyHOVJDRMsqngZjvKkhgkWVbwMR3lSwwSLKl6GozypYYLF3bb+RjySbL32TWVfP8vWEybJ1mvfVPb1xcK3tYAnNSz8UWVfXyx8Wwt4UsPCH1X29cXCt7WAJzUs/FFlX18sfFsLeFLDwh9V9vXFwre1gCc1LPxRZV9fLHxbC3hSw8IfVfb1xcK3tYAnNSz8UWVfXyx8Wwt4UsPCH1X29cXCt7WAJzUs/FFlX18sfFsLeFLDwh9V9vXFwre1gCc1LPxRZV9fLHxbC3hSw8IfVfb1xcK3tYAnNSz8UWVfXyx8Wwt4UsPCH1X29cXCt7WAJzUs/FFlX18sfFsLeFLDwh9V9vXFwre1gCc1LPxRZV9fLHxbC3hSw8IfVfb1xcK3tYAnNSz8UWVfXyx8Wwt4UsPCH1X29cXCt7WAJzUs/FFlX18sfFsLeFLDwh9V9vXFwre1gCc1LPxRZV9fLHxbC3hSw8IfVfb1QCuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbgjXJH18/PTz88vWP2z8gxB9fPz+cHsTv5//rbw+vrfup7cvIk8gbVfQk5vn7l+/D+uXv2z4gxJ+ffrwEP/129gME68X1J5E4quhJTPP3p1fTOredqx8Q4tc3L8Jfz3zEb4L14uqTiBxV9CSm+fx6Wr/c8gEZfn9463/vP+SLYL24+iQiRxU9iVn+OC3my5//fBX/vKE/l39AiNOXBJ/+99c///TXv19rfTv7Ie//p4muPYnMUQnWBJ//e3tz+j+LX5d/QIbTF1ifvr93+fPsm8LTzzR/WvnT2qerTyJzVIJ1v2+vvmR//c/jHxDiy5svBH4790bmlLEvq35Se3XtSYSOKuTMUn++/r9vX898dX71A0K8fVmdfZWdKnb2Vw/jXHsSoaMSrPu9mda5nV39gEzngvUl5pV3zbUnEToqwbrfaVp/ff8Pf535Sv7qB0T6du4t4fPPNH/73y8PD58zXoAXXXsSoaMSrPsJ1k1+O/OT7s8/0/zHy+8v+hT0G7jfufokQkclWPf79OatzZk3Olc/INKn16+4F6efmHn12yGDi3X1SYSOKuPKWg/XpnP1AxJ9PvdFwdvf3f0+aDmuPonQUWVcWUuwbnD6jaOf3n0/yfPv7v56enH+/duncz/HFePqkwgdVcaVtQRruV8vv8/5/j/9+yH5VwyvPonQUWVcWUuwFvs88jNUp6id+/7oPOeeROioMq6sJVhLDfXq+Ze+Pq/x6ezeuScROqqMK2u9+/Wan78F7OoHZPn3j0UZ+BXAmFfhVWeeROioDOJ+fh/WIs+/uWjo9wmsA7QAAAR5SURBVFgJ1ndnnkToqAzifoK1xPNPt7//9cFzBOs7wfrOIO53+i0zP/6E8tMfoXLuewk//IAgzz99NfYzU2e/dyfSuScROirBup8/rWHcc68u/9Lf6b/98YfWvXlsaa4+idBRCdYEr79gfzO04Q9IceG3X/3w5fV//eXjjz22608ic1SCNcGrP/Dxy9k3PFc/IMTzN8h99O02v7/61a4/Yn7l65zrTyJzVII1wfOfrv35x5+u/f6rgqsfEOKXh2vfHvj8XSin5/P38zfTZfy8zFlXn0TmqARrhh9/P9zD658dffU1+/kPSPPHwzv//hf//dPb7/kN+arhrAtPIn1UgjXD67/z8r9fsX+1rfMfkObr9WC9+ZsLk3t16Umkj0qwpvjrxx9d9Orv4H21rfMfkOb1X/x5KVivvrKI/RXCF2efRPqoBGuSP77+s55PX1//TMLrbZ39gDTve3UmWE/ffn/+OeRfc/8srO/OPYn0UQkW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYLGhrw8PD798/w9//vMfHv7a8tNh9wSLLX36p1H/O/PPcJZgsaXfT19V/f38j2++2oKzBItNffknU59P/+ANIQMEi019O3XqtydvCBkiWGzrt1OxvnlDyBDBYmOf/0nVF28IGSJYbOyvhxfeEHKVYLG1//3bK28IuU6w2NwvD94QMkaw2Nwfp179uvVnQQeCxea+ekfIIMFia3/6OXdGCRZb+/Tyq4R+EourBIuNnd4QfvrxHTrwEcFiW89vCP/89unlO3TgI4LFtl6+h/DlO3TgQ4LFpn58D+Hzd+hs/dmwd4LFlv77HsK/T//0+9afDzsnWGzp1R8qc/oOnU/bfjbsnmCxoTd/qMzpO3S+bvnZsH+CBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAb/w8THzzbxUDJMAAAAABJRU5ErkJggg=="
     },
     "metadata": {
      "image/png": {
       "height": 480,
       "width": 600
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f <- function(x){x^2 + 1}\n",
    "\n",
    "a <- 0\n",
    "b <- 5\n",
    "\n",
    "number_of_rectangles <- 2\n",
    "  \n",
    "xgrid <- seq(a, b, length.out=number_of_rectangles+1)\n",
    "xs <- xgrid[1:number_of_rectangles]\n",
    "ys <- f(xs)\n",
    "\n",
    "plot(f,from = -1, to = 6 , col=\"red\", axes=FALSE)\n",
    "axis(1, pos=0, at=xgrid)\n",
    "axis(2, pos=0, labels=FALSE, las=2)\n",
    "abline(v=c(a,b))\n",
    "\n",
    "xlefts <- xgrid[1:number_of_rectangles]\n",
    "xrights <- xgrid[2:(number_of_rectangles+1)]\n",
    "rect(xlefts, 0, xrights, ys, col=\"lightblue\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f95f873-fb35-4b4c-a6a1-19dd2fa9b6a2",
   "metadata": {},
   "source": [
    "You may find that in the case of two rectangles you may find it is easier to do the below calculation by hand rather than by using a computer. However, the beauty of numerical methods is that once you do the code carefully, you can easily customize it to work in situations where doing it by hand will be much less efficient. In what follows we will write the code for computing the area of the two rectangles above which will be easy to customize to use any number of rectangles.\n",
    "\n",
    "To find the area of these rectangles follow the steps below.\n",
    "\n",
    "1. Define a function implementing `f(x)` \n",
    "2. (optionally) Define variables `a` and `b` storing the values of the lower and higher limits of integration.\n",
    "3. (optionally) Define a variable, say `number_of_rectangles` storing the value of the number of rectangles used in the approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a7935f17-ce2c-46b4-9f96-f0a998776fe0",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "f <- function(x){x^2 + 1}\n",
    "\n",
    "a <- 0\n",
    "b <- 5\n",
    "\n",
    "number_of_rectangles <- 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a124272a-0e59-47fb-ada8-02637389a07d",
   "metadata": {},
   "source": [
    "    Now we will need the widths and the heights of the rectangles. \n",
    "\n",
    "4. Proceed by creating a grid, say `xgrid`, of the horizontal endpoints of the rectangles. Since they are of equal width we can use the `seq()` function. Notice that the grid will have `number_of_rectangles+1` elements (for example, in the case of two rectangles the endpoints are 0, 2.5, and 5 as seen above)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5a59f187-c36f-4ba9-bf0f-5aa0573303dd",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>0</li><li>2.5</li><li>5</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0\n",
       "\\item 2.5\n",
       "\\item 5\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0\n",
       "2. 2.5\n",
       "3. 5\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 0.0 2.5 5.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xgrid <- seq(a, b, length.out=number_of_rectangles+1)\n",
    "xgrid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cade1fb-1453-4faa-a52b-b17d7790078a",
   "metadata": {},
   "source": [
    "5. Store the width of the rectangles in a variable, say `dx`. The width of each rectangle is simply the distance between `b` and `a` divided by the number of rectangles. For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a85af519-6eff-43be-90da-54a451cea107",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "2.5"
      ],
      "text/latex": [
       "2.5"
      ],
      "text/markdown": [
       "2.5"
      ],
      "text/plain": [
       "[1] 2.5"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dx <- (b-a)/number_of_rectangles\n",
    "dx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e0a3048-553a-4f1b-aa7b-7eb1fa0fbd10",
   "metadata": {},
   "source": [
    "    so the width of the rectangles when we use two rectangles is 2.5.\n",
    "\n",
    "6. Now we need the height of the rectangles. In general, the height of each rectangle equals the value of the function `f` evaluated at the left endpoint of the rectangle. \n",
    "\n",
    "- Let a vector, say `xs`, store the left endpoints of the rectangles. This vector will consist of the first `number_of_rectangles` elements of `xgrid` above (we need to omit the final element of `xgrid` because no rectangle has it as a left endpoint)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5d90367c-c46b-464a-9df9-6105adf090ba",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>0</li><li>2.5</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0\n",
       "\\item 2.5\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0\n",
       "2. 2.5\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 0.0 2.5"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xs <- xgrid[1:number_of_rectangles]\n",
    "xs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f24cfcbc-d69e-4a8a-b17d-6433726cea5d",
   "metadata": {},
   "source": [
    "so in the case of two rectangles, the left endpoints are 0 and 2.5.\n",
    "    \n",
    "- Next, evaluate the heights of the rectangles and store them in a vector, say `ys`. The heights of the rectangles are simply the value of the function `f` evaluated at the left endpoints of the rectangles (i.e. at `xs`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "405afb65-0f7d-434a-8042-76d616ec37db",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>1</li><li>7.25</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 1\n",
       "\\item 7.25\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 1\n",
       "2. 7.25\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 1.00 7.25"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ys <- f(xs)\n",
    "ys"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "540043d9-fe0a-4f5e-8953-86d9ff1415bd",
   "metadata": {},
   "source": [
    "So in the case of two rectangles,\n",
    "\n",
    "- the height of the first rectangle, `ys[1]`, is 1 ($f(0)=1$)\n",
    "\n",
    "- the height of the second rectangle, `ys[2]`, is 7.25 ($f(2.5)=7.25$)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "875201bf-ee42-4c59-a4db-874bd6879d4f",
   "metadata": {},
   "source": [
    "Finally the sum of the areas of the rectangles is simply \n",
    "\n",
    "$$\n",
    "f(0)\\Delta x + f(2.5)\\Delta x = (f(0)+f(2.5)) \\Delta x\n",
    "$$\n",
    "\n",
    "where $\\Delta x$ is the width of the rectangle. \n",
    "    \n",
    "In terms of the names of the variables we defined in R this is simply `(f(0)+f(2.5))*dx` \n",
    "\n",
    "- or alternatively `(ys[1]+ys[2])*dx` where `ys[i]` denotes the `i`'th element of the vector `ys`\n",
    "\n",
    "- or even more conveniently `sum(ys)*dx`  (using the `sum()` function which returns the sum of all elements of a vector). This option is convenient because it works with the same syntax irrespective of the number of rectangles we have.\n",
    "\n",
    "So the sum of areas of the rectangles, say `sum_areas`, is simply"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9c9eb3a1-dffe-4ee4-bada-6c17e783aee4",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "20.625"
      ],
      "text/latex": [
       "20.625"
      ],
      "text/markdown": [
       "20.625"
      ],
      "text/plain": [
       "[1] 20.625"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sum_areas <- dx*sum(ys)\n",
    "sum_areas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "734bb789-ffe0-4d97-a731-16566fce9622",
   "metadata": {},
   "source": [
    "Now for clarity let's reproduce the whole code we needed for this calculation. Note that it's not a lot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "edfdeec2-b407-43dc-8551-4b19e82b26f4",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "20.625"
      ],
      "text/latex": [
       "20.625"
      ],
      "text/markdown": [
       "20.625"
      ],
      "text/plain": [
       "[1] 20.625"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f <- function(x){x^2 + 1}\n",
    "\n",
    "a <- 0\n",
    "b <- 5\n",
    "\n",
    "number_of_rectangles <- 2\n",
    "xgrid <- seq(a, b, length.out=number_of_rectangles+1)\n",
    "dx <- (b-a)/number_of_rectangles\n",
    "xs <- xgrid[1:number_of_rectangles]\n",
    "ys <- f(xs)\n",
    "sum_areas <- dx*sum(ys)\n",
    "\n",
    "sum_areas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57a3e3b4-b838-47df-abeb-c7b88df4863f",
   "metadata": {},
   "source": [
    "Note that when we approximate the area with the area of two rectangles only we get 20.625 while the true value of the integral is 46.(6). The large approximation error is not surprising given the graph that we saw before - we are clearly omitting quite a bit of the area under the curve. \n",
    "\n",
    "To improve the approximation let's increase the number of rectangles to 5. Graphically, this would amount to approximating the area under the curve by the sum of areas of the five rectangles below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "51d6afd0-88f4-4d34-adc4-2888b7eb204c",
   "metadata": {
    "tags": [
     "hide-input"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABLAAAAPACAMAAAASenCyAAAANlBMVEUAAABNTU1oaGh8fHyMjIyampqnp6et2OaysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD/AAD///+9zLTcAAAACXBIWXMAACToAAAk6AGCYwUcAAAgAElEQVR4nO3d2WIbx5JFUVKU3NfWRP7/z7ZATZyJISOzos5aL+370E4EENhGkQRwdQfQxNXqGwBwLMEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BAtoQLKANwQLaECygDcEC2hAsoA3BKnTl3mW08KXKnr5Y+G5RIXypsqcvFr5bVAhfquzpi4XvFhXClyp7+mLhu0WF8KXKnr5Y+G5RIXypsqcvFr5bVAhfquzpi4XvFhXClyp7+mLhu0WF8KXKnr5Y+G5RIXypsqcvFr5bVAhfquzpi4XvFhXClyp7+mLhu0WF8KXKnr5Y+G5RIXypsqcvFr5bVAhfquzpi4XvFhXClyp7+mLhu0WF8KXKnr5Y+G5RIXypsqcvFr5bVAhfquzpi4XvFhUWLdXt7ZJjn/GUKiRYDCdYVBEshluzVFvplWBVEiyGW7JUt4KVQLAYblGwFhz6Ik+pQoLFcCuWaju9EqxKgsVwC5ZqOxeEglVKsBhuSbCmH/kqT6lCgsVw85dqS70SrEqCxXDTl2pLF4SCVUqwGG5BsCYf+CZPqUKCxXCzl2pbvRKsSoLFcJOXalsXhIJVSrAYbnqwph73Lk+pQoLFcHOXamu9EqxKgsVwU5dqaxeEglVKsBhu5lJtr1eCVUmwGG5usOaddSRPqUKCxXATl2qDvRKsSoLFcPOWaoMXhIJVSrAYbtpSbbJXglVJsBhuYrAmHXQST6lCgsVws5Zqm70SrEqCxXCTlmqbF4SCVUqwGG7OUm21V4JVSbAYblawZpxyBk+pEa5es/qGsTtTlmqzvRKsIQSLWWYs1WYvCAWrlGAx3ISl2nCvBKuSYDFc/VJtuVeCVUmwGK58qTbdK8GqJFgMNyFYxQdcxFOqkGAxXPVSbbtXglVJsBiueKm2fUEoWKUEi+Fql2rrvRKsSoLFcKVLtfleCVYlwWK44mAV/suH8JQqJFgMV7lU2++VYFUSLIYrXKrtXxAKVinBYri6perQK8GqJFgMV7ZULXolWJUEi+EKg1X0Lx7KU6qQYDFc1VL16JVgVRIshitaqh4XhIJVSrAYrmapuvRKsCoJFsOVLFWbXglWJcFiuKJgFfxLS3hKFRIshqtYqj69EqxKgsVwBUvV54JQsEoJFsONX6pOvRKsSoLFcMOXqlWvBKuSYDHc6KXq1SvBqiRYDDc+WGP/fcU8pQoJFsMNXqpmvRKsSoLFcGOXqtkFoWCVEiyGG7pU7XolWJUEi+FGLlW/XglWJcFiuLHBGvfvmsRTqpBgMdzApWrYK8GqJFgMN26pGl4QClYpwWK4YUvVsleCVUmwGG7UUvXslWBVEiyGGxesMf+eyTylCgkWww1aqqa9EqxKgsVwY5aq6QWhYJUSLIYbslRteyVYlQSL4QYFa8C/ZAlPqUKCxXAjlqpvrwSrkmAx3ICl6ntBKFilBIvhLl+qzr0SrEqCxXAjgjXidiziKVVIsBju4qVq3SvBqiRYDHfpUrW+IBSsUoLFcBcuVfNeCVYlwWK4i4M16HYs4ilVSLAY7rKl6t4rwaokWAx30VJ1vyAUrFKCxXCXLFX/XglWJcFiuAuWage9EqxKgsVw5y/VHnolWJUEi+HOXqpd9EqwKgkWw10QrKG3YxFPqUKCxXDnLtU+eiVYlQSL4c5cqn1cEApWKcFiuPOWai+9EqxKgsVwZy3VbnolWJUEi+HOWar99EqwKgkWw52xVDvqlWBVEiyGOytYBbdjEU+pQoLFcKcv1Z56JViVBIvhTl6qPV0QClYpwWK4U5dqX70SrEqCxXAnLtXOeiVYQ1y9ZvUNY3dODlbR7VjEU2oEwWKW05Zqb70SrEqCxXAnLdXeLggFq5RgMdwpS7W/XglWJcFiuBOWaoe9EqxKgsVwJwWr8HYs4ilVSLAY7vil2mOvBKuSYDHc0Uu1xwtCwSolWAx37FLts1eCVUmwGO7IpdpprwSrkmAx3NHBKr4di3hKFRIshjtuqfbaK8GqJFgMd9RS7fWCULBKCRbDHbNU++2VYFUSLIY7Yql23CvBqiRYDPf+Uu25V4JVSbAY7t2l2nWvBKuSYDHce0u1714JViXBYrj3gzXndiziKVVIsBjunaXaea8Eq5JgMdzbS7XzC0LBKiVYDPfmUu2+V4JVSbAY7q2l2n+vBKuSYDHcG0sV0CvBqiRYDPf6UiX0SrAqCRbDvbpUEb0SrEqCxXCvLVVGrwSrkmAx3OvBmns7FvGUKiRYDPfKUoX0SrAqCRbDvbxUIReEglVKsBjuxaWK6ZVgVRIshntpqXJ6JViVBIvhXliqoF4JViXBYrjnS5XUK8GqJFgM91KwVtyORTylCgkWwz1bqqheCVYlwWK4p0sVdUEoWKUEi+GeLFVYrwSrkmAx3OOlSuuVYFUSLIZ7GqxVt2MRT6lCgsVwj5YqrleCVUmwGO7hUsVdEApWKcFiuAdLFdgrwaokWAz3d6kSeyVYlQSL4f4sVWSvBKuSYDHc76XK7JVgVRIshvu1VKG9EqxKgsVwP5cqtVeCVUmwGO5+qWJ7JViVBIvhDkuV2yvBqiRYDPdjqYJ7JViVBIvhrq6SeyVYlQSL4Q7BWn0bFvKUKiRYDPcjWKtvwkqeUoUEi+HClyp7+mLhu0WB2/Clyp6+WPhuMd6tYFEmfLcY7lawVt+APQvfLYa7vU1fquzpi4XvFqMd/qAhfKmypy8WvlsMditYglUpfLcY6+dfuIcvVfb0xcJ3i6F+vSMnfKmypy8WvluM9PsdhOFLlT19sfDdYqA/73gOX6rs6YuF7xYD/XnHc/hSZU9fLHy3GOfvJzSEL1X29KNcvWb1DWMfbgXrl+zpRxEsKj38yL7wpcqevlj4bjHIo48YDV+q7OmLhe8WYzz+SOTwpcqevlj4bjHEk49wD1+q7OmLhe8WIzz9yonwpcqevlj4bjHAs6/ICV+q7OmLhe8Wl3v+lV7hS5U9fbHw3eJiL3wFYfhSZU9fLHy3uNRLX5kavlTZ0xcL3y0u9OJXPIcvVfb0xcJ3i8u8/JX04UuVPX2x8N3iIi/3Kn2psqcvFr5bXOKVXqUvVfb0xcJ3iwu81qv0pcqevlj4bnG+V3uVvlTZ0xcL3y3O92qv0pcqe/pi4bvF2V7vVfpSZU9fLHy3ONfrF4TxS5U9fbHw3eJMb/Uqfamypy8Wvluc581epS9V9vTFwneLs7zdq/Slyp6+WPhucY53epW+VNnTFwvfLc7wXq/Slyp7+mLhu8Xp3u1V+lJlT18sfLc43bu9Sl+q7OmLhe8WJ3u/V+lLlT19sfDd4lTvXxDGL1X29MXCd4sTHdOr9KXKnr5Y+G5xmqN6lb5U2dMXC98tTnJcr9KXKnv6YuG7xSmO7FX6UmVPXyx8tzjBsb1KX6rs6YuF7xZHuz26V+lLlT19sfDd4lgn9Cp9qbKnLxa+WxzphFzFL1X29MXCd4vjnNSr9KXKnr5Y+G5xlNN6lb5U2dMXC98tjnFir9KXKnv6YuG7xRFO7VX6UmVPXyx8t3jfyb1KX6rs6YuF7xbvOr1X6UuVPX2x8N3iXaf3Kn2psqcvFr5bvOeMXqUvVfb0xcJ3i3eccUEYv1TZ0xcL3y3edlav0pcqe/pi4bvFm87rVfpSZU9fLHy3eMuZvUpfquzpi4XvFm84t1fpS5U9fbHw3eJ1Z/cqfamypy8Wvlu85pTPv3oqfKmypy8Wvlu84pJepS9V9vTFwneLl12Sq/ilyp6+WPhu8aLLepW+VNnTFwvfLV5yYa/Slyp7+mLhu8ULLu1V+lJlT18sfLd47uJepS9V9vTFwneLZy7vVfpSZU9fLHy3eGpAr9KXKnv6YuG7xVMDepW+VNnTFwvfLZ4Y0av0pcqevlj4bvHYiAvC+KXKnr5Y+G7xyJhepS9V9vTFwneLhwb1Kn2psqcvFr5bPDCqV+lLlT19sfDd4q9hvUpfquzpR7l6zeobxjaM65VgcTnB4i0DeyVYlAnfLX4Z2av0pcqevlj4bvHT0F6lL1X29MXCd4t7Y3uVvlTZ0xcL3y0OBvcqfameTv/5y/clt2OXwneLu/G9Sl+qp9PfXF0p1ijhu0VBr9KX6un0v38V71fyA7gP043vVfpSCVYh92G2i75/8DXhS/VCsL79+r8Lbs3OuA+jlfQqfameTv/h6urmUCzBGsB9mKwkV/FL9XT6f159l4mEncw9FqyoV+lL9XT6r4I1jnssV1Wv0pfq2fT/CtYw7rFYZb1KX6rn03/5KFiDuMdS1fUqfalem16eBnAfhirsVfpSCVYh92Gmyl6lL5VgFXIfRirtVfpSZU9fLHy3QtX2Kn2psqcvFr5bmYp7lb5U2dMXC9+tSNW9Sl+q7OmLhe9WovJepS9V9vTFwncrUH2v0pcqe/pi4buVZ0Kv0pcqe/pi4buVpubjZJ4KX6rs6YuF71aYOb1KX6rs6YuF71aWObmKX6rs6YuF71aUWb1KX6rs6YuF71aSab1KX6rs6YuF71aQeb1KX6rs6YuF71aOib1KX6rs6YuF71aKSb8e/CV8qbKnLxa+WyHm9ip9qbKnLxa+Wxnm5ip+qbKnLxa+WxFm9yp9qbKnLxa+Wwmm9yp9qbKnLxa+WwHm9yp9qbKnLxa+W/u3oFfpS5U9fbHw3dq9Fb1KX6rs6YuF79bO3d6+8SXplVYPvlT29MXCd2vfVuUqfKmypy8Wvlu7dnh99X8rhC9V9vTFwndrz24Fa43s6YuF79aOHX7cLlgrZE9fLHy39uv+14OCtUL29MXCd2u3fv45g2CtkD19sfDd2qtff34lWCtkT18sfLd26vefiwrWCtnTFwvfrX368+ftgrVC9vTFwndrl/6+HUewVsievlj4bu3Rg7cPCtYK2dMXC9+tHXr4dmfBWiF7+mLhu7U7jz+9XbBWyJ6+WPhu7c2Tb5sQrBWypy8Wvls78/TTrwRrhezpi4Xv1r48+7Q+wVohe/pi4bu1Jy98+aBgrZA9fbHw3dqRl74sVbBWyJ6+WPhu7ceLH94uWCtkT18sfLd24+UvmxCsFbKnLxa+W3vxypfjCNYK2dMXC9+tfXjpx1f3BGuF7OmLhe/WLrzaK8FaInv6YuG7tQdvfFeqYK2QPX2x8N3agbe+21mwVsievlj4bvX35nfRC9YK2dOP4kt6d+nNXgnWEtnTjyJYe/R2rwRriezpi4XvVm+v/3rwF8FaIXv6YuG71dq7vRKsJbKnLxa+W529myvBWiN7+mLhu9XYEb0SrCWypy8Wvlt9HdMrwVoie/pi4bvV1fs/vronWCtkT18sfLeaOrJXgrVE9vTFwnerpyNzJVhrZE9fLHy3Wjq6V4K1RPb0xcJ3q6FjLwcPBGuF7OmLhe9WP6f0SrCWyJ6+WPhutXNKrgRrjezpi4XvVjMnvby6E6w1sqcvFr5bvZzaK8FaInv6YuG71cqpuRKsNbKnLxa+W42c/PLqTrDWyJ6+WPhu9XFOrwRriezpi4XvVhvn5Eqw1sievlj4bjVx1surO8FaI3v6YuG71cO5vRKsJbKnLxa+Wy2cmyvBWiN7+mLhu9XB+b0SrCWypy8Wvlvbd/bl4IFgrZA9fbHw3dq8i3olWEtkT18sfLe27qJcCdYa2dMXC9+tbbvs5dWdYK2RPX2x8N3atIt7JVhLZE9fLHy3tuziXAnWGtnTFwvfre26/OXVnWCtkT19sfDd2qwhvRKsJbKnLxa+W1s1JFeCtUb29MXCd2ujBvVKsJbInr5Y+G5t06heCdYS2dMXC9+tLRrz46t7grVC9vTFwndrgwb2SrCWyJ6+WPhubc/AXAnWGtnTFwvfrc0Z2ivBWiJ7+mLhu7UxIy8HDwRrhezpi4Xv1raM7pVgLZE9fbHw3dqS4bkSrDWypy8WvlsbUtArwVoie/pi4bu1GRW5Eqw1sqcvFr5bW1HTK8FaInv6YuG7tRE1uRKsNbKnLxa+W5tQ9PLqTrDWyJ6+WPhubUFdrwRriezpi4Xv1gbU5Uqw1sievlj4bi1X+PLqTrDWyJ6+WPhurVbbK8FaInv6YuG7tVhtrgRrjezpi4Xv1lLFL6/uBGuN7OmLhe/WSvW9EqwlsqcvFr5bC9XnSrDWyJ6+WPhuLTPh5dWdYK2RPX2x8N1aZU6vBGuJ7OmLhe/WGpNyJVhrZE9fLHy3lpjWK8FaInv6YuG7tcC8XAnWGtnTFwvfrflm9kqwlsievlj4bs02NVeCtUb29MXCd2uyyb0SrCWypy8WvltTzc6VYK2RPX2x8N2aaX6vBGuJ7OmLhe/WPAtyJVhrZE9fLHy3plnSK8FaInv6Ua5es/qGJbi9ffXuLyZY82VPP8qrK736hgVYlivBWiF7+mLhuzXD7eH11ZJw/J9grZA9fbHw3ZrgVrDCZE9fLHy36t3/sF2wkmRPXyx8t6r9+uWgYCXJnr5Y+G4V+/3HDIKVJHv6YuG7Verv314JVpLs6YuF71ah2wd/KypYSbKnLxa+W3Ue/Wm7YCXJnr5Y+G5VefJOHMFKkj19sfDdqnH79J2DgpUke/pi4btV4VmuBCtL9vTFwnerwEufyyBYSbKnLxa+W8O9/DEygpUke/pi4bs12isfeyVYSbKnLxa+W2O9+il9gpUke/pi4bs11OufKipYSbKnLxa+WwO99SHIgpUke/pi4bs1zpsf2i5YSbKnLxa+W8O8/R0TgpUke/pi4bs1yHvfiSNYSbKnLxa+W2O8+x1egpUke/pi4bs1whFfOShYSbKnLxa+W5d74Z2DzwlWkuzpi4Xv1sWO+0ZnwUqSPX2x8N260LFfQC9YSbKnLxa+Wxc56mrwnmAlyZ6+WPhuXeD4XAlWluzpi4Xv1tlOyZVgZcmevlj4bp3ptFwJVpbs6YuF79ZZTs2VYGXJnr5Y+G6d4fRcCVaW7OmLhe/Wyc7JlWBlyZ6+WPhuneqsXAlWluzpi4Xv1mnOzJVgZcmevlj4bp3ivKvBe4KVJHv6YuG7dbwLciVYWbKnLxa+W0e7JFeClSV7+mLhu3Wky3IlWFmypy8WvlvHubRXghUle/pi4bt1jItzJVhZsqcvFr5b77voh+2/CVaS7OmLhe/We4bkSrCyZE9fLHy33jYoV4KVJXv6YuG79ZZhuRKsLNnTFwvfrdcNzJVgZcmevlj4br1maK4EK0v29MXCd+sVY3MlWFmypy8WvlsvGp0rwcqSPX2x8N16weCrwXuClSR7+mLhu/VMRa4EK0v29MXCd+uJmlwJVpbs6YuF79YjVbkSrCzZ0xcL360H6nIlWFmypy8Wvlt/VOZKsLJkT18sfLd+qc2VYGXJnr5Y+G7dq86VYGXJnr5Y+G7dzciVYGXJnr5Y+G5NyZVgZcmevlj0bt3OyZVgZcmevljwbk3LlWBlyZ6+WOxuzavVnWBlyZ6+WOhuTc2VYGXJnr5Y5G5NzpVgZcmevljgbk3PlWBlyZ5+lKvXrL5hk716P9RaEg7BWiJ7+lFefSqtvmEz3U5M1GNLwiFYS2RPXyxpt35cCS57Bi85VrCWyJ6+WMxu/fyrK8Gac+zqR3up7OmLhezW7z8SFaw5x65+vJfKnr5Ywm49+JN2wZpz7OJHfK3s6Yvtf7cevQNHsOYcu/QRXy17+mJ7360nf3QlWHOOXfeAb0D29MV2vVvP394sWHOOXfWIb0L29MV2vFsvfRqDYM05ds0jvhHZ0xfb62698tkxgjXn2AWP+HZkT19sn7v16kddCdacY6c/4luSPX2xPe7WG+9uFqw5x859wDcme/pi+9utNz+MQbDmHDvzAd+c7OmL7W233vnsGMGac+y8B3yDsqcvtqvdev9T2gVrzrGzHvFNyp6+2H526/aYL5UQrDnHznnINyp7+mI72a2janUnWLOOnfCQb1f29MX2sFvH1upOsGYdW/2Qb1r29MX679ZJXy8oWHOOLX3Ety57+mLNd+vUL0MVrDnH1j3iDWRPX6zzbp1wKfibYM05tuohbyF7+mJ9d+usb5oXrDnHljziXWRPX6zpbp1VqzvBmnXs+Ee8kezpi/XbrdvbMy4FfxOsOceOftBbyZ6+WK/dur2oVneCNevYoY96N9nTF2u0Wxe26p5gzTl21GPeUvb0xdrs1uWxOhCsOceOeMTbyp6+WI/dGlOrO8GadeyAh6qv7OmLNditAVeCfwjWnGOHPFhdZU9fbOu7NbJWd4I169hRj1dL2dMX2/RuDa7VnWDNOnbgQ9ZP9vTFNrtbI34n+JxgzTl27KPWTPb0xTa5W5f+tdXrBGvOscMfuE6ypy+2ud2qi9WBYM05tuTB6yJ7+mKb2q3aWB0I1pxjyx7ADrKnL7aV3br4TTfHEaw5x5Y+iFuXPX2xLezWpFgdCNacY8sfyC3Lnr7Y6t2a16p7gjXn2DmP5kZlT19s5W5NjtWBYM05dt4jukHZ0xdbtVu3V4ssegYvOVawlgiZfvJrjV9W7NaPOVflSrDmHDt/qTYkZPoFF0h304P1e8a4Z/CSYwVriaDp50dr4m49/GVg3DN4ybGCtUTY9HOjNWm3ng4V9wxecqxgLRE4/e0jlSdN2K2X5oh7Bi85VrCWCJ1+TrQKd+ut2x/3DF5yrGAtkT198TXi+N26febFY5c8kwRr0rGjl6qV7Ol/Knu1NXa3jn5VGPcMXnKsYC2RPf1fR710OdW43TrpVsU9g5ccK1hLZE//3NBwXb5bZ92SuGfwkmMFa4ns6d8wolwn79bzQ885PO4ZvORYwVoie/ojXNKQY3fr9TPOen0X9wxecqxgLZE9/SneqMprrq5O/f8YckvjnsFLjhWsJbKnv8RRwUp6E7JgTTp29eYvlT19McHa8bGCtUT29MU8g3d8rGAtkT19Mc/gHR8rWEtkT1/MM3jHxwrWEtnTF/MM3vGxgrVE9vTFPIN3fKxgLZE9fTHP4B0fK1hLZE9fzDN4x8cK1hLZ0xfzDN7xsYK1RPb0xTyDd3ysYC2RPX0xz+AdHytYS2RPX8wzeMfHCtYS2dOPsuxNgwRave1LZU8/yuoVJsnqbV8qe/pii3Zr1Uobd8fHbkX29MXCVtq4Oz52K7KnLxa20sbd8bFbkT19sbCVNu6Oj92K7OmLha20cXd87FZkT18sbKWNu+NjtyJ7+mJhK23cHR+7FdnTFwtbaePu+NityJ6+WNhKG3fHx25F9vTFwlbauDs+diuypy8WttLG3fGxW5E9fbGwlTbujo/diuzpi4WttHF3fOxWZE9fLGyljbvjY7cie/piYStt3B0fuxXZ0xcLW2nj7vjYrcievljYSht3x8duRfb0xcJW2rg7PnYrsqcvFrbSxt3xsVuRPX2xsJU27o6P3Yrs6YuFrbRxd3zsVmRPXyxspY2742O3Inv6YmErbdwdH7sV2dMXC1tp4+742K3Inh5oRbCANgQLaEOwgDYEC2hDsIA2BAtoQ7CANgQLaEOwgDYEC2hDsIA2BAtoQ7CANgQLaEOwgDYEC2hDsIA2BAtoQ7CANgQLaEOwgDYEC2hDsIA2BAtoQ7CANgQLaEOwgDYEC2hDsIA2BAtoQ7CANgQLaEOwqnz+dH119eHT5wVHf7ia/rB+/nRzdRj3v7nHfv3nx6zXS+7kH75cTb2jv189NPHgLUmdu9q3D78X68O32Wf/N32dv1z/eR5d/zvx3I+/T72ZficfXM+9o78IlmAV+Xb9YLUmP5m+TX4e3d398+iZ9M+0c2/+Hnq9oFifJnfjX8ESrCI3D1frw9Sjf7Zy5on/XT32v0nnflp3Jx98md2Nj4/u5YkHb0nq3LU+Hzbq45cfS31fri8Tj/712m7iifdXRtf/+3p3+JHS/eHfpxz79XDUpx/Hfrsv1/SfY02/o6+n3bUbJlgVbv5eGR3+s/hp3sm/f5o078T7F1h/rsi+zLso/OfvPXu4kz9OOfSvT7Pv6MPP3K/nHbdRglXg+4NN/j51qw/PopsFVyp/X0P+O+3y7MPfFxzf5l8j/f4J+NwTZ1d5ewSrwJeHr6o+TbwmPLy8urmb/Ox9fNzcQL98G2a4nh6sw38KZv4KdpsEq8Cj1Zq5Zz9/4r32R7KrgjX3p+6H/wx9mDvqx4n/5dsuwSpwWK2vv//H1wxlqPcAAAXQSURBVImv5K8+HI5dGqzvK35j92XaD84eHHj1de4dff8z9+//+5HJm+AXWoJVYFmwfv6d+dJg/Tu7HQc3s39/dn3/WnbqHX3/M/fPv65Erxf9bf96glXg8V9uTu/H0mBdP4z1HPd/OzL1Gfzp58vIqXf04UXdgz9HTi2WYBW4yg3WzfRfZd3/OeWHqY38eUE4+Y5+/Hfu0/+rsBWCVSA3WIe/jZr8Lpn7v+K4mfqm658XhJPv6PswH/5O9u7bv9crflC4DYJVIDZY/yy4WPnw6Z+bua+xPv3OxdQ7+vrBXXv/d2eZvzEUrAKpwbpZ9cOVbzcTn8C/LwhX/5dh+q82NkGwCoQGa1mv7u5ff1xP+j3h7wvClcE6/Or5Zs3RiwlWgWe/JZz7FrA1z6Of77pe9curz1ezPiXizwXh0h8WLv1N8EqZUxdb9ndYPy1Z5vu/EFr350Hfp73kuHpuyrnPb8WCY9fLnLpYYLDuf9y+4lP0fps2tGAtlTl1scOfzPz5Pfvh01fmvpViwTL//NuCyYc+khWsJW+A2gTBKrDs0xp+mv8cuu/V7N9aHc7883P2r9N+UrgqWI/GfbRhUQSrwsMtfrRo00+fYsWfX/288P7z0vWfFc/gqXf0x4f38ccVd/gmCFaFBx8z+nH+tdLsYN2/y23+O0XuP+n0138Kvi5J5tQ7+r8HryE/T3s9uTmCVeH+M91v/nym++Rn0uxgfbha8862+zeoHO7bb/fvs5v/M525d/STcUM/YkawSvz5VsIVz6TJwfq86qfQTw6e/yvKuXf043c/Z/7ZqGAVefidl9N/2T85WI++bWtmsB59H+KKPwGbfEc/HDe1V4JV5Oufjy6a/83Pk59HD78zdm6w7j7/OXvJVz9PvqMfvMbK/A3hgWBV+fzp8P62T/v/D//zXk08/b/DLzU+/LPmkwtmB+vu+/24N/+EfhbWgWABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFgt9urq6+vD7f3z58T+uvq68OWyeYLHS9Y9G/e+Ff4YXCRYr/Xd4VfXt/h8fvdqCFwkWS338kambwz+4IOQIgsVS3w+d+vfOBSFHESzW+vdQrO8uCDmKYLHYzY9UfXRByFEEi8W+Xv3igpB3CRar/e9nr1wQ8j7BYrkPVy4IOY5gsdznQ6/+WX0r6ECwWO6TK0KOJFis9sXP3DmWYLHa9a/fEvohFu8SLBY7XBBe/3mHDrxFsFjr/oLwy/frX+/QgbcIFmv9eg/hr3fowJsEi6X+vIfw/h06q28NWydYrPT3PYTfDv/03+rbw8YJFis9+FCZwzt0rtfeGjZPsFjo0YfKHN6h82nlrWH7BAtoQ7CANgQLaEOwgDYEC2hDsIA2BAtoQ7CANgQLaEOwgDYEC2hDsIA2BAtoQ7CANgQLaEOwgDYEC2hDsIA2BAtoQ7CANgQLaEOwgDYEC2hDsIA2BAtoQ7CANgQLaEOwgDYEC2hDsIA2BAtoQ7CANgQLaEOwgDYEC2hDsIA2BAtoQ7CANgQLaEOwgDYEC2hDsIA2BAtoQ7CANgQLaEOwgDYEC2hDsIA2BAtoQ7CANgQLaEOwgDYEC2hDsIA2BAtoQ7CANgQLaEOwgDYEC2hDsIA2BAtoQ7CANv4fOB6zDlBoxNAAAAAASUVORK5CYII="
     },
     "metadata": {
      "image/png": {
       "height": 480,
       "width": 600
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f <- function(x){x^2 + 1}\n",
    "\n",
    "a <- 0\n",
    "b <- 5\n",
    "\n",
    "number_of_rectangles <- 5\n",
    "  \n",
    "xgrid <- seq(a, b, length.out=number_of_rectangles+1)\n",
    "xs <- xgrid[1:number_of_rectangles]\n",
    "ys <- f(xs)\n",
    "\n",
    "plot(f,from = -1, to = 6 , col=\"red\", axes=FALSE)\n",
    "axis(1, pos=0, at=xgrid)\n",
    "axis(2, pos=0, labels=FALSE, las=2)\n",
    "abline(v=c(a,b))\n",
    "\n",
    "xlefts <- xgrid[1:number_of_rectangles]\n",
    "xrights <- xgrid[2:(number_of_rectangles+1)]\n",
    "rect(xlefts, 0, xrights, ys, col=\"lightblue\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8ca8278-6f12-42d1-a6d3-f2af3f8ab309",
   "metadata": {},
   "source": [
    "Now in order to calculate this area we can use exactly the same code as bofore but we need to set the number of rectangles to 5. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "2d7062a5-2013-4fc9-98ea-fcf2f35b0ecc",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "35"
      ],
      "text/latex": [
       "35"
      ],
      "text/markdown": [
       "35"
      ],
      "text/plain": [
       "[1] 35"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f <- function(x){x^2 + 1}\n",
    "\n",
    "a <- 0\n",
    "b <- 5\n",
    "\n",
    "number_of_rectangles <- 5\n",
    "xgrid <- seq(a, b, length.out=number_of_rectangles+1)\n",
    "dx <- (b-a)/number_of_rectangles\n",
    "xs <- xgrid[1:number_of_rectangles]\n",
    "ys <- f(xs)\n",
    "sum_areas <- dx*sum(ys)\n",
    "\n",
    "sum_areas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09044040-e701-4e57-88ed-613ecf7bd5c8",
   "metadata": {},
   "source": [
    "As we can see the approximation gets better (now the area is 35 rather than 20.625) but we are still far from the true value of 46.667. \n",
    "\n",
    "Note that at this stage it becomes clear that calculating the area of these rectangles by hand will take quite a while, while the R code we use does the calculation in miliseconds.\n",
    "\n",
    "Let's take a big step and use 100 rectangles (if we were to do this by hand it might take weeks, but a computer excels at repetitive tasks)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6074aef0-97bf-4932-a93a-6b5326b9982d",
   "metadata": {
    "tags": [
     "hide-input"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABLAAAAPACAMAAAASenCyAAAANlBMVEUAAABNTU1oaGh8fHyMjIyampqnp6et2OaysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD/AAD///+9zLTcAAAACXBIWXMAACToAAAk6AGCYwUcAAAgAElEQVR4nO3d3XrjxrGGUTl/e2fHSej7v9k9ssf2zJCUCLKruwvfWidxTgYoqfA+AkhRb78ANPG2+gQAHiVYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWIXefHUZLXypsqcvFr5bVAhfquzpi4XvFhXClyp7+mLhu0WF8KXKnr5Y+G5RIXypsqcvFr5bVAhfquzpi4XvFhXClyp7+mLhu0WF8KXKnr5Y+G5RIXypsqcvFr5bVAhfquzpi4XvFhXClyp7+mLhu0WF8KXKnr5Y+G5RIXypsqcvFr5bVAhfquzpi4XvFhXClyp7+mLhu0WF8KXKnr5Y+G5RIXypsqcvFr5bVAhfquzpi4XvFhXClyp7+mLhu0WFRUt1uSw57BWXVCHBYjjBoopgMdyapdqlV4JVSbAYbslSXQQrgWAx3KJgLTjoTS6pQoLFcCuWap9eCVYlwWK4BUu1zw2hYJUSLIZbEqzph7zLJVVIsBhu/lLt1CvBqiRYDDd9qXa6IRSsUoLFcAuCNfmAH3JJFRIshpu9VHv1SrAqCRbDTV6qvW4IBauUYDHc9GBNPdynXFKFBIvh5i7Vbr0SrEqCxXBTl2q3G0LBKiVYDDdzqfbrlWBVEiyGmxusecd6kEuqkGAx3MSl2rBXglVJsBhu3lJteEMoWKUEi+GmLdWWvRKsSoLFcBODNelAh7ikCgkWw81aqj17JViVBIvhJi3VnjeEglVKsBhuzlLt2ivBqiRYDDcrWDOO8gSX1Ahv96w+MU5nylJt2yvBGkKwmGXGUm17QyhYpQSL4SYs1ca9EqxKgsVw9Uu1c68Eq5JgMVz5Um3dK8GqJFgMNyFYxQd4iUuqkGAxXPVS7d0rwaokWAxXvFR73xAKVinBYrjapdq9V4JVSbAYrnSptu+VYFUSLIYrDlbhPz6ES6qQYDFc5VLt3yvBqiRYDFe4VPvfEApWKcFiuLql6tArwaokWAxXtlQteiVYlQSL4QqDVfQPD+WSKiRYDFe1VD16JViVBIvhipaqxw2hYJUSLIarWaouvRKsSoLFcCVL1aZXglVJsBiuKFgF/2gJl1QhwWK4iqXq0yvBqiRYDFewVH1uCAWrlGAx3Pil6tQrwaokWAw3fKla9UqwKgkWw41eql69EqxKgsVw44M19t8r5pIqJFgMN3ipmvVKsCoJFsONXapmN4SCVUqwGG7oUrXrlWBVEiyGG7lU/XolWJUEi+HGBmvcvzWJS6qQYDHcwKVq2CvBqiRYDDduqRreEApWKcFiuGFL1bJXglVJsBhu1FL17JVgVRIshhsXrDH/zmQuqUKCxXCDlqpprwSrkmAx3JilanpDKFilBIvhhixV214JViXBYrhBwRrwjyzhkiokWAw3Yqn69kqwKgkWww1Yqr43hIJVSrAY7vWl6twrwaokWAw3IlgjzmMRl1QhwWK4l5eqda8Eq5JgMdyrS9X6hlCwSgkWw724VM17JViVBIvhXg7WoPNYxCVVSLAY7rWl6t4rwaokWAz30lJ1vyEUrFKCxXCvLFX/XglWJcFiuBeW6gS9EqxKgsVwzy/VGXolWJUEi+GeXqpT9EqwKgkWw70QrKHnsYhLqpBgMdyzS3WOXglWJcFiuCeX6hw3hIJVSrAY7rmlOkuvBKuSYDHcU0t1ml4JViXBYrhnluo8vRKsSoLFcE8s1Yl6JViVBIvhngpWwXks4pIqJFgMd3ypztQrwaokWAx3eKnOdEMoWKUEi+GOLtW5eiVYlQSL4Q4u1cl6JVhDvN2z+sQ4ncPBKjqPRVxSIwgWsxxbqrP1SrAqCRbDHVqqs90QClYpwWK4I0t1vl4JViXBYrgDS3XCXglWJcFiuEPBKjyPRVxShQSL4R5fqjP2SrAqCRbDPbxUZ7whFKxSgsVwjy7VOXslWJUEi+EeXKqT9kqwKgkWwz0crOLzWMQlVUiwGO6xpTprrwSrkmAx3ENLddYbQsEqJVgM98hSnbdXglVJsBjugaU6ca8Eq5JgMdznS3XmXglWJcFiuE+X6tS9EqxKgsVwny3VuXslWJUEi+E+D9ac81jEJVVIsBjuk6U6ea8Eq5JgMdzHS3XyG0LBKiVYDPfhUp2+V4JVSbAY7qOlOn+vBKuSYDHcB0sV0CvBqiRYDHd/qRJ6JViVBIvh7i5VRK8Eq5JgMdy9pcrolWBVEiyGux+sueexiEuqkGAx3J2lCumVYFUSLIa7vVQhN4SCVUqwGO7mUsX0SrAqCRbD3VqqnF4JViXBYrgbSxXUK8GqJFgMd71USb0SrEqCxXC3grXiPBZxSRUSLIa7WqqoXglWJcFiuB+XKuqGULBKCRbD/bBUYb0SrEqCxXDfL1VarwSrkmAx3I/BWnUei7ikCgkWw323VHG9EqxKgsVw3y5V3A2hYJUSLIb7ZqkCeyVYlQSL4f5cqsReCVYlwWK4P5YqsleCVUmwGO73pcrslWBVEiyG+7pUob0SrEqCxXC/LVVqrwSrkmAx3K9LFdsrwaokWAz3vlS5vRKsSoLFcF+WKrhXglVJsBju7S25V4JVSbAY7j1Yq89hIZdUIcFiuC/BWn0KK7mkCgkWw4UvVfb0xcJ3iwKX8KXKnr5Y+G4x3kWwKBO+Wwx3EazVJ3Bm4bvFcJdL+lJlT18sfLcY7f0NDeFLlT19sfDdYrAvt4PvVp/GUtnTFwvfLcZ6exMswaoUvlsM9eXnq//5InypsqcvFr5bjHQRrHfZ0xcL3y0GugjWr7KnLxa+Wwz06xsaBEuwKoXvFuP89rhdsARriLd7Vp8Yp/D7OgmWYA0hWFT6PVWCJViVwneLQS6C9Yfs6YuF7xZjXATrT9nTFwvfLYa4CNY3sqcvFr5bjHD55g0NgiVYlcJ3iwEu374+KFiCVSl8t3jd929oECzBqhS+W7zsh6dXgiVYlcJ3i1f9+LhdsASrUvhu8aKr1wcFS7Aqhe8Wr7l+fVCwBKtS+G7xkhuvDwqWYFUK3y1ecev1QcESrErhu8ULbj69EizBqhS+Wzzv9uN2wRKsSuG7xfME647s6YuF7xZPu/P6oGAJVqXw3eJZ997QIFiCVSl8t3jS3Tc0CJZgVQrfLZ5z/w0NgiVYlcJ3i6d88LhdsASrUvhu8YyPXh8ULMGqFL5bPOHOLzwL1lfZ0xcL3y2O+/j1QcESrErhu8VxHz9uFyzBqhS+Wxz12euDgiVYlcJ3i6M+S5VgCVal8N3ioE8ftwuWYFUK3y2O+fz1QcESrErhu8UhD7yhQbAEq1L4bnHEI4/bBUuwKoXvFo977PVBwRKsSuG7xcMevBkULMGqFL5bPOrRp1eCJViVwneLBz38uF2wBKtS+G7xmEd+f1Cwvsqevlj4bvGQTz5fVLC+kz19sfDd4hFHXh8ULMGqFL5bPODQ0yvBEqxK4bvF5449bhcswaoUvlt86uDrg4IlWJXCd4tPHXx9ULAEq1L4bvGZo4/bBUuwKoXvFh87/vqgYAlWpfDd4mPHUyVYglUpfLf40BOP2wVLsCqF7xYfeeb1QcESrErhu8UHnnpDg2AJVqXw3eK+o7/wLFhfZU9fLHy3uOdy+BeeBeur7OmLhe8Wd1yefUODYAlWpfDd4rbnn14JlmBVCt8tbnrhcbtgCVal8N3illcetwuWYFUK3y1ueOlxu2AJVqXw3eLaa4/bBUuwKoXvFldefHolWIJVKXy3+NGrj9sFS7Aqhe8WPxKsl2VPXyx8t/jBy68PCpZgVQrfLb73+uN2wRKsSuG7xbdGvD4oWIJVKXy3+NaIVAmWYFUK3y2+MeRxu2AJVqXw3eJPY14fFCzBGuLtntUnxh5e/oVnwfoqe/pRBIuPjHrcLliCVSl8t/jNuNcHBUuwKoXvFr8ZlyrBEqxK4bvFrwY+vRIswaoUvlu8G/m4XbCugvWvn/+75DxOKXy3+GXw64OCdRWsv769KdYo4bvF6NcHBesqWL+/FO8l+QF8DcONfn1QsASrkq9htuE3g4J1K1j/+fq/C87mZHwNo41/eiVYV8H6y9vbX9+LJVgD+BomK3jcLlhXwfrH3d8ykbDDfMWCvf73BwXrhh+n/7dgjeMrlmvA3x8UrBuupv+nYA3jKxar5vVBwbrxTvef/yZYg/iKpSp6eiVYd381R54G8DUMVfW4XbAEq5KvYaay1wcFS7Aq+RpGqnt6JVg+raFS+G5lqnzcLliCVSl8tzJVpkqwBKtS+G5FKn16JViCVSl8txLVPm4XLMGqFL5bgYpfHxQswaoUvlt5qh+3C5ZgVQrfrTT1rw8KlmBVCt+tMBNuBgVLsCqF71aWGU+vBEuwKoXvVpQpj9sFS7Aqhe9WkjlPrwRLsCqF71aOWY/bBUuwKoXvVo5ZqRIswaoUvlsp5j29EizBqhS+WyEmPm4XLMGqFL5bGer+Oo5g3ZA9fbHw3YpwuUx83C5YglUpfLcSXKa+PihYglUpfLcCfH18JVjTZE9fLHy3zu8iWLNlT18sfLdO74/n7YI1Tfb0xcJ36+Quf/41esGaJnv6YuG7dW6Xb34fR7CmyZ6+WPhundp37xYVrGmypy8WvltndhGsNbKnLxa+Wyf2w9vbBWua7OmLhe/WeV1+eLeoYE2TPX2x8N06rauPvxKsabKnLxa+W2d1/eEMgjVN9vTFwnfrpG58moxgTZM9fbHw3TqnWx9/JVjTZE9fLHy3Tunmx18J1jTZ0xcL360z+vH1QcGaLHv6YuG7dUJ3/jyOYE2TPX2x8N06nbsf3i5Y02RPXyx8t87m/l+bEKxpsqcvFr5bJ/PBX5sQrGmypy8Wvlvncvtxu2BNlj19sfDdOpPLncftgjVZ9vTFwnfrRD75W6mCNU329MXCd+s8PvvjzoI1Tfb0xcJ36zQ+/ePOgjVN9vTFwnfrLD563C5Yk2VPXyx8t87h48ftgjVZ9vTFwnfrFD55eiVYk2VPXyx8t87gs8ftgjVZ9vTFwnfrBD593C5Yk2VPXyx8t/r7/HG7YE2WPf0ob/esPjFe8cDjdsGaLHv6UQTrjB55eiVYk2VPXyx8t3p77HG7YE2WPX2x8N1q7SJYW8qevlj4bnX26ON2wZose/pi4bvV2MOP2wVrsuzpi4XvVl+P3wwK1mTZ0xcL362ujjy9EqzJsqcvFr5bTR163C5Yk2VPXyx8t3o69rhdsCbLnr5Y+G61dPBxu2BNlj19sfDdaujwzaBgTZY9fbHw3ern+NMrwZose/pi4bvVzoMfJiNYC2VPXyx8t5q5PPO4XbAmy56+WPhu9fLIZ7cL1nLZ0xcL361Wnnx6JViTZU9fLHy3GnnizaKCtUT29MXCd6uP559eCdZk2dMXC9+tNl54eiVYk2VPXyx8t5p46WZQsCbLnr5Y+G718NrTK8GaLHv6YuG71cLTbxYVrCWypy8WvlsdvPq4XbAmy56+WPhu7e+VN4sK1hLZ0xcL363tDXh6JViTZU9fLHy3djficbtgTZY9fbHw3drbZcjjdsGaLHv6YuG7tbUXPptBsBbKnr5Y+G7tbNTjdsGaLHv6YuG7ta9xT68Ea7Ls6YuF79a2XvxsBsFaKHv6YuG7tauRT68Ea7Ls6YuF79amhj69EqzJsqcvFr5bexp7MyhYk2VPXyx8t3Y0+umVYE2WPX2x8N3a0PDH7YI1Wfb0xcJ3az/jH7cL1mTZ0xcL363tFDxuF6zJsqcvFr5bmym5GRSsybKnLxa+W3upeXolWJNlT18sfLd2UvX0SrAmy56+WPhubaTs6ZVgTZY9fbHw3dpG4c2gYE2WPX2x8N3aReXTK8GaLHv6YuG7tYmBnywqWMtlT18sfLe2MPSTRQVruezpi4Xv1g7GfrKoYC2XPX2x8N3aQPnTK8GaLHv6YuG7tVzRrzoL1kLZ0xcL363VZjy9EqzJsqcvFr5bi015eiVYk2VPXyx8t5aadDMoWJNlT18sfLdWmvX0SrAmy56+WPhuLTTtZlCwJsuevlj4bi0z8emVYE2WPX2x8N1aZebNoGBNlj19sfDdWmPu0yvBmix7+mLhu7XE5JtBwZose/pi4bu1wPSnV4I1Wfb0xcJ3a775N4OCNVn29MXCd2u2FU+vBGuy7OmLhe/WZEtuBgVrsuzpi4Xv1lSLnl4J1mTZ0xcL362ZVt0MCtZk2dMXC9+tedY9vRKsybKnLxa+W9NM+6A+wVoue/pR3u5ZfWIJpvyZCcHaRPb0owjWOnP+zIRgbSJ7+mLhuzXD6ptBwZose/pi4bs1wfKfrQRrsuzpi4XvVrkNbgYFa7Ls6YuF71a1HW4GBWuy7OmLhe9WsS2eXgnWZNnTFwvfrVKb3AwK1mTZ0xcL361C2zy9EqzJsqcvFr5bdfa5GRSsybKnLxa+W1X2eO+VYC2RPX2x8N2qsf4XcQRroezpi4XvVoUdfhFHsBbKnr5Y+G4V2O5mULAmy56+WPhuDbf8Y2QEa7ns6YuF79ZoOz69EqzJsqcvFr5bY+359EqwJsuevlj4bg216c2gYE2WPX2x8N0aaNunV4I1Wfb0xcJ3a5x9n14J1mTZ0xcL361hNn56JViTZU9fLHy3Btn6ZlCwJsuevlj4bg2x+c9WgjVZ9vTFwndrgO1vBgVrsuzpi4Xv1uv2vxkUrMmypy8Wvlsva/D0SrAmy56+WPhuvajFzaBgTZY9fbHw3XpJk6dXgjVZ9vTFwnfrBW1uBgVrsuzpi4Xv1tMa/WwlWJNlT18sfLee1OpmULAmy56+WPhuPaXZzaBgTZY9fbHw3XrC5p/KIFjLZU9fLHy3Drvs/qkMgrVc9vTFwnfrqJ0/U1SwNpE9fbHw3Tqm582gYE2WPX2x8N06ouvNoGBNlj19sfDdelzfm0HBmix7+mLhu/WwxjeDgjVZ9vTFwnfrQZfOP1sJ1mTZ0xcL363H9L4ZFKzJsqcvFr5bj/h6N7hBdASrhezpi4Xv1uf+eGP7BtERrBaypy8Wvluf+eaN7RtER7BayJ6+WPhufezy7XsZNoiOYLWQPX2x8N36yA+/5bxBdASrhezpi4Xv1n1Xv+W8QXQEq4Xs6YuF79Y9l+s3tm8QHcFqIXv6YuG7dcet33LeIDqC1UL29MXCd+um27/lvEF0BKuF7OmLhe/WDfd+y3mD6AhWC9nTFwvfrSv3P/Jqg+gIVgvZ0xcL360ffPSRVxtER7BayJ6+WPhufefjj7zaIDqC1UL29MXCd+sbn33+8QbREawWsqcvFr5bf/j88483iI5gtZA9fbHw3frqkc8/3iA6gtVC9vTFwnfrV4/9MZwNoiNYLWRPXyx8t355/I/hbBAdwWohe/pi4bt14I/hbBAdwWohe/pi0bt1OfKXUTeIjmC1kD19seDduvoAmZ2ueMHqK3v6YrG7deMDZHa64gWrr+zpi4Xu1hN/t2uD6AhWC9nTF4vcraf+6vwG0RGsFrKnLxa4W0/+1fkNoiNYLWRPP8rbPatPbLKn/4jzBtERrBaypx9FsH458q6r1Ve8YPWVPX2xpN068q6r1Ve8YPWVPX2xmN06+K6r1Ve8YPWVPX2xkN06/K6r1Ve8YPWVPX2xhN169WcrwRKsI7KnL3b+3Xr9ZyvBEqwjsqcvdvbdeu1Bu2AJ1nHZ0xc79W69/KBdsATruOzpi514twY8aBcswToue/piZ92tkT9bCZZgHZE9fbFz7tbYn60ES7COyJ6+2Bl3a9yDdsESrOOypy92vt0afzMoWIJ1RPb0xc62WxU3g4IlWEdkT1/sVLtVlCrBEqwjsqcvdp7dqvvZSrAE64js6YudZLdKUyVYgnVE9vTFzrBb1T9bCZZgHZE9fbH2uzUhVYIlWEdkT1+s+W5N+dlKsATriOzpizXerWmpEizBOiJ7+mJtd2viz1aCJVhHZE9frOlulfzyzU5XfPPTX70fS2VPX6zfbl2GfsrVrld889NfvSRLZU9frNdu/fEGBsHa+/RXL8pS2dMXa7Rbgz+Qb+8rvvnpr96VpbKnL9ZmtxamSrAE64js6Yv12K21P1sJlmAdkT19sQa7tTxVgiVYR2RPX2z33drhZyvBEqwjsqcvtvVubZIqwRKsI7KnL7btbi1+TXD1Fd/89Fdvz1LZ0xfbcrcWvTV0pyu++emv3qClsqcvtt1uLXxr6E5XfPPTX71FS2VPX2yr3druNlCwBOu47OmL7bJbq35DcNcrvvnpr16npbKnL7bDbq38DcFdr/jmp796pZbKnr7Y6t26bH0bKFiCdVz29MVW7tb+t4GCJVjHZU9fbNVu9bgNFCzBOi5k+su76UddsFvtUiVYgnVEyPSX30096uTdavPEavUV3/z0py7VboKmnx+tibv1+yOr5dduhyu++elPW6odhU0/N1qTduvS8zZQsATruMDpL9+pPFL9bv3wyGqDa7fDFd/89KuXamuh08+JVuFu3X4hcINrt8MV3/z0y5aqg+zpi+8Rx+/W5Zv3gt64A9zg2u1wxTc//dFL1Ur29L8p+2lr7G59mKrm1/0Gp9Dn9EcuVTvZ0//pcmXAPzpit96ubHPhtL3im5/+60vVWPb014aG67XdOpSq5tf9BqfQ5/RfWar2sqf/wIhyHd6tyxOROsV1v8Ep9Dn944t4ItnTP+A6XI8H7LHduo7U02+p2uDa7XDFNz/9V1e6tezpj7hfrrve3n7/30e1unDaXvHNT3/1lbBU9vSveChYEyJ1iut+g1Poc/qrN3+p7OmLFSRppwvHia85/dVrvVT29MVOfuE48TWnv3qtl8qevtjJLxwnvub0V6/1UtnTFzv5hePE15z+6rVeKnv6Yie/cJz4mtNfvdZLZU9f7OQXjhNfc/qr13qp7OmLnfzCceJrTn/1Wi+VPX2xk184TnzN6a9e66Wypy928gvHia85/dVrvVT29MVOfuE48TWnv3qtl8qevtjJLxwnvub0V6/1UtnTFzv5hePE15z+6rVeKmX6xz8wAbhh9SX8m01Oo9jq7zW0t/oi/s0mp1Fs0Vc767DGddh6m5xGsaxvsiv4zMfNOuyVTU6jWNY32RV85uNmHfbKJqdRLOub7Ao+83GzDntlk9MolvVNdgWf+bhZh72yyWkUy/omu4LPfNysw17Z5DSKZX2TXcFnPm7WYa9schrFsr7JruAzHzfrsFc2OY1iWd9kV/CZj5t12CubnEaxrG+yK/jMx8067JVNTqNY1jfZFXzm42Yd9somp1Es65vsCj7zcbMOe2WT0yiW9U12BZ/5uFmHvbLJaRTL+ia7gs983KzDXtnkNIplfZNdwWc+btZhr2xyGsWyvsmu4DMfN+uwVzY5jWJZ32RX8JmPm3XYK5ucRrGsb7Ir+MzHzTrslU1Oo1jWN9kVfObjZh32yianUSzrm+wKPvNxsw57ZZPTKJb1TXYFn/m4WYe9sslpFMv6JruCz3zcrMNe2eQ0imV9k13BZz5u1mGvbHIaxbK+ya7gMx8367BXNjkNgM8JFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgW0IZgARqBES0AAAIYSURBVG0IFtCGYAFtCBbQhmABbQgW0IZgAW0IFtCGYAFtCBbQhmABbQgWC/397e3tL7//n5+//J+3f688HbYnWKz005dG/e+N/4abBIuV/u/9p6r//Pqf3/20BTcJFkv97Uum/vr+H24IeYBgsdR/3zv1z1/cEPIQwWKtf74X679uCHmIYLHYX7+k6m9uCHmIYLHYv9++ckPIpwSL1f73t165IeRzgsVyf3lzQ8hjBIvl/vXeq3+sPgs6ECyW+7s7Qh4kWKz2s2fuPEqwWO2nr68SeojFpwSLxd5vCH/64zd04COCxVq/3hD+/N+fvv6GDnxEsFjr6+8Qfv0NHfiQYLHUH79D+Otv6Kw+G3YnWKz05+8Q/uf9v/5v9fmwOcFipW8+VOb9N3R+Wns2bE+wWOi7D5V5/w2dv688G/YnWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW0IVhAG4IFtCFYQBuCBbQhWEAbggW08f+xnMbaH7J2cQAAAABJRU5ErkJggg=="
     },
     "metadata": {
      "image/png": {
       "height": 480,
       "width": 600
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f <- function(x){x^2 + 1}\n",
    "\n",
    "a <- 0\n",
    "b <- 5\n",
    "\n",
    "number_of_rectangles <- 100\n",
    "  \n",
    "xgrid <- seq(a, b, length.out=number_of_rectangles+1)\n",
    "xs <- xgrid[1:number_of_rectangles]\n",
    "ys <- f(xs)\n",
    "\n",
    "plot(f,from = -1, to = 6 , col=\"red\", axes=FALSE)\n",
    "axis(1, pos=0, labels=FALSE)\n",
    "axis(2, pos=0, labels=FALSE, las=2)\n",
    "abline(v=c(a,b))\n",
    "\n",
    "xlefts <- xgrid[1:number_of_rectangles]\n",
    "xrights <- xgrid[2:(number_of_rectangles+1)]\n",
    "rect(xlefts, 0, xrights, ys, col=\"lightblue\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2aba7633-8198-4f0c-a0b9-79b4a9a7b21d",
   "metadata": {},
   "source": [
    "Use exactly the same code but set `number_of_rectangles <- 100`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "3078e69c-b571-4fee-be99-ebcc0f1eaca9",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "46.04375"
      ],
      "text/latex": [
       "46.04375"
      ],
      "text/markdown": [
       "46.04375"
      ],
      "text/plain": [
       "[1] 46.04375"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f <- function(x){x^2 + 1}\n",
    "\n",
    "a <- 0\n",
    "b <- 5\n",
    "\n",
    "number_of_rectangles <- 100\n",
    "xgrid <- seq(a, b, length.out=number_of_rectangles+1)\n",
    "dx <- (b-a)/number_of_rectangles\n",
    "xs <- xgrid[1:number_of_rectangles]\n",
    "ys <- f(xs)\n",
    "sum_areas <- dx*sum(ys)\n",
    "\n",
    "sum_areas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7d30453-f1bb-429e-a608-fbd0bc1808be",
   "metadata": {},
   "source": [
    "This is much better, with the approximation error less than 1.\n",
    "\n",
    "And for 10000 rectangles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "4d39194e-8181-4a74-a787-bdfa1a9f1630",
   "metadata": {
    "tags": [
     "hide-input"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABLAAAAPACAMAAAASenCyAAAAM1BMVEUAAABNTU1oaGh8fHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD/AAD///89ODILAAAACXBIWXMAACToAAAk6AGCYwUcAAAgAElEQVR4nO3c3WIk13GsUcjyj2zZot7/aa0hKYocDGYawM6Oyoq1bo59o+6sjvqOYJt6+TvAEi/pLwDwKMEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BAtYQLGANwQLWECxgDcEC1hAsYA3BGvTi6XJa+ai6rx9Wvi0mlI+q+/ph5dtiQvmouq8fVr4tJpSPqvv6YeXbYkL5qLqvH1a+LSaUj6r7+mHl22JC+ai6rx9Wvi0mlI+q+/ph5dtiQvmouq8fVr4tJpSPqvv6YeXbYkL5qLqvH1a+LSaUj6r7+mHl22JC+ai6rx9Wvi0mlI+q+/ph5dtiQvmouq8fVr4tJpSPqvv6YeXbYkL5qLqvH1a+LSaUj6r7+mHl22JC+ai6rx9Wvi0mhEb100+Rj33FKzVIsDhOsJgiWByXGdVVeiVYkwSL4yKj+kmwGggWx4WCFfjQb/JKDRIsjkuM6jq9EqxJgsVxgVFd5w9CwRolWBwXCdbTP/JNXqlBgsVxzx/VlXolWJMEi+OePqor/UEoWKMEi+MCwXryB36XV2qQYHHcs0d1rV4J1iTB4rgnj+pafxAK1ijB4rinB+upH/dDXqlBgsVxzx3V1XolWJMEi+OeOqqr/UEoWKMEi+OeOarr9UqwJgkWxz03WM/7rAd5pQYJFsc9cVQX7JVgTRIsjnveqC74B6FgjRIsjnvaqC7ZK8GaJFgc98RgPemD3sUrNUiwOO5Zo7pmrwRrkmBx3JNGdc0/CAVrlGBx3HNGddVeCdYkweK4ZwXrGZ/yAV6pE17ekv5i3M5TRnXZXgnWEYLFszxjVJf9g1CwRgkWxz1hVBfulWBNEiyOmx/VlXslWJMEi+PGR3XpXgnWJMHiuCcEa/gDPsUrNUiwOG56VNfulWBNEiyOGx7Vtf8gFKxRgsVxs6O6eq8Ea5JgcdzoqC7fK8GaJFgcNxyswX/xI7xSgwSL4yZHdf1eCdYkweK4wVFd/w9CwRolWBw3N6oNvRKsSYLFcWOjWtErwZokWBw3GKyhf+GjvFKDBIvjpka1o1eCNUmwOG5oVDv+IBSsUYLFcTOj2tIrwZokWBw3Mqo1vRKsSYLFcUPBGvgXHeGVGiRYHDcxqj29EqxJgsVxA6Pa8wehYI0SLI47P6pNvRKsSYLFccdHtapXgjVJsDju9Kh29UqwJgkWx50P1tl/vWFeqUGCxXGHR7WsV4I1SbA47uyolv1BKFijBIvjjo5qXa8Ea5JgcdzJUe3rlWBNEiyOOxusc/9aT+KVGiRYHHdwVAt7JViTBIvjzo1q4R+EgjVKsDju2KhW9kqwJgkWx50a1c5eCdYkweK4c8E686/zZF6pQYLFcYdGtbRXgjVJsDjuzKiW/kEoWKMEi+OOjGptrwRrkmBx3KFgHfgXifBKDRIsjjsxqr29EqxJgsVxB0a19w9CwRolWBz3+VFt7pVgTRIsjjsRrBPfI8QrNUiwOO7To1rdK8GaJFgc99lRrf6DULBGCRbHfXJUy3slWJMEi+M+HaxD3yPEKzVIsDjuc6Pa3ivBmiRYHPepUW3/g1CwRgkWx31mVPt7JViTBIvjPjGqG/RKsCYJFsd9fFR36JVgTRIsjvvwqG7RK8GaJFgc94lgHf0eIV6pQYLFcR8d1T16JViTBIvjPjiqe/xBKFijBIvjPjaqu/RKsCYJFsd9aFS36ZVgTRIsjvvIqO7TK8GaJFgc94FR3ahXgjVJsDjuQ8Ea+B4hXqlBgsVx7x/VnXolWJMEi+PePao7/UEoWKMEi+PeO6p79UqwJgkWx71zVDfrlWAd8fKW9Bfjdt4drKHvEeKVOkGweJb3jepuvRKsSYLFce8a1d3+IBSsUYLFce8Z1f16JViTBIvj3jGqG/ZKsCYJFse9K1iD3yPEKzVIsDju8VHdsVeCNUmwOO7hUd3xD0LBGiVYHPfoqO7ZK8GaJFgc9+CobtorwZokWBz3cLCGv0eIV2qQYHHcY6O6a68Ea5JgcdxDo7rrH4SCNUqwOO6RUd23V4I1SbA47oFR3bhXgjVJsDjux6O6c68Ea5JgcdwPR3XrXgnWJMHiuB+N6t69EqxJgsVxPw7Wc75HiFdqkGBx3A9GdfNeCdYkweK474/q5n8QCtYoweK4747q9r0SrEmCxXHfG9X9eyVYkwSL474zqoJeCdYkweK4t0fV0CvBmiRYHPfmqCp6JViTBIvj3hpVR68Ea5JgcdzbwXru9wjxSg0SLI57Y1QlvRKsSYLFcd8eVckfhII1SrA47pujqumVYE0SLI771qh6eiVYkwSL474xqqJeCdYkweK416Nq6pVgTRIsjvtWsBLfI8QrNUiwOO7VqKp6JViTBIvjvh5V1R+EgjVKsDjuq1GV9UqwJgkWx/1xVG29EqxJgsVxXwcr9T1CvFKDBIvj/jCqul4J1iTB4rjfj6ruD0LBGiVYHPe7URX2SrAmCRbH/WtUjb0SrEmCxXG/jaqyV4I1SbA47p+j6uyVYE0SLI77dVSlvRKsSYLFcb+MqrVXgjVJsDju51HV9kqwJgkWx30ZVW+vBGuSYHHcP0ZV3CvBmiRYHPfy0twrwZokWBz3JVjp7xDklRokWBz3j2Clv0KSV2qQYHFc+ai6rx9Wvi0GvJSPqvv6YeXb4rwXwWJM+bY47kWw0l/gzsq3xWkvgtV9/bDybXHYi2AJ1qTybXHWi2AJ1qjybXHUi2D9XbBGlW+Lk14E64vu64eVb4uTBOtn3dcPK98WB70I1s+6rx9Wvi3OeRGsX3RfP6x8WxzzIli/6r7+lJe3pL8Yt2BUv+m+/hTBYpBR/Uv39cPKt8UZ/n/B3+m+flj5tjjCv23/ve7rh5VvixP8zxn+oPv6YeXb4gD/g9E/6r5+WPm2+Dz/m5yvdF8/rHxbfJr/1fPXuq8fVr4tPk2wvtZ9/bDybfFZ3+xV96i6rx9Wvi0+6du96h5V9/XDyrfF57zRq+5RdV8/rHxbfMpbveoeVff1w8q3xWe82avuUXVfP6x8W3zC273qHlX39cPKt8XHfadX3aPqvn5Y+bb4sO/1qntU3dcPK98WH/XdXnWPqvv6YeXb4oO+36vuUXVfP6x8W3zMD3rVParu64eVb4sP+VGvukfVff2w8m3xET/sVfeouq8fVr4tPuDHveoeVff1w8q3xfs90KvuUXVfP6x8W7yfYP1A9/XDyrfFuz3Sq+5RdV8/rHxbvNdDveoeVff1w8q3xTs91qvuUXVfP6x8W7zPg73qHlX39cPKt8W7PNqr7lF1Xz+sfFu8x8O96h5V9/XDyrfFOzzeq+5RdV8/rHxbPO4dveoeVff1w8q3xcPe06vuUXVfP6x8WzzqXb3qHlX39cPKt8WD3ter7lF1Xz+sfFs85p296h5V9/XDyrfFQ97bq+5RdV8/rHxbPOLdveoeVff1w8q3xQPe36vuUXVfP6x8W/zYB3rVParu64eVb4sf+kivukfVff2w8m3xIx/qVfeouq8fVr4tfuBjveoeVff1w8q3xfd9sFfdo+q+flj5tviuj/aqe1Td1w8r3xbf8+FedY+q+/ph5dviOz7eq+5RdV8/rHxbvO0TveoeVff1w8q3xZs+06vuUXVfP6x8W7zlU73qHlX39cPKt8UbPter7lF1Xz+sfFt82yd71T2q7uuHlW+Lb/psr7pH1X39sPJt8S2f7lX3qLqvH1a+Lb7h873qHlX39cPKt8VrB3rVParu64eVb4tXTvSqe1Td1w8r3xZfO9Kr7lF1Xz+sfFt85UyvukfVff2w8m3xR4d61T2q7uuHlW+LPzjVq+5RdV8/rHxb/N6xXnWPqvv6YeXb4nfO9ap7VN3XDyvfFv9ysFfdo+q+flj5tvjNyV51j6r7+lNsi+852qvuUXVff4pt8R1ne9U9qu7rh5Vvi18c7lX3qLqvH1a+LX52ulfdo+q+flj5tvjieK+6R9V9/bDybfH3iV51j+rr6//nr3+LfI9bKt8WI73qHtXX1//55UWxTinfFiO96h7V19f/83m0P5cjPMNyI73qHpVgDfIMu830qntU3wjW//36/wa+zc14htWGetU9qq+v/7eXlz9/KVb7cznCM2w21avuUX19/X96Uud4YsXGetU9qq+v/19P6hxPrNdcr7pH9er6v3hSx3hitQZ71T2q19f/9d89qUM8sVKTuSof1VvXtz+XIzzDTrO96h6VYA3yDCsN96p7VII1yDNsNN2r7lF1Xz+sfFudxnvVParu64eVb6vSfK+6R9V9/bDybTV6Qq+6R9V9/bDybRV6Rq+6R9V9/bDybfV5Sq+6R9V9/bDybdV5Tq+6R9V9/bDybbV5Uq+6R9V9/bDybZV5Vq+6R9V9/bDybXV5Wq+6R9V9/bDybVV5Xq+6R9V9/bDybRV5Yq7KR9V9/bDybfV4aq+6R9V9/bDybdV4bq+6R9V9/bDybbV4cq+6R9V9/bDybZV4dq+6R9V9/bDybXV4eq+6R9V9/bDybVV4fq+6R9V9/bDybTUI9Kp7VN3XDyvfVoFEr7pH1X39sPJt3V+kV92j6r5+WPm27i6Tq/JRdV8/rHxbN5fqVfeouq8fVr6te4v1qntU3dcPK9/WreV61T2q7uuHlW/rzoK96h5V9/XDyrd1Y8ledY+q+/ph5du6r2ivukfVff2w8m3dVrZX3aPqvn5Y+bZuKpyr8lF1Xz+sfFv3lM5V+ai6rx9Wvq1bStfqpXxU3dcPK9/WHaVj9UX6GUR1Xz+sfFs3lG7Vz9IPIar7+mHl27qfdKp+kX4KUd3XDyvf1u2kS/Wr9GOI6r5+WPm2bibdqd+kH0RU9/XDyrd1L+lM/Uv6SUR1Xz+sfFu3kq7U76QfRVT39cPKt3Un6Uj9XvpZRHVfP6x8WzeSbtQfpB9GVPf1w8q3dR/pRP1R+mlEdV8/rHxbd5EO1NfSzyOq+/ph5du6iXSfXkk/kKju64eVb+se0nl6Lf1EorqvH1a+rVtI1+kb0o8kqvv6YeXbuoN0nL4l/Uyiuq8fVr6tG0i36ZvSDyWq+/ph5dtaL12mN6QfS1T39afY1h09sUHvkn4uUd3Xn2JbN/TEBL1P+sFEdV8/rHxbu6Wz9Lb0k4nqvn5Y+bZWS1fpO9KPJqr7+mHl21os3aTvSj+cqO7rh5Vva690kr4v/XSiuq8fVr6ttdJF+oH044nqvn5Y+ba2SgfpR9LPJ6r7+mHl21oq3aMfSj+gqO7rh5Vva6V0jR6QfkRR3dcPK9/WRukYPSL9jKK6rx9Wvq2F0i16SPohRXVfP6x8W/ukU/SY9FOK6r5+WPm2tkmH6FHp5xTVff2w8m0tk+7Qw9IPKqr7+mHl29olnaHHpZ9UVPf1w8q3tUq6Qu+QflRR3dcPK9/WIukGvUv6YUV1Xz+sfFt7pBP0PumnFdV9/bDyba2RLtA7pR9XVPf1w8q3tcRP6QC9V/qBRXVfP6x8Wzuk8/N+6ScW1X39sPJtrZCuzwekH1lU9/XDyre1QTo+H5F+ZlHd1w8r39b1pdPzMemnFtV9/bDybV1eujwflH5sUd3XDyvf1tWlw/NR6ecW1X39sPJtXdu6/2uG36SfXFT39cPKt3Vp6ep8QvrRRXVfP6x8W1eWjs5npJ9dVPf1w8q3dV17/xz8Iv30orqvH1a+rctKF+eT0o8vqvv6YeXbuqp0cD4r/fyiuq8fVr6ti0r35tPSDzCq+/ph5du6pHRtDkg/wqju64eVb+uK0rE5If0Mo7qvH1a+rQtKt+aI9EOM6r5+WPm2LiddmkPSjzGq+/ph5du6mnRoTkk/x6ju64eVb+ti0p05Jv0go7qvH1a+rUtJV+ag9KOM6r5+WPm2riQdmZPSzzKq+/ph5du6kHRjjko/zKju64eVb+sy0oU5LP04o7qvH1a+ratIB+a09POM6r5+WPm2LiLdl+PSDzSq+/ph5du6hHRdBqQfaVT39cPKt3UF6bhMSD/TqO7rh5Vv6wLSbRmRfqhR3dcPK99WXLosQ9KPNar7+mHl20pLh2VK+rlGdV8/rHxbWemszEk/2aju64eVbysqXZVB6Ucb1X39sPJtJaWjMin9bKO6rx9Wvq2cdFJmpZ9uVPf1w8q3FZMuyrD0443qvn5Y+bZC0j0Zl37AUd3XDyvfVkY6J/PSTziq+/ph5dtKSMfkGdLPOKr7+mHl2wpIt+Qp0g85qvv6YeXber50Sp4j/ZSjuq8fVr6tZ0uH5FnSzzmq+/ph5dt6snRHnib9oKO6rx9Wvq2nSlfkidKPOqr7+mHl23qmdESeKf2so7qvH1a+redJJ+S50k87qvv6YeXbepp0QZ4s/bijuq8/xbaCnpiKa0g/8Kju60+xrZwnluIi0k88qvv6YeXbeoZ0PBLSzzyq+/ph5dt6gnQ7ItIPPar7+mHl2xqXLkdI+rFHdV8/rHxb09LhSEk/96ju64eVb2tWOhs56Scf1X39sPJtjUpXIyj96KO6rx9Wvq1B6WZEpR9+VPf1w8q3NSedjKz004/qvn5Y+bampIORln7+Ud3XDyvf1ox0LvLSv0BU9/XDyrc14ad0LS4g/RtEdV8/rHxbA9KtuIT0jxDVff2w8m0d599e/Sz9M0R1Xz+sfFunpUNxFenfIar7+mHl2zornYnrSP8SUd3XDyvf1lHpSlxI+qeI6r5+WPm2Dko34lLSP0ZU9/XDyrd1TjoR15L+NaK6rx9Wvq1T0oG4mvTvEdV9/bDybZ2RzsP1pH+RqO7rh5Vv64h0HS4o/ZNEdV8/rHxbB6TbcEnpHyWq+/ph5dv6tHQZLir9s0R1Xz+sfFuflQ7DVaV/l6ju64eVb+tz0lm4rvQvE9V9/bDybX1GOgpXlv5torqvH1a+rU9IN+HS0j9OVPf1w8q39WHpIlxc+ueJ6r5+WPm2Pijdg8tL/0BR3dcPK9/Wh/gP6fuh9E8U1X39sPJtfYBcPSD9I0V1Xz+sfFvvJlcPSf9MUd3XDyvf1nulQ7BF+neK6r5+WPm23iedgT3Sv1RU9/XDyrf1Hv4afFz6t4rqvn5Y+bYeJ1fvkf61orqvH1a+rYelC7BM+ueK6r5+WPm2HpR+/9dJ/2BR3dcPK9/WY9Kv/z7pXyyq+/ph5dt6RPrl3yj9m0V1Xz+sfFs/ln71d0r/alHd1w8r39aPpF/8rdK/W1T39cPKt/V96dd+r/QvF9V9/bDybX1P+qXfLP3bRXVfP6x8W29Lv/K7pX+9qO7rh5Vv6y3pF3679O8X1X39sPJtfVv6dd8v/QtGdV8/rHxb35R+228g/RNGdV8/rHxb35B+128h/SNGdV8/rHxbr6Tf9JtI/4xR3dcPK9/WV9Lv+W2kf8io7uuHlW/rD9Jv+Y2kf8qo7uuHlW/rd9Lv+K2kf8yo7uuHlW/rN+k3/GbSP2dU9/XDyrf1q/T7fTvpHzSq+/ph5dv6Iv1y31H6N43qvn5Y+bbkakb6V43qvn5Y+7YYkf5do7qvH1a9rfRrfV/pXzaq+/phxdtKv9R3lv5to7qvH1a7rfQrfW/pXzeq+/phpdtKv9B3l/59o7qvH1a5rfTrfH/pXziq+/phfdtKv8sV0j9yVPf1p9jWz5740jZL/8xR3defYlt/l6unSf/QUd3XD2vaVvotLpL+qaO6rx9Ws630K9wl/WtHdV8/rGRb6Re4Tfr3juq+fljDtn5Kv7590j95VPf1w+6/rfS7Wyn9o0d1Xz/s7ttKv7ml0j97VPf1w269rfRr2yv9y0d1Xz/sxttKv7TN0r99VPf1w+66rfQbWy7980d1Xz/snttKv6/10gOI6r5+2A23lX5ZueGo3qP7+mG321b6VeWL9Aqiuq8fdq9tpd9TfpUeQlT39cPutK30W8pv0lOI6r5+2G22lX5F+b30GqK6rx92j22l30++kh5EVPf1w+6wrfTbySvpSUR1Xz9s/bbSrybfkl5FVPf1w3ZvK/1e8ob0MKK6rx+2eFvpl5K3pbcR1X39sK3bSr+RfFd6HlHd1w9bua3068iPpBcS1X39sH3bSr+LPCA9kqju64ft2pb/dPYl0kOJ6r5+2KJtqdUe6a1EdV8/bM220q8g75FeS1T39cN2bCv9/vFO6cFEdV8/7PrbSr97fEB6NFHd1w+7+LbSLx4fk95NVPf1wy68rfRLx8eltxPVff2wq24r/cbxKen5RHVfP+yK20q/bXxaekJR3dcPu9y20q8aJ6RXFNV9/bBLbSv9mnFKeklR3dcPu8q2/J+x30p6TlHd1w+7wrbE6nbSk4rqvn5YeltidUvZUYV1Xz8sua30W8WY3KguoPv6Yaltpd8oRmVGdREl1//0xdM/NbCt9MvEvKeP6kpKrv/pn576qU/eVvpF4jmeOqqrKbr++dF63rbS7xBP9KxRXVLZ9c+N1jO2lX57eLr5UV1Y4fU//cHkJw1vK/3ikDE6qqsrvf450RrbVvqVIWloVDt0Xz/8N+L5baXfFS7g9KhW6b7+F2P/buvgttIvCddxbFQbdV//Lz+9cuBf9MS20m8Hl/P5US3Wff1rR8P1iW2lXwqu6xOL3K/7+u84Ua73biv9JrDC+4d4I93XP+B1uB4P2APbSo+ffQ6seq/u69/j7XK96eXlJ/8JL5yWfhOiuq//jIeCBcellx/Vff2w9LK5pfSso7qvH5ZeNreUnnVU9/XD0svmltKzjuq+flh62dxSetZR3dcPSy+bW0rPOqr7+mHpZXNL6VlHdV8/LL1sbik966ju64ell80tpWcd1X39sPSyuaX0rKO6rx+WXja3lJ51VPf1w9LL5pbSs47qvn5YetncUnrWUS3Xp0cGy6Vf4V9c5GsMS//WsF76Jf7FRb7GsNDT7vpY5/rYeRf5GsO6fmRv8J0/t+tjX7nI1xjW9SN7g+/8uV0f+8pFvsawrh/ZG3znz+362Fcu8jWGdf3I3uA7f27Xx75yka8xrOtH9gbf+XO7PvaVi3yNYV0/sjf4zp/b9bGvXORrDOv6kb3Bd/7cro995SJfY1jXj+wNvvPndn3sKxf5GsO6fmRv8J0/t+tjX7nI1xjW9SN7g+/8uV0f+8pFvsawrh/ZG3znz+362Fcu8jWGdf3I3uA7f27Xx75yka8xrOtH9gbf+XO7PvaVi3yNYV0/sjf4zp/b9bGvXORrDOv6kb3Bd/7cro995SJfY1jXj+wNvvPndn3sKxf5GsO6fmRv8J0/t+tjX7nI1xjW9SN7g+/8uV0f+8pFvsawrh/ZG3znz+362Fcu8jWGdf3I3uA7f27Xx75yka8xrOtH9gbf+XO7PvaVi3yNYV0/sjf4zp/b9bGvXORrDOv6kb3Bd/7cro995SJfA+DHBAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEC1hDsIA1BAtYQ7CANQQLWEOwgDUEi6D/eHl5+bd//jd//cd/8/K/ya/D5QkWSX/6R6P+6xv/NXyTYJH031/+XdX//fxf/uHfbcE3CRZR//6PTP35y3/hD0IeIFhE/Wr9irMAAAG2SURBVO1Lp/7yd38Q8hDBIusvX4r1N38Q8hDBIuzP/0jVv/uDkIcIFmH/+/IrfxDyQ4JF2n/90it/EPJjgkXcv734g5DHCBZx//OlV/+Z/hZsIFjE/Ye/CHmQYJH2V/8zdx4lWKT96df/LaH/IRY/JFiEffmD8E+//RM68D2CRdbPfxD+9W9/+vWf0IHvESyyfv1nCH/9J3TguwSLqN/+GcKf/wmd9Lfh6gSLpH/9M4T/9+W/+u/09+HiBIuk3/2Hynz5J3T+lP02XJ5gEfSH/1CZL/+Ezn8kvw3XJ1jAGoIFrCFYwBqCBawhWMAaggWsIVjAGoIFrCFYwBqCBawhWMAaggWsIVjAGoIFrCFYwBqCBawhWMAaggWsIVjAGoIFrCFYwBqCBawhWMAaggWsIVjAGoIFrCFYwBqCBawhWMAaggWsIVjAGoIFrCFYwBqCBawhWMAaggWsIVjAGoIFrCFYwBqCBawhWMAaggWsIVjAGoIFrCFYwBqCBawhWMAaggWsIVjAGoIFrCFYwBqCBawhWMAaggWsIVjAGoIFrCFYwBqCBawhWMAaggWsIVjAGoIFrPH/zMIYwRMl4WUAAAAASUVORK5CYII="
     },
     "metadata": {
      "image/png": {
       "height": 480,
       "width": 600
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f <- function(x){x^2 + 1}\n",
    "\n",
    "a <- 0\n",
    "b <- 5\n",
    "\n",
    "number_of_rectangles <- 10000\n",
    "  \n",
    "xgrid <- seq(a, b, length.out=number_of_rectangles+1)\n",
    "xs <- xgrid[1:number_of_rectangles]\n",
    "ys <- f(xs)\n",
    "\n",
    "plot(f,from = -1, to = 6 , col=\"red\", axes=FALSE)\n",
    "axis(1, pos=0, labels=FALSE)\n",
    "axis(2, pos=0, labels=FALSE, las=2)\n",
    "abline(v=c(a,b))\n",
    "\n",
    "xlefts <- xgrid[1:number_of_rectangles]\n",
    "xrights <- xgrid[2:(number_of_rectangles+1)]\n",
    "rect(xlefts, 0, xrights, ys, col=\"lightblue\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fb5c0c3-3fe6-45bf-ad3f-f52de6f92477",
   "metadata": {},
   "source": [
    "Set `number_of_rectangles <- 10000`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "ce17212c-2a14-4af6-9d26-aadb1a2a176b",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "46.660416875"
      ],
      "text/latex": [
       "46.660416875"
      ],
      "text/markdown": [
       "46.660416875"
      ],
      "text/plain": [
       "[1] 46.66042"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f <- function(x){x^2 + 1}\n",
    "\n",
    "a <- 0\n",
    "b <- 5\n",
    "\n",
    "number_of_rectangles <- 10000\n",
    "xgrid <- seq(a, b, length.out=number_of_rectangles+1)\n",
    "dx <- (b-a)/number_of_rectangles\n",
    "xs <- xgrid[1:number_of_rectangles]\n",
    "ys <- f(xs)\n",
    "sum_areas <- dx*sum(ys)\n",
    "\n",
    "sum_areas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02e6785f-3c00-438c-a42a-f95583783598",
   "metadata": {},
   "source": [
    "And the numerical error is now only about 0.006."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da29042d-a57e-4cca-a255-95974c504bd4",
   "metadata": {},
   "source": [
    "This is the idea of numerical integration. There are different methods used to improve the accuracy beyond using more rectangles (e.g. using trapezoids rather than rectangles) but in all cases numerical integration works by approximating area under a curve by sums of areas of rectangles (or trapezoids, etc) and uses further methods to control the approximation error."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "844bfa41-de0c-41e4-a638-1cdf2bca2250",
   "metadata": {
    "tags": []
   },
   "source": [
    "(section:integrals:integrate)=\n",
    "\n",
    "## Numerical integration in R using `integrate()`\n",
    "\n",
    "I hope that you found the discussion in {ref}`section:integrals:byhand` interesting and useful and now understand what numerical integration means and how it works, as well as you have an intuition about why it works and why is an approximate method. \n",
    "\n",
    "However, in reality we rarely implement numerical integration from first principles and can instead use in-built functions that implement numerical integration efficiently and with high accuracy.\n",
    "\n",
    "The bottom line is that such a function in R is the function `integrate()` with usage\n",
    "\n",
    "`integrate(f,a,b)` \n",
    "\n",
    "where\n",
    "- `f` is the function to be integrated\n",
    "- `a` is the lower limit of integration \n",
    "- `b` is the upper limit of integration\n",
    "\n",
    "To use `integrate()` to compute the definite integral from before just set `f`, `a`, and `b`, as before"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "74c5ee08-baa1-43aa-bcea-2ef51343e61c",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "f <- function(x){x^2 + 1}\n",
    "\n",
    "a <- 0\n",
    "b <- 5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66685a42-4564-447b-8ad5-64dd16e9fcec",
   "metadata": {},
   "source": [
    "and call the function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "2ce6281b-ac85-48bf-bc08-210ebaf01099",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "46.66667 with absolute error < 5.2e-13"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "integrate(f,a,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61529b66-7df3-464b-b22d-87fcf7c98020",
   "metadata": {},
   "source": [
    "As you can see it works very quickly, and produces a result with approximation error of the order 0.0000000000005, which is quite precise.\n",
    "\n",
    "Using the `integrate()` function you are able to quickly evaluate definite integrals for any univariate function that you can define through an expression."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46cb8d04-099e-465a-b573-dfeabc6bf2e2",
   "metadata": {},
   "source": [
    "This completes the discussion of {ref}`chapter:integrals` as well as this set of notes on numerical methods. I hope that you have enjoyed it.\n",
    "\n",
    "Thank you for reading and all the best!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}