A method for computing Laplace and Stokes interactions among N spherical particles arbitrarily placed in a unit cell of a periodic array is described. The method is based on an algorithm by Greengard and Rokhlin [J. Comput. Phys. 73, 325 (1987)] for rapidly summing the Laplace interactions among particles by organizing the particles into a number of different groups of varying sizes. The far-field induced by each group of particles is expressed by a multipole expansion technique into an equivalent field with its singularities at the center of the group. The resulting computational effort increases only linearly with N. The method is applied to a number of problems in suspension mechanics with the goal of assessing the efficiency and the potential usefulness of the method in studying dynamics of large systems. It is shown that reasonably accurate results for the interaction forces are obtained in most cases even with relatively low-order multipole expansions.