$\boldsymbol{\partial^\infty}$-Grid: A Neural Differential Equation Solver with Differentiable Feature Grids
Authors: Navami Kairanda
Shanthika Naik
Marc Habermann
Avinash Sharma
Christian Theobalt
Vladislav Golyanik
Venue: ICLR 2026 Poster (2026)
Keywords: Differentiable Equations; Neural Field and Representations; Feature Grid; RBF Interpolation
Primary area: neurosymbolic & hybrid AI systems (physics-informed, logic & formal reasoning, etc.)
TL;DR: We propose the first feature-grid method for representing signals and solving PDEs, with differentiable interpolation and multiscale design. This enables accurate solutions unlike prior feature grids while outperforming neural solvers in efficiency.

Abstract: We present a novel differentiable grid-based representation for efficiently solving differential equations (DEs). Widely used architectures for neural solvers, such as sinusoidal neural networks, are coordinate-based MLPs that are, both, computationally intensive and slow to train. Although grid-based alternatives for implicit representations (e.g., Instant-NGP and K-Planes) train faster by exploiting signal structure, their reliance on linear interpolation restricts their ability to compute higher-order derivatives, rendering them unsuitable for solving DEs. In contrast, our approach overcomes these limitations by combining the efficiency of feature grids with radial basis function interpolation, which is infinitely often differentiable. To effectively capture high-frequency solutions and enable stable and faster computation of global gradients, we introduce a multi-resolution decomposition with co-located grids. Our proposed representation, $\boldsymbol{\partial^\infty}$-Grid, is trained implicitly using the differential equations as loss functions, enabling accurate modeling of physical fields. We validate $\boldsymbol{\partial^\infty}$-Grid on a variety of tasks, including Poisson equation for image reconstruction, the Helmholtz equation for wave fields, and the Kirchhoff-Love boundary value problem for cloth simulation. Our results demonstrate a 5–20× speed-up over coordinate-based MLP-based methods, solving differential equations in seconds or minutes while maintaining comparable accuracy and compactness.

Published as a conference paper at ICLR 2026
GGRID A NEURAL DIFFERENTIAL EQUATION SOLVER WITH DIFFERENTIABLE FEATURE GRIDS
Navami Kairanda', Shanthika Naik?, Mare Habermann', Avinash Sharma?,
Christian Theobalt', Vladislav Golyanik'
Max Planck Institute for Informatics, Saarland Informatics Campus
Indian Institute Technology Jodhpur
ABSTRACT
We present a novel differentiable grid-based representation for efficiently solv- ing differential equations (DEs). Widely used architectures for neural solvers, such as sinusoidal neural networks, are coordinate-ban based MLPs that are both computationally intensive and slow train. Although grid-based alternatives for implicit representations (e.g., Instant-NGP and K-Planes) train faster by ex- ploiting signal structure, their reliance on linear interpolation restricts their abil- ity to compute higher-order derivatives, rendering them unsuitable for solving DEs. Our approach overcomes these limitations by combining the efficiency of feature grids with radial basis function interpolation, which is infinitely dif- ferentiable. To effectively capture high-frequency solutions and enable stable and faster computation of global gradients, we introduce a multi-resolution de- composition with co-located grids. Our proposed representation, - is trained implicitly using the differential equations accurate modelling of physical fields. We validate -00 -rid on a variety tasks, including the Poisson equation for image reconstruction, the Helmholtz equation for wave fields, and the Kirchhoff--oow boundary value problem for cloth simulation. Our results demonstrate 5-20x speed-up over coordinate- based MLP-based methods, solving differential equations in seconds or minutes while maintaining comparable accuracy and compactness. See our project page: https::/44qgt http://4dvv...ii.i. inf mpg. de/DInffGrid//
loss functions, enabling
INTRODUCTION
Finding representations and parametrisations for fields, which both accurate and fast to evaluate, challenging and has many applications in many machine learning domains (vision, graphics, sim- ulation of physical systems, solving PDEs, engineering) (Xie et al. 2022). We desire representations that model global structure and capture high- figh-frequency as well as local details. Not only do we want represent known signals well, but also recover unknown fields by modelling physical systems governed by differentiable equations. At their core, many of these problems involve mapping spatial or spatio-temporal coordinates x E IRd to signals or fields u(x). E Rm Such systems governed by differential equations can be represented as implicit functions F of the form:
0, where g represents known vectors or functions controlling the system.
F (xx u, Vxu, Vu,...::
Replacing the traditional numerical solvers for equation I with neural networks has gained signifi- cant traction in recent years, particularly in the context of
informeedd machine learning (Kar-) niadakis et al. [2021] Hao al., 2022). The optimisation of network weights can emulate the process of solving physical equations; we refer to such networks as neural solvers for brevity. Current neural solvers are predominantly based on coordinate-basse multi-layer perceptrons (MLPs) or extensions thereof for efficiency, accuracy, etc.; while effective, such solvers exhibit several limitations. For instance, they tend to favour www-frqueecy fitting of signal u, requiring techniques like positional encoding /Vaswani et al. 2017) to handle figg-frreeency signals. Moreover, not all networks suitable for solving DEs, as they are not differentiable beyond the first order with respect to the input coordinates (e.g., those with ReLU activation lead v 0). Siren (Sitzmann al. 2020).

---

Published as a conference paper at ICLR 2026
Helmholtz Equation
Uniform wave velocity
Point wave source
90 seconds
minutes
4 minutes
Initial cloth state
Boundary
Cloth Simulation
Material; Gravity
Inputs
30 seconds
minutes Our solution field (w/ training time)
minute
Reference field
Figure 1: Our proposed method, -Grid, is feature grid-hased representation capable of accurately modelling signals and their higher-order derivatives. We demonstrate its effectiveness in solving complex differential equations, including the Helmholtz equation (top) and neural cloth simulation (bottom) (Kairanda)
et al. 2024) (reference adapted from (Clyde et al. 2017)) achieving both high speed and accuracy.
of these issues by being infinitely differen- tiable and capable of fitting fieghfrre signals, making it preferred choice as a neural field for physical problems in equation []] However, Siren and other neural solvers do not account for the spatial structure of signals and their derivatives in their architecture. Consequently, they are com- putationally expensive, requiring long training times, sometimes up several hours. On the other hand, to overcome the slow training and inference of neural fields (Mildenhall et al. (2020), explicit grid-based or hybrid architectures (Sara Fridovich Keil and Alex Yu at. 2022] Müller al. 2022 Cao & Johnson, (223)) that leverage the spatial locality have been proposed for scene representation. While they offer significantly faster training, they rely on d-linear interpolation, which not differ- entiable the grid boundaries (Vxl) 0) and not infinitely often differentiable in the grid interiors (Vzu. ... = 0). This severely limits their applicability as neural representations for solving DEs.
sinusoidal activation- based architecture, addresses
To address the computational overhead of neural solvers and the failure of grid-based representa- tions (Fridovich- Keil et al., 2023 Müller et al. 2022) in solving diffcrential equations F(.) = we propose novel feature grid representation capable of accurately modeling fields and their deriva- rives uxu, V u,.... We combine the efficiency feature grids with a differentiable interpolation technique supporting higher-order derivatives for effectively solving DEs. Our method leverages radial basis function (RBF) interpolation, which smooth and infinitely differentiable, where the gradients can be computed analytically using automatic differentiation. To represent high-frequency details as well as the low-frequency global structure, we propose a multi-resolution co-located fea- ture grid representation. Our approach offers several advantages. First, the feature grid localises learning, significantly speeding up training compared to global neural solvers. Second, our method provides continuous parameterisation despite using a grid, allowing queries at arbitrary points sim- ilar to neural fields. Consequently, it achieves accuracy comparable to (Sitzmann et al.) (2000) or better than (Raissi et al., 2019). i.e., MLP-based methods, while being significantly (up to 20x) faster. In summary, core technical contributions of our O method are as follows:
A feature grid-based representation for physical fields and a formulation recovering fields from differential equations (Sec. [.1]p
Differentiable RBF interpolation for accurately modelling higher-order derivatives (Sec. 3.2
Multi-resolution grids for faster global gradient flow (Sec. 3.3).
We validate Grid on diverse physical problems; see Fig. where we learn the solutions for
Image with Poisson Equation: A mapping from 2D coordinates to colours for high-resolution images using only gradient or Laplacian supervision;
Helmholtz Equation: Mapping spatial domain complex steady-state wavefields;
2

---

Published as a conference paper at ICLR 2026
Cloth Simulation: Mapping 2D material space to 3D deformations for cloth quasi-staticn mod-
elled as thin-shell boundary- problem:
Eikonal Equation: A signed distance function learned directly from oriented point cloud;
Advection and Heat Equations: Spatio-temporal transport and diffusion of physical quantities.
We prepared the source code for release, which act drop-in model for both representing fields and solving fields from governing equations.
RELATED WORK
Traditionally, differential equations are solved with numerical methods (Butcher 2016). Classi- cal solvers discretise governing equations but struggle with data incorporation, meshing, and high- dimensional costs, motivating neural representations as a promising alternative.
Neural Fields and DE Solvers. Neural fields continuous and non-linear mapping functions parametrised by trainable neural network weights. They have been widely used in representing
radiance fields (Mildenhall et al., 2020), signed distance and occupancy fields (Park et al. 2019
Mescheder et al. 2019; Yang et al.] 2021). scene reconstruction (Li et al., 2023] Wang et al, 2023 Guédon & Lepetit [2024] Yariv et al., (02T). and beyond. Physics-Informed Neural Net- works (PINNs) are neural fields that model solutions to differential equations. solved as optimisa- tion problems supervised by physical laws (Raissi et al. 2019 Hao al. 2022). A wide range of problems, such as solving fluid dynamics (Cai et al. 2, Eikonal equations (Smith et al. 2021), simulation (Bastek & Kochmann| 2023). and many more, have been achieved with neural fields. Early works that used coordinate-based MLPs struggled to model high-frequency details in signals and derivatives. Siren (Sitzmann et al.) 2020) introduces periodic activation functions for MLPs, facilitating gradient supervision with higher-order derivatives, enabling more complex applications (Chen et al. 2023a, Novello et al. 2023; Kairanda et al., 2024). However, Siren and similar global MLP-based methods have shortcomings: they ignore the spatial structure of physical fields, require back-propagating through the full MLP for every coordinate sample, thus leading to slower training times compared to our proposed grid- iad-based representation.
Feature Grids. In contrast to coordinate-basse MLPs, discrete or hybrid representations learn com- pact localised features, reducing computation and training time while improving quality. Various discrete representations, such as grids (Saaa Fridovich- Keil and Alex Yu et al. 2022] Martel et al.
2021. octrees (Takikawa al.) (2021) and feature planes (Fridovich-Ker) et al. [023]] Chan et al. Cao & Johnson 2023). in combination with MLPs, have been introduced for faster learning. Müller et al. (2022) employ a multi-resolltion structure with a hash encoding to provide faster lookup, significantly speeding up computation. Unfortunately. the linear com- ponents (linear interpolation, ReLU) of these feature learning methods do not support higher-order derivative supervision required for solving DEs. While some methods (Li et al. 2023 Huang & Alkhalifah 2024 Wang et al. 2024) support higher-order derivative supervision using finite dif- ference methods, they are known to be approximations and are susceptible to discretisation errors. Chetan et al. (2025) recently introduced a post-hoc differential operator for pre-trained hybrid fields, such Instant-NGP. Due to the reliance on pre-trained representation, they fail solve differen- tial equations. On the other hand, our proposed method combines the advantages of localised grid feature learning with high-order differentiability enable faster training and convergence.
2022 Chen et al. [022]] Cao
Radial Basis Functions. (Mai-Duy & Tran-Cong, 2003) were among the first to employ RBFs as activation functions, with later works such as (Ramasinghe & Lucey| 2022) and GARF (Chng al.)
2022) further demonstrating the utility of Gaussian activations. Recently methods inspired by 3D
Gaussian Splatting (Kerbl al. 2023), including Lagrangian Hashing (Govindarajan et al.. 2025)
and SC-GS (Huang et al. 2024). leverage Gaussian interpolation for feature aggregation. though their focus is on radiance fields rather than solving DEs. The closest to ours is NeuRBF (Chen et al. 2023b), which introduced feature grids with RBF interpolation. However, its objectives and therefore the design differ fundamentally: NeuRBF proposes neural field representations that adapt NeRFs), whereas -c0-Grid -Grid solves PDEs from derivative-only supervision without observing the target field; see App. B for details. Their Gaussian interpolation remains discontinuous across grid cells due to the single-ring neigh- bourhood, and the method requires signal-deennde RBF initialisation with ground-truth signals (or proxy hybrids) in the loss formulation, leading to failure of BFuRBF for our problem setting.
to known target signals (images. SDFs, or Instant-NGP-distiie
3

---

Published as a conference paper at ICLR 2026
RBF Interpolated
Features
V. Fxx,, Van, P
Fo
f(x:
Feature
Decoder
p[||x
x.|
6 (x))
5
u
Van. V.u. U.
F
Input Domain
(f(x))
22
Es
(x)
Multi-resolution Feature Grid
Boundary
Solution Field
Differential Quantities
Figure 2: Illustration of -0 -Grid. Our grid-based representation accurately represents fields 92 R and efficiently solves differential equations F. We first sample query coordinates x from the input domain Ω (here, 2D). We apply infinitely differentiable RBF interpolation yx(x) to extract smooth, query-dependent features from the multi-scale learnable feature grids {F [Fs)ecoo....ss These interpolated features {f, (x)), at different scales
concatenated as f(x) and then passed through decoder (f((x): Q) produce the final signal u(x) (e.g., image cloth deformation). The converged neural field (solution) is obtained by optimising the governing equation J as the loss.
3 METHOD
Our goal to represent and solve for the field subject to physical constraints F, expressed as dif- ferentiable equations, as shown in Eq. (I. The solution must be accurate, and the optimisation needs to be fast. This section presents our feature grid-based representation, which is computationally ef- ficient (Sec. 3.1. 3.11 While grid-based representations have been explored in prior work (Müller al.) 2022 Fridovich- Keil et al. 2023), these approaches often lack differentiability with respect to inputs making them unsuitable for solving differential equations accurately. In contrast, our method sup- ports higher-order derivatives by combining feature grids with the smoothness and diffcrentiability introduce multi-resolution technique effec- tively capture high-frequency signals which ensures stable optimisation with fast global gradient propagation (Sec. 3.3 An overview of our approach provided in Fig.
of RBF interpolation (Sec. 3.2p. Additionally,
3.1 FEATURE GRID FORMULATION
We aim to parameterise the field uxx): Ω + Rm where Ω C IRd is a low-dimensional (d = 2,3) spatio-temporal input domain. Directly processing the query coordinates x with a coordinate-based MLP, as in neural fields (Xie et al.] 2022] Hao et al.] 2022), requires updating all weights during backpropagation, which is computationally expensive and bypasses the opportunity to exploit the local structure of the spatial domain. To address these limitations, we propose novel encoding- decoding framework. The signal is computed as a composition of two functions:
(2) where f a feature encoder parametrised by the feature grid F, and d is a small decoder. This for- mulation leverages the spatial structure in the feature encoder, requiring only small set of trainable parameters of grid F adjaceent to the query coordinate x x-to updated during optimisation.
u(x) = dff(x; F); θ);
Feature Encoder. We adopt grid-based representation, where the domain of interest is discre- tised into a regular grid. Each grid node stores learnable feature vector, which optimised during training. Unlike traditional coordinate-basse MLPs, this representation exploits spatial locality, abling efficient computation. Since feature grids are discrete, interpolation is required to query continuous points. Formally, let F (N+) denote the feature grid, where N is the resolution (assumed identical along all spatial dimensions for simplicity), and F is the feature dimension at each grid node. Given query coordinate E Rd, the corresponding feature vector f(x) E R is obtained via interpolation over the grid. The interpolated feature is computed as:
f(x)
Σ (EA(x)
w(x, x; F()))
(3)
4

---

Published as a conference paper at ICLR 2026
where (xx) is the set of neighbouring grid nodes enclosing x; F(x;) and w(x,x,) x,) are the learnable feature vector and the scalar interpolation weight corresponding to node x,, respectively.
Feature Decoder. To map the interpolated feature f(x) to the target signal u(x), employ feature decoder d(;; Q) If + IRm with weights θ The decoder can be implemented as either a simple linear layer or a small MLP (with smooth activation e,g, tanh), depending on the complexity of the task. This additional decoding step allows the model to map the features to the target flexibly.
Differential Equation Solver. We propose to compute gradients, Jacobians, and Laplacians of the decoded interpolated features with respect to query samples x using automatic differentiation (autograd). These derivatives are essential for solving differential equations. Recalling our goal of optimising for in problems of the form Eq. ( we define the loss function as:
((F / F()
(,,)))) Σ (x,, u(x;), Vx, u(x)), V2 u(xi)..... :x(x)).
(4)
where u(F, +) parameterised in Eq. 2 The loss function represents the PDE residual over the apatio-tempora domain, optionally including boundary constraints or data terms. Our framework accommodates various formulations: strong form (e.g., Poisson), variational form with boundary constraints (e.g.. cloth simulation), and data-driven terms (e.g., SDF). To guarantee uniqueness and stability, we typically impose Dirichlet boundary or initial conditions as hard constraints, which known improve convergence (Kairanda et al., 2024 Lu et al.) [02I]. This is achieved by modulating the predicted field in Eq. (2) with a bwwndaoy-wwweeo....ooaaaa aware function:
u(x) d(f(x: F); B(xx + h(x) (1-B(x)),
(5)
where h denotes the known boundary function and ((x) = 1 - e- xxxl= /a is a distance-based weighting that enforces u(x) h(x) for all x ön. This formulation is equally applicable to initial conditions, where the boundary is defined by t= 0.
During training, optimise both the feature grid F and the decoder θ using gradient- based op- timisation. In practice, the integral in Eq. 4 is estimated by sampling points from the domain Ω. A dataset D {(x;; gxx))); constructed by stratified random sampling of coordinates x; E Ω for PDE evaluation, optionally including data coordinates for loss terms and their associated known quantities g(x)). The dataset 2 is sampled dynamically during training, progressively refining the approximation of L as the number of samples increases.
3.2 RADIAL BASIS FUNCTION INTERPOLATION
In Eq. (3), the interpolation weights w(x, x.) play a crucial role in determining the contribution of neighbouring grid points to the interpolated features. These weights can be computed using various methods. Prior works on scene representations often rely on d-linear interpolation (Müller et al. 2022; Fridovich-Kad Keil et al. 2023), where the weights determined based on the relative distances between the query point x and the nearest grid points along each dimension. Formally, let (I1,, 22,,, ,....,,,a, aaa) represent the query point, and x; (Ia ... . Fea) denote grid node. The weight for dllinear interpolation is given by:
w(x, x;)
(1 A=
II
d
Fh 1ick
(6)
where a is the grid cell size. While linear interpolation is computationally efficient, it is only C*( (continuous but not differentiable) at grid nodes, and (²R( in the grid interior, when differentiating with respect to query coordinates. The lack of support for higher-order derivatives prohibits the computation of differential quantities, such as Jacobians and Laplacians, which are essential for solving DEs as described in Eq. (4). An alternative to linear interpolation is to employ higher-order basis functions, such as quadratic, cubic, or RBFs. We employ RBF interpolation (Wright) 2003). which is smooth across neighbourhoods and infinitely differentiable (C*((R) (R)). The interpolation weights are defined as
p[[|x x+ w w(x, EjeNto p([[x xIll)
5

---

Published as a conference paper at ICLR 2026
where (x) denotes the neighbourhood of x, and the denominator normalises the weights. We use Gaussian RBF kernels, ((()) e e-(rr), where & controls the kernel width. RBFs provide differentiable representation that supports higher-order derivatives when required. For instance, fourth-order PDEs such as the Euler- BeleolBernoulli beam or Kirchhoff Love plate equa- tions (Guo et al. 2022)) demand differentiability beyond what cubic bases can provide, particularly in strong-form settings. Moreover, RBFs offer tunable support via the shape parameter & (which we set hyperparameter), enabling flexible trade-off between locality and accuracy that degreee polynomials cannot achieve.
Efficient RBF Interpolation. In RBF interpolation, neighbour- ing grid points N (x) are weighted according to their distance from the query point As shown in Eqs. (3) and ( Gaussian RBFs (p(r) = e (er)²) are globally supported, requiring the neighbour- hood to be the full grid-Icading to high memory and computa- tional costs. Consequently, standard RBF interpolation becomes impractical for high-resolution grids or large numbers of samples.
Nx x
(a) RBF interpolation ID
F
(x
To match the efficiency of linear interpolation with the smooth- ness of RBFs, we sample the domain in structured way and pre- compute neighbourhood indices. At each training step, we stratify samples over domain St ensuring they are both uniformly dis- tributed and random for more accurate approximation of Eq. (4).
In conjunction with stratified sampling, we propose precompu- tation of the neighbourhood ((x), which can be reused across training iterations. For given query point x, the effective neigh- bourhood (x). where the RBF weights are non-zero, is deter- mined by the shape parameter of the RBF kernel. For example, if = 2, only the 1-ring neighbourhood (1(x) needs to be con- sidered, whereas for e required. Since the sample point may lie anywhere within grid cell, including near boundaries, we include an additional ring in practice. For instance, N2(xx) is used for 2; see Fig. for an illustration in d = 1. This approach ensures smooth interpolation with feature overlap across query samples, ensuring accurate learning of the feature grid.
(b) RBF support region in 2D
Figure 3: Illustration of effective neighbourhoods (x in d 1,2, determined by the shape pa- rameter E. An additional ring ac- counts for stratified sampling.
the 2-ring neighbourhood (x)x)
3.3 MULTI-RESOLUTION FEATURE GRID
We observe that our feature encoder with a single-resolutio grid captures local gradients effectively and requires many time steps to propagate gradient information across the entire domain. This limitation stems from the small interpolation neighbourhoods used in relatively large feature grids. Drawing inspiration from recent advances in scene representations (Müller al. 2022] Fridovich- Keil et al. 2023), propose multi-resolution approach to overcome this challenge. We define a set of multi-scale feature grids, each representing a different resolution. Formally, for given the base grid resolution Nmax and the number of scales S, the feature grids {F, [F=)sc(o.....-1))
learnable parameters: F, E R(N+1)*xF where N = N2ax/2" is the resolution at scale s. With this formulation, the interpolated features from all scales are concatenated to form a multi-scale feature vector f(x) = (fo (x), fi (x)..... .....fs- f.-fsxs)).... (x)) e RS- where f, (x) is computed using the RBF interpolation described in Eqs. (3) and (7).
The multi-resolution feature grid offers several advantages. First, it is computationally efficient, as the hierarchical structure reduces the cost global gradient propagation. Second, it is highly expressive, as the combination of coarse and fine features, along with their gradients, enables the stable and accurate recovery of complex signals without requiring adaptive RBF shapes.
EXPERIMENTS
We experimentally demonstrate the versatility of representation across diverse signal types, in- cluding the Poisson equation (Sec. (.1). the Helmholtz equation (Sec. (42),, and neural cloth sim- ulation (Sec. 4.3 Unlike prior grid-based methods such as NGS (Müller al. 2022

---

Published as a conference paper at ICLR 2026
Supervision with grdd
Supervision with Liupl
x()
Siperision
Au(
Supervisos
Siren
Ours
Siren
Ours
Training Time (mins)
Figure Image Reconstruction. Comparison of 8 -Grid with Siren (Sitzmann et al. 2020) for the recon- struction from gradient and Laplacian fields (Poisson equation). Results are visualised at training time marked by "+" and "v", and the supervision signal is highlighted with green border. Our method achieves higher PSNR with faster convergence for high-resolution images, while Siren struggles, especially for Laplacian.
Table 1: Comparison against Siren (Sitzmann et al. 2020) and K-planes (Fridovich Keil et al. (2023 for image reconstruction (Fig. 4 go--Grid achieves higher PSNR with substantially reduced training time than Siren.
K- K-planes struggle (for gradient supervision) or fail (for Laplacian) due to non-differentiable interpolation.
Model
K-Planes
Siren
Ours
Supervision
Grad.
Lapl.
Grad.
Lapl.
Grad.
Lapl.
PSNR
17.96
32.12
11.82
32.24
12.19
Training time
9.5mins
10mins
1hr56mins
25s
15mins
Parameters
8.4m
329 99k
1.32m
702.94k
700.81k
K-Planes (Fridovich-Ker al.) 2023), and NeuRBF (Chen 2023b our feature grid repre- sentation accurately solves for fields. We show that 000 -Grd achieves significant speed-ups over sinusoidal neural networks (Sitzmann et al., (2020) while maintaining comparable accuracy, and that it can recover high-frequenee solution fields where PINNs struggle (see Sec. [B) Apart from solving DES, we demonstrate that our representation supports direct signal fitting (Sec. We also include ablations the RBF kernel size and the multiscale design of the feature grid (Sec. 4.5) Additional experiments are provided in the appendices, i.e., SDF solutions to the Eikonal equation, the spatio- temporal heat equation, and ID/2D advection with Gaussian waves (with comparison to (Chen al. 2023ah). We further demonstrate highly challenging cases, such as Zalesak's s disk and Helmholtz problems with higher wavenumbers.
4.1 IMAGE RECONSTRUCTION POISSON EQUATION SOLVE
To demonstrate that our feature grid representation can accurately reconstruct a field and its deriva- tives, we first consider simple differentiable equation- Poisson equation. This equation rep- resents class of problems where the objective is to estimate an unknown signal g from discrete samples of its Laplacian Ag = V Vg. As a simpler variant, we additionally estimate the signal from its gradients Vg. which can be formulated as minimising one of the following losses (as instance of Eq. (4)):
Eegaad / [[Vxu(x) Vxg(x)] dx. Lippl = ||An(x) Ag(x))l dx.
(8)
We solve the Poisson equation for images by supervising the model solely with higher-order deriva- tives-either ground-trutt truth gradients, as in Landd or Laplacians Llapt- Our model successfully recon- structs the target image, shown in Figs. [ and 7 For colour images supervise gradients channel-wise and visualise the red-channel field for clarity. Global intensity variations remain due to the ill-posedness of the inverse problem, boundary conditions not explicitly specified.
To assess computational efficiency against coordinate- MLPs, we evaluate high-resolution im- age reconstruction with method and compare it Siren (Sitzmann al.] 2020). While the

---

Published as conference paper at ICLR 2026
Error (Real)
Real
U (x)
Error (Imag)
Imaginary
Ua(x)
K-Planes
Instant-NGP
Siren
Ours
Reference
Training Time (mins)
Figure 5: Helmholtz Equation. Comparison for solving the Helmholtz equation (sccond-order PDE) for single point source at the centre. Prior hybrid architectures, such as K-Planes (Fridovich- Keil et al. 2023 and Instant-NGP (Mäller al. 2022). break down their linear components equate for second-order derivatives. Our method matches reference with significantly faster training time compared to Siren, as shown in the & error plot on the right and e--rrror visualisations in Fig. VI
the full 512 x 512 resolution for a fairer comparison and clearer demonstration training speed. We adopt the smallest Siren ar- chitecture that converges reliably, as reducing capacity can speed up training but could also lead to underfitting. As shown in Fig. (right) our method captures both global structure and fine local detail more accurately- and converges significantly faster (20x). Quantitative comparisons for the Poisson equation are reported in Tab. for our method, Siren, and K K-Planes (Fridovich- Keil et al.) 2023). Notably, existing grid-based approaches such as K-Planes cannot handle higher-order super- vision: with linear interpolation, first derivatives are discontinuous at grid boundaries and second derivatives vanish everywhere, preventing recovery from Laplacian fields.
original Siren paper reports results 256 256 images,
4.2 HELMHOLTZ EQUATION
We further demonstrate the applicability of our method to wave propagation problems governed by complex boundary conditions. Specifically, we consider the Helmholtz equation, second-order partial differential equation commonly used to model steady-state wave behaviour and diffusion phenomena. The Helmholtz equation is given by:
H(m) uxx) -x(),, with H(m) (A + m(x)
where gxx) denotes known source function, u(x) the unknown complex-vaie valued wavefield, w is the wave number, and m(x) = e(x)a represents the squared slowness as a function of wave velocity c(x). To solve for the wavefield, tion. The decoder outputs two channels corresponding to the real and imaginary components of the wavefield. We set up the input domain as {x E R² ||x|| < 1]. and consider spatially uniform wave velocity with single point source, modelled as a Gaussian with variance o² 10-4. The model is trained using loss function derived directly from the Helmholtz formulation:
parameterise u(x) using our proposed 00 -Grd representa-
Hellmoltz 1 xxx) [[(m)) uxx) + s(x) dx,
(10)
where A(x) = l (a tunable hyperparameter) for locations where 9(x) * O, corresponding to the inhomogeneous part of the domain, and X(x) elsewhere. Unlike the image reconstruction task from derivative fields-whert data sampling is required- employ stratified sampling here, since Litelmholtz can be evaluated arbitrary points within the domain. We incorporate perfectly matched layers (PML), which modify the Helmholtz PDE to attenuate outgoing waves the boundary, ensuring uniqueness of the solution (see Sec. for implementation details). All methods, including ours, predict only the interior of the domain, the PML-modified Helmholtz formulation absorbs steady-state wavefields generated by the point source. The reference solution computed in closed form using Hankel functions and therefore extends across the entire domain.
Our representation achieves tigh-fidelity reconstructions the wavefield, as shown in Figs. and for 20, where compare against the closed-form solution. While both methods match the
8

---

Published as a conference paper at ICLR 2026
Table 2: Ablation study for varying RBF shapes and neighbourhood sizes for the high-resolution gradient- based image reconstruction in Fig.
Shape
Ring p
PSNR
SSIM 1
Training +
0.6
28.74
0.87
15 min
4
28.74
0.87
17 min
K-Planes
Instant-NGP
Ours
3
28.73
0.87
12 min
2
28.82
0.86
12 min
Figure 6: Kirchhoff-Love Simulation (Kairanda
1"
27.50
0.82
15 min
etal. 2024). -0- captures fine -innggg cloth wrinkles by solving highly nonlinear problem.
2"
1
17.72
0.74
6 min
reference closely, our approach converges faster, yielding a 4xspeed-up over Siren; see Fig. (right). Additional results for more challenging cases with w 30, are presented in Sec. and Tab. [I] We also compare with prior feature grid methods: (Fig. (left): They fail as effective solvers since they cannot capture the higher-order derivatives required by the Helmholtz equation.
4.3 CLOTH SIMULATION
NeuralClothSim (Kairanda et al., 2024) introduces quasistatic cloth simulation framework that implicit neural deformation field (NDF). Given
models surface deformations as a coordinate- simulation cruup-defined by the cloth's 's initial configuration, material parameters, boundary con- straints and applied external forces-the the method learns the equilibrium deformation field governed by the Kirchhoff-Love thin shell theory. The deformation represented as an NDF u(x) - IR: where Ω C R.2 denotes the cloth's s parametric domain. In the context of volumetric thin shells such as cloth, the Kirchhoff-Love model (Wempner & Talaslidis) (200)) provides a reduced kine- matic formulation, where the shell midsurface fully characterises the strain distribution across the thickness. Following standard thin shell assumptions, the deformation field u is used to compute the Green strain, which is decomposed into membrane strain and bending strain K, capturing in- plane stretching and out-of-plane curvature changes, respectively. The internal hyperelastic energy (ie,, k:; Φ) is then formulated as a functional of these geometric strains and the material properties 4. Under the action of external forces and boundary conditions b, the simulation seeks an equilib- rium configuration by minimising the total potential energy, combining internal energy and external work. This principle is implemented by defining the loss as (see Sec. for further details):
Lalath = / (ee(u(x)), <((());; +b] g(x) u(x) dx, subject to u(x) = b(x) on M.. (11)
Note that Kairanda al. (2024) employ Siren as the representation for NDF u. while we use our -00 gGidd. As a specific example, we consider an analytically defined square piece of cloth of mass density p held with two handles and subject to a gravity force g [0, 0, 9.8p). We incorpo- rate Dirichlet constraints directly into the loss following Eq. (5). and assume linear elastic material model for 4. See Fig. 6 results, where we recover realistic simulations, including folds and wrin- kles, by solving highly non-linear PDE where stretching and bending strains depend non-linearly on the deformation. Additional comparisons with K-Planes and Instant-NGP show their failure, as they are unable to model the membrane and bending strains, which require higher-order derivatives of the deformation field. As multiple equilibria are possible for the cloth deformation (thus, no single reference solution), exclude cloth simulation from the numerical accuracy evaluation.
4.4 SIGNAL REPRESENTATION
The following additional experiment underscores the versatility of our representation (although not being the central focus of our work). Beyond solving DEs, our method can also signals directly by instantiating Eq. (4) with supervision signal, i.e..
Laggnal
|u(x) - g(x)ll dx,
(12)
where g gxx denotes the known target signal such as a 2D image or an SDF. We showcase both cases in Fig. 7 (a,c), where method preserves high-frequency details on par with previous competitive

---

Published as a conference paper at ICLR 2026
Supervision with (n) g(x)
(b) Van(x)
SDF(x)
(d) Lelkonal
u(x)
u(x
Vegix
(D
Gradient Image
Reconstructed
Fitted Surface
Oriented Pointeloud
Reconstructed
Image
Fitted Image
Surface
Image
Surface
Figure 7: Results of signal representation. Our -00-Gria matches high-fregueen frequency details when fitting images and SDFs directly (a.c) and, unlike prior grid methods (Müller al. 2022 Chen et al. (023b) (Figs. | and
Tab. also succeeds with gradient and Eikonal supervision (b,d). Inputs are outlined in green.
gridd based baselines (Müller et al. 2022 Chen et al., 2023b). Our representation further supports indirect supervision: We recover images from gradient fields and SDFs from oriented point clouds by solving the Eikonal equation (Sec. C.O as illustrated in Fig. ().dd.. In these settings, the competing grid-based methods fail (Figs. and and Tab, Tab.[Tin the Appendix). Together, these results show that single differentiable grid accommodates both direct signal fitting and PDE-constrained
reconstruction without architectural changes.
4.5 ABLATIONS AND FURTHER COMPARISONS
Neighbourhood Ring Selection. The choice of RBF shape parameter E, and consequently the effective neighhourhood N, balances computational efficiency and accuracy. To study the trade-off between truncation, accuracy, and speed well potential instabilities we conducted ablations on image recovery from the higher-resolution 768 768 gradient image in Fig. [ First, we vary and adjust the neighbourhood size adaptively, that the RBF weights are non-zero (e.g., P 4 for e = 0.6). Smaller e values with larger rings produce smoother, more accurate reconstructions but at higher computational cost, while overly narrow kernels (e.g., & 2. p 1) introduce discontinuities resembling linear interpolation. Second, we fixed æ = and varied the ring size p. The 1 1-ring setting showed visible artefacts, whereas & with 2- or 3-ring neighbourhood offers stable training and smooth interpolation without artefacts, matching our intuition from Fig.
Multi-resolution Grid. We demonstrate the effectiveness of multi-resolution feature grids in Fig. m In the SDF example, the single-resolution grid restricts the Eikonal loss to local regions, leading to noticeable artefacts. In contrast, the multi-resolution grids enable gradient propagation across the entire domain while incurring minimal additional parameter overhead. Next, as illustrated by the image recovery example (Fig. Mb. this design leads to significantly faster convergence for PSNR.
Comparisons to NeuRBF (Chnn et al. 2023b). We modified NeuRBF's umage-fitting objective to use our gradient-based Poisson loss and retained their gradient-weighted RBF initialisation; the reconstruction failed to converge at PSNR/SSIM 10.85/0.17 10 with visible artefacts (cf. Fig. [b. Ap- plying the same pipeline to Eikonal SDF fitting from oriented point clouds also failed to recover a smooth surface (ef. Fig. V. highlighting the difficulty of training NeuRBF under purely differential supervision. Please refer to App. B for the detailed analysis.
5 CONCLUSION
The experiments demonstrate that 000 gGid enables computation of higher-order derivatives di- rectly from the interpolated features, thanks to smooth and infinitely differentiable RBF interpola- tion weights. We can reliably compute gradients, Jacobians, and Laplacians for solving PDEs. This leads to accurate results for Poisson and Helmholtz equations, allows recovery of high- -reeurryyyncy fields for challenging and non-linear PDEs, such as the Kirchhoff- boundary value problem for cloth simulation and the Eikonal equation. Moreover, we experimentally show faster convergence compared to coordinate-based MLPs. One current limitation is that the interpolation becomes ex- pensive for higher-dimensee dimensional inputs, similar to numerical solvers, and a potential solution for that in the future could be projection to diver-dieeniionl spaces. Moreover, while we observed conver- gence across all experiments with the proposed RBF neighborhood truncation, further investigation needed establish theoretical guarantees for C under truncation.
10

---

Published as a conference paper at ICLR 2026
REFERENCES
Jan-Hendrik Bastek and Dennis M Kochmann. Physics-informed neural networks for shell struc- tures. European Journal of Mechanics-ASooiis, 2023.
John C. Butcher. Numerical Methods for Ordinary Differential Equations. John Wiley & Sons, 2016.
Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Kariiadakis. Physics- informed neural networks (pinns) for fluid mechanics: review. Acta Mechanica Sinica, 2021.
Ang Cao and Justin Johnson. Hexplane: A fast representation for dynamic scenes. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2023.
Eric Chan, Connor Z. Lin, Matthew A. Chan, Koki Nagano, Boxiao Pan, Shalini De Mello, Orazio Gallo, Leonidas J. Guibas, Jonathan Tremblay, Sameh Khamis, Tero Karras, and Gordon Wetzstein. Efficient geometry-awan 3d generative adversarial networks. In Computer Vision and Pattern Recognition (CVPR), 2022.
Anpei Chen, Zexiang Xu, Andreas Geiger, Jingyi Yu, and Hao Su. Tensorf: Tensorial radiance fields. In European Conference on Computer Vision (ECCV), 2022.
Honglin Chen, Rundi Wu, Eitan Grinspun, Changxi Zheng, and Peter Yichen Chen. Implicit neural spatial representations for time-dependent pdes. In International Conference on Machine Learn- ing (ICML), 2023a.
Zhang Chen, Zhong Li, Liangchen Song, Lele Chen, Jingyi Yu, Junsong Yuan, and Yi Xu. Neurbf: A neural fields representation with adaptive radial basis functions. In International Conference on Computer Vision (ICCV), 2023b.
Zhongying Chen, Dongsheng Cheng, Wei Feng, and Tingting Wu. An optimal 9-point finite differ- ence scheme for the helmholtz equation with pml. International Journal of Numerical Analysis & .Modeling, 2013.
Aditya Chetan, Guandao Yang, Zichen Wang, Steve Marschner, and Bharath Hariharan. Accu- rate differential operators for hybrid neural fields. In Computer Vision and Pattern Recognition (CVPR), 2025.
Shin-Fand Chng. Sameera Ramasinghe, Jamie Sherrah, and Simon Lucey. Gaussian activated neural radiance fields for high fidelity reconstruction and pose estimation. In European Conference on Computer Vision (ECCV), 2022.
David Clyde, Joseph Teran, and Rasmus Tamstorf. Modeling and data-driven parameter estimation for woven fabrics. In Proc. ACM SIGGRAPH/ Eurographics Symposium on Computer Animation (SCA), 2017.
Todd F Dupont and Yingjie Liu. Back and forth error compensation and correction methods for removing errors induced by uneven gradients of the level set function. Journal of Computational Physics, 2003.
Ronald Fedkiw, Jos Stam, and Henrik Wann Jensen. Visual simulation of smoke. In Proceedings of the 28th annual conference on Computer graphics and interactive techniques, 2001.
Sara Fridovich-Keie Giacomo Meanti, Frederik Rahbak Warburg, Benjamin Recht, and Angjoo Kanazawa. K-planes: Explicit radiance fields in space, time, and appearance. In Computer Vision and Pattern Recognition (CVPR), 2023.
Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural networks. In International Conference Artificial Intelligence and Statistics. JMLR Workshop and Conference Proceedings, 2010.
Shrisudhan Govindarajan, Zeno Sambugaro, Akhmedkhan Shabanov, Towaki Takikawa, Daniel Re- bain, Weiwei Sun, Nicola Conci, Kwang Moo Yi, and Andrea Tagliasacchi. Lagrangian hash- ing for compressed neural field representations. In European Conference on Computer Vision (ECCV), 2025.
11

---

Published as a conference paper at ICLR 2026
Antoine Guédon and Vincent Lepetit. Sugar: Surface-aligned gaussian splatting for efficient 3d mesh reconstruction and high-quality mesh rendering. Computer Vision and Pattern Recognition (CVPR), 2024.
Hongwei Guo, Xiaoying Zhuang, and Timon Rabczuk. A deep collocation method for the bending analysis of kirchhoff plate. arXiv preprint wrXiv:2102.02611 2021.
Zhongkai Hao, Songming Liu, Yichi Zhang, Chengyang Ying, Yao Feng, Hang Su, and Jun Zhu. Physics-informed machine learning: A survey on problems, methods and applications. arXiv preprint arXiv:22/1 XXi:222/1.00064, 2022.
Xinquan Huang and Tariq Alkhalifah. Efficient physics-informed neural networks using hash en- coding. Journal of Computational Physics, 2024.
Yi-Hua Huang, Yang- Tian Sun, Ziyi Yang, Xiaoyang Lyu, Yan-Pei Cao, and Xiaojuan Qi. Sc-gs: Sparse-controlles gaussian splatting for editable dynamic scenes. In Computer Vision and Pattern Recognition (CVPR), 2024.
Navami Kairanda, Mare Habermann, Christian Theobalt, and Vladislav Golyanik. Neuralelothsim: Neural deformation fields meet the thin shell theory. In Advances Neural Information Process- ing Systems (NeurIPS), 2024.
George Em Karniadakis, Ioannis Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 2021.
Bernhard Kerbl, Georgios Kopanas, Thomas Leimkühler, and George Drettakis. 3d gaussian splat- ting for real-time radiance field rendering. ACM Transactions on Graphics, 2023.
Zhaoshuo Li, Thomas Müller, Alex Evans, Russell H Taylor, Mathias Unberath, Ming- Yu Liu, and Shen-Hsuan Lin. Neuralangelo: High-fidelity neural surface reconstruction. In IEEE Conference Computer Vision and Pattern Recognition (CVPR), 2023.
Lu Lu, Raphael Pestourie, Wenjie Yao, Zhicheng Wang, Francesc Verdugo, and Steven G Johnson. Physics-informed neural networks with hard constraints for inverse design. SIAM Journal on Scientific Computing, 2021.
Nam Mai-Duy and Thanh Tran-Cong. Approximation of function and its derivatives using radial basis function networks. Applied Mathematical Modelling, 2003.
Julien N. P. Martel, David B. Lindell, Connor Z. Lin, Eric R. Chan, Marco Monteiro, and Gordon Wetzstein. Acom: adaptive coordinate networks for neural scene representation. ACM Transac- tions on Graphics, 2021.
Lars Mescheder, Michael Occhsle, Michael Niemeyer, Sebastian Nowozin, and Andreas Geiger. Occupancy networks: Learning 3d reconstruction in function space. In Computer Vision and Pattern Recognition (CVPR), 2019.
Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng. Nerf: Representing scenes as neural radiance fields for view synthesis. In European Conference on Computer Vision (ECCV), 2020.
Thomas Müller, Alex Evans, Christoph Schied, and Alexander Keller. Instant neural graphics prim- itives with multiresolution hash encoding. ACM Transactions on Graphics, 2022.
Tiago Novello, Vinícius da Silva, Guilherme Schardong, Luiz Schirmer, Hélio Lopes, and Luiz Velho. Neural implicit surface evolution. In International Conference on Computer Vision (ICCV), 2023.
Jeong Joon Park, Peter Florence, Julian Straub, Richard Newcombe, and Steven Lovegrove. Deepsdf: Learning continuous signed distance functions for shape representation. In Computer Vision and Pattern Recognition (CVPR), 2019.
12

---

Published as a conference paper at ICLR 2026
F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Pretten- hofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Couraapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 2011.
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, 2019.
Sameera Ramasinghe and Simon Lucey. Beyond periodicity: towards unifying framework for activations coordinate-mlps. In European Conference on Computer Vision (ECCV), 2022.
Sara Fridovich-Keil and Alex Yu, Matthew Tancik, Qinhong Chen, Benjamin Recht, and Angjoo Kanazawa. Plenoxels: Radiance fields without neural networks. In Computer Vision and Pattern Recognition (CVPR), 2022.
Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein. Im- plicit neural representations with periodic activation functions. In Advances Neural Information Processing Systems (NeurIPS), 2020.
Jonathan D. Smith, Kamyar Azizzadenesheli, and Zachary E. Ross. Eikonet: Solving the eikonal equation with deep neural networks. IEEE Transactions on Geoscience and Remote Sensing, 2021.
Towaki Takikawa, Joey Litalien, Kangxue Yin, Karsten Kreis, Charles Loop, Derek Nowrouzezahrai, Alec Jacobson, Morgan McGuire, and Sanja Fidler. Neural geometric level of detail: Real-time rendering with implicit 3D shapes. In Computer Vision and Pattern Recognition (CVPR), 2021.
Jiaxiang Tang. Torch-ngp:
pytorch implementation of instant-ngp. 2022.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Advances in Neural Infor- mation Processing Systems (NeurIPS), 2017.
Haoxiang Wang, Tao Yu, Tianwei Yang, Hui Qiao, and Qionghai Dai. Neural physical simulation with multi-resolution hash grid encoding. AAAI Conference on Artificial Intelligence, 2024.
Yiming Wang, Qin Han, Marc Habermann, Kostas Danilidis,, Christian Theobalt, and Lingjie Liu. Neus2: Fast learning of neural implicit surfaces for multi-view reconstruction. In International Conference on Computer Vision (ICCV), 2023.
Gerald Wempner and Demosthenes Talaslidis. Mechanics of solids and shells. CRC, Boca Raton, 2003.
Grady Barrett Wright. Radial basis function interpolation: numerical and analytical developments. University of Colorado at Boulder, 2003.
Yiheng Xie, Towaki Takikawa, Shunsuke Saito, Or Litany, Shiqin Yan, Numair Khan, Federico Tombari, James Tompkin, Vincent Sitzmann, and Srinath Sridhar. Neural fields in visual comput- ing and beyond. In Computer Graphics Forum (Earographics State the Art Reports), 2022.
Guandao Yang, Serge Belongie, Bharath Hariharan, and Vladlen Koltun. Geometry processing with neural fields. In Advances in Neural Information Processing Systems (NeurIPS), 2021.
Lior Yariy, Jiatao Gu, Yoni Kasten, and Yaron Lipman. Volume rendering of neural implicit surfaces. In Advances Neural Information Processing Systems (NeurIPS), 2021.
13

---

Published as a conference paper at ICLR 2026
APPENDICES
Reproducibility Details
15
Further Comparisons
16
B. Comparison with NeuRBF (Chen et al. (023b)
16
B.2 Comparison with PINNs
17
C Eikonal Equation
18
Implementation Details
18
Comparison with Grid-based Methods
18
D Poisson Equation
19
D.I Implementation Details
19
Helmholtz Equation
19
Perfectly Matched Layers Formulation
20
Implementation Details
20
Higher Wavenumber
Cloth Simulation
21
Implementation Details
Heat Equation
22
H Advection Equation
23
ID Advection
23
H.2 2D Advection
23
H.3 Zalesak's Disk Problem
23
imiitaiion and Future Work
23
14

---

Published as a conference paper at ICLR 2026
g(x)
Vx9(x)
Input for Supervision
Caignal ! ||n(x) o(x)ll dx
Garad / [[xxu(x) Vxo(x)I1 dx
u(x)
Vxu(x)
Fitted Gradient
Fitted Signal
Vxu(x)
uxx
Reconstructed Gradient
Reconstructed
Signal
NeuRBF
Ours
NeuRBF
Ours
Figure I: Comparisons to NeuRBF CCen al. 2023b for fitting 2D images with direct image supervision
(left) and gradient image supervision (right). Note that while NeuRBF fit signals well, it fails compute accurate gradients (bottom left) and, consequently, cannot recover images with gradient supervision.
A REPRODUCIBILITY DETAILS
@ Implementation. Our method is implemented in PyTorch, and all our experiments and the speed comparison for Siren are conducted
single H100 GPU, Following (Müller et al. 2022). each feature grid is initialised with small random values drawn from a uniform distribution, i.e. F, 14(-10- 10 2), whereas the decoder (whether linear MLP) weights are initialised with Xavier (Glorot& & Bengio 2010) and biases are set to zero. In the case of MLP decoder, we employ tanh as the activation. learning rate of 5 10-3 is used with the Adam optimiser. We set 1, which selects the neighbourhood (x) (i.e., p 3) around each query point x. Because ring (x) contains (2p) neighbours per scale, the Gaussian RBF therefore uses 6d neighbours every scale/grid. Neighbourhood indices Nx(x,) are precomputed for all training and test samples }x, Additional hyperparameters, including the scale S, grid resolution N Nuux: and feature dimension F, are experiment specific and detailed in the respective section for reproducibility.
K-Planes Comparison. We provide the settings used reproduce our K-Planes results with the offi- cial implementation (Fridovich- Keil et al.) (203) neural cloth simulation (Sec.4..3) (4.3) and Helmholtz equation (Sec. 4.2p. We extend K-Planes using the same losses as in our method. The architectures
summarised as follows:
NeuralClothSim: single-resolution 2D grid size 64 x 64 with feature dimension 16, followed by an MLP with two 64-wide ReLU layers and linear head; learning rate 10-2 with 0.1 decay every 2000 steps for 10000 epochs.
Helmholtz: single-resolution 2D grid of size 128 x 128 with feature dimension 32, followed by an MLP with two 64-wide ReLU layers and linear head; learning rate 2 x 10-5 for 50000 epochs. We employ single resolution the multi-grid K-Planes struggles converge otherwise.
15

---

Published as a conference paper at ICLR 2026
Image Reconstruction Lgrad
Eikonal Equation
Ground Truch
u(x)
PSNR
S
2306
S2 S-4 #params =213.47k
-0.3
02
S #pparams =526.34k #params 180
S-1
Truining Time (mins)
S-2 #params =009.59k
18
Figure II: Effectiveness of the multi-resolution grids. While the single- resolution feature encoder converges similar the multi-resolution grid in terms of the differential equation loss (as shown in the gradient loss plot), r -Grid multi-resolution architecture enables faster and more accurate field reconstruction (as illustrated in the PSNR plot, achieves higher PSNR in fewer iterations). The image reconstruction results are visualised at =30 seconds of training.
Ni-1
X--1 X;
Xi+1
Neel
(a) Lincar Interpolation
(b) NeuRBF Interpolation
(c) Our RBF interpolation
(Müller Figure III: Comparison interpolation schemes. Linear interpolation in Mrtant-NGPPK-AA
ct al.
2022 Fridovich- Keil et al. (023) is only C grid nodes x., whereas our RBF interpolation is C because centres fixed the grid with overlapping neighbourhoods. NeuRBF (Chen et al. 2023b also
RBFs but places adaptive centres off-grid and restricts each query single neighhourhood, leading discontinuities.
Instant-NGP Comparison. We use the PyTorch implementation (Tang 2022) Instant-NGP (that supports autograd) and apply the same losses in our method. The architecture details as follows:
NeuralClothSim: eight-level hash grid (base resolution 4, final resolution 64) with feature dimen- sion 2, followed by a linear head; learning rate 10 with 0.1 decay every 1000 steps for 10000 epochs.
Helmholtz equation: eight-level hash grid (base resolution 4, final resolution 128) with feature dimension followed by MIP with two 64-wide ReLUU layers and linear head; learning rate
2 x 10-5 for 50000 epochs.
SDF fitting: sixteen-level hash grid (base resolution 16, final resolution 2048) with feature dimen- sion 2, followed by an MLP with two 64-wide ReLU layers and a linear head; learning rate 10-4
for 20 epochs.
B FURTHER COMPARISONS
B.1 COMPARISON WITH NEURBF (CHN ET AL.. 2023B
NeuRBF and - both use RBF-interpolated feature grids, yet they target very different prob- lems. NeuRBF focuses on signal fitting: it reconstructs images and SDFs from ground-truth signals and tackles NeRFs via distillation, always coupled with hybrid backbone such Instant-NGP, K- Planes, or TensoRF. It does not attempt to solve PDEs directly from differential constraints, which precisely our setting.
Consequently, NeuRBF assumes access to the full target signal (or proxy provided by the hybrid model), whereas losses (Eqs. (8) and (1)) only require derivative supervision. We recover solu-
16

---

Published as a conference paper at ICLR 2026
tions from gradients, Laplacians, or oriented point clouds without seeing the ground-truth field. For images and SDFs, NeuRBF augments its adaptive RBFs with Instant-NGP features for C continu- ity; for NeRF, it first warms up K-Planes TensoRF for 1-2k steps and then distils their predictions to initialise the RBFs. Such proxy supervision is infeasible in our PDE setting because Instant- NGP/K-PPlanes themselves fail when trained with purely differential signals.
NeuRBF chooses adaptive RBFs to position and scale kernels based on the known signal, while we adopt fixed-shape Gaussian RBFs to guarantee smooth, differentiable interpolation suitable for PDE losses. Even if one tried to adapt NeuRBF to PDEs, several critical problems remain that we summarise below:
Signal-dependent initialisation. NeuRBF requires hand-crafted initialisation of RBF centres and shapes via weighted k-means. This biases the RBF distribution towards data points with higher weights (making their RBFs adaptive). Weights depend on the GT signal: w; = [xx())] for images, W; = 1/(SSDF (x:) 10-9) for SDFs. For NeRF, NeuRBF first warms up K- Planes or TensoRF for 1-2k steps and distils their predictions to initialise the RBFs as w; Such proxy supervision infeasible for PDEs, because Instant- NGP/K-Planes themselves fail in our settings, and no ground-truth field is available. By contrast, we place fixed-shape Gaussian RBFs on grid nodes without bespoke initialisation, which allows to additionally precompute the effective neighbourhood; see Fig. III for RBF placement of ours vs NeuRBF.
(1- (1-exp( ex.((o
Discontinuities in interpolation. NeuRBF restricts interpolation to a single RBF ring, causing sharp discontinuities whenever the active neighbour set changes; they lean on Instant-NGP to en- force C continuity grid nodes (Fig. MD. which is insufficient for PDE supervision. Expanding the neighbourhood is non-trivial because it would require KD-tree lookups at every query point RBF centres change, and their gegnal-dpeedent initialisation scheme is provided only for the single-ring case. PDE losses such as the Eikonal constraint demand smooth gradients across cell boundaries, so we precompute overlapping N3 neighhourhoods on the grid for arbitrary query points, guaranteeing continuous interpolation and stable higher-order derivatives (Figs. 3 ald III..
Empirical comparison. To check if NeuRBF can handle PDE supervision, we modify the of- ficial NeuRBF implementation to (i) replace its image-fitinng objective with our gradient-based Poisson loss (Eq. (8) and (ii) optimise the Eikonal loss (Fq. (13) for SDF reconstruction from oriented point clouds. We keep their initialisation scheme- weighting k-means by gradient mag- nitudes for images and by inverse SDF magnitudes for point clouds note that no analogous scheme exists for higher-order PDE supervision (e.g., Kirchhoff-Love problem for cloth simulation). Even under this favourable setup. NeuRBF failed to converge: Poisson reconstruction plateaued at PSNR/SSIM 10.85/0.17 versus 29. 44/.87 for ours (Fig. D. and the Eikonal experiment did not recover smooth surface (Fig. V These observations reinforce that NeuRBF's reliance hybrid architecture (Instant-NGP), signal-dependent initialisation and discontinuous interpolation make unsuitable for PDE solving, whereas our differentiable grid is explicitly engineered for derivative- only supervision.
B.2 COMPARISON WITH PINNS
We also compare against neural PDE solvers such PINNs (Raissi et al. 2019). While PINNs perform reasonably on smooth signals, they struggle with high-freyeency content, failing to capture wrinkles cloth simulations or oscillations Helmholtz. For instance, the image reconstruction experiment, GELUU-basd PINN achieves only PSNR/SSIM with our method; results for image reconstruction from gradient field (Sec. 4. and Helmholtz equations (.2c. are shown in Fig. IV This limitation arises from their use of smooth activation functions (e.g., Swish, GELU, tanh), which, though highly differentiable, are known to underfit high- frequency signals. We, therefore, adopt and primarily compare with Siren, which is better suited for such cases and widely used in physics-informat informed neural fields (Chen et al. 2023a Kairanda et al. 2024
compared to 29.44/00 /9.44/0.87
17

---

Published as a conference paper at ICLR 2026
C EIKONAL EQUATION
In addition the experiments presented in the main paper, demonstrate the learning of 3D shapes represented as signed distance func- tions (SDFs) by solving the Eikonal equation.
An SDF defines a mapping from 3D space to the distance of point from watertight face, with the zero level set corresponding to
PANS
the surface itself. Rather than assuming cess to ground-truth SDF values, we explore fitting SDFs directly oriented point clouds. This task corresponds to solving the Eikonal boundary value problem, non-linear PDE, that enforces the constraint [x)]] = al- most everywhere in the domain. Following the training procedure described in Sitzmann al.
Image Reconstruction
Helmholtz Wavefield
Figure Comparison PINN. compare 8o Grid to MLP with GELU activation function, which struggles with the high freeffreueency content.
2020). we formulate the loss function as:
Ceikonal [ 1| dx+
(|u(x)]]
(Vxu(x), n(x))) dx+
((((x)) dx.
(13)
where ! denotes the full domain, and S4 refers to the zero level set of the SDE The regularisation function defined as (+.) exp( |u(-))), with a >> 1, which penalises off-surface points that produce SDF values close to zero. As illustrated in Figs. [T] and (right), - is able to accurately solve the Eikonal equation and recover high-quality SDFs from oriented point clouds.
Although the encoder lives on spatio-temporal grid [-1, 1] we support irregular geometries by sampling PDE residuals and boundary terms on the target (potentially unstructured) domain f4. This arbitrary off- off-grid query sets (points and normals) can be supervised, enabling fits to surfaces (2D manifolds in 3D) from oriented point clouds. The trade-off is efficiency: empty grid points still store feature parameters, making our grid less parameter-efffccent than coordinate- MLPs.
works because uxx) in Eq. (5) is interpolated continuously,
C.I IMPLEMENTATION DETAILS
Dataset. In the SDF experiment, learn mapping u(x) [-1, - IR using the loss function deseribed above. The dataset D {(x;, nxx))] comprises samples on the point cloud x and their ground-truth surface normals n nx).
Architecture. The feature encoder comprises = 4 learnable grids, with the finest grid having a resolution of Nmax 56, and each grid possessing a feature dimension of F 1. Each query point is interpolated using 63 = 216 features per scale. The feature decoder is linear layer that maps the 3D input features to ID SDF.
Training. We employ stratified sampling over the domain, using 56 training samples along each input dimension. Additionally, in a total batch size of 2 x 56 A learning rate of 10-2 is used with the Adam optimiser. The SDF converges within 3000 iterations, requiring approximately six minutes of training time. For visualisation, the predicted SDF is evaluated on uniformly sampled 112 x 112 grid covering the same domain, and meshes are extracted using the Marching Cubes algorithm.
equal number of points sampled the surface, resulting
C.2 COMPARISON WITH GRID-BASED METHODS
Many grid-based methods Müller et al. 2022: Li et al. 2023) represent signed distance functions (SDPs) using linear interpolation. In contrast to solving the Fikonal equation from oriented point clouds as in Eq. (13), these methods usually overfit to ground-truth SDFs computed from meshes, typically using a loss of the form:
Lalr = / |lu(x) SDF(x)] dx.
(14)
18

---

Published as a conference paper at ICLR 2026
L211
Input: Oriented Point cloud
NeuR
Instant-NGP
Ours
Figure V: Solutions to the Eikonal equation from oriented point clouds. For each method, we visualise 2D slices (left) and the mesh extracted from the SDF (right). Instant-NGP and NeuRBF can overfit when ground- truth SDF supervision (Car) available, but they cannot reliably solve the Eikonal equation given partial point- cloud observations (Ceiooaa) because their gradients are not differentiable across grid boundaries (Fig. [. In contrast, our smooth formulation handles these gradients and solves the Eikonal problem comparably to baselines such as sinusoidal MLP.
where the loss penalises deviations from the ground-truth SDE. This approach often used either directly (Müller et al. 2022 Chen et al. 2023b) or indirectly, for example, in multi-view surface reconstruction (Li al. 2023 Wang et al. 2023 (sometimes with additional regulariser enforcing the Eikonal condition). While such methods can produce accurate reconstructions when provided with high-quality ground- truth SDFs, they struggle when solving the Eikonal equation directly using feikonal (Eq. (13).. In particular, the Eikonal loss is not smoothly optimised across grid boundaries when using linear interpolation, which often results in incorrect SDF values (e.g., incorrect sign) on either side of the boundary. This issue is visualised in 2D slices for Instant-NGP in Fig, V In contrast, our method employs differentiable grid interpolation, enabling the recovery of consistent and smooth SDFs directly from oriented point clouds.
D POISSON EQUATION
D.1 IMPLEMENTATION DETAILS
Dataset. We use the Camera image from sci-kit (Pedregosa et al. 2011). Following Siren (Sitz- mann al. 2020), the ground-truth gradient image is computed using the Sobel filter and scaled by a constant factor of 10 for training. The ground-truth Laplacian image is computed using Laplace filter and scaled by factor of 10k.
Architecture. The feature encoder comprises S = 4 learnable grids, with the finest grid of the resolution Nmax = 512, and each grid of the feature dimension F = We vary these hyperparame- ters for gradient and Laplacian supervision; see Tab, for the number of trainable parameters. Each query point is interpolated using 62 36 features per scale. The feature decoder is a linear layer.
Training. We train by evaluating every pixel of the gradient Laplacian image at each it- eration. With regular fixed sampling over the domain, iee,, 512 training samples along each input dimension, results in total batch size of 5122. A learning rate of 10-3 is used with the Adam optimiser. The convergence varies per experiment and is visualised in Fig. and Tab.
E HELMHOLTZ EQUATION
We describe the full details of the Helmholtz experiment, as provided in Sec. 4.2 of the main paper. The Helmholtz equation is second-order partial differential equation commonly used to model steady-state wave behaviour and diffusion phenomena. It is expressed as:
(A + mx) w²) u(x) = -g(x),
(15)
where g(x) denotes known source function, u(x) is the unknown complex-var valued wavefield, and m(x)= = ã represents the squared slowness, defined in terms of the wave velocity c(x). We solve for the wavefield using the -Grdd representation, with results presented in Figs. and 5 and
19

---

Published as a conference paper at ICLR 2026
Real Error
Imaginary Error
3 mins 30 secs
21 mins 42 secs
Converged
mins 30
6 mins 35 Converged
Siren
Ours
Figure er-er for solution to the Helmholtz equation. As presented Fig. while both ours and
Sitzmann et al. 2020) closely match the reference solution, our method trains significantly faster.
error visualisation in Fig. VI In the following, we outline the specific formulation of the Helmholtz equation variant used in our implementation.
E.1 PERFECTLY MATCHED LAYERS FORMULATION
To obtain unique solution to the Helmholtz equation within bounded domain, we adopt per- fectly matched layer (PML) strategy, following the formulation introduced in (Chen et al. 2013] Sitzmann et al., (0200). The PML serves to absorb outgoing waves at the domain boundaries, thereby preventing artificial reflections. This results the following variant of the Helmholtz equation:
a 0x
du(x) dx
Ōu(x)
ez, Ca2
+
Cs, Exa
+ exe Ezgh kack ux(x) = -g(x),
(16)
Ox2
where x (I1, 22) E Ω, the complex coefficients given by Cz2 = 11-1 1-1-jt0 j and k w/c. The damping term oz. along each spatial axis specified as:
anw
Lrg Texa
if r; E an, otherwise.
Ox,
=
(17)
where dn determines the attenuation strength (set to ag 5), and c = is the wave propagation speed. The term 4,, represents the distance to the nearest PML boundary along the axiaxis while LyMi, denotes the total PML width. The PML applied exclusively near the domain boundary, defined as OΩ fx 0.5 < ||x|| < 1). and the original Helmholtz equation remains unchanged elsewhere. We optimise Equation [6] using the loss function introduced in the main text (Eq. 9-
(main)), with the spatial weighting factor set A(x) batch
5x100"
E.2 IMPLEMENTATION DETAILS
Dataset. In the Helmholtz experiment, we learn mapping u(x) [-1, -> R using the loss function described above. The dataset D = ((x; g(x;))) comprises sampled spatial coordinates x; and values of the Gaussian source function, defined g(x) (xx [0, 0), 10 10-4).
Architecture. The feature encoder comprises S 4 learnable grids, with the finest grid having resolution of Nmux 128, and each grid possessing feature dimension of F = 2. Each query point is interpolated using 6 = 36 features per scale. The feature decoder is an MLP that takes the concatenated features input, projects them to hidden layer of width 64 with Swish activation, and produces a 2D output through a final linear layer. The total number of trainable parameters is approximately 45.10 19k:
20

---

Published as a conference paper at ICLR 2026
Wavenumber
Wavenumber w
Re Reference
Real Frror
Imag Error
20
0.0008
0.0007
30
0.0040
0,0040
(x)
[2)
40
0.0250
0.0240
Ours
Table I: Helmholtz equation with higher wavenumbers: We report & eтror with the reference solution for the real and imaginary parts of the wavefield, Parameter count: 45.19k for w 20, 30; 176.7
Real
Imaginary
Real
Imaginary
Figure VII: Results Helmholtz equation with w = 30, 40 compared with closed-form solutions. Our 8 -Grid methad remains accurate overall, with minor errors near the source at 40.
Training. We employ stratified sampling over the domain, with 256 training samples along each input dimension, resulting in a total batch size of 256 learning rate of 10-3 used with the Adam optimiser. The Helmholtz wavefields converge within 2500 iterations, requiring approx- imately 6 minutes of training time. For visualisation, the predicted wavefield is evaluated on uniformly sampled 256 256 grid covering the same domain.
E.3 HIGHER WAVENUMBER
To further evaluate our method under more challenging settings, we reproduce the Helmholtz equa- tion experiments with higher wavenumbers,
While Fig. [ presented results for 20, here we extend the experiments to w = 30 and w = 40. Comparisons with the closed-form reference solutions are shown in Fig. VII and Tab. For w 30, our method closely matches the analyti- cal solution, whereas for = 40 we observe noticeable errors near the wave source but otherwise maintain good accuracy across the domain.
F CLOTH SIMULATION
We describe the further details of the cloth simulation experiment, presented in Section 4.3 of the main paper. Given the rest state the cloth and the material elastic parameters L, quasistatic cloth simulation involves finding the equilibrium deformation field u Ω - R of the midsurface under the influence of external forces g and boundary constraints on 02. The stable equilibrium configuration is obtained by minimising the total potential energy functional subject Dirichlet boundary constraints:
u' arg min fa v(e(u(x)), m(u(x)); -g - dx,
(18)
u x2) b on an,
where the internal elastic energy density, and E, represent the membrane and bending strain components, respectively. We solve for the deformation field using the o representation. Results are shown in Figures and 6 of the main paper.
F.1 IMPLEMENTATION DETAILS
Dataset. In the simulation experiment, we learn mapping u(x) []1 -> R33 using the loss function described above. The dataset D consists of sampled spatial coordinates x and their corre- sponding rest state positions, set as [z), 22, 0J. The material properties are defined by the parameter set Φ, with values:p = 0.144, h = 0.0012, E 5000, and // = 0.25. We impose Dirichlet boundary conditions as hard constraints (following NeuralClothSim (Kairanda al. 2024 by displacing the tep-left and top-right vertices the cloth towards the centre of the input domain by 0.4 units each.
TTe external force is modelled as constant gravitational acceleration [0. -9.8, 0).
21

---

Published as a conference paper at ICLR 2026
(=0s
4s
PDE Residual
MAE
1.0
Data
35-
Ours
Our
0.5
1.3
2.0
Snan
0.0
20
1.5
-1.
-14
Graund
05
-10
-53
0.5
1.0
-1.0
-05
0.0
055 100
Heat Field
Training Time (mins)
Figure VIII: Solutions of the Heat Equation. We solve the spatio-temporal heat equation successfully, where an initial sinusoidal wave diffuses over time (visualised = 48), Compared to a coordinate- bascd MLP (Sitz-) mann et al. 2020), ours converges to the analytical solution more accurately and efficiently, as demonstrated with PDE residual and mean absolute enror of the heat field.
Table II: ID Advection comparison with INSR (Chen et al
2023a and Grid solver (Fedkiw et al. 2001
Method
Error (MAE)
Training Time
Memory
INSR Chen al.
2023a)
0.0030
5.33h
3.53KB
Grid Fedkiw et
2001
0.0029
1.80s
27.35KB
Our Grdd approach
0.0023
1.5mins
16.52KB
Architecture. The feature encoder consists of a single learnable grid (S = 1) with feature dimension of F 4. We use grid resolution of N 32 for the simulation in Fig. 1-(main), and N = 64 for that in Fig. 6-(main). Each query point interpolated using 62 = 36 features per scale. The feature decoder is linear layer that maps the 3-dimensional input features to 3I) deformation field. The total number of trainable parameters is approximately 4.37k for the results shown in Fig. 1-(main), and 16.91 for those in Fig. 6-(main).
Training. Similar to the Helmholtz experiment, we employ stratified sampling over the domain, using 64 and 128 training samples along each input dimension for Figs. and 6-(main), respectively. A learning rate of x 10-2 is used with the Adam optimiser. The cloth simulation converges within 500 iterations, requiring at most 5 minutes of training time. For visualisation, the predicted deformation field is evaluated on a uniformly sampled 64 x 64 grid spanning the same domain.
G
HEAT EQUATION
Here, we additionally demonstrate our result for the ID heat equation. The heat equation is diffu- sion PDE used model temporal smoothing signals:
Ou 8u
(19)
0 0.2
where u(æ, t) represents the temperature at position x and time t, and a is the thermal diffusivity constant. We consider spatial domain x [-1, 1] and temporal domain E [0, 4] with Dirichlet boundary conditions and an initial sinusoidal condition u(I; 0) sin(TT) with 1. It admits an analytical solution:
uxx, at sin(()))
(20)
which we use to quantitatively evaluate reconstruction fidelity. We train with the PDE loss with hard boundary constraints. Fig. VIII summarises the experiment: the sinusoidal wave diffuses as expected and matches the analytical solution with a mean absolute error of 0.004 over uniformly sampled space-time points. Relative to the coordinate- based MLP baseline (Sitzmann et al., 2020), the PDE residual and the heat-field error reduce faster. Thus, - faithfully recovers the spatial and temporal evolution of the diffusion process.
22

---

Published as a conference paper at ICLR 2026
4
F=
/=
u1x 0)
u(x; O)
DFx. ()
01x 01
x1
ID Gaussian Wave
2D Gaussian Wave
2D Zalesak's Disk
Figure IX: solve spatio-temporal advection successfully for various ID and fields, including and 2D Gaussian waves as well as 2D Zalesak's disk.
H ADVECTION EQUATION
We consider the spatio-temporal advection equation
+ vxu=0,
(21)
Øt
which describes the transport of field with advection velocity a. We demonstrate method on ID and 2D cases, including Gaussian wave propagation and the Zalesak benchmark with a challenging discontinuous boundary.
H.1 ID ADVECTION
2023a). where Gaussian-shaped wave, centred at - 1.5 propagates rightwards the spatial domain [- 2, 2]. with advection velocity a = 0.25. In contrast to INSR, impose boundary and initial conditions hard constraints and treat time continuous input to our feature grid, rather than discretising it (they use midpoint integration). We present the results in Tab. M and Fig. [X] Our results closely match both the analytical solution and INSR, with mean absolute error (MAE) reported at = 4 seconds. The INSR and the finite difference-basee grid solver (Fedkiw et al.) 2001) values are taken from Tab. 1 of (Chen al., 2023a). - suffers minimal dissipation; we observe an MAE of 2.3 x 10-3 at 4s, i.e., less than 0.3% amplitude loss compared the analytical solution, while still training 200x faster than INSR. Overall, our approach achieves lower error relative to the ground truth, reduced training time compared to INSR, and lower memory use compared numerical solvers.
We follow the INSR setup (Chen et al.
H.2 2D ADVECTION
For the 2D Gaussian wave, we adopt the domain and initial conditions of Chetan al. 2025 sim- ulating for 4 seconds. Relative to the closed-form solution, our method achieves mean advection enrors of 0.0047 (after minutes), 0.0032 (after minutes), and 0.00246 (after 15 minutes, fully converged); see Fig. (X-ntee))e) for qualitative results. Note that in Chetan et al. (2025) advec- tion example, they employ forward/explicit Fuler (possible with post hoc gradient operator) for timestepping, unlike INSR (Chen et al.) 2023a) (which can support mid-point due to underlying Siren representation) and ours (which supports auto diff and mid-point similar to Siren/INSR).
H.3 ZALESAK'S S DISK PROBLEM
For the Zalesak benchmark, follow the setup in Dupont & Liu (2003). and solve for the signed distance function (SDF) of the slotted disk under advection. The results converge and preserve the overall shape. although some numerical diffusion is visible blurring along the slot edges; Fig. (right). After half a revolution, we obtain an fi-error of 0.065 compared to the closed-form rody rotation solution with constant angular velocity.
I LIMITATIONS AND FUTURE WORK
While method leverages differentiable interpolation to enable accurate representation of signals and their derivatives, exhibits several limitations. First, the interpolated features are sensitive to the choice of the Gaussian RBF shape parameter, which in turn determines the effective neighbourhood.
23

---

Published as a conference paper at ICLR 2026
As result, the representation introduce some degree of smoothing, depending on this config- uration. While we observed convergence across all experiments, further investigation is needed to establish theoretical guarantees for C under neighborhood truncation. A well-studied drawback of RBF methods, particularly near the domain boundaries, is the occurrence of boundary errors (Runge's phenomenon), which negatively impact reconstruction accuracy in those regions. Our grid approach also suffers from the curse of dimensionality; for example, fitting SDFs in 3D is no- tably slower than in 2D. The efficiency gain our method is primarily due to the use of structured Cartesian grids, which is lost when solving PDEs on manifolds. As direction for future work, we suggest exploring planar projection strategies improve scalability and computational efficiency in high-dimensional domains.
Moreover, our novel representation, although it improves speed over neural solvers and introduces differentiabilityy into feature grids, it might not be as efficient traditional numerical solvers, such as multi-grid methods. Since this requires further thorough investigation, in our evaluations, we provided initial experiments and comparisons to numerical methods help guide future work in this direction. In the future, custom CUDA kernel implementation of do--Grid could further accelerate optimisation compared to current Py Torch implementation.
24


Reviews:
[{"review_id": "Lmho9Uxl2n", "reviewer": "Reviewer_ai8J", "summary": "The paper proposes a fast neural PDE solver based on grid-based neural fields like Instant NGP. Traditional grid-based neural fields are neural fields with features on a grid that are modulated by a small MLP. For arbitrary query points, these fields compute the feature by using d-linear interpolation. However, relying on d-linear interpolation means that these networks do not yield meaningful high-order derivatives. The authors propose to mitigate this by proposing a differentiable RBF-based interpolation that can yield accurate higher-order derivatives. The authors use a Gaussian RBF kernel to compute interpolation weights for the grid points surrounding a query point. \n\nThe authors apply their method to solving a variety of PDEs, ranging from the Helmholtz equation, the Advection and Heat Equation, Eikonal Equation, etc, showing increased accuracy over contemporary grid-based neural field architectures like Instant NGP and increased efficiency as compared to MLP-based architectures like SIREN.", "questions": "1. Appreciate the comparison to NeuRBF, but could the authors shed some light on how NeuRBF's formulation differs from their approach, which makes it possible for their approach to be used as an effective neural solver? The current explanation is unclear.\n\t- As per the authors' explanation it seems that NeuRBF's design makes it suitable to fit zeroth-order signals but not gradient-/derivative-based residuals like the losses described by the authors in their experiments which makes it difficult to use NeuRBF as a neural solver, but I am still not convinced by the author's explanation of on what specific design detail of NeuRBF makes it infeasible to use it as a neural solver.\n\t- Currently, the authors point to the requirement of ground-truth signals by NeuRBF as an issue, but that is because NeuRBF has been floated as a signal representation. Even if the authors' proposed architecture were to be used as a signal representation, for instance, to learn the initial conditions of a PDE, they would require the ground-truth signal for training.\n2. Can the author's architecture be used for overfitting to any input signal as well? This could increase the scope of their contribution. \n\t- Asking because Instant NGP reported that moving to higher-order interpolation approaches hurt performance on reconstruction \n\t- If the authors' proposed approach can do better would be a great contribution to the community. Even if not, it would be something useful for the community to know.\n\nI am generally leaning towards accepting this work. If the authors can address the weaknesses and provide an explanation for the above questions, I think that would make this work worthy of being accepted.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The problem of obtaining accurate higher-order derivatives from grid-based neural fields is an important one, and the authors' differentiable interpolation scheme is a great step in this direction.\n2. The authors have done a rigorous evaluation across a broad range of PDEs.", "weaknesses": "Weaknesses\n\n1. Writing can be improved: \n\t- Some of the authors' claims about their contributions seem a bit exaggerated. For instance, multi-resolution grids were an idea that was proposed in both Instant NGP and K-planes. While the authors have applied this idea to a novel setting of PDE-solvers, claiming the idea itself as a contribution seems misleading to me.\n\t- Notational issue: it seems that $r$ is both a parameter (input to the RBF kernel) and a hyper-parameter (extent of neighborhood considered for kernel) - this should be either resolved or the explanation of the $r$ hyper-parameter should be improved.\n\t- For 2D Advection, the authors claim that their approach shows no numerical dissipation. However, the change in the peak value on the colorbar between the initial value and the final state shows that there has been some dissipation. The authors have also not shown quantitative results for 1D and 2D advection\n2. Missing Results?: Could not find results for the heat equation in the supplementary material. (Appendix G)", "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": "eqvpFBR1UY", "reviewer": "Reviewer_qdmK", "summary": "The paper proposes a differentiable grid-based representation for solving differential equations. By combining multi-resolution feature grids with smooth radial basis function (RBF) interpolation, the method enables higher-order differentiability, allowing it to learn continuous fields whose derivatives satisfy physical constraints. The approach is applied to a variety of PDE problems, including the Poisson equation, the Helmholtz equation, and the Kirchhoff–Love model, achieving faster convergence than coordinate-based neural solvers such as SIREN.", "questions": "- How does the accuracy or convergence rate scale? Have the authors attempted any PDEs that require actual 3D fields (as opposed to the 2D-to-3D mapping in the cloth simulation)?\n\n- Can the method handle PDEs with non-smooth solutions or discontinuous coefficients, where RBF smoothness may become a drawback?\n\n- In the Eikonal equation experiment in the Appendix, how is the proposed method now different from NeuRBF? Would the authors be able to supplement a comparison with NeuRBF regarding this specific setup?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 2, "strengths": "- The paper is written with clarity, organization, and completeness. Overall, the derivation is straightforward and easy to follow.\n\n- The paper includes a comprehensive and thorough list of experiments, including Poisson Equation, Helmholtz Equation, Cloth Simulation, Eikonal Equation, and Advection and Heat Equations, and with both qualitative and quantitative results. \n\n- I appreciated the numerous experiments the authors supplement in the Appendix, especially the section on Eikonal Equations.\n\n- The method achieves comparable or better accuracy than MLP-based PINNs and SIREN while providing 5–20× faster training and inference.", "weaknesses": "- The use of RBF interpolation in feature grids is conceptually similar to prior work such as NeuRBF (Chen et al., ICCV 2023). Both frameworks ultimately learn continuous fields under differentiable consistency constraints, and while the authors claim NeuRBF \"adapts to known target signals,\" this is not the case for NeRF-style trainings where the supervision comes from images and the ground truth 3D field is also unknown. I would sincerely suggest that the authors acknowledge this part and limit their contribution to \"applying existing RBF feature grid to PDEs.\"\n\nThis is my main assessment of the concerns of the paper, but from a \"solving PDEs\" perspective, there is still merit in having this paper published.", "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": "vNwtL6M1T9", "reviewer": "Reviewer_mShH", "summary": "The submission introduces ∂∞-Grid, a differentiable grid-based neural representation for solving partial differential equations (PDEs).Unlike coordinate-based multilayer perceptrons (MLPs) such as PINNs, ∂∞-Grid employs Radial Basis Function (RBF) interpolation over learned feature grids, resulting in an infinitely differentiable representation that supports analytic computation of spatial derivatives.\nThe framework further integrates multi-resolution, co-located grids to facilitate gradient propagation across scales.\n\nEmpirical evaluations encompass a range of PDE problems including gradient-based reconstruction, Helmholtz, Kirchhoff–Love plate deformation, Eikonal, and advection equations and demonstrate 5–20× faster convergence than coordinate-based baselines, with comparable or improved accuracy.\nThe approach aims to unify the strengths of grid-based representations, efficient yet non-differentiable, and coordinate-based networks, which are flexible but computationally slower.", "questions": "1.\tCan ∂∞-Grid be extended to non-rectangular or manifold domains?\n2.\tCould the framework be combined with coordinate-based PINNs to handle complex geometries while preserving grid efficiency?\n3.\tCould the authors clarify the implementation details for both K-Planes and Instant-NGP? For example, whether K-Planes is formulated as\n$u(x) = \\mathrm{MLP}(\\mathrm{linear_interpolation}(F, x))$.\n\n5.\tCould the authors evaluate or discuss the sensitivity of the proposed approach to grid resolution, the number of scales, and the dimensionality of the feature vectors?\n\n6.\tFigure 5: Could the authors include an image illustrating the difference between SIREN and the proposed method?\n7.\tWhat specific classes of PDEs can the method address effectively, and under what assumptions regarding smoothness or boundary conditions?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "1.\tThe combination of feature grids with Gaussian RBF interpolation is intersting, technically sound and conceptually straightforward. Since Gaussian RBFs are $C^{\\infty}$, the method inherently supports analytic higher-order derivative computation, which is important for PDE solvers.\n2.\tThe experiments report significant improvements in convergence rate and accuracy, achieving 5–20x faster training compared to SIREN and K-Planes. These speedups represent a substantial computational gain, reducing training time from hours to minutes or even seconds, and demonstrate the method’s scalability and efficiency across diverse PDE formulations.\n3.\tThe presentation is clear, and the methodology is described in a structured and coherent manner.\n4.\tThe experimental evaluation across multiple classes of PDEs suggests that the proposed approach is broadly applicable rather than tailored to a specific problem type.", "weaknesses": "1.\tThe formulation assumes a Cartesian grid. In contrast to coordinate-based PINNs, ∂∞-Grid cannot directly address irregular geometries or complex boundary conditions without additional coordinate transformations or masking procedures, limiting its applicability to structured domains.\n2.\tThe comparisons focus mainly on SIREN and K-Planes. Incorporating recent operator-learning approaches such as the Fourier Neural Operator or DeepONet would strengthen the evaluation, as the proposed framework requires solving each new boundary condition or PDE configuration from scratch, unlike operator-based methods that generalize across parameter variations.\n3.\tThe comparison to K-Planes may not be entirely fair for derivative-based tasks, as K-Planes relies on non-differentiable piecewise linear interpolation. \n4.\tThe paper lacks an analysis of how the solver’s performance varies with grid resolution, the number of multi-resolution levels, and the dimensionality of the feature vectors. It is also unclear how these parameters are selected.\n5.\tWhile the paper presents interesting empirical results, a more formal treatment of aspects such as convergence, stability, or error bounds relative to classical numerical solvers or PINNs could strengthen the contribution. Clarifying the underlying functional space (e.g., RKHS, $C^2$) would also enhance the theoretical grounding.\n\nMinor comments\n\n1.\tThe symbol f is used for different quantities in Eqs. (1) and (2). Although boldface differentiates the latter, this overlap could cause ambiguity.\n2.\tThe “Poisson” experiment corresponds to a gradient or Laplacian-based reconstruction task rather than a direct PDE solution; clearer terminology would help prevent confusion.", "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 (Poster)