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

Sandy Herho, Siti Kaban

► To cite this version:

Sandy Herho, Siti Kaban. Quantitative Performance Analysis of Spring-Mass-Damper Control Systems: A Comparative Implementation in Python and R. 2025. hal-04860538

HAL Id: hal-04860538

https://hal.science/hal-04860538v1

Preprint submitted on 1 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Quantitative Performance Analysis of Spring-Mass-Damper Control Systems: A Comparative Implementation in Python and R

Sandy H. S. Herho<fnref n="1" />* and Siti N. Kaban<fnref n="2" />

## Abstract

The numerical simulation of control of spring-mass-damper (SMD) systems offer critical insights into dynamical systems and computational methodologies. This study provides a comprehensive comparative analysis of implementing SMD systems across two prominent open-source scientific computing platforms: Python and R. By examining both open-loop and closed-loop system configurations, the research investigates the computational performance, numerical accuracy, and implementation characteristics of these platforms. Utilizing an idealized one-dimensional SMD system with a Proportional-Integral-Derivative (PID) controller, the study conducted extensive numerical simulations and statistical performance analyses. Results revealed Python's significant advantages in execution speed, achieving up to 63.57% reduction in runtime for controlled system simulations, while R demonstrated superior consistency in execution and memory usage. The controlled system demonstrated exceptional performance, with a final position error of merely 0.4% and enhanced damping characteristics. This work not only bridges theoretical stability analysis with empirical performance insights but also promotes reproducibility and transparency in computational dynamics research by leveraging open-source platforms.

**Keywords:** Control systems, Cross-platform implementation, Numerical methods, Performance analysis, Spring-mass-damper dynamics

<footnote n="1">¹*Department of Earth and Planetary Sciences, University of California, 900 University Ave., Riverside, 92521, CA, USA.</footnote>

<footnote n="2">²Offshore Engineering Research Group, Bandung Institute of Technology (ITB), Jalan Ganesha 10, Bandung, 40132, West Java, Indonesia.</footnote>

<footnote>*Corresponding author(s). E-mail(s): sandy.herho@email.ucr.edu;</footnote>
# 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.
75
Third, with the growing emphasis on research transparency, these platforms provide
76
an ideal foundation for studying how different numerical schemes behave in both open
77
and closed-loop dynamical systems, contributing to both applied mathematics and
78
computational science.
79

## 2 Methods

### 2.1 Mathematical Formulation

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

<figure><img src="image_2.png" /><figcaption>**Fig. 1** Physical representation of the idealized SMD system showing coordinate system and force components. The origin was set at the spring's equilibrium position, with positive displacement defined in the rightward direction.</figcaption></figure>

The mathematical model emerged directly from Newton's Second Law, which
92
stated that the sum of forces acting on the mass equaled the product of mass and
93
acceleration. In our system, three distinct forces acted on the mass: the spring force
94
following Hooke's law ($F_s = -kx(t)$), the damping force proportional to velocity ($F_d = -d\dot{x}(t)$), and an external force $F(t)$. 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}(t) = F(t) - kx(t) - d\dot{x}(t), \quad (1)$$

which, upon substituting our system parameters, became:

$$100\ddot{x}(t) + 50\dot{x}(t) + 50x(t) = F(t). \quad (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}(t) = [x_1(t), x_2(t)]^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{aligned} \dot{x}_1 &= x_2, \\ \dot{x}_2 &= -0.5x_1 - 0.5x_2 + 0.01F(t), \end{aligned} \quad (3)$$

which could be expressed in the canonical matrix form:

$$\begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} = \underbrace{\begin{bmatrix} 0 & 1 \\ -0.5 & -0.5 \end{bmatrix}}_{\text{A}} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} + \underbrace{\begin{bmatrix} 0 \\ 0.01 \end{bmatrix}}_{\text{B}} F(t). \quad (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$ Ns/m for derivative action. Given a desired position trajectory $x_d(t)$, the control law synthesized the input force as:

$$F(t) = K_p[x_d(t) - x_1(t)] + K_i \int_0^t [x_d(\tau) - x_1(\tau)]d\tau + K_d[x_d(t) - x_2(t)]. \quad (5)$$

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

$$\begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \\ \dot{x}_3 \end{bmatrix} = \begin{bmatrix} 0 & 1 & 0 \\ -2.5 & -1.5 & -0.5 \\ -1 & 0 & 0 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} x_d(t). \quad (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.
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Block diagram representation of the SMD system: Open-loop configuration showing direct force input with system transfer function and closed-loop configuration illustrating PID control architecture with complete signal flow paths and transfer functions.</figcaption></figure>

## 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 Frame 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.
131 The solver required a system model function that returned the state derivatives:

$$ \frac{d}{dt} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = f(\mathbf{x}, t, F), \qquad (7) $$

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

$$ f(\mathbf{x}, t, F) = \begin{bmatrix} x_2 \\ \frac{F-50x_1-50x_2}{100} \end{bmatrix}. \qquad (8) $$

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

$$ \mathbf{x}_{n+1} = \mathbf{x}_n + h \sum_{i=0}^{k} \beta_i \mathbf{f}_{n-i}, \qquad (9) $$

138 where $h$ represented the step size and $\beta_i$ were the method coefficients. For stiff regions, 139 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}, \qquad (10) $$

140 where $\alpha_i$ and $\beta_0$ were the BDF coefficients.

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

$$ t_i = t_0 + ih, \quad i = 0, 1, \dots, N, \qquad (11) $$

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

146 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 147 times $t_1$ and $t_2$, we implemented a linear transition:

$$ F(t) = \begin{cases} F_1, & t < t_1 \\ F_1 + \frac{F_2-F_1}{t_2-t_1}(t-t_1), & t_1 \le t \le t_2, \\ F_2, & t > t_2 \end{cases} \qquad (12) $$

148 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.

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

$$f(\mathbf{x}, t, x_d) = \left[ \frac{x_2}{\frac{200(x_d - x_1) + 50x_3 - 100x_2 - 50x_1 - 50x_2}{100}} \right]. \quad (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(t)$ 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:

$$LTE_n = \|\mathbf{x}_{n+1} - \tilde{\mathbf{x}}_{n+1}\| \leq tol, \quad (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_{\text{new}} = h_{\text{old}} \left( \frac{\text{tol}}{LTE_n} \right)^{1/p}, \quad (15)$$

with $p$ denoting the current method order and 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} = \begin{bmatrix} 0 & 1 \\ -\frac{k}{m} & -\frac{d}{m} \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -0.5 & -0.5 \end{bmatrix}, \quad (16)$$
from which we formulated the characteristic equation through determinant evaluation:

$$ \det(s\mathbf{I} - \mathbf{A}) = \begin{vmatrix} s & -1 \\ 0.5 & s + 0.5 \end{vmatrix} = s^2 + \frac{d}{m}s + \frac{k}{m} = s^2 + 0.5s + 0.5 = 0. \quad (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{\frac{k}{m}} = \sqrt{0.5}, \quad \zeta = \frac{d}{2\sqrt{mk}} = \frac{50}{2\sqrt{100 \cdot 50}} = 0.354. \quad (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(t) = A e^{-\zeta \omega_n t} \cos(\omega_d t + \phi), \quad (19) $$

where $\omega_d = \omega_n \sqrt{1 - \zeta^2}$ represented the damped natural frequency, and A and $\phi$ 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, \quad (20) $$

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

$$ \begin{vmatrix} s^2 & 1 & 0.5 \\ s^1 & 0.5 & 0 \\ s^0 & 0.5 & 0 \end{vmatrix}. \quad (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(s) = \frac{1}{ms^2 + ds + k} = \frac{0.01}{s^2 + 0.5s + 0.5} \quad (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 + \frac{d + K_d}{m}s^2 + \frac{k + K_p}{m}s + \frac{K_i}{m} = 0, \quad (23)$$

necessitating an expanded Routh array for stability analysis:

$$\left| \begin{array}{ccc} s^3 & \frac{k+K_p}{m} & 0 \\ s^2 & \frac{d+K_d}{m} & 0 \\ s^1 & \frac{(d+K_d)(k+K_p)-mK_i}{m(d+K_d)} & 0 \\ s^0 & \frac{K_i}{m} & 0 \end{array} \right. . \quad (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 ($\text{Re}(\lambda) < 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(s) = \frac{K_p s^2 + K_d s + K_i}{s(ms^2 + ds + k)}. \quad (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$ controlled, open and language $j \in$ Python, 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, 20, 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} = \frac{1}{n} \sum_{k=1}^{n} X_{i,j,k}, \qquad (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{\frac{1}{n-1} \sum_{k=1}^{n} (X_{i,j,k} - \mu_{i,j})^2}, \qquad (27) $$

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

$$ \eta_{i,j} = \frac{\chi_{i,j}}{\mu_{i,j}} \times 100\%, \qquad (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 = \frac{1}{n} \sum_{k=1}^{n} \frac{(X_{i,j,k} - \mu_{i,j})^3}{\chi_{i,j}^3}, \qquad (29) $$

where $\gamma_1$ represents the skewness of the distribution.

$$ \gamma_2 = \frac{1}{n} \sum_{k=1}^{n} \frac{(X_{i,j,k} - \mu_{i,j})^4}{\chi_{i,j}^4} - 3, \qquad (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 = \frac{\left(\sum_{i=1}^{n} a_i x_{(i)}\right)^2}{\sum_{i=1}^{n} (x_i - \bar{x})^2}, \qquad (31)$$

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

$$D_n = \sup_x |F_n(x) - F(x)|, \qquad (32)$$

where $D_n$ is the Kolmogorov-Smirnov test statistic, $F_n(x)$ is the empirical distribution function, and $F(x)$ 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} [\operatorname{sgn}(x_{2,r} - x_{1,r}) \cdot R_r], \qquad (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 = \frac{(N-k) \sum_{i=1}^{k} N_i (Z_i - Z_{i..})^2}{(k-1) \sum_{i=1}^{k} \sum_{j=1}^{N_i} (Z_{ij} - Z_{i..})^2}, \qquad (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_{i..}$ 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 = \frac{\mu_{\text{Python}} - \mu_{\text{R}}}{\sqrt{\frac{\chi_{\text{Python}}^2 + \chi_{\text{R}}^2}{2}}}, \qquad (35)$$

where $\delta$ is Cohen's d effect size measure, $\mu_{\text{Python}}$ and $\mu_{\text{R}}$ are the means of Python and R implementations respectively, and $\chi_{\text{Python}}^2$ and $\chi_{\text{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 ($|\delta| < 0.2$), small ($0.2 \le |\delta| < 0.5$), medium ($0.5 \le |\delta| < 0.8$), or large ($|\delta| \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:
$$|\zeta| = \left| \frac{x - \mu}{\chi} \right| > 3, \qquad (36)$$

where $|\zeta|$ 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_{\text{mod}}| = \left| \frac{0.6745(x - \text{median})}{\text{MAD}} \right| > 3.5, \qquad (37)$$

where $|\zeta_{\text{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\text{IQR} \text{ or } x > Q_3 + 1.5\text{IQR}, \qquad (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}_{\text{mod}} \cup \mathcal{A}_{\text{IQR}}, \qquad (39)$$

where $\mathcal{A}_{\zeta}$, $\mathcal{A}_{\text{mod}}$, and $\mathcal{A}_{\text{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.
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Temporal evolution of the open-loop system response to a 50 N step input. From top to bottom: (a) Applied force showing the transition from 0 to 50 N between t = 4 s and t = 5 s, (b) Position response demonstrating underdamped oscillations converging to 1.000 m, (c) Velocity profile reaching a maximum of 0.4 m/s before dampened oscillations to zero, and (d) Acceleration response highlighting the system's dynamic behavior with peak values of ± 0.2 m/s².</figcaption></figure>
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.

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Phase portrait of the open-loop system showing the relationship between position and velocity. The spiral trajectory, beginning from the origin and converging to the equilibrium point (1 m, 0 m/s), illustrates the system's stable focus behavior. The decreasing spiral radius demonstrates the energy dissipation through damping, while the elliptical shape reflects the interchange between potential and kinetic energy.</figcaption></figure>

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 $\pm 0.2$ m/s$^2$. 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.
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Closed-loop system response with PID control. From top to bottom: (a) Position tracking showing desired (dashed red) and actual (solid blue) trajectories, demonstrating effective reference tracking, (b) Control force profile exhibiting initial peak of 259.2 N followed by smooth regulation, (c) Position error evolution showing maximum deviation of 0.107 m and final error of 0.004 m, (d) Velocity response with peak of 1.04 m/s, and (e) Acceleration profile with maximum of 2.14 m/s².</figcaption></figure>

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 char-
acteristics, as evidenced by the increased damping ratio of ζ = 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.

<figure><img src="image_7.png" /><figcaption>Fig. 6 Phase portrait of the closed-loop system illustrating the controlled state-space trajectory. The smooth path from initial conditions to the desired setpoint (1.0 m, 0 m/s) demonstrates the controller's ability to regulate both position and velocity simultaneously. The absence of multiple spirals indicates improved damping compared to the open-loop response.</figcaption></figure>

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. 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 trajec-
tory 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 track-
ing 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 aggres-
sive yet controlled motions when necessary while maintaining system stability. These
performance metrics validate the theoretical design of the PID controller with gains
Kp = 200 N/m, Ki = 50 N/(ms), and Kd = 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 itera-
tions 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.

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Violin plots of execution time distributions for Python and R under controlled and open-loop implementations. The left panel represents controlled implementations, while the right panel shows open-loop implementations.</figcaption></figure>
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.

<figure><img src="image_9.png" /><figcaption>**Fig. 8** Violin plots of memory usage distributions for Python and R under controlled and open-loop implementations. The left panel represents controlled implementations, while the right panel shows open-loop implementations.</figcaption></figure>

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] emphasize 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
452

453

454

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.

455

456

457

**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.

458

459

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

460

## Declarations

461

**Conflict of interest.** The authors declare there is no conflict.

462

463

**Competing interests.** Authors do not have any competing financial interest to declare.

464

## References

465

466

467

468

[1] Irawan, D.E., Pourret, O., Besançon, L., Herho, S.H.S., Ridlo, I.A., Abraham, J.: Post-Publication Review: The Role of Science News Outlets and Social Media. Annals of Library and Information Studies **71**, 465–474 (2024) https://doi.org/10.56042/alis.v71i4.1425

469

470

471

472

[2] Fraser, N., Brierley, L., Dey, G., Polka, J.K., Pálfy, M., Nanni, F., Coates, J.A.: The evolving role of preprints in the dissemination of COVID-19 research and their impact on the science communication landscape. PLoS Biology **19**(4), 3000959 (2021) https://doi.org/10.1371/journal.pbio.3000959

473

474

475

476

[3] Sugimoto, C.R., Work, S., Larivière, V., Haustein, S.: Scholarly use of social media and altmetrics: A review of the literature. Journal of the Association for Information Science and Technology **68**(9), 2037–2062 (2017) https://doi.org/10.1002/asi.23833

477

478

479

480

481

[4] Harris, C., Millman, K., Walt, S., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N., Kern, R., Picus, M., Hoyer, S., Kerkwijk, M., Brett, M., Haldane, A., Río, J., Wiebe, M., Peterson, P., Oliphant, T.: Array Programming with NumPy. Nature **585**(7825), 357–362 (2020) https://doi.org/10.1038/s41586-020-2649-2

482

483

484

485

486

[5] Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Courneau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., Walt, S.J., Brett, M., Wilson, J., Millman, K.J., Mayorov, N., Nelson, A.R.J., Jones, E., Kern, R., Larson, E., Carey, C.J., Polat, , Feng, Y., Moore, E.W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E.A., Harris, C.R.,
487

488

489

Archibald, A.M., Ribeiro, A.H., Pedregosa, F., Mulbregt, P.: SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods **17**(3), 261–272 (2020) https://doi.org/10.1038/s41592-019-0686-2

490

491

492

[6] McKinney, W.: Data Structures for Statistical Computing in Python. In: Walt, S., Millman, J. (eds.) Proceedings of the 9th Python in Science Conference, pp. 51–56 (2010). https://doi.org/10.25080/majora-92bf1922-00a

493

494

495

[7] Soetaert, K., Petzoldt, T., Setzer, R.W.: Solving Differential Equations in R: Package deSolve. Journal of Statistical Software **33**(9), 1–25 (2010) https://doi.org/10.18637/jss.v033.i09

496

497

498

499

500

501

[8] Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L.D., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, T.L., Miller, E., Bache, S.M., Müller, K., Ooms, J., Robinson, D., Seidel, D.P., Spinu, V., Takahashi, K., Vaughan, D., Wilke, C., Woo, K., Yutani, H.: Welcome to the Tidyverse. Journal of Open Source Software **4**(43), 1686 (2019) https://doi.org/10.21105/joss.01686

502

503

[9] Ogata, K.: Modern Control Engineering. Prentice Hall, Upper Saddle River, NJ, USA (2010)

504

505

[10] Khalil, H.K.: Nonlinear Systems. Prentice Hall, Upper Saddle River, New Jersey, USA (2002)

506

507

[11] Dorf, R.C., Bishop, R.H.: Modern Control Systems. Prentice Hall, Upper Saddle River, NJ, USA (2011)

508

509

[12] Nise, N.S.: Control Systems Engineering. John Wiley & Sons, New York City, NY, USA (2020)

510

511

[13] Franklin, G.F., Powell, J.D., Emami-Naeini, A.: Feedback Control of Dynamic Systems. Pearson, Upper Saddle River, NJ, USA (2015)

512

513

[14] Åström, K.J., Murray, R.M.: Feedback Systems: An Introduction for Scientists and Engineers. Princeton University Press, Princeton, NJ, USA (2008)

514

515

[15] Chen, C.-T.: Linear System Theory and Design. Oxford University Press, Oxford, UK (1995)

516

517

[16] Waskom, M.L.: Seaborn: Statistical Data Visualization. Journal of Open Source Software **6**(60), 3021 (2021) https://doi.org/10.21105/joss.03021

518

519

520

[17] Georges, A., Buytaert, D., Eeckhout, L.: Statistically Rigorous Java Performance Evaluation. ACM SIGPLAN Notices **42**(10), 57–76 (2007) https://doi.org/10.1145/1297105.1297033

521

[18] Thompson, S.K.: Sampling. Wiley Series in Probability and Statistics, New York
522 City, NY, USA (2012). https://doi.org/10.1002/9781118162934

523 [19] Herho, S., Anwar, I., Herho, K., Dharma, C., Irawan, D.: COMPARING
524 SCIENTIFIC COMPUTING ENVIRONMENTS FOR SIMULATING 2D NON-
525 BUOYANT FLUID PARCEL TRAJECTORY UNDER INERTIAL OSCIL-
526 LATION: A PRELIMINARY EDUCATIONAL STUDY. Indonesian Physical
527 Review 7(3), 451–468 (2024) https://doi.org/10.29303/ipr.v7i3.335

528 [20] Herho, S., Fajary, F., Herho, K., Anwar, I., Suwarman, R., Irawan, D.: Reap-
529 praising Double Pendulum Dynamics across Multiple Computational Platforms.
530 Preprints (2024) https://doi.org/10.20944/preprints202405.0232.v1

531 [21] Herho, S., Kaban, S.N., Irawan, D.E., Kapid, R.: Efficient 1D Heat Equation
532 Solver: Leveraging Numba in Python. Eksakta : Berkala Ilmiah Bidang MIPA
533 (2024) https://doi.org/10.24036/eksakta/vol25-iss02/487

534 [22] Jain, R.: The Art of Computer Systems Performance Analysis: Techniques for
535 Experimental Design, Measurement, Simulation, and Modeling. John Wiley &
536 Sons, New York, USA (1991)

537 [23] Brown, R.J.C., Brown, R.F.C.: Statistical Analysis of Measurement Data. Royal
538 Society of Chemistry, London, United Kingdom (1998)

539 [24] Hoaglin, D.C., Mosteller, F., Tukey, J.W.: Understanding Robust and Exploratory
540 Data Analysis. Wiley-Interscience, Hoboken, NJ, USA (2000)

541 [25] Mostafavi, S., Hakami, V., Paydar, F.: Performance Evaluation of Software-
542 Defined Networking Controllers: A Comparative Study. Computer and Knowl-
543 edge Engineering 2(2), 63–73 (2020) https://doi.org/10.22067/cke.v2i2.84917

544 [26] Montgomery, D.C.: Design and Analysis of Experiments. John Wiley & Sons,
545 New York City, NY, USA (2017)

546 [27] Bulmer, M.G.: Principles of Statistics. Dover Publications, New York, USA
547 (1979)

548 [28] Razali, N.M., Wah, Y.B.: Power Comparisons of Shapiro-Wilk, Kolmogorov-
549 Smirnov, Lilliefors and Anderson-Darling Tests. Journal of Statistical Modeling
550 and Analytics 2(1), 21–33 (2011)

551 [29] Siegel, S.: Nonparametric Statistics for the Behavioral Sciences. McGraw-Hill,
552 New York City, NY, USA (1956)

553 [30] Wilcoxon, F.: Individual Comparisons by Ranking Methods. Biometrics Bulletin
554 1(6), 80–83 (1945) https://doi.org/10.2307/3001968

555 [31] Levene, H.: Robust Tests for Equality of Variances. In: Olkin, I. (ed.) Contributions to Probability and Statistics, pp. 278–292. Stanford University Press,
557

Stanford, California, USA (1960)

558

[32] Brown, M.B., Forsythe, A.B.: Robust Tests for the Equality of Variances. Journal of the American Statistical Association 69(346), 364–367 (1974) https://doi.org/10.1080/01621459.1974.10482955

559

560

561

[33] Sullivan, G.M., Feinn, R.: Using Effect Size—or Why the P Value is Not Enough. Journal of Graduate Medical Education 4(3), 279–282 (2012) https://doi.org/10.4300/JGME-D-12-00156.1

562

563

564

[34] Cohen, J.: Statistical Power Analysis for the Behavioral Sciences. Routledge, New York City, NY, USA (1988). https://doi.org/10.4324/9780203771587

565

566

[35] Sawilowsky, S.S.: New Effect Size Rules of Thumb. Journal of Modern Applied Statistical Methods 8(2), 597–599 (2009) https://doi.org/10.22237/jmasm/1257035100

567

568

569

[36] Iglewicz, B., Hoaglin, D.C.: Volume 16: How to Detect and Handle Outliers. ASQC Quality Press, Milwaukee, WI, USA (1993)

570

571

[37] Chandola, V., Banerjee, A., Kumar, V.: Anomaly Detection: A Survey. ACM Computing Surveys 41(3), 1–58 (2009) https://doi.org/10.1145/1541880.1541882

572

573

574

[38] Rousseeuw, P.J., Croux, C.: Alternatives to the Median Absolute Deviation. Journal of the American Statistical Association 88(424), 1273–1283 (1993) https://doi.org/10.1080/01621459.1993.10476408

575

576

577

[39] Tukey, J.W.: Exploratory Data Analysis. Addison-Wesley, Reading, MA, USA (1977)

578

579

[40] Aggarwal, C.C.: Outlier Analysis. Springer, New York City, NY, USA (2013). https://doi.org/10.1007/978-1-4614-6396-2

580

581

[41] Mytkowicz, T., Diwan, A., Hauswirth, M., Sweeney, P.F.: Producing Wrong Data Without Doing Anything Obviously Wrong! ACM SIGARCH Computer Architecture News 37(1), 265–276 (2009) https://doi.org/10.1145/2528521.1508275

582

583

584

[42] Altman, D.G., Machin, D., Bryant, T.N., Gardner, M.J.: Statistics with Confidence: Confidence Intervals and Statistical Guidelines. BMJ Books, London, UK (2013)

585

586

587

[43] Masini, S., Bientinesi, P.: High-Performance Parallel Computations Using Python as High-Level Language. In: Guarracino, M.R., Vivien, F., Träff, J.L., Cannatoro, M., Danelutto, M., Hast, A., Perla, F., Knüpfer, A., Di Martino, B., Alexander, M. (eds.) Euro-Par 2010 Parallel Processing Workshops, pp. 541–548. Springer, Berlin, Heidelberg (2011). https://doi.org/10.1007/978-3-642-21878-1_66

588

589

590

591
[44] Raschka, S., Patterson, J., Nolet, C.: Machine Learning in Python: Main Developments and Technology Trends in Data Science, Machine Learning, and Artificial Intelligence. Information **11**(4), 193 (2020) https://doi.org/10.3390/info11040193

[45] Watson, A., Babu, D.S.V., Ray, S.: Sanzu: A data science benchmark. In: 2017 IEEE International Conference on Big Data (Big Data), pp. 263–272 (2017). https://doi.org/10.1109/BigData.2017.8257934