08. Численное интегрирование.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from scipy import integrate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Численное интегрирование\n",
    "\n",
    "Задача численного интегрирования\n",
    "заключается в нахождении некоторого приближения к значению определенного интеграла\n",
    "$$\n",
    "  I = \\int_a^b f(x) dx,\n",
    "$$\n",
    "где $[a, b]$ - отрезок интегрирования, а $f(x)$ - заданная функция.\n",
    "В основе методов численного интегрирования лежит\n",
    "суммирование (с некоторыми весами) значений функции в узлах $x_0, x_1, \\dots, x_n$, \n",
    "выбираемых на отрезке интегрирования:\n",
    "$$\n",
    "  \\int_a^b f(x) dx = \\sum_{k=0}^n A_k f(x_k) + R,\n",
    "$$\n",
    "где $R$ - остаточный член (ошибка). Обозначим\n",
    "$$\n",
    "  h = \\max \\left\\{x_k - x_{k-1}: ~ k = 1,2,\\dots,n\\right\\}.\n",
    "$$\n",
    "\n",
    "## Формула прямоугольников\n",
    "\n",
    "Пусть \n",
    "$$\n",
    "   a = x_0 < x_1 < \\dots < x_{n-1} < x_n = b,  \\tag{*}\n",
    "$$\n",
    "тогда\n",
    "$$\n",
    "  \\int_a^b f(x) dx \\approx \\sum_{k=0}^{n-1} f(x_k) \\cdot (x_k - x_{k-1}).\n",
    "$$\n",
    "<img src='http://aco.ifmo.ru/el_books/numerical_methods/lectures/images/glava2_clip_image_p002.png'>\n",
    "Эта формула называется *формулой левых прямоугольников*.\n",
    "В случае равноотстоящих узлов ($x_k - x_{k-1} = h$) получаем:\n",
    "$$\n",
    "  \\int_a^b f(x) dx \\approx h\\sum_{k=0}^{n-1} f(x_k).\n",
    "$$\n",
    "Она соответствует замене площади криволинейной трапеции\n",
    "площадью ступенчатой фигуры.\n",
    "Аналогична *формула правых прямоугольников*:\n",
    "$$\n",
    "  \\int_a^b f(x) dx \\approx \\sum_{k=1}^n f(x_k) \\cdot (x_k - x_{k-1}).\n",
    "$$\n",
    "<img src='http://aco.ifmo.ru/el_books/numerical_methods/lectures/images/glava2_clip_image_p003.png'>\n",
    "\n",
    "Для равноотстоящих узлов имеем\n",
    "$$\n",
    "  \\int_a^b f(x) dx \\approx h\\sum_{k=1}^n f(x_k).\n",
    "$$\n",
    "\n",
    "Если $f(x)$ - непрерывно дифференцируемая на отрезке $[a, b]$\n",
    "функция, то для остаточного члена формул\n",
    "правых и левых прямоугольников справедливо:\n",
    "$$\n",
    "  R = \\frac{b-a}{4} f(\\xi)\\, h, \\quad \\mbox{где $\\xi \\in [a, b]$}.\n",
    "$$\n",
    "Из формулы для $R$ следует, что если $h\\to 0$, то $R\\to 0$ (для непрерывно дифференцируемой функции),\n",
    "иными словами, чем мельче шаг, тем меньше ошибка.\n",
    "\n",
    "Небольшая модификация позволяет улучшить точность метода.\n",
    "В качестве узлов возьмем точки\n",
    "$$\n",
    "  x_{k + \\frac{1}{2}} = \\frac{x_k - x_{k-1}}{2} \\quad (k=0,1,\\dots,n-1),\n",
    "$$\n",
    "тогда\n",
    "$$\n",
    "  \\int_a^b f(x) dx \\approx \\sum_{k=0}^{n-1} f\\left(x_{k+\\frac{1}{2}}\\right) \\cdot (x_k - x_{k-1}).\n",
    "$$\n",
    "Эта формула называется *формула (средних) прямоугольников*.\n",
    "<img src='http://aco.ifmo.ru/el_books/numerical_methods/lectures/images/glava2_clip_image_p001.png'>\n",
    "В случае равноотстоящих узлов получаем:\n",
    "$$\n",
    "  \\int_a^b f(x) dx \\approx h\\sum_{k=0}^{n-1} f\\left(x_{k+\\frac{1}{2}}\\right).\n",
    "$$\n",
    "\n",
    "Если $f(x)$ - дважды непрерывно дифференцируемая на отрезке $[a, b]$\n",
    "функция, то для остаточного члена формулы прямоугольников справедливо:\n",
    "$$\n",
    "  R = \\frac{b-a}{24} f''(\\xi)\\, h^2, \\quad \\mbox{где $\\xi \\in [a, b]$}.\n",
    "$$\n",
    "Так же, как и для формулы остаточного члена в левых и правых прямоугольниках, \n",
    "получаем, что если $h\\to 0$, то $R\\to 0$ (для дважды непрерывно дифференцируемой функции).\n",
    "Обратите внимание, что $h$ теперь стоит во 2-й степени, а в формулах для остаточного члена правых и левых прямоуголниках был в 1-й степени. Таким образом, *теперь $R$ стремится к нулю быстрее (2-й порядок сходимости, а раньше был $1$-й) и, следовательно, с помощью средних прямоугольников можно получить большую точность при том же количестве узлов*. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вычисления по формуле прямоугольников с равноотстоящими узлами \n",
    "легко реализовать с помощью функции `np.sum()`.\n",
    "Рассмотрим, например, интеграл\n",
    "$$\n",
    "I = \\int_a^b \\frac{x}{\\sin x} dx  \\tag{**}\n",
    "$$\n",
    "\n",
    "График подынтегральной функции:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de3hc9X3n8fd3dLUulmXLlm0sbAzm6oCDvA65OXISCBAW0jZp4WkISdM6zabPNm3yLLTblm7SdLtNmzxpaEJIMCHZFKcbQkMMgbjEKtBwswk2NsbY2AbLli3rrtF9NN/9Y46ILEb2eEaay/Hn9TzzzLn8zpzf1wOfOfrNmXPM3RERkfCK5LoDIiIysxT0IiIhp6AXEQk5Bb2ISMgp6EVEQq441x1Ipq6uzpctW5bWtv39/VRWVk5vh3IgLHWAaslXYaklLHVAZrVs27at3d3nJ1uXl0G/bNkytm7dmta2zc3NNDU1TW+HciAsdYBqyVdhqSUsdUBmtZjZa1Ot09CNiEjIKehFREJOQS8iEnIKehGRkFPQi4iEnIJeRCTkFPQiIiGnoBcRyQObXzrGwwdGZuS18/IHUyIiZ5qf7zrKY6/FZuS1dUQvIpIHugZGqSyxGXltBb2ISB7oGRyhqmRmXltBLyKSB3RELyISct0Do1Qp6EVEwsnd6R4Y0RG9iEhY9Y+MEYs7VaUKehGRUOrqT5w/X6kvY0VEwqlncBQgd2P0ZtZgZlvMbLeZ7TKzPw6WzzWzzWa2N3iunWL7W4I2e83slukuQESk0HUNjB/R527oJgZ8zt0vAq4APmNmFwO3AY+5+wrgsWD+BGY2F7gdeBuwBrh9qg8EEZEzVfdAcESfqzF6d2919+eD6T5gN3AWcANwb9DsXuBDSTb/ALDZ3TvdvQvYDFw9HR0XEQmL7oE8GqM3s2XAW4FngHp3b4XEhwGwIMkmZwGHJsy3BMtERCTwxhH9DA3dpHxRMzOrAu4HPuvuvWYpdShZI5/i9dcD6wHq6+tpbm5OtWsniEajaW+bT8JSB6iWfBWWWsJQx4uvDFNeBEMD/TNSS0pBb2YlJEL+B+7+42DxMTNb5O6tZrYIaEuyaQvQNGF+CdCcbB/ufhdwF8Dq1au9qakpWbNTam5uJt1t80lY6gDVkq/CUksY6niw7QXm9XRSVRWZkVpSOevGgLuB3e7+lYl9A8bPorkF+EmSzR8FrjKz2uBL2KuCZSIiEugeGGVOxQwN0JPaGP07gZuB95rZC8HjWuDvgCvNbC9wZTCPma02s+8AuHsn8EXgueDxhWCZiIgEugdGqK0onbHXP+XQjbs/SfKxdoD3JWm/Ffj9CfMbgA3pdlBEJOy6BkZZPGcWMDgjr69fxoqI5Fh7dJi6qrIZe30FvYhIDo3E4vQNxZhbOXNDNwp6EZEc6gwuaDavSkEvIhJK7dFhAOZVauhGRCSUdEQvIhJyHf3jR/QKehGRUOqIBkf0GroREQmnjv4RiiPG7FkpX3rstCnoRURyqCM6zLyqUlK8UGRaFPQiIjnU2T/C3BkctgEFvYhITrVHR6ibwTNuQEEvIpJTHf3DM/qrWFDQi4jkVGd0ZEbPuAEFvYhIzgyOjNE/MjajP5YCBb2ISM4c70v8WGp+tY7oRURCqa1vCIAFMxz0pzxD38w2ANcBbe6+Mlj2Q+CCoMkcoNvdVyXZ9iDQB4wBMXdfPU39FhEpeMd6E0f0C6rLZ3Q/qfwU67vAHcD3xhe4+++MT5vZPwI9J9l+nbu3p9tBEZGweuOIfnaOj+jd/XEzW5ZsXXDj8N8G3ju93RIRCb+2vmGKI8bcGbxfLIC5+6kbJYJ+0/jQzYTla4GvTDUkY2YHgC7AgW+5+10n2cd6YD1AfX1948aNG1Ms4UTRaJSqqqq0ts0nYakDVEu+CksthVzHt3cMs7tzjK80VQCZ1bJu3bptUw6Pu/spH8AyYGeS5d8EPneS7RYHzwuA7cDaVPbX2Njo6dqyZUva2+aTsNThrlryVVhqKeQ6Pvqdp/36rz/xxnwmtQBbfYpMTfusGzMrBn4T+OFUbdz9SPDcBjwArEl3fyIiYXO8b5j5M/xFLGR2euX7gZfdvSXZSjOrNLPq8WngKmBnBvsTEQmVtr7hGf8iFlIIejO7D3gKuMDMWszsk8GqG4H7JrVdbGYPB7P1wJNmth14FnjI3R+Zvq6LiBSukViczv6RGT+HHlI76+amKZZ/PMmyI8C1wfR+4LIM+yciEkrjNwWf6XPoQb+MFRHJiba+8aDPg6EbERGZfkd7Ej+Wqp+tI3oRkVBq7RkEYNEcBb2ISCgd7RmitDjCvBm+6Qgo6EVEcuJIzxCLaspn9Kbg4xT0IiI50No9yMIsjM+Dgl5EJCdae4ZYPGdWVvaloBcRybKxuHOsNzF0kw0KehGRLGuPDhOLu4JeRCSsWoNz6BfVaOhGRCSUWrsT59Av1BG9iEg4HQmO6PVlrIhISB3uGmRWSRG1FSVZ2Z+CXkQkyw51DbCkdlZWfiwFCnoRkaxr6RqkYW5F1vaXyo1HNphZm5ntnLDsr83ssJm9EDyunWLbq81sj5ntM7PbprPjIiKFyN1p6RygoTY74/OQ2hH9d4Grkyz/qruvCh4PT15pZkXAPwPXABcDN5nZxZl0VkSk0PUMjtI3HGNJbR4d0bv740BnGq+9Btjn7vvdfQTYCNyQxuuIiIRGS1fi1MqGudk7oj/lrQRP4o/M7GPAVuBz7t41af1ZwKEJ8y3A26Z6MTNbD6wHqK+vp7m5Oa1ORaPRtLfNJ2GpA1RLvgpLLYVWx3NHYwAc2/8Sze17Tlg3Y7W4+ykfwDJg54T5eqCIxF8EXwI2JNnmI8B3JszfDHw9lf01NjZ6urZs2ZL2tvkkLHW4q5Z8FZZaCq2Ob/3HPl966ybv7h9507pMagG2+hSZmtZZN+5+zN3H3D0OfJvEMM1kLUDDhPklwJF09iciEhYtXYNUlxdTk6Vz6CHN0yvNbNGE2d8AdiZp9hywwszOMbNS4EbgwXT2JyISFoc6B2jI4hexkMIYvZndBzQBdWbWAtwONJnZKsCBg8CngraLSQzXXOvuMTP7I+BREsM8G9x914xUISJSIA52DHDRouqs7vOUQe/uNyVZfPcUbY8A106Yfxh406mXIiJnothYnEOdA1yzcmFW96tfxoqIZMnh7kFicWdZXWVW96ugFxHJkgPt/QAsm6egFxEJpYPjQV+X3S9jFfQiIllysGOAytIi5leVZXW/CnoRkSw52NHPsrrKrF2eeJyCXkQkSw6292f9i1hQ0IuIZMVwbIxDXYMsV9CLiITTgfZ+xuLOivrs/lgKFPQiIlmx91gUgBULqrK+bwW9iEgW7D3WR8Rg+XwN3YiIhNLetijL5lVSVlyU9X0r6EVEsmBvW5TzcjBsAwp6EZEZNxKLc7C9nxX1CnoRkVA62NFPLO6cn4MzbkBBLyIy48bPuNHQjYhISL1yrA8zOHd+nga9mW0wszYz2zlh2ZfN7GUz22FmD5jZnCm2PWhmL5rZC2a2dTo7LiJSKPa1RTl7bgXlJdk/4wZSO6L/LnD1pGWbgZXufinwCvBnJ9l+nbuvcvfV6XVRRKSw7W7t5YIcjc9DCkHv7o8DnZOW/dzdY8Hs08CSGeibiEjB6x+OcaCjn0sW1+SsD+bup25ktgzY5O4rk6z7KfBDd/+/SdYdALpI3ET8W+5+10n2sR5YD1BfX9+4cePGFEs4UTQapaoqN+Ng0yksdYBqyVdhqSXf69jbNcaXnhnis5eXsWrByW/TnUkt69at2zblyIm7n/IBLAN2Jln+P4EHCD4wkqxfHDwvALYDa1PZX2Njo6dry5YtaW+bT8JSh7tqyVdhqSXf67jnyf2+9NZNfrRn8JRtM6kF2OpTZGraZ92Y2S3AdcDvBjtJ9iFyJHhuCz4Q1qS7PxGRQrTrSC91VaUsqM7uXaUmSivozexq4FbgencfmKJNpZlVj08DVwE7k7UVEQmrnUd6uWRxTdbvKjVRKqdX3gc8BVxgZi1m9kngDqAa2BycOnln0HaxmT0cbFoPPGlm24FngYfc/ZEZqUJEJA8Nx8bYe6yPSxbPzmk/Tv7NAODuNyVZfPcUbY8A1wbT+4HLMuqdiEgBe+VolFjcWXlW7s64Af0yVkRkxuw60gOQ8yN6Bb2IyAx58XAP1eXFNNRW5LQfCnoRkRny/OvdrGqYQySSuy9iQUEvIjIjosMx9hzt5fKza3PdFQW9iMhM2HGom7jD5UsV9CIiobTttS4AVjUkvbhvVinoRURmwPOvd7FiQRU1s0py3RUFvYjIdIvHnedf786L8XlQ0IuITLv97f30DI5y+dLcD9uAgl5EZNo9H4zP64heRCSknt7fwbzK0pzdI3YyBb2IyDRyd375agdXnDsv5z+UGqegFxGZRgfa+znaO8Tbl8/LdVfeoKAXEZlGT+3vAOAd5yroRURC6ZevdlA/u4xz6ipz3ZU3pBT0ZrbBzNrMbOeEZXPNbLOZ7Q2ek369bGa3BG32BrcfFBEJJXfn6Vc7eMe5dTm9o9RkqR7Rfxe4etKy24DH3H0F8FgwfwIzmwvcDryNxP1ib5/qA0FEpNDtOdZHR/8Ib8+jYRtIMejd/XGgc9LiG4B7g+l7gQ8l2fQDwGZ373T3LmAzb/7AEBEJhf/YcxyAd51Xl+OenMjcPbWGZsuATe6+Mpjvdvc5E9Z3uXvtpG0+D5S7+98E838JDLr7PyR5/fXAeoD6+vrGjRs3plVQNBqlqio/zl3NRFjqANWSr8JSSz7V8b+fGWQgBl9856y0ts+klnXr1m1z99XJ1p3ynrEZSjZIlfSTxd3vAu4CWL16tTc1NaW1w+bmZtLdNp+EpQ5QLfkqLLXkSx29Q6Ps+/lm1q9dTlPThWm9xkzVkslZN8fMbBFA8NyWpE0L0DBhfglwJIN9iojkpSf3tjMWd9574YJcd+VNMgn6B4Hxs2huAX6SpM2jwFVmVht8CXtVsExEJFR+8XIbNbNKeGseXH9+slRPr7wPeAq4wMxazOyTwN8BV5rZXuDKYB4zW21m3wFw907gi8BzweMLwTIRkdCIx53mPcdZe/58iovy7+dJKY3Ru/tNU6x6X5K2W4HfnzC/AdiQVu9ERArArw510x4d5n15OGwD+mWsiEjGfvZiK6VFEd57kYJeRCR03J2f7TzKu1fUMbs897cNTEZBLyKSge0tPRzuHuSatyzKdVempKAXEcnAz15spThiXHlRfa67MiUFvYhImtydh15s5R3n1VFTkZ/DNqCgFxFJ29bXumjpGuRDqxbnuisnpaAXEUnT/dtaqCgt4gOXLMx1V05KQS8ikoah0TEe2tHK1SsXUlk205cNy4yCXkQkDZtfOkbfcIwPX74k1105JQW9iEgafrSthUU15VyRRzcBn4qCXkTkNB3qHODxvcf5yOoGIpH8uWXgVBT0IiKn6QfPvE7EjJvWNJy6cR5Q0IuInIbh2Bj/b+sh3nfhAhbVpHcnqWxT0IuInIZHdh6lo3+Ej16xNNddSZmCXkQkRe7OhicPcE5dZd7dAPxk0g56M7vAzF6Y8Og1s89OatNkZj0T2vxV5l0WEcmNZw90sr2lh0++65yC+BJ2XNpn+bv7HmAVgJkVAYeBB5I0fcLdr0t3PyIi+eLbT+xnbmUpv1UA585PNF1DN+8DXnX316bp9URE8sq+tij/vruNm69YyqzSolx357SYu2f+ImYbgOfd/Y5Jy5uA+4EW4AjweXffNcVrrAfWA9TX1zdu3Lgxrb5Eo1GqqqrS2jafhKUOUC35Kiy1ZKuOb+0YYtvRMf7hPRXMLpuZYZtMalm3bt02d1+ddKW7Z/QASoF2oD7JutlAVTB9LbA3lddsbGz0dG3ZsiXtbfNJWOpwVy35Kiy1ZKOOvcd6/ZzbNvnfPvTSjO4nk1qArT5Fpk7H0M01JI7mjyX5EOl192gw/TBQYmaF81W1iAjwtcf2UV5SxPq1y3PdlbRMR9DfBNyXbIWZLTQzC6bXBPvrmIZ9iohkxZ6jfWzacYSPv2MZ86rKct2dtGR0bU0zqwCuBD41YdkfArj7ncCHgU+bWQwYBG4M/sQQESkIX938CpWlxfzBuwvzaB4yDHp3HwDmTVp254TpO4A7Jm8nIlIIntnfwSO7jvLZ96+gtrI0191Jm34ZKyKSxFjc+cKml1hcU86n1p6b6+5kREEvIpLEj7YdYteRXm695sKCO29+MgW9iMgkfUOjfPnRPTQureX6y/L7xt+pUNCLiEzy94/soaN/hNv/68UEJw4WNAW9iMgEzx3s5PtPv8bH37GMS5fMyXV3poWCXkQkMDQ6xq337+CsObP4/FUX5Lo70yaj0ytFRMLk67/Yy/7j/Xzv99ZQWRaeeNQRvYgIiXPmv9n8Kh9pXMLa8+fnujvTSkEvIme8noFR/uSHL3D23Apuv/6SXHdn2oXnbxMRkTS4O3/2wA7a+oa5/9PvoCpEQzbjdEQvIme0e/7zIA+/eJTPf+ACLmsIx1k2kynoReSM9dSrHXzp4d1cdXE96wv4omWnoqAXkTPS4e5BPvMvz7NsXgX/+NuXFdTNvk+Xgl5Ezjg9g6P83j3PMRqLc9fHVlNdXpLrLs2o8H3rICJyEsOxMT71/a3sb4/y3U+s4dz5hX/f3FNR0IvIGSMed/70X7fz9P5OvnbjKt553plxZ9OMh27M7KCZvWhmL5jZ1iTrzcz+ycz2mdkOM7s8032KiJwud+f2B3fx0I5W/vzaC7lh1Vm57lLWTNcR/Tp3b59i3TXAiuDxNuCbwbOISFa4O3/1k118/+nX+NTa5QV9W8B0ZOPL2BuA73nC08AcM1uUhf2KiBCPO3/5k51vhPxt11wYiksPnw7L9F7dZnYA6AIc+Ja73zVp/Sbg79z9yWD+MeBWd986qd16YD1AfX1948aNG9PqTzQapaqq8L9cCUsdoFryVVhqOVkdY3Hn3pdGeLwlxrXnlPCR80vyOuQzeU/WrVu3zd1XJ13p7hk9gMXB8wJgO7B20vqHgHdNmH8MaDzZazY2Nnq6tmzZkva2+SQsdbirlnwVllqmqqN/eNR/755nfemtm/zLj7zs8Xg8ux1LQybvCbDVp8jUjMfo3f1I8NxmZg8Aa4DHJzRpARomzC8BjmS6XxGRqXREh/nkvVvZ0dLNFz+0kpuvWJrrLuVURmP0ZlZpZtXj08BVwM5JzR4EPhacfXMF0OPurZnsV0RkKrtbe/nQN/6T3a29fPOjjWd8yEPmZ93UAw8EY17FwL+4+yNm9ocA7n4n8DBwLbAPGAA+keE+RUSS+un2I/yPH+1g9qxiNq6/greeXZvrLuWFjILe3fcDlyVZfueEaQc+k8l+REROZiQW58uPvsy3nzjA6qW1fOOjl7OgujzX3cob+mWsiBS01mic3/zmf7LzcC8fe/tS/uKDF1NarMt4TaSgF5GC5O5sfO4Qtz81SEVZjG/d3MgHLlmY627lJQW9iBScw92D/MUDL7Jlz3Eunhfhnk+tpX62hmqmoqAXkYIxFnfu/eVB/uHnewD4y+su5pzRgwr5U1DQi0hB+NXrXfz1T19i+6Fumi6Yz998aCVLaitobn4t113Lewp6EclrrT2D/P0je3jgV4eZX13G125cxfWXLc7rSxnkGwW9iOSlnsFR7n5iP3c9sZ+4w2fWncunm86jqkyxdbr0LyYieSU6HOOeJw/w7Sf20zsU44OXLuK2qy+kYW5FrrtWsBT0IpIXeodG+cHTr3PX46/SNTDK+y+q50+uXMEli2ty3bWCp6AXkZw63D3IPU8eYONzh4gOx3jP+fP50yvP57KGObnuWmgo6EUk69yd51/v5ntPHWTTjsQ1Dq+7dBF/8O7lrDxLR/DTTUEvIlnTOzTKv/3qMP/yzOu8fLSPqrJiPvGOZXziXedw1pxZue5eaCnoRWRGxePOswc7eeD5wzy4/QiDo2OsPGs2f/sbb+H6VYt1Fk0W6F9YRGbE7tZe/u2Fw/z0hSMc6RmiorSI6y9bzO9ecTaXLtH4ezYp6EVkWrg7e471sXnXMTbtaGXPsT6KIsbaFXXces2FXHlxPRWlipxcSPtf3cwagO8BC4E4cJe7f21SmybgJ8CBYNGP3f0L6e5TRPJLbCzO1te62PzSMX7+0lEOdQ4C0Li0li/ccAkffMsi5lWV5biXksnHawz4nLs/H9xOcJuZbXb3lya1e8Ldr8tgPyKSR9p6h3hyXztP7G2neU8bXQOjlBZFeOd58/j0e87j/RctYIEuMpZX0g764L6vrcF0n5ntBs4CJge9iBSwwZExnj3YyROvHOfJfe28fLQPgLmVpbzn/PlcdclC1p4/X1+q5jFL3OkvwxcxWwY8Dqx0994Jy5uA+4EW4AjweXffNcVrrAfWA9TX1zdu3Lgxrb5Eo1GqqqrS2jafhKUOUC35aqpaoiPO3u4xXumK80rnGAd744w5FEfg/NoIl8wrYmVdEQ3VESJ5cGGxM+E9ScW6deu2ufvqZOsyDnozqwL+A/iSu/940rrZQNzdo2Z2LfA1d19xqtdcvXq1b926Na3+NDc309TUlNa2+SQsdYBqyVfNzc28e+172H88yvaWHp5/vYutBzt55VgUgNKiCJcuqeG/nDOXK5bPY82yucwqLcpxr98sbO9JurWY2ZRBn9HfWmZWQuKI/QeTQx5g4tG9uz9sZt8wszp3b89kvyJy+tydQ52DbG/pZkdLN4/vHOQzv3iU/pExAKrKimlcWssNq85i9dJaLmuYQ3lJ/gW7nL5Mzrox4G5gt7t/ZYo2C4Fj7u5mtgaIAB3p7lNEUhMdjrHnaB8vH+1NPLcmpnuHYgCUFkdYUgkfblzCW5bM4bIlNSyfX0VRJPdDMTL9MjmifydwM/Cimb0QLPtz4GwAd78T+DDwaTOLAYPAjT4dXwqICAB9Q6McaO9n//F+9rVFeTkI95auwTfaVJUVc8HCaq67bDErF9dw6ZIazq+v5pdPPk5T08oc9l6yJZOzbp4ETvrx7+53AHekuw8RgdGxOC1dg+w/HuVAez+vHu9n//Eo+9v7Od43/Ea7ooixvK6SVQ1zuGnN2VxQX80FC6tZUjtLd2M6w+l8KJEcG4s7bX1DHOoc5FDnAIe6BmjpSky3dA3S2jNIfMLfwbUVJSyfX0XT+fM5Z34ly+uqOHd+JWfPq6CsWGPq8mYKepEZFI87nQMjHO0Z4mjPEK29QxzrGeJob2K+pWuAw92DjI6dOKJZP7uMhtoK1pwzl4baWTTMrWD5/CqW11VSW1mao2qkUCnoRdIwOhans3+E9ugw7dER2vuG6egfpq13+I0QP9o7xLHeoTeFeMRgQXU59TXlXLK4hg+sXEhDbQUNcytoqJ3F4jmzdLaLTCsFvQgwEovTMzhK98AI3YOjdPWPnBjk0WHao8O8fmyAwcd/TtfAaNLXKSuOsKimnPrZ5axeWkt9TTmLZpezsKachTWzWDi7nLqqUoqLIlmuUM5kCnoJjdhYnOhwjL6h8ccovUOxRHgPjNI9GDxPnh4YeeNc8mSqy4qpqy5jXmUpi6oiXHTOIuZVllFXXcb8qlLmVZVRV1VGXVUpVWXF+uJT8o6CXnJqdCzOwMgYAyMxBkbGGBwZo384xsBoYjo6FKN3aJS+oVgQ4r+e7g3CPBoE++Do1GENibNS5swqoaaihNqKUhbOLueChdXMmVVKbUUJcypKqKlITNfMKmFeVSLcJw6jJH65+JaZ/mcRmVYKenmTeNwZGYszPBpnKDbG8Gic4dgYQ8HzcCzO0Gji+Y3lo+PL4wyOnhjch1qH+MaepxIhPhJjcGTsjXCfPH59MpWlRVSXl1BdXkxVeTE1s0pYUjuL6rLixLKyxLpfPxLztRWl1FSUUK2jbTlDKeinSTzujLkzFnfi7sTinlg2/hhfF4dYPE7cnbE4b1o//nipY4zIK8cZHYszOubBc5zYWCKEk02PjsUZmWJ68uuMjDmjsV8H98TwHonFM/q3KC2OUFFaRGVpMbNKixgbcSoM5leXcXZpBZWlRVQE6ypKiqgoK6aitCh4JKZnBfPV5SVUlRVTVVasX22KpClUQf/Bf3qCjp4BKrY24ySu7ZF4BscTz8EB5MR18QnTjLebvH0wTbDuhDB3Z0Z+7/vcsyk3LYoYJUVGSSRCSXGEkiKjOBKhdMJ0SXGE0mB6VmmEkohRXlJEWUmEsuIiyoojifniCGUlEcqLi054LisuojxJ24nblBUXvSmQE8Mdb5/ufx0RSVGogv78+moOx/tZsGA2ZoYBZgTPiXkMDJuwHCKWmGfS8hPbJcJrfHnEoKjIKDKjOGJEIonpSCQxXxQxIpZ4fuMxYT4yab4oAkWRSPAa8OL27TRe/lZKiiKUFEUoLf51WE8O9JJIhIiOdkVkCqEK+q/+zqrg6PHyXHclYyOHili9bG6uuyEiIaCTeUVEQk5BLyIScgp6EZGQU9CLiIScgl5EJOQyCnozu9rM9pjZPjO7Lcn6MjP7YbD+GTNblsn+RETk9KUd9GZWBPwzcA1wMXCTmV08qdkngS53Pw/4KvB/0t2fiIikJ5Mj+jXAPnff7+4jwEbghkltbgDuDaZ/BLzPdLEREZGsyuQHU2cBhybMtwBvm6qNu8fMrAeYB7RPfjEzWw+sD2ajZrYnzX7VJXv9AhSWOkC15Kuw1BKWOiCzWpZOtSKToE92ZD75ii+ptEksdL8LuCuD/iR2aLbV3Vdn+jq5FpY6QLXkq7DUEpY6YOZqyWTopgVomDC/BDgyVRszKwZqgM4M9ikiIqcpk6B/DlhhZueYWSlwI/DgpDYPArcE0x8GfuE+I9d5FBGRKaQ9dBOMuf8R8ChQBGxw911m9gVgq7s/CNwNfN/M9pE4kr9xOjp9ChkP/+SJsNQBqiVfhaWWsNQBM1SL6QBbRCTc9MtYEZGQU9CLiIRcQQZ9mC69kEItHzez42b2QvD4/Vz081TMbIOZtZnZzinWm5n9U1DnDjPL27vDpFBLk5n1THhP/irbfUyFmTWY2RYz221mu8zsj5O0KYj3JcVaCuV9KTezZ81se1DL/0rSZnozzN0L6kHii4ePMfQAAAKjSURBVN9XgeVAKbAduHhSm/8G3BlM3wj8MNf9zqCWjwN35LqvKdSyFrgc2DnF+muBn5H4bcUVwDO57nMGtTQBm3LdzxTqWARcHkxXA68k+e+rIN6XFGsplPfFgKpgugR4BrhiUptpzbBCPKIP06UXUqmlILj745z8NxI3AN/zhKeBOWa2KDu9Oz0p1FIQ3L3V3Z8PpvuA3SR+rT5RQbwvKdZSEIJ/62gwWxI8Jp8VM60ZVohBn+zSC5Pf8BMuvQCMX3oh36RSC8BvBX9W/8jMGpKsLwSp1loo3h786f0zM7sk1505leBP/7eSOHqcqODel5PUAgXyvphZkZm9ALQBm919yvdlOjKsEIN+Wi+9kGOp9POnwDJ3vxT4d379KV9oCuU9ScXzwFJ3vwz4OvBvOe7PSZlZFXA/8Fl37528Oskmefu+nKKWgnlf3H3M3VeRuKLAGjNbOanJtL4vhRj0Ybr0wilrcfcOdx8OZr8NNGapb9MtlfetILh77/if3u7+MFBiZnU57lZSZlZCIhh/4O4/TtKkYN6XU9VSSO/LOHfvBpqBqyetmtYMK8SgD9OlF05Zy6Tx0utJjE0WogeBjwVneVwB9Lh7a647lQ4zWzg+Xmpma0j8f9SR2169WdDHu4Hd7v6VKZoVxPuSSi0F9L7MN7M5wfQs4P3Ay5OaTWuGZXL1ypzw/L30wmlLsZb/bmbXAzEStXw8Zx0+CTO7j8RZD3Vm1gLcTuJLJtz9TuBhEmd47AMGgE/kpqenlkItHwY+bWYxYBC4MU8PJN4J3Ay8GIwHA/w5cDYU3PuSSi2F8r4sAu61xM2bIsC/uvummcwwXQJBRCTkCnHoRkREToOCXkQk5BT0IiIhp6AXEQk5Bb2ISMgp6EVEQk5BLyIScv8fTES5TnYIvRcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xx = np.linspace(1e-3, 3, 500)\n",
    "yy = xx/np.sin(xx)\n",
    "plt.plot(xx, yy)\n",
    "plt.grid()\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Этот интеграл не берется в элементарных функциях.\n",
    "\n",
    "Вычислим его численно при $a=0$, $b=3$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.449457117179284"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = 0\n",
    "b = 3\n",
    "n = 100\n",
    "h = (b - a)/n\n",
    "x = np.arange(a + h/2, b, h)\n",
    "y = x/np.sin(x)\n",
    "I_rect = h*sum(y)\n",
    "I_rect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как оценить точность полученного значения?\n",
    "Давайте удвоим количество узлов и сравним полученные значения. Скорее всего, старшие цифры, в которых ответы совпадают, точные."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.453808551489745"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n2 = 200\n",
    "h = (b - a)/n2\n",
    "x = np.arange(a + h/2, b, h)\n",
    "y = x/np.sin(x)\n",
    "h*sum(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Таким образом, можем утверждать, что $I \\approx 8.45$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "С помощью функции `np.cumsum()` можно вычислять численно интеграл с переменным верхним пределом (первообразную):\n",
    "$$\n",
    "  I(x) = \\int_0^x \\frac{t}{\\sin t} dt\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAfBklEQVR4nO3deXxV5b3v8c+TEchARkKYkjCEeUxAEKvBoVJFrdZavXWqx9LaW7U99nTUY09tb9vT2/bYnvZ4vbU9Vq3YKmqLrTOgIgIJYZJJQgYSEjKHDGTY2c/5I4FazLATsrPW3vm+X6+8Xjvstff+/VjwzZNnPWstY61FRETcK8TpAkREpG8KahERl1NQi4i4nIJaRMTlFNQiIi4X5o83TUpKsunp6QN+XXNzM1FRUUNfkAPUizupF3dSL5CXl1dtrU3u8Ulr7ZB/ZWVl2cHYuHHjoF7nRurFndSLO6kXa4Fc20umaupDRMTlFNQiIi6noBYRcTkFtYiIyymoRURcTkEtIuJyCmoREZfzywkvIiIjzev7T3C0uonpfrh0tEbUIiJDYH1+KU+8V0yIMUP+3gpqEZEhkF9Sz6LJ8X55bwW1iMg5OnGylfKGVhZPjvPL+yuoRUTOUX5JPQCLpiioRURcKf9YHeGhhjmpsX55fwW1iMg52lVSz5zUWEaFh/rl/RXUIiLnwNPpZW9ZA4un+OdAIiioRUTOycGKRlraO1nkpwOJ4GNQG2O+aox53xizzxjztDFmlN8qEhEJIDuKagFYmpHgt8/oN6iNMROBe4Bsa+08IBS40W8ViYgEkB1FtUyMG83EuNF++wxfpz7CgNHGmDBgDHDcbxWJiAQIay3bC+tYmu6/+WkAY304L90Ycy/wA+AU8Kq19rM9bLMWWAuQkpKStW7dugEX09TURHR09IBf50bqxZ3UizsFai8VzV6++fYpbp8bQc7kcGDwvaxatSrPWpvd45O93Uzx9BcQD7wJJAPhwAvAzX29Rje3VS9upV7cKVB7eWZ7iU37xgb7wYmTZ/7MqZvbXgoUWmurrLUdwHrg/AH/uBARCTLbCmtJiIpgWrJ/fxvwJahLgOXGmDHGGANcAhzwa1UiIi5nrWVrQTXLpyZg/HDFvA/rN6ittduAZ4GdwN7u1zzq16pERFyuqKaF4w2trJiW5PfP8unGAdbaB4EH/VyLiEjA2HKkGoCV0xL9/lk6M1FEZBDeLagmdewoMpKi/P5ZCmoRkQHyei3vFtRw/rQkv89Pg4JaRGTA3j9+kvqWDlZO9/+0ByioRUQGbPPhSgA+NiN5WD5PQS0iMkCbDlUxb2IsyTGRw/J5CmoRkQFoaOlgZ0kdOZnjhu0zFdQiIgPw9pEqvBZyZg7PtAcoqEVEBmTToSpiR4X59UYBZ1NQi4j4yNPp5c2DleTMHEdY6PDFp4JaRMRHO4rqqG1uZ/W88cP6uQpqEREfvfJ+BZFhIVyUOXzz06CgFhHxibWWV96v4MLMZKIifbpM0pBRUIuI+GBPaQPlDa2snju80x6goBYR8cnL71cQGmK4ZPbwrZ8+TUEtItIPay0v76tgxdRE4sZEDPvnK6hFRPrxQWUThdXNXD7Mqz1OU1CLiPTjlX0VGAOXz0lx5PMV1CIi/XhpbzmLJ8cxLnaUI5+voBYR6cP+4yc5WNHIJxdPdKwGBbWISB/W7ywlPNSwZsEEx2pQUIuI9MLT6eXF3cdZNXMcCVHDv9rjNAW1iEgv3jlSTVVjG9ctmeRoHQpqEZFerN9ZxtjR4ayaNbzX9jibglpEpAeNrR288n4FVy1MJTIs1NFaFNQiIj34294K2jxex6c9QEEtItKjZ3eWkpEUxeJhvJNLbxTUIiJn+eBEI9sLa/l09iSMMU6Xo6AWETnbk+8VExEawmeyJztdCqCgFhH5B81tHp7bWcaVC1JJjI50uhxAQS0i8g9e2FVGU5uHm5enOV3KGQpqEZFu1lqe2FrMnNRYlkxx/iDiaQpqEZFuecV1HKxo5JYVaa44iHiaglpEpNsT7xUTExnGNYucuwBTTxTUIiJAecMpXtpTzvXZkxgTMbx3Ge+PglpEBHjs7UIscMfKDKdL+QgFtYiMeA0tHTy9vYSrFqQyOWGM0+V8hIJaREa8J7cV09zeydoLpzldSo8U1CIyorV2dPK7LYVclJnMnAmxTpfTI5+C2hgTZ4x51hhz0BhzwBizwt+FiYgMh+d2llLd1M4XL3LnaBrA10ObDwMvW2uvN8ZEAO6bxBERGSBPp5dH3zrKwkljWT41welyetXviNoYEwtcCDwGYK1tt9bW+7swERF/W59fRnFNC/971XRXneByNmOt7XsDYxYBjwL7gYVAHnCvtbb5rO3WAmsBUlJSstatWzfgYpqamoiOjh7w69xIvbiTenEnJ3rxeC3ffPsUMeGGf10xasiCerC9rFq1Ks9am93jk9baPr+AbMADnNf9/cPAQ329Jisryw7Gxo0bB/U6N1Iv7qRe3MmJXp58r8imfWODffPgiSF938H2AuTaXjLVl4OJpUCptXZb9/fPAksG/ONCRMQlWjs6+eUbR1gyJY6cTGdvXOuLfoPaWlsBHDPGzOz+o0vomgYREQlIT28voeJkK/d9fKar56ZP83XVx93AU90rPo4Cn/NfSSIi/tPc5uHXmwo4LyOB86clOl2OT3wKamvtLrrmqkVEAtr/e+soVY1tPHJzVkCMpkFnJorICFLecIpH3ypgzYJUstLinS7HZwpqERkxfvLKIbwWvrF6ltOlDIiCWkRGhL2lDazfWcYdKzNceYW8viioRSToWWv5/kv7SYyK4Eur3HtNj94oqEUk6P1lTznbCmv56mWZxI4Kd7qcAVNQi0hQO9nawUMb9rNg0lhuWjbF6XIGxV03BhMRGWI/e/Uw1U1tPHZbNqEhgbEc72waUYtI0NpX1sDvtxZxy/I0FkyKc7qcQVNQi0hQ6vRavvP8XhKiIrnv4zP7f4GLKahFJCj9bkshu0sbeGDNbMaODrwDiB+moBaRoFNY3cxPXjnEpbPHcfXCCU6Xc84U1CISVLxey9ef3U1kWAg/uHZ+wFzPoy8KahEJKo9vLWJHUR0PXjWXlNhRTpczJBTUIhI0Cqub+feXD7FqZjLXLZnodDlDRkEtIkGh3ePl3nX5RISF8MPrFgTFlMdpOuFFRILCz18/zJ7SBh65eQnjxwbHlMdpGlGLSMB7t6CaRzYXcOPSyayel+p0OUNOQS0iAa2+pZ37/rib9MQoHlgzx+ly/EJTHyISsLxey1ef2UV1UxvP3XU+UZHBGWkaUYtIwPqvzQVsPFTFA2vmBPS1PPqjoBaRgPTukWp++uohrl44gVuWpzldjl8pqEUk4Jw42co96/KZmhzND68LjrMP+xKcEzoiErRaOzpZ+/tcTrV38vTnlwTtvPSHBX+HIhI0rLV8a/1edpc28OgtWcxIiXG6pGGhqQ8RCRiPvnWU5/PLuO+yTD4+d7zT5QwbBbWIBIQ3DpzgRy8f5Mr5qXz54ulOlzOsFNQi4nq7jtXz5T/kM3/iWH7y6eC6jocvFNQi4mpF1c3c8d87SI6J5LHbljImYuQdWlNQi4hr1TS1cdvvtmOt5b8/t5TkmEinS3LEyPvRJCIBoaXdwx2P51LR0MofPr+cqcnRTpfkGI2oRcR1Ojq93P2HfPaW1vPLmxaTlRbvdEmO0ohaRFzF0+nlK8/s4o2DlTz0yXkjahlebzSiFhHX6Lox7R5e2lPOt6+YFfTX8PCVglpEXMFay3de2Mf6/DL++bJM1l44zemSXENBLSKOs9byb3/Zz9PbS/hSzjTuHmEntPRHc9Qi4ihrLesOtfNKURF3rMzgXy6fOeJOaOmPRtQi4hiv13L/C/t4pcjDbSvSeGDNbIV0D3wOamNMqDEm3xizwZ8FicjI4On08rVnd/PUthKuyAjnu1fPVUj3YiBTH/cCB4BYP9UiIiNEu8fLV5/ZxUt7y/nnyzKZH1KqkO6DTyNqY8wk4ErgN/4tR0SCXWtHJ3c9mcdLe8u5/8rZ3HPJDIV0P3wdUf8H8HVgZFylW0T8ora5nTsf30H+sXq+/8l53Kx10j4x1tq+NzBmDXCFtfZLxpgc4GvW2jU9bLcWWAuQkpKStW7dugEX09TURHR0cJzPr17cSb04p6rFy09zW6lutXxhQSRLx/99nBhovfRlsL2sWrUqz1qb3eOT1to+v4AfAqVAEVABtABP9vWarKwsOxgbN24c1OvcSL24k3pxxp5j9Tbrodfsgu++YncU1nzk+UDqpT+D7QXItb1kar9z1Nbab1lrJ1lr04EbgTettTcP+MeFiIxImw5V8plHtxIZFsJzd51PdnqC0yUFHJ3wIiJ+Ya3lt1uK+MFL+5mdGsvvbl/KuNhRTpcVkAYU1NbaTcAmv1QiIkGjzdPJAy/s44+5payeO56f3rCQqEiNCwdLf3MiMqSqm9q468k8dhTVcc/F0/nKpZmEhGj53blQUIvIkNlX1sAXnsijuqmNX960mKsWTnC6pKCgoBaRIfHMjhIeePF9EqMi+NMXV7BgUpzTJQUNBbWInJPWjq756D/llXLB9CQevnERidEj8ya0/qKgFpFBK65p5q4nd7K//CT3XDydey/NJFTz0UNOQS0ig/LirjLuf34fISGG392+lFWzxjldUtBSUIvIgDS1efjXF/exfmcZWWnx/MdnFjE5YYzTZQU1BbWI+GzXsXruXZfPsdoW7r1kBndfPJ2wUN1/xN8U1CLSr06v5ZHNBfz8tcOkxI7imS+sYKlOBR82CmoR6dOx2ha+9qfdbCusZc2CVH5w7XzGjg53uqwRRUEtIj3yei1PbivmR387SIgx/OT6BVyfNUkX+XeAglpEPqKkpoV/ebZrFH1hZjI/vG4+E+NGO13WiKWgFpEzvF7LE+91jaLDQgw//tR8bsierFG0wxTUIgJAQVUT31q/l+2FtVzUPYqeoFG0KyioRUa41o5Ofr3xCI9sPsqo8BD+/VML+HS25qLdREEtMoK9dbiKB17cR3FNC59cNIHvXDmH5Bhdp8NtFNQiI1DlyVa+t2E/G/aUMzUpiqfuPI+V05OcLkt6oaAWGUHaPV5+v7WIh1//gLZOL1+9NJMv5kwlMizU6dKkDwpqkRFi48FKHtqwn6PVzVyUmcx3r55LRlKU02WJDxTUIkHuSGUT339pP5sOVTE1KUpXugtACmqRINXQ0sEv3vyAx98tYnR4KPdfOZtbV6QTEaaLKAUaBbVIkGnt6OTxd4v41cYjNLZ5uHHpFO77eCZJuutKwFJQiwSJTq/l+fwyfvbqIY43tJIzM5lvrJ7F7NRYp0uTc6SgFglw1lo2Harixy8f5GBFIwsnjeX/3rCQ86dpuV2wUFCLBLB3C6r5+WuH2VFUR1riGP7zfy3myvmpOqswyCioRQLQtqM1/HDbKQ69vI2U2EgeumYun1k6RQcKg5SCWiSA5BbV8vPXD7PlSA1jIw0PXjWHm5ZNYVS4TlgJZgpqkQCQX1LHz1//gLcOV5EUHcH9V85mSnsxH1+Z4XRpMgwU1CIuZa3l3YIafr3pCFuO1JAQFcG3PjGLW1akMSYijE2bSpwuUYaJglrEZbxey6v7T/Bfm46wu7SB5JhIvvWJWdy8PI2oSP2XHYm010VcoqPTy4u7jvPI5gKOVDaRljiG/3PtfK5bMlFz0COcglrEYU1tHv6Ue4zfvF1IWf0pZo2P4Rc3LeaKeeMJC9UqDlFQizimrP4Uj79bxNPbS2hs9bA0PZ7vf3IeOTOTtQ5a/oGCWmSY5ZfU8dg7hfxtXwUAn5g3nn+6IIPFU+IdrkzcSkEtMgw8nV5e3X+Cx94pJK+4jphRYfzTBRncdn46E3UDWemHglrEjyobW3lm+zH+sL2E8oZWpiSM4cGr5vDp7MlEawWH+Ej/UkSGmLWWHUV1PPFeMS/vK6ej0/KxGUl89+q5XDo7hdAQzT/LwCioRYZIU5uH5/PLeHJrMYdONBI7KoxbV6Tz2fOmMDU52unyJID1G9TGmMnA74HxgBd41Fr7sL8LEwkUBytO8tR7JazfWUpzeydzJ8Ty40/N5+qFExkdofXPcu58GVF7gPustTuNMTFAnjHmNWvtfj/XJuJaja0d/GV3Oc/sKGF3aQMRYSGsWZDKLcvTWDQ5TsvrZEj1G9TW2nKgvPtxozHmADARUFDLiGKtJbe4jmd2HOOlPeWc6uhkZkoMD6yZw7WLJ5IQFeF0iRKkjLXW942NSQfeAuZZa0+e9dxaYC1ASkpK1rp16wZcTFNTE9HRwTGXp17caTC9NLRZthzv4K1SDxXNllGhcF5qGBdNCiNjbIhjo+eRvl/carC9rFq1Ks9am93Tcz4HtTEmGtgM/MBau76vbbOzs21ubu6AC920aRM5OTkDfp0bqRd38rWXNk8nGw9Wsn5nGW8erMTjtWSnxfOZpZO5ckEqYyKcPw4/EvdLIBhsL8aYXoPap39txphw4Dngqf5CWiRQWWvJK65jfX4ZL+0pp+FUB8kxkdxxQQY3ZE9m+rjgGPFJ4PFl1YcBHgMOWGt/5v+SRIZXUXUz6/PLeCG/jJLaFkaHh3L53BSuXTKJldMSdWEkcZwvI+qVwC3AXmPMru4/+7a19q/+K0vEv+pb2vnLnnKe31nKzpJ6jIGV05K495IZXD5vvM4aFFfxZdXHO4DWGknAa2rzsPW4hycf38Hmw1V0dFpmpsTwzU/M4ppFE0gdq2tuiDtp2CBB7VR7J28erGTDnuO8ebCSNo+X8bEnuW1FOtcumcic1FiteRbXU1BL0GnzdLL5UBUb9pTz+oETtLR3khQdyY1LJzPRe4I7r7mYEF1vQwKIglqCQrvHy5aCajbsLufV9ytobPMQNyacaxZN5KoFqZw3NZHQEMOmTdUKaQk4CmoJWC3tHt46XMXL+yp442Alja0eYkaFcfm88axZkMrK6UmEa8WGBAEFtQSUhlMdvHnwBC/vq2Dz4SpaO7zEjQln9dzxXD53PB/LTCIyTBdCkuCioBbXq25q47X9XeH8bkE1HZ2WlNhIbsiezOq541mWkaC1zhLUFNTiSsdqW7rC+f0Kcotq8VqYkjCGO1ZmcPm88SyaFKe5ZhkxFNTiCp1ey65j9bxx4ARvHKjk0IlGAGaNj+Hui2ewet54Zo2P0VI6GZEU1OKY5jYPb39QxesHKtl4sJKa5nZCQwxL0+O5/8rZXDI7hYykKKfLFHGcglqGVVn9Kd44cILXD1TyXkEN7Z1eYkeFkTNzHJfMHkdO5jjGjgl3ukwRV1FQi195vZbdpfW8caCS1w+c4GBF15RGRlIUt65I45LZKWSnx2sZnUgfFNQy5Kqb2nj7gyo2Hari7Q+qqe2e0shKi+fbV8ziktkpTNPNXkV8pqCWc+bp9LLrWD2bD3eF896yBgASoyLIyUzmopnJXJSZTNwY3apKZDAU1DIoJ062svlwFZsPVfH2B1WcbPUQYmDJlHjuuyyTnJnjmDshVkvoRIaAglp80u7xkldc1z1qrjwz1zwuJpLV88ZzUeY4LpiepAOBIn6goJYeWWspqGpmy5FqXtjZypc3vkZTm4ewEEN2ejzfWD2LnJnJWtssMgwU1HJGVWMbW45U886RarYcqaa8oRWA5NGGqxZO5qLMZFZOTyRmlEbNIsNJQT2CtbR72FZYy5YPusL59HTG2NHhrJyeyN3Tk7lgehJH924nJ2e+w9WKjFwK6hHE0+llb1kD73QH886SOjo6LRGhIWSnx/P11TO5YHoScyeMJfRDBwGPOliziCiog9rpeeatBV3B/G5BDY2tHgDmpMZyx8oMVk5PYml6AqMjdGlQEbdSUAcRay3FNS1sPVrD1oIath6toaqxDYCJcaO5Yl4qF8xI4vxpiSRGRzpcrYj4SkEd4I7VtpwJ5feO1vz9AGBMJCumJrJiWiLLpyaSnjhGqzNEApSCOsAcrz91Jpi3FtRQVn8K6DoLcPnURJZPS2TF1ESmJUcpmEWChILa5U6cbGVrQddoeevRGoprWgCIGxPO8oxE1l44lRXTEpkxLlrBLBKkFNQuU3mylW2FtV1TGQU1HK1uBiB2VBjnTU3k1hXprJiayKzxMTo9W2SEUFA7yFpLad0pthXWsqOwlu1FtRR2B3N0ZBjLMhK4adkUVkxLZHZq7D8smRORkUNBPYy6lss1sa2wlu3dX6cP/sWNCWdpegKfPW8KS9MTmDshVjdsFRFAQe1XnV7LgfKT3cFcw46iOmqb24Guixkty0jgvIwElmV0zTFrKkNEeqKgHkLtHi97y+rZVljLy7mt3L3xVRrbuk4wmZIwhlUzx3UHcwJpWi4nIj5SUJ+DlnYP+SX1Z0bM+SX1tHm8AEyINly9aDLLuoM5dexoh6sVkUCloB6AmqY2covryC2qJbe4jr2lDXi8lhADcybE8tnz0liWkcDS9Hj25m7VhYxEZEgoqHthreVodTN5RXXsKKolr7juzFK5iNAQ5k8ay+cvnMqyjASy0uKJ1aU/RcRPFNTd2jyd7Cs7SV5xLTuK6sgr/vuBv7gx4WSnxXPD0slkp8Uzb+JYRoXrIkYiMjxGbFA3tHSQV9IdykV17Cqtp717fjk9sevA39L0eLLT45mapBUZIuKcERHU1lqO1Z4i98xouZbDJ5oACAsxzJ04lluXp5GdHk9WWgLJMbqynIi4R1AGtafTy/7yk+QW1ZFbXEtuUR2V3Zf7jIkMY0laPFcvnEBWWgKLJsfpWswi4mpBEdSNrR3kl9SfWZGx61g9Le2dQNd1mFdMSyQ7PYHstHgyU2J0KraIBJSADOrj9af+vkyuqI6DFSfxWggxMDs1lhuyJ5OV1jW/rPXLIhLofApqY8xq4GEgFPiNtfZHfq3qQzq9lkMVjf+wGuP0NZjHRISyeEocd188g+z0eBZPiSc6MiB/9oiI9KrfVDPGhAK/Ai4DSoEdxpg/W2v3+6OglnYPu47Vd61fLq4jv7juzGnYKbGRZKclcOfHMshOS2B2aowuXCQiQc+X4ecy4Ii19iiAMWYdcA0wpEHd5unke1tPUfLqq3i8FmMgc1wMVy+aQHZ6PNlpCUyKH63rY4jIiGOstX1vYMz1wGpr7Z3d398CnGet/fJZ260F1gKkpKRkrVu3bsDF/CqviZSYCGbEhzA9LpSo8MAN5aamJqKjo50uY0ioF3dSL+402F5WrVqVZ63N7uk5X0bUPaXlR9LdWvso8ChAdna2zcnJGUiN3TYxuNe5z6ZN6sWN1Is7qZe++TLBWwpM/tD3k4DjQ1qFiIj0ypeg3gHMMMZkGGMigBuBP/u3LBEROa3fqQ9rrccY82XgFbqW5/3WWvu+3ysTERHAx3XU1tq/An/1cy0iItIDLUIWEXE5BbWIiMspqEVEXE5BLSLicv2emTioNzWmCigexEuTgOohLscp6sWd1Is7qRdIs9Ym9/SEX4J6sIwxub2dQhlo1Is7qRd3Ui9909SHiIjLKahFRFzObUH9qNMFDCH14k7qxZ3USx9cNUctIiIf5bYRtYiInEVBLSLico4EtTFmtTHmkDHmiDHmmz08H2mMeab7+W3GmPThr9I3PvRyuzGmyhizq/vrTifq7I8x5rfGmEpjzL5enjfGmF9097nHGLNkuGv0lQ+95BhjGj60T/51uGv0lTFmsjFmozHmgDHmfWPMvT1sExD7xsdeAmLfGGNGGWO2G2N2d/fybz1sM3Q5Zq0d1i+6LpVaAEwFIoDdwJyztvkS8Ej34xuBZ4a7ziHs5XbgP52u1YdeLgSWAPt6ef4K4G903fFnObDN6ZrPoZccYIPTdfrYSyqwpPtxDHC4h39jAbFvfOwlIPZN9991dPfjcGAbsPysbYYsx5wYUZ+5Wa61th04fbPcD7sGeLz78bPAJcadd7X1pZeAYK19C6jtY5NrgN/bLu8BccaY1OGpbmB86CVgWGvLrbU7ux83AgeAiWdtFhD7xsdeAkL333VT97fh3V9nr8wYshxzIqgnAsc+9H0pH91ZZ7ax1nqABiBxWKobGF96AfhU96+kzxpjJvfwfCDwtddAsaL719a/GWPmOl2ML7p/dV5M1+jtwwJu3/TRCwTIvjHGhBpjdgGVwGvW2l73y7nmmBNB7cvNcn26oa4L+FLnX4B0a+0C4HX+/hM20ATKPvHFTrquq7AQ+CXwgsP19MsYEw08B3zFWnvy7Kd7eIlr900/vQTMvrHWdlprF9F1H9llxph5Z20yZPvFiaD25Wa5Z7YxxoQBY3Hnr7L99mKtrbHWtnV/+/+BrGGqbagFzU2OrbUnT//aarvuXhRujElyuKxeGWPC6Qq2p6y163vYJGD2TX+9BNq+AbDW1gObgNVnPTVkOeZEUPtys9w/A7d1P74eeNN2z8i7TL+9nDVXeDVd83KB6M/Ard0rDJYDDdbacqeLGgxjzPjTc4XGmGV0/T+ocbaqnnXX+RhwwFr7s142C4h940svgbJvjDHJxpi47sejgUuBg2dtNmQ55tM9E4eS7eVmucaY7wG51to/07UznzDGHKHrJ9CNw12nL3zs5R5jzNWAh65ebnes4D4YY56m64h7kjGmFHiQrgMkWGsfoeuemVcAR4AW4HPOVNo/H3q5HrjLGOMBTgE3unQgALASuAXY2z0fCvBtYAoE3L7xpZdA2TepwOPGmFC6fpj80Vq7wV85plPIRURcTmcmioi4nIJaRMTlFNQiIi6noBYRcTkFtYiIyymoRURcTkEtIuJy/wMZvSQa/TB/6QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, h*np.cumsum(y))\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Формула трапеций\n",
    "\n",
    "Пусть точки $x_k$ ($k=0,1,\\dots, n$) удовлетворяют условиям ($*$).\n",
    "Тогда\n",
    "$$\n",
    "  \\int_a^b f(x) dx \\approx \\sum_{k=1}^n \\frac{f(x_k) + f(x_k)}{2} \\cdot (x_k - x_{k-1}).\n",
    "$$\n",
    "Эта формула называется *формулой трапеций*.\n",
    "Она соответствует замене площади криволинейной трапеции на сумму площадей прямоугольных \n",
    "трапеций.\n",
    "\n",
    "<img src=http://aco.ifmo.ru/el_books/numerical_methods/lectures/images/glava2_clip_image012.png>\n",
    "\n",
    "Для равноотстоящих узлов имеем\n",
    "$$\n",
    "  \\int_a^b f(x) dx \\approx h \\, \n",
    "     \\left( \\frac{f(x_0) + f(x_n)}{2} + f(x_1) + f(x_2) + \\dots + f(x_{n-1}) \\right).\n",
    "$$\n",
    "Если $f(x)$ дважды непрерывно дифференцируема на отрезке $[a, b]$, \n",
    "то для остаточного члена формулы трапеций справедливо:\n",
    "$$\n",
    "  R = \\frac{b-a}{12} f''(\\xi)\\, h^2, \\quad \\mbox{где $\\xi \\in [a, b]$}.\n",
    "$$\n",
    "*Также имеем 2-й порядок сходимости, как и для формулы центральных прямоугольников.*\n",
    "\n",
    "Метод трапеций реализует функция `trapz()` из модуля `scipy.integrate`.\n",
    "Есть два варианта обращения к ней: с одним входным аргументом\n",
    "и с двумя. Если вектор `x` содержит узлы, а `y` - значения функции в этих узлах, то\n",
    "`trapz(y, x)`\n",
    "численно вычисляет значение определенного интеграла по формуле трапеций.\n",
    "\n",
    "Вычислим интеграл ($**$) по формуле трапеций введем:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.466935060243005"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.finfo(float).eps\n",
    "b = 3\n",
    "n = 100\n",
    "h = (b - a)/n\n",
    "x = np.linspace(a, b, n + 1)\n",
    "y = x/np.sin(x)\n",
    "I_trapz = integrate.trapz(y, x)\n",
    "I_trapz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В точке $0$ функция имеет (устранимую) особенность,\n",
    "поэтому мы отступили от $0$ на машинное эпсилон.\n",
    "\n",
    "Вариант функции `trapz` вида `trapz(y)`\n",
    "предназначен для интегрирования по формуле трапеций с равноотстоящими узлами.\n",
    "Чтобы получить численное значение интеграла, \n",
    "необходимо значение `trapz(y)` домножить на шаг интегрирования $h$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.466935060243006"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate.trapz(y)*h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Правило Симпсона\n",
    "\n",
    "Пусть точки $x_k$ ($k=0,1,\\dots, n$) \n",
    "удовлетворяют условиям ($*$), $n$ четно.\n",
    "Заменяя на каждом отрезке $[x_k, x_{k+2}]$ $(k=0,2,4,\\dots,n)$\n",
    "подынтегральную функцию $f(x)$ ее интерполяционным многочленом $f_k(x)$ $2$-го порядка,\n",
    "построенным по узлам $x_k$, $x_{k+1}$, $x_{k+2}$, \n",
    "и интегрируя эти многочлены, мы приходим\n",
    "к *формуле Симпсона*. В частности, для случая равноотстоящих узлов\n",
    "получим:\n",
    "$$\n",
    "  \\int_a^b f(x) dx \\approx \n",
    "\\frac{h}{3} \n",
    "\\Bigl(\n",
    "  y_0 + y_n\n",
    "  + 2\\bigl(f(x_2) + f(x_4) + \\dots + f(x_{n-2})\\bigr)\n",
    "  + 4\\bigl(f(x_1) + f(x_3) + \\dots + f(x_{n-1})\\bigr)\n",
    "\\Bigr).\n",
    "$$\n",
    "\n",
    "<img src='http://aco.ifmo.ru/el_books/numerical_methods/lectures/images/glava2_clip_image025.png'>\n",
    "\n",
    "Если функция $f(x)$ имеет непрерывную на отрезке $[a, b]$\n",
    "производную $4$-го порядка, то для остаточного члена формулы Симпсона справедливо:\n",
    "$$\n",
    "  R = \\frac{b-a}{180} f^{IV}(\\xi)\\, h^4, \\quad \\mbox{где $\\xi \\in [a, b]$}.\n",
    "$$\n",
    "Обратите внимание, что *поряок сходимости для формулы Симпсона составляет 4, т.е. эта формула намного точнее формул прямоугольников и трапеций*.\n",
    "Метод Симпсона реализован в функции `simps` модуля `scipy.integrate`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.455462165555241"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "I_simps = integrate.simps(y, x)\n",
    "I_simps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.455462165555243"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate.simps(y)*h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Значения, которые мы получили, разбивая отрезок интегрирования на 100 частей. Ожидаемо, правило Симпсона должно давать наиболее точный результат. Сравнивая его с относительно точным значением, полученным адаптивными функциями ниже, установим, что правилосимпсона со 100 узлами дало 4 верные значащие цифры ($8.455$).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rectangaular method  8.449457117179284\n",
      "Trapezoidal  method  8.466935060243005\n",
      "Simpson rule         8.455462165555241\n"
     ]
    }
   ],
   "source": [
    "print('Rectangaular method ', I_rect)\n",
    "print('Trapezoidal  method ', I_trapz)\n",
    "print('Simpson rule        ', I_simps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Автоматический выбор шага\n",
    "\n",
    "\n",
    "Как мы уже отмечали, для оценки точности вычисленного значения интеграла можно увеличивать количество узлов, например, вдвое, и сравнивать полученный ответ с предыдущим. Если они отличаются незначительно, то можно считать, что необходимаю точность получена. Иначе можно еще раз увеличить количество узлов вдвое и т.д.\n",
    "\n",
    "Более интеллектуальная стратегия - увеличивать количество точек не на всем отрезке интегрирования, а только на проблемных участках - такого сорта алгоритмы называются адаптивными, так как они подстариваются (адаптируются) под поведение функции на конкретном участке. Например, применим метод Симпсона к заданной функции для $3$ и $5$ точек и сравним результаты. Если разница превышает некоторый порог, то разобьем отрезок интегрирования на два равных отрезка и применим к каждому из них тот же алгорритм рекурсивно.\n",
    "\n",
    "Функция `quad` из модуля `scipy.integrate` реализует некоторый сложный адаптивный алгоритм численного вычисления интеграла.\n",
    "\n",
    "Форма обращения к функции `quad` следующая:\n",
    "```\n",
    "    I, abserr = quad(f, a, b)\n",
    "```\n",
    "здесь $a$, $b$ - начало и конец отрезка интегрирования (могут быть $\\pm$`np.Inf`), а\n",
    "`f` - подынтегральная функция, `abserr` - полученная абсолютная ошибка.\n",
    "\n",
    "Для функции `quad` можно задавать некоторые другие параметры, в частности,\n",
    "можно задать абсолютную `epsabs` и относительную `epsabs` ошибки для остановки алгоритма\n",
    "(по умолчанию они равны $1.49\\times 10^{-8}$). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8.45527024499398, 2.6426765123466225e-11)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def fun(x):\n",
    "    return x/np.sin(x)\n",
    "\n",
    "I, abserr = integrate.quad(fun, np.finfo(float).eps, 3)\n",
    "I, abserr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Интеграл вычислен с 10 верными значащими цифрами. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8.455270245966526, 0.0005948624018474938)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "I, abserr = integrate.quad(fun, np.finfo(float).eps, 3, epsabs = 1e-3)\n",
    "I, abserr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Выводы\n",
    "\n",
    "* Для достаточно гладких функций \n",
    "    * метод левых и правых прямоугольников имеет 1-й порядок сходимости,\n",
    "    * средних прямоугольников и метод трапеций - 2-й порядок, \n",
    "    * метом Симпсона - самый точный - 4-й порядок сходимости.\n",
    "\n",
    "* Для вычисления интеграла с необхоимой точностью используют адаптивные алгоритмы.\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}