Fast training of accurate physics-informed neural networks without gradient descent
Authors: Chinmay Datar
Taniya Kapoor
Abhishek Chandra
Qing Sun
Erik Lien Bolager
Iryna Burak
Anna Veselovska
Massimo Fornasier
Felix Dietrich
Venue: ICLR 2026 Oral (2026)
Keywords: physics-informed neural networks
extreme learning machines
random features
partial differential equations
optimization
training
causality
neural PDE solvers
optimization
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: Our approach - Frozen-PINNs addresses longstanding training and accuracy bottlenecks of Physics-Informed Neural Networks (PINNs) and makes PINNs highly realize high-precision, temporal causality, and extremely fast training.

Abstract: Solving time-dependent Partial Differential Equations (PDEs) is one of the most critical problems in computational science. While Physics-Informed Neural Networks (PINNs) offer a promising framework for approximating PDE solutions, their accuracy and training speed are limited by two core barriers: gradient-descent-based iterative optimization over complex loss landscapes and non-causal treatment of time as an extra spatial dimension. We present Frozen-PINN, a novel PINN based on the principle of space-time separation that leverages random features instead of training with gradient descent, and incorporates temporal causality by construction. On nine PDE benchmarks, including challenges like extreme advection speeds, shocks, and high-dimensionality, Frozen-PINNs achieve superior training efficiency and accuracy over state-of-the-art PINNs, often by several orders of magnitude. Our work addresses longstanding training and accuracy bottlenecks of PINNs, delivering quickly trainable, highly accurate, and inherently causal PDE solvers, a combination that prior methods could not realize. Our approach challenges the reliance of PINNs on stochastic gradient-descent-based methods and specialized hardware, leading to a paradigm shift in PINN training and providing a challenging benchmark for the community.

Published as a conference paper at ICLR 2026
FAST TRAINING OF ACCURATE PHYSICS-INFORMED NEURAL NETWORKS WITHOUT GRADIENT DESCENT
Chinmay Datar Taniya Kapoor Abhishek Chandra Qing Sun 1.3 Erik Bolager!
Iryna Burak 1.3 Anna Veselovska', Massimo Fornasier Felix Dietrich 1.3,,
TTechiiaal University of Munich (TUM)
2TUM Institute for Advanced Study
Munich Center for Machine Learning
Wageningen University & Research
6 KTH Royal Institute of Technology
6 Munich Data Science Institute
ABSTRACT
Solving time-dependent Partial Differential Equations (PDEs) is one of the most critical problems in computational science. While Physics-Informed Neural Net- works (PINNs) offer a promising framework for approximating PDE solutions, their accuracy and training speed are limited by two core barriers: gradient-descent based iterative optimization over complex loss landscapes and non-causal treatment of time as an extra spatial dimension. We present Frozen- PINN, a novel PINN based on the principle of space-time separation that leverages random features instead of training with gradient descent, and incorporates temporal causality by construction. On nine PDE benchmarks, including challenges like extreme ad- vection speeds, shocks, and high-dimensionality. Frozen-PINNs achieve superior training efficiency and accuracy over state-of-the-ar PINNs, often by several orders of magnitude. Our work addresses longstanding training and accuracy bottlenecks of PINNs, delivering quickly trainable, highly accurate, and inherently causal
PDE solvers, a combination that prior methods could not realize. Our approach challenges the reliance of PINNS on stochastic gradient-descent- based methods and specialized hardware, leading to paradigm shift in PINN training and providing a challenging benchmark for the community.
INTRODUCTION
Partial Differential Equations (PDEs) provide a unifying framework for modeling complex dynamical systems across physics, biology, and engineering, yet developing efficient methods to solve them remains a longstanding challenge (Farlow) 1993). Deep neural networks have recently shown significant promise for approximating solutions of PDEs because of the mesh-free construction of
basis functions, high expressivity of neural networks (Rudi & Rosasco| 2021). their ability to represent functions high dimensions
[
2020 Wu & Long [022 Han et al. 2018), and powerful software
for automatic differentiation (c.g.. Pytorch (Paszke al.
2017). TensorFlow (Abadi et al., 2015).
DeepXDE (Lu et al.) 2021b)). Earlier work on solving PDEs using neural networks (Dissanayake &
Phan-Thien 1994| Lagaris al. 1998) was recently popularized in the form of
neural networks (PINNs) (Raissi et al. 2019 Karniadakis et al. 2021 Sirignaao & Spiliopoulos 2018). PINNs incorporate physical constraints by minimizing a loss function involving the PDE, boundary condition, and initial condition residuals during training. Despite their promise, we identify two root causes limiting the performance of PINNS in terms of accuracy and training time.
1. Inherent challenges posed by the PINN optimization problem: Many studies (Wang et al. 2021 2022) show that even in very simple settings, the PINN loss is quite challenging to minimize using iterative gradient- based optimization methods leveraging the classical nuck-propagation algorithm (Rumelhart et al. 1986). Krishnapriyan et al. (202T) show that incorporating PDE-based soft constraints into the PINN loss function yields highly nontrivial loss landscape, rendering optimization particularly challenging. Wang et al. (2022) analyze PINN training dynamics via the Neural Tangent Kernel (NTK) and highlight issues with spectral bias and different convergence rates across different loss components. Rathore et al. (2024) show that differential operators in the PDE
"Corresponding author, chinmay.datar@uum tumm.d.e.

---

Published as a conference paper at ICLR 2026
residual loss induce "ill-conditionnng.. characterized by steep and shallow gradients in different directions near the optimum, complicating the optimization.
Efforts to improve PINN training, such as balancing loss terms (ao et al. 2023), effective regulariza- tion (Lu et al. 2021c Yu et al. 2022), architectural innovations (Wang et al. (024b), and improved optimizers (Müller & Zeinhofer 2023 Liu et al.] 2024), have been explored. We assert that such approaches address the symptoms rather than the root cause that makes training PINNs extremely challenging: the PINN optimization problem
dihe-dimeniinall (large number of trainable parame- ters), multi-objective (simultaneous minimization of PDE, and initial and boundary condition losses), and con-convex, with inherently conflicting loss terms (Liu et al., 2024) and further complicated by treating time as an additional dimension in space.
Non-causal treatment of time as an extra spatial dimension: The temporal structure of initial value PDEs is inherently Markovian as the solution at each subsequent time step depends solely on the solution at the preceding time step. Most PINN-based approaches fail to incorporate temporal causality explicitly, and time is treated neural bases spanning the entire space-time domain, exacerbating the optimization. Such approaches struggle to capture high-frequency temporal dynamics (Krishnapriyan et al., (021). and solving PDEs over a long time horizon, without resorting to domain decomposition techniques (Meng et al. 2020).
extra dimension in the input layer. This leads to
Previous studies have sought to enforce temporal causality by progressively penalizing residuals in time (Wang et al., 2024d), training distinct models across disjoint intervals with integral form losses within each interval (Jung et al., 2024), applying implicit dife-differencing with transfer learning to sequentially update PINNS on each interval (Li et al.) (002)). Nonetheless, such approaches difficult implement, require precise tuning of temporal windows and weight scheduling, and
remain computationally demanding (Kim & Son 2025 Li et al. 2024 Penwarden et al. (023). 3) See
Appendix A for an extended literature review and Appendix B. for a detailed discussion on PINNs.
(1) Space-lime separition
Dasis functions in space: sampled frozen
Time dependdpendent oulput layer coefficients
Total lrea:
L X Lpos
Ay Lpc
Ag Lae
MV, by
citt
Classical ODE with step size control
Impose hard constraints
Solve lenst squares problem
W0ga + b,)
Eie a aw b.)
i(x..)
(2) Decouple
and optimine each term sepurately
Space
Time
Frozen- PINN solution
Figure 1: Core ideas of Frozen-PINNs: (I) Space-time separation: For x E Rd, spatial bases Φ; = (w; x + b,) with a tanh, E Rd by E IR are sampled and frozen (shown for d = 1); output layer parameters (t) are evolved via ODEs. Each pair (øi, 6;) color-matched. (2) Loss decoupling: PDE, boundary, and initial condition losses LPDE, LBC, Lic are optimized independently.
To address the root causes of accuracy and training bottlenecks of PINNs rather than the symptoms, we investigate: How can the PINN optimization problem be simplified while enforcing temporal causality for time-dependent PDEs? We propose "Frozen-PINN"" based on space-time separation -a novel approach that simplifies the PINN optimization problem and enforces temporal causality by construction. We achieve this by: (a) sampling and freezing space-depender hidden layer parameters to reduce the dimensionality, (b) decoupling the PINN loss and optimizing each term separately, and (c) computing time-dependent output layer parameters using least squares and adaptive Ordinary Differential Equation (ODE) solvers, replacing gradient-descent-e based training (see Figure T). In Figure 8 we contrast Prozen-PINNss with classical PINNs. Our key contributions are:
1. Training algorithm: Frozen-PINNS break the longstanding training and accuracy bottle- necks PINNs, making PINNs rapidly trainable, temporally causal, and highly accurate, a combination realized for the first time, defining new sfate-of-the-art to our knowledge.
2. Extensive empirical evaluation: Across nine challenging PDE benchmarks and rigorous ablation studies, show that Frozen PPNNss achieve up to 4-5 orders of magnitude faster training than state-of-th -hr (SOTA) PINNs, attain high-preeision accuracies that are com- parable efficient besh-based methods in low dimensions, which most SOTA neural PDE
2

---

Published as a conference paper at ICLR 2026
solvers fail to match, and scale efficiently to high-dimensional problems where mesh-based solvers fail.
3. Adaptive solution-driven network parameters: We use solution data from previous time-steps to compute efficient neural network parameters. This extends previous work on random feature methods (Bolager al. (023) for self-supervised PDE learning tasks.
4. Model compression: We introduce SVD layer that reduces the number of neurons in the last hidden layer the network by up to 20 times and speeds up training up to 75 times.
2 SOLVING TIME-DEPENDENT PDES USING FROZEN-PINNS
In this section, we discuss the theoretical details Frozen-PINNS.
2.1 FROZEN-PIN ANSATZ
In this work, we consider time-dependent PDEs on domain C Rd for space dimension d with boundary an, seeking solutions u Ω x R - IR of PDEs defined by linear operators L and B that only involve derivative operators space, forcing f -> R. boundary aS - R, initial condition uo Ω + I, and nonlinear operator > for E R 0 for linear PDEs):
u (xx 1) + Lu(x, 1) + (N(u) f(x), x E Ω E [0, 7],
(1a)
where u, denotes the time derivative u, with boundary and initial conditions given by
Bu(æ,t) 9(x),
iΩ, t [0, and, u(zz 0) un(x), Z
(lb)
respectively. We parameterize the approximation of the solution to the PDE (Equation with Frozen-PINN having a single hidden layer with M neurons and activation function tanh
a(x, t) C(t)
(W + b) co()).
Here, c(t) RIXM and (t) R are time-dependent parameters, W E Mxd and b E RM1 space-dependent parameters, and C' := [c, cul E R (M++1) The activation functions stacked in = [di..... dM). where om (x) = o (w + b. Note that our approach does not require the PDE solution be separable in space and time. We next discuss how sample parameters W and b.
2.2 COMPUTING MIDDEN LAYER PARAMETERS WITHOUU GRADIENT DESCENT
We sample space-dependent hidden layer parameters in Frozen-PINNY using either ELM or SWIM. Hidden layer parameters are frozen (kept independent of time) after sampling (except Section 3.4
the Extreme Learning Machine (ELM) approach (Huang et al., 2006). the weights are sampled from Gaussian distribution, and biases are sampled from uniform distribution r. 27] for each hidden layer, where 1) hyper- parameer
ELM (Data-agnostic):
SWIM (Data-dependent): The Sample Where It Matters (SWIM) approach follows Bolager et al. (2023) and samples weights and biases using data deta-deeenddnt distribution. Each pair (Wm: bm
computed using two collocation points (1)
(0)) (1)) Ω: Wm. 81 Ta. J
bm -w (2))
$2: where 31, $2 E I depend on the activation function. In the unsupervised setting, one can choose pairs of collocation points from a uniform distribution over all possible pairs of collocation points, which the default setting in this paper, as we do not know the solution of the PDE beforehand. In the supervised setting (Section 3.4 Section 3.7 collocation pairs (2x(1) .x(2)) are sampled with density (2)
(1) (1) output is 0.5 (1) and +0.5 at (2), ensuring centers of activations tanh lie inside the domain
xx1)]l, ||. Neuron weights and biases set that the tanh
and are aligned with the direction zx()) ( (?), unlike ELM. The suitability of each of the proposed approaches depends on the true PDE solution's gradient distribution. See Appendix forr details. In Figure 2 illustrate the difference between the basis functions sampled with ELM and SWIM.
3

---

Published as a conference paper at ICLR 2026
drived sampling
Data- Data-agnostic sampling
Target function
SWIM basis functions
Gaussian (u=0, o=1)
FLM basis funclions
R(())
0.25
-2.5
2.5
Figure 2:
Sampling in Frozen-PINNs: (Left): SWIM (data-driven, places bases with steep gradients near regions with shocks) vs. (Right): ELM (data-agnostic, no control over basis placement).
2.3 SOLVING TIME-DEPENDERS PDES USING FROZEN-PINNS BY SEPARATION OF VARIABLES
We now discuss the computation of time-dependent output layer parameters c(I). We insert the ansatz (Equation (2) into the PDE Equation (1a), reformulating it as an ODE for c(I), preserving the inherent causal structure time-dependent PDEs, thereby enforeing temporal causality by design. We assemble No collocation points in X E RN, sample weights and biases of M neurons, compute hidden layer output ((X), and obtain the ODE
Ct(1)= R(X, (X) where R(X, C()) ((00))-((()())))) -XC0)//XX))01000XX00XXXX0X00)010 1J)
(3)
[
where [x(),, 1] E (M and the pseudo-inverse denoted by . The initial condition is computed via a least squares solution: C(0) = u(X, 0) (X), 1] which decouples the initial condition loss from PDE and boundary losses, simplifying the optimization problem. We compute C(t) via ODE solvers with step-size control (e.g.. RK45 (Dormand & Prince 1980), LSODA (Petzold 1983)) instead of gradient descent, and interpolate solutions at test points. See Appendix [B..22 Appendix B.2.3 for detailed derivations of PDE-to-ODE: reformulations for all PDEs considered here.
2.4 APPROACHES FOR SATISFYING BOUNDARY CONDITIONS FOR GROZEN-PINNY
We propose two different strategies to satisfy boundary conditions for Frozen-PINNs: the first utilizes a boundary-complint layer, and the second augments the reformulated ODE.
Boundary-comppe
Certain boundary conditions can be enforced via a linear map A e M, (M. == M) applied after the sampled hidden layer, forming houndary-complicnn layer
layer:
(see Figure Defining & = [AΦ, 1] and C(t) E IRI(MM++11 we rewrite Equation 33 to
(t = R(X, ((t))
where
R(X, C(t)) CC(()) N(((((() + [f()))
Boundary conditions defined by and g determine the construction of A; see Appendix B.2.4 for details. With a decoupling the PINN loss so that the ODE solver minimizes only the PDE residual. The rationale for outer basis functions is discussed in Appendix B.2.I
layer, boundary conditions are satisfied by construction, fully
Augmented ODE:
This strategy eliminates the need for a boundary-compliant layer by augmenting the ODE with correction term enforcing boundary conditions. For Dirichlet boundary condition u(z)= = (z)), we add û,(2)) = -x =(ū(x) - g(x)) for E 0Ω and solve the augmented system:
(t) [R(X. C(t), -(())))))(((()))))))))) (xx)) 9(X6)
(A(XX,
(5)
ERIx(N++M, (,)
(R(Net Ng)
where x > a fixed parameter, X are the N. collocation points and X,, E Nsxd collection of N₂ points on the boundary 8Q. For consistency of notation, we set A = / in Equation [ when using the augmented ODE. In practice, we skip the boundary-complie The intuition behind this technique is that the augmented ODE (Equation (5)) corrects the solution by steering axx t) toward g(z) for x E A at rate n(ū - g), with Æ = 10" as a default value. We empirically investigate the effect of K on the boundary loss and the time solution (see Figure [5)). This still partially decouples the PINN loss, with the initial condition treated separately. Depending on the PDE, domain, and boundary type, either strategy can be applied (see Appendix B.2.5
layer if we adopt this approach.
4

---

Published as a conference paper at ICLR 2026
2.5 SVD LAYER
As the last step in the Frozen-PINN architecture, we add a linear layer to reduce the stiffness of the associated ODE (Equation (4) and the size of the ODE system. To achieve this, we propose orthogonalizing the basis functions using an SVD layer. We compute a truncated singular value decomposition of ( E &M to obtain matrices V,, Ey, and Up with r S Mo such that V,E.U, = (X) O(E++1). We then define A, := V, A and use it instead of the matrix A and C(t) E Ix(++1) This ensures (X) are orthogonal functions on the data X, and the matrix (X) has bounded condition number. The SVD layer accelerates computation by up 75 times while reducing the ODE system dimension 20 times, as validated by an extensive ablation study (see Appendix C.. Figure 3 visualizes the complete Frozen-PINN architecture.
2.6 SUMMARY OF THE TRAINING ALGORITHM FOR PROZEN-PINNS
We summarize training process in Algorithm where CSVD is the SVD threshold that governs the SVD-layer width. See Appendix for for additional methodological details, and Appendix [B2.]] for extended discussion PINN vs. Frozen-PINN training, comparison between sampling strategies, influence of random sampling, rationale for outer bases, and the Kolmogorov n-width barrier.
Boundary- compliant layer
Algorithm Frozen-PINN training algorithm
ELM/
Input layer
SVD layer
Output layer
Input: PDE (Equation [ D. test grid points Xast Tice
Output: PDE Solution the test grid points ((X Tua)
Parameters: No. M,, M, E N, SSVD E
Sample N₂ collocation points: E R. Xu
V. (XX
2: Construct hidden layer params (Mom bm) Mr (SWIM/LM) Section 2.2
3: Compute hidden layer output (1 X.N,
(IXN.)
(Nexd)
((X)
A 440X)
(rxN2)
4: Construct boundary-compllant layer: (+x( E RM
Section 2.4
5: Compute truncated SVD: V.E.UU = ( and SVD layer output
[M, (M,)N,)
V," (x) (XX)
6: Compute neural hases: A. (x) == (A, 1) E R(+
Architecture of Frozen- PINNs trained with a gradient- descent-free training algorithm.
Figure 3:
7: Initialize outpuz-laser params (least-squares); C(O) u(X T (X)
8: Solve ODE for C(t) 3 (++1) using 4 Ar.
Equation
9: Evaluate 6(Xue. Tai) C(Tue)&A (Xuc))
Equation
3 EMPIRICAL RESULTS
In this section, with a comprehensive empirical study across nine challenging low- and high- dimensional PDE benchmarks, we demonstrate that Poozen-PINNs consistently outperform existing state-of-the-are neural PDE solvers with orders-of-magnitude faster training in all cases and higher accuracy in almost all cases without requiring specialized hardware like GPUs. Moreover, our work includes rigorous evaluation against the classical SOTA approaches like IGA-FEM (see (ppendix(B.3) (Hughes et al., 2005 Cottrell et al. 2006 2009) or FEM for low-dimensional PDEs, bridging a gap not suffciintlyy addressed in the literature between neural and mesh-based solvers.
Appendix C contains details of the PDEs, important ablation studies for our experiments (for the SVD layer and the width of the network), metrics used for comparison, train and test data, software and hardware environments, the absolute error plots on test points, and elaborate explanations of results. Figure 12 visually summarizes all the PDE benchmarks used for evaluation, identifies the specific challenges posed by each PDE, and shows true solutions. We perform all experiments with three seeds and report the mean and standard deviation. The code reproduce experiments from the paper, and a refactored version that will be actively maintained are available at:
https://gitlabs com/ fd heasearchsswimpdee http://w...aa.ooooeeeeeeeerrraas
To ensure fair comparisons, we follow the two rules outlined by McGreivy & Hakim 2024): (i) we benchmark at (almost) equal accuracy, defining low-precision (1e-2 to 1e-4) and preciiiiiiioon (le-5 le-10) regimes, configuring Frozen-PINNs to marginally outperform the best PINN baselines the low-pree preciiiioonn regime and aligning FEM/IGA -EM fidelity with Frozen-PINNs in the high-precision regime; (ii) compare against efficient numerical methods, including SOTA
5

---

Published as a conference paper at ICLR 2026
IGA-FEM or classical FEM for low-dimensional PDEs, while highlighting neural solvers' scalability in high-dimensiona benchmarks where FEM and IGA-FEM suffer from the curse of dimensionality.
3.1 HIGH ADVECTION SPEEDS, FAST CONVERGENCE, AND LONG-TIME SIMULATION
We benchmark the linear advection equation to demonstrate how Frozen- PINNs resolve three im- portant well-known challenges for PINNs: (1) handling high advection speeds (Krishnapriyan et al. 202Ip, (2) achieving fast convergence with increasing width (Cuomo et al. 2022). and (3) long-time time
simulations (Lippe et al. 2024
Kapoor et al. 2024a). We describe all details in Appendix C.I
High advection speeds: We solve the advection equation for increasing advection coefficients, denoted by B3. (igure[[(Leff) shows that approaches using basis functions in the entire spatiotemporal domain, such as PINNs, ELM, and SWIM, completely fail as the flow velocity 3 increases beyond 40, In contrast, Frozen-PINN, can accurately solve the PDE, even for extremely high values of 3 (as high as 10t) with relative L errors less than 10-4 Table 1 shows that for B 40, Frozen-PINNs train 45 to 533 times faster than other alternatives at similar accuracy in the low-precision regime. With the exception of Frozen-PINNA none the neural PDE solvers evaluated here attain high sigh-precision accuracy. Frozen-PINNs outperform existing neural PDE solvers by over six orders of magnitude in accuracy and approach the fidelity of IGA-FEM, which unsurprisingly is the most accurate solver.
Fast convergence (error decay with hidden layer width): For low value of advection coefficient = 10, Figure (Middle) shows that errors with classical PINNs do not decay quickly with width, primarily due the difficulties training. In contrast, the relative L error decays exponentially with hidden layer width for Mrozen-PINNs, ultimately plateauing at a value more than four orders of magnitude smaller than that obtained with PINNs.
Long-time simulation: Neural PDE solvers employing joint space-time basis functions, like vanilla PINNs, encounter substantial challenges in accurately approximating dynamics over extended time spans. Here, we consider the advection equation with the advection coefficient / = 1. As shown in Figure F (Right), Frozen-PINNe can simulate the advection equation for 1000 seconds with relative L error under 0.001% in just 0.94 seconds.
PNN
LIM
PrzzenFllleeih Filsssi
PINN
Gausel PAN
SWIM
Relatice crro
10
10-4
10-7
Reaite
30
Frozer-PINN-elm
10-4
Detetive
10-4
10-9
Frozer-PIMM swim
30
10-4
250
500
750
1000
Time
10-2
100
101
107
107
Number functions
Figure 4: Illustration of experimental results for the advection equation: (Leftt): high advection speeds effect advection coefficient / on the test error for different PDE solvers, (Middle): fast convergence with B = 10, Frozen PINNs achieve exponential decay in test error as indicated by the reference dotted line, while standard PINNs display plateaued error decay despite increasing number of basis functions (hidden layer size), (Right): long time simulation Slow error growth with time.
3.2 HIGHER-ORDER DERIVATIVES IN SPACE AND TIME
We consider two variants of the Euler-Bernoulli beam equation -classical Buler-Bernould beam equation and its extension with a Winkler foundation. See Appendix C.2 for details. The main challenge posed by both PDEs for PINNs is the higher-order differential terms (fourth- and second- order derivatives in space and time, respectively). Prozen-PINNM eliminate expensive evaluation of higher-order derivatives via backpropagation, cutting training cost by four orders of magnitude in the low- precision regime, while achieving (GA-FEM-lee accuracy that more than six orders of magnitude accurate compared to other SOTA PINN benchmarks considered here (see Table T).

---

Published as a conference paper at ICLR 2026
3.3 MULTI-SCALE SOLUTIONS
To demonstrate the capability of our method to solve PDEs with multi-scale solutions, we consider a Wave equation benchmark (Hao et al. 2024). We examine two settings: one with two distinct frequencies and another with three well-separated frequencies, which increases the spatial complexity and significantly broadens the range of scales in the solution. For the two-frequency setup, we compare the performance against prior PINN baselines in Table [ and observe that CPU-trained Grozen-PINNY achieve 625 to 5500 times faster training than GPU-trained competing PINN variants, while simultaneously being four to five orders of magnitude more accurate. Frozen-PINNs also solve the wave equation in the three- frreeeecy scenario (illustrated in Figure 19) extremely quickly and with high precision, reinforcing their potential for solving PDEs with complex, multiscale dynamics. Additional implementaiion details and extended results are provided in Appendix C.3
3.4 NON-LINEARITY AND SHOCKS
In this example, we highlight how using pairs of data points to sample neural basis functions using the SWIM algorithm can be leveraged resolve locally steep gradients in the solution of the non-linear viscous Burgers' equation, as shown in Figure Saareee[Lcee (Left). See Appendix C.4 for details.
Prozen-PINN-see swim creates numerous basic functions with steep gradients, accurately placing them near the location of the shock, leveraging the SWIM algorithm and solutions from previous time-steps to fit neural basis functions, given enough collocation points in the domain's center (see Figure 200 (Left)). To concentrate collocation points near the shock in the domain's 's center, we resample them periodically after a set number of time steps. guided by probability distribution that leverages the gradient of the approximate solution (see Figure 5a (Top)). At the resampling time t, E T). approximate the probability density p(x) ~ Vū(c. (. which we then use to re-sample collocaioo points illustrated in Figure 5a (Bottom), placing more collocation points near the shock region.
A. Initial PDE solution (Varies only zz))
D. SWIM basis function (without projection)
C. SWIM besis function (with projection)
Probability distrbutton
10
10-5
sp (2)
(3)
a()
0.0
D.5
Sampled points
Unariented bas fundtion gradient ution gradeet
Wall cor.-reeeddddeeeee basts Function Solution gradien
Solution gratlent
0
(b) Illustration of embedding directional information to orient Frozen-PINN- swim basis functions along the gradient of the initial condition: (Left): A toy initial condition up R - R varying a single direction, (Middle); randomly selected pair of points (zx(1) (²) leads to SWTM hasis function misaligned with the gradient of the initial solution, (Right): A projected pair of points (a()), (()) yields basis function aligned with the gradient of up.
-1
0.0
(a) (Top): Probability distribution, (Bottom): Sampled collocation points.
Figure 5: Constructing useful PiNcenn--INN-swiin bases. (Left): shock-aware sampling (Burgers,
(Section 3.4b) and (Right): direction-aware bases rreactionddiffsiion Section 3.7).
As shown in Table 1 Frozen-PINNs achieve 46 2945 times speedups in training time over other PINN variants in the low-precision regime. Remarkably, even in the figh-precision regime, CPU- trained Frozen-PINNs remain 203 to 535 times faster than state-of-the-
GPU-trained PINNS at comparable accuracy. While optimizers like SSBroyden (Kiyani et al. 2025) can offer higher accuracy, they are extremely slow, resource-intensiee and difficult to implement. Furthermore, Frozen-PINN-swin basis functions handle shocks significantly better than Fourier or Chebyshev bases used in classical spectral methods (see Figure 22 Figure 23 Appendix C.4.1).
3.5 NON-LINEARITY AND COMPLICATED DOMAIN GEOMETRY
In this example, we consider a non- non-linear diffusion equation on a complicated domain geometry. See Appendix C.5 for details. For besh-based methods, meshing can be resource-ittensve and technically
7

---

Published as a conference paper at ICLR 2026
demanding (see Figure 24). unlike neural PDE solvers. As shown in Table Frozen-PINNs 145 to 456 times faster than PINNs and 4.83 times faster than FEM at comparable low-precision accuracy, and can achieve over 1000 times better accuracy than other PINNs. Notably, Frozen-PINNs require only 350 basis functions versus around 2000 finite elements in FEM for similar accuracy (see Table [8), mainly due the global support of neural bases. For fairness, the FEM grid points reused as collocation points for minimizing the PDE residual in Frozen-PINNs.
3.6 CHAOS AND STRONG NON-LINEARITY
We tackle the highly nonlinear Kuramoto Kuramoto-Sivashinskky equation, which models laminar flame-from instabilities that exhibit spatiotemporal chaos. As shown in Figure 6 our Frozen-PIIN captures the characteristic chaotic pattern over a long-time horizon E (0, 5). with an average training time of only 6.9 seconds CPU (averaged over 5 seeds). Further experimental details are provided in Appendix C.7 Since chaotic dynamics amplify small numerical differences, trajectory-level errors not meaningful, and assess performance based on the qualitative spatiotemporal pattern..
3.7
PDES WITH LOW-DIMENSIONAL SOLUTION MANIFOLDS
In this benchmark ZZng et al. 2020). we solve a five-dimensional non-linear reaction-diffusion equation, where the solution only changes two dimensions that are a priori unknown. We construct SWIM basis functions aligned with the two intrinsic dimensions of variation, directly embedding directional information unlike PINNS and ELMs. by using spatial coordinates projected onto the gradient of the initial solution to sample SWIM basis functions, as shown in Figure [b Appendix C.6 for further details.
See
Table 1 shows that Frozen-PINN-seee PINN-swim is over 3400 times faster than other PINNs at comparable low-precision accuracy.
the only method reach the high-precision regime, achieving 2-3 orders of magnitude higher accuracy than other PINN variants and weak adversarial networks (Zang) iet al.) 2020). These results confirm that explicitly embedding informative basis functions yields far more efficient and accurate models than relying on iterative optimization to learn them implicitly.
3.8 HIGI--DIMENSIONALII
High--dimensonnl PDEs, such as the 100-dimensional heat equation, are computationally prohibitive for grid-based methods, which require more than 1030 grid points, considering only points per dimension. The following examples demonstrate Frozen-PINNs' ability to solve such PDEs efficiently established benchmarks: one introduced in Wang & Dong (2024). which addresses the heat equation in up to 10 dimensions on unit hypercube, and another introduced [He al. (2023). which focuses on 100-dimensional variant of the heat equation on a unit ball. We discuss all details in Appendix C.8
and accurately. We evaluate our approach
Miizen-PINNNelmm is consistently 10-1000 times more accurate than classical PINNs for up to 100- dimensional PDEs (esuuree(top) with error decaying rapidly with network width until saturation Figure 7 (bottom). For the 10-d heat equation, Frozen-PINi PINN-elm trains 100 1000 times faster than other PINNs while achieving higher accuracy. For the 100-d heat equation, CPU-trained Frozen-PINN, remain hundreds of times faster than GPU-trained PINNs while delivering an order- of-magnitude better accuracy (Table underscoring both their computational efficiency and high accuracy. Table 2 summarizes the advantages our algorithm over classical mesh-based based and physics-informed methods based on iterative gradient-descent-ban methods.
4 CONCLUSION
Frozen-PINNS directly address the longstanding training and accuracy bottlenecks of PINNs by fundamentally simplifying the optimization problem and enforcing temporal causality by construc- tion, leveraging the idea of space-time separation. Our extensive empirical analysis reveals that Frozen-PINNe consistently realize extremely fast training and high precision (often several orders of magnitude better than SOTA PINNs), and preserve temporal causality broad range of PDEs involving challenges such as extreme flow velocities, long-time simulation, higher-order spatial and temporal derivatives, complicated spatial domains, non-linearities, shocks, and high-dimensionality,
8

---

Published as a conference paper at ICLR 2026
Table 1: Summary of empirical results on eight PDE benchmarks, including results from prior works: dashes denote training times not reported in prior works; Training times labeled with were obtained using GPUs; thus, CPU-based training. as with Frozen-PINNs, would lead to substantially larger values. For each PDE, solvers above/below the horizontal line correspond to low-/high- precision regimes. Normalized training times relative to Frozen-PINNs are computed as the ratio of each method's training time to that of Frozen-PINNS and are computed at similar accuracy.
PDE benchmark
Method
Training time
Normalized training time
Relative L error
Advection (33 40)
PINN (Adam)
Fail for 8=40
SWTM
Fail for 8000
ELM
Fail for 8=40
Causal
357.63
PINN (L- BFGG)
30.5
45.5
6.92e 1 1: 2.96e-2
Krichnuprivan
PINN (seq2seq, L-BFGS)
Krismapryaa
2021
PINN (Curriculam training, L-BFGS)
5.336-2
Frozen-PINN (our)
0.67
4.19e-3 di 2,97e-3
Prozen-PiK N-swis (our) o
0.7
R.42-9 E 1.126-8"
Mesh-based method (IGA)
0.07
D.1
1.17e-10
Fuler-Bernone Bernonlli (classical)
PINN (Adam)
4209.82
84196
3.95e- 4: 1.796-2
Kapoor et al 2023
PINN (L-BFGS)
2303.71
46074
4.21o-3 + 9.56c-
Frozen-PINN-eeln (our)
0.05
2.82e4 2.15e-4
Mesh-bused melhoo (IGA)
0,94
0.13
4.2Te-7
Frozen-PIN Frozen-PINN-elm (our)
6.90
9.33e-9 1: 4.36e-9
Fuler- Bernoulli (Winkler)
PINN (L-BFGS)
1858
37 370
5.33e40
Kapoor et al
2024b
Adaptive PINN
3807 89
76140
5.32c+0
Kapoor
Self-adaptive PINN
4042.57
Kapoor et at 0144
Wavelei PINN
4764.25
95280
4.38e+0
Kapoor et al
M1245
Causal PINN
1873
37460*
3.00e-2
Frozen- P-NNeemm (our)
0.05
4.19e-
Frozaa-PINNN-siiaa (Gar)
3.41
T.47-7 E 1206-7
Mesh-based method (IGA)
1.88
0.44
2.70e-8
Wave Hao
2024
PINN (FBPINN)
3090
5517.9*
5.9ee1 1 4,74e-2
Han et al 2134
PINN (L-BFGS)
350
625
5.88 1 1 9.63e-2
Han et ml M134
PINN (gPINN)
775
1383.9
5.56e + 1.67e-2
Ilu et al
X024
PINN (NTK)
840
1500
9.79e 1 7.72e-3
Frozen-PiNN-emm (our)
0.56
I8īe- Σ 101c-6
Burgers
Causal PINN
1531.79
2945:75
1660 8.97e-3
PINN (L.BFGS)
275.2
$29.2
3.88e 261ee3
Kiyani ct
2025
PINN (BFGS with trust region)
24
46.1
1.l-3
Mozen-PINN-swim (our)
0.52
1.00-3 1.13e-3
Chen et al
2024b
PINN (residual-based sttention)
8.22-4 + 2.33e-
McClenny
Hraga-Neto
2023
Self-adaptive PINN
4.80-4 + lc-4
PINN (balanced residual decay rate)
1338-4 0.85c-4
Kiyani
2025
PINN (RAdam BFGS)
1070
203
6e-6
Urbán et al
2025
PINN (SSBroyden)
2.9--6 + 0.46-
Frozen-PINN. Piwen. (our)
5.25
2.27e-7 1: 6.89e-8
Mesh-bbssed method (IGA)
76.32
14,5
1.126-7
Kiyani ct al. 2025
PINN (Adam SSBroyden)
2812
535
1.62-8
Nonlinear diffusion
PINN (Adam)
81,36
145.2
2.09-22 +: 14e-3
PINN (L-BFGS)
255.9
456.9
1.22e 1: 2.38e
bessh-bsedd metiod (FEM)
2.71
4.83
266e-3
Prozen-PIN--ee (our)
0.56
1
2.60e-3 + 1.61e-3
Prozea-PIN--swiie (our)
423
2.00e-6 + 1.99-66
5-d Reaction diffusion
PINN (Adam)
171.33
3428.6
3.40e-1 + 1.79c-2
PINN (L-BFGS)
183.38
3667.6
3.336-2 1.54c-2
Zang et al 2020
Weak Adversarial Network
2.8e-2
Frozen- Frozen-PINN-swiin P-wN-swim (our)
0.05
1.07e2 4.52e-4
Frozen-PINN.o PINN-swim (our)
12.43
9.99c-5 1 6.21e-9
10-d heat
PINN (Adam)
1002.49
3037.8
1.68-1 3.21e-2
Wang Dong 2024 2024
PINN (L-BFGS)
189.6
574.3
6.06e- + 1.00c-4
benchmark exteded to
FraeenPPINNeelm (our)
0.33
4.35e 5.91e-5
Frozen.PiNN.eim (our)
T68.6
3.6633 2.1e.5
100-d head He et
2023
PINN (Adam)
141
1084.6
0.60e-2
He ct 2023
PINN (no stacked-baokproppgttonn
49.8
383.1
0.63e-2
PINN (Adam+L-BFGS)
26.25*
201.9 1
4.980-3 2.96c4
Frozen-PPP P-NN- (our)
0.13
4.12e-4 +: 1.70e-5

---

Published as a conference paper at ICLR 2026
without requiring specialized hardware like GPUs. Prozen-PINNs maintain high precision over long time spans and capture
temporal dynamics where prior neural PDE solvers fail.
Within the scope of the empirical study in this work, in low dimensions, Frozen-PINNs match the accuracy of classical mesh- solvers while retaining advantages such as mesh-free free basis functions, ease of implementation, the ability to handle complex domains, spectral convergence for PDEs with smooth solutions, and scalability for high-dimensional PDEs where mesh-based approaches struggle.
Kuramoto Swoshinsky countion
PINN
Frozen-PINN-elle
Relative to eerar
10
10-1
10
Figure 6: Frozen-PINN solution of the strongly non- linear and chaotic -iraaoto-Sivashinsky equation.
10-1
100
Dimension on of the PDE
PDE wetting
IGA-FEM/
PINNs
Froee-PINN.
FEM
10
d=
Solutions with shocks
/ (SWIM)
Complex domains
Easy
10
5
High dimensionality
(Cab)
Relatve
d
Performanc//eatures
10
Accumacy/Prcciioo
High Fast
Often low Slow (triining)
lligh
d 10
Spead
10a
T0
100
Temporal causality
x (sofi constraint)
Width (before SVD layer)
4 100
Table 2: Comparison of Frozen- with mesh-based FEM and classical PINNs in different problem settings presented in this paper: The comparison is grounded in results reported in Section [3] for the PDEs and solvers studied. / denotes compatibility, and x denotes either incompatibility or the need for substantial modifications. Curse of Dimensionality is abbreviated as CoD)
Figure 7: High-dimensional heat equation:
(Top): comparison of test errors for vary- ing PDE dimensions (different hatch pat- terns indicate different benchmarks), (Bot- tom): fast decay of test error with network width (dashed:
PINN-swim, solid:
Frozen-PINN-elme
Limitations and future work: Our method assumes knowledge of the PDE, but its speed makes it well-suited for inverse problems via fast forward solves. While Frozen-PINN: efficiently deal with extreme temporal complexity, as shown in the advection equation with extreme flow velocities, dealing with spatial complexity encountered while solving PDEs like Navier-Stokes is exciting next step, where one could leverage domain decomposition to deal with the added complexity (Moseley al., 2023 Howard al. 2024). Finally, universal approximation properties concerning specific PDE settings and understanding the role of re-sampling network parameters in overcoming the Kolmogorov n-width barrier (Peherstorfer 2022) are some the most challenging, yet important theoretical open areas of investigation, beyond the scope of this paper.
Frozen-PINNS take a decisive step toward practical neural PDE solvers through lightweight opti- mization process and extremely fast training without GPUs, promoting low-carbon Al development (Verdecchia et al.) 2023), advancing sfate-of-hheaar performance, and establishing a formidable benchmark for the community to build upon in advancing rapid and accurate neural PDE solvers.
Reproducibility statement:
The source code used reproduce the experimental results, along with comprehensive reproducibility instructions, is included in the supplementary material and publicly available as an open-source repository. All experiments are run with multiple seeds, and the corresponding seed values are stored in the repository to ensure reproducibility.
Ethics statement: Neural networks are inherently dual-use technologies, and ethical considerations are essential for any new machine learning approach. Frozen-PINNs are grounded in classical scientific computing principles, which offer well-understood behavior and interpretability. By bridging neural PDE solvers with classical numerical methods, framework enables clearer analysis of robustness, failure modes, and reproducibility. We helieve this transparency reduces the risk of misuse and enhances controllability, making Frozen-PINNA safe and interpretable. Thus, believe that the benefits of our approach far outweigh the potential downsides of misuse because a system that is better understood can also be controlled more straightffrwardly.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
We are grateful Tim Bürchner for providing an initial version of the IGA code and Jana Huhne for providing an initial version of the Frozen-PINN code. We sincerely thank Dinesh Parthasarathy, Ana Cukarska, and Adwait Datar for the detailed feedback on the manuscript. Their careful review and suggestions greatly strengthened the manuscript. We further express our appreciation to Wil Schilders for engaging discussions and valuable perspectives. F.D., I.B., and E.B. acknowledge funding by the German Research Foundation- project 468830823, and association with DFG-SPP-229 C.D. is partially funded by the Institute for Advanced Study (IAS) at the Technical University of Munich. F.D. and Q.S. are supported by the TUM Georg Nemetschek Institute Artificial Intelligence for the Built World. F.D. and C.D. acknowledge Sølve Eidnes and Martine Mahlum for a great workshop on physics in machine learning in Oslo, 2024, which provided with good feedback the work.
REFERENCES
Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng. TensorFlow: Large-scale machine learning on heterogencous systems, 2015. Software available from tensorflow.org.
Olah,
Joubine Aghili, Joy Zialesi Atokple, Marie Billaud-Friess, Guillaume Garnier, Olga Mula, and Norbert Tognon. A Dynamical Neural Galerkin Scheme for Filtering Problems, January 2024.
Igor A. Baratta, Joseph P. Dean, Jørgen S. Dokken, Michal Habera, Jack S. Hale, Chris N. Richardson, Marie E. Rognes, Matthew W. Scroggs, Nathan Sime, and Garth Wells. DOLFINx: the next generation FEniCS problem solving environment. preprint, 2023.
Cea Basdevant, M Deville, P Haldenwang, JM Lacroix, J Ouazzani, R Peyret, Paolo Orlandi, and AT0612 Patera. Spectral and finite difference solutions of the burgers equation. Computers & fluids, 14(1):23-41, 1986.
Jules Berman and Benjamin Peherstorfer. Randomized sparse neural galerkin schemes for solving evolution equations with deep networks. Advances in Neural Information Processing Systems, 36, 2024.
Jules Berman, Paul Schwerdtner, and Benjamin Peherstorfer. Neural Galerkin schemes for sequential- in-time solving of partial differential equations with deep networks. In Handbook of Numerical Analysis, volume 25, pp. 389-418. Elsevier, 2024. ISBN 978-0-443--2998-77 doi: 10. 10.1016/bs.hma 2024.05.006.
Dimitri Beskos. Boundary element methods in dynamic analysis. Applied Mechanics Reviews, 40 (1):1-23, 1987.
Erik L Bolager, Iryna Burak, Chinmay Datar, Qing Sun, and Felix Dietrich. Sampling weights of deep neural networks. In Advances in Neural Information Processing Systems, volume 36, pp. Curran Associates, Inc., 2023.
Mayur P Bonkile, Ashish Awasthi, Lakshmi. Vijitha Mukundan, and VS Aswin A systematic literature review of burgers' equation with recent advances. Pramana, 90:1-21, 2018.
N. F. Britton. Reaction-diftusio equations and their applications biology. Academic Press, London, 1986. ISBN 9780121351403.
Joan Bruna, Benjamin Peherstorfer, and Eric Vanden-Eijnden. Neural Galerkin schemes with active learning for high-dimensional evolution equations. Journal of Computational Physics, 496:1 12588, January 2024. ISSN 09919991. doi: 11006//jip.222231112588
Hans-Joachin Bungartz and Michael Griebel. Sparse grids. Acta numerica, 13:147-269, 2004.
11

---

Published as a conference paper at ICLR 2026
Francesco Calabro, Gianluca Fabiani, and Constantinos Siettos. Extreme learning machine collocation for the numerical solution of elliptic PDEs with sharp gradients. Computer Methods in Applied Mechanics and Engineering, 387:1 14188, December 2021. ISSN 00457825, doi: 10.1016/j.cma
2021.1 11488..
Jingrun Chen, Weinan and Yifei Sun. Optimization of random feature method in the high-precision regime. Communications on Applied Mathematics and Computation, pp. 1-28, 2024a.
Peter Yichen Chen, Jinxu Xiang, Dong Hcon Cho, Yue Chang, G A Pershing, Henrique Teles Maia, Maurizio M Chiaramonte, Kevin Thomas Carlberg, and Eitan Grinspun. CROM: Continuous reduced-order modeling of PDEs using implicit neural representations. In The Eleventh Interna- tional Conference Learning Representations, 2023. URL
Forum?1d-FUORR1tGG00
Wen Chen, Zhuo-Jia Fu, and C.S. Chen. Recent Advances in Radial Basis Function Collocation Meth- ods. Springer Berlin Heidelberg, 2014. ISBN 9783642395727. doi: 10.
Wenqian Chen, Amanda A Howard, and Panos Stinis. Self-adaptive weights based on balanced residual decay rate for physics-informed informed neural networks and deep operator networks. arXiv preprint arXiv:2400..00 006111,,,1,,,. 2024b.
Pao-Hsiune Chiu, Jian Cheng Wong, Chinchun Ooi, My Ha Dao, and Yew-Soon Ong. Can-pinn: differentiation
A fast physics- phfsies-ifformed neural network based on method. Computer Methods in Applied Mechanics and Engineering, 395:1 14909, 2022.
Junwoo Cho, Scungtae Nam, Hyummo Yang. Seok-Bac Yun, Youngjoon Hong, and Eunbyung Park. Separable physics-infommn informed neural networks. Advances in Neural Information Processing Systems, 36, 2024.
Ronald R Coifman and Stéphane Lafon. Diffusion maps. Applied and computational harmonic analysis, 21(1):5-30, 2006.
J. Austin Cottrell, Thomas J. R. Hughes, and Yuri Bazilevs. Isogeometric Analysis: Toward Integration of CAD and FEA. Wiley Publishing, 1st edition, 2009. ISBN 0470748737.
J.A. Cottrell, A. Reali, Y. Bazilevs, and T.J.R. Hughes. Isogeometric analysis of structural vibrations. Computer Methods in Applied Mechanics and Engineering, 195(41):522) (95(41):5277-5296 2006. ISSN 0045-7825. doi: https://doi.org/10.100 www.:000.6...m22000........ John H. Argyris Memorial Issue. Part II.
M. G. COX. The Numerical Evaluation B-Splines* IMA Journal of Applied Mathematics, 10(2): 134-149, 10 1972. ISSN 0272-4960. doi: 10. 100.109//mamat/10.2.25
Salvatore Cuomo, Vincenzo Schiano Di Cola, Fabio Giampaolo, Gianluigi Rozza, Maziar Raissi, and Francesco Piccialli. Scientific machine learning through physics-informed neural networks: Where we are and what's next. Journal of Scientific Computing. 92(3):88, 2022.
Chinmay Datar, Adwait Datar, Felix Dietrich, and Wil Schilders. Systematic construction of continuous-time neural networks for linear dynamical systems. SIAM Journal on Scientific Computing, 4784):C800 2025.
Carl de Boor. On calculating with b-splines. Journal of Approximation Theory, 6(1):60-66) 1972.
ISSN 0021-9045. doi: https://doi.ons .rg//0../////ooo wttp://i..rr//111111101111000 11011///002-190555229008881
MWMG Dissanayake and Nhan Phan-Thien. serrrll-etwooo--baaaeae approximations for solving partial differential equations. communications in Numerical Methods in Engineering. 10(3): 195-201, 1994.
Suchuan Dong and Zongwei Li. Local extreme learning machines and domain decomposition for solving linear and nonlinear partial differential equations. Computer Methods in Applied Mechanics and Engineering, 387:1 14129, 2021.
12

---

Published as a conference paper at ICLR 2026
Suchuan Dong and Jielin Yang. On computing the hyperparameter of extreme learning machines: Algorithm and application to computational pdes, and comparison with classical and high-order finite elements. Journal of Computational Physics, 463:1 11290. 2022.
John R Dormand and Peter J Prince. A family of embedded runge-kutte formulae. Journal of computational und applied mathematics, 6(1):19-26, 1980.
Gideon Dresdner, Dmitrii Kochkov, Peter Christian Norgaard, Leonardo Zepeda-Nunez, Jamie Smith, Michael Brenner, and Stephan Hoyer. Learning to correct spectral methods for simulating turbulent flows. Transactions on Machine Learning Research, 2022.
Yifan Du and Tamer A Zaki. Evolutional deep neural network. Physical Review E, 104(4):045303, 2021.
Yiheng Du, Nithin Chalapathi, and Aditi Krishnapriyan. Neural spectral methods: Self-supervised learning in the spectral domain. arXiv preprint 1rXiv:2312.05225. 2023.
Vikas Dwivedi and Balaji Srinivasan. Physics informed extreme learning machine (pielm)-a rapid method for the numerical solution of partial differential equations. Neurocomputing, 391:96-118 2020.
Vikas Dwivedi, Nishant Parashar, and Balaji Srinivasan. Distributed learning machines for solving forward and inverse problems in partial differential equations. Neurocomputing, 420:299-316 2021.
Weinan E. Towards a Mathematical Understanding of Neural Network-Based Machine Learning: What We Know and What We Don't. CSIAM Transactions on Applied Mathematics, 1(4):561-615, June 2020. ISSN 2708-0560, 2708-0579,
doi: 10.4208/csiam- am. S0-20000000.
Gianluca Fabiani, Francesco Calabrò, Lucia Russo, and Constantinos Sicttos. Numerical solution and bifurcation analysis of nonlinear partial differential equations with extreme learning machines. Journal of Scientific Computing, 89(2):44, November 2021. ISSN 0885-7474, 1573-7691. doi: 1.1007//10911--21-011666
Stanley J. Farlow. Partial Differential Equations for Scientists and Engineers. Dover Publications, New York, 1993. ISBN 9780486676203.
Mare Finzi, Andres Potapczynski, Matthew Choptuik, and Andrew Gordon Wilson. A stable and scal- able method for solving initial value pdes with neural networks. arXiv preprint arXiv:2304. 14994, 2023.
Evangelos Galaris, Gianluca Fabiani, Ioannis Gallos, Ioannis Kevrekidis, and Constantinos Sicttos. Numerical Bifurcation Analysis of PDEs From Lattice Boltzmann Model Simulations: A Parsi- monious Machine Learning Approach. Journal of Scientific Computing, 92(2):34, August 2022.
ISSN 0885-7474, 1573-7691. doi: 1..1077/119915-02--018333
Christophe Geuzaine and Jean-Françoh, Remacle. Gmsh: A 3-d finite element mesh generator with built-in pre- and Engineering. 79:1309 1331, 09 2009. doi: 10. 1071/mmmm.2.59.
facilities. International Journal for Numerical Methods in
R. A. Gingold and J. Monaghan. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society, 181(3):355-389, 12 1977. ISSN 0035-8711. doi: 10. (0.109//mnrr///81....377
Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural networks. In Proceedings of the thirteenth international conference artificial intelligence and statistics, pp. 249-256. JMLR Workshop and Conference Proceedings, 2010.
David Gottlieb and Chi- Chi-Wang Shu. On the gibbs phenomenon and its resolution. SIAM review, 39(4): 644-668, 1997.
Jiequn Han, Aruuff Jentzen, and Weinan E. Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences, 115(34)850058511) 2018.
13

---

Published as a conference paper at ICLR 2026
Zhongkai Hao, Jiachen Yao, Chang Su, Hang Su, Ziao Wang, Fanzhi Lu, Zeyu Xia, Yichi Zhang, Songming Liu, Lu Lu, and Jun Zhu. Pinnacle: a comprehensive benchmark of physicc- informed neural networks for solving pdes. In Proceedings of the 38th International Conference on Neural Information Processing Systems, NIPS '24, Red Hook, NY, USA, 2024. Curran Associates Inc. ISBN 9798331314385.
Di He, Shanda Li, Wenlei Shi, Xiaotian Gao, Jia Zhang, Jiang Bian, Liwei Wang, and Tie- Yan Liu. Learning physics-informed neural networks without stacked back-propagation. In International conference on artificial intelligence and statistics, pp. 3034-3047. PMLR, 2023.
Amanda A Howard, Bruno Jacob, Sarah H Murphy, Alexander Heinlein, and Panos Stinis. Finite basis kolmogorov-arnold networks: domain decomposition for data-driven and physics-informed problems. arXiv preprint arXiv:2406. 19662, 2024.
Guang-Bin Huang, Qin- Yu Zhu, and applications. Neurocomputing, 70(1-3)4499-011, 2006.
Siew. Extreme learning machine: theory and
Xinquan Huang and Tariq Alkhalifah. Efficient thysics--iffrmed neural networks using hash encoding. Journal of Computational Physics, 501:1 12760, 2024.
T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs. Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194(39)441354419) 2005.
Tanveer ul Islam and Prasanna S Gandhi. Viscous fingering in multiport hele shaw cell for controlled shaping of fluids. Scientific reports, 7(1):16602, 2017.
Ameya D Jagtap, Kenji Kawaguchi, and George Em Karniadakis. Adaptive activation functions accelerate convergence in deep and physics-informee neural networks. Journal of Computational Physics, 404:109136 2020.
Jacques Janssen, Oronzio Manca, and Raimondo Manca. Applied diffusion processes from engineering to financ.. John Wiley & Sons, 2013.
Jeahan Jung, Heechang Kim, Hyomin Shin, and Minseok Choi. Ceens: Causality-enforeed evolutional networks for solving time-dependent partial differential equations. Computer Methods Applied Mechanics and Engineering, 427: 17036. 2024.
Taniya Kapoor, Hongrui Wang, Alfredo Núñez, and Rolf Dollevoet. Physics-informed neural networks for solving forward and inverse problems in complex beam systems. IEEE Transactions on Neural Networks and Learning Systems, 2023.
Taniya Kapoor, Abhishek Chandra, Daniel M Tartakovsky, Hongrui Wang, Alfredo Nunez, and Rolf Dollevoet. Neural oscillators for generalization of physics-informed machine learning. In Proceedings of the AAAI Conference on Artificial Intellience, volume 38, pp. 13059-13067, 2024a.
Taniya Kapoor, Hongrui Wang, Alfredo Núñez, and Rolf Dollevoet. Transfer learning for improved generalizability in causal physics-ifformed neural networks for beam simulations. Engineering Applications of Artificial Intelligence, 133: 108085, 2024b.
George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422-440, 2021.
Mariella Kast and Jan Hesthaven. Positional embeddings for solving pdes with evolutional deep neural networks. Journal of Computational Physics, 1008:12986... 2024.
Jaeseung Kim and Hwijae Son. -ausality-aware training of physics-informed neural networks for solving inverse problems. Mathematics, 13(7):1057, 2025.
Elham Kiyani, Khemraj Shukla, Jorge F Urbán, Jérôme Darbon, and George Em Karniadakis. Which optimizer works best for wwwsscs-iioormnnneeee informed neural networks and kolmogorov- arnold networks? arXiv preprint arXiv:2501.1637/ 2025.
14

---

Published as a conference paper at ICLR 2026
Nikola B. Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhatacharya, Andrew M. Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces. CoRR, abs/2108. 08481, 2021.
Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. Char- acterizing possible failure modes in physics-informed neural networks. Advances in neural information processing systems, 34:26548 -6560,, 2021.
Isaac E Lagaris, Aristidis Likas, and Dimitrios Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE transactions on neural networks, 965)997-10000 1998.
King-Yeun Yeung Lam and Yuan Lou. Introduction to reaction-diffusiion equations: Theory and applica- fions to spatial ecology and evolutionary biology. Springer Nature, 2022.
Peter Lancaster and Kestutis Salkauskas. Surfaces generated by moving least squares methods. Mathematics of Computation, 37:141-158, 1981.
Henning Lange, Steven L Brunton, and J Nathan Kutz. From fourier to koopman: Spectral methods for long-term time series prediction. Journal of Machine Learning Research, 22(41):1-38, 2021.
Huilai Li, Zhuoqun Wu, Jingxue Yin, and Junning Zhao. Nonlinear Diffission Equations. World Scientific, 2001.
Ye Li, Siqi Chen, Bin Shan, and Sheng-Jun Huang. Causality-enhanced discreted physics-informed neural networks for predicting evolutionary equations. In Proceedings of the Thirty-Thin Interna- tional Joint Conference Artificial Intelligence, pp.
2024.
Zongyi Li, Nikola Borislavov Kovachki, Kamyar Arizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, Anima Anandkumar, et al. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, 2020.
Phillip Lippe, Bas Veeling, Paris Perdikaris, Richard Turner, and Johannes Brandstetter. Pde-
refiner: Achieving accurate long rollouts with neural pde solvers. Advances in Neural Information Processing Systems, 36, 2024.
Qiang Liu, Mengyu Chu, and Nils Thuerey. Config: Towards conflict-free training of physics informed neural networks. arXiv preprint arXiv:2408.11100 2024.
Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature machine intelligence, 3(3):218-22) 2021a.
Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Kamiaaakis. DeepXDE: A deep learning library for solving differential equations. SIAM Review, 63(1):208-22) 228, 2021b. doi: 10.1137/ 19M1274067.
Lu Lu, Raphael Pestourie, Wenjie Yao, Zhicheng Wang, Francesc Verdugo, and Steven Johnson. Physies-ifformed neural networks with hard constraints for inverse design. SIAM Journal Scientific Computing, 43(6):B1105 111110--1130 B1132, 2021c.
L. B. Lucy. A numerical approach to the testing of the fission hypothesis. The Astronomical Journal, 82:1013-1024, December 1977. doi: 10.1086/11220 12164.
Levi D McClenny and Ulisses M Braga-Neto. Self-adaptive physics-informed neural networks. Journal of Computational Physics, 474:111722, 2023.
Robert McGraw, Fan Yang, and Laura M Fierce. Preserving tracer correlations in moment-based atmo- spheric transport models. Journal of Advances Modeling Earth Systems, 1665)=200335000361,, 2024.
Nick McGreivy and Ammar Hakim. Weak baselines and reporting biases lead to overoptimism in machine learning for fluid-related partial differential equations. Nature Machine Intelligence, 100:1256-1269 2024.
15

---

Published as a conference paper at ICLR 2026
Xuhui Meng, Zhen Li, Dongkun Zhang, and George Em Karniadakis. Ppinn: Parareal physics- informed neural network for denendeendent pdes. Computer Methods in Applied Mechanics and Engineering. 370:1 13250, 2020.
Brek Meuris, Saad Qadeer, and Panos Stinis. Machine-leaming- based spectral methods for partial differential equations. Scientific Reports, 13(1):1739, 2023.
Ben Moseley, Andrew Markham, and Tarje Nissen-Meyer. Finite basis physics-informan neural networks (fopinns): scalable domain decompasition approach for solving differential equations. Advances in Computational Mathematics, 49(4))62, 2023.
Johannes Müller and Marius Zeinhofer. Achieving high accuracy with pinns via energy natural gradient descent. In International Conference on Machine Learning, pp. 25471-25485. PMIR, 2023.
Masao Nagasawa. Schrödinger equations and diffusion theory. Springer Science & Business Media, 2012.
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang. Zachary Vito,o, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in PyTorch, 2017.
Benjamin Peherstorfer. Breaking the kolmogorov barrier with nonlinear model reduction. Notices of the American Mathematical Sociery, 69(5):725-733, 2022.
Michael Penwarden, Ameya D Jagtap, Shandian Zhe, George Em Kamiaadaii, and Robert M Kirby. A unified scalable framework for causal sweeping strategies for physics-informed neural networks (pinns) and their temporal decompositions. Journal of Computational Physics, 493:112464 12464, 2023.
Linda Petzold. Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations. SIAM journal on scientific and statistical computing, 4(1):136-148, 1983.
Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez, and Peter Battaglia. Learning mesh-based simulation with graph networks. In International Conference on Learning Representations, 2021.
Les Piegl and Wayne Tiller. The NURBS book (2nd ed.). Springer- Opringer-Verlag, Verlag, Berlin, Heidelberg, 1997. ISBN 3540615458.
M M Powell. The Theory of Radial Basis Function Approximation in 1990. In Advances in Numerical Analysis: Wavelets, Subdivision Algorithmss, and Radial Basis Functions. Oxford University Press, 04 1992. ISBN 9780198534396 doi: 10. 9T0//////978009993333996
Serge Prudhomme, Frédéric Pascal, J. J.Tinsley Oden, and Albert Romkes. A priori error estimate for the naumann-oden version of the discontinuous galerkin method. Comptes Rendus de l'Académie des Sciences Series / Mathematics, 332(0)8551-856 2001. ISSN 0764-4442. doi: https: 0001//00000//000044400000000 J019336--X
Atamert Rahma, Chinmay Datar, and Felix Dietrich. Training hamiltonian neural networks without backpropagation. arXiv preprint 11X1v:2411.17511, 2024.
Atamert Rahma, Chinmay Datar, Ana Cukarska, and Felix Dietrich. Rapid training of hamiltonian graph networks without gradient descent. arXiv preprins arXiv:25006065588 2025.
Maziar Raissi, Paris Perdikaris, and George 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:686-707, 2019.
Pratik Rathore, Weimu Lei, Zachary Frangella, Lu Lu, and Madeleine Udell. Challenges in training pinns: A loss landscape perspective. arXiv preprint arXiv:2402.00 (1000200.800 868, 2024.
Richard B Rood. Numerical advection algorithms and their role in atmospheric transport and chemistry models. Reviews of geophysics, 25(1):71-100, 1987.
Alessandro Rudi and Lorenzo Rosasco. Generalization Properties of Learning with Random Features, April 2021.
16

---

Published as a conference paper at ICLR 2026
David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams. Learning representations by back-propagating errors. nature, 523(6088) 1986.
Guillermo Sapiro. Geometric partial differential equations and image analysis. Geometric Partial Differential Equations and Image Analysis, pp. 440, 2001.
M.S. Shadloo, G. Oger, and D. Le Touzé. Smoothed particle hydrodynamics method for fluid flows, towards industrial applications: Motivations, current state, and challenges. Computers Fluids, 136:11-34, 2016. ISSN 0045-7930. doi: https://doi.org/I0.0 www.oo///.........o.oo.oooooo cotp://f..oompuuu222066.........
Yong Shang and Fei Wang. Randomized Neural Networks with Petrov-Galerkin Methods for Solving Linear Elasticity and Navier-Stokes Equations. Journal of Engineering Mechanics, 150 (4):04024010, April 2024. ISSN 0733-9399) 1943-7889. doi: 10.061/JEMMDTEEEENG 7463.
Ramansh Sharma and Varun Shankar. Accelerated training of physics-informed neural networks (pinns) using meshless discretizations. Advances in Neural Information Processing Systems, 35: 1034-1046, 2022.
Donald Shepard. two-dimensional interpolation function for irregularly-spacee data. In Proceedings of the 1968 23rd ACM National Conference, ACM '68, pp. 517-524, New York, NY, USA, 1968. Association for Computing Machinery. ISBN 9781450374866. doi: 10.1145/800186.. 810616
Justin Sirignano and Konstantinos Spiliopoulos. Dgm: A deep learning algorithm for solving partial differential equations. Journal of computational physics, 375:1339- 1364, 2018.
Shashank Subramanian, Robert M Kirby, Michael W Mahoney, and Amir Gholami. Adaptive self- supervision algorithms for physics-informet informed neural networks. In ECAI 2023, pp. 2234-2241. IOS Press, 2023.
Jingbo Sun, Suchuan Dong, and Fei Wang. Local randomized neural networks with discontinu- ous Galerkin methods for partial differential equations. Journal of Computational and Applied Mathematics, 445:115830, August 2024. ISSN 03770427. doi: 00.1066f.cam.2224.11533
Ziya Uddin, Sai Ganga, Rishi Asthana, and Wubshet Ibrahim. Wavelets based physics informed neural networks to solve non-linear differential equations. Scientific Reports, 13(1):2882, 2023.
Jorge F Urbán, Petros Stefanou, and José A Pons. Unveiling the optimization process of physics informed neural networks: How accurate and competitive can pinns be? Journal of Computational Physics, 13656, 2025.
Roberto Verdecchia, June Sallou, and Lufs Cruz. A systematic review of green ai. Wiley Interdisci- plinary Reviews: Data Mining and Knowledge Discovery, 13(4) 2023.
Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C Carey, ilhan Polat, Yu Feng. Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261-22 272, 2020. doi: 10. 10.1038//41 592-0119-22---6688
Haoxiang Wang, Tao Yu, Tianwei Yang. Hui Qiao, and Qionghai Dai. Neural physical simulation with multi-resolution hash grid encoding. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 38, pp. 5410-5418, 2024a.
Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient flow pathologies in physics-informed neural networks. SIAM Journal on Scientific Computing, 43(5):A3355-A300) 2021.
Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why pinns fail to train: A neural tangent kernel perspective. Journal of Computational Physics, 449: 10768, 2022.
17

---

Published as a conference paper at ICLR 2026
Sifan Wang, Bowen Li. Yuhan Chen, and Paris Perdikaris. Piratenets: Physics-infommen deep learning with residual adaptive networks. Journal of Machine Learning Research, 25(402):1-51, 2024b.
Sifan Wang, Shyam Sankaran, and Paris Perdikaris. Respecting causality for training physics- informed neural networks. Computer Methods in Applied Mechanics and Engineering, 421: 116813, 2024c.
Sifan Wang, Shyam Sankaran, and Paris Perdikaris. Respecting causality for training physics- informed neural networks. Computer Methods Applied Mechanies and Engineering. 421: 116813, 2024d.
Yiran Wang and Suchuan Dong. An extreme learning machine- method for computational pdes in higher dimensions. Computer Methods in Applied Mechanics and Engineering, 418:116578 2024.
Lei Wu and Jihao Long. Spectral-Based Analysis of the Separation between Two-Layer Neural Networks and Linear Methods. Journal of Machine Learning Research, 23(1), January 2022. ISSN 1532-4435.
Mingtao Xia, Lucas Böttcher, and Tom Chou. Spectrally adapted physics-informed neural networks for solving unbounded domain problems. Machine Learning: Science and Technology, 4(2): 025024, 2023.
Chenhui Xu, Dancheng Liu, Yuting Hu, Jiajie Li, Ruiyang Qin, Qingxiao Zheng, and Jinjun Xiong. Sub-sequential physics-informed learning with state space model. In Forty-second International Conference on Machine Learning, 2025.
Jiachen Yao, Chang Su, Zhongkai Hao, Songming Liu, Hang Su, and Jun Zhu. Multiadam: Parameter- wise scale-invariant optimizer for multiscale training of physies-informed neural networks. In International Conference Machine Learning, pp. 39702-39721. PMLR, 2023.
Yuan Yin, Matthieu Kirchmeyer, Jean-Yves Franceschi, Alain Rakotomamonjy, and Patrick Gallinari. Continuous pde dynamics forecasting with implicit neural representations. International Conference Learning Representations, 2023. URL htllhttps:// /wtp///peeeevveeeeeee www. forum?
hnttps://openrvvieee
id=B73niNjbPs
Jeremy Yu, Lu Lu, Xuhui Meng, and George Em Karniadakis. Gradient-enhanced
informed neural networks for forward and inverse pde problems. Computer Methods in Applied Mechanies and Engineering, 393:1 14823, 2022.
Yaohua Zang, Gang Bao, Xiaojing Ye. and Haomin Zhou. Weak adversarial networks for high- dimensional partial differential equations. Journal of Computational Physics, 411:109409, 2020.
Zhiyuan Zhao, Xueying Ding, and Aditya Prakash. Pinnsformer: transformer-basee framework for physics-informed neural networks. In The Twelfth International Conference Learning Representations, 2024.
18

---

Published as a conference paper at ICLR 2026
APPENDIX
CONTENTS
Introduction
2
Solving time-dependent PDEs using Frozen-PINNe
2.1
Frozen- N ansatz
2.2
Computing hidden layer parameters without gradient descent
2.3
Solving time-dependent PDEs using PioNee--PNNss by separation of variables
2.4
Approaches for satisfying boundary conditions for FIozes-PINNs
2.5
SVD Layer
2.6
Summary of the training a algorithm for Frozen- PINNS
Empirical results
3.1
High advection speeds, fast convergence, and long-time simulation
3.2
Higher-order derivatives in space and time
3.3
Multi-scale solutions
3.4
Non-linearity and shocks
3,5
Non-lincarity and complicated domain geometry
3.6
Chaos and strong non-linearity
3.7
High-Dimensional PDHs with low-dimensional solution manifolds
3.8
High-dimennionnalty
4
Conclusion
Extended review related work
21
Supplementary methodological details on PDE solvers
22
B.
Physics-Informan Neural Networks
22
B.2
Pivenn-PNNN-wim and Piieen--INNNelm..
................................
23
B.2.1
Extended Discussion on Frozen-PINNS
23
B.2.2
Computing spatial and temporal differential operators in PDES
27
B.2.3
Reformulating PDEs as ODEs using Frozen-PINN ansatz
28
B.2.4
Handling boundary conditions via boundary-compliant layer
31
B.2.5
Handling boundary conditions via augmented ODE
B.3
IGA-FEM
31
Supplementary details on numerical experiments
32
Linear advection equation
33
Buler-Bernould equation)
35
Wave equation
38
19

---

Published as a conference paper at ICLR 2026
C.4
Burgers
42
C.4.1
Comparison with classical spectral methods
44
C.5
Nonlinear diffusson equation
46
Nonlinear reaction-diffusion equation
50
Kuramoto-Sivashinsly equation
54
High-dimensional diftusion equation
54
20

---

Published as a conference paper at ICLR 2026
A EXTENDED REVIEW OF RELATED WORK
In this section, we provide a comprehensive extended review of the literature and highlight how it relates to our work.
Physics-informed neural networks are widely used to solve PDEs with neural networks. In this work, we benchmark our approach against various PINN variants such as adaptive activation
PINNs (Jagtap et al. 2020), self-adaptive PINNs (McClenny & Braga-Neto 2023). wavelet PINNs (Uddin et al. 2023). and causal PINNs (Wang et al. 2024c). among others. For high-frequency temporal variations in the PDE solutions, Krishnapriyan et al. (2021) propose curriculum learning with gradually increasing advection coefficients. Compared to curriculum learning, our approach with space-time separation is much easier implement, computationally efficient, and accurate, as we demonstrate in Section [.1] Subramanian et al. (2023) propose using adaptive self-supervision of PINNs for sampling collocation points using the gradient of the loss function. We instead use the solution gradient to capture locally sharp features in the solution (see Section 3.4). Many specialized
approaches based on PINNs (Cho al. 2024, Meng et al., 2020, Sharma & Shankar 2022] Chiu
et al. 2022), methods based on hash-encoding (Huang & Alkhalifah| 2024 Wang et al. (024a). and
transfer learning (Kapoor et al.) 2024b) have been proposed, but are still based on gradient-bases iterative optimization and back prop-prooag,tion, unlike ours.
Other recent advances of PINNs include methods that model the PDE system as pseudo-sequences. For instance, PINNsFormer employs Transformer-based architecture that constructs pseudo- sequences from spatio-tempora temporal samples and uses self-attention to model long-range temporal dependencies (Zhao et al.] (2024). Another work, PINNMamba, is based State Space Models (SSMs) and sub-sequence alignment, enabling continuous discrete temporal modeling and improved propagation of initial-condition information (Xu et al.] 2025). Although these methods model PDE systems as pseudo-sequences, these architectures often lead to more computational time and out-of-memory issues owing their architecture, as presented by Xu al. (2025).
Physics-informed approaches using randomized neural networks for solving PDEs have mostly been studied by combining Extreme Learning Machines (ELMs) with the self-supervised setting of PINNs (Chen et al. [0244a] Wang & Dong 2024 |Shang & Wang [2024] (Sun et al. (024). For instance, Dwivedi & Srinivasan (2020) propose physics-iffomnet extreme learning machine (PIELM) to efficiently solve linear PDES, while Calabrò et al. (2021); Galaris et al. (2022) employ ELMs to learn invariant manifolds as well as PDEs from data. Dong & Yang (2022) show that given a fixed computational budget, ELMS achieve substantially higher accuracy compared to classical second-order FEM and slightly higher accuracy compared to higher-order FEM. For static, nonlinear PDEs, ELMs can be used together with nonlinear optimization schemes (Fabiani et al.) 2021). On larger spatiotemporal domains, Dong & Li (202T) and Dwivedi et al. (2021) propose using multiple distributed ELMs multiple subdomains. Although the aforementioned methods simplify the optimization problem by randomly sampling hidden layer parameters and fixing them, they treat time as merely another spatial dimension. As result, their neural basis functions span the full spatiotemporal domain, which limits their accuracy on PDEs exhibiting figh-frrquency temporal dynamics, unlike our approach.
2025) discuss While the problem setting is restricted to Hamiltonian systems, Rahma et al. 2024
how to train Hamiltonian neural networks and Hamiltonian graph neural networks using ELM and SWIM approaches, and demonstrate how random sampling can be leveraged to significantly speed up training compared gradient- based iterative optimization. In this work, we show how random sampling can speed up training and resolve optimization challenges of PINNs for time-dependent PDEs.
Neural Galerkin schemes (Finzi et al. 2023, Aghili et al., 2024 Berman et al. 2024 Bruna et al. 2024) offer an alternative to the full spatiotemporal approach of the randomized neural networks and PINNs. These approaches treat all or sparse subsets of network parameters, beyond just the last layer's 's parameters, as dine-dependent. This leads to a much larger system of ODEs compared to our approach. The work neural implicit representations (Chen et al. 2023 Yin et al.) (023)) also uses neural basis functions to represent only the space component, but relies on gradient- based iterative optimization via back propaaapp,,,,,n,, unlike our approach.
21

---

Published as a conference paper at ICLR 2026
Spectral methods for solving PDEs promise fast convergence with much fewer basis functions. Mcuris et al. (2023) present a method to extract hierarchical spatial basis functions from a trained DeepONet and employ it in spectral method to solve the given PDE. Xia et al. (2023) integrate adaptive techniques into PINN-based PDE solvers to obtain numerical solutions of unbounded domain problems that standard PINNs cannot efficiently approximate. Laane et al. (202T) propose spectral methods that fit linear and nonlinear oscillators to data and facilitate long-term forecasting of temporal signals. Dresdner et al. (2022) demonstrate spectral solvers that provide sub-grid corrections to classical spectral methods to improve their accuracy. Du et al. (2023) use fixed orthogonal bases to learn PDE solutions as map between spectral coefficients and introduce training strategy based on spectral loss. These methods differ from ours in problem setting, architecture, and training.
Neural operator frameworks (Lu et al. 2021a Kovachki et al., 2021 Li et al. 2020 Platf et al. 2021) are promising but are typically trained with PDE solutions with different initial conditions, spatial domains (geometries), or parameter settings. |atar et al. (2025) have demonstrated how continuous-time neural networks can be constructed for linear operator approximation for linear and time-invariant systems. Instead, in our setting here, we solve the PDE using given coefficients, domain, and initial conditions without relying on any training data. The ease of implementation, rapid training, and high accuracy of our backpropagation-fie approach can be leveraged to generate PDE solution data for training operator networks.
Mesh-free methods are typically based on radial basis functions (RBFs, (Powell 1992 Chen et al. 2014)) Moving Least Squares (MLS) (Shepard 1968: Lancaster & Salkauskas 1981). These often do not have user-friendly software or are only applicable in specialized settings (c.g., smoothed particle hydrodynamics, (Lucy 1977 Gingold & Monaghan 1977 Shadloo et al. 2016)). Moreover, despite the ease of dealing with complicated geometries, these methods typically suffer from many challenges, such as the choice of kernel, imposing boundary conditions, and convergence issues. These methods are not the focus of this work.
Classical numerical methods such as finite elements, finite volumes, and finite differences have been used to solve PDEs for decades. They often have rich theoretical grounding and high accuracy. Isogeometric analysis (IGA) is one such method, in which spline- -ased basis functions are defined over structured grid (Hughes et al. 2005; Cottrell et al., 2209 2006). Mesh-based methods often entail a time-consuming setup phase, especially when mesh generation is challenging. Methods like sparse grids enable adaptivity through hierarchical bases but pose significant implementation challenges, particularly for irregular domains (Bungartz & Griebel 2004). In this work, we benchmark our results against IGA and finite-element-bassee methods.
B SUPPLEMENTARY METHODOLOGICAL DETAILS ON PDE SOLVERS
B. PHYSICS-INFORMEN NEURAL NETWORKS
This work benchmarks Frozen PINNs against many prominent variants of physics-informed neural networks. While we directly report results from other works for many PINN variants for different PDE benchmarks (see Table D.. we also implement two PINN variants for certain PDEs classical physics- informed neural network (PINN) (Raissi et al. 2019) and causality-respecting physics-informed neural network (causal PINN) (Wang et al. 2024c). We now describe these two variants.
Classical PINNS are feedforward deep neural networks designed to approximate PDE solutions by incorporating physical laws into the learning process. The architecture of a vanilla PINN includes deep neural network that maps inputs (c.g., space and time coordinates) to outputs (c.g.. physical quantities of interest) and is trained to minimize loss function that combines data and physics-basee errors. The data term ensures that the neural network fits the provided data points, while the physics term enforces the PDE constraints with automatic differentiation. The constraints on initial and boundary conditions are satisfied via additional loss terms. The loss function for a classical PINN is:
L(p)= LPD(()))) + )eLic(4) + ABL + A1
(6)
where μ represents the trainable network parameters, and X;, for 1, 2, 3, 4 represent the weighting factors for individual loss terms, which are hyperparameters. In this work, we consider the setting of unsupervised learning and thus neglect the data loss term.
22

---

Published as a conference paper at ICLR 2026
Let N be the total number of training points, which is the sum of the number of interior training points Nin (where the PDE residual evaluated), initial condition training points Nic (where the initial condition is evaluated), and boundary condition training points N, (where the boundary condition is evaluated), We denote the neural network solution at point (x("), {(m) in the computational domain by w" (x(1), ((m)). We consider the generic nonlinear PDE defined by equation defined by
The PDE loss term
Nint Σ n=1
LPDE(A)) = Nint
+ Lue (())) + (x(a")(
s(a())))1P. (7)
The boundary condition loss term is defined as
No.
LBc (1+) = Nb
Σ
Bu Bui ((0)) g(2(0))|1p
(8)
Similarly, the initial condition loss term defined as
1 Nic
Lic(p) = Nic Σ (1(g(z((0)) ==(
(9)
We now describe a Causal PINN, which modifies the PINN loss function impose temporal causality, inherent in time-dependent PDEs, as soft constraint. In conventional PINNs, the loss is computed without prioritizing accuracy at earlier times, which disrupts temporal causality. The Causal PINN remedies this by assigning weights at each time step based on the cumulative loss from previous steps, ensuring that the model concentrates on accurately approximating solutions at earlier times before moving forward. This tries to incorporate the causal structure of the physical problem being solved as a soft constraint. The causal PDE loss term defined by
N
LPDE(H) =
w. W.LpDE(ti, u), where
(10)
W1 1, W2 = Ei EPDE(tu.))
for = 2.
Here N₂ represents the number of time steps into which the computational domain is divided. The causality hyperparameter regulates the steepness of the weights and is incorporated in the loss function, similar to Kapoor et al. (2024b). This modification introduces a weighting factor w; for the loss at each time level t,, with being dependent on the cumulative PDE loss up time f,. The network prioritizes fully resolved solution at earlier time levels by exponentiating the negative of this accumulated loss. Consequently, the modified PDE loss term for causal PINN is expressed as
LPDE(H) N₂
w LoD((Et()) + Eisl Ledz(te.)) LPDE(L...)) LPppEtt))
(11)
B.2 FROZEN-PINTY PINN-SWIM AND FROZEN-PINE PINNN--INN-ELLL ELM
B.2.1 EXTENDED DISCUSSION ON FROZEN-PINNS
Difference compared to training physics-inforee neural networks:
We summarize the difference between training classical physics-informed neural networks and Frozen PINNs in Figure g
Comparison between Frozen-PINN and Frozen-PINN-elmm One of the main factors in- and Frozen-PINN-ellm the underlying solution
fluencing the performance of Frozen-PINNY of the PDE. We explain, with example of the Burgers' equation, how the SWIM sampling can be leveraged when the solution has steep gradients, as one can sample localized basis functions in the part of the domain where the solution has steep gradients. For ELM. the probability of sampling steep basis functions with the vanilla ELM is lower, as illustrated in the Figure 2 Even one uses different distribution sample the network parameters such that more basis functions with steep
23

---

Published as a conference paper at ICLR 2026
Compute derivatives
Basis functions in the space-time domain
Lpdz(0)
Physics/PDE loss
Total loss:
Lpc
1 0) x L.LLpp LFDE(
Boundary condition loss
+ Ay Loc
+ As LIc(
Initial condition loss
Network parameters:
Solve the coupled optimization problem: man L(0)
Algorithm: Gradient-base Iterative optimization (back-propageel
Classical training physics-informet informed neural networks using gradient- based iterative optimization.
Basis functions only the space domain (sampled and frozen)
Time-dependent output layer parameters
Compute derivatives
LPDE (8)
Classical
solvers
step-size
Physics/PDE loss
LBc(#)
control
ú(x, t
OR
Boundary condition loss
Impose hard constraints
Initial condition lass
L
Solve Least
squares
problem
Network parameters:
Total loss: L(0) A LPDe(() Ag LBc(0) + Lic(0)
de-coupled optimization prohlem: min L(
(Minimize the initial condition loss separately)
Algorithm: Sampling, least squares, and classical ODE solvers
Frozeen PINNs: Training physicx-informeed neural networks without gradient-descent.
Figure 8: Comparison of Prozen-PINNY (bottom row) that leverage gradient-descent-Fie training algorithm, with conventional PINNs (top row) that rely on gradient-based iterative optimization: conventional PINNs basis functions in the entire spatio-temporal domain and solve fully coupled optimization problem involving multiple loss terms via gradient-baat based iterative training. In contrast, Frozen-PINNs sample basis functions only in space, make time dependence explicit only in the output layer, decouple initial/boundary conditions, and leverage least squares and adaptive ODE solvers. Parameters dependent on space, time, and both are indicated by blue, orange, and blue-orange colors, respectively, offering direct visual representation of the space-time separation in Frozen-PINNI Notation: The network output a(æ, t, 8) approximates the solution to the PDE. The total loss term (L(O)) sums three loss terms one for the initial condition (Lic(!)), one for the boundary conditions (fB)())), and one for the PDE residual (LPde(f) that together impose physical constraints.
gradients are sampled, placing the basis functions at appropriate spatial locations is another challenge. With ELM, one cannot resample or choose basis functions using data as it is data-agnostic. Thus, especially if the solution has localized steep gradients. the performance of ELM typically worse compared to SWIM. We additionally demonstrate with snapshot of the Burgers' solution that SWIM basis functions exhibit rapid exponential decay of error with increasing network width, where ELM,, Fourier, and Chebyshev basis functions used in classical spectral methods suffer from the Gibbs phenomenon (see Appendix C.4.I) and lead to poor scaling and accuracy (see Figure 23 Figure 22).
If the underlying solution is sufficiently smooth and does not have steep gradients anywhere in the domain, ELM typically performs very well, as seen in the example with the Advection equation (see Section 3.1. Euler Bernoulli equation (see Section 3.2), and sishdimnnssionll diffusion equation
(Section 3.8 where Frozen- ProzzenPPIN--ee PINN-elm performs much better than Frozen-PINN. PINN-swim as shown in
Table Whit While just the vanilla SWIM algorithm in the presented results, easily adapt
24

---

Published as a conference paper at ICLR 2026
the algorithm and, after sampling the network parameters with SWIM, multiply the basis functions with a tunable scaling factor before applying the non-linearity to sample many more basis functions with shallow slopes.
Thus, the choice between the two strategies is particularly governed by the underlying solution of the PDE. Apart from the favorable cases for each method mentioned above, both methods have comparable performance and typically outperform PINNs by several orders of magnitude in speed and time. Thus, the rapid training of our approach could be leveraged to try out both approaches if one has no information about what the solution of the PDE could look like.
Influence of random sampling on the method Similar to the question of how PINNs trained with Adam/SGD perform based on their random network initialization, understanding the influence of weights on the output is a challenge. There are two main differences between (stochastic) gradient- based optimization and our setting. First, after fixing the internal weights, we regularized least-squares (not a stochastic method) to fit the initial condition. Second, we do not use a stochastic method to solve over time. Therefore, even though PINNs can adapt their random initialization over the gradient-based optimization, precisely that optimization also adds stochasticity. If the number of neurons for the model increases, the randomness in case decreases because the regularized least-squares fit to the initial condition (which converges to a single solution in the limit of many neurons), while stochastic gradient descent will only converge a distribution (because of mini-batch optimization). This has been observed for the supervised learning problems in Bolager et al. (2023) particularly in the transfer learning experiments. In Table [ we observe that our model's 's performance is often orders of magnitude better, and the variance is on the same scale as the magnitude.
"data-driven" and "data-agnostic" sampling In this work, we assume that we do not have access to the true solution of the PDE. The term "data-driven sampling" can be misleading for the problem setting of this paper, which concerns unsupervised learning tasks. Thus, here we clarify what we mean by data-driven sampling. Our data are random pairs of collocation points, but we do not have access to the true function values (because, at the initial time point / = O, we have not solved the PDE yet). Thus, even though we do not have access to the true solution of the PDE, we call this "data-driven" sampling because we create the parameters of our basis functions (neurons) so that they are centered strictly within the domain. We achieve this by using data points sampled in the domain, thereby considering the geometry and bounds of the spatial domain. Note that with data-agnostic sampling in ELM, the neurons can easily be centered outside the spatial domain because weights and biases are chosen without considering any information about the geometry and bounds of the spatial domain. To summarize, though our algorithm proposes "data-driven" sampling, we do not start with time-series data and instead work in self-supervised setting.
Rationale for constructing outer basis functions One might reasonably ask that if one knows the outer basis functions analytically, why add another layer just to approximate them with tach basis functions? When analytical basis functions are known, they should be used directly. However, in many cases, such expressions are not readily available. We argue that this idea of a boundary- compliant layer can be quite powerful for PDEs where the basis functions are not known analytically but only through boundary conditions, which we can then incorporate by constructing useful outer basis functions. For instance, to solve the diffusion equation complex geometries, one use the optimal bases consisting of the eigenfunctions of the Laplacian operator computed numerically at discrete points as the outer basis functions (Coilman & Laton 2006). Thus, representing them with tanh basis functions facilitates a straightforward computation of the derivatives needed for solving the PDE.
Kolmogorov n-width barrier Without resampling the internal network parameters, our method faces the Kolmogorov n-width barrier PPeherstorfer (2022); Du & Zaki (2021) ); Berman & Peherstorfer (2024); Kast & Hesthaven (2024) because our basis functions are not time-dependent. However, resampling basis functions at certain time points of the Prozeen-INNN-winn (as done in the Burgers' equation in Section 3.4) results in solution- and time-dependent basis approximation of the solution manifold and, thus, in theory, can break the barrier. PINNS can theoretically break the Kolmogorov n- n-wid barrier as time is treated an extra spatial dimension, and internal network parameters are dependent. However, for PINNs, the optimization issues pose much more severe challenges even on very simple PDEs and in low dimensions (Krishnapriyan et al. 2021 Wang al., 2021 2022).
25

---

Published as a conference paper at ICLR 2026
Rel i, error (space-tine solution)
107
Frozen-PINN PINN-l
Frozen-PINN-die
Frozen-PINN-swim
Frozen-PINN-swim
10"
104
10-2
10-5
10-8
10-M
10-10 10-s 10-2
10-12 10-9
10-6
10-1
10
Fzond
Rel La error (initial condition)
Figure 9: (Advection Equation): Empirical analysis of how initial condition loss affects the perfor- mance of Frozen-PINNY PINNs: (Left): Relative error of the full space- time solution Vs reond, (Right): Relative error of the full space-time solution grows with the relative error of the initial condition (controlled via rcond).
So even though our vanilla Poozen-PINNNssii///oooPPPPNN approach (without periodically resampling hidden layer weights) faces the Kolmogorov n-width barrier, outperform PINNs, typically by several orders of accuracy and time in practice.
Analysis of residual initial condition loss and its impact on model performance For all PDEs considered here, the initial condition is relatively easy to fit, and one can approximate it accurately by sampling enough collocation points at using enough basis functions, and setting a relatively low regularization constant (+e 10 10-12). So. the initial condition loss is not the bottleneck in any of the PDEs we considered here. However, interesting investigate the impact of initial condition loss on the model's performance.
Given N₂ collocation points X E R.N.xd .M neurons, and hidden layer output ((X), the initial condition is computed via least squares solution:
C(0) u) (x))
(12)
where [o(X,, 1] (M++)) and the pseudo-inverse is denoted by .+. We compute this least squares solution C(0) using the Python function np linalg..... lstsq, which takes as an argument reond which is the cut-off ratio for small singular values of [x),), 1). High cut-off ratios reduce the accuracy of the least-squares solution, while very low ratios lead to poorly conditioned systems that introduce numerical errors.
We perform an experiment by progressively increasing the cut-off value to deliberately degrade the mntiaa-condiiton fit and study its impact the overall PDE residual. We solve the advection equation given in Equation (24). using the same hyperparameter settings as in Table 5 fix the advection coefficient 40, and vary rcond from 10 10 10-17 Figure 9 shows that small rcond values yield highly accurate initial cnntal-condition fits and low relative error. As r rcond increases, more dominant singular values are discarded in the least-squares solve, degrading the initial-condttonn representation and leading to larger errors in the full space time solution. Thus, for Frozen-PINNs, maintaining reasonably accurate initial condition fit important, inaccuracies can influence the ODE solve and increase the overall error.
On the choice of the SVD truncation threshold: The SVD truncation threshold crucial hyperparameter for Frozen-PINNs determining the dimensionality of the ODE solver and affecting the speed and accuracy of Frozen-PINNs. We conduct an ablation study on the SVD truncation threshold for Burgers' equation (Equation (8 the nonlinear diffusion equation (Equation (2) ). and the 10-dimensional diffusion equation (Equation (33)) with hyperparameters described in Table [4]
Table 18 Table 29
respectively. We change the width to 200 for the non-linear diffusion equation, and for all PDEs, we only vary the SVD truncation thresholds. For detailed problem setups please refer Appendix C.4 Appendix C.5 Appendix C.8
Figure highhightt the keoln trade-off between accuracy and speed. Retaining fewer singular values reduces the dimensionality and stiffness of the last-layer ODE, yielding faster solutions with
26

---

Published as a conference paper at ICLR 2026
100
10
10a
Releevv E3 errer
10-1
Relative &2 errer Time
Resetvvls f, errer
10"
Relative
102
20-1
10
Tima 64
10-1
Tima (4)
Reeetive (s
10"3
Relative ig artor Time
Time 6O
10
20-1
10
10
-100
10
10-0
10-1
10
10-15
10-39
10-1 10-4
10
10-9
10-4
16
SVD cutoff
SVD cutaff
SVD cutoff
(a) Burgers' equation
(b) Nonlinear diffusion
(c) 10-dimensional diffusion
Figure 10: Impact of the SVD truncation threshold E used for the SVD layer on the time-to- solution and accuracy of Frozen-PINN across three PDEs. Black triangles in the nonlinear diffusion plot indicate solution blow-up at large SVD cutoffs.
similar or slightly lower accuracy. Retaining more singular values increases dimensionality and stiffness, which slows the solver but improves accuracy. Importantly, the performance is robust for SVD truncation thresholds eşvn 10-10
The choice of the optimal SVD truncation threshold depends highly on the application constraints. Higher thresholds (esv D 10-8) are suitable for faster solutions with moderate accuracy, while very low thresholds (ESVD S. 10-13) are preferable when high accuracy is the primary goal. The default value we choose 10
represents a good trade-off between accuracy and speed. We always set the rcond (regularization constant for the initial least squares solve) to the same value as the SVD truncation threshold because: (a) it also represents the cut-off ratio for the SVD of the feature matrix for the initial least squares solve and it does not make sense to solve this with extremely high or low precision when the data has already been passed through the SVD layer, to maintain similar level of truncation the SVD layer. (b) Empirically, we set the regularization constant to be equal to the SVD threshold (see Figure ( where the relative eπror stays the same for the rcond values in the range 10 10-11 when the SVD truncation threshold is set at 10 10-14).
observe very robust performance if
On the choice between the two strategies for enforcing boundary conditions: In practice, the choice between the boundary-complian layer and the augmented ODE follows a simple cost-benefit tradeoff.
We recommend using coundayy-complann layer when a rroleem-specfiic transformation @ is easy to derive (e.g., zero Dirichlet, periodic boundary conditions on simple domains). It enforces boundary conditions (almost) exactly and does not enlarge the ODE system, so typically more efficient. The main limitation is that it requires deriving Φ (X), which may be non-trivial for complex geometries or boundary conditions.
We recommend using the Augmented ODE strategy when boundary geometry or constraints make mapping difficult. This is universally applicable and requires no
an analytic problem-specific engineering since it soft-enforces boundary conditions by augmenting the state, at the cost increasing system dimension and possibly worsening stiffness.
B.2.2 COMPUTING SPATIAL AND TEMPORAL DIFFERENTIAL OPERATORS IN PDES
We use the notation described in Section of the manuscript. We first discuss how to compute different spatial and temporal derivative terms appearing in the PDEs described in this manuscript using the neural network ansatz. We then use these expressions to reformulate the PDEs described in this manuscript as corresponding ODEs. We consider neural networks in the most general setting by considering the outer basis functions and the SVD layer (cf Algorithm[1]
Computing spatial derivatives: We list and describe how to compute the spatial derivatives the approximate PDE solutions:
27

---

Published as a conference paper at ICLR 2026
First-order spatial derivative of the approximate PDE solution is computed as:
ü. (x t) (.))
(13)
=C o e
where @ is the Hadamard product, and
(2(2:) (0=(=) (=(=))===+> :b...., (:():...
RRM.xd
(14)
with (2(2) E RRM, and Ø= is the first derivative of the tanh activation function.
Second-order spatial derivative of the approximate PDE solution is computed as:
Uza (3r, 1) C(
(15)
= =] & W ( E Riddd
where Tz (x) is defined equivalently as o (z) but with dzx being the second-order spatial derivative of the tanh activation function.
The Laplacian of the approximate PDE solution is computed as:
Aa(r,t) (=)) where, 1 E Rdx1 [A, @ W e oz (2), 0]4 R
(16)
Fourth-order spatial derivative of the approximate PDE solution is computed as: Gurrz (ar, = [( (e)za
(17)
= [[ & W WOW @ @zzza( 0] E IR
where Ozzza is the fourth-order spatial derivative of the lanh activation function.
Computing time derivatives: We now list and describe how to compute the time derivatives of the approximate PDE solutions:
First-order time derivative the approximate PDE solution computed as:
(1 t) = A,,
(18)
Second-order time derivative of the approximate PDE solution computed as:
Dez (x. t) = C.. (t) (t).()....
(19)
B.2.3 REFORMULATING PDES AS ODES USING FROZEN- PINN ANSATZ
The partial differential equations considered in this work recast ordinary differential equations for evolving output layer coefficients, making use of the spatial and temporal derivatives derived in Appendix B.2.2 We denote the pseudo-inverse by
Advection equation:
The one-dimensional advection equation is
Ue (z. t) + Bug t) = O,
where 3 is scalar. Approximating the solution with neural network ansatz (Equation 2 and substituting Equation (18) and Equation (13) in the advection equation, we get,
C.) = - C.(t) = (9
The initial condition given by
C(0) u ((x00) (X)) [x]]
28

---

Published as a conference paper at ICLR 2026
Euler-Bernoulli equation: The Euler-Bernouli PDE considered in this manuscript is
Ust + Uyans [ 1).
Approximating the solution with neural network ansatz (Fquation and substituting Equation 17 and Equation (19) in the Beler-Bernoulli equation, we get,
(x)))((X) ((X,)) C(I) manz (X)
We rewrite this second-order ODE as a combination of first-order ODEs given by
(t) D(t), (x())((X) f(x,t) C(t) Pawre (X).
We then reformulate the ODEs as
(C.(L) D.(1) = (C(1) D(1) 0
- - (e(X)annXX(X)) + +(0 x)=x
The initial condition is given by
C(0) - m(X,0) x(X)*
D u, (X,0) (xx)* (xx))
The extension to the Fuler-Bernoulli beam equation on a Winkler foundation straightforward, where the reformulated ODE is written as:
(C,(t) D.(() = (C(t) D(t)) (i - -(((X)xxa) JeeZz (x)))))(((x)) + (0 L)(s(X,()) ((x))
0
Wave equation:
The wave equation considered in this manuscript is
Utt KUl s(X, t)-
Approximating the solution with neural network ansatz (Equation and substituting Equation and Equation (19) in the wave equation, we get,
15
(+(())((X) = ((X,)) - ,x))
We rewrite this second-order ODE as a combination of first-order ODEs given by
(t()) D(t),
(x())((X) = -
We then reformulate the ODEs as
(C.(1) D.(1) (C(1) D(O) (0 -xD(X) xx0(()) +0 [[FXXX))TT((X)"""""
The initial condition given by
C(0) x(x,0)=x(x) D(0) n. (X,0) (x)+
Burgers' equation:
The one-dimensional Burgers' PDE we consider is
Wr + uw, Arll... 0,
where a scalar.
Approximating the solution with neural network ansatz (Equation 2 and 18 Equation (33 and Equation
substituting Equation
(15 in the Burgers equation, we get,
(t A. (X) - ( & (()) +
(t(t) = (C 4 (X) @ (1) (x)))(X) + a (C) { Ar (X))) (X))*
Note that the non-linearity is transferred to the right-hand side of the ODE, The initial condition is given by
C(0) u) (X)+
29

---

Published as a conference paper at ICLR 2026
Nonlinear diffusion equation: The two-dimensional nonlinear diffusion equation consider is
udu - f(r, t). x E S2 C R E [0, 1].
(20)
Approximating the solution with neural network ansatz (Equation and Equation (15) in the nonlinear diffusion equation, we get,
) substituting Equation (8).
(()))((X) = (C(t) @ [C(t) Pe (X)]1) Cx(t) (C (t) (X) @ (t 2 (X) +)]" (
[f(X,t)]* (f(x,))
Note that the non-linearity is transferred to the right-hand side of the ODE. The initial condition is given by
((0) = u(x.0) (((xX)))T(xx'''.
Nonlinear reaction-diffusion equation: The two-dimensional nonlinear diffusion equation consider is
Us Au - f(,,t, t), I TEΩCH R. E [0.
(21)
Approximating the solution with neural network ansatz (Equation substituting Equation (18). and Equation (15) in the nonlinear diffusion equation, get,
(t) [=t (C()) (CC)))))))))))) (())))))()))))))))))))))())))))) (F(X)))))))()))) []],
(t = ([C(t) (C()) & C(t)
+ t)) o(x)+
The non- non-lincarriy u² is transferred to the right- side of the ODE. The initial condition is given by
C(O) u( (xx)
High-dimensional diffusion equation: The d-dimensional diffusion equation we consider
Au f(z, t), E C Rdd t [0.
(22)
Approximating the solution with neural network ansatz (Equation 2 substituting Equation (18) and Equation (15) in the diffusion equation, get,
(x)) = Cx(t) = (IC (t) Pxx (I(XX)))+++[X))X)))
(f(X,t))
[f(X, +) 4(x)*
The initial condition is given by
C(0) = ul u(X,0) (X)+
Kuramoto-Sivashinkt. equation: The governing PDE considered this manuscript is
"4 0ut, Bugz Murzz f(X,t)
Approximating the solution with neural network ansatz (Equation (2) and substituting Equation Equation 15) Equation (17) and Equation (18) in the Kuramoto-Sivashinsk equation, we get,
13
C,(t) = f(X,t) aC(t)((X) @ C(t) (X) CCt (X) - 7C) (X).
We can reformulate as:
Ct(t) = (f(X,t) - aC(t) @ C(t)
- C(t) - (e)) 4(X)+
The initial condition given by
-(0)= u(x 0) +
We use the solve_ive routine of the SciPy package Virtanen et al. 2020)) One can pass test points in time
Note on ODE solvers and interpolation time:
argument to the method solve ivp. One can optionally set the parameter dense dewee_eeueeeeeeeeee output to true, which means that the output of the ODE is function handle that can be evaluated by interpolation at any time point Ω. The method specified dictates the interpolation order. RK23 uses a cubic Hermite polynomial, while DOPRI85 uses seventh-order polynomial.
30

---

First introduced in Hughes et al. al. (2005). Isogeometric analysis (IGA) is a numerical method developed unify the fields of computer-aided design (CAD) and finite element analysis (FEA). The key idea is to represent the solution space for the numerical analysis using the same functions that define the Rational Published as a conference paper at ICLR 2026
B.2.4 HANDLING BOUNDARY CONDITIONS VIA BOUNDARY-COMPLIEE LAYER
To enforce periodic boundary conditions, is sufficient for each basis function satisty the periodic condition individually, as the Frozen-PINN ansatz, which is a linear combination of these functions, will inherently satisfy as well. For instance, for a one- oeedimensional spatial domain, we find A so that (x1))) Ax,x,, where I,, are the left and right boundary points of the domain. In this paper, for certain PDEs (see Appendix [.. for x E Ω and = 1, ........ Mss we approximate (A)) = sin(kx) (for k even) and [A
cos(kar) (for & odd) and set co = for all t. For zero Dirichlet boundary condition given by u(z) = O, we can use the technique described above by choosing basis functions so that Aő(x) 0 for öΩ. For other boundary conditions, propose using the augmented ODE trick to satisfy the boundary conditions.
B.2.5 HANDLING BOUNDARY CONDITIONS VIA AUGMENTED ODE
Our approaches to satisfying the Dirichlet and periodic boundary conditions are already explained in the main text. Here, we explain how we handle deeeddepndennt Dirichlet boundary conditions and Neumann boundary conditions.
Time-dependent Dirichlet boundary conditions: For handling time-dependent Dirichlet bound- ary conditions (u(2, () = g(z, () for E 892), we set A to the identity map and augment the ODE (Equation (3) with an additional equation given by
de (a., 1) ge (1. t) for x E an = Ct = [ C(t)), (Xb, t)] 4 (x((X,x,)'
EREX (
ER(Ne M6X(CM, 11)
we know the solution on the boundary at all time points, which continuously differentiable. If the solution on the boundary points is not available at all time points, one can interpolate and approximate the derivative of the solution on the boundary.
In the example Section 3.5
Neumann boundary conditions: For simple spatial domains, one can choose appropriate outer basis functions as described in Section 2.4 that inherently satisfy the Neumann boundary conditions. For instance, for zero Neumann boundary conditions on a one-dimensional domain, one can choose outer basis functions consisting of cosines of different frequencies scaled to the domain (function value is the boundaries) so that their spatial derivatives, which are the sine functions, are zero on the boundary points.
On complicated domain geometries, to satisfy Neumann boundary conditions (Vu(x, t)- f(x) for x E 8?), we set A to the identity map and augment the ODE (Equation with an additional equation for the boundary points and solve
(1( [ C()), 1 (X),,, (X.)))X)
ERAX(N++M)
CR(No (R(N++M))x(M,,+))
B.3 IGA-FEM
geometry in CAD (Cottrell et 2209), which include the B B-Splines and
B-Splines (NURBS) (Piegl & Tiller 1997).
In this paper, we use B-Splines as the basis functions. The B-Splines are defined using the Cox-de Boor recursion formula (COX 1972 de Boor 1972), i.e.,
[ & < otherwise,
Neo(E)
0
5- N NN.p-1(E) + Eitp+1 - N++1--1((E) Eitp &
NE,(E)
Eip41 E
31

---

Published as a conference paper at ICLR 2026
where & is the ith knot, and is the polynomial degree. The vector E E2.......... + is the knot vector, where n is the number of B-Splines. By specifying the knot vector, we define the basis functions
to solve the PDEs. We use uniform open knot vector, where the first and last knots have multiplicity p+ 1, the inner knots have no multiplicity, and all knots that have different values are uniformly distributed. We refer the knots with different values as "nodes". The intervals between two successive nodes are knot spans, which can be viewed as "elements". The elements form "patch". A domain can be partitioned into subdomains, and each represented by patch. In our work, we use a single patch to represent the entire ID domain. Figure 11a shows an example of such a patch, where the B-Splines are CC"cooninuous within the knot spans and C continuous at the inner knots. In order to address the boundary conditions, we adapt the B-Splines shown in so that the boundary conditions are directly built into the solution space.
Figure
11b
Figure
11c
1.0
1.0
1.0
0.4
0.2
0.2
0.0
0.4
0.0
0.4
0.0 0.400..000
10
Number of basis
func-
(b)
Number of basis
func-
(c) tions
Number of basis func-
tions
Lions
Figure 11: Examples of B-Splines representing the ID domain [0, 1]. Number of nodes 6 and degree of polynomials 2. (Left): The original B-Splines. (Middle): Adapted B-Splines to satisfy the Dirichlet boundary condition. (Right): Adapted B B-Splines satisfy the periodic boundary condition. Note that the first (blue) spline is identical to the second last (brown) one, and the second (orange) spline is identical the last (pink) one, they share the coefficient. The gray dashed lines indicate where the domain starts and ends.
In the following, we refer to the adapted B-Splines as basis functions (xr) Thus, the solutions of PDEs are approximated by
K
u(x, () Σ
(a(D)(()).
We solve the PDEs in the weak formulation. For the linear advection equation (see Equation (4) the weak form of the equation is
Σ (t) / (n(z)u(z)dz
+ + E cn(t) / ()(())(()) = 0,
(23)
where v(x) are the test functions. The test functions are chosen to be the same as the basis functions. The integral of the functions is computed using Gaussian quadrature. Then we solve the linear Ordinary differential equation (ODE)
Mć + Kc 0, O,
where matrix M and matrix contain the integral the B-Splines and their derivatives, and the coefficient 3, which are given. We solve the Suler-Bernoull equation equation 25 and the Burgers' equation equation 28 in a similar way. The boundary condition for the Euler-Bernon equation is, in addition, weakly imposed, is done in Prudhomme et al. (2001).
C SUPPLEMENTARY DETAILS ON NUMERICAL EXPERIMENTS
Here, we discuss additional experimental details for the PDEs considered in this work. We start by listing the details the code repository, FEM software, hardware, metrics, and ablation studies:
32

---

Published as a conference paper at ICLR 2026
Code repository: The source code, along with the instructions reproducing the results, is provided in the supplemental material (zipped file) and is made publicly available (see Section B) The code repository provides Python scripts and notebooks that can be executed and tested readily. Moreover, we make the refactored repository publicly available and intend to actively maintain it.
FEM code: In this paper, we use DOLF N 0.8.0 to solve the nonlinear diffusion equation (see equation 29), DOLFINx (Baratta et al., 2023), which is part the FEniCS project, is C++ and Python library used for solving PDEs with the finite element method (FEM). It provides tools for defining complex geometries, formulating variational problems, and solving them efficiently distributed architectures. We used the software Gmsh (Geuzaine & Remacle 2009) to generate a mesh for this experiment with complicated geometry, as shown in Figure 24a
Hardware details: The computational experiments for Frozen PINNs, FEM, and IGA-FEM were performed with: Ubuntu 20.04.6 LTS, NVIDIA driver 515.105.00 and i7 CPU.
Metries for computing errors: We use the Root Mean Squared Error (RMSE) and the relative L error to quantify errors in all experiments (see Appendix for the definitions). We compute the test error on uniform grid for all PDEs with 256 points in space and 100 points in time, unless otherwise specified. We use float 64 numerical precision in all the experiments.
Let d be the dimension of space and Ω x [0, T] C Rd x be the spatio-temport domain. Given N points in test set X, the error metrics we to compare numerical results Root Mean Squared Error (RMSE) and relative 2 error given by
EBEx (Mtrue (I) Upred (a))²
RMSE
N
and
Exe (Utrus (z) Upred xx))
Relative 2 error
Ezee
(Merue (2:)²
For each experiment, the mean and standard deviation of the RMSE and the relative L error are computed with three seeds.
blation studies for neural architecture and SVD layer: We perform ablation studies whenever necessary for the neural architectures considered in this work. Importantly, also perform an ablation study on the SVD layer. To quantify the compression in width after the SVD layer, we define a compression ratio as C, Ma, where M, is the width of the (sampled) hidden layer before the SVD layer (assuming no-boundary-comple r the width of the SVD layer (see Figure 3 We define a speed-up in computation time as Tooad as the ratio of computational time without the SVD layer to the time required Ted with the SVD layer.
layer), and
We now describe the detailed problem setups, ablation studies, and plots comparing the results of Frozen-PINNS with those other approaches for all PDEs considered here (see Figure [2].
C.I LINEAR ADVECTION EQUATION
The linear advection equation describes the transport of a quantity and is used model many real-life applications. such as simplified traffic flow models, transport of pollutants in rivers or the atmosphere (Rood) 1987; McGraw al. 2024). Here, we consider the linear advection equation with periodic boundary conditions given by
Problem setup:
ue(x,)) t) + Bug (2. t) O,
for x E [0, e [0, 1].
(24a)
u(æ, 0) = sin(z),
for x E [0, 2m],
(24b)
((0,t) u(2n, t), for t [0, 1).
(24c)
The analytical solution of Equation 24 is given by u(x, t) = sin(æ - Bt). We describe detailed hyperparameter settings used for the experiments on: (a) high-advection speeds (how the error grows
33

---

Published as a conference paper at ICLR 2026
Application
Key challenges posed by the PDE
PDE and ground truth solution the PDE
BC Periodic
equation
Extremely high- frequency temporal variations
Rapid ascillations
u(x, t)
1000.
(2): Long time span
Case 2. 8 1.7 1000
Euler-Bernoulli
Higher-order deriva- tives time
PDE (ciossical): Maser = f(x,
beam equation
(winkler foundation):
Mrr Moxse KM f(z t)
BC: (u, Rkgs)
Ceesical
Winkler foundation
Wave
Multi-scale solutions
PDE:
Vag
Mas
0
equation
BC
-.3337
Diffusion equation
Non-linearity complicated domain geomtery
NÂN f(x, t)
Complicated
BC Time- dome-ddeenndntt Dirichett
geometry
t=
1=03
t= na
1.
Burgers' equation
Non lincarity and shock
PDE:
Mt + Ui2,, MA
BC: Dirichlet
Shock
5-dimensional diffusion equation
Non-linearity high dimensionality
PDE: 14 Au u f(x, t)
Timeedeeendent Dirichlet
twn
Twn4 twbs
100 dimensional
High-dimensionallit
PDE: 14: Au (1c,
diffusion equation
BC
Time--dependent Dirichlet
Kuramoto- Sivashinsky equation
Non-linearity and Chaotic system
PDE: Ut callyy + Blur 7 Mynss
Figure 12: Overview of the PDE benchmarks considered in this study, highlighting the core challenges associated with each problem and their corresponding ground truth solutions. Boundary conditions abbreviated as BC.
34

---

Published as a conference paper at ICLR 2026
with the advection coefficient 3), (b) convergence (how the enror decays with the number of basis
functions for a fixed advection coefficient 3 = 10), (c) error for advection coefficient _ = 40 (for a comparison with other PINN-based variants), and (d) long-time simulation for T 1000 seconds for a fixed advection coefficient & = 1 in Table [3] Table Table5] and Table [ respectively. The hidden layer weights for ELM and Prozen-PIIN-es
are sampled from the Gaussian distribution and biases from a uniform distribution in |-4, 4]. For SWIM and ELM, we use 1000 interior points for E {10-2 2, 1001-01-000 10), and we use 8000 interior points for 8 (40, 1003. The code repository contains all the necessary Python notebooks to reproduce results for Frozen-PINNs for all the different cases of the advection equation considered here, including the three key experiments concerning high advection speeds, convergence, and long-time simulation (see Section 3.1
Deeper networks and further optimization experiments: Additional experiments are carried out for baseline PINN on deeper networks with 10 and 20 hidden layers, where each hidden layer has 30 neurons. The experiments are run for 20000 epochs using Adam and L-BFGS optimizers under multiple learning rates for the advection equation with 8 = 10 and / 40. Tables[[lnn] summarize the RMSE, relative L errors, and training times for 3 = 10 and B = 40 cases, respectively, The results show consistency with the known literature of PINNS. First, deeper PINNS do not directly lead to better performance. Increasing depth from 10 to 20 layers often degrades accuracy for both optimizers, reflecting the optimization difficulty of fully-connected PINNs as the models become deeper. This behavior is also discussed previously in the literature, for instance by Wang et al. (2024b), where the authors show that PINN performance is known degrade when larger and deeper neural network architectures are employed. Second, the 3 = 40 case is known to be challenging for PINNs (Krishnapriyan al. 202T) due to high frequency features, and the results presented in Table 7 show failures across depths. learning rates, and optimizers. The results show that even with larger networks, longer training. and different optimizers, standard PINNs face challenges in achieving high accuracy. especially for 3 40, and require more computational time. Thus, matching the accuracy of standard PINNs with the proposed method inherently challenging and computationally expensive.
Ablation studies: For the advection coefficient 3 = 10, the ablation study for Frozen-PIND PINN-swim, Prozen-PINN-elm, and vanilla PINNs is already presented in Figure4[Middie for varying the number of neurons and interior points. The ablation studies for PINNs for the network width and number of interior points are presented in Table 9 and Table I0, respectively. Since the network width is already quite low for optimal parameters, the SVD layer does not further reduce the dimension of the ODE system. Hence, we do not perform ablation studies for the SVD layer in PINNN, as is not used in this case.
Comparison of results: Figure 13 shows the absolute errors obtained with the Frozen-PINN --NN-swim, Frozen-PINN
PINN, Causal PINN, and IGA methods along with the ground truth for 3 40. One can observe that all approaches considered here, besides Frozen-PINNs and IGA-FEM, fail to capture the high-frequency temporal dynamics. Figure 14 shows the true solution at B 1 for the example with long time simulation.
C.2 EULER-BEREE BERNOULLI EQUATION
Problem Setup: The time-dependent Euler-Bernoulli beam equation models the dynamic behavior of beams, including vibrations and transient loads. It is used to model loads on rail tracks, bridges, and aircraft wings, among many other applications (Beskos) 1987). We consider two different types of the Fuler-Bernoulli beam equations in this work. The first is the classical Fuler-Bernouli beam equation that models a simply supported beam with varying transverse force and is described as
Utt Uxxxx = f(x, 4) x E [0, / E [0, 1],
(25a)
where f(r,t) = (1 - 1677) sin (z) cos with initial and boundary conditions
u(x,0) = sin(x), ue 0) = 0, u(0, t) = u(T,t) t) = uxx (0x t) = uxx (
(25b)
(T, t) 0.
(25c)
The forcing function and the analytical solution are taken from Kapoor et al. (2023).
We consider another variant of the Bereollll beam equation, with a Winkler foundation (an elastic, deformable foundation) given by:
Met Mxxxx p(=, f(x, 1) [0, [0, 1].
(26a)
35

---

Published as a conference paper at ICLR 2026
5.0
5.0
5.0
2.5
2.5
2.5
0.0
0.0
0.8
0.0
0.2
0.4
0.6
0.8
1.0
3.3e-16
8.1e-08
1 1.1e-16
4.5e-08
8.0e-11
1.6e-10
(a) Frozen- Frozen-PINN -ee----NN---iii-----
(b) Frozen- Seeenn-INNN-elm
(c) IGA
5.0
5.0
2.5
2.5
2.5
8.0
0.2
0.4
0.6
0.8
1.0
0.0
0.0
0.2
0.4
0.5
0.8
1.0
0.0.
0.2
0.4
0.6
6.0e-08
5.4e-01
1.1e+00
6.3e-07
1.5e+00
3.0e+00 00
-8.0e-01
0.0e+
8.0e-01
(d) PINN
(c) Causal PINN
(f) Ground truth
Figure 13: Advection equation (B 40): absolute error plots and ground truth.
True solution
5
0
0
200
400
600
800
1000
t
-0.50 -0.25
0.00
0.25
0.50
0.75
Figure 14: Advection equation 3 = 1, (long time simulation): Analytical solution u(xr, t)
sin(x 3t).
The forcing term is f(z,t) = - 7²? sin (x) cos((Tt), with the initial and boundary conditions
u(r, 0) = sin( ux( 0) = 0,
(26b)
u(0,t) u(ST = (0, t) = Mxx (8T, t) = 0.
(26c)
The foundation reaction force p(z, t) assumed be proportional to the displacement of the beam and modeled as p(z, t) = czu t), where x is the spring constant and is set to 1 in this case. The forcing function and the analytical solution for the Euler-Bernoulli beam equation with Winkler foundation are taken from Kapoor et al. (2024b).
The ablation studies for the PINN-based variants for classical Euler Bernoulli and the one with the Winkler foundation could be found in Kapoor et al. (2023 2024b). The hyperparameters for various neural PDE solvers used for solving the classical Euler-Bernoult PDE and the one with the Winkler foundation are described in Table I1 and Table 12 respectively. The
Ablation studies:
hidden layer weights for Frozen-PINe -iinn-PINN---- are sampled from the Gaussian distribution and biases from uniform distribution in [-2,2] 2).
For this example, we have employed the augmented ODE strategy to satisfy boundary conditions and obtain the results presented in Table We empirically investigate the effect of the penalty term in the augmented ODE on the performance of Frozen-PINNs, considering both accuracy and computation time. For this experiment, we use the hyperparameters described in Table [1 for Frozen-PINN-elll (high-precision regime) and vary the value x. Figure 15 shows that: (a) for k: > 10S. the boundary loss is negligible (RMSE 10 10-10) and the total loss is very low (RMSE 10 6-10-9) -0), and (b) for extremely large x 2 108, the augmented ODE becomes slightly stiffer, resulting in increased solution time.
36

---

5 on the losses and Published as a conference paper at ICLR 2026
Table Advection equation (high-advection speeds): Network hyper-parameters used for 3 E {10-?, 10 1, 10, 40, 100, 1000, 10000) to study the influence of the advection coefficient on the errors (optimal hyper-patee parammarseeterss in bold) (see Figure Left).
Parameter
Value
Frozen-PINN PINN-swim,
Number of hidden layers
Frozen-INNN-clm
Hidden layer width Outer busis lunctions Activation
1140, 380, 560)
110, 14, 20, tanh
-eguullarization
[10 10 10, 10
12 10 14
Loss
squared erroe
houndury condition stritegy
houndary-complin layer
SWIM, ELM
Numher of hidden layers
SVD curoff
10 10
Hidken layer width Activation 1. Prregularivation
[140, 560] tanh
[10 s 10 10
Loss
mean-squared error
Initial and boundary points
400
IGA
Number of modes
16
Degree of polynomials
Number of basis functions
15
PINN
Number of hidden layers
width
[10, 20,30, 40]
Activation
tanh
Optimizer
LBFGS
Epochs
5000
mean-squared
Learning rate
0.1
Batch size
200
Purametter initialization
Xavier Glarott Bengio 2010
Loss weights, A Ag
1.
Interior points Initial and houndary points
[500, 1000, 1500, 2000)
600
Causal PINN
Number of hiddee layers
Layer width
30
Activation
Optimizer
ADAM followed by LBFGS
Fpochs
LBPGS Epochs
5000
Loss
mean-squared error
Learning
0.1
Ratch size
2000
Parameter initialization
Xavier Glorot Bongio 2010
Loss weights, As Aa
1.1
Interior points
40000
Initial and boundary points
6000
Causality parameter,
10
Total loss
Boundary loss
Time
102
10a
RMSE
10-4
10
Time (s)
10-8
100
10a
10a
10a
107
10a
Augmented ODE penalty term x
Figure 15: The effect of penalty term in the augmented ODE (see Equation time to solution.
Comparison of results: Figure and Figure 17 present the absolute errors for the classical Euler- Bernuuli PDE and its variant with Winkler foundation, respectively, using Frozen-PINN-swinn Frozen-PINN-elm PINN, and IGA methods, along with the true solution. The error plots for the Fuler-Bernoulli beam equation with Winkler foundation for other variants of PINNs, such as
37

---

Published as a conference paper at ICLR 2026
Table 4: Advection equation (convergence for 10): Optimal hyper- saper-parameters in the experiment designed to study how the error decays with the number of basis functions in the neural network (see Figure Middle).
Parameter
Value
Frozen- PINNS (both variants)
Number of hidden layers
Hidden layer width
[2...., 30 tanh
Activation
1.2 regulatization
[10 10-5, 10 1,,1,,0000111111 1,1010,, 100-1,100111 10-10, 10011
Loss
mean-squared error
boundary condition stralegy
boumdary compliant layer
Number of hidden layers
Layer width
[10, 20, 20,30,440000 40] tana
Activation
Optimizer
LBFGS
Epochs
5000
Loss
meon-squared
Learning rate
0.1
Batch
200
Purameter initialization
Xavier Glarot & Bengio 2010 1.
Loss weights, A1. A2
Table 5: Advection equation (for 3 = 40): Hyper-parametter for the results in Table
Parameter
Value
Frozen-PINN-eln (low-precision)
Number hidden layers
Hidden layer width Outer basis functions cutoff Activation
50
[14]
[10-12]
tanh
1. regulaarization
10
ODF solver tolerance
Loss
squared
houndary condition strategy
compllaan layer
Frozen: PiNN-swim (high-pret precision)
Number of hiddee layers
Hidden layer Outer basis functions svd cutoff Activation
[14]
[10-12]
tanh
L--regularization
10
ODE solver tolcrance
10
Loos
mean-squared error
boundary condition stratogy
boundary-compliant layer
hidden layers
Layer width
[10, 20, 40
Activation
tanh
Optimizer
LBFGS
Epochs
5000
mean-squared
Learning rate
0.1
Batch size
Parameter initialization
Xavier Glorot Bengio 2010
Loss weights, A1, Au
Wavelet PINN, causal PINN, adaptive PINN, and self-adaptiv PINNs, can be found in Kapoor et al. (20245). Table ssoows the summary results for the classical the variant considering the Winkler foundation for different methods.
beam equation and
C.3 WAVE EQUATION
The acoustic wave equation models the propagation of sound waves through a medium. It describes how pressure velocity evolve over time. We consider the wave equation Ω = 00, 1] for time / E [0,1] from Hao et al. (0024), given by:
Problem Setup:
Uet 4ugr 0, E t 10,
(27a)
38

---

Published as a conference paper at ICLR 2026
Table 6: Results of PINNs for advection equation (3 10) for dense networks (number of hidden layers = 10 and 20, with each hidden layer having 30 neurons). The experiment studies the perfor- mance of PINNs with L-BFGS and ADAM optimizers under different learning rates. Each case is run for 20000 epochs.
RMSE
Relative L error
Training time (s)
Optimizer Learning rate
Hidden layers
Hidden layers
Hidden layers
10
20
10
20
10
20
0.1
6.02e-4
1.23e-2
8.51e-4
1.74e-2
204.25
336.32
L-BEGS
0.01
1.57e-3
1.30e-2
2.22e-3
1.84e-2
213.26
327.41
0.001
2.26e-2
1.64e-2
3.19e-2
2.32e-2
174.84
301.12
Adam
0.0001
6.72e-3
2.07e-2
9.5e-3
2.92e-2
190.24
317.97
Table Results of PINNs for advection equation (3 = 40) for dense networks (number of hidden
layers = 10 and 20, with each hidden layer having 30 neurons). The experiment studies the perfor-
mance of PINNs with ADAM and L-BFGS optimizers under different learning rates. Each case is run for 20000 epochs.
RMSE
Relative L error
Training time (s)
Optimizer Learning rate
Hidden layers
Hidden layers
Hidden layers
10
20
10
20
10
20
0.1
3.09e-3
1.53c+1
4.37e-3
2.17e+1
207.25
320.43
L-BFGS
0.01
1.46e-2
7.07c-1
2.07e-2
1.00c+0
205.08
323.43
0.001
3.66e-1
7.07e-1
5.18e-1
1.00e+0
168.14 14
331.03
Adam
0.0001
6.67e-1
6.82e-1
9.44e-1
9.65e-1
181.89
316.90
Table 8: Advection equation (long-time simulation for = 1): Optimal syper-parameters for Frozen-PINNs in the experiment used to demonstrate that the errors with Frozen-PINNs stay low for simulations up to 1000 seconds (see Figure (ighht)).
Parameter
Value
Frozen- PINN
Number hidden layers
(both variants)
Hidden layer width basis functions Activation
250
Lanh
L d'-reullarization
[10-19]
mean-squared boundary compliant layer
houndary condition strategy
Table 9: Advection equation (3 10): Ablation study for PINN (LBFGS) with respect the network width. The mean is computed over 3 seeds.
Layer width
Training time (s)
RMSE
Relative L error
10
24.47 + 0.19
1.24e-3 + 2.38e-4
1.76e-3 + 3.37e-4
20
27.46 + 0.08
6.52e-4 = 2.59e-4
9.22e-4 + 3.66e-4
30
30.43 + 0.50
3.69e-4 + 4.33e-5
5.23e-4 + 6.13e-5
40
33.64 + 0.41
3.86e-4 9.37e-5
5.46e-4 + 1.32e-4
with the initial condition
u(x, = sin(TT)) + sin(4 x Ω
(27b)
u1 0) = 0, E St
(27c)
and the boundary condition
m(0,t) - u(1, te [0.1].
(27d)
39

---

Published as a conference paper at ICLR 2026
Table 10: Advection equation (3 10): hyperparameter optimization for PINN (LBFGS) varying the number of interior points. The mean is computed over seeds.
Interior points
Training time (8)
RMSE
Relative I,2 error
500
25.76 + 0.29
4.10e-4 + 7.20e-5
5.80e-4 + 1.01e-4
1000
27.44 + 0.25
3.72e-4 + 4.06e-5
5.27e-4 +: 5.74e-5
1500
29.61 0.16
5.68e-4 1.97e-4
8.03e-4 + 2.79e-4
2000
30.43 + 0.50
3.69e-4 + 4.33e-5
5.23e-4 = 6.13e-5
2.5
x
08.00
0.25
0,50
0.75
1.00
8.00
0.25
0,50
0.75
1.00
0,0
0.2
0.8
2e-08
4e-08
2e-08
5e-08
7e-08
0.0e+00
2.0e-07
4.0e-07
(a) PNniss-m
(h) Frozen-PINN-ee
(c) IGA
2.5
2.5
2.5
8.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.0
0.2
0.5
0.8
1.0
6.0c-08
5.4e-03
1.1c-02
1.8c-06
6.1e-02
1.2c-01
-8.0e-01
0.0e+00
8.0e-01
(d) PINN (LBFGG)
(e) PINN (Adam)
(f) Ground truth
Figure 16: The classical Euler-Bernoude Bernoulli beam equation: absolute error plots and ground truth.
0
10
20
10
20
10
20
2e-06
4e-06
5e-07
1e-06
De+00
2e-08
4e-08
(a) Frozen-PPin -INGnPPIN------ swim
(b) Frozen-PINN-oo
(c) IGA
10
0.5
0.5
Ground truth
Figure 17: The Kuler-Bernoud beam equation on Winkler foundation: absolute error plots and ground truth.
The analytical solution of the problem
u(x, () (1x)))))) c2T() + sin(ATx) cos (8T().
(27e)
40

---

Published as a conference paper at ICLR 2026
Table 11: The classical muler-Bemouulli beam equation for the results in Table hyperparameters.
summary of all
Value Hidden layer width
Parameter
Prozeu-PINN-clm (low-precision)
Number hidden layers
50
SVD-cutolI
10-6
Activation
tamh
1,2 re-rraaaarizatonn
10
Loss
mcan-squared ervor
boundary condition strategy
augmented
Frozen-PINe PINN elm (high precisio))ssim)
Number hidden layers
2
Hidden layer width
100
SVD-cutoff
10
Activation
tamb
I,a regularizzatin
10 10
Loss
mean mean-squared error
boundary condition strategy
sugmented ODE
IGA
Number of modes
27
Degree polynomials
Number of basis functions
33
PINN
Number hidden layors
Layer width
20
Activation
tash
Optimizer
LBFGS (ADAM)
Epochs
15000 (30000)
Loss
mesn-squared error
Learning rate
0.1
Batch
2000
Parameter initialization
Xavier Glorot & Bengio 2010
Loss weights, A1.
Interior points
10000
Initial and boundary points
6000
Table 12: The Euler-Bernoulli beam equation Winkler foundation for the results in Table summary of all hyperparameters.
Parameter
Value
Frozen-PINN-elm (low-precision)
Nomber hidden layers
2
Hidden layer width
200
SVD-cutolf
10 8
Activation
tanh
1.-regularizztton
10
mean-squared error
boundary condition strategy
augmeeted ODE
Frozen PINN swim (high precision)
hickden layers
Hidden layer width
400
SVD-cutoff
10 10
Activation
tanh
T.3 .:-regalarization
10
Loss
mean-squared enorr
boundary condition strategy
augmented ODE
IGA
Number modes
Degres polynomials
Number of basis functions
63
PINN
Number of hidden layers
4
Layer width
200
Activation
tanh
Optimizer
LBFGS
Epochs
10000
Loss
mean- squuree-urrde eroor
Leaniing rate
0.1
Batch size
500
Parameter initialization I Loss weights, A1.
Xavigr Glorot & Bengio 2010
# Interior points
10000
boundary points
1500
41

---

Published as a conference paper at ICLR 2026
In addition, present significantly challenging scenario involving multi-scall solution by employing the initial condition
1 u(z, 0) sin(Tx) sin (ATx) sin(972), S
(271)
Me = 0, E
(27g)
for which the corresponding analytical solution is given by
1 sin sin(4 cos(8 + sin() cos (18mt.
u(x, t) (in( cos (2 + 2
(27h)
4
Comparison of results: Table [ we compare Frozen-PINNy PIinn--PNN---e--A swim result of Equation 27e with other PINN methods from Hao et al. (2024), see Figure [8] The hyperparameters used for this experiment can be found in Table [3] Additionally, we demonstrate the capability of our proposed method by solving the multi-scale problem in Equation (27h) using the same neural network architecture, and the relative r,a error is less than 10 The results are shown in Figure 19
Table 13: Wave equation:
for the result in Table
Parameter
Value
Frozen-PINN.oo PINN-PNNN--IN Number of hidhle layers
Hidden layer width Outer basis functions Activation L -regularization
400
tanh
10
Loss boundary condition strategy
mean-squared error
boundary- layer
True solution
Frozen-PINN-stn PINN-swim solution
Absolute difference
1.0
1.0
1.0
0.8
0.8
0.8
02
0.2
0.2
0.0
0.0
00.5
1.0
-1.0
-0.5
1.0
Figure 18: Wave equation Equation
27cr: Ground truth, Prozen-PINN-swim solution, absolute error.
C.4 BURGERS
Problem Setup: The Burgers' equation different settings used model traffic flows, large- scale structure formation in cosmology, and shock formation in inviscid flows, among other appli- cations (Bonkile et al. 2018). The inviscid Burgers' equation is nonlinear PDE, which can form shock waves. We consider Burgers' equation on Ω = [-1, 1] for time t (0, 1] from Raissi et al. (2019), given by: Wr + uu, - (0.01/7)um =0, Ω i [0, 1],
(28a)
with initial and boundary conditions
w(0, x) sin((Tx)) E Ω
(28b)
utt, -) u(t, 1) - t [0.
(28c)
We consider the analytical solution provided by Basdevant al. (1966..
42

---

Published as a conference paper at ICLR 2026
True solution
Frozen- Paann--iN---iiin solution
Absolute difference
1.0
1.0
1.0
0.8
0.4
0.4
0.4
0.2
0.2
0.0
0.5
0.0
1.0
0.0
0.5
2.0
-1.0
-0.5
0.0
0.5
1.0
-1.0
-0.5
0.0
0.5
1.0
Je-5
(Left)) Ground truth, (Middle): Prozen-PINN-swiig solution, (Right); point-wise absolutc error.
1.345
-1.337
(b) visualization of the multi-scale wave solution.
Figure 19: Wave equation Equation (2hh: Ground truth, Frozen- solution, absolute error, and multi-scale solution visualization.
Why Frozen-PINN-elm can't resolve shocks in PDE solutions? To accurately resolve PDE solutions with sharp gradients, is essential to: (a) construct basis functions with steep gradients, and (b) position them particularly near the shock regions within the domain. Figure 20a (Right) illustrates why solution- or data-agnostic ELM basis functions make it very difficult for Frozen- PINN-elm to capture the sharp features in the solution, particularly at the center of the domain, due to the exponentially small probability of sampling steep basis functions (Huang et al. (0066). While sampling weights from wider uniform distribution, as discussed by Calabró et al. (2021) for linear PDEs, can increase the probability of sampling steeper basis functions, offers no spatial control over their placement.
Ablation studies: We describe additional details in solving the Burgers' equation with various neural PDE solvers in Table [4 and Table 15) The results of the ablation study with the number of neurons in the hidden layer for PINN-wim
presented in Table T6] We observe that starting with width of 1200, the error decreases for a width up to 600 and increases again below 600. We believe that for widths lower than 600, the network capacity seems to be the reason for the loss of accuracy. For very high widths, the regularization constant has to be kept to higher value to avoid overfitting- Otherwise, the ODE system becomes highly stiff. With this high regularization constant, the training becomes stable, but it affects the training accuracy. We do not include results for Frozen-PINN-ellm as it fails on all widths, as it is not able to capture the sharp shocks and exhibits Gibbs phenomenon (Gottlieb & Shu 1997). which is explained in detail in Appendix C.4.1
We also perform an ablation study for the SVD layer for Frozen-PINN PINN-swim. Please refer to Table 17 The ablation study reveals that the SVD layer compresses the number of neurons by a factor of 1.58, which reduces the output computation time by a factor of for almost the same accuracy. This highlights the utility of the SVD layer.
43

---

Published as a conference paper at ICLR 2026
SWIM basis functions
ELM basis functions
1.0
1.0
0.0
-0.5
-1.0
-1.0
-0.5
0.5
1.0
-1.0
-0.5
1,0
Burgers Gradient =882.69
Burgers
Gradient - 40.7
Gradient 150.62
Gradient - 0.94
Gradient 2.09
Gradient 141.73
Gradient 0.66
(a) (Left): Re-sampled SWIM basis functions (with steep gradients centered around the shock) at 0.66, (Right): data-agnostic ELM basis functions.
Ground truth
Prozen-PINN-swim (resampl ing)
Ground truth
Prozen-PIMN-swimm
Absolute error
t-
0.25
t-
0.50
t- 0.75
0.5
0,5
0,5
x
0.0
0.0
0.0
uft ()
-0.5
-0.5
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
-0.8
0.0
67
16:66
(b) (Left): Ground truth, (Middle): Grozen-PINN-swing solution, where black and gray dashed lines mark time snapshots selected for comparison (in (d) the right) and the collocation points resampling times, respectively. (Right): point-wist wise absolute error.
(c) Comparison between Frozen- and
numerical solutions at three time instances.
Figure 20: Illustration of experimental results for the Burgers' equation.
Comparison of results: Figures 20b and 20c present a comparison between the Frozen-PINN-swinn solution and the numerical solution from Raissi al. (201)), validating the ability of Frozen-PINN- swim to resolve shocks with high accuracy. We demonstrate with a snapshot of the Burgers' solution that SWIM basis functions exhibit a rapid exponential decay of eπror with increasing network width, where Fourier and Chebyshev basis functions suffer from the Gibbs phenomenon Gottlieb & Shu (1997) (See Figure 22 Figure 23 Appendix [C..TT) Figure (ouure27shore shows the absolute errors obtained with the PINN, Causal PINN, and IGA methods.
C.4.1 COMPARISON WITH CLASSICAL SPECTRAL METHODS
In this section, we study how the basis functions sampled with SWIM and ELM approaches perform in comparison the basis functions typically employed in traditional spectral methods. We try to approximate single snapshot of the solution to the Burgers' equation, which has a locally steep gradient. If method fails even approximate single snapshot well enough, is highly unlikely to achieve better results in approximating the entire space-time solution of the PDE.
Figure shows the approximation of the Burgers' equation solution at f = 0.99, using SWIM basis functions, ELM basis functions, Fourier series, and Chebyshev polynomials, respectively. The number basis functions is 102 for all methods. Figure 23 shows the approximation error using a different number of basis functions. We can see that for ELM basis functions. Fourier basis functions, and Chebyshev polynomials, there are oscillations near the shock, and the error is large compared to the SWIM basis functions, where sampling appropriate basis functions in order to adapt to the target function well. Note that this experiment, the weights for the ELM basis functions sampled from Gaussian distribution
able to take advantage of resampling data points and
44

---

Published as a conference paper at ICLR 2026
1.0
1.0
0.0
x 0.0
0.0
8.0
-1.8 0.0
0.2
0.4
0.5
0.8
-1 0.0
0.2
0.4
0.6
0.8
1.0
0.0e++00
8.0e. 03
1.6e:
2.0e-02
4.0e-02
0.0e+000
1.0e-03
2.0e-03
(a) PINN
(b) Causal PINN
(c) IGA
1.0
0.0
0.0
0.2
0.4
0.6
0.8
1.0
.8. 3. 01
0.0e+00
8.0e-01
(d) Ground truth
Figure 21: Burgers' equation: absolute error plots and ground truth.
Table 14: Burgers' equation: Summary of hyper-parameters for Frozen-PINy PPize (see Table
Parameter
Value
Frozen-PINN- --ii- (low-precision)
Number of hidden layers
2
Hidden layer width Activation
[300]
tanh
//-regularization
[10
7,
10
10
10
12)
svd cutoff
10-6
Loss
mean-squared
# collocation points (space)
6600)
sampling points
10000
ODE solver tolerance
10-3
time windows for resampling boundary condition strategy
augmented ODE
Frozen-PINN-sww hhgg-procisioo)
Number hidden layors
2
Hidden layer width
[450]
Activation r.2 Pregulaarizatinn
tonh
[10
10 1000-
10 10-10012 12,01,2,000 18)
svd cutolf
Loss
mean squared emror
collocation points (space)
1000]
# sampling points
6000
ODE solver tolerance
10-6
windows resampling
houndary condition strategy
augmented ODE
Prowen-PINN-tie
Number hidden layers
2
Hidden layer width
[2000]
Activation L -egullarrzzio
tanh
[10-6 10-7, 10 10-10, 10 10
Loss
mean-squared error
collocation points (space)
3000]
# sampling points
[6000]
boundary condition strategy
augmented
with standard deviation of 10 in order to increase the number of basis functions. The biases sampled from a uniform distribution [-10, 10]. For the Fourier basis functions and Chebyshev polynomials,
equispaced grid points. We also experimented with quadrature points and placed more points near the steep gradient an attempt to mitigate the oscillations associated with the Gibbs phenomenon and the Runge phenomenon, but it did not lead to any significant improvement in the
45

---

Published as a conference paper at ICLR 2026
Table 15: Burgers' equation (see Table D Network hyper-parameterr used for PINN, Causal PINN, and IGA.
Parameter
Value
Number hiddea layers
Layer width
20
Activation
tanh
Optimizer
LBFGS
Epochs
10000
Loss
mcan-squared
Leurnnng rate
0.1
Batch size
Parameter initialization
Xavier (Glorot Bengio 2010
Loss weights, A1. Aa
1.
Interior points
10000
Initial and boundary points
600
Causal PINN
Number of hidden layers
Layer width
20
Activation
tanh
Optimizer
ADAM followed by LBFGS
ADAM Epochs
5000
LBFGS Epochs
10000
Loss
maannssuarod error
Leurning
0.1
Batch size
200
Parameter initialization
Xavier Gloort Rengio 2010
Loss weights, A1. Az
1. 1
Interior points
40000
Initial and boundary points
600
Causality parameter,
IGA
Number of nodes
750
Degree of polynomials
Number of busis luncciions
756
Table 16: Burgers' equation: ablation study for the network width for Frozen-PINN-swin.
Width
Relative 1 error
240
4.27e-4
550
2.27e-7
800
2.78e-6
1200
1.54e-6
Table 17: Burgers' Equation: Ablation Study for the SVD layer with PoNzen-PINN-swim..
With SVD layer
Without SVD layer
Ratio
Number of neurons
500
316
Width Compression 1.58x
Time (s) Rel. L2 error
141.5
989.84
Speed-up x 7x
3.34e-4
3.28e-4
results. This conclusively demonstrates that SWIM basis functions perform better than traditional bases used in spectral methods in accurately resolving shocks.
C.5 NONLINEAR DIFFUSION EQUATION
Problem Setup:
The non-linear diffusion equation in different forms is used model the spread of populations, bacterial colonies, and forest fires, as well as to model groundwater and ice-sheet flow in glaciers, and mass diffusion in reactive flows (Li al. 2001). We consider two-dimensional nonlinear diffusion equation given by
Ut - uAu = (.,,),t), (a.y) e a te [0,1),
(29a)
with forcing function
f(x, y, 5c 55 "ssinTTT) -3 (-1 é sin( (-12 22)))
(29b)
46

---

Published as a conference paper at ICLR 2026
groud truth
approximation
SWIM basis
ELM basis
Function approximation error
0.5
(ttxx)
u(txx)
10
0.0
-0.5
Relative to error
10-2
SWIM
10-3
Chebyshev
ELM
Fourier function
polynomial
Fourier
10-4
Chebyshev
0.5
u(t.x)
0.0
u(t,x)
10-5
-0.5
10-6
-1
0
1
-1
50
100 150 200
X
Number of basis functions
Figure 22: Approximation of Burgers' equation solution at 0.99 with four types of basis functions. The number of basis functions in all cases 102. Oscillations can be seen near the steep gradient for the methods using ELM ba- sis functions, Fourier functions, and Chebyshev polynomials.
Figure 23: Approximation error for four types of basis functions. Here, we directly fit the Burgers' equation solution at t 0.99. The approxima- tion error decreases as we increase the number of basis functions, and the SWIM basis functions yield the best result among all methods.
on a complicated geometry inspired by tree-like pattern occurring during the controlled shaping of fluids Islam & Gandhi (2017). The initial condition and time-dependent Dirichlet boundary conditions obtained from the constructed solution of the PDE
(x, y, t) 5c 5t (in(( 3 (x. EQ, t [0.
(29c)
The training is performed on 1500 data points in the interior and boundary. We test the neural-PDE solvers with 5000 data points in the interior and on the boundary. The weights of the hidden layer for the Frozen-PINN-ell distribution in [-1, 1]. For our approach to handling time-dependent Dirichlet boundary conditions, please refer to Appendix B.2.5] The hyperparameters for various neural PDE solvers are outlined in Table 18 Figure 24 shows the mesh generated for the FEM and the sampled collocation points for the neural PDE solvers. For the mesh consider for this problem (see Figure 24). we could not improve the accuracy further with FEM by using higher-order polynomial basis functions. While mesh refinement is possible, it's time-consuming, and our method avoids this by working directly with point clouds.
are sampled from the Gaussian distribution and biases from a uniform
The ablation study for the number of neurons in the hidden layer of the network for Frozen-PINN-elm and Prozen-PINN-swin is presented in Table [9 For PINN, the results for the ablation studies for the width of the network and the number of data points are included in Table 20 Table 21] Additionally, we perform an ablation study for the SVD layer demonstrate its impact on the computation time saved in Table 22 Particularly, we observe that with the SVD layer, the number of basis functions (width after the SVD layer) is reduced by up to 22x for Frozen- PINeen- and up and we obtain substantial speed-ups (more than a factor of 50) in the
Ablation studies:
to 1.5x for computation time.
Comparison of results: The comparison of training times and errors is presented Table [ Figure 25 shows the ground truth and Figure shows the error plots with all approaches.
47

---

Published as a conference paper at ICLR 2026
*
#
(a) Generated Mesh: FEM
Sampled collocation points: Neural PDE solvers
Figure 24: Advantages of free methods: (a) For mesh-based methods, complicated mesh must be constructed, whereas (b) for neural PDE solvers, one can easily sample arbitrary points in the interior (blue) and on the boundary (red) of the domain and work directly with point clouds.
Table 18: Non-linear diffusion equation (see Table I Summary of hyper-parameters.
Parameter
Value
Frozcu-PINN-elmn (low-procision)
Number of hidden layers
(scalinear and SVD layer)
Hidden layer wilth Activation L -erggllaaiiaiionn
350
tanh
11
SVD cutoff
x
10-11
ODE solver tolerance Loss
10-8
mean-squared emorr
boundary condition strategy
augmented ODE
Pozzen-PNNN-siim (high-precision)
Number hidden layers
(nonlinear SVD layer)
Hidden layer width Activation
500
tanh
[Frrgulaaizzation
SVD cutoff
10-18
ODE solver tolerance
10-8
Loss
mean-squared error
boundary condition strategy
augmented ODE
FEM
Number cotities
154
Number of nodes
1193
Number of elements
2070
Type elements
Lagrange triangle
Shape of clements
Degree polynomials
Number of basis functions
1193
Solver
Newton solver
Timestep size
0.001
PINN
Number hidlden layers
Layer width
[10, 20, 30, 40) Lanh
Activation
Optimizer
LBFGS ADAM 10000 mean-squared
Epochs
Loss
Leurning rate
0.01
Batch
1000
Parameter initialization
Xavier Gloret & Bengio 2010
Loss weights. Ar, Aa
0.01,
* Interior points
[8790, 1760, 880, 440]
# Initial and boundary points
(11-10, 630, 320, 160)
48

---

Published as a conference paper at ICLR 2026
Ground truth
t 0
t= 0.3
t= 0.6
t 1
0.9
0.8
7.5
0.7
5.0
0.70 0.70.80.9
0.7 0.70.80.9
0.70.80.9
0.70.80.9
Figure 25: Non-linear diffusion equation: ground truth.
Absolute Frozen-PINN-swiim
0.3
0.9
0.8
0.7
2
2.5 5.0
0.5 1.0
1.5
2 3
1e-9
1e-6
1e-5
le-3
(a) Absolute Prozen-PINN-soee
Absolute lellen-
0.9
0.8
0.7
0.7
0.9
C.7
0.7
0.8
0.9
0.7
0.8
0.9
X1
X1
X1
2
0.5 10 1.5
1e-4
le-3
(b) Absolute eтror: Frozen- Frozen-PINN-ooe PINN-elm.
Absolute error: FEM
t=
0.3
= 0.6
t=1
0.9
0.7
0.7
0.8
0.9
0.7
0.8
0.9
0.7
0.9
0.7 0.8
2
0
2
0
2
5
1e-1
1e-1
1e-1
1e-2
(c) Absolute emror: FEM.
Absolute error: PINN
0
t 0.3
0.6
t=1
0.8
0.7
0.8
0.9
0.7
0.8
0.9
0.7
0.8 X
0.9
0.7 0.8 0.9
X
0.5
7.5
2.5 5.0 7.5
1e-1
1e-1
1e-2
1e-2
(d) Absolute emor: PINN (LBFGS).
Figure 26: Non- -ineli diffusion equation: absolute error plots and ground truth at four-time instants.
49

---

Published as a conference paper at ICLR 2026
Table 19: Non-linear diffusion equation: ablation study for the network width for Frozen- and Frozen-PINN-elm. The mean is computed over 3 seeds.
Width
Relative I,2 error (Frozen-PINN-swimm
Relative 1,2 error (Frozen-PINN-elmm)
200
1.34e-4
4.92e-3
300
5,07e-6
3.13e-5
400
2.88e-6
1.02e-5
500
3.02e-7
1.52e-5
Non-linear diffusion equation: hyperparameter optimization for PINN varying layer width. The mean is computed over seeds.
Table 20:
Layer width
Training time (s)
RMSE
Relative L error
10
61.09 + 1.62
4.11e-2 + 2.04e-3
1.50e-2 7.48e-4
20
68.05 + 1.56
3.74e-2 +: 1.04e-3
1.37e-2 + 3.82e-4
30
76.01 + 0.57
3.67e-2 = 1.03e-3
1.34e-2 3.78e-4
40
82.43 + 0.45
3.76e-2 = 1.69e-3
1.37e-2 + 6.21e-4
Table 21: Non-linear diffusion equation: hyperparameter optimization for PINN varying interior
points.
Interior points
Training time (s)
RMSE
Relative L error
600
65.08 + 4.23
3.74e-2 + 1.04e-3
1.37e-2 + 3.82e-4
1200
98.48 3.78
3.51e-2 6.67e-4
1.28e-2 = 2.44e-4
2390
143.31 +: 5.50
3.34e-2 + 6.53e-4
1.22e-2 +: 2.38e-4
Table 22: Non-linear diffusion equation: Ablation study of the SVD layer in Frozen-PINN-swinn and Frozen-PINe PINN-elm. We report for runtimes exceeding 3 hours. Two variants of Frozen-PINe PINN--INNnn--e-- are shown: Prozen-PIIN-elmm eecmuccceette (higher accuracy, longer runtime) and Frozen- PINN-elm (lower runtime, with error comparable or better than PINNs, enabling fair comparison). The ratio of the hidden layer width to the SVD layer width is denoted by Cr-
Method
Quantity
With SVD layer
Without SVD layer
Ratio
Trozen--IINN-elm-ccccraat
Width
62
300
C Rv 22.8x Speed-up =o 52x
Time (s)
60.98
7087.38
Rel. La emrr
6.49c-8
1.02c-6
Frozee-PINN-ellm-fatt
Width
35
300
C, R 8.5x Speed-up no
Time (s)
30.57
Rel. La error
5.12e-5
Frozen-PINN PINN-
Width
316
1.5x Speed-up no
Time (s)
328.00
0c
Rel. L2 error
2c-6
C.6 NONLINEAR REACTION-DIFFUSIN EQUATION
Problem Setup: The non-linear reaction-diffusion equation models biological pattern forma- tion, such as Zebra stripes, fish spots, in myriad chemical reactions, and flame propagation during
combustion (Britton 1986 Lam & Lou 2022).
In this benchmark from Zang et al. (2000) we consider a five-dimensional nonlinear diffusion equation given by
U1
Au u² f(z,t),
rEΩC C Rd
t [0, 1].
(30a)
f(x 1) = -17² 2) sin
π
cos
4sin² 2
cos 52
(30b)
2
50

---

Published as a conference paper at ICLR 2026
on the domain Ω = 1|, The initial condition and time-dependent Dirichlet boundary conditions are obtained from the constructed solution of the PDE
u(x, t) 2in
π
COS
π 22 e
(31)
Note that the solution is independent of three out of five dimensions. The training is performed on 1000 data points in the interior and 1000 data points on the boundary. The test data set is generated the same way as in Zang al. (220)) to evaluate the weak adversarial networks. In particular, compute the error in the 5-dimensiona domain, we use a mesh of size 100 x 100 for the two coordinate directions in which the solution changes (here, (21, z22) and uniformly randomly sample the other coordinates (here, (zz. 24, σ5)) in the domain. The hidden layer weights for the -INzenn sampled from the standard Gaussian distribution and biases from uniform distribution in [-1,1] 1]. Please refer to Appendix B.2.5 for our approach to handling time-dependent Dirichlet boundary conditions.
The sampling strategy for basis functions described in Section 3.7 using projected pairs of data points substantially improves efficiency and accuracy. Our approach requires 20 times fewer training points in the interior compared to Zang et al. (2020) while simultaneously achieving a relative L error more than two orders of magnitude lower.
Details on the sampling well-oriented basis functions: For each pair of collocation points in the spatial domain (1),, (2) E Ωt we project the vector x2(2)_ xx onto the two-dimensional hyper-plane spanned by the gradient of the initial solution at xx2) and use the projected points as the new pair of points a .x(2) E Q. Since f(()) - _i(1) always points in the direction of the gradient of the initial solution, this allows the SWIM algorithm to embed directional information into basis functions,
unlike PINNs and ELMs, which lack this control. This idea is illustrated in Figure 5b
Ablation studies: The ablation study for the number of neurons in the hidden layer of the network for Frozen-PINN-elm and Frozen-PINN-swin is presented in Table 23 We further validate the efficiency of sampling basis functions using projected pairs of data points with the Frozen-PINNN swim approach by performing an ablation study varying the number of internal collocation points in Table 24 Our results show that using just 1,000 data points achieves training errors that are nearly identical to those with 20,000 points. This highlights the effectiveness of the projection trick in reducing the need for excessive collocation points, thereby significantly lowering computational cost without compromising accuracy. Additionally, we perform an ablation study for the SVD layer to demonstrate its impact on the computation time saved Table 25 We observe that with the SVD layer, the number of basis functions (width after the SVD layer) is reduced by up to 1.57x for Prozen-PINN--simm, and we obtain substantial speed-ups by a factor of in the computation time.
Table 23: Non-linear reaction diffusion equation: ablation study for the network width for Frozen- PINN-swim and Frozen-PINN-elm. The mean is computed over 3 seeds.
Width
(with projection)
-INzen-PNNN-wwim
Frozen-PINN-on
100
1.44e-4
7.65e-2
2.08e-1
400
9.99e-5
1.75e-2
6.37e-2
700
9.92e-5
8.72e-3
3.65e-2
1000
9,87e-5
5.70e-3
2.58e-2
2000
9.86e-5
8.62e-3
1.67e-2
4000
9.86e-5
9.98e-3
3.68e-2
Table 24: Non-linear reaction diffusion equation: ablation study for the number of interior collocation points for Frozen Srizen-PINNN-eee PINN-swim. The mean is computed over seeds, and the network width is 400.
Interior points
Training time (s)
RMSE
Relative L error
1000
12.43
3.68e-5 +- 6.78e-12
9.99e-5 +- 1.84e-17
2000
102.74
3.67e-5 +- 3.00e-10
9.90e-5 +- 8. 1.0
20000
689.81
3.63e-5 +- 2.02e-9
9.87e-5 +- 5.49e-09
51

---

Published as a conference paper at ICLR 2026
Table 25: Non-linear reaction diffusion equation: Ablation Study for the SVD layer with Frozen- PwNN-swim.
With SVD layer
Without SVD layer
Ratio
Number of neurons
254
400
Width Compression 1.57x
Time (s) Rel. Lg error
12.86
53.51
Speed-up v 4.1x
9.99e-5
9.99e-5
The exact architectures and comparison training times and errors are presented in Table and Table 27 We observe that Frozen-PINN-swwnn with the projected pairs of points Piizze--INN--wim-pp far outperforms all the other approaches by around orders of magnitude, while simultaneously being 50 times faster. Figure 27 shows the errors with all approaches and the ground truth.
Comparison of results:
Table 26: Non-linear reaction diffusion equation: Summary of hyper-parameters.
Parameter
Value
PINN
Number of hidden layers
Activation
tanh
Optimizer
LBFGS & ADAM
Epochs
10000
squared emorr
Leamning rate
0.001
Batch
Parameter initialization
Xavier (Glorot Bengio 2010
Loss weights, A1, Az Interior points
0.01,
1000
Initial and boundary points
1000
Prozea-PINN-swing
Number of hidden layers
(nonlincar and SVD layer)
Halden layer width Activation 1.2 r.-rrgularization
1000
tanh
10 10
SVD cutoll
10 10
ODE solver tolerance
10
Loss
mean-squared
boundary condition straregy
augmented ODE
Frezea-PINN-elm
Number of hidden layers
(nonlinear and SVD layer)
Halden layer width Activation 1,2 I.-regularization
2000
10 10
SVD cutoff
10 10
ODE solver tolerance
10-
Loss
mean-squared
boundary condition straregy
augmented ODE
Mrozen-PNNN-wwin (with projection)
Number of bidden layers
(nonlinear and SVD layer)
Hidden layer width Activation
700
I.-regularization
10
SVD cutoff
10 10-10
ODE solver tolerance
10-4
Loss
mean-squared
boundary condition strategy
augmented ODE
Table 27: Mon-linear reaction diffusion equation: Summary of results.
Method
Train time (8)
RMSE
Relative eror
architecture
PINN (ADAM)
171.43
1.25e-1 1: 6.60e-3
3.40e-1 1.79
(6. 40. 1)
(LBFGS)
3.33e 1: 1.54e-2
(6. 4x.4,
Prozen- Prozen-PINN-chm
621.2
6.17c-3 + 2.02c-4
1.67e22
5.49c-4
(5. 2000,
Procn-PPINN-swin
117.24
2.09c-3 + 1.9ee5
5.70c-3 5.19e-5
(5. 1000.
Frozen-PIhe PPNN--si (projection)
12.43
3,67e-5 :: 2.28c-9
9,99e-5 6.21e-9
(5. 7007
52

---

Published as a conference paper at ICLR 2026
Ground truth (test)
1
1
0
-1
0
1
1
(a) Ground truth.
Absolute error: Mrozen-PINN--Ce
0.3
O
0
-1
-1
1
X1
X1
X1
2.5
5.0
0.5 1.0
0.5 1.0
2.5 5.0
1e-6 6
1e-4
1e-5
(b) Absolute error: Frozen- PINN-swim (with projected pairs of data points).
Absolu Frozen-PINN-ellm
-1
X1
X1
X1
2.5
5.0
0.5
10
0.5
1.0
2.5
5.0
(c) Absolute emror: Frozen Mrizen-PINN-elm...
Figure 27: Non-linear reaction diffusion equation: absolute error plots and ground truth at four time instants.
53

---

Published as a conference paper at ICLR 2026
Reference
Frozen-PINN-swim
Difference
5.0
5.0
5.0
2.5
2.5
2.5
0.75
1.00
0.75
1.00
t
-2.
0.0
2.5
2e+00
Se+00
Figure 28: The Kuramoto-Sivashiee Sivashinskyy equation: Reference solution, Frozen-PINN solution, and
difference between two solutions.
C.7 KURAMOTO- SIVASIIIS EQUATION
Problem Setup: The
equation is a fourth-order nonlinear PDE, which
models the instabilities in flames and exhibits chaos. We consider the equation from Hao al. 2024 of the form:
+ +owu + Buaz YUzzan 0, E S2 L E,
(32a)
with the parameters
100
100
100
16
3
162
Y
164
(32b)
The domain is Ω x T [0, 2x
x [0, 5] for the experiment in Figure 6 or Ω x T [0, x [0, 1]
for the experiment in Figure 28
The initial condition
u(x, 0) cos(x)(1 + sin(x)),
E ?
(32c)
and we apply the periodic boundary condition.
Comparison of results: In Figure 28 we compare our Frozen-PINN solution with the reference
solution in Hao et al. (2024). The training time of Frozen- is 2.7 seconds on CPU (averaged over seeds). We also solve the PDE for longer time span, as shown in Figure 6 The architecture of the model can be found in Table 28
Table 28: Kuramoto-Sivashinkt equation: Hyper-parameters for the result in Figure Figure 28
Mrozen-PINN-swiid
Number of hidden layers
2
Hilden layer width Outer hasis functions Activation 1,2 regularization
2000
100
tanh
10-12
Loss boundary condition strategy
mean-squared error
boundary-complian layer
C.8 DIGH-DDMENSIONAL DIFFUSION EQUATION
Problem setup: High-dimensional diffusion plays an important role in various fields, including
image processing, finance, and quantum mechanics (Sapiro 200T Janssen al. 2013 Nagasawa
2012).
We consider two benchmarks for the high-dimensional diffusion equation. the first case, following
WWag & Dong 2024), we solve the diffusion equation defined over the domain Z = [-1, 1] and
time interval (0, I[, for dimension d (3.5, 5, 7, 10} given by
Au (2 -1-1) cos
Σ
z 2=11
exp (-t), JEΩ, t [0,
(33)
d
54

---

Published as a conference paper at ICLR 2026
with the exact solution given by
1 Σ
u(x, t) cos
exp (-t).
(34)
The initial and boundary conditions are derived from Equation (34). For this high-dimensional diffusion equation, use 16000 training points in the interior and 4000 points the boundary
randomly sampled using the Latín hypercube strategy. The test data contains 8000 points in the interior and 2000 points on the domain's boundary, which were also sampled with a Latin hypercube strategy.
For the second benchmark from He al. (023). we consider a 100-dimensional heat equation for x E [00, 1). E (0, 1) given by
U = Au, x E B(O, 1). E [0,
(35a)
u(x, 0) = ||^||
xE B(,,
(35b)
2N
1(x,1) = + 2N x E @B(0, 1), [0
(35c)
where the true solution is
||]||
u(a, t) =
(36)
2d
The value of d is 100 and represents the dimension of the PDE. To solve the 100-dimensional heat equation, we generate 1000 interior and 1000 boundary training samples using Latin hypercube sampling. The test dataset comprises 8000 interior points and 2000 boundary points, also selected via Latin hypercube sampling.
Extended discussion on results: Note that, in general, is extremely hard accurately represent arbitrary 100-dimensional functions with a few hundred basis functions unless the solution is already in their span (e-g, approximating a linear solution with linear bases). The 100-dimensional heat equation benchmark from He et al. (2023), indeed, admits true solution with very shallow gradients in space that varies linearly in time (see Equation (36)). Although this benchmark technically admits quadratic analytical solution, rendering second-order polynomials natural fit for this particular example, higher-order approximations, in general, quickly become infeasible due to the curse of dimensionality. For instance, cubic approximation already requires millions of basis functions. A natural alternative could be to use lower-order upproximations like linear regression. However, one linear bases to solve the diffusion equation, they cannot capture temporal dynamics because linear basis functions are harmonic (the Laplacian is zero). By contrast, our approach provides mildly non-linear bases that have non-zero Laplacians, facilitating "almost linear" approximation each point in time.
The fact that Frozen-PINN-ell with single hidden layer yields a significantly accurate and faster ap- proximation compared to PINNs with multiple hidden layers trained with classical pack-propagation reveals an interesting observation that one does not necessarily benefit from using deeper neural net- works. While stochastic and iterative training methods might eventually identify suitable parameters, the highly non-linear, non-convex loss landscape makes such optimization particularly challenging.
Due to the smoothness and lack of steep gradients in the solution of the PDE, Frozen-PINN-eld more suitable for approximating the solution of the chosen PDE and is one three orders of magnitude more accurate than vanilla Prozen-PINN-swim, as one would expect (see Section 2.2
The ablation study with respect to the network width for Frozen-PINN-cin and Prozen-PINN-swir is already presented in Figure 7 where observe rapid exponential decay of error with respect to increasing width of the network (even exponential convergence for the high-dimensional diffusion equation in 3 and 5 dimensions).
Ablation studies:
The hyperparameters for all neural PDE solvers considered in this work for the 10-d heat equation and the 100-d heat equation are presented in Table 29 and Table 30 respectively. The results for up
55

---

Published as a conference paper at ICLR 2026
to 100-dimensional diffusion equations summarized in Table 31 Please refer to He et al. (2023) for details on hyperparameters for PINNs for the 100-dimensional heat equation.
The results of the ablation study for the SVD layer with the high-dimensional diffusion equations for different dimensions are presented in Tables 32 and 33 We observe that for Piitn--PNNN-liii the SVD layer results in substantial speed-ups for 3, and 7 dimensional heat equations by factors of 52, 77, and 21, respectively. We observe that the compression ratios achieved with the SVD layer are also substantial 22.8, 5, and 1.2, for dimensions 3, 5, and 7, respectively. For the 10-dimensional diffusion equation, to cover the high-dimensional space, we observe (relatively lower compared to other dimensions) compression ratio of 1.4, as more basis functions are required to represent functions in high dimensions accurately. Thus, the time required with the SVD layer is around 6 percent less than the time required without the SVD layer. In all the cases, the loss always in the same order as the one without the SVD layer.
Note that in all cases, the extra cost of computing the SVD easily pays off by substantially saving time in the ODE solver for Frozen- feature matrix and the reduction in the size of the ODE system to solved. With Frozen-PINN- swim, the observations are similar but with lower compression ratios and speed-ups. But, for this problem,
This is because of the improved conditioning of the
is not the preferred method, as the underlying solution is smooth, has low-frequeney spatial variations, and does not have steep gradients anywhere in the domain. Thus, SWIM basis functions are not optimal in the vanilla setting. See Appendix 1.2.1 for details on this.
Comparison of results: We demonstrate that Frozen-PINN-ell accurately solves the 10- dimensional and 100-dimensional heat equation by visualizing the time evolution of the solution at some sampled points in space in Figure 29 and Figure 30 in certain dimensions.
Table 29: Summary of hyperparameters for the 10-dimensional diffusion equation.
Parameter
Value
Frozen:
Number hichlen layers
(nonlinear layer)
Hiöden layer width Activation
4900
tanh
1.2 regerualarizttonn
10-10 10
SVD cutoff
10
ODE solver tolerance
10-
Loss
mean-squared ervor
boundary condition strategy
augmented ODE
Mroven-PNNN.oo elm (low- precision)
Number hidden layers
(nonlinear and layer)
Hidden layer width Activation
400
tanh
E. regullaarizatonn
10-5
SVD cutoff
10
ODE solver tolcrance
10-4
parameter range I-Tms Pim Loss
I'm 0.05
squared
boundary condition strategy
augmented ODE
Prozeo-PINN-elmn (high-precision)
Number of hidden layers
(ocolincar and SVD layer)
Hidden layer width Activation L reguullarizttion
4900
tanh
10 10
SVD cutoff
10-10
ODE solver tolerance
10-6
parameter range [- Fmrl
Tun 0.05
Loss
mean-squared crror
boundary condition strategy
augmented ODE
PINN
Number of hidden layers
Layer width Actisation
20
tanh
Optimizer
LBEGS (ADAM)
Epochs
1000 (5000)
Loss
mean-squared orror
Learning rate
0.1
Batch size
4000
Parameter initialization
Xavier
Glorat Benpio 2010
Loss weights, A1, Aa
1. 1
# Interior points
16000
# Initial and boundary points
4000
56

---

Published as a conference paper at ICLR 2026
Ground Truth
Frozen-PINN-eile
Absolute Error
1.0
1.0
1.0
0.00010
0.75
0.5
0.5
0.00005
0.50
0.50
0.0
0.0
0.0
X3
X3
X3
(a) 23-4 plane.
Ground Truth
Trozen-PINN-elme
Absolute Error
1.0
1.0
1.0
0.00015
0.00010
0.75
0.75
0.5
0.5
0.5
0.00005
0.50
0.50
0.0
0.0
0.0
-1
1
0
1
-1
Xs
(b) Zo-t plane.
Ground Truth
Trozen-PINN-elm
Absolute Error
1.0
0.00010
0.75
0.75
0.5
0.5
~ 0.5
0.00005
0.50
0.50
0.0
0.0
0.0
-1
-1
-1
Xy
X7
X7
(c) Ø7-t plane.
Figure 29: 10-dimensional diffusion equation: Ground truth, Frozen-PINN-elm solution, and point- wise absolute eиror at various planes different time points. The rest of the spatial coordinates are set to the center of the spatial-tempont domain.
57

---

Published as a conference paper at ICLR 2026
Ground Truth
Frozen-PINN-ein
Absolute Error
1.0
1.0
1,0
1.0
1.0
0.0055
0.5
0.5
0.5
0.0050
0.0
0.0
0.0
-0.25
0.00
-0.25
0.00
-0.25
0.00
X1
X1
X1
(a) #1-t plane.
Ground Truth
Frozen- PINN-elm
Absolute Error
1.0
1.0
1.0
1.0
1.0
0.0055
0.5
0.5
0.5
0.5
0.5
0.0050
0.0
0.0
0.0
0.00
0.25
0.00
0.25
0.00
0.25
X21
X21
X21
(b) Ø21-t plane.
Ground Truth
Pozzen-PINNeeiee
Absolute Error
1.0
1.0
1.0
0.0055
0.5
0.5
0.5
0.5
0.5
0.0050
0.0
0.0
0.0
0.25
0.00
- 0.25
0.00
0.25
0.00
X100
X100
X100
(c) Fion-t plane.
Figure 30: 100-dimensional diffusion equation: Ground truth, Frozen-PINNNellm solution, and point- wise absolute eтror at various planes different time points. The rest of the spatial coordinates set to the center of the spatial-temporal domain.
58

---

Published as a conference paper at ICLR 2026
Table 30: Summary of hyper- parameters for the 100-dimensional diffusion equation.
Parameter
Value
Froven Proven-PINNN PiNN-swim
Number of hicken layers
(nonlinear and SVD layer)
Hidden layer Activation
200
tanh
L -regularization
10-5
SVD cutoff
10
ODE solver tolerance
10-4
Loss
squard- error
boundary condition stralegy
augmentted ODE
Number of hiclen layers
(nonlinear and SVD layer)
Hidden layer width Activation L -regularization
125
tanh
10
SVD cutoff
10-4
ODB solver tolerance
purameterr range
0.05
Loss
mean squared emor
boundary condition stratogy
augmented ODE
Table 31: Summary of results for high-dimensional diffusion equation. We denote the Frozen- PINN-elm results in the low-precision and high-preciste regimes with Frozen-PINN- and
respectively.
Dimension
Method
Time (s)
RMSE
Relative L error
3-d
PINN (LBFGS)
102.32
2.84c-4
3.73e-5
4.54c + 5.97c-5
Frozen-PINN-swinn (our)
95.73
2.18e-6 1.93e-6
5.37e-6 4.27e-7
Frozen-PINN-elm-te (our)
2.42e-6
1.37e-6
3.90e-6 2.98e-6
Frozen-PINN-cm-accrra (our)
60.98
3.48e-8
2.17e-6
6.49c-8 4.31c-8
5-d
PINN (LBFGS)
133.95
2.91e-4 = 5.34e-5
4.52e-4 8.30e-5
Frozen-PINN-see PINN-swim (our)
129.65
1,03e-4 = 5.94e-5
2.39c-4 8.69e-5
Frozen-PINN-elm- (our)
1.2
1.25e-4 2.42e-5
3.74e-4 5.37ee-5
Frozen-PINN-elm-accurat (our)
102.95
4.71e-7 3.56e-7
7.5e-7 = 33.22-77
7-d
PINN (LBFGS)
163.89
3,05e-4 2.94e-5
4.69e-4 44.5le5
Frozen-PIND PINN-swim- (our)
198.20
3.96e-4 = 1.03e-4
7.8e-4 = 2-500e4
Frozen-PINN-elm-fff (our)
1.05e-5 8.76c-6
2.21e-5 1.01e-5
Frozen-PINNellm-acuraa (our)
176.95
1.19e-6 2.93e-7
2.54e-6 5.10e-7
10-d
PINN (LBFGS)
189.67
3,98e-4 - 6.59e-5
6.06e-4 1.00e-4
Frozen-PINN-swinn (our)
61.07
1.0le-3 3.09e-4
2.31e-3 1 1.03e-3
Crozen-PINN-elmte (our)
2.07
2.89e-4 5.91e-5
4.46e-4 1 9.61e-
Frozen-PINN-ell--cccrattt (our)
182.91
1.4e-5 -3.22-6
2.28e-5 1 5.91e-6
100-d
Vanilla PINN (He al.
2023
141
6.00e-3
PINN (He et al, (2023)
49.8
6.30e-3
Frozen-PINN-swwe (our)
68.39
1.00e-3 = 1.75e-5
1.71e-3 + 3.0le-5
Frozen-PINN-ee (our)
5.24
2.40e-4 = 9.92e-6
4.12e-4 + 1.70e-5
59

---

Published as a conference paper at ICLR 2026
Table 32: High-dimensional diffusion equation: Ablation Study for the SVD layer with Frozen- P-NN- swim.
Dimension
Quantity
With SVD layer
Without SVD layer
Ratio
3-d
Width
1391
4000
Compression s 2.9x
Time (s)
95.73
388.12
Speed-up ~ 4x
Rel. L2 enror
5.29e-6
4.77e-6
5-d
Width
1437
4000
Compression 2.8x Speed-up = 1.5x
Time (s)
129.65
199.92
Rel. L2 error
2.39e-4
2.18e-4
7-d
Width
3114
4000
Compression 2 1.3x
Time (8)
120.32
198.31
Speed-up 1.6x
Rel. 12 error
7.83e-4
7.83e-4
10-d
Width
3100
4000
Compression R 1.3x Speed-up Rs 0.91x
Time (s)
121.93
111.8
Rel. L2 error
2.30e-3
2.30e-3
100-d
Width
200
200
Compression 1x Speed-up Rs 0.97x
Time (s)
5.24
5.13
Rel. Lg error
3.82e-3
3.82e-3
Table 33: High-dimensional diffusion equation: Ablation Study for the SVD layer with Frozen-PINNy
elm.
Dimension
Quantity
With SVD layer
Without SVD layer
Ratio
3-d
Width
175
4000
Compression < 22.8x
Time (s)
60.98
7087.38
Speed-up Rs 52x
Rel. L2 enrrr
6.49e-8
1.02e-6
5-d
Width
794
4000
Compression 5x Speed-up z 77x
Time (s)
89.27
6873.8
Rel. L2 error
7.30e-7
2.19c-6
7-d
Width
3336
4000
Compression ss 1.2x
Time (s)
176.95
3770.09
Speed-up x 21x
Rel. Le error
2.54e-6
4.06e-6
10-d
Width
2856
4000
Compression ~ 1.4x Speed-up Rs 1.06x
Time (s)
119
127
Rel. La error
5.57e-5
4.36e-5
100-d
Width
552
600
Compression & 1.3x Speed-up Rs 0.96x
Time (s)
68.39
71.38
Rel. Le error
1.71e-3
1.71e-3
60


Reviews:
[{"review_id": "5sdu8wL5Pp", "reviewer": "Reviewer_19cn", "summary": "The paper proposes Frozen-PINN, a gradient-descent-free approach for time-dependent PDEs that enforces temporal causality by separating space and time. Hidden-layer spatial bases (weights and biases) are sampled and then frozen (via either data-agnostic ELM sampling or data-dependent SWIM sampling), while the time-dependent output coefficients are evolved by solving an ODE obtained by inserting the frozen spatial basis into the PDE. Initial coefficients are set by least squares and boundary conditions are handled via either a boundary-compliant layer or an augmented ODE. The method further compresses the neural bases with a truncated SVD layer to reduce ODE dimension and accelerate computation. Empirically, the authors evaluate Frozen-PINN on eight PDE benchmarks and show improved accuracy compared to many PINN variants and classical mesh solvers in the studied regimes. The manuscript additionally provides ablation studies (e.g., on the SVD layer and SWIM vs ELM sampling), implementation details, and a reproducibility statement.", "questions": "1. What is the computational cost and memory scaling of SWIM resampling and basis projection in high dimensions (d≥10)? Are there practical heuristics (e.g., maximum number of collocation points, subsampling strategies) you recommend to keep SWIM tractable?\n\n2. Could you provide more explicit guidance (or an automated procedure) for choosing the SVD cutoff/target rank so a new practitioner can reliably obtain the claimed speedups without extensive manual tuning? Please include how SVD choice impacts stability of the ODE solver.\n\n3. For PDEs with truly complex spatial structure (e.g., 2D Navier–Stokes), do you expect Frozen-PINN to require domain decomposition or local basis strategies? Do you have preliminary results or runtimes/cost models for such extensions?\n\n4. It will be interesting to see efficiency comparisons with numerical methods for solving PDEs.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "1. The paper introduces a clear and intuitively appealing space–time separation idea: freezing spatial basis functions and evolving only the time-dependent output layer via ODE solvers, which enforces temporal causality by construction and directly addresses a known PINN weakness (treating time as an extra spatial dimension).\n\n2. The manuscript presents two pragmatic sampling strategies for spatial bases: (i) data-agnostic Extreme Learning Machine (ELM) sampling and (ii) data-dependent SWIM sampling (Sample-Where-It-Matters), where SWIM uses pairs of collocation points (and optional projection onto solution gradients) to place basis functions near shocks or informative regions. This design is well motivated and demonstrated to be useful for handling shocks and directional features.\n\n3. The authors implement an SVD compression layer to orthogonalize and truncate the neural basis, substantially reducing the dimension of the resulting ODE system and lowering computational cost; the paper supplies ablation evidence that the SVD layer can yield very large runtime improvements in some settings (e.g., tens of speedups in several ablations).\n\n4. The experimental evaluation is broad and systematic within the chosen set of PDEs: eight benchmarks spanning advection extremes, shocks, higher-order spatial derivatives, complex geometries, and high dimensionality are considered, and the authors report means and standard deviations over multiple seeds with numerous ablations (resampling, basis widths, SVD thresholds). This breadth strengthens the empirical narrative.\n\n5. The method removes reliance on iterative gradient-descent training for the bulk of parameters, enabling very fast CPU training in many reported cases and promising a low-resource alternative to GPU-centric PINN training; the authors emphasize carbon-efficiency and implementation simplicity as benefits.\n\n6. The paper is reasonably reproducible in spirit: many methodological details, hyperparameter settings, and appendices are provided, and the authors pledge to release code and seeds to reproduce experiments.", "weaknesses": "1. Regarding the experimental claims of extreme speedups and accuracy (e.g., “up to 100,000× faster training” and multiple orders-of-magnitude accuracy gains), it should be carefully noted that these comparisons are across different CPU vs GPU settings (CPUs for Frozen PINNs vs GPUs for PINNs), and different solver regimes (low-precision vs high-precision). The paper does acknowledge benchmarking rules (matching accuracy regimes and noting GPU timings), but the mixed hardware/configurations needs to be specified while making broad claims of speedups and accuracy (which could be attributed to the differences in configurations too).\n\n2. The scope of theoretical guarantees and assumptions is limited in interpretability for practitioners. While the paper references and partially justifies why random-feature style bases and ODE evolution reduce optimization complexity, it lacks clear, practically checkable conditions that explain when SWIM/ELM sampling will produce sufficiently expressive bases (e.g., how many bases are needed as a function of Kolmogorov n-width or solution regularity). The theoretical framing would benefit from tighter bounds or guidance connecting basis counts, SVD truncation thresholds, and expected approximation errors in realistic regimes.\n\n3. The overhead and robustness of repeated resampling in high dimensions or with sparse data should be quantified and discussed more explicitly. Some ablations show SVD reduces runtime dramatically, but the cost of SWIM itself (and its scaling w.r.t. dimension and collocation count) needs clearer accounting.\n\n4. The generalization to fully spatially complex PDEs (e.g., Navier–Stokes turbulence or domains requiring domain decomposition) is asserted but not demonstrated. Claims about matching mesh solvers in low-dimensions and scaling to high dimensions are promising but would be stronger with at least one medium-scale 2D/3D fluid example or a domain-decomposition demonstration.\n\n5. The method’s dependence on choices of basis sampling (ELM vs SWIM), SVD cutoff, and collocation density appears substantial from the ablations, but sensitivity studies across a broader hyperparameter grid and more seeds would increase confidence in robustness and reproducibility. Some ablations show very large speedups only after aggressive SVD compression. This raises questions about stability and about how to choose compression thresholds in new problems.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "0zDfPyGhyT", "reviewer": "Reviewer_6QHM", "summary": "The paper proposes Frozen-PINN, a physics-informed neural network framework that avoids gradient-descent altogether by freezing space-dependent features and evolving only time-dependent coefficients through ODE solvers, thereby enforcing temporal causality by design. It argues that two root causes of slow/inaccurate PINNs are (i) hard, multi-objective, nonconvex optimization over many parameters and (ii) treating time as just another spatial dimension, and shows that space–time separation plus loss decoupling sidestep both. Frozen-PINNs train orders of magnitude faster, while reaching accuracies comparable to or better than SOTA PINNs and even close to mesh-based solvers in low dimensions. The method also introduces an SVD layer to reduce stiffness and further accelerate solving.", "questions": "You propose two strategies: boundary-compliant layer vs. augmented ODE. How should a practitioner decide between the two for a new PDE/domain? Is there a cost/accuracy tradeoff (e.g. augmented ODE enlarges the system and can worsen stiffness; boundary layer needs a problem-specific A)? A decision table or rule-of-thumb would make this much more usable.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper removes SGD/backprop from the loop by freezing spatial features and evolving only time-dependent coefficients with classical ODE solvers and least squares, which is a non-standard but very clean way to “solve” the training bottleneck of PINNs rather than tuning optimizers.\n\n2. By constructing the solution as frozen spatial bases plus time-evolving coefficients, the method enforces the Markov/causal structure of time-dependent PDEs “by construction,” something many domain-decomposition / curriculum PINNs only approximate with scheduling. That’s an insightful modeling decision with high originality.\n\n3. The 10D–100D heat-equation results and the comparison plot (Frozen-PINN-elm vs. PINN) show the method keeps errors small while increasing width, a case where many neural PDE solvers and all classical meshes struggle; this is important significance-wise for scientific ML.\n\n4. Because training is reduced to classical ODE solves and least squares on frozen bases, the method does not lean on big GPUs or specialized accelerators, which broadens the audience and strengthens the paper’s practical significance.", "weaknesses": "1. Missing literatures regarding model-architecture-based temporal PINN methods. In the introduction section, it is claimed that Non-causal treatment of time as an extra spatial dimension. But there are indeed such researches like PINNsFormer[1] and PINNMamba[2] which consider and model the PDE system as a (pseudo-)sequence. Although this doesn't hurt the novelty and significance of this paper, it should be discussed.\n\n[1] Zhao Z, Ding X, Prakash B A. Pinnsformer: A transformer-based framework for physics-informed neural networks. ICLR 2024.\n\n[2] Xu C, Liu D, Hu Y, et al. Sub-sequential physics-informed learning with state space model. ICML 2025.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "eUja4HXrsK", "reviewer": "Reviewer_iLK9", "summary": "The paper proposes Frozen-PINN, a gradient-free Physics-Informed Neural Network that replaces iterative optimization with a space–time separation strategy. Spatial features are randomly sampled and fixed, while time-dependent coefficients are solved using ODE solvers, ensuring temporal causality by design. This approach greatly simplifies training and removes dependence on backpropagation. Across multiple PDE benchmarks, Frozen-PINN achieves much faster convergence and higher accuracy than existing PINNs and even matches mesh-based solvers, establishing a fast, causal, and interpretable alternative to gradient-descent-based PINN training.", "questions": ">Q1. See W1.\n\n\n>Q2. See W2.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": ">S1. \n\nThis paper presents a methodology that fundamentally differs from conventional PINN training approaches. Through a space-time separation framework, the proposed method circumvents the multi-objective optimization problem inherent in PINNs and naturally incorporates time causality, thereby avoiding several fundamental challenges that have long limited PINN performance.\n\n>S2. \n\nThe paper is well-organized with clear motivation. The proposed approach offers a reasonable alternative framework with solid justification for why it sidesteps the challenges faced by traditional PINNs.\n\n>S3.\n\nFrozen PINN demonstrates substantially faster training times and higher accuracy compared to the conventional PINN approach. The experiments are conducted on benchmark problems known to be challenging for PINNs, including cases with extreme advection speeds, shocks, and high-dimensionality. Across these benchmark, Frozen PINN consistently exhibits superior performance.", "weaknesses": "> W1.\n\nFrozen-PINN satisfies the initial condition by solving a least squares problem. While this decouples the initial condition from the PDE and boundary conditions, thereby simplifying the optimization, the least squares solution may not always be perfect. An analysis of the residual initial condition loss and its impact on model performance would be valuable. More broadly, investigating the sources of different loss components and their effects on accuracy would enhance our understanding of Frozen-PINN's behavior.\n\n>W2.\n\nDespite the various challenges in PINN training, increasing the number of collocation points (combined with appropriate sampling strategies and architectures) generally improves performance when sufficient training is applied. In my understanding, Frozen-PINN could also reduce PDE residual loss by investing more computation to solve the ODEs more accurately. However, since the initial condition relies on a predetermined least squares solution, it is unclear whether additional computation can improve the initial condition fitting. Is there a mechanism in Frozen-PINN to leverage increased computational resources for better initial condition satisfaction?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "HA6LJdqUf5", "reviewer": "Reviewer_Ukoa", "summary": "This paper presents Frozen-PINNs, a gradient-free framework for solving time-dependent partial differential equations. The method separates spatial and temporal components of the solution: spatial bases are randomly sampled and fixed, while the temporal coefficients evolve through a least-squares initialization and an adaptive ODE solver, ensuring temporal causality. By avoiding gradient descent, the approach simplifies optimization and achieves remarkable computational speedups and accuracy across various PDE benchmarks, including advection, diffusion, Burgers’, and high-dimensional heat equations.", "questions": "1. If I understand correctly, the authors use neural bases to represent the spatial domain and then integrate the resulting ODE system using a numerical ODE solver. If that is the case, why not simply apply a conventional ODE solver directly? The authors may argue that their approach improves upon standard PINNs, which I acknowledge. However, this direction seems to diverge from the central objective of PINN research in machine learning—to understand how well neural networks can approximate physical dynamics without relying on traditional numerical solvers. Such investigations are valuable because they provide insights transferable to other frameworks such as physics-informed neural ODEs, physics-informed neural operators, and physics-informed diffusion models. How do the authors justify this shift in focus? Do Frozen-PINNs have the potential to achieve higher accuracy than classical ODE solvers? What additional benefits or conceptual insights beyond accuracy does this approach provide?\n\n2. What numerical precision is used during training—float32 or float64? Please clarify, as this can significantly affect both runtime and accuracy.\n\n3. Is the method scalable in practice? It seems to rely on full-batch computation. How many collocation points are used for the 100-dimensional diffusion equation, and how does the cost scale with dimensionality?\n\n\n4. How were the SVD truncation threshold and L2 regularization chosen? Was any ablation study conducted to assess their sensitivity and impact on performance?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. The paper is clearly written and well organized, with informative figures and tables that effectively convey key ideas.\n2. The method delivers strong empirical performance, achieving large speedups and accuracy gains across diverse PDE benchmarks.\n3. The experiments are thorough and well documented, with detailed ablations and transparent reporting that make reproduction straightforward.", "weaknesses": "1. The proposed approach essentially reformulates the PDE into an ODE and integrates it with a standard solver, while the neural network merely provides a fixed spatial basis. As a result, the method’s “learning” component is limited, and the main computational work is performed by the numerical ODE solver. This makes the conceptual novelty questionable and leaves unclear why such a neural formulation is preferable to directly applying classical ODE or PDE solvers.\n\n2. The evaluation focuses mainly on relatively simple or smooth problems. To demonstrate robustness and generality, the method should be tested on more strongly nonlinear or coupled systems—such as the Navier–Stokes equations (e.g., lid-driven cavity flow)—where the dynamics are more complex and multi-scale behavior is prominent.\n\n3. From the hyperparameter tables (Tables 4 and 5, etc.), the proposed model appears considerably larger than the PINN and causal PINN baselines. Moreover, the baseline settings deviate from common practice—the PINN models are quite small, trained for fewer than 10,000 iterations, and in some cases use an unusually large initial learning rate of 0.1 with Adam. While this setup shortens their reported training time, it likely underrepresents their attainable accuracy, making the comparison of efficiency and performance less conclusive.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Oral)