The hybrid simulation method is developed for simulating wave propagation only in a localized heterogeneous media with hybrid inputs obtained once for all from a known reference model. Despite the fact that the hybrid simulation method has a wide range of applications in computational seismology, the associated error control of this method has received relatively little attention in previous research works. We quantitatively discuss the error of the two-step hybrid method in acoustic wave cases and propose a spatial refinement scheme to compute hybrid inputs based on the multi-elements spline interpolation, which is preferable to traditional Lagrange interpolation since it uses more polydirectional interpolated points. This method can also be used for local refinement of wavefield in more general applications, such as saving smooth wavefield in the full-waveform inversion framework. Furthermore, to save memory requirements, hybrid inputs are proposed to be sparsely stored with a high upsampling ratio during the global simulation, and the Fourier interpolation method is introduced to recover them to their original time series. To demonstrate the effect of the proposed methods, we perform several 2D and 3D hybrid wave numerical simulations using the spectral element method. We find that when the global and local meshing differs, the proposed spatial interpolation method can appreciably reduce the error of the hybrid waveforms caused by inaccurate hybrid inputs. We also point out that the Fourier interpolation can efficiently recover the original waveform, allowing hybrid inputs to be stored with time steps toward the Nyquist limit. Our method is expected to become a standard method to reduce the error of hybrid waveforms and save the memory requirements during hybrid simulations and has potential implications for further improving the accuracy of the so-called box tomography.