<p>The generalized wave equation extends traditional elastic wave theory by incorporating couple-stress and strain-gradient effects, thereby providing a more comprehensive description of subsurface mechanical behavior. This formulation is particularly suitable for characterizing rotational motion and scale-dependent phenomena induced by microstructural heterogeneity. In this study, the generalized wave equation with couple-stress and strain-gradient terms is derived and numerically solved using a finite-difference (FD) method because of its simplicity, computational efficiency, and suitability for parallel implementation. However, FD methods often suffer from numerical dispersion. To address this issue, we develop optimized staggered-grid finite-difference coefficients using an improved self-adaptive differential evolution (SADE) algorithm. The proposed algorithm employs a multi-strategy adaptive mutation mechanism and a dynamically updated strategy pool to balance global exploration and local exploitation. The resulting optimized coefficients for first-order spatial derivatives achieve high accuracy and significantly reduce numerical dispersion. Numerical experiments on homogeneous and Marmousi models demonstrate that the proposed method reproduces wavefield phase and amplitude characteristics more accurately than conventional FD coefficients and coefficients optimized using other metaheuristic algorithms. The high-precision finite-difference coefficients derived in this study provide effective support for imaging and inversion based on the generalized wave equation.</p>