Search code examples
pythonmodelrunge-kutta

How can I find out the amount of susceptible, infected and recovered individuals in time = 50, where S(50), I(50), R(50)? (SIR MODEL)

How can I find out the amount of susceptible, infected and recovered individuals in time = 50, where S(50), I(50), R(50)? (SIR MODEL)

# Equações diferenciais e suas condições iniciais
h = 0.05
beta = 0.8
nu = 0.3125

def derivada_S(time,I,S):
    return -beta*I*S

def derivada_I(time,I,S):
    return beta*I*S - nu*I

def derivada_R(time,I):
    return nu*I

S0 = 0.99
I0 = 0.01
R0 = 0.0

time_0 = 0.0
time_k = 100
data = 1000

# vetor representativo do tempo
time = np.linspace(time_0,time_k,data)

S = np.zeros(data)
I = np.zeros(data)
R = np.zeros(data)

S[0] = S0
I[0] = I0
R[0] = R0

for i in range(data-1):
    S_k1 = derivada_S(time[i], I[i], S[i])
    S_k2 = derivada_S(time[i] + (1/2)*h, I[i], S[i] + h + (1/2)*S_k1)
    S_k3 = derivada_S(time[i] + (1/2)*h, I[i], S[i] + h + (1/2)*S_k2)
    S_k4 = derivada_S(time[i] + h,  I[i], S[i] + h + S_k3)
    
    S[i+1] = S[i] + (h/6)*(S_k1 + 2*S_k2 + 2*S_k3 + S_k4)

    I_k1 = derivada_I(time[i], I[i], S[i])
    I_k2 = derivada_I(time[i] + (1/2)*h, I[i], S[i] + h + (1/2)*I_k1)
    I_k3 = derivada_I(time[i] + (1/2)*h, I[i], S[i] + h + (1/2)*I_k2)
    I_k4 = derivada_I(time[i] + h,  I[i], S[i] + h + I_k3)
    
    I[i+1] = I[i] + (h/6)*(I_k1 + 2*I_k2 + 2*I_k3 + I_k4)
    
    R_k1 = derivada_R(time[i], I[i])
    R_k2 = derivada_R(time[i] + (1/2)*h, I[i])
    R_k3 = derivada_R(time[i] + (1/2)*h, I[i])
    R_k4 = derivada_R(time[i] + h, I[i])
    
    R[i+1] = R[i] + (h/6)*(R_k1 + 2*R_k2 + 2*R_k3 + R_k4)

plt.figure(figsize=(8,6))
plt.plot(time,S, label = 'S')
plt.plot(time,I, label = 'I')
plt.plot(time,R, label = 'R')
plt.xlabel('tempo (t)')
plt.ylabel('Susceptível, Infectado e Recuperado')
plt.grid()
plt.legend()
plt.show()

I'm solving an university problem with python applying Runge-Kutta's fourth order, but a I don't know how to collect the data for time = 50.

Solution

  • The easiest way to get a value at time 50 is to compute a value at time 50. As you compute data over 100 days, with about 10 data points per day, reflect this in the time array construction

    time = np.linspace(0,days,10*days+1)

    Note that linspace(a,b,N) produces N nodes that have between them a step of size (b-a)/(N-1).

    Then you get the data for day 50 at time index 500 (and the 9 following).

    For this slow-moving system and this relatively small time step, you will get reasonable accuracy with the implemented order-1 method, but will get better accuracy with a higher-order method like RK4.

    You need to apply associated updates to all components everywhere. This requires to interleave the RK4 steps that you have, as for instance the corrected step

       S_k2 = derivada_S(time[i] + (h/2), I[i] + (h/2)*I_k1, S[i] + (h/2)*S_k1)

    requires that the value I_k1 is previously computed. Note also that h should be a factor to S_k1, it should not be added.

    In total you should get

    for i in range(data-1):
    S_k1 = derivada_S(time[i], I[i], S[i])
    I_k1 = derivada_I(time[i], I[i], S[i])
    R_k1 = derivada_R(time[i], I[i])

    S_k2 = derivada_S(time[i] + (1/2)*h, I[i] + (h/2)*I_k1, S[i] + (h/2)*S_k1)
    I_k2 = derivada_I(time[i] + (1/2)*h, I[i] + (h/2)*I_k1, S[i] + (h/2)*S_k1)
    R_k2 = derivada_R(time[i] + (1/2)*h, I[i] + (h/2)*I_k1)

    S_k3 = derivada_S(time[i] + (h/2), I[i] + (h/2)*I_k2, S[i] + (h/2)*S_k2)
    I_k3 = derivada_I(time[i] + (h/2), I[i] + (h/2)*I_k2, S[i] + (h/2)*S_k2)
    R_k3 = derivada_R(time[i] + (h/2), I[i] + (h/2)*I_k2)

    S_k4 = derivada_S(time[i] + h, I[i] + I_k3, S[i] + S_k3)
    I_k4 = derivada_I(time[i] + h, I[i] + I_k3, S[i] + S_k3)
    R_k4 = derivada_R(time[i] + h, I[i] + I_k3)
    
    S[i+1] = S[i] + (h/6)*(S_k1 + 2*S_k2 + 2*S_k3 + S_k4)
    I[i+1] = I[i] + (h/6)*(I_k1 + 2*I_k2 + 2*I_k3 + I_k4)
    R[i+1] = R[i] + (h/6)*(R_k1 + 2*R_k2 + 2*R_k3 + R_k4)

    Note that h is a factor to I_k1, S_k1 etc. You have a sum there.

    Replacing just this piece of code gives the plot

    enter image description here

    But there is another problem before that. You defined the time step as 0.05 so that t=50 is reached at the last place. As the system is autonomous, the contents of the time array makes no difference, but the labeling of the x axis has to be divided by 2. The values that you want are in fact the last values computed with data = 10*time_k+1.

    S[-1]=0.10483, I[-1]=8.11098e-05, R[-1]=0.89509

    For the previous discussion to remain valid, you could also set h=t[1]-t[0], so that t=50 is reached in the middle at i=500.