  We present new formalism for description of the neutrino oscillations in
matter with varying density. The formalism is based on the Magnus expansion and
has a virtue that the unitarity of the S-matrix is maintained in each order of
perturbation theory. We show that the Magnus expansion provides better
convergence of series: the restoration of unitarity leads to smaller deviations
from the exact results especially in the regions of large transition
probabilities. Various expansions are obtained depending on a basis of neutrino
states and a way one splits the Hamiltonian into the self-commuting and
non-commuting parts. In particular, we develop the Magnus expansion for the
adiabatic perturbation theory which gives the best approximation. We apply the
formalism to the neutrino oscillations in matter of the Earth and show that for
the solar oscillation parameters the second order Magnus adiabatic expansion
has better than 1% accuracy for all energies and trajectories. For the
atmospheric $\Delta m^2$ and small 1-3 mixing the approximation works well ($<
3 %$ accuracy for $\sin^2 \theta_{13} = 0.01$) outside the resonance region
(2.7 - 8) GeV.
