An interval arithmetic method is described for finding the global maxima or minima of multivariable functions. The original domain of variables is divided successively, and the lower and the upper bounds of the interval expression of the function are estimated on each subregion. By discarding subregions where the global solution can not exist, one can always find the solution with rigorous error bounds. The convergence can be made fast by Newton's method after subregions are grouped. Further, constrained optimization can be treated using a special transformation or the Lagrange-multiplier technique.