## 2 Materials and Methods
### Integration of Dynamic Equations
Thus, to calculate the inverse of a diagonal matrix of non-zero elements, it is only necessary to calculate the inverse of the particular element \(m_{i}\), that is, \(1/m_{i}\). Once these positions are calculated, the contact forces are calculated for which the mass matrix is no longer necessary following [PERSON]'s law where the nodes, which are the ends of the contact spring, and the stiffness of the spring are known. In Equation (18), these terms appear together although they can be calculated, as mentioned above, in a decoupled manner. According to reference [23], this approach introduces a small error that can be ignored in structural systems with low damping, as seen in our case. Next, the velocity and acceleration in the coordinates are approximated by the following equations:  
\[ \begin{array}{c}\dot{q}_{n-\frac{1}{2}}=\frac{1}{2\Delta t}(q_{n}-q_{n-1}) \\ \ddot{q}_{n}=\frac{1}{2\Delta t}\left(\dot{q}_{n+\frac{1}{2}}-\dot{q}_{n-\frac {1}{2}}\right)=\frac{1}{\Delta t^{2}}(q_{n+1}-2q_{n}+q_{n-1})\end{array}  \tag{19}\]  
Substituting the previous equations into Equation (18) and solving for the vector \(q_{n+1}\) yields:  
\[\left(\frac{1}{\Delta t^{2}}\right)Mq_{n+1}=R_{n}-K_{n}q_{n}-\varnothing_{n}^{ t}\lambda_{n}+\frac{1}{\Delta t^{2}}M\left(q_{n}+\Delta t\dot{q}_{n-\frac{1}{2}} \right)-C_{n}\dot{q}_{n-\frac{1}{2}} \tag{20}\]  
In this case, it is possible to calculate the vector \(q_{n+1}\) directly without the need to solve a system of equations because the mass matrix is diagonal. However, some variables remain to be determined: the position of the contact elements \((y_{3})_{n+1}\) and \((y_{4})_{n+1}\), and contact forces \((\lambda_{1})_{n+1}\) and \((\lambda_{2})_{n+1}\). The mass of the contact elements is null and Equation (20) cannot be applied; however, the positions of the catenary nodes in the pantograph environment are known at \(t_{n+1}\). Thus, when the contact elements are activated, we use Equation (10) particularized in \(t_{n+1}\) for \((y_{3})_{n+1}\) and a similar equation for \((y_{4})_{n+1}\). When the contact element is deactivated, its position coincides with that of the head mass; in this case, we will use Equation (14) or (15). For the contact forces in \(t_{n+1}\), these forces are calculated according to the deformation in the corresponding contact spring, according to the expressions:  
\[\left(\lambda_{1}\right)_{n+1}=k_{cont}\big{[}\left(y_{3}\right)_{n+1}+v_{1}- \left(y_{2}\right)_{n+1}\big{]},\ \ \left(\lambda_{2}\right)_{n+1}=k_{cont}\big{[}\left(y_{4}\right)_{n+1}+v_{2}- \left(y_{2}\right)_{n+1}\big{]} \tag{21}\]  
As can be seen in the previous equations, an additional correction term is added, given by the variables \(v_{1}\) and \(v_{2}\), to the theoretical deformation of the spring given by the difference in the positions of the contact element and the head mass. This term takes into account the effect of the slope in the sectioning span as a consequence of the previous configuration of the beam. The values of functions \(v_{1}\) and \(v_{2}\) are discussed in the following section. Because the contact force is a compression force, each integration step in Equation (21) must be solved by checking that the forces are negative. In this case, the value of the force is correct and the corresponding stiffness of the contact spring \(k_{3}\) or \(k_{4}\) is the theoretical value \(k_{cont}\) of the stiffness matrix in Equation (7), and this value is maintained for the next integration step. However, if a value of the force is positive, it means that the spring works in traction, and therefore a take-off, has taken place. In this case, the contact force becomes zero and the corresponding rigidity is annulled for the next step. This is equivalent to eliminating the contact catenary-pantograph. Based on the above, the algorithm for detecting take-off at instant \(t_{n+1}\) is:  
\[