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

114 4 SolvingOrdinaryDifferentialEquations 4.2.3 ProgrammingtheNumericalMethod;theSpecialCase The computation of (4.26)–(4.28) can be readily made in a computer program SIR1.py: 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() plt.savefig(’tmp.pdf’); plt.savefig(’tmp.png’) This programwaswritten to investigate the spreading of aflu at thementioned boardingschool, andthereasoningfor thespecificchoicesˇ and goesas follows. At someother schoolwhere thediseasehasalready spread, itwasobserved that in thebeginningofadaytherewere40susceptiblesand8infected,while thenumbers were 30 and18, respectively, 24hours later. Using 1h as time unit,we then have from (4.11) thatˇ D 10=.40 8 24/. Among 15 infected, it was observed that 3 recovered during a day, giving D 3=.15 24/. Applying these parameters to anewcasewherethere isoneinfectedinitiallyand50susceptibles,gives thegraphs inFig.4.9. Thesegraphsare just straight linesbetweenthevaluesat times ti D i t as computedby theprogram.Weobserve thatS reducesas I andRgrows. After about30dayseveryonehasbecomeill and recoveredagain. Wecanexperimentwithˇand toseewhetherwegetanoutbreakofthedisease or not. Imagine that a “wash your hands” campaignwas successful and that the other school in this case experienced a reduction ofˇ by a factor of 5. With this