Motivated by fractional derivative models in viscoelasticity, a class of semilinear stochastic Volterra integro-differential equations, and their deterministic counterparts, are considered. A generalized exponential Euler method, named here as the Mittag-Leffler Euler integrator, is used for the temporal discretization, while the spatial discretization is performed by the spectral Galerkin method. The temporal rate of strong convergence is found to be (almost) twice compared to when the backward Euler method is used together with a convolution quadrature for time discretization. Numerical experiments that validate the theory are presented.
