  In contrast with the diffusion equation which smoothens the initial data to
$C^\infty$ for $t>0$ (away from the corners/edges of the domain), the
subdiffusion equation only exhibits limited spatial regularity. As a result,
one generally cannot expect high-order accuracy in space in solving the
subdiffusion equation with nonsmooth initial data. In this paper, a new
splitting of the solution is constructed for high-order finite element
approximations to the subdiffusion equation with nonsmooth initial data. The
method is constructed by splitting the solution into two parts, i.e., a
time-dependent smooth part and a time-independent nonsmooth part, and then
approximating the two parts via different strategies. The time-dependent smooth
part is approximated by using high-order finite element method in space and
convolution quadrature in time, while the steady nonsmooth part could be
approximated by using smaller mesh size or other methods that could yield
high-order accuracy. Several examples are presented to show how to accurately
approximate the steady nonsmooth part, including piecewise smooth initial data,
Dirac--Delta point initial data, and Dirac measure concentrated on an
interface. The argument could be directly extended to subdiffusion equations
with nonsmooth source data. Extensive numerical experiments are presented to
support the theoretical analysis and to illustrate the performance of the
proposed high-order splitting finite element methods.
