  A class of evolution equations with nonlocal diffusion is considered in this
work. These are integro-differential equations arising as models of propagation
phenomena in continuum media with nonlocal interactions including neural
tissue, porous media flow, peridynamics, models with fractional diffusion, as
well as continuum limits of interacting dynamical systems. The principal
challenge of numerical integration of nonlocal systems stems from the lack of
spatial regularity of the data and solutions intrinsic to nonlocal models. To
overcome this problem we propose a semidiscrete numerical scheme based on the
combination of sparse Monte Carlo and discontinuous Galerkin methods. An
important feature of our method is sparsity. Sparse sampling of points in the
Monte Carlo approximation of the nonlocal term allows to use fewer
discretization points without compromising the accuracy.
  We prove convergence of the numerical method and estimate the rate of
convergence. There are two principal ingredients in the error of the numerical
method related to the use of Monte Calro and Galerkin approximations
respectively. We analyze both errors. Two representative examples of
discontinuous kernels are presented. The first example features a kernel with a
singularity, while the kernel in the second example experiences jump
discontinuity. We show how the information about the singularity in the former
case and the geometry of the discontinuity set in the latter translate into the
rate of convergence of the numerical procedure. In addition, we illustrate the
rate of convergence estimate with a numerical example of an initial value
problem, for which an explicit analytic solution is available. Numerical
results are consistent with analytical estimates.
