Accurately characterizing the forces acting on particles in fluids is of fundamental importance for understanding particle dynamics and binding kinetics. Conventional asymptotic solutions may lead to poor accuracy for neighboring particles. In this paper, we develop an accurate boundary integral method to calculate forces exerted on particles for a given velocity field. We focus our study on the fundamental two-bead oscillating problem in an axisymmetric frame. The idea is to exploit a correspondence principle between the unsteady Stokes and linear viscoelasticity in the Fourier domain such that a unifying boundary integral formulation can be established for the resulting Brinkman equation. In addition to the dimension reduction vested in a boundary integral method, our formulation only requires the evaluation of single-layer integrals, which can be carried out efficiently and accurately by a hybrid numerical integration scheme based on kernel decompositions. Comparison with known analytic solutions and existing asymptotic solutions confirms the uniform third-order accuracy in space of our numerical scheme.