According to Nyquist-Shannon sampling theory, two sampling points of the displacements per minimal wavelength are required to reconstruct a wavefield. However, high-frequency components above the Nyquist frequency can be theoretically recovered by introducing spatial derivatives of the displacements. We have developed a new finite-difference scheme for solving the acoustic equation, named the dispersion-relation preserving stereo-modeling method, which is practical for seismic forward modeling beyond the Nyquist frequency. The key idea in this model is to approximate the high-order spatial derivatives with wavenumber-domain optimization in wavefields represented by the wave displacements and their gradients. We investigated the theoretical properties of this method, including the absolute error of spatial operators, numerical error, stability criterion, and computational efficiency. Our results indicated that the method can effectively suppress the numerical dispersion on extremely coarse grids with less than two grid points per minimal wavelength, surpassing the pseudospectral and high-order Lax-Wendroff correction methods. Numerical experiments for the homogeneous medium, layered medium, and 2D SEG/EAGE salt model indicated that the new method performed well in terms of computational efficiency and simulation accuracy, making it a highly useful tool for seismic modeling.