Page - 119 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

4.2 SpreadingofDiseases 119 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 = dt*N_t # End time U_0 = [50, 1, 0] u, t = ode_FE(f, U_0, dt, T) S = u[:,0] I = u[:,1] R = u[:,2] fig = plt.figure() l1, l2, l3 = plt.plot(t, S, t, I, t, R) fig.legend((l1, l2, l3), (’S’, ’I’, ’R’), ’lower right’) plt.xlabel(’hours’) plt.show() # Consistency check: N = S[0] + I[0] + R[0] eps = 1E-12 # Tolerance for comparing real numbers for n in range(len(S)): SIR_sum = S[n] + I[n] + R[n] if abs(SIR_sum - N) > eps: print ’*** consistency check failed: S+I+R=%g != %g’ %\ (SIR_sum, N) if __name__ == ’__main__’: demo_SIR() Recall that the u returned fromode_FE contains all components (S, I,R) in the solution vector at all time points. We therefore need to extract the S, I, andR values in separatearrays for furtheranalysis andeasyplotting. Another key feature of this higher-quality code is the consistency check. By adding the threedifferential equations in theSIRmodel,we realize thatS0 CI 0 C R0 D0,whichmeansthatSCICRD const.Wecancheckthat this relationholds bycomparingSnCIn CRn to the sumof the initial conditions. Thecheck isnot afull-fledgedverification,but it isamuchbetter thandoingnothingandhopingthat thecomputationiscorrect. Exercise4.5suggestsanothermethodforcontrollingthe qualityof thenumerical solution. 4.2.7 Time-RestrictedImmunity Letusnowassume that immunityafter thediseaseonly lasts for somecertain time period. Thismeans that there is transport fromtheRstate to theSstate: