## 3 Numerical Calculation Method in Different States  
### Numerical Calculation Method in the Slip State  
From hereon, the superscript \(*\) representing dimensionless quantities is omitted for simplicity. The unknown friction force \(f\) can be eliminated by substituting Equation (11) into Equations (14)-(16). The equations for the unknown variables \(\theta\), \(\varphi\), and \(N\) become  
\[A_{1}\ddot{\theta}+B_{1}\dot{\varphi}=C_{1}+NH, \tag{18}\]  
\[A_{2}\ddot{\theta}+B_{2}\ddot{\varphi}=C_{2}+Nl_{1}G_{1}, \tag{19}\]  
\[A_{3}\ddot{\theta}+B_{3}\ddot{\varphi}=C_{3}+Nl_{2}G_{2}, \tag{20}\]  
where \(A_{1}\), \(A_{2}\), \(A_{3}\), \(B_{1}\), \(B_{2}\), \(B_{3}\), \(C_{1}\), \(C_{2}\), \(C_{3}\), \(H\), \(G_{1}\), and \(G_{2}\) are  
\[A_{1}=\frac{-m_{d}l_{1}\sin\theta}{\cos\theta_{0}}+m_{a}l_{g1}\sin(\theta+ \theta_{0})+m_{b}l_{1}\tan\theta_{0}\cos(\theta+\theta_{0}), \tag{21}\]  
\[A_{2}=\frac{-m_{a}l_{1}l_{g1}\sin\theta}{\cos\theta_{0}}\sin(\theta+\theta_{0 })+\frac{m_{b}l_{1}^{2}\cos\theta}{\cos\theta_{0}}\cos(\theta+\theta_{0})+m_ {a}l_{g1}^{2}+I_{1}, \tag{22}\]  
\[A_{3}=\frac{m_{b}l_{1}l_{g2}\cos\varphi}{\cos\theta_{0}}\cos(\theta+\theta_{0 }), \tag{23}\]  
\[B_{1}=\frac{-l_{2}\sin\varphi}{\cos\theta_{0}}+m_{b}l_{g2}\sin(\varphi+\theta _{0}), \tag{24}\]  
\[B_{2}=m_{b}l_{1}l_{g1}\cos(\varphi-\theta)-\frac{m_{e}l_{2}\sin\varphi\sin( \theta+\theta_{0})}{\cos\theta_{0}}, \tag{25}\]  
\[B_{3}=\frac{-m_{b}l_{2}l_{g2}\sin\varphi}{\cos\theta_{0}}\sin(\varphi+\theta_{ 0})+m_{b}l_{g2}^{2}+I_{2}, \tag{26}\]  
\[C_{1}= \dot{\theta}^{2}\bigg{[}\frac{m_{d}l_{1}\cos\theta}{\cos\theta_{0}}-m_ {a}l_{g1}\cos(\theta+\theta_{0})+m_{b}l_{1}\tan\theta_{0}\sin(\theta+\theta_{ 0})\bigg{]}\]  
\[A_{2}\ddot{\theta}+B_{2}\ddot{\varphi}=C_{2}+Nl_{1}G_{1}, \tag{19}\]0  
\[A_{2}\ddot{\theta}+B_{2}\ddot{\varphi}=C_{2}+Nl_{1}G_{1}, \tag{19}\]1 \[A_{2}\ddot{\theta}+B_{2}\ddot{\varphi}=C_{2}+Nl_{1}G_{1}, \tag{19}\]2  
\[A_{2}\ddot{\theta}+B_{2}\ddot{\varphi}=C_{2}+Nl_{1}G_{1}, \tag{19}\]3 \[A_{2}\ddot{\theta}+B_{2}\ddot{\varphi}=C_{2}+Nl_{1}G_{1}, \tag{19}\]4  
\[A_{2}\ddot{\theta}+B_{2}\ddot{\varphi}=C_{2}+Nl_{1}G_{1}, \tag{19}\]5  
\[A_{2}\ddot{\theta}+B_{2}\ddot{\varphi}=C_{2}+Nl_{1}G_{1}, \tag{19}\]6  
\[A_{2}\ddot{\theta}+B_{2}\ddot{\varphi}=C_{2}+Nl_{1}G_{1}, \tag{19}\]7  
Using Equation (20), the normal force \(N\) is expressed in terms of \(\theta\) and \(\varphi\) as  
\[A_{2}\ddot{\theta}+B_{2}\ddot{\varphi}=C_{2}+Nl_{1}G_{1}, \tag{19}\]8Substituting Equation (33) into Equations (18) and (19) yields  
\[A_{2}\ddot{\theta}+B_{2}\ddot{\varphi}=C_{2}+Nl_{1}G_{1}, \tag{19}\]9  
\[A_{3}\ddot{\theta}+B_{3}\ddot{\varphi}=C_{3}+Nl_{2}G_{2}, \tag{20}\]0  
Equations (34) and (35) can also be written in matrix form as follows:  
\[ \begin{bmatrix}a_{11}&a_{12}\\ a_{21}&a_{22}\end{bmatrix}\begin{bmatrix}\ddot{\theta}\\ \dot{\varphi}\end{bmatrix} = \begin{bmatrix}P_{1}\\ P_{2}\end{bmatrix} , \tag{36}\]  
where \(a_{11}\), \(a_{12}\), \(a_{21}\), \(a_{2}\), \(P_{1}\), and \(P_{2}\) are  
\[a_{11}=A_{1}-\frac{A_{3}H}{l_{2}G_{2}}, \tag{37}\]  
\[a_{12}=B_{1}-\frac{B_{3}H}{l_{2}G_{2}}, \tag{38}\]  
\[a_{21}=A_{2}-\frac{l_{1}A_{3}G_{1}}{l_{2}G_{2}}, \tag{39}\]  
\[a_{22}=B_{2}-\frac{l_{1}B_{3}G_{1}}{l_{2}G_{2}}, \tag{40}\]  
\[P_{1}=C_{1}-\frac{C_{3}H}{l_{2}G_{2}}, \tag{41}\]  
\[P_{2}=C_{2}-\frac{l_{1}C_{3}G_{1}}{l_{2}G_{2}}. \tag{42}\]  
The differential equations for \(\theta\) and \(\varphi\) can be rewritten as follows using Equation (36):  
\[\ddot{\theta}=\frac{a_{22}P_{1}-a_{12}P_{2}}{a_{11}a_{22}-a_{12}a_{21}}, \tag{43}\]  
\[\ddot{\varphi}=\frac{-a_{21}P_{1}+a_{11}P_{2}}{a_{11}a_{22}-a_{12}a_{21}}. \tag{44}\]  
The state equation can then be written as Equation (45), and the Runge-Kutta method can be used for the numerical calculation of \(\theta\) and \(\varphi\). \[\frac{\mathrm{d}}{\mathrm{d}t} \begin{bmatrix}\theta\\ \theta\\ \varphi\\ \dot{\varphi}\end{bmatrix} = \begin{bmatrix}\dot{\theta}\\ \dot{\theta}\\ \dot{\varphi}\\ \dot{\varphi}\end{bmatrix} . \tag{45}\]