integrals_oskireeva.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.8.5"
    },
    "colab": {
      "name": "integrals_oskireeva.ipynb",
      "provenance": [],
      "collapsed_sections": []
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "skVfdh-n_F7u"
      },
      "source": [
        "Before you turn this problem in, make sure everything runs as expected. First, **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart) and then **run all cells** (in the menubar, select Cell$\\rightarrow$Run All).\n",
        "\n",
        "Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\", as well as your name and collaborators below:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mu5C-jj4_F70"
      },
      "source": [
        "NAME = \"\"\n",
        "COLLABORATORS = \"\""
      ],
      "execution_count": 39,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uD4-B0SY_F71"
      },
      "source": [
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "76725925c1c2a3fda23efe2691cfa949",
          "grade": false,
          "grade_id": "cell-4d48d7afa1c0c5b4",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "6U0uE7RQ_F72"
      },
      "source": [
        "# Вычисление определенных интегралов"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "9a9c3cddcc11d2a8cc7f46935462ad6d",
          "grade": false,
          "grade_id": "cell-fc379192341494d5",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "oyvFuym9_F72"
      },
      "source": [
        "## I. Двухточечная квадратура Гаусса\n",
        "\n",
        "Постройте квадратурную формулу Гаусса, интегрирующую точно многочлены степеней вплоть до третьей на интервале $[a, b]$. Заметим, что для этого достаточно построить _двухточечную_ квадратуру.\n",
        "\n",
        "Напишите функцию, которая возвращает веса, $w_1$ и $w_2$, и узлы, $x_1$ и $x_2$, двухточечной квадратурной формулы Гаусса с весовой функцией $w(x) = 1$, т.е., интегралы вида\n",
        "\n",
        "$$\n",
        "\\int_a^b\\! f(x)\\, dx \\approx w_1 f(x_1) + w_2 f(x_2)\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "904166366d67e67506e3087dd19ce395",
          "grade": false,
          "grade_id": "cell-cf3c018fa371396b",
          "locked": false,
          "schema_version": 3,
          "solution": true,
          "task": false
        },
        "id": "mmBEbw5t_F73"
      },
      "source": [
        "from math import sqrt\n",
        "def gauss_2(a, b):\n",
        "    r\"\"\"Return nodes and weights for a two-point Gauss quadrature on [a, b].\n",
        "    \n",
        "    Parameters\n",
        "    ----------\n",
        "    a, b : float\n",
        "       Boundaries of the integration interval\n",
        "       \n",
        "    Returns\n",
        "    -------\n",
        "    x1, x2, w1, w2 : float\n",
        "       Nodes and weights of the quadrature.\n",
        "    \"\"\"\n",
        "    w1 = w2 = (b - a) / 2\n",
        "    x1 = x2 = w1\n",
        "    tmp = (b - a) / (2 * sqrt(3))\n",
        "    x1 -= tmp\n",
        "    x2 += tmp\n",
        "    return x1, x2, w1, w2"
      ],
      "execution_count": 40,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "0359e8b293f2fb6d01a58c9fedf8394d",
          "grade": true,
          "grade_id": "cell-b6739a7f6045a9ec",
          "locked": true,
          "points": 5,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "DWAGd-vs_F73"
      },
      "source": [
        "from numpy.testing import assert_allclose\n",
        "\n",
        "x1, x2, w1, w2 = gauss_2(0, 1)\n",
        "\n",
        "def f(x, n): \n",
        "    return x**n\n",
        "\n",
        "for n in [0, 1, 2, 3]:\n",
        "    assert_allclose(w1*f(x1, n=n) + w2*f(x2, n=n),\n",
        "                    1./(n+1), atol=1e-14)"
      ],
      "execution_count": 41,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "f0b361e7c3254cdb24248df735a6a6ef",
          "grade": false,
          "grade_id": "cell-08fc4ae6160d61c1",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "P2SvXv6L_F74"
      },
      "source": [
        "## II. Вычитание интегрируемых сингулярностей\n",
        "\n",
        "\n",
        "Вычислите определённый интеграл методом трапеций с вычитанием сингулярности\n",
        "\n",
        "$$\n",
        "I = \\int_{0}^{1}\\frac{e^x}{\\sqrt{x(1-x)}}dx.\n",
        "$$\n",
        "\n",
        "Вам могут пригодиться значения следующих определенных интегралов:\n",
        "\n",
        "$$\n",
        "\\int_0^1 \\frac{1}{\\sqrt{x (1-x)}} \\, dx=\\pi,\\quad \\int_0^1 \\frac{x}{\\sqrt{x (1-x)}} \\, dx=\\pi/2.\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XEHHV2nW_F74"
      },
      "source": [
        "Преобразуйте данный интеграл, вычитая сингулярности. Выпишите расчетные формулы.\n",
        "(5 баллов, защита)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "85546100a8f30436f5555da4da442cca",
          "grade": false,
          "grade_id": "cell-6c554eec0b99fcb4",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "h8px2y5G_F75"
      },
      "source": [
        "Cоставьте функцию, возвращающую значение интеграла методом трапеций."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "094424dc42e9c4d3807980eaea9631f0",
          "grade": false,
          "grade_id": "cell-7fb345ac8e79dced",
          "locked": false,
          "schema_version": 3,
          "solution": true,
          "task": false
        },
        "id": "flkdRlgX_F76"
      },
      "source": [
        "def integ(npts=10):\n",
        "    \"\"\"Compute the value of the integral above.\n",
        "    \n",
        "    Subtract the singularities and use the trapezoid rule. \n",
        "    \n",
        "    Parameters\n",
        "    ----------\n",
        "    npts : int\n",
        "        The number of points for the trapezoid rule\n",
        "        \n",
        "    Returns\n",
        "    -------\n",
        "    I : float\n",
        "       The computed value of the integral\n",
        "    \"\"\"\n",
        "    h = 1.0/npts\n",
        "    f = lambda x: (1-x)**(5/2)\n",
        "    s = 3*np.pi/2 + 6/13\n",
        "    x = 0\n",
        "    for i in range(npts-1):\n",
        "        x+= h\n",
        "        s+=h*(f(x-h)+f(x))/2\n",
        "    print(s)\n",
        "    return s"
      ],
      "execution_count": 42,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "6d0b77bbebc9c1ead838b263af7ac7ae",
          "grade": true,
          "grade_id": "cell-e24dd8dc92f2c8bd",
          "locked": true,
          "points": 5,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "0NioWcse_F76"
      },
      "source": [
        "# this is a test to check your computed value\n",
        "from numpy.testing import assert_allclose\n"
      ],
      "execution_count": 43,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "ba059a5eaf1c17a3aa3e98590a0a6f9b",
          "grade": false,
          "grade_id": "cell-70721ae1c137e72e",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "pL9tVy-w_F77"
      },
      "source": [
        "## III. Интеграл от осциллирующей функции\n",
        "\n",
        "Рассмотрим определенный интеграл\n",
        "\n",
        "$$\n",
        "I = \\int_0^\\infty\\! \\frac{\\sin(x) \\cos{(\\cos{(x)})}}{x}\\,dx\n",
        "$$\n",
        "\n",
        "Вычислите значение данного интеграла с относительной точностью $10^{-6}$. Для упрощения задачи можете воспользоваться функционалом `scipy.integrate.quad`.\n",
        "\n",
        "Заметим, что \"из коробки\" интеграл вычислить не удается, и нужно что-то придумать."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "11f9ea8995c2fc47c6518788c1e3a53b",
          "grade": false,
          "grade_id": "cell-515f8f19764eb256",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "aoIBtxfm_F77",
        "outputId": "cb1197fe-d2a4-4e82-97dd-cd612f15d6d2"
      },
      "source": [
        "from math import sin, cos\n",
        "\n",
        "from scipy.integrate import quad\n",
        "quad(lambda x: sin(x) * cos(cos(x)) / x, 0, float('inf'))"
      ],
      "execution_count": 44,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:4: IntegrationWarning: The integral is probably divergent, or slowly convergent.\n",
            "  after removing the cwd from sys.path.\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(1.9653912540956746, 4.089174284042278)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 44
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "markdown",
          "checksum": "53e2ba2b52667aada8dc1a7f43391ef4",
          "grade": false,
          "grade_id": "cell-e34178c9823af271",
          "locked": true,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "jenngRMK_F77"
      },
      "source": [
        "Напишите функцию, которая возвращает значение данного интеграла. Только само значение интеграла, без оценки погрешности.\n",
        "(оборачиваем интеграл в функцию только для удобства автопроверки)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "deletable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "f9fa2326b0f5d084ce9385ccd5a2fbd7",
          "grade": false,
          "grade_id": "cell-ddd4e0b6d148fcfd",
          "locked": false,
          "schema_version": 3,
          "solution": true,
          "task": false
        },
        "id": "l3xLRnlL_F78"
      },
      "source": [
        "from scipy.integrate import quad\n",
        "import numpy as np\n",
        "from scipy.special import jv\n",
        "\n",
        "def integ():\n",
        "  \"\"\"Return a single float, the computed value of the integral.\"\"\"\n",
        "  a = jv(0,1)*0.6247132564277136 + 2 * jv(2,1) * 0.4512847288161426 + 2 * jv(4,1) * 0.1493606415273971 + 2 * jv(6,1) * (-0.04766731534829027)\n",
        "  return quad(lambda x: sin(x) * cos(cos(x)) / x, 0, 1)[0] + a"
      ],
      "execution_count": 45,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "deletable": false,
        "editable": false,
        "nbgrader": {
          "cell_type": "code",
          "checksum": "28f736930b29b8de099e68693bd5f160",
          "grade": true,
          "grade_id": "cell-b341c2ce36bc957d",
          "locked": true,
          "points": 10,
          "schema_version": 3,
          "solution": false,
          "task": false
        },
        "id": "VVmem0f5_F78"
      },
      "source": [
        "from numpy.testing import assert_allclose\n",
        "\n",
        "# this is a test cell, keep it intact\n"
      ],
      "execution_count": 46,
      "outputs": []
    }
  ]
}