Numerical computation of shape gradients from Eulerian derivatives is essential to wildly used gradient type methods in shape optimization. Boundary type Eulerian derivatives are popularly used in literature. The volume type Eulerian derivatives hold more generally, but are rarely noticed and used numerically. We investigate thoroughly the accuracy of Galerkin finite element approximations of the two type shape gradients for optimization of elliptic eigenvalues. Under certain regularity assumptions on domains, we show \emph{a priori} error estimates for the two approximate shape gradients. The convergence analysis shows that the volume integral formula converges faster and generally offers better accuracy. Numerical experiments verify theoretical results for the Dirichlet case. For the Neumann case, however, the boundary formulation surprisingly converges as fast as the volume one. Numerical results are presented.