KP-PINNs: Kernel Packet Accelerated Physics Informed Neural Networks

1 Introduction

Differential equations, including ordinary differential equations (ODEs) and partial differential equations (PDEs), are widely applied in engineering for modelling, analyzing, and optimizing dynamic systems [Kreyszig, 2007]. They provide a mathematical framework for describing natural phenomena and solving complex problems. Differential equations are widely used in simulating fluid dynamics [Ragheb, 1976], electromagnetics [Deschamps, 1981], heat transfer [Isachenko et al., 1980] and acoustics [Kaltenbacher, 2018] Differential equations solutions are crucial for accurately predicting and optimizing the behaviour of dynamic systems, ensuring effective design and operation in engineering applications.

˚Corresponding author.

Because analytical solutions of differential equations are often difficult to obtain in practice, researchers seek numerical methods instead. For example, Euler's method and Runge-Kutta method [Atkinson et al., 2009] can be applied to solving ODEs, and finite element method (FEM), finite difference method (FDM), finite volume method (FVM), and spectral methods can be used to solve PDEs. These methods provide a foundation for solving various differential equations encountered in engineering and scientific applications. However, such methods require substantial computational resources and memory, particularly in high-dimensional or complex physical scenarios. Moreover, they lack flexibility, as adapting to different physical problems or parameter changes often needs significant model reconfiguration [Karniadakis et al., 2021].

With the development of deep learning, Physics informed neural networks (PINNs) incorporate physical laws into training neural networks to solve differential equations [Kianiharchegani, 2023], particularly PDEs, by effectively combining deep learning with domain-specific knowledge. PINNs have strong flexibility and can be effectively combined with existing traditional methods. However, PINNs can lead to incorrect solutions for some complex PDEs, affecting the accuracy of the solutions obtained [Zhang et al., 2024]. Several methods are proposed to address this issue. For example, [Haitsiukevich and Ilin, 2022] used ensemble agreement for domain expansion, and the improved algorithm can make the training of PINNs more stable. [Hu et al., 2021] theoretically analyzed the convergence and generalization of extended PINNs. [Yu et al., 2022] proposed gradient-enhanced PINNs, which utilize the gradient information of PDEs residuals and embed the gradient into the loss function to improve the effectiveness of solving PDEs. Other works include [Nguyen et al., 2023; Forootani et al., 2024; Aliakbari, 2023; Yuan et al., 2022; Lin et al., 2022].

However, it has been noticed that the original PINNs are not stable [Wang et al., 2022], due to the default L2 loss function. To address the instability of the original PINNs, [Wang et al., 2022] proposed the L8 loss function and proposed an adversarial training method. However, the computation of L8 norm is complicated, and the adversarial training is timeconsuming. Another class of loss functions is based on the Sobolev norm, which is proposed by [Son et al., 2021], and the corresponding method is named as Sobolev-PINNs. However, the computation of the Sobolev norm in Sobolev-PINNs is complex, and accurate derivative information is hard to get.

In this work, we propose a novel form of loss function using the Reproducing Kernel Hilbert Space (RKHS) norm. Under specific choice of the kernel function, the RKHS norm is equivalent to the (tensored) Sobolev norm, thus addressing the stability issue of the original PINNs. The proposed method can also avoid directly computing the derivative because the RKHS norm of a function can be approximated by the RKHS norm of its interpolation, while the later has a close form. To compute the RKHS norm, we employ Kernel Packet (KP) method [Chen et al., 2022], which can efficiently compute the inverse of the kernel matrix. Therefore, KP accelerated PINNs (KP-PINNs) can solve differential equations effectively and efficiently. We also provide theoretical guarantee for our proposed method.

The rest is organized as follows. Section 2 briefly introduces differential equations, PINNs and KP. Section 3 analyses the proposed KP-PINNs algorithm. Section 4 gives some theoretical results. Section 5 shows four numerical examples to demonstrate the performance of the KP-PINNs algorithm. Section 6 gives this research's discussion and future work.

2 Preliminaries

2.1 Differential Equations

Differential equations are mathematical equations that describe the relationship between a function and its derivatives, representing how a quantity changes over time or space. They are broadly classified into ordinary differential equations (ODEs), which involve functions of a single variable, and partial differential equations (PDEs), which involve functions of multiple variables. Differential equations are foundational in modeling dynamic systems and natural phenomena. By capturing the principles of change and interaction, they provide powerful tools for analyzing and solving complex problems. In general, a differential equation system can be expressed [Evans, 2022] by:

"

where Ω is an open domain with boundary BΩ, Lθ represents a partial differential operator, Bθ denotes the boundary condition, fθ and gθ are two known functions defined on Ω and BΩ, respectively, θ is the parameter, and u is an unknown function. The goal of solving a differential equation system is to find u such that (1) holds. Usually, (1) is called ODE if d " 1, and PDE if d ą 1. An example of PDEs is the second-order linear PDE, with

x " px1, ..., xdqJ P Rd , and the parameters are αi,j, βk, γ for i, j, k P t1, . . . , du.

Besides the linear PDEs, where Lθ is a linear differential operator, nonlinear PDEs, in which the solution or its derivatives appear nonlinearly, make them more complex to solve compared to linear PDEs. Because of the flexibility of nonlinear PDEs, it has been widely applied in fluid dynamics [Yan et al., 2025], optics [Evans and Souganidis, 1989], and biology [Ghergu and Radulescu, 2011]. One particular example is the Hamilton-Jacobi-Bellman (HJB) equations, which are widely applied in electro-hydraulic systems [Guo, 2022], energy systems [Sieniutycz, 2000] and finance [Witte and Reisinger, 2011]. The HJB equations characterise the optimal control strategy for a dynamic system under certain constraints and objectives. Specifically, the form of the HJB equations [Yong and Zhou, 2012] can be expressed by

where upx, tq is the value function, kpx, tq is the given control function, rpx, kq is the cost rate during the process and gpxq is the final cost at the terminal state [Wang et al., 2022]. If the cost rate function is rpx, kq " a1|k1|b1 `¨ ¨ ¨`ad|kd|bd ´ hpx, tq then the HJB equations can be expressed [Yong and Zhou, 2012] as

$'&'%

where Ai " paibiq´ 1 bi´1 ´ aipaibiq´ bi bi´1 ą 0 with ci " bi bi´1 ą 1, L0u :" B Bt u ´ ∆u is a linear operator, and ˜ LHJBu :" L0u ` řdi"1 Ai|Biu|ci is a nonlinear operator.

In the forward problem of differential equations, the goal is to obtain a function u such that (1) is satisfied, where all the operators and functions Lθ, Bθ, hθ and gθ are known. In practice, there are many scenarios where the parameter θ is unknown and needs to be estimated based on the equation form and observed data. This is called inverse problems, where the goal is to recover the parameters. In this work, we consider both forward problems and inverse problems.

2.2 Physics Informed Neural Networks

Physics informed neural networks (PINNs) provide a flexible way for numerically solving the differential equation systems (1).

For notational convenience, let

Clearly, u is the solution to (1) if and only if h1 " h2 " 0. Therefore, for the forward problem, it is natural to minimize the loss function

which was considered by [Raissi et al., 2019]. Choose xi P Ω, i " 1, ...,NL and si P BΩ, i " 1, ...,NB. Then (6) can be approximated by

with

PINNs use a deep neural network, denoted by unn, to minimize the loss function in (7).

For the inverse problems, suppose that we have observed function values at points xi P Ω, i " 1, ...,NL and si P BΩ, i " 1, ...,NB, denoted by udpxi q and udpsi q, respectively. We define the loss function in the inverse problem as

where Losspuq is as in (7), and

with h3 " u ´ ud. The goal of solving an inverse problem is to solve

where Θ is the parameter space. In practice, the minimizers of the loss functions in both forward and inverse problems can be found via optimization algorithms such as Adam [Kingma and Ba, 2014] and L-BFGS [Byrd et al., 1995].

2.3 Kernel Packet

Our proposed method relies on Kernel Packet (KP), which was proposed in [Chen et al., 2022]. In functional analysis, a Hilbert space of functions is termed a Reproducing Kernel Hilbert Space (RKHS) if the evaluation is a continuous linear operator at any point within the space. Assume that Ω Ă Rd is compact with Lipschtiz boundary, and Φ : ΩˆΩ Ñ R is a symmetric positive definite kernel function [Wang and Jing, 2022; Wang et al., 2020]. Define the linear space

and equip this space with the bilinear form

C

where the RKHS NΦpΩq generated by kernel Φ is defined as the closure of FΦpΩq under the inner product x¨, ¨yΦ. The norm of NΦpΩq is ||f||NΦpΩq " bxf, fyNΦpΩq, where x¨, ¨yNΦpΩq is induced by x¨, ¨yΦ.

Let x1, ..., xn P Ω Y BΩ be the points of interest. KP provides an efficient way to compute the inverse of the kernel matrix K :" pΦpxj, xkqqjk P Rnˆn , which is essential in the application of the RKHSs, for example, prediction in kernel ridge regression, and computing the posterior variance in Bayesian optimization. In KP, it is assumed that the kernel function Φ is a Mat´ern kernel function [Abramowitz and Stegun, 1968] as

where ℓ ą 0 is the scale parameter, and ν is the half integer smoothness parameter. For example, when ν is equal to 12 , the corresponding kernel function is Φ12 pxj, xkq " exp p´ℓ|xj ´ xk|q. The Mat´ern kernel function is widely applied in practice; see [Muyskens et al., 2024] for example.

For a moment, let us consider the case d " 1. Without loss of generality, assume x1 ă x2 ă . . . ă xn. A nonzero function ϕ is termed an s degree KP if it can admit the representation as

with the support of ϕ being the interval rx1, xns. In (11), Aj is the coefficients of Φpx, xjq, which can be obtained by solving the following linear systems:

where l " 0, ..., s´3 2 , δ " ˘1, c2 " 2ν ℓ2 (ν and ℓ are the smoothness parameter and scale parameter of Mat´ern kernel function respectively). By (11), we have AK " ϕpxq, and the inverse of the kernel matrix can be computed by

where x " px1, x2, . . . , xnqJ is the input vector, and the pi, jqth entry of ϕpXq and A are ϕjpxiq and the i element of Aj, respectively.

It has been shown in [Chen et al., 2022] that both A and ϕpxq are sparse banded matrices, hence the computation of K´1 can be fast. Specifically, [Chen et al., 2022] showed that computing A and ϕpxq needs only Opp2ν ` 2q3nq computation time and requires storage space of Opp2ν`2qnq, where ν is usually small in most practical applications. LU decomposition [Davis, 2006] can be applied for solving rϕpxqs´1 , and [Chen et al., 2022] showed that the total computation time is only Opn2 q.

If d ą 1, i.e., the input X is multi-dimensional, then special structure of the input points X " px1, x2, ..., xdq is needed. KP requires the input points to be tensor-like, that is, the input points can be described as the Cartesian combination of multiple one-dimensional point collections. Specifically, XFG " ŚNi"1 xpiq , where each xpiq represents a distinct one-dimensional point collection. With tensor-like input points, K can be decomposed into Kronecker products as

where Kxi " pΦpxij , xikqq, 1 ď i ď d, 1 ď j, k ď ni, ni is the number of points for dimension i and thus the dimension of K is n ˆ n, where n " śdi"1 ni. Then the inverse K can be directly computed by

where each K´1 xi can be computed using one-dimensional KP. Putting all things together, we obtain

whose computation time is still Opn2 q.

3 KP Accelerated PINNs

In this section, we introduce the proposed method, KP accelerated PINNs (KP-PINNs). Instead of using L2 loss in (6), which was shown to be unstable in some specific PDE systems, we consider tensor Sobolev norm. The tensor Sobolev space Hν T pΩq [Hochmuth et al., 2000] is defined as the set of functions f : ΩÑR satisfying

where ||f||HνT pΩq is the norm on Hν T pΩq, explicitly defined via the Fourier transform as

where ˜ fpωq represents the Fourier transform of f. We define the loss function with tensor Sobolev norm as

where λ ą 0 is a tuning parameter that balances the weights between the difference between the space and its boundary, and uNN is the predicted value obtained from the neural network. However, directly computing the tensor Sobolev norm may be time-consuming and inaccurate. By the equivalence between the RKHS norm and tensor Sobolev norm, we modify the loss function (19) to

where NΦνT is tensor RKHS, i.e. the RKHS generated by ΦνT pxq " śdi"1 Φν pxiq. Tensor RKHS and Sobolev space are widely considered for complexity reduction in highdimensional spaces [Ding et al., 2020; K¨uhn et al., 2015; D˜ung and Nguyen, 2021].

Computing the loss function in (20) can be efficient via KP. Specifically, let

and we consider ||h1NN||2NΦνT pΩq " ||LuNN ´ f||2NΦνT pΩq first. Based on the Theorem 11.23 in [Wendland, 2004], ||h1NN||2NΦνT pΩq can be approximated by

where yh1 " ph1NNpx1 q, ..., h1NNpxnqq and x1 , ..., xn P Rd means all the training data. K is as in (14). The other term ||h2NN||2NΦνT pBΩq can be approximated in a similar manner.

Algorithm 1 KP-PINNs

Input: PDE systems (1), iterations niter, known points Output: Approximated solution ˆupxq and equations' parameters (if inverse problem)

• 1: Initialize neural network parameters.
• 2: while i ă niter do
• 3: Forward pass and compute the predicted value.
• 4: Compute derivatives using automatic differentiation.
• 5: Compute A and ϕpxq by (12) and (11).
• 6: Compute the loss of KP-PINNs based on (23).
• 7: Update parameters.
• 8: end while

With (13), (16) and approximation (22), the loss function (20) can be approximated by:

By KP method introduced in Section 2.3, LossNΦνT can be easily computed.Some well-established optimization algorithms, such as Adam or L-BFGS, can be used to minimize the loss effectively. The KP-PINNs algorithm is summarized in Algorithm 1. For the inverse problem, comparing (7) and (8), it can be seen that the loss function of the inverse problem has more LossDpuq part. For this part, the computation of the RKHS norm is similar to the analysis above. In this way, the KP-PINNs algorithm can effectively solve both forward and inverse problems.

4 Theoretical Results

In this section, we analyze the stability of PDEs solved with the loss function (19), or equivalently (20). The stability of a PDE system is defined as follows [Wang et al., 2022].

Definition 1 (Stability of PDEs [Wang et al., 2022]). Suppose Z1,Z2, and Z3 are three Banach spaces. For the exact solution u˚ pxq and the predicted solution u, if for ||Lu´f ||Z1 , ||Bu´g||Z2 Ñ 0, it has ||u˚´u||Z3 " Op||Lu´ f||Z1 ` ||Bu ´ g||Z2 q, then the PDEs (1) is pZ1,Z2,Z3qstable.

Stability defined in Definition 1 provides a theoretical guarantee for solving a PDE. Specifically, if the loss function is ||Lu´f||Z1 `||Bu´g||Z2 , thenminimizing the loss function can ensure the convergence of the solution under Z3-norm. In this section, we consider two classes of PDEs: second-order linear elliptic equations and Hamilton-Jacobi-Bellman (HJB) equations. We start with the second-order linear elliptic equations satisfying the following assumption. While in the previous section, we considered the parameters of the PDE as constants, here we extend the analysis to a more general case, where the parameters are functions. This allows us to establish the theoretical results in a broader context.

Assumption 1. Suppose in (1), the operator L " řdi,j"1 ai,jpxqBxi Bxj ` řdi"1 bipxqBxi ` cpxq and B is the identity operator. The functions ai,j P C 2 satisfy the uniformly elliptic condition. In addition, the only solution with zero input data is the zero solution.

Lemma 1 (Theorem 2.1 in [Bramble and Schatz, 1970]). In addition to Assumption (1), assume that L is with C 8 coefficients, defined on C 8 bounded domain Ω on Rd . For any real number l, ||u||Hl pΩq ď C ´||Lu||Hl´2 pΩq ` ||u||Hl´ 1 2 pBΩq¯ for all u P C 8p¯Ωq, where ¯Ω is a closure of Ω and C is independent of u.

Theorem 1 (Stability of second-order linear elliptic equation). Suppose the RKHS NΦ1T coincides with the space H1 T . For any C 8 bounded domain Ω Ď Rd , if Assumption 1 is satisfied and all the coefficient functions are in C 8, then the equation

"

is pNΦ1T pΩq, NΦ1T pBΩq, H1 pΩqq-stable.

Proof. Sets l " 1, by Lemma 1, we have

Since the RKHS NΦ1T coincides with the space H1 T , we have that ||u||H1 pΩq ď C4p||Lu||NΦ1T pΩq ` ||Bu||NΦ1T pBΩqq.

The proof of Theorem 1 is a natural result derived from the classical theory of second-order linear elliptic equations [Shin et al., 2023]. It shows that the numerical solution for the second-order linear elliptic equation is stable. Next, we show the stability of the HJB equation (4).

Lemma 2 (Stability of HJB equation - L8 norm [Wang et al., 2022]). For p, q ě 1, let r0 " pd`2qq d`q . Suppose the subsequent inequalities are valid for p, q and r0:

where ˇ c " max1ďiďd ci in (4). For any r P r1, r0q, any bounded open set Ω Ă Rd ˆ r0, Ts, and for a large R such that Ω Ă BR ˆ r0, T s, where BR is a ball with radius R, (4) maintains pLppΩq, LqpBRq, W 1,r pΩqq-stability when ˇ c ď 2.

Theorem 2 (Stability of HJB equation). For p, q ě 1, let r0 " pd`2qq d`q . Suppose the subsequent inequalities are valid for p, q and r0:

where ˇ c " max1ďiďd ci in (4). Then, for any r P r1, r0q, ν1, ν2 ě 1 and any bounded open set Ω Ă Rd ˆ r0, Ts, (4) maintains pNΦν1 T pRd ˆ r0, Tsq,NΦν2 T pRd q,W 1,r pΩqq-stable when ˇ c ď 2.

Proof. First of all, according to the Sobolev embedding theorem [Adams and Fournier, 2003; Wendland, 2004; Ding et al., 2019], when bounded set Ω Ă Rd , p ă 8, @f P H1 , Dc1, c2, c3, c4, c5 ą 0, it has the norm inequality

If ||LHJBupxq ´ φpxq||Hν1 , ||BHJBupxq ´ gpxq||Hν2 Ñ 0, known by (24), ||LHJBupxq ´ φpxq||Lp , ||BHJBupxq ´ gpxq||Lq Ñ 0. By Lemma 2, (4) is stable. It can be represented as pLppΩq,LqpBRq,W 1,r pωqq-stable, and thus pHν1 pRd ˆ r0, Tsq,Hν2 pRd q,W 1,r pΩqqstable. Additionally, the equation is also pNΦ1T pRd ˆ r0, Tsq,NΦ1T pRd q,W 1,r pΩqq-stable, thus pNΦν1 T pRd ˆ r0, Tsq,NΦν2 T pRd q,W 1,r pΩqq-stable.

Theorem 2 is a modified version of Theorem 4.3 in [Wang et al., 2022]. It demonstrates that when the dimension of the state function d is large, selecting an appropriate RKHS ensures the stability of the equations.

5 Experiments

In this section, we apply KP-PINNs to four representative differential equations: an ODE - Stiff equation, a second-order linear PDE - Helmholtz equation, a HJB equation - LQG equation, and a nonlinear PDE - NS equation. For each differential equation, both the forward and inverse problems are addressed. We compare the KP-PINNs algorithm with several baseline PINNs approaches including L2-PINNs, RKHSPINNs, and Sobolev-PINNs, as summarized in Table 1. The implementation details and source code are available at: https://github.com/SiyuanYang-sy/KP-PINNs.

We consider the accuracy and computational time, where the accuracy is measured by the relative L2 error defined by

where upxiq represents the true solution at xi, ˆupxiq denotes the prediction obtained from the algorithm, and Ntest is the size of the test set. Also, the standard error (SE) is given to reflect the variability across multiple experiments.

5.1 Stiff Equation

The Stiff equation [Wanner and Hairer, 1996] exhibits steep gradients caused by fast dynamics, and is commonly treated with time-scale decomposition to reduce numerical and reaction-related errors. It is defined by

where λ and µ are two parameters. Here we set λ " ´2.0 and µ " 2.0. The analytic solution to (26) is

There is only one initial point so NB " 1, and NL " 50 in (7) in the forward problem. In the inverse problem, we assume that λ is unknown and set NB`NL " 50. In both cases Ntest is equal to 2000. Table 2 provides the average results of three independent runs. KP-PINNs and RKHS-PINNs yield similar results due to the same loss function form, but differ in the computation of K´1 , leading to the run-time gap observed in Table 6, where KP-PINNs demonstrate significantly faster computation across all four equations. The SobolevPINNs algorithm uses automatic differentiation to compute first- and higher-order derivatives, leading to error accumulation and difficulty in obtaining accurate results.

5.2 Helmholtz Equation

The Helmholtz equation is a second-order elliptic PDE commonly encountered in electromagnetism that describes the spatial behavior of electromagnetic wave propagation. It is defined by

where ∆u " B2 u Bx2 ` B2 u By2 , k is the wave number, p is the source term (excitation term) [Li et al., 2022] defined by ppx, yq " px ` yq sinpπxq sinpπyq ´ 2π2 px ` yq sinpπxq sinpπyq ` 2π cospπyq sinpπxq ` 2π cospπxq sinpπyq. Here, we set k " 1.0. The analytic solution for k " 1.0 [Jagtap et al., 2020] is

In this experiment, NB " 100 and NL " 100 ˆ 100 in (7). For the inverse problem, parameter k is unknown and NB`NL " 30ˆ30 in (7). In both casesNtest is equal to 500ˆ 600. Table 3 gives the average results of five times. It can be seen that the first three algorithms can get the solutions well. For inverse problems, only KP-PINNs (ν " 1 2 ) and RKHSPINNs (ν " 12 ) can predict the k relatively correctly, while the L2-PINNs fails to work. Figure 2 (b), Figure 2 (c) and Figure 3 (a) show the inverse problem results of L2-PINNs and KP-PINNs (ν " 12 ). The parameter ν determines the degree of smoothness for the kernel function, and a smaller value of ν is suitable for less smooth conditions. It can be seen that a smaller ν (such as ν " 12 ) can get better results.

5.3 Linear Quadratic Gaussian (LQG) Equation

One application of the HJB equation is LQG control problem. The form of LQG [Han et al., 2018] in d dimensions are

with the solution

where gpxq " ln ´ 1`||x||2 2 ¯ is the terminal cost function.

We set d " 2, µ " 1.0, T " 1.0,NB " 2000 and NL " 30 ˆ 20 ˆ 10 in (7) in the forward problem. For the inverse problem, µ is unknown and NB `NL " 10ˆ10ˆ10. In both cases, Ntest is equal to 100 ˆ 80 ˆ 50. Table 4 gives the average results of three times experiments. For the inverse problem, KP-PINNs with ν " 12 can obtain the closest estimation among all algorithms. Figure 3 (b) gives a one-case inverse problem results of KP-PINNs (ν " 12 ).

5.4 Navier-Stokes (NS) Equation

The widely-used NS equations are fundamental PDEs that describe the motion of viscous and incompressible fluids. In this experiment, the two-dimensional incompressible NS equations are considered to model the evolution of fluid flow. The governing equation for the vorticity ωpt, x, yq is

Table 6: Time of computing K´1 (ν " 12 as example). NB `NL " 2000 for Stiff equation. NB ` NL " 200 ˆ 300 for Helmholtz equation. NB`NL " 50ˆ50ˆ40 for LQG equation. NB`NL " 60 ˆ 40 ˆ 20 for NS equation.

where µ is the viscosity coefficient. In our experiment, we set µ " 0.01, T " 1.9, xstart " 1.0, xend " 8.0, ystart " ´2.0 and yend " 2.0. The initial condition is the situation of t " 0. The boundary conditions are the situation of x " xstart, x " xend, y " ystart and y " yend. The vorticity is defined in terms of the stream function ψpt, x, yq by the Poisson equation ω " ´pψxx ` ψyyq. The velocity components pu, vq can be recovered from the stream function as u " ψy, v " ´ψx. The two components of the velocity field data upt, x, yq and vpt, x, yq are obtained from [Raissi et al., 2019], where the numerical solution is performed. We approximate the stream function ψ using a neural network. The velocity field and vorticity can be derived by automatic differentiation.

We set NB " 5000 and NL " 30ˆ20ˆ10 in (7). For the inverse problem, µ is unknown and NB `NL " 10ˆ10ˆ10. In both cases, Ntest " 98 ˆ 48 ˆ 20. Table 5 shows the averaged results over three independent runs. It can be seen that KP-PINNs, L2-PINNs and RKHS-PINNs can solve the equation in general. Among them, KP-PINNs and RKHS-PINNs consistently demonstrate superior performance in terms of both accuracy and robustness. For the inverse problem, KPPINNs(ν " 1 2 ) can obtain the closest estimation.

6 Conclusion

This work proposes KP-PINNs, an efficient method for solving differential equations, with stability guarantees for second-order linear elliptic and HJB equations under certain conditions. Numerical experiments demonstrate that KPPINNs outperform L2-PINNs and Sobolev-PINNs for both forward and inverse problems, significantly accelerating the computation of the inverse of kernel matrix K through the KP technique. Future directions include extending KP-PINNs to handle sparse or irregular data, optimizing neural network architecture, and refining kernel selection in the loss function to improve accuracy and robustness. These advancements aim to enhance KP-PINNs' versatility in solving complex differential equations.

Acknowledgments

This research is supported by Guangdong Provincial Fund - Special Innovation Project (2024KTSCX038).