A stable numerical method for integral epidemic models with behavioral changes in contact patterns

We propose a non-standard numerical method for the solution of a system of integro-differential equations describing an epidemic of an infectious disease with behavioral changes in contact patterns. The method is constructed in order to preserve the key characteristics of the model, like the positivity of solutions, the existence of equilibria, and asymptotic behavior. We prove that the numerical solution converges to the exact solution as the step size h of the discretization tends to zero. Furthermore, the method is first-order accurate, meaning that the error in the discretization is O(h), it is linearly implicit, and it preserves all the properties of the continuous problem, unconditionally with respect to h. Numerical simulations show all these properties and confirm, also by means of a case-study, that the method provides correct qualitative information at a low computational cost.