This article presents a wavenumber error optimized nearly analytic discrete differentiator (WONAD) using Taylor polynomial constraints and wavenumber domain optimizing in the function space expanded by the displacement and its gradient. We also apply the differentiator in a 2D scalar‐wave equation for seismic‐wave forward modeling in heterogeneous media. Subsequent analysis shows the WONAD generates both low numerical dispersion and numerical error. The maximum phase velocity error of the WONAD is as low as 3.12%, even in an extreme case with only two sampling points per minimum wavelength when the Courant number is 0.5. In our numerical experiments, the maximum relative error of the WONAD in a simple harmonic wavefield is less than 1.12% after a 300 s simulation. On the same grids, both the numerical dispersion and the numerical error of the WONAD are lower than what have been found in cases using traditional high‐order methods, such as the 24th‐order Lax–Wendroff correction (LWC) method and so on. To simulate a wavefield without visible numerical dispersion, the computational efficiency of the WONAD is 341%, 651%, and 316%, compared with that of the fourth‐order staggered grid method and the fourth‐order and 24th‐order LWC methods, while the phase misfit of the WONAD is still lower than the other three methods. The WONAD, showing a promising prospective in large‐scale long‐time seismic modeling, performs well in computational efficiency, simulation accuracy, and numerical stability.