We have developed a finite-difference iterative solver of the Helmholtz equation for seismic modeling and inversion in the frequency domain. The iterative solver involves the shifted Laplacian operator and two-level preconditioners. It is based on the application of the preconditioners to the Krylov subspace stabilized biconjugate gradient method. A critical factor for the iterative solver is the introduction of a new preconditioner into the Krylov subspace iteration method to solve the linear equation system resulting from the discretization of the Helmholtz equation. This new preconditioner is based on a reformulation of an integral equation-based convergent Born series for the Lippmann-Schwinger equation to an equivalent differential equation. We have determined that our iterative solver combined with the novel preconditioner when incorporated with the finite-difference method accelerates the convergence of the Krylov subspace iteration method for frequency-domain seismic wave modeling. A comparison of a direct solver, a one-level Krylov subspace iterative solver, and our two-level iterative solver verified the accuracy and accelerated convergence of the new scheme. Extensive tests in full-waveform inversion demonstrate the solver’s applicability to such problems.