Nonlinear piezoelectric energy harvesters (NPEHs) have attracted much attention in recent years due to their wide harvesting bandwidth and high harvesting efficiency. Due to the complexity of NPEHs, researchers usually use theoretical numerical and experiments to analyze. However, two problems prevent us from more accurate analysis and design. Firstly, the common magnetic model named magnetic dipole model (MDM) will usually cause an unacceptable error in NPEHs; secondly, the harmonic balance method (HBM) which usually only uses the first-order harmonic terms is accurate in the monostable configuration for low-amplitude vibrations, but for the bistable or multi-stable ones, this method has major differences compared to numerical methods such as sweep method (SM) or experiment. To solve the above two problems, we propose a more accurate magnetic model based on the MDM to reduce the error. Meanwhile, we use the swarm intelligent optimization algorithm (SIOA) to solve the high-dimensional nonlinear algebraic equations obtained by high order HBM, which is in agreement with the SM, and find out more nonlinear dynamical phenomena. What’s more, the high-energy orbital solutions which can be easily neglected in the SM and experiment but have significant advantages to energy harvesting have been found.