Variational theory based on self-adjoint equations of motion cannot fully represent the interaction of the earth's seismic free oscillations in the presence of lateral structure, attenuation, and rotation. The more general Galerkin procedure can model correctly the frequencies and attenuation rates of hybrid oscillations. Implementation of either algorithm leads to a generalized matrix eigenvalue problem in which the potential and kinetic energy interactions are separated into distinct matrices. The interaction of the earth's seismic free oscillations due to aspherical structure, attenuation, and rotation is best treated as a matrix eigenvalue problem. The presence of attenuation causes the matrices to be non-Hermitian and requires the use of a general Galerkin procedure. Physical dispersion, represented as a logarithmic function in frequency, must be represented by a truncated Taylor series about a fiducial frequency in order to be incorporated in the Galerkin formalism in a numerically tractable manner. The earth's rotation introduces an interaction matrix distinct from the potential and kinetic energy matrices, leading to a quadratic eigenvalue problem. A simple approximation leads to an eigenvalue problem linear in squared frequency. Tests show that this approximation is accurate for calculations using modes of frequencies f>1 mHz, unless interaction across a wide frequency band is modeled. Hybrid oscillation particle motions are represented by matrix eigenvectors that can be significatnly nonorthogonal. The degrees of freedom in the low-frequency seismic system remain distinct, since source excitation is calculated by using dual eigenvectors. Synthetic seismograms that are constructed from Galerkin coupling calculations without reference to this eigenvector nonorthogonality can be disastrously noncausal.