  The asymptotic stable region and long-time decay rate of solutions to linear
homogeneous Caputo time fractional ordinary differential equations (F-ODEs) are
known to be completely determined by the eigenvalues of the coefficient matrix.
Very different from the exponential decay of solutions to classical ODEs,
solutions of F-ODEs decay only polynomially, leading to the so-called
Mittag-Leffler stability, which was already extended to semi-linear F-ODEs with
small perturbations. This work is mainly devoted to the qualitative analysis of
the long-time behavior of numerical solutions. By applying the singularity
analysis of generating functions developed by Flajolet and Odlyzko (SIAM J.
Disc. Math. 3 (1990), 216-240), we are able to prove that both $\mathcal{L}$1
scheme and strong $A$-stable fractional linear multistep methods (F-LMMs) can
preserve the numerical Mittag-Leffler stability for linear homogeneous F-ODEs
exactly as in the continuous case. Through an improved estimate of the discrete
fractional resolvent operator, we show that strong $A$-stable F-LMMs are also
Mittag-Leffler stable for semi-linear F-ODEs under small perturbations. For the
numerical schemes based on $\alpha$-difference approximation to Caputo
derivative, we establish the Mittag-Leffler stability for semi-linear problems
by making use of properties of the Poisson transformation and the decay rate of
the continuous fractional resolvent operator. Numerical experiments are
presented for several typical time fractional evolutional equations, including
time fractional sub-diffusion equations, fractional linear system and
semi-linear F-ODEs. All the numerical results exhibit the typical long-time
polynomial decay rate, which is fully consistent with our theoretical
predictions.
