High-resolution image and waveform inversion of small-scale targets requires the handling of high-frequency seismic wavefields. However, conventional finite-difference (FD) methods have strong numerical dispersions in the presence of high-frequency components. To reduce these numerical dispersions, we optimized the constant coefficients of the FD operator by maximizing the wavenumber coverage within a given error limitation. We set up three general criteria to enhance the convergence of the algorithm and reduce the optimization effort. We selected the error limitation to be 0.0001, this being the smallest in the literature, which led to perfect agreement between theoretical analyses and numerical experiments. The accuracy of our optimized FD methods can even reach that of much higher order unoptimized FD methods, which means great savings of computational efforts and memory demand. These advantages become even more apparent with 3D modeling, especially for saving memory demand.