KernelOnet：一种基于核函数的可解释神经算子

# KernelOnet: An Interpretable Neural Operator Based on Kernel Functions
Source: [https://arxiv.org/html/2609.35938](https://arxiv.org/html/2609.35938)
Yuan GuoHanshu ChenAffiliation:College of Mechanics and Engineering Science, Hohai University, Nanjing 211100, ChinaQiang XiAffiliation:College of Mechanics and Engineering Science, Hohai University, Nanjing 211100, ChinaTimon RabczukAffiliation:Institute of Structural Mechanics, Bauhaus\-University Weimar, Weimar 99423, GermanyZhuojia Fu††thanks:Corresponding author\.paul212063@hhu\.edu\.cnAffiliation:College of Mechanics and Engineering Science, Hohai University, Nanjing 211100, ChinaAffiliation:Key Laboratory of Ministry of Education for Coastal Disaster and Protection, Hohai University, Nanjing 210098, China

###### Abstract

This paper proposes a new interpretable neural operator framework, termed the Kernel Operator Network \(KernelOnet\), which explicitly incorporates kernel functions into the neural operator architecture so that the operator structure is consistent with the kernel\-expansion form used in boundary\-type kernel\-expansion methods\. Unlike traditional neural operators such as DeepONet, which rely on deep networks to learn basis functions implicitly, KernelOnet replaces the trunk network with explicit kernel functions and provides three complementary kernel construction strategies: the first is a data\-driven learnable kernel, in which a radial basis function is parameterized by a neural network and learned directly from data; for constant\-coefficient linear problems, the learned kernel can be regarded as a non\-singular fundamental solution; the second is a physics\-informed kernel, which explicitly embeds physical information such as analytic fundamental solutions into the network structure, so that the expansion automatically satisfies the governing equation and can be trained without supervision using boundary conditions alone, without any interior solution data; the third is a hybrid kernel, which splits the solution according to the linear principal part of the governing equation into a homogeneous part spanned by analytic fundamental solutions and a source part carried by low\-rank learned correction kernels, thereby balancing physical priors against data\-fitting capability on nonlinear problems for which no analytic fundamental solution exists\. On three benchmark problems and one engineering problem in a shallow\-water waveguide, KernelOnet attains high\-accuracy solutions; on the examples that can be compared directly with DeepONet it achieves higher accuracy with fewer learnable parameters\. Its unsupervised configuration can be trained without any interior solution labels, its per\-query inference cost is far below that of per\-instance solvers, and it offers an effective route to acoustic propagation in unbounded exterior domains, which general\-purpose neural operators struggle to handle\.

*Keywords*Neural Operator⋅\\cdotKernel Function⋅\\cdotDeepONet⋅\\cdotFundamental Solution⋅\\cdotShallow Water Acoustics

## 1Introduction

Partial differential equations \(PDEs\) are ubiquitous in many fields such as physics, engineering and materials science, and their efficient numerical solution has long been one of the core problems of computational mathematics\. Traditional numerical methods, such as finite difference methods[Smith \(1985\)](https://arxiv.org/html/2609.35938#bib.bib57), finite element methods[Reddy \(1984\)](https://arxiv.org/html/2609.35938#bib.bib58)and various meshless methods, generally require discretizing the solution domain and approximating the solution by solving large\-scale linear or nonlinear algebraic systems, and their computational cost grows sharply with the problem size and dimension; moreover, for tasks with a fixed geometry the system must still be reassembled and solved for every new set of boundary conditions or source terms, which makes it difficult to meet the efficiency requirements of emerging applications such as real\-time prediction, inverse problems and uncertainty quantification\.

In recent years, the rapid development of deep learning[LeCun et al\. \(2015\)](https://arxiv.org/html/2609.35938#bib.bib51)has brought a new paradigm to scientific computing, and neural operators, as a new class of learning frameworks, have attracted widespread attention\. Their goal is to learn directly from data the mapping between function spaces, that is, the operator itself, so that given a new input function such as a boundary condition, a source term or an initial condition, the solution function can be produced directly by a single forward pass without solving the PDE again\. Within the family of neural operators, DeepONet[Lu et al\. \(2021a\)](https://arxiv.org/html/2609.35938#bib.bib1)constructs the solution operator on the basis of a branch–trunk architecture together with the universal approximation theorem for operators[Chen and Chen \(1995\)](https://arxiv.org/html/2609.35938#bib.bib48), and subsequent work has further introduced physical information to improve accuracy and generalization[Wang et al\. \(2021\)](https://arxiv.org/html/2609.35938#bib.bib11); this physics\-informed line of thought has since developed further: on the one hand, the weak form of the PDE is incorporated into operator training in a variational manner, giving rise to the variational physics\-informed neural operator \(VINO\)[Eshaghi et al\. \(2025\)](https://arxiv.org/html/2609.35938#bib.bib54); on the other hand, the learned operator is used for the pretraining and warm\-starting of iterative solvers such as finite element methods, so as to accelerate conventional numerical solution procedures[Wang et al\. \(2026\)](https://arxiv.org/html/2609.35938#bib.bib55); the Fourier neural operator \(FNO\)[Li et al\. \(2020a\)](https://arxiv.org/html/2609.35938#bib.bib2)and its generalizations perform global convolution in the spectral domain and have produced variants such as multipole graph structures[Li et al\. \(2020b\)](https://arxiv.org/html/2609.35938#bib.bib3), learnable deformations on general geometries[Li et al\. \(2022\)](https://arxiv.org/html/2609.35938#bib.bib4)and large\-scale weather forecasting[Kurth et al\. \(2023\)](https://arxiv.org/html/2609.35938#bib.bib17); operator learning based on Transformers[Cao \(2021\)](https://arxiv.org/html/2609.35938#bib.bib9);[Hao et al\. \(2023\)](https://arxiv.org/html/2609.35938#bib.bib7)and network layers based on Clifford algebras[Brandstetter et al\. \(2022\)](https://arxiv.org/html/2609.35938#bib.bib10)have also appeared; in addition, tensor\-product methods for operators with multiple inputs[Jin et al\. \(2022\)](https://arxiv.org/html/2609.35938#bib.bib6), multi\-head neural operators for multifield and interface\-dynamics modelling[Eshaghi et al\. \(2026\)](https://arxiv.org/html/2609.35938#bib.bib56), Green’s function learning for nonlinear boundary value problems[Gin et al\. \(2021\)](https://arxiv.org/html/2609.35938#bib.bib16), and physics\-informed neural networks \(PINNs\) that embed PDE information into the loss function[Raissi et al\. \(2019\)](https://arxiv.org/html/2609.35938#bib.bib12)together with their follow\-up work[Karniadakis et al\. \(2021\)](https://arxiv.org/html/2609.35938#bib.bib13);[Cuomo et al\. \(2022\)](https://arxiv.org/html/2609.35938#bib.bib14);[Lu et al\. \(2021b\)](https://arxiv.org/html/2609.35938#bib.bib15)jointly form a flourishing landscape of operator learning and physics\-informed machine learning, for which a review can be found in[Kovachki et al\. \(2021\)](https://arxiv.org/html/2609.35938#bib.bib8)\. However, most of the above methods rely on deep networks to learn basis functions or operator kernels implicitly, and their structure lacks an explicit connection with the specific governing equation, which raises two problems: first, poor interpretability, since it is difficult to understand how the model works from its internal structure; second, insufficient physical consistency, since a large amount of high\-fidelity training data is usually required, and generalization is limited in the small\-sample regime\.

In contrast, kernel methods and numerical methods based on kernel expansions have constructed kernel functions explicitly from the very beginning\. In machine learning, the theory of reproducing kernel Hilbert spaces \(RKHS\)[Aronszajn \(1950\)](https://arxiv.org/html/2609.35938#bib.bib19)provides a rigorous mathematical foundation for kernel functions, and the success of support vector machines[Cortes and Vapnik \(1995\)](https://arxiv.org/html/2609.35938#bib.bib25), Gaussian process regression[Williams and Rasmussen \(1996\)](https://arxiv.org/html/2609.35938#bib.bib29);[Rasmussen and Williams \(2006\)](https://arxiv.org/html/2609.35938#bib.bib22)and general kernel methods[Hofmann et al\. \(2008\)](https://arxiv.org/html/2609.35938#bib.bib26);[Schölkopf and Smola \(2002\)](https://arxiv.org/html/2609.35938#bib.bib27)shows that embedding data into a suitable function space through a kernel function yields good generalization performance; radial basis function \(RBF\) networks[Park and Sandberg \(1991\)](https://arxiv.org/html/2609.35938#bib.bib21)and regularization network theory[Girosi et al\. \(1995\)](https://arxiv.org/html/2609.35938#bib.bib28)have further revealed the intrinsic connection between kernel methods and neural networks\. In numerical methods, radial basis function interpolation theory[Micchelli \(1986\)](https://arxiv.org/html/2609.35938#bib.bib20)provides the convergence foundation for meshless and boundary\-type kernel methods; the method of fundamental solutions \(MFS\)[Golberg \(1995\)](https://arxiv.org/html/2609.35938#bib.bib36);[Chen et al\. \(1998\)](https://arxiv.org/html/2609.35938#bib.bib37);[Fairweather and Karageorghis \(1998\)](https://arxiv.org/html/2609.35938#bib.bib38);[Fairweather et al\. \(2003\)](https://arxiv.org/html/2609.35938#bib.bib39);[Golberg and Chen \(1999\)](https://arxiv.org/html/2609.35938#bib.bib40), the Trefftz method[Kita and Kamiya \(1995\)](https://arxiv.org/html/2609.35938#bib.bib41), the boundary element method \(BEM\)[Brebbia et al\. \(1984\)](https://arxiv.org/html/2609.35938#bib.bib42)and various meshless methods[Belytschko et al\. \(1996\)](https://arxiv.org/html/2609.35938#bib.bib32);[Liu \(2009\)](https://arxiv.org/html/2609.35938#bib.bib33);[Fasshauer \(2007\)](https://arxiv.org/html/2609.35938#bib.bib34);[Monaghan \(2005\)](https://arxiv.org/html/2609.35938#bib.bib35)all employ kernels that contain physical information, such as fundamental solutions, as basis functions in order to achieve high accuracy and high computational efficiency\. Schaback and Wendland[Schaback and Wendland \(2006\)](https://arxiv.org/html/2609.35938#bib.bib30)pointed out that the kernel techniques used in machine learning and the kernel functions used in meshless and boundary\-type kernel methods are in fact the same mathematical object manifested in two fields\. However, traditional kernel methods and kernel\-expansion\-based numerical methods usually rely on a manually chosen kernel function and discretization scheme, such as the source\-point layout of a boundary collocation method or the choice of the virtual boundary for a boundary integral, and are difficult to adjust automatically for complex nonlinear problems; and mainstream neural operators learn basis functions implicitly, sacrificing physical consistency and interpretability\. This observation motivates us to combine the two and to propose an operator learning framework that uses explicit kernel functions as the basis functions of a neural network\.

To address the above issues, this paper proposes the Kernel Operator Network \(KernelOnet\), a direct improvement over DeepONet: it replaces the trunk network with explicit kernel functions, so that the operator output structure is consistent with the kernel\-expansion form used in boundary\-type numerical methods based on kernel expansions, such as the boundary\-type collocation method of fundamental solutions \(MFS\) and the boundary integral boundary element method \(BEM\), while remaining amenable to analysis within the well\-developed theory of kernel methods, such as radial basis function interpolation theory and reproducing kernel Hilbert space theory\. Specifically, this paper provides three complementary kernel construction strategies: \(i\) a data\-driven learnable kernel \(KernelOnet\-RBF\), whose kernel is a neural\-network\-parameterized radial basis function learned directly from data by a shallow network, so that for constant\-coefficient linear problems the learned kernel can be regarded as a non\-singular fundamental solution; \(ii\) a physics\-informed kernel \(KernelOnet\-PIKF\), whose kernel is the analytic fundamental solution of the governing equation, so that the expansion automatically satisfies the governing equation and can therefore be trained without supervision using boundary conditions alone, without any interior solution data; \(iii\) a hybrid kernel \(KernelOnet\-HK\), which splits the solution according to the linear principal part of the governing equation into a homogeneous part spanned by analytic fundamental solutions and a source part spanned byKcK\_\{c\}low\-rank learned correction kernels; the former satisfies the linear principal part exactly while the latter carries the nonlinear source terms that the kernel functions cannot represent, so that when the problem admits no complete fundamental solution, such as a nonlinear equation with nonlinear source terms, a compromise between physical prior and breadth of applicability is still achievable with the help of supervised data\. Numerical experiments on four examples—the Laplace equation on a circular domain, the nonlinear modified Helmholtz equation on a star\-shaped domain, the complex Helmholtz equation in an unbounded exterior domain, and the underwater acoustic radiation and propagation induced by spherical\-shell vibration in a shallow\-water waveguide—show that the proposed KernelOnet achieves high prediction accuracy while significantly improving interpretability, that the unsupervised configuration requires no interior solution labels at all, and that the three variants exhibit complementary advantages on different problems\. It should be noted that both the physics\-informed kernel and the hybrid kernel presuppose the existence of an analytic fundamental solution of the governing equation or of its linear principal part, which is the condition under which the framework takes effect\.

The main contributions of this paper are summarized as follows: \(1\) the KernelOnet framework is proposed, which embeds kernel functions explicitly into a neural operator to replace the implicit trunk network of DeepONet, establishing an intrinsic connection between neural operators and boundary\-type numerical methods based on kernel expansions; \(2\) three complementary kernel construction strategies are given, namely the data\-driven learnable kernel, the physics\-informed kernel and the hybrid kernel, and their physical meaning and parameter efficiency are analysed; \(3\) on three benchmark examples, the advantages of KernelOnet in accuracy, interpretability and unsupervised training are systematically verified, and in particular an effective solution route is provided for unbounded\-domain problems that general\-purpose neural operators struggle to solve; \(4\) with a view to engineering applications, KernelOnet is applied to the underwater acoustic radiation and propagation problem induced by spherical\-shell vibration in a shallow\-water waveguide, where the waveguide Pekeris kernel and the normal\-mode kernel are used to construct the kernel functions, so that the near field and the far field are both solved without supervision and good robustness is exhibited under various sound speed profiles\.

The remainder of this paper is organized as follows\. Section 2 introduces the methodology and theoretical derivation of KernelOnet in detail: it first gives the problem definition of operator learning and its differences from conventional neural networks, then presents the two baseline architectures DeepONet and PI\-DeepONet together with their loss constructions, and afterwards introduces kernel functions and kernel\-expansion\-based boundary\-type solution methods as the theoretical bridge connecting neural operators with PDE solution methods; on this basis, three complementary kernel construction strategies are proposed and compared, and finally the convergence and accuracy of each variant are analysed from the perspective of kernel theory\. Section 3 systematically verifies the proposed method on three benchmark examples and one practical engineering example\. Section 4 concludes the paper and outlines future work\.

## 2Methodology

### 2\.1Operator Learning

Let𝒜\\mathcal\{A\}and𝒰\\mathcal\{U\}denote the input function space and the output solution function space, respectively\. The goal of operator learning is to approximate the solution operator determined by a PDE,

𝒢:𝒜→𝒰,a↦u,\\mathcal\{G\}:\\mathcal\{A\}\\rightarrow\\mathcal\{U\},\\qquad a\\mapsto u,\(1\)where the input functionaa\(for example a boundary condition, a source term or an initial condition\) and the output solutionuusatisfy the governing equationℒ​u=f\\mathcal\{L\}u=ftogether with the corresponding well\-posedness conditions \(boundary and initial conditions\)\. Unlike the traditional "solve problem by problem" paradigm, operator learning simultaneously exploits multiple input–output function pairs

𝒟=\{\(a\(i\),u\(i\)\)\}i=1N\\mathcal\{D\}=\\left\\\{\\left\(a^\{\(i\)\},u^\{\(i\)\}\\right\)\\right\\\}\_\{i=1\}^\{N\}\(2\)to learn a parametric approximation𝒢θ≈𝒢\\mathcal\{G\}\_\{\\theta\}\\approx\\mathcal\{G\}of the operator𝒢\\mathcal\{G\}, so that during inference it directly yields the solutionu=𝒢θ​\(a\)u=\\mathcal\{G\}\_\{\\theta\}\(a\)for an unseen input functionaawithout solving the PDE again\. Representative operator learning methods include DeepONet[Lu et al\. \(2021a\)](https://arxiv.org/html/2609.35938#bib.bib1)and the Fourier neural operator \(FNO\)[Li et al\. \(2020a\)](https://arxiv.org/html/2609.35938#bib.bib2), among others; a review can be found in[Kovachki et al\. \(2021\)](https://arxiv.org/html/2609.35938#bib.bib8)\. Compared with traditional neural network methods, operator learning differs in two essential respects\. First, conventional neural networks learn a mapping between fixed\-dimensional vectors, and both the input and output dimensions are fixed at training time, so that changing the grid or the resolution requires retraining; in operator learning, the inputs and outputs are functions, and the discrete values at a finite number of nodes are merely sampled realizations of those functions, so that in principle the dependence on a fixed discrete layout can be removed\. Whether genuine resolution independence is achieved, however, depends on the specific network architecture: the Fourier neural operator[Li et al\. \(2020a\)](https://arxiv.org/html/2609.35938#bib.bib2)and related methods parameterize functions in the frequency domain and can be evaluated directly on inputs and outputs of different resolutions, and therefore possess resolution invariance; the trunk network of DeepONet can likewise be evaluated at arbitrary output coordinates, but its branch network takes its input from a fixed set of sensors, so that changing the input sampling layout requires interpolation or readjustment of the model\. To address this limitation of the branch network, the resolution\-independent neural operator proposed by Bahmani et al\.[Bahmani et al\. \(2025\)](https://arxiv.org/html/2609.35938#bib.bib5)learns a dictionary of continuous basis functions by means of an implicit neural representation and projects input functions with differing sensor layouts onto a fixed\-dimensional embedding space, thereby obtaining input\-side resolution independence without modifying the DeepONet architecture\. Second, operator learning provides a theoretical guarantee at the level of universal approximation: under certain conditions, operators constructed by neural networks can approximate nonlinear continuous operators to arbitrary accuracy[Chen and Chen \(1995\)](https://arxiv.org/html/2609.35938#bib.bib48), which justifies the feasibility of approximating solution operators with a finite number of parameters\.

### 2\.2DeepONet and PI\-DeepONet

DeepONet builds on the classical results of Chen and Chen[Chen and Chen \(1995\)](https://arxiv.org/html/2609.35938#bib.bib48), Hornik et al\.[Hornik et al\. \(1989\)](https://arxiv.org/html/2609.35938#bib.bib49)and Cybenko[Cybenko \(1989\)](https://arxiv.org/html/2609.35938#bib.bib50)on the universal approximation of neural networks: for any continuous operator𝒢\\mathcal\{G\}and any compact setKK, there exists a finite\-term expansion of the form

𝒢⁡\(a\)​\(𝐱\)≈∑k=1pβk​\(a\)​τk​\(𝐱\)\\mathcal\{G\}\(a\)\(\\mathbf\{x\}\)\\approx\\sum\_\{k=1\}^\{p\}\\beta\_\{k\}\(a\)\\,\\tau\_\{k\}\(\\mathbf\{x\}\)\(3\)that approximates𝒢​\(a\)​\(𝐱\)\\mathcal\{G\}\(a\)\(\\mathbf\{x\}\)uniformly onKK, where the coefficientsβk​\(a\)\\beta\_\{k\}\(a\)depend only on the input functionaaand the basis functionsτk​\(𝐱\)\\tau\_\{k\}\(\\mathbf\{x\}\)depend only on the spatial coordinate𝐱\\mathbf\{x\}\.

DeepONet uses two neural networks to approximate the coefficients and the basis functions of this expansion separately: the branch networkℬ:ℝnb→ℝp\\mathcal\{B\}:\\mathbb\{R\}^\{n\_\{b\}\}\\rightarrow\\mathbb\{R\}^\{p\}maps the discrete values of the input function at a set of sensors to the coefficientsβk\\beta\_\{k\}, and the trunk network𝒯:ℝd→ℝp\\mathcal\{T\}:\\mathbb\{R\}^\{d\}\\rightarrow\\mathbb\{R\}^\{p\}maps the spatial coordinates to the basis functionsτk\\tau\_\{k\}; the final output is the inner product of the two plus a bias:

𝒢θ​\(a\)​\(𝐱\)=∑k=1pbk​\(a\)​tk​\(𝐱\)\+b0,\\mathcal\{G\}\_\{\\theta\}\(a\)\(\\mathbf\{x\}\)=\\sum\_\{k=1\}^\{p\}b\_\{k\}\(a\)\\,t\_\{k\}\(\\mathbf\{x\}\)\+b\_\{0\},\(4\)whereb0b\_\{0\}is a learnable bias term introduced in this paper andppis the number of basis functions \(that is, the output dimension of the trunk network\)\. In practice the trunk output is usually passed through atanh\\tanhactivation to ensure boundedness\. Although DeepONet possesses universal approximation capability, in the fully connected trunk network adopted in this paper the basis functionstk​\(𝐱\)t\_\{k\}\(\\mathbf\{x\}\)are learned entirely implicitly from data, so that neither the physical structure of the governing equation is explicitly exploited nor is a direct physical interpretation available; this is the starting point of the improvement proposed in this paper\. The architecture of DeepONet is shown in Figure[1](https://arxiv.org/html/2609.35938#S2.F1)\.

DeepONet can be trained in a supervised manner; in this paper, the discrete values of the input function at a finite number of boundary points serve as the input of the branch network, and the discrete values of the output solution at interior evaluation points serve as the supervision signal\. Let the training data be

𝒟=\{\(𝐲b\(i\),𝐮\(i\)\)\}i=1N,𝐲b\(i\)∈ℝnb,𝐮\(i\)∈ℝnt,\\mathcal\{D\}=\\left\\\{\\left\(\\mathbf\{y\}\_\{b\}^\{\(i\)\},\\mathbf\{u\}^\{\(i\)\}\\right\)\\right\\\}\_\{i=1\}^\{N\},\\qquad\\mathbf\{y\}\_\{b\}^\{\(i\)\}\\in\\mathbb\{R\}^\{n\_\{b\}\},\\quad\\mathbf\{u\}^\{\(i\)\}\\in\\mathbb\{R\}^\{n\_\{t\}\},\(5\)wherenbn\_\{b\}is the number of boundary points andntn\_\{t\}is the number of interior evaluation points\. The model parametersθ\\thetaare obtained by minimizing the data loss, that is, the pointwise squared error between the output and the reference solution over all interior evaluation points:

𝒥⁡\(θ\)=𝒥data​\(θ\)=1N​nt​∑i=1N∑k=1nt\(𝒢θ​\(𝐲b\(i\)\)​\(𝐱\(k\)\)−uk\(i\)\)2,\\mathcal\{J\}\(\\theta\)=\\mathcal\{J\}\_\{\\mathrm\{data\}\}\(\\theta\)=\\frac\{1\}\{N\\,n\_\{t\}\}\\sum\_\{i=1\}^\{N\}\\sum\_\{k=1\}^\{n\_\{t\}\}\\left\(\\mathcal\{G\}\_\{\\theta\}\\left\(\\mathbf\{y\}\_\{b\}^\{\(i\)\}\\right\)\\left\(\\mathbf\{x\}^\{\(k\)\}\\right\)\-u\_\{k\}^\{\(i\)\}\\right\)^\{2\},\(6\)where𝐱\(k\)\\mathbf\{x\}^\{\(k\)\}is thekk\-th interior evaluation point and𝒢θ​\(𝐲b\(i\)\)​\(𝐱\(k\)\)\\mathcal\{G\}\_\{\\theta\}\\left\(\\mathbf\{y\}\_\{b\}^\{\(i\)\}\\right\)\\left\(\\mathbf\{x\}^\{\(k\)\}\\right\)is the value of the operator output at that point\.

To introduce physical information into operator learning, the physics\-informed DeepONet \(PI\-DeepONet\)[Wang et al\. \(2021\)](https://arxiv.org/html/2609.35938#bib.bib11)was proposed on the basis of DeepONet: instead of relying solely on the purely data\-driven loss above, it adds the residual of the governing equation to the loss function as a soft constraint\. Specifically, in the unsupervised setting where no interior solution labels are used, its loss function consists of two weighted terms:

𝒥⁡\(θ\)=λPDE​𝒥PDE​\(θ\)\+λBC​𝒥BC​\(θ\),\\mathcal\{J\}\(\\theta\)=\\lambda\_\{\\mathrm\{PDE\}\}\\,\\mathcal\{J\}\_\{\\mathrm\{PDE\}\}\(\\theta\)\+\\lambda\_\{\\mathrm\{BC\}\}\\,\\mathcal\{J\}\_\{\\mathrm\{BC\}\}\(\\theta\),\(7\)where𝒥PDE\\mathcal\{J\}\_\{\\mathrm\{PDE\}\}is the PDE residual loss evaluated at interior collocation points,𝒥BC\\mathcal\{J\}\_\{\\mathrm\{BC\}\}is the boundary condition residual loss, andλPDE\\lambda\_\{\\mathrm\{PDE\}\}andλBC\\lambda\_\{\\mathrm\{BC\}\}are the corresponding weight coefficients\. When supervised data are available, a data loss term𝒥data\\mathcal\{J\}\_\{\\mathrm\{data\}\}, that is, the mean squared error loss of DeepONet, can be added to the above expression to further improve accuracy and training stability\. The PDE residual loss is obtained by differentiating the network output to high order, making use of automatic differentiation; for the Laplace equation, for instance,𝒥PDE=1N​∑i∑k\(∇2𝒢θ​\(a\(i\)\)​\(𝐱k\)\)2\\mathcal\{J\}\_\{\\mathrm\{PDE\}\}=\\frac\{1\}\{N\}\\sum\_\{i\}\\sum\_\{k\}\\left\(\\nabla^\{2\}\\mathcal\{G\}\_\{\\theta\}\(a^\{\(i\)\}\)\(\\mathbf\{x\}\_\{k\}\)\\right\)^\{2\}\. This allows the network to be trained with the help of the PDE constraint when interior solution labels are scarce, thereby alleviating to some extent the dependence on large amounts of labelled data\.

However, PI\-DeepONet incorporates physical information into the loss function as a soft constraint, which entails several inherent limitations\. First, the weightsλ\\lambdabetween the PDE residual, the boundary residual and the optional data term must be tuned manually; the optimal weights differ substantially from problem to problem, and imbalanced weights lead to unstable training\. Second, the PDE residual loss requires repeated computation of high\-order derivatives of the network output, which is computationally expensive\. Third, even when training converges, the network output only approximately satisfies the governing equation at the collocation points, so that physical consistency is not rigorously guaranteed\. These limitations show that the soft\-constraint approach can hardly reconcile training stability with rigorous physical consistency\.

Figure 1:Architectures of DeepONet and PI\-DeepONet\. DeepONet learns the basis functionstkt\_\{k\}implicitly through a fully connected trunk network and optimizes the data loss in a supervised manner; on this basis, PI\-DeepONet adds the PDE residual and the boundary condition residual to the loss function as soft constraints, so that in the unsupervised case the loss contains only these two terms, balanced by the weightsλPDE\\lambda\_\{\\mathrm\{PDE\}\}andλBC\\lambda\_\{\\mathrm\{BC\}\}, while the data loss term can be added when supervised data are available to further improve performance\.
### 2\.3Kernel Functions and Kernel\-Expansion Methods for PDEs

A kernel function is a bivariate functionK⁡\(𝐱,𝐲\)K\(\\mathbf\{x\},\\mathbf\{y\}\)defined on a product space\. In function approximation, given a set of centres\{𝐱s\(j\)\}j=1Nc\\left\\\{\\mathbf\{x\}\_\{s\}^\{\(j\)\}\\right\\\}\_\{j=1\}^\{N\_\{c\}\}, a kernel function can be used directly as a basis function to expand a target functionu⁡\(𝐱\)u\(\\mathbf\{x\}\)as a linear combination of kernels,

u⁡\(𝐱\)≈∑j=1Ncβj​K​\(𝐱,𝐱s\(j\)\),u\(\\mathbf\{x\}\)\\approx\\sum\_\{j=1\}^\{N\_\{c\}\}\\beta\_\{j\}\\,K\\left\(\\mathbf\{x\},\\mathbf\{x\}\_\{s\}^\{\(j\)\}\\right\),\(8\)where the coefficientsβj\\beta\_\{j\}are determined by data or by the equation\.

As noted above, kernel methods have been widely used in machine learning, with a theoretical foundation given by the theory of reproducing kernel Hilbert spaces \(RKHS\)[Aronszajn \(1950\)](https://arxiv.org/html/2609.35938#bib.bib19): when a kernel function is symmetric and positive definite, it induces a complete inner product space and becomes the reproducing kernel of that space; the relevant properties will be discussed in Section 2\.5\.

This section focuses on another use of kernel functions: as basis functions for the solutions of partial differential equations\. Since the solution of a PDE is itself a function, if its approximation space is spanned by kernel functions, solving the PDE is reduced to determining the coefficients of the kernel expansion; when the kernel functions and the centres are suitably configured, this expansion can approximate the target function in the corresponding function space\. The specific form of the kernel function determines the approximation capability of the approximation space, so that choosing a suitable kernel for a given problem is essential\. In this paper, general kernel functions are denoted byφ\\varphiand physics\-informed kernels such as fundamental solutions byΦ\\Phi\. The remainder of this section first presents these two classes of kernel functions \(Section 2\.3\.1\) and then summarizes the kernel\-expansion solution methods derived from them \(Section 2\.3\.2\)\.

#### 2\.3\.1General and Physics\-Informed Kernels

Depending on whether they contain physical information of the governing equation, kernel functions can be divided into two broad classes: general kernel functions and physics\-informed kernel functions\. A general kernel function does not depend on a specific equation but is determined only by geometric and smoothness requirements; the most typical example is the radial basis function \(RBF\)[Buhmann \(2003\)](https://arxiv.org/html/2609.35938#bib.bib23);[Wendland \(2005\)](https://arxiv.org/html/2609.35938#bib.bib24)\. An RBF depends only on the radial distance between two points, i\.e\.K⁡\(𝐱,𝐲\)=φ⁡\(\|𝐱−𝐲\|\)K\(\\mathbf\{x\},\\mathbf\{y\}\)=\\varphi\\left\(\\left\|\\mathbf\{x\}\-\\mathbf\{y\}\\right\|\\right\), and its rotational invariance makes it naturally suited to isotropic problems\. Common RBFs include the Gaussian, multiquadric and thin\-plate spline kernels, which are positive definite or conditionally positive definite and possess universal approximation capability[Micchelli \(1986\)](https://arxiv.org/html/2609.35938#bib.bib20);[Park and Sandberg \(1991\)](https://arxiv.org/html/2609.35938#bib.bib21);[Girosi et al\. \(1995\)](https://arxiv.org/html/2609.35938#bib.bib28); precisely because they do not satisfy any governing equation in advance, their form can be chosen with great flexibility, which makes them suitable for complex geometries and nonlinear problems\.

In contrast to general kernel functions, a physics\-informed kernel function \(PIKF\) is a kernel function that contains, wholly or in part, information about the governing equation \(PDE\); it may be chosen as a fundamental solution, a Green’s function, a harmonic function, a radial Trefftz function or even the solution of some linearly simplified PDE[Fu et al\. \(2024\)](https://arxiv.org/html/2609.35938#bib.bib52);[Kita and Kamiya \(1995\)](https://arxiv.org/html/2609.35938#bib.bib41)\. The most basic of these is the fundamental solution of the governing equation\. Letℒ\\mathcal\{L\}be a linear partial differential operator\. If there exists a functionΦ\\Phisatisfying

ℒ​Φ=−δ⁡\(𝐱−𝐱s\),\\mathcal\{L\}\\Phi=\-\\delta\(\\mathbf\{x\}\-\\mathbf\{x\}\_\{s\}\),\(9\)whereδ\\deltais the Dirac distribution, and𝐱s\\mathbf\{x\}\_\{s\}is the source point, thenΦ\\Phiis called the fundamental solution of the operatorℒ\\mathcal\{L\}at𝐱s\\mathbf\{x\}\_\{s\}[Brebbia et al\. \(1984\)](https://arxiv.org/html/2609.35938#bib.bib42); unlike a Green’s function, which further satisfies specific homogeneous boundary conditions, a fundamental solution carries no boundary conditions whatsoever\. Taking the two\-dimensional Laplace operatorℒ=∇2\\mathcal\{L\}=\\nabla^\{2\}as an example, for a radial function the Laplacian in polar coordinates reduces to∇2=∂2∂r2\+1r​∂∂r\\nabla^\{2\}=\\dfrac\{\\partial^\{2\}\}\{\\partial r^\{2\}\}\+\\dfrac\{1\}\{r\}\\dfrac\{\\partial\}\{\\partial r\}; for a rotationally symmetric functionΦ⁡\(r\)\\Phi\(r\)that depends only on the radial distance, the equation∇2Φ=−δ⁡\(𝐱−𝐱s\)\\nabla^\{2\}\\Phi=\-\\delta\(\\mathbf\{x\}\-\\mathbf\{x\}\_\{s\}\)reduces, forr≠0r\\neq 0, to the ordinary differential equation

Φ′′​\(r\)\+1r​Φ′​\(r\)=0⟹Φ⁡\(r\)=C1​ln⁡r\+C2\.\\Phi^\{\\prime\\prime\}\(r\)\+\\frac\{1\}\{r\}\\Phi^\{\\prime\}\(r\)=0\\quad\\Longrightarrow\\quad\\Phi\(r\)=C\_\{1\}\\ln r\+C\_\{2\}\.\(10\)The constants of integration are determined by the flux condition at the source point, which gives the fundamental solution of the two\-dimensional Laplace operator,

Φ⁡\(r\)=−12​π​ln⁡r,r=\|𝐱−𝐱s\|,\\Phi\(r\)=\-\\frac\{1\}\{2\\pi\}\\ln r,\\qquad r=\\left\|\\mathbf\{x\}\-\\mathbf\{x\}\_\{s\}\\right\|,\(11\)satisfying∇2Φ=−δ⁡\(𝐱−𝐱s\)\\nabla^\{2\}\\Phi=\-\\delta\(\\mathbf\{x\}\-\\mathbf\{x\}\_\{s\}\)\. Similarly, the fundamental solution of the Helmholtz operatorℒ=∇2\+k2\\mathcal\{L\}=\\nabla^\{2\}\+k^\{2\}is given by the zeroth\-order Hankel function of the first kind:

Φ⁡\(r\)=i4​H0\(1\)​\(k​r\),\\Phi\(r\)=\\frac\{i\}\{4\}H\_\{0\}^\{\(1\)\}\(kr\),\(12\)which also automatically satisfies the Sommerfeld radiation condition at infinity and is therefore particularly suitable for unbounded exterior domain problems\. Unlike a general kernel function, a PIKF embeds the physical information of the governing equation explicitly into the kernel function, so that the kernel automatically satisfies, or partially satisfies, the governing equation and thus naturally preserves physical consistency in function approximation\.

#### 2\.3\.2Kernel\-Expansion Solution Methods for PDEs

The two classes of kernel functions above correspond precisely to two kernel\-expansion routes: general kernel functions to the domain\-type route and physics\-informed kernel functions to the boundary\-type route\. What these methods have in common is that the solution is written as a linear combination of kernel functions at a number of centres \(source points\); they differ in three respects—whether the kernel function satisfies the governing equation in advance, whether the centres lie in the whole domain or on the boundary, and whether the boundary conditions are enforced pointwise in a strong form or in a weak form through a boundary integral equation\. Along these three dimensions, three basic classes of methods can be distinguished: \(i\) domain\-type collocation, which uses general kernel functions as basis functions with centres distributed throughout the domain; \(ii\) boundary\-type collocation, which uses physics\-informed kernel functions as basis functions with centres \(source points\) placed on the boundary or on a virtual boundary close to it; and \(iii\) boundary integral methods, which use fundamental solutions as integral kernels to construct single\- or double\-layer potentials on the boundary and replace the collocation coefficients by boundary densities\. None of these three classes requires a volume mesh; in particular, \(ii\) and \(iii\) require only a boundary discretization, thereby alleviating the limitations imposed by mesh generation on complex geometries and large\-scale problems\. Beyond this framework, when the equation contains nonlinear terms and the operator no longer admits a usable fundamental solution, one may take the fundamental solution of the linear principal part as the kernel and supplement it with linearization or correction kernels, which constitutes a fourth class of extension \(see item \(4\)\)\.

##### \(1\) Domain\-type collocation

Domain\-type collocation places collocation points both in the interior of the solution domain and on its boundary, the classical representative being Kansa’s method[Kansa \(1990\)](https://arxiv.org/html/2609.35938#bib.bib31)\. This method uses general radial basis functionsφ\\varphias basis functions and expands the solution as a linear combination overNNcentres\{𝐱j\}j=1N\\\{\\mathbf\{x\}\_\{j\}\\\}\_\{j=1\}^\{N\}distributed throughout the domain,

u⁡\(𝐱\)≈∑j=1Nαj​φ​\(\|𝐱−𝐱j\|\),u\(\\mathbf\{x\}\)\\approx\\sum\_\{j=1\}^\{N\}\\alpha\_\{j\}\\,\\varphi\\left\(\\left\|\\mathbf\{x\}\-\\mathbf\{x\}\_\{j\}\\right\|\\right\),\(13\)and then requires the interior collocation points to satisfy the governing equationℒ​u=f\\mathcal\{L\}u=fand the boundary collocation points to satisfy the boundary conditions; the linear system formed by these collocation equations determines the unknown coefficientsαj\\alpha\_\{j\}\. Since RBFs need not satisfy the governing equation in advance, domain\-type collocation can handle complex geometries and nonlinear problems flexibly; the price is that the collocation points must cover the entire domain and that the resulting algebraic system is usually dense, so that the computational cost rises markedly as the number of degrees of freedom increases\. Moreover, the collocation matrix of Kansa’s method is in general non\-symmetric and may be singular or severely ill\-conditioned, its well\-posedness depending on the collocation layout and on the choice of the shape parameter of the basis functions[Fasshauer \(2007\)](https://arxiv.org/html/2609.35938#bib.bib34)\.

##### \(2\) Boundary\-type collocation

Unlike the domain\-type methods, a boundary\-type collocation method places collocation points only on the boundary of the solution domain and selects physics\-informed kernel functions that automatically satisfy the governing equation as its basis functions\. Since the governing equation is already satisfied exactly inside the domain by the basis functions, it suffices to enforce the boundary conditions at the boundary collocation points, and this feature reduces the dimension of the problem by one—a three\-dimensional problem becoming two\-dimensional, for instance—so that the number of collocation points required is greatly reduced\. Its classical representative is the method of fundamental solutions \(MFS\)[Golberg and Chen \(1999\)](https://arxiv.org/html/2609.35938#bib.bib40);[Fairweather and Karageorghis \(1998\)](https://arxiv.org/html/2609.35938#bib.bib38), which uses the fundamental solution of the governing equation as the basis function and enjoys advantages such as exponential convergence and high accuracy[Golberg \(1995\)](https://arxiv.org/html/2609.35938#bib.bib36);[Chen et al\. \(1998\)](https://arxiv.org/html/2609.35938#bib.bib37);[Fairweather et al\. \(2003\)](https://arxiv.org/html/2609.35938#bib.bib39)\. To avoid the singularity of the fundamental solution at the source point, the MFS places the source points on a virtual boundary outside the physical domain\. Let\{𝐱s\(j\)\}j=1Ns\\left\\\{\\mathbf\{x\}\_\{s\}^\{\(j\)\}\\right\\\}\_\{j=1\}^\{N\_\{s\}\}beNsN\_\{s\}source points distributed outside the solution domain; then the solution can be approximated as

u⁡\(𝐱\)≈∑j=1Nsαj​Φ​\(\|𝐱−𝐱s\(j\)\|\),u\(\\mathbf\{x\}\)\\approx\\sum\_\{j=1\}^\{N\_\{s\}\}\\alpha\_\{j\}\\,\\Phi\\left\(\\left\|\\mathbf\{x\}\-\\mathbf\{x\}\_\{s\}^\{\(j\)\}\\right\|\\right\),\(14\)where the coefficientsαj\\alpha\_\{j\}are determined by the boundary conditions; by the denseness of the solution space of elliptic equations \(a Runge\-type approximation theorem\)[Fairweather and Karageorghis \(1998\)](https://arxiv.org/html/2609.35938#bib.bib38), this kernel expansion approaches the true solution as the number of source pointsNsN\_\{s\}increases, when the source points are suitably configured, exhibiting spectral convergence for smooth problems\.

A key difficulty of boundary\-type collocation is the source\-point layout: the source points must avoid the boundary in order to circumvent the singularity of the fundamental solution, yet their positions directly affect the condition number of the coefficient matrix and the accuracy\. The MFS avoids the singularity by moving the source points outward onto a virtual boundary outside the solution domain, but the position of this virtual boundary is usually chosen by experience\. The singular boundary method \(SBM\)[Gu et al\. \(2011\)](https://arxiv.org/html/2609.35938#bib.bib46);[Fu et al\. \(2020\)](https://arxiv.org/html/2609.35938#bib.bib53)instead places the source points directly on the physical boundary nodes and subtracts the singular part of the fundamental solution by means of source intensity factors, so that no virtual boundary is needed; the price is that a source intensity factor has to be estimated for every boundary node\. Both classes of methods show that, although an expansion based on fundamental solutions is physically consistent and highly accurate, the distance between the source points and the boundary is a hyperparameter that must be handled with care\.

This sensitivity to the hyperparameter stems from the analyticity of the fundamental\-solution kernel: the source points of the MFS lie outside the domain, so that its kernel is analytic and smooth on the boundary and the corresponding boundary operator is a smoothing \(infinitely differentiable\) operator, whose singular values decay geometrically and whose condition number grows exponentially with the number of source pointsNsN\_\{s\}\. In other words, the MFS trades spectral accuracy for severe ill\-conditioning, and the choice of the source distance is essentially a compromise between approximation accuracy and matrix conditioning\.

##### \(3\) Boundary integral methods

Besides serving as basis functions for collocation, fundamental solutions can also act as the integral kernel of boundary integral methods, and are used to construct boundary integral equations \(BIEs\)\. The starting point is to represent the solution inside the domain as a potential integral over the boundary; taking the single\-layer potential constructed from a general fundamental solutionΦ\\Phias an example,

u⁡\(𝐱\)=∫ΓΦ⁡\(\|𝐱−𝐲\|\)​σ​\(𝐲\)​d​Γ𝐲,u\(\\mathbf\{x\}\)=\\int\_\{\\Gamma\}\\Phi\\left\(\\left\|\\mathbf\{x\}\-\\mathbf\{y\}\\right\|\\right\)\\sigma\(\\mathbf\{y\}\)\\,\\mathrm\{d\}\\Gamma\_\{\\mathbf\{y\}\},\(15\)whereσ\\sigmais the unknown boundary density\. After discretizing the boundary into elements and approximating the integral numerically, the above expression can formally be written as a linear combination of basis functions,

u⁡\(𝐱\)≈∑j=1Neσj​∫ΓjΦ⁡\(\|𝐱−𝐲\|\)​d​Γ𝐲,u\(\\mathbf\{x\}\)\\approx\\sum\_\{j=1\}^\{N\_\{e\}\}\\sigma\_\{j\}\\int\_\{\\Gamma\_\{j\}\}\\Phi\\left\(\\left\|\\mathbf\{x\}\-\\mathbf\{y\}\\right\|\\right\)\\mathrm\{d\}\\Gamma\_\{\\mathbf\{y\}\},\(16\)whereσj\\sigma\_\{j\}is the constant density on thejj\-th element and the basis functions are the integrals of the fundamental solution over that element \(the self\-influence element contains an integrable weak singularity\)\. This structure is precisely the foundation of the boundary element method \(BEM\)[Brebbia et al\. \(1984\)](https://arxiv.org/html/2609.35938#bib.bib42): it shares the same fundamental\-solution kernel with the MFS, the difference being that a weak\-form integral equation replaces strong\-form pointwise collocation\. Because the integration is performed on the true boundary, the BEM needs to introduce no virtual source distance; although the self\-influence element integral contains a weak \(logarithmic\) singularity, this singularity is integrable and can be handled exactly by analytic integration or by singularity subtraction\.

From the point of view of operator properties, the kernel of the BEM is only weakly singular on the boundary and the operator smooths only to first order, so that its singular values mostly decay algebraically and its condition number grows polynomially with the number of elementsNeN\_\{e\}—usually milder than the exponential ill\-conditioning of the MFS—but near corners the boundary density becomes singular and graded meshes or special elements must be used\. For comparison, the virtual boundary method \(VBM\)[Sun and Yao \(1997\)](https://arxiv.org/html/2609.35938#bib.bib47)places the boundary integral on a virtual boundary enclosing the solution domain, so that no singular integral has to be handled because the integration points are far from the true boundary; its price is the same as that of the MFS—the position of the virtual boundary must be chosen artificially and directly affects the computational accuracy and the condition number of the coefficient matrix\.

##### \(4\) Extension to nonlinear problems

When the governing equation contains nonlinear terms, the principle of superposition no longer holds and the operator has no usable fundamental solution in the usual sense, so that the kernel expansion formed by a linear superposition of fundamental solutions no longer applies\. The typical way to handle such problems is to "approximate the nonlinear by the linear": the nonlinear operator is split into a linear principal partℒ0\\mathcal\{L\}\_\{0\}and a nonlinear remainder, or the nonlinear term is linearized by a Newton or Picard iteration at the current approximate solution, reducing the original problem to a sequence of linear subproblems, each of which is then solved by the kernel expansion described above[Brebbia et al\. \(1984\)](https://arxiv.org/html/2609.35938#bib.bib42);[Wrobel and Brebbia \(1987\)](https://arxiv.org/html/2609.35938#bib.bib44)\.

After linearization, the nonlinear remainder is moved to the right\-hand side of the equation and acts as an equivalent source term, so that each step still requires the solution of a linear equation of the formℒ0​u=g\\mathcal\{L\}\_\{0\}u=g\. Its solution can be written as a homogeneous part spanned by the fundamental solution of the linear principal part plus a particular\-solution kernel for the equivalent source term,

u⁡\(𝐱\)≈∑j=1Nsαj​Φ0​\(\|𝐱−𝐱j\|\)\+∑k=1Kcγk​Ψ0​\(\|𝐱−𝐱k\|\),u\(\\mathbf\{x\}\)\\approx\\sum\_\{j=1\}^\{N\_\{s\}\}\\alpha\_\{j\}\\,\\Phi\_\{0\}\\left\(\\left\|\\mathbf\{x\}\-\\mathbf\{x\}\_\{j\}\\right\|\\right\)\+\\sum\_\{k=1\}^\{K\_\{c\}\}\\gamma\_\{k\}\\,\\Psi\_\{0\}\\left\(\\left\|\\mathbf\{x\}\-\\mathbf\{x\}\_\{k\}\\right\|\\right\),\(17\)whereΦ0\\Phi\_\{0\}is the fundamental solution of the linear principal\-part operatorℒ0\\mathcal\{L\}\_\{0\}andΨ0\\Psi\_\{0\}is the particular\-solution kernel of the equivalent source term, whileαj\\alpha\_\{j\}andγk\\gamma\_\{k\}are the expansion coefficients, the latter being determined by a functional approximation of the equivalent source term\. In the Newton or Picard iteration it suffices to substitute the current residual into the equivalent source term at each step, so that the same set of kernel expansions can be reused: the homogeneous part satisfiesℒ0​Φ0=δ\\mathcal\{L\}\_\{0\}\\Phi\_\{0\}=\\deltaexactly, while all the nonlinear information is carried by the particular\-solution kernel and updated as the iteration proceeds\. This construction—"fundamental solution of the linear principal part plus nonlinear correction", of which the dual reciprocity method is a classical realization[Nardini and Brebbia \(1983\)](https://arxiv.org/html/2609.35938#bib.bib43);[Wrobel and Brebbia \(1987\)](https://arxiv.org/html/2609.35938#bib.bib44);[Partridge et al\. \(1991\)](https://arxiv.org/html/2609.35938#bib.bib45)—gives kernel\-expansion methods the ability to handle nonlinear problems while preserving physical consistency, and it is also the common starting point of all kinds of nonlinear kernel methods, whether domain\-type collocation or boundary\-type expansion\.

In summary, kernel functions occupy a central position in kernel\-expansion\-based PDE solution methods: domain\-type collocation uses general radial basis functions as kernels and is flexible in form but requires a domain\-wide discretization; boundary\-type collocation and boundary integral methods take the fundamental solution \(PIKF\) of the governing equation as their core and embed the physical information of "automatically satisfying the governing equation" explicitly into the basis functions, thereby obtaining higher accuracy with fewer degrees of freedom; and for nonlinear problems one takes the fundamental solution of the linear principal part and supplements it with linearization or correction kernels\. It can thus be seen that domain\-type collocation, boundary\-type collocation and boundary integral methods share one and the same library of kernel functions and differ only in the choice of kernel function and in the way the boundary conditions are imposed, and are therefore unified at the level of the kernel expansion\.

### 2\.4Kernel Operator Network

The kernel\-expansion perspective of Section 2\.3 has a direct implication for the design of neural operators: the trunk network of DeepONet[Lu et al\. \(2021a\)](https://arxiv.org/html/2609.35938#bib.bib1)can be regarded as learning basis functions implicitly in a data\-driven sense, whereas the Kernel Operator Network \(KernelOnet\) proposed in this paper makes the kernel functions explicit, so that the network structure and the kernel expansion are mathematically consistent—the coefficients are given by the branch network and the kernel functions by the trunk network\. Based on this observation, the general form of KernelOnet is

𝒢θ​\(a\)​\(𝐱\)=∑j=1Bbj​\(a\)​ψj​\(𝐱\),\\mathcal\{G\}\_\{\\theta\}\(a\)\(\\mathbf\{x\}\)=\\sum\_\{j=1\}^\{B\}b\_\{j\}\(a\)\\,\\psi\_\{j\}\(\\mathbf\{x\}\),\(18\)where the coefficients are given by the branch network𝐛​\(a\)=ℬθ​\(a\)\\mathbf\{b\}\(a\)=\\mathcal\{B\}\_\{\\theta\}\(a\)and the family of basis functionsψj​\(𝐱\)\\psi\_\{j\}\(\\mathbf\{x\}\)is no longer learned implicitly by a free network but is explicitly constrained to be a kernel function\. Figure[2](https://arxiv.org/html/2609.35938#S2.F2)shows the overall architecture of KernelOnet: the three variants share the same branch network and the same operator expansion structure and differ only in the construction of the kernel trunkψj​\(𝐱\)\\psi\_\{j\}\(\\mathbf\{x\}\)—KernelOnet\-RBF parameterizes a radial basis function by a neural network and learns the kernel shape from data; KernelOnet\-PIKF takes the analytic fundamental solution as its kernel, and admits two equivalent forms, collocation and boundary integral \(the former uses aγ\\gamma\-shifted fundamental solution, the latter a single\-layer potential integral over the true boundary\), so that the expansion satisfies the governing equation exactly; KernelOnet\-HK, on the other hand, mixes analytic fundamental\-solution kernels with low\-rank learned correction kernels, absorbing the nonlinear source terms into a correction branch withKc≪BK\_\{c\}\\ll B\. The three variants are introduced in the following three subsections\.

Figure 2:Architecture of KernelOnet\. The three variants share the same branch network and the same operator expansion structureu⁡\(𝐱\)=∑j=1Bbj​\(a\)​ψj​\(𝐱\)u\(\\mathbf\{x\}\)=\\sum\_\{j=1\}^\{B\}b\_\{j\}\(a\)\\,\\psi\_\{j\}\(\\mathbf\{x\}\), differing only in the kernel trunkψj​\(𝐱\)\\psi\_\{j\}\(\\mathbf\{x\}\): KernelOnet\-RBF uses a neural\-network\-parameterized radial basis function; KernelOnet\-PIKF takes the analytic fundamental solution as its kernel \(the collocation form uses aγ\\gamma\-shifted fundamental solution, the boundary integral form a single\-layer potential kernel on the true boundary\) and its expansion satisfies the governing equation exactly, so that it can be trained without supervision; KernelOnet\-HK mixes analytic fundamental\-solution kernels with low\-rank learned correction kernels,ψ=Φ⁡\(γ​𝐱b\)⊕φθ​\(𝐱−𝐭k\)\\psi=\\Phi\(\\gamma\\,\\mathbf\{x\}\_\{b\}\)\\oplus\\varphi\_\{\\theta\}\(\\mathbf\{x\}\-\\mathbf\{t\}\_\{k\}\), absorbing the nonlinear source terms into a correction branch withKc≪BK\_\{c\}\\ll Band requiring supervised training\.#### 2\.4\.1KernelOnet\-RBF: Data\-Driven Learnable Kernel

KernelOnet\-RBF learns the kernel function directly in a data\-driven manner, using a neural\-network\-parameterized radial basis function as the kernel:

ψj​\(𝐱\)=φθ​\(rj\),rj=\|𝐱−𝐱b\(j\)\|,\\psi\_\{j\}\(\\mathbf\{x\}\)=\\varphi\_\{\\theta\}\\left\(r\_\{j\}\\right\),\\qquad r\_\{j\}=\\left\|\\mathbf\{x\}\-\\mathbf\{x\}\_\{b\}^\{\(j\)\}\\right\|,\(19\)whereφθ\\varphi\_\{\\theta\}is a univariate shallow network taking only the radial distancerras input, and is learned jointly with the branch network during training\. The kernel centres are taken directly as the boundary discretization points\{𝐱b\(j\)\}j=1B\\left\\\{\\mathbf\{x\}\_\{b\}^\{\(j\)\}\\right\\\}\_\{j=1\}^\{B\}, so that every output coefficient of the branch network corresponds exactly to the weight of the basis function at one boundary point, with a clear physical meaning\. This kernel function does not assume a specific form in advance but learns the radial shape adaptively from data\. For constant\-coefficient linear equations, by the approximation theory of Section 2\.5 the solution can be spanned by an expansion of fundamental solutions at the boundary centres, and the fundamental solution is the optimal choice for a single radial kernel; in this case the learned kernel function can be regarded as a regularized approximation of the fundamental solution\. It must be emphasized that this interpretation holds only in that case: once the equation contains nonlinear source terms or variable coefficients, the solution can no longer be spanned by an expansion of a single fundamental solution, and the kernel function has to fit both the homogeneous field and the source term, so that it can no longer be interpreted as a fundamental solution\. At the same time, its inductive bias of depending only on the radial distance guarantees rotational invariance, consistent with the symmetry of isotropic fundamental solutions in constant\-coefficient linear problems\. KernelOnet\-RBF is trained with the data loss of Section 2\.2 and therefore belongs to supervised learning; its range of applicability is also the broadest of the three variants, since it can naturally handle strongly nonlinear problems, variable coefficients and complex geometries for which no analytic fundamental solution exists\. The learned kernel function has the potential to transfer across problems and can still offer a certain degree of interpretability in the form of a readable radial curve\. The price is that the physical guarantee of "the expansion automatically satisfying the governing equation" is lost entirely: the form to which the kernel function converges is determined completely by data, and physical consistency is only soft and indirect\.

#### 2\.4\.2KernelOnet\-PIKF: Physics\-Informed Kernel

KernelOnet\-PIKF embeds physical information explicitly into the kernel function\. It rests on the concept of the physics\-informed kernel function \(PIKF\) introduced in Section 2\.3: taking a PIKF as the basis function of the trunk network embeds PDE information directly into the network structure, without the need to add PDE constraints to the loss function as PINNs do; this idea has been validated by the physics\-informed kernel function neural network \(PIKFNN\) proposed by Fu et al\.[Fu et al\. \(2024\)](https://arxiv.org/html/2609.35938#bib.bib52)\. KernelOnet\-PIKF takes the analytic fundamental solution of the governing equation as its kernel and constructs the following collocation expansion:

𝒢θ​\(a\)​\(𝐱\)=∑j=1Bbj​\(a\)​Φ​\(\|𝐱−γ​𝐱b\(j\)\|\),\\mathcal\{G\}\_\{\\theta\}\(a\)\(\\mathbf\{x\}\)=\\sum\_\{j=1\}^\{B\}b\_\{j\}\(a\)\\,\\Phi\\left\(\\left\|\\mathbf\{x\}\-\\gamma\\,\\mathbf\{x\}\_\{b\}^\{\(j\)\}\\right\|\\right\),\(20\)whereΦ\\Phiis the analytic fundamental solution of the governing equation,γ\\gammais a learnable radial scaling scalar andbj​\(a\)b\_\{j\}\(a\)is given by the branch network; the source points are the boundary discretization points𝐱b\(j\)\\mathbf\{x\}\_\{b\}^\{\(j\)\}scaled byγ\\gamma\. SinceΦ\\Phisatisfiesℒ​Φ=−δ⁡\(𝐱−𝐱s\)\\mathcal\{L\}\\Phi=\-\\delta\(\\mathbf\{x\}\-\\mathbf\{x\}\_\{s\}\)\(whereℒ\\mathcal\{L\}is the very governing equation to be approximated\) and the source points lie outside the solution domain, each basis function satisfies the homogeneous governing equation inside the domain, so that for arbitrary values of the coefficients the expansion automatically satisfies the governing equation inside the domain and the PDE residual vanishes identically\. Consequently KernelOnet\-PIKF requires no interior solution data and no computation of a PDE residual term for training; it suffices to enforce the boundary conditions on the boundary, which realizes unsupervised operator learning in the true sense\.

Within the same boundary\-type framework, KernelOnet\-PIKF can also be written as a kernel expansion in boundary integral \(single\-layer potential\) form, that is, with the fundamental\-solution kernel distributed on the true boundary:

𝒢θ​\(a\)​\(𝐱\)=∫ΓΦ⁡\(\|𝐱−𝐲\|\)​b​\(a\)​\(𝐲\)​d​Γ𝐲≈∑j=1Bbj​\(a\)​∫ΓjΦ⁡\(\|𝐱−𝐲\|\)​d​Γ𝐲,\\mathcal\{G\}\_\{\\theta\}\(a\)\(\\mathbf\{x\}\)=\\int\_\{\\Gamma\}\\Phi\\left\(\\left\|\\mathbf\{x\}\-\\mathbf\{y\}\\right\|\\right\)b\(a\)\(\\mathbf\{y\}\)\\,\\mathrm\{d\}\\Gamma\_\{\\mathbf\{y\}\}\\;\\approx\\;\\sum\_\{j=1\}^\{B\}b\_\{j\}\(a\)\\int\_\{\\Gamma\_\{j\}\}\\Phi\\left\(\\left\|\\mathbf\{x\}\-\\mathbf\{y\}\\right\|\\right\)\\mathrm\{d\}\\Gamma\_\{\\mathbf\{y\}\},\(21\)whereb⁡\(a\)b\(a\)is the boundary density function output by the branch network andΓj\\Gamma\_\{j\}is thejj\-th boundary patch\. Both expressions take the fundamental solutionΦ\\Phias their kernel and therefore belong to the same PIKF framework, corresponding respectively to the boundary collocation and the boundary integral implementation\. The two are not the same limit: in the integral form the source distribution lies on the true boundaryΓ\\Gamma, and although the self\-influence element integral contains a weak singularity \(logarithmic in two dimensions\) it is integrable, so that the value at a boundary collocation point is finite and no source point needs to be moved; the collocation form, by contrast, is a point source of zero measure, and when the source point falls on the boundary the kernel is unbounded there so that the collocation equation cannot hold directly, and the source points must be moved off the boundary\. In other words, the point\-source form corresponds to a kernel expansion on a virtual boundary while the integral form corresponds to a kernel expansion on the true boundary\.

Thus KernelOnet\-PIKF structurally unifies the boundary\-type kernel methods of Section 2\.3: its collocation form agrees with the collocation forms of the MFS[Golberg and Chen \(1999\)](https://arxiv.org/html/2609.35938#bib.bib40), the Trefftz method[Kita and Kamiya \(1995\)](https://arxiv.org/html/2609.35938#bib.bib41)and the SBM[Gu et al\. \(2011\)](https://arxiv.org/html/2609.35938#bib.bib46), while its boundary integral form shares the same integral kernel as the single\-layer potential of the BEM[Brebbia et al\. \(1984\)](https://arxiv.org/html/2609.35938#bib.bib42)\. The difference is that the expansion coefficients are no longer obtained by solving simultaneously for the source\-point layout and the boundary conditions but are given directly by the branch network, and that in the collocation form the source points give a virtual surface automatically through the scaling of the boundary points, which removes the need for the virtual boundary selection of the MFS/VBM and for the source intensity factor estimation of the SBM; the boundary integral form, meanwhile, performs panel integration directly on the true boundary\.

The training loss contains only the boundary residual:

𝒥⁡\(θ\)=1N​B​∑i=1N∑k=1B\|𝒢θ​\(𝐲b\(i\)\)​\(𝐱b\(k\)\)−yb,k\(i\)\|2,\\mathcal\{J\}\(\\theta\)=\\frac\{1\}\{N\\,B\}\\sum\_\{i=1\}^\{N\}\\sum\_\{k=1\}^\{B\}\\left\|\\mathcal\{G\}\_\{\\theta\}\\left\(\\mathbf\{y\}\_\{b\}^\{\(i\)\}\\right\)\\left\(\\mathbf\{x\}\_\{b\}^\{\(k\)\}\\right\)\-y\_\{b,k\}^\{\(i\)\}\\right\|^\{2\},\(22\)where𝐲b\(i\)=\(yb,1\(i\),…,yb,B\(i\)\)\\mathbf\{y\}\_\{b\}^\{\(i\)\}=\\left\(y\_\{b,1\}^\{\(i\)\},\\ldots,y\_\{b,B\}^\{\(i\)\}\\right\)is the boundary condition of theii\-th sample, so that the supervision signal is the input itself\. Unlike physics\-informed neural networks \(PINNs\)[Raissi et al\. \(2019\)](https://arxiv.org/html/2609.35938#bib.bib12)and the physics\-informed DeepONet \(PI\-DeepONet\)[Wang et al\. \(2021\)](https://arxiv.org/html/2609.35938#bib.bib11), which add the PDE residual to the loss as a soft constraint and satisfy the governing equation only approximately through optimization, KernelOnet\-PIKF hard\-codes the physical information into the kernel structure, which both avoids the computational cost of repeatedly evaluating high\-order PDE residuals and eliminates the problem of balancing the weights between the PDE residual and the boundary residual\.

The above expansion and training scheme can be extended directly to the complex case: taking the complex Helmholtz equation in an unbounded exterior domain as an example, one chooses the complex fundamental solutionΦ⁡\(r\)=i4​H0\(1\)​\(k​r\)\\Phi\(r\)=\\frac\{i\}\{4\}H\_\{0\}^\{\(1\)\}\(kr\)as the kernel, and both the kernel function and the expansion coefficients take complex values, so that the form of the expansion and of the loss remains unchanged, while the real and imaginary parts of the output solution are given by the cross inner products of the real and imaginary parts of the coefficients with those of the kernel function\. Since the complex fundamental solution itself satisfies the governing equation and the Sommerfeld radiation condition exactly, the above properties continue to hold in the complex case, as will be verified in Example 3\.

Moreover, because the kernel function is taken directly as the analytic fundamental solution, KernelOnet\-PIKF has the strongest interpretability of the three variants, since after training the physics learned by the network can be read off pointwise; the price is a strict precondition: the whole governing equation must admit an analytic fundamental solution and the solution must satisfy the principle of superposition, so that it applies only to linear problems that themselves possess fundamental solutions, such as Laplace, Helmholtz, modified Helmholtz, biharmonic and convection–diffusion problems\. The fundamental solutions of common operators are given in Appendix[A](https://arxiv.org/html/2609.35938#A1)\. Among them, the fundamental solution automatically satisfies the radiation condition at infinity, which makes it particularly suitable for unbounded exterior domains and wave propagation problems; for nonlinear problems whose governing equation admits no fundamental solution, KernelOnet\-PIKF no longer applies and KernelOnet\-HK, introduced in the next subsection, should be used instead\.

When the family of kernel functions is numerically ill\-conditioned under the given discretization of source and evaluation points, so that using it directly as a kernel would pollute the gradients, the kernel basis must be orthogonalized as a preprocessing step\. In this paper, the kernel function matrixG⁡\(𝐱\)G\(\\mathbf\{x\}\)is orthogonalized by a singular value decomposition \(SVD\) of the kernel basis, taking the firstqqprincipal singular directions to form a well\-conditioned orthogonal kernel basisUq​\(𝐱\)U\_\{q\}\(\\mathbf\{x\}\)\(the criterion for selecting the retained directions is given in Appendix[B](https://arxiv.org/html/2609.35938#A2): usually the principal singular directions of the boundary kernel matrix are taken, but when the magnitude of the evaluation block differs greatly from that of the boundary block the retained directions must instead be selected from the input subspace, which is the case for the far field of Example 4\), so that the operator output in both the collocation and the boundary integral form can be written uniformly as𝒢θ​\(a\)​\(𝐱\)=Uq​\(𝐱\)​𝐚\\mathcal\{G\}\_\{\\theta\}\(a\)\(\\mathbf\{x\}\)=U\_\{q\}\(\\mathbf\{x\}\)\\,\\mathbf\{a\}; this preprocessing does not destroy the physical consistency of the expansion\. Its construction details, truncation error analysis, continuous basis extension and coupling with source\-point learning and offline preprocessing are given in Appendix[B](https://arxiv.org/html/2609.35938#A2)\.

#### 2\.4\.3KernelOnet\-HK: Hybrid Kernel

When the governing equation contains nonlinear terms and its operator admits no usable fundamental solution, the principle of superposition no longer holds and the whole equation cannot be represented by a linear superposition of a single fundamental solution\. If KernelOnet\-PIKF of Section 2\.4\.2 were still used in this case, with the fundamental solution of its linear principal part as the kernel, the expansion would satisfy the homogeneous equationℒ​uh=0\\mathcal\{L\}u\_\{h\}=0identically and be unable to represent the solution driven by the nonlinear remainder, thus losing its approximation capability\. For this reason, this section proposes a third construction strategy for KernelOnet—the hybrid kernel \(HK\)—whose central idea is to decompose the solution according to the linear principal part of the governing equation: the solution is written as the sum of a homogeneous branch and a correction branch,

u=uh\+up,ℒ​uh=0,ℒ​up=−N⁡\[u\],u=u\_\{h\}\+u\_\{p\},\\qquad\\mathcal\{L\}\\,u\_\{h\}=0,\\qquad\\mathcal\{L\}\\,u\_\{p\}=\-N\[u\],\(23\)whereℒ\\mathcal\{L\}is the linear principal part of the governing equation andN⁡\[u\]N\[u\]is the nonlinear remainder\. Hereuhu\_\{h\}is the homogeneous branch, lying in the null space ofℒ\\mathcal\{L\}and representable by an expansion of the fundamental solutionΦ\\Phiofℒ\\mathcal\{L\}; the nonlinearity is concentrated entirely in the correction branchupu\_\{p\}, which is determined by the source term−N⁡\[u\]\-N\[u\]and in general no longer satisfiesℒ​up=0\\mathcal\{L\}u\_\{p\}=0, so that it must be represented by a family of basis functions that do not satisfy this homogeneous equation\. Accordingly, KernelOnet\-HK takes the kernel function to be a mixture of an analytic fundamental\-solution kernel and a low\-rank learned correction kernel\. As in Section 2\.4\.2, the homogeneous branch admits two equivalent forms, collocation \(aγ\\gamma\-shifted fundamental solution\) and boundary integral \(a single\-layer potential kernel on the true boundary\); the collocation form is taken as an example below:

𝒢θ​\(a\)​\(𝐱\)=∑j=1Bbj​\(a\)​Φ​\(\|𝐱−γ​𝐱b\(j\)\|\)⏟analytic fundamental\-solution kernel \(satisfies​ℒ​uh=0​, a physical prior\)\+∑k=1Kcck​\(a\)​φθ​\(\|𝐱−𝐭k\|\)⏟learned source\-term correction kernel \(low\-rank,​Kc≪B​\),\\mathcal\{G\}\_\{\\theta\}\(a\)\(\\mathbf\{x\}\)=\\underbrace\{\\sum\_\{j=1\}^\{B\}b\_\{j\}\(a\)\\,\\Phi\\\!\\left\(\\left\|\\mathbf\{x\}\-\\gamma\\,\\mathbf\{x\}\_\{b\}^\{\(j\)\}\\right\|\\right\)\}\_\{\\text\{analytic fundamental\-solution kernel \(satisfies \}\\mathcal\{L\}u\_\{h\}=0\\text\{, a physical prior\)\}\}\+\\underbrace\{\\sum\_\{k=1\}^\{K\_\{c\}\}c\_\{k\}\(a\)\\,\\varphi\_\{\\theta\}\\\!\\left\(\\left\|\\mathbf\{x\}\-\\mathbf\{t\}\_\{k\}\\right\|\\right\)\}\_\{\\text\{learned source\-term correction kernel \(low\-rank, \}K\_\{c\}\\ll B\\text\{\)\}\},\(24\)whereΦ\\Phiis the analytic fundamental solution of the linear principal partℒ\\mathcal\{L\}\(the definition ofγ\\gammais given in Section 2\.4\.2\);φθ\\varphi\_\{\\theta\}is a shallow network taking only the radial distance as input, belonging to the same family as the kernel function of KernelOnet\-RBF;\{𝐭k\}k=1Kc\\left\\\{\\mathbf\{t\}\_\{k\}\\right\\\}\_\{k=1\}^\{K\_\{c\}\}areKcK\_\{c\}fixed centres of the correction kernel, sampled inside the solution domain by Latin hypercube sampling \(LHS\) \(for unbounded exterior problems the centres are taken in a bounded region inside the obstacle\); andbj​\(a\)b\_\{j\}\(a\)andck​\(a\)c\_\{k\}\(a\)are output at once by the same branch network, whose output dimension isB\+KcB\+K\_\{c\}\. Since the parametersθ\\thetaofφθ\\varphi\_\{\\theta\}are shared across all samples and only the coefficientsck​\(a\)c\_\{k\}\(a\)vary with the input, for any given input the correction term can only lie in the

𝒱Kc=span⁡\{φθ​\(\|𝐱−𝐭1\|\),…,φθ​\(\|𝐱−𝐭Kc\|\)\}\\mathcal\{V\}\_\{K\_\{c\}\}=\\mathrm\{span\}\\left\\\{\\varphi\_\{\\theta\}\\\!\\left\(\\left\|\\mathbf\{x\}\-\\mathbf\{t\}\_\{1\}\\right\|\\right\),\\ldots,\\varphi\_\{\\theta\}\\\!\\left\(\\left\|\\mathbf\{x\}\-\\mathbf\{t\}\_\{K\_\{c\}\}\\right\|\\right\)\\right\\\}\(25\)KcK\_\{c\}\-dimensional subspace\. ThusKc≪BK\_\{c\}\\ll Blimits the dimensionality occupied by the data\-driven part, so that the analytic fundamental\-solution kernel continues to dominate the expansion and the model does not degenerate into the purely data\-driven KernelOnet\-RBF, while the correction branch is guaranteed to have a clear physical meaning rather than being a black box\.

The hybrid kernel also connects the two endpoints of Sections 2\.4\.1 and 2\.4\.2 through a single tunable dimensionKcK\_\{c\}: whenKc=0K\_\{c\}=0the correction branch disappears and the model degenerates to KernelOnet\-PIKF; whenKc=BK\_\{c\}=B, and the centres and basis\-function form of the correction kernels agree with those of KernelOnet\-RBF and the weights of the first branch tend to zero, it degenerates to KernelOnet\-RBF\. The three variants are therefore not independent of one another but form a continuous spectrum from "strong physical prior" to "fully data\-driven": PIKF and RBF are the two endpoints and HK is the intermediate form continuously tuned between them byKcK\_\{c\}\.It should be noted that the RBF limit is idealized: it requires the centres and the basis\-function form of the correction kernel to agree with those of KernelOnet\-RBF\.

Since the correction branch destroys the property that the expansion satisfies the governing equation exactly, training KernelOnet\-HK requires interior solution data for supervised learning, and its data loss is

𝒥data​\(θ\)=1N​nt​∑i=1N∑k=1nt\(upred,k\(i\)−utrue,k\(i\)\)2\.\\mathcal\{J\}\_\{\\mathrm\{data\}\}\(\\theta\)=\\frac\{1\}\{N\\,n\_\{t\}\}\\sum\_\{i=1\}^\{N\}\\sum\_\{k=1\}^\{n\_\{t\}\}\\left\(u\_\{\\mathrm\{pred\},k\}^\{\(i\)\}\-u\_\{\\mathrm\{true\},k\}^\{\(i\)\}\\right\)^\{2\}\.\(26\)The physical prior enters the model structure only through the analytic kernel basis \(as a hard constraint\) and the training loss contains no PDE residual term; compared with purely data\-driven methods that require large numbers of interior solution labels, its training does not rely on interior labels at all\. It should be emphasized that the treatment in this paper introduces no inference\-time iteration: the value of operator learning lies in amortized inference—once training is complete, a single forward pass yields the solution for a new input function\. This construction is consistent with the idea described in Section 2\.3\.2 of "taking the fundamental solution of the linear principal part as the homogeneous kernel and absorbing the nonlinearity into a particular\-solution kernel for the source term", the difference being that here the particular\-solution kernel is obtained from data through the low\-rank learned branch\. The hybrid kernel also brings interpretability: after training, the coefficient norms‖b⁡\(a\)‖\\\|b\(a\)\\\|and‖c⁡\(a\)‖\\\|c\(a\)\\\|of the homogeneous and correction branches can be read off separately, so as to determine quantitatively the proportions of the "homogeneous component" \(ℒ​uh=0\\mathcal\{L\}u\_\{h\}=0\) and the "source component" in the solution; the correction kernel centres\{𝐭k\}\\left\\\{\\mathbf\{t\}\_\{k\}\\right\\\}and the learned radial shapeφθ\\varphi\_\{\\theta\}also provide physical insight into the spatial distribution of the nonlinear source term\. When the governing equation itself admits a fundamental solution,N⁡\[u\]≡0N\[u\]\\equiv 0in Eq\. \([23](https://arxiv.org/html/2609.35938#S2.E23)\), and one may takeup≡0u\_\{p\}\\equiv 0, so that KernelOnet\-HK degenerates to KernelOnet\-PIKF and linear problems are still handled by PIKF in an unsupervised manner\.

#### 2\.4\.4Comparison and Summary

Table[1](https://arxiv.org/html/2609.35938#S2.T1)summarizes the similarities and differences among DeepONet, PI\-DeepONet and the three KernelOnet variants from seven aspects: the form of the basis functions, the way physical information is embedded, PDE consistency, the supervision information, the learnable kernel parameters, and the applicability to nonlinear and unbounded\-domain problems\. As regards the form of the basis functions, DeepONet and PI\-DeepONet learn the basis functions implicitly through a network, whereas KernelOnet makes the basis functions explicit as kernel functions—a neural\-network\-parameterized radial basis function for KernelOnet\-RBF, an analytic fundamental\-solution kernel \(with two equivalent forms, collocation and boundary integral\) for KernelOnet\-PIKF, and a mixture of the analytic fundamental\-solution kernel with low\-rank learned correction kernels for KernelOnet\-HK \(Eq\. \([24](https://arxiv.org/html/2609.35938#S2.E24)\)\)\. As regards the way physical information is embedded and PDE consistency, DeepONet makes no use of physical information at all and cannot satisfy PDE consistency, PI\-DeepONet adds the PDE residual to the loss as a soft constraint and satisfies the equation only approximately at the collocation points, the kernel function of KernelOnet\-RBF does not satisfy the governing equation in advance, KernelOnet\-PIKF hard\-codes the physical information into the kernel function so that the expansion satisfies the governing equation exactly, and in KernelOnet\-HK the homogeneous branch satisfies the linear principal part exactly while the correction branch carries the nonlinear source term, so that the expansion as a whole is satisfied only approximately\. As regards the supervision information, DeepONet, KernelOnet\-RBF and KernelOnet\-HK all require interior solution labels, PI\-DeepONet requires boundary conditions together with the PDE residual, and KernelOnet\-PIKF requires only boundary conditions and can be trained without supervision\. As regards the learnable kernel parameters, the basis functions of the two DeepONet variants are learned by the network as a whole and require no dedicated kernel parameters \(PI\-DeepONet requires tuning of the loss weightsλ\\lambda\), KernelOnet\-RBF learns a shallow kernel function network, KernelOnet\-PIKF learns only one scalarγ\\gammain its collocation form \(the boundary integral form contains no learnable kernel parameter\), and KernelOnet\-HK learns a low\-rank correction kernel networkφθ\\varphi\_\{\\theta\}\(whoseKcK\_\{c\}correction kernel centres are fixed sampling points, withKc≪BK\_\{c\}\\ll B\)—the strong inductive bias brought by explicit kernel functions reduces the number of learnable parameters and yields higher accuracy and better interpretability for the same amount of data\. As regards the range of applicability, for nonlinear problems KernelOnet\-PIKF does not apply, KernelOnet\-HK applies only when the linear principal part admits an analytic fundamental solution, and KernelOnet\-RBF can cover more general problems such as variable coefficients and strong nonlinearities; for unbounded\-domain problems, only KernelOnet\-PIKF and KernelOnet\-HK, which take analytic fundamental solutions as kernels, are applicable\.

In summary, the three kernel construction strategies proposed in this section form a continuous spectrum from a strong physical prior to fully data\-driven learning, complementing rather than replacing one another\. This trade\-off between "strength of the physical prior" and "breadth of applicability" is precisely the design motivation behind the three variants in Figure[2](https://arxiv.org/html/2609.35938#S2.F2)sharing one architecture and differing only in the kernel trunk, and it will be systematically verified in the four examples of Section 3\.

Each of the three variants continues one of the existing kernel\-expansion routes of Section 2\.3: KernelOnet\-RBF continues the domain\-type collocation route and can be regarded as the counterpart in operator learning of Kansa’s method with general radial basis functions as kernels; KernelOnet\-PIKF continues the boundary\-type route, its collocation form being of the same origin as the method of fundamental solutions, the Trefftz method and the singular boundary method, while its boundary integral form shares the same integral kernel as the single\-layer potential of the boundary element method; and KernelOnet\-HK continues the nonlinear extension route, its construction of "fundamental solution of the linear principal part plus source\-term correction" being consistent with Newton collocation and the dual reciprocity method[Nardini and Brebbia \(1983\)](https://arxiv.org/html/2609.35938#bib.bib43);[Partridge et al\. \(1991\)](https://arxiv.org/html/2609.35938#bib.bib45)\. The three thus bring all three technical routes of traditional kernel\-expansion methods into the operator learning framework\.

Table 1:Comparison of DeepONet, PI\-DeepONet, KernelOnet\-RBF, KernelOnet\-PIKF and KernelOnet\-HK\.

### 2\.5Theoretical Analysis: Convergence and Accuracy

This section analyses the convergence and accuracy of the three KernelOnet variants at the level of approximation theory, starting from kernel function theory, in order to clarify the mathematical basis of the trade\-off between "strength of the physical prior" and "breadth of applicability" described in Section 2\.4\. For the sake of uniformity, the operator expansions of the three variants are all regarded as approximations of the solution in an approximation space spanned by kernel functions: givenBBcentres \(source points\)\{𝐱s\(j\)\}j=1B\\left\\\{\\mathbf\{x\}\_\{s\}^\{\(j\)\}\\right\\\}\_\{j=1\}^\{B\}, define the approximation space

𝒱B=span⁡\{ψ1,ψ2,…,ψB\},\\mathcal\{V\}\_\{B\}=\\mathrm\{span\}\\left\\\{\\psi\_\{1\},\\psi\_\{2\},\\ldots,\\psi\_\{B\}\\right\\\},\(27\)so that the prediction𝒢θ​\(a\)\\mathcal\{G\}\_\{\\theta\}\(a\)of KernelOnet is an element of𝒱B\\mathcal\{V\}\_\{B\}whose coefficients are given by the branch network\. For KernelOnet\-HK the approximation space additionally containsKcK\_\{c\}low\-rank correction kernels, that is, the direct sum𝒱B⊕𝒱Kc\\mathcal\{V\}\_\{B\}\\oplus\\mathcal\{V\}\_\{K\_\{c\}\}\(withKc≪BK\_\{c\}\\ll B\), and its error is correspondingly divided into a homogeneous\-branch part and a source\-branch part, corresponding respectively touhu\_\{h\}andupu\_\{p\}in Eq\. \([23](https://arxiv.org/html/2609.35938#S2.E23)\)\. Different from the above division according to the components of the solution \(homogeneous branch and source branch\), the total error can also be decomposed by origin into two parts: first, the "approximation error" of the approximation space𝒱B\\mathcal\{V\}\_\{B\}with respect to the true solution, that is, whether the shape of the kernel function "fits" the true solution; and second, the "estimation error" introduced by the determination of the coefficients and by training\. The differences among the three variants are essentially differences in the former, that is, in the choice of the kernel function space\.

##### \(1\) Error bound in the reproducing kernel space framework—the accuracy basis of KernelOnet\-RBF

KernelOnet\-RBF uses the radial basis functionψ​\(r\)=φθ​\(r\)\\psi\(r\)=\\varphi\_\{\\theta\}\(r\)as its kernel\. In reproducing kernel Hilbert space \(RKHS\) theory[Aronszajn \(1950\)](https://arxiv.org/html/2609.35938#bib.bib19);[Schaback and Wendland \(2006\)](https://arxiv.org/html/2609.35938#bib.bib30), a radial basis function is rigorously regarded as the reproducing kernel of some function space \(its native space\), from which the standard error bound for interpolation approximation follows\. It should be pointed out that this bound presupposes that the kernel function is strictly or conditionally positive definite; theφθ\\varphi\_\{\\theta\}learned by KernelOnet\-RBF does not automatically satisfy this condition, so that the bound holds only when it is \(approximately\) positive definite or is made positive definite by regularization\. LetΩ\\Omegabe the solution domain, letφ\\varphibe a strictly or conditionally positive definite radial basis function with native space𝒩φ​\(Ω\)\\mathcal\{N\}\_\{\\varphi\}\(\\Omega\), and letX=\{𝐱s\(j\)\}X=\\\{\\mathbf\{x\}\_\{s\}^\{\(j\)\}\\\}be the discrete set of centres; then for anyu∈𝒩φ​\(Ω\)u\\in\\mathcal\{N\}\_\{\\varphi\}\(\\Omega\), the RBF interpolation approximantsssatisfies the pointwise bound

\|u⁡\(𝐱\)−s⁡\(𝐱\)\|≤Pφ,X​\(𝐱\)​\|u\|𝒩φ,\|u\(\\mathbf\{x\}\)\-s\(\\mathbf\{x\}\)\|\\leq P\_\{\\varphi,X\}\(\\mathbf\{x\}\)\\,\|u\|\_\{\\mathcal\{N\}\_\{\\varphi\}\},\(28\)wherePφ,X​\(𝐱\)P\_\{\\varphi,X\}\(\\mathbf\{x\}\)is the power function associated with the kernel function and the centre distribution and\|u\|𝒩φ\|u\|\_\{\\mathcal\{N\}\_\{\\varphi\}\}is the semi\-norm of the native space[Wendland \(2005\)](https://arxiv.org/html/2609.35938#bib.bib24);[Buhmann \(2003\)](https://arxiv.org/html/2609.35938#bib.bib23)\. The dependence of the error on the density of the centre distribution is characterized by the fill distancehX,Ω=sup𝐱∈Ωminj⁡\|𝐱−𝐱s\(j\)\|h\_\{X,\\Omega\}=\\sup\_\{\\mathbf\{x\}\\in\\Omega\}\\min\_\{j\}\|\\mathbf\{x\}\-\\mathbf\{x\}\_\{s\}^\{\(j\)\}\|: as the centres become denser,hX,Ω→0h\_\{X,\\Omega\}\\to 0and the approximation error converges at an order determined by the smoothness of the kernel function—algebraic convergence for finitely smooth kernels, and spectral convergence for smooth kernels such as the Gaussian and multiquadric kernels\. At the same time, Schaback and Wendland[Schaback and Wendland \(2006\)](https://arxiv.org/html/2609.35938#bib.bib30)revealed the well\-known "uncertainty principle" of kernel methods: there is an intrinsic trade\-off between the error and the condition number, the two cannot be improved simultaneously, and good accuracy is inevitably accompanied by severe ill\-conditioning\. This result explains theoretically the characteristics of KernelOnet\-RBF—its kernel function needs finite truncation and regularization to maintain numerical well\-posedness, while under data\-driven learning it adaptively learns a kernel shape that "fits" the true solution, thereby gaining approximation room within a broader family of kernel shapes\.

##### \(2\) Spectral convergence of the fundamental\-solution expansion—the accuracy basis of KernelOnet\-PIKF

KernelOnet\-PIKF uses the fundamental solutionΦ\\Phiof the governing equation as its kernel\. Since the source points lie outside the domain, the fundamental solution satisfies the homogeneous governing equation inside the domain, so that the expansion𝒢θ​\(a\)\\mathcal\{G\}\_\{\\theta\}\(a\)automatically satisfies the governing equation and its approximation error is determined entirely by the residual of the boundary condition fitting; this is mathematically the same structure as the method of fundamental solutions \(MFS\)[Golberg and Chen \(1999\)](https://arxiv.org/html/2609.35938#bib.bib40);[Fairweather and Karageorghis \(1998\)](https://arxiv.org/html/2609.35938#bib.bib38)and the Trefftz method[Kita and Kamiya \(1995\)](https://arxiv.org/html/2609.35938#bib.bib41)\. For such boundary\-type methods, when the solution is analytic or sufficiently smooth in the domain, the family of fundamental solutions is inherently of high resolution and the approximation error exhibits "spectral convergence" as the number of basis functionsBBincreases, approaching exponential or geometric decay in the sufficiently smooth case, its convergence rate depending only on the analyticity of the solution and on the configuration of the source pointsγ\\gamma, and not on any mesh size\. In particular, by the denseness of the solution space of elliptic equations \(a Runge\-type approximation theorem\), the space spanned by the family of fundamental solutions with source points distributed outside the domain is dense in the space of analytic solutions of the homogeneous equation, and the approximation theory of the MFS guarantees the spectral accuracy of the interpolation approximation, which provides theoretical support for "obtaining high accuracy with the fewest degrees of freedom"; moreover, onceγ\\gammadeparts from 1 the family of fundamental solutionsΦ⁡\(\|𝐱−γ​𝐱b\|\)\\Phi\(\|\\mathbf\{x\}\-\\gamma\\mathbf\{x\}\_\{b\}\|\)tends to become linearly dependent and its condition number rises as the source distance increases, which is exactly why Section 2\.4 makes the source\-point scaling a learnable scalar adaptedγ\\gamma, so as to compromise between avoiding the singularity and maintaining well\-posedness\. In summary, the core advantage of KernelOnet\-PIKF is that its kernel space𝒱B\\mathcal\{V\}\_\{B\}is constrained within an approximation subspace of the homogeneous solution space, so that the approximation space matches the physical structure of the solution naturally and the approximation error can be made arbitrarily small; this is the root of its ability to achieve high accuracy with very few parameters and in an unsupervised manner\.

##### \(3\) Decomposition error of the source\-term correction—bias analysis of KernelOnet\-HK

KernelOnet\-HK takes the approximation space to be the direct sum of the analytic fundamental\-solution subspace and the low\-rank correction subspace,𝒱B⊕𝒱Kc\\mathcal\{V\}\_\{B\}\\oplus\\mathcal\{V\}\_\{K\_\{c\}\}, and its total error is correspondingly divided into two parts according to the decomposition of Eq\. \([23](https://arxiv.org/html/2609.35938#S2.E23)\):

u−𝒢θ​\(a\)=\(uh−uh,θ\)⏟approximation error of the homogeneous branch\+\(up−up,θ\)⏟approximation error of the source branch,u\-\\mathcal\{G\}\_\{\\theta\}\(a\)=\\underbrace\{\\bigl\(u\_\{h\}\-u\_\{h,\\theta\}\\bigr\)\}\_\{\\text\{approximation error of the homogeneous branch\}\}\+\\underbrace\{\\bigl\(u\_\{p\}\-u\_\{p,\\theta\}\\bigr\)\}\_\{\\text\{approximation error of the source branch\}\},\(29\)whereuh,θ=∑j=1Bbj​\(a\)​Φ​\(\|𝐱−γ​𝐱b\(j\)\|\)u\_\{h,\\theta\}=\\sum\_\{j=1\}^\{B\}b\_\{j\}\(a\)\\,\\Phi\(\|\\mathbf\{x\}\-\\gamma\\mathbf\{x\}\_\{b\}^\{\(j\)\}\|\)is spanned by the analytic fundamental solutions andup,θ=∑k=1Kcck​\(a\)​φθ​\(\|𝐱−𝐭k\|\)u\_\{p,\\theta\}=\\sum\_\{k=1\}^\{K\_\{c\}\}c\_\{k\}\(a\)\\,\\varphi\_\{\\theta\}\(\|\\mathbf\{x\}\-\\mathbf\{t\}\_\{k\}\|\)is spanned by the low\-rank correction kernels\. Note that the decomposition of Eq\. \([23](https://arxiv.org/html/2609.35938#S2.E23)\) is not unique—uhu\_\{h\}andupu\_\{p\}may differ by aℒ\\mathcal\{L\}\-homogeneous solution; this degree of freedom is precisely absorbed by the homogeneous branch \(it suffices to takeuhu\_\{h\}to be the homogeneous part matching the boundary conditions\), so that the error analysis below is insensitive to it\. The error of the homogeneous branch is exactly the same as in the case of KernelOnet\-PIKF in item \(2\): since every fundamental solution satisfiesℒ=0\\mathcal\{L\}=0inside the domain, its approximation error is determined only by the boundary fitting residual and converges spectrally asBBincreases\. The error of the source branch, on the other hand, is determined by the "compressibility" ofupu\_\{p\}in the correction subspace𝒱Kc\\mathcal\{V\}\_\{K\_\{c\}\}: by the best approximation theorem,

‖up−up,θ‖Ω≤\(1\+Λ\)​infv∈𝒱Kc‖up−v‖Ω,\\\|u\_\{p\}\-u\_\{p,\\theta\}\\\|\_\{\\Omega\}\\leq\\bigl\(1\+\\Lambda\\bigr\)\\inf\_\{v\\in\\mathcal\{V\}\_\{K\_\{c\}\}\}\\\|u\_\{p\}\-v\\\|\_\{\\Omega\},\(30\)whereΛ\\Lambdais a stability constant related to sampling and optimization and the second factor on the right\-hand side is the best approximation error ofupu\_\{p\}in𝒱Kc\\mathcal\{V\}\_\{K\_\{c\}\}\. Therefore, as long as the nonlinear source termupu\_\{p\}can be compressed efficiently in theKcK\_\{c\}\-dimensional subspace—for example when its energy is concentrated in a few modes—the error of the source branch decays rapidly withKcK\_\{c\}; otherwise it decays slowly\.

Applying the linear principal part to Eq\. \([24](https://arxiv.org/html/2609.35938#S2.E24)\) gives the PDE residual of the whole expansion\. Since the homogeneous branch satisfiesℒ\\mathcal\{L\}exactly and contributes nothing, we obtain

ℛθ​\(𝐱\)=ℒ​𝒢θ​\(𝐱\)\+N⁡\[𝒢θ\]​\(𝐱\)=∑k=1Kcck​\(a\)​ℒ​φθ​\(\|𝐱−𝐭k\|\)\+N⁡\[𝒢θ\]​\(𝐱\)\.\\mathcal\{R\}\_\{\\theta\}\(\\mathbf\{x\}\)=\\mathcal\{L\}\\mathcal\{G\}\_\{\\theta\}\(\\mathbf\{x\}\)\+N\[\\mathcal\{G\}\_\{\\theta\}\]\(\\mathbf\{x\}\)=\\sum\_\{k=1\}^\{K\_\{c\}\}c\_\{k\}\(a\)\\,\\mathcal\{L\}\\varphi\_\{\\theta\}\\\!\\left\(\\left\|\\mathbf\{x\}\-\\mathbf\{t\}\_\{k\}\\right\|\\right\)\+N\[\\mathcal\{G\}\_\{\\theta\}\]\(\\mathbf\{x\}\)\.\(31\)When theℒ\\mathcal\{L\}\-image ofφθ\\varphi\_\{\\theta\}is made to approximate a set of source\-term basis functionsχk\\chi\_\{k\}, the sum∑kck​χk\\sum\_\{k\}c\_\{k\}\\chi\_\{k\}is the low\-rank expansion of the source term−N⁡\[u\]\-N\[u\], and Eq\. \([31](https://arxiv.org/html/2609.35938#S2.E31)\) is precisely the remainder of that expansion\. The residual is therefore determined entirely by the approximation capability of the correction branch and is independent of the homogeneous branch, with a clear physical origin: the low\-rank source\-term subspace carries the error and can be reduced adaptively withKcK\_\{c\}and the compressibility of the source term , representing−N⁡\[u\]\-N\[u\]spatially adaptively as a function of𝐱\\mathbf\{x\}\.

TheKcK\_\{c\}in Eq\. \([24](https://arxiv.org/html/2609.35938#S2.E24)\) is therefore an explicit "bias–capacity" knob: whenKc=0K\_\{c\}=0we have𝒱Kc=\{0\}\\mathcal\{V\}\_\{K\_\{c\}\}=\\\{0\\\}, KernelOnet\-HK degenerates to KernelOnet\-PIKF and its error reduces to a purely homogeneous approximation error, which for nonlinear problems manifests itself as an irreducible model bias; asKcK\_\{c\}increases, the approximation capability of the source branch is enhanced and the approximation error of Eq\. \([30](https://arxiv.org/html/2609.35938#S2.E30)\) decreases, but the number of coefficients and network parameters to be calibrated grows and the dependence on interior labels increases, so that the estimation error rises accordingly, and the model approaches the fully data\-driven KernelOnet\-RBF asKc→BK\_\{c\}\\to Bwith the centres and basis\-function form of the correction kernel agreeing with those of KernelOnet\-RBF \(an idealized limit\)\. This is the theoretical basis on which KernelOnet\-HK exchanges a controllable physical bias for the ability to handle nonlinearity: its accuracy ceiling is determined by the compressibility of the source term in a low\-dimensional subspace\.

Taken together, the errors of the three variants can be decomposed uniformly into an approximation error and an estimation error, where the approximation error is determined by the choice of the kernel space \(for KernelOnet\-HK by the direct sum𝒱B⊕𝒱Kc\\mathcal\{V\}\_\{B\}\\oplus\\mathcal\{V\}\_\{K\_\{c\}\}\) and the estimation error by the amount of training data and the optimization process\. Table[2](https://arxiv.org/html/2609.35938#S2.T2)compares the three at the theoretical level in terms of the properties of the kernel space, the origin of the error and the convergence\. It can be seen that the differences among the three essentially reduce to a compromise between "the agreement between the kernel space and the physical structure of the true solution" and "the flexibility of the approximation space": KernelOnet\-PIKF exchanges a highly restricted but physically consistent kernel space for a minimal approximation error and unsupervised training; KernelOnet\-RBF exchanges a flexible kernel space for adaptability to problems that differ from the fundamental solution; and KernelOnet\-HK superimposes a low\-rank source\-term correction on the spectral approximation of the homogeneous branch, withKcK\_\{c\}as a knob providing a continuous transition between the two ends\.

Table 2:Comparison of KernelOnet\-RBF, KernelOnet\-PIKF and KernelOnet\-HK from the perspective of kernel function theory\.It should be noted that the above analysis provides, at the level of approximation theory, a "qualitative or semi\-quantitative description of the convergence order and the origin of the error"; the constants and the specific convergence rates also depend on practical factors such as the smoothness of the problem solution, the regularization parameters of the kernel function, the capacity of the branch network and the degree of convergence of the optimization, which are usually difficult to identify analytically in advance\. This is exactly why Section 3 systematically verifies the accuracy of each method through numerical experiments\. In addition, universal approximation at the operator level—that is, the ability of the branch network to approximate the coefficient functions—is guaranteed by the classical universal approximation theorem for neural networks[Chen and Chen \(1995\)](https://arxiv.org/html/2609.35938#bib.bib48);[Park and Sandberg \(1991\)](https://arxiv.org/html/2609.35938#bib.bib21), so that KernelOnet possesses universal approximation capability for the target operator when the number of kernel functionsB→∞B\\to\\inftyand the capacity of the branch network is sufficient \(for KernelOnet\-PIKF this conclusion presupposes that the target solution belongs to, or can be approximated by, the homogeneous solution space\); the actual performance for finiteBBand finite data is then determined jointly by the approximation error and the estimation error discussed above\.

## 3Numerical examples and discussions

In this section, we evaluate the performance of the proposed KernelOnet framework through four representative numerical examples: the first three are benchmark cases that verify the accuracy, interpretability, and physical consistency of the framework on idealized mathematical problems from different perspectives, while the last one is an engineering\-oriented example of underwater acoustic radiation in a shallow\-water waveguide that tests the framework in a realistic complex physical scenario: \(1\) the Laplace equation on a circular domain, which is used to examine whether the data\-driven kernel can recover the fundamental solution from data when the expansion strictly satisfies the governing equation; \(2\) the nonlinear modified Helmholtz equation on a star\-shaped domain, which is used to examine the applicability of the hybrid kernel when no analytic fundamental solution exists; \(3\) the complex Helmholtz equation in an unbounded exterior domain, which is used to examine how the analytic fundamental\-solution kernel handles traveling waves and the radiation condition, and to compare two implementations, namely boundary collocation and boundary integral; \(4\) underwater acoustic radiation and propagation induced by spherical\-shell vibration in a shallow\-water waveguide, which tests the ability of the framework to solve problems in a realistic complex physical scenario\. The model performance is measured by the mean of the sample\-level relativeL2L\_\{2\}error:

ℰ=1N​∑i=1N∑k=1nt\(upred,k\(i\)−utrue,k\(i\)\)2∑k=1nt\(utrue,k\(i\)\)2,\\mathcal\{E\}=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\\frac\{\\sqrt\{\\sum\_\{k=1\}^\{n\_\{t\}\}\\left\(u\_\{\\mathrm\{pred\},k\}^\{\(i\)\}\-u\_\{\\mathrm\{true\},k\}^\{\(i\)\}\\right\)^\{2\}\}\}\{\\sqrt\{\\sum\_\{k=1\}^\{n\_\{t\}\}\\left\(u\_\{\\mathrm\{true\},k\}^\{\(i\)\}\\right\)^\{2\}\}\},\(32\)whereupred,k\(i\)u\_\{\\mathrm\{pred\},k\}^\{\(i\)\}andutrue,k\(i\)u\_\{\\mathrm\{true\},k\}^\{\(i\)\}are the predicted and reference values, respectively, of theii\-th sample at thekk\-th interior evaluation point; for the complex\-valued problems of Case 3 and Case 4, the error is computed separately for the real and imaginary parts and then averaged\.

All experiments were carried out on a Linux server configured with an AMD Ryzen 9 9950X3D processor, 96 GB of memory, an NVIDIA GeForce RTX 5090 D v2 graphics card with 24 GB of GPU memory, and the Ubuntu 24\.04 operating system\. All models are implemented in PyTorch: the branch network of every method is a three\-layer fully connectedBranchNet⁡\[160,160,160\]\\mathrm\{BranchNet\}\[160,160,160\]withtanh\\tanhactivation, trained with the Adam optimizer[Kingma and Ba \(2014\)](https://arxiv.org/html/2609.35938#bib.bib18)at a learning rate of10−410^\{\-4\}\. The input dimension of the branch network equals the number of boundary pointsnbn\_\{b\}of each case: 160 for Case 1 and Case 3, 200 for Case 2, and 50 for Case 4; for Case 4, the complex acoustic pressure is encoded separately in terms of its real and imaginary parts, so the actual input dimension is2×50=1002\\times 50=100\. The training and test sets of all four cases contain 2000 samples, and the number of training epochs is uniformly set to5×1055\\times 10^\{5\}\.

The construction of the trunk network of each method is detailed in Section 2\.4: DeepONet uses a standard fully connected trunkTrunkNet⁡\[160,160,160\]\\mathrm\{TrunkNet\}\[160,160,160\]; KernelOnet\-RBF uses a data\-driven radial basis function networkRBFTrunk⁡\[160,160\]\\mathrm\{RBFTrunk\}\[160,160\]; KernelOnet\-PIKF uses an analytic fundamental\-solution kernel trunk, whose collocation form is theγ\\gamma\-shifted fundamental solution and whose boundary integral form is the single\-layer potential kernel on the true boundary; and KernelOnet\-HK appends a low\-rank radial correction branch to the analytic fundamental\-solution branch, this shallow network taking the radial distance as input, with hidden layers\[160,160\]\[160,160\]andtanh\\tanhactivation, and with fixed sampling points as the correction kernel centers\.

As for data generation, in Case 1 and Case 3 the boundary conditions are sampled from a Gaussian random field and the reference solutions are generated by the method of fundamental solutions, while in Case 2 the boundary conditions are sampled from a Gaussian random field and the reference solution is generated by the finite element method; in Case 4 the acoustic pressure on the shell surface is synthesized by randomly superposing several vibration modes, and the reference solutions are generated separately for the near and far fields, with the Pekeris waveguide Green’s function taken as the kernel in the near field and the normal\-mode Green’s function as the kernel in the far field, the fields being synthesized after the virtual source strengths are determined by least squares from the Dirichlet data on the shell surface, and with a separate dataset generated for each of the three sound speed profiles in the far field\. It should be noted that not all methods require the interior reference solutions provided by the training set: the unsupervised configurations \(PI\-DeepONet and KernelOnet\-PIKF\) take only the boundary data from it, and their training involves no interior solution labels at all\. The reference solutions generated in the dataset are therefore used only to compute the test error and to train the three supervised baselines, namely DeepONet, KernelOnet\-RBF, and KernelOnet\-HK; Case 3 and Case 4 adopt unsupervised configurations throughout, and the datasets of these two cases serve only for error evaluation and take no part in training\. As for the choice of kernel form, Case 1 uses KernelOnet\-PIKF and KernelOnet\-HK, whose analytic fundamental\-solution branches are both of collocation form; Case 2 uses KernelOnet\-HK, whose homogeneous branch is likewise an analytic fundamental solution of collocation form; Case 3 compares both the collocation and the boundary integral forms together with their SVD\-truncated variants; and Case 4 uses the collocation form with SVD truncation applied\. In addition, Case 3 and Case 4 are both unbounded exterior problems, for which only KernelOnet\-PIKF, with the analytic fundamental solution as its kernel, is directly applicable, whereas DeepONet and KernelOnet\-RBF do not apply and KernelOnet\-HK would also introduce an additional estimation error; both cases are therefore solved with KernelOnet\-PIKF\. More detailed computational costs, accuracy comparisons, and additional analyses are given in Appendix[D](https://arxiv.org/html/2609.35938#A4)\.

### 3\.1Case 1: Laplace equation on a circular domain

The first example considers the two\-dimensional Laplace equation on a circular domain of radiusR=0\.5R=0\.5:

∇2u\(x,y\)=0,\(x,y\)∈Ω=\{x2\+y2<0\.25\};u=g,𝐱∈Γ=∂Ω,\\nabla^\{2\}u\(x,y\)=0,\\quad\(x,y\)\\in\\Omega=\\\{x^\{2\}\+y^\{2\}<0\.25\\\};\\qquad u=g,\\quad\\mathbf\{x\}\\in\\Gamma=\\partial\\Omega,\(33\)
Table[3](https://arxiv.org/html/2609.35938#S3.T3)summarizes the training results of the five methods on Case 1, including the number of learnable parameters, the relativeL2L\_\{2\}error, and the average running time per 100 training epochs\. Among them, neither KernelOnet\-PIKF nor PI\-DeepONet uses interior solution labels: the former takes the boundary residual alone as its loss, whereas the latter takes the boundary residual together with the PDE residual at interior collocation points; all the remaining methods are trained in a supervised manner using interior solution labels\. It can be seen that KernelOnet\-PIKF attains the highest accuracy among the five methods with8\.89×10−48\.89\\times 10^\{\-4\}, and has the fewest total learnable parameters, only 103,041, of which the kernel parameter is a single scalarγ\\gamma\. The supervised KernelOnet\-RBF \(1\.29×10−31\.29\\times 10^\{\-3\}\) also outperforms DeepONet \(1\.89×10−31\.89\\times 10^\{\-3\}\) while using fewer parameters\. It is worth noting that this case is a linear homogeneous problem, so the low\-rank correction branch of HK is not theoretically necessary: with the correction kernel centers taken by in\-domain LHS, the trained correction branch is almost zero, itsL2L\_\{2\}energy ratio to the homogeneous branch being only about7×10−67\\times 10^\{\-6\}\. At this point HK and KernelOnet\-PIKF are virtually indistinguishable in terms of function space, and the only substantial difference between them lies in the loss function—the loss of HK contains only the interior data term and no boundary residual, so it does not directly constrain the boundary condition; its accuracy is therefore2\.04×10−32\.04\\times 10^\{\-3\}, lower than the8\.89×10−48\.89\\times 10^\{\-4\}of KernelOnet\-PIKF, which fits the boundary residual directly, and also slightly lower than the1\.29×10−31\.29\\times 10^\{\-3\}of KernelOnet\-RBF, which adaptively learns the kernel shape from data; the gap thus arises from the type of loss rather than from the correction branch itself\. PI\-DeepONet converges to3\.53×10−33\.53\\times 10^\{\-3\}, an accuracy about1/41/4that of KernelOnet\-PIKF, and takes 0\.53 s per 100 epochs, about 9 times the 0\.06 s of the latter\.

Table 3:Case 1: comparison of the training results of DeepONet, PI\-DeepONet, KernelOnet\-RBF, KernelOnet\-PIKF, and KernelOnet\-HK\.Figure[3](https://arxiv.org/html/2609.35938#S3.F3)shows the predicted solutions and absolute error fields of the five methods on the same randomly selected test sample, where Figure[3](https://arxiv.org/html/2609.35938#S3.F3)\(a\) is the MFS reference solution and Figure[3](https://arxiv.org/html/2609.35938#S3.F3)\(b\) displays the distribution of the training collocation points: the boundary points \(red crosses\) are uniformly distributed on the circumference, serving both to encode the Dirichlet boundary condition and, through the radial scalingγ​𝐱b\(j\)\\gamma\\,\\mathbf\{x\}\_\{b\}^\{\(j\)\}, to generate the kernel centers; the interior evaluation points \(blue dots\) are arranged in the radial\-angular direction and are naturally denser near the center of the circle\. In terms of the predicted fields, all five methods qualitatively reproduce the overall morphology of the reference solution, but the differences among the absolute error fields lie mainly in their magnitudes: apart from PI\-DeepONet, whose error is spread widely over the domain, the errors of the remaining methods are all concentrated near the boundary\. The error of DeepONet \(Figure[3](https://arxiv.org/html/2609.35938#S3.F3)\(c\),L2=2\.14×10−3L\_\{2\}=2\.14\\times 10^\{\-3\}\) is concentrated in a ring near the boundary, indicating that a free trunk network without physical priors struggles to capture the rapid variation of the solution in the vicinity of the boundary\. KernelOnet\-PIKF \(Figure[3](https://arxiv.org/html/2609.35938#S3.F3)\(d\),L2=4\.75×10−4L\_\{2\}=4\.75\\times 10^\{\-4\}\) achieves the smallest error, which is strictly confined to a thin boundary layer and is almost zero in the interior region\. Since the PDE residual is identically zero, the error is a harmonic function whose amplitude is controlled by the boundary residual according to the maximum principle; the boundary data of this problem contain high\-frequency components, whose harmonic extension decays rapidly toward the interior, so the error exhibits a thin boundary\-layer structure\. PI\-DeepONet \(Figure[3](https://arxiv.org/html/2609.35938#S3.F3)\(e\),L2=2\.84×10−3L\_\{2\}=2\.84\\times 10^\{\-3\}\) is comparable in magnitude to DeepONet, and its error is more widely distributed over the domain\. The difference between the two lies in the way physical information is embedded: PI\-DeepONet penalizes the PDE residual at interior collocation points as a soft constraint, whereas KernelOnet\-PIKF hard\-codes the physical information into the kernel structure, the consequences of which will be further reflected in the loss convergence in Figure[4](https://arxiv.org/html/2609.35938#S3.F4)\. The errors of KernelOnet\-HK \(Figure[3](https://arxiv.org/html/2609.35938#S3.F3)\(f\),L2=2\.11×10−3L\_\{2\}=2\.11\\times 10^\{\-3\}\) and KernelOnet\-RBF \(Figure[3](https://arxiv.org/html/2609.35938#S3.F3)\(g\),L2=7\.59×10−4L\_\{2\}=7\.59\\times 10^\{\-4\}\) are likewise mainly distributed near the boundary, with magnitudes entering the range10−3∼10−410^\{\-3\}\\sim 10^\{\-4\}; the error structure of HK is similar to that of KernelOnet\-PIKF but somewhat larger in magnitude, for the reason given above\.

![Refer to caption](https://arxiv.org/html/2609.35938v1/fig3.png)Figure 3:Case 1: predicted solutions and absolute error fields for a representative test sample\. \(a\) MFS reference solution\. \(b\) Distribution of the boundary points \(red crosses, used to encode the Dirichlet data and generate the kernel centers\) and the interior evaluation points \(blue dots\)\. \(c\)–\(g\) Predicted solutions and absolute errors of DeepONet, KernelOnet\-PIKF, PI\-DeepONet, KernelOnet\-HK, and KernelOnet\-RBF, with the corresponding relativeL2L\_\{2\}errors annotated\. Note that KernelOnet\-PIKF and PI\-DeepONet did not use any interior solution labels during training\.Figure[4](https://arxiv.org/html/2609.35938#S3.F4)shows the convergence history of the training loss\. Since the loss definitions of the supervised and unsupervised groups of methods differ \(the former being the mean squared error of the interior solution and the latter the boundary residual\), they cannot be compared directly and are therefore plotted in two separate subfigures\. Figure[4](https://arxiv.org/html/2609.35938#S3.F4)\(left\) compares the three supervised methods: the losses of all three keep decreasing during training and eventually fall to about the1×10−61\\times 10^\{\-6\}level—with final\-stage averages of1\.3×10−61\.3\\times 10^\{\-6\}\(KernelOnet\-HK\),1\.7×10−61\.7\\times 10^\{\-6\}\(KernelOnet\-RBF\), and1\.9×10−61\.9\\times 10^\{\-6\}\(DeepONet\), respectively\. This indicates that with sufficient training all three supervised methods can suppress the interior fitting error to a very low level, and that the differences in their final prediction accuracy \(Table[3](https://arxiv.org/html/2609.35938#S3.T3)\) stem mainly from their ability to resolve the error in the boundary neighborhood rather than from the loss values themselves\. Figure[4](https://arxiv.org/html/2609.35938#S3.F4)\(right\) compares the two unsupervised methods\. The boundary residual of KernelOnet\-PIKF decreases monotonically from about1\.5×10−21\.5\\times 10^\{\-2\}to1\.2×10−61\.2\\times 10^\{\-6\}and converges smoothly, corresponding to the highest prediction accuracy in Table[3](https://arxiv.org/html/2609.35938#S3.T3); the weighted loss of PI\-DeepONet decreases from3\.9×10−13\.9\\times 10^\{\-1\}to about5\.7×10−45\.7\\times 10^\{\-4\}\(of which the PDE residual is about1\.9×10−41\.9\\times 10^\{\-4\}and the boundary residual about3\.8×10−53\.8\\times 10^\{\-5\}\); although it does not stagnate, its final accuracy is still lower than that of KernelOnet\-PIKF\. The reason is that PI\-DeepONet must simultaneously optimize two competing objectives, the PDE residual and the boundary residual, and the balance between the weightsλPDE\\lambda\_\{\\mathrm\{PDE\}\}andλBC\\lambda\_\{\\mathrm\{BC\}\}requires manual tuning; whereas for KernelOnet\-PIKF the PDE residual is identically zero by virtue of the kernel structure, so the optimization objective degenerates into a single boundary fitting problem, formally equivalent to a well\-conditioned least\-squares problem that requires no weight tuning and suffers no gradient conflict, and therefore converges rapidly and smoothly\. This is mutually corroborated by the fact in Table[3](https://arxiv.org/html/2609.35938#S3.T3)that PI\-DeepONet takes about nine times as long per epoch as KernelOnet\-PIKF: the hard constraint is superior to the soft constraint in accuracy, stability, and efficiency simultaneously\.

Figure 4:Case 1: convergence history of the training loss\. The left panel shows the supervised methods \(DeepONet, KernelOnet\-RBF, KernelOnet\-HK\), whose loss is the mean squared error relative to the interior reference solution; the right panel shows the unsupervised methods \(KernelOnet\-PIKF and PI\-DeepONet\), whose losses are the boundary residual and the weighted sum of the PDE residual and the boundary residual, respectively\.Unlike the fully connected trunk of DeepONet, the kernel of KernelOnet is a univariate radial function that can be plotted directly and compared pointwise with the analytic fundamental solution, thereby allowing one to "read out" the physics learned by the network\. Figure[5](https://arxiv.org/html/2609.35938#S3.F5)\(a\) is the kernelψ​\(r\)=φθ​\(r\)\\psi\(r\)=\\varphi\_\{\\theta\}\(r\)learned by KernelOnet\-RBF after training \(red solid line\)\. Although this kernel is learned entirely from data without being told anything about the Laplace operator during training, its shape is highly similar to that of the analytic fundamental solutionΦ⁡\(r\)=−12​π​ln⁡r\\Phi\(r\)=\-\\frac\{1\}\{2\\pi\}\\ln r\(blue solid line\): it likewise rises steeply asr→0r\\to 0and decays monotonically asrrincreases\. Quantitatively, fitting it with a two\-parameter affine transformation yieldsψ⁡\(r\)≈1\.05​Φ​\(r\)−0\.40\\psi\(r\)\\approx 1\.05\\,\\Phi\(r\)\-0\.40\(green dashed line\), and this fitted curve almost completely coincides with the learned kernel\. This result shows that the network autonomously "discovers" the fundamental solution of the Laplace operator from the data, differing only by a scale factor and a constant shift—the scale factor can be exactly absorbed by the coefficientsbjb\_\{j\}output by the branch network, while the constant shift corresponds to a common component of the coefficients and can likewise be adjusted by the network, so neither affects the expressive capacity\. The only essential difference between the two lies in the behavior at the origin: the analytic fundamental solution diverges atr=0r=0, whereas the learned kernel takes the finite valueψ⁡\(0\)=0\.477\\psi\(0\)=0\.477, meaning that the network automatically learns a regularized, non\-singular fundamental solution\. This provides direct numerical evidence for the claim in Section 2\.4 that "in constant\-coefficient linear problems, the learned kernel can be viewed as a regularized approximation of the fundamental solution," and also explains why KernelOnet\-RBF can achieve higher accuracy with significantly fewer parameters than DeepONet: its trunk is constrained to a rotation\-invariant function that depends only on the radial distance, so the search space is greatly reduced compared with that of a free trunk, while the function family containing the true solution remains within it\.

To further verify this claim from a functional point of view, we directly freeze the trained network and perform a boundary\-type solve with this kernel as the fundamental solution: the source points are taken as the boundary points of Case 1 \(160160points on the circle of radius0\.50\.5\), the solution domain is the square inscribed in this circle, the target solution is the harmonic functionu∗=x3−3​x​y2u^\{\*\}=x^\{3\}\-3xy^\{2\}, and the coefficients are obtained by least squares from collocation on the boundary of the square according tou⁡\(𝐱\)=∑jcj​ψ​\(\|𝐱−𝐬j\|\)u\(\\mathbf\{x\}\)=\\sum\_\{j\}c\_\{j\}\\psi\(\|\\mathbf\{x\}\-\\mathbf\{s\}\_\{j\}\|\)\. Figure[5](https://arxiv.org/html/2609.35938#S3.F5)\(b\) gives the result of this solve, whose left, middle, and right panels are the analytical solution, the NN\-based MFS solution, and the absolute error field, respectively: the prediction is visually almost indistinguishable from the analytical solution, the relativeL2L\_\{2\}error on the grid is3\.04×10−33\.04\\times 10^\{\-3\}\(maximum absolute error5\.3×10−45\.3\\times 10^\{\-4\}\), and the error is concentrated near the four vertices that coincide with the source points, i\.e\., exactly where the learned kernel deviates from the fundamental solution\. It can be seen that the learned kernel not only approximates the fundamental solution in shape but can also be used directly for boundary\-type solving just like the analytic fundamental solution\.

![Refer to caption](https://arxiv.org/html/2609.35938v1/case1_kernel_mfs.png)Figure 5:Case 1: kernel diagnostics and functional verification of the learned kernel\. \(a\) The data\-driven kernelψ​\(r\)=φθ​\(r\)\\psi\(r\)=\\varphi\_\{\\theta\}\(r\)of KernelOnet\-RBF \(red\), which is well described by the affine fit1\.05​Φ​\(r\)−0\.401\.05\\,\\Phi\(r\)\-0\.40\(green dashed line\), i\.e\., the network recovers the fundamental solutionΦ⁡\(r\)=−12​π​ln⁡r\\Phi\(r\)=\-\\frac\{1\}\{2\\pi\}\\ln rup to a scale and an offset, while remaining finite at the origin withψ⁡\(0\)=0\.477\\psi\(0\)=0\.477\. \(b\) The NN\-based MFS obtained by taking this learned kernel as the fundamental solution and collocating on the boundary of the inscribed square: left, the analytical solutionu∗=x3−3​x​y2u^\{\*\}=x^\{3\}\-3xy^\{2\}; middle, the predicted solution; right, the absolute error field \(the relativeL2L\_\{2\}error is given in the title\); the source points are160160boundary points on the circle of radius0\.50\.5, and the solution domain is the square inscribed in this circle\.
### 3\.2Case 2: nonlinear modified Helmholtz equation on a star\-shaped domain

The previous example considered a linear problem on a simple geometry\. As a further test of the applicability of the framework, this example combines "complex geometry" with "nonlinearity": we solve the modified Helmholtz equation with a cubic nonlinear term on a five\-pointed star domain

∇2u−k2​u\+ε​u3=0,𝐱∈Ω;u=g,𝐱∈Γ,\\nabla^\{2\}u\-k^\{2\}u\+\\varepsilon u^\{3\}=0,\\quad\\mathbf\{x\}\\in\\Omega;\\qquad u=g,\\quad\\mathbf\{x\}\\in\\Gamma,\(34\)whereΩ=\{ρ​R​\(θ\)​\(cos⁡θ,sin⁡θ\):ρ∈\[0,1\]\}\\Omega=\\\{\\rho R\(\\theta\)\(\\cos\\theta,\\sin\\theta\):\\rho\\in\[0,1\]\\\}andR⁡\(θ\)=1\+0\.2​cos⁡\(5​θ\)R\(\\theta\)=1\+0\.2\\cos\(5\\theta\), so that the boundaryΓ\\Gammais aC∞C^\{\\infty\}curve containing no corners; we fixk=2k=2and sweepε=0\.1,0\.5,1,2,3,4\\varepsilon=0\.1,0\.5,1,2,3,4to examine the influence of the nonlinearity strength, withε=4\\varepsilon=4taken as the main case\. The cubic nonlinearity is chosen both to distinguish the problem from the linear one and to make the nonlinear term depend only onuuitself\. Here the nonlinear and linear terms are of the same order: writingρnl=ε​u2/k2\\rho\_\{\\mathrm\{nl\}\}=\\varepsilon u^\{2\}/k^\{2\}, its typical value overu∈\[0,1\]u\\in\[0,1\]is about0\.50\.5, reaching1\.01\.0whenu=1u=1; in the main caseε=4=k2\\varepsilon=4=k^\{2\}, the equation can be rewritten as∇2u=k2​u​\(1−u2\)\\nabla^\{2\}u=k^\{2\}u\\,\(1\-u^\{2\}\), and the solution tends to the saturation valueu→1u\\to 1in the inner region, which is a rather strong nonlinearity\. The reference solution is generated by the finite element method implemented inscikit\-femtogether with Newton iteration; its discretization accuracy is verified against a manufactured solution and lies far below the operator accuracy being compared\.

Figure 6:Case 2: relativeL2L\_\{2\}error as a function of the nonlinearity strengthε\\varepsilon\. HK with the correction branch switched off \(Kc=0K\_\{c\}=0\) deteriorates rapidly withε\\varepsilon, followed by KernelOnet\-RBF, whereas HK with the low\-rank correction branch retained hardly degrades and surpasses RBF forε≥1\\varepsilon\\geq 1\.Figure 7:Case 2: \(a\) comparison of the kernelψ​\(r\)=φθ​\(r\)\\psi\(r\)=\\varphi\_\{\\theta\}\(r\)learned by KernelOnet\-RBF \(red\) with the analytic fundamental solutionΦ⁡\(r\)=12​π​K0​\(k​r\)\\Phi\(r\)=\\frac\{1\}\{2\\pi\}K\_\{0\}\(kr\)\(blue\) atε=0\.1\\varepsilon=0\.1, where the red dot isψ⁡\(0\)\\psi\(0\); \(b\) the learned kernelsφθ​\(r\)\\varphi\_\{\\theta\}\(r\)atε=1,2,3,4\\varepsilon=1,2,3,4, whose shapes are similar while their amplitudes increase with the nonlinearity strength\.![Refer to caption](https://arxiv.org/html/2609.35938v1/star_case2_nl_fields.png)Figure 8:Case 2: \(a\) FEM reference solution; \(b\) finite element mesh \(40×20040\\times 200, central fan triangulation\); \(c\)–\(e\) predicted solutions \(left\) and absolute error fields \(right\) of DeepONet, KernelOnet\-RBF, and KernelOnet\-HK, where the titles of the error plots give the relativeL2L\_\{2\}error for this sample\.Figure[6](https://arxiv.org/html/2609.35938#S3.F6)gives the variation of the relativeL2L\_\{2\}error of each method with the nonlinearity strengthε\\varepsilon, where "HK \(Kc=0K\_\{c\}=0\)" denotes the baseline with the low\-rank correction branch switched off and only the analytic fundamental\-solution branch retained; to use interior solution labels as the other methods do, this baseline is trained in a supervised manner, and its error is precisely the "irreducible model bias" described in Section 2\.5\. Asε\\varepsilonincreases, HK without the correction branch \(Kc=0K\_\{c\}=0\) deteriorates rapidly, its error growing from1\.77×10−31\.77\\times 10^\{\-3\}to2\.08×10−22\.08\\times 10^\{\-2\}\(ε=4\\varepsilon=4, about1212times\); the fully data\-driven KernelOnet\-RBF also degrades markedly \(1\.57→5\.11×10−31\.57\\to 5\.11\\times 10^\{\-3\}, about3\.33\.3times\); whereas KernelOnet\-HK, which retains the low\-rank correction branch, increases only slowly from1\.81×10−31\.81\\times 10^\{\-3\}to2\.75×10−32\.75\\times 10^\{\-3\}\(about1\.51\.5times\)\. Fitting the error as an approximate power law inε\\varepsilon\(err∼εp\\mathrm\{err\}\\sim\\varepsilon^\{p\}\) givesp≈0\.68p\\approx 0\.68\(Kc=0K\_\{c\}=0\),0\.320\.32\(RBF\), and0\.110\.11\(HK\), i\.e\., the error of HK grows most slowly with the nonlinearity\. Correspondingly, the advantage of HK over RBF widens monotonically withε\\varepsilon: RBF is slightly better forε≤0\.5\\varepsilon\\leq 0\.5, the two are close atε=1\\varepsilon=1, after which HK overtakes it and pulls ahead, and atε=4\\varepsilon=4HK is about1\.91\.9times better than RBF and about7\.67\.6times better thanKc=0K\_\{c\}=0; by contrast, the error of DeepONet stays near5×10−35\\times 10^\{\-3\}throughout the sweep and does not improve withε\\varepsilon\. The mechanism behind this trend is that the analytic fundamental\-solution branchΦ=K0​\(k​r\)/\(2​π\)\\Phi=K\_\{0\}\(kr\)/\(2\\pi\)provides the correct structure of the homogeneous field and the near\-boundary behavior for free, so the network only needs to learn a low\-dimensional source\-term correction, whereas RBF must learn both the homogeneous field and the nonlinear source term from data and struggles increasingly as the nonlinearity grows\. The amplitude of the correction branch also grows in step withε\\varepsilon, quantitatively confirming that it carries the nonlinear component: the fraction of the prediction’sL2L\_\{2\}energy that it accounts for increases from about0\.3%0\.3\\%atε=0\.1\\varepsilon=0\.1to about26\.6%26\.6\\%atε=4\\varepsilon=4, and the ratio of coefficient norms‖c‖/‖b‖\\\|c\\\|/\\\|b\\\|increases from0\.110\.11to0\.200\.20; both consistently indicate that the correction branch accommodates the anharmonic component introduced by the cubic nonlinear termε​u3\\varepsilon u^\{3\}\. It should be pointed out thatε≥2\\varepsilon\\geq 2already exceeds the sufficient uniqueness condition given by the maximum principle \(ε<k2/\(3​umax2\)≈1\.33\\varepsilon<k^\{2\}/\(3u\_\{\\max\}^\{2\}\)\\approx 1\.33\); however, using Newton iteration from two different initial guesses, we obtained the same solution forε=2,3,4\\varepsilon=2,3,4\(maximum difference≤10−13\\leq 10^\{\-13\}\), indicating that the solution is unique and can be solved robustly over the range considered\.

Figure[7](https://arxiv.org/html/2609.35938#S3.F7)examines the relationship between the learned kernel and the fundamental solution of this linear principal part\. Figure[7](https://arxiv.org/html/2609.35938#S3.F7)\(a\) compares, atε=0\.1\\varepsilon=0\.1\(weak nonlinearity\), the kernelψ​\(r\)=φθ​\(r\)\\psi\(r\)=\\varphi\_\{\\theta\}\(r\)learned by KernelOnet\-RBF pointwise with the fundamental solutionΦ⁡\(r\)=12​π​K0​\(k​r\)\\Phi\(r\)=\\frac\{1\}\{2\\pi\}K\_\{0\}\(kr\)of the modified Helmholtz equation: the two are qualitatively similar only near the origin—ψ\\psitakes a finite value at the origin \(ψ⁡\(0\)=1\.342\\psi\(0\)=1\.342, whereasΦ\\Phidiverges there\) and decreases monotonically withrr; but asrrincreases the two separate noticeably,ψ\\psidecaying more slowly thanΦ\\Phi, crossing zero atr≈0\.6r\\approx 0\.6and becoming negative, reaching about−0\.96\-0\.96atr=2r=2, wherer=2r=2already exceeds the maximum radius of the star\-shaped domain and corresponds to the value of the kernel function itself, whereasΦ\\Phiis always positive and tends rapidly to zero asrrincreases\. It can thus be seen that in this problem the learned kernel does not coincide with an affine transformation ofΦ\\Phi, in contrast to Case 1 \(the Laplace equation\), where the kernel almost recovers an affine image of the fundamental solution: for problems with a nonlinear source term, the single data\-driven radial kernel carries information about both the homogeneous field and the source term in its fitting\. This contrast also delineates the scope of applicability of the statement that "the learned kernel is the fundamental solution": it stems from the property that in constant\-coefficient linear problems the solution can be spanned by an expansion in fundamental solutions, and is therefore a conclusion for linear problems; for nonlinear problems, no usable fundamental solution exists for the equation, and the learned kernel can no longer be interpreted as a fundamental solution, see Section 2\.4\.1\. Figure[7](https://arxiv.org/html/2609.35938#S3.F7)\(b\) further gives the learned kernels atε=1,2,3,4\\varepsilon=1,2,3,4: the four have similar shapes, but their amplitudes increase systematically with the nonlinearity strength,ψ⁡\(0\)\\psi\(0\)being about1\.341\.34atε=1\\varepsilon=1and increasing to2\.012\.01atε=4\\varepsilon=4, i\.e\., the network adapts the scale of the kernel to cope with a stronger nonlinear source term\. The network evaluation part of HK only needs to be computed forKc=32K\_\{c\}=32correction centers, whereas RBF requires computation over allB=200B=200centers, which is a direct manifestation of the cost difference between the two\.

Figure[8](https://arxiv.org/html/2609.35938#S3.F8)gives the field distribution of the main case \(k=2k=2,ε=4\\varepsilon=4\) on a representative test sample: Figure[8](https://arxiv.org/html/2609.35938#S3.F8)\(a\) is the FEM reference solution, \(b\) is the finite element mesh \(40×20040\\times 200, central fan triangulation\), and \(c\)–\(e\) are in turn the predicted solutions \(left\) and absolute error fields \(right\) of DeepONet, KernelOnet\-RBF, and KernelOnet\-HK\. All three methods reproduce the overall morphology of the reference solution, and the differences in accuracy are mainly reflected in the resolution of the boundary neighborhood: the absolute error of DeepONet is distributed fairly uniformly over the domain, with a magnitude comparable to that of RBF \(L2=4\.77×10−3L\_\{2\}=4\.77\\times 10^\{\-3\}on this sample\), indicating that a fully connected trunk lacking physical priors struggles to capture the rapid variation near the boundary caused by the strong nonlinearity; KernelOnet\-RBF hasL2=4\.80×10−3L\_\{2\}=4\.80\\times 10^\{\-3\}on this sample, but its error is clearly concentrated near the five convex extrema of the star\-shaped domain and is more confined to the thin boundary layer; the error field of KernelOnet\-HK has the lowest overall magnitude \(2\.77×10−32\.77\\times 10^\{\-3\}\), and the peaks at the five convex extrema are also lower than those of RBF, being more uniformly distributed over the whole domain\. In the test\-set average sense, the relativeL2L\_\{2\}error of HK is2\.75×10−32\.75\\times 10^\{\-3\}, better than5\.11×10−35\.11\\times 10^\{\-3\}for RBF and5\.71×10−35\.71\\times 10^\{\-3\}for DeepONet\.

### 3\.3Case 3: complex Helmholtz equation in an unbounded exterior domain

The third case study considers the two\-dimensional complex Helmholtz equation in an unbounded exterior domain, which is used to test the ability of KernelOnet to solve problems on infinite domains and to capture traveling\-wave propagation:

∇2u\+k2​u=0,𝐱∈Ω=\{𝐱∈ℝ2:\|𝐱\|\>R\};u=g,𝐱∈Γ=∂Ω,\\nabla^\{2\}u\+k^\{2\}u=0,\\quad\\mathbf\{x\}\\in\\Omega=\\\{\\mathbf\{x\}\\in\\mathbb\{R\}^\{2\}:\\ \|\\mathbf\{x\}\|\>R\\\};\\qquad u=g,\\quad\\mathbf\{x\}\\in\\Gamma=\\partial\\Omega,\(35\)whereR=0\.5R=0\.5, the wavenumberkkis adjustable, and the boundary condition is of Dirichlet type; the solutionu⁡\(𝐱\)=ure​\(𝐱\)\+i​uim​\(𝐱\)u\(\\mathbf\{x\}\)=u^\{\\mathrm\{re\}\}\(\\mathbf\{x\}\)\+i\\,u^\{\\mathrm\{im\}\}\(\\mathbf\{x\}\)is complex\-valued and satisfies the Sommerfeld radiation condition at infinity \(i\.e\., only outward\-propagating waves exist\)\. This section scansk∈\{5,10,15,20\}k\\in\\\{5,10,15,20\\\}and takesk=20k=20as the main case of interest\. Here only KernelOnet\-PIKF, whose kernel is the analytic fundamental solution, and its boundary integral form KernelOnet\-PIKF\-SL are applicable: its kernelG⁡\(r\)=i4​H0\(1\)​\(k​r\)G\(r\)=\\frac\{i\}\{4\}H\_\{0\}^\{\(1\)\}\(kr\)inherently satisfies the Helmholtz equation and the Sommerfeld radiation condition, whereas DeepONet, PI\-DeepONet, and KernelOnet\-RBF all lack a mechanism to enforce the radiation condition \(see Section 2\.4\.2\)\.

Table[4](https://arxiv.org/html/2609.35938#S3.T4)gives the average relativeL2L\_\{2\}error overk∈\{5,10,15,20\}k\\in\\\{5,10,15,20\\\}for the two kernel\-basis forms, collocation and boundary integral, together with their SVD\-truncated variants\. All four configurations are trained without supervision, with only the boundary residual imposed: the collocation form moves the source points radially inward, i\.e\.,xs=γ​xbx\_\{s\}=\\gamma\\,x\_\{b\}with a learnableγ<1\\gamma<1; the boundary integral form adopts the single\-layer potential kernel basis of Section 2\.4\.2; and the SVD variants apply the orthogonal kernel\-basis preprocessing of Appendix[B](https://arxiv.org/html/2609.35938#A2)to the boundary kernel matrix\. Without truncation, the error of the collocation form varies slowly with the wavenumber at the10−310^\{\-3\}level, whereas the boundary integral form degrades markedly as the wavenumber increases; after SVD truncation, not only does the magnitude of both drop across the board, but their growth with the wavenumber also becomes the gentlest, indicating that the truncated orthogonal kernel basis is more robust as the wavenumber rises\. The truncation order is selected offline according to the reconstructability of the boundary data: for the collocation kernel,q=40q=40suffices at all wavenumbers, whereas the singular values of the single\-layer potential kernel decay slowly and the required order grows with the wavenumber, so that ifq=40q=40were taken uniformly, the integral form would truncate away the high\-order modes carrying the boundary data fork≥15k\\geq 15, degrading the error to the10−110^\{\-1\}level\. This shows that SVD preprocessing is exactly the key to eliminating the numerical ill\-conditioning of the single\-layer potential kernel and bringing the integral form to the same accuracy level as the collocation form\. After training, the source\-point scaling coefficient of the collocation form converges toγ≈0\.50\\gamma\\approx 0\.50atk=20k=20, lying inside the obstacle, so that the fundamental solution strictly satisfies the governing equation and the radiation condition in the exterior domain\.

Table 4:Case 3: average relativeL2L\_\{2\}error \(averaged over the real/imaginary parts\) of the two kernel\-basis forms, collocation and boundary integral, and of their SVD\-truncated variants, on the complex Helmholtz problem in an unbounded exterior domain, as a function of the wavenumberkk\. The SVD truncation order is selected offline according to the boundary reconstructability criterion:q=40q=40for the collocation form at all wavenumbers, andq=40,40,75,110q=40,40,75,110for the boundary integral form atk=5,10,15,20k=5,10,15,20, respectively\.Figure[9](https://arxiv.org/html/2609.35938#S3.F9)gives the predicted solution and absolute error field of KernelOnet\-PIKF \(collocation form\) on a representative test sample atk=20k=20: the top and bottom rows show the real and imaginary parts of the solution, respectively, and the three columns give in turn the MFS reference solution, the predicted solution, and the absolute error field\. The method accurately captures the oscillation and outward decay of the solution in the unbounded exterior domain: the predicted fields of the real and imaginary parts are visually almost indistinguishable from the reference solution, and the error is mainly concentrated in the near\-field region around the circular boundary and decays rapidly outward\. Since the expansion strictly satisfies the governing equation and the Sommerfeld radiation condition, the error field is likewise an outward\-propagating solution satisfying that equation and radiation condition, and its distribution is determined by the boundary fitting residual; it is therefore concentrated in the boundary neighborhood and decays toward the far field\.

![Refer to caption](https://arxiv.org/html/2609.35938v1/fig11.png)Figure 9:Case 3: predicted solution and absolute error field for a representative test sample atk=20k=20\(KernelOnet\-PIKF, collocation form\)\. The top and bottom rows show the real and imaginary parts of the solution, respectively, and the three columns give in turn the MFS reference solution, the KernelOnet\-PIKF predicted solution, and the absolute error field annotated with the relativeL2L\_\{2\}error\.
### 3\.4Case 4: underwater acoustic radiation and propagation caused by spherical\-shell vibration in a shallow\-water waveguide

As an extension from benchmark problems to engineering applications, this case applies KernelOnet to a practically oriented acoustic problem: underwater acoustic radiation and propagation caused by the vibration of a shell structure in a shallow\-water environment\. This problem is widespread in application scenarios such as structural acoustics of underwater vehicles and ocean engineering / ocean ambient noise, and the difficulty of solving it lies in the fact that, as acoustic waves propagate in the shallow\-water waveguide, they undergo multiple reflections from the sea surface and the seafloor, forming a complex waveguide interference structure; therefore, a dedicated Green’s function capable of characterizing the waveguide reflection effects must be used, rather than the free\-space fundamental solution; at the same time, the computational domain is usually unbounded, and both the complex near\-field acoustic field of the structure and the traveling\-wave propagation in the far field must be captured simultaneously\. This scenario is a direct continuation of the unbounded\-exterior\-domain traveling\-wave problem of Case 3, but the physical model is closer to engineering practice\. The problem background and model setup of this case follow the study of Fu et al\. on shell acoustic radiation in a shallow\-water waveguide[Fu et al\. \(2020\)](https://arxiv.org/html/2609.35938#bib.bib53)\.

As shown in Figure[10](https://arxiv.org/html/2609.35938#S3.F10), consider a thin\-walled spherical shell fully immersed in a shallow sea, with radiusR=0\.5R=0\.5m, sphere center at depthh=15h=15m, and sea depthH=25H=25m, seawater densityρ1=1025\\rho\_\{1\}=1025kg/m3, seafloor sediment densityρ2=2600\\rho\_\{2\}=2600kg/m3, and sound speedc2=1620c\_\{2\}=1620m/s\. Rectangular coordinates\(x,y,z\)\(x,y,z\)are adopted, with the origin at the sea surface and thezzaxis pointing vertically upward; the sea surfacez=0z=0is a pressure\-release boundary and the seafloorz=−Hz=\-His a penetrable boundary; the center of the spherical shell is located atz=−hz=\-h\. Since the spherical\-shell structure and its surface vibration are axisymmetric about thezzaxis, the acoustic fieldppis independent of the azimuthal angleη\\eta, and the three\-dimensional model can be reduced to a two\-dimensional axisymmetric problem \(consistent with the axisymmetric treatment of Fu et al\.[Fu et al\. \(2020\)](https://arxiv.org/html/2609.35938#bib.bib53)\)\. Collocation points\{θj\}\\\{\\theta\_\{j\}\\\}are taken uniformly along the two\-dimensional axisymmetric shell boundary \(θj\\theta\_\{j\}being the polar angle relative to thezzaxis\), and each collocation point corresponds to one azimuthal ring on the three\-dimensional shell surface; the source points are arranged along the same meridian\. The acoustic propagation frequency is set tof=240f=240Hz, corresponding to the angular frequencyω=2​π​f\\omega=2\\pi f\. It should be noted that the sound speed in a real ocean is a depth\-dependent sound speed profile; however, since the depth variation of the sound speed has little influence on near\-field underwater acoustic propagation \(negligible compared with the far field\), an isospeed approximation is adopted in the near\-field computation of this case, with a reference sound speedc1=1510c\_\{1\}=1510m/s and the corresponding wavenumberk=ω/c1k=\\omega/c\_\{1\}\. In the far field, the depth\-dependent sound speed profilec1​\(z\)c\_\{1\}\(z\)is instead incorporated, and three profiles are examined together,c1​\(z\)=1507−0\.24​zc\_\{1\}\(z\)=1507\-0\.24z,c1​\(z\)=1510c\_\{1\}\(z\)=1510, andc1​\(z\)=1513\+0\.24​zc\_\{1\}\(z\)=1513\+0\.24z, corresponding respectively to near\-bottom\-accelerating, uniform, and near\-bottom\-decelerating shallow\-water waveguides; under different profiles, the modes and eigenvalues of the normal modes differ, and the required propagating modes and the waveguide interference structure change accordingly\. The acoustic field satisfies the Helmholtz equation in the frequency domain

∇2p\+k2​p=0,\\nabla^\{2\}p\+k^\{2\}p=0,\(36\)whereppis the complex acoustic pressure: its modulus\|p\|\|p\|represents the amplitude of the acoustic\-pressure oscillation at that frequency, and its argumentarg⁡p\\arg prepresents the phase of the oscillation\. Since, in a linear frequency\-domain acoustic field,ppdiffers from the velocity potential only by a constant factor,ppthereby also carries the meaning of the velocity potential\. For ease of presentation, the sound pressure level \(SPL\) is taken asSPL=20​log10⁡\(\|p\|/10−6\)\\mathrm\{SPL\}=20\\log\_\{10\}\(\|p\|/10^\{\-6\}\)dB, i\.e\., with10−610^\{\-6\}Pa as the reference pressure and the decibel as the unit for the acoustic\-pressure amplitude; the acoustic\-pressure distributions in the figures of this section are displayed in this quantity\.

![Refer to caption](https://arxiv.org/html/2609.35938v1/figures/fig15.png)Figure 10:Case 4: schematic of the numerical model for shell structural acoustic radiation in a shallow\-water waveguide \(not to scale\)\. The origin of the rectangular coordinates\(x,y,z\)\(x,y,z\)lies at the sea surface, with thezzaxis pointing vertically upward;Ωo\\Omega\_\{o\}is the seawater domain \(densityρ1\\rho\_\{1\}, sound speedc1c\_\{1\}\); the sea surface is a pressure\-release boundaryΓ1\\Gamma\_\{1\}, and at a depthHHbelow it lies a penetrable boundaryΓ2\\Gamma\_\{2\}, beneath which is the seafloor sediment layer \(densityρ2\\rho\_\{2\}, sound speedc2c\_\{2\}\)\. This case takes a spherical\-shell radiusR=0\.5R=0\.5m, sphere\-center depthh=15h=15m, sea depthH=25H=25m, and frequencyf=240f=240Hz; the sea\-surface reflection coefficient isa2=−1a\_\{2\}=\-1and the seafloor reflection coefficient isa1=0\.4626a\_\{1\}=0\.4626\.![Refer to caption](https://arxiv.org/html/2609.35938v1/fig16.png)Figure 11:Case 4 \(near\-field acoustic radiation in a shallow\-water waveguide,x∈\[2,50\]x\\in\[2,50\]m,z∈\[−24,−1\]z\\in\[\-24,\-1\]m\): reference solution \(left\), unsupervised KernelOnet\-PIKF prediction \(middle\), and relative amplitude error \(right\) for a representative test sample\. The left and middle columns show the sound pressure level SPL=20​log10⁡\(\|p\|/10−6\)=20\\log\_\{10\}\(\|p\|/10^\{\-6\}\)\(dB\), and the right column shows the relative amplitude errorε~=\|\|ppred\|−\|pref\|\|/\|pref\|\\tilde\{\\varepsilon\}=\\bigl\|\|p\_\{\\mathrm\{pred\}\}\|\-\|p\_\{\\mathrm\{ref\}\}\|\\bigr\|/\|p\_\{\\mathrm\{ref\}\}\|\. The shallow\-water waveguide interference structure formed by multiple reflections at the sea surface and the seafloor is clearly visible, and the prediction is visually almost indistinguishable from the reference solution\.The kernel functions of the near and far fields are taken as the shallow\-water waveguide Green’s functionGnG^\{n\}and the normal\-mode Green’s functionGfG^\{f\}, respectively:GnG^\{n\}explicitly incorporates the waveguide reflection effects through an infinite mirror\-image superposition of the sea surface and the seafloor and is applicable to the near field for incidence angles below40∘40^\{\\circ\};GfG^\{f\}is formed by superposing the normal modes determined by the sound speed profile and is applicable to the far field\. Their definitions, the incidence\-angle dependence of the reflection coefficients, and the normal\-mode eigenvalue problem are given in Appendix[C](https://arxiv.org/html/2609.35938#A3)\. This case constructs the kernels separately for the near and far fields and solves them separately: the near\-field kernel is summed overNηN\_\{\\eta\}virtual nodes on the ring in the form of a ring source, whereas the far\-field kernel is evaluated as a point source since the field distance is far larger than the shell size; the reference acoustic fields useGnG^\{n\}andGfG^\{f\}as kernels, respectively, and are synthesized after the virtual source strengths are obtained by least squares from the same Dirichlet data on the shell\-surface collocation points\. The operator takes the complex acoustic pressurep⁡\(θj\)p\(\\theta\_\{j\}\)on the shell\-surface collocation points as input and the acoustic field as output, i\.e\., it is the solution operator from the shell\-surface complex acoustic pressure to the acoustic field\. Since the far\-field kernel is determined by the sound speed profilec1​\(z\)c\_\{1\}\(z\), a dataset is generated for each of the three profiles and one operator model is trained for each, so as to examine the adaptability of the operator to different waveguides\. The source points of both the near\- and far\-field kernels are taken on concentric virtual spheres obtained by contracting the shell\-surface collocation points radially toward the sphere center; the kernel matrices are numerically ill\-conditioned under this source\-point arrangement, so the boundary kernel matrices are preprocessed with the SVD orthogonal kernel basis according to Appendix[B](https://arxiv.org/html/2609.35938#A2), withq=40q=40for the near field; for the far field, because the magnitudes of the kernel values on the evaluation points and of the boundary kernel values differ drastically, the retained directions are instead selected according to the input subspace \(taking the orthogonal basis of the subspace spanned by the shell\-surface vibration modes, of dimensionJ=4J=4; the criterion is given in Appendix[B](https://arxiv.org/html/2609.35938#A2)\), and the complex coefficients are output by the branch network\.

Table[5](https://arxiv.org/html/2609.35938#S3.T5)summarizes the relativeL2L\_\{2\}errors of the various configurations on the test set together with three acoustic quality metrics\. All three metrics are defined in terms of the complex acoustic pressure and are computed sample by sample and then averaged over the20002000test samples: the relative amplitude error is\|\|ppred\|−\|pref\|\|/\|pref\|\\bigl\|\|p\_\{\\mathrm\{pred\}\}\|\-\|p\_\{\\mathrm\{ref\}\}\|\\bigr\|/\|p\_\{\\mathrm\{ref\}\}\|, the phase error is\|arg⁡\(ppred/pref\)\|\\bigl\|\\arg\(p\_\{\\mathrm\{pred\}\}/p\_\{\\mathrm\{ref\}\}\)\\bigr\|, and the transmission loss deviation is the root mean square ofTLpred−TLref\\mathrm\{TL\}\_\{\\mathrm\{pred\}\}\-\\mathrm\{TL\}\_\{\\mathrm\{ref\}\}, whereTL=−20​log10⁡\(\|p⁡\(r,z\)\|/\|p⁡\(r0,z0\)\|\)\\mathrm\{TL\}=\-20\\log\_\{10\}\\bigl\(\|p\(r,z\)\|/\|p\(r\_\{0\},z\_\{0\}\)\|\\bigr\)and the reference point\(r0,z0\)\(r\_\{0\},z\_\{0\}\)is taken at the source depth at the closest distance within the window concerned; the first two are averaged over the evaluation points\. The statistical windows are consistent with the corresponding acoustic\-field figures, i\.e\.,x∈\[2,50\]x\\in\[2,50\]m for the near field andx∈\[950,1000\]x\\in\[950,1000\]m for the far field, and only the interior of the water column is counted, excluding one row each at the sea surface and the seafloor, so as to avoid meaningless relative errors wherep≡0p\\equiv 0on the pressure\-release boundary\. For the near field, the unsupervised KernelOnet\-PIKF \(collocation form,q=40q=40\) attains an accuracy of8\.57×10−48\.57\\times 10^\{\-4\}; as in Cases 1 and 3, this shows that when the kernel function is physically consistent, high accuracy can be obtained without any interior\-solution labels\. For the far field, the normal\-mode kernel is evaluated as a point source and contains only33propagating modes, and its shell\-surface boundary matrix is more ill\-conditioned than the near\-field ring\-source kernel \(with a condition number of101710^\{17\}\); if truncated according to the firstqqsingular directions of the boundary kernel matrix, the projection error on the shell surface is only2\.4×10−42\.4\\times 10^\{\-4\}, but the amplification factor of the discarded directions on the far field,wk=∥Gi​𝐯k∥2/σkw\_\{k\}=\\lVert G\_\{i\}\\mathbf\{v\}\_\{k\}\\rVert\_\{2\}/\\sigma\_\{k\}, is as high as103∼10810^\{3\}\\sim 10^\{8\}, and the far\-field truncation error reaches0\.610\.61, consistent with the truncation floor, indicating that its root cause is the representational capability of the kernel basis rather than insufficient optimization, and that no amount of training can get past it \(measured value0\.6090\.609\)\. The far field therefore instead selects the retained directions according to the input subspace: the input is a complex linear combination ofJ=4J=4shell\-surface vibration modes, whose spanning subspace has orthogonal basisΨd\\Psi\_\{d\}; the relative amplification factor on this subspace is only about22, and the boundary residual again becomes a stable surrogate for the interior field\. In this setting, the unsupervised KernelOnet\-PIKF attains1\.18×10−31\.18\\times 10^\{\-3\},1\.31×10−31\.31\\times 10^\{\-3\}, and1\.18×10−31\.18\\times 10^\{\-3\}under the three sound speed profiles, respectively, without using any interior acoustic\-field labels during training, showing that the operator with the normal\-mode Green’s function as its kernel is applicable to different sound speed profiles\. The three quality metrics corroborate the relativeL2L\_\{2\}error: the relative amplitude error is7\.2×10−4∼1\.0×10−37\.2\\times 10^\{\-4\}\\sim 1\.0\\times 10^\{\-3\}, comparable in magnitude to the relativeL2L\_\{2\}error, indicating that the prediction error comes mainly from the amplitude rather than the phase; the phase error does not exceed4\.1×10−24\.1\\times 10^\{\-2\}degrees, which, converted with the near\-field wavenumberk=ω/c1≈1\.0k=\\omega/c\_\{1\}\\approx 1\.0rad/m, corresponds to an equivalent acoustic path error of less than11mm, far smaller than the wavelength of about6\.36\.3m; and the transmission loss deviations are all below10−210^\{\-2\}dB, far smaller than the variation scale of the transmission loss of the reference solution itself, which decays by23\.023\.0dB within the near\-field window and fluctuates by about2\.52\.5dB within the far\-field window, see Table[D6](https://arxiv.org/html/2609.35938#A4.T6)\(see Appendix[D](https://arxiv.org/html/2609.35938#A4)\)\. In addition, under the three sound speed profiles, the maximum relative differences in the relativeL2L\_\{2\}error, the relative amplitude error, and the phase error all do not exceed20%20\\%, and although the transmission loss deviation varies relatively more, its absolute value is likewise below10−210^\{\-2\}dB, further indicating that the operator is insensitive to changes in the sound speed profile\.

Table 5:Case 4: training results with the unsupervised KernelOnet\-PIKF for both the near field and the far field under the various sound speed profiles\. The near field uses the collocation form with truncation orderq=40q=40, and the far field selects the retained directions of the kernel basis according to the input subspace; the relativeL2L\_\{2\}error is reported separately for the real and imaginary parts and then averaged, and the definitions and statistical conventions of the relative amplitude error, the phase error, and the transmission loss deviation are given in the main text\.Figure[11](https://arxiv.org/html/2609.35938#S3.F11)gives, for a random test sample in the near field, the reference solution, the predicted solution, and the relative amplitude error of the unsupervised KernelOnet\-PIKF; the left and middle columns use the sound pressure level SPL \(dB\) color scale, so as to reveal simultaneously the strong field near the shell surface and the weak field in the far region\. It can be seen that the predicted solution and the reference solution are visually almost indistinguishable: the acoustic field is strongest near the spherical shell, decays rapidly outward along the radial direction, and exhibits oblique modal interference fringes formed by multiple reflections at the sea surface and the seafloor, which is exactly the typical feature distinguishing shallow\-water waveguide propagation from free\-space propagation\. The error field is mainly concentrated in the near\-field region around the spherical shell and decays rapidly outward\. In the error field, the relative deviation is most pronounced at the waveguide interference nodes, i\.e\., where\|p\|\|p\|is small: a node is where the multiple paths of sea\-surface and seafloor reflections cancel each other coherently, and its residual amplitude is determined by the tiny imbalance among the complex phases of the various paths, making it extremely sensitive to the complex\-phase error of the kernel functionGnG^\{n\}—with even a slight deviation in the complex phase, the waves that should cancel cannot cancel strictly, and the node positions shift accordingly\.

![Refer to caption](https://arxiv.org/html/2609.35938v1/fig17.png)Figure 12:Case 4 \(far\-field acoustic propagation using the normal\-mode kernel,x∈\[950,1000\]x\\in\[950,1000\]m,z∈\[−24,−1\]z\\in\[\-24,\-1\]m\): reference solution \(left column\), unsupervised KernelOnet\-PIKF prediction \(middle column\), and relative amplitude error \(right column\) for a representative test sample, corresponding to three shallow\-water sound speed profilesc1​\(z\)c\_\{1\}\(z\)\. Herezzis the vertical coordinate \(sea surfacez=0z=0, seafloorz=−Hz=\-H\); the left and middle columns show the sound pressure level SPL \(dB\), and the right column shows the relative amplitude errorε~=\|\|ppred\|−\|pref\|\|/\|pref\|\\tilde\{\\varepsilon\}=\\bigl\|\|p\_\{\\mathrm\{pred\}\}\|\-\|p\_\{\\mathrm\{ref\}\}\|\\bigr\|/\|p\_\{\\mathrm\{ref\}\}\|\. The three profiles: \(a\)c1​\(z\)=1507−0\.24​zc\_\{1\}\(z\)=1507\-0\.24z, \(b\)c1​\(z\)=1510c\_\{1\}\(z\)=1510, and \(c\)c1​\(z\)=1513\+0\.24​zc\_\{1\}\(z\)=1513\+0\.24z\. Under all three profiles, the operator accurately reproduces the far\-field waveguide modal interference structure \(nodes and antinodes along the depth direction\); the figure shows a single test sample, and the average errors over the test set are given in Table[5](https://arxiv.org/html/2609.35938#S3.T5)\.Figure[12](https://arxiv.org/html/2609.35938#S3.F12)shows the reference solution, the unsupervised KernelOnet\-PIKF predicted solution, and the relative amplitude error in the far field under the three sound speed profiles, at distancesx∈\[950,1000\]x\\in\[950,1000\]m\. Each row corresponds to one sound speed profilec1​\(z\)c\_\{1\}\(z\), wherezzis the vertical coordinate \(sea surfacez=0z=0, seafloorz=−Hz=\-H\), in turnc1​\(z\)=1507−0\.24​zc\_\{1\}\(z\)=1507\-0\.24z,c1​\(z\)=1510c\_\{1\}\(z\)=1510, andc1​\(z\)=1513\+0\.24​zc\_\{1\}\(z\)=1513\+0\.24z\. The far\-field acoustic pressure exhibits a clear normal\-mode standing\-wave structure along the depth direction, and the reference solution and the predicted solution agree closely, indicating that the operator has successfully learned the coherent\-superposition law of the normal modes\. Similar to the near field, the relative amplitude deviation in the error field is most pronounced near the waveguide interference nodes \(where\|p\|\|p\|is small\), which is caused by the coherent cancellation of the multiple paths at a node and by the fact that the residual amplitude is determined by the imbalance among the complex phases of the various paths\. At the same time, comparing the positions of the modal standing waves in the three rows, one can observe that the node/antinode positions shift as the sound speed profile changes, indicating that the far\-field waveguide modal structure is sensitive to the sound speed profile; this observation is consistent with the conclusions of Fu et al\. Although the modal structure changes markedly with the profile, the operator can still stably characterize the far\-field propagation in the three waveguides, further verifying the ability of KernelOnet\-PIKF with the normal\-mode Green’s function as its kernel to characterize the modal propagation properties of the shallow\-water far\-field waveguide\. As a reference, Table[D6](https://arxiv.org/html/2609.35938#A4.T6)\(see Appendix[D](https://arxiv.org/html/2609.35938#A4)\) gives the transmission loss of the reference solution itself: in the near field it decreases by23\.023\.0dB from22m outward at26\.826\.8m, close to the22\.522\.5dB of spherical spreading; within the far\-field window it fluctuates by only about2\.52\.5dB, whereas cylindrical spreading over the same distance brings only about0\.20\.2dB of loss, showing that the far\-field structure is dominated by modal interference\.

## 4Conclusion and Future Work

This paper has proposed and systematically verified the Kernel Operator Network \(KernelOnet\), an interpretable operator learning framework that embeds kernel functions explicitly into a neural operator\. By replacing the implicit trunk network of DeepONet with explicit kernel functions, KernelOnet is mathematically consistent with the kernel expansion of meshless collocation methods, and three complementary kernel construction strategies are given, namely the data\-driven learnable kernel \(KernelOnet\-RBF\), the physics\-informed kernel \(KernelOnet\-PIKF\) and the hybrid kernel \(KernelOnet\-HK\), which together form a spectrum from a strong physical prior to fully data\-driven learning; the accuracy basis of the three can also be given respectively by reproducing kernel space theory and by the spectral convergence of the fundamental\-solution expansion\. In the four examples—the Laplace equation on a circular domain, the nonlinear modified Helmholtz equation on a star\-shaped domain, the complex Helmholtz equation in an unbounded exterior domain, and the underwater acoustic radiation and propagation caused by spherical\-shell vibration in a shallow\-water waveguide—KernelOnet achieved high\-accuracy solutions in all cases; on the two examples that can be compared directly with DeepONet it attains higher accuracy with fewer learnable parameters\. Its advantages can be summarized along three main lines\. The first is physical consistency\. The physics\-informed kernel hard\-codes the governing equation into the network structure, so that the expansion satisfies the governing equation exactly, can therefore be trained without supervision and without interior solution labels, and dispenses with the cost of repeatedly evaluating the PDE residual by automatic differentiation; this is confirmed by all configurations of Examples 3 and 4 as well as by KernelOnet\-PIKF in Example 1\. For unbounded exterior domains and travelling\-wave propagation problems, KernelOnet\-PIKF, which takes the analytic fundamental solution as its kernel, is also the only one of the operator learning methods compared in this paper that is directly applicable\. The second is interpretability\. The analytic fundamental\-solution kernel allows the homogeneous part of the expansion to be compared pointwise with the analytic fundamental solution, while the few coefficients and kernel parameters of the correction branch give physical insight into the nonlinear source term: in Example 1 the data\-driven kernel function is almost recovered as an affine image of the fundamental solution, whereas in Example 2 theL2L\_\{2\}energy proportion of the correction branch grows from0\.3%0\.3\\%to26\.6%26\.6\\%with the nonlinearity, directly quantifying the weight of the nonlinear component in the solution\. The third is computational efficiency and engineering applicability\. For problems that admit an analytic fundamental solution, the physics\-informed kernel and the hybrid kernel attain high accuracy without interior solution labels; for strongly nonlinear and variable\-coefficient problems that admit none, the data\-driven KernelOnet\-RBF compensates for the breadth of applicability by adaptively learning the kernel shape, at the price of a weaker physical prior\. In the shallow\-water waveguide example, both the near\-field ring\-source kernel and the far\-field normal\-mode kernel are numerically ill\-conditioned, and once the retained directions of the kernel basis are selected from the input subspace both the near field and the far field can be trained without supervision and remain robust under three sound speed profiles\. In addition, the inference cost of the operator is far below that of per\-instance solvers, and in nonlinear problems and scenarios with expensive boundary integrals the training cost is recovered after a few hundred to two thousand queries; once training is complete, the operator can be evaluated directly on evaluation grids of arbitrary resolution without retraining\. These results show that combining explicit kernel functions with physical priors achieves a good compromise among accuracy, interpretability, computational efficiency and breadth of applicability\.

Several limitations remain in this work, and they indicate corresponding directions for improvement\. First, the construction of the physics\-informed kernel has not yet been automated: the analytic fundamental solution or Green’s function must be derived by hand for each specific physical setting—Example 4, for instance, required the separate construction of a near\-field mirror\-image superposition kernel and a far\-field normal\-mode kernel for the shallow\-water waveguide, together with a separate delimitation of their respective ranges of validity\. A feasible improvement is to build a library of kernel functions for common operators and media and to automate derivations such as the method of images and separation of variables; for media that cannot be solved analytically, one may precompute numerical Green’s functions, or retain the structure of the hybrid kernel and let the low\-rank learned branch absorb the bias of the analytic kernel\. Second, when extending to high frequencies and large scales the size of the kernel basis grows accordingly: the order of the SVD truncation introduced to suppress ill\-conditioning rises markedly with the wavenumber—in Example 3 the single\-layer potential kernel grows from4040atk=5k=5to110110atk=20k=20, and the far field of Example 4 even requires a switch to a criterion that selects the directions from the input subspace—which shows that the selection of the retained directions still relies on case\-by\-case analysis\. The corresponding improvement is to incorporate the amplification factorwkw\_\{k\}of the evaluation block into an automatic order\-selection criterion, and to combine preconditioning, block low\-rank factorization and multiscale kernel bases layered by frequency band, so as to decouple the size of the kernel basis from the wavenumber\. Third, the theoretical guarantee of the overall approximation is still incomplete: the error bounds of Section 2\.5 presuppose the positive definiteness of the learned kernel, and the spectral convergence of the fundamental\-solution expansion is an asymptotic conclusion; the two have not yet been combined into a convergence rate and generalization error bound that covers both the kernel expansion and the learning of the coefficients by the branch network\. One may, within the native space framework, regard the operator approximation as a restricted approximation problem, decompose the error according to the size of the kernel basis and the number of training samples, and give an explicit rate for the correction branch of the hybrid kernel under a compressibility assumption on the source term\.

## CRediT authorship contribution statement

Yuan Guo: Writing – review & editing, Writing – original draft, Visualization, Validation, Methodology, Investigation, Formal analysis, Data curation, Conceptualization\. Hanshu Chen: Writing – review & editing, Software, Validation\. Qiang Xi: Writing – review & editing, Software, Validation\. Timon Rabczuk: Supervision, Writing – review & editing\. Zhuojia Fu: Writing – review & editing, Methodology, Supervision, Project administration, Funding acquisition\.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper\.

## Data availability

## Acknowledgments

The research was supported by the National Natural Science Foundation of China \(12122205, 12372196\)\.

## Appendix ACommon fundamental solutions of differential operators

The KernelOnet\-PIKF presented in Section 2\.4 of this paper uses the analytic fundamental solution of the governing equation as the kernel function, whose applicability presupposes that the governing equation admits an explicit fundamental solution\. To facilitate the direct construction of KernelOnet\-PIKF for different problems, this appendix collects the fundamental solutions of several common linear differential operators in two dimensions \(2D\) and three dimensions \(3D\); the complete listing of the corresponding harmonic functions and radial Trefftz functions can be found in the appendix of the physics\-informed kernel function neural network \(PIKFNN\)[Fu et al\. \(2024\)](https://arxiv.org/html/2609.35938#bib.bib52)\.

Letℒ\\mathcal\{L\}be a constant\-coefficient linear partial differential operator,𝐱∈ℝd\\mathbf\{x\}\\in\\mathbb\{R\}^\{d\}, the radial distance between the field point𝐱\\mathbf\{x\}and the source point𝐱s\\mathbf\{x\}\_\{s\}ber=\|𝐱−𝐱s\|r=\|\\mathbf\{x\}\-\\mathbf\{x\}\_\{s\}\|, and the fundamental solutionΦ\\Phisatisfy

ℒ​Φ​\(𝐱,𝐱s\)=−δ⁡\(𝐱−𝐱s\)\.\\mathcal\{L\}\\,\\Phi\(\\mathbf\{x\},\\mathbf\{x\}\_\{s\}\)=\-\\delta\(\\mathbf\{x\}\-\\mathbf\{x\}\_\{s\}\)\.\(37\)Table[A1](https://arxiv.org/html/2609.35938#A1.T1)lists the fundamental solutions of several common operatorsℒ\\mathcal\{L\}in 2D and 3D, whereΔ\\Deltais the Laplace operator,kkis the wavenumber \(or decay coefficient\),DDis the diffusion coefficient,𝐯\\mathbf\{v\}is the convection velocity,μ=k2/D\+\|𝐯\|2/\(4​D2\)\\mu=\\sqrt\{k^\{2\}/D\+\|\\mathbf\{v\}\|^\{2\}/\(4D^\{2\}\)\},H0\(1\)H\_\{0\}^\{\(1\)\}is the zero\-order Hankel function of the first kind, andK0K\_\{0\},J0J\_\{0\},I0I\_\{0\}are the zero\-order modified Bessel function of the second kind, the zero\-order Bessel function of the first kind, and the zero\-order modified Bessel function of the first kind, respectively\. These fundamental solutions satisfy the corresponding homogeneous governing equation \(ℒ​Φ=0\\mathcal\{L\}\\Phi=0,r≠0r\\neq 0\), and can therefore be directly used as the kernel backbone of KernelOnet\-PIKF\.

Table A1:Fundamental solutionsΦ\\Phiof common differential operators \(2D and 3D\)\.It should be emphasized that the above fundamental solutions \(together with their corresponding higher\-order variants and the fundamental solutions of time\-dependent operators, such as those of the heat operator∂/∂t−k​Δ\\partial/\\partial t\-k\\Deltaand the wave operator∂2/∂t2−c12Δ\\partial^\{2\}/\\partial t^\{2\}\-c\_\{1\}^\{2\}\\Delta\) can all be constructed as PIKF in the manner described in Section 2\.3 and Section 2\.4\. If the governing equation of a certain class of problems itself admits no analytic fundamental solution \(such as nonlinear equations with nonlinear source terms\), then the strict version of KernelOnet\-PIKF is no longer applicable; in this case, one may construct the hybrid kernel KernelOnet\-HK using the fundamental solution of the linear principal part of the governing equation, or instead fully adopt the data\-driven KernelOnet\-RBF; the specific trade\-off is discussed in Section 2\.4\.

## Appendix BSVD Orthogonalization of the Kernel Basis

When the kernel function family is numerically ill\-conditioned under a given discretization of source points and evaluation points, using it directly as the kernel severely pollutes the gradients\. To this end, this paper applies an orthogonalization preprocessing of the kernel basis, based on the singular value decomposition \(SVD\), to the kernel matrixG⁡\(𝐱\)G\(\\mathbf\{x\}\): for the collocation form \(Eq\. \([20](https://arxiv.org/html/2609.35938#S2.E20)\)\) and the boundary integral form \(Eq\. \([21](https://arxiv.org/html/2609.35938#S2.E21)\)\), the basis function familyψj​\(𝐱\)\\psi\_\{j\}\(\\mathbf\{x\}\)consists of linear combinations \(pointwise or boundary\-integral\) of fundamental solutions, so the preprocessing steps are identical for the two, and only the construction ofG⁡\(𝐱\)G\(\\mathbf\{x\}\)differs\.

For the collocation form \(Eq\. \([20](https://arxiv.org/html/2609.35938#S2.E20)\)\), the basis functions areψj​\(𝐱\)=Φ⁡\(\|𝐱−γ​𝐱b\(j\)\|\)\\psi\_\{j\}\(\\mathbf\{x\}\)=\\Phi\\bigl\(\\left\|\\mathbf\{x\}\-\\gamma\\,\\mathbf\{x\}\_\{b\}^\{\(j\)\}\\right\|\\bigr\)\. The kernel matrixG⁡\(𝐱\)G\(\\mathbf\{x\}\)is constructed from theBBbasis functions as follows: at a single evaluation point𝐱\\mathbf\{x\},𝐠⁡\(𝐱\)=\[ψ1​\(𝐱\),…,ψB​\(𝐱\)\]∈ℂ1×B\\mathbf\{g\}\(\\mathbf\{x\}\)=\\bigl\[\\psi\_\{1\}\(\\mathbf\{x\}\),\\ldots,\\psi\_\{B\}\(\\mathbf\{x\}\)\\bigr\]\\in\\mathbb\{C\}^\{1\\times B\}is the row vector of the kernel matrix at that point, so that thejj\-th column of the kernel matrix is the basis functionψj​\(𝐱\)\\psi\_\{j\}\(\\mathbf\{x\}\); stacking these rows at thentn\_\{t\}evaluation points\{𝐱m\}m=1nt\\\{\\mathbf\{x\}\_\{m\}\\\}\_\{m=1\}^\{n\_\{t\}\}yields the discrete kernel matrixG⁡\(𝐱\)∈ℂnt×BG\(\\mathbf\{x\}\)\\in\\mathbb\{C\}^\{n\_\{t\}\\times B\}with entriesGm​j=ψj​\(𝐱m\)G\_\{mj\}=\\psi\_\{j\}\(\\mathbf\{x\}\_\{m\}\)\. A singular value decomposition is then performed on this kernel matrix

G⁡\(𝐱\)=U⁡\(𝐱\)​Σ​\(𝐱\)​VH​\(𝐱\),Σ⁡\(𝐱\)=diag⁡\(σ1​\(𝐱\),σ2​\(𝐱\),…,σr​\(𝐱\)\),σ1​\(𝐱\)≥σ2​\(𝐱\)≥⋯≥σr​\(𝐱\)\>0,G\(\\mathbf\{x\}\)=U\(\\mathbf\{x\}\)\\,\\Sigma\(\\mathbf\{x\}\)\\,V^\{H\}\(\\mathbf\{x\}\),\\qquad\\Sigma\(\\mathbf\{x\}\)=\\mathrm\{diag\}\(\\sigma\_\{1\}\(\\mathbf\{x\}\),\\sigma\_\{2\}\(\\mathbf\{x\}\),\\ldots,\\sigma\_\{r\}\(\\mathbf\{x\}\)\),\\qquad\\sigma\_\{1\}\(\\mathbf\{x\}\)\\geq\\sigma\_\{2\}\(\\mathbf\{x\}\)\\geq\\cdots\\geq\\sigma\_\{r\}\(\\mathbf\{x\}\)\>0,\(38\)
wherer≤min⁡\(nt,B\)r\\leq\\min\(n\_\{t\},B\)is the rank of the matrix, andU⁡\(𝐱\)=\[𝐮1​\(𝐱\),…,𝐮r​\(𝐱\)\]U\(\\mathbf\{x\}\)=\[\\mathbf\{u\}\_\{1\}\(\\mathbf\{x\}\),\\ldots,\\mathbf\{u\}\_\{r\}\(\\mathbf\{x\}\)\]andV⁡\(𝐱\)=\[𝐯1​\(𝐱\),…,𝐯r​\(𝐱\)\]V\(\\mathbf\{x\}\)=\[\\mathbf\{v\}\_\{1\}\(\\mathbf\{x\}\),\\ldots,\\mathbf\{v\}\_\{r\}\(\\mathbf\{x\}\)\]are the left and right singular vectors, respectively: the left singular vector𝐮k​\(𝐱\)\\mathbf\{u\}\_\{k\}\(\\mathbf\{x\}\)is defined on the evaluation points and is a function of spatial position, while the right singular vector𝐯k​\(𝐱\)\\mathbf\{v\}\_\{k\}\(\\mathbf\{x\}\)is likewise a function of spatial position and varies with the evaluation point set;\(⋅\)H\(\\cdot\)^\{H\}denotes the conjugate transpose\. When the kernel function family is approximately linearly dependent under the source\-point arrangement, the singular values decay sharply from some point onward, the numerical rank is far smaller than the number of collocation points, and the condition numbercond⁡\(G⁡\(𝐱\)\)=σ1​\(𝐱\)/σr​\(𝐱\)\\mathrm\{cond\}\(G\(\\mathbf\{x\}\)\)=\\sigma\_\{1\}\(\\mathbf\{x\}\)/\\sigma\_\{r\}\(\\mathbf\{x\}\)reaches as high as101610^\{16\}, whereupon the small singular directions amplify tiny perturbations of the boundary fitting residual into huge fluctuations of the expansion coefficients and make the gradient computation severely ill\-conditioned\. Taking the firstqqprincipal singular directions forms a well\-conditioned orthogonal kernel basis

Uq\(𝐱\)=U\(𝐱\)\[:,1:q\],U\_\{q\}\(\\mathbf\{x\}\)=U\(\\mathbf\{x\}\)\\,\[:,1:q\],\(39\)
The selection of retained directions cannot rest solely on the magnitude of the singular values of the boundary kernel matrix; the amplification of the evaluation block must also be taken into account\. LetG\(b\)​\(𝐱\)G^\{\(b\)\}\(\\mathbf\{x\}\)denote the boundary kernel matrix andG⁡\(𝐱\)G\(\\mathbf\{x\}\)the kernel matrix at the evaluation points, with𝐯k\\mathbf\{v\}\_\{k\}the right singular vector ofG\(b\)G^\{\(b\)\}; then the amplification factor of thekk\-th direction over the evaluation block iswk=∥G⁡\(𝐱\)​𝐯k∥2/σk\(b\)w\_\{k\}=\\lVert G\(\\mathbf\{x\}\)\\mathbf\{v\}\_\{k\}\\rVert\_\{2\}/\\sigma\_\{k\}^\{\(b\)\}\. Since the interior field error caused by the boundary fitting residual𝜹\\boldsymbol\{\\delta\}is exactly∥G⁡\(𝐱\)​𝜹∥2\\lVert G\(\\mathbf\{x\}\)\\boldsymbol\{\\delta\}\\rVert\_\{2\}, directions for whichσk\(b\)\\sigma\_\{k\}^\{\(b\)\}is very small whilewkw\_\{k\}is very large must be retained or discarded as a whole: truncating merely by the magnitude of the singular values discards their contributions to the interior field along with them\. When the boundary block and the evaluation block are of comparable magnitude \(Cases 1–3\),wkw\_\{k\}is bounded and taking the firstqqprincipal singular directions suffices; when the two differ greatly \(the normal\-mode kernel of the far field in Case 4, wherewkw\_\{k\}reaches103∼10810^\{3\}\\sim 10^\{8\}\), the retained directions should instead be selected according to the input subspace: let the input data span aJJ\-dimensional subspace with complex orthonormal basisΨd\\Psi\_\{d\}, and take

Ψb=Ψd,Ψi=G⁡\(𝐱\)​pinv​\(G\(b\)​\(𝐱\)\)​Ψd,\\Psi\_\{b\}=\\Psi\_\{d\},\\qquad\\Psi\_\{i\}=G\(\\mathbf\{x\}\)\\,\\mathrm\{pinv\}\\bigl\(G^\{\(b\)\}\(\\mathbf\{x\}\)\\bigr\)\\,\\Psi\_\{d\},\(40\)
that is,Ψb\\Psi\_\{b\}andΨi\\Psi\_\{i\}serve as the kernel bases of the boundary block and the evaluation block, respectively\. In this case the boundary residual always lies within the input subspace, its amplification to the interior field is characterized byΨi​Ψb\+\\Psi\_\{i\}\\Psi\_\{b\}^\{\+\}and is bounded in norm, and the boundary residual again becomes a stable surrogate for the interior field\. When the input is a linear combination of known modes,Ψd\\Psi\_\{d\}can be given analytically; otherwise it can be estimated offline by the proper orthogonal decomposition \(POD\) of the training boundary data matrix\. Neither relies on interior solution labels\.

that is, the columns𝐮1​\(𝐱\),…,𝐮q​\(𝐱\)\\mathbf\{u\}\_\{1\}\(\\mathbf\{x\}\),\\ldots,\\mathbf\{u\}\_\{q\}\(\\mathbf\{x\}\)ofUq​\(𝐱\)U\_\{q\}\(\\mathbf\{x\}\)are taken as the kernel basis, and this kernel basis is well\-conditioned\. The operator output is taken as

𝒢θ​\(a\)​\(𝐱\)=Uq​\(𝐱\)​𝐚,𝐚=ℬθ​\(a\)∈ℂq,\\mathcal\{G\}\_\{\\theta\}\(a\)\(\\mathbf\{x\}\)=U\_\{q\}\(\\mathbf\{x\}\)\\,\\mathbf\{a\},\\qquad\\mathbf\{a\}=\\mathcal\{B\}\_\{\\theta\}\(a\)\\in\\mathbb\{C\}^\{q\},\(41\)
where the complex coefficients𝐚\\mathbf\{a\}are learned by the branch network: since the network outputs real numbers only, each complex coefficient is formed by two outputs, its real part and its imaginary part, so the output dimension of the branch network is2​q2q, and the real and imaginary parts of the solution are given by the cross inner products of the real and imaginary parts of the coefficients with the real and imaginary parts ofUq​\(𝐱\)U\_\{q\}\(\\mathbf\{x\}\)\.

For the boundary integral form \(Eq\. \([21](https://arxiv.org/html/2609.35938#S2.E21)\)\), the basis functions are integrals of the fundamental solution over boundary patches,ψj​\(𝐱\)=∫ΓjΦ⁡\(\|𝐱−𝐲\|\)​d​Γ𝐲\\psi\_\{j\}\(\\mathbf\{x\}\)=\\int\_\{\\Gamma\_\{j\}\}\\Phi\\bigl\(\\left\|\\mathbf\{x\}\-\\mathbf\{y\}\\right\|\\bigr\)\\,\\mathrm\{d\}\\Gamma\_\{\\mathbf\{y\}\}, i\.e\., flattened boundary\-element kernels\. The kernel matrixG⁡\(𝐱\)G\(\\mathbf\{x\}\)is constructed in the same way and subjected to the same singular value decomposition and truncation, yielding the kernel basisUq​\(𝐱\)U\_\{q\}\(\\mathbf\{x\}\)and the operator output in the same form\. Hence, under the SVD preprocessing the collocation form and the boundary integral form are unified into the operator structure𝒢θ​\(a\)​\(𝐱\)=Uq​\(𝐱\)​𝐚\\mathcal\{G\}\_\{\\theta\}\(a\)\(\\mathbf\{x\}\)=U\_\{q\}\(\\mathbf\{x\}\)\\,\\mathbf\{a\}, the only difference being the construction of the kernel matrixG⁡\(𝐱\)G\(\\mathbf\{x\}\)\(pointwise fundamental solutions or boundary integral kernels\)\.

The accuracy of this truncation can be quantified by the singular spectrum\. Suppose that the reference solution lies in the reachable space of the kernel function, i\.e\., there exists a coefficient vector𝜶\\boldsymbol\{\\alpha\}such thatp⁡\(𝐱\)=G⁡\(𝐱\)​𝜶p\(\\mathbf\{x\}\)=G\(\\mathbf\{x\}\)\\,\\boldsymbol\{\\alpha\}; writingGq​\(𝐱\)=Uq​\(𝐱\)​Σq​\(𝐱\)​VqH​\(𝐱\)G\_\{q\}\(\\mathbf\{x\}\)=U\_\{q\}\(\\mathbf\{x\}\)\\,\\Sigma\_\{q\}\(\\mathbf\{x\}\)\\,V\_\{q\}^\{H\}\(\\mathbf\{x\}\), the projection error ofpponto the truncated subspacespan​\(Uq​\(𝐱\)\)\\mathrm\{span\}\(U\_\{q\}\(\\mathbf\{x\}\)\)satisfies

‖p−Uq​\(𝐱\)​UqH​\(𝐱\)​p‖2=‖\(I−Uq​\(𝐱\)​UqH​\(𝐱\)\)​\(G⁡\(𝐱\)−Gq​\(𝐱\)\)​𝜶‖2≤σq\+1​\(𝐱\)​‖𝜶‖2,\\left\\\|p\-U\_\{q\}\(\\mathbf\{x\}\)U\_\{q\}^\{H\}\(\\mathbf\{x\}\)p\\right\\\|\_\{2\}=\\left\\\|\(I\-U\_\{q\}\(\\mathbf\{x\}\)U\_\{q\}^\{H\}\(\\mathbf\{x\}\)\)\(G\(\\mathbf\{x\}\)\-G\_\{q\}\(\\mathbf\{x\}\)\)\\boldsymbol\{\\alpha\}\\right\\\|\_\{2\}\\leq\\sigma\_\{q\+1\}\(\\mathbf\{x\}\)\\left\\\|\\boldsymbol\{\\alpha\}\\right\\\|\_\{2\},\(42\)that is, the truncation error is controlled jointly by the largest discarded singular valueσq\+1​\(𝐱\)\\sigma\_\{q\+1\}\(\\mathbf\{x\}\)and the norm of the expansion coefficients\. It should be emphasized that this estimate is in the boundary metric: the truncation error of the evaluation block \(the interior field\) is∥G⁡\(𝐱\)​\(I−Uq​\(𝐱\)​UqH​\(𝐱\)\)​𝜶∥2\\lVert G\(\\mathbf\{x\}\)\(I\-U\_\{q\}\(\\mathbf\{x\}\)U\_\{q\}^\{H\}\(\\mathbf\{x\}\)\)\\boldsymbol\{\\alpha\}\\rVert\_\{2\}, which is controlled by the amplification factorwkw\_\{k\}rather than byσk\\sigma\_\{k\}; a smallσq\+1/σ1\\sigma\_\{q\+1\}/\\sigma\_\{1\}therefore does not imply that the interior field remains accurate after truncation—the far field of Case 4 is a counterexample, where the retained directions must instead be selected according to the input subspace\. The ratio of the sum of squares of the firstqqsingular values to the total sum of squares,∑k=1qσk2​\(𝐱\)/∑k=1rσk2​\(𝐱\)≈100%\\sum\_\{k=1\}^\{q\}\\sigma\_\{k\}^\{2\}\(\\mathbf\{x\}\)\\big/\\sum\_\{k=1\}^\{r\}\\sigma\_\{k\}^\{2\}\(\\mathbf\{x\}\)\\approx 100\\%, is by itself insufficient to assert that the truncation is lossless; the premise that “the reference solution is indeed a reachable vector of this kernel expansion” is also required\. By the Eckart–Young theorem, the error of the optimalqq\-rank approximation ofG⁡\(𝐱\)G\(\\mathbf\{x\}\)in the spectral norm is exactly the largest discarded singular value

‖G⁡\(𝐱\)−Uq​\(𝐱\)​Σq​\(𝐱\)​VqH​\(𝐱\)‖2=σq\+1​\(𝐱\),\\left\\\|G\(\\mathbf\{x\}\)\-U\_\{q\}\(\\mathbf\{x\}\)\\,\\Sigma\_\{q\}\(\\mathbf\{x\}\)\\,V\_\{q\}^\{H\}\(\\mathbf\{x\}\)\\right\\\|\_\{2\}=\\sigma\_\{q\+1\}\(\\mathbf\{x\}\),\(43\)and whenσq\+1​\(𝐱\)\\sigma\_\{q\+1\}\(\\mathbf\{x\}\)has fallen to machine precision relative toσ1​\(𝐱\)\\sigma\_\{1\}\(\\mathbf\{x\}\), the truncation is numerically exact\.

It should be emphasized that the improvement in the condition number comes from switching to the orthogonal basis of left singular vectors rather than from the truncation itself: the columns ofU⁡\(𝐱\)U\(\\mathbf\{x\}\)are themselves mutually orthogonal and of unit norm, socond⁡\(U⁡\(𝐱\)\)=1\\mathrm\{cond\}\(U\(\\mathbf\{x\}\)\)=1; the additional role of the truncation is to constrain the approximation space to the physically meaningful principal singular subspace, removing the spurious modes associated with small singular values that are numerically unreliable upon extrapolation, and thereby eliminating their ill\-conditioned amplification of the gradients\.

When the kernel matrix is fixed during training \(as in the boundary integral form, or when the source\-point positions are fixed\), this preprocessing depends only on the kernel function and the discretization layout and is independent of the sample data; it can be completed offline once before training and introduces no learnable parameters\. If the source\-point positions take part in learning, thenG⁡\(𝐱\)G\(\\mathbf\{x\}\)changes accordingly and the singular value decomposition must be recomputed, so that it is no longer a one\-off offline operation\.

It should be emphasized that the SVD truncation does not destroy the physical consistency of the kernel\. FromUq​\(𝐱\)=G⁡\(𝐱\)​Vq​\(𝐱\)​Σq​\(𝐱\)−1U\_\{q\}\(\\mathbf\{x\}\)=G\(\\mathbf\{x\}\)\\,V\_\{q\}\(\\mathbf\{x\}\)\\,\\Sigma\_\{q\}\(\\mathbf\{x\}\)^\{\-1\}it follows that every column ofUq​\(𝐱\)U\_\{q\}\(\\mathbf\{x\}\)is a column vector of the kernel matrixG⁡\(𝐱\)G\(\\mathbf\{x\}\), i\.e\., a linear combination of the basis functions \(fundamental solutions or their boundary integrals\); since the fundamental solutions satisfy the governing equation exactly and linear superposition does not change the nature of the solution, each column still satisfies the homogeneous governing equation and the corresponding boundary conditions, so that the expansionUq​\(𝐱\)​𝐚U\_\{q\}\(\\mathbf\{x\}\)\\,\\mathbf\{a\}automatically satisfies the governing equation for arbitrary coefficients and the property that the PDE residual vanishes identically is preserved\. Therefore, training remains an unsupervised boundary residual fit, enforcing the boundary conditions only at the boundary collocation points

𝒥⁡\(θ\)=1N​B​∑i=1N∑k=1B\|\(Uq\(b\)​\(𝐱\)​𝐚i\)k−yb,k\(i\)\|2,\\mathcal\{J\}\(\\theta\)=\\frac\{1\}\{N\\,B\}\\sum\_\{i=1\}^\{N\}\\sum\_\{k=1\}^\{B\}\\left\|\\bigl\(U\_\{q\}^\{\(b\)\}\(\\mathbf\{x\}\)\\,\\mathbf\{a\}\_\{i\}\\bigr\)\_\{k\}\-y\_\{b,k\}^\{\(i\)\}\\right\|^\{2\},\(44\)
whereUq\(b\)​\(𝐱\)U\_\{q\}^\{\(b\)\}\(\\mathbf\{x\}\)is the value of the kernel basisUq​\(𝐱\)U\_\{q\}\(\\mathbf\{x\}\)at the boundary collocation points, i\.e\., the rows of the kernel matrixG⁡\(𝐱\)G\(\\mathbf\{x\}\)corresponding to the boundary collocation points after truncation, and𝐚i=ℬθ​\(𝐲b\(i\)\)\\mathbf\{a\}\_\{i\}=\\mathcal\{B\}\_\{\\theta\}\\bigl\(\\mathbf\{y\}\_\{b\}^\{\(i\)\}\\bigr\)is the expansion coefficient output by the branch network for theii\-th sample\.

Finally, the coefficient representationUq​\(𝐱\)=G⁡\(𝐱\)​Vq​\(𝐱\)​Σq​\(𝐱\)−1U\_\{q\}\(\\mathbf\{x\}\)=G\(\\mathbf\{x\}\)\\,V\_\{q\}\(\\mathbf\{x\}\)\\,\\Sigma\_\{q\}\(\\mathbf\{x\}\)^\{\-1\}also reveals the functional nature of the truncated basis and provides the theoretical basis for the discrete invariance discussed above: sinceUq​\(𝐱\)⊂col⁡\(G⁡\(𝐱\)\)U\_\{q\}\(\\mathbf\{x\}\)\\subset\\mathrm\{col\}\(G\(\\mathbf\{x\}\)\), the truncated basis can be extended to a continuous basis through the combination coefficientsC⁡\(𝐱\)=Vq​\(𝐱\)​Σq​\(𝐱\)−1∈ℂB×qC\(\\mathbf\{x\}\)=V\_\{q\}\(\\mathbf\{x\}\)\\,\\Sigma\_\{q\}\(\\mathbf\{x\}\)^\{\-1\}\\in\\mathbb\{C\}^\{B\\times q\}—thekk\-th column of the kernel basis,𝐮k​\(𝐱\)=∑j=1BCj​k​\(𝐱\)​ψj​\(𝐱\)\\mathbf\{u\}\_\{k\}\(\\mathbf\{x\}\)=\\sum\_\{j=1\}^\{B\}C\_\{jk\}\(\\mathbf\{x\}\)\\,\\psi\_\{j\}\(\\mathbf\{x\}\), is a linear combination of the basis functions defined at an arbitrary spatial point, and on the training gridG⁡\(𝐱\)​C​\(𝐱\)=Uq​\(𝐱\)G\(\\mathbf\{x\}\)\\,C\(\\mathbf\{x\}\)=U\_\{q\}\(\\mathbf\{x\}\)holds exactly elementwise, so that the extension is not an approximation but an exact recovery of the functional nature ofUq​\(𝐱\)U\_\{q\}\(\\mathbf\{x\}\)\. Therefore the basis can be evaluated at any new evaluation point viaUq​\(𝐱\)=G⁡\(𝐱\)​C​\(𝐱\)U\_\{q\}\(\\mathbf\{x\}\)=G\(\\mathbf\{x\}\)\\,C\(\\mathbf\{x\}\), and changing the resolution of the evaluation grid requires no retraining, so that the discrete invariance of the operator is preserved\. In Case 4, both the near\-field Pekeris kernel and the far\-field normal\-mode kernel are numerically ill\-conditioned under the source\-point arrangement, and the above preprocessing is therefore adopted: for the near fieldq=40q=40is taken \(determined by the criterionσq\+1​\(𝐱\)/σ1​\(𝐱\)≤5×10−4\\sigma\_\{q\+1\}\(\\mathbf\{x\}\)/\\sigma\_\{1\}\(\\mathbf\{x\}\)\\leq 5\\times 10^\{\-4\}\); although the far\-field normal\-mode kernel is likewise ill\-conditioned, its small singular directions carry the dominant energy of the far field \(σ50/σ1∼10−16\\sigma\_\{50\}/\\sigma\_\{1\}\\sim 10^\{\-16\}whilew50∼108w\_\{50\}\\sim 10^\{8\}\), and truncating to the firstqqsingular directions gives a far\-field error of0\.610\.61, so the retained directions are instead selected according to the input subspace \(the orthonormal basis of the subspace spanned by the shell\-surface vibration modes, of dimensionJ=4J=4\)\. See Section 3\.4 for details\.

## Appendix CGreen’s Functions for the Shallow\-Water Waveguide

This section gives the complete construction of the near\-field and far\-field kernel functions in Case 4\.

### C\.1Near field: Pekeris waveguide Green’s function

Since the field points in the near field are not far from the spherical shell, the acoustic field is determined jointly by the superposition of the reflected waves and the direct wave; the kernel function is therefore taken to be the shallow\-water waveguide Green’s functionGnG^\{n\}, which is generalized from the fundamental solution of Case 3: the standard Helmholtz fundamental solution is replaced by the waveguide Green’s function obtained from the infinite image superposition over the sea surface and the seafloor, so that the reflection effects are explicitly incorporated into the kernel function\. Here the sea surface is a pressure\-release boundary with reflection coefficienta2=−1a\_\{2\}=\-1; the seafloor is a penetrable boundary with reflection coefficienta1=0\.4626a\_\{1\}=0\.4626\. This reflection coefficient is determined by the incidence angleθ\\theta, and is computed as

a1=\{a​cos⁡θ−b2−sin2⁡θa​cos⁡θ\+b2−sin2⁡θ,\|sin⁡θ\|<b,1,\|sin⁡θ\|≥b,a\_\{1\}=\\begin\{cases\}\\dfrac\{a\\cos\\theta\-\\sqrt\{b^\{2\}\-\\sin^\{2\}\\theta\}\}\{a\\cos\\theta\+\\sqrt\{b^\{2\}\-\\sin^\{2\}\\theta\}\},&\|\\sin\\theta\|<b,\\\\\[6\.0pt\] 1,&\|\\sin\\theta\|\\geq b,\\end\{cases\}\(45\)wherea=ρ2/ρ1a=\\rho\_\{2\}/\\rho\_\{1\},b=c1/c2b=c\_\{1\}/c\_\{2\}, andθ\\thetais the incidence angle\. Substituting the parameters of this case givesa=2\.537a=2\.537andb=0\.932b=0\.932, witha1=0\.4626a\_\{1\}=0\.4626atθ=0\\theta=0; this coefficient varies very little as long as the incidence angle is below40∘40^\{\\circ\}, so it can be approximated as the constanta1=0\.4626a\_\{1\}=0\.4626in the near field\. This also delineates the range of validity of the simplified Pekeris waveguide Green’s functionGnG^\{n\}: incidence angles below40∘40^\{\\circ\}, i\.e\., the horizontal distance between the field point and the source point must not be too large, so thatGnG^\{n\}applies only to the near field\. In this way,GnG^\{n\}inherently satisfies the Helmholtz equation and the reflecting boundaries of the shallow\-water waveguide\. The source point is mapped about thezzaxis intoNηN\_\{\\eta\}circumferential virtual nodes and the kernel function is summed along the circumferential direction; its explicit expression is

Gn​\(x,y,z,x0,y0,z0\)=∑ε=0Nη−1∑λ=0∞\(a1​a2\)λ​\(e−i​k​R1\(ε\)R1\(ε\)\+a1​e−i​k​R2\(ε\)R2\(ε\)\+a2​e−i​k​R3\(ε\)R3\(ε\)\+a1​a2​e−i​k​R4\(ε\)R4\(ε\)\),G^\{n\}\(x,y,z;x\_\{0\},y\_\{0\},z\_\{0\}\)=\\sum\_\{\\varepsilon=0\}^\{N\_\{\\eta\}\-1\}\\sum\_\{\\lambda=0\}^\{\\infty\}\(a\_\{1\}a\_\{2\}\)^\{\\lambda\}\\left\(\\frac\{e^\{\-ikR\_\{1\}^\{\(\\varepsilon\)\}\}\}\{R\_\{1\}^\{\(\\varepsilon\)\}\}\+a\_\{1\}\\frac\{e^\{\-ikR\_\{2\}^\{\(\\varepsilon\)\}\}\}\{R\_\{2\}^\{\(\\varepsilon\)\}\}\+a\_\{2\}\\frac\{e^\{\-ikR\_\{3\}^\{\(\\varepsilon\)\}\}\}\{R\_\{3\}^\{\(\\varepsilon\)\}\}\+a\_\{1\}a\_\{2\}\\frac\{e^\{\-ikR\_\{4\}^\{\(\\varepsilon\)\}\}\}\{R\_\{4\}^\{\(\\varepsilon\)\}\}\\right\),\(46\)where the four propagation paths and the circumferential virtual nodes are

\{R1\(ε\)=\(x−xε\)2\+\(y−yε\)2\+\(2​λ​H\+z−zε\)2,R2\(ε\)=\(x−xε\)2\+\(y−yε\)2\+\(2​λ​H\+2​\(H−h\)\+z\+zε\)2,R3\(ε\)=\(x−xε\)2\+\(y−yε\)2\+\(2​λ​H\+2​h−z−zε\)2,R4\(ε\)=\(x−xε\)2\+\(y−yε\)2\+\(2​\(λ\+1\)​H−z\+zε\)2,xε=x02\+y02cosηε,yε=x02\+y02sinηε,zε=z0,ηε=2​π​εNη,\\left\\\{\\begin\{aligned\} R\_\{1\}^\{\(\\varepsilon\)\}&=\\sqrt\{\(x\-x\_\{\\varepsilon\}\)^\{2\}\+\(y\-y\_\{\\varepsilon\}\)^\{2\}\+\(2\\lambda H\+z\-z\_\{\\varepsilon\}\)^\{2\}\},\\\\ R\_\{2\}^\{\(\\varepsilon\)\}&=\\sqrt\{\(x\-x\_\{\\varepsilon\}\)^\{2\}\+\(y\-y\_\{\\varepsilon\}\)^\{2\}\+\(2\\lambda H\+2\(H\-h\)\+z\+z\_\{\\varepsilon\}\)^\{2\}\},\\\\ R\_\{3\}^\{\(\\varepsilon\)\}&=\\sqrt\{\(x\-x\_\{\\varepsilon\}\)^\{2\}\+\(y\-y\_\{\\varepsilon\}\)^\{2\}\+\(2\\lambda H\+2h\-z\-z\_\{\\varepsilon\}\)^\{2\}\},\\\\ R\_\{4\}^\{\(\\varepsilon\)\}&=\\sqrt\{\(x\-x\_\{\\varepsilon\}\)^\{2\}\+\(y\-y\_\{\\varepsilon\}\)^\{2\}\+\(2\(\\lambda\+1\)H\-z\+z\_\{\\varepsilon\}\)^\{2\}\},\\\\ x\_\{\\varepsilon\}&=\\sqrt\{x\_\{0\}^\{2\}\+y\_\{0\}^\{2\}\}\\cos\\eta\_\{\\varepsilon\},\\quad y\_\{\\varepsilon\}=\\sqrt\{x\_\{0\}^\{2\}\+y\_\{0\}^\{2\}\}\\sin\\eta\_\{\\varepsilon\},\\quad z\_\{\\varepsilon\}=z\_\{0\},\\quad\\eta\_\{\\varepsilon\}=\\frac\{2\\pi\\varepsilon\}\{N\_\{\\eta\}\},\\end\{aligned\}\\right\.\(47\)where\(x,y,z\)\(x,y,z\)and\(x0,y0,z0\)\(x\_\{0\},y\_\{0\},z\_\{0\}\)are the coordinates of the field point and the source point, respectively,\(xε,yε,zε\)\(x\_\{\\varepsilon\},y\_\{\\varepsilon\},z\_\{\\varepsilon\}\)is the circumferential virtual node obtained by rotating the source point about thezzaxis through the azimuthal angle2​π​ε/Nη2\\pi\\varepsilon/N\_\{\\eta\},hhandHHare the immersion depth of the spherical shell and the sea depth, respectively,kkis the near\-field wavenumber, andNηN\_\{\\eta\}is the number of circumferential nodes;R1R\_\{1\}is the direct wave,R2R\_\{2\}the single seafloor image,R3R\_\{3\}the single sea\-surface image, andR4R\_\{4\}the double seafloor–sea\-surface image, whose coefficients are11,a1a\_\{1\},a2a\_\{2\}, anda1​a2a\_\{1\}a\_\{2\}, respectively\. Asλ\\lambdaincreases, the images move farther and farther from the field point and their magnitudes decay according to\(a1​a2\)λ\(a\_\{1\}a\_\{2\}\)^\{\\lambda\}, andλ→∞\\lambda\\to\\inftycorresponds to the exact infinite\-image solution; since\|a1​a2\|=0\.4626<1\|a\_\{1\}a\_\{2\}\|=0\.4626<1, the series converges rapidly inλ\\lambda, and taking the highest order of the image superposition to beM=200M=200already makes the truncation error negligible\. KernelOnet\-PIKF usesGnG^\{n\}as its kernel, so the expansion automatically satisfies the governing equation and the waveguide boundaries, and training can therefore be completed unsupervised using only the boundary residuals on the shell\-surface collocation points\.

### C\.2Far field: normal\-mode Green’s function

In the far field the horizontal distance becomes larger, the incidence angle at the seafloor accordingly exceeds40∘40^\{\\circ\}, and the above approximation of treatinga1a\_\{1\}as a constant no longer holds, so the normal\-mode Green’s function is used instead[Fu et al\. \(2020\)](https://arxiv.org/html/2609.35938#bib.bib53)

Gf​\(x,y,z,x0,y0,z0\)=i​πρ1​∑d=1Nmϕd​\(z0\)​ϕd​\(z\)​H0\(1\)​\(μd​\(x−x0\)2\+\(y−y0\)2\),G^\{f\}\(x,y,z;x\_\{0\},y\_\{0\},z\_\{0\}\)=\\frac\{i\\pi\}\{\\rho\_\{1\}\}\\sum\_\{d=1\}^\{N\_\{m\}\}\\phi\_\{d\}\(z\_\{0\}\)\\,\\phi\_\{d\}\(z\)\\,H\_\{0\}^\{\(1\)\}\\\!\\Big\(\\mu\_\{d\}\\sqrt\{\(x\-x\_\{0\}\)^\{2\}\+\(y\-y\_\{0\}\)^\{2\}\}\\Big\),\(48\)whereH0\(1\)H\_\{0\}^\{\(1\)\}is the zeroth\-order Hankel function of the first kind andNmN\_\{m\}is the number of propagating modes; under the strictly penetrable \(Robin\) seafloor condition, all three sound speed profiles of this case haveNm=3N\_\{m\}=3propagating modes \(becausec2\>c1c\_\{2\}\>c\_\{1\}, trapped modes exist only in the narrow bandk2<μ<k1k\_\{2\}<\\mu<k\_\{1\}\), andμd\\mu\_\{d\}andϕd\\phi\_\{d\}are the eigenvalue and the normalized mode of thedd\-th normal mode, respectively\. The modesϕd\\phi\_\{d\}and the eigenvaluesμd\\mu\_\{d\}are obtained by solving the Sturm\-Liouville eigenvalue problem determined by the shallow\-water sound speed profilec1​\(z\)c\_\{1\}\(z\):

d2​ϕd​\(z\)d​z2\+\[ω2c12​\(z\)−μd2\]​ϕd​\(z\)=0,−H<z<0,\\frac\{d^\{2\}\\phi\_\{d\}\(z\)\}\{dz^\{2\}\}\+\\left\[\\frac\{\\omega^\{2\}\}\{c\_\{1\}^\{2\}\(z\)\}\-\\mu\_\{d\}^\{2\}\\right\]\\phi\_\{d\}\(z\)=0,\\qquad\-H<z<0,\(49\)and satisfy the pressure\-release boundary condition at the sea surface and the penetrable boundary condition at the seafloor

ϕd​\(0\)=0,ϕd​\(−H\)\+ρ2ρ1​1μd2−\(ω/c2\)2​d​ϕd​\(−H\)d​z=0\.\\phi\_\{d\}\(0\)=0,\\qquad\\phi\_\{d\}\(\-H\)\+\\frac\{\\rho\_\{2\}\}\{\\rho\_\{1\}\}\\frac\{1\}\{\\sqrt\{\\mu\_\{d\}^\{2\}\-\(\\omega/c\_\{2\}\)^\{2\}\}\}\\,\\frac\{d\\phi\_\{d\}\(\-H\)\}\{dz\}=0\.\(50\)whereκ=μd2−\(ω/c2\)2\\kappa=\\sqrt\{\\mu\_\{d\}^\{2\}\-\(\\omega/c\_\{2\}\)^\{2\}\}is the vertical decay wavenumber of the normal mode in the sediment layer; this condition is derived from the continuity of pressure and of normal displacement at the seafloor interface, andμd2−\(ω/c2\)2\\sqrt\{\\mu\_\{d\}^\{2\}\-\(\\omega/c\_\{2\}\)^\{2\}\}should appear in the denominator to keep the dimensions consistent\.

## Appendix DSupplementary Results for the Numerical Examples

This appendix reports supplementary results that are not expanded in the main text: the accuracy levels and H1 semi\-norm errors of the various methods, the computational cost of training and inference, the break\-even comparison with classical per\-instance solvers, the batch\-size ablation, the resolution invariance with respect to the evaluation grid, and the transmission loss of the reference solution\. All times are measured on the machine used in this paper, and the training and inference times are synchronized bytorch\.cuda\.synchronize\(\)to reflect the true GPU wall\-clock time; except in the batch\-size ablation subsection, all methods are trained with the full batch, i\.e\., every epoch traverses all20002000training samples\.

### D\.1Relative H1 semi\-norm error and comparison at equal accuracy

The relativeH1H^\{1\}semi\-norm error measures the gradient difference between the predicted and the reference solution, and is computed at the sample level:

ℰH1=1N​∑i=1N∥∇\(upred\(i\)−utrue\(i\)\)∥∥∇utrue\(i\)∥,\\mathcal\{E\}\_\{H^\{1\}\}=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\\frac\{\\bigl\\lVert\\nabla\\bigl\(u\_\{\\mathrm\{pred\}\}^\{\(i\)\}\-u\_\{\\mathrm\{true\}\}^\{\(i\)\}\\bigr\)\\bigr\\rVert\}\{\\lVert\\nabla u\_\{\\mathrm\{true\}\}^\{\(i\)\}\\rVert\},\(51\)where∥⋅∥\\lVert\\cdot\\rVertis the discreteL2L\_\{2\}norm on the evaluation grid: the gradient is approximated by second\-order central differences, the star\-shaped\-domain grid of Case 2 uses body\-fitted curvilinear coordinates and the gradient is transformed to physical space through the metric tensor of that coordinate system; the discrete summation is weighted by the area element of each evaluation grid, namelyr​d​r​d​θr\\,\{\\rm d\}r\\,\{\\rm d\}\\thetafor the polar grids of Cases 1 and 3, the analytic area elementρ​R​\(θ\)2​d​ρ​d​θ\\rho R\(\\theta\)^\{2\}\\,\{\\rm d\}\\rho\\,\{\\rm d\}\\thetafor Case 2, and a constant area element for the uniform meridional\-plane grid of Case 4, which is equivalent to equal\-weight summation\. For the complex\-valued fields of Cases 3 and 4, the gradient modulus is given by the sum of the squares of the real and imaginary parts, consistent with the treatment of the relativeL2L\_\{2\}error\.

The equal\-accuracy comparison in Table[D1](https://arxiv.org/html/2609.35938#A4.T1)shows that the kernel\-based methods converge markedly faster than DeepONet: on Cases 1 and 2 they reach10−210^\{\-2\}in less than one third of the epochs required by the latter, neither DeepONet nor PI\-DeepONet ever drops to10−310^\{\-3\}within the whole budget, whereas KernelOnet\-PIKF reaches that accuracy on Case 1 at3\.8×1053\.8\\times 10^\{5\}epochs\. The H1 semi\-norm errors give exactly the same ordering as theL2L\_\{2\}errors, which shows that the high accuracy of the kernel\-based methods is not limited to theL2L\_\{2\}norm\. The ratio of the two also varies with the type of problem: the solutions of Cases 1 and 2 are smooth and the H1 error is about6∼246\\sim 24times theL2L\_\{2\}error, indicating that the gradient is harder to fit than the function values themselves; the acoustic fields of Cases 3 and 4 are oscillatory, and the ratio of the gradient to the function value is governed by the wavenumber content of the field itself, so the two are close, namely0\.93∼1\.040\.93\\sim 1\.04for Case 3 and0\.999∼1\.0010\.999\\sim 1\.001for Case 4\. The near equality in Case 4 arises because its far\-field reference solution and the network kernel are both normal\-mode expansions and the error likewise lies in the space spanned by the same set of propagating modes, so that the gradient norm in the numerator and the denominator are scaled by the same wavenumber factor and cancel\. Moreover, several configurations in the table reach10−210^\{\-2\}already at the10410^\{4\}\-epoch level, with the corresponding training wall\-clock times given in Table[D2](https://arxiv.org/html/2609.35938#A4.T2), only seconds to tens of seconds, which shows that the5×1055\\times 10^\{5\}\-epoch budget mainly serves to give the comparison methods ample opportunity to converge rather than being necessary for the method of this paper\.

Table D1:RelativeL2L\_\{2\}and relativeH1H^\{1\}semi\-norm errors of each case and each method, together with the numbers of training epochs required for the validation error to first fall to given accuracy levels, all computed over20002000test samples; “–” means that the level was not reached within the5×1055\\times 10^\{5\}\-epoch budget\.Case 2 takesε=4\\varepsilon=4and also lists theKc=0K\_\{c\}=0baseline with the correction branch disabled; Case 3 takes the main casek=20k=20; for Case 4 the three sound speed profiles are listed for the far field\.CaseMethodrel\.L2L\_\{2\}rel\. H1to5×10−25\\times 10^\{\-2\}to10−210^\{\-2\}to5×10−35\\times 10^\{\-3\}to10−310^\{\-3\}1DeepONet1\.89×10−31\.89\\times 10^\{\-3\}3\.47×10−23\.47\\times 10^\{\-2\}9,500105,000182,500–1PI\-DeepONet3\.53×10−33\.53\\times 10^\{\-3\}4\.17×10−24\.17\\times 10^\{\-2\}21,500225,000389,000–1KernelOnet\-RBF1\.29×10−31\.29\\times 10^\{\-3\}3\.10×10−23\.10\\times 10^\{\-2\}3,00014,50040,000–1KernelOnet\-PIKF8\.89×10−48\.89\\times 10^\{\-4\}6\.30×10−36\.30\\times 10^\{\-3\}1,50010,50020,000379,5001KernelOnet\-HK2\.04×10−32\.04\\times 10^\{\-3\}2\.19×10−22\.19\\times 10^\{\-2\}3,00029,00093,000–2 \(ε=4\\varepsilon=4\)DeepONet5\.71×10−35\.71\\times 10^\{\-3\}7\.06×10−27\.06\\times 10^\{\-2\}19,500275,000––2 \(ε=4\\varepsilon=4\)KernelOnet\-RBF5\.11×10−35\.11\\times 10^\{\-3\}4\.76×10−24\.76\\times 10^\{\-2\}7,50044,000––2 \(ε=4\\varepsilon=4\)KernelOnet\-HK2\.75×10−32\.75\\times 10^\{\-3\}2\.75×10−22\.75\\times 10^\{\-2\}6,00056,000159,500–2 \(ε=4\\varepsilon=4\)KernelOnet\-HK \(Kc=0K\_\{c\}=0baseline\)2\.08×10−22\.08\\times 10^\{\-2\}1\.28×10−11\.28\\times 10^\{\-1\}8,000–––3 \(k=20k=20\)KernelOnet\-PIKF \(collocation\)1\.90×10−31\.90\\times 10^\{\-3\}1\.77×10−31\.77\\times 10^\{\-3\}2,00014,50045,000–3 \(k=20k=20\)KernelOnet\-PIKF\-SVD \(collocation\+\+SVD\)1\.36×10−31\.36\\times 10^\{\-3\}1\.37×10−31\.37\\times 10^\{\-3\}2,50010,50018,500–3 \(k=20k=20\)KernelOnet\-PIKF\-SL \(boundary integral\)5\.15×10−35\.15\\times 10^\{\-3\}5\.20×10−35\.20\\times 10^\{\-3\}19,000106,000––3 \(k=20k=20\)KernelOnet\-PIKF\-SL\-SVD \(boundary integral\+\+SVD\)1\.85×10−31\.85\\times 10^\{\-3\}1\.92×10−31\.92\\times 10^\{\-3\}2,0009,50017,500–4 \(near field\)KernelOnet\-PIKF8\.57×10−48\.57\\times 10^\{\-4\}8\.56×10−48\.56\\times 10^\{\-4\}3,50013,00021,000284,0004 \(far field\)KernelOnet\-PIKF \(c1=1507−0\.24​zc\_\{1\}\{=\}1507\{\-\}0\.24z\)1\.18×10−31\.18\\times 10^\{\-3\}1\.18×10−31\.18\\times 10^\{\-3\}3,5006,50012,500–4 \(far field\)KernelOnet\-PIKF \(c1=1510c\_\{1\}\{=\}1510\)1\.31×10−31\.31\\times 10^\{\-3\}1\.31×10−31\.31\\times 10^\{\-3\}3,5006,50012,500–4 \(far field\)KernelOnet\-PIKF \(c1=1513\+0\.24​zc\_\{1\}\{=\}1513\{\+\}0\.24z\)1\.18×10−31\.18\\times 10^\{\-3\}1\.18×10−31\.18\\times 10^\{\-3\}3,5006,50013,000–

### D\.2Model size, training cost, and GPU memory usage

Table[D2](https://arxiv.org/html/2609.35938#A4.T2)summarizes the learnable parameters, the training cost, and the GPU memory usage of the various methods\. As for the learnable parameters, the two DeepONets consist of two fully connected networks, a branch network and a trunk network, and have the largest parameter counts, about1\.8×1051\.8\\times 10^\{5\}; the kernel of KernelOnet\-PIKF has only11scalar parameter, giving the smallest parameter count,1\.03×1051\.03\\times 10^\{5\}in Case 1; KernelOnet\-RBF and KernelOnet\-HK additionally contain a shallow kernel network and a low\-rank correction network, respectively, and their parameter counts lie in between, about1\.3∼1\.5×1051\.3\\sim 1\.5\\times 10^\{5\}\. As for the training speed, the physics\-informed kernel, which needs to be evaluated only at the boundary collocation points, costs the least per epoch, namely0\.0610\.061s/100 epochs for KernelOnet\-PIKF in Case 1, slightly below the0\.0760\.076s/100 epochs of DeepONet; KernelOnet\-RBF needs to be evaluated at all160160boundary source points and reaches0\.4310\.431s/100 epochs, about seven times the former; PI\-DeepONet must evaluate the PDE residual at the interior collocation points and back\-propagate it in every epoch, costing0\.5340\.534s/100 epochs, comparable to RBF\. The peak training GPU memory is consistent with the scale of the trunk evaluation: in Case 1 KernelOnet\-PIKF uses only4848MiB, whereas KernelOnet\-RBF and PI\-DeepONet reach689689and687687MiB, respectively; apart from the332332MiB of KernelOnet\-RBF, the peak inference memory does not exceed7171MiB, the latter again arising from the need to evaluate all source points\. The differences in Case 3 come from the size of the kernel basis: the untruncated collocation form costs the most per epoch,0\.1050\.105s/100 epochs, while its boundary integral form and the SVD\-truncated variants drop to0\.063∼0\.0690\.063\\sim 0\.069s/100 epochs\. The complete training wall\-clock time and per\-sample inference time of each configuration are given in Table[D3](https://arxiv.org/html/2609.35938#A4.T3)\.

Table D2:Learnable parameters, training cost, and GPU memory usage of each case and each method\. The per\-sample training time is obtained by dividing the time per100100epochs by100×2000100\\times 2000\. Case 2 takes the nonlinear main caseε=4\\varepsilon=4and also lists theKc=0K\_\{c\}=0baseline with the correction branch disabled; Case 3 takes the main casek=20k=20; both the near field and the far field of Case 4 are trained unsupervised, and the three sound speed profiles are listed for the far field\.CaseMethodLearnable parameterss/100 ep\.Training per sample \(μ\\mus\)Memory, training \(MiB\)Memory, inference \(MiB\)1DeepONet180,8010\.0760\.38128211PI\-DeepONet180,8010\.5342\.67687211KernelOnet\-RBF129,2810\.4312\.156893321KernelOnet\-PIKF103,0410\.0610\.3048201KernelOnet\-HK132,5030\.1020\.51169592 \(ε=4\\varepsilon=4\)DeepONet187,2010\.0740\.37103202 \(ε=4\\varepsilon=4\)KernelOnet\-RBF142,1210\.4072\.046413122 \(ε=4\\varepsilon=4\)KernelOnet\-HK147,2750\.1190\.59184712 \(ε=4\\varepsilon=4\)KernelOnet\-HK \(Kc=0K\_\{c\}=0baseline\)142,1230\.0980\.49112343 \(k=20k=20\)KernelOnet\-PIKF \(collocation\)128,8010\.1050\.52104563 \(k=20k=20\)KernelOnet\-PIKF\-SVD \(collocation\+\+SVD\)90,1600\.0650\.3397433 \(k=20k=20\)KernelOnet\-PIKF\-SL \(boundary integral\)128,8000\.0690\.34158543 \(k=20k=20\)KernelOnet\-PIKF\-SL\-SVD \(boundary integral\+\+SVD\)112,7000\.0630\.3299444 \(near field\)KernelOnet\-PIKF80,5600\.0630\.3162324 \(far field\)KernelOnet\-PIKF \(c1=1507−0\.24​zc\_\{1\}\{=\}1507\{\-\}0\.24z\)68,9680\.0620\.3137244 \(far field\)KernelOnet\-PIKF \(c1=1510c\_\{1\}\{=\}1510\)68,9680\.0620\.3137244 \(far field\)KernelOnet\-PIKF \(c1=1513\+0\.24​zc\_\{1\}\{=\}1513\{\+\}0\.24z\)68,9680\.0620\.313724

### D\.3Reference\-solution generation cost and break\-even comparison

This subsection compares the cost of a classical per\-instance solver with the inference cost of the neural operator on a common basis\. The per\-instance solve is timed in two ways:*pre\-factorized*means that the instance\-independent matrix and its factorization are constructed only once, so that each instance needs only a coefficient solve and an evaluation, which is feasible for linear problems with fixed geometry;*from scratch*means that the matrix is reconstructed and factorized for every instance, in which case the geometry, the coefficients, and the source term may all vary from instance to instance\.

Table[D3](https://arxiv.org/html/2609.35938#A4.T3)shows that the benefit of break\-even depends on how expensive the classical solver is\. Case 2 is a nonlinear problem in which every instance must be reassembled and solved through Newton iterations, costing about1\.01\.0s per instance, so that the operator needs only a few hundred to two thousand queries to amortize the training cost, namely3\.6×1023\.6\\times 10^\{2\}queries for DeepONet and5\.8×1025\.8\\times 10^\{2\}for KernelOnet\-HK\. The near field of Case 4 requires the solution of a ring\-source boundary integral for the source strengths, and if the matrix is allowed to be reconstructed for every instance, the break\-even point is about6\.3×1026\.3\\times 10^\{2\}queries\. Conversely, Cases 1 and 3 and the far field of Case 4 are linear problems with fixed geometry, for which the pre\-factorized method of fundamental solutions costs only0\.005∼0\.30\.005\\sim 0\.3ms per instance, below the operator inference cost, so that no break\-even point exists in the table; here the value of the operator lies not in replacing a single solve but in obtaining the solution for arbitrary new boundary data in a single forward pass\. The stratification of the inference times is precisely what causes the differences in the break\-even points: the configurations trained directly on the boundary residual all cost0\.063∼0\.2520\.063\\sim 0\.252ms/sample, whereas KernelOnet\-RBF must evaluate all source points and reaches1\.4∼1\.51\.4\\sim 1\.5ms/sample, and hence has the highest break\-even point on Case 2\. The reference\-solution generation cost differs enormously among the cases: the finite\-element solution of Case 2 requires about4141minutes per dataset, whereas the method of fundamental solutions of Cases 1 and 3 requires only0\.150\.15and0\.340\.34s, respectively; since the unsupervised configurations need only boundary data during training, this cost does not enter their training cost but is used for evaluation and for the supervised baselines\. As an order\-of\-magnitude reference for GPU\-native solvers, the GPU\-based Galerkin finite\-element solver TensorMesh[Wen et al\. \(2026\)](https://arxiv.org/html/2609.35938#bib.bib59)takes about10∼7710\\sim 77s to solve a Helmholtz problem with on the order of one million nodes on an A100, and the cost grows with the wavenumber; such solvers are far faster than the CPU implementation of this paper for a single large\-scale solve, but the system must still be reassembled and re\-solved whenever the boundary conditions or the source term change\.

Table D3:Reference\-solution generation cost, the inference cost of the classical per\-instance solver and of the neural operator, and the break\-even number of queriesN⋆N^\{\\star\}for each case\. All times are measured on the machine used in this paper, averaged over100100instances\. Solvers: the method of fundamental solutions \(MFS\) for Cases 1 and 3, finite elements\+\+Newton for Case 2, a ring\-source boundary integral for the near field and normal modes for the far field of Case 4\. The per\-instance solve andN⋆N^\{\\star\}are both given in the pre\-factorized and the from\-scratch settings; Case 2 is a nonlinear problem whose Jacobian changes with the solution and which has no pre\-factorized setting, so only a single value is listed\. “–” means that the single\-instance cost of the classical solver is below the operator inference cost, in which case no break\-even point exists;†denotes unsupervised training that uses no interior solution labels, whose training does not depend on the reference\-solution generation cost\.
### D\.4Batch\-size ablation

Table[D4](https://arxiv.org/html/2609.35938#A4.T4)examines the influence of the batch size on the three data\-driven models in Case 2\. Four batch sizes are used,500500,10001000,15001500, and the full batch, and all training settings other than the batch size are identical to those of this case in Table[D1](https://arxiv.org/html/2609.35938#A4.T1); when the batch size is smaller than the training\-set size, every epoch traverses all training samples batch by batch, and the time per epoch increases as the batch size decreases\. All three models benefit from mini\-batches, but to different degrees: the relativeL2L\_\{2\}error of DeepONet drops from5\.71×10−35\.71\\times 10^\{\-3\}with the full batch to2\.15×10−32\.15\\times 10^\{\-3\}at a batch size of500500, an improvement of a factor of2\.662\.66, making it the most sensitive to the batch size; KernelOnet\-RBF drops from5\.11×10−35\.11\\times 10^\{\-3\}to3\.03×10−33\.03\\times 10^\{\-3\}, an improvement of1\.691\.69times; and KernelOnet\-HK drops from2\.75×10−32\.75\\times 10^\{\-3\}to1\.75×10−31\.75\\times 10^\{\-3\}, an improvement of1\.571\.57times, with the ordering of the H1 semi\-norm errors exactly the same\. Consequently, the relative merits of DeepONet and KernelOnet\-RBF change with the batch size: with the full batch RBF is slightly better, while with mini\-batches DeepONet overtakes it; KernelOnet\-HK, by contrast, remains the best at all four batch sizes\. The gain in accuracy is paid for in training time: the time per100100epochs increases by a factor of about3\.5∼3\.93\.5\\sim 3\.9, and the total wall\-clock time accordingly grows from0\.10∼0\.570\.10\\sim 0\.57h to0\.36∼2\.200\.36\\sim 2\.20h; the peak GPU memory, in contrast, decreases slightly as the batch size decreases, from103103to6161MiB for DeepONet, whereas the two kernel models change only marginally, their memory being dominated by the evaluation of the kernel function at all source points\. It can be seen that the full batch is not the training setting with the best accuracy, but since the main text adopts the full batch uniformly for the three models in every case, this does not affect the comparison among the methods; if accuracy is the sole objective, a batch size of500500can further reduce the error at about four times the training wall\-clock time\.

Table D4:Batch\-size ablation for the three data\-driven models on Case 2 \(ε=4\\varepsilon=4\): relativeL2L\_\{2\}/H1 semi\-norm errors, training time per100100epochs, total training wall\-clock time, and peak training GPU memory\. The full batch means that the batch size equals the training\-set size, and its values agree with the corresponding rows of this case in Table[D2](https://arxiv.org/html/2609.35938#A4.T2)\.
### D\.5Resolution invariance

Table[D5](https://arxiv.org/html/2609.35938#A4.T5)shows that the relativeL2L\_\{2\}errors of the three models all remain at the10−310^\{\-3\}level throughout the refinement, and that those of DeepONet and KernelOnet\-HK vary by less than10%10\\%: the solution of the proposed operator is determined by the data on the training grid, so that once training is complete it can be evaluated directly on an evaluation grid of arbitrary resolution without retraining\. The error of KernelOnet\-RBF rises to4\.3×10−34\.3\\times 10^\{\-3\}when the evaluation grid is refined to twice the training grid or more, and stays within1\.5×10−31\.5\\times 10^\{\-3\}on all the other grids\.

Table D5:Resolution invariance of the three data\-driven models with respect to the evaluation grid on Case 1: relativeL2L\_\{2\}error as a function of the evaluation gridnr×nθn\_\{r\}\\times n\_\{\\theta\}\. All models are trained on the native40×4040\\times 40grid and then evaluated directly on each grid with frozen weights; the reference solution is sampled from the same fundamental solution on each grid, and the average over20002000test samples is reported\. The number of evaluation\-grid points increases from100100to25,60025\{,\}600, a span of256256times\.
### D\.6Transmission loss of the reference solution in Case 4

Table[D6](https://arxiv.org/html/2609.35938#A4.T6)gives the transmission loss of the reference acoustic field\. The reference solutions are synthesized waveguide Green’s functions in which the modal excitation amplitudes are random for each sample, so that the absolute sound level has no physical meaning, and the transmission loss in the table therefore represents the relative sound\-level change within the window\. The near field decays rapidly outward from22m, dropping by23\.023\.0dB at26\.826\.8m, in agreement with the20​log10⁡\(r/r0\)≈22\.520\\log\_\{10\}\(r/r\_\{0\}\)\\approx 22\.5dB of spherical spreading, which indicates that the acoustic field at this distance is still dominated by spherical spreading; within the far\-field window the fluctuation is only about2\.52\.5dB, whereas cylindrical spreading yields only about0\.20\.2dB of loss over the same5050m distance, so that the far\-field structure is dominated by modal interference\. The far\-field reference point happens to lie at a relative minimum of the modal interference,9∼139\\sim 13dB below the peak of the window, so that the transmission loss at most points in the window is negative, i\.e\., stronger than at the reference point\.

Table D6:Transmission loss of the reference solution of Case 4 in each evaluation window, over20002000test samples\.TL=−20​log10⁡\(\|p⁡\(r,z\)\|/\|p⁡\(r0,z0\)\|\)\\mathrm\{TL\}=\-20\\log\_\{10\}\\bigl\(\|p\(r,z\)\|/\|p\(r\_\{0\},z\_\{0\}\)\|\\bigr\), where the reference point\(r0,z0\)\(r\_\{0\},z\_\{0\}\)is taken at the source depth at the nearest distance within the window; the TL at the source depth is the sample mean, and the TL within the window is listed as the minimum, median, and maximum over the pooled set of samples and window points\.

## References

- Aronszajn \(1950\)N\. AronszajnTheory of reproducing kernels\.Transactions of the American Mathematical Society68\(3\),pp\. 337–404\.External Links:[Document](https://dx.doi.org/10.1090/s0002-9947-1950-0051437-7)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1),[§2\.3](https://arxiv.org/html/2609.35938#S2.SS3.p2.1),[§2\.5](https://arxiv.org/html/2609.35938#S2.SS5.SSS0.Px1.p1.1)\.
- Bahmaniet al\.\(2025\)B\. Bahmani, S\. Goswami, I\. G\. Kevrekidis, and M\. D\. ShieldsA resolution independent neural operator\.Computer Methods in Applied Mechanics and Engineering444,pp\. 118113\.External Links:[Document](https://dx.doi.org/10.1016/j.cma.2025.118113)Cited by:[§2\.1](https://arxiv.org/html/2609.35938#S2.SS1.p1.3)\.
- Belytschkoet al\.\(1996\)T\. Belytschko, Y\. Krongauz, D\. Organ, M\. Fleming, and P\. KryslMeshless methods: an overview and recent developments\.Computer Methods in Applied Mechanics and Engineering139\(1\-4\),pp\. 3–47\.External Links:[Document](https://dx.doi.org/10.1016/s0045-7825%2896%2901078-x)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1)\.
- Brandstetteret al\.\(2022\)J\. Brandstetter, R\. van den Berg, M\. Welling, and J\. K\. GuptaClifford neural layers for PDE modeling\.External Links:2209\.04934Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- Brebbiaet al\.\(1984\)C\. A\. Brebbia, J\. C\. F\. Telles, and L\. C\. WrobelBoundary element techniques: theory and applications in engineering\.Springer Berlin Heidelberg,Berlin\.External Links:ISBN 9783642488627,[Document](https://dx.doi.org/10.1007/978-3-642-48860-3)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1),[§2\.3\.1](https://arxiv.org/html/2609.35938#S2.SS3.SSS1.p2.2),[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px3.p3.1),[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px4.p1.1),[§2\.4\.2](https://arxiv.org/html/2609.35938#S2.SS4.SSS2.p4.1)\.
- Buhmann \(2003\)M\. D\. BuhmannRadial basis functions: theory and implementations\.Cambridge University Press\.External Links:[Document](https://dx.doi.org/10.1017/cbo9780511543241)Cited by:[§2\.3\.1](https://arxiv.org/html/2609.35938#S2.SS3.SSS1.p1.1),[§2\.5](https://arxiv.org/html/2609.35938#S2.SS5.SSS0.Px1.p1.2)\.
- Cao \(2021\)S\. CaoChoose a transformer: Fourier or Galerkin\.External Links:2105\.14995Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- Chenet al\.\(1998\)C\. S\. Chen, M\. A\. Golberg, and Y\. C\. HonNumerical justification of fundamental solutions and the quasi\-Monte Carlo method for Poisson\-type equations\.Engineering Analysis with Boundary Elements22\(1\),pp\. 61–69\.External Links:[Document](https://dx.doi.org/10.1016/s0955-7997%2898%2900036-8)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1),[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px2.p1.1)\.
- Chen and Chen \(1995\)T\. Chen and H\. ChenUniversal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems\.IEEE Transactions on Neural Networks6\(4\),pp\. 911–917\.External Links:[Document](https://dx.doi.org/10.1109/72.392253)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1),[§2\.1](https://arxiv.org/html/2609.35938#S2.SS1.p1.3),[§2\.2](https://arxiv.org/html/2609.35938#S2.SS2.p1.1),[§2\.5](https://arxiv.org/html/2609.35938#S2.SS5.SSS0.Px3.p5.1)\.
- Cortes and Vapnik \(1995\)C\. Cortes and V\. VapnikSupport\-vector networks\.Machine Learning20\(3\),pp\. 273–297\.External Links:[Document](https://dx.doi.org/10.1007/bf00994018)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1)\.
- Cuomoet al\.\(2022\)S\. Cuomo, V\. S\. Di Cola, F\. Giampaolo, G\. Rozza, M\. Raissi, and F\. PiccialliScientific machine learning through physics\-informed neural networks: where we are and what’s next\.Journal of Scientific Computing92\(3\),pp\. 88\.External Links:[Document](https://dx.doi.org/10.1007/s10915-022-01939-z)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- Cybenko \(1989\)G\. CybenkoApproximation by superpositions of a sigmoidal function\.Mathematics of Control, Signals, and Systems2\(4\),pp\. 303–314\.External Links:[Document](https://dx.doi.org/10.1007/bf02551274)Cited by:[§2\.2](https://arxiv.org/html/2609.35938#S2.SS2.p1.1)\.
- Eshaghiet al\.\(2025\)M\. S\. Eshaghi, C\. Anitescu, M\. Thombre, Y\. Wang, X\. Zhuang, and T\. RabczukVariational physics\-informed neural operator \(VINO\) for solving partial differential equations\.Computer Methods in Applied Mechanics and Engineering437,pp\. 117785\.External Links:[Document](https://dx.doi.org/10.1016/j.cma.2025.117785)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- Eshaghiet al\.\(2026\)M\. S\. Eshaghi, N\. Valizadeh, C\. Anitescu, Y\. Wang, X\. Zhuang, and T\. RabczukMulti\-head neural operator for modeling interfacial dynamics\.International Journal of Mechanical Sciences317,pp\. 111363\.External Links:[Document](https://dx.doi.org/10.1016/j.ijmecsci.2026.111363)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- Fairweatheret al\.\(2003\)G\. Fairweather, A\. Karageorghis, and P\. A\. MartinThe method of fundamental solutions for scattering and radiation problems\.Engineering Analysis with Boundary Elements27\(7\),pp\. 759–769\.External Links:[Document](https://dx.doi.org/10.1016/s0955-7997%2803%2900017-1)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1),[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px2.p1.1)\.
- Fairweather and Karageorghis \(1998\)G\. Fairweather and A\. KarageorghisThe method of fundamental solutions for elliptic boundary value problems\.Advances in Computational Mathematics9\(1\-2\),pp\. 69–95\.External Links:[Document](https://dx.doi.org/10.1023/a%3A1018981221740)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1),[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px2.p1.1),[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px2.p1.2),[§2\.5](https://arxiv.org/html/2609.35938#S2.SS5.SSS0.Px2.p1.1)\.
- Fasshauer \(2007\)G\. E\. FasshauerMeshfree approximation methods with MATLAB\.Interdisciplinary Mathematical Sciences, Vol\.6,World Scientific\.External Links:[Document](https://dx.doi.org/10.1142/6437)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1),[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px1.p1.2)\.
- Fuet al\.\(2020\)Z\. Fu, Q\. Xi, Y\. Li, H\. Huang, and T\. RabczukHybrid fem\-sbm solver for structural vibration induced underwater acoustic radiation in shallow marine environment\.Computer Methods in Applied Mechanics and Engineering369,pp\. 113236\.External Links:[Document](https://dx.doi.org/10.1016/j.cma.2020.113236)Cited by:[§C\.2](https://arxiv.org/html/2609.35938#A3.SS2.p1.1),[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px2.p2.1),[§3\.4](https://arxiv.org/html/2609.35938#S3.SS4.p1.1),[§3\.4](https://arxiv.org/html/2609.35938#S3.SS4.p2.1)\.
- Fuet al\.\(2024\)Z\. Fu, W\. Xu, and S\. LiuPhysics\-informed kernel function neural networks for solving partial differential equations\.Neural Networks172,pp\. 106098\.External Links:[Document](https://dx.doi.org/10.1016/j.neunet.2024.106098)Cited by:[Appendix A](https://arxiv.org/html/2609.35938#A1.p1.1),[§2\.3\.1](https://arxiv.org/html/2609.35938#S2.SS3.SSS1.p2.1),[§2\.4\.2](https://arxiv.org/html/2609.35938#S2.SS4.SSS2.p1.1)\.
- Ginet al\.\(2021\)C\. R\. Gin, D\. E\. Shea, S\. L\. Brunton, and J\. N\. KutzDeepGreen: deep learning of Green’s functions for nonlinear boundary value problems\.Scientific Reports11\(1\),pp\. 21610\.External Links:[Document](https://dx.doi.org/10.1038/s41598-021-00773-x)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- Girosiet al\.\(1995\)F\. Girosi, M\. Jones, and T\. PoggioRegularization theory and neural networks architectures\.Neural Computation7\(2\),pp\. 219–269\.External Links:[Document](https://dx.doi.org/10.1162/neco.1995.7.2.219)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1),[§2\.3\.1](https://arxiv.org/html/2609.35938#S2.SS3.SSS1.p1.1)\.
- Golberg and Chen \(1999\)M\. A\. Golberg and C\. S\. ChenThe method of fundamental solutions for potential, Helmholtz and diffusion problems\.InBoundary Integral Methods: Numerical and Mathematical Aspects,M\. A\. Golberg \(Ed\.\),Computational Engineering, Vol\.1,pp\. 103–176\.External Links:ISBN 9781853125294Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1),[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px2.p1.1),[§2\.4\.2](https://arxiv.org/html/2609.35938#S2.SS4.SSS2.p4.1),[§2\.5](https://arxiv.org/html/2609.35938#S2.SS5.SSS0.Px2.p1.1)\.
- Golberg \(1995\)M\. A\. GolbergThe method of fundamental solutions for Poisson’s equation\.Engineering Analysis with Boundary Elements16\(3\),pp\. 205–213\.External Links:[Document](https://dx.doi.org/10.1016/0955-7997%2895%2900062-3)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1),[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px2.p1.1)\.
- Guet al\.\(2011\)Y\. Gu, W\. Chen, and C\. ZhangSingular boundary method for solving plane strain elastostatic problems\.International Journal of Solids and Structures48\(18\),pp\. 2549–2556\.External Links:[Document](https://dx.doi.org/10.1016/j.ijsolstr.2011.05.007)Cited by:[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px2.p2.1),[§2\.4\.2](https://arxiv.org/html/2609.35938#S2.SS4.SSS2.p4.1)\.
- Haoet al\.\(2023\)Z\. Hao, Z\. Wang, H\. Su, C\. Ying, Y\. Dong, S\. Liu, Z\. Cheng, J\. Song, and J\. ZhuGNOT: a general neural operator transformer for operator learning\.External Links:2302\.14376Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- Hofmannet al\.\(2008\)T\. Hofmann, B\. Schölkopf, and A\. J\. SmolaKernel methods in machine learning\.The Annals of Statistics36\(3\),pp\. 1171–1220\.External Links:[Document](https://dx.doi.org/10.1214/009053607000000677)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1)\.
- Horniket al\.\(1989\)K\. Hornik, M\. Stinchcombe, and H\. WhiteMultilayer feedforward networks are universal approximators\.Neural Networks2\(5\),pp\. 359–366\.External Links:[Document](https://dx.doi.org/10.1016/0893-6080%2889%2990020-8)Cited by:[§2\.2](https://arxiv.org/html/2609.35938#S2.SS2.p1.1)\.
- Jinet al\.\(2022\)P\. Jin, S\. Meng, and L\. LuMIONet: learning multiple\-input operators via tensor product\.SIAM Journal on Scientific Computing44\(6\),pp\. A3490–A3514\.External Links:[Document](https://dx.doi.org/10.1137/22m1477751)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- Kansa \(1990\)E\. J\. KansaMultiquadrics—a scattered data approximation scheme with applications to computational fluid\-dynamics—I surface approximations and partial derivative estimates\.Computers & Mathematics with Applications19\(8\-9\),pp\. 127–145\.External Links:[Document](https://dx.doi.org/10.1016/0898-1221%2890%2990270-t)Cited by:[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px1.p1.1)\.
- Karniadakiset al\.\(2021\)G\. E\. Karniadakis, I\. G\. Kevrekidis, L\. Lu, P\. Perdikaris, S\. Wang, and L\. YangPhysics\-informed machine learning\.Nature Reviews Physics3\(6\),pp\. 422–440\.External Links:[Document](https://dx.doi.org/10.1038/s42254-021-00314-5)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- Kingma and Ba \(2014\)D\. P\. Kingma and J\. BaAdam: a method for stochastic optimization\.External Links:1412\.6980Cited by:[§3](https://arxiv.org/html/2609.35938#S3.p2.1)\.
- Kita and Kamiya \(1995\)E\. Kita and N\. KamiyaTrefftz method: an overview\.Advances in Engineering Software24\(1\-3\),pp\. 3–12\.External Links:[Document](https://dx.doi.org/10.1016/0965-9978%2895%2900067-4)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1),[§2\.3\.1](https://arxiv.org/html/2609.35938#S2.SS3.SSS1.p2.1),[§2\.4\.2](https://arxiv.org/html/2609.35938#S2.SS4.SSS2.p4.1),[§2\.5](https://arxiv.org/html/2609.35938#S2.SS5.SSS0.Px2.p1.1)\.
- Kovachkiet al\.\(2021\)N\. Kovachki, Z\. Li, B\. Liu, K\. Azizzadenesheli, K\. Bhattacharya, A\. Stuart, and A\. AnandkumarNeural operator: learning maps between function spaces\.External Links:2108\.08481Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1),[§2\.1](https://arxiv.org/html/2609.35938#S2.SS1.p1.3)\.
- Kurthet al\.\(2023\)T\. Kurth, S\. Subramanian, P\. Harrington, J\. Pathak, M\. Mardani, D\. Hall, A\. Miele, K\. Kashinath, and A\. AnandkumarFourCastNet: accelerating global high\-resolution weather forecasting using adaptive Fourier neural operators\.InProceedings of the Platform for Advanced Scientific Computing Conference,pp\. 1–11\.External Links:[Document](https://dx.doi.org/10.1145/3592979.3593412)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- LeCunet al\.\(2015\)Y\. LeCun, Y\. Bengio, and G\. HintonDeep learning\.Nature521\(7553\),pp\. 436–444\.External Links:[Document](https://dx.doi.org/10.1038/nature14539)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- Liet al\.\(2022\)Z\. Li, D\. Z\. Huang, B\. Liu, and A\. AnandkumarFourier neural operator with learned deformations for PDEs on general geometries\.External Links:2207\.05209Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- Liet al\.\(2020a\)Z\. Li, N\. Kovachki, K\. Azizzadenesheli, B\. Liu, K\. Bhattacharya, A\. Stuart, and A\. AnandkumarFourier neural operator for parametric partial differential equations\.External Links:2010\.08895Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1),[§2\.1](https://arxiv.org/html/2609.35938#S2.SS1.p1.3)\.
- Liet al\.\(2020b\)Z\. Li, N\. Kovachki, K\. Azizzadenesheli, B\. Liu, K\. Bhattacharya, A\. Stuart, and A\. AnandkumarMultipole graph neural operator for parametric partial differential equations\.External Links:2006\.09535Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- Liu \(2009\)G\. LiuMeshfree methods: moving beyond the finite element method\.CRC Press\.External Links:[Document](https://dx.doi.org/10.1201/9781420082104)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1)\.
- Luet al\.\(2021a\)L\. Lu, P\. Jin, G\. Pang, Z\. Zhang, and G\. E\. KarniadakisLearning nonlinear operators via DeepONet based on the universal approximation theorem of operators\.Nature Machine Intelligence3\(3\),pp\. 218–229\.External Links:[Document](https://dx.doi.org/10.1038/s42256-021-00302-5)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1),[§2\.1](https://arxiv.org/html/2609.35938#S2.SS1.p1.3),[§2\.4](https://arxiv.org/html/2609.35938#S2.SS4.p1.1)\.
- Luet al\.\(2021b\)L\. Lu, X\. Meng, Z\. Mao, and G\. E\. KarniadakisDeepXDE: a deep learning library for solving differential equations\.SIAM Review63\(1\),pp\. 208–228\.External Links:[Document](https://dx.doi.org/10.1137/19m1274067)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- Micchelli \(1986\)C\. A\. MicchelliInterpolation of scattered data: distance matrices and conditionally positive definite functions\.Constructive Approximation2\(1\),pp\. 11–22\.External Links:[Document](https://dx.doi.org/10.1007/bf01893414)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1),[§2\.3\.1](https://arxiv.org/html/2609.35938#S2.SS3.SSS1.p1.1)\.
- Monaghan \(2005\)J\. J\. MonaghanSmoothed particle hydrodynamics\.Reports on Progress in Physics68\(8\),pp\. 1703–1759\.External Links:[Document](https://dx.doi.org/10.1088/0034-4885/68/8/r01)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1)\.
- Nardini and Brebbia \(1983\)D\. Nardini and C\. A\. BrebbiaA new approach to free vibration analysis using boundary elements\.Applied Mathematical Modelling7\(3\),pp\. 157–162\.External Links:[Document](https://dx.doi.org/10.1016/0307-904X%2883%2990003-3)Cited by:[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px4.p2.2),[§2\.4\.4](https://arxiv.org/html/2609.35938#S2.SS4.SSS4.p3.1)\.
- Park and Sandberg \(1991\)J\. Park and I\. W\. SandbergUniversal approximation using radial\-basis\-function networks\.Neural Computation3\(2\),pp\. 246–257\.External Links:[Document](https://dx.doi.org/10.1162/neco.1991.3.2.246)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1),[§2\.3\.1](https://arxiv.org/html/2609.35938#S2.SS3.SSS1.p1.1),[§2\.5](https://arxiv.org/html/2609.35938#S2.SS5.SSS0.Px3.p5.1)\.
- Partridgeet al\.\(1991\)P\. W\. Partridge, C\. A\. Brebbia, and L\. C\. WrobelThe dual reciprocity boundary element method\.Springer Netherlands,Dordrecht\.External Links:ISBN 9789401136907,[Document](https://dx.doi.org/10.1007/978-94-011-3690-7)Cited by:[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px4.p2.2),[§2\.4\.4](https://arxiv.org/html/2609.35938#S2.SS4.SSS4.p3.1)\.
- Raissiet al\.\(2019\)M\. Raissi, P\. Perdikaris, and G\. E\. KarniadakisPhysics\-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations\.Journal of Computational Physics378,pp\. 686–707\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2018.10.045)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1),[§2\.4\.2](https://arxiv.org/html/2609.35938#S2.SS4.SSS2.p5.2)\.
- Rasmussen and Williams \(2006\)C\. E\. Rasmussen and C\. K\. I\. WilliamsGaussian processes for machine learning\.MIT Press\.External Links:[Document](https://dx.doi.org/10.7551/mitpress/3206.001.0001)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1)\.
- Reddy \(1984\)J\. N\. ReddyAn introduction to the finite element method\.McGraw\-Hill Book Company,London\.Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p1.1)\.
- Schaback and Wendland \(2006\)R\. Schaback and H\. WendlandKernel techniques: from machine learning to meshless methods\.Acta Numerica15,pp\. 543–639\.External Links:[Document](https://dx.doi.org/10.1017/s0962492906270016)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1),[§2\.5](https://arxiv.org/html/2609.35938#S2.SS5.SSS0.Px1.p1.1),[§2\.5](https://arxiv.org/html/2609.35938#S2.SS5.SSS0.Px1.p1.2)\.
- Schölkopf and Smola \(2002\)B\. Schölkopf and A\. J\. SmolaLearning with kernels: support vector machines, regularization, optimization, and beyond\.MIT Press\.External Links:[Document](https://dx.doi.org/10.7551/mitpress/4175.001.0001)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1)\.
- Smith \(1985\)G\. D\. SmithNumerical solution of partial differential equations: finite difference methods\.3rd edition,Oxford University Press,Oxford\.Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p1.1)\.
- Sun and Yao \(1997\)H\. Sun and W\. YaoVirtual boundary element\-linear complementary equations for solving the elastic obstacle problems of thin plate\.Finite Elements in Analysis and Design27\(2\),pp\. 153–161\.External Links:[Document](https://dx.doi.org/10.1016/S0168-874X%2896%2900087-X)Cited by:[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px3.p4.1)\.
- Wanget al\.\(2021\)S\. Wang, H\. Wang, and P\. PerdikarisLearning the solution operator of parametric partial differential equations with physics\-informed DeepONets\.Science Advances7\(40\),pp\. eabi8605\.External Links:[Document](https://dx.doi.org/10.1126/sciadv.abi8605)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1),[§2\.2](https://arxiv.org/html/2609.35938#S2.SS2.p4.1),[§2\.4\.2](https://arxiv.org/html/2609.35938#S2.SS4.SSS2.p5.2)\.
- Wanget al\.\(2026\)Y\. Wang, Z\. Hao, M\. S\. Eshaghi, C\. Anitescu, X\. Zhuang, T\. Rabczuk, and Y\. LiuPretrain finite element method: a pretraining and warm\-start framework for PDEs via physics\-informed neural operators\.Journal of the Mechanics and Physics of Solids214,pp\. 106682\.External Links:[Document](https://dx.doi.org/10.1016/j.jmps.2026.106682)Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p2.1)\.
- Wenet al\.\(2026\)S\. Wen, M\. Chi, T\. Yu, B\. Moseley, M\. Y\. Michelis, P\. Ren, H\. Sun, and S\. MishraLearning, solving and optimizing PDEs with TensorGalerkin: an efficient high\-performance Galerkin assembly algorithm\.External Links:2602\.05052Cited by:[§D\.3](https://arxiv.org/html/2609.35938#A4.SS3.p2.1)\.
- Wendland \(2005\)H\. WendlandScattered data approximation\.Cambridge University Press\.External Links:[Document](https://dx.doi.org/10.1017/cbo9780511617539)Cited by:[§2\.3\.1](https://arxiv.org/html/2609.35938#S2.SS3.SSS1.p1.1),[§2\.5](https://arxiv.org/html/2609.35938#S2.SS5.SSS0.Px1.p1.2)\.
- Williams and Rasmussen \(1996\)C\. K\. I\. Williams and C\. E\. RasmussenGaussian processes for regression\.InAdvances in Neural Information Processing Systems 8 \(NIPS 1995\),pp\. 514–520\.Cited by:[§1](https://arxiv.org/html/2609.35938#S1.p3.1)\.
- Wrobel and Brebbia \(1987\)L\. C\. Wrobel and C\. A\. BrebbiaThe dual reciprocity boundary element formulation for nonlinear diffusion problems\.Computer Methods in Applied Mechanics and Engineering65\(2\),pp\. 147–164\.External Links:[Document](https://dx.doi.org/10.1016/0045-7825%2887%2990010-7)Cited by:[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px4.p1.1),[§2\.3\.2](https://arxiv.org/html/2609.35938#S2.SS3.SSS2.Px4.p2.2)\.