A notable feature of the anomalous sub-solution equation is its solution’s algebraic decay over extended time periods, a phenomenon commonly associated with Mittag-Leffler type stability. For power-nonlinear sub-diffusion models with variable coefficients, we prove Mittag-Leffler stability under natural decay assumptions on the source functions, with decay rate \(\Vert u(t)\Vert _{L^{s}(\Omega )}=O( t^{-(\alpha +\beta )/\gamma } )\) as \(t\rightarrow \infty \) , where \(\alpha \) , \(\gamma \) are positive constants, \(\beta \in (-\alpha ,\infty )\) and \(s\in (1,\infty )\) . We then develop a structure-preserving algorithm for these models. For the complete monotonicity-preserving ( \(\mathcal{C}\mathcal{M}\) -preserving) schemes developed by Li and Wang (Commun. Math. Sci., 19(5):1301-1336, 2021), we show they satisfy the discrete comparison principle for time-fractional differential equations with variable coefficients. By carefully constructing the fine the discrete supersolutions and subsolutions, we obtain the numerical solution’s long-time optimal decay rate \(\Vert u_{n}\Vert _{L^{s}(\Omega )}=O( t_n^{-(\alpha +\beta )/\gamma } )\) as \(t_{n}\rightarrow \infty \) , which aligns perfectly with the theoretical decay rate. Finally, we validate our analysis through numerical experiments.
