ParticleWNN: a Novel Neural Networks Framework for Solving Partial
  Differential Equations
Deep neural networks (DNNs) have been widely used to solve partial
differential equations (PDEs) in recent years. In this work, a novel deep
learning-based framework named Particle Weak-form based Neural Networks
(ParticleWNN) is developed for solving PDEs in the weak form. In this
framework, the trial space is defined as the space of DNNs, while the test
space consists of functions compactly supported in extremely small regions,
centered around particles. To facilitate the training of neural networks, an
R-adaptive strategy is designed to adaptively modify the radius of regions
during training. The ParticleWNN inherits the benefits of weak/variational
formulation, requiring less regularity of the solution and a small number of
quadrature points for computing integrals. Additionally, due to the special
construction of the test functions, ParticleWNN enables parallel implementation
and integral calculations only in extremely small regions. This framework is
particularly desirable for solving problems with high-dimensional and complex
domains. The efficiency and accuracy of ParticleWNN are demonstrated through
several numerical examples, showcasing its superiority over state-of-the-art
methods. The source code for the numerical examples presented in this paper is
available at https://github.com/yaohua32/ParticleWNN.
Yaohua Zang, Gang Bao
# ParticleWNN: a Novel Neural Networks Framework for Solving Partial Differential Equations

###### Abstract

Deep neural networks (DNNs) have been widely used to solve partial differential equations (PDEs) in recent years. In this work, a novel deep learning-based framework named Particle Weak-form based Neural Networks (ParticleWNN) is developed for solving PDEs in the weak form. In this framework, the trial space is chosen as the space of DNNs, and the test space is constructed by functions compactly supported in extremely small regions whose centers are particles. To train the neural networks, an R-adaptive strategy is designed to adaptively modify the radius of regions during training. The ParticleWNN inherits the advantages of weak/variational formulation, such as requiring less regularity of the solution and a small number of quadrature points for computing the integrals. Moreover, due to the special construction of the test functions, the ParticleWNN allows local training of networks, parallel implementation, and integral calculations only in extremely small regions. The framework is particularly desirable for solving problems with high-dimensional and complex domains. The efficiency and accuracy of the ParticleWNN are demonstrated with several numerical examples. The numerical results show clear advantages of the ParticleWNN over the state-of-the-art methods.

## 1 Introduction

Due to the powerful approximation ability of DNNs, the DNN-based methods for solving PDEs and related inverse problems have sprung up in recent years. Generally, the PDE-based forward and inverse problems are governed by the following PDE

\[\mathcal{A}[u(t,\mathbf{x});\mathbf{\gamma}] =0,\quad\text{in }(0,T]\times\Omega, \tag{1a}\] \[\mathcal{B}[u(t,\mathbf{x});\mathbf{\gamma}] =g(t,\mathbf{x}),\quad\text{on }[0,T]\times\partial\Omega,\] (1b) \[\mathcal{I}[u(0,\mathbf{x});\mathbf{\gamma}] =h(\mathbf{x}),\quad\text{on }\Omega, \tag{1c}\]

where \(\Omega\) is a bounded domain in \(\mathbb{R}^{d}\) with boundary \(\partial\Omega\), \(T>0\), and \(u(t,\mathbf{x})\) is the solution of the PDE. Here, \(\mathcal{A}\) is the forward operator, which can be of parabolic, hyperbolic, or elliptic type. \(\mathcal{B}\) and \(\mathcal{I}\) represent the boundary condition operator and the initial condition operator, respectively. The function \(\mathbf{\gamma}\) denotes parameters in the PDE. The forward problem is to solve the PDE problem given the function \(\gamma\) along with appropriate boundary conditions (1b) and/or initial conditions (1c). The inverse problem, on the other hand, is to determine the function \(\gamma\) from additional boundary or interior measurements of the solution.

Among the existing DNN-based PDE solvers, one popular framework is based on the strong form of PDEs. The most notable one is the physics-informed neural networks (PINN) [30], which approximates the PDE solutions with DNNs and formulates the loss function as a combination of the strong-form PDE residuals and data mismatch. However, the strong-form methods usually require a massive amount of collocation points, leading to a high training cost [11]. Moreover, the
non-convex optimization problem obtained by strong-form methods may result in many local minima, which makes learning complex physical phenomena more challenging[12; 18]. Different from the strong-form methods, the weak-form methods formulate the loss based on the weak/variational form of the PDEs, e.g., [41; 23; 19; 17; 16; 43; 2; 15]. Therefore, these methods have the advantage of requiring less smoothness of the solution, a smaller number of quadrature points for computing the integrals, and allowing local learning through domain decomposition [15]. Usually, these methods use DNNs to parameterize the trial functions. The difference lies in choosing different types of test functions, including global or piecewise polynomials [15; 17], localized non-overlapping high order polynomials [16], and deep neural networks [43; 2]. However, existing weak-form methods define the test functions in the whole domain or sub-domain, which requires accurate integration methods and a large number of integration points to reduce integral approximation errors. Although domain decomposition [19; 16] can alleviate the problem in low-dimensional cases, high-dimensionality makes domain decomposition much more difficult.

Main contributions.In this paper, a novel deep learning-based framework, ParticleWNN is developed for solving PDEs based on the weak form. The trial space is chosen as the space of DNNs, and the test space is constructed by functions compactly supported in extremely small regions whose centers are particles. There are several advantages of constructing test functions this way. First, the locally defined test functions allow the local training of network parameters [16] and enable parallel implementation. Second, since the particles are randomly sampled in the domain and the radius can be extremely small, the proposed method is desirable for high-dimensional and complex domain problems, and avoids calculating integrals over large domains. An R-adaptive strategy is then introduced to train the neural network, i.e., the upper bound of the radius decreases adaptively with iterations. Finally, the effectiveness of the ParticleWNN is confirmed with several numerical examples. Our primary contributions are listed below:

* We propose a novel weak-form framework for solving PDEs by using DNNs, where the test space is constructed by functions compactly supported in extremely small regions whose centers are particles.
* We develop several training techniques for the proposed framework to improve its convergence and accuracy.
* We demonstrate the efficiency of the proposed framework with several numerical examples and show its superiority over the state-of-the-art methods in dealing with difficult problems, such as PDE problems with complicated solutions and inverse problems with noisy data.

## 2 Related Works

According to different frameworks, deep learning-based PDE methods can be roughly divided into four categories: the strong-form methods, the weak-form methods, the hybrid methods of traditional methods and deep learning, and methods based on neural operators. Typically, the strong-form methods utilize DNNs to approximate the solution of the equation. Then, collocation points are generated on the domain and the domain boundary to evaluate the strong-from residuals of the governing equations in the domain and the mismatches of solutions on the boundary, respectively. Finally, the network parameters are learned by minimizing the weighted combination of residuals and mismatches. The most representative methods of this type include the deep Galerkin method (DGM) [33] and the physical-informed neural networks (PINN) method [30]. Based on PINN, more methods are proposed to solve fractional advection-diffusion equations (fPINN) [28], Navier-Stokes equations (NSFnets) [13], PDE-based inverse problems (B-PINN) [39], and PDEs in irregular domain (PhyGeoNet) [10]; or to improve the PINN from different aspects, including adaptive collocation points strategy [1; 25], domain decomposition [12; 11], loss modification [29; 42], and the new training strategy [18]. For additional information on the strong-from methods, refer to [14].

The weak-form methods, on the other hand, formulate the loss based on the weak or variational form of the PDEs. Based on the variational principle, the Deep Ritz method (DeepRitz) [41] solves PDEs by minimizing an energy function, which was further extended to PDEs with essential boundary conditions in the deep Nitsche method [23]. A deep domain decomposition method (D3M) was proposed based on the variational principle in [19]. In [43; 2], the weak adversarial network (WAN) was proposed based on the weak form to convert the problem into an operator norm minimization
problem. The variational neural networks (VarNet) method [17] is based on the Petrov-Galerkin method, in which the trial function was approximated by DNN and the test functions were a set of piecewise polynomials that have compact support in the domain. The variational formulation of PINN (VPINN) [15] has a similar formulation, except that the test functions were chosen to be polynomials with globally compact support over the domain. Combined with domain decomposition, the hp-VPINN [16] was developed based on the sub-domain Petrov-Galerkin method. It divides the domain into multiple non-overlapping subdomains (h-refinement) and sets test functions as localized non-overlapping high-order polynomials (p-refinement).

The hybrid methods combine traditional methods with DNN-based methods to overcome the shortcomings of pure DNN-based methods. The classical Galerkin method was combined with neural networks in [34] to solve initial boundary value problems. In [31], a DiscretizationNet method was developed to combine the finite volume discretization with a generative CNN-based encoder-decoder for solving the Navier-Stokes equations. A coupled automatic-numerical differentiation PINN (CANN PINN) was developed in [7], where the numerical differentiation-inspired method was coupled with the automatic differentiation to define the loss function. DNN-based methods were also combined with the finite difference method in [32; 35] and the finite element method in [40]. There are also recent methods based on neural operators [21; 24; 20; 22; 6], which focus on learning mappings between function spaces.

## 3 Methodology

### The ParticleWNN framework

To illustrate the proposed method, we consider the classical Poisson equation with the Dirichlet boundary condition

\[\begin{cases}-\Delta u(\mathbf{x})=f(\mathbf{x}),\quad\text{in}\;\Omega\subset\mathbb{ R}^{d},\\ u(\mathbf{x})=g(\mathbf{x}),\quad\text{on}\;\partial\Omega.\end{cases} \tag{2}\]

The weak formulation of Poisson's equation (2) involves finding a function in \(\{u\in H^{1}(\Omega)|u|_{\partial\Omega}=g\}\) such that, for all test functions \(v\in H^{1}_{0}(\Omega)\), the following equation holds

\[\int_{\Omega}\nabla u\cdot\nabla v\;d\mathbf{x}=\int_{\Omega}fv\;d\mathbf{x}, \tag{3}\]

where \(H^{1}(\Omega)\) denotes the Sobolev space of functions with square-integrable derivatives and \(H^{1}_{0}(\Omega)\) contains \(H^{1}(\Omega)\) functions with zero boundary conditions. Under appropriate conditions for \(f\) and \(g\), the weak form (3) allows a unique solution \(u\), which is called the _weak solution_[8]. Generally, the weak-form DNN-based methods approximate the function \(u\) with a neural network \(u_{NN}(\mathbf{x};\theta)\), which is usually comprised of \(l\) hidden layers with \(\mathcal{N}_{i}\) neurons in each layer and activation function \(\sigma(\cdot)\) that takes the following form

\[u_{NN}(\mathbf{x};\theta)=T^{(l+1)}\circ T^{l}\circ T^{(l-1)}\circ\cdots\circ T^{ (1)}(\mathbf{x}). \tag{4}\]

Here, the linear mapping \(T^{(l+1)}:\mathbb{R}^{\mathcal{N}_{l}}\rightarrow\mathbb{R}\) indicates the output layer, and \(T^{(i)}(\cdot)=\sigma(\mathbf{W}_{i}\cdot+\mathbf{b}_{i}),\;i=1,\cdots l\) are nonlinear mappings with weights \(\mathbf{W}_{i}\in\mathbb{R}^{\mathcal{N}_{i}\times\mathcal{N}_{i-1}}\) and biases \(\mathbf{b}_{i}\in\mathbb{R}^{\mathcal{N}_{i}}\), and \(\theta=\left\{\mathbf{W}_{i},\mathbf{b}_{i}\right\}_{i=1}^{l+1}\) collects the network parameters. Then, these methods train the network by minimizing a loss function that is usually defined as the root mean square (RMS) or mean squared error (MSE) of the weak-form residuals along with some penalty terms. Following [15; 16], we denote \(\mathcal{R}\) the weak-form residual in (3)

\[\mathcal{R}(u_{NN};v)=\int_{\Omega}\nabla u_{NN}\cdot\nabla v\;d\mathbf{x}-\int_ {\Omega}fv\;d\mathbf{x}. \tag{5}\]

Different choices of the test functions vary in different weak-form methods. For example, the VarNet [17] defines the test functions as a set of piecewise polynomials that have compact support over the entire domain, the VPINN [15] selects the test functions to be polynomials with globally compact support, the hp-VPINN [16] uses localized non-overlapping high-order polynomials as test functions, and the WAN [43] represents the test space with DNNs. The special feature of our work is to choose test functions to be compactly supported functions defined in small neighborhoods \(B(\mathbf{x}^{c},R)\subset\Omega\), where \(\mathbf{x}^{c}\) is a particle in \(\Omega\) and \(R\) is the radius of the neighborhood. Advantages of
defining the test function in this way include allowing the DNN to focus on extremely local regions, avoiding integrating over the entire domain, being able to parallelize, etc. Specifically, we choose the compactly supported radial basis functions (CSRBFs) as test functions in this work. Usually, the CSRBFs defined in \(B(\mathbf{x}^{c},R)\) have the following form

\[v(r)=\begin{cases}v_{+}(r),&r(\mathbf{x})\leq 1,\\ 0,&r(\mathbf{x})>1,\end{cases} \tag{6}\]

where \(r(\mathbf{x})=\frac{\|\mathbf{x}-\mathbf{x}^{c}\|}{R}\). In fact, any functions in \(H^{1}(\Omega)\) that is compactly supported in \(B(\mathbf{x}^{c},R)\) can be used as test functions. To improve the training efficiency, we use multiple test functions to formulate the loss function. That is, we generate \(N_{p}\) particles \(\left\{\mathbf{x}_{i}^{c}\right\}_{i}^{N_{p}}\) and the corresponding \(\left\{R_{i}\right\}_{i}^{N_{p}}\) randomly or with some rules in the domain 1, and then define \(N_{p}\) CSRBFs \(\left\{v_{i}\right\}_{i}^{N_{p}}\) in each small neighbourhood \(B(\mathbf{x}_{i}^{c},R_{i})\). Therefore, we obtain the MSE of the weak-form residuals

Footnote 1: To ensure that \(B(\mathbf{x}^{c},R)\subset\Omega\), we generate \(R\) first, and then sample \(\mathbf{x}^{c}\) in \(\tilde{\Omega}=\{\mathbf{x}\in\Omega|\text{dist}(\mathbf{x},\partial\Omega)\geq R\}\).

\[\mathcal{L}_{\mathcal{R}}=\frac{1}{N_{p}}\sum_{i=1}^{N_{p}}|\mathcal{R}(u_{NN} ;v_{i})|^{2}. \tag{7}\]

For the boundary condition (and/or initial condition), we can treat it as a penalty term

\[\mathcal{L}_{\mathcal{B}}=\frac{1}{N_{bd}}\sum_{j=1}^{N_{bd}}|\mathcal{B}[u_{ NN}(\mathbf{x}_{j})]-g(\mathbf{x}_{j})|^{2}, \tag{8}\]

where \(\{\mathbf{x}_{j}\}_{j=1}^{N_{bd}}\) are sampled points on the \(\partial\Omega\). Finally, we formulate our loss function as

\[\mathcal{L}(\theta)=\lambda_{\mathcal{R}}\mathcal{L}_{\mathcal{R}}+\lambda_{ \mathcal{B}}\mathcal{L}_{\mathcal{B}}, \tag{9}\]

where \(\lambda_{\mathcal{R}}\) and \(\lambda_{\mathcal{B}}\) are weight coefficients in the loss function.

### Calculation of the loss

To evaluate the loss in (9), we need to calculate \(N_{p}\) integrals that are defined in \(B(\mathbf{x}_{i}^{c},R_{i}),\;i=1,\cdots,N_{p}\). A straightforward way to evaluate integrals is to use the Monte Carlo integration. Unfortunately, it requires an immense sample size to ensure admissible integration errors. An alternative is the quadrature rule method. This method works efficiently in low-dimensional case or when the integrand is simple. However, in the case of high dimensional problems and complicated integrands, one needs to further increase the quadrature points, thus greatly increasing the computational cost. Other numerical techniques, such as sparse grids [27] and quasi-Monte Carlo integration [26], can also be employed. For our framework, thanks to the special construction of the test functions, we only need to evaluate integrals on the small region \(B(\mathbf{x}_{i}^{c},R_{i})\) rather than on the entire domain. Through a simple coordinate transformation, we can convert the calculation of integrals on \(N_{p}\) small regions into a standard region \(B(\mathbf{0},1)\). In fact, the coordinate transformation \(\mathbf{x}=\mathbf{s}R_{i}+\mathbf{x}_{i}^{c}\) can be applied in (5) to get

\[\mathcal{R}(u_{NN};v_{i})=R_{i}^{d}\bigg{(}\int_{B(\mathbf{0},1)}\nabla_{\mathbf{x}}u_ {NN}(\mathbf{x})\cdot\nabla_{r}v_{i}\cdot\nabla_{\mathbf{x}}r\;d\mathbf{s}-\int_{B(\mathbf{0},1)}f(\mathbf{x})v_{i}\;d\mathbf{s}\bigg{)}, \tag{10}\]

where we use \(\nabla_{\mathbf{x}}v_{i}=\nabla_{r}v_{i}\cdot\nabla_{\mathbf{x}}r\). Therefore, we only need to generate one set of integration points in \(B(\mathbf{0},1)\) to calculate \(N_{p}\) integrals. For example, assume that \(\left\{\mathbf{s}_{k}\right\}_{k=1}^{m}\) are \(K_{int}\) integration points uniformly sampled from \(B(\mathbf{0},1)\), then we denote \(\mathbf{x}_{k}^{(i)}:=\mathbf{s}_{k}*R_{i}+\mathbf{x}_{i}^{c}\) and approximate \(\mathcal{R}(u_{NN};v_{i})\) by

\[\mathcal{R}(u_{NN};v_{i})\approx\frac{R_{i}^{d}\mathcal{V}_{B(\mathbf{0},1)}}{K_{ int}}\sum_{k=1}^{K_{int}}\bigg{(}\nabla_{\mathbf{x}}u_{NN}(\mathbf{x}_{k}^{(i)})\cdot \nabla_{r}v_{i}(\mathbf{x}_{k}^{(i)})\cdot\nabla_{\mathbf{x}}r(\mathbf{x}_{k}^{(i)})-f( \mathbf{x}_{k}^{(i)})v_{i}(\mathbf{x}_{k}^{(i)})\bigg{)}, \tag{11}\]

where \(\mathcal{V}_{B(\mathbf{0},1)}\) indicates the volume of \(B(\mathbf{0},1)\). Consequently, the loss (7) can be approximated by

\[\mathcal{L}(\mathcal{R})\approx\frac{\mathcal{V}_{B(\mathbf{0},1)}^{2}}{N_{p}K_{ int}^{2}}\sum_{i=1}^{N_{p}}R_{i}^{2d}\bigg{(}\sum_{k=1}^{K_{int}}\bigg{(}\nabla_{ \mathbf{x}}u_{NN}(\mathbf{x}_{k}^{(i)})\cdot\nabla_{r}v_{i}(\mathbf{x}_{k}^{(i)})\cdot \nabla_{\mathbf{x}}r(\mathbf{x}_{k}^{(i)})-f(\mathbf{x}_{k}^{(i)})v_{i}(\mathbf{x}_{k}^{(i)}) \bigg{)}\bigg{)}^{2}. \tag{12}\]
From (12), we can see that the value of \(\mathcal{L}(\mathcal{R})\) depends on \(R_{i}^{2d}\). To avoid the training failure, we remove \(R_{i}^{2d}\) in (12). We also remove the fixed term \(\mathcal{V}_{B(\mathbf{0},1)}^{2}\), which does not affect training. We denote \(\tilde{\mathcal{L}}(\mathcal{R})\) as the approximation of \(\mathcal{L}(\mathcal{R})\) with \(R_{i}^{2d}\) and \(\mathcal{V}_{B(\mathbf{0},1)}^{2}\) removed. Then, the loss function (9) is modified and approximated as

\[\begin{split}\tilde{\mathcal{L}}(\theta)&=\frac{ \lambda_{\mathcal{R}}}{N_{p}K_{int}^{2}}\sum_{i=1}^{N_{p}}\bigg{(}\sum_{k=1}^{K _{int}}\bigg{(}\nabla_{\mathbf{x}}u_{NN}(\mathbf{x}_{k}^{(i)})\cdot\nabla_{r}v_{i}(\bm {x}_{k}^{(i)})\cdot\nabla_{\mathbf{x}}r(\mathbf{x}_{k}^{(i)})-f(\mathbf{x}_{k}^{(i)})v_{i}( \mathbf{x}_{k}^{(i)})\bigg{)}\bigg{)}^{2}\\ &\quad+\frac{\lambda_{\mathcal{B}}}{N_{bd}}\sum_{j=1}^{N_{bd}} \bigg{(}\mathcal{B}[u_{NN}(\mathbf{x}_{j})]-g(\mathbf{x}_{j})\bigg{)}^{2}.\end{split} \tag{13}\]

## 4 The Implementation Details

In this section, we discuss some training techniques that will improve the performance of the proposed ParticleWNN.

Selection of test functions.In our framework, the choice of test functions can be very diverse. For different integration methods, different test functions will have different effects on the integration error, and will eventually affect the accuracy of the method. For numerical examples in this paper, we choose meshgrid integration points and select CSRBFs as test functions. There are many types of CSRBFs to choose from, such as the Bump function [9], Wendland's function [37], Wu's function [38], Buhmann's function [5], and so on. One can also refer to [44] for how to construct it. Typically, we consider the following Wendland's type CSRBFs

\[\phi_{d,2}(r)=\begin{cases}\frac{(1-r)^{l+2}}{3}[(l^{2}+4l+3)r^{2}+(3l+6)r+3],\quad r<1,\\ 0,\quad r\geq 1,\end{cases} \tag{14}\]

where \(l=\lfloor d/2\rfloor+3\), \(d\) is the dimension of the domain, and \(\lfloor\cdot\rfloor\) indicates the flooring function.

R-adaptive strategy.In the ParticleWNN framework, \(R\) is a very important parameter, which determines the size of the interest area of the DNN. If \(R\) is too large, there will be more overlaps between compact supports, so the loss will become insensitive to the change of particles. Moreover, a large \(R\) usually means a complicated integrand, which results in a larger approximation error of the integral. However, \(R\) should not also be too small since it will cause the precision loss of floating-point numbers. In this work, we generate \(R_{i}\) for each particle \(\mathbf{x}_{i}^{c}\) from a range \([R_{min},R_{max}]\) where \(R_{min}\) is a small fixed number and \(R_{max}(\geq R_{min})\) is another small number that varies with iterations. Generally, there are three common ways for \(R_{max}\) to change with iterations, that is, **R-ascending**: \(R_{max}\) gradually increases with iterations until it reaches a given upper bound, **R-fixed**: \(R_{max}\) remains unchanged, and **R-descending**: \(R_{max}\) gradually decreases with iterations until a given lower bound is encountered. In the experiment, we found that the **R-descending** strategy has the best performance (see details in the Supplementary Material). Therefore, we adopt the **R-descending** strategy to train the ParticleWNN.

Adaptive selection of particles.A common way to improve the training of ParticleWNN is to select the particles with some smart rules. In this paper, we found that the \(topK\) technique can improve the performance of the proposed ParticleWNN (see the Supplementary Material for details). The implementation of the \(topK\) technique is as follows: first, we randomly sample \(N_{p}\) particles in the domain; then, we evaluate the square residuals for each particle with (11); finally, we select the largest \(topK(\leq N_{p})\) residuals to calculate the loss.

Finally, we summarize the algorithm for the ParticleWNN in 1.

## 5 Experiments

In this section, we show the efficiency of the ParticleWNN in both PDE problems and inverse problems. For comparison, we choose the DeepRitz method and the vanilla PINN as baselines.
Experimental setups.Without specific clarifications, we choose the trial function \(u_{NN}\) as a ResNet with \(4\) hidden layers and \(50\) nodes in each layer, and give the activations in specific problems. During the training process, we randomly sample \(N_{p}\) particles in the domain and \(N_{bd}=2d\tilde{N}_{bd}\) points (\(\tilde{N}_{bd}\) for each of \(2d\) sides) on the boundary, and sample radius \(R_{i}\) from \([R_{min},R_{max}]\) for each region. We set \(R_{min}=10^{-6}\), \(R_{max}=10^{-4}\), \(\lambda_{\mathcal{R}}=1\), \(\lambda_{\mathcal{B}}=5\), \(maxIter=20000\), and generate \(K_{int}\) integration points in \(B(0,1)\) to evaluate integrals. We adopt the Adam optimizer with \(\tau_{\theta}=0.001\) and apply the StepLR scheduler with a step size of 1 epoch and a \(\gamma=1-2./maxIter\). We use the Relative error \(\|u_{NN}-u\|_{2}/\|u\|_{2}\) and the maximum absolute error (MAE) \(\max|u_{NN}-u|\) as evaluation metrics and run each example with 5 random seeds to obtain the mean and the standard deviation. The experiments in Section 5.2 were executed using Kaggle Notebook with a Tesla P100 GPU, and all other experiments were executed using Google Colab [4] with a Tesla T4 GPU.

### One Dimensional Poisson's Equation

We first consider the following one dimensional Poisson's equation

\[\begin{cases}-\frac{\partial^{2}u}{\partial x^{2}}(x)=f(x),\quad x\in\Omega= (-1,1),\\ u(x)=g(x),\quad x\in\partial\Omega.\end{cases} \tag{15}\]

We construct a solution \(u(x)=x\cos(\omega x)\) for (15) and use it to calculate the forcing term \(f(x)\) and boundary condition \(g(x)\). It has been shown in [3] that, given a large frequency \(\omega\) (e.g., \(15\pi\)), the vanilla PINN failed to solve the problem with a Fully Connected Feedforward Neural Network with the activation function being the Tanh. Here, we use the Sine function as the activation in the ResNet. We set \(N_{p}=200\), \(topK=150\), \(K_{int}=49\), and \(\tilde{N}_{bd}=1\) for the proposed method. For comparison, we generate \(9800\) integration points in \(\Omega\) for the DeepRitz method and \(9800\) (\(topK=7350\)) collocation points for the vanilla PINN. For other parameters, we keep their settings consistent. We show the experiment results in Figure 1. It can be seen from Figure 1(c) that, although the DeepRitz method converges fastest, its accuracy is greatly limited by the number of integral sample points. The vanilla PINN does not have the problem of being limited by the integral points. However, it converges slowly and takes a long computation time (almost twice that of the other two methods). The ParticleWNN method converges quickly and obtains the smallest Relative error and MAE, see Figure 1(a) and 1(b). The point-wise error in Figure 1(d) indicates that the solution obtained by the ParticleWNN is the closest to the ground truth.

### Allen-Cahn equation

We then consider the following nonlinear Allen-Cahn equation

\[\begin{cases}u_{t}-\lambda u_{xx}+5u^{3}-5u=0,\quad t\in(0,1],\ x\in[-1,1],\\ u(x,0)=x^{2}\cos(\pi x),\\ u(t,-1)=u(t,1),\ u_{x}(t,-1)=u_{x}(t,1),\end{cases} \tag{16}\]

where \(\lambda=0.0001\). The equation is very difficult to solve as it has a sharp solution (see Figure 2(c)). To show the efficiency of the proposed method, we solve the problem with algorithm 1 by choosing
the test function \(v(\mathbf{x})\) only depends on \(\mathbf{x}\) and multiplying it with both sides of the equation (16). Then, we get the weak formulation through variation

\[\int_{\Omega}(u_{t}+5u^{3}-5u)v\;d\mathbf{x}+\lambda\int_{\Omega}\nabla_{\mathbf{x}}u \cdot\nabla_{\mathbf{x}}v\;d\mathbf{x}=0,\quad t\in(0,1]. \tag{17}\]

We approach the initial value conditions in the same way as we handled the boundary conditions and sample \(N_{init}=200\) points in \(\Omega\) to evaluate the mismatch of initial conditions. To evaluate the weak-form residuals, we sample \(N_{t}=100\) time points in \((0,1]\). At each time point, we randomly sample \(N_{p}=50\) particles in \(\Omega\) and set \(K_{int}=25\) to calculate the integrals in (17). For other parameters, we set \(N_{bd}=200\), \(topK=4000\), \(\lambda_{\mathcal{R}}=100\), \(maxIter=50000\), and increase the number of neurons per layer in the ResNet to \(100\) and use the composite function of the Tanh and the Sine as the activation. For comparison, we generate \(125000\) (\(topK=100000\)) collocation points for the vanilla PINN to evaluate the strong residuals and keep other settings consistent with the ParticleWNN. We show the experiment results in Figure 2, from which we can see that the ParticleWNN achieves higher accuracy and converges much faster than the vanilla PINN. The Figure 2(b) shows that the ParticleWNN takes almost half the computation time of the vanilla PINN.

### Inverse Problems

The particleWNN can also be applied in solving inverse problems. Typically, we consider the following coefficient identification problem which is often used to verify the effectiveness of DNN-based methods

\[-\nabla\cdot(a(x,y)\Delta u)=f(x,y),\quad(x,y)\in\Omega=[-1,1]\times[-1,1], \tag{18}\]

where \(u\) is the solution of the equation, \(f\) indicates the source term, and \(a\in\{a\in L^{\infty}(\Omega):0<\underline{a}\leq a(x,y)\leq\overline{a}\;a.e. \;in\;\Omega\}\) represents the coefficient. Given the source term \(f\), the inverse problem is to

Figure 1: Performance of the ParticleWNN, the vanilla PINN, and the DeepRitz in solving the problem (15). (a) Average Relative errors obtained by the ParticleWNN (\(0.0031\pm 0.0018\)), the vanilla PINN (\(0.0061\pm 0.0043\)), and the DeepRitz (\(0.1049\pm 0.0064\)); (b) Average MAEs obtained by the ParticleWNN (\(0.0032\pm 0.0019\)), the vanilla PINN (\(0.0045\pm 0.0025\)), and the DeepRitz (\(0.0893\pm 0.0035\)); (c) Average Relative errors vs. Average computation times. (d) Average pointwise errors.

Figure 2: Performance of the ParticleWNN and the vanilla PINN in solving the problem (16). (a) Average Relative errors obtained by the ParticleWNN (\(0.073\pm 0.015\)) and the vanilla PINN (\(0.323\pm 0.167\)); (b) Average Relative errors vs. average computation times. (c) The exact \(u\) and the average point-wise errors obtained by the ParticleWNN and the vanilla PINN.
identify the coefficient \(a\) with inexact measurements \(u^{\delta}\) where \(\delta\) indicates Gaussian noise. We set the exact coefficient as

\[a(x,y)=0.1+\exp^{-\frac{(x-x_{1})^{2}}{\sigma_{1}}+\frac{(y-y_{1})^{2}}{\sigma_{2 }}}, \tag{19}\]

where \(\sigma_{1},\sigma_{2}\) are sampled from \([0.01,0.5]\), and \(x_{1},y_{1}\) are sampled from \([-0.5,0.5]\), see Figure 3(a). For the sake of convenience, we manufacture a solution \(u(x,y)=\sin(\pi x)\sin(\pi y)\) and use the exact coefficient and solution to compute \(f\). We assume that the measurements \(u^{\delta}\) can be obtained from 100 sensors inside the domain and \(\tilde{N}_{bd}=25\) equally distributed sensors at each boundary for the Dirichlet boundary condition, see Figure 3(b). To solve the problem, we use two independent DNNs with the same structure to parameterize \(a\) and \(u\), respectively. We set \(N_{p}=200\), \(topK=150\), \(K_{int}=60\), and use the Tanh as the activation for the ParticleWNN method. For a fair comparison, we use \(12000\) (\(topK=9000\)) collocation points for the vanilla PINN and keep other settings consistent with the ParticleWNN. We consider two different noise levels \(\delta\sim\mathcal{N}(0,0.01^{2}),\;\delta\sim\mathcal{N}(0,0.1^{2})\), and show the experiment results in Table 1. From Table 1, we can see that the vanilla PINN is very sensitive to noise, while the ParticleWNN achieves good results at both noise levels. Moreover, the ParticleWNN converges much faster and takes less than half of the computation time of the vanilla PINN (See Figure 3(c)-3(d)).

### Additional ablations

In the following, we study the effect of \(N_{p}\), \(K_{int}\) on the ParticleWNN. For this purpose, we consider solving Poisson's equation (15) by the proposed ParticleWNN approach with different combinations of \(N_{p}\) and \(K_{int}\). We show experimental results in Figure 4 (Relative errors, MAEs, and computation times are recorded in the Supplementary Material). From Figures 4(a)-4(b), we can see that the algorithm performs worst when \(K_{int}=5\), and increasing \(N_{p}\) does not improve the algorithm much, indicating that the ParticleWNN still needs a sufficient number of integration points to obtain a considerable precision of integral approximation. When there are enough integration points, the number of particles will have a dominant influence on the algorithm. However, the consuming time increases linearly with \(N_{p}K_{int}\) (see Figure 4(c)), so it is necessary to select an appropriate combination to obtain a balance between accuracy and time consumption.

\begin{table}
\begin{tabular}{c|c c|c c} \hline \hline  & \multicolumn{2}{c|}{noise level 0.01} & \multicolumn{2}{c}{noise level 0.1} \\ \cline{2-5}  & ParticleWNN & vanilla PINN & ParticleWNN & vanilla PINN \\ \hline Relative error for \(a\) & \(\mathbf{0.036}\pm 0.005\) & \(0.755\pm 1.014\) & \(\mathbf{0.061}\pm 0.006\) & \(0.255\pm 0.393\) \\ MAE for \(a\) & \(\mathbf{0.071}\pm 0.009\) & \(0.613\pm 0.697\) & \(\mathbf{0.087}\pm 0.004\) & \(0.277\pm 0.408\) \\ Relative error for \(u\) & \(\mathbf{0.013}\pm 0.002\) & \(0.272\pm 0.318\) & \(\mathbf{0.022}\pm 0.004\) & \(0.139\pm 0.232\) \\ MAE for \(u\) & \(\mathbf{0.039}\pm 0.013\) & \(0.502\pm 0.577\) & \(\mathbf{0.058}\pm 0.009\) & \(0.277\pm 0.449\) \\ Time (s) & \(\mathbf{610.31}\pm 0.98\) & \(1441.65\pm 1.98\) & \(\mathbf{607.60}\pm 0.24\) & \(1436.66\pm 2.64\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Performance (avg\(\pm\)std) of the ParticleWNN and the vanilla PINN in solving the inverse problem (18) with different noise levels.

Figure 3: Performance of the ParticleWNN and the vanilla PINN in solving the inverse problem (18) with different noise levels. (a) The exact \(a\); (b) The exact \(u\) and sensors (black dots); (c) Average Relative error vs. average computation time at noise level 0.01; (d) Average Relative error vs. average computation time at noise level 0.1.
## 6 Conclusion

In this work, we propose a novel framework, ParticleWNN, for solving PDEs in weak-form with DNNs. The novelty of the framework is to employ test functions compactly supported in extremely small regions. The sizes and locations of these support regions can be selected arbitrarily. In this way, the method not only inherits the advantages of weak-form methods such as requiring less regularity of the solution and a small number of quadrature points for computing the integrals, but also outperforms them in solving complicated solution problems. In addition, the flexible definition of the test function makes ParticleWNN easily applicable for solving high-dimensional problems and complex domain problems. To improve the training of the ParticleWNN, several training strategies are developed, such as the R-descending strategy and adaptive selection of particles. Finally, we demonstrate the advantages of the ParticleWNN over other state-of-the-art methods in solving PDE problems and inverse problems.

Although the ParticleWNN avoids integrals over the entire domain or subdomains, it still suffers from the limitation imposed by integral calculations, such as errors caused by integral approximations. Therefore, an important future direction is to develop more efficient integral calculation techniques. Another future direction is to investigate new techniques to improve the ParticleWNN method, such as smarter particle selection rules and more efficient training strategies,including causal training technique [36] and seq2seq learning technique [3].

The potential risk is that the appropriate conditions for the algorithm to converge are still unclear due to the lack of theoretical analysis of DNNs, especially in the case of nonlinear PDEs. However, this is a problem common to all DNN-based methods, which in particular is known to researchers in this area. Overall, we believe our new framework is likely to have a positive impact on the society. It brings a novel framework for numerical solution of PDEs with diverse applications in science and technology.

## Acknowledgments and Disclosure of Funding

We are grateful to Qian Huang for providing valuable feedbacks on the draft. The work of G. Bao was supported in part by the National Natural Science Foundation of China (Grant No. 11621101, U21A20425) and by the Key Laboratory of Zhejiang Province.
