Paper deep dive
One-Step Evolution for Long-Time Extrapolation: An Error-Bound-Informed and Prior-Guided Neural Residual Framework for Autonomous PDEs
Maqun Zhang, Feng Gao, Wankun Chen, Hui Yu, Yanhai Gan, Junyu Dong
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%
Last extracted: 8/25/2026, 7:55:19 AM
Summary
The paper proposes a neural residual framework for long-time extrapolation of autonomous partial differential equations (PDEs) without ground-truth trajectory supervision. The method combines a low-cost numerical prior to generate a baseline state prediction with a convolutional neural network that learns a residual correction. This correction is constrained by a weak-form PDE residual loss, which serves as a proxy for the one-step approximation error. The approach is validated on five benchmark cases across four PDE classes, demonstrating superior long-time simulation accuracy compared to ten existing physics-informed learning methods.
Entities (8)
Relation Signals (7)
Neural Residual Framework → solves → Autonomous PDEs
confidence 98% · For a well-posed autonomous PDE... we propose a numerical-prior-guided, physics-constrained method
Neural Residual Framework → uses → Weak-Form PDE Residual
confidence 96% · a weak-form PDE residual provides a computable proxy for the one-step error term
Neural Residual Framework → uses → Numerical Prior
confidence 95% · a low-cost numerical prior reduces the difficulty of approximating the one-step evolution operator
Weak-Form PDE Residual → proxies → One-Step Error
confidence 94% · provides a computable proxy for the one-step error term in the error-propagation bound
Residual Compensation Network → corrects → Numerical Prior
confidence 93% · a correction network corrects errors in the prior update
Neural Residual Framework → outperforms → Physics-Informed Neural Networks (PINNs)
confidence 90% · outperforms the best competing baseline in each case
Neural Residual Framework → outperforms → Neural Operators
confidence 90% · outperforms the best competing baseline in each case
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Accurate simulation of the long-time evolution of systems governed by partial differential equations (PDEs) is central to scientific computing. Among existing deep learning?based approaches for solving PDEs, neural operators typically rely on extensive trajectory data, whereas physics-informed meth?ods often exhibit limited stability during long-time extrapolation. For a well-posed autonomous PDE, long-time trajectories can be generated by repeated composition of a fixed-step evolution operator; hence, long-time extrapolation depends on controlling the approximation error of this operator and the propagation of that error under recursive composition. Accordingly, we propose a numerical-prior-guided, physics-constrained method trained without ground-truth trajectory supervision: a low-cost numerical prior reduces the difficulty of approximating the one?step evolution operator, while a weak-form PDE residual provides a computable proxy for the one-step error term in the error?propagation bound. We validate the method on five benchmark cases spanning four PDE classes and compare it with ten physics?informed learning methods under a unified protocol that excludes ground-truth trajectories from training and model selection. The results indicate that, in all five cases, the proposed method reduces long-time extrapolation error relative to the numerical prior and outperforms the best competing baseline in each case, thereby improving long-time simulation accuracy across different PDEs without ground-truth trajectory supervision. The source code developed for this paper will be made publicly available upon acceptance of the manuscript.
Tags
Links
- Source: https://arxiv.org/abs/2608.22026v1
- Canonical: https://arxiv.org/abs/2608.22026v1
Trouble viewing inline? Open PDF directly →
Full Text
55,582 characters extracted from source content.
Expand or collapse full text
One-Step Evolution for Long-Time Extrapolation: An Error-Bound-Informed and Prior-Guided Neural Residual Framework for Autonomous PDEs Maqun Zhang Feng Gao Wankun Chen Hui Yu Yanhai Gan Junyu Dong Thanks: This work was supported by the Leverhulme Trust through Project VP1-2020-044; the National Natural Science Foundation of China under Grant 42406192; the Fundamental Research Funds for the Central Universities under Grants 202413040 and 202572015; the National Science and Technology Major Project of China under Grant 2022ZD0117201; and the Postdoctoral Project of Qingdao under Grant QDBSH20240102021. Thanks: Corresponding authors: Yanhai Gan (ganyanhai@ouc.edu.cn) and Junyu Dong (dongjunyu@ouc.edu.cn). Thanks: Maqun Zhang, Feng Gao, Wankun Chen, Yanhai Gan, and Junyu Dong are with the State Key Laboratory of Physical Oceanography and Faculty of Information Science and Engineering, Ocean University of China, Qingdao 266100, China. Hui Yu is with the School of Psychology and Neuroscience, University of Glasgow, Glasgow G12 8Q, U.K. Abstract Accurate simulation of the long-time evolution of systems governed by partial differential equations (PDEs) is central to scientific computing. Among existing deep learning-based approaches for solving PDEs, neural operators typically rely on extensive trajectory data, whereas physics-informed methods often exhibit limited stability during long-time extrapolation. For a well-posed autonomous PDE, long-time trajectories can be generated by repeated composition of a fixed-step evolution operator; hence, long-time extrapolation depends on controlling the approximation error of this operator and the propagation of that error under recursive composition. Accordingly, we propose a numerical-prior-guided, physics-constrained method trained without ground-truth trajectory supervision: a low-cost numerical prior reduces the difficulty of approximating the one-step evolution operator, while a weak-form PDE residual provides a computable proxy for the one-step error term in the error-propagation bound. We validate the method on five benchmark cases spanning four PDE classes and compare it with ten physics-informed learning methods under a unified protocol that excludes ground-truth trajectories from training and model selection. The results indicate that, in all five cases, the proposed method reduces long-time extrapolation error relative to the numerical prior and outperforms the best competing baseline in each case, thereby improving long-time simulation accuracy across different PDEs without ground-truth trajectory supervision. The source code developed for this paper will be made publicly available upon acceptance of the manuscript. Index Terms: autonomous PDEs, long-time extrapolation, numerical prior, physics-constrained learning without ground-truth trajectories, weak-form loss. I INTRODUCTION Partial differential equations (PDEs) are widely used to describe spatiotemporal processes, including fluid motion, heat transfer, wave propagation, and phase-field evolution. Accurate simulation of the long-time evolution of these processes from prescribed initial conditions is fundamental to scientific computing and prediction in engineering [1, 2]. For a well-posed autonomous PDE with specified physical parameters and boundary conditions, the state transition over a fixed time interval can be represented by a well-defined evolution operator [3, 4]. Deep neural networks can approximate complex nonlinear mappings and continuous operators, thereby providing a new means of learning this evolution operator [5, 6]. Deep learning-based PDE solvers primarily follow two paradigms: data-driven neural operators and physics-informed methods [7, 8, 9]. Neural operators such as DeepONet and the Fourier neural operator (FNO) learn complex function-to-function mappings but typically require paired samples generated by high-fidelity solvers. For time-dependent problems, these samples take the form of trajectories, and generalization error depends on the input-function distribution and training-sample coverage [7, 8, 10]. Physics-informed neural networks (PINNs) and their variants reduce reliance on ground-truth data by incorporating governing equations and initial and boundary conditions [9]. However, conventional PINNs directly optimize the full spatiotemporal solution and remain subject to loss ill-conditioning, gradient imbalance, and temporal-causality difficulties for complex equations and long time intervals [11, 12, 13]. Physics-constrained autoregressive networks and PhyCRNet can recursively solve time-dependent PDEs without ground-truth trajectory supervision [14, 15], but one-step errors accumulate during rollout, limiting the transfer of training-window accuracy to long-time extrapolation [16, 17]. Numerical methods can provide physically meaningful state updates for deep learning models [18, 19, 20], but finite resolution still introduces discretization errors [21, 22]. Existing numerical-learning hybrid methods often rely on high-fidelity trajectory supervision or corrections tailored to specific numerical schemes and are therefore not readily applicable to long-time extrapolation without ground-truth trajectories [18, 19, 20]. Accordingly, this work investigates how a numerical prior can reduce the difficulty of approximating the evolution operator and constrain its approximation error without ground-truth trajectory supervision, allowing a model trained within a short time window to be applied to long-time extrapolation. For an autonomous PDE that is well-posed over the target time horizon, the solution flow forms a time-homogeneous semigroup; consequently, long-time evolution can be represented by repeated composition of the same fixed-step state-transition operator [3, 4]. This structure allows the one-step dynamics learned within a short time window to be reused at later times, but errors introduced by repeated application of the approximate operator continue to propagate; therefore, the recursive structure of the exact evolution does not automatically translate into long-time model accuracy [16, 17]. Over the set of states reachable during extrapolation, long-time error is jointly controlled by the uniform approximation error of the one-step operator and system stability. Accordingly, using an existing physical approximation and correcting its one-step error, rather than learning the full evolution operator, provides a natural means of reducing approximation difficulty and improving long-time extrapolation. To realize this one-step error correction, the numerical prior first produces a state prediction, while a fully convolutional network learns a correction; the next state is obtained through additive fusion and hard enforcement of physical constraints. Because the ideal correction is unavailable, we construct a physics-based loss from the weak-form PDE residual [23] and, based on Eq. (7), use it as a conditional a posteriori proxy for the one-step approximation error of the evolution operator, thereby training the correction network without ground-truth trajectory supervision. The model is trained only within a short time window and then applied recursively beyond the training interval. The main contributions are summarized as follows: 1. An operator-approximation framework for long-time extrapolation with neural PDE solvers. Starting from the existence and recursive structure of the fixed-step evolution operator for a well-posed autonomous PDE, we apply the one-step error-propagation relation to neural approximation models, identify the one-step approximation error as the key quantity linking short-time training to long-time extrapolation, and establish the design principles for numerical-prior guidance and physics-based constraints without ground-truth trajectory supervision. 2. A prior-correction and error-constraint mechanism without ground-truth trajectory supervision. A low-cost numerical prior provides the baseline evolution update, a correction network corrects errors in the prior update, and the weak-form PDE residual serves as a conditional a posteriori proxy for the one-step error; hard constraints preserve the initial and boundary conditions. 3. A unified evaluation across PDEs and solution structures. Under a unified protocol that excludes ground-truth trajectories from training and model selection, we compare the proposed method with ten physics-informed learning methods on five benchmark cases spanning four PDE classes and systematically evaluate its long-time extrapolation performance relative to the numerical prior and the competing methods. I RELATED WORK I-A Data-Driven Neural Operators Neural operators directly learn function-space mappings from initial and boundary conditions, equation coefficients, or source terms to PDE solutions. DeepONet, FNO, and neural operator theory form the foundations of this line of research [6, 7, 8]. Recent surveys have organized the broader landscape of deep neural networks for PDEs, including physics-informed solvers and neural operators [42]. Subsequent studies incorporated PDE constraints through PINO [24] and extended neural operators to complex geometries [25], cross-equation pretraining [26, 27], and time-invariant evolution modeling [28]. Recent graph neural operator architectures further exploit multiscale spatial and frequency-domain features for PDE solving and retain strong performance with limited training samples and low-resolution data [43]. Although these methods improve operator representation and transferability, most still rely on paired solution data or large-scale pretraining trajectories, limiting their long-time extrapolation by data costs and training-distribution coverage. I-B Physics-Informed PDE Solvers PINNs incorporate governing equations and initial and boundary conditions into the loss function, thereby directly approximating PDE solutions without ground-truth data [9]; however, multi-term losses involving differential operators often lead to gradient imbalance and ill-conditioned optimization [11, 12]. To address these issues, VPINN [29] and hp-VPINN [23] introduce weak-form constraints; SA-PINN [30], gPINN [31], and Causal PINN [13] improve training through loss weighting, gradient enhancement, and temporal causality, respectively; and PINNsFormer [32] and RoPINN [33] further improve network architectures and sampling strategies. Benchmark studies indicate that these variants still exhibit inconsistent performance on complex PDEs [34]; moreover, most methods directly optimize the full solution over a prescribed spatiotemporal domain, so extrapolation beyond the temporal domain is not a natural output of their formulations. I-C Autoregressive Neural PDE Solvers Autoregressive methods learn fixed-step state-transition mappings and generate long-time trajectories by repeatedly applying the same model. PDE-Net was among the early methods to employ constrained convolutions for dynamical time stepping [35]. AR-DenseED and PhyCRNet further combined physical constraints with convolutional recurrent architectures to enable recursive PDE solution without ground-truth trajectory supervision [14, 15], while graph-network methods extended this paradigm to irregular discretizations [16, 17]. The main challenges are train–inference distribution shift and the accumulation of one-step errors; PDE-Refiner, APEBench, and recurrent neural operators investigate long-rollout stability through prediction refinement, unified benchmarking, and recurrent training, respectively [36, 37, 38]. However, most methods still rely on ground-truth trajectories, while the long-time accuracy of physics-constrained models trained without ground-truth trajectory supervision remains controlled by the one-step approximation error and its propagation. I-D Hybrid Numerical–Learning PDE Solvers Hybrid numerical–learning methods retain parts of the discrete structure or mechanistic model, while neural networks learn unknown dynamics, closure terms, or discretization errors. Interpretable spatiotemporal neural models have also been developed for PDE-governed distributed parameter systems [44]. PDE information has further been incorporated as a physics-informed prior in Bayesian neural networks to obtain physically consistent forecasts from limited observations [45]. Representative approaches include Universal Differential Equations, which combine known equations with learnable terms [39]; methods that learn corrections to coarse-grid discretizations [18, 40]; and Solver-in-the-Loop [20], machine-learning-accelerated computational fluid dynamics [19], and Differentiable Turbulence [41], which incorporate differentiable solvers into training. These studies demonstrate that numerical structures can effectively constrain the learning process, but many rely on high-fidelity trajectories or are tailored to specific solvers. In contrast, the proposed method obtains the baseline evolution through low-cost numerical time stepping and learns an error correction from weak-form residuals without ground-truth trajectory supervision. I PRIOR-GUIDED PHYSICS-CONSTRAINED EXTRAPOLATION WITHOUT GROUND-TRUTH TRAJECTORY SUPERVISION I-A Evolution-Operator Foundations and Methodological Motivation Fig. 1: Overall architecture and training workflow of the proposed method. The theoretical relations at the top connect recursive fixed-step evolution, one-step error propagation, and weak-form error constraints, with colored terms corresponding to their modules. Low-Cost Numerical Prior obtains the prior prediction through restriction, prior stepping, and prolongation. Residual Compensation Network concatenates the current state and prior prediction at C and outputs a correction. Residual Fusion & Hard Physics scales the correction by s, adds it to the prior prediction, and applies Πℳ _M to produce the next state for recursive rollout. Weak-Form Physics-Unsupervised Optimization constructs ℒweakL_weak from the PDE residuals of short predicted trajectories through test-function projection and weighted aggregation, enabling training without ground-truth trajectory supervision. As shown in Fig. 1, the proposed framework integrates a low-cost numerical prior, a residual compensation network, hard physical projection, and weak-form physics-unsupervised optimization for recursive long-time rollout. This section starts from the temporal evolution law of autonomous PDEs and establishes a fixed-step evolution operator under the assumption that the corresponding initial-boundary value problem is well posed. It then analyzes the error between the exact evolution operator and its approximation, together with its propagation, thereby clarifying the roles of the numerical prior, neural correction, and weak-form constraints in improving operator approximation and long-time extrapolation. I-A1 Fixed-Step Evolution Operator for Autonomous PDEs This work considers long-time extrapolation for autonomous PDEs. Let u(t)u(t) denote the system state. The governing equations are written in the general form ∂tu=(u,μ),ℬ(u)=b,u(0)=u0. _tu=G(u;μ), (u)=b, u(0)=u_0. (1) Here, G denotes the evolution law defined by the PDE, μ denotes the fixed physical parameters, and ℬ(u)=bB(u)=b represents time-invariant boundary conditions. Autonomy means that G has no explicit dependence on absolute time: given the same current state, the system follows the same evolution over an equal subsequent time interval. Assuming that the corresponding initial-boundary value problem is well posed, the exact evolution from the initial state to time t can be represented by the solution flow t S_t. Autonomy further yields the semigroup relation u(t)=t(u0),t+s=t∘s,0=ℐ.u(t)= S_t(u_0), S_t+s= S_t S_s, S_0=I. (2) For a fixed time interval τ, define the exact one-step evolution operator as τ=τ T_τ= S_τ. Then, un+1=τ(un),un=τn(u0).u_n+1= T_τ(u_n), u_n= T_τ^\,n(u_0). (3) As shown in (3), exact long-time evolution can be generated by repeated composition of the same one-step operator. We therefore reformulate the learning target from the full spatiotemporal solution to a one-step evolution operator that can be applied recursively. However, this result applies only to the exact operator and does not imply that its learned approximation ^τ,θ T_τ,θ also supports stable extrapolation. The approximate model introduces an error at each step, which continues to propagate under recursive application. Therefore, connecting the recursive structure of the exact operator to the long-time extrapolation capability of the model requires a relation between the one-step approximation error and its propagation. Building on this relation, the next subsection introduces the numerical prior, neural correction, and physics-based loss. I-A2 Approximation of the Exact Evolution Operator and the Proposed Method The exact evolution operator can be repeatedly composed, but a practical model can only approximate it. Let K denote the set of states reachable during extrapolation. We denote the learnable one-step evolution operator by ^τ,θ T_τ,θ and define its uniform one-step error as u^n+1=^τ,θ(u^n),εθ=supu∈‖^τ,θ(u)−τ(u)‖X. u_n+1= T_τ,θ( u_n), _θ= _u \| T_τ,θ(u)- T_τ(u) \|_X. (4) If the exact evolution operator has a local stability constant LτL_τ on K, the predicted states remain in K, and both trajectories share the same initial state, then the error en=‖u^n−un‖Xe_n=\| u_n-u_n\|_X at step n satisfies en+1≤Lτen+εθ,en≤εθ∑j=0n−1Lτj.e_n+1≤ L_τe_n+ _θ, e_n≤ _θ _j=0^n-1L_τ^j. (5) Thus, the recursive structure of the exact operator only provides the basis for extrapolation; long-time model accuracy depends on the one-step operator error and its propagation. To reduce the difficulty of approximating the full operator τ T_τ, we introduce a low-cost numerical prior τ,hP_τ,h and let the network learn the ideal correction to the prior error: τ,h(u) _τ,h(u) =τ(u)−τ,h(u),^τ,θ(u) = T_τ(u)-P_τ,h(u), T_τ,θ(u) (6) =Πℳ[τ,h(u)+θ(u,τ,h(u))]. = _M [P_τ,h(u)+C_θ (u,P_τ,h(u) ) ]. Here, θC_θ approximates the ideal correction operator τ,hD_τ,h, while Πℳ _M enforces the necessary physical constraints. The numerical prior provides a baseline approximation, and the convolutional neural network (CNN) learns only the correction to the prior error, thereby narrowing the operator-approximation range when the prior captures the dominant dynamics. However, neither the exact operator nor the ideal correction is available during training, so Eq. (4) cannot be optimized directly. Both strong-form and weak-form residuals can provide physics-based criteria without ground truth; the weak form is adopted because stability estimates for evolution problems typically control state error through the dual norm of the spatiotemporal residual, which weak testing can discretize directly. If the PDE is stable and the test space is sufficiently rich, ‖v(τ)−τ(u)‖X \|v(τ)- T_τ(u) \|_X ≤Cstab‖ℛ(v,u)‖′ ≤ C_stab \|R(v;u) \|_Y (7) ≤Cstab(CΦ1/2ℒweak1/2(v,u)+δΦ). ≤ C_stab (C_ ^1/2L_weak^1/2(v;u)+ _ ). Here, ′Y denotes the dual space of the spatiotemporal residual over one time step. The weak-form loss therefore converts the unavailable one-step operator error into a computable residual proxy that constrains how accurately ^τ,θ T_τ,θ approximates the one-step action of τ T_τ. The conditions and derivation of Eq. (7) are provided in Appendix A. During training, the same approximate operator is applied recursively over a short time window, while the weak-form residual constrains the entire predicted trajectory. Parameter sharing ensures that each time step follows the same evolution rule, and Eq. (5) shows that reducing the one-step approximation error correspondingly lowers its accumulated upper bound under recursive application. In summary, the numerical prior provides the baseline approximation, the CNN learns the prior correction, and the weak-form loss constrains the operator error without ground truth; together, these components form the proposed pathway to long-time extrapolation. I-B Prior-Guided Residual Correction Operator This section presents the computational implementation of the approximate evolution operator ^τ,θ T_τ,θ. Given the current state unu_n, a fixed numerical prior first produces a baseline prediction, after which the Residual CNN estimates a correction to the prior error. The next state is then obtained through residual fusion and the necessary physical constraints. A forward pass is summarized as un⟶un+1P⟶Δun+1⟶u^n+1.u_n u_n+1^P u_n+1 u_n+1. Neither the numerical prior nor the physical constraints contain trainable parameters; all network parameters reside in the correction operator θC_θ. I-B1 Cross-Resolution Prior Embedding Let hfh_f and hch_c denote the spatial scales of the target and prior grids, respectively, with hc≥hfh_c≥ h_f. We represent the numerical prior as a fixed one-step time-marching operator τ,hcP_τ,h_c. The cross-resolution computation is un+1P=ℐhc→hf[τ,hc(ℛhf→hc(un))].u_n+1^P=I_h_c→ h_f [P_τ,h_c (R_h_f→ h_c(u_n) ) ]. (8) Here, ℛhf→hcR_h_f→ h_c is the restriction operator from the target grid to the prior grid, and ℐhc→hfI_h_c→ h_f is the prolongation operator. The prior module only needs to receive the current state, advance it by a fixed time step τ, and return a prediction containing the same physical variables as the original state; it is not tied to a specific numerical scheme. Fig. 2: Full-state RMSE for TGV over steps 0–200. Training ends at step 20. Fig. 3: TGV vorticity field at the representative extrapolation step 100. Because the analytical solution decays rapidly over time, the field at step 200 is provided in Fig. 19. Fig. 4: RMSE for the smooth heat equation over steps 0–200. Training ends at step 20. The prior is evaluated on a coarser grid to reduce the cost of baseline time stepping. If the number of nodes along each spatial direction decreases from NfN_f to NcN_c, the prior degrees of freedom for a local grid-based method in d dimensions decrease to approximately (Nc/Nf)d(N_c/N_f)^d of those on the target grid. This ratio applies only to the prior computation; prolongation and the CNN remain on the target grid. Errors introduced by coarse-grid discretization, restriction, and prolongation jointly constitute the prior error to be corrected by the network. The specific prior schemes, grid resolutions, and time steps are given in Appendix C and Appendix F-B. I-B2 Prior-Error Correction with a Residual CNN The correction network approximates the error unresolved by the prior, and its input must contain the complete state required to determine the next-step evolution. We concatenate the current state and the upsampled prior prediction along the channel dimension and output a correction of the same dimensionality: zn=Concat(un,un+1P),Δun+1=θ(zn).z_n=Concat (u_n,u_n+1^P ), u_n+1=C_θ(z_n). (9) The network architecture must balance representation capacity, computational cost, and physical identifiability. Its effective receptive field should cover the dominant spatial scale of the prior error over one time step, without requiring global coupling among all locations. Because the network is repeatedly evaluated at every extrapolation step, we adopt a lightweight, moderate-capacity, fully convolutional architecture to limit both computational overhead and degrees of freedom left unconstrained by weak-form physical testing. Convolutional weights are shared over the full spatial field, with padding consistent with the boundary conditions. If tiling is necessary, overlapping regions are retained to avoid boundary artifacts from independent nonoverlapping patches. A Residual CNN is adopted as the standard implementation. An initial 3×33× 3 convolution maps znz_n to C hidden channels, followed by B residual blocks. Each block contains two 3×33× 3 convolutions, GroupNorm, and SiLU activations, and updates the features through a skip connection: F(ℓ+1)=F(ℓ)+ℋℓ(F(ℓ)),ℓ=0,…,B−1.F^( +1)=F^( )+H_ (F^( ) ), =0,…,B-1. (10) The network width C, the number of residual blocks B, and any required dilation rates jointly determine its capacity and receptive field; their values are selected according to the spatial scale of the prior error and the target resolution. The output convolution maps the hidden features to correction terms for the physical variables. Its weights and biases are initialized to zero, so at the start of training, θ0(zn)=0,u^n+1=Πℳ(un+1P).C_ _0(z_n)=0, u_n+1= _M (u_n+1^P ). (11) If the prior prediction already satisfies the physical constraints, the initial model is exactly equal to the prior. Zero initialization and residual output jointly limit the initial perturbation, allowing the network to correct the prior progressively rather than relearn the full evolution. The network size and boundary treatment are specified in Appendix F-A; temporal consistency is maintained by the parameter-shared one-step operator and the physics-based loss introduced in Section I-C. I-B3 Residual Fusion and Physical Constraints The network output represents an additive correction to the prior prediction. To account for scale differences among physical variables, the correction is first scaled channel-wise and then fused with the prior prediction: u~n+1=un+1P+s⊙Δun+1. u_n+1=u_n+1^P+s u_n+1. (12) Here, s is the variable-scale vector, and ⊙ denotes channel-wise multiplication. The scales are determined by nondimensionalization or prescribed physical magnitudes and do not depend on ground-truth supervision. The fused result is then mapped to the physically admissible set to obtain the final state: u^n+1=Πℳ(u~n+1)=^τ,θ(un). u_n+1= _M ( u_n+1 )= T_τ,θ(u_n). (13) The operator Πℳ _M enforces the hard constraints required by the specific PDE, such as periodic or fixed boundary conditions, incompressibility, or variable positivity. If explicit projection is unnecessary for a given case, Πℳ=ℐ _M=I. The numerical prior, CNN correction, and physical constraints thus constitute a complete one-step evolution. Parameter optimization and recursive application are introduced in the next section. Fig. 5: Temperature field for the smooth heat equation at the representative extrapolation step 120. The field at step 200 is provided in Fig. 20. Fig. 6: Full-state RMSE for DSL over steps 0–500. Both the raw curve for Ours and its 15-step centered moving average are shown, with the latter included only for readability. All statistics are computed from the unsmoothed data. Fig. 7: DSL vorticity field at the final extrapolation step 500. Fig. 8: Full-state RMSE for Gaussian Mound over steps 0–200. Fig. 9: Free-surface perturbation and velocity field for Gaussian Mound at step 200. I-C Physics-Constrained Optimization without Ground-Truth Trajectory Supervision and Recursive Extrapolation This section describes how the composite operator is optimized without ground-truth trajectory supervision. The weak-form loss constrains short-window predictions to satisfy the governing equations, while necessary PDE-specific physical regularization further restricts inadmissible states. After training, the parameters are frozen, and the same operator is repeatedly applied for long-time extrapolation. I-C1 Weak-Form Physics Constraint for One-Step Evolution Starting from state unu_n, the model generates a predicted sequence of length K: u^n+k=^τ,θk(un),k=1,…,K. u_n+k= T_τ,θ^\,k(u_n), k=1,…,K. (14) Because the exact evolution operator is unavailable during training, the predicted sequence is evaluated through the governing equation. For the abstract equation ∂tu=(u) _tu=G(u), define the residual as ℛ(u^)=∂tu^−(u^).R( u)= _t u-G( u). (15) Select test functions ϕmm=1M\ _m\_m=1^M over the training spatiotemporal window Q, and define rm=∫Qℛ(u^)ϕmt,ℒweak=1M∑m=1Mωm|rm|2.r_m= _QR( u)\, _m\,d x\,dt, _weak= 1M _m=1^M _m|r_m|^2. (16) The conditions under which this loss serves as a proxy for the one-step error are given in Appendix A; its additional numerical and optimization advantages are discussed in Appendix B, and the test-function sets and quadrature rules used for each case are reported in Appendix F-C. Rigorous a posteriori error bounds for PDE-defined PINNs have likewise shown that computable prediction-error estimates can be obtained without access to the exact solution [46]. All initial and boundary conditions are enforced as hard constraints. The initial condition directly initializes the recursion, while the boundary conditions are preserved through boundary-consistent numerical operations and the physical projection Πℳ _M. Therefore, no initial-condition or boundary-condition losses are added to Eq. (16). I-C2 Overall Objective and Training without Ground-Truth Trajectory Supervision The overall training objective is ℒ(θ)=ℒweak+∑r=1Rλrℒphys(r).L(θ)=L_weak+ _r=1^R _rL_phys^(r). (17) Here, ℒphys(r)L_phys^(r) denotes PDE-specific physical regularization, such as constraints on divergence, conserved quantities, equilibrium states, or correction magnitude. These terms exclude nonphysical solutions that are not sufficiently constrained by the weak-form tests. Different physical variables are nondimensionalized or scale-normalized before evaluation to prevent variables with larger numerical magnitudes from dominating the optimization. Because the initial and boundary conditions are enforced by construction, the objective contains no initial-condition or boundary-condition losses and no ground-truth fitting terms. Training is performed only over a prescribed short time window. For each window, the current model recursively generates a predicted sequence, after which Eq. (17) is evaluated and the parameters are updated. Model selection is based on physical residuals evaluated over held-out time windows and independent test functions. Reference numerical or analytical solutions are excluded from both the loss and model selection and are used only for error evaluation. The training windows, loss weights, and optimizer settings are given in Appendix F-D. I-C3 Recursive Extrapolation with Frozen Parameters After training, the parameters θ∗θ^* are frozen, and the same composite operator is repeatedly applied from the prescribed initial state: u^0=u0,u^n+1=^τ,θ∗(u^n),n=0,1,…,N−1. u_0=u_0, u_n+1= T_τ,θ^*( u_n), n=0,1,…,N-1. (18) The parameters θ∗θ^* are shared across all extrapolation steps. The state is not reset at the end of the training window, and neither ground-truth correction nor further network updates are introduced. Consequently, the one-step operator learned over the short training window is applied directly beyond the training interval, and its long-time performance depends entirely on one-step physical consistency and stability under recursive composition. IV NUMERICAL EXPERIMENTS AND RESULTS This section evaluates the long-time extrapolation capability of the proposed method without ground-truth trajectory supervision using five PDE cases with distinct dynamics. It first presents the unified experimental setup and evaluation protocol, then compares the errors and field evolution of different methods for each case, and finally analyzes the sustained correction provided by the network relative to the numerical prior across all cases. IV-A Experimental Setup and Evaluation Protocol We select five benchmark cases with distinct dynamics to evaluate long-time extrapolation after training over a short time window. All methods roll out continuously from step 0 to the test endpoint. Table I summarizes the case settings and validation targets; the complete problem definitions, reference-solution sources, and numerical priors are provided in Appendix C, while the accuracy characterization and information-isolation protocol for the reference solutions are specified in Appendix C-F. TABLE I: Benchmark Settings and Validation Targets Case Governing equation Training/test steps Evaluation grid Validation target Taylor–Green vortex (TGV) Incompressible Navier–Stokes equations 0–20 / 0–200 1282128^2 Smooth analytical flow Smooth Heat Equation Heat equation 0–20 / 0–200 1282128^2 Dissipation and near-zero stability Double shear layer (DSL) Incompressible Navier–Stokes equations 0–150 / 0–500 1282128^2 Nonlinear instability and vortex structures Gaussian Mound Shallow-water equations 0–20 / 0–200 2562256^2 Wave propagation and variable coupling Cahn–Hilliard (CH) Two-Droplet Cahn–Hilliard equation 0–20 / 0–260 2562256^2 Fourth-order conservation and interface evolution Training and test steps denote the physics-based optimization window and the full rollout interval, respectively. We compare the proposed method with PINN [9], PINNsFormer [32], RoPINN [33], gPINN [31], Causal PINN [13], hp-VPINN [23], SA-PINN [30], PINO trained without ground-truth supervision [24], AR-DenseED [14], and PhyCRNet [15]. AR-DenseED and PhyCRNet represent a physics-constrained autoregressive convolutional network and a convolutional recurrent PDE solver, respectively. Each method retains the core architecture and training mechanism of the original work and adopts the same governing equations, initial and boundary conditions, training window, and evaluation grid for each case. The adaptation and model-selection protocols are detailed in Appendix E. All state variables are evaluated using per-step mean absolute error (MAE) and root mean square error (RMSE). The main text reports RMSE curves and representative extrapolated fields and analyzes the error reduction of the proposed method relative to the numerical prior in Section IV-C. Complete MAE results and longest-horizon fields are provided in Appendix D, the evaluation protocol in Appendix E, and the full implementation and randomness settings of the proposed method in Appendix F. IV-B Long-Time Extrapolation Results for Different Evolution Equations This section compares the long-time extrapolation performance of the proposed method and ten physics-informed learning methods across five benchmark cases. The main text reports per-step RMSE and representative extrapolated fields; complete problem definitions are provided in Appendix C, and the correction relative to the numerical prior is analyzed across all cases in Section IV-C. The vertical dashed lines in the curves indicate the end of training. IV-B1 TGV: Smooth Analytical Flow Fig. 2 presents the full-state RMSE over the rollout, and Fig. 3 shows the vorticity field at step 100. This case evaluates extrapolation accuracy for a smooth analytical flow. Over steps 21–200, the average RMSEs of the proposed method and the best-performing baseline, PINN, are 3.86×10−53.86× 10^-5 and 9.47×10−29.47× 10^-2, respectively. At step 100, the proposed method retains a vortex structure consistent with the reference solution, whereas the competing methods exhibit clear deviations. These results show that the proposed method maintains high long-time accuracy for smooth evolution. IV-B2 Smooth Heat Equation: Linear Dissipative Evolution Fig. 4 presents the RMSE over the rollout, and Fig. 5 shows the temperature field at step 120. This case evaluates stability during long-time dissipation and the near-zero regime. Over steps 21–200, the average RMSEs of the proposed method and the best-performing baseline, PINNsFormer, are 9.36×10−59.36× 10^-5 and 1.88×10−21.88× 10^-2, respectively. At step 120, the proposed method preserves a smooth temperature distribution and the correct decay trend without evident rebound or nonphysical oscillations. These results show that the proposed method stably extrapolates linear dissipative evolution. IV-B3 DSL: Nonlinear Shear Instability Fig. 6 presents the full-state RMSE over the rollout, and Fig. 7 shows the vorticity field at step 500. This case evaluates error growth and vortex-structure preservation after nonlinear instability. Over steps 151–500, the average RMSEs of the proposed method and the best-performing baseline, Causal PINN, are 4.77×10−24.77× 10^-2 and 5.04×10−15.04× 10^-1, respectively. At step 500, the competing methods generally exhibit excessive smoothing or structural displacement, whereas the proposed method retains the primary vortex rolls and their surrounding structures. These results show that the proposed method better preserves spatial organization after nonlinear instability. IV-B4 Gaussian Mound: Smooth Hyperbolic Wave Propagation Fig. 8 presents the full-state RMSE over the rollout, and Fig. 9 shows the free-surface perturbation and velocity field at step 200. This case evaluates the extrapolation of coupled variables during hyperbolic wave propagation. Over steps 21–200, the average RMSEs of the proposed method and the best-performing baseline, gPINN, are 1.71×10−31.71× 10^-3 and 3.99×10−33.99× 10^-3, respectively. At step 200, the proposed method accurately preserves the propagation position, symmetry, and velocity direction of the annular wave. These results show that the proposed method reduces multivariable coupling errors in smooth hyperbolic systems. IV-B5 CH Two-Droplet: Fourth-Order Conservative Phase-Field Evolution Fig. 10 presents the RMSE over the rollout, and Fig. 11 shows the phase-field distribution at step 200. Fig. 10: RMSE for CH Two-Droplet over steps 0–260. Fig. 11: Phase-field distribution for CH Two-Droplet at the representative extrapolation step 200. This case evaluates interface and topology evolution in a fourth-order conservative system. Over steps 21–260, the average RMSEs of the proposed method and the best-performing baseline, gPINN, are 4.29×10−24.29× 10^-2 and 6.40×10−26.40× 10^-2, respectively. At step 200, the proposed method preserves clear phase-separation interfaces and the main droplet structures. These results show that the proposed method applies not only to second-order transport and diffusion problems but also to fourth-order conservative phase-field evolution. IV-C Cross-Case Analysis of Prior Correction To isolate the contribution of the learned correction in the full model, we compare its per-step error with that of the numerical prior and examine whether it consistently improves recursive predictions beyond the training interval. For each case, the full model and the numerical prior are rolled out continuously from step 0, and their RMSEs are computed. Let Ec,m(n)E_c,m^(n) denote the RMSE of method m for case c at step n, where m∈prior,fullm∈\prior,full\ denotes the numerical prior and the full model, respectively. The normalized rollout progress and peak-normalized error reduction are defined as sc(n)=nNc,Ic(n)=100%Ec,prior(n)−Ec,full(n)max0≤j≤NcEc,prior(j).s_c^(n)= nN_c, I_c^(n)=100\% E_c,prior^(n)-E_c,full^(n) _0≤ j≤ N_cE_c,prior^(j). (19) Here, NcN_c is the final test step, and Ic(n)>0I_c^(n)>0 indicates that the full model has a lower error than the numerical prior. Normalization by the peak prior RMSE aligns the error scales across cases and avoids excessive amplification of pointwise relative ratios when the errors approach zero. Fig. 12: Per-step RMSE reduction of the proposed method relative to the numerical prior. The circles mark the end of training for each case, with the training interval on the left and the extrapolation interval on the right; positive values indicate that Ours outperforms the numerical prior. The corresponding MAE results are provided in Appendix D. As shown in Fig. 12, all curves start from zero gain at the common initial state, and all five cases maintain positive gains at every extrapolation step. TGV and Heat are dominated by solution-amplitude decay, and their gains gradually decrease in the middle and late stages. DSL and CH are dominated by vortex-structure and interface evolution over the test interval, and their gains generally increase during rollout. Gaussian Mound maintains a relatively stable positive gain during wave propagation. These differences show that the correction magnitude does not vary monotonically with equation complexity but depends on the accuracy of the numerical prior and how its error accumulates during evolution. The MAE results in Appendix D follow the same trend as the RMSE results, showing that the learned correction consistently reduces the prior error over the extrapolation intervals of different PDEs and that the performance improvement does not arise solely from the numerical prior. V CONCLUSION This study investigates long-time extrapolation without ground-truth trajectory supervision for deep learning-based PDE solvers. To address the dependence of neural operators on trajectory data and the tendency of physics-informed methods to accumulate errors during recursive extrapolation, we start from the existence of the discrete evolution operator for an autonomous PDE and the propagation of its approximation error, interpret long-time rollout as repeated approximation of the same exact evolution mapping, and characterize extrapolation stability in terms of controlling the one-step approximation error and its recursive propagation. Based on this understanding, we develop a numerical-prior-guided, physics-constrained framework for extrapolation without ground-truth trajectory supervision. To reduce the difficulty of approximating the full evolution operator, we introduce a low-cost numerical prior to provide the baseline evolution; to correct the error unresolved by the prior, we further construct a correction network; and to constrain the one-step error without ground-truth trajectory supervision, we employ the weak-form PDE residual as a computable proxy for its upper bound, following Eq. (7). Consequently, after training, the model can be applied recursively from the initial state, yielding a unified evolution approximation for long-time extrapolation. The numerical experiments include five benchmark cases—TGV, the smooth heat equation, DSL, Gaussian Mound, and Cahn–Hilliard Two-Droplet—spanning smooth analytical flow, linear dissipation, nonlinear shear instability, hyperbolic wave propagation, and fourth-order conservative phase-field evolution. The proposed method outperforms the corresponding numerical prior in every case and achieves more stable long-time extrapolation than multiple physics-informed learning baselines. A cross-case analysis of prior correction further shows that the learned correction consistently reduces the prior error throughout the extrapolation interval, indicating that the performance gains do not arise solely from the numerical prior. The primary role of the method is not to replace numerical evolution rules, but to build on an existing physical prior to construct a more tractable and tightly constrained approximation of the evolution operator, thereby mitigating error accumulation during long-time rollout. This study has several limitations. Current experiments focus mainly on two-dimensional regular grids with fixed initial and boundary conditions, and applicability to complex geometries, multiphysics coupling, high-dimensional systems, and parametric PDE families requires further validation. Moreover, the test-function set currently requires PDE-specific adjustment; future work will investigate automatic selection algorithms and their effects on the test-space error term in Eq. (7). Future work will also analyze the relationship between numerical-prior accuracy and long-time error propagation and extend the framework to more complex scientific computing settings. References [1] L. C. Evans, Partial Differential Equations, 2nd ed. Providence, RI, USA: American Mathematical Society, 2010, doi: 10.1090/gsm/019. [2] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, “Physics-informed machine learning,” Nat. Rev. Phys., vol. 3, p. 422–440, 2021, doi: 10.1038/s42254-021-00314-5. [3] A. Pazy, Semigroups of Linear Operators and Applications to Partial Differential Equations. New York, NY, USA: Springer, 1983, doi: 10.1007/978-1-4612-5561-1. [4] M. G. Crandall and T. M. Liggett, “Generation of semi-groups of nonlinear transformations on general Banach spaces,” Amer. J. Math., vol. 93, no. 2, p. 265–298, 1971, doi: 10.2307/2373376. [5] T. Chen and H. Chen, “Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems,” IEEE Trans. Neural Netw., vol. 6, no. 4, p. 911–917, 1995, doi: 10.1109/72.392253. [6] N. Kovachki et al., “Neural operator: Learning maps between function spaces with applications to PDEs,” J. Mach. Learn. Res., vol. 24, no. 89, p. 1–97, 2023. [7] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis, “Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators,” Nat. Mach. Intell., vol. 3, p. 218–229, 2021, doi: 10.1038/s42256-021-00302-5. [8] Z. Li et al., “Fourier neural operator for parametric partial differential equations,” in Proc. Int. Conf. Learn. Represent. (ICLR), 2021. [9] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” J. Comput. Phys., vol. 378, p. 686–707, 2019, doi: 10.1016/j.jcp.2018.10.045. [10] S. Lanthaler, S. Mishra, and G. E. Karniadakis, “Error estimates for DeepONets: A deep learning framework in infinite dimensions,” Trans. Math. Appl., vol. 6, no. 1, Art. no. tnac001, 2022, doi: 10.1093/imatrm/tnac001. [11] A. S. Krishnapriyan, A. Gholami, S. Zhe, R. M. Kirby, and M. W. Mahoney, “Characterizing possible failure modes in physics-informed neural networks,” in Adv. Neural Inf. Process. Syst., vol. 34, p. 26548–26560, 2021. [12] S. Wang, Y. Teng, and P. Perdikaris, “Understanding and mitigating gradient flow pathologies in physics-informed neural networks,” SIAM J. Sci. Comput., vol. 43, no. 5, p. A3055–A3081, 2021, doi: 10.1137/20M1318043. [13] S. Wang, S. Sankaran, and P. Perdikaris, “Respecting causality for training physics-informed neural networks,” Comput. Methods Appl. Mech. Eng., vol. 421, Art. no. 116813, 2024, doi: 10.1016/j.cma.2024.116813. [14] N. Geneva and N. Zabaras, “Modeling the dynamics of PDE systems with physics-constrained deep auto-regressive networks,” J. Comput. Phys., vol. 403, Art. no. 109056, 2020, doi: 10.1016/j.jcp.2019.109056. [15] P. Ren, C. Rao, Y. Liu, J. Wang, and H. Sun, “PhyCRNet: Physics-informed convolutional-recurrent network for solving spatiotemporal PDEs,” Comput. Methods Appl. Mech. Eng., vol. 389, Art. no. 114399, 2022, doi: 10.1016/j.cma.2021.114399. [16] A. Sanchez-Gonzalez, J. Godwin, T. Pfaff, R. Ying, J. Leskovec, and P. W. Battaglia, “Learning to simulate complex physics with graph networks,” in Proc. 37th Int. Conf. Mach. Learn. (ICML), vol. 119, p. 8459–8468, 2020. [17] J. Brandstetter, D. Worrall, and M. Welling, “Message passing neural PDE solvers,” in Proc. Int. Conf. Learn. Represent. (ICLR), 2022. [18] Y. Bar-Sinai, S. Hoyer, J. Hickey, and M. P. Brenner, “Learning data-driven discretizations for partial differential equations,” Proc. Natl. Acad. Sci. USA, vol. 116, no. 31, p. 15344–15349, 2019, doi: 10.1073/pnas.1814058116. [19] D. Kochkov, J. A. Smith, A. Alieva, Q. Wang, M. P. Brenner, and S. Hoyer, “Machine learning–accelerated computational fluid dynamics,” Proc. Natl. Acad. Sci. USA, vol. 118, no. 21, Art. no. e2101784118, 2021, doi: 10.1073/pnas.2101784118. [20] K. Um, Y. Fei, P. Holl, R. Brand, and N. Thuerey, “Solver-in-the-loop: Learning from differentiable physics to interact with iterative PDE-solvers,” in Adv. Neural Inf. Process. Syst., vol. 33, p. 6111–6122, 2020. [21] P. D. Lax and R. D. Richtmyer, “Survey of the stability of linear finite difference equations,” Commun. Pure Appl. Math., vol. 9, no. 2, p. 267–293, 1956, doi: 10.1002/cpa.3160090206. [22] R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems. Philadelphia, PA, USA: Society for Industrial and Applied Mathematics, 2007, doi: 10.1137/1.9780898717839. [23] E. Kharazmi, Z. Zhang, and G. E. Karniadakis, “hp-VPINNs: Variational physics-informed neural networks with domain decomposition,” Comput. Methods Appl. Mech. Eng., vol. 374, Art. no. 113547, 2021, doi: 10.1016/j.cma.2020.113547. [24] Z. Li et al., “Physics-informed neural operator for learning partial differential equations,” in Proc. Int. Conf. Learn. Represent. (ICLR), 2022. [25] Z. Li, D. Z. Huang, B. Liu, and A. Anandkumar, “Fourier neural operator with learned deformations for PDEs on general geometries,” J. Mach. Learn. Res., vol. 24, no. 388, p. 1–26, 2023. [26] Z. Hao et al., “DPOT: Auto-regressive denoising operator transformer for large-scale PDE pre-training,” in Proc. 41st Int. Conf. Mach. Learn. (ICML), vol. 235, p. 17616–17635, 2024. [27] M. Herde et al., “Poseidon: Efficient foundation models for PDEs,” in Adv. Neural Inf. Process. Syst., vol. 37, p. 72525–72624, 2024, doi: 10.52202/079017-2311. [28] Z. Zhou, W. Zhang, and L. Liu, “Time-invariant neural operators with applications in solving time-dependent PDEs,” arXiv preprint arXiv:2607.06188, 2026. [29] E. Kharazmi, Z. Zhang, and G. E. Karniadakis, “Variational physics-informed neural networks for solving partial differential equations,” arXiv preprint arXiv:1912.00873, 2019. [30] L. D. McClenny and U. M. Braga-Neto, “Self-adaptive physics-informed neural networks,” J. Comput. Phys., vol. 474, Art. no. 111722, 2023, doi: 10.1016/j.jcp.2022.111722. [31] J. Yu, L. Lu, X. Meng, and G. E. Karniadakis, “Gradient-enhanced physics-informed neural networks for forward and inverse PDE problems,” Comput. Methods Appl. Mech. Eng., vol. 393, Art. no. 114823, 2022, doi: 10.1016/j.cma.2022.114823. [32] Z. Zhao, X. Ding, and B. A. Prakash, “PINNsFormer: A transformer-based framework for physics-informed neural networks,” in Proc. Int. Conf. Learn. Represent. (ICLR), 2024. [33] H. Wu, H. Luo, Y. Ma, J. Wang, and M. Long, “RoPINN: Region optimized physics-informed neural networks,” in Adv. Neural Inf. Process. Syst., vol. 37, p. 110494–110532, 2024, doi: 10.52202/079017-3508. [34] Z. Hao et al., “PINNacle: A comprehensive benchmark of physics-informed neural networks for solving PDEs,” in Adv. Neural Inf. Process. Syst., vol. 37, p. 76721–76774, 2024. [35] Z. Long, Y. Lu, X. Ma, and B. Dong, “PDE-Net: Learning PDEs from data,” in Proc. 35th Int. Conf. Mach. Learn. (ICML), vol. 80, p. 3208–3216, 2018. [36] P. Lippe, B. S. Veeling, P. Perdikaris, R. E. Turner, and J. Brandstetter, “PDE-Refiner: Achieving accurate long rollouts with neural PDE solvers,” in Adv. Neural Inf. Process. Syst., vol. 36, p. 67398–67433, 2023, doi: 10.52202/075280-2946. [37] F. Koehler, S. Niedermayr, R. Westermann, and N. Thuerey, “APEBench: A benchmark for autoregressive neural emulators of PDEs,” in Adv. Neural Inf. Process. Syst., vol. 37, p. 120252–120310, 2024, doi: 10.52202/079017-3822. [38] Z. Ye, C. S. Zhang, and W. Wang, “Recurrent neural operators: Stable long-term PDE prediction,” arXiv preprint arXiv:2505.20721, 2025. [39] C. Rackauckas et al., “Universal differential equations for scientific machine learning,” arXiv preprint arXiv:2001.04385, 2020. [40] J. Zhuang, D. Kochkov, Y. Bar-Sinai, M. P. Brenner, and S. Hoyer, “Learned discretizations for passive scalar advection in a two-dimensional turbulent flow,” Phys. Rev. Fluids, vol. 6, no. 6, Art. no. 064605, 2021, doi: 10.1103/PhysRevFluids.6.064605. [41] V. Shankar, D. Chakraborty, V. Viswanathan, and R. Maulik, “Differentiable turbulence: Closure as a partial differential equation constrained optimization,” Phys. Rev. Fluids, vol. 10, no. 2, Art. no. 024605, 2025, doi: 10.1103/PhysRevFluids.10.024605. [42] S. Huang, W. Feng, C. Tang, Z. He, C. Yu, and J. Lv, “Partial differential equations meet deep neural networks: A survey,” IEEE Trans. Neural Netw. Learn. Syst., vol. 36, no. 8, p. 13649–13669, Aug. 2025, doi: 10.1109/TNNLS.2025.3545967. [43] P. Bie, N. Song, N. Zhang, J. Nie, M. Ye, and X. Liang, “Space–frequency cross-attention node feature optimization graph neural operator for partial differential equations,” IEEE Trans. Neural Netw. Learn. Syst., vol. 37, no. 6, p. 2983–2993, Jun. 2026, doi: 10.1109/TNNLS.2025.3643632. [44] P. Wei and H.-X. Li, “Spatiotemporal transformation-based neural network with interpretable structure for modeling distributed parameter systems,” IEEE Trans. Neural Netw. Learn. Syst., vol. 36, no. 1, p. 729–737, Jan. 2025, doi: 10.1109/TNNLS.2023.3334764. [45] L. Menicali, D. H. Richter, and S. Castruccio, “Bayesian neural networks with physics-informed priors with application to boundary layer velocity,” IEEE Trans. Neural Netw. Learn. Syst., vol. 36, no. 10, p. 18048–18061, Oct. 2025, doi: 10.1109/TNNLS.2025.3577508. [46] B. Hillebrecht and B. Unger, “Rigorous a posteriori error bounds for PDE-defined PINNs,” IEEE Trans. Neural Netw. Learn. Syst., vol. 36, no. 1, p. 1583–1593, Jan. 2025, doi: 10.1109/TNNLS.2023.3335837.