Trapeioidal rule for a multiple integral over a hyperquadrilateral is devised. The N-dimensional hyperquadrilateral is partitioned into 2(N)N! hypertriangles, in each of which the integrand is interpolated by a linear function of the N coordinate variables. The resulting N-dimensional trapezoidal rule is a useful quadrature formula for approximating a multiple integral by a weighted sum of the values of the integrand at the nodes of a hyperquadrilateral lattice. For multiple integrals, likely the trapezoidal rule gives better approximation than higher-degree rules when the dimensionality is high. The latter interpolatory rules have the shortcoming that some of the basis polynomials may be not everywhere nonnegative, incurring the possibility of rendering the associated nodal weights negative. As an application, the trapezoidal rule is applied to a surface integral to obtain finite-sum expressions for partial derivatives of a function of three variables in non-orthogonal coordinates. The obtained finite-sum approximation has better accuracy than corresponding finite-difference approximation by accounting for the coupling effect of multiple dimensionality. (C) Elsevier Science Inc., 1997.