A direct variational method based on an energy principle is applied to obtain approximate magnetohydrodynamic equilibria for tokamak plasmas. The geometry of the nested magnetic flux surfaces is specified by a model that includes displacement, elongation and triangularity effects. The radial dependence in flux coordinates is described by a consistent set of trial functions which allows analytical calculation of the flux-surface averaged internal energy density of the plasma. Approximate solutions of the variational problem are obtained for arbitrary aspect-ratio tokamaks using a one-parameter optimization procedure.