Nonlinear dynamic analysis of the assembled structures involves the complex nonlinearity of the joint interfaces. By combining the multi-harmonic balance method with the Newton–Raphson iteration process, a high-efficiency nonlinear dynamic analysis method is proposed to analyze the steady-state dynamic responses influenced by the hysteresis nonlinearity of the joint interfaces. The proposed method adopts a gradient vector based on the chain derivation of the continuous function to iterate the steady-state dynamic responses instead of the Jacobian matrix in each iteration process, and a 2-order interpolation function is chosen to compensate the continuity of the discrete harmonic coefficients. The proposed method is verified by a comparison with the Jacobian matrix-based Newton iteration method and the direct numerical integration by using a single degree-of-freedom (DOF) system, and the effects of the excitation amplitudes are also investigated. Based on the generalized mode superposition method, the complex lap-type bolted joint beam structure is reduced into a series of single DOF oscillators with coupling nonlinear dynamic behaviors of the joint interfaces. The proposed method is applied to calculate the frequency response function (FRF) and verified by a comparison with the prediction of the Jacobian matrix-based iteration method. The results show that the solutions of the proposed method agree well with the results predicted by the direct numerical integration and Jacobian matrix-based iteration in literature. Larger excitation amplitude will induce the loss of the joint stiffness and degrade the nonlinearity effects. The predicted FRFs of the bolted joint beam system agree well with the results iterated by the Jacobian matrix-based method, and indicate a better convergence performance and a higher efficiency with much less computational cost.