The characteristics of the partial differential equations describing three‐dimensional toroidal magnetohydrodynamic equilibria with nested flux surfaces in inverse flux coordinates are derived and examined. The equilibrium equations are then variationally reduced to a truncated set of ordinary differential equations by decomposing the flux surface geometry into a spectral representation. The magnetic field lines on the flux surfaces are given in terms of a variable stream function to allow optimum choice of the angle coordinates over the flux surfaces and to simplify the treatment in the vicinity of a rational magnetic surface. Analytic properties of the spectral representation and moment equations are considered. Comparative calculations are performed numerically. The results agree well with those calculated using a standard three‐dimensional equilibrium code, but the variational spectral method is substantially faster computationally. The Mercier stability criterion is given.