There is substantial promise in the use of piezoelectric energy harvesting systems to power electronics and devices through the conversion of ambient vibrations into electrical energy. Accurately predicting the energy harvesting capability in these systems is of paramount importance as the harvesting structure and the associated circuit must be ‘tuned’ in concert to the signature of the ambient vibrations. Here, we develop a numerical simulation methodology for piezoelectric harvesting systems that allows for high-fidelity finite element models of the harvester electromechanics as well as realistic non-linear circuits models. Our approach synthesizes conventional harmonic finite element analysis and a harmonic balance approach that accounts efficiently for the non-linear circuit behavior, bypassing the need for computationally expensive time marching schemes. The proposed methodology solves, via an iterative residual reduction scheme, for the coupled electrical output and structural dynamics of the piezoelectric harvesting structure. The proposed method is well suited for integration into model-predictive design methods. By numerical examples we illustrate the capability of the proposed approach to predict previously unreported coupling phenomena between the non-linear circuit response and piezoelectric structural behavior. The increased accuracy in the circuit model and flexibility of this method over current models outweighs the increased numerical costs.