Summary: This paper develops a stabilization procedure in piezoelectric media to ensure the temporal stability of node-based smoothed finite element method (NS-FEM), and applies it to obtain certified solution with bounds in both static and frequency analyses of piezoelectric structures using three-node triangular elements. For such stabilized NS-FEM, two stabilization terms corresponding to squared-residuals of two equilibrium equations, i.e., mechanical stress equilibrium and electric displacement equilibrium, are added into the smoothed potential energy functional of the original NS-FEM. A gradient smoothing operation is then performed on second-order derivatives of shape functions to achieve the stabilization terms. Due to the use of divergence theory, the smoothing operation relaxes the requirement of shape functions, so that the square-residuals can be evaluated using linear elements. The effectiveness of the present stabilized NS-FEM is demonstrated via numerical examples.