  We consider three types of subdiffusion models, namely single-term,
multi-term and distributed order fractional diffusion equations, for which the
maximum-principle holds and which, in particular, preserve nonnegativity. Hence
the solution is nonnegative for nonnegative initial data. Following earlier
work on the heat equation, our purpose is to study whether this property is
inherited by certain spatially semidiscrete and fully discrete piecewise linear
finite element methods, including the standard Galerkin method, the lumped mass
method and the finite volume element method. It is shown that, as for the heat
equation, when the mass matrix is nondiagonal, nonnegativity is not preserved
for small time or time-step, but may reappear after a positivity threshold. For
the lumped mass method nonnegativity is preserved if and only if the
triangulation in the finite element space is of Delaunay type. Numerical
experiments illustrate and complement the theoretical results.
