We present several methods, which utilize symplectic integration techniques based on two and three part operator splitting, for numerically solving the equations of motion of the disordered, discrete nonlinear Schrödinger (DDNLS) equation, and compare their efficiency. Our results suggest that the most suitable methods for the very long time integration of this one-dimensional Hamiltonian lattice model with many degrees of freedom (of the order of a few hundreds) are the ones based on three part splits of the system’s Hamiltonian. Two part split techniques can be preferred for relatively small lattices having up to N ≈ 70 sites. An advantage of the latter methods is the better conservation of the system’s second integral, i.e. the wave packet’s norm.