AbstractThe peridynamic theory brings advantages in dealing with discontinuities, dynamic loading, and nonlocality. The integro‐differential formulation of peridynamics poses challenges to numerical solutions of complicated and practical problems. Some important issues attract much attention, such as the computation of infinite domains, the treatment of softening of boundaries due to an incomplete horizon, and time accumulation error in dynamic processes. In this work, we develop the boundary element method of peridynamics (PD‐BEM). The numerical examples demonstrate that the PD‐BEM exhibits several features. First, for nondestructive cases, the PD‐BEM can be one to two orders of magnitude faster than the meshless particle method of peridynamics (PD‐MPM) that directly discretizes the computational domains; second, it eliminates the time accumulation error, and thus conserves the total energy much better than the PD‐MPM; third, it does not exhibit spurious boundary softening phenomena. For destructive cases where new boundaries emerge during the loading process, we propose a coupling scheme where the PD‐MPM is applied to the cracked region and the PD‐BEM is applied to the uncracked region such that the time of computation can be significantly reduced.