## 2 System Model
### Dynamic Modeling of the Coupled System
\begin{array}{l}\frac{\partial T}{\partial\theta}=\Big{\{}2J+\Big{[}\rho A \int_{0}^{l}x^{2}dx\Big{]}+\Big{[}\rho A\int_{0}^{l}\varphi^{2}(x)dx\Big{]}p^{ 2}\Big{\}}\dot{\theta}+\Big{[}\rho A\int_{0}^{l}x\varphi(x)dx\Big{]}\dot{p}\\ \frac{d}{dt}\left(\frac{\partial T}{\partial\theta}\right)=\Big{\{}2J+\Big{[} \rho A\int_{0}^{l}x^{2}dx\Big{]}+\Big{[}\rho A\int_{0}^{l}\varphi^{2}(x)dx \Big{]}p^{2}\Big{\}}\ddot{\theta}+2\Big{[}\rho A\int_{0}^{l}\varphi^{2}(x)dx \Big{]}\dot{\theta}\dot{p}p\\ \qquad\qquad+\Big{[}\rho A\int_{0}^{l}x\varphi(x)dx\Big{]}\ddot{p}\\ \frac{\partial T}{\partial\theta}=0,\\ \frac{\partial V}{\partial\theta}=2k\theta,\\ \frac{\partial D}{\partial\theta}=2c\dot{\theta}\end{array}  \tag{22}\]  
Then, we obtain  
\[ \begin{array}{l}\Big{\{}2J+\Big{[}\rho A\int_{0}^{l}x^{2}dx\Big{]}+\Big{[} \rho A\int_{0}^{l}\varphi^{2}(x)dx\Big{]}p^{2}\Big{\}}\ddot{\theta}+\Big{[} \rho A\int_{0}^{l}x\varphi(x)dx\Big{]}\ddot{p}\\ +2k\theta+2c\dot{\theta}+2\Big{[}\rho A\int_{0}^{l}\varphi^{2}(x)dx\Big{]} \dot{\theta}\dot{p}p=Q_{\theta}\end{array}  \tag{23}\]  
For generalized coordinate \(p\), based on the Lagrange equation Equations (3) and (4), we have:  
\[\frac{d}{dt}\left(\frac{\partial T}{\partial\dot{p}}\right)-\frac{\partial T}{ \partial p}+\frac{\partial V}{\partial p}+\frac{\partial D}{\partial\dot{p}}= Q_{p} \tag{24}\]  
where  
\[ \begin{array}{l}\frac{\partial T}{\partial\dot{p}}=\Big{[}\rho A\int_{0}^{l} x\varphi(x)dx\Big{]}\dot{\theta}+\Big{[}\rho A\int_{0}^{l}\varphi^{2}(x)dx \Big{]}\dot{p}\\ \frac{d}{dt}\left(\frac{\partial T}{\partial p}\right)=\Big{[}\rho A\int_{0}^{l }x\varphi(x)dx\Big{]}\ddot{\theta}+\Big{[}\rho A\int_{0}^{l}\varphi^{2}(x)dx \Big{]}\dddot{p}\\ \frac{\partial T}{\partial p}=\Big{[}\rho A\int_{0}^{l}\varphi^{2}(x)dx\Big{]} \dot{\theta}^{2}p\\ \frac{\partial V}{\partial p}=\Big{[}\rho A\int_{0}^{l}\varphi^{2}(x)dx\Big{]} \omega^{2}p\\ \frac{\partial D}{\partial\dot{p}}=0\end{array}  \tag{25}\]  
Then, we obtain \[ \begin{split}&\left[\rho A\int_{0}^{l}x\varphi(x)dx\right]\ddot{ \theta}+\left[\rho A\int_{0}^{l}\varphi^{2}(x)dx\right]\ddot{p}+\left[\rho A \int_{0}^{l}\varphi^{2}(x)dx\right]\omega^{2}p\\ &-\left[\rho A\int_{0}^{l}\varphi^{2}(x)dx\right]\dot{\theta}^{2} p=Q_{p}\end{split}  \tag{26}\]  
Equations (23) and (26) can be rewritten in matrix form:  
\[ \begin{bmatrix}M_{11}&M_{12}\\ M_{12}&M_{22}\end{bmatrix}\begin{bmatrix}\ddot{\theta}\\ \ddot{p}\end{bmatrix} + \begin{bmatrix}2k&0\\ 0&(2B_{2}\omega^{2})\end{bmatrix}\begin{bmatrix}\theta\\ p\end{bmatrix} + \begin{bmatrix}2c&0\\ 0&0\end{bmatrix}\begin{bmatrix}\dot{\theta}\\ \dot{p}\end{bmatrix} +\begin{bmatrix}H_{1}\\ H_{2  
Figure 4 presents the frequency response of the system under varying current intensities. In Table 2, when the electric current increases from 1 to 6 A, the joint torsional stiffness increases from 0.38 to 2.31 N \(\cdot\) m/rad, and the first resonant frequency increases from 0.21 to 0.48 Hz. Moreover, the resonance peak decreases from 54.27\({}^{\circ}\) to 23.74\({}^{\circ}\). In Figure 5, the comparison of the system frequency-shift effect is presented under different conditions: single joint and dual joints. Overall, the system natural frequency with dual joints is less than that with one joint. The reason is that although the torsional stiffness increases for dual joints, the load also increases. Therefore, the frequency shift of the variable stiffness caused by the increased load is much greater than that caused by the increased torsional stiffness.