A nonlinear Fredholm integral equation of the first kind for the perturbation potential can be derived by interpreting the acoustic velocity as a perturbation of a reference velocity. We present an accurate and efficient method to formulate and numerically solve this equation with no restriction on the size of the perturbation by introducing a velocity‐weighted wavefield function (i.e., the scalar product of the potential function with the perturbed velocity function inside a scatterer), which is an intermediate function associating the observed wavefield with the perturbed velocity function. We start with a singularity analysis of the Fredholm integral equation when the observation point coincides with the scattering point to establish a linear singular integral equation with respect to the velocity‐weighted wavefield function. Then we formulate the relation between the velocity‐weighted wavefield and scattered field in the f‐k domain. A numerical scheme is developed to solve the forward problem for the velocity‐weighted wavefield inside a scatterer. Then the observed wavefield can be calculated through back substitutions of the Fredholm integral equation. The total wavefield and scattered field from a triangular object in a host medium observed from four sides are computed and compared with those given by the boundary‐element (BE) method and the approach using Born approximation. A second model computes both the common shot‐point (CSP) and common depth‐point (CMP) records of three velocity layers in a host medium. In addition, theoretical analyses and numerical experiments show that we can also use the singular integral equation to decompose accurately the velocity‐weighted wavefield given by inverse algorithms to recover the perturbed velocity function.