#### Runge Kutta Methods

H
{
"nbformat": 4,
"nbformat_minor": 0,
"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",
"\\n",
"\\frac{d S(t)}{d t}=A-\\mu S(t)-\\frac{r I}{N_0} S(t)\n",
"\\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:"
],
"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)"
],
"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)"
],
"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)"
],
"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."
],
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "snFbzheU5k8S",
},
"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 = []"
],
"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**"
],
"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]"
],
"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)) "
],
"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**"
],
"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."
],
"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))"
],
"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)"
],
"id": "qrXWMyh7Wy03"
},
"execution_count": 27,
"outputs": []
},
{
"cell_type": "code",
"source": [
"# print(I_data)"
],
"id": "zSjyo5IDkMeT"
},
"execution_count": 28,
"outputs": []
},
{
"cell_type": "code",
"source": [
"# print(R_data)"
],
"id": "J8oCbC9FkNr8"
},
"execution_count": 29,
"outputs": []
},
{
"cell_type": "code",
"source": [
"# print(Time)"
],
"id": "rGsiEb3zkUpd"
},
"execution_count": 30,
"outputs": []
},
{
"cell_type": "code",
"source": [
"# print(N_data)"
],
"id": "59T-ht2x-9h5"
},
"execution_count": 31,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"**Variation of S, I, R and N with time.**"
],
"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)"
],
"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>"
]
},
"execution_count": 32
},
{
"output_type": "display_data",
"data": {
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
],
},
"needs_background": "light"
}
}
]
},
{
"cell_type": "code",
"source": [
"plt.scatter(Time, I_data)"
],
"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>"
]
},
"execution_count": 33
},
{
"output_type": "display_data",
"data": {
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
],
},
"needs_background": "light"
}
}
]
},
{
"cell_type": "code",
"source": [
"plt.scatter(Time, R_data)"
],
"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>"
]
},
"execution_count": 34
},
{
"output_type": "display_data",
"data": {
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
],
},
"needs_background": "light"
}
}
]
},
{
"cell_type": "code",
"source": [
"plt.scatter(Time, N_data)"
],
"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>"
]
},
"execution_count": 35
},
{
"output_type": "display_data",
"data": {
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
],
},
"needs_background": "light"
}
}
]
}
]
}


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.

$$\frac{d S(t)}{d t}=A-\mu S(t)-\frac{r I}{N_0} S(t)$$
$$\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;