Paper deep dive
A neural operator framework for data-driven discovery of stability and receptivity in physical systems
Chengyun Wang, Liwei Chen, Nils Thuerey
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 97%
Last extracted: 4/26/2026, 11:22:08 PM
Summary
The paper introduces a data-driven neural operator framework designed to discover the stability and receptivity properties of complex physical systems without requiring explicit governing equations. The method involves training a neural network (NN) as a dynamics emulator and using automatic differentiation to extract its Jacobian. This Jacobian serves as a proxy for the system's local linear operator, enabling both linear stability analysis (via eigenvalue decomposition) and resolvent analysis (via weighted SVD in a reduced subspace). The framework is demonstrated on the Lorenz-63 chaotic model and high-dimensional fluid flows, showing that it can effectively capture dominant instability modes and optimal forcing responses even in strongly nonlinear regimes, outperforming traditional data-driven methods like Dynamic Mode Decomposition (DMD) in nonlinear settings.
Entities (8)
Relation Signals (5)
Lorenz-63 → istestedon → Neural Operator Framework
confidence 100% · We first consider the canonical Lorenz system as a simple model...
Neural Operator Framework → performs → Linear Stability Analysis
confidence 100% · ...to perform both linear stability analysis via an eigenvalue decomposition to find unstable modes...
Neural Operator Framework → performs → Resolvent Analysis
confidence 100% · ...and for resolvent analysis via a weighted SVD on the projected resolvent...
Neural Network (NN) Emulator → provides → Jacobian
confidence 100% · The Jacobian of the trained emulator is then extracted via automatic differentiation...
Neural Operator Framework → uses → Neural Network (NN) Emulator
confidence 100% · By training a neural network as a dynamics emulator...
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Understanding how complex systems respond to perturbations, such as whether they will remain stable or what their most sensitive patterns are, is a fundamental challenge across science and engineering. Traditional stability and receptivity (resolvent) analyses are powerful but rely on known equations and linearization, limiting their use in nonlinear or poorly modeled systems. Here, we introduce a data-driven framework that automatically identifies stability properties and optimal forcing responses from observation data alone, without requiring governing equations. By training a neural network as a dynamics emulator and using automatic differentiation to extract its Jacobian, we can compute eigenmodes and resolvent modes directly from data. We demonstrate the method on both canonical chaotic models and high-dimensional fluid flows, successfully identifying dominant instability modes and input-output structures even in strongly nonlinear regimes. By leveraging a neural network-based emulator, we readily obtain a nonlinear representation of system dynamics while additionally retrieving intricate dynamical patterns that were previously difficult to resolve. This equation-free methodology establishes a broadly applicable tool for analyzing complex, high-dimensional datasets, with immediate relevance to grand challenges in fields such as climate science, neuroscience, and fluid engineering.
Tags
Links
- Source: https://arxiv.org/abs/2604.19465v1
- Canonical: https://arxiv.org/abs/2604.19465v1
Trouble viewing inline? Open PDF directly →
Full Text
75,214 characters extracted from source content.
Expand or collapse full text
A neural operator framework for data-driven discovery of stability and receptivity in physical systems Chengyun Wang1, Liwei Chen2 & Nils Thuerey1,∗ 1School of Computation, Information and Technology, Technical University of Munich, Boltzmannstrasse 3, D-85748 Garching bei München, Germany 2Beijing Institute of Astronautical Systems Engineering, Beijing, China ∗Corresponding author: nils.thuerey@tum.de Abstract Understanding how complex systems respond to perturbations, such as whether they will remain stable or what their most sensitive patterns are, is a fundamental challenge across science and engineering. Traditional stability and receptivity (resolvent) analyses are powerful but rely on known equations and linearization, limiting their use in nonlinear or poorly modeled systems. Here, we introduce a data-driven framework that automatically identifies stability properties and optimal forcing responses from observation data alone, without requiring governing equations. By training a neural network as a dynamics emulator and using automatic differentiation to extract its Jacobian, we can compute eigenmodes and resolvent modes directly from data. We demonstrate the method on both canonical chaotic models and high-dimensional fluid flows, successfully identifying dominant instability modes and input–output structures even in strongly nonlinear regimes. By leveraging a neural network–based emulator, we readily obtain a nonlinear representation of system dynamics while additionally retrieving intricate dynamical patterns that were previously difficult to resolve. This equation-free methodology establishes a broadly applicable tool for analyzing complex, high-dimensional datasets, with immediate relevance to grand challenges in fields such as climate science, neuroscience, and fluid engineering. 1 Introduction Real-world systems are inherently nonlinear and exhibit rich, multi-scale behavior in both space and time. Dynamical systems theory offers a unified mathematical framework to describe, analyze, and predict the evolution of such systems, capturing the complex interactions among quantities that evolve over time. From its foundations in the seminal work of Poincaré on celestial mechanics, this field has expanded to encompass a vast range of phenomena across science (Jendrejack et al., 2003; Sun and Wu, 2015) and engineering (Holmes, 2012; Battat et al., 2022). Within this universal framework, a fundamental challenge persists: to understand and predict how a system responds to perturbations. Answering this involves determining if a system’s state is stable, identifying the conditions that trigger transitions to instability, and discovering which external inputs will most effectively steer the system’s evolution. Given that nonlinearity remains a primary challenge in analyzing and controlling dynamical systems, a leading perspective is to investigate the system’s local behavior through linearization around a fixed point to extract the dominant modes of perturbation evolution. This strategy has given rise to a powerful suite of modal analysis techniques, many of which were developed within the hydrodynamics and cybernetics community (Simon and Mitter, 1968; Proctor et al., 2016; Rowley and Dawson, 2017; Sun et al., 2020). These techniques can be broadly divided into two categories (Taira et al., 2017). The first consists of operator-based methods, such as linear stability analysis (Theofilis, 2011) and resolvent analysis (McKeon and Sharma, 2010), which are built upon the linearized operator derived from the system’s governing equations. Their profound insights, however, are contingent on this explicit, model-based knowledge. To address systems not amenable to first principles, a second category of data-driven methods has emerged. These approaches, including dynamic mode decomposition (DMD) (Schmid, 2010), proper orthogonal decomposition (POD) (Holmes, 2012), and sparse identification of nonlinear dynamics (SINDy)(Brunton et al., 2016b), rely purely on observational data to uncover underlying dynamical structures, offering a pathway to analyze systems where governing equations are unknown. Linear stability analysis and resolvent analysis, as two representative operator-based methods, have been demonstrated to be highly effective in providing physical insights for all kinds of systems (Bertozzi and Brenner, 1997; Kolokolnikov et al., 2011; Symon et al., 2018). The former linearizes the governing equations around the base state and assumes exponential perturbation growth. This assumption elegantly converts the linear initial-value problem into an eigenvalue problem for the Jacobian, whose eigenvalues can indicate whether the base state is stable. However, this framework is limited to describing the asymptotic perturbation evolution and does not account for short-term transient behavior (Trefethen et al., 1993; Schmid, 2007). To remedy this shortcoming, resolvent analysis was developed to investigate the linearized system’s receptivity to external forcings. Instead of focusing on internal instabilities, it can identify the optimal forcing patterns that elicit the most amplified responses, which usually represents the natural starting point for control design (Luhar et al., 2014; Yeh and Taira, 2019). As an input-output response framework, it provides not only a transfer-function viewpoint of analyzing the system’s receptivity characteristics, but also reveals the transient growth of perturbations, thereby explaining the transient energy amplification that arises in linearly stable systems. Therefore, linear stability analysis and resolvent analysis complement each other to provide a comprehensive picture of the system’s linearized dynamics. Figure 1: Schematic of our data-driven algorithm. (1) Time-resolved state snapshots collected from a dynamical system are used. (2) An N emulator is trained using a rollout strategy to accurately learn the system’s evolution map. (3) The Jacobian of the trained emulator is then extracted via automatic differentiation and used to approximate the local Jacobian. This data-driven operator is the key ingredient for (4) N-based modal analysis to perform both linear stability analysis via an eigenvalue decomposition to find unstable modes, and for resolvent analysis via a weighted SVD on the projected resolvent to identify optimal forcing modes, response modes, and the gain distribution. Despite the effectiveness of operator-based methods, their practical application to many dynamical systems of interest remains very challenging. In numerous scientific frontiers, such as neuroscience, epidemiology, and ecology, there is a fundamental lack of known physical laws from which to derive the governing equations of motion. Even in systems where such equations are well-established, such as turbulence and combustion, the construction of the linearized operator often requires specialized numerical solvers or adjoint simulations (Rolandi et al., 2024). Meanwhile, the rapid progress in experimental measurements and numerical simulations has led to data generation at an unprecedented scale. This data deluge is driving a paradigm shift to data-driven methods that can extract meaningful information directly from datasets without prior knowledge. Among which, DMD has been widely applied beyond its original domain of fluid dynamics (Schmid, 2010), in large part due to its simple framing in terms of linear regression and strong connections to nonlinear dynamical systems via Koopman spectral theory (Rowley et al., 2009). In addition, these two classical categories are being complemented by a promising hybrid perspective that seeks to perform operator-based analysis in a purely data-driven manner. A leading example is a DMD-based approach to resolvent analysis for linearly stable systems, in which DMD is applied to time-resolved flow snapshots to approximate the system’s eigenmodes. The resolvent operator is then constructed within the reduced-order subspace spanned by these DMD modes (Herrmann et al., 2021). This idea is further refined by restricting the admissible DMD model space to certain matrix manifolds that preserve the desired physical properties (Baddoo et al., 2023). However, the underlying assumption of DMD that the system evolution can be approximated by a global linear operator leads to a severely contaminated approximation of the eigenmodes when handling nonlinear dynamical systems. It thus limits the applicability to linear or weakly nonlinear regimes. Neural networks (NNs) offer a pathway beyond these limitations. As universal function approximators, N emulators can learn complex nonlinear dynamics from spatiotemporal measurements of the system(Colen et al., 2021; Floryan and Graham, 2022; Yu and Wang, 2024; Page et al., 2024). They can be trained either in a purely data-driven manner (Morton et al., 2018; Thuerey et al., 2020; Chen and Thuerey, 2023) or augmented with physical information (Bar-Sinai et al., 2019; Raissi et al., 2019; Kochkov et al., 2021; List et al., 2022). Recent works also demonstrated that the emulators trained on airfoil datasets can capture intrinsic aerodynamic characteristics, such as Kelvin–Helmholtz instabilities and shock wave motion (Chen and Thuerey, 2024). Crucially, these emulators, which can be viewed as learned approximations of classical numerical solvers, are fully differentiable. This property is the key to our approach: the local Jacobian of the trained N can be computed at any point in the state space via automatic differentiation, opening the door to an equation-free, data-driven modal analysis for nonlinear systems. To this end, we propose a data-driven approach to learn a N emulator from data generated by truly nonlinear systems, with a schematic illustration presented in Fig. 1. The Jacobian of the learned N emulator serves as an approximation of the system’s local linear operator, which is subsequently used to perform operator-based modal analysis. To demonstrate the power and generality of this framework, we diagnose four representative dynamical systems of increasing complexity, successfully extracting their essential stability characteristics and input-output resolvent modes. Furthermore, by comparing with traditional data-driven methods like DMD, we provide insights on how to select the appropriate tool for data-driven analysis. 2 Results Operator-based Jacobian of the system Considering a forced nonlinear dynamical system (temporally continuous but spatially discretized into N degrees of freedom), the governing equation can be expressed in terms of the state variable ∈ℂN q ^N in the compact operator form as ddt=()+, d qdt=N( q)+ f, (1) where N is the nonlinear operator collecting all the terms except for the time derivative and ∈ℂN f ^N is an external forcing term. To linearize this ordinary differential equation, the state variable can be decomposed as a sum of base state b∈ℂN q_b ^N and infinitesimal perturbation ′∈ℂN q ^N, which yields dbdt+d′dt=(b+′)+=(b)+∂|b′+O(|′|2)+. d q_bdt+ d q dt=N( q_b+ q )+ f=N( q_b)+ ∂ q |_ q_b q +O(| q |^2)+ f. (2) While the base state could be chosen as a time-averaged solution, we focus on the stability around steady-state solutions, and thus b q_b is an equilibrium point of Eq. 1 satisfying db/dt=(b)=0d q_b/dt=N( q_b)=0. After neglecting the higher-order infinitesimals in Eq. 2, we obtain the linearized equation of the perturbation ′ q as d′dt=∂|b′+=′+, d q dt= ∂ q |_ q_b q + f=A q + f, (3) where ∈ℂN×NA ^N× N represents the discretized linear operator of this continuous dynamical system at base state b q_b, i.e., the local Jacobian. Based on this Jacobian A, operator-based modal analyses such as linear stability and resolvent analysis can be performed; for details see Sec. 4. Data-driven neural network-based Jacobian Our data-driven approach adopts a time-stepping framework using a N emulator to evolve the original unforced dynamical system of Eq. 1. The state n q_n at time t is mapped to the subsequent state n+1 q_n+1 as: n+1=fθ(n), q_n+1=f_θ( q_n), (4) where fθf_θ is an arbitrary N architecture with weights θ. Let P represent a classical numerical solver that also advances the system state with a time step Δt t. Then the learning task is to find the optimal weights θopt _opt of the N to best approximate the numerical solver, i.e., fθ≈Pf_θ≈ P. When trained on a sufficiently diverse dataset, the N emulator can capture the underlying latent dynamics and generalize effectively to previously unseen system states. Since the N emulator operates in a temporally discrete manner, whereas operator-based analysis relies on temporally continuous equations, an essential step is to bridge the gap between these two representations of the dynamical system. We can also linearize the temporally discrete dynamical system (Eq. 4) around the equilibrium points b=fθ(b) q_b=f_θ( q_b) based on the base-perturbation decomposition: n+1′=∂fθ∂|bn′=n′, q _n+1= ∂ f_θ∂ q |_ q_b q _n=N q _n, (5) where the Jacobian N of our learned neural operator fθf_θ can be calculated via automatic differentiation evaluated at b q_b (Part 3, Fig. 1). The differentiability of the N emulator actually allows for the use of automatic differentiation to evaluate N at any input state. In the absence of external forcing, the general solution of Eq. 3 takes an exponential form, and thus two successive perturbation states ′ q should satisfy: n+1′=exp(Δt)n′ q _n+1= (A t) q _n. Hence, we can derive the relationship between the N-based Jacobian N and the operator-based Jacobian A as =exp(Δt).N= (A t). (6) Benefiting from this relationship, if a learned surrogate model can accurately capture the underlying system dynamics, it can then be utilized to evaluate the stability of the original system without knowing its governing equations. N-based linear stability and resolvent analysis According to Eq. 6, the eigenvalues of the N-based and operator-based Jacobians, denoted by the diagonal matrices _N and _A, respectively, should satisfy the exponential relationship: =exp(Δt) _N= ( _A t). The associated eigenvectors corresponding to each eigenvalue should be identical, i.e., =V_N=V_A. Provided the learned model accurately approximates the original system, the extracted eigenvalues and eigenvectors from N serve as good approximations of those from A, whose spectral characteristics can thereby be inferred indirectly (Part 4a, Fig. 1). This constitutes the core procedure of N-based linear stability analysis. In essence, the entire workflow can be interpreted as a data-driven eigendecomposition, enabling stability assessment and extraction of dominant unstable modes without requiring explicit governing equations. As for resolvent analysis, a natural idea is to replace the Jacobian A by log()/Δt (N)/ t to compute the resolvent operator (ω)=(−iω−)−1H(ω)=(-i -A)^-1. However, it is not advisable to directly take the logarithm of the Jacobian N, as the logarithmic function is a multivalued function in the complex domain, and as such, is strongly underdetermined. N typically corresponds to a large number of A that satisfy Eq. 6. Instead, we can leverage the results of the linear stability analysis to implement the dimensionality reduction (Reddy and Henningson, 1993). In a nutshell, we can project the forced system Eq. 3 into the subspace expanded by a truncated r eigenvectors (Schmid and Henningson, 2012; Herrmann et al., 2021). Supposing the original state variable ′ q and f in Eq. 3 can be expressed in the subspace spanned by the first r linearly independent eigenvectors ,rV_A,r, we have d,rdt=,r+,r, dV_A,r xdt=AV_A,r x+V_A,r y, (7) where =[x1,x2,⋯,xr]T∈ℂr x=[x_1,x_2,·s,x_r]^T ^r and =[y1,y2,⋯,yr]T∈ℂr y=[y_1,y_2,·s,y_r]^T ^r are the vectors of expansion coefficients in the eigenvector coordinate. Based on the definition of eigenvalue ,r=,r,rAV_A,r=V_A,r _A,r, we can obtain the new governing equation in the eigenvector coordinate: ddt=,r+. d xdt= _A,r x+ y. (8) Given then weighting matrix Q used to define the state norm in the original coordinate as ‖^‖2=^∗^|| q||^2_Q= q^*Q q, the corresponding weighting matrix in the new coordinate becomes ~=,r∗,r Q=V_A,r^*QV_A,r, which can be factorized as ~=~∗~ Q= F^* F via Cholesky factorization. Under this new coordinate system, we can carry out a weighted SVD of the projected resolvent: ~(−iω−,r)−1~−1=~(ω)(ω)~∗(ω), F(-i - _A,r)^-1 F^-1= _ F(ω) (ω) ^*_ F(ω), (9) where the diagonal entries of the singular value matrix Σ(ω) (ω) are the resolvent gains and the column vectors in ~(ω) _ F(ω) and ~(ω) _ F(ω) are the corresponding response and forcing modes at frequency ω. The conversion of resolvent modes back to the original space is then performed via =,r~−1~(ω) =V_A,r F^-1 _ F(ω) and =,r~−1~(ω) =V_A,r F^-1 _ F(ω). In practice, the derivation of ,r _A,r and ,rV_A,r requires the knowledge of the corresponding governing equation, so ,r _N,r and ,rV_N,r acquired from data are employed instead (Part 4b, Fig. 1). However, it is important to avoid indiscriminately increasing the dimensionality of the subspaces (Herrmann et al., 2021). A reliable projection of a high-dimensional system onto a reduced subspace requires that the chosen eigenvectors be linearly independent. But especially in non-normal systems, where eigenvectors are generally non-orthogonal, care must be taken not to select an overly large number of modes r. An excessively large r may lead to a redundant coordinate system, resulting in a non-unique representation of the system state in terms of eigenvector coordinates. This will deteriorate the validity of resolvent operator estimation instead of enriching the basis. Example 1: Lorenz-63 We first consider the canonical Lorenz system as a simple model based on ordinary differential equations (ODEs). Although it is a chaotic system, NNs have been proven capable of characterizing its temporal dynamics, including accurately advancing the solution in time for unseen initial conditions (Brunton and Kutz, 2022). The state of the Lorenz system given by =[x,y,z]T q=[x,y,z]^T is governed by dxdt dxdt =σ(y−x), =σ (y-x ), (10) dydt dydt =x(ρ−z)−y, =x (ρ-z )-y, dzdt dzdt =xy−βz. =xy-β z. Here, the standard parameter values of σ=10σ=10, ρ=28ρ=28, and β=8/3β=8/3 are chosen, and thus the system exhibits chaotic motion confined to a strange attractor. This system has three equilibrium points, one at the origin ([0,0,0]T[0,0,0]^T) and two symmetric off-origin points for ρ>1ρ>1 ([±β(ρ−1),±β(ρ−1),ρ−1][± β(ρ-1),± β(ρ-1),ρ-1]). The analytical local Jacobian around a base state b=[xb,yb,zb]T q_b=[x_b,y_b,z_b]^T is: =∂|b=[−σ0ρ−zb−1−xbybxb−β].A= ∂ q |_ q_b= bmatrix-σ&σ&0\\ ρ-z_b&-1&-x_b\\ y_b&x_b&-β bmatrix. (11) Figure 2: Time evolution of the true Lorenz system (red solid) compared with the N emulator prediction (green dotted) from an unseen initial point (blue circle). (a) Trajectories in phase space; (b) Component-wise time evolution series. Even when trained on the trajectory dataset with a one-step rollout training strategy (Training setup in Sec. 4), it can be used to autoregressively generate new trajectories as far into the future as desired, even for initial conditions not present in the training dataset. Fig. 2(a) presents the trajectory comparison between the N emulator’s prediction and ground truth starting from a randomly selected initial condition over an interval of 8 time units. Fig. 2(b) offers a clearer view of the approximation’s performance by showing the time evolution of three components separately. The close agreement between both simulations in the early stage suggests that the trained N emulator captures the system dynamics with high fidelity, effectively mimicking the update rule of a classical Runge–Kutta scheme. Despite a noticeable deviation observed for the N emulator at the end of the simulation, the performance remains commendable given the lack of rollout steps at training time. Base State b q_b N-based Jacobian N Analytical Results exp(Δt) (A t) equilibrium points [000] bmatrix0\\ 0\\ 0 bmatrix [0.91850.0947−0.00010.26571.00360.00050.0002−0.00010.9735] bmatrix0.9185&0.0947&-0.0001\\ 0.2657&1.0036&0.0005\\ 0.0002&-0.0001&0.9735 bmatrix [0.91790.095100.26631.00350000.9737] bmatrix0.9179&0.0951&0\\ 0.2663&1.0035&0\\ 0&0&0.9737 bmatrix [626227] bmatrix6 2\\ 6 2\\ 27 bmatrix [0.90540.0945−0.00410.00740.9861−0.08360.07890.08800.9705] bmatrix0.9054&0.0945&-0.0041\\ 0.0074&0.9861&-0.0836\\ 0.0789&0.0880&0.9705 bmatrix [0.90520.0946−0.00410.00600.9869−0.08320.08000.08730.9700] bmatrix0.9052&0.0946&-0.0041\\ 0.0060&0.9869&-0.0832\\ 0.0800&0.0873&0.9700 bmatrix [−62−6227] bmatrix-6 2\\ -6 2\\ 27 bmatrix [0.90500.09480.00390.00710.98650.0835−0.0812−0.08690.9697] bmatrix0.9050&0.0948&0.0039\\ 0.0071&0.9865&0.0835\\ -0.0812&-0.0869&0.9697 bmatrix [0.90520.09460.00410.00600.98690.0832−0.0800−0.08730.9700] bmatrix0.9052&0.0946&0.0041\\ 0.0060&0.9869&0.0832\\ -0.0800&-0.0873&0.9700 bmatrix non-equilibrium points [1.46898.7045−7.6905] bmatrix1.4689\\ 8.7045\\ -7.6905 bmatrix [0.92210.09410.00010.33881.0075−0.01800.08550.02380.9734] bmatrix0.9221&0.0941&0.0001\\ 0.3388&1.0075&-0.0180\\ 0.0855&0.0238&0.9734 bmatrix [0.92150.0952−0.00070.33931.0071−0.01450.08470.01870.9736] bmatrix0.9215&0.0952&-0.0007\\ 0.3393&1.0071&-0.0145\\ 0.0847&0.0187&0.9736 bmatrix Table 1: Comparison of the N-based Jacobian matrix N and the operator-based Jacobian matrix A around three equilibrium points and one non-equilibrium point (the initial point in Fig. 2). Meanwhile, we leverage automatic differentiation to compute the N-based Jacobians around the three equilibrium points and compare them with the analytical Jacobians derived from Eq. 11. As shown in Tab. 1, the N-based and analytical Jacobians exhibit an almost perfect element-wise match. This suggests that the following modal analysis performed using the N-based Jacobian can reproduce theoretical results, despite the fact that the N itself is trained purely on data and is agnostic to the underlying governing equations. This paves the way for data-driven modal analysis approaches for both linear and nonlinear scenarios in higher-dimensional settings. Furthermore, we found that for this simple ODE system, the Jacobian agreement is not limited to equilibrium points but instead extends to arbitrary points along the trajectory (see the last row of Tab. 1). This indicates that the trained neural operator can also faithfully represent local dynamical behavior at any points and may facilitate tasks such as computing finite-time Lyapunov exponents (FTLE) (Shadden et al., 2005). This is not the primary focus of our investigation. Example 2: Complex Ginzburg-Landau equation The second example considered is the one-dimensional complex Ginzburg–Landau equation (Bagheri et al., 2009). It is widely recognized as a reduced-order model for analyzing convective and absolute instabilities (Huerre and Monkewitz, 1985; Sreenivasan et al., 1989; Kyle and Sreenivasan, 1993; Utada et al., 2008) in fluid systems under varying parameter configurations. As its linearized operator can be obtained analytically through the manipulation of the differential matrices, it has become a benchmark for the evaluation and development of stability analysis techniques (Herrmann et al., 2021; Martini et al., 2021). Figure 3: Sample trajectories from datasets of (a) linear and (b) nonlinear complex Ginzburg-Landau systems. The nonlinear complex Ginzburg–Landau equation defined on an infinite interval x∈(−∞,∞)x∈(-∞,∞) is written as ∂t=(−ν∂x+γ∂2∂x2+μ(x))−a||2=−a||2, ∂ q∂ t= (-ν ∂ x+γ ∂^2∂ x^2+μ(x) ) q-a| q|^2 q\\ =A q-a| q|^2 q, (12) where −a||2-a| q|^2 q is the cubic nonlinear term. Since =0 q=0 is an equilibrium point, the linearized system around this point simply involves omitting the cubic term and thus A is the linearized operator. The advection and dispersion properties are controlled by the complex terms ν=U+2icuν=U+2ic_u and γ=1+icdγ=1+ic_d, respectively. The real-valued term μ(x)=μ0−cu2+μ2x2/2μ(x)= _0-c_u^2+ _2x^2/2 is defined as a quadratic function to model exponential instabilities (Bagheri et al., 2009). A detailed description of the parameters can be found in Tab. 2. The homogeneous boundary conditions are considered at x→±∞x→±∞, so that perturbations are allowed to grow and decay throughout the entire domain. Variable Description Value U Mean advection velocity 2.0 cuc_u Most unstable wavenumber 0.2 cdc_d Dispersion parameter -1.0 μ0 _0 Bifurcation parameter 0.23(0.41) μ2 _2 Degree of non-parallelism -0.01 a Nonlinearity 0(1) Table 2: Parameter descriptions and settings for the complex Ginzburg–Landau system. Values corresponding to the nonlinear scenario are indicated in parentheses. Figure 4: Comparison of the real part of the Jacobians between the N-based operator and the operator-based ground truth. The first and second rows correspond to the linear and nonlinear scenarios, respectively. We aim to demonstrate the applicability of the N emulator approach to both linear and nonlinear regimes even in complex number space. For this purpose, we consider two datasets derived respectively from the linear (Herrmann et al., 2021) and nonlinear (Bagheri et al., 2009) forms of this system. The sample trajectories and the specific parameter settings used in each scenario are presented in Fig. 3 and Tab. 2. A truncation level of r=24r=24 is used for the eigenbasis in this example. As shown in Fig. 4, we compare the Jacobian heatmaps of the learned emulator with those of the analytical operator for both linear and nonlinear scenarios. We can find that the network perfectly recovers the operator in the linear scenario, whereas for the nonlinear scenario, it successfully captures the main features but exhibits minor discrepancies in finer details. These deviations are most likely caused by nonlinear mappings being inherently more difficult to approximate than linear ones. Figure 5: Linear stability analysis results of the linear and nonlinear Ginzburg-Landau systems. (a)(b) Eigenspectra obtained from the N-based Jacobian (∘ ), DMD-based Jacobian (++) and operator-based ground truth (∙ ). The dashed line denotes the stability boundary in the complex plane. (c)(d) The leading three eigenmodes of the N-based Jacobians in both scenarios, where solid and dashed lines show the real part and magnitude of the modes. The thick gray lines in the background show operator-based ground-truth for comparison. In the context of linear stability analysis, Fig. 5(a) and (b) compare the eigenspectra of the N-based and operator-based Jacobians for linear and nonlinear scenarios, respectively. For reference, the eigenspectra derived from the DMD model are also included. Fig. 5(c) and (d) show the leading three eigenmodes identified by the N emulator for both scenarios. For the linear scenario, it is anticipated that the DMD-based results show excellent agreement with the analytical solution, as DMD is fundamentally designed to approximate linear operators. Even so, the N-based results also show excellent alignment of the leading eigenvalues and eigenvectors with the analytical reference. For the nonlinear scenario, DMD-based approach totally breaks down, as demonstrated by the identification of more than one unstable mode located in the right half of the complex plane (see Fig. 5(b)). In contrast, the N-based Jacobian still yields an accurate leading unstable eigenvalue, even though it differs in finer structures according to Fig. 4. This suggests that the N emulator can effectively identify the most dominant mode regardless of the dynamical properties of datasets. While the remaining higher-order modes deviate more noticeably, they decay rapidly and contribute minimally to the long-term behavior of perturbations. Figure 6: Resolvent analysis results of the complex Ginzburg–Landau systems. The first and the second row correspond to the linear and nonlinear scenarios, respectively. (a)(c) Resolvent gain distribution for the first three modes with respect to the frequency. (b)(d) The first three forcing and response modes at the peak gain frequency. The thick gray lines in the background show operator-based ground-truth for comparison. In the context of resolvent analysis, Fig. 6(a) and (c) show the gain curves of the first three resolvent modes as a function of the forcing frequency for the two scenarios. The corresponding dominant response and forcing modes at the peak gain frequency are presented in Fig. 6(b) and (d). For the linear scenario, our N-based approach can yield highly accurate gain curves and mode shapes. For the nonlinear scenario, it also succeeds in predicting the gain curve near the dominant frequency, and accurately reconstructs the response and forcing modes corresponding to this frequency. Though applicable to both linear and nonlinear regimes, the higher-order mode shapes predicted by the N emulator become increasingly influenced by noise and deviate from ground truth more noticeably, as shown in Fig. 6(c) and (d). There are two important aspects to note for the gain distributions of Fig. 6(c). First, the N emulators have worse performance for gain prediction at frequencies far from the peak gain. This phenomenon can be interpreted via dyad expansion of the resolvent operator (Symon et al., 2018). The resolvent operator H(ω)H(ω) can be expressed as a sum of dyadic products of left eigenvectors ~j g_j and right eigenvectors ~j h_j, each weighted by the distance between the eigenvalue and the imaginary axis 1/(iω−λj)1/(iω- _j) as H(ω)=(iω−)−1=∑j=1n1iω−λj~j~j∗≈1iω−λ1~1~1∗,H(ω)=(iω I- Q)^-1= _j=1^n 1iω- _j g_j h_j^*≈ 1iω- _1 g_1 h_1^*, (13) where we can see that near the frequency of maximal gain, i.e., the forcing frequency ω is close to the imaginary part of the dominant eigenvalue, the resolvent operator admits a good rank-1 approximation. Since the leading eigenmode is accurately recovered (see Fig. 5), the corresponding resolvent operator is also highly accurate for those frequency intervals satisfying the rank-1 approximation. Second, the N-based method tends to exhibit more noise in the predicted forcing modes compared to the response modes. However, this behavior is in line with other data-driven approaches and can be attributed to the same cause: in non-normal systems, the spatial structure of the direct eigenvectors resembles that of the response modes, but can differ significantly for the forcing modes (Herrmann et al., 2021). Example 3: 2D transitional channel flow The third example is the two-dimensional plane Poiseuille flow with Re = 2000 (based on the centerline velocity and the channel half-height) governed by the incompressible Navier-Stokes equations. The analytical derivation of the linearized operator for this example can be obtained by considering perturbations in the wavenumber space, which can be reduced to the Orr–Sommerfeld (OS) equation (Schmid and Henningson, 2012). The dynamics of the velocity perturbations with a particular streamwise wavenumber kxk_x is governed by the OS operator, which we use for the ground truth reference (Details in Supplementary Materials). Figure 7: Snapshots of the wall-normal velocity perturbation field during the transient evolution stage of two-dimensional channel flow. The initial perturbation profiles are the same but with different magnitude scaling factors ϵε: (a) weakly nonlinear scenario with ϵ=10−6ε=10^-6; (b) nonlinear scenario with ϵ=10−1ε=10^-1. To simplify the analysis, the flow perturbations with a single wavenumber kx=1k_x=1 are advanced in time. The initial velocity fields are of the form u(x,y) u(x,y) =ϵ⋅ℜ[dϕ(y)dy⋅eikxx], =ε· [ dφ(y)dy· e^ik_xx ], (14) v(x,y) v(x,y) =−ϵ⋅ℜ[i⋅ϕ(y)⋅eikxx] =-ε· [i·φ(y)· e^ik_xx ] where ϕ(y)φ(y) is the initial profile constructed as a superposition of multiple random Gaussian profiles satisfying specific boundary conditions. ϵε serves as a controllable scaling factor for the initial perturbation magnitude, determining the strength of the nonlinear effects on the fluid behavior during the transient phase (Fig. 7). We consider a weakly nonlinear dataset with ϵ=10−6ε=10^-6 and a strongly nonlinear dataset with ϵ=10−1ε=10^-1. As shown in Fig. 7, the flow field snapshots sampled from the respective datasets illustrate that, the decaying perturbations in the weakly nonlinear dataset preserve symmetric structures during their evolution, whereas those in the strongly nonlinear dataset exhibit symmetry-breaking phenomena. A truncation level of r=30r=30 is used for the eigenbasis of the projection subspace. Figure 8: Comparison of the eigenspectra obtained from the N-based Jacobian (∘ ), DMD-based Jacobian (++) and operator-based ground truth (∙ ) for the 2D channel flow system. The dashed line denotes the stability boundary in the complex plane. (a) Weakly nonlinear dataset. (b) Nonlinear dataset. Figure 9: The spatial structures (eigenvectors) of the leading A and P modes of the N-based Jacobian obtained from weakly nonlinear and nonlinear datasets. The first column is the reference results obtained from the OS operator. For the chosen Reynolds number of 2000, the channel system is stable (Schmid and Henningson, 2012), which can be verified by the eigenspectrum of the OS operator from Fig. 8. All operator-based eigenvalues fall in the stable complex half-plane and they are located on three branches forming a Y-shaped pattern. These branches have been labeled A (cr→0c_r→ 0), P (cr→1c_r→ 1), and S (cr≈2/3c_r≈ 2/3) by (Mack, 1976). The branches describe distinct features of the perturbation dynamics: P modes (center modes) on the P-branch with much higher phase speeds capture the fast dynamics in the center of the channel, and A modes (wall modes) with rather small phase velocities capture the slow dynamics near the channel wall. The S modes approaching a phase speed 2/32/3 are highly damped. In Fig. 8, we note that the first leading eigenvalues of the A-branch and P-branch agree well between the results of the OS operator and the N emulator trained on both weakly nonlinear and nonlinear datasets. This indicates that the N emulator can successfully identify the two most dominant perturbation structures near the wall and centerline, i.e., the least-damped A and P modes (see Fig. 9), despite the strong nonlinearities of the dataset. While the remaining higher-order A and P modes appear to be spurious, and the S modes with high decay rates were not recovered, their contribution to the system’s dynamics is negligible. We attribute this to limitations in the numerical precision of the dataset, as well as the inherent noise in the trained emulator, which potentially obscures the evolution of these highly damped modes. Additionally, the randomly dispersed DMD-based eigenvalues in Fig. 8(b) again highlight its inherent flaws when handling nonlinear datasets. It is noteworthy that another spurious eigenvalue appears on the imaginary axis close to the stability boundary for the case of the nonlinear dataset (Fig. 8(b)), which is absent for the operator-based ground truth and the case of the weakly nonlinear dataset. This discrepancy may be attributed to the preprocessing applied to the training data. Given the negligible magnitude of streamwise perturbations, the raw data were preprocessed by subtracting the base flow and subsequently normalized. While such treatment is reasonable in the linear or weakly nonlinear regime, it can cause distortions in nonlinear scenarios where the linear superposition principle does not hold. Consequently, zero-frequency modes associated with non-oscillatory components were introduced in the dataset. These modes are subsequently identified by the N emulator during feature extraction, as evidenced by their appearance along the imaginary axis in Fig. 8(b). Figure 10: Resolvent analysis results of the 2D channel flow system. (a) Resolvent gain curves of the analytical operator and the N-based resolvent operators derived from weakly nonlinear and nonlinear datasets. (b) Forcing and response modes at the frequency with maximum gain. Results are shown from top to bottom: operator-based, N-based (weakly nonlinear), N-based (nonlinear). As for the resolvent analysis, Fig. 10(a) displays the gain curves of three distinct operators in response to harmonic forcings with frequencies in the range of [0,0.5][0,0.5]. The peak frequency at ω=0.31ω=0.31 is accurately predicted by the N emulator despite the nonlinear nature of this scenario. We note that the corresponding resolvent gain curve exhibits an extra peak at zero frequency for the nonlinear scenario. This phenomenon aligns with the presence of non-oscillatory spurious modes discussed above, and the additional gain peak comes from resonating with them. Nevertheless, the N emulator accurately predicts the physically relevant peak frequency of ω=0.31ω=0.31 originating from a near-resonance with the first A-branch eigenmode for both scenarios. The resolvent modes of streamwise and wall-normal velocity perturbations at ω=0.31ω=0.31 are illustrated in Fig. 10(b), revealing that the N emulator detects the area that is most sensitive for perturbation amplification for channel flow at the current Reynolds number. Since the leading eigenvectors of the N-based Jacobian align closely with those of the true system (Fig. 9), it naturally leads to accurate resolvent mode predictions. Similar to the observations in the Ginzburg–Landau case, the reconstructed response mode shows near-exact agreement with the ground truth, whereas the reconstructed forcing mode yields a less accurate approximation, which is also attributed to the non-normal nature of the channel flow. This behavior is more prominent in the context of nonlinear datasets, indicating that nonlinearity may exacerbate the challenges associated with training under the current neural network architecture. Example 4: 2D cylinder flow As the last example for the N-based modal analysis, we consider a flow with an immersed obstacle, the two-dimensional cylinder flow at Re = 100 (based on the cylinder diameter and inflow velocity). This Reynolds number setting is greater than the critical threshold, and the cylinder wake flow experiences a Hopf bifurcation and evolves to be a time-periodic non-symmetric state. In contrast to channel flow, the linearized operator for such an intricate configuration cannot be derived analytically and must instead be computed numerically. The projection subspace is truncated to r=30r=30 modes. Figure 11: Linear stability analysis results of the 2D cylinder flow system. (a) Eigenspectra of the numerical direct stability analysis (operator-based) and N-based Jacobian for unsteady two-dimensional cylinder flow. (b) Eigenmodes of the unstable eigenvalue from the operator-based and N-based Jacobians. From top to bottom: pressure, streamwise velocity and transverse velocity. Figure 12: Resolvent analysis results of the 2D cylinder flow system. (a) Resolvent gain as a function of frequency in the interval [0, 3]. (b) Optimal pressure forcing and response modes evaluated at the frequency corresponding to the peak gain. The eigenspectra obtained from numerical direct stability analysis and N-based Jacobian are presented in Fig. 11(a), along with the leading eigenmodes shown in Fig. 11(b). As expected, the cylinder flow at Re=100Re=100 exhibits one unstable leading eigenvalue located in the upper half of the complex plane. The leading eigenvalue predicted by the N emulator at frequency f=0.12f=0.12 aligns well with the numerical ground truth, and its eigenmode in the downstream area of the cylinder also shows an excellent agreement. This observation again demonstrates that the N emulator can not only advance the flow field in an autoregressive fashion, but also capture spatial perturbation structures consistent with those in the numerical modal analysis, thereby reflecting the underlying dynamics of the system. Note that the frequency of the leading mode differs from the vortex shedding frequency, i.e., Strouhal number (0.1660.166) at Re=100Re=100. Previous literature (Barkley, 2006) has reported that only linear stability analysis around the mean flow can predict the vortex shedding frequency with good accuracy, but we reiterate that our current method is not applicable to time-averaged states, as this would invalidate the aforementioned derivation. Since the steady-state solution of cylinder flow at Re=100Re=100 represents an artificial pseudo-equilibrium that is not physically realizable, any perturbation introduced around this base state will grow exponentially, ultimately leading to vortex shedding. Consequently, there is no reference for resolvent analysis based on this steady-state solution. In prior resolvent studies of cylinder wakes, the time-averaged flow is typically adopted as the base state, with the unknown nonlinear terms treated as effective forcing to uncover coherent structures. Nevertheless, we can follow the procedure of previous examples to conduct the resolvent analysis around the steady-state flow based on the obtained N-based Jacobian. The resulting gain curves and resolvent pressure modes at the optimal forcing frequency are presented in Fig. 12. Despite the lack of a ground truth solution, the results offer some clear physical insights: the optimal forcing frequency again coincides with the resonance frequency in Fig. 11, and the corresponding response mode forms symmetric lobe-like patterns downstream of the cylinder, similar to the leading eigenmode. The forcing mode also displays lobe-like structures in the downstream region much closer to the cylinder. This is likely associated with the wavemaker region, the area immediately downstream of the cylinder that is most sensitive to external forcing (Giannetti and Luchini, 2007). 3 Discussion In the context of modal analysis, we have demonstrated the promising capabilities of neural networks, as nonlinear function approximators, to serve as a robust and encouraging avenue for extending modal analysis to general nonlinear regimes. We have proposed an algorithm for linear stability analysis and resolvent analysis by leveraging the representational capabilities of N emulators. The method leverages the Jacobian of trained N emulators to approximate the eigenvalues and eigenvectors of the system’s local linearization, and to construct the resolvent operator in this eigenbasis. Its effectiveness is demonstrated on four canonical systems of increasing complexity and nonlinearity, underscoring its potential for analyzing complex dynamics. The fact that such results are achieved purely through the representational capacity of primitive network architectures is highly encouraging. The N emulator and the widely used DMD approach share a close connection, as both are equation‑free methods that rely solely on data rather than explicit governing equations. They can extrapolate future flow evolution from arbitrary initial snapshots at a fraction of the computational cost of traditional solvers, and thus serve as reduced‑order models. Moreover, both approaches are capable of revealing spectral information within their respective applicable regimes, contributing to the broader concept of data‑driven modal analysis. However, the applicability of DMD is limited by its assumption that system evolution can be approximated by a linear mapping between successive snapshots. While this works well for linear cases where DMD modes coincide with those from linear stability analysis, the extracted modes only reflect tangent linear structures rather than the true stability modes for strongly nonlinear cases (Schmid, 2010). By contrast, the N emulator, particularly the time‑stepper model adopted here, is designed to learn nonlinear mappings and therefore possesses broader applicability. In essence, the N emulator can serve as a nonlinear alternative to DMD and effectively replace it in eigenbasis approximation for data‑driven resolvent analysis. As a novel data-driven approach, previous extensions and modifications for traditional methods such as DMD may offer valuable insights for improving our approach. One promising direction lies in preprocessing the data, for instance by projecting it into a reduced-order subspace using autoencoder networks (Wiewel et al., 2019), which can help alleviate challenges associated with increasing dimensionality. Moreover, since DMD can handle non-uniformly sampled datasets to circumvent the Nyquist criterion (Guéniat et al., 2015), our N emulator could also be extended to such data by adopting higher-order integration schemes rather than the direct mapping used in this work. As with other data-driven approaches, N-based modal analysis has the potential to significantly lower the barrier to entry for conducting operator-based modal analysis across diverse fields such as finance (Mann and Kutz, 2016), video processing (Grosek and Kutz, 2014) and neuroscience (Brunton et al., 2016a). N-based modal analysis serves a dual role: First, it acts as a diagnostic tool, revealing physically interpretable structures of a dynamical system. Second, it provides a benchmark for assessing the fidelity of a trained N emulator in capturing the system’s true dynamics, especially when reference modal results are available. This synergy between machine learning and classical modal analysis opens a promising avenue toward interpretable, physics-informed modeling and advances the integration of data-driven and theoretical approaches. 4 Methods 4.1 Linear Stability and resolvent analysis When there is no exogenous forcing applied to the original dynamical system, the governing equation of the perturbation ′ q (Eq. 3) degrades to a homogeneous system: d′dt=′, d q dt=A q , (15) where we can assume the perturbation in the form of normal mode ′=^eσt q = qe^σ t with σ=λ+iωσ=λ+iω. The real part λ and the imaginary part ω are the growth rate and characteristic frequency, respectively, and q describes the mode shape. By substituting the normal mode into Eq. 15, we finally obtain an eigenvalue problem as follows: ^=σ^.A q=σ q. (16) Since the Jacobian matrix A depends on the base state b q_b, it indicates that for a dynamic system with multiple fixed points, they may have different stability properties. It depends on whether each corresponding A has at least one eigenvalue with a positive real part λ. However, linear stability analysis provides only a partial description of the system dynamics because it solely focuses on unstable modes. In practice, a system can remain asymptotically stable with all eigenvalues residing in the stable region, yet still exhibit significant transient growth of perturbations. This motivates the use of such tools as resolvent analysis, in order to obtain a input-output viewpoint of the perturbation evolution. Resolvent analysis aims to examine the response outputs of perturbation ′ q for different forcing inputs and identify the optimal input mode which yields the most amplified energy by internal dynamics. Therefore, we assume the response of perturbations ′ q and the external forcing f to be harmonic, which can be expressed as ′(t)=^e−iωt+c.c. q (t)= qe^-iω t+c.c. and (t)=^e−iωt+c.c. f(t)= fe^-iω t+c.c., where c.c.c.c. represents the complex conjugate and ω is the angular driving frequency. Based on this assumption, Eq. 3 can be rewritten as ^=(−iω−)−1^=(ω)^, q=(-i -A)^-1 f=H(ω) f, (17) where (ω)∈ℂN×NH(ω) ^N× N is defined as the resolvent operator, i.e., transfer function matrix between the forcing inputs and response outputs at frequency ω. Since the goal is to obtain the optimal energy amplification over all possible forcing modes q, a physically meaningful norm for energy measurement needs to be selected. A state variable norm can be defined as ‖^‖2=^∗^|| q||^2_Q= q^*Q q, where ()∗()^* denotes the Hermitian transpose and =∗Q=F^*F is a positive-definite weighting matrix that may account for grid discretization size (Yuan et al., 2023), compressible flow energy (Chu, 1965) and other factors. The Q norm is related to the standard Euclidean norm as ‖^‖2=‖^‖22|| q||^2_Q=||F q||^2_2. The largest input-output gain can be quantified via the ratio of their respective norms σ2(ω)=max^‖′^‖2‖^‖2=max^‖(ω)^‖2‖^‖2=‖(ω)−1‖22.σ^2(ω)= _ f || q ||^2_Q|| f||^2_Q= _ f ||H(ω) f||^2_Q|| f||^2_Q=||FH(ω)F^-1||^2_2. (18) Eq. 18 is an optimization problem, whose solution is given by the weighted SVD of the resolvent operator (ω)−1=(ω)(ω)∗(ω),FH(ω)F^-1= _F(ω) (ω) ^*_F(ω), (19) where ∈ℝn ^n is the singular value matrix whose diagonal entries are the optimal gains associated with the leading n resolvent modes, and the column vectors in −1(ω)F^-1 _F(ω) and −1(ω)F^-1 _F(ω) are the corresponding response and forcing modes. 4.2 Dataset description The dataset of the Lorenz system is constructed from high-fidelity simulations via the fourth-order Runge–Kutta method. 100100 different initial conditions are advanced in time for a total of m=800m=800 snapshots with a fixed 0.010.01 time units. The dataset of the complex Ginzburg-Landau equation is created using the spectral collocation method. To discretize the Ginzburg-Landau equation, we create a domain discretization using Hermite functions (Bagheri et al., 2009). A total of N=220N=220 spectral collocation points are given by the roots of Hn(χx)H_n(χ x), where HnH_n is the n-th Hermite polynomial and χ=(−μ2/(2γ))1/4χ=(- _2/(2γ))^1/4 is a scaling factor. With this choice of N, the discretized domain stretches from x1=−84.99x_1=-84.99 to xN=84.99x_N=84.99. Then 3030 independent simulations are conducted with distinct initial conditions derived from the superposition of multiple random Gaussian profiles. For each simulation, m=100m=100 snapshots are collected at regular time intervals of 0.50.5. The corresponding Hermite differentiation matrices provided by (Weideman and Reddy, 2000) are adopted to formulate the first derivative ∂/∂x∂/∂ x and second derivative ∂2/∂x2∂^2/∂ x^2 to construct the discrete operator matrix A of size N×N× N. The dataset of the channel flow simulations is also created using the spectral method, which consists of snapshots of the velocity perturbation field in two spatial directions collected during the transient evolution stage. The N emulator aims to model the temporal evolution of infinitesimal perturbations around the base parabolic velocity profile U(y)=1−y2U(y)=1-y^2. The computational domain is 2π×22π× 2 in the streamwise and wall-normal directions, where periodicity is imposed along the streamwise direction. The spatial domain is discretized by Ny=64N_y=64 Chebyshev collocation points in the wall-normal direction and Nx=32N_x=32 Fourier collocation points in the streamwise direction. A total of 6060 trajectories are simulated with m=400m=400 snapshots saved every 0.50.5 time units. The Orr–Sommerfeld operator is formulated using Chebyshev differentiation matrices at the same wall-normal collocation points. The dataset of the cylinder flow simulations is generated by an open-source spectral element solver (Fischer et al., 2007), which can also be used for steady-state solution and numerical stability analysis. The cylinder with a diameter D is placed at (0,0)(0,0) in a grid extending from −15D-15D to 35D35D and −15D-15D to 15D15D in x and y directions, respectively. We employ the same boundary conditions as in (Barkley, 2006) for linearized and nonlinear NS equations, respectively. The simulations are initialized with the superposition of the base flow and different infinitesimal perturbations satisfying specific boundary conditions for faster saturation. The dataset comprises 30 trajectories, each sampled for m=200m=200 snapshots every 0.30.3 time units. The flow snapshots are uniformly resampled within [−8,12]×[−3,3][-8,12]×[-3,3] using an 80×2480× 24 Cartesian grid. A masking approach is applied by assigning zero values to the network inputs and outputs inside the cylinder. For each dynamical system, we allocate 80%80\% of the data to the training set, 10%10\% to the validation set, and 10%10\% to the test set. 4.3 Network architecture We employ two types of neural network architectures depending on the dynamical system under consideration. Schematics of the four employed neural network architectures are shown in Supplementary Materials. For the Lorenz system and the complex Ginzburg–Landau (CGL) equation, we use multilayer perceptrons (MLPs) due to the relatively low-dimensional nature of their latent dynamics. Specifically, the Lorenz system is modeled using a three-layer feedforward neural network with 10 neurons per hidden layer, following the design proposed in (Brunton and Kutz, 2022). For the CGL equation, we adopt a deeper, complex-valued MLP consisting of six hidden layers with 512 neurons each, allowing the network to capture the richer spatial–temporal features and complex-valued field interactions. For the channel flow and cylinder wake problems, we employ convolutional neural network (CNN) architectures based on the U-Net design (Ronneberger et al., 2015). In the channel flow case, the model uses four encoding and four decoding stages, with 32 channels at the top layer to balance computational cost and predictive accuracy. The cylinder wake network shares the same U-Net architecture but differs in its input representation, which includes three physical quantities: pressure, streamwise velocity, and transverse velocity—reflecting the multi-field nature of the wake dynamics. 4.4 Rollout Training Since the network is trained by minimizing the difference between predicted and true trajectories, their interplay during the learning process is crucial. A straightforward way of training is one-step training, which only minimizes the one-step prediction error (Tran et al., 2021). By contrast, rollout training by autoregressively applying fθf_θ for multiple time steps can learn long-term dependencies and correct for compounding errors (Um et al., 2020). This can significantly improve the long-term prediction performance, especially when dealing with high-dimensional nonlinear dynamical systems. Suppose the ground truth trajectories are 0(j),1(j),⋯,m(j) q^(j)_0, q^(j)_1,·s, q^(j)_m, where the superscript j∈1,⋯,nj∈ \1,·s,n \ indexes the trajectories corresponding to n distinct initial conditions, and m is the maximum sampling time step per trajectory (Part 1, Fig. 1). Given an arbitrary state i(j) q^(j)_i extracted from the ground truth datasets, the trajectories predicted by iteratively applying the learned model fθf_θ for t steps are denoted as i(j),fθ1(i(j)),⋯,fθt(i(j)) q^(j)_i,f^1_θ( q^(j)_i),·s,f^t_θ( q^(j)_i). The total loss function for training with t-step rollouts averaged on n trajectories is: ℒ(θ)=1n∑j=1n∑i=0m−t∑k=1t[ℒ2(fθk(i(j)),i+k(j))],L(θ)= 1n _j=1^n _i=0^m-t _k=1^t [L_2 (f_θ^k( q^(j)_i), q^(j)_i+k ) ], (20) where ℒ2L_2 represents L2-norm loss function. Full gradients can be computed owing to the differentiability of the N emulator, which enables backpropagation through t rollout steps and subsequent updates of the network weights θ (Part 2, Fig. 1). For the Lorenz and complex Ginzburg–Landau systems, the relative simplicity of the dynamics allows the use of a one-step rollout training strategy. In contrast, the channel flow problem is trained with a ten-step rollout to better capture the transient evolution. For the cylinder wake flow, a twenty-step unrolling strategy is employed to ensure coverage of a full period of vortex shedding. 4.5 Ensemble-Based Jacobian Estimation To ensure a robust Jacobian estimation, we adopt a model ensemble strategy (Chua et al., 2018; Allen et al., 2022) by averaging the Jacobians predicted from multiple trained models initialized with different random seeds. This can not only mitigate the effects of weight initialization diversity but also lead to a more representative average gradient. Five independent networks are trained for each case unless otherwise specified. Data Availability The source code for the current study is available at Github (https://github.com/tum-pbs/NonlinearRA). Data available upon request from chengyun.wang@tum.de. Previously published platform (Koehler et al., 2024) used for this work is available at GitHub (https://github.com/tum-pbs/apebench). Acknowledgements C.W. is supported by the China Scholarship Council (No. 202406340048) and L.C. is supported by the National Natural Science Foundation of China (Grant No. 92470120). We are grateful for many fruitful discussions with Jean-Christophe Loiseau and Eduardo Martini. We would also like to thank Benjamin Herrmann for kindly providing the dataset used in his paper. Competing Interests The authors declare no competing interests. References [1] K. Allen, T. Lopez-Guevara, K. L. Stachenfeld, A. Sanchez Gonzalez, P. Battaglia, J. B. Hamrick, and T. Pfaff (2022) Inverse design for fluid-structure interactions using graph network simulators. Advances in Neural Information Processing Systems 35, p. 13759–13774. Cited by: §4.5. [2] P. J. Baddoo, B. Herrmann, B. J. McKeon, J. Nathan Kutz, and S. L. Brunton (2023) Physics-informed dynamic mode decomposition. Proceedings of the Royal Society A 479 (2271), p. 20220576. Cited by: §1. [3] S. Bagheri, D. Henningson, J. Hœpffner, and P. J. Schmid (2009) Input-output analysis and control design applied to a linear model of spatially developing flows. Applied Mechanics Reviews 62 (2), p. 020803. Cited by: §2, §2, §2, §4.2. [4] Y. Bar-Sinai, S. Hoyer, J. Hickey, and M. P. Brenner (2019) Learning data-driven discretizations for partial differential equations. Proceedings of the National Academy of Sciences 116 (31), p. 15344–15349. Cited by: §1. [5] D. Barkley (2006) Linear analysis of the cylinder wake mean flow. Europhysics Letters 75 (5), p. 750. Cited by: §2, §4.2. [6] S. Battat, D. A. Weitz, and G. M. Whitesides (2022) Nonlinear phenomena in microfluidics. Chemical Reviews 122 (7), p. 6921–6937. Cited by: §1. [7] A. L. Bertozzi and M. P. Brenner (1997) Linear stability and transient growth in driven contact lines. Physics of Fluids 9 (3), p. 530–539. Cited by: §1. [8] B. W. Brunton, L. A. Johnson, J. G. Ojemann, and J. N. Kutz (2016) Extracting spatial–temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition. Journal of neuroscience methods 258, p. 1–15. Cited by: §3. [9] S. L. Brunton and J. N. Kutz (2022) Data-driven science and engineering: machine learning, dynamical systems, and control. Cambridge University Press. Cited by: §2, §4.3. [10] S. L. Brunton, J. L. Proctor, and J. N. Kutz (2016) Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences 113 (15), p. 3932–3937. Cited by: §1. [11] L. Chen and N. Thuerey (2023) Towards high-accuracy deep learning inference of compressible flows over aerofoils. Computers & Fluids 250, p. 105707. Cited by: §1. [12] L. Chen and N. Thuerey (2024) Deep learning-based predictive modeling of transonic flow over an airfoil. Physics of Fluids 36 (12). Cited by: §1. [13] B. Chu (1965) On the energy transfer to small disturbances in fluid flow (part i). Acta Mechanica 1 (3), p. 215–234. Cited by: §4.1. [14] K. Chua, R. Calandra, R. McAllister, and S. Levine (2018) Deep reinforcement learning in a handful of trials using probabilistic dynamics models. Advances in neural information processing systems 31. Cited by: §4.5. [15] J. Colen, M. Han, R. Zhang, S. A. Redford, L. M. Lemma, L. Morgan, P. V. Ruijgrok, R. Adkins, Z. Bryant, Z. Dogic, et al. (2021) Machine learning active-nematic hydrodynamics. Proceedings of the National Academy of Sciences 118 (10), p. e2016708118. Cited by: §1. [16] P. Fischer, J. Lottes, and H. Tufo (2007) Nek5000. Technical report Argonne National Laboratory (ANL), Argonne, IL (United States). Cited by: §4.2. [17] D. Floryan and M. D. Graham (2022) Data-driven discovery of intrinsic dynamics. Nature Machine Intelligence 4 (12), p. 1113–1120. Cited by: §1. [18] F. Giannetti and P. Luchini (2007) Structural sensitivity of the first instability of the cylinder wake. Journal of Fluid Mechanics 581, p. 167–197. Cited by: §2. [19] J. Grosek and J. N. Kutz (2014) Dynamic mode decomposition for real-time background/foreground separation in video. arXiv preprint arXiv:1404.7592. Cited by: §3. [20] F. Guéniat, L. Mathelin, and L. R. Pastur (2015) A dynamic mode decomposition approach for large and arbitrarily sampled systems. Physics of Fluids 27 (2). Cited by: §3. [21] B. Herrmann, P. J. Baddoo, R. Semaan, S. L. Brunton, and B. J. McKeon (2021) Data-driven resolvent analysis. Journal of Fluid Mechanics 918, p. A10. Cited by: §1, §2, §2, §2, §2, §2. [22] P. Holmes (2012) Turbulence, coherent structures, dynamical systems and symmetry. Cambridge university press. Cited by: §1, §1. [23] P. Huerre and P. A. Monkewitz (1985) Absolute and convective instabilities in free shear layers. Journal of Fluid Mechanics 159, p. 151–168. Cited by: §2. [24] R. M. Jendrejack, E. T. Dimalanta, D. C. Schwartz, M. D. Graham, and J. J. de Pablo (2003) DNA dynamics in a microchannel. Physical review letters 91 (3), p. 038102. Cited by: §1. [25] D. Kochkov, J. A. Smith, A. Alieva, Q. Wang, M. P. Brenner, and S. Hoyer (2021) Machine learning–accelerated computational fluid dynamics. Proceedings of the National Academy of Sciences 118 (21), p. e2101784118. Cited by: §1. [26] F. Koehler, S. Niedermayr, N. Thuerey, et al. (2024) APEBench: a benchmark for autoregressive neural emulators of pdes. Advances in Neural Information Processing Systems 37, p. 120252–120310. Cited by: Data Availability. [27] T. Kolokolnikov, H. Sun, D. Uminsky, and A. L. Bertozzi (2011-07) Stability of ring patterns arising from two-dimensional particle interactions. Phys. Rev. E 84, p. 015203. External Links: Document Cited by: §1. [28] D. Kyle and K. Sreenivasan (1993) The instability and breakdown of a round variable-density jet. Journal of Fluid Mechanics 249, p. 619–664. Cited by: §2. [29] B. List, L. Chen, and N. Thuerey (2022) Learned turbulence modelling with differentiable fluid solvers: physics-based loss functions and optimisation horizons. Journal of Fluid Mechanics 949, p. A25. Cited by: §1. [30] M. Luhar, A. S. Sharma, and B. J. McKeon (2014) Opposition control within the resolvent analysis framework. Journal of Fluid Mechanics 749, p. 597–626. Cited by: §1. [31] L. M. Mack (1976) A numerical study of the temporal eigenvalue spectrum of the blasius boundary layer. Journal of Fluid Mechanics 73 (3), p. 497–520. Cited by: §2. [32] J. Mann and J. N. Kutz (2016) Dynamic mode decomposition for financial trading strategies. Quantitative Finance 16 (11), p. 1643–1655. Cited by: §3. [33] E. Martini, D. Rodríguez, A. Towne, and A. V. Cavalieri (2021) Efficient computation of global resolvent modes. Journal of Fluid Mechanics 919, p. A3. Cited by: §2. [34] B. J. McKeon and A. S. Sharma (2010) A critical-layer framework for turbulent pipe flow. Journal of Fluid Mechanics 658, p. 336–382. Cited by: §1. [35] J. Morton, A. Jameson, M. J. Kochenderfer, and F. Witherden (2018) Deep dynamical modeling and control of unsteady fluid flows. In Advances in Neural Information Processing Systems, Cited by: §1. [36] J. Page, P. Norgaard, M. P. Brenner, and R. R. Kerswell (2024) Recurrent flow patterns as a basis for two-dimensional turbulence: predicting statistics from structures. Proceedings of the National Academy of Sciences 121 (23), p. e2320007121. Cited by: §1. [37] J. L. Proctor, S. L. Brunton, and J. N. Kutz (2016) Dynamic mode decomposition with control. SIAM Journal on Applied Dynamical Systems 15 (1), p. 142–161. Cited by: §1. [38] M. Raissi, P. Perdikaris, and G. E. Karniadakis (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics 378, p. 686–707. Cited by: §1. [39] S. C. Reddy and D. S. Henningson (1993) Energy growth in viscous channel flows. Journal of Fluid Mechanics 252, p. 209–238. Cited by: §2. [40] L. V. Rolandi, J. H. M. Ribeiro, C. Yeh, and K. Taira (2024) An invitation to resolvent analysis. Theoretical and Computational Fluid Dynamics 38 (5), p. 603–639. Cited by: §1. [41] O. Ronneberger, P. Fischer, and T. Brox (2015) U-net: convolutional networks for biomedical image segmentation. In International Conference on Medical image computing and computer-assisted intervention, p. 234–241. Cited by: §4.3. [42] C. W. Rowley and S. T. Dawson (2017) Model reduction for flow analysis and control. Annual Review of Fluid Mechanics 49 (1), p. 387–417. Cited by: §1. [43] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S. Henningson (2009) Spectral analysis of nonlinear flows. Journal of fluid mechanics 641, p. 115–127. Cited by: §1. [44] P. J. Schmid and D. S. Henningson (2012) Stability and transition in shear flows. Vol. 142, Springer Science & Business Media. Cited by: §2, §2, §2. [45] P. J. Schmid (2007) Nonmodal stability theory. Annu. Rev. Fluid Mech. 39 (1), p. 129–162. Cited by: §1. [46] P. J. Schmid (2010) Dynamic mode decomposition of numerical and experimental data. Journal of fluid mechanics 656, p. 5–28. Cited by: §1, §1, §3. [47] S. C. Shadden, F. Lekien, and J. E. Marsden (2005) Definition and properties of lagrangian coherent structures from finite-time lyapunov exponents in two-dimensional aperiodic flows. Physica D: Nonlinear Phenomena 212 (3-4), p. 271–304. Cited by: §2. [48] J. D. Simon and S. K. Mitter (1968) A theory of modal control. Information and control 13 (4), p. 316–353. Cited by: §1. [49] K. Sreenivasan, S. Raghu, and D. Kyle (1989) Absolute instability in variable density round jets. Experiments in Fluids 7 (5), p. 309–317. Cited by: §2. [50] L. Sun and R. Wu (2015) Mapping complex traits as a dynamic system. Physics of life reviews 13, p. 155–185. Cited by: §1. [51] Y. Sun, Q. Liu, L. N. Cattafesta I, L. S. Ukeiley, and K. Taira (2020) Resolvent analysis of compressible laminar and turbulent cavity flows. AIAA journal 58 (3), p. 1046–1055. Cited by: §1. [52] S. Symon, K. Rosenberg, S. T. Dawson, and B. J. McKeon (2018) Non-normality and classification of amplification mechanisms in stability and resolvent analysis. Physical Review Fluids 3 (5), p. 053902. Cited by: §1, §2. [53] K. Taira, S. L. Brunton, S. T. Dawson, C. W. Rowley, T. Colonius, B. J. McKeon, O. T. Schmidt, S. Gordeyev, V. Theofilis, and L. S. Ukeiley (2017) Modal analysis of fluid flows: an overview. AIAA journal 55 (12), p. 4013–4041. Cited by: §1. [54] V. Theofilis (2011) Global linear instability. Annual Review of Fluid Mechanics 43 (1), p. 319–352. Cited by: §1. [55] N. Thuerey, K. Weißenow, L. Prantl, and X. Hu (2020) Deep learning methods for reynolds-averaged navier–stokes simulations of airfoil flows. AIAA Journal 58 (1), p. 25–36. Cited by: §1. [56] A. Tran, A. Mathews, L. Xie, and C. S. Ong (2021) Factorized fourier neural operators. arXiv preprint arXiv:2111.13802. Cited by: §4.4. [57] L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A. Driscoll (1993) Hydrodynamic stability without eigenvalues. Science 261 (5121), p. 578–584. Cited by: §1. [58] K. Um, R. Brand, Y. R. Fei, P. Holl, and N. Thuerey (2020) Solver-in-the-loop: learning from differentiable physics to interact with iterative pde-solvers. Advances in neural information processing systems 33, p. 6111–6122. Cited by: §4.4. [59] A. S. Utada, A. Fernandez-Nieves, J. M. Gordillo, and D. A. Weitz (2008) Absolute instability of a liquid jet in a coflowing stream. Physical review letters 100 (1), p. 014502. Cited by: §2. [60] J. A. Weideman and S. C. Reddy (2000) A matlab differentiation matrix suite. ACM transactions on mathematical software (TOMS) 26 (4), p. 465–519. Cited by: §4.2. [61] S. Wiewel, M. Becher, and N. Thuerey (2019) Latent space physics: towards learning the temporal evolution of fluid flow. In Computer graphics forum, Vol. 38, p. 71–82. Cited by: §3. [62] C. Yeh and K. Taira (2019) Resolvent-analysis-based design of airfoil separation control. Journal of Fluid Mechanics 867, p. 572–610. Cited by: §1. [63] R. Yu and R. Wang (2024) Learning dynamical systems from data: an introduction to physics-guided deep learning. Proceedings of the National Academy of Sciences 121 (27), p. e2311808121. Cited by: §1. [64] H. Yuan, J. Kou, C. Gao, and W. Zhang (2023) Resolvent and dynamic mode analysis of flow past a square cylinder at subcritical reynolds numbers. Physics of Fluids 35 (7). Cited by: §4.1.