Runge Kutta Methods

H
{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "Runga Kutta methods are used to approximate the value to simultaneous non-linear differential equations. <br>\n",
        "More about Runga Kutta Methods [here](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods).\n",
        "\n",
        "Here we will try to solve the following three equations with unknown variables S, I and R.\n",
        "\n",
        "\\begin{equation}\n",
        "\\frac{d S(t)}{d t}=A-\\mu S(t)-\\frac{r I}{N_0} S(t)\n",
        "\\end{equation}\n",
        "<br>\n",
        "$$\n",
        "\\frac{d I(t)}{d t}=\\frac{r I(t)}{N_0} S(t)-(\\mu+a) I(t)\n",
        "$$\n",
        "<br>\n",
        "$$\n",
        "\\frac{d R(t)}{d t}=a I(t)-\\mu R(t)\n",
        "$$\n",
        "<br>\n",
        "&emsp;&emsp;&emsp;&emsp;&emsp;&emsp;&emsp;&emsp;&emsp;\n",
        "where, $\\mathrm{S}(\\mathrm{t})+\\mathrm{I}(\\mathrm{t})+\\mathrm{R}(\\mathrm{t})=\\mathrm{N}(\\mathrm{t})$\n",
        "<br><br>\n",
        "$$\n",
        "=>\\frac{N(t)}{d t}=A-\\mu N\n",
        "$$\n",
        "\n",
        "These equations have a physical significance also, they represent the SIR model for the transmission of tuberclosis as documented [here](https://drive.google.com/file/d/1t2Rgh1jdEmT_aJ7RKAkDf0ZH_lhFM2sp/view?usp=sharing). \n",
        "\n",
        "S, I and R represent the following three kinds of people:\n",
        "\n",
        "\n",
        "1.   S: The first category is individuals who have not been infected with the TB virus, sometimes known as the suspicion group.  \n",
        "2.   I: The second category consists of people who have been infected with the TB virus.\n",
        "3.   R: The last category consists of people who recovered or died\n",
        "after being infected with the TB virus.\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "Other than S, I and R, all other terms are constants, and the equations can be represented as follows:"
      ],
      "metadata": {
        "id": "JCXtWuN-BYzF"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def dS_dt(S, I, R, t):\n",
        "    return (1449401 -  0.001167*S - (0.5/1449401)*I*S)"
      ],
      "metadata": {
        "id": "NqEw5s-KgR1g"
      },
      "execution_count": 18,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def dI_dt(S, I, R, t):\n",
        "    return ((0.5/1449401)*I*S - (0.001167 + 0.111111)*I)"
      ],
      "metadata": {
        "id": "TH77cRchhc6-"
      },
      "execution_count": 19,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def dR_dt(S, I, R, t):\n",
        "    return  (0.111111*I - 0.001167*R)"
      ],
      "metadata": {
        "id": "GxXtTA92kJRB"
      },
      "execution_count": 20,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "t = int(input())  # Enter a required value of time where you want the value of S, I and R."
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "snFbzheU5k8S",
        "outputId": "855bae9e-af27-48ad-b2e3-327277207cf8"
      },
      "execution_count": 22,
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "10\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "S_data = []\n",
        "I_data = []\n",
        "R_data = []\n",
        "SIR_data = []\n",
        "N_data = []\n",
        "Time = []"
      ],
      "metadata": {
        "id": "KVl7jTB_Q6nq"
      },
      "execution_count": 23,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Runge-Kutta method implementation with S,I,R as dependent variables and t as independent variable**"
      ],
      "metadata": {
        "id": "uWlRx0N7sdwt"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def rungeKutta(S0, I0, R0, t, h):\n",
        "\n",
        "  # No of iterations = n\n",
        "  # Step size = h\n",
        "    n = (int)((t)/h)\n",
        "  \n",
        "  # setting initial values\n",
        "    S = S0\n",
        "    I = I0\n",
        "    R = R0\n",
        "\n",
        "  # implementation of runge-kutta\n",
        "    for i in range(1, n + 1):\n",
        "        kS1 = dS_dt(S, I, R, t)\n",
        "        kI1 = dI_dt(S, I, R, t)\n",
        "        kR1 = dR_dt(S, I, R, t)\n",
        "        \n",
        "        kS2 = dS_dt(S + 0.5 * h * kS1, I + 0.5 * h * kI1, R + 0.5 * h * kR1, t + 0.5*h)\n",
        "        kI2 = dI_dt(S + 0.5 * h * kS1, I + 0.5 * h * kI1, R + 0.5 * h * kR1, t + 0.5*h)\n",
        "        kR2 = dR_dt(S + 0.5 * h * kS1, I + 0.5 * h * kI1, R + 0.5 * h * kR1, t + 0.5*h)\n",
        "\n",
        "        kS3 = dS_dt(S + 0.5 * h * kS2, I + 0.5 * h * kI2, R + 0.5 * h * kR2, t + 0.5*h)\n",
        "        kI3 = dI_dt(S + 0.5 * h * kS2, I + 0.5 * h * kI2, R + 0.5 * h * kR2, t + 0.5*h)\n",
        "        kR3 = dR_dt(S + 0.5 * h * kS2, I + 0.5 * h * kI2, R + 0.5 * h * kR2, t + 0.5*h)\n",
        "\n",
        "        kS4 = dS_dt(S + kS3 * h, I + kI3 * h, R + kR3 * h, t + h)\n",
        "        kI4 = dI_dt(S + kS3 * h, I + kI3 * h, R + kR3 * h, t + h)\n",
        "        kR4 = dR_dt(S + kS3 * h, I + kI3 * h, R + kR3 * h, t + h)\n",
        " \n",
        "        # Updating S, I, R\n",
        "        S = S + (1.0 / 6.0)*(kS1 + 2 * kS2 + 2 * kS3 + kS4)*h\n",
        "        I = I + (1.0 / 6.0)*(kI1 + 2 * kI2 + 2 * kI3 + kI4)*h\n",
        "        R = R + (1.0 / 6.0)*(kR1 + 2 * kR2 + 2 * kR3 + kR4)*h\n",
        "\n",
        "        S_data.append(S)\n",
        "        I_data.append(I)\n",
        "        R_data.append(R)\n",
        "        SIR_data.append(S+I+R)\n",
        "        Time.append(t)\n",
        "\n",
        "        # Updating value of t\n",
        "        t += h\n",
        "    \n",
        "    # printing N(t) at desired t (N(t)=S(t)+I(t)+R(t))\n",
        "    print(\"N(t) = \" ,S+I+R)\n",
        "    return [S,I,R]"
      ],
      "metadata": {
        "id": "-zl-oP5SOxbO"
      },
      "execution_count": 24,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Driver Code\n",
        "S0 = 1446093\n",
        "I0 = 1885\n",
        "R0 = 1423\n",
        "h = 0.01\n",
        "print ('The value of S(t), I(t) and R(t) is:', rungeKutta(S0, I0, R0, t, h)) "
      ],
      "metadata": {
        "id": "KmEB8JRcIlCw",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "e8af4878-e13a-4ebe-e820-59c5e12918d5"
      },
      "execution_count": 25,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "N(t) =  15842350.284962129\n",
            "The value of S(t), I(t) and R(t) is: [410394.8072337986, 10323985.625711355, 5107969.852016973]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Runge-Kutta method implementation with N as dependent variable and t as independent variable**"
      ],
      "metadata": {
        "id": "YaT1SXaUsvc5"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We implemented it just to prove that the equation N(t) = S(t)+I(t)+R(t) is consistent with the values of S,I,R we got above for the same value of t."
      ],
      "metadata": {
        "id": "vPhSPAG_s1z5"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def dN_dt(t, N):\n",
        "\treturn (1449401-0.001167*N)\n",
        "\n",
        "# implementing runge kutta method for ODE: dN/dt = A - m*N where m is death rate\n",
        "def rungeKutta(t0, A, t, h):\n",
        "\n",
        "\t# no. of iterations = n\n",
        "\t# step size = h\n",
        "\tn = (int)((t - t0)/h)\n",
        "\t\n",
        "\tN = A\n",
        "\tfor i in range(1, n + 1):\n",
        "\t\tk1 = dN_dt(t0, N)\n",
        "\t\tk2 = dN_dt(t0 + 0.5 * h, N + 0.5 * h * k1)\n",
        "\t\tk3 = dN_dt(t0 + 0.5 * h, N + 0.5 * h * k2)\n",
        "\t\tk4 = dN_dt(t0 + h, N + h * k3)\n",
        "\n",
        "\t\t# Update value of y\n",
        "\t\tN = N + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4)*h\n",
        "\t\tN_data.append(N)\n",
        "\n",
        "\t\t# Update next value of t0\n",
        "\t\tt0 += h\n",
        "\treturn N\n",
        "\n",
        "# Driver Code\n",
        "t0 = 0\n",
        "A = 1449401\n",
        "h = 0.01\n",
        "print ('The value of N(t) at required t is:', rungeKutta(t0, A, t, h))"
      ],
      "metadata": {
        "id": "ydrVeL51kVuu",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "928c42e2-27e7-4567-8f2b-5b143d0391ac"
      },
      "execution_count": 26,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "The value of N(t) at required t is: 15842350.284962097\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# print(S_data)"
      ],
      "metadata": {
        "id": "qrXWMyh7Wy03"
      },
      "execution_count": 27,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# print(I_data)"
      ],
      "metadata": {
        "id": "zSjyo5IDkMeT"
      },
      "execution_count": 28,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# print(R_data)"
      ],
      "metadata": {
        "id": "J8oCbC9FkNr8"
      },
      "execution_count": 29,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# print(Time)"
      ],
      "metadata": {
        "id": "rGsiEb3zkUpd"
      },
      "execution_count": 30,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# print(N_data)"
      ],
      "metadata": {
        "id": "59T-ht2x-9h5"
      },
      "execution_count": 31,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Variation of S, I, R and N with time.**"
      ],
      "metadata": {
        "id": "HoWaiRbBK5C3"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import math\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "t = np.linspace(20,38, 1000)\n",
        "plt.scatter(Time, S_data)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 293
        },
        "id": "RouYbro0S9Cl",
        "outputId": "2c037811-2167-4a4b-f88f-8422f84d9a11"
      },
      "execution_count": 32,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.collections.PathCollection at 0x7ff40cbf3590>"
            ]
          },
          "metadata": {},
          "execution_count": 32
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEDCAYAAAAcI05xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAYQ0lEQVR4nO3df5BddXnH8c+z+YEagxjZpE4gjaYBBs0v3ZrQWIs4a1GsrOMgOptOO2XIjOM4xigKZkd0SloUZpv+0R+TFYsdUgoI3DJi0dSypWaW6Ib8WELQNDYiK5LVQMFISSBP/7gXumz2nvO9u+f3fb9mMmz2nmyeQzafnPuc5/s95u4CABRXR94FAACiEdQAUHAENQAUHEENAAVHUANAwRHUAFBwqQW1mX3dzI6Y2cOBx3/EzB4xs/1m9k9p1QUAZWNpzVGb2bsk/VrSP7r7W2OOXSrpdkkXuftTZjbf3Y+kUhgAlExqV9Tu/oCko+M/Z2ZLzOw+M9tlZv9pZuc1XrpS0t+4+1ONX0tIA0BD1j3qrZI+6e5vl/RZSX/b+Pw5ks4xsx1m9qCZXZxxXQBQWDOz+o3M7LWSfk/SHWb20qdPG1fHUkkXSjpL0gNmtszdn86qPgAoqsyCWvWr96fdfeUkrz0uaae7n5D032b2Y9WD+4cZ1gcAhZRZ68Pdn1E9hC+TJKtb0Xi5pvrVtMzsTNVbIT/JqjYAKLI0x/NulTQk6Vwze9zMrpDUK+kKM9srab+kSxuHf0fSr8zsEUn3S7rK3X+VVm0AUCapjecBAJLBykQAKLhUbiaeeeaZvnjx4jS+NABU0q5du37p7p2TvRYb1GZ2rqTbxn3qzZK+6O5bmv2axYsXa3h4uOVCAaBdmdlPm70WG9Tu/iNJKxtfaIakUUl3J1YdACBSqz3q90g65O5Nkx8AkKxWg/qjkm6d7AUzW29mw2Y2PDY2Nv3KAACSWghqM5st6YOS7pjsdXff6u5d7t7V2TlpPxwAMAWtXFG/T9JD7v5kWsUAAE7Vynjex9Sk7QG0qndgSDsOHZ30tZkdphsvW6GeVQszrgoopqCViWY2R9Jjkt7s7v8Td3xXV5cznoeJ+mojuuXBx1r6NWuXzNO2Ky9IqSKgOMxsl7t3TfZaUOvD3Y+5+xtCQhqYqK82osVX39tySEvSjkNH9aar71Vt92gKlQHlwBJypKa2e3TKAT2eS9pw2x71DgwlUxhQMgQ1UtE7MKQNt+1J9GvuOHRU3f2DiX5NoAwIaiSuu3+w6Y3C6Tp45BhX1mg7BDUStXrzdh08cizV32PHoaP0rNFWCGok5rxN39aTzx5v+dctnT9H69YsaunXbEy4rQIUWSoPDmA8r/2ct+nb+t8Xw7+XOkzq/8jKU2ala7tHg3vbS+fP0faNF7ZSJlBY0x7PA5p5abIjNKQ7TNpy+Ur95C8vmXRBS8+qhTp8/SVau2Re7Nc6eOSY+mojLdcMlA1BjSnrq420NNmxdP6cpgE90bYrLwgK6+mO/gFlQFBjSmq7R1sKybVL5rXcpth25QVaOn9O7HFcVaPqCGpMyadbuJLecvnKKS8D377xQs3ssMhjuKpG1RHUaNnya+9TSEf6VTNMh68Pa3VEufGyFbHHcFWNKiOo0ZLl196nZ55/Mfa400+boUc3vz+R37Nn1cLY8T2uqlFlBDWCrd68PSikF8ydrX1fvjjR3/u6nmWxx3BVjaoiqBGkd2AoaDHLgrmztXNTdyo1cFWNdkVQI1ZfbSRo747TT5uRWkhL9avquBuLLC1HFRHUiBS62f/pp81IvN0xmbgbi9fctS/1GoCsEdRoqpVZ6SxCWqrfWDxtZvNv2+dOnMykDiBLBDWa2nh72Kz0lstXplzJK33lw8sjX6f9gaohqDGp7v5BnQwYll63ZlHmD6GN+/2+dM/+jCoBskFQ4xTd/YNBe0qvXTIvaGwuDXNmz2j62tPPnciwEiB9BDVeoa82EhTSS+fPyfXp4Js/FP0PBDPVqBKCGq8QcvNwwdzZue8DHdf+YKYaVRIU1GZ2hpl908weNbMDZpbfpRRSs3rz9qDj0pyVbsXrXzMr8nVuKqIqQq+o/1rSfe5+nqQVkg6kVxLy0N0/GLTyMOsJjyjX/tFbIl/npiKqIjaozex1kt4l6SZJcvfj7v502oUhO70DQ0F96TwmPKLEzVRzUxFVEXJF/SZJY5L+wcx2m9nXzOyU3dzNbL2ZDZvZ8NjYWOKFIh213aNBy8PznPCIwkw12kFIUM+U9DZJf+fuqyQdk3T1xIPcfau7d7l7V2dnZ8JlIi0hi1oWzJ2d64RHFGaq0Q5CgvpxSY+7+87Gz7+penCj5EIXtRTl5mEzUTcVaX+gCmKD2t1/IelnZnZu41PvkfRIqlUhdaF96SLdPGwm7qYiUHahUx+flLTNzPZJWinpL9IrCWkL3ba0aDcPm4mrkT41yi4oqN19T6P/vNzde9z9qbQLQzpCd8RbOn9OIW8eTgV9apQdKxPbzOfvDNuvOe+Vh62iT40qI6jbSF9tRM+/EL9fcxn60hPRp0aVEdRtIrTlUZa+9ET0qVFlBHWbCJmXLuqillBR7Q/61CgzgroNhMxLd0iFXdQSKqr9QZ8aZUZQV1xt92jQvHR/CfvSE9H+QFUR1BUX2vIoY1+6VbQ/UFYEdYWt3rw9tuWR95NaksaYHqqIoK6o3oGhoP2lyzYvHYcxPVQRQV1BoVuXlnFeOk47tHDQfgjqCmr3vnSHNX+NG4ooI4K6YnoHhmL70kXeXzoJUefPDUWUEUFdIaG74hV9f+npWnjGq5u+xg1FlBFBXRGhS8Sr2Jee6Ko/PDf+IKBECOqKCNkVb+n8OZXtS4/HwhdUDUFdAaG74lVtFC8K+36gSgjqkqPlMTn2/UCVENQl99k79sYeU+VRvGba7XxRbQR1iXX3D+qFmFm8KuyKlwb61CgTgrqk+mojbbMr3lTRp0ZVENQlFdKXbseWx3j0qVEVBHUJdfcPxh5TtV3xpqKd/5FCtQQFtZkdNrMRM9tjZsNpF4XmegeGYlseHWqvUbwo7PuBKmjlivrd7r7S3btSqwaRQnfFa+e+9ETs+4EqoPVRIiGrD9u9Lz0R+36gCkKD2iV918x2mdn6yQ4ws/VmNmxmw2NjY8lVCElhqw8ZxTsV+36gCkKD+p3u/jZJ75P0CTN718QD3H2ru3e5e1dnZ2eiRba70NWHtDxOxbsLVEFQULv7aOO/RyTdLekdaRaFVwpZfbhuzSJCCaio2KA2szlmNveljyW9V9LDaReGupDVh7M6pOt6lmVUUbUw+YEyCLmiXiDp+2a2V9IPJN3r7velWxak8NWHN1xGyyMKKxRRdrFB7e4/cfcVjR9vcffNWRQGVh8mhRWKKDvG8wqqd2Ao9hhWH4bhHzKUHUFdQCHPPmT1YXLoU6PoCOqCYRQvHfSpUWYEdcGw+jAd9KlRZgR1gbD6MD38w4YyI6gLgpZH+thJD2VFUBcEzz5MX9S6oRu+86PsCgFaRFAXQO/AEM8+zEDUTnqjTz+XYSVAawjqnLHHdHaidtKL6IoAuSOoc7bx9j2xx9DySEbU/8Po9zNAvgjqHPUODEX2TSVWHwIgqHMT0vJg9WHyoiY/gKIiqHMS0vKgL528qHcwjOihqAjqHHT3Dwa1POhLJy9q8oOl5Cgqgjpjtd2jQXtM0/JIR9TkB0vJUVQEdcZCFrZsoeWRGt6loIwI6gyFPFaLUbz0sZQcZUNQZ6R3YCio5cEoXvpYSo6yIagzELr6kJZHNlhKjrIhqDPAHtPFwlJylA1BnTL2mC4elpKjbAjqFLHHdDlxQxFFExzUZjbDzHab2bfSLKhKQkbx1q1ZRMsjBzxDEWXSyhX1pyQdSKuQqumrjcSO4s3qkK7rWZZRRRiPZyiiTIKC2szOknSJpK+lW051hLQ8briMlkdeeBeDMgm9ot4i6XOSmt4VM7P1ZjZsZsNjY2OJFFdWqzdvjz2GKY/8sfAFZREb1Gb2AUlH3H1X1HHuvtXdu9y9q7OzM7ECy6Z3YEhPPns88hj2mC4GFr6gLEKuqNdK+qCZHZb0z5IuMrNbUq2qpEIXtrDhUjGw8AVlERvU7n6Nu5/l7oslfVTSv7v7utQrK6GQPabXrVmUQSUIwcIXlAVz1AkJ2WN6wdzZTHkUCAtfUBYzWznY3QclDaZSSYmF7jG9c1N3BtUAqBquqBPAHtPlxTMUUQYE9TSxx3S58QxFlAFBPQ3sMV1+PEMRZUBQTxF7TFcDz1BEGRDUU8Qe09XAnw/KgKCeAvaYrhaWkqPoCOoWscd09bCUHEVHULeIPaarh6XkKDqCugXsMV1NLCVH0RHULWCP6WpiKTmKjqAOxB7T7YsbisgbQR2APaarj2coosgI6hjsMd0eeIYiioygjsEe0+2BlhWKjKCOwB7TAIqAoG6CPaYBFAVB3QR7TGM8Jj+QJ4J6Euwx3Z6Y/EBREdQT9NVG2GO6TTH5gaIiqCcIWX1Iy6OaeIeEoiKox+kdGIo9hpZHtbHlKYqIoG7oq43ELmxhj+nqY8tTFFFsUJvZq8zsB2a218z2m9mXsygsS+wxjZew5SmKKOSK+nlJF7n7CkkrJV1sZmvSLStbIaN4tDzaQ9SWp0BeYoPa637d+Omsxo/K7P4Yssc0LY/2EfePMX1q5CGoR21mM8xsj6Qjkra7+85JjllvZsNmNjw2NpZ0namh5YFW0KdGHoKC2t1fdPeVks6S9A4ze+skx2x19y537+rs7Ey6zlR09w/GHkPLo/1ELXyhT408tDT14e5PS7pf0sXplJOd3oGh2IUt7DHdnqIWvvBoLuQhZOqj08zOaHz8akndkh5Nu7A0scc0ovBoLhRNyBX1GyXdb2b7JP1Q9R71t9ItK12fv3Nf7DHsMY1muKGIrM2MO8Dd90lalUEtmeirjej5F05GHtMhniTe7l7/mll66jeT7+/xpXv2c98CmWqrlYksbEEoNmhCkbRVUIcsbFm3ZhFXS+B7AIXSNkHdOzAUu7BlVgctDwDF0xZBHTrlccNltDzw/9hJD0XRFkEd8iRxFrZgoqg3YDzxBVmqfFD3DgzFPkmchS2YTNROetxQRJYqHdQhLY8OsbAFk2MnPRRFpYM6pOXBKB6a6Vm1MHLJOH1qZKWyQd3dPxjU8qAvjShR30L0qZGVSgZ1bfdo0JPEaXkgDn1qFEElgzpkYQtPEkcI+tQogsoFdcjCFkbxEIonvqAIKhXUIU8Sl3isFloT9SAB+tTIQmWCOnTDJVoeaBUbNCFvlQnqkD2maXlgKvieQd4qEdS13aNBe0zT8kAa6FMjbZUIaha2IE88mRxpK31Qr968PXZhCy0PTBdPJkeeSh3UfbURPfns8chjaHkgCVE3FIG0lTqoeawWssI8NfJU2qDu7h+MPYaWB7LCPDXSVMqg7h0Yit3LY8Hc2bQ8kKioPjXz1EhTbFCb2dlmdr+ZPWJm+83sU1kU1kzoY7V2burOoBq0k7g+Ne0PpCXkivoFSZ9x9/MlrZH0CTM7P92ymgt9kjiQtLg2Gu0PpCU2qN39CXd/qPHxs5IOSMql8Ruy4VKHeJI40kP7A3loqUdtZoslrZK0c5LX1pvZsJkNj42NJVPdOKEtD6Y8kCbaH8hDcFCb2Wsl3Slpg7s/M/F1d9/q7l3u3tXZ2ZlkjZLCWx5MeSBNtD+Qh6CgNrNZqof0Nne/K92SThW6xzQtD2SB9geyFjL1YZJuknTA3fvTL+mVQlsejOIhK6xSRNZCrqjXSvpjSReZ2Z7Gj/enXNfLQjZcYsoDWWKVIrI2M+4Ad/++JMugllOEPkmclgeK5Ev37OdeCRJV2JWJIU8S7xBPEkc+6FMjS4UN6pApD0bxkBfG9JClQgZ1d/8gTxJHocV9711zV/yj4YBQhQvqvtpIbMtDYsoD+Ytqfzx34iRX1UhM4YKaJ4mjLOLaHyx+QVIKFdS9A0Oxx9DyQFHEfR9yUxFJKUxQ99VGYhe28FgtFE1U+0PipiKSUYigru0e5bFaKKW49gc3FZGEQgT1prtHYo+h5YEi6lm1UHNmz2j6+nMnTmZYDaqqEEF97PiLka/T8kCRbf5Q9MpY2h+YrkIEdRxaHiiyuHd6V90Rv18NEKUQQR21kQgtD5RBVPvjxEmuqjE9hQjq3ia73y2dP4eWB0ohrv3BVTWmoxBBfV3PMq1bs0gzrH5tPcNM69YsYsMllEbPqoWR7wy5qsZ0mHvMPqJT0NXV5cPDw4l/XaDI+mojkWOmr57VoQN//r4MK0KZmNkud++a7LVCXFEDVRC3LzqjepgqghpIUNzThkK2SQAmIqiBBMVdVe84dJReNVpGUAMJixrVk5gAQesIaiBhcaN6J07WbzwCoQhqIGE9qxbqtJnRf7VuefAxWiAIFhvUZvZ1MztiZg9nURBQBV/58PLYY2iBIFTIFfXNki5OuQ6gUnpWLdTaJfMij2ERDELFBrW7PyApekd/AKcI2f5gw21cVSMePWogRXFz1ZK0evP2DCpBmSUW1Ga23syGzWx4bGwsqS8LlNp1Pcs0syNqFxDpyWePsxAGkRILanff6u5d7t7V2dmZ1JcFSu/Gy1bEHrPj0FFG9tAUrQ8gZSE3FqX6yB5hjcmEjOfdKmlI0rlm9riZXZF+WUC1bLvyAi2YOzv2uFsefIw2CE4RMvXxMXd/o7vPcvez3P2mLAoDqmbnpu6g43YcOqru/sF0i0Gp0PoAMhQyBSJJB48cYxoELyOogQxd17NMS+fPCTr2yWePa/HV97IoBgQ1kLXtGy8MDmupvijmd77wbQK7jRHUQA5aDesXTro23LZHb76GK+x2RFADOWk1rCXppNevsBdffS/TIW2Eh9sCOesdGNKOQ9PbTmftknlBe4uguKIebktQAwVQ2z2a6AZNS+fP0faNFyb29ZA+ghooieXX3qdnnn8x1d+DEC8mghookb7aiG558LG8y5gULZb0ENRACRU5sMvitJkd+sqHl6tn1cK8S4lFUAMlRmCXz1TeeUQFNeN5QMFd17NMh6+/RFsuX6lZ/I0thR2HjiY6Pjkzsa8EIFU9qxa+4i18d/+gDh45lmNFiDLdkcvxCGqgpCZObiQxj41iIqiBiojqiRLi5UZQA21gqiN1td2juuqOPTpxMuGC2kDIU31CEdQAmprYFy+DIkzJJD1vzngeABQA43kAUGIENQAUHEENAAVHUANAwRHUAFBwqUx9mNmYpJ9O8ZefKemXCZZTBpxze+Ccq2865/vb7t452QupBPV0mNlwsxGVquKc2wPnXH1pnS+tDwAoOIIaAAquiEG9Ne8CcsA5twfOufpSOd/C9agBAK9UxCtqAMA4BDUAFFyuQW1mXzezI2b28LjPzTOz7WZ2sPHf1+dZY5KanO8NZvaome0zs7vN7Iw8a0zaZOc87rXPmJmb2Zl51JaWZudsZp9s/FnvN7Ov5lVfGpp8b680swfNbI+ZDZvZO/KsMWlmdraZ3W9mjzT+TD/V+HziGZb3FfXNki6e8LmrJX3P3ZdK+l7j51Vxs0493+2S3uruyyX9WNI1WReVspt16jnLzM6W9F5JVXy89s2acM5m9m5Jl0pa4e5vkXRjDnWl6Wad+uf8VUlfdveVkr7Y+HmVvCDpM+5+vqQ1kj5hZucrhQzLNajd/QFJE58PdKmkbzQ+/oaknkyLStFk5+vu33X3Fxo/fVDSWZkXlqImf8aS9FeSPiepcnezm5zzxyVd7+7PN445knlhKWpyzi7p9MbHr5P080yLSpm7P+HuDzU+flbSAUkLlUKG5X1FPZkF7v5E4+NfSFqQZzEZ+zNJ/5p3EWkzs0sljbr73rxrydA5kn7fzHaa2X+Y2e/mXVAGNki6wcx+pvo7iKq9W3yZmS2WtErSTqWQYUUM6pd5fXawcldckzGzTaq/ldqWdy1pMrPXSPqC6m+F28lMSfNUf4t8laTbzczyLSl1H5f0aXc/W9KnJd2Ucz2pMLPXSrpT0gZ3f2b8a0llWBGD+kkze6MkNf5bqbeIkzGzP5X0AUm9Xv3B9iWS3iRpr5kdVr3V85CZ/VauVaXvcUl3ed0PJJ1UfQOfKvsTSXc1Pr5DUqVuJkqSmc1SPaS3uftL55p4hhUxqO9R/Q9Yjf/+S461pM7MLla9V/tBd/9N3vWkzd1H3H2+uy9298WqB9jb3P0XOZeWtpqkd0uSmZ0jabaqv6vczyX9QePjiyQdzLGWxDXeEd0k6YC79497KfkMc/fcfki6VdITkk6o/hf2CklvUP1O6UFJ/yZpXp41ZnC+/yXpZ5L2NH78fd51pn3OE14/LOnMvOvM4M95tqRbJD0s6SFJF+VdZwbn/E5JuyTtVb13+/a860z4nN+peltj37i/v+9PI8NYQg4ABVfE1gcAYByCGgAKjqAGgIIjqAGg4AhqACg4ghoACo6gBoCC+z/OQAC01naN2AAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "plt.scatter(Time, I_data)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 293
        },
        "id": "nR549ji3ljC-",
        "outputId": "eef6d6cc-9fb1-4521-ba58-e250a72de7b4"
      },
      "execution_count": 33,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.collections.PathCollection at 0x7ff40cbf3e10>"
            ]
          },
          "metadata": {},
          "execution_count": 33
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEDCAYAAAAlRP8qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAVUUlEQVR4nO3df7Dd9V3n8ec7BJDtL4wJVUPSMJh2xUpB71IUf2AtNu04hLFLCQszdmWaWad0rKmsuDDYoo62OKnOiNsN2sEKAqG11+xsNLK1WKeTdAm9ARqQErFCLlViW1q3rUDK2z/OST29ufd8v/fe7znfH+f5mLlzz/l+P/ec93du8prPfX8/3++JzESS1H4r6i5AklQNA12SOsJAl6SOMNAlqSMMdEnqCANdkjqi1kCPiA9GxNMR8ZkSY98fEQf6X5+NiGfGUaMktUXUuQ49In4M+P/AhzLz1Yv4uXcA52bmz42sOElqmVpn6Jn5CeCLg9si4syI+IuIuD8i/iYi/uM8P3o5cMdYipSkllhZdwHz2AH8t8x8LCJeC/w+8LpjOyPiFcAZwF/VVJ8kNVKjAj0iXgz8MHB3RBzbfPKcYVuAD2fmN8ZZmyQ1XaMCnV4L6JnMPGfImC3A28dUjyS1RqOWLWbmV4C/j4hLAaLnNcf29/vp3w7sralESWqsupct3kEvnF8VEYcj4irgCuCqiHgAOAhsHviRLcCd6S0iJek4tS5blCRVp1EtF0nS0tV2UnT16tW5YcOGut5eklrp/vvv/+fMXDPfvtoCfcOGDezfv7+ut5ekVoqIf1hony0XSeoIA12SOsJAl6SOMNAlqSMMdEnqiKbdy0WSOmd6ZpZr7j7A8y/8+7aTV67gvW8+m0vOXVvZ+xjoklSh6ZlZtt11gBcKxj179AW27TwAUFmoG+iStETXTz/EbfueWPLPv5Bw055HDXRJGqeyM+/FeuqZr1f2Wga6JM0xX897VL771FMqey0DXdJEG2d4z7Ui4Jo3vKqy1ysM9Ij4IPDTwNOZ+ep59gfwu8CbgK8Bb83MT1dWoSRV6Ipb9vLJv/ti8cARC2D7W84Z+yqXW4HfAz60wP43Ahv7X68F/mf/uyTVqs7Z9zAXnLmK29/2Q5W/bmGgZ+YnImLDkCGbgQ/1P0VoX0ScGhHflZmfr6hGSSqlKbPvQVeev55fv+T7x/JeVfTQ1wJPDjw/3N92XKBHxFZgK8D69esreGtJk2q5SwartnJF8NuXvqbSFsqiaxjnm2XmDmAHwNTUlJ99J6m0pgX4OGfeZVUR6LPAuoHnp/e3SdKSNal90sTwnk8Vgb4LuDoi7qR3MvTL9s8lLUaTTl6O6oTlOJRZtngHcCGwOiIOA78KnAiQmR8AdtNbsniI3rLF/zqqYiV1Q1MCfBQ3yKpTmVUulxfsT+DtlVUkqXOaEuBtaZ0slVeKSqpcEwJ842kv4p5tF9ZXQA0MdEmVuGj7vTz29Fdre/+uz77LMNAlLUmdq1C61vuuioEuqbS6ZuFtXnkyTga6pAXVNQu3fbI0Brqkb1FHiDsDr4aBLmnsIW6Aj4aBLk2ocYa4AT4eBro0Qca1PtwAr4eBLk2AUa9OcRlhMxjoUkeNejbuLLx5DHSpY6ZnZtl21wGqznFn4c1noEsdMYogN8TbxUCXOqDKHvkk3tSqKwx0qcWqWnpoiHeDgS610PTMLO+868CyX8dL7LvFQJdaZrmzclendJeBLrXIa3/jHv7pX55b0s8a5N1noEstcfav/gVfefYbi/452yqTw0CXGm6p/XJn5JPHQJcabCn9coN8chnoUkMtNsxfevIJPPieTSOsSE23ou4CJB3v+umHFhXmF5y5yjCXgS41zfXTD3HbvidKj7/y/PW2WATYcpEaZXpmtnSYB/D+y87xPiv6JgNdapBtO8utZnn5S07iU9ddNOJq1Da2XKSGuGj7vbyQxeM2nvYiw1zzMtClBrjilr2l7pboTbQ0jIEu1azsihbDXEUMdKlGZU+CvvwlJxnmKlQq0CNiU0Q8GhGHIuLaefavj4iPR8RMRDwYEW+qvlSpe375Iw+WGmfPXGUUBnpEnADcDLwROAu4PCLOmjPsemBnZp4LbAF+v+pCpa6Znpnl2aPFHxj3O5edM4Zq1AVlZujnAYcy8/HMfA64E9g8Z0wCL+0/fhnwVHUlSt30S3c/UDjmyvPXu85cpZUJ9LXAkwPPD/e3DXo3cGVEHAZ2A++Y74UiYmtE7I+I/UeOHFlCuVI3XHHLXo4WrFG84MxV3vZWi1LVSdHLgVsz83TgTcAfR8Rxr52ZOzJzKjOn1qxZU9FbS+0yPTNbuKplBXg5vxatTKDPAusGnp/e3zboKmAnQGbuBb4NWF1FgVLXlDkRut2+uZagTKDfB2yMiDMi4iR6Jz13zRnzBPCTABHxvfQC3Z6KNEeZE6EXnLnKvrmWpDDQM/MocDWwB3iE3mqWgxFxY0Rc3B/2LuBtEfEAcAfw1swscRGzNFmu++hDQ/fbatFylLo5V2bupneyc3DbDQOPHwYuqLY0qVumZ2b56nPDPxPUVouWwytFpTEp6p2fuAJbLVoWA10agzK985sudXau5THQpTEomp2fcuIKZ+daNgNdGrEys/Pf/Jmzx1SNusxAl0asaGWLyxRVFQNdGrGilS0uU1RVDHRphK6fHj47v/L89WOqRJPAQJdG6PaCD6/w5luqkoEujcj0zCzDLpc+5UT/+6la/ouSRqRoqaIrW1Q1A10agaKlil4VqlEw0KUReM//Pjh0v1eFahQMdGkEvvS154fud3auUTDQpYpNz8z9/Jdv5VJFjYqBLlWsqN3iUkWNioEuVayo3SKNioEuVaio3XLqKSeOqRJNIgNdqlDRjbjeffH3jakSTSIDXapI0UfMec9zjZqBLlWk6GSoV4Zq1Ax0qSKuPVfdDHSpAp4MVRMY6FIFitotngzVOBjoUgWGtVs8GapxMdClEfNkqMbFQJeWqah/7uxc42KgS8tU1D+XxsVAl5ZpWP/c1S0aJwNdWoaidourWzROpQI9IjZFxKMRcSgirl1gzFsi4uGIOBgRf1JtmVIzFbVb7J9rnFYWDYiIE4CbgYuAw8B9EbErMx8eGLMR+BXggsz8UkScNqqCpSax3aImKTNDPw84lJmPZ+ZzwJ3A5jlj3gbcnJlfAsjMp6stU2oe2y1qmjKBvhZ4cuD54f62Qa8EXhkRn4yIfRGxab4XioitEbE/IvYfOXJkaRVLDWG7RU1T1UnRlcBG4ELgcuCWiDh17qDM3JGZU5k5tWbNmoreWqqH7RY1TZlAnwXWDTw/vb9t0GFgV2Y+n5l/D3yWXsBLE8l2i+pQJtDvAzZGxBkRcRKwBdg1Z8w0vdk5EbGaXgvm8QrrlBrFq0PVRIWBnplHgauBPcAjwM7MPBgRN0bExf1he4AvRMTDwMeBazLzC6MqWqqbV4eqiQqXLQJk5m5g95xtNww8TmBb/0vqPPvnaiKvFJUWyeWKaioDXVoklyuqqQx0aZFst6ipDHRpEWy3qMkMdGkRbLeoyQx0aRFst6jJDHSpIrZbVDcDXSrJq0PVdAa6VJJXh6rpDHSppGH987WnnjLGSqT5GehSBa55w6vqLkEy0KUy7J+rDQx0qQT752oDA10qwfXnagMDXSrg5f5qCwNdKnDTnkeH7rd/rqYw0KUCs898fcF9tlvUJAa6VCCG7LPdoiYx0KUhpmdmySH7bbeoSQx0aQiXK6pNDHRpCJcrqk0MdGmJ7J+raQx0aQFe7q+2MdClBRStP5eaxkCXFuD6c7WNgS4twPXnahsDXZqH68/VRga6NA/Xn6uNDHRpHq4/VxuVCvSI2BQRj0bEoYi4dsi4N0dERsRUdSVK4+XtctVWhYEeEScANwNvBM4CLo+Is+YZ9xLgF4BPVV2kNE5F7Rb752qqMjP084BDmfl4Zj4H3AlsnmfcrwHvBf61wvqksbPdorYqE+hrgScHnh/ub/umiPgBYF1m/p9hLxQRWyNif0TsP3LkyKKLlepmu0VNtuyTohGxAtgOvKtobGbuyMypzJxas2bNct9aqpyX+6vNygT6LLBu4Pnp/W3HvAR4NXBvRHwOOB/Y5YlRtdGw/rntFjVdmUC/D9gYEWdExEnAFmDXsZ2Z+eXMXJ2ZGzJzA7APuDgz94+kYmmEhvXPbbeo6QoDPTOPAlcDe4BHgJ2ZeTAiboyIi0ddoNQUtlvUdCvLDMrM3cDuOdtuWGDshcsvSxq/ov651HReKSr1ebm/2s5Al/pcf662M9AlvNxf3WCgS3i5v7rBQJew3aJuMNA18Wy3qCsMdE082y3qCgNdE892i7rCQNdEs92iLjHQNdFst6hLDHRNNNst6hIDXVqA7Ra1jYEuLcB2i9rGQNfEun76obpLkCploGti3b7vibpLkCploGsiTc/MkkP2e0JUbWSgayIVLVf0hKjayEDXRBq2XPHEFZ4QVTsZ6Jo4RVeH3nTpOWOqRKqWga6J49Wh6ioDXRPHq0PVVQa6NMCToWozA10Tpah/brtFbWaga6IU9c+lNjPQNVHsn6vLDHRNjKJ7t9g/V9sZ6JoYRfdusX+utjPQNRG8d4smgYGuiXDdR223qPtKBXpEbIqIRyPiUERcO8/+bRHxcEQ8GBEfi4hXVF+qtDTTM7N89blvLLjfe7eoKwoDPSJOAG4G3gicBVweEWfNGTYDTGXm2cCHgfdVXai0VEVLFb13i7qizAz9POBQZj6emc8BdwKbBwdk5scz82v9p/uA06stU1q6YUsVwdm5uqNMoK8Fnhx4fri/bSFXAX++nKKkqhRdGerJUHXJyipfLCKuBKaAH19g/1ZgK8D69eurfGtpXn6QhSZJmRn6LLBu4Pnp/W3fIiJeD1wHXJyZz873Qpm5IzOnMnNqzZo1S6lXWhTbLZokZQL9PmBjRJwREScBW4BdgwMi4lzgf9EL86erL1NavKIrQ223qGsKAz0zjwJXA3uAR4CdmXkwIm6MiIv7w24CXgzcHREHImLXAi8njc1tBVeG2m5R15TqoWfmbmD3nG03DDx+fcV1SctSNDs/5cQVtlvUOV4pqk4qum/Lb/7M2WOqRBofA12dU3TfFq8MVVcZ6OqcX/7Ig0P3e2WouspAV6dMz8zy7NEXho5xdq6uMtDVKb909wND97tUUV1moKszrrhlL0dfGNY9d6mius1AVydMz8zyyb/74tAxF5y5ynaLOs1AVycUnQgFuP1tPzSGSqT6GOhqvTInQq8835vBqfsMdLVe0YnQFcCvX/L94ylGqpGBrlYrcyJ0+2WuO9dkMNDVWmVOhHpVqCaJga7W2rbzQOEYrwrVJDHQ1UoXbb+Xgk6LyxQ1cQx0tc5F2+/lsae/OnTMClymqMljoKtVyoQ5eCJUk8lAV2tcccveUmG+8bQX2WrRRDLQ1QrXTz9UuKLlmHu2XTjaYqSGKvURdFKdrrhlb+kw/x1bLZpgBroarWzPHHqX99tq0SQz0NVI0zOzvPOu4nXmx1x5/nov79fEM9DVOIuZlUNvvblhLhnoapDFzsqht6LF9eZSj4Gu2i0lyKEX5q5okf6dga7aLGb1ylyGuXQ8A11jdf30Q9y274llvcYFZ66yzSLNw0DXSE3PzHLN3Qd4fvgHCpWyImD7W85xaaK0AANdlakyvOdyVi4VM9BV2vTMLNvuOsAI8npBzsql8gz0DhvljHkcvFhIWpxSgR4Rm4DfBU4A/iAzf2vO/pOBDwE/CHwBuCwzP1dtqe0PKJVjkEtLUxjoEXECcDNwEXAYuC8idmXmwwPDrgK+lJnfExFbgPcCl1VZ6FLXKqsdXIYoLV+ZGfp5wKHMfBwgIu4ENgODgb4ZeHf/8YeB34uIyMyCDwkr76Y9j1b1UmqIk1eu4L1vPtv+uFSRMoG+Fnhy4Plh4LULjcnMoxHxZeA7gH8eHBQRW4GtAOvXr19UoU898/VFjVcz2U6RRmesJ0UzcwewA2BqampRs/fvPvUUZg31VnGpoTReZQJ9Flg38Pz0/rb5xhyOiJXAy+idHK3MNW94lT30hnG2LTVLmUC/D9gYEWfQC+4twH+ZM2YX8LPAXuA/A39VZf8c+Gaf1VUuS+eMWeq2wkDv98SvBvbQW7b4wcw8GBE3Avszcxfwh8AfR8Qh4Iv0Qr9yl5y71hNokrSAUj30zNwN7J6z7YaBx/8KXFptaZKkxVhRdwGSpGoY6JLUEQa6JHWEgS5JHREVry4s/8YRR4B/WOKPr2bOVagTwGOeDB7zZFjOMb8iM9fMt6O2QF+OiNifmVN11zFOHvNk8Jgnw6iO2ZaLJHWEgS5JHdHWQN9RdwE18Jgng8c8GUZyzK3soUuSjtfWGbokaQ4DXZI6ovGBHhEfjIinI+IzA9tWRcQ9EfFY//u311lj1RY45psi4m8j4sGI+GhEnFpnjVWb75gH9r0rIjIiVtdR26gsdMwR8Y7+7/pgRLyvrvpGYYF/2+dExL6IOBAR+yPivDprrFJErIuIj0fEw/3f5y/0t48kwxof6MCtwKY5264FPpaZG4GP9Z93ya0cf8z3AK/OzLOBzwK/Mu6iRuxWjj9mImId8FPAE+MuaAxuZc4xR8RP0PuM3tdk5vcBv11DXaN0K8f/nt8HvCczzwFu6D/viqPAuzLzLOB84O0RcRYjyrDGB3pmfoLePdYHbQb+qP/4j4BLxlrUiM13zJn5l5l5tP90H71PjuqMBX7PAO8H/jvQubP3CxzzzwO/lZnP9sc8PfbCRmiBY07gpf3HLwOeGmtRI5SZn8/MT/cf/wvwCL3PYB5JhjU+0Bfw8sz8fP/xPwIvr7OYGvwc8Od1FzFqEbEZmM3MB+quZYxeCfxoRHwqIv46Iv5T3QWNwTuBmyLiSXp/kXTtr08AImIDcC7wKUaUYW0N9G/qf9Rd52ZvC4mI6+j9GXd73bWMUkT8B+B/0PsTfJKsBFbR+/P8GmBnRES9JY3czwO/mJnrgF+k9wlonRIRLwY+ArwzM78yuK/KDGtroP9TRHwXQP97p/4sXUhEvBX4aeCKqj+ztYHOBM4AHoiIz9FrMX06Ir6z1qpG7zDwp9nz/4AX6N3Iqct+FvjT/uO7gc6cFAWIiBPphfntmXnsOEeSYW0N9GMfSk3/+5/VWMtYRMQmer3kizPza3XXM2qZ+VBmnpaZGzJzA72g+4HM/MeaSxu1aeAnACLilcBJdP9OhE8BP95//DrgsRprqVT/r6s/BB7JzO0Du0aTYZnZ6C/gDuDzwPP0/lNfBXwHvTPDjwH/F1hVd51jOOZDwJPAgf7XB+quc9THPGf/54DVddc5ht/zScBtwGeATwOvq7vOMRzzjwD3Aw/Q6y//YN11Vni8P0KvnfLgwP/dN40qw7z0X5I6oq0tF0nSHAa6JHWEgS5JHWGgS1JHGOiS1BEGuiR1hIEuSR3xbyDZvnQz/eSHAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "plt.scatter(Time, R_data)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 293
        },
        "id": "EmJtZ1aLljYH",
        "outputId": "e12b863b-2f6a-43fe-abdd-c858ee83d008"
      },
      "execution_count": 34,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.collections.PathCollection at 0x7ff40c9e3710>"
            ]
          },
          "metadata": {},
          "execution_count": 34
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEDCAYAAAAcI05xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAATRElEQVR4nO3dfWydd3nG8etK0hQW0rEsTsZos6AqBTEoLXhNWBmDSmYdQ7TSxGBzpk2riIQYovQFxhKtQiISo5Pp/mCaElqVifLSqNRDCAoZq9eB3IDTJE1fGGWsQE1pzEvVwlhLknt/nOPKTY7P8zv2836+H8mK7fOLfT9yc/Xx/fxeHBECANTXiqoLAAD0R1ADQM0R1ABQcwQ1ANQcQQ0ANUdQA0DNFRbUtm+yfcz2fYnj/8T2A7bvt/3JouoCgKZxUfOobb9W0s8k/UtEvCxj7BZJt0q6JCJ+antDRBwrpDAAaJjC7qgj4i5JP1n4Odvn2r7D9kHb/2n7Jd2X3i7poxHx0+7fJaQBoKvsHvUeSe+KiFdJukbSP3U/f56k82x/zfbdti8tuS4AqK1VZX0j28+T9LuS9tme//SZC+rYIul1ks6WdJftl0fE42XVBwB1VVpQq3P3/nhEXNDjtUckHYiIX0r6H9vfUie4v1FifQBQS6W1PiLiCXVC+C2S5I5XdF+eVOduWrbXq9MK+U5ZtQFAnRU5Pe9TkqYlvdj2I7avkDQu6QrbRyTdL+my7vAvSfqx7Qck3Snp2oj4cVG1AUCTFDY9DwCQD1YmAkDNFfIwcf369bF58+YivjQAtNLBgwd/FBEjvV4rJKg3b96smZmZIr40ALSS7e8u9hqtDwCoOYIaAGqOoAaAmiOoAaDmCGoAqLky9/oAgFaaPDSr67/0X/rB47/Qbz7/ubr2D16syy98YW5fn6AGgGWYPDSrq/cd0YmTnVXes4//QlfvOyJJuYU1rQ8AWIarbj38TEjPO3EytPP2o7l9D4IaAJZo6+79OrnIdkk/f/pEbt8nqfVh+2FJT0o6Iel4RIzmVgEANND43mk99uTTpXyvQXrUr4+IHxVWCQA0xOShWX3tv3/Sd4z7vjoYWh8AMKBrug8L+xnftim375ca1CHpy93Tw3f0GmB7h+0Z2zNzc3O5FQgAdTK+d1rHF2tMd21cu1ofvPzluX3P1KB+TUS8UtIfSnqn7deeOiAi9kTEaESMjoz03KkPABotpeUhSQd2juX6fZOCOiJmu38ek3S7pItyrQIAGuCqWw9njrnhrb3O716ezKC2vcb22vn3Jb1B0n25VwIANdZvKt68M1bkt8hloZRZHxsl3W57fvwnI+KO3CsBgJraNXk0aSre9W/J/25aSgjqiPiOpFcU8t0BoAE+cff3Msds37apkLtpiel5ANDX2MRU5piLz12X6yyPUxHUALCI8b3TeujYz/uO2bh2tW55+6sLrYOgBoAeqpqK1wtBDQA9pKw+3J7j6sN+CGoAOEXK6sMVUqF96VO/FwCgK7XlMVHAwpbFENQAsEBqy6OoqXi9ENQA0DU2MZXZ8ih6Kl4vBDUAqLP6MGsqnqTCp+L1QlADgNJXH1aBoAYw9FJWH27ZsKb0lsc8ghrAUEtZfbhC0v6rXldKPYt9fwAYSnWcitcLQQ1gaKVMxbv43HWlTsXrhaAGMJRSVx9WMcujVx0AMFSa0vKYR1ADGDp1XH3YD0ENYKjUdfVhPwQ1gKFR59WH/RDUAIZGnVcf9kNQAxgKdV992A9BDaD1mrD6sB+CGkCrNW0qXi8ENYBWu+rWw5lj6rD6sB+CGkBrjU1MKWMmXm1WH/ZDUANopclDs0lT8erc8phHUANopZSWR51WH/ZDUANona2792e2POq2+rAfghpAq+yaPKrHnnw6c1zd+9ILEdQAWiVl9eENDehLL5Qc1LZX2j5k+/NFFgQAS7V19/7MMXWfitfLIHfU75b0YFGFAMByjO+dzmx5bFy7ulEtj3lJQW37bEl/JOljxZYDAIPbNXk0afXhgZ1jJVSTv9Q76hskvVfSycUG2N5he8b2zNzcXC7FAUCWyUOzjd0VL1VmUNt+k6RjEXGw37iI2BMRoxExOjIykluBANDP+267N3PMxrWrGzMVr5eUO+qLJb3Z9sOSPi3pEtufKLQqAEiwa/Konjq+6C/6z2hqy2NeZlBHxPsj4uyI2CzpbZL+PSK2F14ZAPSR2vJo2lS8XphHDaCRUg6obeJUvF5WDTI4IqYkTRVSCQAkGt87nXlAbRN2xUvFHTWARmnDQQCDIqgBNMowtTzmEdQAGmNsYiqz5bFlw5rWtDzmEdQAGiHlgFqpvgfULgdBDaD2UvvSTV592A9BDaD2UvrSWzasafTqw34IagC1ljoVr40tj3kENYDaGsapeL0Q1ABqq00H1C4HQQ2glsYmpjIPqD1jhVrbl16IoAZQO7smjyZNxbv+Le1uecwjqAHUyiAHAbS95TGPoAZQK6lLxIeh5TGPoAZQGylLxKX27IqXiqAGUAupfek2HAQwKIIaQC2k9KXbtiteKoIaQOXGJqYyx7RxV7xUBDWASqXsitf2JeJZCGoAldk1eZQl4gkIagCVSJ0vPax96YUIagCVeN9t92aOadMBtctBUAMo3a7Jo3rq+MnMccPe8phHUAMoFUvEB0dQAyhVytalw7ZEPAtBDaA0KVuX0pc+HUENoBSTh2aTlojTlz4dQQ2gFKktD/rSpyOoARRu6+79mS2PYV4inoWgBlCosYkpPfbk05njhnmJeBaCGkBh2Lo0H5lBbfs5tr9u+4jt+21/oIzCADQfS8TzsSphzFOSLomIn9k+Q9JXbX8xIu4uuDYADZaydenGtavpSyfIDOqICEk/6354Rvct+6wcAEMrZetSSTqwc6yEapovqUdte6Xtw5KOSdofEQd6jNlhe8b2zNzcXN51AmiIyUOzSVuX0pdOlxTUEXEiIi6QdLaki2y/rMeYPRExGhGjIyMjedcJoCGYL52/gWZ9RMTjku6UdGkx5QBosvG90ywRL0DKrI8R28/vvv9cSWOSvll0YQCahdNaipMy6+MFkj5ue6U6wX5rRHy+2LIANAlblxYrZdbHvZIuLKEWAA2V0pfesmENW5cuESsTASxLytalEkvEl4OgBrBkLBEvB0ENYEk4Rbw8BDWAJblm35HMMWxdmg+CGsDAxvdO63hCY5q+dD4IagADYYl4+QhqAANhiXj5CGoAyThSqxoENYAkHKlVHYIaQKbU+dLbt20qoZrhQ1ADyJQyX5ol4sUhqAH0tXX3/swxG9eupuVRIIIawKJS+9IcqVUsghpAT+zjUR8ENYCe2MejPghqAKcZm5jKHLNx7WrmS5eEoAbwLON7p5NaHvSly0NQA3hG6rmH9KXLRVADkMT+0nVGUAOQlLa/9AqJvnQFCGoAyftLT9DyqARBDQy51L709m2baHlUhKAGhtggfWn28agOQQ0MMfrSzUBQA0OKvnRzENTAEKIv3SwENTBk6Es3D0ENDJmUw2npS9cLQQ0MkbGJqczDaSX60nWTGdS2z7F9p+0HbN9v+91lFAYgX4Oce0hful5WJYw5LunqiLjH9lpJB23vj4gHCq4NQE7oSzdb5h11RDwaEfd0339S0oOS+N8t0CApfektG9bQl66pgXrUtjdLulDSgSKKAZC/rbv3J/WlOZy2vpKD2vbzJN0m6cqIeKLH6ztsz9iemZuby7NGAEs0vnc66XBa9peut6Sgtn2GOiF9S0R8tteYiNgTEaMRMToyMpJnjQCWYPLQbNKiFvaXrr+UWR+WdKOkByNioviSAOSBvnR7pNxRXyzpzyVdYvtw9+2NBdcFYBnOv+6OzL70xrWr6Us3ROb0vIj4qiSXUAuAHIxNTOmJp05kjuNw2uZgZSLQIoMsakFzENRAS6QuatmyYQ2LWhqGoAZaIuXhIX3pZiKogRZIXdRCX7qZCGqg4VjU0n4ENdBgLGoZDgQ10GAsahkOBDXQUCxqGR4ENdBALGoZLgQ10DAsahk+BDXQICxqGU4ENdAgLGoZTgQ10BApDw8l+tJtRFADDZD68JBFLe1EUAM1l/rwkEUt7UVQAzU2yMNDFrW0F0EN1BgPDyER1EBtsSMe5hHUQA2xIx4WIqiBmtk1eZQd8fAsBDVQI7smj/LwEKchqIGaSJ3hIYmHh0OGoAZq4pp9R5LG0ZcePgQ1UANjE1M6njDFY/u2TfSlhxBBDVRsfO908spDdsQbTgQ1UKHUMw95eDjcCGqgQu/5TPbKwxXi4eGwI6iBipx/3R1KWHioCR4eDj2CGqjA1t37k7Yt5eEhJIIaKN3YxFTS8vDt2zbx8BCSEoLa9k22j9m+r4yCgDZLneHBmYdYKOWO+mZJlxZcB9B6qXt4nHXmSh4e4lkygzoi7pKU/V8XgEWlLg8/68yVuvcD3Bfh2XLrUdveYXvG9szc3FxeXxZohZQDACQR0ugpt6COiD0RMRoRoyMjI3l9WaDxUg8AYA8PLIZZH0CBUmd4sLc0+iGogYKMTUwlz/BgeTj6SZme9ylJ05JebPsR21cUXxbQbIOENDM8kGVV1oCI+NMyCgHaInWutMQeHkhD6wPIUepcaYmHh0hHUAM5ST3vUGIPDwyGoAZyMMh5h+zhgUER1EAOUvaVlghpLA1BDSxT6r7SHKWFpSKogWU4/7o7kvaVZq40liNzeh6A3l6y8wv6vxPZ99LMlcZycUcNLEFqSG9cu5qQxrIR1MCAzr/ujqSQPuvMlTqwc6yEitB2tD6ARJOHZnVl4uwO9pVGnrijBhIQ0qgSQQ1kGCSkn7PShDRyR+sD6GOQZeHPWWl9c/cbC64Iw4g7amARg4T0WWeuJKRRGO6ogR7G904n74JHTxpFI6iBU2zdvT/p+CyJkEY5CGqga5CHhhIhjfIQ1IAGa3VInRWHLGZBWQhqDL3UjZXmEdIoG0GNoTXIrI55hDSqQFBj6Azai5538bnr2KoUlSCoMVTGJqaSTwhf6Ia3XsAZh6gMQY2hsNSAZmYH6oCgRmtNHprVtfsO65cnl/b3aXWgLghqtM5SHhKeilYH6oSgRivkEc4Sx2ahnghqNNLkoVld9ZnDWmJX4zSW9BHuolFTBDVqb7m95iz0olF3BDVqYdAl3HkgoNEUBDWWrIpwzQMBjaZJCmrbl0r6R0krJX0sIj6UdyFN/UeP5ti+bZM+ePnLqy4DGFhmUNteKemjksYkPSLpG7Y/FxEP5FUEIY2inLlqhf7+j8/nISEaLeWO+iJJ346I70iS7U9LukxSbkFNSCNPhDPaJiWoXyjp+ws+fkTS1lMH2d4haYckbdq0KZfigFS0NdBmuT1MjIg9kvZI0ujoaOT1dYFTEcoYNilBPSvpnAUfn939XG4uPncd7Q88g9WBwLOlBPU3JG2x/SJ1Avptkv4szyJuefureaDYcPSFgeJkBnVEHLf915K+pM70vJsi4v68C2FeKwD0ltSjjogvSPpCwbUAAHpYUXUBAID+CGoAqDmCGgBqjqAGgJpzRP5rU2zPSfruEv/6ekk/yrGcJuCahwPX3H7Lud7fioiRXi8UEtTLYXsmIkarrqNMXPNw4Jrbr6jrpfUBADVHUANAzdUxqPdUXUAFuObhwDW3XyHXW7seNQDg2ep4Rw0AWICgBoCaqzSobd9k+5jt+xZ8bp3t/bYf6v75a1XWmKdFrvd629+0fa/t220/v8oa89brmhe8drXtsL2+itqKstg1235X92d9v+0PV1VfERb5b/sC23fbPmx7xvZFVdaYN9vn2L7T9gPdn+m7u5/PPcOqvqO+WdKlp3zubyR9JSK2SPpK9+O2uFmnX+9+SS+LiPMlfUvS+8suqmA36/Rrlu1zJL1B0vfKLqgEN+uUa7b9enXOGn1FRPy2pH+ooK4i3azTf84flvSBiLhA0t91P26T45KujoiXStom6Z22X6oCMqzSoI6IuySdelrAZZI+3n3/45IuL7WoAvW63oj4ckQc7354tzon6LTGIj9jSfqIpPdKat3T7EWu+R2SPhQRT3XHHCu9sAItcs0h6azu+78q6QelFlWwiHg0Iu7pvv+kpAfVOWM29wyr+o66l40R8Wj3/R9K2lhlMSX7K0lfrLqIotm+TNJsRBypupYSnSfp92wfsP0ftn+n6oJKcKWk621/X53fINr22+IzbG+WdKGkAyogw+oY1M+IztzB1t1x9WJ7pzq/St1SdS1Fsv0rkv5WnV+Fh8kqSevU+RX5Wkm32na1JRXuHZLeExHnSHqPpBsrrqcQtp8n6TZJV0bEEwtfyyvD6hjUj9l+gSR1/2zVr4i92P5LSW+SNB7tn9h+rqQXSTpi+2F1Wj332P6NSqsq3iOSPhsdX5d0Up0NfNrsLyR9tvv+PkmtepgoSbbPUCekb4mI+WvNPcPqGNSfU+cHrO6f/1phLYWzfak6vdo3R8T/Vl1P0SLiaERsiIjNEbFZnQB7ZUT8sOLSijYp6fWSZPs8SavV/l3lfiDp97vvXyLpoQpryV33N6IbJT0YERMLXso/wyKisjdJn5L0qKRfqvMP9gpJv67Ok9KHJP2bpHVV1ljC9X5b0vclHe6+/XPVdRZ9zae8/rCk9VXXWcLPebWkT0i6T9I9ki6pus4Srvk1kg5KOqJO7/ZVVdeZ8zW/Rp22xr0L/v2+sYgMYwk5ANRcHVsfAIAFCGoAqDmCGgBqjqAGgJojqAGg5ghqAKg5ghoAau7/ASwqmADs4MP+AAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "plt.scatter(Time, N_data)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 293
        },
        "id": "uN90E-c4ljkh",
        "outputId": "7d10ea76-834d-425f-d9ca-fff71cd1ff02"
      },
      "execution_count": 35,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.collections.PathCollection at 0x7ff40c9c4dd0>"
            ]
          },
          "metadata": {},
          "execution_count": 35
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEDCAYAAAAlRP8qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAWK0lEQVR4nO3df7DddZ3f8ec74ceuCItshN2F3L0ME2wVEXevBIqt0Zo1utbQ6cjvTqlIpo44q1i6sVgiVmdc2QLOlF16dbPRSsHVpXczBcF0uywdC2mCgWhiBcRduBeXuCoLVQcJvPvHObHHm3tzz/me7/d8z/me52Mmk3vO9/O95/2dm7zmfb/fz+f7jcxEkjT6ltVdgCSpHAa6JDWEgS5JDWGgS1JDGOiS1BAGuiQ1RK2BHhGbI2JfRHyji7E3RMSD7T8PR8TTg6hRkkZF1DkPPSL+EfB/gc9l5mk97Pc+4LWZ+a7KipOkEVNrh56Z9wI/6HwvIk6JiLsi4oGI+J8R8fcW2PVC4NaBFClJI+KwugtYwDTwrzLzkYhYDfwB8KYDGyPi14GTgf9RU32SNJSGKtAj4qXAPwC+GBEH3j5y3rALgC9l5guDrE2Sht1QBTqtU0BPZ+YZhxhzAfDeAdUjSSNjqKYtZuYzwHci4p0A0fKaA9vb59NfBtxXU4mSNLTqnrZ4K61wfkVEzEbEZcDFwGUR8RCwB1jfscsFwG3pLSIl6SC1TluUJJVnqE65SJKKq+2i6IoVK3JycrKuj5ekkfTAAw/8bWa+fKFttQX65OQkO3furOvjJWkkRcRfL7bNUy6S1BAGuiQ1hIEuSQ2xZKB3c4vbiFjTvq3tnoj4y3JLlCR1o5sOfQuwbrGNEXEsrRtovSMzXwW8s5zSJEm9WHKWS2beGxGThxhyEXB7Zj7eHr+vnNIkqTlmds3xka17ePonzwPwspcczqZ/8irOfe2JpX1GGdMWTwUOj4h7gKOBT2Xm5xYaGBEbgA0AExMTJXy0JA2/iz99H1/99s89+oEf/vh5rvrSQwClhXoZF0UPA34T+G3gLcC/i4hTFxqYmdOZOZWZUy9/+YLz4iWpMWZ2zTG58Y6DwvyA519Irrv7W6V9Xhkd+izw/cz8EfCjiLgXeA3wcAnfW5JG0kJd+UKefPonpX1mGR36nwGvj4jDIuIlwGrgmyV8X0kaOUt15fP92rG/WNpnL9mht29xuwZYERGzwCbgcIDMvDkzvxkRdwG7gReBz2TmolMcJamp1l5/D4/s+1FP+1z1lleU9vndzHK5sIsx1wHXlVKRJI2YmV1zvP8LD/a836rjjxq6WS6SNLaKdOXQCvNtV64ptRYDXZIKKNqVB3DD+WeU2pkfYKBLUg9mds1x5Rce5MUC+1bRlXcy0CWpS91ORZxvWcD151XTlXcy0CVpCUVPrwCcc8px3HL52SVXtDADXZIOYdi78k4GuiQtYFS68k4GuiTNU3QqYh1deScDXZLa+unKLzlrgo+d++qSK+qNgS5JDNcCoaIMdEljbRgXCBVloEsaW03oyjsZ6JLGTtGuvO6Lnksx0CWNlaJdeV1TEXthoEsaCx+e+Tqfv//xnvcb9q68k4EuqdFGfSpiLwx0SY1V9PTKMUcuZ/e16yqoqFoGuqTGGaeuvFM3zxTdDLwd2JeZpx1i3OuA+4ALMvNL5ZUoSd0bt66807IuxmwBDnmUEbEc+D3gKyXUJEk9m9k1x+TGOwqF+SVnTYx8mEN3D4m+NyImlxj2PuBPgdeVUJMk9aRoV37C0Uew/eq1FVRUj77PoUfEicA/Bd6IgS5pgMb1XPliyrgoeiPwu5n5YkQccmBEbAA2AExMTJTw0ZLGVZMXCBVVRqBPAbe1w3wF8LaI2J+ZM/MHZuY0MA0wNTWVJXy2pDEzDguEiuo70DPz5ANfR8QW4L8tFOaS1K/VH9/GU8/+tOf9mnh6ZSHdTFu8FVgDrIiIWWATcDhAZt5caXWSRPGuvAlTEXvRzSyXC7v9Zpl5aV/VSNI8duXdc6WopKFkV947A13SUHEqYnEGuqSh4QKh/hjokmpnV14OA11SrezKy2OgS6pF0YueYFe+GANd0sAVnYq46vij2HblmvILaggDXdLAFO3KA7jh/GYv2y+DgS5pIE7fdBfPPPdCz/t5eqV7BrqkSrlAaHAMdEmVsSsfLANdUumKduVOReyPgS6pNC4QqpeBLqkUF3/6Pr767R/0vJ9TEctjoEvqSz9d+Y1ORSyVgS6psKLL9u3Kq2GgS+pZP8v27cqrY6BL6knRqYh25dUz0CV1xWX7w6+bh0RvBt4O7MvM0xbYfjHwu7R+bs8C78nMh8ouVFJ97MpHQzcd+hbgPwKfW2T7d4A3ZOYPI+KtwDSwupzyJNXJrny0LBnomXlvREweYvv/6nh5P3BS/2VJqlM/UxHPOeU4brn87JIrUjfKPod+GfDlxTZGxAZgA8DExETJHy2pDC4QGl2lBXpEvJFWoL9+sTGZOU3rlAxTU1NZ1mdL6p8LhEZfKYEeEacDnwHempnfL+N7Shocu/Jm6DvQI2ICuB3455n5cP8lSRoUu/Jm6Wba4q3AGmBFRMwCm4DDATLzZuAa4JeBP4gIgP2ZOVVVwZLK4bL95ulmlsuFS2x/N/Du0iqSVKmiXblTEYefK0WlMWJX3mwGujQG7MrHg4EuNdjMrjmu/MKDvFhgX7vy0WOgSw1VdCrisoDrz7MrH0UGutQwLtsfXwa61CBFL3ralTeDgS41gF25wECXRp5duQ4w0KUR1U9XfslZE3zs3FeXXJHqZqBLI8gFQlqIgS6NkKJduadXxoOBLo2Iol25Fz3Hh4EuDTm7cnXLQJeGmF25emGgS0PowzNf5/P3P97zfnbl481Al4bM6o9v46lnf9rzfk5FlIEuDYmiXfkxRy5n97XrKqhIo8ZAl4aAXbnK0M0zRTcDbwf2ZeZpC2wP4FPA24AfA5dm5tfKLlRqoqK3uLUr10KWdTFmC3CofzlvBVa1/2wA/rD/sqRmm9k1x+TGOwqF+SVnTRjmWlA3D4m+NyImDzFkPfC5zEzg/og4NiJ+NTO/W1KNUqMUnYp4wtFHsP3qtRVUpKYo4xz6icATHa9n2+8Z6FIHb6alqg30omhEbKB1WoaJiYlBfrRUK7tyDUIZgT4HrOx4fVL7vYNk5jQwDTA1NZUlfLY01IpORQS7cvWujEDfClwREbcBq4G/8/y5VHwqosv2VVQ30xZvBdYAKyJiFtgEHA6QmTcDd9KasvgorWmL/7KqYqVR4AIh1aWbWS4XLrE9gfeWVpE0wlwgpDq5UlQqgV25hoGBLvXp9E138cxzL/S8n125ymagSwUV7cqdiqiqGOhSAXblGkYGutQDu3INMwNd6oLL9jUKDHRpCUVvcbvq+KPYduWa8guSFmGgS4vopyu/8Xyf66nBM9ClBRS9mZZduepkoEsd+rmZll256magS21FpyLalWtYGOgae0W78gBusCvXEDHQNbb6uehpV65hZKBrLBWdirgs4Prz7Mo1nAx0jZV+unIfPKFhZ6BrbBTtyg9bFvz+O19jV66hZ6Cr8VwgpHFhoKvRXCCkcWKgq5GKduVORdQo6yrQI2Id8ClgOfCZzPzEvO0TwGeBY9tjNmbmnSXXKnXFrlzjaslAj4jlwE3AWmAW2BERWzNzb8ewDwN/kpl/GBGvBO4EJiuoV1qUXbnGXTcd+pnAo5n5GEBE3AasBzoDPYFj2l//EvBkmUVKS7Erl7oL9BOBJzpezwKr5435CPCViHgfcBTw5oW+UURsADYATExM9FqrdJCiXbkLhNREZV0UvRDYkpn/ISLOBv5zRJyWmS92DsrMaWAaYGpqKkv6bI2pol25C4TUVN0E+hywsuP1Se33Ol0GrAPIzPsi4heAFcC+MoqUOtmVSwvrJtB3AKsi4mRaQX4BcNG8MY8D/xjYEhF/H/gF4HtlFiqBXbl0KEsGembuj4grgLtpTUncnJl7IuKjwM7M3Ap8EPh0RHyA1gXSSzPTUyoqTdFb3NqVa5x0dQ69Paf8znnvXdPx9V7gnHJLk1pWf3wbTz370573u+SsCT527qsrqEgaTq4U1dAq2pUfc+Rydl+7roKKpOFmoGvo9HMzLbtyjTMDXUOl6EVPu3LJQNeQsCuX+megq3ZFu/ITjj6C7VevraAiaTQZ6KpN0YueYFcuLcRAVy2KTkW0K5cWZ6BroOzKpeoY6BoYFwhJ1TLQVTkXCEmDYaCrUqdvuotnnnuh5/3syqXeGeiqRNGu3IueUnEGukpnVy7Vw0BXaezKpXoZ6Oqby/al4WCgqy8Xf/o+vvrtH/S836rjj2LblWvKL0gaYwa6CumnK7/xfJ8gJFXBQFfP7Mql4bSsm0ERsS4ivhURj0bExkXGnBcReyNiT0T8l3LL1DCY2TXH5MY7CoX5jeefYZhLFVuyQ4+I5cBNwFpgFtgREVvbzxE9MGYV8CHgnMz8YUQcX1XBqkfRW9zalUuD080plzOBRzPzMYCIuA1YD+ztGHM5cFNm/hAgM/eVXajqUfRceQA3eK5cGqhuAv1E4ImO17PA6nljTgWIiK8Cy4GPZOZdpVSo2tiVS6OlrIuihwGrgDXAScC9EfHqzHy6c1BEbAA2AExMTJT00SqbXbk0mroJ9DlgZcfrk9rvdZoFtmfm88B3IuJhWgG/o3NQZk4D0wBTU1NZtGhVY2bXHFd+4UFeLLCvXblUv24CfQewKiJOphXkFwAXzRszA1wI/HFErKB1CuaxMgtVtYpORVwWcP15duXSMFgy0DNzf0RcAdxN6/z45szcExEfBXZm5tb2tt+KiL3AC8BVmfn9KgtXOfpZIHTOKcdxy+Vnl1yRpKIis54zH1NTU7lz585aPlstduXS6ImIBzJzaqFtrhQdQ3blUjMZ6GOm6FREu3Jp+BnoY8Jb3ErNZ6CPARcISePBQG8wFwhJ48VAbyi7cmn8GOgNU7Qr96KnNPoM9AYp2pU7FVFqBgO9AT4883U+f//jPe9nVy41i4E+4lZ/fBtPPfvTnvdzKqLUPAb6iCralR9z5HJ2X7uugook1c1AH0F25ZIWYqCPELtySYdioI+I0zfdxTPPvdDzfnbl0vgw0Idc0a78hKOPYPvVayuoSNKwMtCHmF25pF4Y6EPIrlxSEQb6EPEWt5L6YaAPiaKPg/NmWpIOWNbNoIhYFxHfiohHI2LjIcb9s4jIiFjweXc62MyuOSY33lEozG88/wzDXNLPLNmhR8Ry4CZgLTAL7IiIrZm5d964o4HfAbZXUWgT2ZVLKlM3p1zOBB7NzMcAIuI2YD2wd964fw/8HnBVqRU2UD/nym/0wROSFtFNoJ8IPNHxehZY3TkgIn4DWJmZd0TEooEeERuADQATExO9V9sARZft25VLWkrfF0UjYhlwPXDpUmMzcxqYBpiamsp+P3uUFJ2K6OPgJHWrm0CfA1Z2vD6p/d4BRwOnAfdEBMCvAFsj4h2ZubOsQkdZ0QVCduWSetFNoO8AVkXEybSC/ALgogMbM/PvgBUHXkfEPcC/NsztyiUN1pKBnpn7I+IK4G5gObA5M/dExEeBnZm5teoiR00/Fz3tyiUV1dU59My8E7hz3nvXLDJ2Tf9lja6iUxF9HJykfrlStCQu25dUNwO9BC4QkjQMDPQ+uEBI0jAx0Atae/09PLLvRz3vZ1cuqSoGeo+KduVORZRUNQO9B3blkoaZgd4Fu3JJo8BAX4JduaRRYaAvomhX7gIhSXUx0BdQtCs/55TjuOXysyuoSJKWZqB3KHozLbtyScPAQG8r+uAJl+1LGhZjH+hFu/JjjlzO7mvXVVCRJBUz1oFuVy6pScYy0O3KJTXR2AV60cfB2ZVLGnZjE+hFu/ITjj6C7VevraAiSSrXWAS6XbmkcdBVoEfEOuBTtJ4p+pnM/MS87VcC7wb2A98D3pWZf11yrT2zK5c0TpYM9IhYDtwErAVmgR0RsTUz93YM2wVMZeaPI+I9wCeB86souBs+Dk7SOOqmQz8TeDQzHwOIiNuA9cDPAj0z/6Jj/P3AJWUW2QsfBydpXHUT6CcCT3S8ngVWH2L8ZcCX+ymqCB8HJ2nclXpRNCIuAaaANyyyfQOwAWBiYqK0z7Url6TuAn0OWNnx+qT2ez8nIt4MXA28ITOfW+gbZeY0MA0wNTWVPVc7j125JP1/3QT6DmBVRJxMK8gvAC7qHBARrwX+E7AuM/eVXuUCfPCEJP28JQM9M/dHxBXA3bSmLW7OzD0R8VFgZ2ZuBa4DXgp8MSIAHs/Md1RRsI+Dk6SFdXUOPTPvBO6c9941HV+/ueS6FlR0XrlduaRxMDIrRWd2zfUc5nblksbJyAT6dXd/q6fxduWSxs3IBPqTT/+kq3E+Dk7SuBqZQP+1Y3+RuSVC3Yc0Sxpny+ouoFtXveUVixa7LFrzyg1zSeNsZDr0A6dQPnT7bn7y/ItAK8gvWu3NtCQJRijQoRXqnhuXpIWNzCkXSdKhGeiS1BAGuiQ1hIEuSQ1hoEtSQ0Rm37clL/bBEd8Dij5IegXwtyWWMwo85vHgMY+Hfo751zPz5QttqC3Q+xEROzNzqu46BsljHg8e83io6pg95SJJDWGgS1JDjGqgT9ddQA085vHgMY+HSo55JM+hS5IONqoduiRpHgNdkhpi6AM9IjZHxL6I+EbHe8dFxLaIeKT998vqrLFsixzzdRHxfyJid0T814g4ts4ay7bQMXds+2BEZESsqKO2qix2zBHxvvbPek9EfLKu+qqwyL/tMyLi/oh4MCJ2RsSZddZYtohYGRF/ERF72z/T32m/X3qODX2gA1uAdfPe2wj8eWauAv68/bpJtnDwMW8DTsvM04GHgQ8NuqiKbeHgYyYiVgK/BfT2hPDRsIV5xxwRbwTWA6/JzFcBv19DXVXawsE/508C12bmGcA17ddNsh/4YGa+EjgLeG9EvJIKcmzoAz0z7wV+MO/t9cBn219/Fjh3oEVVbKFjzsyvZOb+9sv7gZMGXliFFvk5A9wA/BugcVfvFznm9wCfyMzn2mP2DbywCi1yzAkc0/76l4AnB1pUxTLzu5n5tfbXzwLfBE6kghwb+kBfxAmZ+d32138DnFBnMTV4F/DluouoWkSsB+Yy86G6axmgU4F/GBHbI+IvI+J1dRc0AO8HrouIJ2j9RtK03z5/JiImgdcC26kgx0Y10H8mW/MuG9e9LSYirqb1K9wtdddSpYh4CfBvaf0KPk4OA46j9av5VcCfRETUW1Ll3gN8IDNXAh8A/qjmeioRES8F/hR4f2Y+07mtrBwb1UB/KiJ+FaD9d6N+LV1MRFwKvB24OJu/gOAU4GTgoYj4K1qnmL4WEb9Sa1XVmwVuz5b/DbxI60ZOTfYvgNvbX38RaNRFUYCIOJxWmN+SmQeOtfQcG9VA30rrHwHtv/+sxloGIiLW0TqX/I7M/HHd9VQtM7+emcdn5mRmTtIKut/IzL+pubSqzQBvBIiIU4EjaP6dCJ8E3tD++k3AIzXWUrr2b1h/BHwzM6/v2FR+jmXmUP8BbgW+CzxP6z/1ZcAv07oq/Ajw34Hj6q5zAMf8KPAE8GD7z81111n1Mc/b/lfAirrrHMDP+Qjg88A3gK8Bb6q7zgEc8+uBB4CHaJ1b/s266yz5mF9P63TK7o7/v2+rIsdc+i9JDTGqp1wkSfMY6JLUEAa6JDWEgS5JDWGgS1JDGOiS1BAGuiQ1xP8DGgmlQkBInMsAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    }
  ]
}
About this Algorithm

Runga Kutta methods are used to approximate the value to simultaneous non-linear differential equations.
More about Runga Kutta Methods here.

Here we will try to solve the following three equations with unknown variables S, I and R.

\begin{equation} \frac{d S(t)}{d t}=A-\mu S(t)-\frac{r I}{N_0} S(t) \end{equation}
$$ \frac{d I(t)}{d t}=\frac{r I(t)}{N_0} S(t)-(\mu+a) I(t) $$
$$ \frac{d R(t)}{d t}=a I(t)-\mu R(t) $$
          where, $\mathrm{S}(\mathrm{t})+\mathrm{I}(\mathrm{t})+\mathrm{R}(\mathrm{t})=\mathrm{N}(\mathrm{t})$

$$ =>\frac{N(t)}{d t}=A-\mu N $$

These equations have a physical significance also, they represent the SIR model for the transmission of tuberclosis as documented here.

S, I and R represent the following three kinds of people:

  1. S: The first category is individuals who have not been infected with the TB virus, sometimes known as the suspicion group.
  2. I: The second category consists of people who have been infected with the TB virus.
  3. R: The last category consists of people who recovered or died after being infected with the TB virus.

Other than S, I and R, all other terms are constants, and the equations can be represented as follows:

def dS_dt(S, I, R, t):
    return (1449401 -  0.001167*S - (0.5/1449401)*I*S)
def dI_dt(S, I, R, t):
    return ((0.5/1449401)*I*S - (0.001167 + 0.111111)*I)
def dR_dt(S, I, R, t):
    return  (0.111111*I - 0.001167*R)
t = int(input())  # Enter a required value of time where you want the value of S, I and R.
10
S_data = []
I_data = []
R_data = []
SIR_data = []
N_data = []
Time = []

Runge-Kutta method implementation with S,I,R as dependent variables and t as independent variable

def rungeKutta(S0, I0, R0, t, h):

  # No of iterations = n
  # Step size = h
    n = (int)((t)/h)
  
  # setting initial values
    S = S0
    I = I0
    R = R0

  # implementation of runge-kutta
    for i in range(1, n + 1):
        kS1 = dS_dt(S, I, R, t)
        kI1 = dI_dt(S, I, R, t)
        kR1 = dR_dt(S, I, R, t)
        
        kS2 = dS_dt(S + 0.5 * h * kS1, I + 0.5 * h * kI1, R + 0.5 * h * kR1, t + 0.5*h)
        kI2 = dI_dt(S + 0.5 * h * kS1, I + 0.5 * h * kI1, R + 0.5 * h * kR1, t + 0.5*h)
        kR2 = dR_dt(S + 0.5 * h * kS1, I + 0.5 * h * kI1, R + 0.5 * h * kR1, t + 0.5*h)

        kS3 = dS_dt(S + 0.5 * h * kS2, I + 0.5 * h * kI2, R + 0.5 * h * kR2, t + 0.5*h)
        kI3 = dI_dt(S + 0.5 * h * kS2, I + 0.5 * h * kI2, R + 0.5 * h * kR2, t + 0.5*h)
        kR3 = dR_dt(S + 0.5 * h * kS2, I + 0.5 * h * kI2, R + 0.5 * h * kR2, t + 0.5*h)

        kS4 = dS_dt(S + kS3 * h, I + kI3 * h, R + kR3 * h, t + h)
        kI4 = dI_dt(S + kS3 * h, I + kI3 * h, R + kR3 * h, t + h)
        kR4 = dR_dt(S + kS3 * h, I + kI3 * h, R + kR3 * h, t + h)
 
        # Updating S, I, R
        S = S + (1.0 / 6.0)*(kS1 + 2 * kS2 + 2 * kS3 + kS4)*h
        I = I + (1.0 / 6.0)*(kI1 + 2 * kI2 + 2 * kI3 + kI4)*h
        R = R + (1.0 / 6.0)*(kR1 + 2 * kR2 + 2 * kR3 + kR4)*h

        S_data.append(S)
        I_data.append(I)
        R_data.append(R)
        SIR_data.append(S+I+R)
        Time.append(t)

        # Updating value of t
        t += h
    
    # printing N(t) at desired t (N(t)=S(t)+I(t)+R(t))
    print("N(t) = " ,S+I+R)
    return [S,I,R]
# Driver Code
S0 = 1446093
I0 = 1885
R0 = 1423
h = 0.01
print ('The value of S(t), I(t) and R(t) is:', rungeKutta(S0, I0, R0, t, h)) 
N(t) =  15842350.284962129
The value of S(t), I(t) and R(t) is: [410394.8072337986, 10323985.625711355, 5107969.852016973]

Runge-Kutta method implementation with N as dependent variable and t as independent variable

We implemented it just to prove that the equation N(t) = S(t)+I(t)+R(t) is consistent with the values of S,I,R we got above for the same value of t.

def dN_dt(t, N):
	return (1449401-0.001167*N)

# implementing runge kutta method for ODE: dN/dt = A - m*N where m is death rate
def rungeKutta(t0, A, t, h):

	# no. of iterations = n
	# step size = h
	n = (int)((t - t0)/h)
	
	N = A
	for i in range(1, n + 1):
		k1 = dN_dt(t0, N)
		k2 = dN_dt(t0 + 0.5 * h, N + 0.5 * h * k1)
		k3 = dN_dt(t0 + 0.5 * h, N + 0.5 * h * k2)
		k4 = dN_dt(t0 + h, N + h * k3)

		# Update value of y
		N = N + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4)*h
		N_data.append(N)

		# Update next value of t0
		t0 += h
	return N

# Driver Code
t0 = 0
A = 1449401
h = 0.01
print ('The value of N(t) at required t is:', rungeKutta(t0, A, t, h))
The value of N(t) at required t is: 15842350.284962097
# print(S_data)
# print(I_data)
# print(R_data)
# print(Time)
# print(N_data)

Variation of S, I, R and N with time.

import math
import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(20,38, 1000)
plt.scatter(Time, S_data)
&lt;matplotlib.collections.PathCollection at 0x7ff40cbf3590&gt;
plt.scatter(Time, I_data)
&lt;matplotlib.collections.PathCollection at 0x7ff40cbf3e10&gt;
plt.scatter(Time, R_data)
&lt;matplotlib.collections.PathCollection at 0x7ff40c9e3710&gt;
plt.scatter(Time, N_data)
&lt;matplotlib.collections.PathCollection at 0x7ff40c9c4dd0&gt;