The shape derivative of Laplacian eigenvalues with respect to domain deformations was theoretically investigated by Hadamard in the early 20th century. However, the rigorous computation of these derivatives is not an easy task, since there exists the singularity for repeated eigenvalues. In this study, we propose a verified computation method for the shape derivative of Laplacian eigenvalues using guaranteed computation of both eigenvalues and eigenfunctions.