  We study quench dynamics of the Bose-Hubbard model by exact diagonalization.
Initially the system is at thermal equilibrium and of a finite temperature. The
system is then quenched by changing the on-site interaction strength $U$
suddenly. Both the single-quench and double-quench scenarios are considered. In
the former case, the time-averaged density matrix and the real-time evolution
are investigated. It is found that though the system thermalizes only in a very
narrow range of the quenched value of $U$, it does equilibrate or relax well in
a much larger range. Most importantly, it is proven that this is guaranteed for
some typical observables in the thermodynamic limit. In order to test whether
it is possible to distinguish the unitarily evolving density matrix from the
time-averaged (thus time-independent), fully decoherenced density matrix, a
second quench is considered. It turns out that the answer is affirmative or
negative according to the intermediate value of $U$ is zero or not.
