In this work, we study the long time behaviors, including asymptotic contractivity and dissipativity, of the solutions to several numerical methods for fractional ordinary differential equations (F-ODEs). The existing algebraic contractivity and dissipativity rates of the solutions to the scalar F-ODEs are first improved. In order to study the long time behavior of numerical solutions to fractional backward differential formulas (F-BDFs), two crucial analytical techniques are developed, with the first one for the discrete version of the fractional generalization of the traditional Leibniz rule, and the other for the algebraic decay rate of the solution to a linear Volterra difference equation. By mens of these auxiliary tools and some natural conditions, the solutions to F-BDFs are shown to be contractive and dissipative, and also preserve the exact contractivity rate of the continuous solutions. Two typical F-BDFs, based on the Grünwald-Letnikov formula and L1 method respectively, are studied. For high order F-BDFs, including some second order F-BDFs and 3-α order method, their numerical contractivity and dissipativity are also developed under some slightly stronger conditions. Numerical experiments are presented to validate the long time qualitative characteristics of the solutions to F-BDFs, revealing very different decay rates of the numerical solutions in terms of the the initial values between F-ODEs and integer ODEs and demonstrating the superiority of the structure-preserving numerical methods. .hk).1 As noted in [22], a numerical method preserves the contractivity in (4) is sometimes referred to be dissipative in the numerical literature, but this conflicts with the corresponding terminology in dynamical systems. At the same time, the authors in [22] give an accurate definition of dissipativity for ODEs, which mainly emphasizes the existence of a global attracting set. This definition was later widely accepted in the numerical literature, and we also follow this definition.