The fluid is assumed to be inviscid and incompressible and the flow irrotational. The three_dimensional numerical model is established to simulate the interaction of multiple bubbles and the fast Fourier transform on multipoles (FFTM) method is combined with the higher order boundary method (BEM) to study the physics of multiple bubble dynamics. FFTM method is employed to speedup the solution of the boundary integral equation while achieving the same order of accuracy, enabling to simulate the dynamics of multiple bubbles in a reasonable time. Elastic mesh technique (EMT), which is a new mesh regulation technique, is applied to maintain the regularity of the triangular element mesh used to discretize the dynamic boundary surface during the evolution of bubbles. All these measures make the present approach viable and robust, which is validated by computations of several bubble dynamics problems. Numerical analysi was carried out for the interaction of multiple bubbles and the bubble dynamics. Some physical behaviors of the multiple bubbles are presented in this work. The factors affecting the expansion, collapse, moving of multiple bubbles and the jet formation are also discussed.