We present the Boltzmann moment equation approach for the dynamics of stars (beads-2d), which is a finite-difference Eulerian numerical code designed for the modelling of anisotropic and non-axisymmetric flat stellar discs. The beads-2d code solves the Boltzmann moment equations up to second order in the thin-disc approximation. This allows us to obtain the anisotropy of the velocity ellipsoid and the vertex deviation in the plane of the disc. We study the time-dependent evolution of exponential stellar discs in the linear regime and beyond. The discs are initially characterized by different values of the Toomre parameter Qs and are embedded in a dark matter halo, yielding a rotation curve composed of a rigid central part and a flat outer region. Starting from a near equilibrium state, several unstable modes develop in the disc. In the early linear phase, the very centre and the large scales are characterized by growing one-armed and bisymmetric positive density perturbations, respectively. This is in agreement with expectations from the swing amplification mechanism of short-wavelength trailing disturbances, propagating through the disc centre. In the late linear phase, the overall appearance is dominated by a two-armed spiral structure localized within the outer Lindblad resonance (OLR). During the non-linear evolutionary phase, radial mass redistribution due to the gravitational torques of spiral arms produces an outflow of mass, which forms a ring at the OLR, and an inflow of mass, which forms a transient central bar. This process of mass redistribution is self-regulatory and it terminates when spiral arms diminish due to a shortage of matter. Finally, a compact central disc and a diffuse ring at the OLR are formed. An increase in Qs stabilizes the discs at Qs≈ 3.15, in agreement with the theoretical predictions. Considerable vertex deviations are found in regions with strongly perturbed mass distributions, that is, near the spiral arms. The vertex deviations are especially large at the convex edge of the spiral arms, whereas they are small at the concave edge. The mean vertex deviations correlate well with the global Fourier amplitudes, reaching mean values of about 12° in the saturation stage. Local values of the vertex deviation can reach up to almost 90°. Near the convex edge of the spiral arms, the ratio of radial to azimuthal components of the velocity ellipsoid can deviate considerably from the values predicted from the epicycle approximation.