Quantitative Performance Analysis of Spring-Mass-Damper Control Systems: A Comparative Implementation in Python and R

1. Introduction

The numerical analysis and computational implementation of spring-mass-damper (SMD) systems present rich opportunities for investigating fundamental questions in dynamical systems and numerical methods. These systems, while conceptually straightforward, embody essential characteristics of more complex dynamical systems and serve as valuable benchmarks for evaluating numerical methods and computational frameworks. With the rise of open science practices and increasing demands for research transparency [1], the choice of computational platforms becomes particularly significant in ensuring both mathematical rigor and reproducibility.

The mathematical structure of SMD systems manifests across multiple scales of computational physics and engineering analysis, from simple mechanical oscillators to sophisticated control systems. Our investigation encompasses both open-loop and closed-loop configurations, allowing us to examine how different numerical methods handle varying degrees of system complexity. In the current landscape, where open science platforms are reshaping scientific communication [2], implementing these systems in accessible, transparent environments becomes crucial for validating both theoretical predictions and numerical approximations. Two leading open-source platforms, Python and R, have emerged as dominant choices for scientific computing, each offering distinct approaches to numerical computation and algorithm implementation.

Contemporary scientific computing increasingly relies on these open-source platforms, which align with broader movements toward research transparency and reproducibility. Python’s scientific computing ecosystem, built around NumPy and SciPy, exemplifies how open-source tools can enable sophisticated numerical computations while maintaining accessibility. The platform’s implementation of variable-step, variable-order methods provides particular insight into handling stiff differential equations and adaptive error control. Similarly, R’s statistical computing framework, particularly through packages like deSolve, demonstrates how community-driven development can produce robust tools for complex mathematical modeling, offering unique perspectives on numerical stability and error propagation. The computational challenges become particularly evident when implementing variable-step, variable-order methods for these systems, especially in the context of PID control where system stiffness can vary significantly. As scientific practices evolve toward greater openness and accessibility [3], the ability to examine and validate implementation details becomes crucial. Both Python and R provide transparent implementations of numerical methods, allowing researchers to understand and verify the underlying algorithms - a key advantage for studying numerical stability, convergence properties, and error accumulation in long-time integration.

The significance of this comparative study extends beyond mere platform evaluation, addressing fundamental questions in computational mathematics and numerical analysis. In an era where research integrity increasingly depends on computational reproducibility [2], understanding the relative strengths and limitations of different open-source implementations becomes crucial. First, the choice between Python and R affects not only computational performance but also numerical accuracy and stability, particularly in handling stiff systems and adaptive step size control. Second, as open-source platforms continue to evolve, systematic comparisons help inform both theoretical understanding of numerical methods and their practical implementation. Third, with the growing emphasis on research transparency, these platforms provide an ideal foundation for studying how different numerical schemes behave in both open and closed-loop dynamical systems, contributing to both applied mathematics and computational science.

2. Methods

2.1. Mathematical Formulation

Our analysis began with an idealized one-dimensional SMD system, consisting of a point mass $m=100$ kg coupled to both an ideal linear spring of stiffness $k=50$ N/m and a viscous damper with coefficient $d=50$ Ns/m, as illustrated in Figure 1. The system’s idealization encompassed several key simplifications: the mass behaved as a point particle devoid of rotational dynamics, the spring exhibited perfectly linear behavior following Hooke’s law with negligible mass, the damper provided purely viscous damping, and the system operated in the absence of friction or additional constraints. All connections between components were considered rigid and massless, ensuring that the system’s behavior was governed solely by the interplay of inertial, elastic, and dissipative forces.

The mathematical model emerged directly from Newton’s Second Law, which stated that the sum of forces acting on the mass equaled the product of mass and acceleration. In our system, three distinct forces acted on the mass: the spring force following Hooke’s law ($F_s=-kx\left(t\right)$), the damping force proportional to velocity ($F_d=-d\dot{x}\left(t\right)$), and an external force $F\left(t\right)$. The negative signs in the spring and damping forces indicated their opposition to displacement and motion, respectively. Applying Newton’s Second Law and substituting these forces yielded the fundamental equation of motion:

$$m\ddot{x}\left(t\right)=F\left(t\right)-kx\left(t\right)-d\dot{x}\left(t\right),$$

(1)

which, upon substituting our system parameters, became:

$$100\ddot{x}\left(t\right)+50\dot{x}\left(t\right)+50x\left(t\right)=F\left(t\right).$$

(2)

To facilitate modern control analysis and implementation, we transformed this second-order differential equation into state-space form. The transformation introduced a state vector $\mathbf{x}\left(t\right)={[x_1\left(t\right),x_2\left(t\right)]}^T$, where $x_1$ represented position and $x_2$ represented velocity. This choice of states provided clear physical insight: $x_1$ captured the system’s configuration while $x_2$ described its instantaneous motion. The resulting state equations took the form:

$$\begin{array}{cc}\hfill {\dot{x}}_1& =x_2,\hfill \\ \hfill {\dot{x}}_2& =-0.5x_1-0.5x_2+0.01F\left(t\right),\hfill \end{array}$$

(3)

which could be expressed in the canonical matrix form:

$$\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\end{array}\right]=\underset{\mathbf{A}}{\underbrace{\left[\begin{array}{cc}0& 1\\ -0.5& -0.5\end{array}\right]}}\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]+\underset{\mathbf{B}}{\underbrace{\left[\begin{array}{c}0\\ 0.01\end{array}\right]}}F\left(t\right).$$

(4)

For precise position control, we implemented a Proportional-Integral-Derivative (PID) controller with carefully tuned gains: $K_p=200$ N/m for proportional action, $K_i=50$ N/(ms) for integral action, and $K_d=100$ N·s/m for derivative action. Given a desired position trajectory $x_d\left(t\right)$, the control law synthesized the input force as:

$$F\left(t\right)=K_p[x_d\left(t\right)-x_1\left(t\right)]+K_i{\int }_0^t[x_d\left(\tau \right)-x_1\left(\tau \right)]d\tau +K_d[x_d\left(t\right)-x_2\left(t\right)].$$

(5)

The integral action in the controller necessitated augmenting our state space with an additional state variable $x_3={\int }_0^t[x_d\left(\tau \right)-x_1\left(\tau \right)]d\tau$ representing the accumulated position error. This augmentation led to the complete closed-loop dynamics:

$$\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\\ {\dot{x}}_3\end{array}\right]=\left[\begin{array}{ccc}0& 1& 0\\ -2.5& -1.5& -0.5\\ -1& 0& 0\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]x_d\left(t\right).$$

(6)

The complete system architecture, including both open-loop and closed-loop configurations, was illustrated in Figure 2. This representation explicitly showed the signal flow and system interconnections, highlighting the fundamental difference between direct force input and feedback control strategies.

2.2. Numerical Implementation

The numerical solution of our SMD system required careful consideration of both accuracy and stability, particularly given the potential stiffness introduced by the control parameters. Our implementation leveraged the Python scientific computing ecosystem, specifically employing NumPy [4] for array operations and the Livermore Solver for Ordinary Differential Equations with Automatic method switching (LSODA) through SciPy’s integrate module [5]. For data manipulation and storage of results, we utilized the pandas library [6], which provided efficient DataFrame structures for handling time series data and simulation outputs.

The core of our numerical approach lay in the transformation of our continuous-time system into a form suitable for computational solution. For the open-loop system, we began with the state-space representation derived in our mathematical formulation. The solver required a system model function that returned the state derivatives:

$${\displaystyle \frac{d}{dt}}\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=f(\mathbf{x},t,F),$$

(7)

where $\mathbf{x}={[x_1,x_2]}^T$ represented our state vector, and F was the input force. This formulation expanded to our specific parameter values:

$$f(\mathbf{x},t,F)=\left[\begin{array}{c}{x}_2\\ {\displaystyle \frac{F-50x_1-50x_2}{100}}\end{array}\right].$$

(8)

The LSODA algorithm employed an adaptive strategy, automatically switching between methods optimal for stiff and non-stiff regions of the solution. For non-stiff regions, it utilized Adams-Moulton methods of orders 1-12, which took the general form:

$${\mathbf{x}}_{n+1}={\mathbf{x}}_n+h\sum _{i=0}^k{\beta }_i{\mathbf{f}}_{n-i},$$

(9)

where h represented the step size and ${\beta }_i$ were the method coefficients. For stiff regions, the solver switched to backward differentiation formulas (BDF):

$$\sum _{i=0}^k{\alpha }_i{\mathbf{x}}_{n+1-i}=h{\beta }_0{\mathbf{f}}_{n+1},$$

(10)

where ${\alpha }_i$ and ${\beta }_0$ were the BDF coefficients.

The simulation time domain was discretized with careful attention to numerical stability and accuracy requirements. We employed a base step size of $h=0.01$ seconds over a simulation interval $[0,100]$ seconds for the open-loop system and $[0,20]$ seconds for the controlled system. This discretization resulted in:

$$t_i=t_0+ih,i=0,1,\dots ,N,$$

(11)

where $N=(t_f-t_0)/h$ represented the total number of time steps.

The input force profile was generated through a carefully designed transition function to avoid numerical artifacts. For a step input from $F_1$ to $F_2$ occurring between times $t_1$ and $t_2$, we implemented a linear transition:

$$F\left(t\right)=\left\{\begin{array}{cc}{F}_1,\hfill & t<t_1\hfill \\ F_1+{\displaystyle \frac{F_2-F_1}{t_2-t_1}}(t-t_1),\hfill & t_1\le t\le t_2\hfill \\ F_2,\hfill & t>t_2\hfill \end{array}\right.,$$

(12)

with specific values $F_1=0$ N, $F_2=50$ N, $t_1=4$ s, and $t_2=5$ s for the open-loop case.

For the closed-loop system, the numerical implementation became more intricate due to the PID control law. The augmented state vector $\mathbf{x}={[x_1,x_2,x_3]}^T$ required modification of our system model to incorporate the integral term:

$$f(\mathbf{x},t,x_d)=\left[\begin{array}{c}{x}_2\\ {\displaystyle \frac{200(x_d-x_1)+50x_3-100x_2-50x_1-50x_2}{100}}\\ x_d-x_1\end{array}\right].$$

(13)

The numerical solver handled this augmented system with the same adaptive strategy, but special care was taken in the initialization and update of the integral term to prevent numerical drift. The desired position trajectory $x_d\left(t\right)$ followed a similar smooth transition profile, but with values scaled to position units: $x_1=0$ m to $x_2=1$ m between $t_1=2$ s and $t_2=2.5$ s.

Error control in our numerical solution was maintained through careful monitoring of local truncation error at each step:

$${\mathrm{LTE}}_n=\parallel {\mathbf{x}}_{n+1}-{\tilde{\mathbf{x}}}_{n+1}\parallel \le \mathrm{tol},$$

(14)

where ${\tilde{\mathbf{x}}}_{n+1}$ represented a higher-order approximation used for error estimation. The step size was adaptively adjusted based on this error estimate:

$$h_{\mathrm{new}}=h_{\mathrm{old}}{\left({\displaystyle \frac{\mathrm{tol}}{{\mathrm{LTE}}_n}}\right)}^{1/p},$$

(15)

with p denoting the current method order and $\mathrm{tol}={10}^{-8}$ as our specified tolerance.

The complete implementation was organized into modular Python classes, following object-oriented programming principles, which facilitated easy modification and extension of the system model and control algorithms. All numerical computations were performed using Python, with NumPy, SciPy, and pandas. For comparative validation and robustness testing, parallel implementations were developed in R using the deSolve package [7] for differential equation solving and the tidyverse ecosystem [8] for data manipulation and visualization. The results from both implementations showed excellent agreement, providing confidence in the numerical accuracy of our solutions.

2.3. Stability Analysis

Our stability analysis methodology incorporated multiple theoretical frameworks to thoroughly assess both open-loop and closed-loop system behavior. Following [9], we employed three complementary approaches of eigenvalue analysis, Routh-Hurwitz criterion, and Bounded-Input, Bounded-Output (BIBO) stability assessment, with each method providing unique insights into the system’s dynamic characteristics.

The first stage of our analysis centered on eigenvalue decomposition. For a linear time-invariant system, stability required that all eigenvalues had negative real parts [10]. The analysis began with our open-loop system matrix:

$$\mathbf{A}=\left[\begin{array}{cc}0& 1\\ -{\displaystyle \frac{k}{m}}& -{\displaystyle \frac{d}{m}}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -0.5& -0.5\end{array}\right],$$

(16)

from which we formulated the characteristic equation through determinant evaluation:

$$det(s\mathbf{I}-\mathbf{A})=\left|\begin{array}{cc}s& -1\\ 0.5& s+0.5\end{array}\right|=s^2+{\displaystyle \frac{d}{m}}s+{\displaystyle \frac{k}{m}}=s^2+0.5s+0.5=0.$$

(17)

Following [11], we determined two fundamental parameters characterizing the system’s behavior. The natural frequency ${\omega }_n$ represented the system’s undamped oscillation rate, while the damping ratio $\zeta$ characterized the decay rate of oscillations:

$${\omega }_n=\sqrt{{\displaystyle \frac{k}{m}}}=\sqrt{0.5},\zeta ={\displaystyle \frac{d}{2\sqrt{mk}}}={\displaystyle \frac{50}{2\sqrt{100·50}}}=0.354.$$

(18)

The nature of system damping played a crucial role in determining the transient response characteristics. As derived by [9], the damping ratio $\zeta =0.354$ indicated an underdamped system, as it fell within the range $0<\zeta <1$. In such cases, the system’s free response exhibited oscillatory behavior with exponential decay. The complete solution took the form:

$$x\left(t\right)=A{e}^{-\zeta {\omega }_{n}t}cos({\omega }_{d}t+\varphi ),$$

(19)

where ${\omega }_d={\omega }_n\sqrt{1-{\zeta }^2}$ represented the damped natural frequency, and A and $\varphi$ were determined by initial conditions. This underdamped characteristic arose from our physical parameters, with the mass-spring-damper combination producing complex conjugate eigenvalues ${\lambda }_{1,2}=-\zeta {\omega }_n\pm i{\omega }_d=-0.25\pm 0.661i$.

Asymptotic stability analysis followed the theoretical framework of [10], which required examination of the system’s long-term behavior as $t\to \infty$. For our linear system, asymptotic stability demanded that all eigenvalues had strictly negative real parts. The real component of our eigenvalues, $-\zeta {\omega }_n=-0.25$, being negative, ensured exponential decay of any initial conditions or disturbances. This stability manifested through the exponential term $e^{-\zeta {\omega }_{n}t}$ in the system response, which approached zero as time increased.

The Routh-Hurwitz stability criterion [12] provided an algebraic method for determining stability without explicit eigenvalue calculation. For our second-order characteristic equation:

$$a_2{s}^2+a_{1}s+a_0=s^2+0.5s+0.5=0,$$

(20)

we constructed the Routh array and examined its first column entries:

$$\begin{array}{ccc}{s}^2& 1& 0.5\\ s^1& 0.5& 0\\ s^0& 0.5\end{array}.$$

(21)

The stability conditions derived from the Routh array required all elements in its first column to be positive. Following [13], this translated to the requirements $a_2>0$, $a_1>0$, and $a_0>0$, which were all satisfied in our system, confirming stability from an algebraic perspective.

The BIBO stability assessment followed [14] methodology, focusing on the system’s transfer function:

$$G\left(s\right)={\displaystyle \frac{1}{m{s}^2+ds+k}}={\displaystyle \frac{0.01}{s^2+0.5s+0.5}}.$$

(22)

This approach required examination of both the positivity of denominator coefficients and the Hurwitz stability of the denominator polynomial. Our transfer function satisfied both conditions, as all denominator coefficients were positive and matched the previously analyzed characteristic equation.

For the closed-loop analysis with PID control, we extended these methodologies to the augmented third-order system. The characteristic equation expanded to include controller parameters:

$$s^3+{\displaystyle \frac{d+K_d}{m}}{s}^2+{\displaystyle \frac{k+K_p}{m}}s+{\displaystyle \frac{K_i}{m}}=0,$$

(23)

necessitating an expanded Routh array for stability analysis:

$$\begin{array}{ccc}{s}^3& 1& {\displaystyle \frac{k+K_p}{m}}\\ s^2& {\displaystyle \frac{d+K_d}{m}}& {\displaystyle \frac{K_i}{m}}\\ s^1& {\displaystyle \frac{(d+K_d)(k+K_p)-m{K}_i}{m(d+K_d)}}& 0\\ s^0& {\displaystyle \frac{K_i}{m}}\end{array}.$$

(24)

The relationship between damping and asymptotic stability was further illuminated through the root locus analysis detailed by [13]. The characteristic equation roots, or system poles, lying in the left half of the complex plane (Re$\left(\lambda \right)<0$) guaranteed that any oscillations would eventually decay. The distance of these poles from the imaginary axis, determined by $\zeta {\omega }_n$, dictated the decay rate, while their imaginary components $\pm {\omega }_d$ determined the oscillation frequency.

The robustness of our closed-loop system was examined through the method of [15], which analyzed the loop transfer function:

$$L\left(s\right)={\displaystyle \frac{K_p{s}^2+K_{d}s+K_i}{s(m{s}^2+ds+k)}}.$$

(25)

This transfer function formed the basis for determining critical frequency-domain stability margins, providing insight into the system’s robustness against parameter variations and modeling uncertainties. The stability margins derived from this transfer function offered quantitative measures of how much gain variation and phase shift the system could tolerate while maintaining stability.

2.4. Statistical Framework for Performance Analysis

Building upon the system dynamics simulation framework described before, we developed and implemented a comprehensive statistical analysis approach to compare the performance characteristics of Python and R implementations. All experiments were conducted on a ThinkPad P52s laptop running Fedora Linux 39 (Budgie) x86_64, equipped with an Intel i7-8550U processor (8 threads, 4.000GHz). Given that our numerical simulations primarily involved CPU-bound computations without parallel processing requirements, this standard laptop configuration provided an adequate testing environment representative of typical scientific computing setups in practice.

Our analysis pipeline utilized Python’s scientific computing stack, including NumPy for numerical computations [4], SciPy for statistical testing [5], pandas for data management [6], and Seaborn for statistical visualization [16]. Following established practices in performance evaluation [17] and statistical sampling theory [18], we collected $n=1000$ measurements for each implementation type $i\in \mathrm{controlled},\mathrm{open}$ and language $j\in \mathrm{Python},\mathrm{R}$ after discarding $k=5$ warm-up runs. These warm-up runs were essential to mitigate the effects of just-in-time compilation and cache warming, as demonstrated by Herho et al. [19], Herho et al. [20], Herho et al. [21].

Our analysis began with fundamental statistical measures for each implementation-language pair $(i,j)$. Following the recommendations of Jain [22] for computer systems performance analysis, we computed the central tendency and dispersion metrics:

$${\mu }_{i,j}={\displaystyle \frac{1}{n}}\sum _{k=1}^n{X}_{i,j,k},$$

(26)

where ${\mu }_{i,j}$ represents the mean performance measure, n is the number of observations, and $X_{i,j,k}$ is the k-th performance measurement for implementation i and language j.

$${\chi }_{i,j}=\sqrt{{\displaystyle \frac{1}{n-1}}\sum _{k=1}^n{(X_{i,j,k}-{\mu }_{i,j})}^2},$$

(27)

where ${\chi }_{i,j}$ denotes the standard deviation of the performance measurements.

$${\eta }_{i,j}={\displaystyle \frac{{\chi }_{i,j}}{{\mu }_{i,j}}}\times 100\%,$$

(28)

where ${\eta }_{i,j}$ represents the coefficient of variation.

The coefficient of variation in Equation (28) proved particularly valuable for performance comparisons, as it provided a dimensionless measure of variability [23]. Following [24] and recent performance analysis guidelines [25], we interpreted $\eta$ values below 5% as very consistent performance, 5–10% as good consistency, 10–15% as moderate variability, and above 15% as high variability.

Distribution shape analysis was crucial for understanding performance characteristics and selecting appropriate statistical tests [26]. We computed skewness (${\gamma }_1$) and kurtosis (${\gamma }_2$):

$${\gamma }_1={\displaystyle \frac{1}{n}}\sum _{k=1}^n{\displaystyle \frac{{(X_{i,j,k}-{\mu }_{i,j})}^3}{{\chi }_{i,j}^3}},$$

(29)

where ${\gamma }_1$ represents the skewness of the distribution.

$${\gamma }_2={\displaystyle \frac{1}{n}}\sum _{k=1}^n{\displaystyle \frac{{(X_{i,j,k}-{\mu }_{i,j})}^4}{{\chi }_{i,j}^4}}-3,$$

(30)

where ${\gamma }_2$ represents the excess kurtosis of the distribution.

These shape parameters were essential for detecting potential performance anomalies and understanding the underlying performance distribution patterns [27].

The normality of distributions was assessed using both Shapiro-Wilk and Kolmogorov-Smirnov tests, as recommended by Razali and Wah [28] for their complementary strengths in different sample sizes and distribution types:

$$W={\displaystyle \frac{{\left({\sum }_{i=1}^n{a}_i{x}_{\left(i\right)}\right)}^2}{{\sum }_{i=1}^n{(x_i-\overline{x})}^2}},$$

(31)

where W is the Shapiro-Wilk test statistic, $x_{\left(i\right)}$ represents ordered sample values, $a_i$ are constants derived from the covariance matrix of order statistics, and $\overline{x}$ is the sample mean.

$$D_n=\underset{x}{sup}|F_n\left(x\right)-F\left(x\right)|,$$

(32)

where $D_n$ is the Kolmogorov-Smirnov test statistic, $F_n\left(x\right)$ is the empirical distribution function, and $F\left(x\right)$ is the theoretical normal distribution.

For comparing Python and R implementations, we employed non-parametric tests due to their robustness against violations of normality assumptions [29]. The Wilcoxon signed-rank test [30] was used for paired comparisons:

$$W=\sum _{r=1}^{n_r}[\mathrm{sgn}(x_{2,r}-x_{1,r})·R_r],$$

(33)

where W is the test statistic, $n_r$ is the number of non-zero differences, sgn is the sign function, and $R_r$ is the rank of the absolute difference.

This was complemented by Levene’s test [31] for assessing variance homogeneity, chosen for its robustness against non-normality [32]:

$$\Lambda ={\displaystyle \frac{(N-k)}{(k-1)}}{\displaystyle \frac{{\sum }_{i=1}^k{N}_i{(Z_{i.}-Z_{..})}^2}{{\sum }_{i=1}^k{\sum }_{j=1}^{N_i}{(Z_{ij}-Z_{i.})}^2}},$$

(34)

where $\Lambda$ is the test statistic, N is the total sample size, k is the number of groups, $N_i$ is the sample size of the i-th group, $Z_{ij}=|X_{ij}-{\tilde{X}}_i|$ with ${\tilde{X}}_i$ being the group median, $Z_{i.}$ represents group means of the $Z_{ij}$, and $Z_{..}$ is the overall mean of $Z_{ij}$.

Following modern statistical practice [33], we complemented significance testing with effect size analysis using Cohen’s d [34]:

$$\delta ={\displaystyle \frac{{\mu }_{\mathrm{Python}}-{\mu }_{\mathrm{R}}}{\sqrt{\frac{{\chi }_{\mathrm{Python}}^2+{\chi }_{\mathrm{R}}^2}{2}}}},$$

(35)

where $\delta$ is Cohen’s d effect size measure, ${\mu }_{\mathrm{Python}}$ and ${\mu }_{\mathrm{R}}$ are the means of Python and R implementations respectively, and ${\chi }_{\mathrm{Python}}^2$ and ${\chi }_{\mathrm{R}}^2$ are their respective variances.

This measure provided a standardized assessment of practical significance, with effect sizes interpreted according to established guidelines [35] as negligible ($\left|\delta \right|<0.2$), small ($0.2\le \left|\delta \right|<0.5$), medium ($0.5\le \left|\delta \right|<0.8$), or large ($\left|\delta \right|\ge 0.8$).

Performance anomalies were identified using a robust multi-method approach [36], combining three complementary detection methods as recommended by Chandola et al. [37]. The traditional Z-score method:

$$\left|\zeta \right|=\left|{\displaystyle \frac{x-\mu }{\chi }}\right|>3,$$

(36)

where $\left|\zeta \right|$ is the absolute Z-score, x is the observation value, $\mu$ is the mean, and $\chi$ is the standard deviation.

This was supplemented by the modified Z-score using median absolute deviation, which offered greater robustness against multiple outliers [38]:

$$|{\zeta }_{\mathrm{mod}}|=\left|{\displaystyle \frac{0.6745(x-\mathrm{median})}{\mathrm{MAD}}}\right|>3.5,$$

(37)

where $|{\zeta }_{\mathrm{mod}}|$ is the modified Z-score, MAD is the median absolute deviation, and $0.6745$ is the normalization constant.

The interquartile range method provided a distribution-free approach [39]:

$$x<Q_1-1.5\mathrm{IQR}\mathrm{or}x>Q_3+1.5\mathrm{IQR},$$

(38)

where $Q_1$ and $Q_3$ are the first and third quartiles respectively, and IQR is the interquartile range.

The final set of anomalies $\mathcal{A}$ was determined through consensus voting [40]:

$$\mathcal{A}={\mathcal{A}}_{\zeta }\cup {\mathcal{A}}_{\mathrm{mod}}\cup {\mathcal{A}}_{\mathrm{IQR}},$$

(39)

where ${\mathcal{A}}_{\zeta }$, ${\mathcal{A}}_{\mathrm{mod}}$, and ${\mathcal{A}}_{\mathrm{IQR}}$ represent the sets of anomalies identified by the Z-score, modified Z-score, and IQR methods respectively.

This statistical framework described in Equations (26)–(39) provided a rigorous foundation for comparing performance characteristics between implementations, adhering to established practices in both statistical methodology and performance analysis [41]. All numerical results were standardized to three decimal places for consistency in technical reporting, as recommended by Altman et al. [42] for reproducible research practices.

3. Results and Discussion

3.1. System Dynamics

The SMD system exhibited complex dynamic behavior characteristic of underdamped second-order mechanical systems, as evidenced by both mathematical stability analysis and numerical simulations. Figure 3 illustrates the system’s temporal evolution under a step force input, while Figure 4 provides insight into the system’s state-space behavior through its phase portrait.

The eigenvalue analysis revealed complex conjugate poles at $-0.25\pm 0.661i$, mathematically confirming the system’s underdamped nature. This finding aligns with the calculated damping ratio $\zeta =0.354$, indicating oscillatory behavior with gradual energy dissipation. The characteristic equation coefficients $[1,0.5,0.5]$ satisfied stability criteria across multiple analytical frameworks, including eigenvalue, Routh-Hurwitz, and BIBO stability analyses.

As shown in Figure 3, the system response to the step input force demonstrates characteristic underdamped behavior. The position response (Figure 3b) shows an initial overshoot followed by decaying oscillations before settling to its final value of 1 m. The velocity profile (Figure 3c) exhibits a maximum speed of 0.4 m/s during the initial response phase, with subsequent oscillations diminishing as the system approaches equilibrium. The acceleration response (Figure 3d) shows sharp transitions corresponding to the force application, with peak values of approximately ±0.2 m/s². The phase portrait in Figure 4 provides additional insight into the system’s dynamic behavior. The spiral trajectory illustrates the system’s evolution from its initial state at the origin to its final equilibrium point at (1 m, 0 m/s). The decreasing radius of the spiral confirms the system’s stability and energy dissipation through damping, while the elliptical shape of the trajectory reflects the continuous exchange between potential and kinetic energy characteristic of spring-mass systems. The simulation results confirm the theoretical stability predictions, with the system achieving its expected final state: position at 1 m, velocity at 0 m/s, and acceleration at 0 m/s².

The implementation of PID control significantly enhanced the system’s performance characteristics, introducing more sophisticated dynamic behavior while maintaining stability. Figure 5 presents the comprehensive closed-loop system response, while Figure 6 illustrates the controlled system’s state-space behavior through its phase portrait.

The closed-loop stability analysis revealed a richer pole constellation compared to the open-loop system, with eigenvalues at $-0.5$ and $-0.25\pm 0.661i$. The additional real pole at $-0.5$, introduced by the controller, enhanced the system’s damping characteristics, as evidenced by the increased damping ratio of $\zeta =0.474$ compared to the open-loop value of 0.354. The characteristic equation coefficients $[1,1.5,2.5,0.5]$ maintained stability across all analytical criteria while enabling improved transient response.

The controlled system demonstrated superior performance metrics, as shown in Figure 5. The position tracking (Figure 5a) achieved a steady-state value of 0.996 m, resulting in a final position error of only 0.004 m (0.4%). The maximum overshoot of 0.107 m (10.7%) occurred during the initial transient response, while the settling time of 12.8 seconds represented a significant improvement over the open-loop behavior. The control force profile (Fig. Figure 5b) showed an initial peak of 259.2 N, necessary for rapid response, before settling to a steady-state value that maintained the desired position.

The phase portrait in Figure 6 reveals the effectiveness of the control strategy in managing the system’s state trajectory. Unlike the open-loop case, the controlled system exhibits a more direct path to the desired equilibrium point, with the controller actively shaping the system’s dynamic behavior. The single, well-defined trajectory indicates effective energy management by the controller, avoiding the multiple oscillations characteristic of the uncontrolled response.

The RMS position error of 0.124 m provides a quantitative measure of tracking performance throughout the entire operation. The maximum speed of 1.04 m/s and acceleration of 2.14 m/s² demonstrate the controller’s ability to generate aggressive yet controlled motions when necessary while maintaining system stability. These performance metrics validate the theoretical design of the PID controller with gains $K_p=200$ N/m, $K_i=50$ N/(ms), and $K_d=100$ Ns/m.

3.2. Performance Analysis

Our statistical analysis revealed nuanced and significant differences between Python and R implementations across all metrics and scenarios, evaluated through 1,000 iterations per test. The results show a complex interplay of execution time, memory usage, and stability across controlled and open-loop systems, highlighting the strengths and weaknesses of each platform in scientific computing workflows.

Python demonstrated substantial advantages in execution time performance. In controlled systems, Python achieved a mean execution time of 2.037 seconds ($\sigma =0.337$), significantly outperforming R, which recorded a mean of 5.591 seconds ($\sigma =0.514$). This reduction of 63.566% was statistically significant, confirmed by the Wilcoxon signed-rank test ($p<0.001$, $d=-8.175$). However, Python’s performance came at the cost of higher variability (CV = 16.553%, categorized as "High variability"), while R showed more consistent execution times (CV = 9.201%, "Good consistency"). Figure 7, left panel, illustrates this contrast with Python exhibiting a broader distribution than R.

In open-loop implementations, Python retained its performance edge with a mean execution time of 3.103 seconds ($\sigma =0.473$) compared to R’s 4.223 seconds ($\sigma =0.667$). Although the margin of improvement was reduced to 26.517%, the difference remained statistically significant ($p<0.001$, $d=-1.939$). Both platforms exhibited high variability in this scenario (Python CV = 15.227%, R CV = 15.784%). Levene’s test further confirmed significant variance differences ($F=13.157$, $p<0.001$). The right panel of Figure 7 highlights these trends, showing Python’s faster, albeit more variable, performance relative to R.

Memory usage patterns presented a contrasting narrative. In controlled systems, Python consumed slightly more memory, with a mean of 125.461 MB ($\sigma =0.217$) compared to R’s 123.583 MB ($\sigma =0.139$), marking a 1.519% increase. This difference was statistically significant ($p<0.001$) with a large effect size ($d=10.305$). However, both platforms demonstrated exceptional stability, as reflected by CV values of 0.173% for Python and 0.112% for R. These results suggest that Python’s higher baseline interpreter memory footprint accounts for the slight difference.

In open-loop systems, memory usage between the platforms converged. Python recorded a mean of 132.965 MB ($\sigma =0.224$), while R reached 133.020 MB ($\sigma =0.112$). The negligible difference of -0.042% was statistically significant ($p<0.001$), though with a small effect size ($d=-0.312$), indicating practical equivalence. Both implementations maintained highly consistent memory utilization, with CV values below 0.2%. Figure 8, right panel, illustrates this near-parity in memory usage distributions.

Distribution shape analysis provided additional insights. Execution time data showed significant right skewness in both platforms, with Python exhibiting skewness $\ge 1.094$ and R exceeding 3.369 in controlled implementations. Non-normality tests (Shapiro-Wilk $p<0.001$) confirmed deviations from Gaussian behavior, as shown by broader tails in Figure 7. Memory usage data, on the other hand, displayed symmetric or moderately skewed distributions depending on the scenario, reflecting greater consistency in resource allocation.

Anomaly counts varied significantly across configurations. For controlled execution times, Python recorded a high anomaly rate (33.0% of data) compared to R’s 10.6%, reflecting the variability induced by Python’s dynamic optimization strategies. Conversely, in open-loop memory usage, R exhibited more anomalies (1.8%) compared to Python’s 0.4%, underscoring R’s sensitivity to subtle system conditions in memory-intensive scenarios.

Overall, Python’s execution speed makes it ideal for computationally intensive tasks, especially in controlled systems where performance is critical. However, its variability and higher anomaly rates suggest careful consideration for applications requiring extreme precision and stability. R’s strength lies in its consistent memory usage and lower variability, making it a better choice for scenarios where stability and resource predictability are paramount. These results underscore the importance of aligning platform selection with specific computational requirements to optimize performance and efficiency.

For example, Masini and Bientinesi [43] emphasizes Python’s broader ecosystem and optimization strategies as critical for scalable applications, while Raschka et al. [44] highlight R’s robustness for statistical workflows. These findings align with earlier research by Watson et al. [45], suggesting platform-specific trade-offs in computational reliability.

4. Conclusion

This study provides a comprehensive evaluation of Python and R for simulating and controlling SMD systems, with a focus on computational performance and numerical accuracy. Python exhibited a clear advantage in execution speed, achieving reductions of up to 63.57% in mean runtime compared to R in controlled system simulations. However, this efficiency came with trade-offs, including higher variability in execution times and a greater anomaly rate, reflecting the platform’s dynamic optimization strategies. Conversely, R demonstrated superior consistency in execution and memory usage, making it more suitable for scenarios prioritizing stability and resource predictability. These findings underscore the distinct strengths of each platform and the importance of selecting the right tool for specific computational needs.

The implementation of PID control significantly enhanced system dynamics, yielding improved position tracking and stability. With a final position error of just 0.4% and enhanced damping characteristics, the controlled system achieved superior performance metrics over the open-loop configuration. The controlled system’s steady-state precision and reduced settling time validate the PID controller’s design and its effectiveness in managing transient responses. Overall, this work bridges theoretical stability analysis with empirical performance insights, providing a valuable resource for researchers in computational dynamics and control systems. By leveraging the capabilities of open-source platforms, this research promotes reproducibility and transparency in the numerical study of dynamic systems.

Code and Data Availability. All code used for numerical simulations, statistical analyses, and the complete dataset generated during this study are openly available in the GitHub repository at https://github.com/sandyherho/smdCompare.

Funding. This work was supported by the Dean’s Distinguished Fellowship, University of California, Riverside, in 2023.