```python
import sympy as sym
from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt
```


```python
S, I, R = sym.Function("S"), sym.Function("I"), sym.Function("V")
N, mu, alpha, beta, t = sym.symbols("N, mu, alpha, beta, t")
```


```python
eq1 = sym.Derivative(S(t), t) - (- alpha * S(t) * I(t) - mu * R(t))
eq2 = sym.Derivative(I(t), t) - (alpha * I(t) * S(t) / N - beta * I(t))
eq3 = sym.Derivative(R(t), t) - (beta * I(t) + mu * R(t))
```


```python
sym.dsolve((eq1, eq2, eq3))
```

Further investigation shows that an exact solution to this system of differential equations is difficult. Let us do this numerically:


```python
def dx(x, t, alpha, beta, mu):

    return (- alpha * x[1] * x[0] - mu * x[0],
            alpha * x[1] * x[0]  - beta * x[1],
            beta * x[1] + mu * x[0])
```


```python
alpha = 1 / 1000  # Every 1000 interactions leads to infection
beta = 1 / 5  # take 5 time units to recover
N = 10 ** 4  # Population of 10 thousand people
mu = 0  # 0 vaccination percentage

ts = np.linspace(0, 10, 5000)
xs = integrate.odeint(func=dx, y0=np.array([N - 1, 1, 0]), t=ts, args=(alpha, beta, mu))
S, I, R = xs.T
plt.figure()
plt.plot(ts, S, label="Susceptibles")
plt.plot(ts, I, label="Infected")
plt.plot(ts, R, label="Recovered")
plt.legend()
plt.title(f"$\max(I)={round(max(I))}$ ($\\alpha={alpha}$, $\\beta={beta}$, $\mu={mu}$)")
plt.savefig("base_scenario.pdf");
```


```python
mu = 1 / 2  # Vaccinate half the population
ts = np.linspace(0, 10, 5000)
xs = integrate.odeint(func=dx, y0=np.array([N - 1, 1, 0]), t=ts, args=(alpha, beta, mu))
S, I, R = xs.T
plt.figure()
plt.plot(ts, S, label="Susceptibles")
plt.plot(ts, I, label="Infected")
plt.plot(ts, R, label="Recovered")
plt.legend()
plt.title(f"$\max(I)={round(max(I))}$ ($\\alpha={alpha}$, $\\beta={beta}$, $\mu={mu}$)")
plt.savefig("moderate_vaccination_rate.pdf");
```


```python
mu = 99 / 100  # Vaccinate 99% of the population
ts = np.linspace(0, 10, 5000)
xs = integrate.odeint(func=dx, y0=np.array([N - 1, 1, 0]), t=ts, args=(alpha, beta, mu))
S, I, R = xs.T
plt.figure()
plt.plot(ts, S, label="Susceptibles")
plt.plot(ts, I, label="Infected")
plt.plot(ts, R, label="Recovered")
plt.legend()
plt.title(f"$\max(I)={round(max(I))}$ ($\\alpha={alpha}$, $\\beta={beta}$, $\mu={mu}$)")
plt.savefig("high_vaccination_rate.pdf");
```


```python
vaccination_rates = np.linspace(0, 1, 500)
max_percent_of_infected = []
for mu in vaccination_rates:
    xs = integrate.odeint(func=dx, y0=np.array([N - 1, 1, 0]), t=ts, args=(alpha, beta, mu))
    S, I, R = xs.T
    max_percent_of_infected.append(max(I) / N)
plt.figure()
plt.plot(vaccination_rates, max_percent_of_infected)
plt.xlabel("$\mu$")
plt.ylabel("% of population infected")
plt.savefig("effect_of_vaccination_rate.pdf");
```

77b20038a74e92a0044d418e4df7516270f31fd5
ipynb
Jupyter Notebook
assets/rsc/sir/main.ipynb
geraintpalmer/cfm
fa3f98cf45b225015f28be461e8ae661fa966b61
MIT
assets/rsc/sir/main.ipynb
geraintpalmer/cfm
fa3f98cf45b225015f28be461e8ae661fa966b61
MIT
assets/rsc/sir/main.ipynb
geraintpalmer/cfm
fa3f98cf45b225015f28be461e8ae661fa966b61
MIT
Qwen/Qwen-72B
1. YES
2. YES