Three-dimensional (3D) flutter theory is crucial in bridge flutter research. The existing multi-modal time-domain method suffers from low computational efficiency, and there's no clear approach to determine the initial vibration parameters for 3D multi-modal systems. Based on the theoretical formulas of flutter characteristics and the implicit-iteration algorithm, a 3D multi-modal time-domain analytical method (3DMTM-IA) is established. Meanwhile, a calculation method for determining the initial vibration parameters of the multi-modal system is also determined. Based on 3DMTM-IA, a reduced-order two-dimensional (2D) two-degree-of-freedom time-domain method (2DTM-IA) is obtained. A 2D nonlinear system and a long-span bridge are used to verify the applicability and calculation efficiency of 2DTM-IA and 3DMTM-IA, respectively. The results show that both implicit methods effectively study soft flutter responses. For the two cases here, the torsional mode mainly drives the soft flutter response, while vertical bending modes provide additional aerodynamic damping, affecting the torsional amplitude. Compared to previous time-domain methods, the implicit methods offer better computational efficiency in nonlinear flutter calculations. In addition, accurate initial vibration parameters are crucial for improving computational efficiency and avoiding divergence when using multi-modal time-domain analytical methods.