Sparse grids provide efficient approximations of smooth functions. More precisely, they are a device to describe a function up to a prescribed accuracy with very few degrees of freedom. If conventional Sobolev-norms are used to measure accuracy, sparse grids can be shown to be optimal or near optimal. A tremendous reduction of the amount of data is achieved compared to standard schemes of multivariate approximations that rely on low-order polynomials. It turns out that sparse grids are a priori adaptive grids. The crucial idea is to drop certain insignificant contributions of hierarchical representations of functions. So far, sparse grids schemes have been based on linear and higher order Lagrange polynomials for H1 conforming finite elements. We generalize the discretization on sparse grids to discrete differential forms (Whitney forms). The extension to general l-forms in d dimensions includes mixed elements for and can be used for the Galerkin discretization of mixed boundary value problems. The tensor product structure and the hierarchical multilevel principle provide the hierarchical decomposition of the Whitney spaces. We define the sparse grid interpolation operator relying on the hierarchical basis. The interpolation estimates generalize the known results for Lagrangian finite elements. Approximate interpolation is needed for the Galerkin method for boundary value problems on sparse grids. The combination technique and a two point quadrature rule ensure that similar error estimate as for the exact interpolation hold. Some inherent limitations of the proposed method are given by the fact that only tensor product domains have been successfully investigated. Further, no techniques for boundary value problems with variable coefficients are available yet. However, we will discuss adaptively refined sparse grids, too.