We investigate a variably distributed-order time-fractional wave partial differential equation, which could accurately model, e.g., the viscoelastic behavior in vibrations in complex surroundings with uncertainties or strong heterogeneity in the data. A standard composite rectangle formula of mesh size $\sigma$ is firstly used to discretize the variably distributed-order integral and then the L-1 formula of degree of freedom $N$ is applied for the resulting fractional derivatives. Optimal error estimates of the corresponding fully-discrete finite element method are proved based only on the smoothness assumptions of the data.
To maintain the accuracy, setting $\sigma=N^{-1}$ leads to $O(N^3)$ operations of evaluating the temporal discretization coefficients and $O(N^2)$ memory. To improve the computational efficiency, we develop a novel time-stepping scheme by expanding the fractional kernel at a fixed fractional order to decouple the fractional operator from the variably distributed-order integral. Only $O(\log N)$ terms are needed for the expansion without loss of accuracy, which consequently reduce the computational cost of coefficients from $O(N^3)$ to $O(N^2\log N)$ and the corresponding memory from $O(N^2)$ to $O(N\log N)$. Optimal-order error estimates of this time-stepping scheme are rigorously proved via novel and different techniques from the standard analysis procedure of the L-1 methods. Numerical experiments are presented to substantiate the theoretical results.