CODEWITHSUNDEEP
pythonpython-3.xodeodeint
ParaH2
ParaH2

Reputation: 517

Why these two ODE solving program don't match?

I don't understand why this code (reference):

from numpy import zeros, linspace
import matplotlib.pyplot as plt

# Time unit: 1 h
beta = 10./(40*8*24)
gamma = 3./(15*24)
dt = 0.1             # 6 min
D = 30               # Simulate for D days
N_t = int(D*24/dt)   # Corresponding no of hours

t = linspace(0, N_t*dt, N_t+1)
S = zeros(N_t+1)
I = zeros(N_t+1)
R = zeros(N_t+1)

# Initial condition
S[0] = 50
I[0] = 1
R[0] = 0

# Step equations forward in time
for n in range(N_t):
    S[n+1] = S[n] - dt*beta*S[n]*I[n]
    I[n+1] = I[n] + dt*beta*S[n]*I[n] - dt*gamma*I[n]
    R[n+1] = R[n] + dt*gamma*I[n]

fig = plt.figure()
l1, l2, l3 = plt.plot(t, S, t, I, t, R)
fig.legend((l1, l2, l3), ('S', 'I', 'R'), 'upper left')
plt.xlabel('hours')
plt.show()

doesn't produce the same results as this code I made:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt


beta = 10. / (40 * 8 * 24)
gamma = 3. / (15 * 24)


def flu(y, t):
    S = y[0]
    P = y[1]
    R = y[2]

    S = - beta * S * P
    P = beta * S * P - gamma * P
    R = gamma * P

    return [S, P, R]


C_I = [50, 1, 0]
t = np.linspace(0, 1000, 1000)

y = odeint(flu, C_I, t)

S = y[:, 0]
P = y[:, 1]
R = y[:, 2]

fig, ax = plt.subplots()
ax.plot(t, S, 'b--', label='S')
ax.plot(t, P, 'r--', label='I')
ax.plot(t, R, 'g--', label='R')
legend = ax.legend(loc='upper right', shadow=True, fontsize='x-large')
legend.get_frame().set_facecolor('#FFFCCC')
plt.show()

I used P instead of I to avoid confusion.

The equations solved with odeint should be the same as the ones provided in the reference link above. And if the equations I use are correct, which I am convinced they are, I don't understand where the mistake(s) lie(s).

Thank you for your help

Upvotes: 0

Views: 59

Answers (1)

Matthew Flamm
Matthew Flamm

Reputation: 595

You set S=y[0] then you set S=- beta * S * P. This overwrites y[0]!!! Similar problems for P and R

Try this:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt


beta = 10. / (40 * 8 * 24)
gamma = 3. / (15 * 24)


def flu(y, t):
    S = y[0]
    P = y[1]
    R = y[2]

    dS = - beta * S * P
    dP = beta * S * P - gamma * P
    dR = gamma * P

    return [dS, dP, dR]


C_I = [50, 1, 0]
t = np.linspace(0, 1000, 1000)

y = odeint(flu, C_I, t)

S = y[:, 0]
P = y[:, 1]
R = y[:, 2]

fig, ax = plt.subplots()
ax.plot(t, S, 'b--', label='S')
ax.plot(t, P, 'r--', label='I')
ax.plot(t, R, 'g--', label='R')
legend = ax.legend(loc='upper right', shadow=True, fontsize='x-large')
legend.get_frame().set_facecolor('#FFFCCC')
plt.show()

Upvotes: 1

Related Questions