## 3 Methodology
The central difference formula, which yields better results, is commonly used. In contrast to the forward and backward schemes, where the error decreases linearly, the error in the central difference method decreases quadratically as the step size decreases [26; 27]. In the FD approximation, the values of the dependent variable at a point in the solution region are related to the values at neighboring points. The solution region is divided into \(n\) subintervals, each of length \(\Delta x\). For an accurate approximation, the step size should be sufficiently small. In practice, large step sizes can accelerate simulations, but they may also lead to instability. To determine the appropriate step size for the central FD method and demonstrate its applicability, the forced vibration of a damped one-degree-of-freedom system is analyzed using the piecewise linear interpolation method. A semi-sinusoidal load, which is a piecewise function and shown in Figure 3, is applied to the system as an external load. Figure 3: A semi-sinusoidal, piecewise loading function. The load applied to the system is given by the following form:  
\[P(t)=\left\{ \begin{aligned} 0,& t<0\\ 10 sin\bigg{(}\frac{\pi t}{0.6}\bigg{)},& 0\leq t\leq 0.6\\ 0,& t>0.6\end{aligned} \right. \tag{17}\]  
The equation of motion for the forced vibration of the damped, one-degree-of-freedom system is as follows:  
\[m+ci\dot{u}+ku=P(t) \tag{18}\]  
The recurrence formulas for evaluating the dynamic state (\(u_{n+1}\), \(i_{n+1}\)) at time \(t_{n+1}\), given the state (\(u_{n}\), \(i_{n}\)) at time \(t_{n}\), are expressed as follows:  
\[u_{n+1}=Au_{n}+Bi\dot{u}_{n}+C(P_{n-1}\Delta t)+D(P_{n}\Delta t) \tag{19}\]  
\[\dot{u}_{n+1}=A^{\prime}u_{n}+B^{\prime}\dot{u}_{n}+C^{\prime}(P_{n-1}\Delta t )+D^{\prime}(P_{n}\Delta t) \tag{20}\]  
The coefficients \(A\) through \(D^{\prime}\) for the underdamped case (\(\zeta\) = 0.05) are as follows:  
\[A=e^{-\zeta\omega_{n}\Delta t}\Bigg{(}\frac{\zeta}{\sqrt{1-\zeta^{2}}}sin \omega_{D}\Delta t+cos\omega_{D}\Delta t\Bigg{)} \tag{21}\]  
\[B=e^{-\zeta\omega_{n}\Delta t}\Bigg{(}\frac{1}{\omega_{D}}  
The displacement vs. time graphs, obtained using the piecewise linear interpolation and central FD methods for time steps \(\Delta t\) = 0.1 and 0.01 s, are compared in Figures 4 and 5, respectively. As clearly shown in Figures 4 and 5, as the time step decreases, the results obtained from both methods converge. For \(\Delta t\) = 0.01 s, the convergence of the results is satisfactory for the central FD method. The subsequent analysis is performed with \(\Delta t\) = 0.002 s, starting with vanishing initial conditions, to obtain an accurate solution. By reducing \(\Delta t\) from 0.01 to 0.002 s, the number of time steps increases by a factor of 5, which also raises the computational time required for each iteration. However, earthquake data are provided at a finer time resolution (\(\Delta t\) = 0.002 s). Therefore, using a coarser time step (\(\Delta t\) = 0.01 s) could lead to aliasing or inaccuracies in how the earthquake data are applied. Using \(\Delta t\) = 0.002 s ensures that the model accurately captures the variations in the earthquake input.