The theory of adiabatic invariants has a long history and important applications in physics but is rarely rigorous. Here we treat exactly the general time-dependent 1-D harmonic oscillator, $\ddot{q} + \omega^2(t) q=0$ which cannot be solved in general. We follow the time-evolution of an initial ensemble of phase points with sharply defined energy $E_0$ and calculate rigorously the distribution of energy $E_1$ after time $T$, and all its moments, especially its average value $\bar{E_1}$ and variance $\mu^2$. Using our exact WKB-theory to all orders we get the exact result for the leading asymptotic behaviour of $\mu^2$.
