A variational approach to the analysis of the natural decay rates and eigenmodes of cavity-enclosed diffusive fields in general anisotropic and heterogeneous materials is presented. In the bulk material, diffusivity and volume relaxivity are accounted for. The interaction of the cavity's medium with the embedding material is modeled via a surface relaxivity on the boundary surface. The pertaining eigenmodes are proven to be orthogonal and to form a complete expansion of an initially excited diffusive field. In view of the variational approach, a finite-element type of computation presents itself as the natural tool for numerics. The resulting implementation on a simplicial mesh allows for the modeling of cavities of arbitrary shape. To investigate the feasibility of using the approach in the inverse problem of reconstructing the shape and size of cavities from measured values of the natural decay rates of the eigenmodes, we carry out a number of numerical experiments on the forward problem. They demonstrate the method to be simple and robust, both in 2D and 3D complex geometries. For a benchmark problem with a known analytic solution, error estimates are presented. Applications are found in, for example, nuclear magnetic resonance imaging of subsurface rock pore geometry, biological cell structure and the analysis of neurological defects in medical diagnostics.