arXiv CSopen access

Residual-loss Anomaly Analysis of Physics-Informed Neural Networks: An Inverse Method for Change-point Detection in Nonlinear Dynamical Systems with Regime Switching

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

machine learning, deep learning, neural networks

Residual-loss Anomaly Analysis of Physics-Informed Neural Networks: An Inverse Method for Change-point Detection in Nonlinear Dynamical Systems with Regime Switching Yuhe Baia , Chengli Tanb , Jiaqi Lia , Xiangjun Wanga , Zhikun Zhangb,˚ a School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan, 430074, China

arXiv:2604.25655v1 [stat.ML] 28 Apr 2026

b School of Mathematics and Statistics, Northwestern Polytechnical University, Xi’an, 710072, China

Abstract Nonlinear dynamical systems with regime transitions are typically described by ordinary differential equations with jumping parameters parameters. Traditional methods often treat change-point detection and parameter estimation as separate tasks, ignoring the inherent coupling between them. To address this, we propose residual-loss anomaly analysis of physics-informed neural networks, a unified framework that leverages dynamical consistency within the physics-informed learning paradigm. This approach jointly infers piecewise parameters and transition points under a single set of constraints. The method follows a two-stage strategy: First, local physical residuals are analyzed through overlapping subinterval decomposition. When a subinterval spans a true transition point, the residual exhibits a distinct structural elevation in noise-free conditions, which has a non-zero lower bound, enabling effective localization of potential transition intervals. Second, within our framework, change-point locations and piecewise parameters are integrated into a unified physical loss function for joint optimization, enabling simultaneous identification. Experiments on benchmark nonlinear dynamical systems, including Malthusian and logistic growth models, Van der Pol oscillator, Lotka-Volterra model and Lorenz system, demonstrate that the proposed method outperforms traditional decoupled approaches in both change-point localization and parameter estimation accuracy. This study provides an efficient, unified solution for structurally coupled inverse problems in nonlinear dynamical systems with regime switching. Keywords: Nonlinear Dynamical System; Parameter Estimation; Inverse Problem; Physics-Informed Neural Network; Overlapping-Domain Decomposition; Change-Point Detection.

˚ Corresponding author

Email address: [email protected] (Zhikun Zhang )

Preprint submitted to Elsevier

April 29, 2026

Contents

1 Introduction

2 Nonlinear Dynamical Systems with Regime Switching

2.1

Mathematical Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.2

Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3 Methodology

3.1

Physics-Informed Machine Learning Framework . . . . . . . . . . . . . . . . . . . . . . . .

3.2

Stage I: Coarse Localization via Residual Signature . . . . . . . . . . . . . . . . . . . . . .

3.3

Stage II: Joint Inference within Candidate Interval . . . . . . . . . . . . . . . . . . . . . .

3.4

Theoretical Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4 Numerical Experiment

4.1

Malthus and Logistic Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.2

Van der Pol Oscillator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.3

Lotka-Volterra Model

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.4

Lorenz System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5 Parallel Overlapping Domain Strategy

6 Comparative Experiment

6.1

Traditional Statistical Method

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.2

Existing Neural Network Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7 Conclusion and Discussion

1. Introduction The evolutionary processes of complex systems are typically not governed by a single, stable dynamical mechanism. Instead, these systems often undergo abrupt regime transitions and dynamical reconfigurations in response to external environmental changes, internal structural reorganizations, or shifts in regulatory mechanisms [1, 2, 3]. These regime transitions are prevalent in both natural and engineered systems [4], such as critical transitions in ecosystems, functional remodeling in biological networks, mode switching in engineering system operations, and non-stationary evolutionary processes in climate and environmental systems [5]. Identifying when and how the dynamics of a system change is critical for understanding the mechanisms of complex systems and for improving the reliability of prediction and control [6, 7]. These systems are typically described by ordinary differential equations (ODEs) that govern their state evolution. Transitions in the dynamical mechanism often manifest as discontinuous jumps in model parameters over time [8, 9], which can be modeled as jump coefficients. Compared to stationary models, parameter discontinuities can alter local dynamics and cause qualitative changes in global trajectories, resulting in piecewise regime behavior [10]. However, these changes are not directly observable and can only be inferred from limited system output trajectories, making it a challenging task to infer both the system dynamics and the timing of these transitions from observational data, particularly in the absence of prior structural assumptions [11, 12, 13]. For inverse problems involving dynamical systems with regime transitions, the object to be inferred is no longer a finite-dimensional static parameter, but a time-dependent parameter function with an unknown piecewise structure, where both the segmentation locations and their corresponding values jointly determine the system’s dynamical behavior [14, 15]. This structural complexity leads to an inherent coupling between parameter estimation and transition time identification [10]. If the abrupt transition times are not accounted for, parameter estimation becomes distorted due to the confounding of dynamics across different regimes. Conversely, without an accurate parametric model, distinguishing regime switches from stochastic fluctuations is difficult, leading to unreliable change-point detection [16]. Therefore, treating structural identification and parameter inference sequentially or independently often fails to capture their interdependence [17]. This circular dependency highlights the need for integrated methods capable of simultaneously inferring change-point locations and jump parameters [18]. Such integration is both a theoretical necessity and a foundation for understanding transition mechanisms in systems and achieving accurate prediction and intervention. Regime-transition problems in dynamical systems are often formalized as structural break models or piecewise parameter models, where the unknown change-point locations are treated as discrete structural variables, and the piecewise dynamical parameters are assumed to be locally stationary or constant within each subinterval [19, 20, 21]. Classical methods estimate change-point locations and piecewise parameters 3

jointly or alternately using techniques such as generalized likelihood ratio statistics, penalized maximum likelihood functions, or sparsity-regularized objective functions, with consistency and convergence rate analyses under certain regularity conditions [22, 23]. Bayesian methods incorporate the number and locations of change points into a unified posterior framework and use methods such as reversible jump Markov chain Monte Carlo or sequential Monte Carlo for joint inference of model dimensionality and parameter space [24, 25, 26, 27]. However, these methods often rely on pre-specified piecewise forms or searchable structural spaces, decomposing change-point localization and parameter estimation into relatively independent or sequentially optimized steps [28, 29]. In recent years, data-driven methods have provided new approaches for modeling dynamical systems through high-dimensional function approximation and end-to-end differentiable optimization frameworks [30, 31]. These approaches represent system evolution as nonlinear mappings or implicit probabilistic models and leverage automatic differentiation and gradient-based optimization for unified training [32]. Building on this, physics-informed neural networks (PINNs) embed differential equation residuals into the loss function, so that dynamical constraints serve as differentiable regularization terms in the optimization process, enabling simultaneous state reconstruction and parameter inversion [33]. However, in systems with regime transitions, these methods generally assume parameter continuity over the training interval or a known piecewise structure, or they rely on external rules to pre-partition the time domain [34]. This results in a relatively decoupled workflow between change-point detection and parameter learning [35]. While statistical inference and data-driven methods differ significantly in their theoretical tools and numerical implementations, they often share an implicit assumption when addressing regime transition problems: treating ”structural identification” and ”parameter inference” as independent tasks that can be optimized separately [36]. However, when transition times and parameter values are structurally coupled, such a decoupled approach can undermine both consistency and stability [37]. Segmentation without dynamical constraints struggles to accurately capture regime differences, while parameter inversion that ignores structural uncertainty is highly sensitive to misspecified intervals [38]. Therefore, it is essential to re-examine the interplay between structural changes and parameter learning under a unified dynamical constraint, overcoming the limitations of decoupled paradigms [39, 40]. The above analysis indicates that the core challenge lies not in designing more sophisticated segmentation strategies or more refined parameter estimators, but in uncovering the intrinsic relationship between structural changes and parameter evolution under a unified dynamical constraint [41]. When a single dynamical model is used to fit a time interval that spans a true transition point, the structural mismatch caused by regime differences will inevitably manifest as a systematic deviation in the differential equation residuals [42]. In an ideal noise-free scenario, this mismatch typically corresponds to an irreducible non-zero residual level. Therefore, this residual reflects structural model misspecification rather than stochastic noise [43, 44].

Based on the above analysis, this paper proposes a unified inference framework that uses dynamical consistency as an intrinsic signal to identify regime transitions, simultaneously capturing structural changes and parameter evolution within the same dynamical system. Specifically, the method adopts a two-stage strategy: In the first stage, a local dynamical consistency measure is constructed using overlapping subintervals and steady-state statistics to robustly localize potential transition intervals. In the second stage, a differentiable change-point representation is introduced within candidate intervals, which makes the change-point locations optimizable variables to be solved jointly with piecewise parameter inversion under a unified loss function. This design leverages the discriminative power of structural residual signals while maintaining the physical consistency of parameter inversion, enabling accurate and simultaneous identification of both transition locations and jump parameters. The framework provides a unified approach to characterizing regime reconfigurations in complex systems and offers new insights for solving structurally coupled inverse problems. The framework proposed in this paper offers a novel approach for identifying regime transitions in dynamical systems and demonstrates broad application potential, particularly in modeling and predicting complex systems such as biological regulation, financial markets, and climate change. This method holds significant theoretical and practical value. By using this framework, we can more accurately capture abrupt behavioral changes in systems, providing a powerful tool for interdisciplinary scientific research. As our understanding of complex systems deepens, the proposed method is expected to lay the foundation for further research and applications, driving theoretical innovation and technological advancements in these fields. The following sections detail the implementation of the framework, including the specific steps of the two-stage change-point detection, its validation, and its application across multiple dynamical systems. This paper is organized as follows: Section 2 formulates nonlinear dynamical systems with jump coefficients and establishes the well-posedness of the inverse problem. Section 3 develops the unified, physics-consistent inference framework, introducing the two-stage strategy that leverages residual-induced structural mismatch for joint change-point localization and parameter estimation. Section 4 evaluates the proposed framework on representative dynamical systems and presents a detailed analysis of the results. Section 5 introduces a parallel overlapping-domain implementation for scalable computation. Section 6 compares the proposed method with traditional statistical approaches and existing PINNsbased techniques. Finally, Section 7 concludes with a discussion of future directions.

2. Nonlinear Dynamical Systems with Regime Switching In this section, nonlinear dynamical systems with jump coefficients are formulated, and the theoretical foundation of the associated inverse problem is established. A parametric ordinary differential equation with piecewise constant parameters governed by a discrete switching process is introduced. The existence 5

and uniqueness of solutions are then proved, providing a rigorous basis for subsequent joint parameter estimation and change-point detection.

2.1. Mathematical Setup In this section, the basic definition of nonlinear dynamical systems with jump coefficients is introduced. A mathematical framework for general parametric ODEs is first established. Let T ą 0 denote the terminal time. Consider a function x : r0, T s Ñ Rn satisfying the following parametric ordinary differential equation: dxptq “ f pt, xptqqdt,

(2.1)

where f : r0, T s ˆ Rn Ñ Rn , together with the initial condition pt0 , x0 q P G Ă Rn`1 . Within the PINNs framework, parameter identification and change-point detection problems can be formulated as inverse problems. We study statistical inference for system parameters based on observations of the system output. Suppose that (2.1) depends on a finite number of unknown parameters, and rewrite it in the parametric form dxptq “ f pt, xptq; θq dt,

(2.2)

where θ “ pθ1 , θ2 , . . . , θd q is a d-dimensional parameter vector. To model parameter switching over time, we introduce a discrete state variable rptq P I :“ t1, 2, . . . , Ku. Let Θ “ tθ 1 , θ 2 , . . . , θ K u be the associated set of parameter regimes, and define the active parameter by θptq “ θ rptq . Assume that rptq is a rightcontinuous piecewise constant function on rt0 , T q. Given a partition t0 ă t1 ă ¨ ¨ ¨ ă tn ă tn`1 “ T , we write rptq “

n ÿ

rk 1rtk ,tk`1 q ptq,

(2.3)

k“0

where 1rtk ,tk`1 q ptq is the indicator function of rtk , tk`1 q. We then define the vector field with jumping coefficients as follows. f r pt, xptq, rptq; Θq “ f pt, xptq; θ rptq q,

(2.4)

so that the nonlinear dynamical system with jumping coefficients takes the form dxptq “ f r pt, xptq, rptq; Θq dt,

(2.5)

where f r : rt0 , T q ˆ Rn ˆ I Ñ Rn , and pt0 , x0 q P G Ă Rn`1 . Before performing parameter inversion and change-point detection for the proposed nonlinear dynamical system with jumping coefficients, we establish its well-posedness, namely, the existence and uniqueness

of solutions. To this end, we apply the classical Picard-Lindelöf theorem [45, 46] to prove the following existence and uniqueness result for the proposed system. Theorem 2.1 (Existence and Uniqueness Theorem for Nonlinear Dynamical Systems). Consider the nonlinear dynamical system defined by (2.5). If the following conditions hold. 1. (Lipschitz Condition): There exists a constant L such that for all t P rt0 , T s, r P I and any x1 , x2 P Rn , the following holds: }f r pt, x1 , rptq; Θq ´ f r pt, x2 , rptq; Θq} ď L}x1 ´ x2 }.

(2.6)

2. (Linear Growth Condition): There exists a constant K such that for all t P rt0 , T s, r P I and x P Rn , the following holds. }f r pt, x, rptq; Θq} ď Kp1 ` }x}q.

(2.7)

Then, in a neighborhood of t0 , there exists a unique piecewise continuously differentiable solution x : pt0 ´ ϵ, t0 ` ϵq Ñ Rn , where ϵ ą 0, that satisfies equation(2.5) and the initial condition xpt0 q “ x0 , $ &dxptq “ f r pt, xptq, rptq; Θqdt, %

(2.8)

xpt0 q “ x0 .

The solution depends continuously on t and the parameter rptq, and the solution is locally unique around t0 for any rptq that satisfies the conditions. Proof. Since rptq takes values in a finite set I and is right-continuous, it is a piecewise constant function on rt0 , T s with finitely many jumps. Hence, there exist switching points tτi uki“1 with k ă 8 such that t0 “ τ0 ă τ1 ă ¨ ¨ ¨ ă τk ă τk`1 “ T,

(2.9)

and rptq “ ri ,

t P rτi , τi`1 q,

i “ 0, 1, . . . , k,

(2.10)

where ri P I. For each i “ 0, 1, . . . , k, on the interval rτi , τi`1 q, equation (2.5) reduces to dxptq “ f ri pt, xptq, ri ; Θq dt “ f pt, xptq; θ ri q dt.

(2.11)

By the Lipschitz condition and linear growth condition, the classical Picard-Lindelöf theorem implies

that, for each i, there exists a unique continuously differentiable solution xri : rτi , τi`1 q Ñ Rn ,

(2.12)

satisfying (2.11) with initial value xpτi q. Starting from the initial condition xpt0 q “ x0 , we construct the solution successively: for i “ 0, we obtain xr0 on rτ0 , τ1 q. For i ě 1, we take xpτi q as the initial value and obtain xri on rτi , τi`1 q. Define xptq “ xri ptq,

t P rτi , τi`1 q,

i “ 0, 1, . . . , k.

(2.13)

Then xptq is well-defined on rt0 , T s and satisfies equation (2.5). Uniqueness follows from the uniqueness of solutions to (2.11) on each subinterval and the consistency of initial conditions at τi . Therefore, equation (2.5) admits a unique solution on rt0 , T s. Remark. Theorem 2.1 establishes the well-posedness of the dynamical system under piecewise constant switching coefficients. In particular, the existence and uniqueness of a global solution ensure that the induced piecewise dynamics are mathematically consistent. This result is essential for the subsequent analysis of parameter estimation and change-point detection.

2.2. Problem Statement Building upon the formulation of nonlinear dynamical systems with jump coefficients, we consider the following inverse problem. Let xptq denote the system state governed by dxptq “ f pt, xptq; θptqqdt,

(2.14)

where the parameter θptq is assumed to follow an unknown piecewise constant structure, characterized pkq K`1 by a set of change-point locations tτk uK uk“1 . k“1 and the corresponding regime-specific parameters tθ

Given partial observations txd pti quN i“1 of the system trajectory, the objective is to jointly recover the change-point locations and the associated parameters, namely τ̂k and θ̂

pkq

The intrinsic entanglement between change-point detection and parameter estimation poses a fundamental challenge in this setting. Unknown change points distort parameter inference across regimes, while misestimated parameters can in turn obscure the true transitions. We refer to this coupled inverse task as change-point detection in nonlinear dynamical systems (NDS-CPD). To address it, we propose a unified physics-informed framework that treats change points and regime-specific parameters as coupled variables, and exploits residual-driven structural inconsistencies as endogenous signals for joint and data-efficient inference.

3. Methodology In this section, we introduce the residual-loss anomaly analysis framework of physics-informed neural networks (RAA-PINNs), our inverse framework for change-point detection and parameter estimation in nonlinear dynamical systems with regime switching. The main objective is to recover both the changepoint locations and the regime-specific parameters from observational data in a unified manner. To present the main idea clearly, we first consider the case of a single change-point. Let θptq denote the time-dependent parameter vector. Assume that there exists an unknown τ P p0, T q such that $ ’ &θ ´ ,

t ă τ,

’ %θ ` ,

t ě τ,

θptq “

θ´ ‰ θ` ,

(3.1)

where θ ´ and θ ` denote the parameter vectors before and after the transition, respectively. The inverse task is therefore to infer the piecewise parameters pθ ´ , θ ` q together with the unknown change-point τ . The key idea of RAA-PINNs is that change points manifest as anomalies in the physics loss. On intervals contained in a single regime, locally constant parameters allow accurate fitting of the dynamics. On intervals that cross a true change point, the same assumption creates an unavoidable residual mismatch. We therefore use the residual-loss anomaly as the basic signal for change-point detection. Accordingly, RAA-PINNs follows a two-stage strategy. Stage I uses overlapping subintervals and local physics-constrained training to identify a candidate interval from the residual-loss profile. Stage II then refines the result by jointly optimizing the differentiable change-point variable and the piecewise parameters θ ´ and θ ` . This combines coarse screening with local refinement in a unified inverse framework. Figure 1 provides a schematic overview of the proposed two-stage workflow. Using the Lorenz system as an illustrative example, the figure shows how overlapping-domain decomposition is first used in Stage I to detect high-residual subintervals and how the resulting candidate interval is then refined in Stage II through joint optimization of the change-point location and the piecewise parameter vectors. 3.1. Physics-Informed Machine Learning Framework We first introduce the local physics-informed inverse formulation that serves as the basic building block of RAA-PINNs. On a given time interval I “ r0, T s, we consider the nonlinear dynamical system in (2.2), where x : I Ñ Rn denotes the unknown state trajectory and θ P Θ Ă Rm is the unknown constant parameter vector on that interval, with Θ denoting the admissible parameter set. We approximate x by a neural network xϕ : I Ñ Rn , where ϕ P RM collects the trainable weights and biases. The network parameters ϕ and the system parameter vector θ are optimized jointly. d Let Id “ ttj uN j“1 Ă I denote the set of observation times, and let xd ptj q be the corresponding

observations. For any sufficiently smooth function u : I Ñ Rn and any θ P Θ, we define the physics 9

Figure 1: Schematic illustration of the proposed RAA-PINNs framework for change-point detection and parameter estimation via overlapping-domain decomposition.

residual and the data residual by Rint ru, θsptq “

` ˘ du ptq ´ f t, uptq; θ , dt

Rdata rusptj q “ uptj q ´ xd ptj q,

t P I, (3.2)

tj P Id .

These residuals quantify how well the pair pu, θq satisfies the governing equation and fits the observed data. In particular, if x is the exact solution corresponding to the true parameter vector θ, then Rint rx, θsptq “ 0 for t P I and Rdata rxsptj q “ 0 for tj P Id . The PINN estimator is obtained by jointly optimizing over ϕ and θ so as to balance fidelity to the governing dynamics and consistency with the observations. In an abstract form, this can be written as pϕ̂, θ̂q “ arg min p}Rint rxϕ , θs}Y ` }Rdata rxϕ s}Z q ,

(3.3)

pϕ,θq

where Y and Z are normed spaces for the physics and data residuals, respectively. In this work, we take Y “ Lpy pIq and Z “ ℓpz on the observation set Id , where 1 ď py , pz ă 8. This yields the continuous objective ˜ż pϕ̂, θ̂q “ arg min pϕ,θq

Nd ÿ › › › › ›Rint rxϕ , θsptq›py dt ` ›Rdata rxϕ sptj q›pz

¸ .

(3.4)

j“1

int In practice, the continuous physics residual term is approximated by quadrature. Let Sint “ tsi uN i“1 Ă I

int denote the collocation set, with positive quadrature weights twi uN i“1 . For the data term, we introduce d positive weights tvj uN j“1 . The resulting discrete loss is

Jpϕ, θq “

Nd ÿ

N int ÿ › ›pz › ›py vj ›Rdata rxϕ sptj q› ` λ wi ›Rint rxϕ , θspsi q› ,

j“1

i“1

(3.5)

where λ ą 0 balances the data-misfit and physics-residual terms. To further stabilize training, we may optionally include a regularization term on the network parameters as ´ ¯ pϕ̂, θ̂q “ arg min Jpϕ, θq ` λreg Jreg pϕq ,

(3.6)

pϕ,θq

where Jreg : RM Ñ R is a regularization functional and λreg ě 0 is the corresponding penalty parameter. A common choice is Jreg pϕq “ }ϕ}qq , where q “ 2 gives standard L2 regularization, while q “ 1 can be used to promote sparsity.

3.2. Stage I: Coarse Localization via Residual Signature To obtain a coarse localization of the change-point, we partition the temporal domain into overlapping subintervals. Let 0 “ t0 ă t1 ă ¨ ¨ ¨ ă tK “ T,

(3.7)

be a partition of I “ r0, T s, and let δ ą 0 denote the overlap width. For each k “ 1, . . . , K, we define Ik “ rtk´1 ´ δ, tk ` δs X r0, T s.

(3.8)

By construction, every change-point τ P p0, T q is contained in at least one such subinterval. At this stage, the subinterval problems are treated independently. No interface or continuity conditions are imposed across overlapping regions, because the goal of Stage I is change-point screening rather than global reconstruction. On each Ik , we train a local PINN xϕk : Ik Ñ Rn together with a locally constant parameter vector θ k P Θ. Let pkq N r

k Skr “ tsi ui“1 Ă Ik ,

pkq N d

k Skd “ ttj uj“1 Ă Ik X Id ,

(3.9) pkq N r

pkq N d

pkq › pkq ›py wi ›Rint rxϕk , θ k spsi q› .

(3.10)

k k which denote the collocation set and the observation set on Ik , respectively. Let twi ui“1 and tvj uj“1

be the associated positive weights. The local training objective is Nkd

Jk pϕk , θ k q “

pkq › pkq ›pz vj ›Rdata rxϕk sptj q› ` λ

j“1

Nkr

i“1

Since the local problems are uncoupled, all subinterval PINNs can be trained in parallel. The overlap is introduced to improve the robustness of change-point localization near subinterval boundaries. After training on Ik , we compute the physics residual energy Nkr

Ek “

pkq › pkq ›py wi ›Rint rxϕk , θ k spsi q› .

(3.11)

i“1

To make these energies comparable across subintervals, we normalize them by the total collocation weight and define Ek Ēk “ řN r pkq . k i“1 wi

(3.12)

Because PINN training may exhibit transient oscillations or occasional spikes, we summarize the terminal behavior of the normalized residual over the last M iterations. Let Ēk,ℓ denote the value of (3.12) at iteration ℓ, and let W “ tL ´ M ` 1, . . . , Lu,

(3.13)

where L is the total number of training iterations. We then define Sk “ mediantĒk,ℓ : ℓ P Wu.

(3.14)

This score provides a robust summary of the terminal residual level on subinterval Ik . The key observation underlying RAA-PINNs is that Sk tends to be elevated when Ik contains the change-point. If Ik lies entirely on one side of τ , then the dynamics on Ik are governed by a single parameter regime, and a locally constant parameter vector can fit the subinterval well, leading to a relatively small residual level. By contrast, if τ P Ik , then the dynamics on Ik involve two distinct parameter regimes associated with θ ´ and θ ` . Any single constant parameter vector θ k must then compromise between the two regimes, producing a residual-loss anomaly that cannot be removed by training. To formalize this residual-loss anomaly mechanism, we introduce the following assumptions. These conditions are tailored to the present regime-switching setting and are closely related to standard identifiability assumptions in parameter estimation and consistency requirements in collocation and physicsinformed formulations [47, 48]. Assumption 1 (Parameter-affine structure and local identifiability). Let x‹ : I Ñ Rn denote the true trajectory. Assume that the vector field f : I ˆ Rn ˆ Θ Ñ Rn is affine in the parameter vector θ, namely, f pt, x; θq “ Gpt, xqθ ` bpt, xq,

(3.15)

where G : I ˆ Rn Ñ Rnˆm ,

b : I ˆ Rn Ñ Rn .

(3.16)

Moreover, assume that there exists a constant α ą 0 such that for every subinterval J Ă I with positive length, the associated information matrix ż

` ˘J ` ˘ G t, x‹ ptq G t, x‹ ptq dt,

M pJq “

(3.17)

satisfies ` ˘ λmin M pJq ě α|J|,

(3.18)

where |J| denotes the length of J and λmin p¨q denotes the smallest eigenvalue. Assumption 2 (Quadrature consistency of the interior residual). For each subinterval Ik , let Skr “ pkq N r

pkq N r

k k tsi ui“1 Ă Ik be the collocation set, where Nkr is the number of collocation points in Ik , and let twi ui“1

be the associated positive weights. Assume that the corresponding quadrature rule is consistent on Ik , in the sense that for any continuous function h : Ik Ñ R, Nkr

lim r

Nk Ñ8

ż pkq

` pkq ˘ h si “

hptq dt.

(3.19)

i“1

Assumption 1 excludes degenerate regimes in which the dynamics are locally insensitive to the unknown parameters. Assumption 2 ensures that the discrete PINN residual provides a faithful approximation to its continuous-time counterpart. Under these conditions, the next theorem explains why the interior residual remains systematically larger on subintervals that cross the change-point. We denote by Rk the corresponding continuous-time residual energy on Ik . Theorem 3.1 (Residual lower bound on change-point subintervals). Under Assumption 1, let Ik be any Stage I subinterval and define Ik´ “ Ik X r0, τ q and Ik` “ Ik X rτ, T s. For any constant parameter vector θ P Θ, define the continuous-time residual energy › ‹ › › dx ` ‹ ˘›2 › dt, › ptq ´ f t, x ptq; θ › dt › Ik 2

ż Rk pθq “

(3.20)

and let Rk “ inf Rk pθq. θPΘ

(3.21)

If τ R Ik , then either Ik´ “ H or Ik` “ H, and hence Rk “ 0 in the absence of observational noise. If τ P Ik and θ ´ ‰ θ ` , then Rk ě α

|Ik´ | |Ik` | }θ ´ ´ θ ` }22 . |Ik´ | ` |Ik` |

(3.22)

In particular, α mint|Ik´ |, |Ik` |u }θ ´ ´ θ ` }22 . 2

Rk ě

(3.23)

Proof. We consider separately the cases τ R Ik and τ P Ik . First, suppose that τ R Ik . If Ik lies entirely on one side of the change-point, then the true dynamics on Ik are governed by a single constant parameter vector θ true P tθ ´ , θ ` u, which means ` ˘ dx‹ ptq “ f t, x‹ ptq; θ true , dt

@ t P Ik .

(3.24)

In the noise-free setting, choosing θ “ θ true in (3.20) yields ż

› ‹ ˘› ` ›x9 ptq ´ f t, x‹ ptq; θ true ›2 dt “ 0. 2

(3.25)

Rk “ inf Rk pθq “ 0.

(3.26)

Rk pθ true q “ Ik

Hence, θPΘ

We now turn to the case τ P Ik and θ ´ ‰ θ ` . When the subinterval Ik crosses the change-point, the true trajectory satisfies the piecewise dynamics

x9 ‹ ptq “

$ ˘ ` ’ &f t, x‹ ptq; θ ´ ,

t P Ik´ ,

’ %f `t, x‹ ptq; θ ` ˘,

t P Ik` .

(3.27)

For any candidate constant parameter vector θ P Θ, the corresponding residual energy admits the decomposition ż Rk pθq “ I

› ‹ ` ˘› ›x9 ptq ´ f t, x‹ ptq; θ ›2 dt 2

żk “

Ik´

› ` ‹ ˘ ` ˘› ›f t, x ptq; θ ´ ´ f t, x‹ ptq; θ ›2 dt ` 2

ż Ik`

› ` ‹ ˘ ` ˘› ›f t, x ptq; θ ` ´ f t, x‹ ptq; θ ›2 dt.

(3.28)

By Assumption 1, the vector field is affine in the parameter vector and can be written as f pt, x; θq “ Gpt, xqθ ` bpt, xq.

(3.29)

Substituting this representation into (3.28) and canceling the common offset term bpt, x‹ ptqq yields ż Rk pθq “

Ik´

› › ` ‹ ˘ ´ ›G t, x ptq pθ ´ θq›2 dt ` 2

ż Ik`

› › ` ‹ ˘ ` ›G t, x ptq pθ ´ θq›2 dt 2

“ pθ ´ ´ θqJ M pIk´ qpθ ´ ´ θq ` pθ ` ´ θqJ M pIk` qpθ ` ´ θq,

(3.30)

where, for any interval J Ă I, ż M pJq “

` ˘J ` ˘ G t, x‹ ptq G t, x‹ ptq dt

(3.31)

denotes the associated information matrix. Invoking Assumption 1, there exists α ą 0 such that for any relevant interval J Ă I, ` ˘ λmin M pJq ě α|J|.

(3.32)

Rk pθq ě α|Ik´ | }θ ´ ´ θ}22 ` α|Ik` | }θ ` ´ θ}22 .

(3.33)

Applying this bound to (3.30) gives

The right-hand side of (3.33) is a strictly convex quadratic function of θ. Define Qpθq “ |Ik´ | }θ ´ ´ θ}22 ` |Ik` | }θ ` ´ θ}22 .

(3.34)

A direct calculation shows that the unique minimizer of Q over Rm is θ̄ k “

|Ik´ |θ ´ ` |Ik` |θ ` . |Ik´ | ` |Ik` |

(3.35)

|Ik´ | |Ik` | }θ ´ ´ θ ` }22 . |Ik´ | ` |Ik` |

(3.36)

Substituting θ̄ k back into Q yields infm Qpθq “

θPR

Since Θ Ă Rm , it follows that inf Qpθq ě infm Qpθq.

θPΘ

θPR

(3.37)

Combining this inequality with (3.33) and taking the infimum over θ P Θ gives Rk “ inf Rk pθq ě α θPΘ

|Ik´ | |Ik` | }θ ´ ´ θ ` }22 , |Ik´ | ` |Ik` |

(3.38)

which proves (3.22). Finally, since 1 ab ě minta, bu, a`b 2

a, b ą 0,

(3.39)

we also obtain (3.23). This completes the proof. Remark. Theorem 3.1 provides the theoretical basis for Stage I. It shows that any subinterval containing

the true change-point must exhibit a strictly positive residual mismatch under a single constant-parameter fit. Hence, the elevated residual used in the screening step is not merely an empirical phenomenon, but a structural consequence of regime switching. Unlike purely statistical change-point criteria, the detection signal here arises directly from incompatibility with the governing dynamics. Theorem 3.1 concerns an idealized continuous-time quantity defined along the true trajectory. In practice, the score Sk in (3.14) is the corresponding learned discrete counterpart and is also affected by approximation error, measurement noise, and quadrature error. Under Assumption 2, together with sufficient network capacity and adequate training, these effects are reduced, so that the relative ordering of tSk uK k“1 remains informative for identifying subintervals that contain the change-point. To reduce sensitivity to global scaling and outliers across subintervals, we further standardize the scores tSk uK k“1 using median absolute deviation (Mad) normalization. Let MedpSq “ mediantS1 , . . . , SK u,

(3.40)

␣ ( MadpSq “ median |Sk ´ MedpSq| : k “ 1, . . . , K .

(3.41)

and define

We then set Zk “

Sk ´ MedpSq , MadpSq ` ε

(3.42)

where ε ą 0 is a small numerical stabilization parameter. We define the candidate index set by Kc “ tk : Zk ě γu,

(3.43)

k ‹ “ arg max Sk .

(3.44)

Ic “ Ik‹ ´1 Y Ik‹ Y Ik‹ `1 ,

(3.45)

and select kPKc

Finally, we define the candidate interval by

with the obvious boundary truncation when k ‹ P t1, Ku. By construction, Ic is a narrow region that is expected to contain the true change-point τ with high probability and is therefore used as the search region in Stage II.

3.3. Stage II: Joint Inference within Candidate Interval Stage I yields a candidate interval Ic “ rtL , tR s rather than a point estimate of the change-point. In Stage II of RAA-PINNs, we refine this coarse localization by treating the change-point as a trainable variable within a local PINN defined on Ic . To enforce the constraint τ P ptL , tR q during optimization, we introduce an unconstrained scalar variable η P R and define τ pηq “ tL ` ptR ´ tL q σpηq,

σpzq “

1 . 1 ` e´z

(3.46)

This reparameterization ensures that the change-point remains in the interior of Ic throughout the optimization process. To enable gradient-based optimization with respect to τ , we replace the discontinuous step profile with the smooth gating function ` ˘ gpt; τ, κq “ σ κpt ´ τ q ,

(3.47)

where κ ą 0 is a prescribed sharpness parameter. Larger values of κ yield a closer approximation to a hard change-point. Using this gate, we introduce the time-dependent parameter vector ` ˘ θpt; θ ´ , θ ` , ηq “ θ ´ ` pθ ` ´ θ ´ q g t; τ pηq, κ ,

(3.48)

which serves as a differentiable surrogate for the piecewise-constant parameterization in (3.1). In the limit κ Ñ 8, this representation approaches a sharp transition. We then train a new local PINN xϕ : Ic Ñ Rn jointly with the variables pθ ´ , θ ` , ηq in order to refine both the change-point location and the regime-specific parameters. Let pcq N r

c Scr “ tsi ui“1 Ă Ic ,

pcq N d

c Scd “ ttj uj“1 Ă Ic X Id

(3.49) pcq N r

pcq N d

c c denote the collocation set and the observation set on Ic , respectively. Let twi ui“1 and tvj uj“1 be the

associated positive weights. The physics residual and the data residual under the gated parameterization are defined by Rint rxϕ , θp¨qsptq “

` ˘ dxϕ ptq ´ f t, xϕ ptq; θpt; θ ´ , θ ` , ηq , dt

t P Ic ,

(3.50)

and Rdata rxϕ sptj q “ xϕ ptj q ´ xd ptj q,

tj P Id X Ic .

(3.51)

We then minimize the local objective Ncd

Jc pϕ, θ , θ

Ncr ÿ pcq › ›pz pcq › pcq pcq ›py , ηq “ vj ›Rdata rxϕ sptj q› ` λ wi ›Rint rxϕ , θp¨qspsi q› , j“1 i“1

(3.52)

where λ ą 0 plays the same balancing role as in (3.5) and (3.10). ´

Let pϕ̂, θ̂ , θ̂ , η̂q be any minimizer of (3.52). The refined change-point estimate is then given by τ̂ “ τ pη̂q, while θ̂

and θ̂

are the corresponding parameter estimates on the two sides of the transition.

By construction, Stage II converts the interval-level output of Stage I into a pointwise change-point estimate through a differentiable parameterization on the reduced temporal region Ic , while simultaneously recovering the parameter vectors before and after the transition.

3.4. Theoretical Error Analysis In this section, we analyze the approximation error of the PINN-based inverse solution for the ` ˘ ODE system dxptq “ f t, xptq; θ dt in the presence of measurement data. This analysis provides a generalization-error interpretation for the local physics-informed inverse module underlying RAAPINNs. Our argument follows the general philosophy of physics-informed inverse learning and conditionalstability-based generalization estimates [49]. Let x̂ :“ xϕ denote the learned PINN approximation of the state trajectory. Let I 1 Ă I and let d r 1 ď px , py , pz ă 8. Let S r “ tsi uN i“1 Ă I be the collocation set for the physics residual, and let S “ 1 d ttj uN j“1 Ă I be the data set associated with the measurement locations. We denote the corresponding Nd r positive quadrature weights by twi uN i“1 and tvj uj“1 , respectively.

For any interval E Ă I, we define the generalization error by EG pEq “ EG pE; ϕ, S r , S d q “ }x̂ ´ x}Lpx pEq .

(3.53)

As indicated above, the generalization error depends on the training sets S r , S d , as well as on the network parameter ϕ. We next define the training residuals as ˜ r

Er,T pθ, S q “

Nr ÿ

¸1{py › ›py wi ›Rint rx̂, θspsi q›

˜ ,

Ed,T pS q “

i“1

Nd ÿ

¸1{pz › ›pz vj ›Rdata rx̂sptj q›

(3.54)

j“1

where Rint and Rdata are defined in (3.2). The quantity Er,T measures the discrete physics residual, while Ed,T measures the discrepancy with the observations. Assumption 3 (Conditional stability). Let X̂ Ă X ˚ Ă X “ Lpx pIq be Banach spaces. For any u, v P X̂, 18

assume that the differential operator D and the restriction operator L satisfy ¯ ` ˘´ ρ }u ´ v}Lpx pEq ď Cpd }u}X̂ , }v}X̂ }Dpu, θq ´ Dpv, θ 1 q}Yp ` }Lpuq ´ Lpvq}ρZd ,

(3.55)

for some 0 ă ρp , ρd ď 1 and any subset I 1 Ă E Ă I, where Y “ Lpy pIq and Z “ Lpz pI 1 q. Here Dpx, θqptq “

` ˘ dx ptq ´ f t, xptq; θ , dt

ˇ Lpxq “ xˇI 1 .

(3.56)

r Let S r “ tsi uN i“1 Ă I be quadrature points with weights wi P R` . For a function h : I Ñ R, define

the quadrature functional QN r phq “

Nr ÿ

wi hpsi q.

(3.57)

i“1

Assume that the corresponding quadrature error satisfies ˇ ˇż ˇ ˇ ` ˘ ´α ˇ ˇ hptq dt ´ QN r phqˇ ď Cq }h}Y,d¯ Nr , ˇ

for some α ą 0,

(3.58)

where d¯ denotes the regularity index required by the quadrature rule. 1 d Similarly, let S d “ ttj uN j“1 Ă I be quadrature points for the data term with weights vj P R` . For a

function hd : I 1 Ñ R, define QN d phd q “

Nd ÿ

vj hd ptj q.

(3.59)

j“1

Assume that the corresponding quadrature error satisfies ˇ ˇż ˇ ˇ ` ˘ ´αd ˇ ˇ hd ptq dt ´ QN , d phd qˇ ď Cqd }hd }Z,d¯ Nd ˇ 1

for some αd ą 0.

(3.60)

Then the following theorem bounds the generalization error in terms of the discrete training residuals and the quadrature errors, in the spirit of the conditional-stability-based PINN analysis in [50]. Theorem 3.2 (PINN generalization error estimate for ODE inverse problems). Let x P X̂ Ă X ˚ Ă X be the solution of the inverse problem associated with (2.2), and assume that the stability estimate (3.55) holds for any I 1 Ă E Ă I. Let x̂ P X̂ be a PINN approximation generated from the training sets S r and S d . Assume further that }Rint rx̂, θ̂sp¨q}py P Y and }Rdata rx̂sp¨q}pz P Z, and that the quadrature errors satisfy (3.58) and (3.60). Then ´ ¯ ´α ρ {p EG pEq ď Cpd Er,T pθ̂, S r qρp ` Ed,T pS d qρd ` Cq Nr´αρp {py ` Cqd Nd d d z ,

(3.61)

where ` ˘ Cpd “ Cpd }x}X̂ , }x̂}X̂ ,

´› › ¯ Cq “ Cq ›}Rint rx̂, θ̂sp¨q}py ›Y , 19

´› › ¯ Cqd “ Cqd ›}Rdata rx̂sp¨q}pz ›Z .

(3.62)

Proof. For notational simplicity, define R “ Dpx̂, θ̂q,

ˇ Rd “ Lpx̂q ´ xˇI 1 .

(3.63)

Since px, θq solves (2.2), we have Dpx, θq ” 0, and therefore R “ Dpx̂, θ̂q ´ Dpx, θq.

(3.64)

Similarly, because Lpxq “ x|I 1 , it follows that Rd “ Lpx̂q ´ Lpxq.

(3.65)

Applying the conditional stability estimate (3.55) with u “ x̂ and v “ x gives ¯ ´ ρ EG pEq “ }x̂ ´ x}Lpx pEq ď Cpd }Dpx̂, θ̂q ´ Dpx, θq}Yp ` }Lpx̂q ´ Lpxq}ρZd ´ ¯ ρ “ Cpd }R}Yp ` }Rd }ρZd .

(3.66)

We next bound the two residual norms by their discrete training counterparts. Since Y “ Lpy pIq, applying the quadrature estimate (3.58) to hptq “ }Rptq}py yields ż

}R}Yy “

}Rptq}py dt ď I

Nr ÿ

wi }Rpsi q}py ` Cq Nr´α “ Er,T pθ̂, S r qpy ` Cq Nr´α .

(3.67)

i“1

Since 0 ă ρp ď 1, we use the subadditivity of a ÞÑ aρp {py on R` to obtain ρ

}R}Yp ď Er,T pθ̂, S r qρp ` Cqρp {py Nr´αρp {py .

(3.68)

Similarly, since Z “ Lpz pI 1 q, applying (3.60) to hd ptq “ }Rd ptq}pz gives ż }Rd }pZz “

}Rd ptq}pz dt ď I1

Nd ÿ

vj }Rd ptj q}pz ` Cqd Nd´αd “ Ed,T pS d qpz ` Cqd Nd´αd .

(3.69)

j“1

Again, since 0 ă ρd ď 1, ρ {pz

}Rd }ρZd ď Ed,T pS d qρd ` Cqdd

´αd ρd {pz

(3.70)

Substituting (3.68) and (3.70) into (3.66) yields (3.61). This bound also suggests a natural notion of a well-trained PINN, namely the regime in which the discrete training residuals are of the same order as, or smaller than, the quadrature-induced generalization

gap as ! ) ρ {p ´α ρ {p max Er,T pθ̂, S r qρp , Ed,T pS d qρd ď Cqρp {py Nr´αρp {py ` Cqdd z Nd d d z .

(3.71)

In this regime, the dominant contribution to the generalization error comes from discretization and quadrature rather than incomplete optimization. Remark. Estimate (3.61) decomposes the generalization error into contributions from training residuals and quadrature approximation under a conditional-stability framework. In particular, it links the PINN error explicitly to both optimization accuracy and sampling complexity. In the context of RAA-PINNs, this result provides theoretical support for the reliability of the local physics-informed inverse modules used in residual-based screening and local refinement. We further show that a parameter jump induces a non-vanishing post-change physics residual when the post-change dynamics are evaluated with a mismatched constant parameter in following theorem. Theorem 3.3 (Post-change residual lower bound under parameter mismatch). Suppose that the true r r parameter satisfies (3.1). Let E` Ă I X rτ, T s, S` “ S r X E` , and Nr,` “ |S` |. For any θ̃ P Θ, assume

that there exists γ` ą 0 such that › › › › ›f p¨, xp¨q; θ ` q ´ f p¨, xp¨q; θ̃q› py L

pE` q

ě γ` }θ ` ´ θ̃}.

(3.72)

Define 9̂ ´ x} 9 Lpy pE` q ` L}x̂ ´ x}Lpy pE` q . ε` “ }x

(3.73)

› › › › ›Dpx̂, θ̃q› py

(3.74)

Then L

pE` q

´ ¯ ě γ` }θ ` ´ θ̃} ´ ε` . `

Moreover, if (3.58) holds on E` with constant Cq,` , then r qě Er,T pθ̃, S`

„´

γ` }θ ` ´ θ̃} ´ ε`

¯py `

ȷ1{py ´α ´ Cq,` Nr,`

(3.75)

In particular, for θ̃ “ θ ´ , the post-change residual is bounded from below by the jump size }θ ` ´ θ ´ } up to approximation and quadrature errors. Proof. On E` , the true trajectory satisfies 9 xptq “ f pt, xptq; θ ` q.

(3.76)

Hence, for any θ̃ P Θ, 9̂ ´ f pt, x̂; θ̃q Dpx̂, θ̃q “ x 9̂ ´ x9 ` f pt, x; θ ` q ´ f pt, x; θ̃q ` f pt, x; θ̃q ´ f pt, x̂; θ̃q. “x

(3.77)

Using the reverse triangle inequality, the Lipschitz condition (2.6), and (3.72), we obtain › › › › ›Dpx̂, θ̃q› py L

pE` q

› › › › ě ›f p¨, x; θ ` q ´ f p¨, x; θ̃q› p

L y pE` q

9̂ ´ x} 9 Lpy pE` q ´ L}x̂ ´ x}Lpy pE` q ´ }x

ě γ` }θ ` ´ θ̃} ´ ε` .

(3.78)

Taking the positive part gives (3.74). It remains to pass from the continuous residual to the discrete training residual. Applying the quadrature estimate on E` to hptq “ }Dpx̂, θ̃qptq}py gives › ›py › › r py Er,T pθ̃, S` q ě ›Dpx̂, θ̃q› p

L y pE` q

´α ´ Cq,` Nr,` .

(3.79)

Combining this inequality with (3.74) yields (3.75). Remark. Theorem 3.3 states that, after a parameter jump, using a pre-change or mismatched parameter necessarily generates a positive post-change physics residual whenever the jump magnitude dominates approximation and quadrature errors. This provides a theoretical justification for locating change points by detecting elevated physics residual losses after candidate transition times. In summary, this section establishes the RAA-PINNs framework for change-point detection and parameter recovery in nonlinear dynamical systems with regime switching. The proposed method combines local physics-informed inverse learning, residual-loss-anomaly-based coarse localization, and differentiable local refinement within a unified optimization framework. The theoretical analysis further shows that PINN generalization error can be controlled by training residuals and quadrature errors, while parameter jumps induce non-vanishing physics residuals under mismatched regimes. These results provide the mathematical foundation for using residual loss as an intrinsic signal of change points. In the next section, we validate the proposed framework on representative nonlinear dynamical systems and examine its performance in both change-point localization and parameter estimation.

4. Numerical Experiment In this section, we evaluate the performance of the proposed RAA-PINNs method for change-point detection and parameter estimation in nonlinear dynamical systems. We consider five representative

systems: the Malthus model, the logistic model, the Van der Pol oscillator, the Lotka-Volterra model, and the Lorenz system. For all experiments, the time domain is divided into overlapping subintervals, and local PINNs are trained independently in each subinterval. Each network is a fully connected feedforward network, taking time as input and system state variables as output. The baseline network has 4 hidden layers with 64 neurons per layer. In the refinement stage, the width is increased to 80. All networks use the tanh activation function. The loss function includes both a data fitting term and a physics-informed residual term with equal weights. Training is performed using the Adam optimizer, with learning rates of 5 ˆ 10´3 for network parameters and 10´3 for physical parameters. In the refinement stage, a three-step training procedure is used: network pre-training, optimization of change-point positions, and joint optimization of parameters and change points. Change points are parameterized using a sigmoid function and linked with smooth transitions to ensure differentiability and numerical stability during optimization.

4.1. Malthus and Logistic Model Population growth is a classical topic in dynamical systems, describing the temporal evolution of biological populations. Such models are widely used in ecology, demography, and resource management. Two representative models are Malthus and logistic models. In this work, we employ the RAA-PINNs to jointly estimate time-varying parameters and detect change points, illustrating the application of our framework to NDS-CPD. The Malthusian model describes exponential population growth assuming unlimited resources. It is expressed as a one-dimensional ordinary differential equation: dPm “ rm Pm , dt

(4.1)

where Pm denotes the population at time t, and rm is the population growth rate. In the context of a single-regime ODE, we denote xptq “ Pm ptq and θ “ rm following the traditional nonmixture model formulation (2.1). We extend this model to a time-varying setting with multiple change points, using RAA-PINNs to jointly estimate the piecewise parameter rm ptq and detect change points. To illustrate, we introduce a discrete state variable rptq to represent the regime switching over the temporal domain r0, T s with T “ 100. Let I “ t1, 2u, and define

rptq “

$ ’ &1, 0 ď t ă 40, ’ %2, 40 ď t ď 100.

(4.2)

#103

10 2 r(t)

Pm (t)

8 6 4 1 2 0

100

#10!4

0.15 PDE loss Search interval

Predicted solution Reference solution

True change-point Predicted change-point

0.1 rm

PDE loss

60 t

0.05

0 0

100

Figure 2: Malthus model. Top row: the output of Pm ptq over r0, 100s (left) and the corresponding sample path (right). Bottom row: Stage I identifies the candidate change-point interval determined from the physical loss in the overlapping domain (left), and Stage II results for refined parameter estimation and change-point localization (right).

Accordingly, the temporal evolution of rptq follows rptq : 1 Ñ 2. The piecewise population growth rate is $ ’ &0.1,

rptq “ 1,

’ %0.05,

rptq “ 2.

θptq “

(4.3)

In Stage I, the temporal domain is divided into overlapping subintervals for coarse change-point localization. The candidate change-point intervals identified by the elevated physical loss are r39, 41s, which fully cover the true change point at t “ 40. Stage II refines this interval to jointly estimate the precise change-point location and the piecewise parameter values. The system output, sample path, and residual results are shown in Figure 2. The logistic model incorporates environmental carrying capacity to capture the transition from rapid growth to a stable population. It is described by ˆ ˙ dPl Pl “ rl Pl 1 ´ , dt Q

(4.4)

where Pl denotes the population at time t, rl is the growth rate, and Q is the carrying capacity. In a

#102 2

3 r(t)

Pl (t)

2 1

100

t 12

#10!4

100

Predicted solution Reference solution

True change-point Predicted change-point

0.1

PDE loss

0.15

PDE loss Search interval

60 t

0.05

2 0

0 0

100

Figure 3: Logistic model. Top row: the output of Pl ptq over r0, 100s (left) and the corresponding sample path (right). Bottom row: Stage I candidate change-point interval (left) and Stage II refined parameter estimation and change-point localization (right).

single-regime ODE, xptq “ Pl ptq and θ “ rl . We consider a time-varying logistic equation with a change point. Introducing a discrete state variable rptq over r0, T s with T “ 100: rptq “

$ ’ &1, 0 ď t ă 60,

(4.5)

’ %2, 60 ď t ď 100. The piecewise parameter vector is $ ’ &0.1,

rptq “ 1,

’ %0.05,

rptq “ 2.

θptq “

(4.6)

In Stage I, the candidate change-point interval identified by the elevated physical loss is r59, 61s, covering the true change point at t “ 60. Stage II jointly refines this interval to estimate the precise change-point location and parameter values. The system output, sample path, and residuals are shown in Figure 3. For both models, the time interval is r0, 100s with window length 2 and step size 1. All observations within each overlapping domain are used to estimate parameters independently with RAA-PINNs. El-

Table 1: Parameter estimation results for different numerical examples.

Numerical Example

Equation Coefficient

Time

True Value

Parameter Estimation

Squared Error

Malthus

r0, 40s r40, 100s

0.1 0.05

0.0972 0.0487

7.840 ˆ 10´6 1.690 ˆ 10´6

Logistic

r0, 60s r60, 100s

0.1 0.05

0.0947 0.0435

2.809 ˆ 10´5 4.225 ˆ 10´5

r0, 40s r40, 80s r80, 100s

1 0.1 0.5

1.0238 0.0912 0.4878

5.664 ˆ 10´4 7.744 ˆ 10´5 1.488 ˆ 10´4

Van der Pol

evated physical loss in domains containing change points determines Stage I candidate intervals, which Stage II refines to obtain precise change-point locations and parameter estimates. Detailed statistical inference and mean square errors are reported in Tables 1 and 4. This example demonstrates that RAA-PINNs can recover time-varying parameters in simple exponential growth systems, illustrating the method’s effectiveness in identifying abrupt transitions in linear population dynamics.

4.2. Van der Pol Oscillator The Van der Pol oscillator is a canonical self-excited nonlinear oscillatory system, originally developed to model electronic circuits. Due to its simple structure yet rich nonlinear behavior, it has become a standard benchmark for studying nonlinear oscillations and stable limit-cycle dynamics, and is widely used in circuits, mechanical systems, and biological rhythms. In this work, we apply RAA-PINNs to estimate time-varying parameters and detect change points, illustrating its application to NDS-CPD. The oscillator is governed by the second-order differential equation: ` ˘ dv d2 v ´ µ 1 ´ v2 ` v “ 0. 2 dt dt

(4.7)

This equation can be equivalently expressed as a first-order system: $ dM ’ & “ N, dt ’ % dN “ µp1 ´ M 2 qN ´ M, dt

(4.8)

where xptq “ pM ptq, N ptqqJ and θptq “ µptq in the traditional nonmixture ODE model (2.1). We consider a Van der Pol oscillator with a piecewise time-varying parameter µptq exhibiting multiple change points. A discrete state variable rptq with state space I is introduced to characterize the parameter

4 3

2 r(t)

(M (t); N (t))>

M(t) N(t)

-2 1 -4 0

100

t 20

#10!3

100

PDE loss Search interval

Predicted solution Reference solution

1.5

True change-point Predicted change-point

1 10

PDE loss

60 t

0.5 5

-0.5 0

100

Figure 4: Van der Pol oscillator model. Top row: the output of pM ptq, N ptqq over r0, 100s (left) and the corresponding sample path (right). Bottom row: Stage I results show candidate change-point intervals determined from the physical loss in the overlapping domain (left), and Stage II results for parameter estimation and change-point localization (right).

switching. Let I “ t1, 2, 3u and T “ 100, and define $ ’ ’ 1, 0 ď t ă 40, ’ ’ & rptq “ 2, 40 ď t ă 80, ’ ’ ’ ’ %3, 80 ď t ď 100.

(4.9)

Accordingly, the temporal evolution of rptq is 1 Ñ 2 Ñ 3, and the piecewise parameter θptq is $ ’ ’ 1, ’ ’ & θptq “ 0.1, ’ ’ ’ ’ %0.5,

rptq “ 1, rptq “ 2,

(4.10)

rptq “ 3.

The system output and sample paths are shown in Figure 4. In this model, the time interval is set to r0, 100s with a window length of 2 and a step size of 1. In Stage I, each window is independently trained for 30,000 iterations to estimate the local constant

parameter µ, and candidate transition intervals are identified based on the physical residual. Specifically, the Stage I candidate change-point intervals are r39, 41s and r79, 81s, fully covering the true change points. In Stage II, these candidate intervals are refined to jointly optimize the transition locations and the piecewise parameters, demonstrating that RAA-PINNs can robustly recover parameter variations and accurately identify change points in nonlinear dynamical systems. Within each overlapping domain, all available observation data are used to estimate parameters independently. Elevated physical loss in domains containing change points is averaged over the last 100 training iterations to determine the Stage I candidate intervals. The Stage II optimization then produces the final parameter estimates and refined change-point locations. Detailed prediction errors, statistical inference results, and mean square errors for parameter recovery and change-point detection are reported in Tables 1 and 4. This case illustrates that RAA-PINNs can recover parameters in self-excited nonlinear oscillatory systems, demonstrating the framework’s robustness in identifying transitions in oscillatory dynamics with pronounced nonlinearity. 4.3. Lotka-Volterra Model The Lotka-Volterra model, commonly known as the predator-prey model, describes species interactions in an ecological system. Its key assumptions include closed population dynamics, a closed system with no migration, and a linear functional response in predator-prey interactions. The dynamics of a prey species S and a predator species W are governed by the following two-dimensional coupled differential equations: $ dS ’ & “ Spα ´ βW q, dt ’ % dW “ ´W pγ ´ δSq, dt

(4.11)

where xptq “ pSptq, W ptqqJ and θ “ pα, β, γ, δq in the traditional nonmixture ODE model (2.1). In this work, we employ the RAA-PINNs to jointly estimate the time-varying parameters and detect change points, illustrating the application to NDS-CPD. We consider a Lotka-Volterra system with multiple time-varying parameters and multiple change points. A discrete state variable rptq with state space I is introduced to characterize the temporal variation of system parameters. Let I “ t1, 2, 3u and T “ 100, and define $ ’ ’1, ’ ’ & rptq “ 2, ’ ’ ’ ’ % 3,

0 ď t ă 20, 80 ď t ď 100, 20 ď t ă 40, 60 ď t ă 80,

(4.12)

40 ď t ă 60.

Thus, the temporal evolution of rptq is 1 Ñ 2 Ñ 3 Ñ 2 Ñ 1, and the corresponding piecewise parameters 28

Table 2: Parameter estimation for the Lotka-Volterra model with time-varying parameters.

Time

Equation Coefficient

True Value

Parameter Estimation

Squared Error of Parameter

r0, 20s

α β γ δ

2 1 2 1

2.0238 0.9863 1.9879 0.9878

5.664 ˆ 10´4 1.877 ˆ 10´4 1.464 ˆ 10´4 1.488 ˆ 10´4

r20, 40s

α β γ δ

4 2 3 4

3.9912 1.9865 3.0365 3.9876

7.774 ˆ 10´5 1.823 ˆ 10´4 1.332 ˆ 10´3 1.538 ˆ 10´4

r40, 60s

α β γ δ

3 4 1 2

2.9878 3.9912 0.9873 1.9684

1.488 ˆ 10´4 7.744 ˆ 10´5 1.613 ˆ 10´4 9.985 ˆ 10´4

r60, 80s

α β γ δ

4 2 3 4

3.9865 1.9932 2.9834 4.0534

1.823 ˆ 10´4 4.624 ˆ 10´5 2.756 ˆ 10´4 2.852 ˆ 10´3

r80, 100s

α β γ δ

2 1 2 1

1.9845 0.9941 2.0523 1.0376

2.403 ˆ 10´4 3.481 ˆ 10´5 2.735 ˆ 10´3 1.414 ˆ 10´3

are $ ’ ’ p2, 1, 2, 1q, ’ ’ & θptq “ p4, 2, 3, 4q, ’ ’ ’ ’ %p3, 4, 1, 2q,

rptq “ 1, rptq “ 2,

(4.13)

rptq “ 3.

The system output and sample paths are shown in Figure 5. In the Lotka-Volterra model, the time interval is r0, 100s with a window length of 2 and a step size of 1. This experiment validates the effectiveness and scalability of RAA-PINNs for multi-parameter coupled nonlinear systems. In each overlapping domain, all observation data are used to estimate the parameters pα, β, γ, δq independently. Elevated physical loss in domains containing change points is averaged over the last 100 training iterations to determine the Stage I candidate intervals, which are r19, 21s, r39, 41s, r59, 61s, and r79, 81s. In Stage II, the candidate intervals are refined to jointly optimize the parameter values and change-point locations. The time-varying parameter estimates obtained using overlapping-domain RAA-PINNs and the corresponding change-point detection results are shown in Figure 5. Prediction errors, statistical inference results, and mean square errors for parameter recovery and change-point detection are reported in Tables 2 and 4. Overall, the results of overlapping-domain RAA-PINNs agree well with the reference solution, successfully capturing all four change points. This demonstrates that RAA-PINNs can robustly handle multi-parameter coupled nonlinear systems, recovering simultaneous changes in interacting species.

S(t) W(t)

(S(t); W (t))>

20 15 10 5 0 0

50 t

100

#10!3

PDE loss Search interval

4 PDE loss

r(t)

3 2 1

0 0

100

60 t

6 Predicted solution Reference solution

True change-point Predicted change-point

5 4

4 -

Predicted solution Reference solution

3 2 1

0 0

100

5 Predicted solution Reference solution

True change-point Predicted change-point

Predicted solution Reference solution

6 5

4 /

True change-point Predicted change-point

3 2 2 1

1 0

100

Figure 5: Lotka-Volterra model. First row: trajectories of pSptq, W ptqq over r0, 100s. Second row: sample path (left) and Stage I candidate change-point intervals determined from the physical loss in the overlapping domain (right). Third and fourth rows: Stage II results for parameter estimation and change-point localization.

4.4. Lorenz System The Lorenz system is a canonical model in nonlinear dynamics, widely used to study chaotic behavior and complex dynamical phenomena. It was originally derived from a truncated Fourier expansion of the Navier-Stokes and heat equations describing fluid convection, and has played a fundamental role in revealing sensitive dependence on initial conditions and long-term unpredictability in deterministic systems. In this work, we employ the RAA-PINNs to jointly estimate time-varying parameters and detect change points, illustrating its application to NDS-CPD. The system is described by the following three-dimensional coupled differential equations for the convective intensity U , the horizontal temperature gradient V and the vertical temperature W : $ dU ’ ’ “ σpV ´ U q, ’ ’ ’ & dt dV “ rU ´ V ´ U W, ’ dt ’ ’ ’ ’ % dW “ U V ´ bW, dt

(4.14)

where xptq “ pU ptq, V ptq, W ptqqJ and θ “ pσ, r, bq in the traditional nonmixture ODE model (2.1). We consider a time-varying Lorenz system with multiple change points. A discrete state variable rptq with state space I is introduced to characterize the temporal evolution of system parameters. Let I “ t1, 2, 3u and T “ 20, and define $ ’ ’ 1, 0 ď t ă 8, ’ ’ & rptq “ 2, 8 ď t ă 16, ’ ’ ’ ’ %3, 16 ď t ď 20.

(4.15)

Accordingly, the temporal evolution of rptq is 1 Ñ 2 Ñ 3, and the piecewise parameter vector is $ ’ ’ r10, 28, 8{3s, rptq “ 1, ’ ’ & θptq “ r12, 16, 10{3s, rptq “ 2, ’ ’ ’ ’ % r14, 20, 3s, rptq “ 3.

(4.16)

The system output and sample paths are shown in Figure 6. Given the sensitive dependence on initial conditions inherent in chaotic systems, trajectories diverge rapidly over longer time spans. Consequently, the simulation interval is restricted to r0, 20s to maintain numerical stability and allow RAA-PINNs to accurately recover the time-varying parameters and detect change points. The time interval is divided into 100 equal-length subintervals, with a window length of 0.2 and a step size of 0.1. In Stage I, each window is independently trained for 30,000 iterations to 31

r(t)

1 0

0.1

10 t

16 PDE loss Search interval

0.08

Predicted solution Reference solution

True change-point Predicted change-point

0.06

13 <

PDE loss

0.04

11 0.02

9 0

10 t

5 Predicted solution Reference solution

True change-point Predicted change-point

Predicted solution Reference solution

4.5

True change-point Predicted change-point

30 25

4 3.5 3

2.5 15

2 0

10 t

Figure 6: Lorenz model. Top row: trajectories of pU ptq, V ptq, W ptqq over r0, 20s (left) and the corresponding sample path (right). Second row: Stage I candidate change-point intervals determined from the physical loss in the overlapping domain (left), and Stage II results for parameter estimation and change-point localization (right). Other rows: Stage II parameter estimation results.

Table 3: Parameter estimation and change-point detection of Lorenz system with time-varying parameters.

Time

Equation Coefficient

True Value

Parameter Estimation

Squared Error of Parameter

r0, 8s

σ r b

10 28 8/3

10.0238 27.9863 2.6877

5.664 ˆ 10´4 1.877 ˆ 10´4 4.424 ˆ 10´4

r8, 16s

σ r b

12 16 10/3

11.9912 15.9865 3.3787

7.774 ˆ 10´5 1.823 ˆ 10´4 2.058 ˆ 10´3

r16, 20s

σ r b

14 20 3

13.9878 19.9912 2.9873

1.488 ˆ 10´4 7.724 ˆ 10´5 1.613 ˆ 10´4

Table 4: Estimates of change points for different numerical examples with time-varying parameters.

Numerical Example

Change point interval

True Value

Change point Estimation

Squared Error

MSE of xptq

Malthus

r39, 41s

40.0093

8.649 ˆ 10´5

3.854 ˆ 10´4

Logistic

r59, 61s

60.0037

1.369 ˆ 10´5

7.297 ˆ 10´4

Van der Pol

r39, 41s r79, 81s

40 80

40.0013 79.9987

1.690 ˆ 10´6 1.690 ˆ 10´6

4.372 ˆ 10´4 8.633 ˆ 10´5

Lotka-Volterra

r19, 21s r39, 41s r59, 61s r79, 81s

20 40 60 80

20.0263 40.0139 60.0237 79.9897

6.917 ˆ 10´4 1.932 ˆ 10´4 5.617 ˆ 10´4 1.061 ˆ 10´4

2.538 ˆ 10´4 7.348 ˆ 10´4 3.283 ˆ 10´4 8.240 ˆ 10´5

Lorenz

r7.9, 8.1s r15.9, 16.1s

8 16

8.0237 16.0128

5.617 ˆ 10´4 1.638 ˆ 10´4

3.184 ˆ 10´4 7.216 ˆ 10´4

estimate local constant parameters, and candidate change-point intervals are determined from elevated physical residuals. The Stage I candidate intervals are r7.9, 8.1s and r15.9, 16.1s, covering the true change points. In Stage II, these intervals are refined to jointly optimize parameter values and change-point locations.

5. Parallel Overlapping Domain Strategy Within each overlapping domain, all available observation data are used to estimate σ, r, and b independently. The final time-varying parameter estimates and change-point locations obtained via RAA-PINNs are shown in Figure 6. Prediction errors, statistical inference results, and mean square errors for parameter recovery and change-point detection are reported in Table 3 and Table 4. This demonstrates that RAA-PINNs can robustly recover parameters and accurately detect change points even in highly sensitive chaotic systems. Overall, this example demonstrates that the framework can robustly recover parameters and accurately identify transitions even in highly sensitive chaotic systems, highlighting its applicability to complex chaotic dynamics. 33

For change-point detection in nonlinear dynamical systems with regime switching, the proposed RAAPINNs framework combines residual-loss anomaly analysis with overlapping-domain training to achieve stable coarse localization and accurate local parameter inversion. However, the computational cost of overlapping-domain decomposition can become a major bottleneck in practical implementations, especially when the temporal resolution is high or the observation horizon is long. This issue is most pronounced in Stage I, where independent PINNs must be trained on many overlapping subintervals in order to obtain reliable residual-loss statistics and identify candidate change-point regions. Since the number of subintervals grows approximately linearly with the temporal resolution and the coverage of the time domain, a serial implementation leads to an almost linear increase in total training time, which substantially limits the scalability of the method. A key observation is that, in Stage I, there are no interface coupling constraints between overlapping subintervals, and the corresponding training tasks are fully independent. Therefore, this stage has a natural parallel structure and can be directly mapped to multi-core CPUs, multi-GPU platforms, or distributed computing environments, leading to a substantial reduction in wall-clock time. Let the temporal domain be I “ r0, T s. In Stage I, we construct a coarse partition 0 “ t0 ă t1 ă ¨ ¨ ¨ ă tK “ T,

(5.1)

and define K overlapping subintervals with overlap width δ ą 0 as Ik “ rtk´1 ´ δ, tk ` δs X r0, T s,

k “ 1, . . . , K.

(5.2)

On each Ik , we train an independent PINN xϕk : Ik Ñ Rn together with a local constant parameter θ k P Θ by minimizing the subinterval loss Jk pϕk , θ k q “

ÿ › ›2 › ›p d › r › wi,k Rint rxϕk , θ k spsi q›p , wj,k Rdata rxϕk , θ k sptj q›2 ` λ

pkq

where Sd

(5.3)

pkq

si PSint

tj PSd

pkq

Ă Ik denotes the measurement set, Sint Ă Ik denotes the interior collocation set, and λ ą 0

balances the data mismatch and physics residual. Crucially, the optimization problems tJk uK k“1 are mutually independent.

The absence of cross-

subinterval coupling enables direct parallelization at the subinterval level, with each pair pϕk , θ k q optimized on a separate computational unit. Let tsub denote the average training time for one PINN on a single Stage I subinterval, and let K be pIq

the number of overlapping subintervals. Then the total runtime of serial Stage I, denoted by Ts , can be

Table 5: Runtime comparison of the proposed framework under serial and parallel implementations. Here P denotes the number of parallel workers, S “ Ts {Tp is the speedup, and E “ S{P is the parallel efficiency.

Example

Stage

Ts psq

Tp psq

Van der Pol

I II Total

100 2

128 128 128

53625.7 167.9 53793.6

425.1 92.6 517.7

126.152 1.811 103.915

0.986 0.014 0.812

Lotka-Volterra

I II Total

100 4

128 128 128

55554.8 401.2 55956.0

594.5 108.2 702.7

93.465 3.771 79.525

0.730 0.029 0.621

Lorenz

I II Total

100 2

128 128 128

72812.3 252.6 73064.9

718.0 119.5 837.5

101.339 2.141 87.288

0.792 0.016 0.682

approximated as TspIq “ K tsub ` c1 ,

(5.4)

where c1 collects additional overhead that is typically much smaller than K tsub and therefore does not affect the overall scaling behavior. Given load-balancing effects and communication overhead, the runtime pIq

with P parallel workers can be approximated by R TppIq “

V K tsub ` c2 , P

(5.5)

where c2 denotes parallelization overhead and is typically much smaller than rK{P s tsub . Hence, the theoretical speedup is given by S pIq “ TspIq {TppIq .

(5.6)

When K " P , we obtain S pIq « P , indicating near-linear scaling. Moreover, Stage II is carried out only on a narrow candidate interval and typically requires training one or only a few local PINNs. As a result, the total runtime of the two-stage RAA-PINNs framework is dominated by Stage I. The total runtime under serial and parallel implementations can therefore be written as Ts “ TspIq ` T pIIq ,

Tp “ TppIq ` T pIIq .

(5.7)

In our experiments, each GPU with 24 GB of memory can execute 16 tasks concurrently. Using 8 GPUs yields a total of 8 ˆ 16 “ 128 concurrent tasks. We evaluate the computational benefit of this parallel strategy on three representative nonlinear dynamical systems with discontinuous parameters: the Van der Pol oscillator, the Lotka-Volterra model, and the Lorenz system. These examples cover oscillatory, coupled population, and chaotic dynamics, respectively, and therefore provide a representative testbed for assessing the scalability of RAA-PINNs. The results are presented in Table 5. Notably, the parallel strategy changes only the scheduling of independent subinterval optimization

tasks. It does not alter the mathematical form of the subinterval losses or the Stage II refinement procedure. Therefore, the parallel and serial implementations produce consistent change-point localization and parameter reconstruction results, while the wall-clock time is substantially reduced. This confirms that the proposed parallel overlapping-domain strategy improves the computational scalability of RAA-PINNs without sacrificing inversion accuracy.

6. Comparative Experiment In this section, we compare the proposed method with both traditional statistical baselines and existing neural-network-based approaches. The goal is to evaluate its performance on two core tasks in nonlinear dynamical systems with regime switching: parameter estimation and change-point localization. The comparative study is carried out on representative systems to highlight the differences in modeling strategy, reconstruction accuracy, and localization precision.

6.1. Traditional Statistical Method For ordinary differential equation systems with discontinuous parameters, we adopt a two-stage statistical pipeline as the baseline for comparison. Specifically, the baseline combines adaptive gradient matching (AGM)[51] for ODE parameter inference with Pruned Exact Linear Time (PELT)[52] for penalized change-point detection, hereafter referred to as AGM-PELT. Under this baseline, parameter estimation is carried out within a probabilistic inference framework using AGM. In this approach, Gaussian processes are used to model the latent system trajectory, while the governing ODE is incorporated as a soft constraint rather than being enforced exactly. The resulting joint inference problem can be formulated as px̂, θ̂q “ arg max ppy | xq ppx9 | f px; θqq , x,θ

(6.1)

which yields a time-dependent parameter sequence tθ t uTt“1 . Within this statistical formulation, the parameters are treated as stochastic quantities, and the governing equations are incorporated through probabilistic consistency between the inferred trajectory and the ODE dynamics. Given the estimated parameter sequence tθ t u, change-point detection is then formulated as a penalized global segmentation problem: # tτi u “ arg min

m,tτi u

m`1 ÿ

+ ´ C

i tθ t uτt“τ i´1

¯ ` ψm ,

(6.2)

i“1

where Cp¨q denotes a segment-wise cost function and ψ is a penalty parameter controlling the number of detected change points. This optimization is solved using the PELT algorithm, which combines dynamic programming with pruning to obtain an exact solution for the chosen penalized cost. Under suitable

8 RAA-PINNs solution AGM solution Reference solution

True change-point Predicted change-point PELT change-point

RAA-PINNs solution AGM solution Reference solution

True change-point Predicted change-point PELT change-point

5 4

4 3

2 1

0 0

100

8 RAA-PINNs solution AGM solution Reference solution

True change-point Predicted change-point PELT change-point

RAA-PINNs solution AGM solution Reference solution

True change-point Predicted change-point PELT change-point

4 4

3 2

1 0

0 0

100

Figure 7: Learned parameter trajectories for the Lotka-Volterra model obtained by different methods, where the parameters are θ “ pα, β, γ, δq.

conditions, its expected computational cost can scale linearly with the sequence length. Accordingly, the detection performance depends on both the choice of the segment-wise cost function and the specification of the penalty parameter. In contrast, the proposed method embeds parameter inference directly into a physics-informed inverse framework. Instead of first estimating a time-varying parameter trajectory and then segmenting it, the governing dynamics are enforced explicitly over multiple overlapping subintervals, so that parameter estimation remains tightly coupled with the underlying dynamical system. This design reduces the potential mismatch between statistically smoothed trajectories and physically admissible dynamics. The proposed two-stage strategy adopts a different perspective on change-point detection. In Stage I, neural networks trained on overlapping subintervals exhibit clear residual-loss anomalies in regions where parameter discontinuities occur, thereby naturally identifying candidate change-point intervals. Then in Stage II, refined local modeling is carried out only within these intervals to jointly infer the piecewise parameters and the change-point locations. In this way, change points are not obtained through explicit segmentation of a pre-estimated parameter sequence, but instead emerge from local breakdowns of physical consistency. Comparative experiments between this statistical baseline and the proposed method are conducted on the Lotka-Volterra system, and the results are reported in Figure 7 and Table 6.

Overall, the results show that the proposed method provides a more faithful reconstruction of discontinuous parameters than the AGM-PELT baseline. For all four coefficients, the recovered trajectories are closer to the reference piecewise-constant structure, exhibit smaller within-regime dispersion, and display sharper transitions at the discontinuities. In contrast, the AGM-PELT baseline shows noticeable scatter around the true parameter levels, and its representation of regime-wise dynamics is less accurate. This qualitative difference is also supported by the mean squared errors reported in Table 6, where the proposed method consistently achieves lower errors than the baseline.

6.2. Existing Neural Network Method For ordinary differential equation systems with jump parameters, the PINNs combined with statistical learning algorithm named expectation-maximization for Gaussian mixture models (PINNs-EM-GMM) framework is a recently proposed hybrid approach for joint parameter estimation and change-point detection [35]. It combines PINNs-based parameter inversion with statistical mechanism identification. Let the system state be observed as txpti quN i“1 and governed by the ODE as ` ˘ 9 xptq “ f t, xptq; θptq ,

θptq P Rd ,

(6.3)

where θptq represents the time-varying parameter vector. The PINNs-EM-GMM framework consists of two stages. In the first stage, sliding-window PINNs are employed to reconstruct the state function xϕ ptq and the parameter function θ ξ ptq locally over each time window rts , te s by minimizing the composite loss J “

}xϕ pti q ´ xpti q}22 ` λ

ti PTobs

ÿ › ` ˘› ›x9 ϕ ptj q ´ f tj , xϕ ptj q; θ ξ ptj q ›2 , 2

(6.4)

tj PTres

where the first term penalizes data mismatch and the second enforces physical consistency, with λ ą 0 controlling the weight of the residual. By leveraging overlapping windows and adaptive residual weighting, the network can capture rapid changes near transition points while mitigating cross-regime compromise, p yielding a smooth continuous parameter trajectory θptq. In the second stage, mechanism discretization and change-point detection are performed. The continuous trajectory is treated as a sample sequence generated from a finite set of latent regimes. For each parameter component θpi ptq, a K-component Gaussian mixture model is employed to characterize its distribution:

K ÿ ˘ ` ˘ ` 2 ηik N θpi ptq | µik , σik , p θpi ptq “

(6.5)

k“1 2 where ηik , µik , and σik denote the mixture weight, mean, and variance of the k-th cluster, respectively.

The EM algorithm estimates these parameters and computes the posterior cluster probability for each

Reference solution PINNs-EM-GMM solution RAA-PINNs solution

6 5 4 -

3 3 2 2

1 0

100

t Reference solution PINNs-EM-GMM solution RAA-PINNs solution

100

t Reference solution PINNs-EM-GMM solution RAA-PINNs solution

6 5 4 /

3 2 2 1

1 0

100

Change Probability

1 Change probability True change-point Predicted change-point

0.8 0.6 0.4 0.2 0 0

50 t

Figure 8: First and second rows: learned parameter trajectories θ “ pα, β, γ, δq of the Lotka-Volterra model obtained using different PINNs methods. Third row: change-point detection results based on the EM-GMM approach.

time point as ` ˘ γik ptq “ P zi ptq “ k | θpi ptq .

(6.6)

A change-point probability is then constructed via a three-point local consistency measure

pi ptq “ 1 ´

K ÿ

γik pt ´ ∆tq γik ptq γik pt ` ∆tq,

(6.7)

k“1

and the probabilities across all parameters are averaged to obtain a combined change-point probability d

pptq “

1ÿ pi ptq. d i“1

(6.8)

The top Ncp change points are selected from the peaks of pptq, and the corresponding credible intervals are formed as ‰ “ Ij “ min τpij , max τpij . i

(6.9)

This framework integrates continuous parameter inversion, statistical regime discretization, and the posterior-based change-point detection into a unified pipeline. As a result, it provides an interpretable statistical representation of both parameter evolution and mechanism transitions, and can exhibit robustness in noisy settings. The comparative results show that the PINNs-EM-GMM method achieves high accuracy in reconstructing parameter trajectories. As shown in Figure 8, it effectively captures the piecewise evolution of the system parameters effectively. In addition, Table 6 reports the mean squared errors of the parameters pα, β, γ, δq, indicating accurate parameter estimation across all four dimensions. Figure 8 also presents the change-point probability together with the true and predicted change-point locations. Pronounced probability peaks appear near the true change points, demonstrating the method’s ability to identify candidate transition regions reliably. The detected change-point intervals, namely r18.6093, 21.2106s, r38.8694, 41.1029s, r59.5799, 61.2806s, and r78.8894, 81.3907s, all contain the true change-point locations, indicating strong interval-level localization performance. However, because it relies on statistical discretization and local posterior-consistency measures, the PINNs-EM-GMM framework intrinsically produces interval estimates rather than direct point estimates of change-point locations. By contrast, the proposed two-stage method first identifies candidate transition regions through residual-loss anomalies observed in overlapping subintervals, and then introduces a differentiable parameterization of the change point within each candidate region so that the transition location and the piecewise parameters can be optimized jointly under a unified physics-informed objective. Consequently, while maintaining comparable accuracy in parameter estimation, the proposed method further improves change-point localization by refining interval-level detection into direct pointwise estimation.

Table 6: Comparison of squared parameter estimation errors for the Lotka-Volterra model with time-varying parameters.

Method

Time

Squared Error of α

Squared Error of β

Squared Error of γ

Squared Error of δ

RAA-PINNs

r0, 20s r20, 40s r40, 60s r60, 80s r80, 100s

5.664 ˆ 10´4 7.774 ˆ 10´5 1.488 ˆ 10´4 1.823 ˆ 10´4 2.403 ˆ 10´4

1.877 ˆ 10´4 1.823 ˆ 10´4 7.744 ˆ 10´5 4.624 ˆ 10´5 3.481 ˆ 10´5

1.464 ˆ 10´4 1.332 ˆ 10´3 1.613 ˆ 10´4 2.756 ˆ 10´4 2.735 ˆ 10´3

1.488 ˆ 10´4 1.538 ˆ 10´4 9.985 ˆ 10´4 2.852 ˆ 10´3 1.414 ˆ 10´3

AGM-PELT

r0, 20s r20, 40s r40, 60s r60, 80s r80, 100s

7.079 ˆ 10´2 1.042 ˆ 10´1 1.818 ˆ 10´2 2.338 ˆ 10´2 2.706 ˆ 10´1

2.778 ˆ 10´2 2.522 ˆ 10´2 1.014 ˆ 10´3 6.539 ˆ 10´2 5.150 ˆ 10´2

2.039 ˆ 10´3 1.476 ˆ 10´1 2.053 ˆ 10´2 3.252 ˆ 10´2 4.065 ˆ 10´1

1.993 ˆ 10´2 2.288 ˆ 10´2 1.215 ˆ 10´1 3.724 ˆ 10´3 2.013 ˆ 10´1

PINNs-EM-GMM

r0, 20s r20, 40s r40, 60s r60, 80s r80, 100s

4.885 ˆ 10´3 8.808 ˆ 10´3 1.555 ˆ 10´4 1.890 ˆ 10´4 2.215 ˆ 10´4

1.619 ˆ 10´3 1.613 ˆ 10´3 6.627 ˆ 10´4 3.785 ˆ 10´3 2.921 ˆ 10´4

1.205 ˆ 10´2 1.162 ˆ 10´4 1.479 ˆ 10´4 2.875 ˆ 10´4 2.937 ˆ 10´3

1.706 ˆ 10´3 1.343 ˆ 10´3 9.451 ˆ 10´4 2.476 ˆ 10´3 1.380 ˆ 10´3

7. Conclusion and Discussion This work presents a residual-loss anomaly-based inverse framework for change-point detection and parameter identification in nonlinear dynamical systems with regime switching. The core idea is to analyze the anomalous behavior of the physics-based residual loss, which is highly sensitive to local inconsistencies caused by parameter discontinuities. Exploiting this mechanism allows the identification of candidate change-point intervals and the recovery of piecewise system parameters. The method proceeds in two stages: in the first stage, overlapping subintervals are analyzed, where parameter changes induce localized violations of physical consistency, enabling effective localization of candidate intervals. In the second stage, these intervals are refined to jointly infer the change-point locations and the corresponding parameters. Systematic experiments on representative systems, including simple growth models, nonlinear oscillators, coupled multi-parameter systems, and chaotic systems, demonstrate that this approach accurately identifies change points, reconstructs parameters, and captures latent state transitions across a wide range of nonlinear dynamics. Each example highlights a specific capability, abrupt transitions in simple growth, nonlinear saturation effects, self-excited oscillatory dynamics, coupled multi-parameter interactions, and robustness under chaos with strong sensitivity to initial conditions. Compared with traditional statistical approaches, such as Gaussian-process-based methods, this framework provides greater modeling flexibility and physical consistency for strongly nonlinear systems with pronounced parameter discontinuities. It maintains robust performance under limited observational data, accurately recovers parameters, and reliably characterizes states near change points, without requiring prior knowledge of transition times. While supported by theoretical analysis and numerical evidence, further studies of its generalization under noisy observations, model misspecification, and uncertainty

quantification would enhance its reliability. Future extensions may include high-dimensional coupled systems, stochastic dynamics, and partial differential equation models with regime switching. Potential applications span practical state monitoring and transition analysis in real-world systems, such as modeling the growth dynamics of hemangiomas or analyzing phase-wise vascular aging from time-series observations derived via computer vision, highlighting the framework’s broad applicability in nonlinear dynamical system identification and change-point analysis.

References [1] J. Bian, T. Huang, X. Zhang, C. Wang, Y. Zhang, C. Zeng, Fluctuation-variable correlation as early warning signals of non-equilibrium critical transitions, Physica A: Statistical Mechanics and Its Applications 661 (2025) 130401. [2] K. Evers, D. Borsboom, E. I. Fried, F. Hasselman, L. Waldorp, Early warning signals of complex critical transitions in deterministic dynamics, Nonlinear Dynamics 112 (21) (2024) 19071–19094. [3] N. Masuda, K. Aihara, N. G. MacLaren, Anticipating regime shifts by mixing early warning signals from different nodes, Nature Communications 15 (1) (2024) 1086. [4] A. Dmitriev, V. Kornilov, V. Dmitriev, N. Abbas, Early warning signals for critical transitions in sandpile cellular automata, Frontiers in Physics 10 (2022) 839383. [5] V. Dakos, C. A. Boulton, J. E. Buxton, J. F. Abrams, B. Arellano-Nava, D. I. Armstrong McKay, S. Bathiany, L. Blaschke, N. Boers, D. Dylewsky, et al., Tipping point detection and early warnings in climate, ecological, and human systems, Earth System Dynamics 15 (4) (2024) 1117–1135. [6] E. Southall, T. S. Brett, M. J. Tildesley, L. Dyson, Early warning signals of infectious disease transitions: a review, Journal of the Royal Society Interface 18 (182) (2021) 20210555. [7] Z. Zhang, Q. Shen, X. Wang, Parameter identification framework of nonlinear dynamical systems with markovian switching, Chaos: An Interdisciplinary Journal of Nonlinear Science 33 (12) (2023) 123117. [8] L. Bettini, M. Cenedese, G. Haller, Model reduction to spectral submanifolds in piecewise smooth dynamical systems, International Journal of Non-Linear Mechanics 163 (2024) 104753. [9] X. Pan, H. Shu, L. Wang, X.-S. Wang, J. Yu, On the periodic solutions of switching scalar dynamical systems, Journal of Differential Equations 415 (2025) 365–382. [10] J. Lemus, B. Herrmann, Multi-objective sindy for parameterized model discovery from single transient trajectory data, Nonlinear Dynamics 113 (10) (2025) 10911–10927. [11] K. Egan, W. Li, R. Carvalho, Automatically discovering ordinary differential equations from data with sparse regression, Communications Physics 7 (1) (2024) 20. [12] O. Strebel, Preprocessing algorithms for the estimation of ordinary differential equation models with polynomial nonlinearities, Nonlinear Dynamics 111 (8) (2023) 7495–7510. [13] W. Zhai, D. Tao, Y. Bao, Parameter estimation and modeling of nonlinear dynamical systems based on runge–kutta physics-informed neural network, Nonlinear Dynamics 111 (22) (2023) 21117–21130. [14] B. Peherstorfer, K. Willcox, Data-driven operator inference for nonintrusive projection-based model reduction, Computer Methods in Applied Mechanics and Engineering 306 (2016) 196–215. [15] S. Panahi, Y.-C. Lai, Global phase-space approach to rate-induced tipping: A brief review, Chaos: An Interdisciplinary Journal of Nonlinear Science 35 (4) (2025) 043139. 42

[16] Z. Liu, X. Zhang, X. Ru, T.-T. Gao, J. M. Moore, G. Yan, Early predictor for the onset of critical transitions in networked dynamical systems, Physical Review X 14 (3) (2024) 031009. [17] S. L. Brunton, J. L. Proctor, J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 113 (15) (2016) 3932–3937. [18] S. H. Rudy, S. L. Brunton, J. L. Proctor, J. N. Kutz, Data-driven discovery of partial differential equations, Science Advances 3 (4) (2017) e1602614. [19] C. Truong, L. Oudre, N. Vayatis, Selective review of offline change point detection methods, Signal Processing 167 (2020) 107299. [20] P. Fryzlewicz, Wild binary segmentation for multiple change-point detection, The Annals of Statistics 42 (6) (2014). [21] H. Dette, W. Wu, Z. Zhou, Change point analysis of correlation in non-stationary time series, Statistica Sinica 29 (2) (2019) 611–643. [22] R. Baranowski, Y. Chen, P. Fryzlewicz, Narrowest-over-threshold detection of multiple change points and change-point-like features, Journal of the Royal Statistical Society Series B: Statistical Methodology 81 (3) (2019) 649–672. [23] H. Dette, J. Gösmann, A likelihood ratio approach to sequential change point detection for a general class of parameters, Journal of the American Statistical Association 115 (531) (2020) 1361–1377. [24] P. J. Green, D. I. Hastie, Reversible jump mcmc, Genetics 155 (3) (2009) 1391–1403. [25] J. Wang, E. Zivot, A bayesian time series model of multiple structural changes in level, trend, and variance, Journal of Business & Economic Statistics 18 (3) (2000) 374–386. [26] A. T. Levin, J. Piger, Bayesian model selection for structural break models, Available at SSRN 1132463 (2008). [27] N. A. Heard, M. J. Turcotte, Adaptive sequential monte carlo for multiple changepoint analysis, Journal of Computational and Graphical Statistics 26 (2) (2017) 414–423. [28] X. Xiu, Y. Yang, L. Kong, W. Liu, Laplacian regularized robust principal component analysis for process monitoring, Journal of Process Control 92 (2020) 212–219. [29] L. Birgé, Model selection via testing: an alternative to (penalized) maximum likelihood estimators, Annales de l’IHP Probabilités et Statistiques 42 (3) (2006) 273–325. [30] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, L. Yang, Physics-informed machine learning, Nature Reviews Physics 3 (6) (2021) 422–440. [31] N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, A. Anandkumar, Neural operator: Learning maps between function spaces with applications to pdes, Journal of Machine Learning Research 24 (89) (2023) 1–97. [32] S. Cuomo, V. S. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, F. Piccialli, Scientific machine learning through physics–informed neural networks: Where we are and what’s next, Journal of Scientific Computing 92 (3) (2022) 88. [33] G. Kissas, Y. Yang, E. Hwuang, W. R. Witschey, J. A. Detre, P. Perdikaris, Machine learning in cardiovascular flows modeling: Predicting arterial blood pressure from non-invasive 4d flow mri data using physics-informed neural networks, Computer Methods in Applied Mechanics and Engineering 358 (2020) 112623. [34] E. O. Oluwasakin, A. Q. Khaliq, Optimizing physics-informed neural network in dynamic system simulation and learning of parameters, Algorithms 16 (12) (2023) 547. 43

[35] G. Zhang, Y. Duan, G. Pan, Q. Chen, H. Yang, Z. Zhang, Data-driven discovery of state-changes in underlying system from hidden change-points in partial differential equations with spatiotemporal varying coefficients, Journal of Computational and Applied Mathematics (2025) 116962. [36] T. Lux, Estimation of regime-switching diffusions via fourier transforms, Statistics and Computing 34 (2) (2024) 88. [37] N. Niknejad, H. Modares, Physics-informed data-driven safe and optimal control design, IEEE Control Systems Letters 8 (2023) 285–290. [38] Y.-J. Huang, C.-W. Chang, C.-h. Hsieh, Detecting shifts in nonlinear dynamics using empirical dynamic modeling with nested-library analysis, PLOS Computational Biology 20 (1) (2024) e1011759. [39] M. Raissi, G. E. Karniadakis, Hidden physics models: Machine learning of nonlinear partial differential equations, Journal of Computational Physics 357 (2018) 125–141. [40] Y. Zhang, Y. Duan, X. Wang, Z. Zhang, Solving fokker–planck-kolmogorov equation by distribution self-adaptation normalized physics-informed neural networks, Physica A: Statistical Mechanics and Its Applications (2026) 131251. [41] M. Farid, Unsupervised data-driven response regime exploration and identification for dynamical systems, Chaos: An Interdisciplinary Journal of Nonlinear Science 34 (12) (2024). [42] J. Stiasny, G. S. Misyris, S. Chatzivasileiadis, Physics-informed neural networks for non-linear system identification for power system dynamics, in: 2021 IEEE Madrid PowerTech, IEEE, 2021, pp. 1–6. [43] S. Mishra, R. Molinaro, Estimates on the generalization error of physics-informed neural networks for approximating pdes, IMA Journal of Numerical Analysis 43 (1) (2023) 1–43. [44] Z. Lin, Y. Li, F. Yin, J. Maroñas, A. H. Thiéry, Efficient transformed gaussian process statespace models for non-stationary high-dimensional dynamical systems, IEEE Transactions on Signal Processing 73 (2025) 5229–5243. [45] G. Teschl, Ordinary differential equations and dynamical systems, Vol. 140, American Mathematical Soc., 2012. [46] J. J. Huerta y Munive, G. Struth, Predicate transformer semantics for hybrid systems: Verification components for isabelle/hol, Journal of Automated Reasoning 66 (1) (2022) 93–139. [47] A. Quarteroni, R. Sacco, F. Saleri, Numerical mathematics, Vol. 37, Springer Science & Business Media, 2006. [48] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics 378 (2019) 686–707. [49] S. Mishra, R. Molinaro, Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for pdes, IMA Journal of Numerical Analysis 42 (2) (2022) 981–1022. [50] Y. Qian, Y. Zhang, Y. Huang, S. Dong, Physics-informed neural networks for approximating dynamic (hyperbolic) pdes of second order in time: Error analysis and algorithms, Journal of Computational Physics 495 (2023) 112527. [51] F. Dondelinger, D. Husmeier, S. Rogers, M. Filippone, Ode parameter inference using adaptive gradient matching with gaussian processes, in: Artificial Intelligence and Statistics, PMLR, 2013, pp. 216–228. [52] R. Killick, P. Fearnhead, I. A. Eckley, Optimal detection of changepoints with a linear computational cost, Journal of the American Statistical Association 107 (500) (2012) 1590–1598.

Record · ID 141485 · SHA-256 5a9c9a5f10df4f11
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.