  Gaussian beams are asymptotically valid high frequency solutions to
hyperbolic partial differential equations, concentrated on a single curve
through the physical domain. They can also be extended to some dispersive wave
equations, such as the Schr\"odinger equation. Superpositions of Gaussian beams
provide a powerful tool to generate more general high frequency solutions that
are not necessarily concentrated on a single curve. This work is concerned with
the accuracy of Gaussian beam superpositions in terms of the wavelength
$\epsilon$. We present a systematic construction of Gaussian beam
superpositions for all strictly hyperbolic and Schr\"odinger equations subject
to highly oscillatory initial data of the form $Ae^{i\Phi/\epsilon}$. Through a
careful estimate of an oscillatory integral operator, we prove that the $k$-th
order Gaussian beam superposition converges to the original wave field at a
rate proportional to $\epsilon^{k/2}$ in the appropriate norm dictated by the
well-posedness estimate. In particular, we prove that the Gaussian beam
superposition converges at this rate for the acoustic wave equation in the
standard, $\epsilon$-scaled, energy norm and for the Schr\"odinger equation in
the $L^2$ norm. The obtained results are valid for any number of spatial
dimensions and are unaffected by the presence of caustics. We present a
numerical study of convergence for the constant coefficient acoustic wave
equation in $\Real^2$ to analyze the sharpness of the theoretical results.
