  Boltzmann solvers are an important tool for the computation of cosmological
observables in the linear regime. They involve solving the Boltzmann equation,
followed by an integration in momentum space, to arrive at the desired fluid
properties. This is a cumbersome, computationally expensive procedure. In this
work we introduce the so-called generalized Boltzmann hierarchy (GBH) for
massive neutrinos in cosmology, a simpler alternative to the usual Boltzmann
hierarchy, where the momentum dependence is integrated out leaving us with a
two-parameter infinite set of ordinary differential equations. Along with the
usual expansion in multipoles, there is now also an expansion in higher
velocity weight integrals of the distribution function. We show that the GBH
produces the density contrast neutrino transfer function to a per mille level
accuracy at both large and intermediate scales compared to the neutrino
free-streaming scale. Furthermore, by introducing a switch to a viscous fluid
approximation after horizon crossing, we show that the GBH can achieve over all
scales the same accuracy as the standard CLASS approach in its default
precision settings. The GBH is then a powerful tool to include neutrino
anisotropies in the computation of cosmological observables in linear theory,
with integration being simpler and potentially faster than standard methods.
