## 2 Materials and Methods
### Integration of Dynamic Equations  
To proceed with the integration of the differential Equation (1), we applied the explicit method of the central differences. Explicit methods are more suitable than implicit methods for dealing with non-linearities, in the particular case of the pantograph-catenary interaction [23], because they correctly take into account the variation of the stiffness matrix. Non-linearities occur due to debonding and coupling. Due to the obvious stability problems of explicit methods, the integration step must be carefully adjusted, which is undertaken in this development. The critical integration step that conditions the stability can be estimated as the element length divided by the wave propagation speed, as shown in Equations (13.10)-(16) (p. 401) of [23]. Regarding the element length, in this work all elements have the same length. If this is not the case, the shortest element would have to be taken for the calculation of the step. With respect to the propagation speed, it depends on the mechanical characteristics of the beam section. However, the authors have deduced the critical integration step by trial and error so that, when the integration step is above the critical time, the algorithm degenerates. Otherwise, the algorithm gives reasonable results. Therefore, different increments of \(t\) have tried to determine the critical value below when the algorithm is stable, and the results are feasible. Thus, for this proposal, the speed and acceleration in the coordinates are initially approximated for a moment of time \(t_{n}\), and an interval \(\Delta t\), in accordance with:  
\[ \begin{array}{l}\dot{q}_{n}=\frac{1}{2\Delta t}(q_{n+1}-q_{n-1})\\ \ddot{q}_{n}=\frac{1}{\Delta t^{2}}(q_{n+1}-2q_{n}+q_{n-1})\end{array}  \tag{16}\]  
Figure 6: Pantograph with different activation cases in the contact elements. Substituting these equations into Equation (1) yields the following expression:  
\[\left(\frac{1}{\Delta t^{2}}M+\frac{1}{2\Delta t}C_{n}\right)q_{n+1}=R_{n}-K_{n}q _{n}-\varnothing_{n}^{t}\lambda_{n}+\frac{1}{\Delta t^{2}}M(2q_{n}-q_{n-1})+ \frac{1}{2\Delta t}C_{n}q_{n-1} \tag{17}\]  
Equation (17) allows the coordinates \(q_{n+1}\) at time \(t_{n+1}\) to be obtained from the elements of the equation and the coordinates at the previous instant \(t_{n}\). However, for this calculation, it is necessary to solve a system of linear equations because, although the mass matrix \(M\) is diagonal, the same does not occur with the damping matrix \(C_{n}\). After all, we assumed a Rayleigh damping model for the catenary according to Equation (4). Although it is possible to assume a diagonal matrix for \(C_{n}\), we preferred to use the Rayleigh model because we think it is a more realistic approach. The same problem exists for the pantograph-damping matrix. One of the main advantages of the explicit integration methods is that the integration variables can be obtained directly without having to solve a system of equations, resulting in a more computationally efficient integration. However, for this purpose, it is necessary to slightly modify the dynamic equation by correcting the velocity vector in the coordinates by half a step, which results in:  
\[\left( \begin{array}{cc}M&0\\ 0&0\end{array} \right)\left( \begin{array}{c}\ddot{q}_{n}\\ \dot{\lambda}_{n}\end{array} \right)+\left( \begin{array}{cc}C_{n}&0\\ 0&0\end{array} \right)\left( \begin{array}{c}\dot{q}_{n-\frac{1}{2}}\\ \dot{\lambda}_{n-\frac{1}{2}}\end{array} \right)+\left( \begin{array}{cc}K_{n }&\varnothing_{n}^{t}\\ \varnothing_{n}&0\end{array} \right)\left( \begin{array}{c}q_{n}\\ \lambda_{n}\end{array} \right)=\left( \begin{array}{c}R_{n}\\ 0\end{array} \right) \tag{18}\]  
Equation (18) involves the calculation of the inverse of a rank-deficient mass matrix. However, the algorithm actually works by firstly calculating the position of the nodes, for which the values of the mass matrix are non-zero following the method of central differences. This method appears in [23].