High order discretizations for spatial dependent SIR models
HHigh order discretization methods for spatial dependent SIRmodels
Bálint Takács ∗ Yiannis Hadjimichael † July 13, 2020
Abstract
In this paper, an SIR model with spatial dependence is studied and results regarding itsstability and numerical approximation are presented. We consider a generalization of theoriginal Kermack and McKendrick model in which the size of the populations differs in space.The use of local spatial dependence yields a system of integro-differential equations. Theuniqueness and qualitative properties of the continuous model are analyzed. Furthermore,different choices of spatial and temporal discretizations are employed, and step-size restric-tions for population conservation, positivity and monotonicity preservation of the discretemodel are investigated. We provide sufficient conditions under which high order numericalschemes preserve the discrete properties of the model. Computational experiments verifythe convergence and accuracy of the numerical methods.
During the millenia of the history of mankind, many epidemics have ravaged the population.Since the plague of Athens in 430 BC described by historian Thucydides (one of the earliestdescription of such epidemics), researchers tried to model and describe the outbreak of illnesses.More recently, the outbreak of COVID-19 pandemic revealed the importance of epidemic re-search and the development of models to describe the public health, social, and economic impactof major virus diseases.Nowadays many of the models used in science are derived from the original ideas of Kermackand McKendrick [26] in 1927, who constructed a compartment model to study the process ofepidemic propagation. In their model the population is split into three classes: S being thegroup of healthy individuals, who are susceptible to infection; I is the compartment of theill species, who can infect other individuals; and R being the class of recovered or immuneindividuals. ∗ Applied Analysis and Computational Mathematics, MTA-ELTE Numerical Analysis and Large Net-works Research Group, Eötvös Loránd University, Pázmány P. s. 1/C, Budapest 1117, Hungary, (email: [email protected] ). † MTA-ELTE Numerical Analysis and Large Networks Research Group, Eötvös Loránd University, PázmányP. s. 1/C, Budapest 1117, Hungary, (email: [email protected] ). a r X i v : . [ m a t h . NA ] J u l he original model of Kermack and McKendrick took into account constant rates of changeand neglected any natural deaths and births or vaccination. In this work, we also considerconstant rates of change and moreover we include a term, c S ( t ) , that describes immunizationeffects through vaccination. The SIR model takes the form ddt S ( t ) = − a S ( t ) I ( t ) − c S ( t ) , ddt I ( t ) = a S ( t ) I ( t ) − b I ( t ) , ddt R ( t ) = b I ( t ) + c S ( t ) , (1.1)where the positive constant parameters a , b and c respectively correspond to the rate of infection,recovery and vaccination.Since the introduction of the model (1.1) in 1927, numerous extensions were constructed todescribe biological processes more efficiently and realistically. A natural extension is to takeinto account the heterogeneity of our domain in a way that we examine not only the change ofthe populations in time, but also we observe the spatial movements. Kendall introduced suchmodels that transformed the system of ordinary differential equations (1.1) into a system ofpartial differential equations [24, 25].The time-dependent functions in (1.1) represent the number of individuals in each class, butcontain no information about their spatial distribution. Instead, one can replace these concen-tration functions with spatial-dependent functions describing the density of healthy, infectiousand recovered species over some domain in R d [35]. In this paper we consider a bounded domainin R , hence the system (1.1) can be recast as ∂∂t S ( t , x , y ) = − a S ( t , x , y ) I ( t , x , y ) − c S ( t , x , y ) , ∂∂t I ( t , x , y ) = a S ( t , x , y ) I ( t , x , y ) − b I ( t , x , y ) , ∂∂t R ( t , x , y ) = b I ( t , x , y ) + c S ( t , x , y ) . (1.2)However, the model (1.2) is still insufficient as it does not allow the disease to spread in thedomain but only accounts for a point-wise infection. Spatial points do not interact with eachother but infect species only at their location. In order to allow a realistic propagation of theinfection, we assume that an infected individual can spread the disease on susceptible speciesin a certain area around its location. Let us define a non-negative function G ( x , y , r , θ ) = g ( r ) g ( θ ) , if ( ¯ x ( r , θ ) , ¯ y ( r , θ )) ∈ B δ (cid:16) x , y (cid:17) ,0, otherwise, (1.3)that describes the effect of a single point ( x , y ) in a δ -radius neighborhood B δ (cid:16) x , y (cid:17) , and set ¯ x ( r , θ ) = x + r cos ( θ ) and ¯ y ( r , θ ) = y + r sin ( θ ) . The function G ( x , y , r , θ ) demonstrates howhealthy individuals at points ( ¯ x ( r , θ ) , ¯ y ( r , θ )) are infected by the center point ( x , y ) , where r ∈ [ δ ] is the distance from the center and θ ∈ [
0, 2 π ) is the angle. We assume that the2ight-hand-side of (1.3) is separable. The effect of the point ( x , y ) depending on the distancefrom the center is described by g ( r ) ; a decreasing, non-negative function that is equal to zerofor values r ≥ δ (since there is no effect outside B δ (cid:16) x , y (cid:17) ). Function g ( θ ) characterizes thepart of the effect depending on the angle, i.e., the direction in which the center is compared topoint ( ¯ x ( r , θ ) , ¯ y ( r , θ )) . The case of constant function g ( θ ) is widely studied in [12] and [13],while a non-constant such function may be useful in the case of modeling the spread of diseasesin a forest with a constant wind blowing in one direction which was described in [35]. In bothcases it is supposed that the function is periodic in the sense that g ( ) = lim θ → π g ( θ ) .The nonlinear terms of the right-hand side of (1.2) describe the interaction of susceptibleand infected species. We can now utilize (1.3) and replace the density of infected species inthese nonlinear terms by Z δ Z π G ( x , y , r , θ ) I ( t , ¯ x ( r , θ ) , ¯ y ( r , θ )) r dθdr ,where we used the fact that G ( x , y , r , θ ) = B δ (cid:16) x , y (cid:17) . Therefore, the model(1.2) can be expressed as a system of integro-differential equations ∂S ( t , x , y ) ∂t = − S ( t , x , y ) Z δ Z π g ( r ) g ( θ ) I ( t , ¯ x ( r , θ ) , ¯ y ( r , θ )) r dθdr − cS ( t , x , y ) , ∂I ( t , x , y ) ∂t = S ( t , x , y ) Z δ Z π g ( r ) g ( θ ) I ( t , ¯ x ( r , θ ) , ¯ y ( r , θ )) r dθdr − bI ( t , x , y ) , ∂R ( t , x , y ) ∂t = bI ( t , x , y ) + cS ( t , x , y ) . (1.4) The aim of this paper is twofold. First, in Section 2 we analyze the stability of the continuousmodel (1.1) and prove that a unique solution exists under some Lipschitz continuity and bound-edness assumptions. Secondly, in Sections 3 and 4 we seek numerical schemes that approximatethe solution of (1.1) and maintain its qualitative properties.We verify that the analytic solution satisfies biologically reasonable properties; however, asshown in Section 2.1 the solution can be only expressed implicitly in terms of S , I , and R and thus it is not directly applicable. A numerical approximation is presented in Section 2.2that provably satisfies the solution’s properties. The first order accuracy of this approxima-tion motivates the search for suitable high order numerical methods that preserve a discreteanalogue of the properties of the continuous model. In Section 3 we use cubature formulas toreduce the integro-differential system (1.4) to an ODE system. We study the accuracy of dif-ferent cubatures and interpolation techniques for approximating the multiple integrals in (1.1).Futhermore, the employment of time integration methods yields an algebraic system to solvenumerically. Section 4 shows that a time-step restriction is sufficient and necessary such thatthe forward Euler method maintains the stability properties of the ODE system. We prove thathigh order strong-stability-preserving (SSP) Runge–Kutta methods can be used under appropri-ate restrictions; thus, we can obtain a high order stable scheme both in space and time. Finallyin Section 5 we demonstrate the theoretical results by conducting numerical experiments.3 Stability of the analytic solution
Analytic results for deterministic reaction epidemic models have been studied by several authors,see for example, [4, 25, 36]. Such models lie in the larger class of reaction-diffusion problems andtherefore one can obtain theoretical results by studying the more general problem. We provethe uniqueness of the solution for system (1.4) by following the work of Capasso and Fortunato[6]. We consider the following semilinear autonomous evolution problem ∂u∂t ( t ) = − Au ( t ) + F ( u ( t )) u = u ( ) ∈ D ( A ) , (2.1)where A is a self-adjoint and positive-definite operator in a real Hilbert space E with domain D ( A ) . Define λ = inf σ ( A ) , where σ ( A ) denotes the spectrum of A . Let us choose E : = L ( Ω ) × L ( Ω ) , where Ω is a bounded domain in R , with a norm defined by (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) u u !(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) = (cid:16) k u k L + k u k L (cid:17) . (2.2)Here u = ( u , u ) (cid:124) ∈ C (cid:16) [ t f ) , D ( A ) (cid:17) , for some final time t f . We also equip D ( A ) with thenorm k u k A = k Au k , u ∈ D ( A ) .Note that it is sufficient to consider only the first two equations in (1.4), since R ( t , x , y ) canbe obtained by using that the sum S ( t , x , y ) + I ( t , x , y ) + R ( t , x , y ) is constant in time for everypoint ( x , y ) . Hence, in view of problem (1.4), the linear operator A is defined as A u u ! : = c b ! u u ! , (2.3)and D ( A ) = E . Because b and c are positive constants, it is easy to see that A is a self-adjointand positive-definite operator. Similarly, F ( u ) consists of the nonlinear terms, and is definedas F u u ! : = − u F ( u ) u F ( u ) ! . (2.4)The function F : L ( Ω ) → L ( Ω ) contains the integral part of (1.4) and is given by F ( t , x , y ) : = F ( I ( t , x , y )) = Z δ Z π g ( r ) g ( θ ) I ( t , ¯ x ( r , θ ) , ¯ y ( r , θ )) r dθdr , (2.5)where I ( t , x , y ) can be viewed as a map I ( t , x , y ) : [ t f ) I t ( x , y ) ∈ L ( Ω ) .The main result of this section is Theorem 2.1 stating that a unique solution of system (1.4)exists. Theorem 2.1 considers the system (2.1) as a generalization of (1.4) and its proof relieson the fact that the function F in (2.4) is Lipschitz-continuous and bounded in k·k A . Therefore,we define the following conditions [6]: 4 A ) F is locally Lipschitz-continuous from D ( A ) to D ( A ) , i.e., k F ( u ) − F ( v ) k A ≤ ζ ( d ) k u − v k A for all u , v ∈ D ( A ) such that d ≥
0, and k u k A ≤ d , k v k A ≤ d .( A ) F is bounded, i.e., there exists ν ≥
0, and γ ≥ k F ( u ) k A ≤ ν k u k + γA , ∀ u ∈ D ( A ) .We also let κ = max r ∈ ( δ ) { g ( r ) } , κ = max θ ∈ ( π ) { g ( θ ) } , and µ ( Ω ) be the Lebesguemeasure of Ω . Theorem 2.1.
Consider the problem (1.4) and assume that conditions ( A ) and ( A ) hold.Then, a unique strong solution solution of system (1.4) exists on some interval [ t f [ . Moreover,if any initial condition u belongs in the set K = ( u ∈ E (cid:12)(cid:12)(cid:12)(cid:12) k u k A < min ( b , c ) √ κ κ µ ( Ω ) ) , then the zero solution is the unique equilibrium solution of (1.4) . The proof of Theorem 2.1 is a direct consequence of two main results by Capasso andFortunato [6]. For clarity, we state these two theorems below.
Theorem 2.2. [6, Theorem 1.1] If assumption ( A ) holds, a unique strong solution in D ( A ) of problem (2.1) exists in some interval [ t f [ . Theorem 2.3. [6, Theorem 1.3] Let us assume that ( A ) and ( A ) hold. Then for any u ∈ f K aglobal strong solution in D ( A ) , u ( t ) , of (2.1) exists. Moreover the zero solution is asymptoticallystable in f K . Here f K = ( { u ∈ D ( A ) (cid:12)(cid:12)(cid:12) k u k A < ( λ / ν ) / γ } , if γ > D ( A ) , if γ = and λ > ν .In the rest of the section we show that the function F as defined in (2.4) satisfies conditions( A ) and ( A ). First, to prove that ( A ) holds, we make use of some auxiliary lemmas; theirproofs appear in Appendix A. Lemma 2.1.
The norms k·k and k·k A are equivalent, i.e., k u k ≤ max (cid:18) c , 1 b (cid:19) k u k A and k u k A ≤ max ( c , b ) k u k . Lemma 2.2.
Let F be given by (2.5) . Then, it holds that kF ( u ) k L ≤ ν F k u k L , where ν F = κ κ µ ( Ω ) . orollary 2.1. Consider F given by (2.4) . Then, the condition ( A ) holds with γ = and ν = √ ( c , b ) max (cid:18) c , 1 b (cid:19) κ κ µ ( Ω ) .Proof. Because of Lemma 2.1, it is enough to prove k F ( u ) k ≤ ˜ ν k u k , (2.6)since from this we get k F ( u ) k A ≤ max ( c , b ) k F ( u ) k ≤ max ( c , b ) ˜ ν k u k ≤ max ( c , b ) max (cid:18) c , 1 b (cid:19) ˜ ν k u k A .Now to prove inequality (2.6), consider that k F ( u ) k = (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) − u F ( u ) u F ( u ) !(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) = ( k u F ( u ) k L + k u F ( u ) k L ) / ≤ √ k u k L kF ( u ) k L .Observe that Lemma 2.2 can be used to bound kF ( u ) k L from above, yielding k F ( u ) k ≤ √ ν F k u k L k u k L ,where ν F is defined in Lemma 2.2. Finally, we have that k u k L k u k L ≤ k u k L + k u k L = k u k ,and thus inequality (2.6) holds with ˜ ν = √ ν F = √ κ κ µ ( Ω ) .The following lemma facilitates Lemma 2.4 and its proof can be found in Appendix A. Lemma 2.3.
Consider F given by (2.5) . Then, the following inequality holds: kF ( u ) − F ( v ) k L ≤ C F k u − v k L , where C F = κ κ µ ( Ω ) . Lemma 2.4.
Consider F given by (2.4) . Then, condition ( A ) holds.Proof. Because of Lemma 2.1, it is enough to prove k F ( u ) − F ( v ) k ≤ ν k u − v k . (2.7)If (2.7) holds, then k F ( u ) − F ( v ) k A ≤ max ( c , b ) k F ( u ) − F ( v ) k≤ max ( c , b ) ν k u − v k≤ max ( c , b ) max (cid:18) c , 1 b (cid:19) ν k u − v k A .6o show that inequality (2.7) holds, first consider that k F ( u ) − F ( v ) k = (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) − u F ( u ) + v F ( v ) u F ( u ) − v F ( v ) !(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ≤ √ k u F ( u ) − v F ( v ) k L .We can further bound the right-hand-side of the above inequality, yielding k u F ( u ) − v F ( v ) k L = k u F ( u ) − v F ( u ) + v F ( u ) − v F ( v ) k L ≤ k u F ( u ) − v F ( u ) k L + k v F ( u ) − v F ( v ) k L = kF ( u ) k L k u − v k L + kF ( u ) − F ( v ) k L k v k L . (2.8)Then, by Lemma 2.2 we have kF ( u ) k L k u − v k L ≤ ν F k u k L k u − v k L ,and also from Lemma 2.3 we get kF ( u ) − F ( v ) k L k v k L ≤ C F k u − v k L k v k L . (2.9)Assume that there exists d ≥
0, such that k u k ≤ d and k v k ≤ d . Then, by definition (2.2) wehave that k v k ≤ d and k u k ≤ d . Putting all together, we get k F ( u ) − F ( v ) k ≤ √ k u F ( u ) − v F ( v ) k L ≤ √ (cid:16) d ν F k u − v k L + d C F k u − v k L (cid:17) / ≤ √ ( dC F , d ν F ) k u − v k .Therefore, condition ( A ) holds with Lipschitz constant c ( d ) = √ ( c , b ) max (cid:18) c , 1 b (cid:19) max ( dC F , d ν F )= √ κ κ µ ( Ω ) max ( c , b ) max (cid:18) c , 1 b (cid:19) max ( d , d ) ,where we used that ν F = C F = κ κ µ ( Ω ) .Corollary 2.1 and Lemma 2.4 show that function (2.4) satisfies conditions ( A ) and ( A ).We know from Lemma 2.1 that γ =
1, so the set K in Theorem 2.3 can be computed by usingthat D ( A ) = E and (cid:18) λ ν (cid:19) / γ = min ( b , c ) √ κ κ µ ( Ω ) ,where b and c are the diagonal elements of matrix A in (2.3), and λ = inf σ ( A ) . Finally, it isevident that Theorem 2.1 follows from Theorems 2.2 and 2.3.7 .1 Qualitative behavior of the model When deriving a mathematical model to describe the spread of an epidemic in both space andtime, it is essential that the real-life processes are being represented as accurately as possible.More precisely, numerical discretizations applied to such models should preserve the qualitativeproperties of the original epidemic model.The first, and perhaps most natural property is that the number of each species is non-negative at every time and point of the domain. Next, assuming that the births and naturaldeaths are the same (vital dynamics have no effect on the process), the total number of speciesof all classes should be conserved. Another property concerns the number of susceptible species.Since an individual gets to the recovered class after the infection, the number of susceptiblescannot increase in time. Similarly, the number of recovered species cannot decrease in time.These properties can be expressed as follows: C : The densities X ( t , x , y ) , X ∈ { S , I , R } , are non-negative at every point ( x , y ) ∈ Ω . C : The sum S ( t , x , y ) + I ( t , x , y ) + R ( t , x , y ) is constant in time for all points ( x , y ) ∈ Ω ,yielding Z Ω (cid:16) S ( t , x , y ) + I ( t , x , y ) + R ( t , x , y ) (cid:17) d x d y = const., ∀ t . C : Function S ( t , x , y ) is non-increasing in time at every ( x , y ) ∈ Ω . C : Function R ( t , x , y ) is increasing in time at every ( x , y ) ∈ Ω .As in the previous section, instead of proving the preservation of properties C – C for the par-ticular model (1.4), we can establish theoretical results for a more general system of equations.First, we state the following lemma, whose proof can be found in Appendix A. Lemma 2.5.
The solution of (1.4) depends continuously on the right hand side of the systemof equations.
Let us define the operator J x , y ( t ) : [ t f ] → L ( Ω ) as J x , y ( t ) : = { I ( t , ¯ x , ¯ y ) : ( ¯ x , ¯ y ) ∈ B δ ( x , y ) } ,consisting of the infectious densities at points ( ¯ x , ¯ y ) lying in the δ -radius ball centered at point ( x , y ) at time t . The next theorem considers a generalization of system (1.4) and shows that itssolution satisfies properties C – C . Theorem 2.4.
Consider the following system of equations ∂S ( t , x , y ) ∂t = − S ( t , x , y ) H (cid:16) J x , y ( t ) (cid:17) − c S ( t , x , y ) , ∂I ( t , x , y ) ∂t = S ( t , x , y ) H (cid:16) J x , y ( t ) (cid:17) − b I ( t , x , y ) , ∂R ( t , x , y ) ∂t = b I ( t , x , y ) , (2.10)8 here H is a continuous functional mapping operators J x , y ( t ) to R . Suppose that H is non-negative in the sense that if φ ( t ) = φ t ∈ L ( Ω ) and φ t ( x , y ) ≥ ∀ ( x , y ) ∈ Ω , t ∈ [ t f ] ,then H ( φ t ) ≥ for all t ∈ [ t f ] . Also, suppose that the initial conditions of the system arenon-negative, i.e. X ( x , y ) ≥ , ∀ ( x , y ) ∈ Ω , X ∈ { S , I , R } . In such case, the properties C – C hold without any restriction on the time interval t ∈ [ t f ] .Proof. The proof consists of two parts: first we prove the required properties for a modifiedversion of (2.10), and then by using Lemma 2.5 we derive the statement of the theorem.Let us consider a modified version of (2.10) ∂S ε ( t , x , y ) ∂t = − S ε ( t , x , y ) H (cid:16) J x , y , ε ( t ) (cid:17) − c S ε ( t , x , y ) , ∂I ε ( t , x , y ) ∂t = S ε ( t , x , y ) H (cid:16) J x , y , ε ( t ) (cid:17) − b I ε ( t , x , y ) + ε , ∂R ε ( t , x , y ) ∂t = b I ε ( t , x , y ) , (2.11)where ε : R → R is a constant positive function, and J x , y , ε ( t ) is defined as J x , y , ε ( t ) : = { I ε ( t , ¯ x , ¯ y ) : ( ¯ x , ¯ y ) ∈ B δ ( x , y ) } .We also suppose that all non-negative initial conditions are assigned to the equation. First, wewould like to prove the non-negativity of I ε ( t , x , y ) . Proceeding towards contradiction, let ussuppose that the function takes negative values for some time t at some point ( x , y ) ∈ Ω . Letus define by t the first moment in time after which I ε ( t , x , y ) takes negative values, i.e., t : = inf { t | ∃ ( x , y ) ∈ Ω : I ε ( t , x , y ) < } .This t exists because I ε is continuous and the initial conditions are not negative, i.e., I ε ( x , y ) ≥
0. Because of the continuity of I ε and the definition of t , there is a point ( x , y ) for which I ε ( t , x , y ) =
0, and ∂I ε ( t , x , y ) ∂t ≤
0. (2.12)We know that all the values of I ε at t inside B δ ( x , y ) are non-negative by the definition of t ,and H is a non-negative operator in the sense defined before, so H (cid:16) J x , y , ε ( t ) (cid:17) ≥ ( t , x , y ) , we can see thatthe term − b I ε ( t , x , y ) is zero, so the term S ε ( t , x , y ) H (cid:16) J x , y , ε ( t ) (cid:17) must be negative forcondition (2.12) to hold (since ε is positive). We have already concluded that H (cid:16) J x , y , ε ( t ) (cid:17) ≥
0, so we need that S ε ( t , x , y ) < S ε , and integrate it with respect to time t from 0 to t . In this way we obtainlog ( S ε ( t , x , y )) − log ( S ε ( x , y )) = − Z t (cid:16) H (cid:16) J x , y , ε ( t ) (cid:17) − ct (cid:17) d t .9y reformulating, we get that for ( x , y ) = ( x , y ) : S ε ( t , x , y ) = S ε ( x , y ) exp (cid:18) − Z t (cid:16) H (cid:16) J x , y , ε ( t ) (cid:17) − ct (cid:17) d t (cid:19) . (2.13)Therefore S ε ( t , x , y ) is non-negative, so we get a contradiction.In this way I ε ( t , x , y ) ≥ t ∈ [ t f ] and ( x , y ) ∈ Ω . Consequently, since R ε ( x , y ) is non-negative, we get that R ε ( t , x , y ) is a non-decreasing and a non-negative function. Notealso that the calculations resulting in the formula (2.13) are also true for any time t and point ( x , y ) ∈ Ω , meaning that S ε is also non-negative, and since H (cid:16) J x , y , ε ( t ) (cid:17) is non-negative, wealso get the non-increasing property from the first equation of (2.10). Hence, we proved thatthe solution of (2.11) satisfies C – C .We also know that because of continuous dependence by Lemma 2.5,lim ε → X ε ( t , x , y ) (cid:12)(cid:12)(cid:12) t ∈ [ t f ] − X ( t , x , y ) (cid:12)(cid:12)(cid:12) t ∈ [ t f ] = X ∈ { S , I , R } . Therefore, properties C – C are satisfied by the solution ofsystem (2.10).Note that in the previous theorem it might happen that the functional J x , y ( t ) does notdepend on all of the values of function I ( t , x , y ) for ( x , y ) ∈ B δ ( x , y ) , but only on some of them.This special case will be useful in Section 3: see Remark 3.1.Due to the complicated form of the equations in (1.4) one can suspect that no analyticalsolution can be derived for this system. Because of this, we are going to use numerical methodsto approximate the solution of these equations.However, the analytic solution of the original SIR model (1.1) has been recently describedin the papers by Harko et al. [22] and Miller [29, 30]. Thus, we can get similar results applyingtheir observations to our modified model (1.4). The analytic solution of system (1.4) can bewritten as S ( t , x , y ) = S ( x , y ) e − φ ( t , x , y ) − ct , I ( t , x , y ) = M ( x , y ) − S ( t , x , y ) − R ( t , x , y ) , R ( t , x , y ) = R ( x , y ) + b Z t I ( s , x , y ) ds + c Z t S ( s , x , y ) ds , (2.14)where we use the notations M ( x , y ) : = S ( x , y ) + I ( x , y ) + R ( x , y ) , φ ( t , x , y ) : = Z t F ( I ( s , x , y )) ds ,and F is given by (2.5).It is evident that in (2.14), the values of the functions at a given time t ∗ can only becomputed if the values in the interval [ t ∗ ) are known. Consequently, these formulas are notuseful in practice, since (2.14) is an implicit system in the solutions S ( t , x , y ) , I ( t , x , y ) and R ( t , x , y ) . Later (see Table 5.3 in Section 5.2), an approximation of the solution of (2.14) willbe compared to the numerical solution of a first order scheme.10ince the values of the functions in (2.14) cannot be calculated directly, numerical methodsare needed to approximate them. We can take two possible paths:1. approximate the values of φ ( t , x , y ) by numerical integration; or2. approximate the solution of the original equation (1.4) by a numerical method.The first approach is discussed in Section 2.2, while the rest of the paper considers the secondcase. We focus on the order and convergence rate of our numerical methods, and ensure thatqualitative properties C – C of the analytic solution are preserved by the numerical method.For that a discrete analogue of conditions C – C is required; see section Section 4. As noted before, if we would like to use the solution (2.14) then we have to approximate theinvolved integrals. This can be achieved by partitioning the time interval [ t f ] into uniformspaced sections by using a constant time step τ . With this approach, the integrals can beapproximated by a left Riemann sum, and thus consider the values of densities X ( t , x , y ) , X ∈ { S , I , R } , at the left endpoint of each section. Therefore, for any integer 1 ≤ n ≤ N suchthat t f = τ N , the integral of X ( t , x , y ) can be approximated by Z nτ X ( s , x , y ) ds ≈ τ n − X k = X ( kτ , x , y ) .An important observation is that the integral equations (2.14) can be rewritten in a recursiveform S ( nτ , x , y ) = S (( n − ) τ , x , y ) exp − Z nτ ( n − ) τ F ( I ( s , x , y )) ds − cτ ! , R ( nτ , x , y ) = R (( n − ) τ , x , y ) + b Z nτ ( n − ) τ I ( s , x , y ) ds + c Z nτ ( n − ) τ S ( s , x , y ) ds , I ( nτ , x , y ) = M ( x , y ) − S ( nτ , x , y ) − R ( nτ , x , y ) . (2.15)Let X n ( x , y ) ≈ X ( nτ , x , y ) , X ∈ { S , I , R } , and define F n : = F ( I n ) . Using the approximations τ F n − ≈ Z nτ ( n − ) τ F ( I ( s , x , y )) ds , τ I n − ≈ Z nτ ( n − ) τ I ( s , x , y ) ds ,and choosing to approximate R nτ ( n − ) τ S ( s , x , y ) ds by τ S n yields an approximating scheme for(2.14), given by S n = S n − e − τ F n − − cτ , R n = R n − + bτ I n − + cτ S n , I n = ( S n − + I n − + R n − ) − S n − R n . (2.16a)(2.16b)(2.16c)Note that in this case, the order of the equations in (2.16) is important as estimates at time t n = nτ are used to update the rest of solution’s components.11 heorem 2.5. Consider the solution X n ( x , y ) , X ∈ { S , I , R } of scheme (2.16) on the timeinterval [ t f ] , where ≤ n ≤ N . Let N be the total number of steps such that t f = τ N where τ denotes the time step. If the step-size restriction < τ ≤ / b holds, then the solution of (2.16) satisfies properties C – C at times t n = nτ , ≤ n ≤ N .Proof. Consider the system (2.16) at an arbitrary step n and assume that the properties C – C hold at the n − C is satisfiedby (2.16c). Moreover, by assumption S n − , I n − , and R n − are non-negative and hence bydefinition F n − is also non-negative. As a result, e − τ ( F n − + c ) <
1, and therefore S n is non-negative and monotonically decreasing. Similarly, the right hand side terms of (2.16b) are alsonon-negative, thus R n is non-negative and monotonically increasing.To show that I n is non-negative, we substitute (2.16a) and (2.16b) into (2.16c) to get I n = S n − (cid:16) − ( + cτ ) e − τ ( F n − + c ) (cid:17) + I n − ( − bτ ) .Since by assumption S n − and I n − are non-negative, a sufficient condition for I n to be non-negative is 1 − ( + cτ ) e − τ ( F n − + c ) ≥ − bτ ≥ x − ln ( + x ) ≥ x > −
1. Since c ≥ F n − is non-negativewe then have τ F n − + cτ − ln ( + cτ ) ≥
0. Rearranging the inequality givesln (cid:18) + cτ (cid:19) ≥ − τ ( F n + c ) ,hence 1 − ( + cτ ) e − τ ( F n + c ) ≥ τ >
0. As a result, a sufficient condition for I n toremain non-negative is 0 < τ ≤ / b .Note that by using the same arguments as above we can show that conditions C – C hold atthe first step, i.e., n =
1, provided that the initial conditions are non-negative. This completesthe proof.
Remark 2.1.
Using left Riemann sums to approximate the integrals in (2.15) results in localerrors of order O ( τ ) . Therefore, the solution of (2.15) can only be first order accurate.In the next two sections, we discretize (1.4) by first using a numerical approximation of theintegral on the right hand side of the system, and then applying a time integration method.This approach results in numerical schemes that are high order accurate, both in space andtime. It is evident that the key point of the numerical solution of problem (1.4) is the approximationof F ( t , x , y ) . This can be done in two different ways. The first approach is to approximate thefunction I ( t , ¯ x ( r , θ ) , ¯ y ( r , θ )) by a Taylor expansion, and then proceed further. This method is12tudied in [12] and [13], but is not efficient in the case of non-constant function g ( θ ) as shownin [35]. The other approach is to use a combination of interpolation and numerical integration(by using cubature formulas) to obtain an approximation of F ( t , x , y ) .We consider two-dimensional cubature formulas on the disc of radius δ with positive coeffi-cients. Denote by Q ( x , y ) the set of cubature points in the disk B δ (cid:16) x , y (cid:17) parametrized by polarcoordinates (see [35]), i.e., Q ( x , y ) : = { ( x ij , y ij ) = ( x + r i cos ( θ j ) , y + r i sin ( θ j )) ∈ B δ (cid:16) x , y (cid:17) , i ∈ I , j ∈ J } ,where r i denotes the distance from center point ( x , y ) , θ j is the angle, and I and J are the setof indices of cubature points. Using numerical integration, we get the system ∂S ( t , x , y ) ∂t = − S ( t , x , y ) T ( t , Q ( x , y )) − cS ( t , x , y ) , ∂I ( t , x , y ) ∂t = S ( t , x , y ) T ( t , Q ( x , y )) − bI ( t , x , y ) , ∂R ( t , x , y ) ∂t = bI ( t , x , y ) + cS ( t , x , y ) , (3.1)where T ( t , Q ( x , y )) = X ( x ij , y ij ) ∈Q ( x , y ) w i , j g ( r i ) g ( θ j ) I ( t , x + r i cos ( θ j ) , y + r i sin ( θ j )) ,and w i , j > Remark 3.1.
Note that Theorem 2.4 can be applied to system (3.1); hence, the properties C – C hold without any restrictions for the solution of this system. Moreover, it can be easilyshown that T ( t , Q ( x , y )) satisfies properties ( A ) and ( A ), by following the proofs of Lemma 2.2and Lemma 2.3. As a result system (3.1) admits a unique strong solution. Now we would like to solve (3.1) numerically. The first step is to discretize the problem inspace. Let us suppose that we would like to solve our problem on a rectangle-shaped domain,namely Ω : = [ L ] × [ L ] . For our numerical solutions we will discretize this domain byusing a spatial grid G : = { ( x k , y l ) ∈ Ω | ≤ k ≤ P , 1 ≤ l ≤ P } ,which consists of P × P points with spatial step sizes h and h , and approximate the contin-uous solutions by a vector of the values at the gridpoints.After this semi-discretization, we get the following set of equations dS k , l ( t ) dt = − S k , l ( t ) T k , l ( t , Q ( x k , y l )) − cS k , l ( t ) , dI k , l ( t ) dt = S k , l ( t ) T k , l ( t , Q ( x k , y l )) − bI k , l ( t ) , dR k , l ( t ) dt = bI k , l ( t ) + cS k , l ( t ) , (3.2)13here X k , l ( t ) , X ∈ { S , I , R } , denotes the approximation of the function at gridpoint ( x k , y l ) .The approximation of F ( · , x k , y l ) is denoted by T k , l ( t , Q ( x k , y l )) and defined as T k , l ( t , Q ( x k , y l )) : = X ( ¯ x k , ¯ y l ) ∈Q ( x k , y l ) w i , j g ( r i ) g ( θ j ) ˜ I ( t , ¯ x k , ¯ y l ) , (3.3)where ¯ x k = x k + r i cos ( θ j ) and ¯ y l = y l + r i sin ( θ j ) . Note that the points ( ¯ x k , ¯ y l ) might not beincluded in G ; in such case there are no I k , l values assigned to them. Because of this, we ap-proximate I ( t , ¯ x k , ¯ y l ) by a using positivity preserving interpolation (e.g. bilinear interpolation)with the nearest known I k , l values and positive coefficients. This is the reason why ˜ I is used in(3.3) instead of I . Theorem 3.1.
A unique strong solution for system (3.2) exists, for which properties C – C hold locally at a given point ( x k , y l ) .Proof. The proof of existence and uniqueness comes from the Lipschitz continuity and boundnessof the right hand side, which can be proved similarly as in Corollary 2.1 and Lemma 2.4.Properties C – C can be proved in a similar manner as in Theorem 2.4.The next theorem characterizes the accuracy of interpolation and cubature techniques ofsystem (3.2). Theorem 3.2.
Suppose that a cubature rule approximates the integral (2.5) to order p , i.e., kF ( I ( t , x , y )) − T ( t , Q ( x , y )) k L = O ( δ p ) , (3.4) where δ is the radius of the disk in which the integration takes place. Let us suppose that the(positivity preserving) spatial interpolation ˜ I approximates the values of I to order q , i.e., (cid:13)(cid:13) I ( t , x , y ) − ˜ I ( t , x , y ) (cid:13)(cid:13) L = O ( h q ) (3.5) where h is the spatial step size. Then if ˜ u is the solution of (1.4) evaluated at the gridpoints of G and ˜ v is the solution of (3.2) , it follows that k ˜ u − ˜ v k L = O ( δ p ) + O ( h q ) . Proof.
Let us proceed in the following way: prove that if w is the solution of (3.1) evaluated atthe gridpoints of G , then k ˜ u − w k L = O ( δ p ) and k w − ˜ v k L = O ( h q ) L .It is easy to see that the theorem follows from the above statement.We prove both estimates similarly to Lemma 2.5, namely by applying the formula of constantvariations. Hence, the solutions can be expressed as ˜ u ( t ) = S ( t ) ˜ u + Z t S ( t − s ) F ( s ) ds , (3.6) w ( t ) = S ( t ) ˜ u + Z t S ( t − s ) ( F ( s ) + O ( δ p )) ds , (3.7) ˜ v ( t ) = S ( t ) ˜ u + Z t S ( t − s ) F T ( s ) ds , (3.8)14here {S ( t ) , t ∈ [ ∞ ) } is the analytic semigroup associated with − A in (2.3) [14, p. 101], and F ( s ) : = F ( u ( s )) . We also use the notation F T ( s ) : = F T ( I ( s , · , · )) , where F T is the operatorthat maps ˜ I ( s , ., . ) to T G ( s ) : = ( T k , l ( t , Q ( x k , y l ))) ( x k , y l ) ∈G , i.e., F T ( ˜ I ( s , · , · )) = T G ( s ) , and weuse the assumption (3.4). Then it follows that k ˜ u − w k L = (cid:13)(cid:13)(cid:13)(cid:13)Z t S ( t − s ) O ( δ p ) ds (cid:13)(cid:13)(cid:13)(cid:13) L = O ( δ p ) .For the second one we use the assumption (3.5), and consequently the fact that w ( t ) = S ( t ) ˜ u + Z t S ( t − s )( F T ( s ) + O ( h q )) ds .Thus, k w − ˜ v k L = (cid:13)(cid:13)(cid:13)(cid:13)Z t S ( t − s ) O ( h q ) ds (cid:13)(cid:13)(cid:13)(cid:13) L = O ( h q ) ,which gives the second estimates, and consequently the prove of the theorem is complete.A natural question arises: what is the best type of cubature and interpolation for solvingthe system (3.2)? In the rest of the section we describe two numerical integration proceduresand also discuss suitable interpolation techniques. One can use a direct cubature rule on the general disk, see for example [9, 34]. In such case theintegral of a function f ( x , y ) over the disk with radius δ can be approximated by Q ( f ) = πδ N r · N θ X i = w i f ( x i , y j ) = πδ N r X i = N θ X j = e w i f ( r i cos ( θ j ) , r i sin ( θ j )) , (3.9)where N r is the number of radial points, N θ is the number of equally spaced angles, and w i and e w i are weights in the [
0, 1 ] interval. We use N θ = N r to have a cubature rule that is equallypowerful in both r and θ . The weights and cubature points are calculated by a modificationof the Elhay-Kautsky Legendre quadrature method [11, 23, 28]. The top panel of Figure 3.1shows the distribution of cubature abscissas for N r ∈ {
3, 6, 13 } . The Elhay-Kautsky cubatureresults in points that are evenly spaced in the θ direction. Alternatively, we can transform the disk into a square, and then use a one-dimensional Gauss-Legendre rule to approximate the integral. First, we transform the disk with radius δ to therectangle [ δ ] × [
0, 2 π ] in the r − θ plane. Next, the rectangle [ δ ] × [
0, 2 π ] is mapped to [
0, 1 ] × [
0, 1 ] on the ξ − η plane by using the linear transformation r = δξ , θ = πη ,15hat has a Jacobian 2 πδ . Using these transformations, the original integral Z δ Z π f ( r cos ( θ ) , r sin ( θ )) rdθdr takes the form Z Z f ( δξ cos ( πη ) , δξ sin ( πη )) δξ πδ dη dξ . (3.10)There are several approaches for computing multiple integrals based on numerical integrationof one-dimensional integrals. In this paper, we use the Gauss-Legendre quadrature rule on theunit interval [37]; other options include generalized Gaussian quadrature rules as described in[27]. The integral (3.10) can be approximated by Q ( f ) = N ξ X i = N η X j = w i w j πδ ξ i f ( δξ i cos ( πη j ) , δξ i sin ( πη j )) = N ξ · N η X m = e w m f ( x m , y m ) , (3.11)where ξ i and η i are the i th abscissas corresponding to the Gauss-Legendre quadrature withweights w i . The number of cubature points in the ξ and η direction are denoted by N ξ and N η , respectively, and we let x m = δξ i cos ( πη j ) , y m = δξ i sin ( πη j ) and e w m = w i w j πδ ξ i .The distribution of the abscissas in the unit disk is not uniform as with the Elhay-Kautskycubature and can be seen in the bottom panel of Figure 3.1. For a fair comparison we use N η = N ξ . Experimental results reveal that the Elhay-Kautsky cubature (3.9) performs betterin cases the interpolated function f ( x , y ) is a bivariate polynomial, whereas the Gauss-Legendrequadrature (3.11) or the generalized Gaussian quadrature rule (see [27]) when f ( x , y ) is anarbitrary nonlinear function.In order to determine which cubature rule performs better for the system (3.2), we per-form a convergence test by applying the cubature formulas (3.9) and (3.11) to the function g ( r ) g ( θ ) I r , where g ( r ) = ( − r + δ ) , g ( θ ) = sin ( θ ) ,and I ( r , θ ) = πσ exp − r s ! (3.12)is a Gaussian distribution with deviation σ and centered at zero. This resembles the initialconditions for I at the origin, as we will see later in Section 5. The exact solution of the integralover a disk of radius δ is given by Z δ Z π g ( r ) g ( θ ) I r dθdr = (cid:18) δ − √ π σ erf (cid:18) δ √ σ (cid:19)(cid:19) .where erf ( x ) is the Gauss error function [3, 21]. Figure 3.2 shows the convergence of the twocubature rules over the disk of radius δ , when σ = /
10, as δ goes to zero. We observe that16he Gauss-Legendre quadrature (3.11) gives much smaller errors (close to machine precision)when more than 13 ×
26 points are used, compared to the Elhay-Kautsky cubature (3.9) whichis third-order accurate.The performance of the cubature formulas depends also on the choice and accuracy of in-terpolation. As mentioned before, bilinear interpolation can be used since it preserves thenon-negativity of the interpolant. One possibility is to use higher order interpolations, like cu-bic or spline, but in these cases the preservation of the required properties cannot be guaranteed.However, numerical experiments show that piecewise cubic spline interpolation results in a pos-itive interpolant for sufficiently fine spatial grid. A better choice is the use of a shape-preservinginterpolation, to ensure that negative values are not generated and the interpolant of I ( t , ¯ x k , ¯ y l ) is bounded by max k , l { S k , l + I k , l + R k , l } for every point x k , y l . This can be accomplished by amonotone interpolation that uses piecewise cubic Hermite interpolating polynomials [10, 15]. In MATLAB the relevant function is called pchip but is only available for one-dimensional problems.Extensions to bivariate shape-preserving interpolation have been studied in [7, 8, 16]; how- -1 -0.5 0 0.5 1-1-0.500.51 (a) 3 × -1 -0.5 0 0.5 1-1-0.500.51 (b) 6 ×
12 points -1 -0.5 0 0.5 1-1-0.500.51 (c) 13 ×
26 points -1 -0.5 0 0.5 1-1-0.500.51 (d) 3 × -1 -0.5 0 0.5 1-1-0.500.51 (e) 6 ×
12 points -1 -0.5 0 0.5 1-1-0.500.51 (f) 13 ×
26 points
Figure 3.1:
Top panel : The distribution of abscissas ( N r × N θ ) in the unit disk using the Elhay-Kautsky cubature rule, bottom panel : The distribution of abscissas ( N x i × N η ) in the unit diskusing the Gauss-Legendre quadrature rule. 17 .0031 0.0063 0.0125 0.025 0.05 0.1 0.210 -5 (a) Elhay-Kautsky cubature (3.9) -20 -15 -10 -5 (b) Gauss-Legendre quadrature (3.11) Figure 3.2: Errors by using different number of points for cubature formulas (3.9) and (3.11).ever, this topic goes beyond the purposes of this paper. Another choice is the modified Akimapiecewise cubic Hermite interpolation , makima . Numerical experiments demonstrate good per-formance as it avoids overshoots when more than two consecutive nodes are constant [1, 2], andhence preserves non-negativity in areas where I ( t , ¯ x k , ¯ y l ) is close to zero. The next step is to use time integration methods to solve the system of ordinary differentialequations (3.2). First we study sufficient and necessary time-step restrictions such that theforward Euler method satisfies a discrete analogue of properties C – C , denoted below by D – D . Then, we discuss how high order SSP Runge–Kutta methods can be applied to (3.2).Let X n = { X nk , l } , X ∈ { S , I , R } , be the numerical approximation of X k , l ( t n ) for all 1 ≤ k ≤ P , 1 ≤ l ≤ P , and 0 ≤ n ≤ N , where N is the total number of steps. The numerical solutionshould satisfy the following properties: D : The densities { X nk , l } , X ∈ { S , I , R } , are non-negative for every 1 ≤ k ≤ P , 1 ≤ l ≤ P ,and for all 0 ≤ n ≤ N . D : The sum S nk , l + I nk , l + R nk , l is constant for all 0 ≤ n ≤ N and for every 1 ≤ k ≤ P ,1 ≤ l ≤ P . D : The density S nk , l is non-increasing, i.e., S nk , l ≤ S n − k , l for every 1 ≤ k ≤ P , 1 ≤ l ≤ P , andfor all 1 ≤ n ≤ N . D : The density R nk , l is increasing i.e., R nk , l ≥ R n − k , l for every 1 ≤ k ≤ P , 1 ≤ l ≤ P , and forall 1 ≤ n ≤ N . 18 .1 Explicit Euler scheme and qualitative properties Let us apply the explicit Euler method to the system (3.2) on the interval [ t f ] , and choosean adaptive time step τ n > t n = t n − + τ n , n ≥
1. After the full discretization weget the set of algebraic equations S n = S n − − τ n S n − ◦ T n − − cτ n S n − , I n = I n − + τ n S n − ◦ T n − − bτ n I n − , R n = R n − + bτ n I n − + cτ n S n − , (4.1a)(4.1b)(4.1c)Here, the operator ◦ denotes the element-by-element or Hadamard product of these matrices.Now we examine the bounds of time step τ n such that the method (4.1) gives solutionswhich are qualitatively adequate and satisfy conditions D – D . Theorem 4.1.
Consider the numerical solution (4.1) obtained by forward Euler method appliedto (3.2) with non-negative initial data. Then, the solution satisfies property D without any step-size restrictions. Moreover, properties D , D and D hold if the time step satisfies τ n ≤ min ( min k , l T n − k , l + c , 1 b ) , (4.2) where T n − k , l = X ( ¯ x k , ¯ y l ) ∈Q ( x k , y l ) w i , j g ( r i ) g ( θ j ) ˜ I n − ( ¯ x k , ¯ y l ) (4.3) is an approximation of (3.3) at point ( x k , y l ) ∈ G .Proof. The proof is similar to the one of [35, Theorem 2]. We prove the statement by inductionon the number of steps.First, assume that the properties D – D hold up to step n −
1; we will prove that they alsohold true for step n . Property D can be easily verified by adding all equations in (4.1). Toshow the monotonicity and non-negativity of S n , consider (4.1a) at point ( x k , y l ) ∈ G S nk , l = ( − τ n ( T n − k , l + c )) S n − k , l .By our assumption I n − k , l ≥
0, and a positivity-preserving interpolation guarantees that the inter-polated values ˜ I n − ( ¯ x k , ¯ y l ) = ˜ I n − ( x k + r i cos ( θ j ) , y l + r i sin ( θ j )) are non-negative. Therefore,by (4.3) we get T n − k , l ≥ ≤ k ≤ P , 1 ≤ l ≤ P since the weights w i , j are positive,and functions g and g are non-negative. As a result, τ n ( T n − k , l + c ) ≥ S nk , l ≤ S n − k , l .Moreover, if τ n ≤ / ( T n − k , l + c ) then S nk , l remains non-negative. Equation (4.1b) yields I nk , l = ( − bτ n ) I n − k , l + τ n S n − k , l T n − k , l ,and hence I n is non-negative if τ n ≤ / b . Finally from (4.1c) we have R nk , l = R n − k , l + bτ n I n − k , l + cτ n S n − k , l ,19herefore R n is non-negative R n ≥ R n − . Putting all together we conclude that properties D – D are satisfied if the time step is bounded by (4.2). By using the above argument it canbe shown that D – D also hold at the first step, n =
1, if the initial data are non-negative andthe time step satisfies (4.2).A drawback of time-step restriction (4.2) is that it depends on the solution at the previousstep. This has important complications for higher order methods as we will see in Section 4.2.For any multistage method the adaptive time step bound (4.2) depends not only on the previoussolution, but also on the internal stage approximations. Consequently, an adaptive time-steprestriction based on (4.2) cannot be the same for all stages of a Runge–Kutta method; insteadit needs to be recalculated at every stage to guarantee that conditions D – D hold. Therefore,such bound has no practical use because it is prone to rejected steps and will likely tend to zero.A remedy is to use a constant time step that is less strict than (4.2), but still guarantee that τ ≤ / ( c + T n − k , l ) holds for all 1 ≤ k ≤ P , 1 ≤ l ≤ P and at every step n . We construct amatrix similar to T k , l from (4.3), namely e T : = X ( ¯ x k , ¯ y l ) ∈Q ( x k , y l ) w i , j g ( r i ) g ( θ j ) M , (4.4)where M is a matrix with all entries equal to e m = max ( x k , y l ) ∈G { S ( x k , y l ) + I ( x k , y l ) + R ( x k , y l ) } .It is evident that if τ = min (cid:26) e T + c , 1 b (cid:27) , (4.5)then the condition τ ≤ T nk , l + c still holds. Moreover, e T ≤ e w κ e mN , where κ = max { max r ∈ ( δ ) { g ( r ) } , max θ ∈ [ π ) { g ( θ ) }} , e w = max ij ( w i , j ) , and N is the number of the points in Q ( x , y ) . Hence, the time step (4.5) islarger than the rather pessimistic time step ˜ τ : = min (cid:26) e w κ e mN + c , 1 b (cid:27) , (4.6)proposed in [35, Theorem 2]. Numerical experiments show that (4.5) is pretty close to thetheoretical bound (4.2), and thus a relatively small increase of time step (4.5) may producequalitatively bad solutions which violate one of the conditions D – D (see Section 5.1). Forward Euler method is only first-order accurate; hence, we would like to obtain time-steprestrictions for higher order Runge–Kutta methods. Note that the spatial discretizations dis-cussed in Section 3 can be chosen so that errors from cubature formulas and interpolation arevery small; therefore, it is substantial to have a high order method in time.20onsider a Runge–Kutta method in the Butcher form [5] with coefficients ( a ij ) ∈ R m × m and b ∈ R m . Let K be the matrix given by K = " ( a ij ) b (cid:124) ,and denote by I the ( m + ) -dimensional identity matrix. If there exists r > ( I + r K ) is invertible, then the Runge–Kutta method can be expressed in the canonical Shu–Osher form Q ( i ) = v i Q n − + m X j = α ij (cid:18) Q ( j ) + τr F (cid:16) Q ( j ) (cid:17)(cid:19) , 1 ≤ i ≤ m + Q n = Q ( m + ) , (4.7)where the coefficient arrays ( α ij ) and ( v i ) have non-negative components. Such methods arecalled strong-stability preserving (SSP) Runge–Kutta methods and have been introduced byShu as total-variation diminishing (TVD) discretizations [31], and by Shu and Osher in relationto high order spatial discretizations [32, 33]. The choice of parameter r gives rise to differentShu–Osher representations; thus we denote the Shu–Osher coefficients of (4.7) by α r = ( α ij ) and v r = ( v i ) to emphasize the dependence on the parameter r . The Shu–Osher representationwith the largest value of r such that ( I + r K ) − exists and α r , v r have non-negative componentsis called optimal and attains the SSP coefficient C = max { r ≥ : ∃ ( I + r K ) − and α r ≥ v r ≥ } .The interested reader may consult [17, 19, 20], as well as the monograph [18] and the referenceswithin, for a throughout review of SSP methods.We would like to investigate time-step restrictions such that the numerical solution obtainedby applying method (4.7) to the problem (3.2) satisfies properties D – D . The following theoremprovides the theoretical upper bound for the time step such that these properties are satisfied. Theorem 4.2.
Consider the numerical solution obtained by applying an explicit Runge–Kuttamethod (4.7) with SSP coefficient C > to the semi-discrete problem (3.2) with non-negativeinitial data. Then property D holds without any time-step restrictions. Moreover, the properties D , D and D hold if the time step satisfies τ ≤ C min (cid:26) e T + c , 1 b (cid:27) , (4.8) where e T is given by (4.4) .Proof. Consider an arbitrary stage i , 1 ≤ i ≤ m +
1, of a Runge–Kutta method (4.7) with21on-negative coefficients and SSP coefficient C >
0. Applying the method to (3.2) we get S ( i ) = v i S n − + i − X j = α ij (cid:18) S ( j ) − τ C (cid:16) S ( j ) ◦ T ( j ) − cS ( j ) (cid:17)(cid:19) , (4.9a) I ( i ) = v i I n − + i − X j = α ij (cid:18) I ( j ) + τ C (cid:16) S ( j ) ◦ T ( j ) − bI ( j ) (cid:17)(cid:19) , (4.9b) R ( i ) = v i R n − + i − X j = α ij (cid:18) R ( j ) + τ C (cid:16) bI ( j ) + cS ( j ) (cid:17)(cid:19) . (4.9c)Since all Runge–Kutta methods preserve linear invariants the property D , i.e., S n + I n + R n = S n − + I n − + R n − , ∀ n is trivially satisfied.The remainder of the proof deals with properties D , D and D . We show that all quantities S n , I n , R n remain non-negative, while S n is non-increasing and R n is increasing. From (4.9a)and (4.9b) we have, respectively, S ( i ) = v i S n − + i − X j = α ij S ( j ) ◦ (cid:18) − τ C (cid:16) T ( j ) + c (cid:17)(cid:19) , I ( i ) = v i I n − + τr i − X j = α ij S ( j ) ◦ T ( j ) + (cid:18) − τr b (cid:19) i − X j = α ij I ( j ) ,where is the P × P all-ones matrix.By definition, T ( i ) k , l = X ( ¯ x k , ¯ y l ) ∈Q ( x k , y l ) w i , j g ( r i ) g ( θ j ) ˜ I ( i ) ( ¯ x k , ¯ y l ) , 1 ≤ i ≤ m + ˜ I ( i ) are interpolated values. Since the initial data are non-negative and the choseninterpolation is positivity-preserving, we have that S ( ) = S n − , I ( ) = I n − and T ( ) are allnon-negative. If 0 ≤ − τr b , and 0 ≤ − τ C (cid:16) T ( j ) + c (cid:17) for 1 ≤ j ≤ i −
1, (4.10)then the explicit Runge–Kutta method inductively results in non-negative T ( i ) , S ( i ) , and I ( i ) for each 2 ≤ i ≤ m +
1. Note that by (4.4) it holds that T ( i ) ≤ e T , for all 1 ≤ i ≤ m +
1, (4.11)because ˜ I ( i ) is an interpolated value of I ( i ) and hence bounded by M .Moreover, the non-negativity of T ( i ) implies that − τ C (cid:16) T ( i ) + c (cid:17) ≤ ≤ i ≤ m + S ( i ) ≤ v i S n − + i − s X j = α ij S ( j ) .Consistency requires that v i + P i − j = α ij = ≤ i ≤ m + S ( i ) ≤ ( − i − X j = α ij ) S n − + i − X j = α ij S ( j ) ≤ S n − i − X j = α ij (cid:16) S n − − S ( j ) (cid:17) . (4.12)Let 1 ≤ q ≤ m + S ( i ) ≤ S ( q ) for all 1 ≤ i ≤ m +
1. Then, taking i = q in (4.12) yields S ( q ) ≤ v i S n − + i − X j = α qj S ( q ) − i − X j = α qj S ( q ) ≤ − i − X j = α qj S n − S ( q ) ≤ S n − .Therefore, S ( i ) ≤ S n − for all 1 ≤ i ≤ m +
1. In particular for i = m + S n = S ( m + ) ≤ S n − .Finally, the non-negativity of initial data, S ( j ) and I ( j ) imply that from (4.9c) we have R ( i ) ≥ R n − for all 1 ≤ i ≤ m +
1, and hence R n = R ( m + ) ≥ R n − .Combining (4.10) and (4.11) we conclude that the step-size restriction (4.8) is sufficient forsatisfying properties D – D . In this section we confirm the results proved in the previous sections by using several numericalexperiments. Computational tests are defined in a bounded domain and thus the choice ofboundary conditions is important. Because we have no diffusion in our problem, we considerhomogeneous Dirichlet conditions and we assume that there is no susceptible population outsideof our domain. This means that we are going to assign a zero value to any point which liesoutside of our rectangle in which the problem is defined. In most cases the abscissas of thecubatures rules (3.9) and (3.11) do not belong to the spatial grid. Special attention must begiven to the corners and boundaries of the domain where cubature points assigned to gridpointsnear the boundary lie outside of the domain. To be able to handle them, we are using ghostcells which are set to zero. This enables us to calculate the values corresponding to the cubaturepoints lying outside of the domain. 23or the numerical experiments we are choosing the following functions. Let g ( r ) be alinearly decreasing function, which takes its maximum at r = r = δ , i.e., g ( r ) : = a ( − r + δ ) ,where a is a parameter as in (1.1). Also, we are going to use a non-constant g ( θ ) function,which is going to be symmetrical in a way that g ( θ ) : = β sin ( θ + α ) + β .From now on, we are using the choices of α = β =
1, that is, assuming a northern windon the domain. In our numerical experiments we are also using the parameter values: a = b = c = δ = S k , l = S =
20 for all k , l , and there are no recovered at the beginning,i.e. R k , l = k , l . For the infected species, we use a Gaussian distribution concentrated atthe middle point ( L / L / ) of the domain, and has standard deviation σ = min ( L , L ) / I k , l = πσ exp − h ( k − ) − L σ + h ( l − ) − L σ ,where L = ( P − ) h and L = ( P − ) h as mentioned in (3.12). We also set L = L = t f =
80. As we can see, the number of susceptibles is decreased, andthe number of infected moves towards the boundaries, while forming a wave. Both densities S and I tend to the zero function, which confirms that the zero solution is indeed an asymptoticallystable equilibrium for the first two equations of (1.4), as it was proved in Section 2.In the next section we demonstrate how close the bound of Theorem 4.1 is to the real bound,i.e., the largest time step under which the numerical scheme preserves the desired qualitativeproperties. As we saw in Section 4.1, the bound (4.5) is larger than ˜ τ in (4.6). Thus step-size restriction(4.5) is closer to the best theoretically reasonable bound (4.2) that guarantees the preservationof properties D – D . A natural question is how close is the bound (4.5) to the adaptive step-size restriction, when compared with the pessimistic bound (4.6). In Table 5.1 we have testedseveral different values of a , for which both the bounds ˜ τ from (4.6) and τ from (4.5) werecomputed, and also the best possible time step choice (i.e. the largest possible time step underthe properties D – D hold), denoted by τ e , is calculated by numerical experiments using abisection method. As we can see, the use of the bound (4.5) results in an increase of about 40%.In Table 5.2, a similar experiment was performed, but with varying δ instead of parameter a .In this case, although the bounds differ considerably, their ratios remain similar, resulting in anapproximate 40% difference between the improved bound (4.5) and (4.6).24igure 5.1: The number of susceptibles S (left), infected I (middle) and recovered R (right)at time t f =
80. The SSPRK104 method has been used with the Gauss-Legendre quadrature(3.11) (15 ×
15 points) and spline interpolation. We used 30 grid points in each direction and δ = τ , the methods still behave as expected. However, for values of τ bigger than (4.5), there is noguarantee that properties D – D will be satisfied. Figure 5.2 shows that the values of S becomenegative when a bigger time step than (4.5) is used. Since we cannot compute the exact solution accurately, we are going to compute the errorsusing a reference solution. The reference solution is computed by using the same parametersand method, but with a very small step size.We first observe how well the different cubatures behave. As seen in Section 3, using more a ˜ τ ˜ τ / τ e τ τ / τ e τ e
50 2.4340 0.4870 4.2000 0.8403 4.9980100 1.2320 0.5171 2.1450 0.9003 2.3826150 0.8247 0.5273 1.4403 0.9209 1.5640200 0.6198 0.5348 1.0841 0.9354 1.1590
Table 5.1: The bounds for the forward Euler method for different values of a , with choices P = P = b = c = δ = -point nonuniform quadrature, bilinearinterpolation and final time t f =
80. 25 ˜ τ ˜ τ / τ e τ τ / τ e τ e Table 5.2: The bounds for the forward Euler method for different values of δ , with choices P = P = a = b = c = ×
20 Gauss-Legendre quadraturerule and bilinear interpolation.Figure 5.2: The values of S for parameters a = b = c = δ = ×
40 grid, 20 ×
20 Gauss-Legendre quadrature rule, and bilinear interpolation with timesteps τ = τ = L -norm errors can beseen in Figure 5.3. It is clear that for 25 or 100 points there is no remarkable difference betweenthe interpolations, but for more cubature points cubic and spline interpolation perform better.Bilinear interpolation results in similar errors for both cubatures (3.9) and (3.11). As it can beseen, cubic and spline interpolation perform the same way for cubature (3.9) and smaller errorsare observed with spline interpolation and cubature (3.11).Another important question is the order of the different time integration methods, and howthe forward Euler compares to the first-order direct method described in Section 2.2. Numericalexperiments show that the higher order schemes work as expected, namely that by using enoughcubature points and grid points, a reasonably small error can be achieved with the desiredaccuracy order.Table 5.3 shows the convergence rates for various SSP Runge–Kutta methods and the integral26
10 15 20 25 30 3510 -7 -6 -5 -4 -3 -2 -1 bilinear (uniform)cubic (uniform)spline (uniform)bilinear (nonuniform)cubic (nonuniform)spline (nonuniform) Figure 5.3: L -norm errors using cubatures formulas (3.9) and (3.11) with n points, n ∈{
5, 10, 15, 20, 25, 30, 35 } and different interpolations. The final time is chosen to be t f = h = /
30 and we used the ten-stage, fourth-order SSP Runge–Kutta method(SSPRK104).solution used in Section 2.2. We start with a reasonable time step 0.429, which is below thebound (4.5) and then successively divide by 2. For the reference solution we use a time step thatis the half of the smallest time step in our computations. It is evident that using higher ordermethods is better than solving the integral equation (2.14) numerically, and even the forwardEuler method produces smaller errors than the scheme (2.16).
In this paper the SIR model for epidemic propagation is extended to include spatial dependence.The existence and uniqueness of the continuous solution were proved, along with properties cor-responding to biological observations. For the numerical solution, different choices of cubature,interpolation and time integration methods are studied. It is shown that for a sufficient smalltime-step restriction, the numerical solution preserves a discrete analogue of the properties ofthe original continuous system. This bound was also made sharper, and an adaptive step-sizetechnique is also suggested for the explicit Euler method in Section 4.1. These theorems wereconfirmed by numerical experiments, while the errors of cubature formulas and the order ofaccuracy of the time discretization methods are also discussed.A possible extension of the work presented in this paper is to study diffusion spatial-dependent SIR systems, including the effect of fractional diffusion. Results for the preservationof qualitative properties can be potentially obtained. Moreover the inclusion the births andnatural deaths in the system and dropping the conservation property makes the model morerealistic. It would be interesting to study the influence of such modification in the behavior ofthe continuous and also the numerical solution.27
FE IM0.4000 9.4 × − × − × − × − × − × − × − × − × − × − τ FE SSPRK22 SSPRK33 SSPRK1040.4290 9.9 × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − Table 5.3: L -norm errors and convergence orders of different time integration methods (“IM"denotes the method (2.16)), with final time t f =
80, spacial grid cells P = P =
20, bilinearinterpolation and the cubature rule (3.11) with 10 ×
10 points.
A Proofs of Lemmata in Section 2
In this section we present the proofs of some technical lemmata that were omitted in the previoussections.
Proof of Lemma 2.1.
The proof is the following straightforward calculation: k u k = c c k u k + b b k u k ≤ max (cid:18) c , 1 b (cid:19) ( c k u k + b k u k )= max (cid:18) c , 1 b (cid:19) k u k A ,and k u k A = c k u k + b k u k ≤ max ( c , b )( k u k + k u k ) = max ( c , b ) k u k . Proof of Lemma 2.2.
The term we are going to give an upper bound to is kF ( I ) k L = Z Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Z δ Z π g ( r ) g ( θ ) I ( t , ¯ x ( r , θ ) , ¯ y ( r , θ )) r dθdr (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dxdy = Z Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Z B δ ( x ) g ( r ) g ( θ ) I ( t , ˜ x ) d ˜ x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d x ,where we used the notation x : = ( x , y ) and B δ ( x ) is the ball with radius δ around x . By thedefinition of g and g : kF ( I ) k L = Z Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Z Ω g ( r ) g ( θ ) I ( t , ˜ x ) d ˜ x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d x .28y definition, we know that g and g are bounded. Let us use the notations κ = max r ∈ ( δ ) { g ( r ) } and κ = max θ ∈ ( π ) { g ( θ ) } . Then, kF ( I ) k L ≤ κ κ Z Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Z Ω I ( t , ˜ x ) d ˜ x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d x = κ κ Z Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Z Ω · I ( t , ˜ x ) d ˜ x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d x ≤ κ κ Z Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)sZ Ω d ˜ x sZ Ω ( I ( t , ˜ x )) d ˜ x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d x ≤ κ κ µ ( Ω ) Z Ω Z Ω | I ( t , ˜ x ) | d ˜ x d x ,where we used the Cauchy–Schwarz inequality, and µ ( Ω ) is the Lebesgue measure of Ω . Itholds that Z Ω Z Ω | I ( t , ˜ x ) | d ˜ x d x = Z Ω k I k L d x = µ ( Ω ) k I k L .Consequently, kF ( I ) k L ≤ κ κ µ ( Ω ) k I k L ,from which we get that ν F = κ κ µ ( Ω ) . Proof of Lemma 2.3.
We would like to bound the following expression: kF ( I ) − F ( I ) k L = Z Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Z δ Z π g ( r ) g ( θ ) ( I ( t , ¯ x ( r , θ ) , ¯ y ( r , θ )) − I ( t , ¯ x ( r , θ ) , ¯ y ( r , θ ))) r dθdr (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dxdy .We can proceed similarly as in the proof of Lemma 2.2: kF ( I ) − F ( I ) k L = Z Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Z B δ ( x ) g ( r ) g ( θ ) ( I ( t , ˜ x ) − I ( t , ˜ x )) d ˜ x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d x = Z Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Z Ω g ( r ) g ( θ ) ( I ( t , ˜ x ) − I ( t , ˜ x )) d ˜ x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d x ≤ κ κ Z Ω (cid:12)(cid:12)(cid:12)(cid:12)Z Ω ( I ( t , ˜ x ) − I ( t , ˜ x )) d ˜ x (cid:12)(cid:12)(cid:12)(cid:12) d x ≤ κ κ µ ( Ω ) Z Ω Z Ω | I ( t , ˜ x ) − I ( t , ˜ x ) | d ˜ x d x = κ κ µ ( Ω ) Z Ω k I − I k L d x = κ κ µ ( Ω ) k I − I k L .which completes the proof with C F = κ κ µ ( Ω ) .29 roof of Lemma 2.5. The proof uses the method of variation of constants. Consider the non-homogeneous semilinear equation u ( t ) = Au ( t ) + F ( u ( t )) , (A.1)where A is a linear bounded operator and F is Hölder continuous. Then, the solution corre-sponds of the solution of the following integral equation: u ( t ) = S ( t ) u + Z t S ( t − s ) F ( s ) ds , (A.2)where {S ( t ) , t ∈ [ ∞ ) } is the analytic semigroup associated with the infinitesimal generator A [14, p. 101], and we use the notation F ( s ) : = F ( u ( s )) .For the system (1.4) we use similar choices as in the beginning of this section, namely A isgiven by (2.3) and F ε u u ! = F u u ! + ε ! ,where F is given by (2.4) and ε (cid:28)
1. Note that we do not consider the third equation of (1.4),since it can be omitted as noted before in Section 2.It is clear that A generates an analytical semigroup, and we also know that F is Lipschitz-continuous because of Lemma 2.4; hence, the method of variation of constants is applicable.Consequently, if u ε and u ε are solutions of (A.2), then k u ε − u ε k = (cid:13)(cid:13)(cid:13)(cid:13) S ( t ) u + Z t S ( t − s ) F ε ( s ) ds − S ( t ) u − Z t S ( t − s ) F ε ( s ) ds (cid:13)(cid:13)(cid:13)(cid:13) .Let ˜ ε i = ( ε i ) (cid:124) , i =
1, 2, then k u ε − u ε k = (cid:13)(cid:13)(cid:13)(cid:13)Z t S ( t − s )( F ( s ) + ˜ ε − F ( s ) − ˜ ε ) ds (cid:13)(cid:13)(cid:13)(cid:13) = (cid:13)(cid:13)(cid:13)(cid:13) ( ˜ ε − ˜ ε ) Z t S ( t − s ) ds (cid:13)(cid:13)(cid:13)(cid:13) .As both ε and ε tend to zero, then the above expression tends to zero too, and this completesthe proof. Acknowledgments
The project has been supported by the European Union, co-financed by the European SocialFund (EFOP-3.6.3-VEKOP-16-2017-00002). The authors would like to thank Lajos Lóczi forhis overall support, and Inmaculada Higueras and David Ketcheson for their comments.
References [1] Akima, H., A New Method of Interpolation and Smooth Curve Fitting Based on LocalProcedures,
Communications of the ACM , 17 (1970), pp. 589–602.302] Akima, H., A Method of Bivariate Interpolation and Smooth Surface Fitting Based onLocal Procedures,
Communications of the ACM , 17 (1974), pp. 18–20.[3] Andrews, L. C., Special functions of mathematics for engineers, second, SPIE OpticalEngineering Press, Bellingham, WA; Oxford University Press, Oxford, 1998, pp. xx+480.[4] Aronson, D. G., The asymptotic speed of propagation of a simple epidemic,
Nonlin-ear diffusion (NSF-CBMS Regional Conf. Nonlinear Diffusion Equations, Univ. Houston,Houston, Tex., 1976) , vol. 14, Research Notes in Mathematics, Pitman London, 1977,pp. 1–23.[5] Butcher, J. C., Numerical methods for ordinary differential equations, Third, With aforeword by J. M. Sanz-Serna, John Wiley & Sons, Ltd., Chichester, 2016, pp. xxiii+513.[6] Capasso, V. and Fortunato, D., Stability results for semilinear evolution equations andtheir application to some reaction-diffusion problems,
SIAM Journal on Applied Mathe-matics , 39 (1980), pp. 37–47.[7] Carlson, R. E. and Fritsch, F. N., Monotone piecewise bicubic interpolation,
SIAM Jour-nal on Numerical Analysis , 22 (1985), pp. 386–400.[8] Carlson, R. E. and Fritsch, F. N., An algorithm for monotone piecewise bicubic interpo-lation,
SIAM Journal on Numerical Analysis , 26 (1989), pp. 230–238.[9] Davis, P. J. and Rabinowitz, P., Methods of numerical integration, Corrected reprint ofthe second (1984) edition, Dover Publications, Inc., Mineola, NY, 2007, pp. xii+612.[10] Dougherty, R. L., Edelman, A. S., and Hyman, J. M., Nonnegativity-, monotonicity-, orconvexity-preserving cubic and quintic Hermite interpolation,
Mathematics of Computa-tion , 52 (1989), pp. 471–494.[11] Elhay, S. and Kautsky, J., Algorithm 655: IQPACK: FORTRAN subroutines for theweights of interpolatory quadratures,
Association for Computing Machinery. Transactionson Mathematical Software , 13 (1987), pp. 399–415.[12] Faragó, I. and Horváth, R., On some qualitatively adequate discrete space-time modelsof epidemic propagation,
Journal of Computational and Applied Mathematics , 293 (2016),pp. 45–54.[13] Faragó, I. and Horváth, R., Qualitative properties of some discrete models of diseasepropagation,
Journal of Computational and Applied Mathematics , 340 (2018), pp. 486–500.[14] Friedman, A., Partial differential equations, Holt, Rinehart and Winston, Inc., NewYork-Montreal, Que.-London, 1969, pp. vi+262.[15] Fritsch, F. N. and Carlson, R. E., Monotone piecewise cubic interpolation,
SIAM Journalon Numerical Analysis , 17 (1980), pp. 238–246.[16] Fritsch, F. N. and Carlson, R. E., Monotonicity preserving bicubic interpolation: aprogress report, vol. 2, Surfaces in CAGD ’84 (Oberwolfach, 1984), 1985, pp. 117–121.[17] Gottlieb, S., Ketcheson, D. I., and Shu, C.-W., High order strong stability preservingtime discretizations,
Journal of Scientific Computing , 38 (2009), pp. 251–289.3118] Gottlieb, S., Ketcheson, D. I., and Shu, C.-W., Strong stability preserving Runge-Kuttaand multistep time discretizations, World Scientific Publishing Co. Pte. Ltd., Hackensack,NJ, 2011, pp. xii+176.[19] Gottlieb, S. and Shu, C.-W., Total variation diminishing Runge-Kutta schemes,
Mathe-matics of Computation , 67 (1998), pp. 73–85.[20] Gottlieb, S., Shu, C.-W., and Tadmor, E., Strong stability-preserving high-order timediscretization methods,
SIAM Review , 43 (2001), pp. 89–112.[21] Greene, W. H., Econometric Analysis, Prentice Hall, 2002.[22] Harko, T., Lobo, F. S. N., and Mak, M. K., Exact analytical solutions of the susceptible-infected-recovered (SIR) epidemic model and of the SIR model with equal death and birthrates,
Applied Mathematics and Computation , 236 (2014), pp. 184–194.[23] Kautsky, J. and Elhay, S., Calculation of the weights of interpolatory quadratures,
Nu-merische Mathematik , 40 (1982), pp. 407–422.[24] Kendall, D. G., in discussion with Bartlett, M. S.: Measles periodicity and communitysize,
Journal of the Royal Statistical Society. Series A , 120 (1957), in discussion withBartlett, Maurice S., pp. 48–70.[25] Kendall, D. G., Mathematical models of the spread of infection,
Mathematics andcomputer science in biology and medicine , (1965), pp. 213–225.[26] Kermack, W. O. and McKendrick, A. G., A contribution to the mathematical theory ofepidemics,
Proceedings of the Royal Society of London. Series A , 115 (1927), pp. 700–721.[27] Ma, J.-H., Rokhlin, V. V., and Wandzura, S. M., Generalized Gaussian quadraturerules for systems of arbitrary functions,
SIAM Journal on Numerical Analysis , 33 (1996),pp. 971–996.[28] Martin, R. S. and Wilkinson, J. H., Handbook Series Linear Algebra: The implicit QL algorithm, Numerische Mathematik , 12 (1968), pp. 377–383.[29] Miller, J. C., A note on the derivation of epidemic final sizes,
Bulletin of MathematicalBiology , 74 (2012), pp. 2125–2141.[30] Miller, J. C., Mathematical models of SIR disease spread with combined non-sexual andsexual transmission routes,
Infectious Disease Modelling , 2 (2017), pp. 35–55.[31] Shu, C.-W., Total-variation-diminishing time discretizations,
SIAM Journal on Scientificand Statistical Computing , 9 (1988), pp. 1073–1084.[32] Shu, C.-W., Essentially non-oscillatory and weighted essentially non-oscillatory schemesfor hyperbolic conservation laws,
Advanced numerical approximation of nonlinear hyper-bolic equations , vol. 1697, Lecture Notes in Mathematics, Berlin: Springer, 1998, pp. 325–432.[33] Shu, C.-W. and Osher, S., Efficient implementation of essentially nonoscillatory shock-capturing schemes,
Journal of Computational Physics , 77 (1988), pp. 439–471.3234] Stroud, A. H., Approximate calculation of multiple integrals, Prentice-Hall Series inAutomatic Computation, Prentice-Hall, Inc., Englewood Cliffs, N.J., 1971, pp. xiii+431.[35] Takács, B., Horváth, R., and Faragó, I., Space dependent models for studying the spreadof some diseases,
Computers and Mathematics with Applications , (2019), (in press).[36] Thieme, H. R., A model for the spatial spread of an epidemic,
Journal of MathematicalBiology , 4 (1977), pp. 337–351.[37] Trefethen, L. N., Is Gauss quadrature better than Clenshaw-Curtis?,