Ensemble learning for Physics Informed Neural Networks: a Gradient
  Boosting approach
While the popularity of physics-informed neural networks (PINNs) is steadily
rising, to this date, PINNs have not been successful in simulating multi-scale
and singular perturbation problems. In this work, we present a new training
paradigm referred to as "gradient boosting" (GB), which significantly enhances
the performance of physics informed neural networks (PINNs). Rather than
learning the solution of a given PDE using a single neural network directly,
our algorithm employs a sequence of neural networks to achieve a superior
outcome. This approach allows us to solve problems presenting great challenges
for traditional PINNs. Our numerical experiments demonstrate the effectiveness
of our algorithm through various benchmarks, including comparisons with finite
element methods and PINNs. Furthermore, this work also unlocks the door to
employing ensemble learning techniques in PINNs, providing opportunities for
further improvement in solving PDEs.
Zhiwei Fang, Sifan Wang, Paris Perdikaris
# Ensemble learning for Physics Informed Neural Networks: a Gradient Boosting approach

###### Abstract

While the popularity of physics-informed neural networks (PINNs) is steadily rising, to this date, PINNs have not been successful in simulating multi-scale and singular perturbation problems. In this work, we present a new training paradigm referred to as "gradient boosting" (GB), which significantly enhances the performance of physics informed neural networks (PINNs). Rather than learning the solution of a given PDE using a single neural network directly, our algorithm employs a sequence of neural networks to achieve a superior outcome. This approach allows us to solve problems presenting great challenges for traditional PINNs. Our numerical experiments demonstrate the effectiveness of our algorithm through various benchmarks, including comparisons with finite element methods and PINNs. Furthermore, this work also unclocks the door to employing ensemble learning techniques in PINNs, providing opportunities for further improvement in solving PDEs.

PDE Physics-informed neural networks Gradient boosting Ensemble learning

## 1 Introduction

Physics informed neural networks have recently emerged as an alternative to traditional numerical solvers for simulations in fluids mechanics [1, 2], bio-engineering [3, 4], meta-material design [5, 6], and other areas in science and engineering [7, 8]. However, PINNs using fully connected, or some variants architectures such as Fourier feature networks [9], fail to accomplish stable training and produce accurate predictions at whites, especially when the underlying PDE solutions contain high-frequencies or multi-scale features [10, 11]. To mitigate this pathology, Krishnapriyan _et. al._[12] proposed a sequence-to-sequence learning method for time-dependent problems, which divide the time domain into sub-intervals and solve the problem progressively on each them. This method avoids the pollution of the underlying solution due to the temporal error accumulation. Wang _et. al._[13] elaborated the reason that the PINNs fail to train from a neural tangent kernel perspective, and proposed an adaptive training strategy to improve the PINNs' performance. An empirical learning-rate annealing scheme has been proposed in Wang _et. al._[14], which utilizes the back-propagated gradient statistics during training to adaptively assign importance weights to different terms in a PINNs loss function, with the goal of balancing the magnitudes of the gradients in backward propagation. Although all of these works were demonstrated to produce significant and consistent improvements in the stability and accuracy of PINNs, the fundamental reasons behind the practical difficulties of training fully-connected PINNs still remain unclear [10].

Besides PINNs, many other machine learning tasks suffer from the same issues, and some of these issues have been resolved by gradient boosting method. The idea of gradient boosting is blending several weak learners into a fortified one that gives better predictive performance than could be obtained from any of the constituent learners alone [15]. For example, Zhang and Haghani [16] proposes a gradient-boosting tree-based travel time prediction method, driven by the successful application of random forest in traffic parameter prediction, to uncover hidden patterns in travel time data to enhance the accuracy and interpretability of the model. Callens _et. al._[17] used gradient boosting trees to improve wave forecast at a specific location whose RMSE values in average \(8\%\) to \(11\%\) lower for the correction of significant wave height and peak wave period. Recently, many researchers have contributed to gradient boosting
method and further improved its performance. Friedman _et. al._[18] shows that both the approximation accuracy and execution speed of gradient boosting can be substantially improved by incorporating randomization into the procedure, and this randomized approach also increases robustness against overcapacity of the base learner. Ke [19] found that the efficiency and scalability of Gradient Boosting Decision Tree (GBDT) are unsatisfactory when the feature dimension is high and data size is large and a greedy algorithm has been used to effectively reduce the number of features without hurting the accuracy of split point determination by much and thus solve the issue.

Inspired by the above-mentioned literature review, we arrive at our proposed method in this paper. In this work, we present a gradient boosting physics informed neural networks (GB PINNs), which adopts a gradient boosting idea to approximate the underlying solution by a sequence of neural networks and train the PINNs progressively. Specifically, our main contributions can be summarized into the following points:

1. A simple implementation of the gradient boosting method, which can easily be integrated into existing PINNs code with minimal modifications.
2. The assembly of several weak PDE predictors to form a strong predictor, resulting in increased flexibility for solving intractable problems.
3. Low sensitivity to the choice of neural networks and their arrangement, resulting in fewer efforts required for fine-tuning hyperparameters.
4. The flexibility to combine with other techniques, such as Fourier features, making it a versatile approach for PINNs.

We introduce some preliminaries for key ingredients of our algorithm in section 2. Then we present our algorithm with motives in section 3. Numerical experiments are shown in section 4 to verify our algorithm. We discuss our algorithm and conclude the paper in section 5.

## 2 Preliminaries

In this section, we will provide a brief overview of the related topics that are relevant to the proposed algorithm in this paper. For a more in-depth understanding of these topics, we encourage readers to refer to the original papers cited below.

### Physics informed neural networks

PINNs are a method for inferring a continuous latent function \(u(\mathbf{x})\) that serves as the solution to a nonlinear PDE of the form:

\[\mathcal{N}[u](\mathbf{x}) =0,\quad\text{ in }\Omega, \tag{1}\] \[\mathcal{B}[u](\mathbf{x}) =0,\quad\text{ on }\partial\Omega, \tag{2}\]

where \(\Omega\) is an open, bounded set in \(\mathbb{R}^{d}\) with a piecewise smooth boundary \(\partial\Omega\), \(\mathbf{x}\in\mathbb{R}^{d}\), and \(\mathcal{N}\) and \(\mathcal{B}\) are nonlinear differential and boundary condition operators, respectively.

The solution to the PDE is approximated by a deep neural network, \(u_{\theta}\), which is parameterized by \(\theta\). The loss function for the network is defined as:

\[L(u;\theta)=\frac{\omega_{e}}{N_{p}}\sum_{i=1}^{N_{p}}|\mathcal{N}[u_{\theta} ](\mathbf{x}_{i}^{p})|^{2}+\frac{\omega_{b}}{N_{b}}\sum_{i=1}^{N_{b}}|\mathcal{B} [u_{\theta}](\mathbf{x}_{i}^{b})|^{2}, \tag{3}\]

where \(\{\mathbf{x}_{i}^{p}\}_{i=1}^{N_{p}}\) and \(\{\mathbf{x}_{i}^{b}\}_{i=1}^{N_{b}}\) are the sets of points for the PDE residual and boundary residual, respectively, and \(\omega_{e}\) and \(\omega_{b}\) are the weights for the PDE residual loss and boundary loss, respectively. The neural network \(u_{\theta}\) takes the coordinate \(\mathbf{x}\) as input and outputs the corresponding solution value at that location. The partial derivatives of the \(u_{\theta}\) with respect to the coordinates at \(\mathcal{N}\) in (3) can be readily computed to machine precision using reverse mode differentiation [20].

The loss function \(L(u;\theta)\) is typically minimized using a stochastic gradient descent algorithm, such as Adam, with a batch of interior and boundary points generated to feed the loss function. The goal of this process is to find a set of neural network parameters \(\theta\) that minimize the loss function as much as possible.

It is worth noting that the abstract PDE problem in (1)-(2) can easily be extended to time-dependent cases by considering one component of \(\mathbf{x}\) as a temporal variable. In this case, one or more initial conditions should be included in the PDE system and additional initial condition constraints should be added to the loss function (3).
### Gradient boosting machines

Gradient Boosting (GB) is a powerful machine learning technique that is commonly used in regression and classification tasks. It is an additive ensemble of weak prediction models, similar to AdaBoost, but with a key difference - unlike other ensemble algorithms, GB does not have trainable weights, and the sub-models are trained sequentially instead of in parallel. For the sake of simplicity, in the rest of the paper, we will use \(f(x;\theta)\) to denote a general neural network \(f\) with input \(x\) and parameterized by \(\theta\).

Given a neural network \(f(x;\theta)\) and a training dataset, the loss function \(L(f;\theta)\) is defined as the sum of the individual losses for each sample, as follows:

\[L(f;\theta)=\sum_{i=1}^{N}L(y_{i},f(x_{i};\theta)),\]

where \(N\) is the total number of samples in the dataset, \(y_{i}\) is the true label for sample \(x_{i}\), and \(f(x_{i};\theta)\) is the predicted label for sample \(x_{i}\).

To minimize this loss function, a common approach is to use the stochastic gradient descent algorithm. This algorithm updates the network's parameters, \(\theta\), iteratively using the following update rule:

\[\theta\leftarrow\theta-\gamma\frac{\partial}{\partial\theta}L(f;\theta), \tag{4}\]

where \(\gamma\) is a user-specified learning rate that controls the step size of the updates.

The goal of the gradient boosting (GB) method is to minimize the loss function \(L(f;\theta)\) with respect to the neural network function \(f\). GB method assumes that the surrogate model can be represented in the following iterative form:

\[f_{m}(\mathbf{x};\Theta_{m})=f_{m-1}(\mathbf{x};\Theta_{m-1})+\rho_{m}h_{m}(\mathbf{x}; \theta_{m}),\quad\text{ for }m=1,2,3,\cdots, \tag{5}\]

where \(f_{0}(\mathbf{x};\theta_{0})\) is a pre-selected baseline neural network, \(\rho_{m}\) is the learning rate, \(f_{m}(\mathbf{x};\Theta_{m-1})\) is parameterized by \(\Theta_{m}=\bigcup_{i=0}^{m}\theta_{m}\), and \(h_{m}(\mathbf{x};\theta_{m})\) is a neural network designed to enhance the accuracy of the predictor \(f_{m-1}(\mathbf{x};\Theta_{m-1})\). The gradient descent algorithm is used to choose \(h_{m}(\mathbf{x};\theta_{m})\), which is defined as:

\[h_{m}(\mathbf{x};\theta_{m})=-\frac{\partial}{\partial f_{m-1}(\mathbf{x};\Theta_{m- 1})}L(f_{m-1};\Theta_{m-1}). \tag{6}\]

Therefore, the model update rule is defined as:

\[f_{m}(\mathbf{x};\Theta_{m})=f_{m-1}(\mathbf{x};\Theta_{m-1})-\rho_{m}\frac{\partial }{\partial f_{m-1}(\mathbf{x};\Theta_{m-1})}L(f_{m-1};\Theta_{m-1}). \tag{7}\]

In this fashion, the corresponding loss at the \(m\)-th step reads

\[L(f_{m};\Theta_{m})=L(f_{m-1}+\rho_{m}h_{m};\Theta_{m}). \tag{8}\]

The technique outlined in this construction is commonly referred to as a GB method. The update function, \(h_{m}(\mathbf{x};\theta_{m})\) in equation (6), is similar in nature to the gradient vector in equation (4), however, GB operates by taking the gradient with respect to the function, rather than the parameter vector as traditional gradient descent does. This distinction is the reason why we refer to GB as a method that descends the gradient in function space. For further information on gradient boosting methods, please refer to the reference [21]. Furthermore, it is worth noting that in the context of PINNs, this method has been adapted to a simpler form that is easily implementable.

## 3 Gradient boosting physics informed neural networks

Despite a series of promising results in the literature [22; 4; 11; 1; 2], the original formulation of PINNs proposed by Raissi _et. al._[23] has been found to struggle in constructing an accurate approximation of the exact latent solution. The reasons for this remain poorly understood. However, some observations in the literature can be used to infer potential solutions to this issue. One such observation is that the prediction error in PINNs is often of high frequency, with small and intricate structures, as seen in figures 4(b) and 6(a) and (b) of [13]. As demonstrated in [9], high-frequency functions can be learned relatively easily using Fourier features. Based on these findings, it is natural to consider using a multi-layer perceptrons (MLPs) as a baseline structure in PINNs, followed by a Fourier feature network, to further minimize the error. This idea led to the development of GB PINNs.

The proposed method, referred to as GB PINNs, utilizes a sequence of neural networks in an iterative update procedure to gradually minimize the loss function. As shown in equation (6), the update model \(h_{m}(\mathbf{x};\theta_{m})\) is defined by the
gradient of the loss with respect to the previous output \(f_{m-1}\). However, in the context of PINNs, the PDE residual loss in (3) typically includes gradients of the outputs with respect to the inputs. This necessitates the computation of twisted gradients, which is a unique characteristic of this approach. For example

\[\frac{\partial}{\partial f(\mathbf{x};\theta)}\left[\left(\frac{\partial f(\mathbf{x}; \theta)}{\partial\mathbf{x}}\right)^{2}\right],\]

which is definitely not elementary and should be avoided. Despite the mathematical validity of the gradient

\[\frac{\partial}{\partial f_{m-1}(\mathbf{x};\Theta_{m-1})}L(f_{m-1};\Theta_{m-1}),\]

it can be challenging to compute it using AD due to the fact that \(L(f_{m-1};\theta)\) is typically a leaf node in the computational graph.

Fortunately, we can still utilize the formulation in equation (8) to establish an appropriate GB algorithm for PINNs. After the training of the \((m-1)\)-th step is completed, we add an additional pre-selected neural network \(\rho_{m}h_{m}(\mathbf{x};\theta_{m})\) to the previous predictor \(f_{m-1}(\mathbf{x};\Theta_{m-1})\) and subsequently apply gradient descent with respect to the parameters \(\theta_{m}\). This iterative procedure allows us to gradually minimize the loss and improve the accuracy of the predicted solution.

It is important to note that the neural networks utilized in the proposed GB PINNs algorithm do not need to possess a consistent structure. In fact, they can be composed of a variety of surrogate models, as long as they have the same input and output dimensions. Examples of such models include MLPs, Fourier feature networks, radial basis neural networks, and even finite element representations. This flexibility allows for a more versatile approach to minimizing the loss and improving the accuracy of the approximation to the exact latent solution.

As highlighted in [13], the neural tangent kernel of wide and deep neural networks remains relatively stable during the training process, which can impede the ability to learn solutions with sharp gradients. However, it has been proposed that by decreasing the size of the network, such issues can be addressed. This gives rise to the idea of utilizing a small network for the initial approximation of the solution, and then progressively refining it using larger networks. 1

Footnote 1: It is worth noting that the proposed organization of neural networks is merely a rough idea and may not always yield the best performance. The choice and arrangement of neural networks may vary depending on the specific problem at hand. The numerical experiments will provide further insights into the optimal configuration of neural networks for a given problem.

In GB training, the term \(\rho_{m}\) in equation (5) serves as a learning rate in the function space gradient descent, and it also adjusts the magnitude of \(h_{m}(\mathbf{x};\theta_{m})\). In PINNs scenarios, once the training of \(f_{m-1}(\mathbf{x};\Theta_{m-1})\) is complete, we already have a decent predictor, which implies that the relative error between the current predictor and the ground truth (for example, relative \(l^{2}\) error defined below) is assumed to be small. The additive model \(h_{m}(\mathbf{x};\theta_{m})\) is then used to handle this error, thus it is reasonable to assume that the magnitude of \(h_{m}(\mathbf{x};\theta_{m})\) decreases over iteration step \(m\). As a result, it is also reasonable to assume that the \(\rho_{m}\) gradually decreases over \(m\). In the experiments below, we will assume that the \(\rho_{m}\) exponentially decays, which is similar to traditional gradient descent methods.

The proposed algorithm can be summarized as follows:

```
0: A baseline neural network \(f_{0}(\mathbf{x};\theta_{0})\) and an ordered neural network set \(\{h_{i}(\mathbf{x};\theta_{i})\}_{i=1}^{m}\) that contains models going to be trained in sequence; A set of learning rate \(\{\rho_{i}\}_{i=0}^{m}\) that correspond to \(\{f_{0}(\mathbf{x};\theta_{0})\}\cup\{h_{i}(\mathbf{x};\theta_{i})\}_{i=1}^{m}\). Usually, \(\rho_{0}=1\) and \(\rho_{i}\) is decreasing in \(i\); Set \(f_{i}(\mathbf{x};\Theta_{i})=f_{i-1}(\mathbf{x};\Theta_{i-1})+\rho_{i}h_{i}(\mathbf{x}; \theta_{i}),\quad\text{ for }i=1,2,3,\cdots,m\). Given PDEs problem (1)-(2), establish the corresponding loss (3).
0:
1: Train \(f_{0}(\mathbf{x};\theta_{0})=\rho_{0}f_{0}(\mathbf{x};\theta_{0})\) to minimize loss (3).
2:for\(i=1\) to m do
3: In \(f_{i}(\mathbf{x};\Theta_{i})=f_{i-1}(\mathbf{x};\Theta_{i-1})+\rho_{i}h_{i}(\mathbf{x}; \theta_{i})\), set trainable parameters as \(\theta_{i}\). Train \(f_{i}(\mathbf{x};\Theta_{i})\) to minimize loss (3).
4:endfor
5:return\(f_{m}(\mathbf{x};\Theta_{m})\) as a predictor of the solution of (1)-(2) for any point in \(\overline{\Omega}\).
```

**Algorithm 1** Gradient boosting physics informed neural network.

The proposed algorithm can be summarized as follows:

```
0: A baseline neural network \(f_{0}(\mathbf{x};\theta_{0})\) and an ordered neural network set \(\{h_{i}(\mathbf{x};\theta_{i})\}_{i=1}^{m}\) that contains models going to be trained in sequence; A set of learning rate \(\{\rho_{i}\}_{i=0}^{m}\) that correspond to \(\{f_{0}(\mathbf{x};\theta_{0})\}\cup\{h_{i}(\mathbf{x};\theta_{i})\}_{i=1}^{m}\). Usually, \(\rho_{0}=1\) and \(\rho_{i}\) is decreasing in \(i\); Set \(f_{i}(\mathbf{x};\Theta_{i})=f_{i-1}(\mathbf{x};\Theta_{i-1})+\rho_{i}h_{i}(\mathbf{x}; \theta_{i}),\quad\text{ for }i=1,2,3,\cdots,m\). Given PDEs problem (1)-(2), establish the corresponding loss (3).
0:
1:for\(i=1\) to m do
2: In \(f_{i}(\mathbf{x};\Theta_{i})=f_{i-1}(\mathbf{x};\Theta_{i-1})+\rho_{i}h_{i}(\mathbf{x}; \theta_{i})\), set trainable parameters as \(\theta_{i}\). Train \(f_{i}(\mathbf{x};\Theta_{i})\) to minimize loss (3).
3:endfor
5:return\(f_{m}(\mathbf{x};\Theta_{m})\) as a predictor of the solution of (1)-(2) for any point in \(\overline{\Omega}\).
```

**Algorithm 2** Gradient boosting physics informed neural network.

The proposed algorithm, described in Algorithm 1, utilizes a sequence of neural networks and an iterative update procedure to minimize the loss gradually. At each iteration step \(i\), the forward prediction relies on the union parameter
set \(\Theta_{i}\), while the backward gradient propagation is only performed on \(\theta_{i}\). This results in a mild increase in computational cost during the training of GB iteration. The simplicity of this algorithm allows practitioners to easily transfer their PINNs codes to GB PINNs'. In the following section, we will demonstrate that this simple technique can enable PINNs to solve many problems that were previously intractable using the original formulation of Raissi _et. al._[23].

Additionally, the proposed GB PINNs algorithm also introduces another dimension of flexibility in terms of network architecture design, namely the combination of different neural networks. This opens up new opportunities for fine-tuning the architecture to minimize PDE residual losses and improve overall predictive accuracy. As shown in the following section, the performance of GB PINNs is relatively insensitive to the specific choice and arrangement of networks, as long as their capacity is sufficient.

## 4 Numerical Experiments

In this section, we will demonstrate the effectiveness of the proposed GB PINNs algorithm through a comprehensive set of numerical experiments. To simplify the notation, we use a tuple of numbers to denote the neural network architecture, where the tuple represents the depth and width of the layers. For example, a neural network with a two-dimensional input and a one-dimensional output, as well as two hidden layers with width \(100\) is represented as \((2,100,100,1)\). Our default experimental setup is summarized in Table 4, and will be used in all experiments unless otherwise specified.

To quantify the model's accuracy, we use the relative \(l^{2}\) error over a set of points \(\{x_{i}\}_{i=1}^{N}\):

\[\text{Error}=\frac{\sum_{i=1}^{N}|u_{pred}(x_{i})-u_{true}(x_{i})|^{2}}{\sum _{i=1}^{N}|u_{true}(x_{i})|^{2}}.\]

### 1D singular perturbation

In this first example, we utilize GB PINNs to solve the following 1D singular perturbation problem.

\[-\varepsilon^{2}u^{\prime\prime}(x)+u(x) =1,\qquad\text{ for }x\in(0,1),\] \[u(0) =u(1) =0.\]

The perturbation parameter, \(0<\varepsilon\ll 1\), is set to \(10^{-4}\) in this case. The exact solution to this problem is given by

\[u(x)=1-\frac{e^{-x/\varepsilon}+e^{(x-1)/\varepsilon}}{1+e^{-1/\varepsilon}}.\]

Despite the boundedness of the solution, it develops boundary layers at \(x=0\) and \(x=1\) for small values of \(\varepsilon\), a scenario in which traditional PINNs have been known to perform poorly.

To evaluate the performance of GB PINNs, we used a series of fully connected network structures sequentially \((1,50,1)\), \((1,100,1)\), \((1,100,100,1)\), \((1,100,100,1)\) for the baseline and update neural networks, followed by a Fourier feature neural network \((1,50,50,1)\) with frequencies ranging from \(1\) to \(10\). The details of the Fourier feature method used in this study can be found in the appendix 6.1. The step size \(\rho_{m}\) in (5) was set to \(0.5^{i}\), where \(i=0,1,\cdots,5\) is the model index.

For each GB iteration, we train \(10,000\) steps using a dataset of \(10,000\) uniform random points in \((0,1)\). The weights in the loss function (3) are set to \(\omega_{e}=1\) and \(\omega_{b}=10\), respectively, and the batch size for PDE residual computation is \(10,000\).

The output of GB PINNs is shown in Figure 1, where the relative \(l^{2}\) error is found to be \(0.43\%\). The boundary layers at \(x=0\) and \(x=1\) are clearly visible in the solution, which is a result of the thinness of the layers and the almost right

\begin{table}
\begin{tabular}{l l} \hline \hline Name & Value \\ \hline Activation function & GeLU \\ Method to initialize the neural network & Xavier \\ Optimizer & Adam \\ learning rate & \(10^{-3}\) \\ learning rate decay period & \(10,000\) \\ learning rate decay rate & \(0.95\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Default Experiment Set up
angle curvature of the solution at these points. Despite the singularity present in the solution, GB PINNs were able to provide an accurate solution for this problem. To further highlight the contribution of GB PINNs in this example, an ablation study was conducted. A vanilla PINNs approach, using a network structure of \((1,100,100,100,1)\) and \(20,000\) training steps, was used to solve the same problem. The resulting relative \(l^{2}\) error was found to be \(11\%\), as shown in Figure 2. Additional experiments including ablation studies and comparisons can be found in appendix 6.2.

Furthermore, we demonstrate the robustness of our algorithm against the choice of network structure and arrangement through this example. We solve the problem using a variety of networks, including \((1,50,1)\), \((1,100,1)\), \((1,100,100,1)\), followed by a Fourier feature network \((1,100,100,100,1)\) with frequencies ranging from \(1\) to \(10\). The resulting relative \(l^{2}\) error is around \(1.1\%\), which is comparable to the previously mentioned result. It is possible that there exists a specific set of hyperparameters and configurations that would allow a single neural network to perfectly solve this problem. After all, by the universal approximation theorem [24], even a neural network with a simple structure possesses the ability to approximate a complicated function. However, the fine-tuning of hyperparameters is a common challenge in machine learning tasks and can consume significant computational resources. GB PINNs, on the other hand, do not require such fine-tuning as the multiple networks adjust the results automatically, ultimately saving effort and resources.

### 2D singular perturbation with boundary layers

In this example, we aim to solve the Eriksson-Johnson problem, which is a 2D convection-dominated diffusion equation. As previously noted in the literature, such as in [25], this problem necessitates the use of specialized finite element techniques in order to obtain accurate solutions, such as the Discontinuous Petrov Galerkin (DPG) finite element method.

Let \(\Omega=(0,1)^{2}\). The model problem is

\[-\varepsilon\Delta u+\frac{\partial u}{\partial x} =0 \text{in }\Omega,\] \[u =u_{0} \text{on }\partial\Omega.\]

The manufactured solution is

\[u(x,y)=\frac{e^{r_{1}(x-1)}-e^{r_{2}(x-1)}}{e^{-r_{1}}-e^{-r_{2}}}\sin(\pi y )\quad\text{ with }\quad r_{1,2}=\frac{-1\pm\sqrt{1+4\varepsilon^{2}\pi^{2}}}{-2\varepsilon}.\]

In this example, we set \(\varepsilon=10^{-3}\). To resolve this problem, we employ various neural network architectures, including \((2,50,1)\), \((2,100,1)\), \((2,100,100,1)\), \((2,100,100,1)\), \((2,100,100,1)\), and a Fourier feature network \((1,100,100,1)\) with frequencies ranging from \(1\) to \(50\). For each iteration of our GB algorithm, we train for \(20,000\) steps. We set the weights in (3) as \(\omega_{e}=1\) and \(\omega_{b}=10,000\), respectively. The batch sizes for PDE residuals and boundaries are set at \(10,000\) and \(50,00\), respectively. The predicted solution is visualized in Figure 3. We can see that

Figure 1: Prediction of singular perturbation problem by GB PINNs, \(\varepsilon=10^{-4}\). Left: predicted solution (black) v.s. ground truth (red). Right: pointwise error.
our model prediction is in good agreement with the ground truth, with a relative \(l^{2}\) error of \(1.03\%\). Notably, there is a boundary layer present on the right side of the boundary (\(x=1\)), which is not easily recognizable to the naked eye due to its thinness. However, GB PINNs are able to provide a reasonable degree of predictive accuracy even in this challenging scenario.

To further demonstrate the efficacy of our proposed method, we also attempted to solve this problem using a single fully connected neural network of architecture \((2,100,100,10)\), under the same hyperparameter settings as before. However, the resulting relative \(l^{2}\) error was \(57.7\%\). As can be seen in Figure 4, there is a significant discrepancy between the predicted solution and the reference solution. Additional experimental results, including an ablation study and comparisons, can be found in Appendix 6.2.

### 2D singular perturbation with an interior boundary layer

In this example, we address a 2D convection-dominated diffusion problem featuring curved streamlines and an interior boundary layer. The model problem is

\[-\varepsilon\Delta u+\beta\cdot\nabla u =f\qquad\text{ in }\Omega,\] \[u =u_{0}\qquad\text{ on }\partial\Omega,\]

Figure 3: Prediction of 2D singular perturbation with boundary problem by GB PINNs, \(\varepsilon=10^{-3}\). Left: predicted solution. Middle: ground truth. Right: pointwise error.

Figure 2: Prediction of singular perturbation problem by PINNs for ablation study, \(\varepsilon=10^{-4}\). Left: predicted solution (black) v.s. ground truth (red). Right: pointwise error.
with \(\beta=e^{x}(\sin(y),\cos(y))\) and \(f\), \(u_{0}\) are defined such that

\[u(x,y)=\arctan\left(\frac{1-\sqrt{x^{2}+y^{2}}}{\varepsilon}\right).\]

This example has been solved by DPG finite element method in [25]. A specific value of \(\epsilon=10^{-4}\) was chosen for the purpose of this study. The neural network architectures employed include \((2,200,200,200,1)\), \((2,100,100,1)\), \((2,100,100,1)\), and a Fourier feature network \((2,50,50,1)\) with frequency ranging from \(1\) to \(5\). The weights for the loss function in (3) were set as \(\omega_{e}=1\) and \(\omega_{b}=10,000\), respectively. The batch size for the PDE residual and boundary were set to \(10,000\) and \(5,000\), respectively. The results of this study are shown in Figure 5 and exhibit a relative \(l^{2}\) error of \(3.37\%\).

As a part of an ablation study, we resolve this problem using a fully connected neural network architecture of \((2,200,200,200,1)\), while maintaining the other configurations as same as the previous experiment. The relative \(l^{2}\) error obtained in this case is \(16\%\). The predictions and the corresponding errors are depicted in Figure 6. Additional experiments pertaining to the ablation study and comparisons can be found in the appendix, section 6.2.

### 2D nonlinear reaction-diffusion equation

In this example, we investigate the solution of a time-dependent nonlinear reaction-diffusion equation. As demonstrated in [12], conventional PINNs have been shown to be inadequate in accurately learning the solution of such equations.

Let \(\Omega=(0,2\pi)\). The model problem is

\[\frac{\partial u}{\partial t}-10\frac{\partial^{2}u}{\partial x^{2 }}-6u(1-u) =0,\quad x\in\Omega,t\in(0,1],\] \[u(x,0) =h(x)\quad x\in\Omega,\]

Figure 4: Prediction of 2D singular perturbation with boundary problem by PINNs, \(\varepsilon=10^{-3}\). Left: predicted solution. Middle: ground truth. Right: pointwise error.

Figure 5: Prediction of 2D singular perturbation with interior boundary problem by GB PINNs, \(\varepsilon=10^{-4}\). Left: predicted solution. Middle: ground truth. Right: pointwise error.
with periodic boundary conditions, where

\[h(x)=e^{-\frac{(\varepsilon-\pi)^{2}}{2(\nu/4)^{2}}}.\]

In order to impose an exact periodic boundary condition, we use \((\sin(x),\cos(x))\) as the spatial input instead of \(x\), while maintaining the temporal input unchanged. This eliminates the need for boundary loss. Additionally, we include an additional loss term for the initial condition in the loss function, equation (3). The neural network architecture utilized for this problem is \((2,200,200,200,1)\), \((2,100,100,100,1)\), \((2,100,100,1)\). The weights for the PDE residual and initial condition loss are set to \(1\) and \(1,000\), respectively. The batch sizes for the PDE residual and initial condition loss are \(20,000\) and \(1,000\), respectively. We present our results in Figure 7. The relative \(l^{2}\) error is \(0.58\%\). As shown in [12], the relative error for traditional PINNs with \(\rho=\nu=5\) is \(50\%\). A comparison between the exact solution and the PINNs' prediction can also be found in the figure.

In the aforementioned study, Krishnapriyan _et. al._[12] proposed a sequence-to-sequence learning approach to address this problem, achieving a relative \(l^{2}\) error of \(2.36\%\) for \(\rho=\nu=5\). The method required training the neural network for \(20\) steps. In contrast, our approach, which was implemented with \(\nu=10\) and \(\rho=6\), a more challenging scenario, required training only three networks. Despite this, our method achieved an error rate that was approximately four times lower than that of the previous study. This represents a significant improvement.

## 5 Conclusion

In this paper, we propose a GB PINNs algorithm, which utilizes multiple neural networks in sequence to predict solutions of PDEs. The algorithm is straightforward to implement and does not require extensive fine-tuning of hyperparameters. Additionally, the method is flexible and can be easily integrated with other PINNs techniques. Our experimental results demonstrate its effectiveness in solving a wide range of intractable PDE problems.

Figure 6: Prediction of 2D singular perturbation with interior boundary problem by PINNs, \(\varepsilon=10^{-4}\). Left: predicted solution. Middle: ground truth. Right: pointwise error.

Figure 7: Prediction of nonlinear reaction-diffusion equation by GB PINNs. Left: predicted solution. Middle: ground truth. Right: pointwise error.
However, it should be noted that the algorithm has some limitations. Firstly, it is not suitable for solving conservation laws with derivative blow-ups, such as the inviscid Burgers' equation and the Sod shock tube problem. This is due to the lack of sensitivity of these equations' solutions to PDE loss. The addition of more neural networks alone cannot overcome this issue. Secondly, the optimal combination of neural networks is not always clear, and the current experimental selection is mostly based on experience and prior estimation of the PDE problem. Further research into the theoretical and quantitative analysis of this method is an interesting direction for future work.

## Acknowledgments

This was supported in part by the US Department of Energy (grant DE-SC0019116) and the Air Force Office of Scientific Research (grant FA9550-20-1-0060), and DOE-ARPA grant DE-AR0001201.
