  This article deals with the efficient and accurate computation of the
electrostatic forces between charged, spherical dielectric particles undergoing
mutual polarisation. We use the spectral Galerkin boundary integral equation
framework developed by Lindgren et al. (J. Comput. Phys. 371 (2018): 712-731)
and subsequently analysed in two earlier contributions of the authors to
propose a linear scaling in cost algorithm for the computation of the
approximate forces. We establish exponential convergence of the method and
derive error estimates for the approximate forces that do not explicitly depend
on the number of dielectric particles $N$. Consequently, the proposed method
requires only $\mathcal{O}(N)$ operations to compute the electrostatic forces
acting on $N$ dielectric particles up to any given and fixed relative error.
