## 3 Computational Methods and Validation
### Harmonic Balance Method with Alternating Frequency/Time-Domain Technique
Near turning points, conventional numerical iteration schemes may exhibit poor convergence characteristics or computational divergence due to substantial variations in the solution curve's derivative structure. To overcome this turning point issue, an arc-length continuation algorithm incorporating predictor-corrector methodology [30] is adopted, which can more efficiently and accurately capture turning points by combining prediction and correction steps. By introducing a new parameter \(\lambda\), the solution curve of the equation is set as follows:  
\[\mathbf{Z}=\left\{\mathbf{X};\lambda\right\}_{(N_{H}+1)\times 1} \tag{19}\]  
where \(N_{H}\) is the total number of harmonic balance equations, and \(\lambda\) is the continuation parameter. In the analysis of steady-state response characteristics, the rotational speed \(\Omega\) is selected as the continuation parameter. Figure 3 illustrates a schematic representation of the predictor-corrector methodology, with the computational procedure detailed below. To avoid derivative singularities at turning points that may compromise numerical convergence, a tangent predictor method introduced by [PERSON] [28] is implemented. The prediction step length \(\Delta\mathbf{Z}_{i}\) is as follows:  
\[f\Big{(}\mathbf{Z}_{i}^{k}\Big{)}\Delta\mathbf{Z}_{i}=\mathbf{e}_{N_{H}+1}\Delta\text{s} \tag{20}\]  
where \(\mathbf{e}_{N_{H}+1}\) represents a \(N_{H}+1\) dimensional unit vector with the last component equal to \(1\). \(\Delta\text{s}\) is the defined length step size along the solution path. The Jacobian matrix for the prediction step \(f\Big{(}\mathbf{Z}_{i}^{k}\Big{)}\) is given by the following:  
\[f\Big{(}\mathbf{Z}_{i}^{k}\Big{)}= \begin{bmatrix}\frac{\partial R(\mathbf{X}_{i}^{k})} {\partial\mathbf{X}}&\frac{\partial R(\mathbf{X}_{i}^{k})}{\partial\lambda}\\ (\mathbf{e}_{N_{H}})^{\text{T}}&0\end{bmatrix}  \tag{21}\]  
where \(\mathbf{e}_{N_{H}}\) is the transpose of an \(N_{H}\)-dimensional unit vector, characterized by a single unity component strategically selected to ensure the full rank of \(f\Big{(}\mathbf{Z}_{i}^{k}\Big{)}\). This rank completeness guarantees the existence of a tangent vector at point \((\mathbf{X}_{i},\lambda)\), enabling the computation of the predictor step. The state update equation for the i-th prediction step is expressed as follows:  
\[\mathbf{Z}_{i+1}=\mathbf{Z}_{i}+\Delta\mathbf{Z}_{i} \tag{22}\]  
The numerical continuation employs an arc-length correction methodology. Given the computational complexity of exact arc-length determination, an approximation approach is implemented. The arc-length constraint equation, perpendicular to the tangent vector through point \(\mathbf{Z}_{i-1}^{k}\), takes the form:  
\[p_{arc}\Big{(}\mathbf{Z}_{i}^{k}\Big{)}=\Big{(}\mathbf{Z}_{i}^{k}-\mathbf{Z}_{i-1}^{k} \Big{)}^{\text{T}}\Big{(}\mathbf{Z}_{i}^{k}-\mathbf{Z}_{i-1}^{k}\Big{)}-\Delta\text{s} ^{2}=0 \tag{23}\]  
where \(i\) denotes the prediction step number, \(k\) represents the correction step number. The arc-length constraint equation combines with the original residual \(\mathbf{R}(\mathbf{X})\) to form a new residual function \(\mathbf{R}_{ext}\Big{(}\mathbf{Z}_{i}^{k}\Big{)}\):  
\[\mathbf{R}_{ext}\Big{(}\mathbf{Z}_{i}^{k}\Big{)}=\Big{[}\mathbf{R}(\mathbf{X}_{i}^{k});\;\;\;p _{arc}(\mathbf{Z}_{i}^{k})\Big{]} \tag{24}\]  
Figure 3: Diagram of the predictor–corrector method. Using the Newton-Raphson iteration method to solve Equation (24), the correction step length \(\Delta\mathbf{Z}^{k}\) is expressed as follows:  
\[\Delta\mathbf{Z}^{k}=-J_{\text{ext}}\Big{(}\mathbf{Z}_{i}^{k}\Big{)}^{-1}R_{ext} \Big{(}\mathbf{Z}_{i}^{k}\Big{)},\quad J_{ext}\Big{(}\mathbf{Z}_{i}^{k}\Big{)}= \begin{bmatrix}\frac{\partial\mathbf{R}(\mathbf{X}_{i}^{k})}{\partial\mathbf{X }}&\frac{\partial\mathbf{R}(\mathbf{X}_{i}^{k})}{\partial\lambda}\\ \frac{d_{ppx}}{\partial\mathbf{X}}&\frac{d_{ppx}}{\partial\lambda}\end{bmatrix}