  In this paper, we will show the $L^p$-resolvent estimate for the finite
element approximation of the Stokes operator for $p \in \left( \frac{2N}{N+2},
\frac{2N}{N-2} \right)$, where $N \ge 2$ is the dimension of the domain. It is
expected that this estimate can be applied to error estimates for finite
element approximation of the non-stationary Navier--Stokes equations, since
studies in this direction are successful in numerical analysis of nonlinear
parabolic equations. To derive the resolvent estimate, we introduce the
solution of the Stokes resolvent problem with a discrete external force. We
then obtain local energy error estimate according to a novel localization
technique and establish global $L^p$-type error estimates. The restriction for
$p$ is caused by the treatment of lower-order terms appearing in the local
energy error estimate. Our result may be a breakthrough in the $L^p$-theory of
finite element methods for the non-stationary Navier--Stokes equations.
