In comparison with numerical methods, theoretical characterizations of ion motion in the nonlinear Paul traps always suffer from low accuracy and little applicability. To overcome the difficulties, the theoretical harmonic balance (HB) method was developed, and was validated by the numerical fourth-order Runge-Kutta (4th RK) method. Using the HB method, analytical ion trajectory and ion motion frequency in the superimposed octopole field, ε, were obtained by solving the nonlinear Mathieu equation (NME). The obtained accuracy of the HB method was comparable with that of the 4th RK method at the Mathieu parameter, q = 0.6, and the applicable q values could be extended to the entire first stability region with satisfactory accuracy. Two sorts of nonlinear effects of ion motion were studied, including ion frequency shift, Δβ, and ion amplitude variation, Δ(C(2n)/C0) (n ≠ 0). New phenomena regarding Δβ were observed, although extensive studies have been performed based on the pseudo-potential well (PW) model. For instance, the |Δβ| at ε = 0.1 and ε = -0.1 were found to be different, but they were the same in the PW model. This is the first time the nonlinear effects regarding Δ(C(2n)/C0) (n ≠ 0) are studied, and the associated study has been a challenge for both theoretical and numerical methods. The nonlinear effects of Δ(C(2n)/C0) (n ≠ 0) and Δβ were found to share some similarities at q < 0.6: both of them were proportional to ε, and the square of the initial ion displacement, z(0)(2).