Two fully discrete methods are investigated for simulating the distributed-order sub-diffusion equation in Caputo's form. The fractional Caputo derivative is approximated by the Caputo's BDF1 (called L1 early) and BDF2 (or L1-2 when it was first introduced) approximations, which are constructed by piecewise linear and quadratic interpolating polynomials, respectively. It is shown that the first scheme, using the BDF1 formula, possesses the discrete minimum-maximum principle and nonnegativity preservation property such that it is stable and convergent in the maximum norm. The method using the BDF2 formula is shown to be stable and convergent in the discrete H (1) norm by using the discrete energy method. For problems of distributed order within a certain region, the method is also proven to preserve the discrete maximum principle and nonnegativity property. Extensive numerical experiments are provided to show the effectiveness of numerical schemes, and to examine the initial singularity of the solution. The applicability of our numerical algorithms to a problem with solution which lacks the smoothness near the initial time is examined by employing a class of power-type nonuniform meshes.