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

HAL
open science

# SINDy-Fish: Sparse, interpretable, and parsimonious reduced-order modeling of hydrodynamic forces in undulating swimmers

Ahmad Saeed, Hamid Saeed

► To cite this version:

Ahmad Saeed, Hamid Saeed. SINDy-Fish: Sparse, interpretable, and parsimonious reduced-order modeling of hydrodynamic forces in undulating swimmers. Physics of Fluids, In press, 37 (12), pp.1904. 10.1063/5.0305318. hal-05475410

HAL Id: hal-05475410

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

Submitted on 24 Jan 2026

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# SINDy-Fish: Sparse, interpretable, and parsimonious reduced-order modelling of hydrodynamic forces in undulating swimmers

Ahmad Saeed,<fnref n="1" />, a) and Hamid Saeed<fnref n="2" />, b)

(Dated: 2 October 2025)

Modeling the complex fluid-structure interactions of fish-like propulsion is a significant challenge, with high-fidelity simulations being computationally prohibitive for design optimization or real-time control. To address this, we introduce SINDy-Fish, a sparse, interpretable reduced-order modeling framework for predicting hydrodynamic forces on undulating swimmers. Leveraging sparse regression, specifically, the mixed-integer optimization sparse regression on high-fidelity computational fluid dynamics data, SINDy-Fish discovers the underlying governing equations from data. We apply this framework to an anguilliform (eel-like) swimmer at a Reynolds number of 5000 and across biologically relevant, thrust-producing Strouhal numbers ($St = 0.3–0.5$). Model parsimony and accuracy are systematically balanced using the Akaike Information Criterion. The results show exceptional predictive accuracy, with $R^2$ scores exceeding 0.94 for lift and 0.995 for drag on unseen test data. The identified model for lift is a self-excited, self-limiting nonlinear oscillator with cubic nonlinearity, which successfully recovers a stable limit cycle, while drag is modeled by a concise algebraic expression of the lift dynamics. By generating explicit, computationally efficient equations, SINDy-Fish provides a powerful and scalable tool for the design, optimization, and real-time control of energy-efficient bioinspired underwater robots.

## I. INTRODUCTION

Time-varying flow-induced forces, such as lift and drag, play a central role in the performance and control of underwater robotic systems. The hydrodynamics of natural swimming have been extensively investigated, with researchers examining how aquatic organisms exploit different body and/or caudal fin (BCF) propulsion strategies—namely anguilliform, subcarangiform, carangiform, and thunniform—to navigate and propel themselves efficiently<fnref n="1" />. Among these, carangiform and anguilliform motions have received considerable attention due to their prevalence in nature, with prior studies offering in-depth evaluations of their associated hydrodynamic performance<fnref n="2" />⁻⁶. These unsteady hydrodynamic forces are the result of complex vortex-body interactions and directly impact propulsion efficiency, maneuverability, and energy harvesting capabilities⁷. Accurate and efficient prediction of such forces is especially critical for real-time feedback control in dynamic flow environments, where rapid changes in force response can significantly affect system behavior. The study of the hydrodynamics of aquatic animal locomotion has been revolutionized in recent years by the advent of particle image velocimetry (PIV). This technique enables detailed visualization of water flow both in the wake and around the bodies of swimming fish, offering unprecedented insight into the vortical structures and force generation mechanisms that underlie efficient propulsion⁸⁻¹⁰.

Understanding the basic principles of aquatic propulsion requires more than descriptive studies of live animals. To probe cause–effect relationships we must be able to manipulate morphological and kinematic parameters, alter material properties

such as stiffness, and impose motions that do not occur in nature. Robotic and mechanical models are particularly valuable in this context: they afford precise, repeatable control over geometry, actuation, and timing, and thus enable rigorous tests of mechanistic hypotheses regarding thrust production, efficiency, and maneuverability. By complementing biological observation with engineered models, researchers can assess whether observed fish behaviors correspond to performance optima or reflect ecological and developmental constraints<fnref n="2" />,¹¹⁻¹³.

However, replicating the efficiency, agility, and adaptability of swimming fish remains a central challenge in the development of advanced underwater robots. Colgate and Lynch ¹⁴ present a foundational review of swimming mechanics and control, emphasizing how fish physiology and motion strategies can inform biomimetic designs capable of precise maneuvering, station keeping, and energy-efficient propulsion. Their analysis frames the core challenges—hydrodynamic modeling, stability, maneuverability, actuation, sensing, and trajectory planning—providing a critical bridge between biological performance and engineering implementation.

Accurately modeling the hydrodynamic forces generated by undulating swimmers is essential for both understanding biological locomotion and designing efficient bioinspired underwater vehicles. The fluid-structure interactions involved are inherently complex and nonlinear, making prediction and interpretation challenging. Although high-fidelity computational fluid dynamics (CFD) simulations can resolve these details, they are computationally intensive, making them impractical for real-time control applications, rapid design iterations, or extensive parametric studies, especially when exploring a wide range of kinematic behaviors or environmental conditions ¹⁵,¹⁶.

Reduced-order models (ROMs) offer an attractive alternative to costly CFD simulations. A variety of approaches exist for constructing ROMs, ranging from physics-based formulations

<footnote n="1">¹) *National University of Sciences and Technology, Islamabad, Pakistan*</footnote>

<footnote n="2">²) *University of Engineering and Technology, Lahore, Punjab, Pakistan*</footnote>

<footnote>a) corresponding author: ahmadsaeed1994@gmail.com</footnote>

<footnote>b) 2019mstfe13@student.uet.edu.pk</footnote>
to purely data-driven techniques. Among the most widely used are proper orthogonal decomposition (POD)¹⁷ and dynamic mode decomposition (DMD)¹⁸, which extract the most energetic modes from fluid flow data to yield a low-rank representation of the dynamics. The Navier–Stokes equations can then be projected onto these modes, producing approximate low-dimensional models expressed in terms of mode coefficients. Such models have been successfully applied to canonical problems such as cylinder wakes¹⁹⁻²¹. However, their direct extension to flows around flexible, moving boundaries remains challenging²². While certain POD–Galerkin-based formulations have shown promise for rigid-body motions after suitable modifications²³,²⁴, further advances are needed to extend these approaches to bioinspired propulsion systems with complex fluid–structure interactions.

However, recent developments in computational resources and data-driven methods have enabled the use of non-intrusive modeling and machine learning techniques to construct efficient surrogate models for undulating swimmers²⁵. Reinforcement learning (RL) methods for mobile robotics, and particularly for swimming robots operating in complex, unmodeled fluid environments, require large datasets that span a substantial portion of the relevant state space. Generating such data is often both computationally expensive and time-consuming, whether through high-fidelity simulations or controlled experiments. This challenge highlights the sample efficiency of many RL algorithms when applied to real-world robotic systems. Furthermore, the gap between simulation and reality complicates direct policy transfer, as simplified models cannot fully capture the rich nonlinear dynamics of aquatic locomotion. To mitigate these issues, recent research has explored the use of reduced-order modeling, physics-informed learning, and surrogate data generation techniques to accelerate training while preserving the essential hydrodynamic features necessary for robust control and adaptation in real environments²⁶,²⁷. Rodwell and Tallapragada (2023)²⁸ introduced a physics-informed reinforcement learning framework for motion control of fish-like swimming robots, leveraging curriculum learning and transfer learning with a deterministic policy gradient algorithm. Key results demonstrate accelerated training convergence, emergent "burst-and-coast" swimming strategies mimicking biological fish, and superior performance over direct RL on fluid simulations, advancing biomimetic robotics and efficient control of underactuated aquatic systems.

Brunton et al. (2016)²⁹ introduced the sparse identification of nonlinear dynamics (SINDy) framework, a powerful machine learning approach for discovering governing equations of dynamical systems directly from measurement data. By leveraging sparse regression techniques, the method identifies parsimonious models that balance accuracy and complexity, assuming that the underlying dynamics are sparse in a high-dimensional space of candidate nonlinear functions. Demonstrated on diverse systems, including the chaotic Lorenz equations, fluid vortex shedding behind a cylinder, and parameterized bifurcations like the Hopf normal form, SINDy effectively reconstructs low-dimensional representations from high-dimensional data, offering robustness to noise and ap-

plicability to various fields, including fluid dynamics.

Bertsimas and Gurnee (2023)³⁰ proposed the mixed-integer optimization sparse regression (MIOSR) method, which reformulates the SINDy problem as an exact mixed-integer optimization task to achieve provably optimal sparse regression solutions. MIOSR enhances model discovery accuracy, sample efficiency, and robustness to noise while enabling the seamless incorporation of physical constraints³¹. Evaluated on canonical ordinary and partial differential equations, the approach outperforms heuristic sparse regression methods³²,³³, providing a flexible and interpretable tool for data-driven dynamical system identification in scientific computing.

This study presents a sparse and parsimonious identification of nonlinear dynamics to overcome existing challenges in modeling the nonlinear hydrodynamic force of fish swimming, named SINDy-Fish. The SINDy-Fish framework uses the SINDy with the MIOSR approach specifically for the hydrodynamic forces in undulatory swimming. It seeks to develop sparse, nonlinear, and human-interpretable reduced-order models derived right from simulation data. Its core strength lies in generating explicit mathematical equations that describe the relationship between the swimmer's kinematics and the hydrodynamic forces, moving beyond black-box predictions to reveal governing physical principles. In contrast to prior SINDy studies on stationary bluff bodies³⁴,³⁵, the present work advances to actively undulating swimmers where body kinematics directly drive the flow. By modeling both lift and thrust/drag across a range of Strouhal numbers, we establish SINDy-Fish as the first reduced-order framework tailored to fish propulsion, providing a predictive foundation for control and design of bio-inspired underwater robots.

Specifically, in this study, we develop SINDy-Fish for an-guilliform (eel-like) swimmers. Recognizing the critical role of kinematics, particularly the Strouhal number (St), in determining the swimming regime, we construct separate, param-eterized SINDy-Fish models for different values of St (0.3 to 0.5) with an increase of 0.05. This approach explicitly targets the thrust-producing regimes at Re=5000. The models are rigorously trained and validated against high-fidelity CFD data. The resulting SINDy-Fish ROMs provide high accuracy in force prediction while offering unprecedented interpretability through their sparse, explicit equations. By distilling the complex fluid dynamics into computationally efficient mod-els, SINDy-Fish bridges the gap between high-cost CFD and practical applications. It enables a deeper understanding of aquatic biomechanics and provides a powerful tool for the pre-pictive design and real-time control of energy-efficient bioin-spired underwater vehicles.

The remainder of this manuscript is organized as follows. Section II provides an overview of the numerical simulation framework employed for modeling undulating swimmers. The proposed SINDy-Fish framework for predicting lift and drag coefficients is introduced in Section III. Results and discussion, covering both the numerical simulations and the SINDy-Fish approach, are presented in Section IV. Finally, Section IV B summarizes the key findings and outlines possible future directions for extending this work.
<figure><img src="image_3.png" /><figcaption>FIG. 1: Schematic overview of the proposed SINDy-Fish framework. A sparse dynamical model is identified for the lift coefficient, while an algebraic model is employed for the drag coefficient.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 2: Fish motion and amplitude envelope employed for numerical simulations.</figcaption></figure>
## II. NUMERICAL SIMULATIONS

### A. Geometry and Kinematics

The swimmer is modeled using a NACA0012 airfoil of thickness $D = 0.12L$, where $L$ denotes the chord length. The foil undergoes undulatory motion in a uniform free stream $U_\infty$, with the chord line representing the backbone of a swimmer at equilibrium. The lateral displacement of a point along the chord is denoted by $\Delta y(s,t)$, where $s \in [0, L]$ is the chordwise coordinate and $t$ is time. The nondimensional coordinate is defined as $X = s/L$, with $X = 0$ and $X = 1$ corresponding to the leading edge (LE) and trailing edge (TE), respectively.

The backbone motion is expressed as

$$ \Delta y(s,t) = \alpha(s) \cos [ks - \Omega t], \quad (1) $$

where $k = 2\pi/\lambda$ is the wavenumber and $\lambda$ is the wavelength, $\Omega = 2\pi f$ the angular frequency, and $\alpha(s)$ the amplitude envelope, which increases monotonically from LE to TE.

The nondimensional form of the lateral displacement reads

$$ \Delta Y(X, \tau) = A(X) \cos \left[ 2\pi \left( \frac{X}{\lambda^*} - \frac{St}{2A_{\max}} \tau \right) \right], \quad (2) $$

where $\Delta Y = \Delta y/L$, $A(X) = \alpha(s)/L$, $\tau = tU_\infty/L$, $St = fL/U_\infty$ is the Strouhal number, and $\lambda^* = \lambda/L$ is the nondimensional wavelength.

In this study, an anguilliform amplitude profile with maximum TE amplitude $A_{\max} = 0.1$ is adopted$^{36,37}$, given by

$$ A(X) = 0.0367 + 0.0323X + 0.0310X^2, \quad 0 < X < 1. \quad (3) $$

Figures in 2 illustrate the backbone kinematics and amplitude envelope, respectively.

### B. Flow configuration

In this work, we solve two-dimensional flow past an undulating swimmer using Navier-Stokes equations.

The incompressible Navier-Stokes equations are given as

$$ \nabla \cdot \mathbf{u} = 0, \quad (4) $$

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\nabla p + \frac{1}{Re} \nabla^2 \mathbf{u}, \quad (5) $$

where $\mathbf{u} = (u_1, u_2)^T$ is the velocity vector, $p$ is the pressure field, and $Re = U_\infty L/\nu$ denotes the Reynolds number. The first equation 4 enforces the incompressibility condition, while the second equation 5 represents the conservation of momentum, including the effects of pressure and viscous diffusion.

These equations are solved on a uniform Cartesian grid to simulate the flow around two-dimensional undulating swimmers for different Strouhal numbers. The CFD solver employed in this study is capable of handling dynamic boundaries and resolving unsteady wakes across a wide range of Reynolds numbers$^{38,39}$.

The simulation domain is prescribed with a uniform inflow condition at the inlet, a Neumann boundary condition at the outflow, and free-slip conditions along the upper and lower boundaries. A no-slip condition is enforced on the undulating body surface to capture the interaction between the moving swimmer and the surrounding fluid. The spatial resolution is expressed in terms of grid points per chord length $L$, with a uniform mesh spacing of $\Delta x = \Delta y = L/128$, based on the results of a convergence study. The details of the numerical simulation, including validation and the convergence study, are provided in our previous work$^{25}$. In the present study, we employed the same numerical parameters and simulation setup to investigate the flow past undulating swimmers, thereby ensuring consistency and reliability of the results.

Fishes are typically observed to swim within a narrow Strouhal number range of $0.2 < St < 0.5$, as reported by Triantafyllou et al. (2004)$^{36}$. This regime has been identified as optimal for efficient propulsion, balancing thrust generation with minimal energy expenditure. The consistency of this behavior across diverse species highlights the Strouhal number as a fundamental parameter governing efficient swimming dynamics. Within the SINDy-Fish framework, this provides a natural benchmark: the discovered reduced-order models can be evaluated against this biologically established regime to ensure that the extracted dynamics are not only mathematically parsimonious but also physically consistent with efficient swimming strategies.

The drag coefficient $C_D$ can be obtained by simply normalizing the total force applied on the x-axis ($F_x$) by the dynamic pressure of the free stream $U_\infty$. Similarly, the lift coefficient $C_L$ is obtained by normalizing the total force in the lateral direction ($F_y$) by the dynamic pressure of the free stream $U_\infty$. The expressions for the $C_L$ and the $C_D$ are given as:

$$ C_D = - \frac{F_x}{0.5 \rho U_\infty^2 L} \quad (6a) $$

$$ C_L = \frac{F_y}{0.5 \rho U_\infty^2 L} \quad (6b) $$

In addition, $C_T = -C_D$ where $C_T$ represents the thrust coefficient, and the time-averaged quantities are presented as $\tilde{C}_D$, $\tilde{C}_T$, $\tilde{C}_L$ to distinguish them from their instantaneous counterparts.

### C. Data collection

We performed simulations of the undulating swimmer in the Strouhal number range $St = 0.3-0.5$, at a Reynolds number of $Re = 5000$ and a non-dimensional wavelength of $\lambda^* = 1.5$. For each case, high-fidelity numerical simulations were conducted, and the time histories of the lift and drag coefficients were recorded for the construction of reduced-order models. After the flow reached a statistically steady state, lift coefficient ($C_L$) data spanning 50 undulation cycles, each discretized into 64 time steps, were collected to ensure both tem-
poral accuracy and statistical convergence. Similarly, drag co-efficient ($C_D$) data over 90 undulation cycles were gathered for sparse modeling. This carefully curated dataset provided the foundation for training, validation, and testing of the SINDy-Fish framework, thereby demonstrating its capability to extract governing equations of force dynamics with both efficiency and fidelity.

### III. SINDY-FISH

Here, we provide a general representation for the proposed method of reduced-order modeling for the lift and drag coefficients of the undulating swimmers. Deriving ROMs from empirical data is a core challenge in mathematical physics, with deep roots in fluid dynamics. Traditional approaches constrain dynamics using a priori knowledge (e.g., Galerkin projection) or heuristically select model structures before parameter optimization. Simultaneous identification of model structure and parameters is combinatorially challenging due to the vast space of candidate functions. The SINDy method addresses this by exploiting sparsity: it efficiently identifies parsimonious dynamical systems from data, avoiding exhaustive search over all possible functional forms.

The general representation of the ROM is given by:

$$ \frac{d}{dt} q = f(q) \qquad (7) $$

where $q = [q_1, q_2]^T \in \mathbb{R}^2$. SINDy enables the non-intrusive learning of a model $f$ without requiring the governing equations. We learn ROMs for hydrodynamic forces with the SINDy, which uses either heuristic algorithms such as convex $l_1$-regularized regression (LASSO)$^{32}$, sequentially thresholded least-squares (STLSQ)$^{29}$, sparse relaxed regularized regression (SR3)$^{33}$ or non-heuristic method such as MIOSR$^{30}$, to learn the model $f$ with the fewest terms possible from a library of candidate functions that might describe the dynamics.

Initially, data is assembled in a matrix as:

$$ Q = [q(t_1) \ q(t_2) \ \dots \ q(t_M)]^T \qquad (8) $$

Similarly, a matrix of derivatives is formed as:

$$ \dot{Q} = \left[ \frac{dq}{dt}(t_1) \quad \frac{dq}{dt}(t_2) \quad \dots \quad \frac{dq}{dt}(t_M) \right]^T . \qquad (9) $$

Here, $T$ denotes the transpose of the matrix. These derivatives are calculated using fifth-order spline derivatives, which are also suitable for noisy measurements$^{40}$. However, other alternatives are also available, the total-variation regularized derivatives$^{41}$. Based on the data $Q$, a library of candidate functions is formed as:

$$ \Theta(Q) = [1 \ Q \ P_2(Q) \ \dots \ P_d(Q) \ \dots \ \sin(Q) \ \dots]. \qquad (10) $$

Here, the matrix $P_d(Q)$ represents a matrix with column vectors for sparse regression, constructed by concatenating all $d$-th degree polynomial terms of the state vector $q(t) \in \mathbb{R}^n$ sampled at $M$ time points. Any basis functions may be used to generate this library, although polynomials work well for fluid dynamics.

The dynamical system in Eq. 7 may now be represented in terms of the data matrices in Eqs. 8 and 9 as

$$ \dot{Q} = \Theta(Q)\Xi. \qquad (11) $$

The coefficient matrix $\Xi$ contains column vectors $\xi_i$ that specify the active terms in each governing equation of the dynamical system 7. A parsimonious representation requires these coefficient vectors to yield accurate dynamical reconstructions in Eq. 11 while maintaining maximal sparsity in $\Xi$. This structure is typically identified through sparse regression techniques that simultaneously minimize prediction error and the number of non-zero parameters.

A sparse optimization problem is now defined to solve the sparse matrix of coefficients $\Xi$ that selects the terms from the library that are active in dynamics.

$$ \operatorname{argmin}_{\Xi} \|\dot{Q} - \Theta(Q)\Xi\|_2^2 + \lambda R(\Xi), \qquad (12) $$

where $R(\Xi)$ is the regularization term; $\lambda$ is a hyperparameter for the regularization.

The optimization problem in Eq. 12 is then solved by the MISOR algorithm using the open-source PySINDy code$^{40,42}$ to obtain the dynamics of $q(t)$.

SINDy can also be applied to find the algebraic model from the data. The general representation is given as:

$$ a = h(q) \qquad (13) $$

where $h$ can also be learned non-intrusively similar to $f$ using the candidate library functions in Eq. 10 and sparse optimization in Eq 12.

### A. Model selection

A key aspect of developing parsimonious models is selecting hyperparameters that effectively balance model fit, number of terms, and sparsity, since increasing degrees of freedom consistently reduces in-sample error. We achieve this equilibrium by selecting parameters that minimize the Akaike information criterion (AIC) metric. The AIC is widely used for model selection by balancing model accuracy and complexity$^{43-45}$. For a candidate model $i$, the AIC is defined as

$$ AIC_i = 2k - 2 \ln(E(\mathbf{q}, \gamma)), \qquad (14) $$

where $E(\mathbf{q}, \gamma) = P(\mathbf{q}|\gamma)$ is the likelihood of data $\mathbf{q}$ under the parameter estimate $\gamma$, and $k$ is the number of non-zero terms.
For finite datasets, a corrected form is applied:

$$AIC_c = AIC + \frac{2(k+1)(k+2)}{m-k-2}, \qquad (15)$$

The residual sum squared (RSS) is defined as :

$$RSS = \sum (\dot{Q} - \Theta(Q)\hat{\Xi})^2, \quad AIC = m \ln(RSS/m) + 2k. \quad (16)$$

m is the total number of measurements. Also, this RSS/ΣQ² normalization of RSS renders the error measure scale-invariant, ensuring that states with larger magnitudes do not dominate the criterion. Consequently, the AIC is based on a relative error metric, which provides a more balanced assessment across different state variables.

Relative support among competing models is measured by relative AIC (AIC_r), and it is given as

$$AIC_r = AIC_c - (AIC_c)_{\min}, \qquad (17)$$

The corrected AIC score for each model is rescaled by subtracting the minimum value, yielding the relative score (AIC_r)⁴⁵. AIC_r ≤ 2 indicates strong support, 4 ≤ AIC_r ≤ 7 indicates weak support, and AIC_r ≥ 10 suggests no support.

In this study, AIC_r is employed within the proposed SINDy-Fish framework to objectively rank candidate reduced-order models of hydrodynamic forces, ensuring a balance between predictive accuracy and model parsimony. The sampled trajectory is split into training, validation, and testing intervals, typically using the 80%, 10%, and 10% of the data, respectively. For each algorithm and hyperparameter configuration, models are trained on the training split, and the corrected Akaike information criterion is computed using the validation data. The final model is selected as the configuration that minimizes the AIC_c metric on the validation set.

## B. Model evaluation metrics

To evaluate the performance of the SINDy-Fish, two statistical metrics are employed, including the coefficient of determination/R-squared (R²) and the mean squared error (MSE). These metrics provide complementary perspectives on the accuracy and robustness of the proposed approach.

### 1. Coefficient of determination

The R² score measures the proportion of variance in the observed data that is explained by the model and is defined as

$$R^2 = 1 - \frac{\sum_{i=1}^{M} (q_i - \hat{q}_i)^2}{\sum_{i=1}^{M} (q_i - \bar{q})^2}, \qquad (18)$$

where q_i and ̂q_i denote the observed and predicted values, respectively, q̄ is the mean of the observed values, and M is the total number of samples. An R² value close to 1 indicates a strong agreement between the predictions and the reference data.

### 2. Mean squared error

The mean squared error quantifies the average squared difference between the observed and predicted values:

$$MSE = \frac{1}{M} \sum_{i=1}^{M} (q_i - \hat{q}_i)^2. \qquad (19)$$

Lower values of MSE indicate improved predictive accuracy, with an MSE of zero representing a perfect fit.

## IV. RESULTS AND DISCUSSION

We performed numerical simulations at Re = 5000 and St = 0.30 – 0.50, λ* = 1.5, and data was collected. The instantaneous vorticity snapshots are shown in Fig. 3. The temporal profiles of the lateral force coefficient C_L and drag coefficient C_D exhibit periodic oscillations, with C_D displaying higher-frequency harmonics compared to C_L. Figures 5, 6 show the Fourier spectral analysis. Spectral analysis identifies the vortex-shedding frequency (f_s) as the dominant mode in C_L, while 2f_s prevails in C_D for all St numbers, corroborating earlier findings for cylinders⁴⁶. In particular, odd harmonics 3f_s emerge in C_L, whereas even harmonics 4f_s dominate in C_D, reflecting the nonlinear coupling between body kinematics and vortex dynamics. These results underscore the universality of wake-body interactions across St and highlight the role of St number in governing thrust efficiency of undulating swimmers.

To further characterize the system dynamics, we construct phase maps using C_L and C_D as state variables, illustrated in Fig. 4 for St = 0.30 and St = 0.4. The analysis employs steady-state time histories that span thirty oscillation cycles, ensuring periodicity. The phase portrait reveals two distinct loops, demonstrating that C_D oscillates at twice the frequency of C_L. This frequency doubling aligns with spectral observations (where C_D peaks at 2f_s) while C_L dominates at f_s and reflects the nonlinear coupling between lateral and drag/thrust forces. Such behavior mirrors phase-space patterns observed in vortex-driven systems, reinforcing the role of harmonic interactions in undulatory propulsion.

### A. Lift Modeling

In this section, a low-dimensional model of the lift coefficient for undulating swimmers is presented. The primary objective is to capture the essential unsteady features of the lift dynamics using an ordinary differential equation (ODE) based model. Such a model not only provides physical insight into the mechanisms responsible for lift generation but also enables efficient predictions without the need for computationally expensive simulations. The following subsection describes the identified model and validation against high-fidelity CFD data.

It is well known that the two-dimensional cylinder flow behaves as a self-excited, self-limiting, and harmonic nonlinear
<figure><img src="image_5.png" /><figcaption>FIG. 3: Instantaneous vorticity ($\Gamma_z L/U_\infty$) in range (-5,5) at normalized time ($\tau/T = 0$) (a) St=0.3, (b) St=0.35, (c) St=0.4, (d) St=0.45, (e) St=0.5</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 4: Trajectory of the lift and drag coefficient in the plane $(C_L, C_D)$. (a) St=0.3, (b) St=0.4</figcaption></figure>

oscillator. This behaviour is also visible in the time evolution of the instantaneous lift coefficients of undulating swimmers for the considered $St = 0.3 - 0.5$ numbers and $Re = 5000$ depicted in Figures 5.

Two models often employed for analyzing the lift force on a cylinder in uniform flow are the Rayleigh and van der Pol oscillators. Both are nonlinear, self-excited, and self-limiting oscillators featuring a cubic nonlinearity$^{47-49}$. Some research has incorporated additional terms into these models, such as the Duffing $C_L^3$ term introduced by Griffin$^{50}$ and Skop and Griffin$^{51}$, while Landl$^{52}$ examined a modified van der Pol oscillator incorporating an extra nonlinear aerodynamic damping term of fifth order ($C_L^4 \dot{C}_L$). However, this model involved many constants to be determined.

Nayfeh et al.$^{46}$ provided a rationale for selecting the van der Pol oscillator over the Rayleigh model to represent the lift coefficient in a stationary cylinder due to the wake. They highlighted that the van der Pol oscillator generates a phase of the third harmonic at $3St$ relative to the fundamental at $St$,

which is 90°, as opposed to 270° in the Rayleigh oscillator scenario. Their numerical simulations, involving the solution of the Reynolds-averaged Navier-Stokes equations at various $Re$ values, confirmed a phase close to 90°. Marzouk et al.$^{53}$ incorporated a Duffing component into the van der Pol oscillator, justifying this addition by the model's enhanced ability to precisely capture the relative phase of the third harmonic of $C_L$, which would be constrained to 90° in the absence of the Duffing term. Furthermore, they offered formulas that link the model parameters with the characteristics of $C_L$, such as the frequency and amplitude of the fundamental harmonic.

It is important to note that only odd nonlinearities can be analyzed in the wake oscillator for $C_L$, as its spectral content is predominantly characterized by a central peak at the Strouhal frequency (i.e., $St = f_s$), accompanied by harmonics at odd integer multiples of $St$. This structural insight can be leveraged within the sparse identification of nonlinear dynamics for the undulating swimmers SINDy-Fish framework, where the governing ODEs are modeled using candidate library functions that implicitly account for the lift coefficient dynamics. By restricting the model to physics-aligned nonlinearities, the approach retains interpretability and parsimony while ensuring computational efficiency in reconstructing the nonlinear lift coefficient behavior.

The sparse identification of the non-linear dynamics algorithm helps determine the equations that govern the dynamics of lift. A repository of polynomial functions is utilized, comprising terms that are linear, cubic, and constant, with a maximum polynomial degree of three. This selection of functions, which can be readily expanded if necessary, is an intuitive option to identify nonlinear oscillators based on the above discussion. The library of functions is given as:

$$ \Theta(q_1, q_2) = [1 \quad q_1 \quad q_2 \quad q_1^3 \quad q_1^2 q_2 \quad q_1 q_2^2 \quad q_2^3]. \quad (20) $$

Also, $q_1$ is normalized as $q_1/max|q_1|$ before creating the library. Such normalization ensured $-1 \le q_1 \le 1$, which sim-
<figure><img src="image_7.png" /><figcaption>FIG. 5: Time histories of the lift coefficient $C_L$ with the Fourier spectra for different cases. (a, b) at St=0.3, (c, d) at St=0.35, (e, f) at St=0.4, (g, h) at St=0.45, (i, j) at St=0.5</figcaption></figure>

plifies the sparse optimization problem in system identification. Given our polynomial library, we employ the procedure
<figure><img src="image_8.png" /><figcaption>FIG. 6: Time histories ($C_D$) with the Fourier spectra. (a,b) at St=0.3, (c,d) at St=0.35, (e,f) at St=0.4, (g,h) at St=0.45, (i,j) at St=0.5</figcaption></figure>

proposed by Mangan et al. (2017)⁴⁵, which uses the AIC$_r$ to select the model structure that optimally balances predictive
<figure><img src="image_9.png" /><figcaption>FIG. 7: AIC<sub>r</sub> criteria for models found by SINDy for different Strouhal numbers; Red color indicates AIC<sub>r</sub> = 0
(a)St = 0.3, (b)St = 0.35, (c)St = 0.4, (d)St = 0.45, (e)St = 0.5</figcaption></figure>

accuracy with model complexity. In this framework, models
with lower AIC<sub>r</sub> scores are considered more parsimonious,
providing a compact yet accurate representation of the dy-
namics. This methodology enables the systematic ranking of
candidate models along the pareto front.

Figure 7 shows the distribution of the identified models
in the complexity–AIC plane. It is important to note that
the model with AIC<sub>r</sub> = 0 is not necessarily the globally op-
timal model, but rather the best among those considered in the

present study. For this case, the model that achieves the most
favorable trade-off between accuracy and complexity is rep-
resented by a cubic polynomial. We get the following general
model for all the Strouhal numbers.
<figure><img src="image_10.png" /><figcaption>FIG. 8: Comparison between the true and predicted time histories of lift coefficients ($q_1$); (a) $St = 0.3$, (b) $St = 0.35$, (c) $St = 0.4$, (d) $St = 0.45$, (e) $St = 0.5$</figcaption></figure>

$$ \begin{cases} \frac{dq_1}{dt} = 1.0q_2, \\ \frac{dq_2}{dt} = -\omega^2 q_1 + \mu q_2 - \alpha q_1^2 q_2 - \beta q_1^3. \end{cases} \qquad (21) $$

Surprisingly, the equations governing the dynamics of lift force are identical to the circular cylinder detailed in Ref.⁵³. Table I summarizes the parameters of the reduced-order models identified for different Strouhal numbers. The table reports the time interval $τ$ used for model identification, the shedding frequency $ω_s$ identified from the direct numerical simulations, the shedding frequency predicted from the low-order model $ω ≈ ω_s$, and the linear and nonlinear damping coefficients $μ$ and $α$, along with the nonlinear Duffing parameter $β$. Predictive performance is quantified using the coefficient of determination ($R^2$) and the MSE for both training and test datasets. In all cases, the identified models exhibit high accuracy, with $R^2$ scores greater than 0.94 and low MSE values. Increasing $St$ is accompanied by higher values of $μ$ and $α$, indicating stronger nonlinear effects in the lift dynamics and presence of a limit cycle. Moreover, the noted increases in $μ$ and $α$

with the Strouhal number align with results from previous research on undulating swimmers. For example, Khalid et al.⁵⁴ report that as $St$ increases, wakes become more unsteady, vortex interference and wake deflection intensify, and hydrodynamic efficiency decreases. Similarly, Borazjani & Sotiropoulos (2009)⁵⁵ find that the form drag component peaks and then diminishes, whereas the frictional drag increases continuously with $St$; at higher $St$ and $Re$, the total drag saturates toward a constant value. These phenomena suggest that increasing $St$ drives increased energy dissipation via stronger vortex shedding, which is captured in our SINDY-Fish model by increasing linear and nonlinear damping coefficients, i.e., $μ$ and $α$ with the Strouhal number.

Also, $ω$ is very close to $ω_s$, the shedding frequency, with the mean difference being only 1% across all the $St$ numbers. We get both parameters $μ$ and $α$ as positive, such that the linear damping term is negative while the nonlinear damping term is positive. Consequently, small oscillations are amplified, and large oscillations are suppressed, leading to the emergence of stable limit cycles⁴⁶.

For the SINDY-Fish model, the corresponding scaled amplitude, defined in physical coordinates is given as
<figure><img src="image_11.png" /><figcaption>FIG. 9: Trajectory of the true and the identified systems in the phase plane $(q_1, q_2)$. For CFD, $q_2$ is normalized with $\omega_s$ as $(q_2/\omega_s)$, and for SINDy, it is normalized by $\omega$ as $(q_2/\omega)$. (a) $St = 0.3$, (b) $St = 0.35$, (c) $St = 0.4$, (d) $St = 0.45$, (e) $St = 0.5$; Red: SINDy, Black: CFD</figcaption></figure>

$$ \|q\| = \sqrt{q_1^2(t) + \left(\frac{q_2(t)}{\omega}\right)^2}, \qquad (22) $$

This definition corresponds to the Euclidean norm of the state vector expressed in the scaled coordinates $(q_1, q_2/\omega)$. In this representation, $q_1$ serves as the displacement-like coordinate, while $q_2/\omega$ acts as a velocity-like coordinate normalized by the shedding frequency $\omega$. The resulting amplitude measure provides a convenient way to monitor the periodic modulation of the trajectory and enables a direct comparison of oscillations across different systems on a common scale within the SINDy-Fish framework.

This behavior is closely aligned with the low-dimensional model of the circular cylinder based on Galerkin projection of the Navier-Stokes equations onto the span of the first two POD modes and using the mean flow as the reference state described in Ref$^{56}$. The amplitude of the oscillations in the low-dimensional model of the circular cylinder is perfectly normalized to unity, corresponding to a unit-radius limit cycle in the reduced-order modal space. This exact normalization enforces constant amplitude, and the system reduces to the simple harmonic oscillator.

In case of an identified model of the undulating swimmers, the average amplitude of oscillations yields $\|q\| \approx 1.02$ across all Strouhal numbers considered. Table I provides the $\|q\|$ for each $St$ case. Also, it can be observed from the state space plots in Figure 9. The slight deviation from unity reflects the residual influence of model truncation and numerical fitting errors, yet the near-constant amplitude demonstrates that the learned dynamics successfully recover a stable limit

cycle. Also, normalization is not imposed by construction, but rather emerges from the learned dynamics. Thus, while the cylinder model achieves exact unit normalization by construction, the SINDy-Fish model achieves an approximate but highly consistent normalization, indicating robustness of the reduced-order description across different kinematic regimes based on Strouhal numbers.

Next, we examine the validity of the lift model obtained from SINDy in predicting the lift coefficient. Figures 8, 10 depict the time histories obtained from direct numerical simulation and predicted by the low-order SINDy model for all the $St$ numbers. The SINDy-Fish model for the lift coefficient matches the CFD simulation results in amplitude and phase with excellent accuracy for all the $St$ numbers.

## B. Drag modeling

The identified dynamical system effectively predicts the underlying high-dimensional dynamics using only the instantaneous lift coefficient. However, it does not directly provide the corresponding drag force, which is often critical in undulating swimmers and aerodynamic applications.

To incorporate drag into the reduced-order framework, one possible approach is to augment the dynamical system with a nonlinear algebraic measurement equation identified via sparse regression, similar to the methodology employed for the circular cylinder in Ref$^{34}$. Several studies have modeled the oscillatory drag coefficient as a harmonic function at twice the Strouhal frequency (2$St$), algebraically linked to the lift coefficient. For instance, Nayfeh et al.$^{46}$ proposed a quadratic
<figure><img src="image_12.png" /><figcaption>FIG. 10: Trajectories of true and identified systems at different initial conditions; (a, b, c) represent $St = 0.3$ and (d, e, f) represent $St = 0.5$.</figcaption></figure>

term ($C_L \dot{C}_L$), which constrains the drag-lift phase relationship to 270°. Qin⁵⁷ introduced an alternative formulation based on a quadratic $C_L^2$ term coupled with a linear contribution in $C_L$. Marzouk et al.⁵³ combined these formulations, employing two algebraic terms to allow arbitrary phase adjustment, with coefficients derived analytically from the relative drag phase.

An alternative approach is to construct a new system of coupled ODEs to describe the dynamics of undulating swimmers. For instance, Kim and Perkins⁵⁸ proposed a coupled lift-drag model for elastic cables, where oscillatory drag was represented as the product of a spatial function (dependent on the static longitudinal coordinates) and a van der Pol oscillator at twice the Strouhal frequency (2St). The lift dynamics were modeled using a similar oscillator at St. Their formulation incorporated seven quadratic coupling terms with seven parameters grouped into four manually tuned functions, which were calibrated via ad hoc fitting to experimental data. This heuristic construction, lacking a systematic parameter identification methodology, introduced an excessive number of tunable parameters and ultimately failed to reproduce the correct mean

drag coefficient.

Given its simplicity and accuracy, the former approach is preferred. We model an algebraic expression for the drag coefficient of undulating swimmers in this work. Once the lift coefficient and its time derivative are identified, nonlinear relationships governing the drag coefficient can be determined via SINDy. Using a quadratic library with a bias term constructed from $q_1$ and $q_2$, least-squares minimization yields the following general expression for the drag coefficient a for all Strouhal numbers.

$$a = \eta_0 + \eta_1 q_1^2 + \eta_2 q_1 q_2 + \eta_3 q_2^2 \quad (23)$$

The drag model consists of four contributions: the $\eta_0$ term, which is very close to the mean of the drag coefficient of CFD data, followed by the lift-dependent terms. Among these terms, $\eta_0$ and the first two lift-dependent contributions dominate the dynamics, whereas the fourth term with coefficient $\eta_3$ provides only a minor effect.
<figure><img src="image_13.png" /><figcaption>FIG. 11: True and predicted time histories of drag coefficients; (a) $St = 0.3$, (b) $St = 0.35$, (c) $St = 0.4$, (d) $St = 0.45$, (e) $St = 0.5$</figcaption></figure>

Table II summarizes the identified drag model parameters for different Strouhal numbers. The table reports the time interval $\tau$ used for training, validation, and testing, the predicted $\eta_0 \approx \bar{C}_D$, the polynomial coefficients $\eta_1, \eta_2$, and $\eta_3$ obtained through sparse regression. Model performance is quantified using the coefficient of determination and the mean squared error for both training and test datasets. All cases show consistently high predictive accuracy, with $R^2 \ge 0.995$ and very low MSE values. Also, as $St$ increases, $\eta_0$ becomes more negative and the higher-order coefficients grow in magnitude, highlighting stronger nonlinear effects in the drag dynamics.

Furthermore, the first term $\eta_0$ is not a direct function of the lift dynamics and is very close to the mean drag. This information can be leveraged to incorporate additional kinematic parameters into the drag model by utilizing scaling laws, which relate the mean thrust or drag to fundamental kinematic parameters. This can be a promising direction for future research to investigate how such scaling laws can be directly incorporated into the SINDy-Fish framework, thereby enhancing the predictive capability of the identified models for different kinematic regimes beyond the Strouhal number. For instance, Thekkethil et al.⁵⁹ demonstrated that the mean thrust can be modeled as a function of the squared Strouhal number, $St^2$, and the dimensionless wavelength, $\lambda^*$, at $Re = 5000$. Similar scaling relations for mean thrust and drag have also been reported in other studies (see, e.g., Refs.⁶⁰⁻⁶²), providing a physically grounded means of prescribing the mean drag/thrust within the identified dynamical framework.

Figures 11 present the predicted drag coefficients from the identified SINDy models compared with the true drag coefficients obtained from direct numerical simulations in the testing range for all Strouhal numbers considered in this study. It is important to emphasize that the testing data were not used

during the training phase. The comparison demonstrates that the SINDy-based predictions closely match the DNS results, confirming the model's ability to generalize beyond the training data and accurately capture the underlying drag dynamics.

## CONCLUSION

In this study, we introduced SINDy-Fish, a novel sparse and interpretable reduced-order modeling framework designed to predict the hydrodynamic lift and drag forces on undulating swimmers. Addressing the substantial challenges posed by costly experimental data and computationally intensive high-fidelity simulations such as CFD, SINDy-Fish leverages the sparse identification of nonlinear dynamics framework and sparse regression techniques to uncover governing equations directly from CFD datasets. This approach aims to provide real-time force prediction at a mere fraction of the computational expense, effectively bridging the gap between high-cost simulations and practical applications.

The core strength of SINDy-Fish lies in its ability to generate explicit mathematical equations that delineate the intricate relationship between a swimmer's kinematics and the resulting hydrodynamic forces, thereby transcending "black-box" predictions to unveil underlying physical principles. Specifically developed for anguilliform (eel-like) swimmers operating at biologically relevant and thrust-producing Strouhal numbers (0.3–0.5) and a Reynolds number of 5000, the framework constructs separate, parameterized models for various Strouhal number values. The use of advanced algorithms, such as Mixed-Integer optimization and sparse identification of nonlinear dynamics, ensures provably optimal sparse re-
<table><caption>TABLE I: Lift parameters and model evaluation metrics for different *St* cases.</caption><thead><tr><th>St</th><th>time interval</th><th>ω<sub>s</sub></th><th>ω</th><th>μ</th><th>α</th><th>β</th><th>||q||</th><th>R<sup>2</sup>(train/test)</th><th>MSE (train/test)</th></tr></thead><tbody><tr><td>0.3</td><td>τ ∈ [100.6, 134]</td><td>9.359</td><td>9.212</td><td>0.143</td><td>0.546</td><td>5.096</td><td>1.01</td><td>0.943/0.935</td><td>2.8 × 10<sup>-2</sup> / 3.3 × 10<sup>-2</sup></td></tr><tr><td>0.35</td><td>τ ∈ [86.2, 114.8]</td><td>10.9188</td><td>10.817</td><td>0.182</td><td>0.727</td><td>5.165</td><td>1.01</td><td>0.998/0.998</td><td>6.4 × 10<sup>-4</sup> / 7.0 × 10<sup>-4</sup></td></tr><tr><td>0.4</td><td>τ ∈ [75.5, 100.5]</td><td>12.478</td><td>12.351</td><td>0.222</td><td>0.858</td><td>6.837</td><td>1.01</td><td>0.977/0.962</td><td>1.1 × 10<sup>-2</sup> / 1.9 × 10<sup>-2</sup></td></tr><tr><td>0.45</td><td>τ ∈ [67.1, 89.3]</td><td>14.038</td><td>13.918</td><td>0.251</td><td>0.969</td><td>7.762</td><td>1.01</td><td>0.968/0.943</td><td>1.6 × 10<sup>-2</sup> / 2.8 × 10<sup>-2</sup></td></tr><tr><td>0.5</td><td>τ ∈ [60.4, 80.4]</td><td>15.598</td><td>15.495</td><td>0.320</td><td>1.247</td><td>8.576</td><td>1.06</td><td>0.974/0.965</td><td>1.2 × 10<sup>-2</sup> / 1.7 × 10<sup>-2</sup></td></tr></tbody></table>

<table><caption>TABLE II: Drag parameters and model evaluation metrics for different *St* cases.</caption><thead><tr><th>St</th><th>time interval</th><th>C&#773;<sub>D</sub></th><th>η<sub>0</sub></th><th>η<sub>1</sub></th><th>η<sub>2</sub></th><th>η<sub>3</sub></th><th>R<sup>2</sup> (train/test)</th><th>MSE (train/test)</th></tr></thead><tbody><tr><td>0.3</td><td>τ ∈ [74.005, 134]</td><td>-0.108</td><td>-0.104</td><td>-0.054</td><td>0.018</td><td>4.99 × 10<sup>-4</sup></td><td>0.996/0.995</td><td>1.8 × 10<sup>-5</sup> / 2.2 × 10<sup>-5</sup></td></tr><tr><td>0.35</td><td>τ ∈ [63.4, 114.8]</td><td>-0.182</td><td>-0.161</td><td>-0.087</td><td>0.023</td><td>3.77 × 10<sup>-4</sup></td><td>0.996/0.996</td><td>3.7 × 10<sup>-5</sup> / 3.8 × 10<sup>-5</sup></td></tr><tr><td>0.4</td><td>τ ∈ [55.5, 100.5]</td><td>-0.267</td><td>-0.253</td><td>-0.094</td><td>0.027</td><td>4.20 × 10<sup>-4</sup></td><td>0.995/0.995</td><td>8.9 × 10<sup>-5</sup> / 8.7 × 10<sup>-5</sup></td></tr><tr><td>0.45</td><td>τ ∈ [49.3, 89.3]</td><td>-0.367</td><td>-0.368</td><td>-0.099</td><td>0.032</td><td>4.99 × 10<sup>-4</sup></td><td>0.995/0.995</td><td>1.4 × 10<sup>-4</sup> / 1.5 × 10<sup>-4</sup></td></tr><tr><td>0.5</td><td>τ ∈ [44.4, 84.4]</td><td>-0.486</td><td>-0.480</td><td>-0.127</td><td>0.038</td><td>4.69 × 10<sup>-4</sup></td><td>0.995/0.995</td><td>2.3 × 10<sup>-4</sup> / 2.5 × 10<sup>-4</sup></td></tr></tbody></table>

Regression solutions, enhancing accuracy, sample efficiency, and robustness against noise. Model selection is rigorously performed using the Akaike Information Criterion (AIC$_c$/AIC$_r$), balancing model fit, complexity, and sparsity.

Our results conclusively demonstrate that SINDy-Fish reliably captures force dynamics across a spectrum of kinematic conditions. The identified ROMs exhibit high predictive accuracy, with $R^2$ scores exceeding 0.94 for lift and 0.995 for drag, coupled with very low MSE values. For lift, the discovered governing equations strikingly resemble those of a circular cylinder, behaving as a self-excited, self-limiting nonlinear oscillator with cubic nonlinearity. These models successfully recover a stable limit cycle with an average oscillation amplitude very close to unity (≈ 1.02), indicating robust and consistent dynamics, even if not imposed by construction as in some other models. The SINDy lift models exhibit excellent agreement with CFD simulations in both amplitude and phase across all Strouhal numbers. For drag, an algebraic expression was identified, showcasing excellent accuracy to CFD results and robust generalization beyond the training data.

Despite its robust performance, a noteworthy aspect of the drag model is that its $\eta_0$ component is not directly a function of lift dynamics and is $\eta_0 \approx C_D$. This aspect should be extended to incorporate more kinematic parameters, such as $\lambda^*$, in this drag model. For that purpose, $\eta_0$ must be specified a priori, even though it is learned. Future research can address this by investigating the direct incorporation of established scaling laws for mean thrust or drag, which relate these forces to fundamental kinematic parameters like the squared Strouhal number and dimensionless wavelength $\lambda^*$ directly into the SINDy-Fish framework. This integration would further enhance the predictive capability of the identified models across various kinematic regimes beyond the Strouhal number investigated in this study. Furthermore, a limitation of the present study is that the SINDy-Fish framework has been developed

and validated entirely on two-dimensional simulation data. In realistic swimming and propulsion scenarios, three-dimensional flow structures such as tip vortices, spanwise transport, and oblique vortex shedding may play a significant role in shaping the hydrodynamic forces. These effects can influence both the amplitude and the spectral content of lift and drag signals, potentially altering the reduced-order dynamics. Future extensions of the present framework could incorporate three-dimensional datasets, either from high-fidelity simulations or from 3D particle image velocimetry. Such an extension would enable the discovery of low-dimensional models that account for spanwise coupling and tip effects, thereby broadening the applicability of the SINDy-Fish approach to fully 3D swimmers and robotic systems.

The interpretability, computational efficiency, and high accuracy of SINDy-Fish make it a scalable and powerful tool for advancing the design, optimization, and real-time control of energy-efficient bioinspired underwater robots, while also contributing to a deeper understanding of fundamental aquatic biomechanics.

## AUTHOR CONTRIBUTIONS

**Ahmad Saeed:** Conceptualization (equal); Data cura-
tion (equal); Formal analysis (equal); Investigation (equal);
Methodology (equal); Resources (equal); Software (equal);
Validation (equal); Visualization (equal); Writing—original
draft (equal); Writing—review and editing (equal). **Hamid
Saeed:** Investigation (equal); Validation (equal); Writ-
ing—review and editing (equal).
## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## CONFLICT OF INTEREST

The authors have no conflicts of interest to disclose.

## Appendix A: Polynomial models for SINDy-Fish coefficients as a function of Strouhal number

In order to construct reduced-order models of the lift and drag coefficients, the dependence of the parameters $\alpha, \beta, \mu$, and $\omega$ for lift model and $\eta_0, \eta_1, \eta_2$ and $\eta_3$ for drag model on the Strouhal number $St$ was analyzed over the range $0.3 \le St \le 0.5$. Polynomial regression was employed to provide compact algebraic representations of these trends.

Tables III, IV, summarize the polynomial models developed for the parameters $\alpha, \beta, \mu$, and $\omega$ and Tables V, VI, summaries the polynomial models for the parameters $\eta_0, \eta_1$, and $\eta_2$ as functions of the Strouhal number $St$ in the range $0.30 \le St \le 0.50$. Both cubic (third-order) and quartic (fourth-order) polynomial fits were considered. The cubic models provide compact functional forms with high accuracy, making them suitable for reduced-order modeling where model simplicity and generalization are prioritized. The quartic models, on the other hand, interpolate the five data points exactly ($R^2 = 1.0$), thereby recovering the measured values without error in the specified interval. However, quartic fits introduce more complex coefficients and may exhibit oscillatory behavior outside the calibration range. For example, the plots in Figures 12 for parameters of the lift model illustrate this tradeoff: cubic fits yield smoother approximations while quartic fits reproduce the data points exactly. For practical purposes, cubic models are recommended for predictive reduced-order modeling, while quartic models may be used when exact interpolation of the available dataset is required.

## Appendix B: Verification of SINDy on the van der Pol Oscillator

Before applying the SINDy framework to hydrodynamic data in the SINDy-Fish methodology, we first validate the approach on a well-known nonlinear dynamical system: the Van der Pol oscillator. This serves as a benchmark to demonstrate that the chosen sparse regression method (MIOSR optimizer) can recover the governing equations of a nonlinear oscillator from time-series data.

The van der Pol system is defined as

$$ 
\begin{aligned}
\dot{x} &= y, \\
\dot{y} &= \mu(1 - x^2)y - x,
\end{aligned} 
\qquad (B1) $$

with parameter $\mu = 2$ in this study. The system was integrated using an ODE solver over $t \in [0, 25]$ with a timestep $\Delta t = 0.001$. Initial conditions $[x(0), y(0)] = [2, 0]$ were chosen. The resulting trajectory $\{x(t), y(t)\}$ serves as the input data for SINDy. We employ the MIOSR algorithm as the optimizer within the PySINDy framework. A third-order polynomial feature library was used. The recovered model is given, which yields equations closely matching the analytical van der Pol dynamics:

$$ 
\begin{aligned}
\dot{x} &= 1.000y, \\
\dot{y} &= -1.000x + 2.000y - 2.000x^2y
\end{aligned} 
\qquad (B2) $$

To validate, the identified SINDy model was simulated from new initial conditions $[x(0), y(0)] = [1, -1]$. Figure 13 compares the phase portrait and time series of the true Van der Pol system with those predicted by the SINDy model. The agreement confirms that SINDy with MIOSR successfully captures the nonlinear dynamics.

## REFERENCES

¹J. S. Nelson, T. C. Grande, and M. V. Wilson, *Fishes of the World* (John Wiley & Sons, 2016).

²G. V. Lauder, “Fish locomotion: recent advances and new directions,” Annual review of marine science **7**, 521–545 (2015).

³M. Sfakiotakis, D. M. Lane, and J. B. C. Davies, “Review of fish swimming modes for aquatic locomotion,” *IEEE Journal of oceanic engineering* **24**, 237–252 (1999).

⁴M. S. Triantafyllou, G. S. Triantafyllou, and D. K. P. Yue, “Hydrodynamics of fishlike swimming,” *Annual Review of Fluid Mechanics* **32**, 33–53 (2000).

⁵J. Gray, “Directional control of fish movement,” *Proceedings of the Royal Society of London. Series B, Containing Papers of a Biological Character* **113**, 115–125 (1933).

⁶M. J. Lighthill, “Hydromechanics of aquatic animal propulsion,” *Annual Review of Fluid Mechanics* **1**, 413–446 (1969).

⁷K. Menon and R. Mittal, “Quantitative analysis of the kinematics and induced aerodynamic loading of individual vortices in vortex-dominated flows: a computation and data-driven approach,” *Journal of computational physics* **443**, 110515 (2021).

⁸G. V. Lauder and P. G. Madden, “Advances in comparative physiology from high-speed imaging of animal and fluid motion,” *Annu. Rev. Physiol.* **70**, 143–163 (2008).

⁹E. Standen, “Muscle activity and hydrodynamic function of pelvic fins in trout (oncorhynchus mykiss),” *Journal of Experimental Biology* **213**, 831–841 (2010).

¹⁰G. V. Lauder and E. D. Tytell, “Hydrodynamics of undulatory propulsion,” *Fish physiology* **23**, 425–468 (2005).

¹¹S. Alben, C. Witt, T. V. Baker, E. Anderson, and G. V. Lauder, “Dynamics of freely swimming flexible foils,” *Physics of Fluids* **24** (2012).

¹²F. E. Fish, “Darwin’s devices: What evolving robots can teach us about the history of life and the future of technology. john long.” (2012).

¹³L. Wen and G. Lauder, “Understanding undulatory locomotion in fishes using an inertia-compensated flapping foil robotic device,” *Bioinspiration & biomimetics* **8**, 046013 (2013).

¹⁴J. E. Colgate and K. M. Lynch, “Mechanics and control of swimming: A review,” *IEEE journal of oceanic engineering* **29**, 660–673 (2004).

¹⁵H.-B. Deng, Y.-Q. Xu, D.-D. Chen, H. Dai, J. Wu, and F.-B. Tian, “On numerical modeling of animal swimming and flight,” *Computational Mechanics* **52**, 1221–1242 (2013).

¹⁶A. P. Maertens, A. Gao, and M. S. Triantafyllou, “Optimal undulatory swimming for a single fish-like body and for a pair of interacting swimmers,” *Journal of Fluid Mechanics* **813**, 301–345 (2017).
<figure><img src="image_14.png" /><figcaption>FIG. 12: Polynomial fits of the model coefficients as functions of Strouhal number $St$ in the range $0.3 \le St \le 0.5$. Scatter points represent simulation data, dashed curves indicate cubic fits, and dotted curves indicate quartic fits (for $\alpha, \beta$, and $\mu$). For $\omega$, a cubic polynomial was found sufficient. The figure highlights the excellent agreement of the polynomial models with the data, consistent with the $R^2$ values reported in Tables III IV.</figcaption></figure>

<table><caption>TABLE III: Cubic polynomial models for $\alpha, \beta, \mu$, and $\omega$ as a function of Strouhal number $St$ (valid for $0.3 \le St \le 0.5$).</caption><thead><tr><th>Cubic Models</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>α(St) = 144.667St<sup>3</sup> - 168.629St<sup>2</sup> + 67.521St - 8.442</td><td>0.9987</td></tr><tr><td>β(St) = -1142.667St<sup>3</sup> + 1392.429St<sup>2</sup> - 536.636St + 71.578</td><td>0.9867</td></tr><tr><td>μ(St) = -26.000St<sup>3</sup> + 29.800St<sup>2</sup> - 11.985St + 1.473</td><td>0.9969</td></tr><tr><td>ω(St) = 54.288St<sup>3</sup> - 65.796St<sup>2</sup> + 57.452St - 3.565</td><td>0.9999</td></tr></tbody></table>

<footnote n="17">¹⁷L. Sirovich, “Turbulence and the dynamics of coherent structures. i. coherent structures,” Quarterly of applied mathematics **45**, 561–571 (1987).</footnote>

<footnote n="18">¹⁸P. J. Schmid, “Dynamic mode decomposition of numerical and experimental data,” Journal of fluid mechanics **656**, 5–28 (2010).</footnote>

<footnote n="19">¹⁹B. R. Noack, M. Morzynski, and G. Tadmor, *Reduced-order modelling for flow control*, Vol. 528 (Springer Science & Business Media, 2011).</footnote>

<footnote n="20">²⁰K. Carlberg, R. Tuminaro, and P. Boggs, “Preserving lagrangian structure in nonlinear model reduction with application to structural dynamics,” SIAM Journal on Scientific Computing **37**, B153–B184 (2015).</footnote>

<footnote n="21">²¹K. Carlberg, M. Barone, and H. Antil, “Galerkin v. least-squares petrov-galerkin projection in nonlinear model reduction,” Journal of Computational Physics **330**, 693–734 (2017).</footnote>

<footnote n="22">²²K. Menon and R. Mittal, “Dynamic mode decomposition based analysis of flow over a sinusoidally pitching airfoil,” Journal of Fluids and Structures</footnote>

<footnote>94, 102886 (2020).</footnote>

<footnote n="23">²³H. Gao and M. Wei, “Domain decomposition in pod-galerkin projection for flows with moving boundary,” in *54th AIAA Aerospace Sciences Meeting* (2016) p. 1102.</footnote>

<footnote n="24">²⁴B. Xu, H. Gao, M. Wei, and J. Hrynuk, “Global pod-galerkin roms for fluid flows with moving solid structures,” AIAA Journal **60**, 1400–1414 (2022).</footnote>

<footnote n="25">²⁵A. Saeed, H. Saeed, and J. Li, “Data-driven reduced-order modeling for undulating swimmers,” Physics of Fluids **37** (2025).</footnote>

<footnote n="26">²⁶R. Tong, Y. Feng, J. Wang, Z. Wu, M. Tan, and J. Yu, “A survey on reinforcement learning methods in bionic underwater robots,” Biomimetics **8**, 168 (2023).</footnote>

<footnote n="27">²⁷P. Gunnarson, I. Mandralis, G. Novati, P. Koumoutsakos, and J. O. Dabiri, “Learning efficient navigation in vortical flow fields,” Nature Communications **12**, 7143 (2021).</footnote>
<table><caption>TABLE IV: Quartic polynomial models for $\alpha, \beta, \mu$, and $\omega$ function of Strouhal number $St$ (valid for $0.3 \le St \le 0.5$).</caption><thead><tr><th>Quartic Models</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>α(St) = 1046.667St<sup>4</sup> - 1530.000St<sup>3</sup> + 824.583St<sup>2</sup> - 191.155St + 16.512</td><td>1.00</td></tr><tr><td>β(St) = 19906.667St<sup>4</sup> - 32993.333St<sup>3</sup> + 20282.433St<sup>2</sup> - 5452.426St + 546.181</td><td>1.00</td></tr><tr><td>μ(St) = -420.000St<sup>4</sup> + 646.000St<sup>3</sup> - 368.750St<sup>2</sup> + 91.815St - 8.540</td><td>1.00</td></tr><tr><td>ω(St) = -832.0610St<sup>4</sup> + 1385.5857St<sup>3</sup> - 855.3624St<sup>2</sup> + 263.0899St - 23.4028</td><td>1.00</td></tr></tbody></table>

<table><caption>TABLE V: Cubic polynomial models for $\eta_0, \eta_1$, and $\eta_2$ as a function of $St$.</caption><thead><tr><th>Cubic Models</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>η<sub>0</sub>(St) = -1.6302 + 13.0667St - 34.2000St<sup>2</sup> + 25.3333St<sup>3</sup></td><td>1.000</td></tr><tr><td>η<sub>1</sub>(St) = 2.0669 - 15.9926St + 39.5429St<sup>2</sup> - 32.6667St<sup>3</sup></td><td>0.999</td></tr><tr><td>η<sub>2</sub>(St) = -0.0791 + 0.6581St - 1.5143St<sup>2</sup> + 1.3333St<sup>3</sup></td><td>0.999</td></tr></tbody></table>

<footnote n="28"><sup>28</sup>C. Rodwell and P. Tallapragada, “Physics-informed reinforcement learning for motion control of a fish-like swimming robot,” Scientific reports 13, 10754 (2023).</footnote>

<footnote n="29"><sup>29</sup>S. Brunton, “Discovering governing equations from data by sparse identification of nonlinear dynamics,” in APS March Meeting Abstracts, Vol. 2017 (2017) pp. X49–004.</footnote>

<footnote n="30"><sup>30</sup>D. Bertsimas and W. Gurnee, “Learning sparse nonlinear dynamics via mixed-integer optimization,” Nonlinear Dynamics 111, 6585–6604 (2023).</footnote>

<footnote n="31"><sup>31</sup>J.-C. Loiseau and S. L. Brunton, “Constrained sparse galerkin regression,” Journal of Fluid Mechanics 838, 42–67 (2018).</footnote>

<footnote n="32"><sup>32</sup>R. Tibshirani, “Regression shrinkage and selection via the lasso,” Journal of the Royal Statistical Society Series B: Statistical Methodology 58, 267–288 (1996).</footnote>

<footnote n="33"><sup>33</sup>P. Zheng, T. Askham, S. L. Brunton, J. N. Kutz, and A. Y. Aravkin, “A unified framework for sparse relaxed regularized regression: Sr3,” IEEE Access 7, 1404–1423 (2018).</footnote>

<footnote n="34"><sup>34</sup>J.-C. Loiseau, B. R. Noack, and S. L. Brunton, “Sparse reduced-order modelling: sensor-based dynamics to full-state estimation,” Journal of Fluid Mechanics 844, 459–490 (2018).</footnote>

<footnote n="35"><sup>35</sup>K. Fukami, T. Murata, K. Zhang, and K. Fukagata, “Sparse identification of nonlinear dynamics with low-dimensionalized flow representations,” Journal of Fluid Mechanics 926, A10 (2021).</footnote>

<footnote n="36"><sup>36</sup>E. D. Tytell and G. V. Lauder, “The hydrodynamics of eel swimming: I. wake structure,” Journal of Experimental Biology 207, 1825–1841 (2004).</footnote>

<footnote n="37"><sup>37</sup>W. M. Van Rees, M. Gazzola, and P. Koumoutsakos, “Optimal shapes for anguilliform swimmers at intermediate reynolds numbers,” Journal of Fluid Mechanics 722, R3 (2013).</footnote>

<footnote n="38"><sup>38</sup>G. D. Weymouth and B. Font, “Waterlily. jl: A differentiable fluid simulator in julia with fast heterogeneous execution,” arXiv preprint arXiv:2304.08159 (2023).</footnote>

<footnote n="39"><sup>39</sup>G. D. Weymouth and D. K. Yue, “Boundary data immersion method for cartesian-grid simulations of fluid-body interaction problems,” Journal of Computational Physics 230, 6233–6247 (2011).</footnote>

<footnote n="40"><sup>40</sup>A. A. Kaptanoglu, B. M. de Silva, U. Fasel, K. Kaheman, A. J. Goldschmidt, J. Callaham, C. B. Delahunt, Z. G. Nicolaou, K. Champion, J.-C. Loiseau, J. N. Kutz, and S. L. Brunton, “Pysindy: A comprehensive python package for robust sparse system identification,” Journal of Open Source Software 7, 3994 (2022).</footnote>

<footnote n="41"><sup>41</sup>R. Chartrand, “Numerical differentiation of noisy, nonsmooth data,” International Scholarly Research Notices 2011, 164564 (2011).</footnote>

<footnote n="42"><sup>42</sup>B. de Silva, K. Champion, M. Quade, J.-C. Loiseau, J. Kutz, and S. Brunton, “Pysindy: A python package for the sparse identification of nonlinear dynamical systems from data,” Journal of Open Source Software 5, 2104 (2020).</footnote>

<footnote n="43"><sup>43</sup>G. Schwarz, “Estimating the dimension of a model,” The annals of statistics, 461–464 (1978).</footnote>

<footnote n="44"><sup>44</sup>H. Akaike, “Information theory and an extension of the maximum likelihood principle,” in Selected papers of hirotugu akaike (Springer, 1998) pp. 199–213.</footnote>

<footnote n="45"><sup>45</sup>N. M. Mangan, J. N. Kutz, S. L. Brunton, and J. L. Proctor, “Model selection for dynamical systems via sparse regression and information criteria,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 473, 20170009 (2017).</footnote>

<footnote n="46"><sup>46</sup>A. H. Nayfeh, F. Owis, and M. R. Hajj, “A model for the coupled lift and drag on a circular cylinder,” in International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, Vol. 37033 (2003) pp. 1289–1296.</footnote>

<footnote n="47"><sup>47</sup>R. T. Hartlen and I. G. Currie, “Lift-oscillator model of vortex-induced vibration,” Journal of the Engineering Mechanics Division 96, 577–591 (1970).</footnote>

<footnote n="48"><sup>48</sup>W. Iwan and R. Blevins, “A model for vortex induced oscillation of structures,” Journal of applied mechanics 41, 581–586 (1974).</footnote>

<footnote n="49"><sup>49</sup>N. Mureithi, S. Goda, H. Kanki, and T. Nakamura, “A nonlinear dynamics analysis of vortex-structure interaction models,” J. Pressure Vessel Technol. 123, 475–479 (2001).</footnote>

<footnote n="50"><sup>50</sup>R. Skop, O. Griffin, and G. Koopman, “The vortex-excited resonant vibrations of circular cylinders,” Journal of Sound and Vibration 31, 235–249 (1973).</footnote>

<footnote n="51"><sup>51</sup>R. Skop and O. Griffin, “On a theory for the vortex-excited oscillations of flexible cylindrical structures,” Journal of Sound and Vibration 41, 263–274 (1975).</footnote>

<footnote n="52"><sup>52</sup>R. Landl, “A mathematical model for vortex-excited vibrations of bluff bodies,” Journal of Sound and Vibration 42, 219–234 (1975).</footnote>

<footnote n="53"><sup>53</sup>O. Marzouk, A. H. Nayfeh, I. Akhtar, and H. N. Arafat, “Modeling steady-state and transient forces on a cylinder,” Journal of Vibration and Control 13, 1065–1091 (2007).</footnote>

<footnote n="54"><sup>54</sup>M. S. U. Khalid, J. Wang, H. Dong, and M. Liu, “Flow transitions and mapping for undulating swimmers,” Physical Review Fluids 5, 063104 (2020).</footnote>

<footnote n="55"><sup>55</sup>I. Borazjani and F. Sotiropoulos, “Numerical investigation of the hydrodynamics of anguilliform swimming in the transitional and inertial flow regimes,” Journal of Experimental Biology 212, 576–592 (2009).</footnote>

<footnote n="56"><sup>56</sup>B. R. Noack, K. Afanasiev, M. Morzyński, G. Tadmor, and F. Thiele, “A hierarchy of low-dimensional models for the transient and post-transient cylinder wake,” Journal of Fluid Mechanics 497, 335–363 (2003).</footnote>

<footnote n="57"><sup>57</sup>L. Qin, “Development of reduced-order models for lift and drag on oscillating cylinders with higher-order spectral moments,” (2004).</footnote>

<footnote n="58"><sup>58</sup>W.-J. Kim and N. Perkins, “Two-dimensional vortex-induced vibration of cable suspensions,” Journal of Fluids and Structures 16, 229–245 (2002).</footnote>

<footnote n="59"><sup>59</sup>N. Thekkethil, A. Sharma, and A. Agrawal, “Unified hydrodynamics study for various types of fishes-like undulating rigid hydrofoil in a free stream flow,” Physics of Fluids 30 (2018).</footnote>

<footnote n="60"><sup>60</sup>Y.-L. Yu and K.-J. Huang, “Scaling law of fish undulatory propulsion,” Physics of Fluids 33 (2021).</footnote>

<footnote n="61"><sup>61</sup>M. I. Uddin, G. A. Garcia, and O. M. Curet, “Force scaling and efficiency of elongated median fin propulsion,” Bioinspiration & Biomimetics 17, 046004 (2022).</footnote>

<footnote n="62"><sup>62</sup>D. Yang, J. Wu, K. Khedkar, L.-M. Chao, and A. P. S. Bhalla, “Hydrodynamics and scaling laws for intermittent s-start swimming,” Journal of Fluid Mechanics 984, A2 (2024).</footnote>
<table><caption>TABLE VI: Quartic polynomial models for $\eta_0$, $\eta_1$, and $\eta_2$ as a function of $St$.</caption><thead><tr><th>Quartic Models</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>$\eta_0(St) = 0.5950 - 10.0000St + 54.3667St^2 - 124.0000St^3 + 93.3333St^4$</td><td>1.00</td></tr><tr><td>$\eta_1(St) = 1.9080 - 14.3450St + 33.2167St^2 - 22.0000St^3 - 6.6667St^4$</td><td>1.00</td></tr><tr><td>$\eta_2(St) = -0.3970 + 3.9533St - 14.1667St^2 + 22.6667St^3 - 13.3333St^4$</td><td>1.00</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>FIG. 13: Verification of SINDy on the van der Pol oscillator. Left: Phase portrait of the true system vs. the SINDy model. Right: Time series of $x(t)$.</figcaption></figure>