Modelling and Simulation: Theory and Examples

3. Differentiation and Integration

3.1. Using Python to find derivatives and integrals.

Example 1

Use the Python sympy package to solve the following

(a) \(\frac{dy}{dx} = x^2\)
(b) \(\frac{dy}{dx} = cos(x)\)
(c) \(\frac{dy}{dx} = e^{-x^2}\)
(d) \(\frac{d^3y}{dx^3} = x^4\)
(e) \(\int cos(x)\)
(f) \(\int_0^\infty e^{-x}\)
(g) \(\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} e^{- x^{2} - y^{2}}\, dx\, dy\)

Use the Python scipy package to solve and plot the following

(a) \(\frac{dy(t)}{dt} = -k \ y(t)\) \(k = 0.3; y_0 = 5\)
(b) \(\frac{dy(t)}{dt} = ke^t + y\)

Solution:

sympy

 1from sympy import *
 2# create a "symbol" called x
 3x = Symbol('x')
 4#Define function
 5f = x**2
 6#Calculating Derivative
 7derivative_f = f.diff(x)
 8print(derivative_f)
 9x,y,z = symbols('x y z')
10#sympy.init_printing()
11print(diff(cos(x), x))
12print(diff(exp(x**2), x))
13# third derivative of x^4
14print(diff(x**4, x, x, x))
15print(diff(x**4, x, 3))
16
17print(integrate(cos(x), x))
18print(integrate(exp(-x), (x, 0, oo)))
19print(integrate(exp(-x**2 - y**2), (x, -oo, oo), (y, -oo, oo)))

2*x
-sin(x)
2*x*exp(x**2)
24*x
24*x
sin(x)
1
pi

scipy

\[\begin{align*} \frac{dy(t)}{dt} &= -k \; y(t) \\ k &= 0.3 \\ y_0 &= 5 \\ \frac{dy(t)}{dt} &= ke^t + y \end{align*}\]

 1import numpy as np
 2import matplotlib.pyplot as plt
 3from scipy.integrate import odeint
 4
 5def model(y,t):
 6    k = 0.3
 7    dydt = -k * y 
 8    return dydt 
 9
10y01 = 5
11t1 = np.linspace(0,20)
12y1 = odeint(model, y01, t1)
13
14def model2(y, t, a):
15    dydt = a*np.exp(t) + y 
16    return dydt 
17
18y02 = 1 
19t2 = np.linspace(0, 5) 
20
21y2= odeint(model2, y02, t2, args=(13,)) 
22
23fig, ax = plt.subplots(1,2, figsize=(10,4))
24
25ax[0].plot(t1,y1)
26ax[0].set_xlabel('t')
27ax[0].set_ylabel('y')
28ax[0].set_title(r'$\frac{dy(t)}{dt} = -k \; y(t)$')
29
30ax[1].plot(t2,y2)
31ax[1].set_xlabel('t')
32ax[1].set_title(r'$\frac{dy(t)}{dt} = ke^t + y$')
33#plt.savefig(r'../images/scipy_ode.png')
34plt.show()

Figure 2. Scipy ode calculation

3.2. Spread of virus using the SIR Model.

Example 2

A new infectious flu virus is discovered to have infected \(10\%\) of a community. The rate of infection is \(0.35\) and that of recovery is \(0.1\). Given that no one has recovered so far, use the SIR Model to simulate the spread of this virus.

Solution:

 1import numpy as np
 2import matplotlib.pyplot as plt 
 3from scipy.integrate import odeint 
 4
 5# model 
 6def SIR_model(y, t, betta, gamma):
 7    S, I, R = y
 8    dSdt = -beta*S*I 
 9    dIdt = beta*S*I - gamma*I 
10    dRdt = gamma*I 
11    return ([dSdt, dIdt, dRdt]) 
12
13# initial conditions 
14S0 = 0.9 # 90 % 
15I0 = 0.1 
16R0 = 0.0 
17beta = 0.35 
18gamma = 0.1 
19t = np.linspace(0, 100, 10000) # create the simul time 
20
21# result 
22
23sol = odeint(SIR_model, [S0, I0, R0], t, args=(beta, gamma)) 
24sol = np.array(sol) 
25
26plt.figure(figsize=[6,4]) 
27plt.plot(t, sol[:, 0], 'b', label=r"$\frac{dS}{dt}$") 
28plt.plot(t, sol[:, 1], 'r', label=r"$\frac{dI}{dt}$") 
29plt.plot(t, sol[:, 2], 'lime', label=r"$\frac{dR}{dt}$") 
30plt.xlabel("Unit Time"); plt.ylabel("Percentage")
31plt.legend()
32plt.savefig("../images/sir_model.png")
33plt.show()

Figure 3. SIR model output.

Add Comment

Comments

No comments yet. Be the first!