11footnotetext: Corresponding author. _Email address:_<EMAIL_ADDRESS>(Z. Shi)

# BINet: Learning to Solve Partial Differential Equations with Boundary
Integral Networks

Guochang Lin1, Pipi Hu1,3, Fukai Chen2, Xiang Chen4,
Junqing Chen2, Jun Wang5, Zuoqiang Shi2,3∗
1 Yau Mathematical Sciences Center, Tsinghua University, Beijing, China
2 Department of Mathematical Sciences, Tsinghua University, Beijing, China
3 Yanqi Lake Beijing Institute of Mathematical Sciences and Applications,
Beijing, China
4 Noah’s Ark Lab, Huawei, China
5 University College London, London, United Kingdom

Abstract: We propose a method combining boundary integral equations and neural
networks (BINet) to solve partial differential equations (PDEs) in both
bounded and unbounded domains. Unlike existing solutions that directly operate
over original PDEs, BINet learns to solve, as a proxy, associated boundary
integral equations using neural networks. The benefits are three-fold.
Firstly, only the boundary conditions need to be fitted since the PDE can be
automatically satisfied with single or double layer representations according
to the potential theory. Secondly, the dimension of the boundary integral
equations is typically smaller, and as such, the sample complexity can be
reduced significantly. Lastly, in the proposed method, all differential
operators of the original PDEs have been removed, hence the numerical
efficiency and stability are improved. Adopting neural tangent kernel (NTK)
techniques, we provide proof of the convergence of BINets in the limit that
the width of the neural network goes to infinity. Extensive numerical
experiments show that, without calculating high-order derivatives, BINet is
much easier to train and usually gives more accurate solutions, especially in
the cases that the boundary conditions are not smooth enough. Further, BINet
outperforms strong baselines for both one single PDE and parameterized PDEs in
the bounded and unbounded domains.

##  1 Introduction

Partial differential equations (PDEs) have been widely used in scientific
fields and engineering applications, such as Maxwell’s equations in optics and
electromagnetism [1], Navier–Stokes equations in fluid dynamics [2], the
Schrödinger equations in the quantum physics [3], and Black-Scholes equations
for call option pricing in finance [4]. Therefore, finding the solution to
PDEs has been a critical topic in research over the years. However, in most
cases, the analytical solution of PDEs is infeasible to obtain, such that
numerical methods become the major bridge between PDE models and practical
applications.

In the past decade, deep learning has achieved great success in computer
vision, natural language processing, and many other topics [5]. It is found
that deep neural networks (DNNs) have the attractive capability in
approximating functions, especially in high dimensional space. Therefore, DNNs
hold great potential in solving PDEs with the promise of providing a good
_ansatz_ to represent the solution, where the parameters can be obtained by
training DNNs with proper loss functions.

In the literature, many efforts have been devoted to developing DNN-based
methods for solving different kinds of PDEs, such as DGM [6], Deep-Ritz [7],
and PINN [8]. The main idea of these methods is to use a neural network to
approximate the solution of the PDE directly. The loss function is designed by
either incorporating the PDE residual and the boundary or initial conditions,
or the energy functional derived from the variational form of the PDE.

However, two important issues are not fully considered in most existing works.
First, PDEs are merely utilized to construct the loss function, and the
essence behind PDEs may be further explored to design a new network structure
to cater the need for solving differential equations. Second, when it comes to
complex problems, such as PDEs with oscillatory or even singular solutions,
failure of the aforementioned methods is frequently reported [9] due to high
order differentiation of the neural networks with respect to the inputs. The
appearance of high-order derivatives may lead to instability in training [10]
(for example, amplified oscillation or singularities) such that the network
can not find the exact solution.

To address the two issues above, in this paper, we propose a novel method,
named BINet, combining boundary integral equations and deep learning to solve
PDEs. Utilizing fundamental solutions of PDE and Green’s formula [11], the
solution to PDE can be expressed in the form of a boundary integral, where the
explicit fundamental solution of PDE serves as the integral kernel. A new
network structure is then designed based on this integral expression of the
solution such that the output of our network can satisfy the PDE
automatically. Since the PDE has been satisfied, we only need to take the
boundary condition as the supervisory signal for the loss. In BINet, the prior
information provided by the PDE is fully integrated into the network.
Moreover, the differential operator is substituted by an integral operator,
which avoids the extra differential operations of the neural networks. The
main advantages of BINet are summarized below:

First, BINet adopts an explicit integral representation of the solution such
that the output of BINet satisfies the original PDE automatically. This means
that the training of BINet is naturally confined in the solution space of PDE.
Since BINet is defined in a much smaller space, i.e., the solution function
space, the training of BINet is faster and more stable than the general neural
network. Another advantage of integral representation is that all differential
operators are removed in BINet. Then the regularity requirement of BINet is
relaxed significantly which enables BINet to approximate the solutions with
poor regularity. Moreover, BINet has good theoretical properties. Using neural
tangent kernel (NTK) techniques [12], BINet can be proved to converge as the
width of the neural network goes to infinity.

Second, since the PDE has been satisfied automatically with the integral
representation in BINet, the residual of the boundary condition is the only
component of the loss function. There is no need to balance the residual of
the PDE and the boundary condition, BINet thus fits the boundary condition
better with less parameter tuning.

Third, BINet can solve PDEs in the unbounded domain since the integral
representation holds for both bounded and unbounded domains. For some problems
such as electromagnetic wave propagation, solving PDEs in an unbounded domain
is critical and complicated using traditional methods. Moreover, existing
deep-learning-based models also suffer from the difficulty of sampling in
unbounded domains. Therefore, BINet provides a good choice to solve this kind
of problem.

Fourth, BINet is also capable to learn a solution operator mapping a
parameterized PDE to its solution by feeding the parameters to the network as
input. Note that in the integral representation of the solution, the integral
kernel, i.e., the fundamental solution to the original PDE, has an explicit
form dependent on the differential operator of the PDE. Moreover, the integral
is conducted exactly on the boundary of the domain on which the PDE is solved.
Therefore, BINet has great advantages in learning the solution operator
mapping differential operator or computational domain to the corresponding
solution.

At last, the boundary integral is defined on the boundary whose dimension is
less by 1 than the original computational domain. Lower dimension leads to
fewer sample points which will reduce the computational cost.

The rest of this paper is organized as follows. An overview of related work on
solving PDEs using deep learning approaches is given in Section 2. The
boundary integral method and BINet are introduced in Section 3. In Section 4,
we analyze the convergence of BINet using the NTK techniques. Extensive
numerical experiments are shown in Section 5. At last, conclusion remarks are
made in Section 6.

##  2 Related Work

Solving PDEs with neural network can be traced back to 1990s [13, 14, 15].
Together with the deep learning revolution, solving PDEs with neural networks
also enter a period of prosperity. In a neural network-based PDE solver, the
loss function and network structure are two key ingredients.

Regarding the loss function, one natural choice is the residual of PDE. In [8,
6], $L_{2}$ norm of the residual is used as the loss function. For elliptic
equations, the variation form provides another choice of the loss function. Yu
and E proposed to use Ritz variational form as the loss function in [7] and
Galerkin variational form was formulated as an adversarial problem in [16]. In
[17, 18], to avoid high order derivatives in the loss function, high order
PDEs are first transformed to first-order PDEs system by introducing auxiliary
variables. For the first-order system, we only need to compute first-order
derivatives in the loss function. To solve PDEs, boundary condition has to be
imposed properly. One simple way to enforce the boundary condition is to add
it to the loss function as a penalty term. In this approach, we must tune a
weight to balance the PDEs’ residual and boundary conditions. Usually, this
weight is crucial and subtle to get good results. The other way is to impose
the boundary condition explicitly by introducing a distance function of the
boundary [19]. Regarding that network structure, there are also many works
recently. A fully connected neural network (FCN) is one of the most frequently
used networks. In [7], it is found out that residual neural network (ResNet)
gives better results. For PDEs with multiscale structure, a multiscale neural
network was designedspecifically by introducing multiscale structure in the
network [20]. Activation function is another important part of neural
networks. Choice of the activation function is closely related to the
smoothness of the neural network. To compute high-order derivatives, smooth
activation functions, sigmoid, tanh, etc., are often used in PDE solvers. ReLU
activation which is used most often in machine learning is hardly used due to
its poor regularity. For special PDEs, other activation functions are also
used, such as sReLU [20], sine function[21]. By constrast, our BINet adopts an
explicit integral representation of the solution, therefore the output
satisfies the original PDE automatically.

Another related research is to learn the solution operator, i.e. map from the
parameter space to the solution space. Both the parameter space and solution
space may be infinitely dimensional. Therefore, learning solution operator is
more challenging than solving a single PDE. Solution operators may be
complicated also, and network architecture becomes more important. In [22,
23], while Green’s function and Fourier transform are used respectively to
design good network architecture, the purpose is not for solving PDEs and thus
different from ours. The network solving the single PDE can also be
generalized to learn solution operators [24, 25, 26]. In [27, 28, 29], a
neural network is used to solve the PDEs with uncertainty.

##  3 Boundary Integral Network (BINet)

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded domain, $\overline{\Omega}$ be
the closure of $\Omega$ and
$\Omega^{c}=\mathbb{R}^{d}\backslash\overline{\Omega}$. We consider the PDE in
the following form,

$\displaystyle\mathcal{L}u(x)=0.$ (3.1)

In this paper, $\mathcal{L}$ is chosen to be Laplace operator
$-\Delta=-\sum_{i=1}^{d}\frac{\partial^{2}}{\partial x_{i}^{2}}$ or Helmholtz
operator $-\Delta-k^{2}$. But in general, BINet can be applied as long as the
fundamental solution of $\mathcal{L}$ in $\mathbb{R}^{d}$ can be obtained. We
list more options of $\mathcal{L}$ in Appendix. We consider both interior
problems and exterior problems. In interior and exterior problems, the PDE
$\mathcal{L}u(x)=0$ is defined in $\Omega$ and $\Omega^{c}$ respectively.

In this paper, we consider the Dirichlet type of boundary condition
$u|_{\partial\Omega}=g(x).$ Other types of boundary conditions can be easily
handled in BINet with a small modification of the boundary integral equation.

###  3.1 Potential Theory

In this subsection, we briefly introduce the basics of the potential theory
which provides the theoretical foundation of BINet. We recall an important
theorem in potential theory [11].

###### Theorem 3.1

For any continuous function $h$ defined on $\partial\Omega$, the single layer
potential is defined as

$\mathcal{S}[h](x):=-\int_{\partial\Omega}G(x,y)h(y)ds_{\bm{y}},$ (3.2)

and the double layer potential is defined as

$\mathcal{D}[h](x):=-\int_{\partial\Omega}\frac{\partial
G(x,y)}{\partial\bm{n}_{\bm{y}}}h(y)ds_{\bm{y}}.$ (3.3)

with $\bm{n}_{\bm{y}}$ denotes out normal of $\partial\Omega$ at $\bm{y}$,
$G(x,y)$ is the fundamental solution of equation (3.1). Then, both single
layer potential and double layer potential satisfy (3.1). And for all
$x_{0}\in\partial\Omega$, we have

$\displaystyle\lim_{x\to x_{0}}\mathcal{S}[h](x)$
$\displaystyle=\mathcal{S}[h](x_{0}),$ (3.4) $\displaystyle\lim_{x\to
x_{0}^{\pm}}\mathcal{D}[h](x)$
$\displaystyle=\mathcal{D}[h](x_{0})\mp\frac{1}{2}h(x),$

where $x\to x_{0}^{-}$ and $x\to x_{0}^{+}$ mean converging in $\Omega$ and
$\Omega^{c}$ respectively.

For many important PDEs, fundamental solutions can be written explicitly. For
the Laplace equation $-\Delta u(x)=0$ in $\mathbb{R}^{2}$, the fundamental
solution is $G(x,y)=-\frac{1}{2\pi}\text{ln}|x-y|$, while the fundamental
solution for the Helmholtz equation $-\Delta u(x)-k^{2}u(x)=0$ in
$\mathbb{R}^{2}$ is $G(x,y)=\frac{i}{4}H_{0}^{1}(k|x-y|)$ where $H^{1}_{0}$ is
the Hankel function. For the Laplace equation and the Helmholtz equation in
the high dimensional case and more equations, please refer to Appendix.

Based on Theorem 3.1, the single/double layer potential (3.2) (3.3) give
explicit integral representations for the solution of the PDE. Using these
integral representations, we can construct a network such that the output of
the network solves the PDE automatically even with random initialization. This
is also the main observation in BINet.

###  3.2 The Structure of BINet

In this subsection, we will explain how to use the boundary integral form
$\mathcal{S}[h](x,\theta)$ or $\mathcal{D}[h](x,\theta)$ to construct the
structure of BINet. As shown in Fig. 1, BINet consists of three components:
_input, approximation, integration_.

  * •

From the integral formula of the single/double layer potential, it is clear
that BINet has three inputs: point in the computational domain $x\in\Omega$,
differential operator $\mathcal{L}_{\alpha}$, and domain boundary
$\partial\Omega_{\beta}$. Differential operator$\mathcal{L}_{\alpha}$
determines the fundamental solution $G$ and domain boundary
$\partial\Omega_{\beta}$ gives the domain of the integral.

  * •

In the single/double layer potential, only a density function $h$ is unknown.
In BINet, the density function $h$ is approximated using a multilayer
perceptron (MLP) (or a residual network, a.k.a, ResNet) denoted as
$h(y,\theta)$ with the learning parameter $\theta$. Note that $h$ is defined
on the boundary only.

  * •

Compute single or double layer potential in Theorem 3.1 by kernel integration
of the density function on the boundary $\partial\Omega_{\beta}$ where the
kernel is given by the explicit fundamental solution $G$. The integration can
be done numerically by the methods shown in [30, 31].

Figure 1: In BINet, fundamental solution $G(x,y,\alpha)$ explicitly depends on
the operator $\mathcal{L}_{\alpha}$ of the equation, while the density
function $h(y,\alpha,\beta,\theta)$ is implicitly dependent on
$\mathcal{L}_{\alpha}$ and the boundary $\partial\Omega_{\beta}$. This
implicit dependence is approximated by a neural network whose input is
$\mathcal{L}_{\alpha}$ and $\partial\Omega_{\beta}$. These two parts are
multiplied together and integrated on the boundary $\partial\Omega_{\beta}$,
giving the output of BINet. The boundary condition is taken as the supervisory
signal for the loss.

To train BINet, the loss function is given by (3.4) in Theorem 3.1.

$L(\theta)=\begin{cases}\|\mathcal{S}[h(\ \cdot\
;\theta)](x)-g(x)\|^{2}_{\partial\Omega},&\text{ single layer
potential}\\\\[5.0pt] \|(\frac{1}{2}\mathcal{I}+\mathcal{D})[h(\ \cdot\
;\theta)](x)-g(x)\|^{2}_{\partial\Omega},&\text{ double layer potential
(Interior problem)}\\\\[5.0pt] \|(-\frac{1}{2}\mathcal{I}+\mathcal{D})[h(\
\cdot\ ;\theta)](x)-g(x)\|^{2}_{\partial\Omega},&\text{ double layer potential
(Exterior problem)}\end{cases}$ (3.5)

where $\mathcal{S}$ and $\mathcal{D}$ are the potential operators defined in
Theorem 3.1, and $\mathcal{I}$ is the identity operator.

In BINet, the differential operator $\mathcal{L}_{\alpha}$ and the
computational domain boundary $\partial\Omega_{\beta}$ are naturally
incorporated, which means that BINet has the capability to learn the map from
the differential operator and computational domain to solutions.

##  4 Convergence Analysis of BINet

In recent years, many efforts have been devoted to the development of the
convergence theory for the over-parameterized neural networks. In [12], a
neural tangent kernel (NTK) is proposed to prove the convergence, and this
tangent kernel is also implicit in these works [32, 33, 34]. Later, a non-
asymptotic proof using NTK is given in [35]. It is shown that a sufficiently
wide network that has been fully trained is indeed equivalent to a kernel
regression predictor. In this work, we give a non-asymptotic proof of the
convergence for our BINet.

In BINet, the density function in the boundary integral form is approximated
by a neural network as $h(y,\theta)$. And a boundary integral operator is
performed on the density function, giving the output of BINet on the boundary
as $v(x)=\mathcal{A}[h](x,\theta),x\in\partial\Omega$. Here
$\mathcal{A}=\mathcal{S}$ for the single layer potential and
$\mathcal{A}=\pm\mathcal{I}/2+\mathcal{D}$ for the double layer potential of
the interior problem or the exterior problem. For simplicity, we denote
$\mathcal{A}[f](x)=\int_{\partial\Omega}\tilde{G}(x,y)f(y)dy,x\in\partial\Omega$
as the output of BINet limited on the boundary. And the loss is given by the
difference between the output and the boundary values, see Section 3 for
detail. Due to the operator $\mathcal{A}$, the convergence analysis of this
structure is non-trivial.

In the learning process, the evolution of the difference between the output
and the boundary value obeys the following ordinary differential equation

$\begin{split}\frac{d}{dt}(v(x,\theta(t))-\tilde{v}(x))&=-\int_{\partial\Omega}(v(x^{\prime},\theta(t))-\tilde{v}(x^{\prime}))\mathcal{N}_{t}(x,x^{\prime})dx^{\prime}\end{split}$
(4.1)

where $v(x,\theta)=\mathcal{A}[h](x,\theta)$ is the output of BINet limited on
the boundary and $\tilde{v}(x)$ is the boundary value, i.e., the label
function. For a detailed derivation of (4.1), see Appindex. Here
$\mathcal{N}_{t}(x,x^{\prime})=\sum_{\theta_{p}}\mathcal{A}[\frac{\partial
h}{\partial\theta_{p}(t)}](x)\mathcal{A}[\frac{\partial
h}{\partial\theta_{p}(t)}](x^{\prime})$ is the kernel at training-step index
$t$, with an admissible operator $\mathcal{A}$, see Appendix for detail.

In the following two theorems, we would show that the kernel in (4.1)
converges to a constant kernel independent of $t$ when the width of the layers
goes to infinity. And the proof sketch is listed in the Appendix based on the
works in [35].

###### Theorem 4.1

(Convergence result of kernel at initialization) Fix $\epsilon>0$ and
$\delta\in(0,1)$. Suppose the activation nonlinear function
$\sigma(\cdot)=\max\\{\cdot,0\\}$ is ReLU, the minimum width of the hidden
layer
$\text{min}_{l\in[L]}d_{l}\geq\Omega(\frac{L^{6}}{\epsilon^{4}}\log(L/\delta))$,
and the operator $\mathcal{A}$ is bounded with $\|\mathcal{A}\|_{\infty}\leq
A$. Then for the normalized data $x$ and $x^{\prime}$ where $\|x\|\leq 1$ and
$\|x^{\prime}\|\leq 1$, with probability at least $1-\delta$ we have

$|\mathcal{N}(x,x^{\prime})-[\mathcal{A}\Theta^{(L)}\mathcal{A}](x,x^{\prime})|\leq(L+1)A^{2}\epsilon.$

Here $[\mathcal{A}\Theta\mathcal{A}](x,x^{\prime})$ is the constant kernel of
BINet given by the neural-network kernel $\Theta(y,y^{\prime})$. The front and
the back operator means the operations are performed with the respect to the
first and the second variable of the neural-network kernel.

###### Theorem 4.2

(Convergence result of kernel during training) Fix
$\omega\leq\text{poly}(\frac{1}{L},\frac{1}{n},1/\log(1/\delta),\lambda_{0})$
and $\delta\in(0,1)$. Suppose that
$\text{min}_{l}\\{d_{l}\\}\geq\text{poly}(1/\omega)$, and the operator
$\mathcal{A}$ is bounded with $\|\mathcal{A}\|_{\infty}\leq A$. Then with
probability at least $1-\delta$ over Gaussian random initialization, we have
for all $t\geq 0$,

$|\mathcal{N}_{t}(x,x^{\prime})-\mathcal{N}_{0}(x,x^{\prime})|\leq
A^{2}\omega,$

where $\mathcal{N}_{t}(x,x^{\prime})$ is the kernel along time $t$ and
$\mathcal{N}_{0}(x,x^{\prime})\overset{\Delta}{=}\mathcal{N}(x,x^{\prime})$ is
the kernel when initialization over random Gaussian denoted in Theorem 4.1 to
distinguish with the training process.

Further, we have the following lemma for the positive definiteness of the new
constant kernel.

###### Lemma 4.1

$[\mathcal{A}\Theta^{(L)}\mathcal{A}](x,x^{\prime})$ is positive definite for
double layer potential in BINet. For single layer potential, the positive
definiteness depends on the $C^{\infty}$ compactness of the boundary
$\partial\Omega$.

The proof of Lemma 4.1 is given in Appendix. And the invertibility of the
operator $\mathcal{A}$ is utilized to complete the proof. [36, 37]

By Lemma 4.1, equation (4.1), Theorem 4.1 and 4.2, the error in BINet thus
vanishes for double layer potential after fully training ($t\to\infty$) under
the assumption that the the width of the neural network goes to infinity. And
for single layer potential, the convergence results depend on the boundary,
i.e., $\partial\Omega$ is $C^{\infty}$ compact. The proof of the convergence
results is in the real space, however with the complex form of the kernel
$[\mathcal{A}^{*}\Theta^{(L)}\mathcal{A}](x,x^{\prime})$ with
$\mathcal{A}^{*}[f](x)=\int_{\partial\Omega}\tilde{G}^{*}(x,y)f(y)dy$ for
Helmholtz equations, the results still hold with inner product defined in the
complex space.

##  5 Experiments

We use BINet to compute a series of examples including solving a single PDE,
where differential operator and domain geometry are fixed, and learning
solution operators. PDEs defined on both bounded and unbounded domains will be
considered. In order to estimate the accuracy of the numerical solution $u$,
the relative $L^{2}$ error $\|u-u^{\ast}\|_{2}/\|u^{\ast}\|_{2}$ is used,
where $u^{\ast}$ is the exact solution. We compare our method with two state-
of-the-art methods, the Deep Ritz method and PINN only for interior problems,
since as we claimed before, other deep-learning-based PDE solvers are not able
to handle exterior problems.

In BINet, the fully connected neural network (MLP) or residual neural network
(ResNet) are used to approximate the density function. Since there is no
regularity requirement on density function, we can use any activation
functions including ReLU. For the Laplace equation, the network only has one
output, i.e., the approximation of density $h$, while for the Helmholtz
equation, because its solution is complex, the network has two outputs, i.e.,
the real part and the imaginary part of density $h$. In the experiments, we
choose the Adam optimizer to minimize the loss function and all experiments
are run on a single GPU of GeForce RTX 2080 Ti.

###  5.1 Experimental Results on Solving One Single PDE

Laplace Equation with Smooth Boundary Condition. First, we consider a
Laplacian equation in the bounded domain,

$\displaystyle-\Delta u(x,y)=0,\ (x,y)\in\Omega,$ (5.1) $\displaystyle
u(x,y)=e^{ax}\sin(ay),\ (x,y)\in\partial\Omega,$

where $a$ is a fixed constant. We will compare the results of PINN, Deep-Ritz
method, and BINet for different $a$. For simplicity, we choose
$\Omega=[-1,1]\times[-1,1]$. In this example, we will use a residual neural
network introduced in [7]. We follow [7] to choose $\text{ReLU}^{3}$ as the
activation function in Deep Ritz method and PINN. In BINet, we use ReLU as the
activation function since BINet has less regularity requirement.

When $a=4$, for these three methods, we all selected 800 equidistant sample
points on $\partial\Omega$, and for PINN and Deep-Ritz method, we randomly
selected 1600 sample points in $\Omega$. We all use residual neural networks
with 40 neurons per layer and six blocks.

When $a=8$, for the BINet method, we selected 2000 equidistant sample points
on the boundary. For PINN and Deep-Ritz method, we randomly selected 4000
sample points in $\Omega$ and randomly selected 800 sample points on
$\partial\Omega$. We also use residual neural networks with 100 neurons per
layer and six blocks. But if we look at the solutions on
$[-0.1,0.1]\times[-1,1]$, we find that the solutions of PINN and Deep Ritz
method are quite different from the exact solution, but BINet method still
captures the subtle structure of the exact solution. The results of different
methods including PINN, Deep-Ritz method and BINet for $a=8$ are shown in
Figure 2.

Table 1: Relative $L^{2}$ error of equation (5.1) with different methods. | PINN | Deep Ritz | BINet
---|---|---|---
$a=4$ | 0.0140 | 0.0952 | 0.0031
$a=8$ | 0.0262 | 0.2194 | 0.0002

After training for 20000 epochs, the relative $L^{2}$ error of these methods
is shown in the table 1. In this example, with the same number of layers and
neurons, BINet is always better than the other two methods no matter what the
value of $a$ is. When $a$ increases, unlike other methods, the result of the
BINet does not get worse.

Figure 2: The solutions of Laplace equation (5.1) with smooth boundary
conditions for $a=8$ by PINN, Deep-Ritz method and BINet. The first row shows
the exact solutions and the numerical solutions obtained by the three methods
in $\Omega=[-1,1]\times[-1,1]$. The second row is the zoom-in of the above
figures in the subdomain $[-0.1,0,1]\times[-1,1]$. Only BINet captures the
subtle structure of the exact solution successfully.

Laplace Equation with Non-smooth Boundary Condition. Next, let’s consider a
Laplace equation with a nonsmooth boundary condition. We also assume the
domain $\Omega=[-1,1]\times[-1,1]$ and the boundary value problem is

$\displaystyle-\Delta u(x,y)=0,\ (x,y)\in\Omega,$ (5.2) $\displaystyle
u(x,y)=2-|x|-|y|,\ (x,y)\in\partial\Omega.$

In problem (5.2), the boundary condition is not smooth. In this example, we
also used the ResNet with six blocks and 40 neurons per layer for three
methods. We selected equidistant 800 sample points on $\partial\Omega$ for
three methods, and for PINN and Deep-Ritz method, we randomly selected 1000
sample points in $\Omega$. Figure 3 shows the results of different methods. In
this example, we take the result of the finite difference method with high
precision mesh as the exact solution.

Figure 3: The solutions of Laplace equation (5.2) with non-smooth boundary
conditions by high-precision finite difference as the ground truth, PINN, deep
Ritz and BINet. The red box is the zoom-in of the vicinity of the non-smooth
point (0,1) in each figure with the same scale. Only BINet learns the
singularity on the boundary successfully.

From Figure 3, we can find that for PINN and Deep Ritz methods, the solutions
on the boundary are smooth, which are different from the boundary condition.
However, the boundary condition is well approximated by the solution of the
BINet method. The reason is that, to satisfy the interior smoothness of the
solution, the neural network of the PINN and Deep Ritz methods have to be a
smooth function. So the solutions are still smooth even if they are close
enough to the unsmooth boundary points.

Helmholtz Equation with Different Wavenumbers. In this experiment, we consider
an interior Helmholtz equation

$\displaystyle-\Delta u(x,y)-k^{2}u(x,y)=0,\ (x,y)\in\Omega,$ (5.3)
$\displaystyle u(x,y)=e^{i(k_{1}x+k_{2}y)},\ (x,y)\in\partial\Omega,$

where $(k_{1},k_{2})=(k\cos\frac{\pi}{7},k\sin\frac{\pi}{7})$, and
$\partial\Omega=\\{(\frac{9}{20}\cos(t)-\frac{\cos(5t)}{9}\cos(t),\frac{9}{20}\sin(t)-\frac{\cos(5t)}{9}\sin(t))|t\in[0,2\pi]\\}$.
The Deep-Ritz method can not solve the Helmholtz equation. Hence, we will
compare the BINet method and PINN method for different $k$. We choose a fully
connected neural network with 4 hidden layers with Sigmoid activation function
and 40 neurons per layer. and we choose 800 points on the boundary for BINet
and PINN. In addition, we also randomly selected 2400 sample points in
$\Omega$. For k = 1 and 4, we use the PINN type method and BINet method to
solve the equation respectively. The loss function and results are shown in
Figure 4. We can see the loss function of BINet descends faster, and for
$k=4$, the loss of the PINN method does not converge. In contrast, the loss of
the BINet is always convergent no matter the value of $k$ is. The second and
the third figures also show the result of BINet is much better than PINN.

Figure 4: The experiments of solving the Helmholtz equation (5.3) by PINN and
BINet. Deep Ritz is unable to solve the Helmholtz equation thus not exhibited.
The first figure shows loss functions of different cases. Since for $k\geq 4$,
PINN also fails, we show the errors when $k=1$. The errors of solutions solved
by PINN method and BINet method are shown in the second and third figures,
respectively.

###  5.2 Experimental Results on Solution Operators

The Operator from Equation Parameters to Solutions. In this example, we
consider the Helmholtz equations with variable wavenumber $k$.

$\displaystyle-\Delta u(x,y)-k^{2}u(x,y)=0,\ (x,y)\in\Omega,$ (5.4)
$\displaystyle u(x,y)=H_{0}^{1}(k\sqrt{x^{2}+y^{2}}),\
(x,y)\in\partial\Omega,$

In the training phase, we set $k\in[2,3.5]\cup[4.5,6]$. We also use double
layer potential to construct the loss function, and after 5000 training
epochs, we show the relative $L^{2}$ error versus the wavenumber $k$ in Figure
5. From the first figure, the relative $L^{2}$ error is about $10^{-3}$ or
$10^{-2}$. Compared with solving a single equation, the relative $L^{2}$ error
is still small. The relative error increases slightly with the increase of
$k$, which is because the Helmholtz equation becomes more difficult to solve
when $k$ increases. This means that we have successfully learned the operator
mapping of exterior parametric PDE problems on an unbounded domain. Most
importantly, although $k$ is not selected between $[3.5,4.5]$ during training,
the relative error is still small on the test when we take values in the
interval $[3.5,4.5]$. This shows that our method has good generalization
ability.

Figure 5: The results of the learning operators mapping the equation
parameters $k$ (wavenumbers) to the solutions of Helmholtz equations (5.4).
The first figure shows the relative $L^{2}$ error of the solution with
different wavenumber $k$, which shows BINet has successfully learned the
solution operator and has generalization capability. The last three figures
show the absolute error of the solutions mapping from different wavenumber
$k$, i.e., $k=3,4,6.2$, respectively.

The Operator from Boundary Geometry to Solutions. In this example, we consider
a Laplace equation with parametric boundaries. The problem is

$\displaystyle-\Delta u(x,y)$ $\displaystyle=0,\ $
$\displaystyle(x,y)\in\Omega_{\beta},$ (5.5) $\displaystyle u(x,y)$
$\displaystyle=g(x,y;\beta),\ $ $\displaystyle(x,y)\in\partial\Omega_{\beta},$

where the boundary condition
$g(x,y;\beta)=(x-x_{bc})(y-y_{bc})+(x-x_{bc})+(y-y_{bc})+1$, and
$(x_{bc},y_{bc})$ is the barycenter of the $\Omega_{\beta}$. We assume that
$\Omega_{\beta}$ can take any triangle in a domain. For simplicity, We can fix
one vertex at the origin and one edge on the positive half x-axis, while the
third vertex is in the first quadrant by translation and rotation. Then we can
assume the vertex is $(0,0),(a,0),\text{and}(b,c)$. In this example, we assume
$a,b,c$ can take any value in interval $[0,1]$. In this example, we choose a
ResNet with eight blocks, and 100 neurons per layer. Single potential layer is
used to calculate the boundary integral. We randomly selected 80 triangles to
calculate the loss function, and after every 500 epochs, triangles will be
randomly selected again. After training for 5000 epochs, we randomly choose
two triangles, and the solutions of the each triangle by BINet method has
shown in figure 6. The relative $L^{2}$ error is about $10^{-3}$. From this,
we can see BINet has successfully learned the operator from boundary geometry
to solution.

Figure 6: The results of learning operators from boundary geometry to
solutions of the equation (5.5). This histogram of relative $L^{2}$ error of
the solutions mapping from different triangle boundaries are shown in the
first figure. One hundred triangle boundaries are randomly selected to give
the distribution of the relative error. It is easy to check that more than
half of the errors are less than $1\%$. For a typical triangle domain given by
vertices $(0,0)$, $(0.3,0)$, $(0.6,0.4)$, the solution and its absolute error
are shown in the right two figures respectively.

##  6 Conclusion

We have developed a new neural network method called BINet to solve PDEs. In
BINet, the solution of PDE is represented by boundary integral composed of an
explicit kernel and an unknown density which is approximated by a neural
network. Then the PDE is solved by learning the boundary integral
representation to fit the boundary condition. Since the loss function measures
only the misfit between the integral representation and the boundary
condition, BINet has less hyper-parameters and lower sampling dimensions than
many other neural network-based PDE solvers. Because the boundary integral
satisfies PDE automatically in the interior and exterior of the boundary,
BINet can solve bounded and unbounded PDEs. Furthermore, BINet can learn
operators from PDE parameters including coefficients and boundary geometry to
solutions. Besides, using the NTK technique, we prove that BINet converges as
the width of the network goes to infinity. We test BINet with the Laplace
equation and Helmholtz equation in extensive settings. The numerical
experiments show that BINet works effectively for many cases such as interior
problems, exterior problems, high wavenumber problems. The experiments also
illustrate the capability of BINet in learning solution operators. All the
experiments verify the advantages of BINet numerically. Although our method
exhibits competitive performance against the PINN method and DeepRitz method
in many situations, the requirement of high-precision boundary integration
limits further applications in higher-dimensional problems. This will be the
direction of improving BINet in the future.

## References

  * [1] David J Griffiths.  Introduction to electrodynamics, 2005.
  * [2] Roger Temam.  Navier-Stokes equations: theory and numerical analysis, volume 343\.  American Mathematical Soc., 2001.
  * [3] Erwin Schrödinger.  An undulatory theory of the mechanics of atoms and molecules.  Physical review, 28(6):1049, 1926.
  * [4] James D MacBeth and Larry J Merville.  An empirical examination of the black-scholes call option pricing model.  The journal of finance, 34(5):1173–1186, 1979.
  * [5] Ian Goodfellow, Yoshua Bengio, Aaron Courville, and Yoshua Bengio.  Deep learning, volume 1.  MIT press Cambridge, 2016.
  * [6] Justin Sirignano and Konstantinos Spiliopoulos.  Dgm: A deep learning algorithm for solving partial differential equations.  Journal of computational physics, 375:1339–1364, 2018.
  * [7] E Weinan and Bing Yu.  The deep ritz method: a deep learning-based numerical algorithm for solving variational problems.  Communications in Mathematics and Statistics, 6(1):1–12, 2018.
  * [8] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [9] Sifan Wang, Xinling Yu, and Paris Perdikaris.  When and why pinns fail to train: A neural tangent kernel perspective.  arXiv preprint arXiv:2007.14527, 2020.
  * [10] Quanhui Zhu and Jiang Yang.  A local deep learning method for solving high order partial differential equations.  arXiv preprint arXiv:2103.08915, 2021.
  * [11] Oliver Dimon Kellogg.  Foundations of potential theory, volume 31.  Courier Corporation, 1953.
  * [12] Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  arXiv preprint arXiv:1806.07572, 2018.
  * [13] MWMG Dissanayake and Nhan Phan-Thien.  Neural-network-based approximations for solving partial differential equations.  communications in Numerical Methods in Engineering, 10(3):195–201, 1994.
  * [14] Isaac E Lagaris, Aristidis Likas, and Dimitrios I Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE transactions on neural networks, 9(5):987–1000, 1998.
  * [15] Isaac E Lagaris, Aristidis C Likas, and Dimitris G Papageorgiou.  Neural-network methods for boundary value problems with irregular boundaries.  IEEE Transactions on Neural Networks, 11(5):1041–1049, 2000.
  * [16] Yaohua Zang, Gang Bao, Xiaojing Ye, and Haomin Zhou.  Weak adversarial networks for high-dimensional partial differential equations.  Journal of Computational Physics, 411:109409, 2020.
  * [17] Zhiqiang Cai, Jingshuang Chen, Min Liu, and Xinyu Liu.  Deep least-squares methods: An unsupervised learning-based numerical method for solving elliptic pdes.  Journal of Computational Physics, 420:109707, 2020.
  * [18] Liyao Lyu, Zhen Zhang, Minxin Chen, and Jingrun Chen.  Mim: A deep mixed residual method for solving high-order partial differential equations.  arXiv preprint arXiv:2006.04146, 2020.
  * [19] Jens Berg and Kaj Nyström.  A unified deep artificial neural network approach to partial differential equations in complex geometries.  Neurocomputing, 317:28–41, 2018.
  * [20] Wei Cai and Zhi-Qin John Xu.  Multi-scale deep neural networks for solving high dimensional pdes.  arXiv preprint arXiv:1910.11710, 2019.
  * [21] Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  Advances in Neural Information Processing Systems, 33, 2020.
  * [22] Craig R Gin, Daniel E Shea, Steven L Brunton, and J Nathan Kutz.  Deepgreen: Deep learning of green’s functions for nonlinear boundary value problems.  arXiv preprint arXiv:2101.07206, 2020.
  * [23] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Fourier neural operator for parametric partial differential equations.  arXiv preprint arXiv:2010.08895, 2020.
  * [24] Shengze Cai, Zhicheng Wang, Lu Lu, Tamer A Zaki, and George Em Karniadakis.  Deepm&mnet: Inferring the electroconvection multiphysics fields based on operator approximation by neural networks.  Journal of Computational Physics, 436:110296, 2021.
  * [25] Han Gao, Luning Sun, and Jian-Xun Wang.  Phygeonet: Physics-informed geometry-adaptive convolutional neural networks for solving parametric pdes on irregular domain.  arXiv preprint arXiv:2004.13145, 2020.
  * [26] Lu Lu, Pengzhan Jin, and George Em Karniadakis.  Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators.  arXiv preprint arXiv:1910.03193, 2019.
  * [27] Yuehaw Khoo, Jianfeng Lu, and Lexing Ying.  Solving parametric pde problems with artificial neural networks.  arXiv preprint arXiv:1707.03351, 2017.
  * [28] Yibo Yang and Paris Perdikaris.  Physics-informed deep generative models.  arXiv preprint arXiv:1812.03511, 2018.
  * [29] Yibo Yang and Paris Perdikaris.  Adversarial uncertainty quantification in physics-informed neural networks.  Journal of Computational Physics, 394:136–152, 2019.
  * [30] Bradley K Alpert.  Hybrid gauss-trapezoidal quadrature rules.  SIAM Journal on Scientific Computing, 20(5):1551–1584, 1999.
  * [31] Sharad Kapur and Vladimir Rokhlin.  High-order corrected trapezoidal quadrature rules for singular functions.  SIAM Journal on Numerical Analysis, 34(4):1331–1356, 1997.
  * [32] Simon Du, Jason Lee, Haochuan Li, Liwei Wang, and Xiyu Zhai.  Gradient descent finds global minima of deep neural networks.  In International Conference on Machine Learning, pages 1675–1685. PMLR, 2019.
  * [33] Simon S Du, Xiyu Zhai, Barnabas Poczos, and Aarti Singh.  Gradient descent provably optimizes over-parameterized neural networks.  arXiv preprint arXiv:1810.02054, 2018.
  * [34] Yuanzhi Li and Yingyu Liang.  Learning overparameterized neural networks via stochastic gradient descent on structured data.  arXiv preprint arXiv:1808.01204, 2018.
  * [35] Sanjeev Arora, Simon S Du, Wei Hu, Zhiyuan Li, Ruslan Salakhutdinov, and Ruosong Wang.  On exact computation with an infinitely wide neural net.  arXiv preprint arXiv:1904.11955, 2019.
  * [36] Wenjie Gao.  Layer potentials and boundary value problems for elliptic systems in lipschitz domains.  Journal of Functional Analysis, 95(2):377–399, 1991.
  * [37] Gregory Verchota.  Layer potentials and regularity for the dirichlet problem for laplace’s equation in lipschitz domains.  Journal of functional analysis, 59(3):572–611, 1984.
  * [38] George C Hsiao and Wolfgang L Wendland.  Boundary integral equations.  Springer, 2008.
  * [39] Lexing Ying.  Fast algorithms for boundary integral equations.  Multiscale Modeling and Simulation in Science, pages 139–193, 2009\.
  * [40] Zeyuan Allen-Zhu, Yuanzhi Li, and Yingyu Liang.  Learning and generalization in overparameterized neural networks, going beyond two layers.  arXiv preprint arXiv:1811.04918, 2018.
  * [41] Jiaoyang Huang and Horng-Tzer Yau.  Dynamics of deep neural networks and neural tangent hierarchy.  In International Conference on Machine Learning, pages 4542–4551. PMLR, 2020.
  * [42] Zeyuan Allen-Zhu, Yuanzhi Li, and Zhao Song.  A convergence theory for deep learning via over-parameterization.  In International Conference on Machine Learning, pages 242–252. PMLR, 2019.
  * [43] Yann A LeCun, Léon Bottou, Genevieve B Orr, and Klaus-Robert Müller.  Efficient backprop.  In Neural networks: Tricks of the trade, pages 9–48. Springer, 2012\.
  * [44] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In Proceedings of the IEEE international conference on computer vision, pages 1026–1034, 2015.

## Appendix

##  Appendix A A review of the PINN and Deep Ritz method

###  A.1 PINN method

To solve the linear PDE

$\displaystyle\mathcal{L}u(x)=0,\ x\in\Omega,$ (A.1) $\displaystyle
u(x)=g(x),\ x\in\partial\Omega,$

the main idea of the PINN[8] method is to use a neural network $u(x;\theta)$
as an ansatz to approximate the solution $u(x)$, where $\theta$ represents the
trainable parameters in the neural network. There was other work of the
similar idea such as [19, 15, 6]. Then we can use the automatic
differentiation tool to calculate the derivative $\mathcal{L}u(x;\theta)$ and
define the loss function

$L_{1}(\theta)=\|\mathcal{L}u(x;\theta)\|^{2}_{\Omega}.$

For the boundary conditions, we can define the loss function

$L_{2}(\theta)=\|u(x;\theta)-g(x)\|^{2}_{\partial\Omega}.$

Finally, we can combine the loss function $L_{1}$ and $L_{2}$ with a hyper-
parameter $\beta$ to get loss function,

$L(\theta)=L_{1}+\beta
L_{2}=\|\mathcal{L}u(x;\theta)\|^{2}_{\Omega}+\beta\|u(x;\theta)-g(x)\|^{2}_{\partial\Omega}.$

By minimizing the loss function $L$, PINN will get the approximation solution
of the PDE (A.1).

###  A.2 Deep Ritz method

For the specific PDE problems in equation (A.1), we can change the equation
into a Ritz variational form. This is the main idea of the Deep Ritz
method[7]. For instance, if we consider a Laplace equation,

$\displaystyle-\Delta u(x)=0,\ x\in\Omega,$ (A.2) $\displaystyle u(x)=g(x),\
x\in\partial\Omega,$

we can solve the equation equivalently by minimizing the following Ritz
variational problem

$\displaystyle\int_{\Omega}\frac{1}{2}|\nabla u(x)|^{2}dx.$ (A.3)

We also use a neural network $u(x;\theta)$ to approximate the solution of the
PDE, and we can use the automatic differentiation tool to calculate the
gradient $\nabla u(x;\theta)$ of the neural network. So the variation (A.3)
can naturally be used as a loss function, defined as

$L_{1}(\theta)=\int_{\Omega}\frac{1}{2}|\nabla u(x;\theta)|^{2}dx.$

For the boundary condition, the loss function $L_{2}$ also can be defined as

$L_{2}(\theta)=\|u(x;\theta)-g(x)\|^{2}_{\partial\Omega}.$

Finally, the loss funtion can be defined as

$L(\theta)=L_{1}+\beta L_{2}=\int_{\Omega}\frac{1}{2}|\nabla
u(x;\theta)|^{2}dx+\beta\|u(x;\theta)-g(x)\|^{2}_{\partial\Omega},$ (A.4)

where $\beta$ is also the hyper-parameter. By minimize the loss function
$L(\theta)$, Deep Ritz method will get the solution of the PDE (A.2).

This paper introduce a residual network as the anastz to approximate the
solution. The residual network is also used in our work to approximate the
density function in the boundary integral form. The architecture of the
residual network is shown in Figure 7.

Figure 7: The architecture of the residual network.

##  Appendix B The fundamental solution of different equations

In this section, we make some supplementary introductions to the fundamental
solution. Before defining the fundamental solution, we first introduce the
$\delta$-function,

###### Definition B.1

A function $\delta(x)$ called n-dimensional $\delta$-function if
$\delta(x)=\left\\{\begin{array}[]{cc}0,&x\neq 0\\\
\infty,&x=0,\end{array}\right.$, and for all functions $f$ that is continuous
at $a$ we have $\int_{\mathbb{R}^{n}}f(x)\delta(x-a)dx=f(a).$

Then, for the PDE

$\mathcal{L}u(x)=0,\ x\in\Omega,$ (B.1)

we can define the corresponding fundamental solution $G$ of the equations
(B.1).

###### Definition B.2

A function $G(x,y)$ is called the fundamental solution corresponding to
equations (B.1) if $G(x,y)$ is symmetric about $x$ and $y$ and $G(x,y)$
satisfy

$\mathcal{L}_{y}G(x,y)=\delta(|y-x|),$

where $(x,y)\in\mathbb{R}^{n}\times\mathbb{R}^{n}$ and $L_{y}$ is the
differential operator $L$ which acts on component $y$.

Limited by the length of the article, although we only introduce the
fundamental solutions of Laplace equations and Helmholtz equations in two-
dimensional cases in detail, in general, BINet can be applied as long as the
fundamental solution of $\mathcal{L}$ in $\mathbb{R}^{d}$ can be obtained.
Let’s give a few more examples. More details can be found in [38, 11, 39].

###  B.1 The Laplace Equations

If we consider a Laplace equation

$-\Delta u(x)=0,\ x\in\Omega,$ (B.2)

the fundamental solution has the following form

$G(x,y)=\left\\{\begin{array}[]{cc}-\frac{1}{2\pi}ln|x-y|&n=2\\\\[5.0pt]
\frac{1}{(n-2)w_{n}}\frac{1}{|x-y|^{n-2}}&n\geq 3,\end{array}\right.$

where $n$ is the dimension of the equation and $w_{n}$ is the volume of the
n-dimensional unit sphere. Then the fundamental solution $G$ satisfies

$-\Delta_{y}G(x,y)=\delta(|x-y|).$

###  B.2 The Helmholtz Equations

The Helmholtz equation has the following form

$-\Delta u(x)-k^{2}u(x)=0,\ x\in\Omega,$ (B.3)

where $k$ is a real number. The fundamental solution of the Helmholtz equation
has the following form

$G(x,y)=\left\\{\begin{array}[]{cc}\frac{i}{4}H_{0}^{1}(k|x-y|)&n=2\\\\[1.99997pt]
\frac{1}{(n-2)w_{n}}\frac{e^{ik|x-y|}}{|x-y|^{n-2}}&n\geq
3,\end{array}\right.$

where $n$ is the dimension of the equation and $w_{n}$ is also the volume of
the n-dimensional unit sphere. Then the fundamental solution $G$ satisfies

$-\Delta_{y}G(x,y)-k^{2}G(x,y)=\delta(|x-y|).$

###  B.3 The Navier’s Equations

We consider Navier’s equations (also called Lam$\acute{e}$ system). These are
famous equations in linear elasticity for isotropic materials, and the
governing equations are

$\displaystyle-\mu\Delta u(x)-(\lambda+\mu)\nabla(\nabla\cdot u(x))=0,$ (B.4)

where $\lambda,\mu>0$ are the Lam$\acute{e}$ constants of the elastic
material, and $u(x)\in\mathbb{R}^{n}$ is the displacement vector. The
fundamental solution of the equation (B.4) is

$G(x,y)=\left\\{\begin{aligned}
\frac{\lambda+3\mu}{8\pi\mu(\lambda+2\mu)}\left[\frac{1}{|x-y|}\mathbf{I}_{3}\right.\\\
\left.+\frac{\lambda+\mu}{\lambda+3\mu}\frac{1}{|x-y|^{2}}(x-y)(x-y)^{\mathrm{T}}\right],&n=3,\\\
\frac{\lambda+3\mu}{4\pi\mu(\lambda+2\mu)}\left[\ln\frac{1}{|x-y|}\mathbf{I}_{2}\right.\\\
\left.+\frac{\lambda+\mu}{\lambda+3\mu}\frac{1}{|x-y|^{2}}(x-y)(x-y)^{\mathrm{T}}\right],&n=2.\end{aligned}\right.$
(B.5)

It means that G(x,y) defined by the (B.5) satisfies the following equation,

$-\mu\Delta_{y}G(x,y)-(\lambda+\mu)\nabla_{y}(\nabla_{y}\cdot
G(x,y))=\delta(|x-y|)\mathbf{I}_{n},$ (B.6)

where $\mathbf{I}_{n}$ is the n-order identity matrix.

###  B.4 The Stokes Equations

Stokes equations are well known in the incompressible viscous fluid model. The
general form of the Stokes equations is

$\displaystyle-\mu\Delta u(x)+\nabla p(x)$ $\displaystyle=f(x),$ (B.7)
$\displaystyle\nabla\cdot u(x)$ $\displaystyle=0,\
x\in\Omega\subset\mathbb{R}^{n},$

where $u$ and $p$ are the velocity and pressure of the fluid flow,
respectively, and $\mu$ and $f$ are the given dynamic viscosity of the fluid
and forcing term, respectively.

For n=2 the fundamental solutions of the (B.10) are

$\displaystyle\bm{v}^{k}(x,y)$
$\displaystyle=\frac{1}{4\pi\mu}\left\\{\log\frac{1}{|x-y|}\bm{e}^{k}+\sum_{j=1}^{2}\frac{\left(x_{k}-y_{k}\right)\left(x_{j}-y_{j}\right)\bm{e}^{j}}{|x-y|^{2}}\right\\}$
(B.8) $\displaystyle q^{k}(x,y)$ $\displaystyle=\frac{\partial}{\partial
x_{k}}\left\\{-\frac{1}{2\pi}\log\frac{1}{|x-y|}\right\\},$

and for n=3 the fundamental solutions of the (B.10) are

$\displaystyle\bm{v}^{k}(x,y)$
$\displaystyle=\frac{1}{8\pi\mu}\left\\{\frac{1}{|x-y|}\bm{e}^{k}+\sum_{j=1}^{3}\frac{\left(x_{k}-y_{k}\right)\left(x_{j}-y_{j}\right)\bm{e}^{j}}{|x-y|^{3}}\right\\}$
(B.9) $\displaystyle q^{k}(x,y)$ $\displaystyle=\frac{\partial}{\partial
x_{k}}\left\\{-\frac{1}{4\pi}\frac{1}{|x-y|}\right\\},$

where $k=1,\cdots,n$ and $e^{k}$ denotes the unit vector along the
$x_{k}$-axis. $v^{k}$ and $q^{k}$ satisfy

$\displaystyle-\mu\Delta_{x}v^{k}(x,y)+\nabla_{x}q^{k}(x,y)$
$\displaystyle=\delta(|x-y|)e^{k},$ (B.10) $\displaystyle\nabla_{x}\cdot
v^{k}(x,y)$ $\displaystyle=0,$

where $x,y\in\mathbb{R}^{n}$.

###  B.5 The Biharmonic Equation

The Biharmonic Equation is a single scalar 4th-order equation, which can be
reduced from plane elasticity and plane Stokes flow. We consider a two
dimensional Biharmonic equation,

$\Delta^{2}u(x)=0,\ x\in\Omega\subset\mathbb{R}^{2}.$ (B.11)

The fundamental solution of the equation (B.11) is

$\displaystyle G(x,y)=\frac{1}{8\pi}|x-y|^{2}\log|x-y|,\
x,y\in\mathbb{R}^{2},$ (B.12)

where $G(x,y)$ satisfies

$\Delta_{y}^{2}G(x,y)=\delta(|x-y|).$

##  Appendix C The convergence analysis of BINet

###  C.1 The structure for solving PDEs using neural networks

BINet consists of a neural network such as MLP and an integral operator
performed on the output of the neural network. Thus, the output of BINet reads

$v(x,\theta)=\mathcal{A}[h](x,\theta),$

where $h=h(y,\theta)\in\mathbb{R}$ is the neural network approximating the
density function in the boundary integral form, $y\in\mathbb{R}^{d}$ is the
$d$-dimensional variable. The operator $\mathcal{A}$ is performed on the
output of the neural network which completes the whole architecture. And the
loss function is

$\mathcal{L}=\|v(x,\theta)-\tilde{v}(x)\|^{2}_{2},$

with label function $\tilde{v}(x)$.

For a more general setup, the operator $\mathcal{A}$ has different forms. For
PINN/DGM method, the operator is directly the partial differential operator,
implying

$\mathcal{A}[u](x,\theta)=-\Delta u(x,\theta),$

where $u(x,\theta)$ is the approximation of the solution. The Deep-Ritz method
for solving the Laplace equation $-\Delta u(x)=0$ is to minimize the
optimization problem $\min_{\theta}L(\theta)$ where part of the loss reads

$\displaystyle L_{1}(\theta)$ $\displaystyle=\int_{\Omega}\frac{1}{2}|\nabla
u(x,\theta)|^{2}dx.$

It follows that the corresponding operator has the following form

$\mathcal{A}[u](x)=\frac{1}{\sqrt{2}}\nabla u(x,\theta).$

Therefore in the view of the operator applied on the neural network, different
from the integral type operator of BINet, PINN and Deep Ritz methods have
extra differential operators although the Deep Ritz method decreases the order
from the second to the first.

###### Definition C.1

The operator $\mathcal{A}$ is admissible if the following conditions hold:

  1. 1.

$\mathcal{A}(\lambda_{1}f_{1}+\lambda_{2}f_{2})=\lambda_{1}\mathcal{A}f_{1}+\lambda_{2}\mathcal{A}f_{2}$
(linear property);

  2. 2.

$[\mathcal{A}_{x},\mathcal{A}_{x^{\prime}}]:=\mathcal{A}_{x}\mathcal{A}_{x^{\prime}}-\mathcal{A}_{x^{\prime}}\mathcal{A}_{x}=0$
(commutative property);

  3. 3.

$[\mathcal{A},\frac{\partial}{\partial\theta}]=0$ (parameter variant).

It is easy to check that the operators of PINN, Deep-Ritz, and our BINet all
satisfy the admissible property. And the admissibility is crucial in the
following proof.

The different design of the neural network and the operator makes the network
different. Here, we adopt the typical settings of the neural network as an
MLP. As the integral operator $\mathcal{A}$ is bounded, thus the convergence
results can be obtained in our BINet and the proof is shown in Appendix 3.3.
Here the structure of the neural network is introduced first for the
derivation of the NTK form.

The $L$-hidden layer MLP is defined as

input layer: $\displaystyle\quad g^{(0)}=y,$ (C.1) hidden layer:
$\displaystyle\quad
g^{(l)}(y)\overset{\Delta}{=}\sqrt{\frac{c_{\sigma}}{d_{l}}}\sigma(f^{(l)}(y)),\;\;f^{(l)}(y)=W^{(l)}g^{(l-1)}(y),\;l\in[L]$
(C.2) output layer: $\displaystyle\quad
h(y,\theta)=f^{(L+1)}=W^{(L+1)}g^{(L)}(y),$ (C.3)

where $l\in[L]\overset{\Delta}{=}\\{1,2,\cdots,L\\}$ is the hidden layer,
$\theta$ is the trainable parameters which is the standard representation for
the weights $W^{(l)}$, $d_{l}$ is the width of the $l$-th layer and
$c_{\sigma}={\bf E}_{u\sim N(0,1)}\sigma(u)$.

###  C.2 The dynamic neural tangent kernel

We have chosen the MLP as the neural network in the analysis for simplicity. A
similar analysis can also be done for other structures as the convergence
results of such neural networks are reported in the literature [40, 41].
Applying the integral operator on the neural networks should also give similar
convergence results. Thus different schemes here imply different forms of the
operator $\mathcal{A}$, see Appendix C.1 for detail.

The training process of the neural ODE is basically to minimize the loss by
the method based on the gradient. One typical scheme is the gradient descent
method which has the form

$\theta_{n+1}=\theta_{n}-\alpha\frac{\partial L}{\partial\theta}.$ (C.4)

When the learning rate $\alpha\to 0$, we have the limiting gradient flow

$\frac{d\theta}{dt}=-\frac{\partial L}{\partial\theta},$ (C.5)

where $t$ is the continuous version of index of the learning steps in the
training process. More precisely, for the weight matrix $W^{(l)}$, we have the
evolution $\frac{dW^{(l)}}{dt}=-\frac{\partial L}{\partial W^{(l)}}.$

Hence the evolution of the prediction satisfies the following form

$\displaystyle\begin{split}\frac{dv(x,\theta)}{dt}&=\sum_{\theta_{p}}\frac{\partial
v}{\partial\theta_{p}}\frac{\partial\theta_{p}}{\partial t}\\\
&=\sum_{\theta_{p}}-\mathcal{A}[\frac{\partial
h}{\partial\theta_{p}}](x)\frac{\partial L}{\partial\theta_{p}}\\\
&=\sum_{\theta_{p}}-\mathcal{A}[\frac{\partial
h}{\partial\theta_{p}}](x)\langle\frac{\partial L}{\partial
v}(x^{\prime}),\frac{\partial v}{\partial\theta_{p}}(x^{\prime})\rangle\\\
&=\sum_{\theta_{p}}-\mathcal{A}[\frac{\partial
h}{\partial\theta_{p}}](x)\langle\zeta(x^{\prime}),\mathcal{A}[\frac{\partial
h}{\partial\theta_{p}}](x^{\prime})\rangle\\\
&=-\langle\zeta(x^{\prime}),\sum_{\theta_{p}}\mathcal{A}[\frac{\partial
h}{\partial\theta_{p}}](x)\mathcal{A}[\frac{\partial
h}{\partial\theta_{p}}](x^{\prime})\rangle,\end{split}$ (C.6)

where $\zeta(x^{\prime})=\frac{\partial L}{\partial
v}(x^{\prime})=v(x^{\prime})-\tilde{v}(x^{\prime})$ is the vector of loss, and
$p$ denotes the index of the learning parameter. We denote the dynamic Neural
Tangent Kernel (DNTK) for the PDE-based neural network as

$\mathcal{N}(x,x^{\prime})=\sum_{\theta_{p}}\mathcal{A}[\frac{\partial
f}{\partial\theta_{p}}](x)\mathcal{A}[\frac{\partial
f}{\partial\theta_{p}}](x^{\prime})=\sum_{l}\langle\mathcal{A}[\frac{\partial
f}{\partial W^{(l)}}](x),\mathcal{A}[\frac{\partial f}{\partial
W^{(l)}}](x^{\prime})\rangle_{W},$ (C.7)

where $\langle\cdot,\cdot\rangle_{W}$ is defined as the summation over each
component index of $W^{(l)}$.

Next, we would give the explicit form of the DNTK for further analysis. Recall
that the output of the MLP in PDE-based neural network has the following form

$\displaystyle\begin{split}h(x,\theta)&=f^{(L+1)}(x)=W^{(L+1)}g^{(L)}(x)\\\
&=W^{(L+1)}\sqrt{\frac{c_{\sigma}}{d_{L}}}\sigma\biggl{(}\cdots
W^{(l+1)}\sqrt{\frac{c_{\sigma}}{d_{l}}}\sigma\bigl{(}W^{(l)}\sqrt{\frac{c_{\sigma}}{d_{l-1}}}\sigma(\cdots)\bigr{)}\cdots\biggr{)},\end{split}$
(C.8)

where we have omitted the explicit dependence of $\theta$ in the formula for
simplicity. And thus

$h(x,\theta)=W^{(L+1)}\sqrt{\frac{c_{\sigma}}{d_{L}}}\sigma\biggl{(}\cdots
f^{(l+1)}\cdots\biggr{)},$ (C.9)

where

$\displaystyle\begin{split}f^{(l+1)}&=W^{(l+1)}g^{(l)},\\\
g^{(l)}&=\sqrt{\frac{c_{\sigma}}{d_{l}}}\sigma\bigl{(}f^{(l)}\bigr{)}.\end{split}$
(C.10)

To give the form of the DNTK, the key is to give the form of $\frac{\partial
h}{\partial W^{(l)}}$. From above forms, we can obtain

$\displaystyle\begin{split}\frac{\partial h}{\partial W^{(l)}}&=\frac{\partial
h}{\partial f^{(l)}}\frac{\partial f^{(l)}}{\partial W^{(l)}}\\\
&=\frac{\partial h}{\partial f^{(l+1)}}\frac{\partial f^{(l+1)}}{\partial
g^{(l)}}\frac{\partial g^{(l)}}{\partial f^{(l)}}\frac{\partial
f^{(l)}}{\partial W^{(l)}}\\\ &=\biggl{(}\frac{\partial h}{\partial
f^{(l+1)}}W^{(l+1)}\sqrt{\frac{c_{\sigma}}{d_{l}}}S^{(l)}\biggr{)}^{T}\bigl{(}g^{(l-1)}\bigr{)}^{T}\\\
&=\sqrt{\frac{c_{\sigma}}{d_{l}}}S^{(l)}\bigl{(}W^{(l+1)}\bigr{)}^{T}\bigl{(}\frac{\partial
h}{\partial f^{(l+1)}}\bigr{)}^{T}\bigl{(}g^{(l-1)}\bigr{)}^{T},\end{split}$
(C.11)

where we have used the denotation

$\bigl{[}S^{(l)}\bigr{]}_{ij}=\bigl{[}\dot{\sigma}(f^{(l)})\bigr{]}_{i}\delta_{ij}.$
(C.12)

By defining $b^{(l)}=\bigl{(}\frac{\partial h}{\partial f^{(l)}}\bigr{)}^{T}$
and we obtain

$\frac{\partial h}{\partial
W^{(l)}}=b^{(l)}\bigl{(}g^{(l-1)}\bigr{)}^{T},l\in[L]$ (C.13)

where $b^{(l)}$ satisfies the induction relation
$b^{(l)}=\sqrt{\frac{c_{\sigma}}{d_{l}}}S^{(l)}\bigl{(}W^{(l+1)}\bigr{)}^{T}b^{(l+1)}$.

With the admissible property of $\mathcal{A}$ in the sense of Definition C.1,
we have the DNTK as

$\displaystyle\begin{split}\mathcal{N}(x,x^{\prime})&=\sum_{l}\bigl{\langle}\mathcal{A}[\frac{\partial
h}{\partial W^{(l)}}](x),\mathcal{A}[\frac{\partial h}{\partial
W^{(l)}}](x^{\prime})\bigr{\rangle}_{W}\\\
&=\sum_{l}\bigl{\langle}\mathcal{A}[b^{(l)}\bigl{(}g^{(l-1)}\bigr{)}^{T}](x),\mathcal{A}[b^{(l)}\bigl{(}g^{(l-1)}\bigr{)}^{T}](x^{\prime})\bigl{\rangle}\\\
&=[\mathcal{A}\mathcal{K}\mathcal{A}](x,x^{\prime}).\end{split}$ (C.14)

where $\mathcal{K}(y,y^{\prime})$ is the dynamic neural tangent kernel of the
MLP with the following form [42]

$\mathcal{K}(y,y^{\prime})=\sum_{l}\bigl{\langle}b^{(l)}(y)\bigl{(}g^{(l-1)}(y)\bigr{)}^{T},b^{(l)}(y^{\prime})\bigl{(}g^{(l-1)}(y^{\prime})\bigr{)}^{T}\bigl{\rangle},$
(C.15)

and $[\mathcal{A}\mathcal{K}\mathcal{A}](x,x^{\prime})$ is a function given by
the kernel $\mathcal{K}(y,y^{\prime})$ operated by $\mathcal{A}$ on its head
and the tail for performing with respect to the former variable $y$ and latter
variable $y^{\prime}$.

Denote the constant neural tangent kernel in [35] as

$\Theta^{(L)}(y,y^{\prime})=\sum_{l=1}^{L+1}\biggl{(}\Sigma^{(l-1)}(y,y^{\prime})\prod_{l^{\prime}=l}^{L+1}\dot{\Sigma}^{(l^{\prime})}(y,y^{\prime})\biggr{)},$
(C.16)

where $\Theta^{(L)}(y,y^{\prime})$ is given by a reduction form

$\displaystyle\Sigma^{(0)}(y,y^{\prime})$ $\displaystyle=y^{T}y^{\prime},$
(C.17) $\displaystyle\Lambda^{(l)}(y,y^{\prime})$
$\displaystyle=\begin{pmatrix}\Sigma^{(l-1)}(y,y)&\Sigma^{(l-1)}(y,y^{\prime})\\\
\Sigma^{(l-1)}(y^{\prime},y)&\Sigma^{(l-1)}(y^{\prime},y^{\prime})\\\
\end{pmatrix},$ (C.18) $\displaystyle\Sigma^{(l)}(y,y^{\prime})$
$\displaystyle=c_{\sigma}\mathbb{E}_{(u,v)\sim\mathcal{N}(0,\Lambda^{(l)})}[\sigma(u)\sigma(v)],$
(C.19) $\displaystyle\dot{\Sigma}^{(l)}(y,y^{\prime})$
$\displaystyle=c_{\sigma}\mathbb{E}_{(u,v)\sim\mathcal{N}(0,\Lambda^{(l)})}[\dot{\sigma}(u)\dot{\sigma}(v)],$
(C.20)

and $\dot{\Sigma}^{(L+1)}(y,y^{\prime})=1$.

The convergence results of both initialization and during training depend on
the Gaussian random initialization. The parameters $W^{(l)},l=0,1,\cdots,L+1$
are initialized by random Gaussian, i.e.,

$W^{(l)}\sim\mathcal{G}(0,I),$ (C.22)

where $\mathcal{G}$ is the Gaussian distribution with the identity covariance
matrix $I$. Such initialization agrees with the so-called “LeCun”[43]
and“Kaiming”[44] initialization with only a constant difference.

The convergence results of the initialization case and during-training case
are summarized in Theorem 4.1 and 4.2 respectively. And the main proof is
given in the Appendix C.4 using the results in [35] and the bounded properties
of the operator $\mathcal{A}$.

###  C.3 Proof of Lemma 1

######  Proof.

The result can be obtained using the positive definiteness of the kernel
$\Theta^{(L)}$ and the invertibility of the operator $\mathcal{A}$. For any
$f\neq 0\in L^{2}(\partial\Omega)$, we have

$\displaystyle\int\int
f(x)[\mathcal{A}\theta^{(L)}\mathcal{A}](x,x^{\prime})f(x^{\prime})dxdx^{\prime}$
$\displaystyle=$ $\displaystyle\int\int
f(x)\int\int\tilde{G}(x,y)\Theta^{(L)}(y,y^{\prime})\tilde{G}(x^{\prime},y^{\prime})dydy^{\prime}f(x^{\prime})dxdx^{\prime}$
$\displaystyle=$ $\displaystyle\int\int dydy^{\prime}\int
f(x)\tilde{G}(x,y)dx\,\Theta^{(L)}(y,y^{\prime})\int\tilde{G}(x^{\prime},y^{\prime})f(x^{\prime})dx^{\prime}$
$\displaystyle=$ $\displaystyle\int\int
l(y)\Theta^{(L)}(y,y^{\prime})l(y^{\prime})dydy^{\prime},$

where $l(y)=\int f(x)\tilde{G}(x,y)dx$ and we have omitted the integral domain
$\partial\Omega$ for simplicity. For double layer potential,
$\mathcal{A}:L^{2}(\partial\Omega)\to L^{2}(\partial\Omega)$ is invertible.
[37] For single layer potential, there is an invertible theorem but with some
constraints, i.e., $\mathcal{A}:L^{2}(\partial\Omega)\to
L_{1}^{2}(\partial\Omega)$ is invertible if $\partial\Omega$ is a $C^{\infty}$
compact domain. [36] Therefore $l\neq 0\in L^{2}(\partial\Omega)$ holds for
double layer potential but holds for single layer potential with constraints.
As the constant kernel $\Theta^{(L)}(y,y^{\prime})$ is positive definite
proven in [12], the new kernel
$[\mathcal{A}\Theta^{(L)}\mathcal{A}](x,x^{\prime})$ is thus positive
definite.

###  C.4 Proof of Theorem 4.1 and Theorem 4.2

In this part, we would give proof of Theorem 4.1 and 4.2. The proof is based
on the result of [35].

######  Proof.

(Proof of Theorem 4.1)

###### Lemma C.1

(Adopted from Theorem 3.1 from [35]) Fix $\epsilon>0$ and $\delta\in(0,1)$.
Suppose the activation nonlinear function $\sigma(\cdot)=max\\{\cdot,0\\}$ is
ReLU, the minimum width of the hidden layer
$\text{min}_{l\in[L]}d_{l}\geq\Omega(\frac{L^{6}}{\epsilon^{4}}\log(L/\delta))$.
Then for the normalized data $y$ and $y^{\prime}$ where $\|y\|\leq 1$ and
$\|y^{\prime}\|\leq 1$, with probability at least $1-\delta$ we have

$|\mathcal{K}_{0}(y,y^{\prime})-\Theta^{(L)}(y,y^{\prime})|\leq(L+1)\epsilon,$

where
$\mathcal{K}_{0}(y,y^{\prime})=\sum_{l}\bigl{\langle}b^{(l)}(y,0)\bigl{(}g^{(l-1)}(y,0)\bigr{)}^{T},b^{(l)}(y^{\prime},0)\bigl{(}g^{(l-1)}(y^{\prime},0)\bigr{)}^{T}\bigl{\rangle}$
is the neural tangent kernel at initialization, i.e., $t=0$.

By Lemma C.1, the result of Theorem 4.1 is directly obtained.

######  Proof.

(Proof of Theorem 4.2)

Denote $W^{(l)}$ the learning parameter and $\tilde{W}^{(l)}$ the perturbation
parameter with the perturbation matrices $\Delta W^{(l)}$, i.e.,
$\tilde{W}^{(l)}=W^{(l)}+\Delta W^{(l)}$, where $\|\Delta W^{(l)}\|_{F}$ is
bounded. Presume that the input data satisfy a distribution $\mathcal{P}$ with
the measure $P_{in}$. Let $\int P_{in}(dx)=1$ and we denote $dx=P_{in}(dx)$
for simplicity.

Let $\tilde{\mathcal{N}}(x,x^{\prime})$ denote the perturbation of
$\mathcal{N}(x,x^{\prime})$ by the perturbation matrices $\Delta W^{(l)}$,
i.e.,

$\displaystyle\mathcal{N}(x,x^{\prime})$
$\displaystyle=\sum_{l}\langle\mathcal{A}[\frac{\partial f}{\partial
W^{(l)}}](x),\mathcal{A}[\frac{\partial f}{\partial
W^{(l)}}](x^{\prime})\rangle_{W}$ (C.23)
$\displaystyle\tilde{\mathcal{N}}(x,x^{\prime})$
$\displaystyle=\sum_{l}\langle\mathcal{A}[\frac{\partial
f}{\partial\tilde{W}^{(l)}}](x),\mathcal{A}[\frac{\partial
f}{\partial\tilde{W}^{(l)}}](x^{\prime})\rangle_{W}$ (C.24)

###### Lemma C.2

(The reduction of the kernel of BINet)

$|\mathcal{N}_{t}(x,x^{\prime})-\mathcal{N}_{0}(x,x^{\prime})|\leq
A^{2}\|\mathcal{K}_{t}-\mathcal{K}_{0}\|_{\infty},$ (C.26)

where
$\|\mathcal{K}\|_{\infty}=\sup_{y,y^{\prime}}|\mathcal{K}(y,y^{\prime})|$.

######  Proof.

$\displaystyle\begin{split}|\tilde{\mathcal{N}}(x,x^{\prime})-\mathcal{N}(x,x^{\prime})|&=|\mathcal{A}\sum_{l}(\langle\frac{\partial
f}{\partial\tilde{W}^{(l)}},\frac{\partial
f}{\partial\tilde{W}^{(l)}}\rangle-\langle\frac{\partial f}{\partial
W^{(l)}},\frac{\partial f}{\partial W^{(l)}}\rangle)\mathcal{A}|\\\
&\leq\|\mathcal{A}\|_{\infty}\|\mathcal{A}\|_{\infty}\|\widetilde{\mathcal{K}}-\mathcal{K}\|_{\infty}\\\
&\leq A^{2}\|\widetilde{\mathcal{K}}-\mathcal{K}\|_{\infty}.\end{split}$
(C.27)

We have shown that the perturbation of the new kernel can be controlled by the
perturbation of the kernel of the neural network. Thus, the perturbation
during training of the new kernel of BINet can be controlled by that of the
kernel of the neural network. We complete the proof.

###### Lemma C.3

(Adopted from [35]) We have for all $t\geq 0$,

$|\mathcal{K}_{t}(y,y^{\prime})-\mathcal{K}_{0}(y,y^{\prime})|\leq\omega,$

if the following holds

$\|W^{(l)}(t)-W^{(l)}(0)\|_{F}=\mathcal{O}(A\frac{\sqrt{n}}{\lambda_{0}})\leq\omega\sqrt{m}$
(C.28)

for any $l$ and $t$. Here $\mathcal{K}_{t}(y,y^{\prime})$ is the kernel along
training step $t$ and $\mathcal{K}_{0}(y,y^{\prime})$ is the kernel.

Lemma C.3 is a trivial property of the neural network from [35]. But the
existence of its condition in BINet is nontrivial because of the operator
$\mathcal{A}$. Luckily in our BINet, the condition can be proven to hold since
our integral operator is bounded. Specifically, from [35] if the following
lemma holds for our BINet, the condition of Lemma C.3 exists. Thus the only
thing left is to verify if the following lemma still holds for our BINet.

###### Lemma C.4

(Adopted from Lemma F.7 in [35]) Fix
$\omega\leq\text{poly}(1/L,1/n,1/\log(1/\delta),\lambda_{0})$ and
$\delta\in(0,1)$. Suppose that
$\text{min}_{l}{d_{l}}\geq\text{poly}(1/\omega)$. Fixed $l^{\prime}\in[L+1]$,
we have with probability at least $1-\delta$ over random initialization, for
all $t\geq 0$

$\|W^{(l)}(t)-W^{(l)}(0)\|_{F}=\mathcal{O}(A\frac{\sqrt{n}}{\lambda_{0}})\leq\omega\sqrt{m},$
(C.29)

if following inequalities hold for $\forall
l\in[L+1]\backslash\\{l^{\prime}\\}$:

$\displaystyle\begin{split}&\|{\bf\tilde{v}}_{nn}(t)-{\bf
v}\|_{2}\leq\exp(-\frac{1}{2}\kappa^{2}\lambda_{0}t)\|{\bf\tilde{v}}_{nn}(0)-{\bf
v}\|_{2},\\\ &\|W^{(l)}(t)-W^{(l)}(0)\|_{F}\leq\omega\sqrt{m}.\end{split}$
(C.30)

######  Proof.

$\displaystyle\begin{split}\|W^{(l)}(t)-W^{(l)}(0)\|_{F}&=\|\int_{0}^{t}\frac{dW^{(l)}(\tau)}{d\tau}d\tau\|_{F}\\\
&=\|\int_{0}^{t}\frac{\partial L(\tau)}{\partial W^{(l)}(\tau)}d\tau\|_{F}\\\
&=\|\int_{0}^{t}\frac{1}{n}\sum_{i=1}^{n}(\tilde{v}_{i}(\tau)-v_{i})\mathcal{A}[\frac{\partial
h(\theta(\tau),y_{i})}{\partial W^{(l)}}](x_{i})d\tau\|_{F}\\\ &\leq
C\frac{\|A\|_{\infty}}{n}\max_{0\leq\tau\leq t}\sum_{i=1}^{n}\|\frac{\partial
h_{nn}(\theta(\tau),y_{i})}{\partial
W^{(l)}}\|_{F}\int_{0}^{t}\|{\bf\tilde{v}}_{nn}(\tau)-{\bf
v}\|_{2}d\tau.\end{split}$ (C.31)

Thus by Lemma F.7 in [35], we can prove for our BINet, for any $t$ and $l$,
the following also holds

$\|W^{(l)}(t)-W^{(l)}(0)\|_{F}=\mathcal{O}(\sqrt{\frac{n}{\lambda_{0}}})\leq\omega\sqrt{m}.$
(C.32)

Till now we have verified the lazy properties are satisfied during training,
by Lemma C.3 and C.2, we complete the proof of Theorem 4.2.

##  Appendix D Repeated Experiments with Random Initialization

With the limit of the number of pages, only some experimental results have
been shown in the paper. To better exhibit the accuracy of the experiments, we
repeated each experiment of all methods 5 times and show the all relative
$L^{2}$ error. For each experiment, in addition to the initialization of
network parameters, other conditions like the optimizer and learning rate are
consistent.

Lapalace Equation with Smooth Boundary Condition. In the first experiment, the
PINN, Deep Ritz and BINet methods are used respectively to solve the following
equaiton

$\displaystyle-\Delta u(x,y)=0,\ (x,y)\in\Omega,$ (D.1) $\displaystyle
u(x,y)=e^{ax}\sin(ay),\ (x,y)\in\partial\Omega.$

where $a$ is taken as $4$ and $8$ for the low and high scales. The experiments
are repeated five times at each setup, i.e., we have done totally $5\times
3\times 2$ times training. Figure 8 shows all the results of relative $L^{2}$
error. It is easily to conclude that our BINet is the best and outperforms
related methods by a significant margin.

(a) a=4.

(b) a=8.

Figure 8: The relative $L^{2}$ error equation of (D.1) by using PINN, Deep
Ritz, and BINet five repeated experiments. The left figure shows the results
of $a=4$, and the right one shows the results of $a=8$.

Laplace Equation with Non-smooth Boundary Condition. In this experiment, we
consider a Laplace equation problem

$\displaystyle-\Delta u(x,y)=0,\ (x,y)\in\Omega,$ (D.2) $\displaystyle
u(x,y)=2-|x|-|y|,\ (x,y)\in\partial\Omega.$

with the non-smooth boundary condition. The performance of the PINN, Deep
Ritz, and BINet has been shown in the main part of the paper. In this
appendix, the numerical experiments using each method are run five times and
the results of all relative $L^{2}$ error are shown in Figure 9. It is easy to
find that the results using BINet have much higher accuracy than others with
more than 100x decrease of the $L^{2}$ error at most. Moreover, the numerical
stability of BINet is also much better.

Figure 9: The relative $L^{2}$ error equation of (D.2) by using PINN, Deep
Ritz, and BINet of 5 repeated experiments.

Helmholtz Equation with Different Wave numbers. In the main part of the paper,
we have considered a Helmholtz equation with different wave numbers. The
Helmholtz equation problem has following form

$\displaystyle-\Delta u(x,y)-k^{2}u(x,y)=0,\ (x,y)\in\Omega,$ (D.3)
$\displaystyle u(x,y)=e^{i(k_{1}x+k_{2}y)},\ (x,y)\in\partial\Omega,$

where $(k_{1},k_{2})=(k\cos\frac{\pi}{7},k\sin\frac{\pi}{7})$, and the
boundary has the parametric representation
$\partial\Omega=\\{(\frac{9}{20}\cos(t)-\frac{\cos(5t)}{9}\cos(t),\frac{9}{20}\sin(t)-\frac{\cos(5t)}{9}\sin(t))|t\in[0,2\pi]\\}$.
In the experiment, we assume $k=1$ and $4$. We also repeated the experiments 5
times for both PINN and BINet methods. The relative $L^{2}$ error of each
experiment has shown in Figure 10.

(a) k=1.

(b) k=4.

Figure 10: The relative $L^{2}$ error equation of (D.3) by using PINN and
BINet of 5 repeated experiments. The left figure shows the results of $k=1$,
and the other shows the result of $k=4$.

In Figure 10 we can see the results of repeated experiments are similar. This
also confirms the stability of our method. We know the difficulty of the
Helmholtz equation will increase with the increase of the wave number. PINN
method failed when the wave number $k=4$, but the relative $L^{2}$ error of
the BINet method has no obvious difference between $k=4$ and $k=1$.

The Operator from Equation Parameters to Solutions. In the main part of the
paper, we have verified the ability of BINet to learn operators and solve PDEs
on the unbounded domain. We also tested the generalization ability of BINet.
We consider the Helmholtz equations with various wave numbers.

$\displaystyle-\Delta u(x,y)-k^{2}u(x,y)=0,\ (x,y)\in\Omega,$ (D.4)
$\displaystyle u(x,y)=H_{0}^{1}(k\sqrt{x^{2}+y^{2}}),\
(x,y)\in\partial\Omega,$

In the training phase, we set $k\in[2,3.5]\cup[4.5,6]$. Figure 11 shows the
relative $L^{2}$ error of different wave numbers. Five lines of different
colors represent the results of five repeated experiments. We can see the five
results are similar, too. When BINet learns an operator to solve a parametric
PDE, this method still has good numerical stability.

Figure 11: The relative $L^{2}$ error equation of (D.4) with different wave
numbers by using and BINet of five repeated experiments.

The Operator from Boundary Geometry to Solutions. In the final experiment of
the main paper, we used the BINet method to learn a operator from boundary
geometry to solutions. We consider a Laplace equation with parametric
boundaries. The problem is

$\displaystyle-\Delta u(x,y)$ $\displaystyle=0,\ $
$\displaystyle(x,y)\in\Omega_{\beta},$ (D.5) $\displaystyle u(x,y)$
$\displaystyle=g(x,y;\beta),\ $ $\displaystyle(x,y)\in\partial\Omega_{\beta},$

where the boundary condition
$g(x,y;\beta)=(x-x_{bc})(y-y_{bc})+(x-x_{bc})+(y-y_{bc})+1$, and
$(x_{bc},y_{bc})$ is the barycenter of the $\Omega_{\beta}$. For details of
the parametric boundary, please refer to the main paper. We repeated the
experiment five times by using the BINet method. After each experiment, we
selected 100 triangles and calculate the relative $L^{2}$ error for each
triangle. Figure 12 shows the average relative $L^{2}$ error of the 100
triangles of each experiment.

Figure 12: The relative $L^{2}$ error equation of (D.5) by using BINet of five
repeated experiments. The relative $L^{2}$ error of each experiment is the
average over 100 triangles.

##  Appendix E More Experiments

Comparison of single layer potential and double layer potential. In this
experiment, we use the single layer potential and double layer potential to
construct BINet, respectively. We also consider this Helmholtz equation with
wave number $k=4$,

$\displaystyle-\Delta u(x,y)-16u(x,y)=0,\ (x,y)\in\Omega,$ (E.1)
$\displaystyle u(x,y)=e^{i(k_{1}x+k_{2}y)},\ (x,y)\in\partial\Omega,$

where $\partial\Omega$ has the same parametric representation as the previous
boundary of the problem (D.3), and in this example, we also assume
$(k_{1},k_{2})=(4\cos\frac{\pi}{7},4\sin\frac{\pi}{7})$. We can know the exact
solution of this problem is

$u(x,y)=e^{i(k_{1}x+k_{2}y)},\ (x,y)\in\Omega.$

We select 600 sample points on the boundary to construct the loss function,
and we use the Adam optimizer with learning rate 0.0001. After 40000 training
epochs, we random 1000 points in $\Omega$. By comparing with the exact
solution, the relative $L^{2}$ error is approximate 0.0086 for single layer
potential and 0.0016 for double layer potential, and we have the following
results of the loss function and error.

Figure 13: The results of the equation (E.1). The first figure is the loss
function of BINet with single and double layer potentials respectively. The
two figures on the right are the absolute error between the numerical
solutions and the exact solution with single and double layer potentials
respectively.

From Figure 13, we can see the solution of the double layer potential is
better than the solution of the single layer potential. In fact, because of
the jump of the double layer potential at the boundary, the condition number
of the double layer potential is better than the single, and the BINet with
double layer potential converges faster.

Laplace equation on unbounded domain. We assume the boundary of $\Omega$ is
same as the previous experiment. We will use BINet to solve the following
equation

$\displaystyle-\Delta u(x,y)=0,\ (x,y)\in\Omega^{c},$ (E.2) $\displaystyle
u(x,y)=\frac{x}{x^{2}+y^{2}},\ (x,y)\in\partial\Omega,$

The exact solution of this example is

$u(x,y)=\frac{x}{x^{2}+y^{2}},\ (x,y)\in\Omega^{c}.$

In this experiment, the double layer potential is used to construct the loss
function. Adam optimizer is used in the optimization with the learning rate of
0.001. We selected 600 sample points on the boundary. After 20000 training
epochs, the relative $L^{2}$ error is about 0.0030, and the error map is shown
in Figure 14. This experiment verifies the feasibility of BINet to solve the
exterior Laplace equations on the unbounded domain.

Figure 14: The result of the equation (E.2). The figure shows the error
between the solution of BINet and the exact solution on the domain
$\\{(x,y)|x^{2}+y^{2}<1.6^{2}\\}\cap\Omega^{c}$.

Helmholtz equations with high wave numbers. In this experiment, we consider an
interior problem of Helmholtz equation on a bounded domain with wave number
$k=10$. Let us consider this Helmholtz equation in an interior domain

$\displaystyle-\Delta u(x,y)-100u(x,y)=0,\ (x,y)\in\Omega_{1},$ (E.3)
$\displaystyle u(x,y)=e^{i(k_{1}x+k_{2}y)},\ (x,y)\in\partial\Omega_{1},$

where $(k_{1},k_{2})=(10\cos\frac{\pi}{5},10\sin\frac{\pi}{5})$, and in an
exterior domain with wave number $k=8$

$\displaystyle-\Delta u(x,y)-64u(x,y)=0,\ (x,y)\in\Omega_{2}^{c},$ (E.4)
$\displaystyle u(x,y)=H_{0}^{1}(8\sqrt{(x-0.5)^{2}+y^{2}}),\
(x,y)\in\partial\Omega_{2},$

where $H_{0}^{1}$ is the first kind Hankel function. The boundaries
$\partial\Omega_{1}$ and $\partial\Omega_{2}$ are shown in Figure 15. By
BINet, the relative $L^{2}$ errors of solutions of the equations (E.3) and
(E.4) are $1.5\%$ and $2.8\%$, respectively.

Figure 15: The results of equations (E.3) and (E.4). The left figure is the
absolute error of the equation (E.3) and the other is the absolute error of
the equation (E.4).