  Solutions to many important partial differential equations satisfy bounds
constraints, but approximations computed by finite element or finite difference
methods typically fail to respect the same conditions. Chang and Nakshatrala
enforce such bounds in finite element methods through the solution of
variational inequalities rather than linear variational problems. Here, we
provide a theoretical justification for this method, including higher-order
discretizations. We prove an abstract best approximation result for the linear
variational inequality and estimates showing that bounds-constrained
polynomials provide comparable approximation power to standard spaces. For any
unconstrained approximation to a function, there exists a constrained
approximation which is comparable in the $W^{1,p}$ norm. In practice, one
cannot efficiently represent and manipulate the entire family of
bounds-constrained polynomials, but applying bounds constraints to the
coefficients of a polynomial in the Bernstein basis guarantees those
constraints on the polynomial. Although our theoretical results do not
guaruntee high accuracy for this subset of bounds-constrained polynomials,
numerical results indicate optimal orders of accuracy for smooth solutions and
sharp resolution of features in convection-diffusion problems, all subject to
bounds constraints.
