## 2 Methods
### Pseudo Arc-Length Continuation Method  
For tracing the solution curve of nonlinear equations like Equation (3), the Newton-Raphson method is usually applied. However, the Newton iterative method usually fails when buckling or bifurcation happens, so the pseudo arc-length continuation method [30] is used to solve Equation (3). The detailed procedure is described in the following paragraphs. In the pseudo arc-length continuation method, an extra parameter s is introduced to control the solution searching process. Equation (3) can be rewritten as the following expression:  
\[\psi(\mathbf{Q},\omega)=\mathbf{Z}(\omega)\mathbf{Q}+\mathbf{F}_{\mathbf{n}}( \mathbf{Q})+\mathbf{F}_{\mathbf{u}}=0 \tag{6}\]  
where \(\psi\) is represented to be related to \(\omega\) because \(\omega\) also varies during the iterative process. At the first iterative step, two initial solutions away from strongly nonlinear region is obtained using Newton method based on Equation (3), i.e.,:  
\[\mathbf{Q}^{(m+1)}=\mathbf{Q}^{(m)}-\left(\mathbf{Z}(\omega)+\frac{\partial \mathbf{F}_{n}}{\partial\mathbf{Q}^{(m)}}\right)\mathbf{R}\Big{(}\mathbf{Q}^{ (m)}\Big{)} \tag{7}\]  
where \(m\) is the number of previous iterative step. Then the first arc length \(s_{1}\) for outer iteration is calculated in the following form:  
\[s_{1}=\sqrt{\zeta(\omega_{1}-\omega_{2})^{2}+\left(\mathbf{Q}_{1}-\mathbf{Q} _{2}\right)^{2}} \tag{8}\]  
where \(\omega_{1}\), \(\mathbf{Q}_{1}\) and \(\omega_{2}\), \(\mathbf{Q}_{2}\) are the results obtained using Newton method; \(\zeta\) is the weight coefficient. Then the inner iteration could be executed but an additional vector is required firstly. It is given as:  
\[\mathbf{v}(\mathbf{Q},\zeta)=\left[v_{1},v_{2},\cdots,v_{N},v_{N+1}\right] \tag{9}\]  
where \(v_{i}=(-1)^{i}\text{det}\Big{[}\frac{\partial\mathbf{b}}{\partial\mathbf{Q} _{1}^{\prime}},\frac{\partial\mathbf{b}}{\partial\mathbf{Q}_{2}^{\prime}}, \cdots,\frac{\partial\mathbf{b}}{\partial\mathbf{Q}_{N}^{\prime}},\frac{ \partial\mathbf{b}}{\partial\mathbf{Q}_{N}^{\prime}}\Big{]}\) and the superscript \(\hat{\ }\) denotes elimination of the corresponding column. Obviously, vector \(\mathbf{v}(\mathbf{Q},\omega)\) is the tangential vector of the tracing curve with respect to \(\mathbf{Q}\) and \(\omega\). The unit tangential vector of solution curve is derived as \(\mathbf{\tau}(\mathbf{Q},\omega)=\frac{\mathbf{v}}{\|\mathbf{v}\|}\) where \(\|\mathbf{v}\|\) represents the Euclidean norm of vector \(\mathbf{v}\). For the \(m\)th outer iteration, the constraint equation of \(s_{m}\) can be expressed as follows:  
\[\mathbf{B}(\mathbf{Q},\omega,s)=\mathbf{\tau}^{T}(s)\bigg{\{} \begin{array}{ c}\Delta\mathbf{Q}(s)\\ \Delta\omega(s)\end{array} \bigg{\}}-\Delta s=0 \tag{10}\]  
where \(\Delta\mathbf{Q}(s)=\mathbf{Q}(s)-\mathbf{Q}_{m}\) and \(\Delta\omega(s)=\omega(s)-\omega_{m}\); \(\Delta s\) is related to the arc length for next step of outer iteration, which can be expressed as \(\Delta s=s_{m+1}-s_{m}\). Then the extra constraint equation Equation (10) can be reorganized into Equation (6) and new equations take the form:  
\[\psi^{*}(\mathbf{Q  
After all those derivations, the process of inner iteration takes the form:  
\[\left\{ \begin{array}{c}\Delta\left[\begin{array}{c}\mathbf{Q}\\ \omega\end{array}\right]^{1}=\Delta\left[\begin{array}{c}\mathbf{Q}\\ \omega\end{array}\right]^{E}\\ \Delta\left[\begin{array}{c}\mathbf{Q}\\ \omega\end{array}\right]^{n+1}=\Delta\left[\begin{array}{c}\mathbf{Q}\\ \omega\end{array}\right]^{n}+\mathbf{J}^{*-1}\left[\begin{array}{c}- \mathbf{\Psi}(\left[\begin{array}{c}\mathbf{Q}\\ \omega\end{array}\right]_{m}+\Delta\left[\begin{array}{c}\mathbf{Q}\\ \omega\end{array}\right]^{n})\\ 0\end{array}\right]\end{array} \right\} \tag{13}\]  
And the initial solution searching point for next outer iteration is obtained by:  
\[\left[ \begin{array}{c}\mathbf{Q}\\ \omega\end{array} \right]^{E}=\left[ \begin{array}{c}\mathbf{Q}\\ \omega\end{array} \right]_{m}+\tau(\left[ \begin{array}{c}\mathbf{Q}\\ \omega\end{array}