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

4.2 SpreadingofDiseases 117 andone for the right-handside functions, f.u;t/D . ˇSI;ˇSI I; I/: Theequationu0 Df.u;t/meanssetting thetwovectorsequal, i.e., thecomponents mustbepairwiseequal. Sinceu0 D .S0;I 0;R0/,weget thatu0 Df implies S0 D ˇSI; I 0 DˇSI I; R0 D I : The generalized short notation u0 D f.u;t/ is very handy since we can derive numericalmethods and implement software for this abstract system and in a par- ticular application just identify the formulas in thef vector, implement these, and call functionality that solves thedifferential equationsystem. 4.2.6 ProgrammingtheNumericalMethod;theGeneralCase InPythoncode, theForwardEuler step unC1 DunC tf.un;tn/; beingascalaror avectorequation, canbecodedas u[n+1] = u[n] + dt*f(u[n], t[n]) both in the scalar and vector case. In the vector case, u[n] is a one-dimensional numpyarrayof lengthmC1holding themathematicalquantityun, and thePython functionfmust returnanumpyarrayof lengthmC1. Then theexpressionu[n] + dt*f(u[n], t[n]) is anarrayplusa scalar timesanarray. For all this to work, the complete numerical solutionmust be represented by a two-dimensional array, createdbyu = zeros((N_t+1, m+1)). Thefirst index counts the timepoints and the second thecomponentsof the solutionvector at one timepoint. That is,u[n,i]correspondsto themathematicalquantityuni .Whenwe useonlyoneindex,asinu[n], this is thesameasu[n,:]andpicksoutall thecom- ponents in thesolutionat the timepointwith indexn. Then theassignmentu[n+1] = ... becomes correct because it is actually an in-place assignment u[n+1, :] = .... Thenice featureof these facts is that the samepieceofPythoncodeworks forbotha scalarODEandasystemofODEs! Theode_FE function for thevectorODEisplaced in thefileode_system_FE. pyandwaswrittenas follows: from numpy import linspace, zeros, asarray import matplotlib.pyplot as plt def ode_FE(f, U_0, dt, T): N_t = int(round(float(T)/dt))