A modified higher-order shear deformation plate theory (MHSDT) is proposed to predict the resonance frequencies and related mode shapes of piezoelectric disks of various thicknesses. The theory explicitly accounts for the effects of the electrical potential on the displacement of the disk. The equations of motion for the disk are derived using Hamilton's principle in conjunction with a variational approach. It is shown analytically that the distribution of the electrical field on the surface of the disk is quantitatively similar to the distribution of the sum of the total stresses within the disk. Thus, the deformation of the disk can be approximated directly from the stress distribution on the disk, thereby greatly simplifying the solution process. The validity of the proposed MHSDT is confirmed by comparing the solutions for the natural frequencies, displacement fields, and electric field distributions of the disk for various radius-to-thickness ratios and boundary conditions with those obtained using Reddy’s third-order plate theory and COMSOL finite element simulations, respectively. The results show that the solutions obtained from the MHSDT method are consistently closer to the exact solutions than those determined using Reddy’s third-order method. Finally, the linear relationship between the electric field and the flexural displacement of the disk is leveraged to design surface electrodes capable of exciting specific flexural mode shapes of the disk under a free boundary condition. The efficiency of the designed electrodes is experimentally demonstrated using an amplitude-fluctuation electronic speckle pattern interferometer (AF-ESPI).