Graphical Abstract One-shot learning for the complex dynamical behaviors of weakly nonlinear forced oscillators
arXiv:2604.15181v1 [cs.LG] 16 Apr 2026
Teng Ma, Luca Rosafalco, Wei Cui, Lin Zhao, Attilio Frangi
Highlights One-shot learning for the complex dynamical behaviors of weakly nonlinear forced oscillators Teng Ma, Luca Rosafalco, Wei Cui, Lin Zhao, Attilio Frangi • A one-shot learning method identifies global nonlinear frequency-response curves of weakly nonlinear forced oscillators from a single excitation time history by inferring the governing equations. • The equation learning framework for weakly nonlinear systems is upgraded from autonomous single-frequency dynamics to non-autonomous multi-frequency dynamics. This advancement is achieved by incorporating the generalized harmonic balance theory.
One-shot learning for the complex dynamical behaviors of weakly nonlinear forced oscillators Teng Maa,b , Luca Rosafalcob,∗, Wei Cuia , Lin Zhaoa,c , Attilio Frangib a
State Key Lab of Disaster Reduction in Civil Engineering, Tongji University, Shanghai, 200092, Shanghai, P. R. China b Department of Civil and Environmental Engineering, Politecnico di Milano, Milan, 20133, Lombardia, Italy c Guangxi Laboratory of Whole Life Safety for Land-Sea Corridor Engineering, Guangxi University, Nanning, 530004, Guangxi, P. R. China
Abstract Extrapolative prediction of complex nonlinear dynamics remains a central challenge in engineering. This study proposes a one-shot learning method to identify global frequencyresponse curves from a single excitation time history by learning governing equations. We introduce MEv-SINDy (Multi-frequency Evolutionary Sparse Identification of Nonlinear Dynamics) to infer the governing equations of non-autonomous and multi-frequency systems. The methodology leverages the Generalized Harmonic Balance (GHB) method to decompose complex forced responses into a set of slow-varying evolution equations. We validated the capabilities of MEv-SINDy on two critical Micro-Electro-Mechanical Systems (MEMS). These applications include a nonlinear beam resonator and a MEMS micromirror. Our results show that the model trained on a single point accurately predicts softening/hardening effects and jump phenomena across a wide range of excitation levels. This approach significantly reduces the data acquisition burden for the characterization and design of nonlinear microsystems. Keywords: System identification, Nonlinear dynamics, Weak nonlinearity, Generalized harmonic balance, One-shot learning
∗
email: [email protected]
1. Introduction The rapid advancement of computational resources has enabled the high-fidelity simulation of complex physical phenomena, even within intricate multi-physics and multi-scale frameworks. However, characterizing the complex dynamical behaviors of physical systems remains a formidable challenge. In practical engineering applications, such as MEMS resonators and sensors, understanding the nonlinear response—including frequency shifts, hardening/softening effects, and jump phenomena—requires extensive frequency and amplitude sweeps. When relying on full-order models (FOMs), such as the finite element method (FEM), each individual simulation in the time domain or parametric space can be computationally very demanding. Consequently, exploring a broad range of initial conditions and parameter combinations becomes a computationally prohibitive task or even unfeasible for real-time monitoring and design optimization. Unlike simple output estimation, capturing the evolution of a high-dimensional nonlinear field is intrinsically difficult due to its time-dependent nature and the sensitivity of nonlinear regimes. These challenges drive the urgent need for efficient yet accurate reduced-order models (ROMs) that can bypass the heavy burden of FOMs while preserving the essential physics of the dynamical system. Efficient ROM construction begins with dimensionality reduction (Van Der Maaten et al., 2009). This process projects high-dimensional FOM states onto a low-dimensional manifold from data previously collected from full-order simulations or experiments (Huang et al., 2017). Historically, the reduced basis method has been the standard for building these reduced spaces (Maday and Rønquist, 2002; Quarteroni et al., 2015; Benner et al., 2017). It often utilizes Proper Orthogonal Decomposition (POD) (Amsallem et al., 2012; Pagani et al., 2018; Gobat et al., 2022) to approximate solutions. However, POD-based techniques are typically intrusive. They require direct access to the internal operators of the FOM. This requirement makes it difficult to use them with proprietary or commercial simulation software. Machine learning has emerged as a powerful non-intrusive alternative for ROM construction. These methods reduce dimensionality directly from data streams. They do not require access to the FOM operators. Dynamic Mode Decomposition (DMD) (Schmid et al., 2011; Le Clainche and
2
Vega, 2017) is a prominent linear framework in this category. It extracts dominant frequencies and spatial modes from snapshot data. However, DMD is limited by its linear assumption. Autoencoder (AE) networks (Zhai et al., 2018; Romor et al., 2023; Gonzalez and Balajewicz, 2018) provide a nonlinear alternative. They are particularly effective at uncovering latent features in complex systems. These networks offer better nonlinear compression than linear methods like POD or DMD. Advanced frameworks like deep learning-based ROMs (DLROMs) (Fresca et al., 2022; Franco et al., 2023) and their enhanced version through POD (POD-DL-ROMs) (Fresca and Manzoni, 2022) have further improved this process. They perform dimensionality reduction and learn the parameter-to-solution map at the same time. Data-driven ROMs face significant challenges when predicting scenarios outside the training set. These methods are essentially interpolation tools. They provide high accuracy only within the observed parameter ranges. When these models encounter unseen parameter values, their performance often degrades. This lack of generalization limits their use in exploring unknown physical regimes. Furthermore, constructing reliable non-intrusive ROMs requires a large volume of training data. Generating these datasets through FOM simulations is computationally expensive. This high cost reduces the overall efficiency of the modeling process. There is a clear need for a framework that offers better extrapolation. Such a framework should also require less training data to maintain physical consistency. Equation discovery offers a powerful alternative to black-box models. This approach infers governing equations directly from data. There are two main methodological paradigms in this field. The first paradigm is symbolic discovery (Schmidt and Lipson, 2009; Ruan et al., 2025). It seeks to recover closed-form representations of physical laws. The second paradigm is sparse regression (Brunton et al., 2016; Gao and Yan, 2022). This method identifies parsimonious models from a predefined library of candidate functions. These techniques were first developed for ordinary differential equations. Researchers later extended them to partial differential equations (Rudy et al., 2017; Rao et al., 2023). Recent advances in machine learning have improved these methods. They are now robust against sparse or noisy data (Chen et al., 2021; Hirsh et al., 2022; Niven et al., 2024; Fung et al., 2025). This makes them viable for real-world scientific discovery. Equation discovery has catalyzed the field of scientific 3
machine learning. It helps uncover interpretable laws in cyber-physical systems (Yuan et al., 2019), fluid dynamics (Ma et al., 2023; Loiseau et al., 2018), biological and chemical networks (Mangan et al., 2016; Hoffmann et al., 2019) and complex networks (Gao et al., 2024; Yu et al., 2025). The main appeal of this approach is the generation of an explicit ROM. The resulting model takes the form of a system of ordinary differential equations. Sparsity constraints ensure that the latent dynamics remain interpretable. Users can then analyze the identified system using standard numerical tools for time integration. Most physical and engineering systems exhibit weakly nonlinear oscillatory behaviors with weak linear dissipation. These include applications in structural engineering and microelectro-mechanical systems (MEMS). Weak nonlinearities often arise from geometric effects. They can lead to pronounced impacts on dynamical performance. However, existing methods often overlook these challenges. They focus on extracting only the most dominant governing equations. This process often discards terms that seem insignificant but are physically important. Recent research has introduced evolutionary-based learning to address this issue. A prominent example is EvLOWN (Ma et al., 2026), which stands for Encoding Cumulation to Learn Perturbative Nonlinear Oscillatory Dynamics. This method reformulates optimization objectives using the method of averaging. It bridges the gap between weak nonlinearities and slow-varying evolutionary variables. Only the terms most informative about amplitude and phase evolutions are selected. However, EvLOWN has specific limitations. It can only identify autonomous single-frequency responses but cannot handle multi-frequency oscillations or non-autonomous dynamics. In this work, we propose a novel data-driven approach. We call it Multi-frequency Evolutionary Sparse Identification of Nonlinear Dynamical systems (MEvSINDy). This framework is designed to learn the governing equations of multi-frequency forced weakly nonlinear oscillators. It extends the advantages of evolutionary learning to more complex forced and multi-frequency regimes. The capability to identify governing equations from multi-frequency data opens the door to a more ambitious paradigm. This paradigm is known as one-shot learning (Fei-Fei et al., 2006). It involves training models to recognize global patterns based on only a single data point. The concept originated in computer vision but remains rare in scientific machine 4
learning (SciML) (Ivanov et al., 2020; Darcy et al., 2023; Jiao et al., 2025). Most data-driven ROMs still require massive datasets to characterize physical systems (Conti et al., 2023; Gobat et al., 2022; Franco et al., 2023). However, generating such snapshots for complex nonlinear oscillators is computationally expensive. We leverage MEv-SINDy to achieve oneshot learning for forced nonlinear dynamics. In our framework, the identified governing equations represent the “general knowledge” of the system. We extract this knowledge from a single excitation configuration. This analytical representation allows the model to predict dynamical behaviors under entirely different parameters. Consequently, MEv-SINDy can predict the complete frequency-response curve from one observation. This strategy eliminates the need for exhaustive training data. It also provides a robust tool for exploring unknown physical regimes with minimal cost. To our knowledge, this is one of the first applications of one-shot learning for multi-frequency forced oscillators. Schematic diagram of this work is shown in Figure 1. The paper is structured as follows. Section 2 details the methodology of the proposed MEv-SINDy framework. This section introduces the Generalized Harmonic Balance (GHB) method as an upgrade to the averaging method. The evolutionary variables are expanded from simple amplitude and phase to general harmonic coefficients. Based on GHB, we reformulate the sparse regression optimization targets. We also describe a joint sparse regression method designed for multiple equations. Section 3 presents numerical validations using a 1-DOF oscillator with quadratic and cubic nonlinearities. These tests validate the effectiveness of each step in the MEv-SINDy process. Section 4 considers two high-dimensional applications. These include a Beam MEMS resonator and a MEMS micromirror. We demonstrate that our approach achieves one-shot learning for the frequency-response curve. To further analyze the robustness of the method, we examine the influence of training point selection in Section 4. We describe a multi-point training strategy to improve predictive performance. Finally, Section 5 provides concluding remarks and summarizes the findings.
5
Figure 1: One-Shot Learning framework of weakly nonlinear forced oscillator: The starting point of our framework is the dynamical response at a single excitation configuration, which is also the training data of MEv-SINDy. During the training stage, the multi-frequency dynamics are separated into distinct harmonic components. We use the generalized harmonic balance method to reformulate the sparse regression problem into an evolution space. A two-stage sparse regression process then discovers the final coefficient vector. The primary output of MEv-SINDy is the explicit governing equation of the weakly nonlinear forced oscillator. Finally, we apply a continuation method to the inferred equations to compute the global frequency-response curves.
2. Methodology 2.1. Problem setup The starting point is to consider the general form of weakly nonlinear oscillators with external forcing: ẍi (t) + ωi2 xi (t) + ϵfi (x(t), ẋ(t)) = Fi (t),
(1)
for i = 1, . . . , k,where: x(t) = (x1 (t), · · · , xk (t)) ∈ Rk is the displacement vector at time t; ωi is the frequency of oscillator i; fi (x(t), ẋ(t)) is a smooth function of the velocity ẋ and the displacement x representing the weak nonlinearity and weak dissipation, which are from 6
geometric, internal, material, contact and coupling nonlinearities; Fi (t) is the external actuation that depends on time t; ϵ is a small parameter quantifying the strength of nonlinearity compared to dominating linear spring effects ωi2 xi . Eq. (1) is well suited for modeling various classes of systems, whose dimensionality k represents the number of degrees of freedom arising, for example, from spatial discretization techniques (e.g. finite elements or finite volumes) applied to a system of partial differential equations (PDEs) governing the underlying physical problem. Due to the weak nonlinearity fi , the steady-state periodic response of the system governed by Eq. (1) under different harmonic actuation Fi exhibits distinct characteristics in terms of amplitude and phase, collectively referred to as the frequency-response curves (FRCs), which characterize the dynamic properties of the system. Typically, the actuation can be parameterized by two features, frequency and intensity, both of which strongly affect the steady-state periodic response. The goal of this study is to utilize data snapshots from a single time history and a given actuation to accurately predict responses under any other actuation conditions, thus reconstructing all the possible FRCs of the system. The proposed methodology employs only few set of time history of displacement x and its time-derivative ẋ. In general, ẋ can be either computed directly or approximated numerically from x. 2.2. Sparse Identification of Nonlinear Dynamics: SINDy To identify the system in Eq. (1), it is possible to rely on a linear combination of a set of predetermined functions collected in a library, as done in the SINDy method, (Brunton et al., 2016). To apply SINDy, snapshots of x at different time instants are collected in the following matrix:
x (t ) x2 (t1 ) 1 1 x1 (t2 ) x2 (t2 ) X= .. .. . . x1 (tT ) x2 (tT )
··· ··· .. . ···
xn (t1 )
xn (t2 ) .. . . xn (tT )
(2)
where: xi (tj ) is the i-th entry of x at the j-th time instant; T is the number of time instants for each snapshot. Similarly, a matrix Ẋ is constructed by collecting ẋ. 7
A library Θ(x) = [θ1 (x), . . . , θp (x)] ∈ Rp of p candidate functions is selected. In Sec. 2.7, it will be discussed how to retain only the most relevant functions among this set of candidates. The matrix Θ(X) ∈ RT ×p is thus constructed by applying Θ to the rows of X, for example, as in the following: | | | | Θ(X) = 1 X XP2 XP3 · · · | | | |
|
|
sin(X) cos(X) · · · . | |
(3)
Any nonlinear function, such as polynomial and trigonometric, cna be included in the library. As mentioned, the model to be identified is a linear combination of the candidate functions, with combination coefficients stored in Ξ = [ξ1 , . . . , ξk ] with ξk ∈ Rp . Combination coefficients are determined by solving a regression problem with a sparsity promoting term, as in the following: arg min ∥ẍ − Θ(x, ẋ, t)Ξ∥22 + λ∥Ξ∥1 Ξ
(4)
where || · ||22 is the least squares term; || · ||1 is the L1 norm of the coefficient vector Ξ, which is the sum of absolute values of the coefficients; λ is the regularization parameter that controls the strength of the penalty. In Sec. 2.7, Eq. (4) will be modified for a better coupling with the proposed algorithm. 2.3. Weakly nonlinear oscillator inference: EvLOWN It is worth noting that most existing model identification frameworks such SINDy do not target the specific difficulties introduced by weak nonlinearities in practical systems. Indeed, these methods identify parsimonious governing equations by neglecting library terms weighted by very small terms. While this approach is meaningful in general, it may become a major drawback in the context of WNOs, where the very necessary components for accurately capturing the dynamics may be discarded as the dynamics is governed by the linear term ωi2 xi . 8
To address this gap, Ma et al. (2026) proposed a data-driven approach called EvLOWN (Evolutionary Learning Oscillator with Weak Nonlinearity) for discovering weakly nonlinear oscillator. EvLOWN reformulates the original system of ODEs into an evolutionary-variable representation via averaging theory, naturally separating terms by their order of magnitude. The goal of EvLOWN is to discover the autonomous oscillatory dynamical system: ẍi (t) + ωi2 xi (t) + ϵfi (x(t), ẋ(t)) = 0
(5)
by learning the evolution of amplitude and phase, rather than directly learning the state variables x. Leveraging averaging theory (Sanders et al., 2007), also known as the Krylov-BogoliubovMitropolsky method (Volosov, 1962), we approximate the response of the oscillator (5) when the nonlinear effect is weak (i.e., ϵ is small) as: xi (t) = Ai (t) sin(ω̂i t + ϕi (t))
(6)
where: ω̂i is the oscillating frequency; Ai and ϕi are the amplitude and phase, respectively. In WNOs, the amplitude Ai (t) and phase ϕi (t) evolve much more slowly over time compared to the oscillatory term, changing little over a single period 2π/ω̂i Nayfeh and Balachandran (2008). This separation of timescales allows us to average the system dynamics over one oscillation cycle, yielding explicit relationships between the weak nonlinearities and the time evolution of the amplitude Ai (t) and phase ϕi (t): Z t+π/ω̂i ϵi Ai (t + π/ω̂i ) − Ai (t − π/ω̂i ) ′ Ai (t) = fi (x(s), ẋ(s)) cos(ω̂i s + ϕi (s))ds, = 2π/ω̂i 2π t−π/ω̂i Z t+π/ω̂i ϕi (t + π/ω̂i ) − ϕi (t − π/ω̂i ) ϵi ′ ϕi (t) = =− fi (x(s), ẋ(s)) sin(ω̂i s + ϕi (s))ds, 2π/ω̂i 2Ai (t)π t−π/ω̂i (7) where (·)′ denotes the rate of change over one period. To identify the first-order ODEs governing the evolution of amplitude and phase, we integrate the basis terms in the original library Θ over a single period, which quantifies their contributions to the averaged dynamics. This transformation results in two new libraries: ΘA for amplitude evolution and Θϕ for phase evolution. The discovery task now becomes 9
identifying a sparse subset of basis terms from these libraries that best approximate A′i (t) and ϕ′i (t). Critically, this averaging process naturally separates terms by order of magnitude, allowing weak nonlinear contributions to be isolated from dominant dynamics. 2.4. Generalized Harmonic Balance method As the magnitude of external forcing increases in the forced oscillator Eq.
(1), the
contribution of nonlinearities in the system, such as quadratic and cubic terms, becomes significant. Consequently, the assumption of a single-frequency solution, see Eq. (8) becomes inadequate for complex forced dynamics. To address this issue, we employ the generalized harmonic balance method to replace the averaging method, thus extending the EvLOWN approach to account for multiple-frequency responses. The solution of the weakly nonlinear forcing oscillator can be approximated as follows (Luo and Huang, 2012): (0) x∗i (t) = ai (t) +
N X (n) (n) {bi (t) cos(nω̂i t) + ci (t) sin(nω̂i t)}
(8)
n=1
where the superscript (·)(n) indicates the order of the harmonic frequency, corresponding to (n)
(0)
(n)
the n-multiple of the fundamental frequency ω̂i . The terms ai (t), bi (t), and ci (t) are the slowly varying parameters associated with each harmonic component nω̂i . The first- and second-order time derivatives of x∗ (t) are: (0) ẋ∗i (t) = ȧi (t) +
N X
(n)
(n)
(n)
{[ḃi (t) + nω̂i ci (t)] cos(nω̂i t) + [ċi
(n)
− nω̂i bi (t)] sin(nω̂i t)}
(9)
n=1
(0) ẍ∗i (t) = äi (t) +
N X
(n)
(n)
(n)
{[b̈i (t) + 2nω̂i ċi (t) − (nω̂i )2 bi (t)] cos(nω̂i t)
n=1 (n)
(n)
(10)
(n)
+[c̈i (t) − 2nω̂i ḃi (t) − (nω̂i )2 ci (t)] sin(nω̂i t)} (0)
(n)
(n)
Supposing that ai (t), bi (t), and ci (t) vary slowly with time, we substitute Eqs. (8)(10) into Eq. (1), and then we average for each harmonic terms of cos(nωi t) and sin(nωi t) obtaining for n = 1, 2, · · · , N :
10
ω̂i (0) äi (t) + ωi2 a0 (t) = −
Z t+π/ω̂i
ω̂i π
Z t+π/ω̂i
2π
(n)
(n)
(n)
b̈i (t) + 2ω̂i nċi (t) − ω̂i2 (n2 − 1)bi (t) = − (n)
(n)
(n)
c̈i (t) − 2ω̂i nḃi (t) − ω̂i2 (n2 − 1)ci (t) = −
ω̂i π
[fi∗ (x, ẋ, s)]ds
t−π/ω̂i
t−π/ω̂i Z t+π/ω̂i
[fi∗ (x, ẋ, s)] cos(nω̂i s)ds
(11)
[fi∗ (x, ẋ, s)] sin(nω̂i s)ds,
t−π/ω̂i
where: (12)
fi∗ (x, ẋ, s) = ϵfi (x∗ (s), ẋ∗ (s)) − Fi (x∗ (s), ẋ∗ (s), s)
It is worth stressing that the only assumption is related to the slow temporal variation (0)
(n)
(n)
of ai (t), bi (t), and ci (t), while no limitations on the type of response (steady state or transient) is instead required. In the next section, we will show how to exploit the evolutionary formulation in Eq. (11) to identify the system described, in its original formulation, by Eq. (1). 2.5. Multiple slow-varying evolutions evaluations (0)
(n)
We now describe how to extract the slowly varying evolutionary variables ai (t), bi (t), (n)
and ci (t) (n = 1, 2, . . .) from the observed signal x(t). These variables are subsequently used to perform sparse regression on the evolutionary formulation reported in Eq. (11) in order to identify the unknown components of the governing equations, namely ωi , fi∗ . Given an input signal x(t) = [x1 (t), . . . , xd (t)]T , the oscillation frequency in the i-th dimension can be directly computed from the Fourier transform F(xi (t)) : R → R, with: Z +∞ ω̂i = arg max F(xi (t)) = arg max ω̂i
ω̂i
xi (t)e−iω̂i t dt
(13)
−∞
It is worth noting that the oscillation circular frequency ω̂i differs from the natural circular frequency ωi . This deviation arises from the presence of weakly nonlinear terms and external forcing. To illustrate this, we consider a simple example of a linearly damped oscillator described by ẍ + 2β ẋ + ω 2 x = F cos(Ωt). Under the assumption of light damping (β ≪ 1), the solution of the linear damped oscillator is given by: p f x(t) = x0 e−βt cos( ω 2 − β 2 t + ϕ) + p cos(Ωt − φ) (ω 2 − Ω2 )2 + 4β 2 ω 2 11
Therefore, the oscillation circular frequency ω̂ obtained from Eq. (13) lies within the range p ( ω 2 − β 2 , Ω), and the value of ω̂ depends on the selected portion of the data x(t). If the p transient response dominates, ω̂ will tend to ω 2 − β 2 ; otherwise, it will be close to Ω. (0)
(n)
(n)
To obtain the multiple slow-varying evolutionary variables ai (t), bi (t), and ci (t) (n = 1, 2, . . .) from the given time history xi (t), we proposed a novel identification method involving the Hilbert transform. Unlike what was proposed in (Ma et al., 2026), the Hilbert transform cannot be directly applied to extract the slowly varying amplitude and phase, because the dynamic response of a weakly nonlinear forced oscillator contains multiple frequency components. To address this, we define N frequency bands to address the limitations of standard Fourier analysis. In classical spectral analysis, isolated peaks represent constant coefficients. However, our evolutionary coefficients depend on time. This time-dependency smears the frequency peaks in the spectrum. In particular, we determine the value of N depends on the number of distinct frequency peaks observed in the frequency spectrum. When, as in Figure 2a, x(t) consists of the zero, first and second harmonic components, we set N = 3. Each frequency band to be centered at nω̂ and has a bandwidth of ω̂/4 (Figure 2b). By passing the signal through the corresponding bandpass filter, each harmonic component is isolated, as shown in Figure 2c. In practice, Finite Impulse Response (FIR) filters are designed using the window method. This procedure is implemented in readily available packages such as scipy.signal. (0)
(n)
(n)
After having obtained the slow-varying evolutionary variables ai (t), bi (t), and ci (t) through the Hilbert transform, we can elaborate Eq. (8) as: xi (t) =
N X
(n)
(14)
xi (t),
i=0
with: (n)
xi (t) =
a(0) (t)
n=0
i
(15)
b(n) (t) cos(nω̂t) + c(n) (t) sin(nω̂t) n = 1, 2, · · · i i
The zero-order harmonic component directly corresponds to a(0) (t). While for n > 0, the slow-varying evolutionary variables are calculated from amplitude and phase. The Hilbert (n)
(n)
(n)
transform H is used to compute the analytic signal x̂i (t) = xi (t) + iH(xi (t)). The 12
Figure 2:
Identification process of the slowly varying evolutionary variables a(0) (t), b(n) (t), and c(n) (t) by
combining a bandpass filter and the Hilbert transform. (a) Input signal x(t), representing the time-domain response of a weakly nonlinear forced oscillator; (b) x(t) consists of multiple narrowband harmonic components centered at nω̂, as shown in the frequency domain. Several frequency bands are defined to isolate each n-th harmonic component; (c) Time-domain representation and (d) corresponding frequency-domain representation of the filtered narrowband signals; (e) Hilbert transform applied to each n-th harmonic component to extract the instantaneous amplitude A(n) (t) and phase β (n) (t); (f) The slowly varying coefficients b(n) (t) and c(n) (t) are calculated from the extracted amplitude and phase, while the zero-order term directly corresponds to a(0) (t). (n)
(n)
evolutionary variables Ai (t) and βi (t) of time t can be considered as the averaging value of instantaneous evolutionary variables in one period:
13
R t+π/ω̂i (n) Ai (t) =
(n)
βi (t) =
t−π/ω̂i
(n)
|x̂i (s)|ds
2π/ω̂i R t+π/ω̂i (n) Arg(x̂i (s))ds t−π/ω̂i
(16)
2π/ω̂i (n)
(n)
As Figure 2f shows, the evolutionary variable bi (t) and ci (t) of higher harmonic components can be calculated as Eq. (17): (n)
(n)
(n)
(n)
(n)
(n)
bi (t) = Ai (t) sin(βi (t))
(17)
ci (t) = Ai (t) cos(βi (t)) 2.6. Reformulation of the Regression Problem in the Evolutionary Domain In order to accurately identify the governing dynamics of weakly nonlinear forced oscillators, it is essential to reformulate the regression problem in a form that enhances the visibility of nonlinear effects. Sharing the same starting point of EvLOWN, we transform the problem into an evolutionary domain. In the proposed MEv-SINDy approach, the system behavior is represented by slowly varying evolutionary variables based on generalized harmonic balance method. This reformulation allows the nonlinear interactions to be more effectively captured and facilitates the application of sparse regression for discovering the underlying equations. We thus exploit SINDy to approximate the evolutionary formulation reported in Eq. (11), instead of the original formulation (Eq. (1)). This is equivalent to perform regularized regression reported in Eq. (4), but here involving the evolutionary instead of the original variables to identify the weakly nonlinear coefficients. (0)
(n)
(n)
According to Sec. 2.5, the slow-varying evolutionary variables ai (t),bi (t) and ci (t) (n = 1, 2, . . .), together with the oscillation frequency ω̂i , can be computed from the given time history xi (t). Eq. (11) can then be reorganized such that known quantities are collected in the left hand side, while in the terms to be inferred are assembled in the right hand side.
14
ω̂i (0) äi (t) + ω̂i2 a0 (t) = −
2π
(n)
(n)
(n)
b̈i (t) + 2ω̂i nċi (t) − ω̂i2 (n2 − 1)bi (t) = − (n)
(n)
(n)
c̈i (t) − 2ω̂i nḃi (t) − ω̂i2 (n2 − 1)ci (t) = −
ω̂i π ω̂i π
Z t+π/ω̂i t−π/ω̂i Z t+π/ω̂i
fˆi∗ (x, ẋ, s)ds fˆi∗ (x, ẋ, s) cos(nω̂i s)ds
t−π/ω̂i
Z t+π/ω̂i
fˆi∗ (x, ẋ, s) sin(nω̂i s)ds
t−π/ω̂i
fˆi∗ (x, ẋ, s) = (ωi2 − ω̂i2 )x∗i (s) + ϵfi (x∗ (s), ẋ∗ (s)) − Fi (x∗ (s), ẋ∗ (s), s) n = 1, 2, · · · , N (18) where fˆi∗ (x, ẋ, t) = Θ∗ (x, ẋ, t)Ξ∗ are unknown functions expressed as a sparse combination of a set of linear and nonlinear candidate basis function. Thanks to this rearrangement, the system is expressed in a standard regression form suitable for sparse identification. Unlike conventional sparse regression, our objective is to find a sparse coefficient vector Ξ∗ that simultaneously satisfies all equations in Eq. (18), rather than a single one. This gives rise to the concept of joint sparse regression, where the goal is to estimate a shared sparse representation across multiple regression targets. Mathematically, given M = 2N + 1 equations reported in the following: yi (t) = Θi (t)Ξ∗ , (0) (0) yi (t) äi (t) + ω̂i2 a0 (t) (1,cos) (1) (1) yi (t) b̈i (t) + 2ω̂i ċi (t) (1,sin) (1) (1) yi (t) c̈i (t) − 2ω̂i ḃi (t) , yi (t) = (2,cos) = (2) (2) (2) yi (t) b̈i (t) + 4ω̂i ċi (t) − 3ω̂i2 bi (t) (2,sin) (2) (2) (2) 2 yi (t) c̈i (t) − 4ω̂i ḃi (t) − 3ω̂i ci (t) .. .. . .
15
(19)
(20)
(0) Θi (t)
t+π/ω̂i ω̂i − 2π Θ∗ (s)ds t−π/ω̂i ω̂i R t+π/ω̂i ∗ − π t−π/ω̂i Θ (s) cos(ω̂i s)ds
R
(1,cos) Θi (t) (1,sin) ω̂i R t+π/ω̂i ∗ Θi (t) − π t−π/ω̂i Θ (s) sin(ω̂i s)ds , Θi (t) = = ω̂i R t+π/ω̂i ∗ (2,cos) Θi (t) − π t−π/ω̂i Θ (s) cos(2ω̂i s)ds (2,sin) ω̂i R t+π/ω̂i ∗ Θi (t) − π t−π/ω̂i Θ (s) sin(2ω̂i s)ds .. .. . .
(21)
the objective is to determine a common coefficient vector Ξ∗ that minimizes the overall residual while maintaining sparsity. Only when this condition is met, the reconstructed system can accurately recover the original formulation expressed in Eq. (1). In Eqs. (20) and (35), the superscript (·)(n,t) identifies the harmonic term in Eq. (18). Specifically, n ∈ R+ represents the order of the harmonic frequency, corresponding to the n-fold multiple of the fundamental frequency ω̂i , and {cos, sin} indicates the cosine and sine components. 2.7. Joint sparse regression across multiple equations To infer the sparse coefficient vector Ξ∗ that best satisfies all evolutionary equations, we propose a joint sparse regression method across multiple equations all evolutionary equations Eq. (19), rather than fitting each equation independently. The problem is formulated as minimizing the total residual across the 2N +1 evolutionary equations derived from the weakly nonlinear oscillator, each corresponding to a harmonic component over N + 1 frequencies, with an L1 regularization term to ensure sparsity. A key algorithmic choice is setting which harmonics are considered in the evolutionary representation among yf ull (t) = [y (0) (t), y (1,cos) (t), y (1,sin) (t), y (2,cos) (t), y (2,sin) (t), · · · ]. Once set, goal of the optimization is to identify a sparse coefficient vector Ξ∗ that simultaneously satisfies all evolutionary equations derived from weakly nonlinear forcing oscillator as shown in the following:
! Ξ∗ = arg min Ξ∗
X
∥y (m) − Θ(m) Ξ∗ ∥22 + λ∥Ξ∗ ∥1
(22)
m∈IS
where: IS collects all the index of the effective harmonic components (e.g, 0, (1, cos), (1, sin), · · · ); y (m) and Θ(m) denote the response and library matrices for the m-th equation; λ > 0 con16
trols the sparsity level of the solution. We introduce a two-phase sparse regression method to solve Eq. (27). The first phase utilizes a convex optimization algorithm to select an initial subset of the library. In this work, we employ the CVXOPT solver (Vandenberghe, 2010) to handle the L1 regularization term. This step provides a preliminary sparse structure for the coefficient vector. The second phase involves a sparse fine-tuning process that removes terms with negligible contributions, thus ensuring a concise and generalized model. We apply an approach inspired by orthogonal matching pursuit to narrow down the model space. We select only the most relevant basis functions from the initial subset. The contribution cij of each basis function j for the i-th dimension is calculated using Eq. (23): 2 Z T X Θ̂mij (t)Ξ∗ ij cij = dt (yim (t))2 0 m∈IS
(23)
These values are then normalized to a range between 0 and 1. Finally, we fine-tune the model by applying a cut-off threshold. Any basis function with a contribution value lower than this threshold is removed from the library. This two-phase strategy ensures the identification of a parsimonious dynamical system. It maintains high predictive accuracy while preventing overfitting to numerical noise. 2.8. Parameter continuation method In nonlinear dynamical systems, periodic solutions often depend on one or more system parameters, such as the excitation frequency and intensity and can be computed using several classical numerical techniques, such as the harmonic balance method (Cheung et al., 1990), the shooting method (Osborne, 1969), and collocation-based approaches (Carrington Jr, 2021). These methods transform the governing differential equations into nonlinear algebraic systems by enforcing periodicity conditions, and the steady-state responses are then obtained through iterative solvers such as the Newton–Raphson method (Ypma, 1995). However, these traditional approaches are often computationally expensive because they require repeated integrations or nonlinear solves for each excitation parameter. In addition, they usually depend on initial guesses that are close to stable solutions, which makes them incapable of accurately tracing unstable branches or capturing turning points in the frequency–response 17
curves. As a result, their applicability is limited when a complete description of the system dynamics, including both stable and unstable periodic responses, is required. The continuation method provides a systematic procedure to trace the evolution of these solutions as parameters vary, allowing for the identification of critical phenomena such as bifurcations or stability transitions. Consider a general system of nonlinear algebraic equations derived from the steady or periodic conditions of a dynamical system:
G(u, ζ) = 0
(24)
where u denotes the unknown state vector (e.g., displacement and velocity), and ζ represents the control parameter to be varied (e.g., frequency and amplitude of external forcing). A straightforward approach would be to increment ζ in small steps and solve for u at each step using the previous solution as the initial guess ζk+1 = ζk + ∆ζ. However, this simple “parameter sweep” method fails e.g. near turning points or folds in the solution curve, where the slope ∂u/∂ζ becomes infinite and the solution path cannot be followed uniquely. To overcome this limitation, arc-length continuation (also known as pseudo–arc-length continuation) is commonly employed (Krauskopf et al., 2007). In this method, the solution path is parameterized by an artificial variables, representing the arc length along the curve in the (u, ζ) space. The continuation step is then constrained by: (∆u)T (∆u) + (∆ζ)2 = (∆s)2
(25)
where ∆s is the prescribed step size along the solution trajectory. By solving the extended system composed of Eq. (24) together with the arc-length constraint, the continuation algorithm can robustly track both stable and unstable branches of the solution family, even across fold and bifurcation points. This technique provides a powerful and general framework for the numerical exploration of nonlinear response characteristics. These methods are implemented in many ready-to-use packages like Auto07p (Doedel and Oldeman, 1998), that implements collocation methods in FORTRAN to perform numerical continuation and bifurcation analysis; Manlab, a Matlab tool that uses HB methods 18
and Asymptotic Numerical Method (Guillot et al., 2019); Nlvib, that also exploits HB methods (Krack and Gross, 2019), and many others among which we mention COCO (Dankowicz and Schilder, 2013). The numerical examples addressed in this work will be solved using the Matlab MATCONT package (Dhooge et al., 2006), that exploits collocation methods to perform the continuation of periodic orbits. 2.9. Mean chamfer distance of frequency response curves (MCDRC) Frequency response curves (FRCs) are the primary goals for one-shot learning and prediction, describing complex dynamic behaviors of weakly nonlinear oscillators. FRCs characterize a steady-state output quantity of interest, such as maximum deflection, across various load multipliers β for all excitation frequency values Ω within a specified range. To quantitatively assess the predictive ability, an error metric named Mean Chamfer Distance of frequency response curves (MCDRC) between two FRCs is defined herein. It is worth mentioning that FRCs are generally not single-valued functions of the excitation frequency. For a given frequency, multiple periodic solutions may exist due to nonlinear effects such as jump phenomena, hysteresis, and multi-stability. Furthermore, continuationbased computations generally produce differently spaced samples along solution branches, and distinct response branches may be specified using varying frequency grids. These properties violate the assumptions required by most pointwise or function-based error measures. To address this limitation, we represent FRCs as point sets in the frequency-response space. This allows to naturally support multi-valued responses, discontinuities, and nonuniform sampling, hence offering a consistent and stable foundation for quantitative comparison of nonlinear frequency response curves. We then exploit the Chamfer distance dch (C → Cref ) (Butt and Maragos, 1998) to establish a novel metric termed Mean Chamfer Distance of Frequency Response Curves (MCDRC). The chamfer distance dch (C → Cref ) is defined as the average distance from each point on the predicted curve C = {(xi , yi )}N i=1 to its nearest neighbor on the reference curve Cref = {(xj , yj )}M j=1 . It can be computed as: dch (C → Cref ) =
1 X min ∥p − q∥ N p∈C q∈Cref 19
(26)
By selecting the oscillation amplitude (A) and phase (P ) as primary steady-state descriptors, P N the FRC is discretized into two point sets: C A = {(Ωi , Ai )}N i=1 and C = {(Ωi , Pi )}i=1 . Here, A
and P denote the amplitude and the phase lag relative to the external excitation, respectively. The MCDRC is defined as follows: (27)
P A ) ) + dch (C P → Cref MCDRF = dch (C A → Cref
Prior to calculating the distance between points, a max-min scaling normalization preprocessing step is introduced to ensure that the scaling disparity between the axes does not spoil the representativeness of MCDRC. Figure 3 presents an illustrative example of MCDRC comparing two predicted FRCs against a reference FRC. The MCDRC for FRC1 is 0.24, indicating poor predictive accuracy, whereas the value for FRC2 is 0.01, representing excellent agreement with the reference. This comparison provides an intuitive demonstration of how MCDRC quantifies the similarity between FRCs. ) 5 &
) 5 &
3 K D V H P
$ P S O L W X G H A
) 5 & U H I H U H Q F H
) U H T X H Q F \
) U H T X H Q F \
Figure 3: Two examples of mean chamfer distance of frequency response functions (MCDRC): FRC1 (blue solid line), whose MCDRC is 0.24, exhibits a noticeable deviation from the reference (gray dashed line); FRC2 (red solid line) has a MCDRC of 0.01 shows a much better correlation with the reference.
20
3. Numerical Validations To demonstrate the effectiveness of the proposed framework, we first consider a benchmark one-dimensional (1D) forced oscillator characterized by linear dissipation, weak quadratic and cubic nonlinearities. This model is widely used in nonlinear dynamics to represent various physical systems (Kudryashov, 2021). The objective is to evaluate the framework capability in formulating the governing equations and to numerically verify the equivalence between the original formulation (Eq. (1)) and the generalized evolutionary formulation (Eq. (11)). The equation of motion is expressed as: ẍ + ω 2 x + cẋ + α1 x2 + α2 x3 = β cos(Ωt)
(28)
where x,ẋ and ẍ denote the displacement, velocity and acceleration of the oscillator, respectively; ω is the natural frequency, α1 , α2 ∈ R are parameters characterizing the weak nonlinearities, c ∈ R is the linear damping coefficient, β and Ω are the amplitude and frequency of the external forcing. To generate the training trajectory, the oscillator described by Eq. (28) is numerically integrated using the following model parameters: c = 1 × 10−2 , α1 = 1 × 10−2 , α2 = 1 × 10−4 , ω = 2, Ω = 1.999, β = 0.5, with initial states x0 = [0, 0]. Indeed, the magnitudes of α1 , α2 and c are orders of magnitude smaller than ω. The resulting time-domain and frequency-domain responses are illustrated in Figure 4. The frequency spectrum clearly indicates that the dominant dynamical response occurs near the fundamental frequency of approximately 2. Additionally, two secondary peaks are observed at the frequencies 0 and 4, which correspond to the DC component and second-order harmonic, respectively. These multi-frequency features underline the necessity of a higherorder analytical treatment. It is worth mentioning that the previously proposed EvLOWN approach (Ma et al., 2026) is unable to handle the presence of higher-order harmonics. Through this numerical example, we aim to demonstrate the efficacy of the proposed method regarding three key aspects: • First, to justify the necessity of incorporating multi-frequency responses in forced nonlinear oscillators; 21
x
Figure 4:
7 L P H V
) U H T X H Q F \ + ]
Benchmark 1D forced oscillator. Training data of numerical quadratic and cubic nonlinear
oscillator in time and frequency domain.
• Second, to validate the evaluation framework for multiple slowly varying evolutionary variables and the generalized harmonic balance (GHB) approach; • Third, to showcase the effectiveness of the method in both governing equation inference and prediction of complex nonlinear dynamical behaviors. 3.1. Multi-frequency response under external forcing The emergence of multi-frequency response stems from the system nonlinearities, specifically where the nonlinear restoring force exhibits a nonlinear dependence on the displacement x. As the amplitude of displacement increases, the contributions of higher-order nonlinear terms may grow order of times faster (e.g. x2 and x3 ), leading to the emergence of pronounced harmonic components in the response spectrum. For forced system, the external excitation may drive the system into a regime where higher-order nonlinearities dominate, resulting in a rich multi-frequency spectrum. Conversely, in the absence of external forcing, a free oscillator typically exhibits diminished vibration amplitudes, effectively oscillating at a single frequency (i.e., its natural frequency), as other harmonic components are negligible . To quantitatively illustrate this phenomenon, we perform a numerical test using the nonlinear oscillator ruled by Eq. (28). All parameters remain constant except for the amplitude of the external forcing β, which varies from 0 to 1. The additional harmonic components, specifically the DC offset (zero-order) and second-order harmonic responses, are evaluated under different values of β. The initial conditions are adjusted to x0 = [1, 0]. As shown in 22
Figure 5, when β approaches zero, the magnitudes of these extra harmonic components vanish, indicating a nearly quasi-linear single-frequency response. However, as β increases, these higher-order harmonics exhibit a steady growth. This trend confirms that intensified external excitation amplifies nonlinear effects, thereby underscoring the necessity of incorporating multi-frequency responses when analyzing forced nonlinear oscillators.
' & F R P S R Q H Q W 2 Q G R U G H U K D U P R Q L F F R P S R Q H Q W
, Q W H Q V L W \ R I G L I I H U H Q W R U G H U K D U P R Q L F F R P S R Q H Q W V
Figure 5:
$ P S O L W X G H R I H [ W H U Q D O I R U F L Q J
Benchmark 1D forced oscillator. Influence of external forcing amplitude β on the intensities of
the DC offset (zero-order) and second-order harmonic components.
3.2. Numerical validation of generalized harmonic balance We evaluate the theoretical foundation of the proposed method, specifically focusing on the validity of slow-varying evolutionary variables and the generalized harmonic balance formulation. The fundamental assumption is that the response x(t) of a weakly nonlinear forced oscillator can be accurately represented as a superposition of multiple harmonic components. Figure 6(a) presents a comparison between the ground-truth response x(t), obtained via numerical integration, and the approximated response x∗ (t) reconstructed using Eq. (8). Based on the spectral characteristics identified in Figure 4, only the zero-order and first-order harmonic components are selected for the oscillator in Eq. (28). We neglect the second-order harmonic in this specific case. It is important to note that including insignificant harmonic components can degrade the inference results. The slowly varying evolutionary variables a(0) (t), b(1) (t), and c(1) (t) are extracted following the procedure described in Sec. 2.5. As 23
illustrated in the zoomed-in view in Figure 6(b), the reconstructed signal x∗ (t) exhibits excellent agreement with the reference response x(t), thereby confirming the robustness of the proposed approximation.
D
E
7 L P H > V @ $ S S U R [ L P D W H x *
7 U X H x
( Y R O X W L R Q D U \ H T X D W L R Q V
x
x
F
y(0) 1
cos) y(1, 1
(0) 1 t
sin) y(1, 1
(1, cos) t 1
(1, sin) t 1
7 L P H > V @
7 L P H > V @
Figure 6: Benchmark 1D forced oscillator. Numerical validation of the Generalized Harmonic Balance (GHB) method: (a) Comparison between the ground-truth displacement x (numerical integration) and the reconstructed x∗ obtained via harmonic superposition, see Eq. (8); (b) zoomed-in view of the transient response highlighting the reconstruction accuracy; (c) validation of the evolutionary terms, comparing the analytical transformations with their numerical counterparts for various harmonic components.
Furthermore, the validation of the GHB-based reformulation (Eq. (18)) is presented in Figure 6(c), comparing the terms derived from the evolutionary variables, with those computed by numerically integrating the weak nonlinearities f (x, ẋ) over a single period. The close match verifies the accuracy of the GHB method and establishes the equivalence between the original physical formulation (Eq. (1)) and the transformed evolutionary representation(Eq. (11)). 3.3. Inferred equations The proposed method is applied to the benchmark 1D oscillator to discover the governing equations and reconstruct its complete dynamic portrait, including frequency response curves (FRCs). To capture potential nonlinear behaviors, a candidate library Θ is constructed, consisting of ten pure polynomial terms up to the fifth order and trigonometric functions 24
Table 1: Benchmark 1D forced oscillator. The "Reference" columns report the parameters used for the training ; in "MEv-SINDy" columns, we list the parameter values identified by our method.
Coefficients
Relevant term
Reference
MEv-SINDy
ω2
x
4.0000
4.0012
c
ẋ
1.0000 × 10−2
1.0495 × 10−2
α1
ẋ2
1.0000 × 10−2
1.0071 × 10−2
α2
x3
1.0000 × 10−4
9.6691 × 10−5
β
cos(Ωt)
5.0000 × 10−1
5.0748 × 10−1
(cos(Ωt) and sin(Ωt)). Notably, cross-coupling terms between displacement and velocity are omitted to maintain model parsimony. The algorithm is implemented with the following hyper-parameters: the effective harmonic orders are selected as [0, 1], the residual tolerance is set to 1 × 10−1 , and the lowest contribution ratio is 5 × 10−2 . A significant advantage of the proposed approach is its data efficiency: the exact equation structure is successfully identified using data from only a single excitation case (indicated by the red point in Figure 7). The estimated coefficients, reported in Table 1, exhibit high precision relative to the reference values. The proposed method successfully identifies the correct equation structure from the library by only one excitation case (red Point in Figure 7). Estimated coefficients are reported in Table 1. These results demonstrate that the proposed approach achieves accurate and parsimonious discovery of the weakly nonlinear dynamics. After the governing equation is identified, the system frequency response function (FRF) can be computed either by time integration or by a parameter continuation method. The FRF here describes the stationary oscillation amplitude of the state x as a function of the excitation frequency Ω. As shown in Figure 7, the predicted FRF closely matches the groundtruth results for various excitation frequencies and amplitudes. It is worth emphasizing, although the governing equation is inferred using only a single excitation condition, the entire frequency-response space under arbitrary forcing conditions can subsequently be predicted with high accuracy. 25
7 U D L Q L Q J 3 R L Q W
* U R X Q G 7 U X W K
D
3 U H G L F W L R Q
E
) R U F L Q J L Q W H Q V L W \
3 K D V H P > U D G @
$ P S O L W X G H A
) U H T X H Q F \
) U H T X H Q F \
Figure 7: Benchmark 1D forced oscillator. Comparison between predicted and reference frequency response functions under different excitation amplitudes. The proposed method accurately captures the hardening dynamical performance of oscillator.
4. Applications To highlight the capabilities of the proposed approach, we investigate high-dimensional structural dynamical systems exhibiting large-amplitude vibrations and geometric nonlinearities. Micro-Electro-Mechanical Systems (MEMS) are typically actuated near their resonance frequencies and consequently undergo substantial deformations. These effects are further intensified by the fact that MEMS are monolithic devices often encapsulated in near-vacuum environments, where energy dissipation is minimal. As a result, the system dynamics are strongly nonlinear and characterized by complex coupled behaviors, showing highly nonlinear dynamical features that are rarely observed at the macro scale. Therefore, MEMS provide a representative and valuable platform for applying the proposed approach to evaluate its capability in predicting dynamical performance through the inferred governing equations. In this section, we demonstrate the application of the proposed method to high-dimensional structural dynamical systems using synthetic data. The discussion begins with the formulation of the full-order model (FOM) based on the finite element discretization of the governing 26
equations. Subsequently, the system is reduced through the Proper Orthogonal Decomposition (POD) technique to obtain a reduced-order representation that retains the dominant dynamical features. To illustrate the procedure and evaluate the performance of the method, two case studies are presented: a clamped–clamped beam, and a more complex MEMS micromirror device that demonstrates the scalability of the proposed approach. 4.1. Formulation of full-order model The framework of structures subjected to large displacements and small strains is considered for full-order model simulations (Holzapfel, 2002). This is the operating range of most microsystems since they are often actuated at resonance and large aspect ratios allow reaching displacements within the linear elastic range of the material. In this framework, the Saint Venant-Kirchhoff constitutive model is the most appropriate choice, and is given by:
S = A : E,
1 E = (∇d + ∇T d + ∇T d · ∇d) 2
(29)
where S is the second Piola-Kirchhoff strain tensor, A is the fourth-order elasticity tensor and E is the Green-Lagrangian Strain tensor. Here we denote by d the displacement field and by ∇(·) the (material) gradient defined with respect to the reference configuration. The weak form of the linear momentum conservation law is:
Z
Z
Z
P[d] : ∇ wdΩ0 =
ρ0 d̈·wdΩ0 + Ω0
T
Ω0
Z f ·wdS0 ,
ρ0 F·wdΩ0 + Ω0
∀w ∈ H01 (Ω0 ) (30)
ST 0
where the integrals are expressed in the reference configuration Ω0 . Here, ρ0 denotes the material density, P[d] = (1 + ∇d) · S is the first Piola-Kirchhoff stress tensor, F is the body force per unit mass, f is the surface tractions prescribed on the surface ST 0 and w is the test velocity selected in H01 (Ω0 ), i.e. the space of functions with finite energy that vanish on the portion SU ⊂ ∂Ω0 where Dirichlet boundary conditions are prescribed. By applying spatial discretization to Eq. (30), for instance through the finite element method, and incorporating a Rayleigh damping model, the system can be transformed into
27
a set of coupled nonlinear ordinary differential equations of the following form: MD̈ + CḊ + KD + G(D, D) + H(D, D, D) = B(D, β, Ω, t),
t ∈ (0, T )
(31)
where the vector D ∈ Rk collects all the unknown displacement nodal values, M ∈ Rn×n is the mass matrix, C = (ω0 /Q)M is the Rayleigh model mass proportional damping matrix, defined with respect to a reference eigenfrequency ω0 and a quality factor Q. B ∈ Rn is the nodal force vector which depends on the actuation intensity parameters β, the angular frequency of actuation Ω and in general also on D. The internal force vector has been exactly decomposed in linear, quadratic, and cubic power terms of the displacement: K ∈ Rn×n is the linear stiffness matrix, while G ∈ Rn and H ∈ Rn are vectors given by monomials of second and third order, respectively. The components of these vectors can be expressed using P P an indicial notation: Gi = nj,k=1 Gijk Dj Dk , Hi = nj,k,l=1 Hijkl Dj Dk Dl , i = 1, · · · , n. Eq. (31) represents the high-fidelity full-order model (FOM) which depends on the input parameters β and Ω. It should be recalled that in resonating MEMS FRCs are an important design tool in which a selected quantity, like the midspan deflection of a beam or the rotation of a micromirror, is plotted versus Ω for different β. Indeed, resonators operate close to a reference frequency where the behavior should be predictable. However, approximating and reconstructing different FRCs for high dimensional dynamical system is computationally demanding, making the use of full-order models for this task largely impractical 4.2. Proper Orthogonal Decomposition To substantially accelerate the computationally expensive simulations and obtain a more interpretable and parsimonious set of governing equations, we propose a reduction technique that simultaneously decreases the problem dimensionality and introduces a new set of coordinates better adapted to capture the system dynamics within only a few dominant modes. The proposed reduction technique is based on Proper Orthogonal Decomposition (POD), a model reduction method that constructs an ordered set of basis functions according to their energy content. By projecting the full-order data onto a limited number of dominant high-energy modes, a reduced-order representation of the system can be efficiently obtained (Rowley et al., 28
2004). The snapshots of the FOM solutions are collected to generate a matrix X ∈ Rp×k , whose p denotes the time stamps of solutions. Next, the Singular Value Decomposition (SVD) of the matrix X is computed: X = UΣVT , where the columns of the orthonormal matrix U ∈ Rk×k are left singular vectors, often called Proper Orthogonal Modes (POMs) (Abbott and Kepler, 2005; Lu et al., 2019) or spatial modes, the columns of V ∈ Rp×p give the temporally evolving coefficients, and Σ ∈ Rk×p is a diagonal matrix representing the singular values of the matrix X ordered from the largest to the smallest conventionally. In particular, the rank of X is equal to the number of nonzero singular values, and the optimal rank-k̂ approximation X̂ of X is given by the P rank-k̂ SVD truncation X̂ = k̂i=1 σi Ui ViT , in which σi is the ith singular value contained in the diagonal of Σ, and Ui , Vi are the ith column of U and V, respectively. In this work, k̂ = 3, which gives a reduced-order model of dimension 3. The reduced order dynamical coordinates Λi = σi ViT ,
i = 1, 2, 3.
4.3. High-dimensional case 1: Beam MEMS resonator The doubly clamped beam is investigated as a representative example of a high-dimensional nonlinear structural system, where a straight beam-type MEMS resonator is excited near its resonance frequency. This configuration captures the essential dynamics of micro-resonators (Zega et al., 2020) in which geometric nonlinearities play a significant role in shaping the system’s response. The doubly clamped beam in this case has length L = 1000 µm with rectangular crosssection of dimensions 10µm × 24µm, as shown in Figure 8a. The beam is fabricated from isotropic polysilicon (Corigliano et al., 2004), with density ρ = 2330 kg/m3 , Young modulus E = 167 GPa and Poisson coefficient ν = 0.22. Dirichlet boundary conditions are enforced on the two opposite sides of the beam. The quality factor is assumed equal to Q = 50. The device is excited by a body load F = Mϕ1 β cos(Ωt) proportional to the first eigenmode ϕ1 , with M mass matrix and β load multiplier. The device vibrates according to its first bending mode at ω0 = 0.5475rad/µs. In order to keep the FOM computational time at a reasonable level, a rather coarse mesh with 2607 nodes has been employed. To effectively 29
infer the underlying dynamics of beam MEMS resonator, the POD modes are built retaining in the linear trial space the first three most energetic POD modes, depicted in Figure 8b-d, respectively. This number of POD modes represents the minimum required to achieve a good accuracy while maintaining a computationally efficient model (Conti et al., 2023).
Figure 8: Schematic representation of the beam MEMS resonator. (a) Geometry and mesh of doubly clamped beam MEMS resonator, (b-d) First three most energetic POD modes obtained from SVD
Figure 9: Beam MEMS resonator training data in time, frequency and evolutionary domain. The model is trained using a single parameter set (β = 0.25, Ω = 0.551). (a) Time series of the POD modal coordinates; (b) Frequency spectra, where red boxes indicate the identified effective harmonic orders; (c) Comparison of slowly varying variables for the first three POD modes, showing excellent agreement between the ground truth (solid) and prediction by Eq.(32)(dashed).
To infer the underlying dynamics of the MEMS resonator, training data are collected from a single simulation case with a load multiplier β = 0.25 and excitation frequency Ω = 0.551. Figure 9 illustrates the characteristics of the training set in the frequency, time 30
and evolutionary domains. Specifically, Figure 9a displays the time-history responses of the first three POD modes, while Figure 9b shows their corresponding frequency spectra. The red boxes in the spectra highlight the effective harmonic orders, confirming that the system response is predominantly composed of multiple discrete frequencies. Using the proposed filtering procedure (Sec 2.4), the slow varying evolutionary variables for the POD modes are extracted, as shown in Figure 9c. Building upon these extracted variables, our approach successfully identifies a parsimonious set of governing equations that govern the reducedorder space from a candidate library Θ consisting of ten pure polynomial terms up to the fifth order and trigonometric functions (cos(Ωt) and sin(Ωt)), as explicitly formulated in Eq. (32). The identified model is capable of accurately reconstructing the slow-varying evolutionary variables (Figure 9c dashed line), confirming that the identified nonlinear terms correctly capture the weakly nonlinear effects across different harmonic orders. Specifically, these identified equations capture the linear restorative forces alongside the geometric nonlinearities inherent in the large-deflection regime of the MEMS beam.
= 1.902β cos(Ωt) ü +2.998×10−1 u1 +2.730 × 10−6 u31 + 1.100 × 10−2 u̇1 1 ü2 +3.006×10−1 u2 −1.468×10−6 u1 +5.453×10−9 u31 +2.450×10−5 u̇1 = 4.202×10−3 β cos(Ωt) ü +3.006×10−1 u +2.402×10−6 u2 − 1.574×10−5 u̇2 =0 3
3
1
1
(32) The resulting model provides a compact yet precise representation of the resonator complex dynamics, enabling efficient prediction of its complete dynamic portrait while significantly reducing the computational burden compared to the full-order finite element model. Specifically, to evaluate the generalization capability of the identified equations, we predict the global FRCs of the beam MEMS resonator under various excitation intensities. As illustrated in Figure 10, the predicted FRCs (dashed line) demonstrate remarkable agreement with the reference FOM results (solid line) across entire investigated range of forcing amplitudes (β ∈ [0.125, 0.750]µN), with the left and right panels presenting the amplitude-frequency and frequency-phase curves, respectively. In particular, the identified model successfully captures the nonlinear hardening behavior, characterized by the rightward tilting of the resonance 31
7 U D L Q L Q J 3 R L Q W
* U R X Q G 7 U X W K
1
xmid > m @
1
1
) U H T X H Q F \ > U D G s @
1
1
3 K D V H > U D G @
) R U F L Q J L Q W H Q V L W \
1
) U H T X H Q F \ > U D G s @
3 U H G L F W L R Q
Figure 10: FRCs of beam MEMS resonator. Comparison between predicted and reference frequency response functions under different excitation amplitudes (left: amplitude; right: phase). The proposed method accurately captures the hardening dynamical performance of beam MEMS resonator.
peaks as the excitation amplitude increases. Notably, while the governing equations were inferred using data from only a single training point (the red dot at β = 0.25, Ω = 0.551), the model accurately extrapolates the complex bifurcation characteristics and phase transitions for forcing levels significantly higher than the training condition. This high-fidelity match underscores the physical consistency and robustness of the proposed identification framework for high-dimensional nonlinear structural systems. 4.4. High-dimensional case 2: MEMS micromirror Scanning micromirrors have experienced rapid growth in recent years, driven by their successful implementation in a wide range of applications, from pico-projectors for Augmented Reality (AR) displays to three-dimensional (3D) scanning systems for Light Detection and Ranging (LiDAR) technologies. Because of the inertial and geometrical effects triggered by large rotations, micromirrors are intrinsically nonlinear and the prediction of their dynamic behavior is essential to guarantee a proper design and control of the mirror during the online stage. Specifically, maintaining motion stability is critical for line scanning performance, making the correct characterization of hardening and softening effects a priority. 32
Figure 11: Schematic representation of the MEMS micromirror: (a) Optical picture of the Micromirror; (b) Schematic view of the layout with few details on the device components; (c) Third (torsional) eigenmode that is actuated during operations and the boundary conditions used in the simulation of half of the device; (d-f) First three most energetic POD modes obtained from SVD.
The mirror considered, fabricated by STMicroelectronics, is illustrated in Figure 11. The mirror plate is suspended to a gimbal connected with a torsional beam along the rotation axis and two suspension beams on each side. The mirror is assumed to be made of isotropic polysilicon (Corigliano et al., 2004), with density ρ = 2330kg/m3 , Young modulus E = 167GP a and Poisson coefficient ν = 0.22. Exploiting structural symmetry„ only half of the micromirror is modeled with the FEM and a total of 9723 dofs. The Dirichelet boundary conditions are imposed in the substrate (the dark green areas in Figure 11b) and on the symmetry plane to enforce the symmetry conditions. Given the high stiffness of the central mirror plate, its rotation angle is adopted as the primary state variable for both the timedomain response and the frequency response function (FRF) analysis. To characterize the complex nonlinearities of the scanning micromirror, the training set of the MEv-SINDy framework is derived from a single simulation trajectory with a load multiplier β = 3 and an excitation frequency Ω = 0.183841. Figure 12 illustrates the sampled data across the time, frequency and evolutionary domains. Unlike the beam case, the torsional tilting of the mirror triggers a richer multi-harmonic spectrum, as evidenced by the discrete frequency peaks in Figure 12b. The filtering procedure effectively isolates the evolutionary variables for the underlying manifold (Figure 12c), providing a clean foundation for
33
Figure 12:
MEMS micromirror training data in time, frequency and evolutionary domains: The model is
trained using a single parameter set (β = 3, Ω = 0.183841). (a) Time series of the POD modal coordinates; (b) Frequency spectra, where red boxes indicate the identified effective harmonic orders; (c) Comparison of slowly varying variables for the first three POD modes, showing excellent agreement between the ground truth (solid) and prediction by Eq. (33)(dashed).
sparse regression. The algorithm yields a compact set of reduced-order equations (Eq. (33)) that faithfully reproduce the training signals (see the comparison between continuous and dashed lines in Figure 12c). Notably, the identified model successfully decouples the inertial nonlinearities and cross-modal interactions inherent in the gimbal-suspension assembly.
ü + 3.383 × 10−2 u1 − 1.180 × 10−11 u31 + 1.865 × 10−4 u̇1 = 1.076 × 10−1 β cos(Ωt) 1 ü2 + 3.374 × 10−2 u2 − 1.103 × 10−7 u21 + 3.912 × 10−5 u̇21 = 0 ü + 3.374 × 10−2 u − 9.127 × 10−7 u2 + 2.392 × 10−5 u̇2 = 0 3 3 1 1 (33) The robustness of the identified model is further validated by predicting the global Frequency Response Functions (FRCs). As shown in Figure 13, the MEv-SINDy framework accurately extrapolates the system behavior across a wide range of forcing intensities. Of particular interest is the model ability to capture the softening behaviors, which is a critical factor for ensuring scanning stability in LiDAR and AR applications. The seamless agreement between the ROM predictions and the FOM ground truth—despite the model being trained on a single localized point—underscores the physical consistency of the discovered 34
dynamics.
Figure 13: FRCs of MEMS micromirror. Comparison between predicted and reference frequency response functions under different excitation amplitudes (left: amplitude; right: phase). The proposed method accurately captures the softening dynamical performance of MEMS micromirror.
5. Robustness experiments While the preceding cases demonstrate the high fidelity of the MEv-SINDy framework in identifying complex structural dynamics, a systematic assessment of its robustness is essential to define its practical operational boundaries. This section investigates the sensitivity and reliability of the proposed method under varying data conditions, focusing on two critical aspects. First, we examine the influence of training point selection on global predictive accuracy, specifically investigating whether the location of the training sample affects the model ability to extrapolate the system global response. Second, we evaluate the scaling of prediction precision with the number of training points. By comparing the performance of models trained on sparse versus enriched datasets, we highlight our approach primary advantage: the ability to achieve high-dimensional ROM discovery with minimal data requirements, thereby significantly reducing the computational overhead inherent in traditional data-driven modeling. 35
5.1. Influence of training point selection To provide a comprehensive assessment of MEv-SINDy’s reliability, we conduct extensive testing across two benchmark datasets: the beam MEMS resonator (56 distinct operating points, featuring specific forcing frequency and amplitude) and the MEMS micromirror (22 distinct operating points). For each individual point in these datasets, the governing equations are first identified via the proposed framework. Subsequently, the full-range frequency responses are reconstructed using MATCONT (as detailed in Sec. 2.7), a sophisticated parameter continuation solver. To quantify the fidelity of the identified models, we utilize the Mean Chamfer Distance of Frequency Response Curves (MCDRC) as our primary metric (as detailed in Sec. 2.8). Specifically, a score below MCDRC < 0.1 indicates that the identified model has successfully captured the essential nonlinear characteristics across the entire dynamical regime.
Figure 14: Robustness evaluation of proposed approach across various training conditions. The predictive performance is assessed by identifying the model using a single training point and evaluating its global extrapolation capability. Blue circles indicate a successful inference (MCDRF < 0.1), while red crosses denote a failed inference. (a) Beam MEMS resonator: results for 56 different training points, covering two excitation amplitudes β ∈ {0.125, 0.250} and 28 frequencies Ω ranging from 0.526 to 0.564 rad/µs. A total of 51 training points yield successful global predictions. (b) MEMS micromirror: results for 22 distinct training points, covering two excitation amplitudes β ∈ {0.5, 3.0} and 11 frequencies Ω ranging from 0.1837 to 0.1841 rad/µs. All 22 training points achieve successful inference.
In particular, the reliability of the MEv-SINDy framework is systematically evaluated by 36
analyzing the influence of training point selection on global predictive accuracy. As illustrated in Figure 13, the framework exhibits remarkable robustness across diverse operating regimes. For the beam MEMS resonator, 51 out of 56 candidate training points yield successful global inferences (Figure 14a), while for the MEMS micromirror, the framework achieves a 100% success rate across all 22 tested points (Figure 14b). Notably, accurate global ROMs are identified when the training data is sampled from the strongly nonlinear regimes characterized by significant resonance tilting.
Figure 15: Quantitative and qualitative error analysis of FRF predictions. (a) MCDRF distribution across all successful inferences; violin plots show high error concentration near zero for both devices. (b–c) Qualitative comparison between Ground Truth (dashed) and predictions (solid) for the beam resonator, and (d–e) for the micromirror. The consistent overlap across various training conditions highlights the framework’s minimal sensitivity to training data selection.
This consistency is further quantified in Figure 15a, where the MCDRC distribution for all successful cases is shown to be highly concentrated near zero. The narrow spread of these violin plots underscores that the identification fidelity is largely invariant to the specific choice of training conditions. Qualitatively, Figures 15b–e confirm this high fidelity; the FRC curves 37
reconstructed from various localized training points demonstrate a nearly perfect overlap with the ground-truth simulations. The framework faithfully reproduces complex nonlinear phenomena, including the softening behavior of the beam and the dual hardening-softening transitions in the micromirror. Collectively, these results demonstrate that MEv-SINDy effectively extracts the underlying physical operators rather than merely overfitting localized data, enabling reliable global extrapolation from minimal, localized measurements. To further enhance the reliability of the MEv-SINDy framework within a one-shot learning context, it is imperative to investigate the specific conditions under which identification fails. The five failure cases observed in the beam MEMS resonator (Figure 15a) can be categorized into two distinct modes: 1. Limited nonlinear effects: Several failed points are located at frequencies far from the system natural frequency. In these near-linear regimes, the nonlinear contributions to the dynamical response are extremely small. Consequently, the sparse regression algorithm only discovers the linear damping for the first POD mode. While the resulting model can predict a precise dynamical response within the local linear regime, it lacks the cubic stiffness nonlinearity required for global extrapolation. This underscores that for robust one-shot learning, the training point must be selected from a region where nonlinear effects are sufficiently prominent to be accounted by the sparse identification process. 2. Incorrect identified dynamics: Occasionally, identification failure occurs because the discovered model converges to an incorrect dynamical equation, leading to an erroneous global FRF prediction. Fortunately, such discrepancies are readily detectable during the validation phase; by comparing the model’s predicted time-history against the original sampled trajectory, one can easily identify mismatches in time domain. To rectify this, two primary strategies can be employed: fine-tuning the hyperparameters (such as the sparsity threshold) to re-infer the true underlying dynamics or re-selecting a more informative training point that better characterizes the system’s nonlinear landscape. 38
5.2. Multi-point training and frequency normalization While the MEv-SINDy framework demonstrates high reliability in one-shot learning, incorporating data from multiple excitation conditions can further stabilize the identification process. However, a fundamental mathematical discrepancy arises when attempting to directly merge datasets from different operating points. In the physical domain, the underlying system parameters (Eq. (1)) are intrinsic constants across different excitation cases. The linear frequency term requires special consideration. Its physical coefficient ω02 is split into two distinct parts. The first part is the observed frequency ω̂ obtained via FFT (Eq. (13)). The second part is a small fine-tuning coefficient identified through the regression process (Eq. (19)). So the governing equations in our identification framework (Eq. (19)) are constructed based on the observed response frequency ω̂, which is determined via FFT for each specific excitation case. Because ω̂ varies across different forcing conditions, the mathematical representation of the linear frequency in Θ(t) becomes coupled with the local operating point. Consequently, Eq. (19) yields inconsistent coefficients for the linear frequency term across datasets, preventing a straightforward least-squares concatenation. To resolve this, we introduce a frequency-normalization step after reformulate the evolutionary library. By selecting a reference frequency ω̄ (e.g., the mean of all ω̂), we apply a compensation term to the target vector y(t). This ensures the linear frequency term is mapped to a consistent basis, effectively decoupling the physical parameters from the local observation frequency. The corrected target vector ycorr (t) replaces the original vector y in Eq. (19). This substitution ensures that the regression targets are consistent across all training sets. ycorr (t) = y(t) − (ω̂ 2 − ω̄ 2 )Θx (t)
39
(34)
where Θx stands for the candidate function x in the evolutionary library: R t+π/ω̂i (0) ω̂i − 2π Θx (t) −a (t) x(s)ds t−π/ω̂i 0 (1,cos) ω̂i R t+π/ω̂i Θx (t) − π t−π/ω̂i x(s) cos(ω̂i s)ds −b1 (t) (1,sin) ω̂i R t+π/ω̂i Θx (t) − π t−π/ω̂i x(s) sin(ω̂i s)ds −c1 (t) , = Θx (t) = (2,cos) = ω̂i R t+π/ω̂i Θx (t) x(s) cos(2ω̂i s)ds −b (t) − π Rt−π/ω̂i 2 (2,sin) ω̂i t+π/ω̂i Θx (t) − π t−π/ω̂i x(s) sin(2ω̂i s)ds −c2 (t) .. .. .. . . .
(35)
The efficacy of this multi-point approach is demonstrated in Figure 16. A comparison of the single-point training results reveals a trade-off in predictive performance: the model trained on Point 1 provides a better fit for low-amplitude excitation cases, whereas the model trained on Point 2 exhibits higher fidelity for high-amplitude regimes. This suggests that localized data, while accurate within its own regime, may not fully capture the system sensitivity across the entire forcing range. By implementing the frequency-normalization terms defined in Eq. (34), we can combine the complementary advantages of both points. The composite model (Training on Points 1 and 2) effectively fuses the “low-amplitude” and “high-amplitude” information into a single, unified ROM. As a result, the identified model achieves a superior and more consistent match with the ground-truth FRC under all tested conditions. This shows that our framework can successfully integrate disparate datasets to overcome the limitations of localized sampling, ensuring stable and high-fidelity global predictions. 6. Conclusion Learning the complete dynamic portrait of weakly forced oscillators typically requires large training datasets. This is true especially when increasing popular deep learning are applied. In this study, we proposed a one-shot method to learn frequency-response curves from a single excitation configuration or, in other words, from a single traning point. To infer weakly nonlinear effects, we modified the Sparse Identification of Nonlinear Dynamic (SINDy) algorithm accounting for an evolution-based learning framework (MEv-SINDy). 40
Figure 16: Comparison of Frequency Response Curves (FRCs) identified from different training sets. (left) Amplitude-frequency response; the forcing amplitudes increase from bottom to top. (right) Phase-frequency response; the forcing amplitudes decrease from bottom to top.
We shifted the focus from autonomous single-frequency dynamics to non-autonomous multifrequency dynamics. This advancement was achieved through the Generalized Harmonic Balance (GHB) method. The core of our approach is the identification of explicit governing equations. Once these equations are learned, they provide strong physical extrapolation capabilities. This analytical representation allows the model to predict complex nonlinear responses far beyond the training point. We can then use standard numerical tools to perform rapid parametric sweeps. This strategy effectively bridges the gap between limited data observations and global nonlinear response prediction, addressing the major challenge of data acquisition in real-world MEMS applications. We applied the proposed MEv-SINDy framework considering two high-dimensional MEMS applications. These cases included a beam resonator and a MEMS micromirror. In both examples, our method achieved highly accurate one-shot learning by predicting complex nonlinear responses from only a single training point. To ensure the reliability of the approach, we performed extensive robustness tests. We investigated the influence of training point selection
41
on the identification accuracy. This analysis helped us define effective strategies for choosing optimal data points. Furthermore, we extended the framework to include multi-point training thanks to a frequency normalization preprocessing. The results demonstrate the use of multiple training points can refine the MEv-SINDy outcomes across different excitation levels and frequency ranges. These findings confirm the practical utility of our method for real-world engineering design. The framework significantly reduces the computational burden while preserving the complex physics of forced oscillators.
Acknowledgments: L.R. is supported by the Joint Research Platform “Sensor sysTEms and Advanced Materials” (STEAM) between Politecnico di Milano and STMicroelectronics. T.M, W.C., and L.Z. are supported by the National Key Research and Development Program of China (grant no. 2022YFC3005301), the National Natural Science Foundation of China (grant no. 52478552 and no. 52378527), the Natural Science Foundation of Shanghai (grant no. 23ZR1464900), the Fundamental Research Funds for the Central Universities (grant no. 22120240363) and China Scholarship Council (grant no. 20230626015) Contributions: Conceptualization: T.M., L.R., A.F; Methodology: T.M.; Experiment: T.M., L.R., A.F.; Visualization: T.M.; Funding acquisition: W.C., L.Z.; Results Discussion: T.M., L.R, W.C., L.Z., A.F. Project administration: A,F.; Supervision: A.F.; Writing – original draft: T.M.; Writing – review & editing: T.M., L.R., A.F. Competing interests: Authors declare that they have no competing interests. Code availability: All data are provided in the main text or the supplementary materials. The codes are available in the public GitHub repository: https://github.com/ TengMa25/MEv-SINDy.git References Abbott, L.F., Kepler, T.B., 2005. Model neurons: from hodgkin-huxley to hopfield, in: Statistical Mechanics of Neural Networks: Proceedings of the Xlth Sitges Conference Sitges, Barcelona, Spain, 3–7 June 1990, Springer. pp. 5–18. 42
Amsallem, D., Zahr, M.J., Farhat, C., 2012. Nonlinear model order reduction based on local reduced-order bases. International Journal for Numerical Methods in Engineering 92, 891–916. Benner, P., Ohlberger, M., Cohen, A., Willcox, K., 2017. Model reduction and approximation: theory and algorithms. SIAM. Brunton, S.L., Proctor, J.L., Kutz, J.N., 2016. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc. Natl. Acad. Sci. 113, 3932–3937. Butt, M.A., Maragos, P., 1998. Optimum design of chamfer distance transforms. IEEE Transactions on Image Processing 7, 1477–1484. Carrington Jr, T., 2021. Using collocation to study the vibrational dynamics of molecules. Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 248, 119158. Chen, Z., Liu, Y., Sun, H., 2021. Physics-informed learning of governing equations from scarce data. Nat. Commun. 12, 6136. Cheung, Y., Chen, S., Lau, S., 1990. Application of the incremental harmonic balance method to cubic non-linearity systems. Journal of Sound and Vibration 140, 273–286. Conti, P., Gobat, G., Fresca, S., Manzoni, A., Frangi, A., 2023. Reduced order modeling of parametrized systems through autoencoders and sindy approach: continuation of periodic solutions. Computer Methods in Applied Mechanics and Engineering 411, 116072. Corigliano, A., De Masi, B., Frangi, A., Comi, C., Villa, A., Marchi, M., 2004. Mechanical characterization of polysilicon through on-chip tensile tests. Journal of Microelectromechanical Systems 13, 200–219. Dankowicz, H., Schilder, F., 2013. Recipes for continuation. SIAM. Darcy, M., Hamzi, B., Livieri, G., Owhadi, H., Tavallali, P., 2023. One-shot learning of stochastic differential equations with data adapted kernels. Physica D: Nonlinear Phenomena 444, 133583. 43
Dhooge, A., Govaerts, W., Kuznetsov, Y.A., Mestrom, W., Riet, A., Sautois, B., 2006. Matcont and cl matcont: Continuation toolboxes in matlab. Universiteit Gent, Belgium and Utrecht University, The Netherlands . Doedel, E.J., Oldeman, B., 1998. Auto-07p: continuation and bifurcation software. Montreal, QC: Concordia University Canada . Fei-Fei, L., Fergus, R., Perona, P., 2006. One-shot learning of object categories. IEEE Transactions on Pattern Analysis and Machine Intelligence 28, 594–611. doi:10.1109/ TPAMI.2006.79. Franco, N., Manzoni, A., Zunino, P., 2023. A deep learning approach to reduced order modelling of parameter dependent partial differential equations. Mathematics of Computation 92, 483–524. Fresca, S., Gobat, G., Fedeli, P., Frangi, A., Manzoni, A., 2022. Deep learning-based reduced order models for the real-time simulation of the nonlinear dynamics of microstructures. International Journal for Numerical Methods in Engineering 123, 4749–4777. Fresca, S., Manzoni, A., 2022. Pod-dl-rom: Enhancing deep learning-based reduced order models for nonlinear parametrized pdes by proper orthogonal decomposition. Computer Methods in Applied Mechanics and Engineering 388, 114181. Fung, L., Fasel, U., Juniper, M., 2025. Rapid bayesian identification of sparse nonlinear dynamics from scarce and noisy data, in: Proceedings A, The Royal Society. p. 20240200. Gao, T.T., Barzel, B., Yan, G., 2024. Learning interpretable dynamics of stochastic complex systems from experimental data. Nat. Commun. 15, 6029. Gao, T.T., Yan, G., 2022. Autonomous inference of complex network dynamics from incomplete and noisy data. Nat. Comput. Sci. 2, 160–168. Gobat, G., Opreni, A., Fresca, S., Manzoni, A., Frangi, A., 2022. Reduced order modeling of
44
nonlinear microstructures through proper orthogonal decomposition. Mechanical Systems and Signal Processing 171, 108864. Gonzalez, F.J., Balajewicz, M., 2018. Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems. arXiv preprint arXiv:1808.01346 . Guillot, L., Cochelin, B., Vergez, C., 2019. A taylor series-based continuation method for solutions of dynamical systems. Nonlinear Dynamics 98, 2827–2845. Hirsh, S.M., Barajas-Solano, D.A., Kutz, J.N., 2022. Sparsifying priors for bayesian uncertainty quantification in model discovery. Royal Society open science 9, 211823. Hoffmann, M., Fröhner, C., Noé, F., 2019. Reactive sindy: Discovering governing reactions from concentration data. J. Chem. Phys. 150. Holzapfel, G.A., 2002. Nonlinear solid mechanics: a continuum approach for engineering science. Huang, D., Abdel-Khalik, H., Rabiti, C., Gleicher, F., 2017. Dimensionality reducibility for multi-physics reduced order modeling. Annals of Nuclear Energy 110, 526–540. Ivanov, A., Iben, U., Golovkina, A., 2020. Physics-based polynomial neural networks for one-shot learning of dynamical systems from one or a few samples. ArXiv abs/2005.11699. URL: https://api.semanticscholar.org/CorpusID:218870280. Jiao, A., He, H., Ranade, R., et al., 2025. One-shot learning for solution operators of partial differential equations. Nature Communications 16, 8386. URL: https://doi.org/10. 1038/s41467-025-63076-z, doi:10.1038/s41467-025-63076-z. Krack, M., Gross, J., 2019. Harmonic balance for nonlinear vibration problems . Krauskopf, B., Osinga, H.M., Galán-Vioque, J., et al., 2007. Numerical continuation methods for dynamical systems. volume 2. Springer.
45
Kudryashov, N.A., 2021. The generalized duffing oscillator. Communications in Nonlinear Science and Numerical Simulation 93, 105526. Le Clainche, S., Vega, J.M., 2017. Higher order dynamic mode decomposition. SIAM Journal on Applied Dynamical Systems 16, 882–925. Loiseau, J.C., Noack, B.R., Brunton, S.L., 2018. Sparse reduced-order modelling: sensorbased dynamics to full-state estimation. J. Fluid Mech. 844, 459–490. Lu, K., Jin, Y., Chen, Y., Yang, Y., Hou, L., Zhang, Z., Li, Z., Fu, C., 2019. Review for order reduction based on proper orthogonal decomposition and outlooks of applications in mechanical systems. Mechanical Systems and Signal Processing 123, 264–297. Luo, A.C., Huang, J., 2012. Approximate solutions of periodic motions in nonlinear systems via a generalized harmonic balance. Journal of Vibration and Control 18, 1661–1674. Ma, T., Cui, W., Gao, T., Liu, S., Zhao, L., Ge, Y., 2023. Data-based autonomously discovering method for nonlinear aerodynamic force of quasi-flat plate. Physics of Fluids 35. Ma, T., Gao, T.T., Cui, W., Frangi, A., Yan, G., Zhao, L., 2026. Encoding cumulation to learn perturbative nonlinear oscillatory dynamics. Advanced Science n/a, e19707. doi:https://doi.org/10.1002/advs.202519707. Maday, Y., Rønquist, E.M., 2002. A reduced-basis element method. Journal of scientific computing 17, 447–459. Mangan, N.M., Brunton, S.L., Proctor, J.L., Kutz, J.N., 2016. Inferring biological networks by sparse identification of nonlinear dynamics. IEEE Trans. Mol. Biol. Multi-Scale Commun. 2, 52–63. Nayfeh, A.H., Balachandran, B., 2008. Applied nonlinear dynamics: analytical, computational, and experimental methods. John Wiley & Sons.
46
Niven, R.K., Cordier, L., Mohammad-Djafari, A., Abel, M., Quade, M., 2024. Dynamical system identification, model selection, and model uncertainty quantification by bayesian inference. Chaos: An Interdisciplinary Journal of Nonlinear Science 34. Osborne, M.R., 1969. On shooting methods for boundary value problems. Journal of mathematical analysis and applications 27, 417–433. Pagani, S., Manzoni, A., Quarteroni, A., 2018. Numerical approximation of parametrized problems in cardiac electrophysiology by a local reduced basis method. Computer Methods in Applied Mechanics and Engineering 340, 530–558. Quarteroni, A., Manzoni, A., Negri, F., 2015. Reduced basis methods for partial differential equations: an introduction. volume 92. Springer. Rao, C., Ren, P., Wang, Q., Buyukozturk, O., Sun, H., Liu, Y., 2023. Encoding physics to learn reaction–diffusion processes. Nat. Mach. Intell. 5, 765–779. Romor, F., Stabile, G., Rozza, G., 2023. Non-linear manifold reduced-order models with convolutional autoencoders and reduced over-collocation method. Journal of Scientific Computing 94, 74. Rowley, C.W., Colonius, T., Murray, R.M., 2004. Model reduction for compressible flows using pod and galerkin projection. Physica D: Nonlinear Phenomena 189, 115–129. Ruan, K., Xu, Y., Gao, Z.F., Liu, Y., Guo, Y., Wen, J.R., Sun, H., 2025. Discovering physical laws with parallel symbolic enumeration. Nature Computational Science , 1–14. Rudy, S.H., Brunton, S.L., Proctor, J.L., Kutz, J.N., 2017. Data-driven discovery of partial differential equations. Sci. Adv. 3, e1602614. Sanders, J.A., Verhulst, F., Murdock, J., 2007. Averaging methods in nonlinear dynamical systems. volume 59. Springer. Schmid, P.J., Li, L., Juniper, M.P., Pust, O., 2011. Applications of the dynamic mode decomposition. Theoretical and computational fluid dynamics 25, 249–259. 47
Schmidt, M., Lipson, H., 2009. Distilling free-form natural laws from experimental data. Science 324, 81–85. Van Der Maaten, L., Postma, E.O., Van Den Herik, H.J., et al., 2009. Dimensionality reduction: A comparative review. Journal of Machine Learning Research 10, 1–41. Vandenberghe, L., 2010. The cvxopt linear and quadratic cone program solvers. Online: http://cvxopt. org/documentation/coneprog. pdf 53. Volosov, V.M., 1962. Averaging in systems of ordinary differential equations. Russ. Math. Surv. 17, 1. Ypma, T.J., 1995. Historical development of the newton–raphson method. SIAM review 37, 531–551. Yu, Z., Ding, J., Li, Y., 2025. Discover network dynamics with neural symbolic regression. Nat. Comput. Sci. , 1–13. Yuan, Y., Tang, X., Zhou, W., Pan, W., Li, X., Zhang, H.T., Ding, H., Goncalves, J., 2019. Data driven discovery of cyber physical systems. Nat. Commun. 10, 4894. Zega, V., Gattere, G., Koppaka, S., Alter, A., Vukasin, G.D., Frangi, A., Kenny, T.W., 2020. Numerical modelling of non-linearities in mems resonators. Journal of Microelectromechanical Systems 29, 1443–1454. Zhai, J., Zhang, S., Chen, J., He, Q., 2018. Autoencoder and its various variants, in: 2018 IEEE international conference on systems, man, and cybernetics (SMC), IEEE. pp. 415– 419.
48