Context. Simulations of asteroid binaries commonly use mutual gravitational potentials approximated by series expansions, leading to truncation errors, and also preventing correct computations of force and torque for certain configurations where the bodies have overlapping bounding spheres, such as in the rotational fission model for creating asteroid binaries and pairs. Aims. We address errors encountered when potentials truncated at order two and four are used in simulations of binaries, as well as other errors related to configurations with overlapping bounding spheres where the series diverge. Methods. For this we utilized a recently developed method where the gravitational interaction between two triaxial ellipsoids can be calculated without approximations for any configuration. The method utilizes surface integration for both force and torque calculations, and it is exact for ellipsoidal shapes. We also computed approximate solutions using potentials truncated at second and fourth order, and we compare these with the solutions obtained with the surface integral method. The approximate solutions were generated with the “General Use Binary Asteroid Simulator” (GUBAS). Results. If the secondary is located with its centroid in the equatorial plane of the primary, the error in the force increases as the secondary is moved closer to the primary, but is still relatively small for both second and fourth order potentials. For torque calculations, the errors become more significant, especially if the other body is located close to one of the extended principal axes. On the axes themselves, the second order series approximation fails by 100%. For dynamical simulations of components separated a few primary radii apart, the fourth order approximation is significantly more accurate than the second order. Furthermore, because of larger errors in the torque calculations, the rotational motion is subject to greater inaccuracies than the translational motion. For configurations resembling contact binaries where the bounding spheres overlap, the errors in both force and torque in the initial stages of the simulation are considerable, regardless of the approximation order, because the series diverge. A comparison of the computational efficiency of the force and torque calculations shows that the surface integration method is approximately 82 times and four times slower than the second and fourth order potentials, respectively, but approximately 16 times faster than the order eight potential. Comparing the computation efficiency of full simulations, including the calculations of the equations of motion, shows that the surface integration scheme is comparable with GUBAS when an order four potential is used. Conclusions. The errors generated when mutual gravitational potentials are truncated at second or fourth order lead to larger errors in the rotational than in the translational motion. Using a mathematically exact method for computing forces and torques becomes important when the bodies are initially close and the bounding spheres overlap, in which case both the translational and rotational motion of the bodies have large errors associated with them. For simulations with two triaxial ellipsoids, the computational efficiency of the surface integral method is comparable to fourth order approximations with GUBAS, and superior to eight order or higher.