Long-Time Trajectory Approximation via SA-NODEs: Model Predictive and Floquet Strategies

LONG-TIME TRAJECTORY APPROXIMATION VIA SA-NODES: MODEL PREDICTIVE AND FLOQUET STRATEGIES

arXiv:2608.10738v1 [cs.LG] 11 Aug 2026

ZIQIAN LI1 AND NIKOLAOS M. MATZAKOS21 Abstract. We study the approximation of dynamical systems by semi-autonomous neural ordinary differential equations (SA-NODEs) over long time horizons. For a single network trained on the whole horizon, the available error bound deteriorates double exponentially in the horizon length. We develop two training strategies that avoid this barrier, each built on a reset of the state. The model predictive strategy partitions the horizon adaptively and restarts every window from observed data: when training meets a prescribed tolerance on every window, the composite model meets it uniformly in time, with a parameter budget linear in the horizon for targets with a bounded, uniformly regular reachable tube. The Floquet strategy addresses autonomous targets with a stable limit cycle and uses no data at deployment: a certified contraction of the learned return map confines the error to linear growth in the number of elapsed periods. For the time-periodic architecture we deploy, the scalar certificate degenerates; we prove instead a uniform-in-time orbital guarantee whose hypotheses are measured on the trained model, and an obstruction showing that, for an exactly periodic learned field, small one-period error and a contracting stroboscopic map cannot hold at once. Numerical experiments on four benchmarks confirm the predicted error laws and measure the hypotheses of every guarantee.

1. Introduction 1.1. Problem statement and motivation. A data-driven model of a dynamical system is useful only over the time horizon on which its predictions stay accurate. Neural ordinary differential equations (neural ODEs) parameterize the right-hand side of a differential equation by a neural network, viewing a residual architecture as the discretization of a continuous-time flow [7]. When the training data are sampled from trajectories of an underlying dynamical system ẋ = f (t, x), the trained network approximates the flow itself. This paper quantifies how the error of such models grows with the horizon and develops two training strategies that control the growth. The architecture we work with is the semi-autonomous neural ODE (SA-NODE) [16, eq. (2.1)]. For a horizon T > 0 and an initial condition x0 ∈ Rd , with d the dimension of the state space, it evolves a state x(t) ∈ Rd by  P X 
  ẋ(t) = Wi ◦ σ A1i x(t) + A2i t + Bi , t ∈ (0, T ], (1.1) i=1    x(0) = x0 . We write fΘ (x, t) for the right-hand side of (1.1). Its defining feature, and the reason for the name, is that the time t enters only through the linear bias A2i t + Bi . The trainable parameters Θ = (Wi , A1i , A2i , Bi )Pi=1 are time-independent, with Wi , A2i , Bi ∈ Rd and A1i ∈ Rd×d ; ◦ is the componentwise product, σ the componentwise activation, and P ∈ N the hidden width. On a 2020 Mathematics Subject Classification. 68T07, 93B45, 34C25, 65L05, 41A30. Key words and phrases. Neural ODEs, long-time approximation, model predictive control, Floquet theory, limit cycles. 1 Chair for Dynamics, Control, Machine Learning and Numerics, Department of Mathematics, FriedrichAlexander-Universität Erlangen-Nürnberg, 91058 Erlangen, Germany (ziqian.li@fau.de). 2 School of Pedagogical & Technological Education (ASPETE), Marousi Attikis, 151 22 Athens, Greece (nikmatz@aspete.gr). 1

2

Z. LI AND N. M. MATZAKOS

compact set K ⊂ Rd of initial conditions, the universal approximation theorem [16, Thm. 2.3] provides parameters ΘP with sup x0 ∈K, t∈[0,T ]

b Θ (t; x0 ) Φ(t; x0 ) − Φ P

≤ CT,K,f P −1/2 ,

(1.2)

b Θ is the SA-NODE flow, and the subscripts of the constant record its where Φ is the target flow, Φ P dependence on the horizon, the set of initial conditions, and the field. The setting is approximation, not identification. Estimate (1.2) provides short-horizon closeness of the learned flow. We take that closeness as given and study how the error propagates over long horizons. The experiments train the models from sampled trajectories, but we do not address system identification from finite or noisy data. The question this paper addresses is how the constant in (1.2) grows with T , and how to avoid that growth. The available upper bound on that constant carries the Grönwall factor eLT , with L the Lipschitz constant of f in x, and the explicit constant of [16, Remark 2.6] deteriorates double exponentially in T ; Section 2 states it and shows that both mechanisms behind it are individually sharp. We call this the time horizon barrier. It is a property of the available certificate: no matching lower bound over Lipschitz targets is known. According to this guarantee, a single network trained once on the whole horizon, which we call the monolithic model, is useful only for horizons of a few multiples of 1/L. Whether the barrier is intrinsic is open; Lemma 2.2 reduces the question to a width lower bound for the target field, stated as Problem 6.1. The two strategies below avoid the barrier by changing the learning problem rather than enlarging the network. The barrier was identified in [16, Remarks 2.6 and 2.9], together with a proposal, not developed there, to control it by model predictive ideas. The factor eLT results from integrating error over the whole of [0, T ] without any correction. Structured targets can escape it, but through learning strategies that exploit that structure. In computational terms, we replace one network fitted to the whole horizon by several short-horizon networks restarted from data, or by a single periodic network with a certified stable orbit. If the true state is available at intermediate times, the horizon can be partitioned into windows of length at most τmax ≪ T , one network per window, with the initial condition reset to the true state at each switch time; the Grönwall factor then only reaches eLτmax . This is model predictive control (MPC) applied to approximation rather than control, and we call the injection of the true state a data reset. The name records the mechanics: adaptive multiple shooting with the node values pinned to data, not an online feedback law. If instead the target is autonomous and its trajectories are attracted to a stable limit cycle, the orbit itself corrects errors: every perturbation transverse to the cycle is contracted at each return to a Poincaré section, at the rate given by the Floquet multiplier. A learned model that inherits this contraction needs no data at deployment. We call this a dynamical reset. The second half of this paper asks whether a trained SA-NODE can be made to inherit it, with a certificate. The two resets are complementary. The data reset restarts the error at every switch time; it applies to any Lipschitz target with a bounded, uniformly regular reachable tube, at the cost of the true state at each switch. The dynamical reset contracts transverse errors at every return and needs no data at deployment, but it requires an autonomous target settling onto a limit cycle, and the contraction of the learned return map must be trained in and then certified. Targets of the latter kind arise in chemical oscillators, vortex shedding, structural vibrations, and circadian or cardiac rhythms, where the horizon of interest spans many periods. The main purpose of this paper is to turn each reset into a quantitative long-horizon approximation guarantee for trained SA-NODEs, and to measure the hypotheses of these guarantees on the trained models themselves.

LONG-TIME APPROXIMATION VIA SA-NODES

3

1.2. Main results. The main results are as follows. In the Floquet statements below, and throughout Section 4, states are written unbold, x0 , as points of Rd , and the parameter vector Θ is abbreviated Θ. • Composite data-assisted approximation. We design an adaptive receding-horizon algorithm (Algorithm 3.1) that grows a composite SA-NODE one window at a time, opening a new window when the supervised error first exceeds a tolerance ε. The guarantee is conditional on training: the tolerance must be realized on every window and the run must cover the horizon (Assumption 3.4). In data-IC mode, in which the true state is supplied as the b then satisfies initial condition at each of the N switch times, the composite flow Φ b x0 ) Φ(t; x0 ) − Φ(t;

sup

≤ ε,

t∈[0,T ], x0 ∈K

Ptotal ≲ N Cτ2max ,K∞ ,f ε−2 ,

(1.3)

where the window constant is controlled by the data of f on the closure K∞ of the forward reachable tube; it is independent of T whenever that tube is bounded and its time-shifted Sobolev data are uniformly bounded (Assumption 3.3, Theorem 3.5). Every admissible partition has N ≥ ⌈T /τmax ⌉; the uniform partition attains it, so the budget is then linear in T . For the adaptive partition, the minimum window length of Algorithm 3.1 already forces N = O(T ), with a constant fixed by the resolution of the training data rather than by the target. Whether linear growth survives the removal of that safeguard is open (Problem 6.2), and the measurements of Section 5 are consistent with it. The reset confines the Grönwall factor of (2.1) to a single window; the growth of the reachable tube persists through the domain of the window constant. • Deployment on novel initial conditions. Without state resets the per-window errors chain multiplicatively. We call predicted-IC the deployment in which each window starts from the previous window’s prediction instead of the true state. For the predicted-IC composite flow and every x0 ∈ K we prove the partition-adapted bound sup t∈[0,T ]

b pred (t; x0 ) Φ(t; x0 ) − Φ

≤ ε

N X k=1

eL̄(T −τk ) ≤ ε

eN L̄τmax − 1 , eL̄τmax − 1

(1.4)

where τ1 < · · · < τN are the switch times and L̄ bounds the Lipschitz constants of the learned fields (Proposition 3.7). The bound remains exponential; the gain over the monolithic model lies in the prefactor. • Certified Floquet approximation. For an autonomous target with a hyperbolic stable limit cycle Γ of period T̂ , we train a single SA-NODE with a differentiable Floquet loss that penalizes a surrogate ρ̃T (Θ) of the learned Poincaré multiplier; the subscript of the surrogate marks the transverse direction, not the horizon. Let PΘ denote the first-return map of the learned flow to a section through a base point p of the cycle, and let ρ∗ < 1 be a prescribed contraction threshold. If the certificate ρ(DPΘ (p)) ≤ ρ∗ holds for the autonomous return map, then for every initial condition near Γ and every t ≥ 0,  t ∥ΦΘ (t; x0 ) − Φf (t; x0 )∥ ≤ C ε 1 + , (1.5) T̂ where ε now denotes the one-period accuracy of the learned flow, the analogue for this strategy of the tolerance above (Theorem 4.13). The bound is linear in the number of elapsed periods, with C independent of t. Transverse errors are contracted at every return, and only the phase drift between the two clocks persists. The periodically encoded architecture reduces to the autonomous case by Lemma 4.15, at the cost of adding to ε the measurable oscillation η of the learned field about its time average. That reduction is effective only while η is small, and the trained models of Section 5.4 are not in that regime; this caveat governs how the experiments are read.

4

Z. LI AND N. M. MATZAKOS

• Certification routes. The contraction hypothesis is verifiable. In dimension two the identity of Liouville–Abel ties the surrogate to the autonomous return-map multiplier: exactly when the trajectory loss vanishes, and to O(ε) in general (Corollary 4.18). Writing LF for the Floquet loss built on the surrogate and ρΘ for the spectral radius of the learned return map at its fixed point, for ε small enough and a relaxation returning a closed orbit, as quantified in Corollary 4.18, the chain LF (Θ) = 0 =⇒ ρΘ < 1 =⇒

e(k T̂ ) ≤ Cε(1 + k),

(1.6)

where e(t) := ΦΘ (t; x0 ) − Φf (t; x0 ) is the trajectory error and k counts the elapsed periods, links the training objective directly to the long-horizon guarantee (Corollary 4.18; see also Proposition 4.16). This chain is proved for an autonomous planar field. For the deployed time-periodic architecture we prove a separate orbital guarantee (Theorem 4.22) and show that the trajectory-wise one is quantitatively incompatible with the certified contraction. One-period closeness and a contracting stroboscopic map are quantitatively incompatible, under a contraction hypothesis on a neighborhood of the cycle (Proposition 4.20) or under the local spectral hypothesis the experiments measure (Proposition 4.21). For the encoded architecture the divergence integral controls only the determinant of the monodromy, so we compute the full spectrum instead (Remark 4.19). In higher dimension certification proceeds through C 1 -closeness and a Bauer–Fike margin. • Numerical experiments. Section 5 illustrates each result with a dedicated experiment. On a forced Duffing benchmark the composite meets its tolerance up to a quantified safeguard excess, while a monolithic baseline under matched training cost exceeds it about 4.0-fold. On the pendulum the monolithic error grows exponentially along the horizon, at a measured rate e0.39 t , and the window count grows with T at fixed ε, consistent with a linear law. A two-by-two ablation on two limit-cycle benchmarks measures the spectral radius of the learned one-period monodromy, together with the two hypotheses of Theorem 4.22 and the lower bound that Proposition 4.20 predicts for the one-period error. An autonomous arm shows the same protocol failing without the encoding: no cold-started run produces a certifiable closed orbit. A warm-started autonomous arm repairs the failure: fine-tuning from the converged encoded models yields closed, transversally stable learned orbits whose surrogate–multiplier gap is 65 to 320 times below the measured ε, consistent with Corollary 4.18(ii). A direct comparison on the van der Pol oscillator places the two reset mechanisms side by side under one protocol. 1.3. Novelties. The novelties of this paper concern two audiences: numerical analysis, and the practice of learning dynamical models from data. For numerical analysis, the paper turns two classical mechanisms into approximation guarantees for learned flows. To our knowledge, the MPC strategy gives the first uniform-in-time, datarestarted approximation bound for SA-NODEs with an explicit summed-width budget. In multipleshooting practice the partition serves identification, control, or discretization-error control, not an approximation guarantee for a learned flow. The Floquet strategy establishes a propagation law for learned periodic dynamics: a certified contraction of the trained return map converts one-period closeness into the global-in-time linear bound of Theorem 4.13. The certificate is exact in the plane for an autonomous learned field up to an O(ε) period correction (Corollary 4.18) and perturbative in higher dimension (Proposition 4.16). For the time-periodic architecture we deploy, neither form is available; there the full monodromy spectrum supplies the measured hypothesis of the orbital bound of Theorem 4.22. For machine learning and data science, the results change what a long-horizon model of dynamics should be asked to demonstrate. A small loss on short trajectories does not by itself guarantee long-horizon accuracy: long-time accuracy is governed by the reset mechanism, not by the fit. Data

LONG-TIME APPROXIMATION VIA SA-NODES

5

play two roles: they train the model, and they reset it at deployment. The uniform guarantee of the MPC strategy rests on the second. Without deployment data, stability must be trained in and then verified. The Floquet loss makes a differentiable surrogate of the certificate trainable. For the time-periodic architecture that surrogate is the determinant of the monodromy, which does not control its spectral radius, so what the spectrum then supports is not (4.5) but the orbital guarantee of Theorem 4.22 (Section 4.4). 1.4. Related works. Our work draws on four lines of research: neural ODEs, multiple shooting and receding-horizon methods, stability-regularized learning of dynamics, and the error-growth theory of numerical integration near periodic orbits. Neural ODEs were introduced in [7], who interpret a residual network as the discretization of a continuous-time flow and train the vector field by differentiating through the solver. Beyond supervised learning, they have become a model class for dynamics, and extended to operator learning [17]. The architecture of the present paper is the SA-NODE of [16]. Its universal approximation theorem, recalled as Theorem 2.1, supplies the certificate whose horizon dependence is the object of this paper. Section 3 turns the model predictive remark of [16, Rem. 2.9] into an algorithm and a guarantee (Algorithm 3.1, Theorem 3.5). Partitioning a horizon and training one model per window, with the true state injected at the switch times, has the same mechanics as multiple shooting for trajectory fitting, with the node values pinned to the data. This is the limiting case of the classical multiple-shooting method of Bock for parameter identification [2] and optimal control [3], in which the node values are decision variables coupled by matching constraints. Its neural-ODE incarnations [26] and segment-wise training are by now common in scientific machine learning. Classical MPC theory is surveyed in [5]; boundary-value and shooting methods in [14]; receding-horizon numerics for parameter identification appear in [27], and turnpike phenomena for neural ODEs in [11, 9, 24, 18]. On the Floquet side, stability-regularized neural ODEs have been studied through Lyapunov functions [23, 13, 15], which target equilibria rather than periodic orbits. Limit-cycle learning appears in [22, 28], orbital stability has been trained through transverse contraction criteria [29], and limit cycles with prescribed trajectories and phase response have been designed through Floquet theory directly [21]. These embed stability into training or design in various forms; what is new here is not a stability-constrained objective but the propagation law that a certified return map yields, together with the obstruction that separates orbital from trajectory accuracy for the periodic architecture. Closest to our infinite-horizon question, [25] proves universal approximation of Morse–Smale systems on [0, ∞) in an ε-δ sense that excludes an initial-condition set of small measure; that result is existential, whereas our bounds are conditional on quantities measured on the trained model. On the multiple-shooting side, differentiable shooting layers [19] use the same partition mechanics for parallelism rather than for an approximation guarantee. Matzakos [20] proves that activation saturation obstructs the transverse contraction required by dissipative limit cycles in standard neural ODEs; our architecture counters it with a periodic time encoding, and Section 5.4 measures how much each ingredient contributes. The linear-in-k structure of Theorem 4.13 parallels the error-growth theory of Cano and SanzSerna [6] for one-step integration of attracting periodic orbits; see also [4, 12]. The parallel is structural only. In that theory the per-period error is a discretization error that vanishes with the step size, and the contraction of the numerical return map is inherited from a convergent scheme applied to the true field. Here the field is learned: the per-period error is an approximation error that no step refinement removes, and the return map of the learned flow need not contract at all. Contraction is therefore not automatic. It must be certified on the trained model, or trained in directly by the Floquet loss (Assumption 4.6, Section 4.3). As in that theory, the learned system has its own periodic orbit, its own return time, and a neutral phase direction that no transverse contraction controls. Separating the contracted transverse error from the accumulating phase

6

Z. LI AND N. M. MATZAKOS

mismatch is the technical core of Theorem 4.13. Floquet theory itself is classical [10]; a modern account is [8, §2.4]. 1.5. Organization. Section 2 fixes the SA-NODE class and quantifies the time horizon barrier. Section 3 develops the model predictive strategy: the composite flow, the adaptive algorithm, and the error analysis. Section 4 develops the Floquet strategy: the geometric setting, the main theorem, and the two certification routes. Section 5 reports the numerical experiments, one per theoretical claim. Section 6 concludes and states three open problems. The proofs are deferred to Appendix A, one part per result: Appendix A.1 proves Lemma 2.2, Appendix A.2 the results of Section 3, and Appendices A.3–A.7 the results of Section 4. 2. Preliminaries: the SA-NODE class and the time horizon barrier This section recalls the SA-NODE class and its approximation theorem, and isolates how the approximation constant grows with the horizon. Section 2.1 states the theorem together with the explicit constant (2.1). Section 2.2 shows that the two mechanisms behind that constant are sharp. They are the obstruction that Sections 3 and 4 remove. 2.1. The SA-NODE class and its approximation theorem. We recall the approximation theorem whose constant is the object of the paper, in the notation used throughout. Fix a compact K ⊂ Rd , a horizon T > 0, a hidden width P ∈ N, and a non-polynomial activation σ ∈ C 0 (R) applied componentwise as σ. The SA-NODE (1.1) then has P d(d + 3) trainable parameters; we b Θ (t; x0 ) for its flow, the solution x(t) of (1.1) from x(0) = x0 , and write Φ KT := {Φ(t; x0 ) : x0 ∈ K, t ∈ [0, T ]} for the reachable tube of the target flow. The following is Theorem 2.3 of [16], reproduced for completeness; the constant is the object of interest. Theorem 2.1 (universal approximation, [16, Theorem 2.3]). Let f ∈ C 0 (Rd × [0, T ]) be Lipschitz s (Rd × [0, T ]) with s > (d + 1)/2 + 2; let σ in x uniformly in t, with constant L, and belong to Hloc be the ReLU activation. For every P ≥ 3 there exist parameters ΘP such that (1.2) holds with a constant CT,K,f independent of P and of x0 . The smoothness threshold s > (d + 1)/2 + 2 is the C 2 -embedding threshold in d + 1 variables, matching the fact that (1.1) is a shallow network in the joint input (x, t). The activation is fixed to ReLU because the uniform rate rests on the L∞ approximation estimate for ReLU networks. For activations that are twice weakly differentiable and satisfy Z σ ′′ (x) (1 + |x|) dx < ∞, R

the sigmoid among them, the same rate holds with the supremum over initial conditions weakened to an L2 average in x0 [16, Remark 3.7]. The horizon enters through the constant, which [16, Remark 2.6] makes explicit in the nord/2+3 malization K = [−1, 1]d and T ≥ 1, for f uniformly L-Lipschitz in x and in Hloc (Rd+1 ; Rd ). The spacetime reachable set is then contained in the box ΩL,T := [−T eLT , T eLT ]d+1 , under the normalization on ∥f (0, ·)∥ implicit in the growth lemma of [16, §3.3], and   √ 3 CT,K,f ≤ Cd T ∥f ∥∗ exp 52 LT + d L + Cd e 2 LT ∥f ∥∗ , ∥f ∥∗ := ∥f ∥H d/2+3 (ΩL,T ) , (2.1) with Cd depending only on d. This bound is double exponential in T : the Sobolev norm is taken on a domain that already grows like eLT , and it re-enters the Grönwall exponent through the Lipschitz constant of the learned field.

LONG-TIME APPROXIMATION VIA SA-NODES

7

2.2. Sharpness of the factor eLT . Two mechanisms drive the growth in (2.1), and both are already visible at the level of the flow. We isolate them here because the strategies of Sections 3 and 4 act on the first, which is Grönwall amplification: if ∥fΘ − f ∥ ≤ δ on a neighborhood of KT that contains both flows on [0, T ], with L the Lipschitz constant of f there, then for every x0 ∈ K sup

b Θ (t; x0 ) − Φ(t; x0 ) Φ

≤ δ

t∈[0,T ]

eLT − 1 ≤ δ T eLT , L

(2.2)

and the first inequality is an equality for ẋ = Lx with a constant field perturbation, so the exponential bound is attained. For L = 0 the quotient is read as its limit T . The second mechanism is the growth of the reachable tube: two solutions started diam(K) apart may separate by the factor eLT , and bounded forcing contributes a term of the same exponential order, so
diam(KT ) ≤ diam(K) + C eLT (2.3) for fields with supt ∥f (t, 0)∥ < ∞, with C depending on this bound, on L, and on supx∈K ∥x∥; this order is attained for ẋ = Lx. Neither mechanism is an artifact of a loose estimate, since each is attained on that example. Attainment concerns the two steps of the bound, not the approximation problem: on ẋ = Lx the class (1.1) is exact, since a ReLU network represents a linear field with Lx = Lσ(x) − Lσ(−x), so the best SA-NODE error there is zero while (2.1) is large. The example shows the estimate is tight step by step, not that the approximation error grows. The two do not merely multiply in (2.1). The Grönwall factor sits in the exponent, and the tube growth enters the same exponent through the Sobolev norm on ΩL,T , which is what makes that bound double exponential. The width that (2.1) certifies to be sufficient for accuracy ε is 2 ε−2 . That prescribed width therefore deteriorates at the same rate, the Grönwall of order CT,K,f factor alone contributing e2LT . Both statements concern the bound and the width it prescribes, not the smallest width that suffices, which is unknown. Whether every SA-NODE must incur the factor is a separate question, and no matching lower bound over Lipschitz targets is known. The next lemma says what such a bound would have to control. It converts closeness of flows into closeness of fields, at the cost of a square root, and so reduces a lower bound for the flow problem to one for approximation of the target field on the spacetime reachable set by the class (1.1). Lemma 2.2 (flow closeness forces field closeness). Write ΨΘ (t; s, x) and Ψ(t; s, x) for the solutions b Θ (t; x) = ΨΘ (t; 0, x) and at time t of (1.1) and of the target, launched from x at time s, so that Φ 2 Φ(t; x) = Ψ(t; 0, x). Let T ≥ 1 and r > 0, and let f and fΘ be C with ∥f ∥C 2 , ∥fΘ ∥C 2 ≤ B on Nr × [0, T ], where Nr is the closed r-neighborhood of KT and the norms are taken jointly in (x, t); set C2 := 2B(B + 1). If ε ≤ C2 e−B /4, if (1 + eB )ε ≤ r/2, and if sup

b Θ (t; x0 ) − Φ(t; x0 ) Φ

≤ ε,

x0 ∈K, t∈[0,T ]

then for every t0 ∈ [0, T − 1] and every x ∈ Φ(t0 ; K), p ∥fΘ (x, t0 ) − f (t0 , x)∥ ≤ 2 ε eB C2 .

(2.4)

The proof is in Appendix A.1. By (2.4), a width lower bound for the flow problem follows from one for uniform approximation of f on the spacetime reachable set by the fields (1.1), which are shallow networks in the joint variable (x, t). Problem 6.1 states the question in that form. The linear field ẋ = Lx cannot witness it: with σ the ReLU, Lx = Lσ(x) − Lσ(−x) componentwise, so it is represented exactly at width 2, one unit for each sign and its approximation error is zero at every horizon. A witness must therefore be a target whose field is hard to approximate on a domain that its own flow makes large. Any algorithm relying on (2.1) for a single approximant trained once on [0, T ] must therefore allow its width to grow exponentially. The remedy is to change the learning problem: partition

8

Z. LI AND N. M. MATZAKOS

the horizon, or exploit the orbit structure of the target. These are the two routes of Sections 3 and 4. 3. The model predictive strategy The Grönwall mechanism behind (2.1) is a consequence of global error accumulation, not of the SA-NODE class. This section replaces the single model on [0, T ] by a composite of short-horizon models, each trained against the true trajectory at its left endpoint. The data reset severs the error chain at every switch time, and the Grönwall factor is confined to a single window. The main results are the uniform data-assisted bound of Theorem 3.5 and the predicted-IC bound of Proposition 3.7, which shows what is lost when the resets are withheld. 3.1. The composite flow and the data reset. This subsection defines the composite model and its two evaluation modes, restarted from data or from its own predictions. The distinction between the two carries the whole analysis. Let 0 = τ0 < τ1 < · · · < τN = T be a partition of [0, T ], and write τmax for a bound on its mesh, maxk (τk+1 − τk ) ≤ τmax . The analysis treats τmax as a fixed quantity, independent of the horizon T ; that the partitions produced in Section 5 respect a horizon-independent mesh is reported there. On each window Ik := [τk , τk+1 ), closed at the right for k = N − 1, we place one SA-NODE fΘk driven by the local time s := t − τk ,
ẋ(t) = fΘk x(t), t − τk , t ∈ Ik . (3.1) Local time keeps the time input s of the bias A2i s + Bi in [0, τmax ] rather than [0, T ], and it makes the warm start Θk ← Θk−1 shape-preserving. The composite flow admits two evaluation modes, and the analysis rests on this distinction. In data-IC mode the true state is injected at each switch time,
b Θ t − τk ; Φ(τk ; x0 ) , b x0 ) = Φ Φ(t; t ∈ Ik . (3.2) k

Window k then sees no upstream error: its accuracy depends on its own training only. In predictedIC mode the composite prediction is propagated,
b pred (t; x0 ) = Φ b Θ t − τk ; Φ b pred (τk ; x0 ) , Φ t ∈ Ik , (3.3) k

which is the classical single-shooting chain. Mode (3.2) is the training-time and validation-time regime; mode (3.3) is deployment on a novel initial condition. Theorem 3.5 concerns the first, Proposition 3.7 the second. Figure 3.1 shows the three ingredients: the colored bands carry the per-window fields fΘk , the curved arrows copy parameters between adjacent blocks (the warm start), and the vertical arrows inject the true state at each switch time (the data reset). 3.2. The adaptive algorithm. This subsection specifies how the partition is chosen at run time, through the stopping rule (3.4), and states the training procedure as Algorithm 3.1. The partition is not fixed in advance. After window k is trained, the algorithm simulates it on the remaining horizon and declares τk+1 to be the first time the supervised error crosses the tolerance, n o b Θ (t − τk ; Φ(τk ; x0 )) > ε , τk+1 = inf t ∈ (τk , T ] : sup Φ(t; x0 ) − Φ (3.4) k x0 ∈K

with inf ∅ = T . A minimum window length smin ∆t, with ∆t the sampling step of the training data and smin an integer, excludes degenerate windows, and a window cap S bounds the total count. Training window k on the full remaining interval would waste effort on times covered by a later window, so the loss is evaluated only on [τk , min(τk + H, T )], where H > 0 is a training horizon playing the role of the MPC prediction horizon. The stopping rule (3.4) is still evaluated on the full remaining interval. Truncating the training interval at H changes the constant of a window but not the logic of the reset, so the analysis below is unaffected by the choice of H. The stopping rule and the reset together produce a sawtooth error profile, measured in Figure 5.2(a): a monolithic model accumulates error over the whole horizon, while the composite

LONG-TIME APPROXIMATION VIA SA-NODES copy

Θ0

I0

0

copy

Θ1

τ1

9

I1

Θ2

τ2

I2

warm start Θk−1 → Θk

τ3 = T

t

data reset Φ(τk ; x0 )

Φ(τ0 ; x0 )

Φ(τ1 ; x0 )

Φ(τ2 ; x0 )

Figure 3.1. The MPC–SA-NODE composite: one SA-NODE Θk per band Ik , warm-start copies (curved arrows), and data resets at the switch times (dashed arrows, mode (3.2)); at deployment the dashed arrows are replaced by (3.3). Algorithm 3.1 Adaptive MPC–SA-NODE training Require: sampled trajectory data (t, x) on [0, T ], tolerance ε, width P , epochs Nep , min window smin , window cap S, training horizon H, warm-start flag 1: k ← 0; tleft ← 0; Θprev ← ∅ 2: while tleft < T and k < S do 3: open window Wk on [tleft , min(tleft + H, T )]; copy weights from Θprev if warm start 4: train Wk for Nep epochs on the supervised loss with a regularizer as in [16, §4] (we use the extended Barron norm), with the true state at tleft as initial condition 5: evaluate Wk on [tleft , T ]; find τk+1 by (3.4), enforce τk+1 − tleft ≥ smin ∆t 6: shrink Wk to [tleft , τk+1 ]; tleft ← τk+1 ; Θprev ← Θk ; k ← k + 1 7: end while b of k windows covering [0, tleft ]; tleft = T unless the window cap binds Ensure: composite Φ (Section 3.3)

restarts from the data at each τk and its error re-enters the tolerance band. The experiments of Sections 5.2 and 5.3 measure this contrast. 3.3. Error analysis. The analysis rests on three structural assumptions on the target and one admissibility condition on the partition the algorithm produces. They lead to Theorem 3.5, the uniform data-assisted bound. Assumption 3.1 (Lipschitz target). f is L-Lipschitz in x, uniformly in t ≥ 0, on a neighborhood S of the forward reachable tube t≥0 Φ(t; K). Since Theorem 2.1 asks for a globally Lipschitz field, we tacitly replace f outside that neighborhood by a globally Lipschitz extension in the same Sobolev class; the flow from K is unchanged, and the constants below depend on the extension only through its bounds. s (Rd × [0, ∞)) with s > (d + 1)/2 + 2, so that Assumption 3.2 (approximation class). f ∈ Hloc Theorem 2.1 applies on every window; throughout this section the activation is ReLU, as that theorem requires.

10

Z. LI AND N. M. MATZAKOS

Assumption 3.3 (bounded reachable tube). The forward reachable tube K∞ := compact, and the target has bounded data on it in the sense that

S

t≥0 Φ(t; K) is

ϱ MH s := sup ∥f (·, · + τ )∥H s (K∞ ×[0,τmax ]) < ∞,

τ ≥0 ϱ where K∞ is the closed ϱ-neighborhood of K∞ , for a fixed ϱ > 0 whose choice the constants below

absorb. Here s is additionally taken at least d/2 + 3, the regularity at which the explicit constant (2.1) is available; the τ -uniformity of the window constant below is derived from that explicit form. Assumption 3.3 holds for dissipative or positively invariant targets whose data are uniformly bounded in time; for autonomous and time-periodic fields the time-uniformity of the data is automatic, so there the assumption reduces to boundedness of the reachable tube. Boundedness of the trajectories, uniform over x0 ∈ K and t ≥ 0, alone makes K∞ compact, but does not by itself bound MH s for a genuinely time-dependent target. The assumption fails for the barrier example ẋ = Lx, whose reachable tube grows without bound. Write Cτmax ,K∞ ,f := sup Cτmax , Φ(τ ;K), f (·,·+τ )

(3.5)

τ ≥0

for the uniform window constant, the supremum of the constants of Theorem 2.1 over all windows of length τmax started along the flow. The assumption makes this constant finite and independent of T , for two reasons. Every window trajectory is a true trajectory of the target, so each window’s ϱ tube lies in the fixed compact K∞ and its Sobolev data are measured on K∞ × [0, τmax ]; and MH s bounds those data for every time shift at once. The form inherited from (2.1) is exponential in ecLτmax MH s , hence double exponential in τmax , because the Sobolev norm of the field enters the Grönwall exponent through the Lipschitz constant of the learned field. What matters here is that the bound involves only τmax , MH s and L, never the horizon T . Assumption 3.4 (ε-admissible partition). Algorithm 3.1 terminates with τN = T , its partition has mesh maxk (τk+1 − τk ) ≤ τmax , and every window satisfies εwin := k

sup t∈Ik , x0 ∈K

b Θ (t − τk ; Φ(τk ; x0 )) Φ(t; x0 ) − Φ k

≤ ε.

(3.6)

Both conditions of Assumption 3.4 are properties of the partition that the algorithm produces, and both are verified on the output rather than imposed during the run; Section 5 reports the realized mesh and the realized per-window error. The stopping rule (3.4) enforces the error condition when it alone determines the partition, up to the boundary case of a first crossing exactly at t = T , which perturbs only the endpoint value of the last window. Two effects perturb it. The minimum window length smin ∆t may push a window past its crossing; when it binds, Assumption 3.4 holds with ε replaced by maxk εwin k , and every statement below holds verbatim with that value. The supremum in (3.4) is evaluated on a finite grid of initial conditions and times. We do not bound the gap to the continuous supremum: that would require a fill distance for the grid together with a Lipschitz bound on the error, which we do not establish. Assumption 3.4 is therefore verified on the grid, and every experimental statement in Section 5 is to be read grid-wise. The window cap S behaves differently. If the algorithm reaches k = S with τS < T , it terminates before covering the horizon. The final stretch [τS , T ] then carries no trained window, so Assumption 3.4 fails outright rather than by a bounded factor, and Theorem 3.5 does not apply. A cap that binds indicates that the per-window width is too small for the tolerance; it is not a bounded excess to be absorbed. Section 5.3 exhibits this regime; the safeguard excess, when the cap does not bind, is quantified in Sections 5.2–5.3. Theorem 3.5 (data-assisted composite approximation). Under Assumptions 3.1–3.4, the data-IC composite flow (3.2) satisfies sup t∈[0,T ], x0 ∈K

b x0 ) Φ(t; x0 ) − Φ(t;

≤ ε,

(3.7)

LONG-TIME APPROXIMATION VIA SA-NODES

11

and a total width sufficient to achieve εwin ≤ ε on all N windows obeys k Ptotal =

N −1 X

Pk ≲ N Cτ2max ,K∞ ,f ε−2 ,

(3.8)

k=0

with Cτmax ,K∞ ,f independent of T . The proof is in Appendix A.2. Remark 3.6 (how the budget depends on T ). The budget (3.8) is proportional to the window count, and the window count is bounded below: window lengths are at most τmax and sum to T , so every admissible partition has N ≥ ⌈T /τmax ⌉. (3.9) The uniform partition of mesh τmax attains (3.9), and for it the budget is exactly linear in T , replacing the Grönwall factor e2LT of the monolithic width by T /τmax . The data reset thus confines the Grönwall factor to a single window; the growth of the reachable tube remains. When Assumption 3.3 fails, the uniform bound on the per-window constant is lost, because the constant may grow with T through its domain Φ(τk ; K). For the spreading barrier example ẋ = Lx the certified constant inherited from (2.1) does grow, so the certified budget is super-linear in T , though the Grönwall factor stays capped at eLτmax ; the smallest sufficient width does not grow there, since the linear field is exactly representable (Section 2.2), so the statement concerns the certificate and not the class. (Unboundedness of the tube alone does not force even this: a pure drift ẋ = c has an unbounded tube but translation-invariant window data, and its certified budget stays linear.) The adaptive partition of Algorithm 3.1 is not uniform, and Theorem 3.5 imposes no upper bound on its N . The algorithm itself does impose one, N ≤ T /(smin ∆t) + 1 from the minimum window length, but its constant is the data resolution, not the target. Whether a target-controlled linear law survives without the safeguard is Problem 6.2; Section 5.3 measures the realized counts. 3.4. Deployment on novel initial conditions. Without resets the composite reduces to a single-shooting chain, and upstream errors are amplified by each learned flow. The resulting bound is adapted to the partition. Proposition 3.7 (predicted-IC bound). Let n o L̄ := max L, max Lipx fΘk , k

the Lipschitz constants being taken on a compact neighborhood of the reachable tube large enough to contain the predicted trajectories and the data-restarted window trajectories. This constant is finite, since each fΘk is globally Lipschitz, with an explicit constant in terms of its weights. It is an a posteriori constant of the trained model and the evaluation set, not one controlled a priori by the assumptions. Under Assumptions 3.1–3.4, for every x0 ∈ K the predicted-IC flow (3.3) satisfies sup

b pred (t; x0 ) Φ(t; x0 ) − Φ

≤ ε

t∈[0,T ]

N X k=1

eL̄ (T −τk ) ≤ ε

eN L̄τmax − 1 , eL̄τmax − 1

(3.10)

the quotient being read as its limit N when L̄ = 0. The proof is in Appendix A.2. Remark 3.8 (reading the bound). The bound (3.10) has three regimes. For the uniform partition, N τmax = T and the quotient equals eL̄T − 1 ε L̄τ ; e max − 1

12

Z. LI AND N. M. MATZAKOS

for non-uniform partitions N τmax ≥ T , that substitution is not available, and only the N -dependent quotient of (3.10) is a valid upper bound. As L̄ → 0 both bounds tend to N ε, so the perwindow errors simply add. For L̄τmax ≫ 1 the bound behaves as ε e(N −1)L̄τmax . Deployment therefore remains exponential in the horizon; the gain over the monolithic model is the prefactor (eL̄τmax − 1)−1 , not the exponent. A continuity penalty at the switch times should shrink the constant of the chain; we leave its analysis open. 3.5. Scope. Theorem 3.5 is a statement about the data-assisted regime: it requires the true state at the N switch times, which is available during training and validation but not in autonomous deployment. The scale of τmax is expected to be set by the local Lipschitz behavior of f along the trajectory, and the experiments are consistent with this; a rigorous link to Lyapunov exponents is open (Problem 6.2). Beyond a Lipschitz target with a bounded reachable tube (Assumptions 3.1– 3.3) the strategy makes no structural assumption: within that class it applies to non-autonomous and non-periodic systems alike, and the data dependence provides this generality. The tube hypothesis confines the parameter budget to grow only with the window count; it is the point where the target’s global behavior re-enters the estimate. 4. The Floquet strategy The MPC strategy needs the true state at every switch time. When the target is autonomous and its trajectories are attracted to a stable limit cycle, the dynamics itself can play the role of the reset: transverse perturbations are contracted at every return to a Poincaré section. This section shows how a trained SA-NODE can be made to inherit that contraction, with a certificate. Standard neural ODEs are obstructed here: activation saturation caps the contraction their flows can express [20]. The periodic encoding below is designed against this obstruction. The Floquet loss supplies a differentiable surrogate that drives the learned return map towards contraction, and a partial repair when the encoding is absent (Section 5.4). What that surrogate does and does not certify is the subject of Section 4.3. The main results are Theorem 4.13, for an autonomous learned field, and Theorem 4.22, for the periodic encoding we deploy. In the first, certified contraction of the learned return map confines the deployment error to linear growth in the number of periods; in the second, the measured monodromy spectrum yields a uniform-in-time orbital bound. 4.1. Setting and assumptions. This subsection fixes the geometry of the problem and states the standing assumptions. The essential one is the certified contraction of Assumption 4.6; the others make it meaningful. The target is an autonomous system ż = f (z),

f ∈ C 2 (Rd ; Rd ),

d ≥ 2,

(4.1)

with flow Φf and a T̂ -periodic orbit γ̂ with image Γ := γ̂([0, T̂ ]). The following notation is used throughout this section and in Appendices A.3–A.7: p := γ̂(0), Nδ0 := {x : dist(x, Γ) ≤ δ0 },

v := ∥f (p)∥ ,

vmin := min ∥f (γ̂(s))∥ > 0,

L := sup ∥Df ∥ ,

Σr := {p + w : w ⊥ f (p), ∥w∥ ≤ r}.

s

Nδ0

As announced in Section 1.2, the parameter vector Θ of (1.1) is abbreviated Θ. Here δ0 > 0 is fixed small enough that the compact tube Nδ0 lies in the basin of attraction of Γ and that the normal tubular chart (s, u) 7→ γ̂(s) + U (s)u is injective on it. Such a δ0 exists because Γ is a compact embedded C 2 curve, and it depends only on f and γ̂. The speed vmin is the lower bound used for the winding time in Lemma 4.10, and v the speed at the base point, used for the return-time estimates. The disc Σr is the transverse disc of radius r at the base point p, inside the hyperplane Σ through p orthogonal to f (p), and distH is the Hausdorff distance between compact sets. Five spectral quantities appear in this section and are pairwise distinct: the transverse Floquet radius ρT (f ) of the target orbit; the prescribed contraction threshold ρ∗ ; the spectral radius ρΘ of the

LONG-TIME APPROXIMATION VIA SA-NODES

13

learned return map at its fixed point; the divergence surrogate ρ̃T of (4.8); and the spectral radius ρ(M ) of the one-period stroboscopic monodromy of the encoded flow (Remark 4.19). They coincide only where explicitly stated. We take the constant B of Assumption 4.4 to bound both fields,  max ∥f ∥C 2 (Nδ ) , ∥fΘ ∥C 2 (Nδ ×R) ≤ B, 0

0

which costs nothing, since the constants below are already allowed to depend on f , and which lets every statement about a general field be applied to f itself. All constants K1 , K2 , . . . and C below depend only on (f, γ̂, δ0 , B, ρ∗ , ρT (f )), with ρ∗ the contraction threshold of Assumption 4.6, and never on ε, t, or the initial condition. Assumption 4.1 (hyperbolic stable orbit). The orbit γ̂ is hyperbolic and orbitally asymptotically stable; its transverse Floquet spectral radius satisfies ρT (f ) < 1. Throughout, T̂ denotes the minimal period of γ̂; the statements that let the target’s time-T̂ map fix every point of Γ use this convention. Here and below the subscript in ρT (·) and in the surrogate ρ̃T marks the transverse direction and is unrelated to the horizon T of Section 3. Floquet theory originates in [10]; asymptotic stability of a periodic orbit whose nontrivial multipliers lie inside the unit circle is [8, Theorem 2.82]. The learned model must have a well-defined return map, so its field should be T̂ -periodic in time. We enforce this in the architecture: the scalar time input of (1.1) is replaced by the encoding
t 7−→ sin(2πt/T̂ ), cos(2πt/T̂ ) , (4.2) so that fΘ (·, t + T̂ ) = fΘ (·, t) exactly, and the model is deployed by continuous integration with no time resets. With smooth activations, fΘ ∈ C ∞ . Explicitly, the deployed field is fΘ (x, t) =

P X

 
⊤ Wi ◦ σ A1i x + Ci sin(2πt/T̂ ), cos(2πt/T̂ ) + Bi ,

Ci ∈ Rd×2 ,

(4.3)

i=1

with the remaining parameters as in (1.1): the scalar-time column A2i is replaced by the two-column block Ci , so the encoded model carries P d(d + 4) trainable parameters instead of P d(d + 3). It is a different architecture class from the one covered by Theorem 2.1, which is stated for the scalar-time ReLU model and is not invoked for it; one-period accuracy of (4.3) is the measured hypothesis of Assumption 4.4. Training minimizes the one-period trajectory loss Z T̂ Ltraj (Θ) := ∥ΦΘ (t; p) − γ̂(t)∥2 dt, 0

whose empirical form in Algorithm 4.1 is the mean squared error on the training grid. Assumption 4.2 (near-autonomy on the tube). The learned field is T̂ -periodic in time, fΘ (·, t + T̂ ) = fΘ (·, t), as the encoding (4.2) enforces. Write Z 1 T̂ fΘ (x, s) ds f¯Θ (x) := T̂ 0 for the time-averaged field, and   η := sup fΘ (x, t) − f¯Θ (x) + Dx fΘ (x, t) − Dx f¯Θ (x) (4.4) x∈Nδ0 , t∈R

for its C 1 oscillation on the tube. The case η = 0 is the autonomous case, in which we write fΘ (x). Remark 4.3 (role of the oscillation). The analysis of Section 4.2 is carried out in the autonomous case η = 0, where the certificates are exact. Lemma 4.15 removes the restriction, at the cost of replacing ε by ε + Cη in conclusions (ii)–(iii) of Theorem 4.13.

14

Z. LI AND N. M. MATZAKOS

Assumption 4.4 (flow closeness over one period). fΘ ∈ C 2 (Nδ0 ×R; Rd ) with ∥fΘ ∥C 2 (Nδ ×R) ≤ B. 0

We extend fΘ to Rd ×R by a C 2 cutoff outside Nδ0 , so that ΦΘ is globally defined and the supremum below is taken over trajectories that exist for all time. The extension multiplies the C 2 bound by a factor depending only on δ0 , which we absorb into B, and it is inactive on Nδ0 , where every trajectory the analysis uses is shown to remain. With r1 = r1 (f, γ̂, δ0 ) the confinement radius of Lemma 4.7, which is built from the true flow alone and does not depend on Θ, the one-period flow closeness is  ε := sup ∥ΦΘ (t; x) − Φf (t; x)∥ : x ∈ N2r1 , t ∈ [0, T̂ + 2] , and we write δ := supx∈Γ ∥fΘ (x) − f (x)∥ for the orbit mismatch. For η > 0 the supremum in δ is taken over x ∈ Γ and t ∈ R, and ΦΘ denotes the flow launched at phase 0. Assumption 4.4 is a hypothesis on the trained model, not a consequence of Section 2: Theorem 2.1 is stated for the ReLU activation, which is not C 2 , so no result of this paper asserts that the C 2 architectures deployed in this section achieve a prescribed ε. The smooth-activation variant of the approximation theorem recalled after Theorem 2.1 grounds achievability for twice weakly differentiable activations, and the experiments measure ε directly; the analysis below takes ε as given. The two small quantities of Assumption 4.4 are not independent. Flow closeness over one period controls the orbit mismatch, at the cost of a square root. This is what lets every smallness condition below be imposed on ε alone. Lemma 4.5 (the orbit mismatch is controlled by the flow closeness). 4.4 hold
Let Assumption √ and set C2 := B 2 + B + L supNδ ∥f ∥. If ε ≤ min δ0 /4, C2 (T̂ + 2)2 /2 , then δ ≤ 2C2 ε + 2η, and 0 √ δ ≤ 2C2 ε in the autonomous case η = 0. The proof is in Appendix A.3. Assumption 4.6 (certified transverse contraction). With PΘ : Σr2 → Σ the first-return map of fΘ furnished by Lemma 4.8 for ε ≤ ε2 , which is proved without the present assumption,
ρ DPΘ (p) ≤ ρ∗ < 1. (4.5) We write ρΘ for the spectral radius of the learned return map at its fixed point. By Lemma 4.9 that fixed point lies within O(ε) of p, and DPΘ is Lipschitz, so ρΘ ≤ ρ∗ + ω(Cε), with ω the modulus of continuity of the spectral radius on the norm ball {A : ∥A∥ ≤ CM } that contains all return-map differentials arising here (Lemma 4.8). That modulus is linear when DPΘ (p) is diagonalizable and only of order ε1/m at a Jordan block of size m. The two conditions are therefore equivalent up to that margin rather than interchangeable. For η > 0 the condition is posed for the return map of the averaged field f¯Θ (Lemma 4.15). Assumption 4.6 is the key hypothesis, and it is verifiable: Section 4.3 gives two sufficient conditions computable from a trained model. Theorem 4.13 should be read as a conditional stability result; we do not claim that training always achieves (4.5), and the experiments quantify how often it does. Assumption 4.4 deserves the same warning, and it is the more demanding of the two. It asks for one-period flow closeness at every initial state of the tube N2r1 , whereas the trajectory loss Ltraj is evaluated on the single base orbit through p. A small training loss therefore does not establish Assumption 4.4, and neither does a small Floquet loss, which concerns stability rather than accuracy. Uniform closeness on the tube is a separate approximation hypothesis, of the kind Theorem 2.1 supplies for a compact set of initial conditions, and the two certification routes of Section 4.3 presuppose it rather than prove it.

LONG-TIME APPROXIMATION VIA SA-NODES

15

4.2. Main theorem. The analysis tracks three quantities: the transverse distance to Γ, contracted at rate ρ∗ at each return; the distance between the two return maps, of order ε; and the cumulative phase drift between the two clocks, which the reset cannot remove. Six lemmas make this precise; their proofs are in Appendix A.3. Lemma 4.7 (one-period confinement). There exists r1 ∈ (0, δ0 /4], depending only on (f, γ̂, δ0 ), such that for every ε ≤ δ0 /4, all x ∈ N2r1 and t ∈ [0, T̂ + 2]: Φf (t; x), ΦΘ (t; x) ∈ Nδ0 and ∥ΦΘ (t; x) − Φf (t; x)∥ ≤ ε. Lemma 4.8 (return maps and return times). There exist r2 ∈ (0, r1 ], ε2 > 0 and K1 := 1 + 2B/v such that for ε ≤ ε2 the first-return maps Pf , PΘ : Σr2 → Σ and return times τf , τΘ are well defined and C 1 , with DPf , DPΘ Lipschitz, Pf (p) = p, τf (p) = T̂ , and uniformly on Σr2 |τΘ − τf | ≤ v2 ε,

∥PΘ − Pf ∥ ≤ K1 ε,

τf , τΘ ∈ [T̂ /2, T̂ + 1].

Lemma 4.9 (persistence of the orbit). Let Assumptions 4.1, 4.4 and 4.6 hold with η = 0. Set ρ̄ := (1+ρ∗ )/2 and ρ̄f := (1+ρT (f ))/2. There are adapted norms on Tp Σ with equivalence constants c1 , c2 and c′1 , c′2 , and radii r3 , r3′ ∈ (0, r2 ], such that ∥DPΘ ∥∗ ≤ ρ̄ on Σr3 and ∥DPf ∥∗f ≤ ρ̄f on Σr3′ . There is ε3 ∈ (0, ε2 ] such that for ε ≤ ε3 the map PΘ has a unique fixed point x∗Θ ∈ Σr3 with ∥x∗Θ − p∥ ≤ K2 ε, K2 := c2 K1 /(c1 (1 − ρ̄)); the orbit γΘ := {ΦΘ (t; x∗Θ ) : 0 ≤ t ≤ TΘ } with period TΘ := τΘ (x∗Θ ) is periodic, orbitally exponentially stable, and distH (γΘ , Γ) ≤ K3 ε,

TΘ − T̂ ≤ K4 ε.

If Ltraj (Θ) = 0 then x∗Θ = p, γΘ = Γ, TΘ = T̂ . Lemma 4.10 (monotone winding and uniqueness). There exist r5′ ∈ (0, δ0 /4], δ1 > 0 and CU ≥ 1, depending only on (f, γ̂, δ0 , B), with the following properties. (i) Winding and confinement. Let g ∈ C 2 (Nδ0 × R; Rd ) satisfy ∥g∥C 2 ≤ B and supΓ×R ∥g − f ∥ ≤ δ ≤ δ1 . Then for every r0 ≤ r5′ and every g-trajectory x(·) starting in Nr0 , one has dist(x(t), Γ) ≤ CU (r0 + δ) for all t ∈ [0, T̂ + 1], so the trajectory remains in the tube Nδ0 on that interval, and within that same time it crosses Σ transversally, at a point at distance at most CU (r0 + δ) from p. Nothing is claimed beyond the first winding, and nothing can be: the hypotheses bound g − f on Γ only, so the transverse dynamics of g may expand, and re-entering the estimate at the crossing replaces r0 by CU (r0 + δ) ≥ r0 . Confinement for all t ≥ 0 is recovered only through the contraction of Assumption 4.6, one winding at a time. This part uses no certificate: it holds for any field within δ1 of f on Γ, in particular for f itself, with
δ = 0. (ii) Uniqueness. Set r5 := min r5′ , r3 /(2CU ) , which depends on ρ∗ through r3 but on nothing else beyond (f, γ̂, δ0 , B), and suppose in addition that Assumption 4.6 holds and that δ ≤ r3 /(2CU ), which by Lemma 4.5 follows from a smallness condition on ε alone. Then CU (r0 + δ) ≤ r3 for every r0 ≤ r5 , and every periodic orbit of fΘ contained in Nr0 equals γΘ ; if moreover K3 ε ≤ r0 , then γΘ ⊂ Nr0 by Lemma 4.9 and γΘ is that unique orbit. In particular r5 does not depend on Θ or on ε. Lemma 4.11 (two contractions track each other). Set K5 := (c′2 /c′1 ) · 2K1 /(1 − ρ̄f ). There exist r4 ∈ (0, min(r3 , r3′ , r5 )] and ε4 ∈ (0, ε3 ] such that for ε ≤ ε4 and all yf , yΘ ∈ Σr4 with ∥yf − yΘ ∥ ≤ K1 ε, the iterates remain in Σr3 ∩ Σr3′ and Pfk (yf ) − PΘk (yΘ )

≤ K5 ε

for all k ≥ 0.

Lemma 4.12 (clock comparison). In the setting of Lemma 4.11, define the cumulative return times tf0 := 0, tfn+1 := tfn + τf (Pfn yf ) and analogously tΘ n . Then for all n ≥ 0, f tΘ n − tn ≤ K6 n ε,

K6 := v2 + Lip(τf ) K5 .

16

Z. LI AND N. M. MATZAKOS

Γ

learned flow ΦΘ

Φ Θ (kT̂ ; x0 )

drift ∼ k |TΘ − T̂ | Φ f (kT̂ ; x0 )

ρΘ < 1 Σ PΘk y → xΘ∗

(a) transverse contraction

(b) phase drift

Figure 4.1. The Poincaré reset mechanism. (a) Certified contraction ρΘ < 1: successive returns of the learned flow to the section Σ converge geometrically to the fixed point, so transverse error cannot accumulate. (b) The non-resettable error: the learned clock runs at period TΘ ̸= T̂ , so the phase drifts by k TΘ − T̂ after k periods, the linear term in (4.6). Theorem 4.13 (linear long-horizon bound for autonomous targets with a stable periodic orbit). Let Assumptions 4.1–4.6 hold with η = 0 in Assumption 4.2. There exist ε0 , r0 > 0 and C ≥ 1, depending only on (f, γ̂, δ0 , B, ρ∗ , ρT (f )) and never on ε, t, or x0 , such that for all ε ≤ ε0 : (i) Orbit persistence. fΘ has a unique periodic orbit γΘ ⊂ Nr0 , orbitally exponentially stable, with distH (γΘ , Γ) ≤ Cε and TΘ − T̂ ≤ Cε. (ii) Transverse contraction. For every x0 with dist(x0 , Γ) ≤ r0 and every t ≥ 0,  
dist ΦΘ (t; x0 ), Γ ≤ C ρ̄ t/(T̂ +1) dist(x0 , Γ) + ε . (iii) Linear bound. For every such x0 and t ≥ 0,  t . (4.6) ∥ΦΘ (t; x0 ) − Φf (t; x0 )∥ ≤ C ε 1 + T̂ The proof is in Appendix A.4; we summarize the mechanism. Grönwall is applied only over single periods, contributing the bounded factor eL(T̂ +1) . Transverse errors do not accumulate: Lemma 4.11 caps them at K5 ε uniformly in the number of returns. The accumulating quantity is the phase: each period the two clocks drift apart by O(ε) (Lemma 4.12), and after n ≤ 2t/T̂ + 1 returns the drift is O(nε), whence (4.6). Figure 4.1 displays the two mechanisms. Remark 4.14 (when the linear term is active). The linear term has two sources, which saturate differently. Near the threshold ρΘ → 1− the transverse geometric sums lose their gap and the transverse error itself grows linearly, so the linear term is expected to be attained there, and a pure period mismatch attains it. In polar coordinates on R2 , let the target be ṙ = −(r − 1), θ̇ = 2π/T̂ , whose unit circle Γ is a hyperbolic limit cycle with transverse multiplier e−T̂ . Let the learned field share the radial part, with θ̇ = 2π/((1 + δ)T̂ ) for small δ > 0. Both fields are autonomous and C 2 near Γ, and both return maps on the section {θ = 0} contract at rate at most e−T̂ , so Assumption 4.6 holds. By Grönwall, the one-period flow closeness on a tube around Γ is

LONG-TIME APPROXIMATION VIA SA-NODES

17

ε ≤ c1 δ, with c1 depending only on T̂ . For x0 ∈ Γ both solutions remain on Γ and differ only in phase. After k periods the angular mismatch is 2πkδ/(1 + δ), so ΦΘ (k T̂ ; x0 ) − Φf (k T̂ ; x0 ) =
2 sin πkδ/(1 + δ) ≥ 4kδ/(1 + δ) as long as kδ/(1 + δ) ≤ 21 , by sin(πs) ≥ 2s on [0, 12 ]. Since ε ≤ c1 δ, the stroboscopic error is at least (2/c1 ) ε k until the phase wraps: the envelope Cε(1 + k) of (4.6) is attained up to constants. Well inside the certified regime the transverse error saturates at O(ε/(1 − ρΘ )) and the phase drift dominates. The error at period k then grows like k TΘ − T̂ until the phase wraps around the orbit, after which it oscillates at the level of the orbit diameter, still below the envelope (4.6). When the fitted period matches T̂ closely, the drift per period is negligible and the error saturates far below the envelope; the deployments of Section 5.4 lie outside the scope of this remark, being governed by Theorem 4.22 rather than by (4.6). The encoded architecture is T̂ -periodic, not autonomous, so Theorem 4.13 does not apply to it verbatim. The following lemma closes the gap. Its proof, in Appendix A.6, has two steps. First, the averaged field f¯Θ is autonomous and satisfies Assumption 4.4 with ε + Cη in place of ε, the constant being made explicit in the lemma. Once its return map is certified, which is the hypothesis of part (ii) of Lemma 4.15 below, Theorem 4.13 applies to it. Second, the proof of part (iii) of Theorem 4.13 reruns with f¯Θ as reference and the encoded field as the perturbed system, once the semigroup identities of the perturbed flow are replaced by evolution-family identities for the (s) phase-indexed flows ΦΘ , whose closeness part (i) controls uniformly in the launch phase. Lemma 4.15 (reduction of the periodic encoding). Let Assumptions 4.1, 4.2 and 4.4 hold, and set
CB := (T̂ + 2) 2 + B(T̂ + 2) e2B(T̂ +2) . (s)

Assume that ε + CB η ≤ δ0 /4 and that δ + η ≤ vmin /(2(T̂ + 1)). For s ∈ R write ΦΘ (t; x) for the solution of ẏ = fΘ (y, s + t), y(0) = x: the flow launched at phase s. (i) For every s ∈ R, x ∈ N2r1 and t ∈ [0, T̂ + 2], (s)

(s)

ΦΘ (t; x) − Φf¯Θ (t; x) + Dx ΦΘ (t; x) − Dx Φf¯Θ (t; x)

≤ CB η.

In particular f¯Θ satisfies Assumption 4.4 with closeness ε + CB η and orbit mismatch δ + η. (ii) Suppose in addition that the averaged field is certified, ρ(DP̄ (p)) ≤ ρ∗ < 1 with P̄ its firstreturn map (the proof uses this hypothesis only through the contraction at rate ρ̄ = (1+ρ∗ )/2 that it induces, via the adapted norm of Lemma 4.9, on a neighborhood of the fixed point) and that ε + η ≤ ε6 for a threshold ε6 > 0 depending only on (f, γ̂, δ0 , B, ρ∗ , ρT (f )). Then three conclusions hold. Conclusion (i) of Theorem 4.13 holds for f¯Θ , with ε+CB η in place of
(0) ε. Every trajectory of the encoded flow with dist(x0 , Γ) ≤ r0 satisfies dist ΦΘ (t; x0 ), γ̄Θ ≤
C ρ̄ t/(T̂ +1) dist(x0 , γ̄Θ ) + η for all t ≥ 0, where γ̄Θ is the orbit of f¯Θ ; in particular the distance is eventually of order η. Conclusions (ii)–(iii) hold for the encoded flow itself: for every x0 with dist(x0 , Γ) ≤ r0 and t ≥ 0,  t . (4.7) ∥ΦΘ (t; x0 ) − Φf (t; x0 )∥ ≤ C (ε + η) 1 + T̂ The lemma trades the autonomy of the field for a measurable quantity: η is a supremum of the trained network’s deviation from its time average, evaluated on a grid over the tube by forward passes. It should be read as an asymptotic bridge between the two theorems, valid only for small η. Section 5.4 shows that the trained models lie far outside that regime, so the deployed architecture is governed by Theorem 4.22 alone.

18

Z. LI AND N. M. MATZAKOS

4.3. Certification. Theorem 4.13 is conditional on (4.5), which can be verified along two routes: an exact scalar identity in the autonomous plane, and a perturbative C 1 route in any dimension. Both certify the return map of an autonomous field, which is what Assumption 4.6 asks for and which, when η > 0, is the return map of the average f¯Θ . For the time-periodic architecture we deploy we compute instead the full spectrum of the stroboscopic monodromy. That spectrum does not verify (4.5), and Section 4.4 shows what it does deliver. Proposition 4.16 (C 1 -closeness certification). Write g := fΘ in the autonomous case η = 0 and g := f¯Θ for η > 0, and let ε1 := ∥fΘ − f ∥C 1 (U ) on a tubular neighborhood U of Γ (for η > 0 the norm is taken on U × R/T̂ Z, and averaging gives ∥g − f ∥C 1 (U ) ≤ ε1 , so the conclusions concern the field that Assumption 4.6 requires). There exist C > 0 and ε∗ > 0, depending only on (f, Γ, B), such that for ε1 < ε∗ the following hold. The bound B enters through the uniqueness radius of part (i): a C 1 distance alone gives no modulus of continuity for DPΘ , and the C 2 bound of Assumption 4.4 supplies one. (i) g has a unique periodic orbit ΓΘ near Γ, and distH (ΓΘ , Γ) ≤ Cε1 . (ii) The transverse monodromy matrices satisfy MΘ⊥ − Mf⊥ ≤ Cε1 . (iii) ρΘ ≤ ρ̄f < 1, with ρ̄f = (1 + ρT (f ))/2 as in Lemma 4.9. In particular, after one further shrinking of ε∗ , Assumption 4.6 holds with ρ∗ := (1 + ρ̄f )/2. The step from (iii) to (4.5) is not immediate, since (iii) bounds the spectral radius at the fixed point x∗Θ of the learned return map while (4.5) is posed at the base point p. By (i) the two points lie within Cε1 of each other and DPΘ is Lipschitz on the relevant disc, so ρ(DPΘ (p)) ≤ ρ̄f + ω(C ′ ε1 ), with ω the modulus of continuity of the spectral radius on the norm ball {A : ∥A∥ ≤ CM } of the proof of Lemma 4.9; as ρ̄f < 1, this falls below (1 + ρ̄f )/2 once ε∗ is small enough. The behavior of ω is as described after (4.5), so the smallness demanded of ε∗ depends on f through it. The conclusions concern an autonomous field: when η = 0 they apply to fΘ itself, and when η > 0 to the average f¯Θ , which is what Assumption 4.6 requires. For η > 0 nothing is asserted here about a periodic orbit or a monodromy of the time-periodic field fΘ ; that object is treated in Section 4.4. The proof, a C 1 perturbation argument along the reference orbit, is in Appendix A.5. The second route is computational and integrates with training. In dimension two the Liouville–Abel identity gives the derivative of the return map at its fixed point exactly (Lemma A.2, Appendix A.3): Z T Θ  DPΘ (x∗Θ ) = exp div fΘ (γΘ (t)) dt . 0

Definition 4.17 (Floquet loss). The surrogate evaluates the divergence integral along the learned orbit over one nominal period,  Z T̂ ρ̃T (Θ) := exp divx fΘ (γ̂Θ (t), t) dt , (4.8) 0

where γ̂Θ is the curve obtained by relaxing the learned flow onto its attractor and sampling one nominal period. Two provisos attach to this definition. First, it presupposes that the relaxation returns a closed fΘ -orbit contained in Nr5′ , with r5′ the certificate-free radius of Lemma 4.10(i). Relaxation onto an attractor does not guarantee this by itself. The implementation tests closure, rejecting the model when the endpoint gap exceeds a quarter of the orbit diameter, and rejects trajectories that leave a fixed sanity ball. It does not test membership of the tube Nr5′ around Γ, so the locality half of this proviso is assumed rather than verified (Section 5.4). The radius is stated in the certificate-free form on purpose, since the surrogate is what the certificate will be derived from. Second, even for a genuine orbit the surrogate is not the exact multiplier. The Liouville–Abel identity integrates the divergence over the learned period TΘ , whereas (4.8) integrates over the nominal period T̂ . The two agree when Ltraj (Θ) = 0, in which case TΘ = T̂ ; in general they differ

LONG-TIME APPROXIMATION VIA SA-NODES

19

by O(ε), quantified in the proof of Corollary 4.18. The surrogate is differentiated in Θ through the analytic Jacobian of (1.1); the implicit dependence of γ̂Θ on Θ is not differentiated through. The Floquet loss is
2 LF (Θ) := µF max 0, ρ̃T (Θ) − ρ∗ , µF > 0, ρ∗ ∈ (0, 1). (4.9) Corollary 4.18 (Floquet-loss certification, d = 2). Let d = 2 and Assumptions 4.1–4.4 hold with η = 0 in Assumption 4.2. Assume that the relaxation defining γ̂Θ in (4.8) returns a closed fΘ -orbit contained in Nr , for a radius r ≤ r5′ depending only on (f, γ̂, δ0 , B, ρT (f )) and fixed in the proof, and let ε be small enough that δ ≤ δ1 , which Lemma 4.5 guarantees. If the trained model satisfies ρ̃T (Θ) ≤ ρ∗ , equivalently LF (Θ) = 0, then: (i) if Ltraj (Θ) = 0 and the relaxed orbit passes through the base point, so that γ̂Θ = Γ, then ρ(DPΘ (p)) = ρ̃T (Θ) ≤ ρ∗ exactly; (ii) in general the learned orbit has a return-map fixed point x∗Θ with ∥x∗Θ − p∥ = O(ε) and |ρ(DPΘ (x∗Θ )) − ρ̃T (Θ)| ≤ Cε, so Assumption 4.6 holds with contraction constant ρ∗ + Cε < 1 whenever ε ≤ (1 − ρ∗ )/(2C). In either case Theorem 4.13 applies, and the chain (1.6) holds. The proof, in Appendix A.5, constructs x∗Θ from Assumptions 4.1–4.4 alone, by a fixed-point argument for the perturbed return map, so the contraction is derived and not assumed. Under near-autonomy the surrogate transfers as well. The divergences of fΘ (·, t) and of f¯Θ differ by at most dη on the tube, so along the learned orbit the one-period divergence integrals differ by at most dT̂ η and the surrogates by the factor edT̂ η = 1 + O(η). The O(η) displacement of the averaged field’s orbit contributes at the same order, and is absorbed in the margin of Corollary 4.18(ii). The transfer is effective only while dT̂ η ≪ 1. For the models of Section 5.4, with the measured C 1 oscillations, this product is about 6 on Stuart–Landau and 35 on van der Pol, so there the surrogate must be replaced by the spectral computation of Remark 4.19. Remark 4.19 (what the surrogate does and does not control when η > 0). For a deployed encoded model the field fΘ (·, t) is genuinely T̂ -periodic. Let M denote the one-period stroboscopic monodromy of the encoded flow along the locked orbit, with Floquet multipliers λ1 , λ2 . By Liouville’s formula the divergence integral (4.8) returns the determinant, ρ̃T (Θ) = det M = λ1 λ2 ,

(4.10)

whereas the stability of the locked orbit is governed by the spectral radius ρ(M ) := maxi |λi |. A periodic non-autonomous system has no trivial unit multiplier, so ρ̃T ≤ ρ∗ bounds the product of the multipliers and not their largest modulus. Only in the autonomous limit η = 0 does the surrogate measure a single multiplier. The tangent multiplier is then pinned at one and excluded from the stability question, and det M equals the transverse multiplier itself (Lemma A.2); this is the regime that Lemma 4.15 certifies. When η is not small the surrogate can pass a model with λ1 > 1 > λ1 λ2 , so a small ρ̃T is not by itself a proof of stability. The reliable diagnostic is to compute both eigenvalues of M , one 2 × 2 variational solve along the orbit; Section 5.4 scores every trained model this way. Beyond the autonomous plane no scalar identity controls the spectral radius. This is the case for d > 2, where the transverse monodromy MΘ⊥ ∈ R(d−1)×(d−1) is a matrix, and already for d = 2 when the field is periodic in time (Remark 4.19). Certification then proceeds by the C 1 route. Let  δ∗ := sup δ > 0 : ρ(Mf⊥ + ∆) < 1 for every ∥∆∥ ≤ δ be the stability margin of the true transverse monodromy. For diagonalizable Mf⊥ = V ΛV −1 the Bauer–Fike theorem [1] gives δ∗ ≥ (1 − ρT (f ))/κ(V ), with κ(V ) := ∥V ∥ V −1 the eigenvector condition number. The bound of Proposition 4.16(ii) then yields ρΘ < 1 whenever Cε1 < δ∗ (Problem 6.3). Both routes are post hoc certificates for an autonomous field, and neither asserts that training achieves a small LF or a small ε1 . For the time-periodic architecture that we deploy the situation

20

Z. LI AND N. M. MATZAKOS

Algorithm 4.1 Floquet-regularized SA-NODE training, with post hoc monodromy diagnostics Require: period T̂ , one-period trajectory data from a base point on Γ, threshold ρ∗ ∈ (0, 1), weight µF > 0, epochs Nep , orbit-refresh interval R 1: build the SA-NODE with the periodic encoding (4.2); initialize Θ randomly 2: for epoch = 1, . . . , Nep do 3: Ltraj (Θ) ← MSE between the model trajectory over [0, T̂ ] and the data 4: every R epochs: relax the learned flow onto its attractor and resample the learned orbit γ̂Θ (detached)
R T̂ 5: ρ̃T (Θ) ← exp 0 divx fΘ (γ̂Θ (t), t) dt by quadrature, differentiable through divx fΘ 6: update Θ by Adam on Ltraj (Θ) + µF max(0, ρ̃T (Θ) − ρ∗ )2 7: end for Ensure: trained model; the surrogate ρ̃T (Θ), the quantity trained against; and the monodromy spectrum from one variational solve along the converged orbit. The model meets hypothesis (B1) of Theorem 4.22 if ρ(M ) = maxi |λi | ≤ ρ∗ , which ρ̃T (Θ) ≤ ρ∗ does not imply (Remark 4.19); this is not the certificate (4.5)

is weaker still: a small LF controls det M and not ρ(M ), so it does not by itself certify stability (Remark 4.19). For such models we compute the full monodromy spectrum instead. It is not the hypothesis (4.5), but it is exactly hypothesis (B1) of Theorem 4.22, which is where the deployed architecture is treated. Algorithm 4.1 therefore reports both quantities at convergence, and how often the certificate is met is an empirical question, answered in Section 5.4. The oscillation-penalized runs of Section 5.4 add λosc η̂ 2 to the loss of line 6, with η̂ 2 the mean over a fixed grid of tube points of the variance of fΘ (x, ·) across a uniform grid of phases in [0, T̂ ). It is a C 0 surrogate for (4.4), which is a C 1 supremum, so it penalizes the oscillation without bounding it. The grid is evaluated, matching Definition 4.17. While the early learned flow is still unstable, the true training orbit serves as a fallback anchor set; the implementation records how many epochs the fallback stays active. Figure 4.2 places the two losses in the training loop: the trajectory loss Ltraj pins the orbit to the data, the Floquet loss LF pushes the surrogate ρ̃T (Θ) below ρ∗ , and the gradient of their sum, taken by backpropagation through the solver, updates Θ until the surrogate is below the threshold; the spectrum is computed afterwards. 4.4. The deployed architecture: an orbital guarantee. Theorem 4.13 governs an autonomous learned field, and Lemma 4.15 carries it to the encoded field only while the oscillation η is small. The models we deploy are not in that regime, as Section 5.4 reports. This subsection settles what the encoded architecture does obey. We first show that the obstruction is structural rather than a defect of training. For a field that is exactly T̂ -periodic, one-period closeness on a set containing Γ and a contraction of the stroboscopic map on that same set are quantitatively incompatible, in the norm in which the contraction holds (Proposition 4.20). We then prove the guarantee that survives, which is orbital rather than trajectory-wise and whose hypotheses are exactly the two quantities Algorithm 4.1 already returns (Theorem 4.22). Throughout this subsection fΘ ∈ C 2 (Rd × R; Rd ) satisfies ∥fΘ ∥C 2 ≤ B and is exactly T̂ -periodic in time, fΘ (·, t + T̂ ) = fΘ (·, t), as the encoding (4.2) enforces; ΦΘ is the flow launched at phase 0. Write SΘ := ΦΘ (T̂ ; ·) for the stroboscopic map at phase 0. Periodicity gives the semigroup identity
ΦΘ (nT̂ + s; x) = ΦΘ s; SΘn x , n ∈ N, s ≥ 0,

(4.11)

(4.12)

LONG-TIME APPROXIMATION VIA SA-NODES

21

base point p ∈ Γ integrate fΘ over [0, T̂ ] true orbit γ̂|[0,T̂ ] returns PΘk p

Ltraj

Γ

LF : ρ̃T (Θ) update Θ

Σ

Ltraj + LF

p

backpropagation ∇Θ L Θ ← Θ − α∇Θ L true orbit γ̂ SA-NODE fΘ

training target: ρ̃T (Θ) ≤ ρ∗

(a) geometric setting

(b) optimization loop

Figure 4.2. Floquet-regularized training. (a) Geometric setting: the base-point orbit is the training target, and the returns PΘk p to Σ illustrate the contraction ρΘ < 1. (b) Optimization loop: Ltraj pins the base-point orbit to the data, LF penalizes ρ̃T (Θ) > ρ∗ , and the monodromy spectrum is computed post hoc (Remark 4.19). proved by induction on n: for n = 1 the function z(s) := ΦΘ (T̂ + s; x) solves ż = fΘ (z, T̂ + s) = fΘ (z, s) with z(0) = ΦΘ (T̂ ; x) = SΘ x, so uniqueness gives ΦΘ (T̂ + s; x) = ΦΘ (s; SΘ x); evaluating this at s = T̂ advances the induction. Note that Φf (T̂ ; ·) is the identity on Γ, because Γ is a periodic orbit of f of period exactly T̂ . The target’s stroboscopic map therefore fixes a whole curve, while a contracting SΘ fixes one point. That contrast is the content of the next proposition. Proposition 4.20 (contraction obstructs one-period closeness). Let Γ be a T̂ -periodic orbit of (4.1) and let fΘ be exactly T̂ -periodic. Suppose there are a norm ∥·∥⋄ , a set D ⊇ Γ and a constant q < 1 with ∥SΘ x − SΘ y∥⋄ ≤ q ∥x − y∥⋄ for all x, y ∈ D, and suppose SΘ has a fixed point x†Θ ∈ D. Then sup ΦΘ (T̂ ; x) − Φf (T̂ ; x) ≥ 12 (1 − q) diam⋄ (Γ). (4.13) x∈Γ

⋄

In particular, since Γ ⊆ N2r1 and T̂ ∈ [0, T̂ + 2], if c1 ∥·∥ ≤ ∥·∥⋄ ≤ c2 ∥·∥ then the constant ε of Assumption 4.4 satisfies ε ≥ 21 (c1 /c2 )(1 − q) diam(Γ), and ε ≥ 12 (1 − q) diam(Γ) when the contraction holds in the Euclidean norm. Proof. The proof is three lines: the target’s stroboscopic map fixes every point of Γ, so one-period closeness forces SΘ to move every point of Γ by little, while a contraction pulls every point of Γ to within σ/(1 − q) of its fixed point. Write σ := supx∈Γ ∥SΘ x − x∥⋄ . Since Φf (T̂ ; x) = x for x ∈ Γ, the left side of (4.13) is exactly σ. For x ∈ Γ ⊆ D, x − x†Θ so

x − x†Θ

⋄

⋄

≤ ∥x − SΘ x∥⋄ + SΘ x − SΘ x†Θ

⋄

≤ σ + q x − x†Θ

⋄

,

≤ σ/(1 − q) for every x ∈ Γ, whence diam⋄ (Γ) ≤ 2σ/(1 − q). This is (4.13).

For the last claim, Assumption 4.4 bounds ∥ΦΘ (t; x) − Φf (t; x)∥ by ε for x ∈ N2r1 ⊇ Γ and t ∈ [0, T̂ + 2] ∋ T̂ , and ∥·∥⋄ ≤ c2 ∥·∥ converts (4.13) into a bound on ε after dividing by c2 and using diam⋄ ≥ c1 diam. □

22

Z. LI AND N. M. MATZAKOS

Proposition 4.20 assumes a contraction on a set containing Γ, which is more than the ablation of Section 5.4 measures: there only the spectral radius at the fixed point is computed. The next proposition removes that gap. Its hypothesis is exactly the measured quantity, and its conclusion is the same incompatibility, at the cost of a square root and a d-th root in the constants. Proposition 4.21 (local form of the obstruction). Let Assumption 4.1 hold, so that the confinement radius r1 of Lemma 4.7 is defined, and let ∥f ∥C 2 (Nδ ) ≤ B; let fΘ be exactly T̂ -periodic with 0

∥fΘ ∥C 2 ≤ B, and let ε be the constant of Assumption 4.4. Let x♯ satisfy dist(x♯ , Γ) ≤ ϱ and put M := DSΘ (x♯ ). There are Λ1 , Λ2 ≥ 1, depending only on (d, T̂ , B), such that if ε ≤ r12 Λ2 /2 then h p

d−1 i1/d ρ(M ) ≥ 1 − Λ1 ϱ + 2εΛ2 1 + ∥M ∥ . (4.14) √ Consequently, for every ρ∗ ∈ [0, 1), ρ(M ) ≤ ρ∗ forces Λ1 ϱ + 2εΛ2 ≥ (1 − ρ∗ )d (1 + ∥M ∥)1−d . The proof is in Appendix A.7. Both forms say the same thing, and neither dominates the other. Proposition 4.20 is quantitatively sharp but assumes a contraction on a set containing Γ, which we do not verify. Proposition 4.21 assumes only the spectral radius at one point, which is what Algorithm 4.1 computes, but pays for it with Λ1 and Λ2 , which are governed by the C 2 size of the trained field and are not controlled during training. In the norm in which the contraction holds, (4.13) carries no constants: the tighter the contraction, the larger the one-period error on Γ must be, and the bound degenerates only as q → 1. Passing to the Euclidean ε of Assumption 4.4 costs the factor c1 /c2 , which we do not estimate; the two quantities the paper would like to be small at once are nonetheless in direct conflict. Section 5.4 measures both sides on every trained model. What survives is a guarantee about the orbit, and it needs neither Assumption 4.4 nor a section, a return time or a clock comparison. Theorem 4.22 (orbital guarantee for the deployed architecture). Let Assumption 4.1 hold and let fΘ be as above. Assume (B1) SΘ has a fixed point x†Θ whose monodromy M := DSΘ (x†Θ ) satisfies ρ(M ) ≤ ρ∗ < 1, where ρ∗ is a free parameter of this theorem and not the training threshold of (4.9); † † , Γ) ≤ εΓ . := {ΦΘ (s; x†Θ ) : s ∈ [0, T̂ ]} satisfies distH (γΘ (B2) the closed orbit γΘ Set ρ̄ := (1 + ρ∗ )/2 < 1. Then there exist r0 > 0 and C ≥ 1, depending only on (d, T̂ , B, ρ∗ ), such that for every x0 with x0 − x†Θ ≤ r0 and every t ≥ 0,
dist ΦΘ (t; x0 ), Γ ≤ εΓ + C ρ̄ ⌊t/T̂ ⌋ x0 − x†Θ . (4.15) Both hypotheses are measurable rather than structural: (B1) is the spectral radius that Algorithm 4.1 returns from one variational solve, and (B2) is a Hausdorff distance between two sampled curves. What the implementation returns are floating-point eigenvalues of a numerically integrated monodromy and a distance between finite samples, without residual or quadrature bounds, so the reported values are measurements of the hypotheses and not certificates of them. The radius r0 is likewise not evaluated for the initial conditions used in Section 5.4. The standard route to close this gap is validated integration of the monodromy together with interval bounds on the sampled distances, which would upgrade both measurements to certificates; we do not pursue it here. Neither involves ε, so Proposition 4.20 does not obstruct them, and Section 5.4 reports both. The bound is uniform in time with no linear term: the clock mismatch of Lemma 4.12, which produces the factor 1 + t/T̂ in (4.6), is absent here because both the target and the learned orbit close in exactly T̂ . Remark 4.23 (what the encoding buys, and what it gives up). Theorem 4.22 bounds the distance to Γ, not the distance to the target trajectory through the same initial state, and Proposition 4.20

LONG-TIME APPROXIMATION VIA SA-NODES

23

shows that the stronger conclusion is incompatible with a contraction of SΘ on a set containing Γ. The reason is visible in (4.12): Φf (T̂ ; ·) fixes every point of Γ, so the target retains the phase of its † initial condition forever, whereas a contracting SΘ entrains every nearby state to the single orbit γΘ and forgets that phase. The periodic encoding buys orbital stability with no data at deployment, and pays for it in phase fidelity. Section 5.4 exhibits both halves by varying the launch phase. Both terms of (4.15) are moreover of the correct order. In the polar setting of Remark 4.14, let
the learned field be ṙ = −(r − 1 − εΓ ), θ̇ = 2π/T̂ + κ sin 2πt/T̂ − θ with κ > 0: an exactly T̂ -periodic field whose flow locks to the circle of radius 1 + εΓ traversed in phase with the drive. In the co-rotating variable ψ := θ − 2πt/T̂ the system reads ṙ = −(r − 1 − εΓ ), ψ̇ = −κ sin ψ, so the locked point is a fixed point of SΘ with monodromy spectrum {e−T̂ , e−κT̂ }. Hypothesis (B1) therefore holds with ρ∗ = max(e−T̂ , e−κT̂ ), and (B2) holds with Hausdorff distance exactly εΓ . Every nearby trajectory converges to the locked circle, so its distance to Γ tends to εΓ : the time-independent floor in (4.15) is attained in the limit and cannot be removed, while the transient term decays geometrically at the rate the monodromy prescribes. 5. Numerical experiments This section tests each theoretical claim on a dedicated experiment, including a width sweep for the budget of Theorem 3.5 and an autonomous arm for the regime of Theorem 4.13. The cold-started autonomous arm is a negative result: training there produces no certifiable closed orbit. A warm-started arm, fine-tuned from the converged encoded models, repairs it (Section 5.4). Section 5.2 treats the data-assisted composite bound of Theorem 3.5 on a forced Duffing equation. Section 5.3 works on the pendulum: it exhibits the exponential barrier, shows the drift of Proposition 3.7 when the resets are withheld, and measures the growth of the window count in the horizon and the tolerance. Section 5.4 measures the monodromy spectrum of the deployed architecture, which is hypothesis (B1) of Theorem 4.22 and not a verification of (4.5), through an ablation over the time encoding and the Floquet loss, and Section 5.5 compares the two strategies on the van der Pol oscillator. The implementation, the configuration and archived results of every run, the trained models, and the scripts generating every figure and table are available at https://github.com/DCN-FAU-AvH/SA-NODE-MPC-Floquet.git. 5.1. Setup. We record the shared configuration once, together with the points at which the computations depart from the hypotheses of the theory. The experiments are proof-of-concept by design: low-dimensional canonical benchmarks on which each theoretical mechanism can be tested in isolation and every reported quantity can be recomputed independently; they are not application-scale validation. All models are trained in PyTorch with the Adam optimizer. The MPC experiments of Sections 5.2–5.3 use ReLU activation, per-window width P = 128, and forward Euler integration at the data resolution ∆t = 0.1; these choices follow the numerical protocol of [16, §4]. Theorem 2.1 is existential in the parameters and bounds the exact SA-NODE flow rather than its forward Euler discretization, and the Duffing field of (5.1) is only locally Lipschitz, so the per-window budget is read as a scaling guide rather than a guaranteed count. The Floquet experiments use tanh activation, width 64, and fixed-step RK4, consistent with the C 2 hypotheses of Section 4. Reference trajectories are computed by an adaptive Runge–Kutta method at integrator tolerances several orders of magnitude below the errors reported. Monolithic baselines use the width of a single window and the same total number of gradient steps as the composite, so the comparisons are matched in training cost rather than capacity. Theorem 2.1 and the width scaling of Section 2.2 certify a sufficient width for an equal-parameter monolithic model that still carries the factor e2LT ; they do not assert that the model must incur it. 5.2. Composite approximation of a forced Duffing equation. The first experiment instantiates the data-assisted bound (1.3) of Theorem 3.5: it runs Algorithm 3.1 on the forced Duffing

Z. LI AND N. M. MATZAKOS

2.0

2.00

1.5

1.75

1.0

1.50

0.5

1.25

sup ke(t)k x

0.0

0

x1 (t)

24

1.00

0.5

0.75

1.0

0.50

1.5

0.25

2.0

0.00

0

2

4

6 true

t

8

10

12

14

0

2

4

composite

(a) sample trajectories

6

10

12

monolithic SA-NODE composite

ε

t

8

14

(b) error over time

Figure 5.1. Forced Duffing equation (5.1), T = 15, ε = 0.5. (a) Sample trajectories, true versus composite; the dotted vertical line marks the switch time. (b) Error over time, supx0 ∥e(t)∥ on the initial-condition grid; the sawtooth peaks are left limits at the switch times. equation ẋ1 = x2 ,

ẋ2 = x1 − x31 + 0.1 cos(πt),

(5.1)

a non-autonomous target that activates the time bias of (1.1), on [0, 15] with a 5 × 5 grid of initial conditions in [−1.5, 1.5]2 , tolerance ε = 0.5, and training horizon H = 10. The algorithm selects 2 windows, of lengths 10.4 and 4.6, and realizes the composite error 0.519, which exceeds the tolerance by 3.8 %: the safeguards bind, so Assumption 3.4, checked on the training grid, holds in place of ε; the sampled analogue of (3.7) then holds with that value. Meanwhile with maxk εwin k the monolithic SA-NODE baseline, trained with the same number of gradient steps, reaches 1.98, about 4.0 times the tolerance. Figure 5.1(a) displays sample trajectories of the composite against the truth, and Figure 5.1(b) the error over time. The data reset holds the composite near the tolerance on every window, while the monolithic error exceeds it from early times onward by a factor of about 4.0 and oscillates at the scale of the reachable tube. At each switch time the reported error is the left limit: the composite is reset to the data there, so the sawtooth peaks just before the switch. On a bounded target the barrier appears as a persistent excess rather than unbounded growth; the growth mechanism itself is exhibited on the pendulum below. 5.3. Window count, tolerance, and the exponential barrier. This experiment has three purposes: to exhibit the exponential error growth of Section 2.2 on a monolithic model, to show the drift that Proposition 3.7 predicts when the data resets are withheld, and to measure how the window count of Theorem 3.5 grows with the horizon and the tolerance. The benchmark is the pendulum ẋ1 = x2 , ẋ2 = − sin x1 , run with the protocol of Section 5.2 on a 5 × 5 grid of initial conditions in [−2, 2]2 ; the horizon sweep uses training horizon H = 3, and the tolerance sweep trains each window on the full remaining horizon. On R2 eight of these initial conditions carry energy above the separatrix value 2, so their trajectories rotate and the infinite-time reachable set is unbounded there. Read on the cylinder the reachable set is compact, and all horizons run here are finite; Assumption 3.4 is in any case verified on the output, so the sweep instantiates the reset mechanism rather than the literal infinite-time tube hypothesis.

LONG-TIME APPROXIMATION VIA SA-NODES

25

Table 5.1. Window count and composite error versus horizon, pendulum, ε = 0.5. T

5

10

15

20

N errL∞

3 8 15 23 0.58 0.58 0.72 0.72

Figure 5.2(a) shows the error growth at T = 10. The monolithic neural ODE error grows exponentially along the horizon, at the measured rate e0.39 t (correlation r = 0.88 over the growth phase), before saturating at the scale of the reachable tube. The composite is reset at each switch time and stays at the tolerance. Its excess over ε (realized error 0.583, about 16.6 % above the tolerance) is the safeguard excess anticipated after Assumption 3.4. The same figure isolates the value of the reset. We take the eight trained windows of this run and chain them without data: each window starts from the state the previous one predicted, which is the predicted-IC deployment of Proposition 3.7. No retraining is involved, so the two composites differ only in what happens at the switch times. The error now accumulates across windows and reaches 3.74, that is 7.5 times the tolerance, while the data-IC evaluation of the same models stays at 0.58. The bound of Proposition 3.7 is honored but far from sharp: the weight-norm Lipschitz bound of the trained ReLU windows is L̄ = 29, which puts the envelope near 1089 . The contrast between the two curves is the content of the proposition: the resets, not the trained windows, carry the uniform guarantee. The partitions produced here satisfy the mesh condition of Assumption 3.4 with a horizonindependent bound: across the sweep the realized mesh is 3.00 at every horizon, so τmax can be taken independent of T as the theory requires. Figure 5.2(b) and Table 5.1 report the window count against the horizon T ∈ {5, 10, 15, 20} at ε = 0.5. The least-squares fit N ≈ 1.34 T − 4.50 has correlation r = 0.995, consistent with the linear growth allowed by Theorem 3.5, and far below the exponential width that (2.1) certifies as sufficient for a monolithic model. Two caveats temper this reading. The increments of N are 5, 7, and 8 and mildly convex. The realized error also drifts from 0.58 to 0.72 across the sweep as the safeguards bind more often, so four points do not separate linear from mildly super-linear growth. Every entry exceeds ε = 0.5: the tolerance is met per window only up to that excess. Problem 6.2 records the general question. The width side of the budget is probed separately: a sweep over widths 16 to 256 on a single fixed window, three seeds each, gives a realized error decaying with slope −0.75 against the width on a log-log scale (|r| = 0.97), consistent with, indeed faster than, the P −1/2 rate that Theorem 2.1 certifies. Figure 5.2(c) reports the tolerance sweep at T = 10: the count N = (4, 8, 12, 15, 15) for ε ∈ {1.0, 0.7, 0.5, 0.35, 0.25}, with each window trained on the full remaining horizon and window cap S = 15. At ε = 0.35 and 0.25 the cap binds, producing the admissibility failure anticipated after Assumption 3.4. The crossing fires so early that the minimum window length dominates, and the horizon is left uncovered from t = 9.1 and t = 9.0 onward in the two runs, where the composite is not defined and no error can be reported. Assumption 3.4 fails outright, and Theorem 3.5 does not apply. In this regime the remedy is to add width per window rather than more windows. 5.4. Certification: time encoding and Floquet loss. Corollary 4.18 links the trained surrogate to the certificate (4.5) only for an autonomous planar field. The architecture ablated here is periodic in time, and by Remark 4.19 the surrogate then controls det M rather than ρ(M ), so the experiment is scored on the latter. It measures which ingredients produce contraction. The ablation is two-by-two: the periodic encoding (4.2) against the raw linear time bias of (1.1), and training with LF against training without it. A third, autonomous arm removes the time features altogether, which is the regime of Theorem 4.13 and Corollary 4.18; it is reported at the end

0

100

10 1 0

2

4

t

monolithic SA-NODE monolithic NODE composite, predicted-IC

6

8

10

composite, data-IC ∝ e 0.39 t

ε

(a) error over time, T = 10

window count N

window count N

Z. LI AND N. M. MATZAKOS

sup ke(t)k x

26

20 15 10 5 5

10

15

horizon T

20

15 10 5 1.0

0.6

tolerance ε

0.4

adaptive N(ε) window cap S = 15 cap binds: horizon uncovered

adaptive N(T) fit N ≈ 1.34 T − 4.50 (b) window count vs. horizon

0.8

(c) window count vs. tolerance

Figure 5.2. Pendulum sweeps. (a) Error over time at T = 10, logarithmic scale; the dashed guide is the least-squares fit e0.39 t to the neural ODE baseline. (b) Window count versus horizon at ε = 0.5, with least-squares fit. (c) Window count versus tolerance at T = 10; the dotted line is the window cap S, which binds at the two smallest tolerances and leaves the horizon uncovered from t ≈ 9: an admissibility failure, not a budget limit. of this subsection. The two benchmarks have true multipliers an order of magnitude apart: the Stuart–Landau system ṙ = r (0.1 − r2 ), θ̇ = 1, (5.2) with ρtrue = e−0.4π ≈ 0.285 and threshold ρ∗ = 0.30, and the van der Pol oscillator ẋ1 = x2 ,

ẋ2 = µ (1 − x21 ) x2 − x1 ,

µ = 0.5,

(5.3)

with period T̂ ≈ 6.3807, ρtrue ≈ 0.039, and the strict threshold ρ∗ = 0.05. For each system, encoding, and loss condition (Ltraj only, condition A, or Ltraj + LF , condition B), we train 15 models, one per seed, and run each of them autonomously for 300 periods from an initial condition at distance 0.2 from the orbit. The linear-encoding models reset their clock at each period, the deployment protocol of the original architecture. The reset makes the field piecewise-defined in time and in general discontinuous at the period boundary, hence outside the C 2 -periodic class of Theorem 4.22; this arm is an empirical ablation only. Every model is assessed by the spectral radius ρ(M ) of the full one-period stroboscopic monodromy along its locked orbit, obtained from one variational solve; for a T̂ -periodic field this is what governs the linear stability of that orbit (Remark 4.19). It is neither the return-map radius ρΘ of Assumption 4.6 nor the surrogate ρ̃T = det M of (4.10), which is reported separately. Table 5.3 reports the outcome, Figure 5.3(a) the per-seed values, and Figure 5.3(b) the resulting long-horizon errors. Three readings must be kept apart. The first is whether the locked orbit is linearly stable, that is whether ρ(M ) < 1, which is hypothesis (B1) of Theorem 4.22. It does, in 119 of the 120 runs. The single exception is a linear-encoding van der Pol model trained without LF , with ρ(M ) = 1.116, and it is also the only divergent deployment, reaching error 51 after 300 periods. Linear stability of the locked orbit is therefore the rule rather than the exception. This measurement supplies hypothesis (B1) and no more: it is not the certificate (4.5), which concerns the averaged field and is out of reach at the measured oscillation η, and it leaves Assumption 4.4 unverified, since no training loss controls closeness on the whole tube. The second is the contraction rate, and there the encoding is decisive. With the periodic encoding the median spectral radius is 0.044 on Stuart–Landau and 0.101 on van der Pol; with the

LONG-TIME APPROXIMATION VIA SA-NODES

27

raw linear bias it is 0.532 and 0.695. The encoding thus improves the rate by roughly an order of magnitude on the first system and a factor of about five to seven, depending on the loss condition, on the second. The Floquet loss tightens the rate further on the unprotected architecture, from 0.695 to 0.387 on van der Pol, and it removes the one unstable run. The third is whether a model reaches the threshold ρ∗ used in the loss; the per-cell counts are in Table 5.3. Every Stuart–Landau model with the periodic encoding meets ρ∗ = 0.30, while on van der Pol few models reach ρ∗ = 0.05: their spectral radii sit at two to three times the target’s own multiplier ρtrue ≈ 0.039, from which the threshold was set. The learned orbits contract, but by less than the true orbit does. Scored on det M instead, every periodic-encoding van der Pol seed would pass while few pass on ρ(M ); the Stuart–Landau linear cells show the same gap (Remark 4.19). We do not adjust ρ∗ after the fact; it is a training target inherited from the surrogate. The orbital reading is the one Theorem 4.22 governs, and its two hypotheses are measured on all 60 periodic-encoding models. Hypothesis (B2), the Hausdorff distance εΓ between the locked orbit and the target cycle, has median 7.9 × 10−4 on Stuart–Landau and 1.5 × 10−2 on van der Pol; hypothesis (B1) is the spectral radius already reported. The constant of Assumption 4.4, evaluated over the tube rather than at the base point, is of another order: median 0.88 and 4.84, that is 1124 and 328 times εΓ , above the diameters 0.63 and 4.74 of the target cycles. This is not a training failure but the obstruction at work. Proposition 4.21 explains the mechanism under exactly the measured hypotheses, ρ(M ) at the fixed point and ϱ ≤ εΓ . For the magnitude, the operator norm of the Jacobian of the stroboscopic map over the 0.2-tube, sampled on a grid for the two representative models, is at most 0.052 on Stuart–Landau and 0.186 on van der Pol; at these contraction factors Proposition 4.20 gives the lower bounds ε ≥ 0.30 and ε ≥ 1.93, and the measured tube errors lie above both. The grid values are a finite sampling, not a certificate, so this is a consistency check: on these models one-period closeness on the tube, and with it Theorem 4.13 and Lemma 4.15, is not observed at any tuning. The launch phase separates the two readings. We deploy the same threshold-passing model from the same state at distance 0.2 from the cycle, moving the launch phase from the training phase to half a period later. The supremum of the trajectory error over eight periods rises from 0.21 to 0.63 on Stuart–Landau, and from 0.26 to 4.90 on van der Pol. The distance to the target cycle stays at the εΓ level: 5.5 × 10−4 against 3.0 × 10−4 , and 1.0 × 10−2 against 3.3 × 10−3 , at the end of the run. The trained models track the cycle and not the trajectory, exactly as Remark 4.23 describes, and the earlier protocol did not see it because it launched only at the training phase. The transient of (4.15) is also visible: before saturating at the εΓ floor, the distance to the cycle contracts per period by a measured factor 0.036 on Stuart–Landau and 0.044 on van der Pol, of the order of the measured ρ(M ) and well below ρ̄. The autonomous arm fails, in the way the encoding was designed to prevent. In none of the 10 runs does training produce a certifiable closed orbit: every model relaxes to a stable equilibrium instead of a periodic attractor, and with no nondegenerate closed orbit no periodic-orbit certificate can be defined. The one-period tube errors range from 0.43 to 5.1. This is consistent with the saturation obstruction [20] that motivated the periodic encoding, and it leaves the chain of Corollary 4.18 without a trained positive instance under cold-started autonomous training. A warm-started arm repairs this failure. Transferring the weights of the converged periodicencoding models to the autonomous architecture (the shared matrices verbatim, the constant bias as the phase average of the periodic bias branch) and fine-tuning under the same condition-B loss produces closed learned orbits where cold-started training produced none. On no seed does the transferred field carry a closed orbit before fine-tuning, so the repair comes from the warm-started optimization, not from the initialization alone. On van der Pol, four of five seeds yield a closed, transversally stable orbit with ρ̃T between 0.14 and 0.19 and transverse multiplier ρ(DPΘ (x∗Θ )) between 0.13 and 0.18; the surrogate–multiplier gap of Corollary 4.18(ii) is 0.002–0.011 against a one-period tube error ε of 0.61–0.78, smaller by factors of 65 to 320. On Stuart–Landau every seed

28

Z. LI AND N. M. MATZAKOS

Table 5.2. Warm-started autonomous arm, per seed: divergence surrogate ρ̃T on the relaxed orbit, transverse multiplier ρ(DPΘ ), the absolute gap |ρ̃T − ρ(DPΘ )| (the quantity of Corollary 4.18(ii), computed from unrounded values), one-period tube error ε, and fitted stroboscopic growth per period. Van der Pol seed 2 relaxes to an equilibrium and is excluded; on Stuart–Landau the relaxation closure test does not pass, so ρ̃T is not evaluated there. System

seed

ρ̃T

ρ(DPΘ )

gap

ε

slope

Stuart–Landau 0 1 2 3 4

— — — — —

0.660 0.654 0.656 0.660 0.654

— — — — —

0.128 0.136 0.132 0.128 0.136

0.108 0.107 0.107 0.107 0.107

van der Pol

0.144 0.151 0.168 0.193

0.134 0.143 0.170 0.183

0.010 0.008 0.002 0.011

0.782 0.606 0.623 0.692

0.103 0.089 0.157 0.108

0 1 3 4

yields a closed, transversally stable orbit with ρ(DPΘ ) ≈ 0.66 and tube error ≈ 0.13; the closure test of the surrogate’s relaxation does not pass there, so ρ̃T is not evaluated. An independent recomputation with an adaptive integrator reproduces every reported multiplier to within 10−3 and closes the learned orbits to residuals below 10−10 on van der Pol and 3 × 10−5 on Stuart– Landau: the orbits are genuine, and the failed closure test reflects the strictness of the check, not the absence of a cycle. At deployment the stroboscopic error grows at 0.09–0.16 per period before the phase wraps, consistent with the envelope Cε(1 + k) of Theorem 4.13(iii). The arm is a measured instantiation of the certification mechanism of Corollary 4.18, not a certified instance: the smallness conditions and constants are not evaluated, and the contraction is read from the measurement rather than from the training target ρ∗ . Two further measurements complete the picture. The first is the oscillation quantified by Lemma 4.15. On the tube of radius 0.2 the representative condition-B models have C 1 oscillation 0.46 (Stuart–Landau) and 2.75 (van der Pol). The lemma requires ε + CB η ≤ δ0 /4; already for B = 1 the constant CB is of order 109 , so the admissible η is below 10−10 and the measured values miss the regime by ten orders of magnitude. Deployment succeeds because the trajectories stay locked to the periodic drive, not because the field is near-autonomous. An oscillation penalty on the loss reduces η by a factor of five to seven at the cost of a coarser one-period fit, and cannot close a gap of ten orders. The second is the reference computation behind Remark 4.19, carried out on the representative models rather than on the whole ablation. Without the oscillation penalty, the Stuart–Landau model reports ρ̃T = 3.4 × 10−4 while its multipliers are (0.048, 0.007); the van der Pol spectral radius 0.118 already exceeds ρ∗ = 0.05 although its determinant 0.014 passes. Under the penalty the van der Pol determinant 0.155 stays below one while one multiplier lies above it, at 1.256. That is the failure Remark 4.19 warns of, and it is the reason the ablation above is scored on ρ(M ) rather than on ρ̃T . 5.5. Comparison of the two strategies on the van der Pol oscillator. The final experiment places both strategies on the van der Pol oscillator (5.3) over eight periods under one protocol. Their guarantees are not directly comparable (data-assisted against autonomous); the comparison is of operating characteristics, not of sample efficiency. Each method receives one training trajectory, but not the same one. MPC trains on the eight-period trajectory from x0 , at distance 0.2

LONG-TIME APPROXIMATION VIA SA-NODES

29

Table 5.3. Spectral radius ρ(M ) = maxi |λi | of the one-period stroboscopic monodromy of the encoded flow (Remark 4.19), over 15 seeds per cell of the encodingby-loss ablation. The column <1 counts the models whose locked orbit is linearly stable; the column ≤ρ∗ counts those reaching the threshold used in the Floquet loss. Medians are over the seeds of the cell. System

encoding

without LF

ρ∗

med. ρ(M ) Stuart–Landau periodic linear Van der Pol periodic linear

0.30 0.05

0.044 0.532 0.101 0.695

<1

with LF ≤ρ∗

med. ρ(M )

15/15 15/15 15/15 0/15 15/15 1/15 14/15 0/15

0.040 0.483 0.080 0.387

<1

≤ρ∗

15/15 15/15 15/15 0/15 15/15 4/15 15/15 0/15

102 101 100 ke(kT̂ )k

ρ(M) = max |λi | i

100

10 1

10 1

10 2 10 3

10 2

10 4 SL periodic

SL linear no LF with LF

VdP periodic ρ=1

(a) per-seed values

ρ∗

VdP linear

0

50

100

150

period k

periodic, with LF , median periodic, with LF , worst linear, no LF , median

200

250

300

linear, no LF , worst linear, with LF , median linear, with LF , worst

(b) stroboscopic error, van der Pol

Figure 5.3. Certification ablation. (a) Per-seed monodromy spectral radius ρ(M ), by system, encoding, and loss; dashed lines are the thresholds ρ∗ , the dotted line is ρ = 1. (b) Stroboscopic error across seeds on the van der Pol oscillator: the periodic encoding stays bounded, the linear-encoding worst case diverges. from the cycle; Floquet trains on a single period from the base point p on the cycle. MPC uses tanh activation and width 64 (Algorithm 3.1, H = 5, ε = 0.5) and is evaluated in data-IC mode; Floquet uses condition B and runs autonomously. Table 5.4 and Figure 5.4 summarize the comparison. MPC selects 9 windows and reaches supt ∥e(t)∥ = 0.537, just above the tolerance ε = 0.5, the excess being the safeguard excess of Section 5.2. This run uses tanh rather than the ReLU that Assumption 3.2 fixes, so only the error half of Theorem 3.5, verified through Assumption 3.4 on the output, is claimed for it. The sawtooth of Figure 5.4(a) reflects the external reset: the error restarts near zero at each switch time and ends the run at 8.8 × 10−3 , the final window’s own error. The locked orbit of the Floquet model is linearly stable, with monodromy spectral radius ρ(M ) = 0.106, and the model runs open loop on the periodic clock, with no state observations. As in Section 5.4 it does not reach the threshold ρ∗ = 0.05 inherited from the surrogate. Its two multipliers form a complex conjugate pair of equal modulus, so the surrogate returns ρ̃T = det M = ρ(M )2 = 0.011. This is the gap of

30

Z. LI AND N. M. MATZAKOS

Table 5.4. Comparison on the van der Pol oscillator (µ = 0.5, eight periods, width 64, tanh).

models trained evaluation supt ∥e(t)∥ ∥e(8T̂ )∥ ρ(M ) ρ̃T = det M

MPC–SA-NODE

Floquet–SA-NODE

9 data-IC resets 0.537 8.8 × 10−3 – –

1 observation-free, time-periodic 0.269 0.169 0.106 0.011

100

2

100

10 2

1

ke(t)k

10 2

10 6 0

10

20

t

30

40

50

MPC (9 windows, data-IC) Floquet (ρ̃T = 0.011, autonomous) ε

(a) continuous error

x2

ke(kTLC )k

10 4

0

10 4

1

10 6

2 0

2

4

period k

MPC

6

8

Floquet

(b) stroboscopic error

2

1 Γ

true

0

x1

1

MPC Floquet

2 x0

(c) phase portrait

Figure 5.4. Comparison on the van der Pol oscillator. (a) Continuous error: the MPC sawtooth against the Floquet transient followed by saturation. (b) Stroboscopic error at the period marks. (c) Phase portrait: both learned trajectories shadow the limit cycle Γ. Remark 4.19 in its plainest form: the surrogate is the square of the quantity that governs stability, and reading 0.011 as a contraction rate would understate it by an order of magnitude. Its error peaks at 0.269 during the initial transient and settles to 0.169 at the period marks (Figure 5.4(b)). The deployed encoded model is governed by the uniform bound (4.15), and the internal reset corrects toward the orbit rather than toward the reference trajectory. Both trajectories shadow the limit cycle (Figure 5.4(c)). For horizons of K periods the MPC strategy trains one network per window, so its budget grows with the realized window count N , about one window per period in the comparison experiment, while the Floquet cost is that of a single training run; the dataversus-autonomy trade-off is discussed in Section 3.5. 6. Conclusions and perspectives We first summarize what has been established and what its limits are, and then state three open problems. 6.1. Conclusions. This paper studied how the certified error bound of a learned flow deteriorates with the horizon T and how that deterioration can be stopped. The barrier is a property of the explicit constant (2.1) for a single network trained once on [0, T ]: Grönwall amplification of the field error and growth of the reachable tube, both attained on the linear system ẋ = Lx (Section 2.2). The model predictive strategy resets the state with data. Provided training realizes the tolerance ε on every window, the run covers the horizon, and the reachable tube is bounded with

LONG-TIME APPROXIMATION VIA SA-NODES

31

uniformly bounded data (Assumption 3.3), the composite keeps the error at ε uniformly in time, up to the quantified safeguard excess, with a parameter budget proportional to the window count (Theorem 3.5). The reset confines the Grönwall amplification to a single window but does not remove the growth of the tube. Run from its own predictions instead of the data, the same composite obeys only the exponential bound of Proposition 3.7: the resets, not the trained windows, carry the uniform guarantee. The Floquet strategy replaces data by stability. For an autonomous learned field whose return map is certified to contract (Assumption 4.6), one-period accuracy ε on a tube around the cycle propagates to the global bound Cε(1 + t/T̂ ): transverse errors are contracted at every return, and only the phase drift accumulates (Theorem 4.13, Remark 4.14). For the time-periodic architecture we deploy, the scalar certificate degenerates to the determinant of the monodromy (Remark 4.19), and strong stroboscopic contraction even forces a lower bound on the one-period error (Propositions 4.20 and 4.21). What survives is the orbital guarantee of Theorem 4.22, uniform in time, whose two hypotheses are measured on the trained model: the periodic encoding buys orbital stability and gives up phase fidelity. The experiments confirm this picture. The composite meets its tolerance up to the safeguard excess while a monolithic baseline of matched training cost exceeds it about fourfold, the window count grows linearly with the horizon, and withholding the resets produces the predicted drift. In the certification ablation the periodic encoding is the dominant factor setting the contraction rate, and both hypotheses of Theorem 4.22 are measured on every periodic-encoding model. Coldstarted autonomous training yields no certifiable orbit; warm-starting from the encoded models repairs it, with a surrogate–multiplier gap far below the measured ε (Corollary 4.18(ii)). Longhorizon accuracy is thus governed by the mechanism that removes accumulated error, not by network size: the width certified by (2.1) grows double exponentially in T , while the composite adds fixed-width networks one window at a time and the Floquet model runs at a single fixed width. 6.2. Perspectives. Three problems mark the most important continuations: one on the barrier itself, one on the window count of the model predictive strategy, and one on certification in higher dimension. Problem 6.1 (a lower bound for the barrier). Exhibit a Lipschitz constant L, a compact K, and a family of targets f , L-Lipschitz in x uniformly in t and with a fixed C 2 bound on the reachable set, such that every SA-NODE (1.1) of width P whose flow satisfies sup

b Θ (t; x0 ) − Φ(t; x0 ) Φ

≤ ε

x0 ∈K, t∈[0,T ]

obeys P ≥ ϕ(T, ε) with ϕ super-linear in T at fixed ε; or prove that no such family exists. By 2 Lemma 2.2 it suffices, within the subclass whose p fields obey a fixed C bound B, to bound below the width needed to approximate f within 2 εeB C2 uniformly on the spacetime reachable set by a shallow network in (x, t). The restriction to a fixed B is not cosmetic: for unrestricted parameters B may grow with the width, and the reduction then says nothing. A super-linear answer would make the barrier intrinsic and a decomposition of the horizon necessary rather than merely convenient; a linear answer would make the composite of Theorem 3.5 optimal in T up to constants and confine the exponential growth to the certificate. Problem 6.2 (window count). Prove, or refute, that the partition generated by the unsafeguarded crossing rule (3.4), with smin = 0 and S = ∞, satisfies N ≤ cT + C at fixed ε for Lipschitz targets, with c controlled by the local Lyapunov exponents along the trajectory rather than by the sampling step. The problem is posed with a fixed bound on the window weights: without one, arbitrarily poor windows cross arbitrarily fast, and already for f = 0 infinitely many windows fit in a finite

32

Z. LI AND N. M. MATZAKOS

horizon. The safeguarded algorithm gives N = Θ(T ) trivially, with a constant set by the data resolution; Section 5.3 supports the linear law on the pendulum. Problem 6.3 (higher dimension). The Liouville–Abel identity certifies only in the autonomous plane. In state dimension d > 2 the transverse monodromy MΘ⊥ of the learned orbit, the linearization of its return map, is a matrix, and already for d = 2 a time-periodic field escapes the identity, since the divergence integral then returns the determinant rather than the spectral radius. Construct differentiable surrogates for the spectral radius ρ(MΘ⊥ ) that are tight enough to train against, beyond the Bauer–Fike route of Section 4.3. Appendix A. Proofs This appendix collects all proofs, one part per result. Appendix A.1 proves Lemma 2.2; Appendix A.2 proves the results of Section 3; Appendices A.3–A.7 prove the results of Section 4. A.1. Proof of Lemma 2.2. This part proves the reduction of the width lower bound from flows to fields. Proof of Lemma 2.2. Write x = Φ(t0 ; x0 ) with x0 ∈ K, and set g(h) := ΨΘ (t0 + h; t0 , x) − Ψ(t0 + h; t0 , x),

h ∈ [0, 1],

so that g(0) = 0 and ġ(0) = fΘ (x, t0 ) − f (t0 , x). Differentiating once more along the two flows gives g̈ = Dx fΘ fΘ + ∂t fΘ − Dx f f − ∂t f , whence ∥g̈∥ ≤ 2(B 2 + B) = C2 . b Θ (t0 ; x0 ), so that ∥x − ŷ∥ ≤ ε by hypothesis. The semigroup For the size of g, put ŷ := Φ b Θ (t0 + h; x0 ). Hence property gives Ψ(t0 + h; t0 , x) = Φ(t0 + h; x0 ) and ΨΘ (t0 + h; t0 , ŷ) = Φ ∥g(h)∥ ≤ ∥ΨΘ (t0 + h; t0 , x) − ΨΘ (t0 + h; t0 , ŷ)∥ +

b Θ (t0 + h; x0 ) − Φ(t0 + h; x0 ) Φ

≤ εeBh + ε ≤ 2εeB , the first term by Grönwall for the single field fΘ , whose Lipschitz constant in x is at most B, and the second by hypothesis. Grönwall is legitimate here by a continuity argument. The trajectory b Θ (·; x0 ) lies within ε of it by hypothesis; only ΨΘ (·; t0 , x) Ψ(·; t0 , x) lies on KT , and ΨΘ (·; t0 , ŷ) = Φ needs the argument. On the maximal interval where it stays in the r-neighborhood of KT , the display bounds its distance to KT by (1 + eB )ε ≤ r/2 < r, so it does not reach the boundary before h = 1 and the interval is all of [0, 1]. Taylor expansion at h = 0 Rwith the integral remainder, which is the form valid for a vectora valued g, gives aġ(0) = g(a) − 0 (a − s)g̈(s) ds, hence ∥ġ(0)∥ ≤ 2εeB /a + aC2 /2 for every a ∈ (0, 1]. p which lies in (0, 1] precisely under the standing The right side is minimized at a∗ = 2 εeB /C2 , p −B bound ε ≤ C2 e /4, and the minimum value is 2 εeB C2 . □ A.2. Proofs of Theorem 3.5 and Proposition 3.7. This part proves the two results of Section 3: the data-assisted composite bound and the predicted-IC bound. Proof of Theorem 3.5. The two conclusions have different status, and the proof makes it explicit. Two things are assumed and not proved: that the optimizer realizes the tolerance on each window, and that Algorithm 3.1 terminates having covered the horizon; they are the content of Assumption 3.4, verified a posteriori on each run and reported in Section 5. The theorem itself assumes the continuous supremum (3.6). A further discretization gap attaches only to the experimental verification: the suprema are checked on finite grids of initial conditions and times (the standing gap of Section 5.1), so what the runs verify is the sampled analogue of Assumption 3.4, not the continuous supremum the theorem assumes. Estimate (3.7) unpacks Assumption 3.4: the content of the data reset is that the global error is a maximum over windows and not a sum, so once each window meets the tolerance from the true state, so does the composite. The substance is

LONG-TIME APPROXIMATION VIA SA-NODES

33

Table A.1. Dependency chain of the radii and smallness thresholds of Section 4. Each row may shrink quantities from earlier rows; no entry depends on Θ, ε, t, or the initial condition. The constraint on ε3 anticipates the certificate-free constants δ1 and CU of Lemma 4.10(i). Quantity

Introduced in

Constraint

Depends on

r1 r2 , ε2 r3 , r3′ , ε3 r5′ , δ1 r5 r4 , ε4 r6 , ε5 ε0 , r0

Lemma 4.7 Lemma 4.8 Lemma 4.9 Lemma 4.10(i) Lemma 4.10(ii) Lemma 4.11 Lemma A.1 Theorem 4.13

(0, δ0 /4] r2 ∈ (0, r1 ] √ r3 , r3′ ∈ (0, r2 ], ε3 ≤ ε2 , 2C2 ε3 ≤ min(δ1 , r3 /(2CU )) ′ r5 ∈ (0, δ0 /4] min(r5′ , r3 /(2CU )) r4 ∈ (0, min(r3 , r3′ , r5 )], ε4 ≤ ε3 r6 ∈ (0, r5 /2], ε5 ≤ ε4 r0 := r6 (Nr0 ⊆ N2r1 ), ε0 := min{ε5 , r0 /K3 , v/2}

(f, γ̂, δ0 ) (f, γ̂, δ0 , B) above; ρ∗ via r3 , ρT (f ) via r3′ (f, γ̂, δ0 , B), certificate-free ρ∗ through r3 only all of the above all of the above (f, γ̂, δ0 , B, ρ∗ , ρT (f ))

(3.8), which uses no part of Assumption 3.4 beyond the mesh, and which is where Assumption 3.3 b Θ evaluated from the true state, so no upstream error is does its work. On Ik the composite is Φ k amplified and (3.7) is the maximum over k of the per-window bounds (3.6). For (3.8), Theorem 2.1 on window k has initial-condition set Φ(τk ; K) ⊂ K∞ and horizon at most τmax , and its constant obeys Cτmax ,Φ(τk ;K),f (·,·+τk ) ≤ Cτmax ,K∞ ,f , because window k is the term τ = τk of the supremum (3.5) and Assumption 3.3 makes the right side finite and independent of k and T . Hence Pk ≲ Cτ2max ,K∞ ,f ε−2 on every window, with a constant independent of k and T ; summing over k concludes. □ Proof of Proposition 3.7. We derive a one-step recursion for the error accumulated at the switch times and unroll it. Let Ek denote the accumulated error at the switch time τk and ℓk := τk+1 − τk . On Ik the upstream error Ek enters the learned flow and is amplified by at most eL̄ℓk ; the window adds its own training error at most ε, by (3.6) of Assumption 3.4. Hence Ek+1 ≤ eL̄ℓk Ek + ε,

E0 = 0.

Unrolling the recursion gives EN ≤ ε

N X k=1

eL̄(τN −τk ) = ε

N X

eL̄(T −τk ) ,

k=1

where the k = N term is the last window’s own ε; this is the first bound in (3.10), and between switch times the same amplification argument applies. For the second bound, with L̄ > 0, we use P −1 L̄mτmax T − τk ≤ (N − k)τmax termwise, and the geometric sum N evaluates to the stated m=0 e quotient. □ A.3. Supporting lemmas for the Floquet strategy. This part and the four that follow prove the results of Section 4. Here we prove Lemma 4.5, the six lemmas preceding Theorem 4.13, an auxiliary hitting estimate (Lemma A.1), and the scalar multiplier identity (Lemma A.2). Throughout Appendices A.3–A.5, Assumptions 4.1–4.6 are in force with η = 0 in Assumption 4.2, and we write fΘ (x); the proof of Corollary 4.18 uses only Assumptions 4.1–4.4, as noted there. Lemma 4.15, proved in Appendix A.6, removes the η = 0 restriction, and Appendix A.7 stands outside this convention altogether. All constants depend only on (f, γ̂, δ0 , B, ρ∗ , ρT (f )). Table A.1 records the dependency chain of the radii and smallness thresholds. Proof of Lemma 4.5. The bound is a second-order Taylor expansion of the flow difference, optimized in the expansion time. Since ε ≤ δ0 /4, Lemma 4.7 keeps both trajectories in Nδ0 , where

34

Z. LI AND N. M. MATZAKOS

the bounds B and L apply. Fix x ∈ Γ and set g(t) := ΦΘ (t; x) − Φf (t; x), so that g(0) = 0 and ġ(0) = fΘ (x, 0) − f (x). Differentiating once more along the two flows, g̈ = Dx fΘ fΘ + ∂t fΘ − Df f,

so

∥g̈∥ ≤ C2

by Assumption 4.4 and the C 1 bound on f . Taylor expansion at 0 gives, for every h ∈ (0, T̂ + 2], ∥ġ(0)∥ h ≤ ∥g(h)∥ + 12 C2 h2 ≤ ε + 21 C2 h2 .

p √ Choosing h = 2ε/C2 , which is admissible under the stated bound on ε, yields ∥ġ(0)∥ ≤ 2C2 ε. This bounds the mismatch at phase 0. For η > 0 the supremum defining δ also runs over t ∈ R, and (4.4) gives ∥fΘ (x, t) − fΘ (x, 0)∥ ≤ 2η for every t; adding this to the phase-0 bound and taking the supremum over x ∈ Γ gives the claim. □ Proof of Lemma 4.7. Both claims come from the orbital stability of Γ and the definition of ε. By orbital asymptotic stability and compactness there is r1 ∈ (0, δ0 /4] with Φf (t; x) ∈ Nδ0 /2 for all x ∈ N2r1 , t ∈ [0, T̂ + 2]. Since ∥ΦΘ (t; x) − Φf (t; x)∥ ≤ ε ≤ δ0 /4 (Assumption 4.4), also ΦΘ (t; x) ∈ N3δ0 /4 ⊂ Nδ0 . The flow closeness is the definition of ε. The length T̂ + 2 accommodates the time shifts of size at most 1 + 2ε/v arising in the proof of Theorem 4.13(iii). There a time argument in [−2, T̂ + 3] occurs, outside the interval on which ε is defined. No estimate is needed outside it. The point in question lies on the already constructed forward trajectory, and reindexing it at as many preceding or following crossings as needed exhibits it as ΦΘ (t′′ ; z) with z ∈ Σr2 ⊂ N2r1 and t′′ ∈ [0, T̂ + 2]: each step changes the time argument by a return time in [T̂ /2, T̂ + 1] (Lemma 4.8), less than the length of the target interval, so some index lands in it. Only the confinement claim is used there; the closeness bound is invoked only for t ∈ [0, T̂ + 1]. □ Proof of Lemma 4.8. The return maps come from the implicit function theorem, applied to a transversality function whose time-slope is bounded below uniformly on the section. We construct the true map first and treat the learned one as its perturbation. The true map. Let n := f (p)/v and hf (t, x) := n⊤ (Φf (t; x) − p). At (T̂ , p) we have hf = 0 and ∂t hf = n⊤ f (p) = v > 0, so the implicit function theorem yields C 1 maps τf and Pf (x) = Φf (τf (x); x) on a disc Σr2′ with τf (p) = T̂ , Pf (p) = p; since f ∈ C 2 the flow is C 2 and Dτf , DPf are Lipschitz. It remains to exclude earlier returns. Fix ς > 0 so small that ∥γ̂(t) − p∥ ≤ r for t ∈ [0, ς] ∪ [T̂ − ς, T̂ ], where r > 0 is chosen with n⊤ f (y) ≥ v/2 whenever ∥y − p∥ ≤ 2r; this is possible by continuity of f and f (p) ̸= 0. On the two end intervals hf (·, x) is then strictly increasing for x close to p, so it has no zero there other than the one near T̂ . On the middle interval set d1 := min ∥γ̂(t) − p∥ > 0, t∈[ς,T̂ −ς]

positive because Γ is a simple closed curve and the minimum is over a compact set. Trajectories starting in Σr2′ spread by at most eLT̂ over one period, so ∥Φf (t; x) − γ̂(t)∥ ≤ eLT̂ r2′ there. Imposing r2′ (1 + eLT̂ ) < d1 /2 gives ∥Φf (t; x) − p∥ ≥ d1 − eLT̂ r2′ > r2′ on the middle interval, so the true trajectory does not meet the disc there. Its first return is therefore the crossing near time T̂ and no other. The learned map. Let x ∈ Σr2 with r2 ≤ min(r2′ /2, r1 ) and hΘ (t, x) := n⊤ (ΦΘ (t; x) − p). By Lemma 4.7, |hΘ (τf (x), x)| ≤ ε. For the time derivative, decompose through the base point p ∈ Γ:
∂t hΘ = n⊤ fΘ ΦΘ (t; x) ≥ n⊤ f (p) − ∥fΘ (p) − f (p)∥ − ∥fΘ (ΦΘ (t; x)) − fΘ (p)∥ . | {z } | {z } | {z } =v

√

≤δ

≤ B ∥ΦΘ (t;x)−p∥

Lemma 4.5 gives δ ≤ 2C2 ε, so after shrinking ε2 we have δ ≤ vmin /(2(T̂ + 1)) ≤ v/(2(T̂ + 1)), a fixed margin below v/2. For (t, x) with |t − τf (x)|, ∥x − p∥ and ε small, ∥ΦΘ (t; x) − p∥ ≤

LONG-TIME APPROXIMATION VIA SA-NODES

35

ε + ∥Φf (t; x) − p∥ is as small as desired by continuity of the true flow; after one shrinking of r2 , ε2 and of the time interval, B ∥ΦΘ (t; x) − p∥ ≤ v/2 − δ, hence ∂t hΘ ≥ v/2. A scalar function with derivative at least v/2 and value of modulus at most ε at τf (x) has a unique nearby zero τΘ (x), with |τΘ (x) − τf (x)| ≤ 2ε/v. No earlier crossing of the disc occurs. On the middle interval the fΘ -trajectory stays within ε of the f -trajectory, hence at distance at least d1 /2 − ε ≥ d1 /4 > r2′ from p once ε2 ≤ d1 /4; note r2′ < d1 /4 from the choice r2′ (1 + eLT̂ ) < d1 /2. On the two end intervals the trajectory remains in the 2r-ball around p; after shrinking r so that in addition B(2r + ε) ≤ v/2 − δ, the slope bound ∂t hΘ ≥ v/2 holds there, so hΘ is strictly increasing on both end intervals. In particular the zero at t = 0, where the trajectory starts on the section, is left immediately, and the only zero with a crossing of the disc is the one near τf (x). Smoothness of τΘ , PΘ and the Lipschitz bound on DPΘ follow from the implicit function theorem and the C 2 bounds on ΦΘ from the variational equations, which depend only on (B, T̂ ). Finally, ∥PΘ (x) − Pf (x)∥ ≤ ∥ΦΘ (τΘ ; x) − ΦΘ (τf ; x)∥ + ∥ΦΘ (τf ; x) − Φf (τf ; x)∥ ≤ B · v2 ε + ε = K1 ε, and τf , τΘ ∈ [T̂ /2, T̂ + 1] after one more shrinking.

□

Proof of Lemma 4.9. The proof has two parts. An adapted norm turns the spectral bound of Assumption 4.6 into a genuine contraction, uniformly over the matrices that can arise; Banach’s fixed point theorem then produces the orbit and its bounds. Adapted norms. For a matrix A with ρ(A) ≤ ρ∗ and ς > 0 there is a norm with ∥A∥∗ ≤ ρ∗ + ς ([8, Proposition 2.84]; a Jordan basis with scaled nilpotent part). Apply this to A = DPΘ (p) with ς = (ρ̄ − ρ∗ )/2; by the Lipschitz continuity of DPΘ (constant Λ), for ∥x − p∥ ≤ r3 := min(r2 , c1 (ρ̄ − ρ∗ )/(2c2 Λ)) we get ∥DPΘ (x)∥∗ ≤ ρ̄, and the mean value inequality makes PΘ a ρ̄contraction on Σr3 in the adapted norm. The same construction for f (whose DPf (p) has spectral radius ρT (f ) < 1 by Assumption 4.1) gives ∥·∥∗f , c′1 , c′2 , r3′ . For d = 2 the section is one-dimensional and one may take the absolute value, with c1 = c2 = 1. For d ≥ 3 the equivalence constants have to be uniform over the family  M := A ∈ R(d−1)×(d−1) : ρ(A) ≤ ρ∗ , ∥A∥ ≤ CM , CM := C(f, B, T̂ ), which is compact because the spectral radius is continuous. We claim that one exponent serves the whole family. Write θ := ρ∗ + ς. For each A0 ∈ M, Gelfand’s formula gives an exponent j0 1/j0

1/j

1/j

0 0 with Aj00 ≤ ρ∗ + ς/2; the map A 7→ Aj0 is continuous, so Aj0 ≤ θ on an open ball around A0 . Finitely many such balls cover M; let j1 , . . . , jm be the associated exponents and put j∗ := lcm(j1 , . . . , jm ). A finite subcover alone would only give finitely many exponents; the passage to a common one uses submultiplicativity. Indeed, if A ∈ M lies in the i-th ball and j∗ = q ji , then q Aj∗ ≤ Aji ≤ θ ji q = θ j∗ .

For such an A define the adapted norm ∥x∥∗ :=

jX ∗ −1

θ−i Ai x .

i=0

Summing the same series shifted by one index and using Aj∗ x ≤ θ j∗ ∥x∥ gives ∥Ax∥∗ ≤ θ ∥x∥∗ , i gives ∥x∥ ≤ c ∥x∥ with that is ∥A∥∗ ≤ θ. The i = 0 term gives ∥x∥ ≤ ∥x∥∗ , and Ai ≤ CM 2 ∗ Pj∗ −1 i c2 := i=0 (CM /θ) . Both c1 = 1 and c2 depend only on (j∗ , CM , ρ∗ , ς), hence are uniform over M. Fixed point. By Lemma 4.8, ∥PΘ (p) − p∥∗ ≤ c2 K1 ε. On the ball of radius R := c2 K1 ε/(1 − ρ̄) around p, ∥PΘ (x) − p∥∗ ≤ ρ̄ ∥x − p∥∗ + ∥PΘ (p) − p∥∗ ≤ ρ̄R + c2 K1 ε = R,

36

Z. LI AND N. M. MATZAKOS

√ provided R ≤ c1 r3 ; together with the requirement 2C2 ε ≤ min(δ1 , r3 /(2CU )), where δ1 and CU are the certificate-free constants of Lemma 4.10(i), used in the stability step below and in Lemma 4.10(ii), this defines ε3 . Banach’s fixed point theorem gives the unique x∗Θ with ∥x∗Θ − p∥ ≤ c−1 1 R = K2 ε and geometric convergence of the iterates. For orbital exponential stability, apply Lemma 4.10(i), which uses neither Assumption 4.6 nor the present lemma: it sends every point of a neighborhood of γΘ to Σr3 within time T̂ + 1. The section distance to x∗Θ is then contracted by ρ̄ at each return, while deviations between consecutive crossings are amplified by at most eB(T̂ +1) . Together these convert the geometric decay on the section into exponential decay of the distance to γΘ . For the Hausdorff bound, for t ∈ [0, T̂ + 1], ∥ΦΘ (t; x∗Θ ) − γ̂(t)∥ ≤ ∥ΦΘ (t; x∗Θ ) − Φf (t; x∗Θ )∥ + ∥Φf (t; x∗Θ ) − Φf (t; p)∥ ≤ ε + eL(T̂ +1) K2 ε =: K3 ε, and every point of γΘ and of Γ is of this form. For the period, TΘ − T̂ ≤ |τΘ (x∗Θ ) − τΘ (p)| + |τΘ (p) − τf (p)| ≤ Lip(τΘ )K2 ε + v2 ε =: K4 ε. If Ltraj (Θ) = 0, then fΘ = f on Γ, so γ̂ solves ẋ = fΘ (x) and PΘ (p) = p; uniqueness gives x∗Θ = p, γΘ = Γ, TΘ = T̂ . □ Proof of Lemma 4.10. Write g for the field of part (i); the case of interest is g = fΘ . The strategy is to control the phase and the transverse coordinate of a g-trajectory simultaneously. The two estimates are coupled, because the bound on the phase speed needs the trajectory to stay near Γ while the transverse bound needs the phase speed to be of order one, so we derive both on a maximal interval and close the loop by continuation. Parametrize the tube by (s, u) 7→ γ̂(s) + U (s)u, with U (s) an orthonormal frame of the normal bundle along γ̂. For a small enough radius this is a diffeomorphism onto its image, with distortion constants depending only on (f, γ̂, δ0 ). Differentiating x(t) = γ̂(s(t)) + U (s(t))u(t) along a gtrajectory in the tube and projecting on the tangent and the normal directions gives ˙
  ⟨γ̂(s), g(x, t)⟩ ˙ 1 + O(|u|) , u̇ = U (s)⊤ g(x, t) − γ̂(s) ṡ − U ′ (s)u ṡ . (A.1) ṡ = 2 ˙ γ̂(s) Step 1: the continuation interval. Fix ϱ ≤ min(δ0 , r2 ) small enough that the phase-speed estimate (A.2) below holds on Nϱ and that ϱ < d1 /2, with d1 the middle-interval distance from the proof of Lemma 4.8; shrinking ϱ only strengthens the estimates below, and require r5′ ≤ ϱ/2 from the outset. Set  t∗ := sup t ≥ 0 : |u(τ )| ≤ ϱ for all τ ∈ [0, t] , which is positive because |u(0)| ≤ r0 ≤ r5′ ≤ δ0 /4. All estimates of Steps 2 and 3 are derived on [0, t∗ ), where |u| ≤ ϱ is available, and Step 4 shows t∗ > T̂ + 1. ˙ Step 2: the phase advances. At u = 0 one has γ̂(s) = f (γ̂(s)), so the quotient in (A.1) equals 2 ⟨f, g⟩/ ∥f ∥ evaluated at γ̂(s). For the true field this is exactly one. For g it is 1 + O(δ/vmin ), because ∥g − f ∥ ≤ δ on Γ and ∥f ∥ ≥ vmin there; the two agree only when g = f on Γ. Off the orbit we use ∥g(x, t) − f (γ̂(s))∥ ≤ δ + (B+2L) |u|, which combines ∥g − f ∥ ≤ δ + (B+L) dist(x, Γ) with the L-Lipschitz continuity of f , and obtain on [0, t∗ ) ṡ ≥ 1 −

T̂ δ + (B+2L)ϱ − Cϱ ≥ , vmin T̂ + 1

ṡ ≤ 2,

(A.2)

after shrinking ϱ, and with it r5′ and δ1 , once; the constant C absorbs the O(|u|) distortion of the chart. In particular ṡ > 0, so the phase increases strictly and the trajectory cannot stall. Step 3: the transverse coordinate. Positivity of ṡ constrains the phase only, and by itself does not prevent the trajectory from leaving the tube through its lateral boundary. The transverse coordinate must be estimated separately, and its equation carries the orbit mismatch as a forcing term. Insert g = f + (g − f ) and f (x) = f (γ̂(s)) + Df (γ̂(s))U (s)u + O(|u|2 ) into the second

LONG-TIME APPROXIMATION VIA SA-NODES

37

˙ ˙ equation of (A.1). The term proportional to f (γ̂(s)) = γ̂(s) is tangent, so U (s)⊤ γ̂(s) = 0 and that whole contribution, including the part carrying 1 − ṡ, drops out under the projection. The remaining terms are bounded on [0, t∗ ), using dist(x, Γ) ≤ |u| ≤ ϱ and the bound ṡ ≤ 2 of (A.2), by |u̇| ≤ C1 |u| + δ,

C1 := L + Cϱ + (B+L) + 2 sup U ′ (s) .

(A.3)

s

Grönwall’s inequality applied to (A.3) gives, for t ∈ [0, t∗ ) with t ≤ T̂ + 1,
|u(t)| ≤ eC1 (T̂ +1) r0 + δ(T̂ + 1) ≤ CU (r0 + δ), CU := (T̂ + 1) eC1 (T̂ +1) .

(A.4)

The additive term δ is the forcing produced by g − f and cannot be removed. For the true field δ = 0 and the estimate reduces to |u(t)| ≤ eC1 (T̂ +1) r0 , the form used in Lemma A.1. Step 4: closing the continuation. Shrink r5′ and δ1 once more so that CU (r5′ + δ1 ) < ϱ. Then (A.4) gives |u(t)| < ϱ strictly, for every t < min(t∗ , T̂ + 1). Were t∗ ≤ T̂ + 1, continuity would force |u(t∗ )| = ϱ, contradicting that strict inequality; hence t∗ > T̂ + 1 and both estimates hold on all of [0, T̂ + 1]. In particular the trajectory stays in Nϱ ⊂ Nδ0 and never reaches the lateral boundary, on [0, T̂ + 1] and not beyond. Step 5: arrival at the section. By (A.2) the phase increases at rate at least T̂ /(T̂ + 1), so s traverses a full period, and the trajectory therefore meets Σ, within time T̂ + 1. The crossing is transversal because ṡ > 0 there. At the crossing s ≡ 0 modulo T̂ , so the point is p + U (0)u ˙ with U (0)u ⊥ γ̂(0) = f (p); it therefore lies in the hyperplane Σ, and its distance to p equals |u|, bounded by CU (r0 + δ) through (A.4). This proves part (i). Step 6: uniqueness. Assume now Assumption 4.6 and CU (r0 + δ) ≤ r3 , and let γ ⊂ Nr0 be a periodic orbit of fΘ . By part (i) it crosses Σ at some x̄ with ∥x̄ − p∥ ≤ CU (r0 + δ) ≤ r3 , so x̄ ∈ Σr3 , and periodicity gives PΘm (x̄) = x̄ for some m ≥ 1. Every iterate PΘj (x̄) is again a crossing of the same orbit γ ⊂ Nr0 , so part (i), applied from the preceding crossing, keeps it in Σr3 as well. Each step then contracts in the adapted norm of Lemma 4.9 along the segment to x∗Θ , which stays in the convex disc Σr3 , so ∥x̄ − x∗Θ ∥∗ = ∥PΘm (x̄) − PΘm (x∗Θ )∥∗ ≤ ρ̄ m ∥x̄ − x∗Θ ∥∗ , which forces x̄ = x∗Θ and γ = γΘ . □ Proof of Lemma 4.11. The two orbits are compared through the contraction of the true return map, with the difference of the maps as a forcing term. We first record that the iterates stay where the contraction is available. Write r7 := min(r3 , r3′ ), so that Σr7 = Σr3 ∩ Σr3′ , and recall ∥x∗Θ − p∥ ≤ K2 ε from Lemma 4.9. Both balls below are calibrated to Σr7 , and the one centred at x∗Θ carries the offset K2 ε so that it lies in Σr7 and not merely within r7 of x∗Θ . In the adapted norms of Lemma 4.9 the balls B∗ := {y : ∥y − x∗Θ ∥∗ ≤ c1 (r7 − K2 ε)},

B∗f := {y : ∥y − p∥∗f ≤ c′1 r7 }

both lie in Σr7 , and each is invariant under its own map. For the inclusions, y ∈ B∗f gives ∥y − p∥ ≤ (c′1 )−1 ∥y − p∥∗f ≤ r7 , while y ∈ B∗ gives ∥y − p∥ ≤ ∥y − x∗Θ ∥ + ∥x∗Θ − p∥ ≤ (r7 − K2 ε) + K2 ε = r7 , the offset being what the radius of B∗ was reduced by. Each ball is convex and centred at the fixed point of its own map, and Σr7 ⊆ Σr3 ∩ Σr3′ is where the two contraction bounds of Lemma 4.9 hold, so the mean value inequality along the segment to the centre gives invariance. Choose r4 ≤ 21 min{c1 /c2 , c′1 /c′2 } r7 ,

r4 ≤ r5 ,

and ε4 so small that K2 ε ≤ r67 min{c1 /c2 , c′1 /c′2 }. The disc Σr4 , enlarged by K2 ε, then lies in B∗ ∩ B∗f : this inclusion runs from the Euclidean norm into the adapted ones and so costs c2 and c′2 , requiring c′2 (r4 +K2 ε) ≤ c′1 r7 and c2 (r4 +2K2 ε) ≤ c1 (r7 −K2 ε). The first holds since c′2 r4 ≤ 21 c′1 r7 and c′2 K2 ε ≤ 16 c′1 r7 ; the second since c2 r4 ≤ 21 c1 r7 , 2c2 K2 ε ≤ 31 c1 r7 , and c1 (r7 − K2 ε) ≥ 65 c1 r7 .

38

Z. LI AND N. M. MATZAKOS

Then every iterate of either sequence remains in Σr7 = Σr3 ∩ Σr3′ , which is what the estimate below uses. Set d˜k := P k (yf ) − P k (yΘ ) . Then f

Θ

∗f

d˜k+1 ≤ Pf (Pfk yf ) − Pf (PΘk yΘ )

∗f

+ (Pf − PΘ )(PΘk yΘ )

∗f

≤ ρ̄f d˜k + c′2 K1 ε,

2K1 by Lemma 4.9 and Lemma 4.8. With d˜0 ≤ c′2 K1 ε, induction gives d˜k ≤ c′2 1−ρ̄ ε for all k, and f norm equivalence concludes. □

Lemma A.1 (hitting the section). There are r6 ∈ (0, r5 /2] and ε5 ∈ (0, ε4 ] such that for ε ≤ ε5 and every x0 ∈ Nr6 : the f - and fΘ -trajectories from x0 first cross Σr4 transversally at times σf , σΘ ≤ T̂ + 1 and points yf , yΘ ∈ Σr4 , with |σΘ − σf | ≤ v2 ε and ∥yΘ − yf ∥ ≤ K1 ε. Proof. The estimate is Lemma 4.10(i) for the true field, sharpened so that the crossing lands in Σr4 , together with a simple-zero argument for the learned crossing. Lemma 4.10(i) applied to f itself, the case δ = ε = 0, yields the first f -crossing within time T̂ + 1, transversal, at yf with ∥yf − p∥ ≤ CU r6 ≤ r4 /2 after shrinking r6 ; then ∥yΘ − p∥ ≤ r4 /2 + K1 ε ≤ r4 for small ε. The crossing is a simple zero with slope at least v/2 of hf (t) := n⊤ (Φf (t; x0 ) − p); by Lemma 4.7, |hΘ (σf )| ≤ ε and ∂t hΘ ≥ v/2 near σf as in the proof of Lemma 4.8, so hΘ has a unique nearby zero σΘ with |σΘ − σf | ≤ 2ε/v. That this is the first fΘ -crossing of the disc uses more than monotonicity of the phase, since the trajectory also meets the hyperplane at phases away from 0. Since ϱ < d1 /2 by construction in the proof of Lemma 4.10, shrink r4 so that r4 < d1 /2, whence r4 + ϱ < d1 , with d1 the middle-interval distance of the proof of Lemma 4.8 and ϱ the tube radius of Lemma 4.10. Any point of the tube within r4 of p then has phase within ς of 0 modulo T̂ , so the hyperplane crossings at intermediate phases fall outside Σr4 , and monotonicity of the phase gives the first crossing. The point estimate follows as in Lemma 4.8. □ Proof of Lemma 4.12. Each return adds a clock difference of order ε, and summing gives the bound. Each increment satisfies τΘ (PΘk yΘ ) − τf (Pfk yf ) ≤ (τΘ − τf )(PΘk yΘ ) + τf (PΘk yΘ ) − τf (Pfk yf ) ≤ v2 ε + Lip(τf )K5 ε, by Lemmas 4.8 and 4.11; summing over k < n gives the claim.

□

Lemma A.2 (scalar multiplier in dimension two). Let d = 2 and let Assumptions 4.1–4.4 hold with η = 0 and ε ≤ ε2 . If x̄ ∈ Σr2 is a fixed point of the return map PΘ , with orbit γ̄ := {ΦΘ (t; x̄) : 0 ≤ t ≤ τ̄ } of period τ̄ := τΘ (x̄), then Z τ̄  div fΘ (γ̄(t)) dt > 0. DPΘ (x̄) = exp 0

No contraction hypothesis is used; the identity holds at any periodic orbit of the autonomous planar field fΘ . Proof. The identity is Liouville’s formula for the variational equation along the periodic orbit, read on the section. By Lemma 4.8 the return map PΘ is C 1 at x̄ and fΘ (x̄) ̸= 0, so γ̄ is a genuine periodic orbit of the autonomous field fΘ . Let X(t) solve the variational equation Ẋ = DfΘ (γ̄(t))X, X(0) = I, and M := X(τ̄ ). Since fΘ is autonomous, t 7→ fΘ (γ̄(t)) solves the variational equation and is τ̄ -periodic, so M fΘ (x̄) = fΘ (x̄): the tangent direction is a Floquet direction with multiplier one. In the basis {fΘ (x̄), w} with w spanning Tx̄ Σ, the matrix M is block triangular with
and m = DPΘ (x̄) [8, Proposition 2.80]. Liouville’s formula gives R τ̄ diagonal (1, m) det M = exp 0 tr DfΘ (γ̄) dt , and det M = 1 · m; positivity follows from the exponential. □

LONG-TIME APPROXIMATION VIA SA-NODES

39

A.4. Proof of Theorem 4.13. We now prove the linear long-horizon bound, assembling the lemmas of Appendix A.3. Proof of Theorem 4.13. The proof assembles the lemmas. Transverse error is contracted at every return (Lemmas 4.9 and 4.11), the two clocks drift apart by O(ε) per return (Lemma 4.12), and Grönwall estimates are applied over single periods only, through Lemma 4.7. We reduce to the section, then prove the three parts in order. Reduction to the section. Set r0 := r6 , with r6 shrunk once more if necessary so that Nr6 ⊆ N2r1 , and ε0 := min{ε5 , r0 /K3 , v/2}, so that the hypothesis K3 ε ≤ r0 of Lemma 4.10(ii) and the bound 2ε/v ≤ 1 used in part (iii) hold for every ε ≤ ε0 ; the δ-smallness conditions of Lemmas 4.9 and 4.10(ii) are folded into ε3 , and every statement used below holds a fortiori for smaller radii and thresholds. For t ≤ T̂ + 1, parts (ii)–(iii) follow from Lemma 4.7, part (ii) together with the comparison dist(Φf (t; x0 ), Γ) ≤ eL(T̂ +1) dist(x0 , Γ), absorbed into C. For t > T̂ + 1, Lemma A.1 provides crossing data σf , σΘ ≤ T̂ + 1, |σΘ − σf | ≤ 2ε/v, and yf , yΘ ∈ Σr4 with ∥yf − yΘ ∥ ≤ K1 ε, such that Φf (t; x0 ) = Φf (t − σf ; yf ) and ΦΘ (t; x0 ) = ΦΘ (t − σΘ ; yΘ ). Write t′ := t − σf . Below we evaluate both flows at the common time t′ ; for the learned flow this replaces t − σΘ by t′ , which costs ∥ΦΘ (t − σΘ ; yΘ ) − ΦΘ (t′ ; yΘ )∥ ≤ B |σΘ − σf | ≤ 2Bε/v, an O(ε) term absorbed into the constants of parts (ii)–(iii) without further mention. Part (i). This is Lemmas 4.9 and 4.10. Part (ii). Let n be the number of fΘ -crossings in (0, t′ ], so t′ = t̃Θ n + s with s ∈ [0, T̂ + 1), where Θ t̃n are the fΘ -crossing times from yΘ . Then ΦΘ (t′ ; yΘ ) = ΦΘ (s; PΘn yΘ ) and
dist ΦΘ (t′ ; yΘ ), γΘ ≤ ∥ΦΘ (s; PΘn yΘ ) − ΦΘ (s; x∗Θ )∥ ≤ eB(T̂ +1) ∥PΘn yΘ − x∗Θ ∥ ≤ C ρ̄ n ∥yΘ − x∗Θ ∥ . Since return times are at most T̂ + 1, n ≥ t′ /(T̂ + 1) − 1; together with distH (γΘ , Γ) ≤ K3 ε and ∥yΘ − x∗Θ ∥ ≤ C(dist(x0 , Γ) + ε) this proves (ii); the shifts from t′ to t and from σf to σΘ are absorbed in C. Part (iii). Let n be the number of f -crossings in (0, t′ ]: t′ = tfn + s, s ∈ [0, T̂ + 1), and n Φf (t′ ; yf ) = Φf (s; Pfn yf ). Set δn := (tfn − tΘ n ) + (σf − σΘ ), so that ΦΘ (t; x0 ) = ΦΘ (s + δn ; PΘ yΘ ); by Lemma 4.12, |δn | ≤ K6 nε + 2ε/v. Assume first K6 nε ≤ 1, so that |δn | ≤ 1 + 2ε/v ≤ 2. Then s+δn ∈ [−2, T̂ +3]. The point ΦΘ (s+δn ; PΘn yΘ ) lies on the already constructed forward trajectory, and re-indexed at finitely many preceding or following crossings, until the shifted time first lands in [0, T̂ + 2], which is possible because consecutive crossings are at most T̂ + 1 apart, it is ΦΘ (t′′ ; z) with z ∈ Σr2 and t′′ ∈ [0, T̂ + 2], so the confinement of Lemma 4.7 applies to it and keeps all points below in Nδ0 . The closeness bound of that lemma is used only at the time s ∈ [0, T̂ + 1); hence ∥ΦΘ (t; x0 ) − Φf (t; x0 )∥ ≤ ∥ΦΘ (s + δn ; PΘn yΘ ) − ΦΘ (s; PΘn yΘ )∥ + ∥ΦΘ (s; PΘn yΘ ) − Φf (s; PΘn yΘ )∥ + Φf (s; PΘn yΘ ) − Φf (s; Pfn yf )  t ≤ B |δn | + ε + eL(T̂ +1) K5 ε ≤ Cε 1 + , T̂ using Lemma 4.7 for the middle term (the first term needs no confinement, since it uses only the global bound ∥fΘ ∥ ≤ B supplied by the cutoff extension), Lemma 4.11 for the last, and n ≤ 2t/T̂ +1 from the lower bound T̂ /2 on return times; one may take C = 1 + eL(T̂ +1) K5 + 2B(K6 + 1/v). If K6 nε > 1, both trajectories lie in Nδ0 by confinement and (ii), so ∥e(t)∥ ≤ diam(Nδ0 ) ≤ diam(Nδ0 )K6 nε ≤ Cε(1 + t/T̂ ) upon enlarging C. □ A.5. Proofs of Proposition 4.16 and Corollary 4.18. This part proves the two certification results of Section 4.3: the C 1 route in any dimension and the Floquet-loss route in the autonomous plane.

40

Z. LI AND N. M. MATZAKOS

Proof of Proposition 4.16. We prove the monodromy comparison (ii) first, then the orbit persistence (i), and then the spectral consequence (iii); (ii) comes first because the comparison along the reference orbit is what the other two parts perturb. Part (ii). Let XΘ (t), X(t) be the fundamental matrices of the variational equations of fΘ , f along γ̂, with XΘ (0) = X(0) = I. Writing the difference by variation of constants and using ∥fΘ − f ∥C 1 (U ) ≤ ε1 and the C 1 -bound of f on U gives ∥XΘ (t) − X(t)∥ ≤ T̂ eC T̂ ε1 for t ∈ [0, T̂ ]. Restricting to the section Σ ⊥ f (p) and projecting out the tangent direction f (p), which is the orthogonal projection onto Σ = f (p)⊥ and so has norm one, bounds the transverse mon⊥ − M ⊥ ≤ Cε , where M ⊥ is the projected fundamental matrix of the odromies along γ̂, MΘ,γ̂ 1 f Θ,γ̂ ⊥ fΘ -variational equation along γ̂; the transition to the monodromy MΘ⊥ := MΘ,Γ along the learned Θ orbit itself, which is the object in (ii), is made at the end of part (i). Part (i). The true return map satisfies Pf (p) = p with DPf (p) = Mf⊥ having no eigenvalue 1 (hyperbolicity), so I − DPf (p) is invertible. Since fΘ is C 1 -close to f , its return map PΘ is a Cε1 perturbation of Pf in C 1 . The implicit function theorem applied to x 7→ PΘ (x) − x then yields a unique fixed point x∗Θ with ∥x∗Θ − p∥ ≤ Cε1 . Its orbit ΓΘ satisfies distH (ΓΘ , Γ) ≤ Cε1 . It remains to compare the monodromy of fΘ along its own orbit ΓΘ with the one along γ̂. Boundedness of DfΘ does not suffice, because the two variational equations are posed along different curves and over different periods. What is needed is the Lipschitz continuity of the coefficient, and this follows from f ∈ C 2 alone: writing Λf := Lip(Df ) on the compact tube, a quantity determined by f ,

∥DfΘ (y) − DfΘ (z)∥ ≤ 2ε1 + Λf ∥y − z∥

for y, z ∈ U,

since ∥DfΘ − Df ∥ ≤ ε1 at both points. No bound on ∥fΘ ∥C 2 is used, so the constants below still depend only on (f, Γ). Parametrize ΓΘ by γΘ (t) = ΦΘ (t; x∗Θ ). On U the learned field is (L + ε1 )-Lipschitz, so ∥γΘ (t) − γ̂(t)∥ ≤ ∥ΦΘ (t; x∗Θ ) − ΦΘ (t; p)∥ + ∥ΦΘ (t; p) − Φf (t; p)∥ ≤ e(L+ε1 )t ∥x∗Θ − p∥ + t e(L+ε1 )t ε1 ≤ Cε1 for t ∈ [0, T̂ + 1]. Let Y and XΘ solve the variational equations of fΘ along γΘ and along γ̂, both from the identity. By the display above their coefficients differ by at most 2ε1 + Λf Cε1 , so variation of constants and Grönwall give ∥Y (t) − XΘ (t)∥ ≤ Cε1 on [0, T̂ + 1]. The two periods differ by TΘ − T̂ ≤ Cε1 , and Ẏ ≤ (L + ε1 )e(L+ε1 )(T̂ +1) , so evaluating at TΘ rather than T̂ costs a further Cε1 . One step remains, because the transverse monodromies are read off through different projections: the return map at x∗Θ projects along fΘ (x∗Θ ) and the one at p along f (p), onto the same hyperplane Σ. The two projectors differ by at most C ∥fΘ (x∗Θ ) − f (p)∥ ≤ Cε1 . ⊥ Combining the four estimates with (ii) yields MΘ,Γ − Mf⊥ ≤ Cε1 . Uniqueness of ΓΘ among Θ the periodic orbits of a tube around Γ is deferred, because it is not a consequence of the uniqueness of the fixed point on Σr : an orbit might cross the section several times, or outside that disc. It follows from Lemma 4.10(ii), whose hypotheses are δ ≤ δ1 , satisfied here since δ ≤ ε1 < ε∗ , and Assumption 4.6. The latter is not available yet, so we return to it after part (iii). Part (iii). The spectral radius is continuous at Mf⊥ , so with ρ̄f = (1 + ρT (f ))/2, the constant
⊥ of Lemma 4.9, there is ε∗ > 0, depending only on (f, Γ), such that ρΘ = ρ MΘ,Γ ≤ ρ̄f < 1 for Θ ε1 < ε∗ . This is a bound at the fixed point of the learned return map, whereas (4.5) is posed at the base point p, so one further step is needed. By part (i) the two points lie within Cε1 of each other and DPΘ is Lipschitz on the disc, so ρ(DPΘ (p)) ≤ ρ̄f + ω(C ′ ε1 ), with ω the modulus of continuity of the spectral radius on the norm ball {A : ∥A∥ ≤ CM } of the proof of Lemma 4.9. As ρ̄f < 1, shrinking ε∗ once more makes the right side at most (1 + ρ̄f )/2 < 1, so Assumption 4.6 holds with ρ∗ := (1 + ρ̄f )/2. Lemma 4.10(ii) now applies and completes the uniqueness claim

LONG-TIME APPROXIMATION VIA SA-NODES

41

of part (i); the constant it contributes depends on ρ̄f , hence on ρT (f ), in addition to (f, Γ). In dimension two the transverse monodromy is a scalar and the rate is linear, ρΘ ≤ ρT (f ) + Cε1 with ε∗ = (1 − ρT (f ))/(2C). For d > 2 no linear rate holds in general: a defective Mf⊥ can move its 1/k

spectral radius by O(ε1 ) under an O(ε1 ) perturbation. When Mf⊥ = V ΛV −1 is diagonalizable, the Bauer–Fike theorem [1] restores the quantitative threshold: ρΘ ≤ ρT (f ) + κ(V ) Cε1 ,

κ(V ) := ∥V ∥ V −1 ,

with ε∗ = (1 − ρT (f ))/(2κ(V )C). This is the route used for certification in Section 4.3.

□

Proof of Corollary 4.18. Throughout d = 2 and η = 0, so fΘ is autonomous; only Assumptions 4.1– 4.4 are used, never the certificate being derived. Step 1: the return map, without a certificate. By Lemma 4.8, PΘ : Σr2 → Σ is well defined and C 1 with ∥PΘ − Pf ∥C 0 (Σr ) ≤ K1 ε; this step uses only Assumptions 4.1 and 4.4. The construction of the 2 adapted norm for the true map in the proof of Lemma 4.9 uses only Assumption 4.1 and provides ∥·∥∗f with equivalence constants c′1 , c′2 and a radius r3′ on which ∥DPf ∥∗f ≤ ρ̄f := (1+ρT (f ))/2 < 1. Step 2: existence and location of the learned orbit. Since Pf (p) = p, the closed ∥·∥∗f -ball B̄ := {y : ∥y − p∥∗f ≤ R} with R := c′2 K1 ε/(1 − ρ̄f ) lies inside Σr3′ once ε is small, and there Pf (B̄) ⊆ { ∥· − p∥∗f ≤ ρ̄f R } by the mean value inequality; ∥PΘ − Pf ∥∗f ≤ c′2 K1 ε = (1 − ρ̄f )R then gives PΘ (B̄) ⊆ B̄. In dimension two Σ is a line and B̄ a compact interval, so the continuous self-map PΘ has a fixed point in B̄ by the intermediate value theorem. Moreover every fixed point x̄ ∈ Σr3′ obeys ∥x̄ − p∥∗f = ∥PΘ (x̄) − Pf (p)∥∗f ≤ c′2 K1 ε + ρ̄f ∥x̄ − p∥∗f , the segment [p, x̄] staying ′ in Σr3′ , hence ∥x̄ − p∥ ≤ c′−1 1 R =: K2 ε. The relaxed orbit γ̂Θ of Definition 4.17, once it closes into a periodic orbit inside the tube Nr5′ (the implementation tests closure; membership of the tube is assumed, as discussed after (4.8)), crosses Σ within distance CU (r5′ + δ) of p by Lemma 4.10(i), which uses neither Assumption 4.6 nor an adapted radius. Both terms are at our disposal: r5′ may √ be shrunk, and δ ≤ 2C2 ε by Lemma 4.5, so CU (r5′ + δ) ≤ r3′ once r5′ and ε are small. A closed orbit may a priori meet the section several times before closing, so that crossing is only a periodic point, PΘm (x̄) = x̄ for some m ≥ 1. In dimension two we can rule out m ≥ 2 without a certificate. By Step 2 of the proof of Lemma 4.10 the phase increases strictly along every trajectory in the tube, so γ̂Θ winds monotonically around the annulus Nr5′ . Being a periodic orbit of the autonomous field fΘ , it is a simple closed curve, by uniqueness of solutions; a simple closed curve in the annulus has winding number zero or ±1 about Γ, and strict phase monotonicity excludes zero. Its winding number is therefore one, so it meets the transverse section exactly once per period, whence m = 1. The crossing is therefore a fixed point of PΘ in Σr3′ , and within K2′ ε of p by the a-priori bound. We take x∗Θ to be this crossing and γΘ := {ΦΘ (t; x∗Θ )} = γ̂Θ its orbit, of period TΘ := τΘ (x∗Θ ) with TΘ − T̂ ≤ K4′ ε, where K4′ := Lip(τf )K2′ + 2/v is built from the certificate-free constant K2′ of this step and Lemma 4.8. When the relaxation degenerates instead, the implementation declares the surrogate unevaluable and no certificate is claimed; the intermediate-value fixed point above still guarantees that the periodic orbit exists. Step 3: the scalar multiplier. Lemma A.2 applied to x∗Θ gives Z T Θ  ∗ DPΘ (xΘ ) = exp div fΘ (γΘ ) > 0, 0

which in dimension two equals ρ(DPΘ (x∗Θ )).
R T̂ Step 4: surrogate versus multiplier. The surrogate is ρ̃T (Θ) = exp 0 divx fΘ (γ̂Θ ) , the same integrand along the same orbit γ̂Θ = γΘ , re-parametrized to start at x∗Θ , a phase shift that changes neither integral by periodicity of the integrand; the two differ only in the integration length T̂ versus

42

Z. LI AND N. M. MATZAKOS

TΘ . With |div fΘ | ≤ dB on the tube, log ρ̃T (Θ) − log ρ(DPΘ (x∗Θ ))

Z T̂ =

div fΘ (γΘ ) ≤ dB K4′ ε,

TΘ

so, using |a − b| ≤ max(a, b) (eκ − 1) with κ := dBK4′ ε and max(a, b) ≤ edB(T̂ +1) , we get |ρ̃T (Θ) − ρ(DPΘ (x∗Θ ))| ≤ Cε with C := dB K4′ edB(T̂ +2) , valid once κ ≤ dB; the constant involves no certificate. Step 5: conclusion. If Ltraj (Θ) = 0 then ΦΘ (·; p) is T̂ -periodic, so p is a fixed point of PΘ and Γ is an orbit of fΘ . This alone does not identify p with x∗Θ : the map PΘ is not yet known to be a contraction, and a second cycle of fΘ could carry the relaxation elsewhere. The hypothesis γ̂Θ = Γ of (i) supplies the identification, and the implementation checks it by testing that the relaxed curve passes near p. Under it x∗Θ = p and TΘ = T̂ , so Step 4 collapses to equality and ρ(DPΘ (p)) = ρ̃T (Θ) ≤ ρ∗ , proving (i). In general, Step 4 and ρ̃T (Θ) ≤ ρ∗ give ρ(DPΘ (x∗Θ )) ≤ ρ∗ + Cε < 1 once ε ≤ (1 − ρ∗ )/(2C), which is Assumption 4.6 at x∗Θ ; since ∥x∗Θ − p∥ ≤ K2′ ε and DPΘ is Lipschitz (Lemma 4.8), the same holds at p with C ′ := C + Lip(DPΘ )K2′ , so ρ(DPΘ (p)) ≤ ρ∗ + C ′ ε < 1 once ε ≤ (1 − ρ∗ )/(2C ′ ), proving (ii). Theorem 4.13 then applies and the chain (1.6) holds. □ A.6. Proof of Lemma 4.15. This part proves the reduction of the periodic encoding to the autonomous case. The convention η = 0 of Appendix A.3 is dropped here. Proof of Lemma 4.15. The proof has two parts. Part (i) compares the encoded flow with the flow of the averaged field, in C 0 and in C 1 , by a Grönwall estimate whose forcing is the oscillation η. Part (ii) applies Theorem 4.13 to the averaged field and transfers its conclusions back to the encoded flow. Part (i): the averaged flow is CB η-close. Averaging preserves the C 2 bound, so f¯Θ C 2 (N ) ≤ B δ0

(s) and Df¯Θ is B-Lipschitz. Fix s ∈ R and x ∈ N2r1 , and set w(t) := ΦΘ (t; x) − Φf¯Θ (t; x). While both trajectories remain in Nδ0 ,    (s) (s)  (s) ẇ = fΘ (ΦΘ , s + t) − f¯Θ (ΦΘ ) + f¯Θ (ΦΘ ) − f¯Θ (Φf¯Θ ) ,

so ∥ẇ∥ ≤ η+B ∥w∥, and Grönwall with w(0) = 0 gives ∥w(t)∥ ≤ η t eBt on [0, T̂ +2]; confinement of both trajectories follows as in Lemma 4.7 once ε + CB η ≤ δ0 /4. For the derivative, the variational (s) solutions V (s) := Dx ΦΘ and V̄ := Dx Φf¯Θ satisfy V (s) , V̄ ≤ eBt and

(s) d − V̄ ≤ η + B ∥w∥ eBt + B V (s) − V̄ , dt V using Dx fΘ (·, s + t) − Df¯Θ ≤ η and the B-Lipschitz continuity of Df¯Θ ; Grönwall again gives V (s) (t) − V̄ (t) ≤ η t (1 + Bt) e2Bt . The sum of the two bounds is at most CB η on [0, T̂ + 2]. The closeness of f¯Θ to f follows from the case s = 0 and the triangle inequality; on Γ, f¯Θ − f ≤ δ + η since f¯Θ − fΘ (·, t) has norm at most η. (ii) Step 1: the averaged field. By (i), f¯Θ is an autonomous C 2 field satisfying Assumptions 4.4– 4.6 with closeness ε + CB η, orbit mismatch δ + η, and contraction constant ρ∗ ; after decreasing ε6 these place it in the scope of Theorem 4.13. Its conclusions hold for Φf¯Θ : a unique orbit γ̄Θ with the persistence bounds of part (i), and  t Φf¯Θ (t; x0 ) − Φf (t; x0 ) ≤ C(ε + CB η) 1 + . T̂ Step 2: the encoded field against the averaged field. Rerun the proof of Theorem 4.13(iii) with reference f¯Θ , whose orbit is hyperbolic and stable with ρ(DP̄ ) ≤ ρ∗ , and with the encoded flow as the perturbed system, one-period closeness CB η uniform in the launch phase by (i). The barred constants K̄1 and v̄ below are the analogues of K1 and v for the averaged field f¯Θ ; by part (i)

LONG-TIME APPROXIMATION VIA SA-NODES

43

√ and Lemma 4.5 they are finite and O( ε + η)-close to the originals, while ρ̄ = (1 + ρ∗ )/2 keeps its value. Four replacements are needed; no other step of that proof uses the autonomy of the (s) perturbed field. First, for each phase s the argument of Lemma 4.8 applied to the pair (f¯Θ , ΦΘ ) (s)

(s)

yields a return map PΘ with PΘ − P̄

C 0 (Σr2 )

≤ K̄1 CB η uniformly in s; the k-th crossing of an

(s )

encoded trajectory is zk+1 = PΘ k (zk ) with sk the accumulated crossing time modulo T̂ . Second, in the proof of Lemma 4.11 the recursion becomes d˜k+1 ≤

P̄ (P̄ k ȳ) − P̄ (zk )

(s )

∗

+ (P̄ − PΘ k )(zk )

∗

≤ ρ̄ d˜k + c̄2 K̄1 CB η,

which uses only the contraction of P̄ in its adapted norm, with equivalence constant c̄2 , and the uniform C 0 bound; the confinement of the crossings zk in Σr3 follows from d˜k ≤ K5′ η by the same induction. Third, the initial crossing data is supplied by Lemma A.1 for the pair, whose winding estimate for the encoded flow holds with the time-uniform mismatch δ of Assumption 4.4, the slope bound being pointwise in time. Fourth, in Lemma 4.12 each increment obeys (s ) τΘ k (zk ) − τf¯Θ (P̄ k ȳ) ≤ (2/v̄)CB η + Lip(τf¯Θ )d˜k , and the sums remain linear in n. All constants of this step depend only on (f, γ̂, δ0 , B, ρ∗√ , ρT (f )) for ε + η ≤ ε6 , since the data of f¯Θ (base point, speed, tube, Lipschitz constants) are O( ε + η)-perturbations of those of f : the speed through the orbit mismatch δ + η of Lemma 4.5, the base point through persistence at rate O(ε + η), and the tube and Lipschitz constants only through the uniform bound B. Only their smallness enters the constants. The assembly of Theorem 4.13(iii) then gives  t , ΦΘ (t; x0 ) − Φf¯Θ (t; x0 ) ≤ Cη 1 + T̂ and with Step 1 and the triangle inequality this proves (4.7). The attraction statement does not follow from the tracking bound, whose linear growth in t would spoil uniformity; it follows from the recursion displayed above. Iterating d˜k+1 ≤ ρ̄ d˜k + c̄2 K̄1 CB η gives d˜k ≤ ρ̄ k d˜0 + c̄2 K̄1 CB η/(1 − ρ̄), uniformly in k: this is the cap recorded above. Between consecutive crossings the trajectory is flowed for a time at most T̂ + 1 by a field of spatial Lipschitz constant at most B, so a single-period Grönwall factor eB(T̂ +1) , together with the norm equivalences, converts the crossing bound into

(0) dist ΦΘ (t; x0 ), γ̄Θ ≤ C ρ̄ t/(T̂ +1) dist(x0 , γ̄Θ ) + η for all t ≥ 0: the trajectory enters and remains in a Cη-neighborhood of γ̄Θ . The distance conclusion for the encoded flow follows by adding distH (γ̄Θ , Γ); the tracking bound and the triangle inequality enter only in the trajectory estimate (4.7), which is where the linear term originates. □ A.7. Proofs of Proposition 4.21 and Theorem 4.22. This part proves the local obstruction and the orbital guarantee. Both proofs stand outside the convention of Appendix A.3: they use the standing hypotheses of Section 4.4, namely Assumption 4.1, exact T̂ -periodicity of fΘ in t, and ∥fΘ ∥C 2 (Rd ×R) ≤ B, together with the constant ε of Assumption 4.4 where stated, and their constants depend only on (d, T̂ , B, ρ∗ ). Proof of Proposition 4.21. The proof has four steps. The tangent to Γ is an exact eigenvector of the target’s stroboscopic map, with eigenvalue one. An interpolation inequality upgrades the C 0 closeness of Assumption 4.4 to C 1 closeness, at the cost of a square root. The tangent is therefore an approximate eigenvector of M with approximate eigenvalue one, and Hadamard’s inequality converts that into a lower bound on the spectral radius. Step 1: the unit eigenvalue of the target. Since Γ has period exactly T̂ , Φf (T̂ ; γ̂(s)) = γ̂(s) for ˙ ˙ every s. Differentiating in s gives DSf (γ̂(s)) γ̂(s) = γ̂(s), where Sf := Φf (T̂ ; ·). Fix x ∈ Γ with ♯ x − x ≤ ϱ and let τ be the unit tangent to Γ at x, which is well defined because ∥f ∥ ≥ vmin > 0 on Γ. Then DSf (x)τ = τ .

44

Z. LI AND N. M. MATZAKOS

Step 2: from C 0 to C 1 closeness. Put g := SΘ − Sf on N2r1 . Assumption 4.4 at t = T̂ gives ∥g∥C 0 (N2r ) ≤ ε. Both maps are C 2 : for SΘ this uses only the global bound ∥fΘ ∥C 2 (Rd ×R) ≤ B 1 supplied by the cutoff extension of Assumption 4.4, with no confinement and no condition on ε; for Sf , the confinement of the true flow, built from f alone in the proof of Lemma 4.7 and likewise free of any condition on ε, keeps every f -trajectory from N2r1 in Nδ0 for t ∈ [0, T̂ ], where ∥f ∥C 2 (Nδ ) ≤ B; the first and second variational equations then bound D2 Sf and D2 SΘ on 0

N2r1 by a constant depending only on (d, T̂ , B). Write Λ2 ≥ 1 for the resulting bound on D2 g there. For unit vectors u, w set φ(h) := ⟨w, g(x + hu)⟩, defined for |h| ≤ r1 since x ∈ Γ. Taylor expansion at 0 in both directions and subtraction give, for 0 < a ≤ r1 , ε aΛ2 φ′ (0) ≤ + , a 2 p √ and the choice a = 2ε/Λ2 , admissible under the stated bound on ε, yields |φ′ (0)| ≤ 2εΛ2 . √ Taking the supremum over u and w gives ∥Dg(x)∥ ≤ 2εΛ2 . Step 3: the tangent is an approximate eigenvector of M . With Λ1 a Lipschitz constant for DSΘ , again from the variational equations, DSΘ (x♯ ) − DSΘ (x) + ∥DSΘ (x) − DSf (x)∥ + ∥DSf (x)τ − τ ∥ p ≤ Λ1 ϱ + 2εΛ2 =: ∆,

∥(M − I)τ ∥ ≤

the last term vanishing by Step 1. Step 4: Hadamard. Complete τ to an orthonormal basis of Rd . Hadamard’s inequality bounds the determinant by the product of the column norms in that basis, so |det(M − I)| ≤ ∥(M − I)τ ∥ ∥M − I∥ d−1 ≤ ∆ (1 + ∥M ∥)d−1 . Q Since det(M − I) = i (λi − 1) over the eigenvalues of M , some eigenvalue satisfies |λ − 1| ≤ [∆(1 + ∥M ∥)d−1 ]1/d , and ρ(M ) ≥ |λ| ≥ 1 − |λ − 1| gives (4.14). □ Proof of Theorem 4.22. The proof has three steps. Step 1 turns the spectral hypothesis (B1) into a genuine contraction of SΘ on a ball around its fixed point, in an adapted norm. Step 2 uses the periodicity identity (4.12) to pass from the iterates of SΘ to the flow at intermediate times. Step 3 is a triangle inequality against Γ through hypothesis (B2). No section, return time or clock comparison appears, and Assumption 4.4 is never used. Step 1: contraction of the stroboscopic map. Since ∥fΘ ∥C 2 ≤ B on all of Rd × R, the variational equation along any trajectory gives, for every x,
∥DSΘ (x)∥ ≤ eB T̂ , Lip DSΘ ≤ ΛB := T̂ B e2B T̂ , (A.5) the second by differentiating the variational equation once more and applying Grönwall on [0, T̂ ]; both constants depend only on (T̂ , B). Put θ := (ρ∗ + ρ̄)/2 ∈ (ρ∗ , ρ̄). The matrix M lies in the compact family {A : ρ(A) ≤ ρ∗ , ∥A∥ ≤ eB T̂ }, so the adapted-norm construction in the proof of Lemma 4.9 applies. That construction is pure linear algebra, and works for any compact family of square matrices of spectral radius at most ρ∗ and norm at most a fixed bound. There it is applied to the (d−1)×(d−1) section matrices, here to the full d×d monodromy. In particular it is applied for d = 2, where the one-dimensional shortcut of that proof does not. It supplies a norm ∥·∥⋄ with P −i ∥M ∥⋄ ≤ θ and ∥·∥ ≤ ∥·∥⋄ ≤ c⋄ ∥·∥, the lower constant being 1 because the sum i θ M ix starts at i = 0. We write c⋄ rather than c2 , which is bound elsewhere in this appendix to the section norm; c⋄ is uniform over the family and depends only on (d, T̂ , B, ρ∗ ). For a matrix A this gives ∥A∥⋄ ≤ c⋄ ∥A∥ in the induced norms. Set R :=

ρ̄ − θ , c⋄ ΛB

B⋄ := {y : y − x†Θ

⋄

≤ R}.

LONG-TIME APPROXIMATION VIA SA-NODES

For y ∈ B⋄ we have y − x†Θ ≤ y − x†Θ

⋄

45

≤ R, so by (A.5) and the choice of R,

∥DSΘ (y)∥⋄ ≤ ∥M ∥⋄ + c⋄ ∥DSΘ (y) − M ∥ ≤ θ + c⋄ ΛB R ≤ ρ̄. The set B⋄ is convex and SΘ x†Θ = x†Θ , so the mean value inequality along the segment from x†Θ to y, which lies in B⋄ , gives SΘ y − x†Θ ≤ ρ̄ y − x†Θ . Hence B⋄ is invariant under SΘ and, by ⋄ ⋄ induction, SΘn y − x†Θ

⋄

≤ ρ̄ n y − x†Θ

⋄

,

y ∈ B⋄ , n ∈ N.

(A.6)

≤ c⋄ r0 = R, so x0 ∈ B⋄ ; this is the Define r0 := R/c⋄ . If x0 − x†Θ ≤ r0 then x0 − x†Θ ⋄ direction in which the equivalence is used, and it costs c⋄ . Step 2: from the iterates to the flow. Fix t ≥ 0 and write t = nT̂ + s with n = ⌊t/T̂ ⌋ and † s ∈ [0, T̂ ). By (4.12), ΦΘ (t; x0 ) = ΦΘ (s; SΘn x0 ) and ΦΘ (s; x†Θ ) ∈ γΘ . Both points are flowed for the same time s ≤ T̂ by the same field, whose spatial Lipschitz constant is at most B, so Grönwall gives ΦΘ (t; x0 ) − ΦΘ (s; x†Θ )

≤ eB T̂ SΘn x0 − x†Θ

≤ eB T̂ ρ̄ n x0 − x†Θ

⋄

≤ c⋄ eB T̂ ρ̄ n x0 − x†Θ ,

using (A.6) and then the two equivalences in the directions in which they are needed. † lies within εΓ of Γ, in particular Step 3: comparison with Γ. By (B2) every point of γΘ † dist(ΦΘ (s; xΘ ), Γ) ≤ εΓ . The triangle inequality for the distance to the compact set Γ then gives

dist ΦΘ (t; x0 ), Γ ≤ ΦΘ (t; x0 ) − ΦΘ (s; x†Θ ) + dist ΦΘ (s; x†Θ ), Γ , which with Step 2 is (4.15) for C := c⋄ eB T̂ ≥ 1. All constants depend only on (d, T̂ , B, ρ∗ ).

□

Acknowledgments The authors sincerely thank Enrique Zuazua for suggesting this line of research and for his detailed feedback on earlier versions of the manuscript, which shaped the revision. Ziqian Li is supported by the Alexander von Humboldt Professorship program. References [1] F. L. Bauer and C. T. Fike, Norms and exclusion theorems, Numer. Math., 2 (1960), pp. 137–141. [2] H. G. Bock, Numerical treatment of inverse problems in chemical reaction kinetics, in Modelling of Chemical Reaction Systems: Proceedings of an International Workshop, Heidelberg, Fed. Rep. of Germany, September 1–5, 1980, Springer, 1981, pp. 102–125. [3] H. G. Bock and K.-J. Plitt, A multiple shooting algorithm for direct solution of optimal control problems, IFAC Proceedings Volumes, 17 (1984), pp. 1603–1608. [4] M. Calvo, M. P. Laburta, J. I. Montijano, and L. Rández, Error growth in the numerical integration of periodic orbits, Math. Comput. Simulation, 81 (2011), pp. 2646–2661. [5] E. F. Camacho, C. Bordons, and J. M. Maestre, Model predictive control, Advanced Textbooks in Control and Signal Processing, Springer, Cham, third ed., 2026. [6] B. Cano and J. M. Sanz-Serna, Error growth in the numerical integration of periodic orbits, with application to Hamiltonian and reversible systems, SIAM J. Numer. Anal., 34 (1997), pp. 1391–1417. [7] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, Neural ordinary differential equations, in Advances in Neural Information Processing Systems, vol. 31, 2018. [8] C. Chicone, Ordinary differential equations with applications, vol. 34 of Texts in Applied Mathematics, Springer, Cham, third ed., 2024. [9] C. Esteve-Yagüe and B. Geshkovski, Sparsity in long-time control of neural ODEs, Systems Control Lett., 172 (2023), pp. Paper No. 105452, 14. [10] G. Floquet, Sur les équations différentielles linéaires à coefficients périodiques, Ann. Sci. École Norm. Sup. (2), 12 (1883), pp. 47–88.

46

Z. LI AND N. M. MATZAKOS

[11] B. Geshkovski and E. Zuazua, Turnpike in optimal control of PDEs, ResNets, and beyond, Acta Numer., 31 (2022), pp. 135–263. [12] E. Hairer, C. Lubich, and G. Wanner, Geometric numerical integration, vol. 31 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, second ed., 2006. Structure-preserving algorithms for ordinary differential equations. [13] Q. Kang, Y. Song, Q. Ding, and W. P. Tay, Stable neural ode with lyapunov-stable equilibrium points for defending against adversarial attacks, in Advances in Neural Information Processing Systems, M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, eds., vol. 34, Curran Associates, Inc., 2021, pp. 14925–14937. [14] H. B. Keller, Numerical methods for two-point boundary-value problems, Blaisdell Publishing Co. [Ginn and Co.], Waltham, Mass.-Toronto, Ont.-London, 1968. [15] J. Z. Kolter and G. Manek, Learning stable deep dynamics models, in Advances in Neural Information Processing Systems, H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, eds., vol. 32, Curran Associates, Inc., 2019. [16] Z. Li, K. Liu, L. Liverani, and E. Zuazua, Universal approximation of dynamical systems by semiautonomous neural ODEs and applications, SIAM J. Numer. Anal., 64 (2026), pp. 193–223. [17] Z. Li, K. Liu, Y. Song, H. Yue, and E. Zuazua, Deep neural ODE operator networks for PDEs, Math. Models Methods Appl. Sci., 36 (2026), pp. 1739–1782. [18] K. Liu and E. Zuazua, Representation and regression problems in neural networks: relaxation, generalization, and numerics, Math. Models Methods Appl. Sci., 35 (2025), pp. 1471–1521. [19] S. Massaroli, M. Poli, S. Sonoda, T. Suzuki, J. Park, A. Yamashita, and H. Asama, Differentiable multiple shooting layers, in Advances in Neural Information Processing Systems, M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, eds., vol. 34, Curran Associates, Inc., 2021, pp. 16532–16544. [20] N. M. Matzakos, Activation saturation and Floquet spectrum collapse in neural ODEs, arXiv preprint arXiv:2604.00543, (2026). [21] N. Namura, T. Ishii, and H. Nakao, Designing 2-D limit-cycle oscillators with prescribed trajectories and phase-response characteristics, IEEE Trans. Automat. Control, 69 (2024), pp. 3144–3157. [22] F. Nawaz, T. Li, N. Matni, and N. Figueroa, Learning complex motion plans using neural odes with safety and stability guarantees, in 2024 IEEE International Conference on Robotics and Automation (ICRA), 2024, pp. 17216–17222. [23] I. D. J. Rodriguez, A. Ames, and Y. Yue, LyaNet: A Lyapunov framework for training neural ODEs, in Proceedings of the 39th International Conference on Machine Learning, K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu, and S. Sabato, eds., vol. 162 of Proceedings of Machine Learning Research, PMLR, 17–23 Jul 2022, pp. 18687–18703. [24] D. Ruiz-Balet and E. Zuazua, Neural ODE control for classification, approximation, and transport, SIAM Rev., 65 (2023), pp. 735–773. [25] Á. Ságodi and I. M. Park, Universal approximation theorems for dynamical systems with infinite-time horizon guarantees, 2026. arXiv preprint arXiv:2602.08640. [26] E. M. Turan and J. Jäschke, Multiple shooting for training neural differential equations on time series, IEEE Control Systems Letters, 6 (2022), pp. 1897–1902. [27] D. W. M. Veldman, A. Borkowski, and E. Zuazua, Stability and convergence of a randomized model predictive control strategy, IEEE Trans. Automat. Control, 69 (2024), pp. 6253–6260. [28] D. Wilson, Data-driven model identification using forcing-induced limit cycles, Phys. D, 459 (2024), pp. Paper No. 134013, 15. [29] H. Zhang, L. Cheng, Z. Liu, and Y. Zhang, Learning orbitally stable dynamics via transverse contraction criteria for modeling periodic tasks, Int. J. Robotics Research, 45 (2026).

Related documents

• DeepSeekMath Meets Order Book: Group-Aware Policy Optimization for High-Frequency Directional Trading#630019
• DeepSeekMath Meets Order Book: Group-Aware Policy Optimization for High-Frequency Directional Trading#628157
• pANO-F12: An atomic natural orbital-inspired route to more compact basis sets for F12 explicitly correlated methods#1077893
• Predict Before Replay: Joint FEC and Flight Control for Reliable Scale-Up Links#463189
• Network-Assisted Dynamic Adaptation (NADA): A Unified Congestion Control Scheme for Real-Time Media#502336