# Finite Element Operator Network for Solving Parametric PDEs

Jae Yong Lee
Center for Artificial Intelligence and Natural Sciences
Korea Institute for Advanced Study
Seoul, Republic of Korea
<EMAIL_ADDRESS>
&Seungchan Ko
Department of Mathematics
Inha University
Incheon, Republic of Korea
<EMAIL_ADDRESS>
Youngjoon Hong11footnotemark: 1
Department of Mathematical Sciences
KAIST
Daejeon, Republic of Korea
<EMAIL_ADDRESS>
Corresponding authors.

###### Abstract

Partial differential equations (PDEs) underlie our understanding and
prediction of natural phenomena across numerous fields, including physics,
engineering, and finance. However, solving parametric PDEs is a complex task
that necessitates efficient numerical methods. In this paper, we propose a
novel approach for solving parametric PDEs using a Finite Element Operator
Network (FEONet). Our proposed method leverages the power of deep learning in
conjunction with traditional numerical methods, specifically the finite
element method, to solve parametric PDEs in the absence of any paired input-
output training data. We demonstrate the effectiveness of our approach on
several benchmark problems and show that it outperforms existing state-of-the-
art methods in terms of accuracy, generalization, and computational
flexibility. Our FEONet framework shows potential for application in various
fields where PDEs play a crucial role in modeling complex domains with diverse
boundary conditions and singular behavior. Furthermore, we provide theoretical
convergence analysis to support our approach, utilizing finite element
approximation in numerical analysis.

##  1 Introduction

Solving partial differential equations (PDEs) is vital as they serve as the
foundation for understanding and predicting the behavior of a range of natural
phenomena [1; 2]. From fluid dynamics, heat transfer, to electromagnetic
fields, PDEs provide a mathematical framework that allows us to comprehend and
model these complex systems [3; 4]. Their significance extends beyond the
realm of physics and finds applications in fields such as finance, economics,
and computer graphics [5; 6]. In essence, PDEs play a crucial role in
advancing our understanding of the world around us and have a wide-reaching
impact in various industries and fields of study.

Numerical methods play a crucial role in approximating solutions to PDEs when
exact solutions are not feasible. This is particularly important for complex
systems that cannot be easily solved using traditional methods [7; 8]. Various
techniques, such as finite difference, finite element, and finite volume
methods, are utilized to develop these numerical methods [9; 10; 11]. In
particular, the finite element method (FEM) is widely employed in engineering
and physics, where the domain is divided into small elements and the solution
within each element is approximated using polynomial functions. The
mathematical analysis of the FEM has undergone significant development,
resulting in a remarkable enhancement in the method’s reliability. The FEM is
particularly advantageous in handling irregular geometries and complex
boundary conditions, making it a valuable tool for analyzing structures, heat
transfer, fluid dynamics, and electromagnetic problems [12; 13]. However, it
is important to note that the implementation of numerical methods often comes
with a significant computational cost.

The use of machine learning in solving PDEs has gained significant traction in
recent years. By integrating deep neural networks and statistical learning
theory into numerical PDEs, a new field known as scientific machine learning
has emerged, presenting novel research opportunities. The application of
machine learning to PDEs can be traced back to the 1990s when neural networks
were employed to approximate solutions [14]. More recently, the field of
physics-informed neural networks (PINNs) has been developed, where neural
networks are trained to learn the underlying physics of a system, enabling us
to solve PDEs and other physics-based problems [15; 16; 17; 18; 19]. Despite
their advantages, these approaches come with limitations. One prominent
drawback of PINN methods is that they are trained on a single instance of
input data, encompassing initial conditions, boundary conditions, and external
force terms. As a consequence, if the input data changes, the entire training
process must be repeated, posing challenges in making real-time predictions
for varying input data. This limitation restricts the applicability of PINNs
in dynamic and adaptive systems where the input data is subject to change.

Operator learning, a recent research area, tackles the challenge by using
data-driven methods to learn mathematical operators governing a physical
system for solving parametric PDEs [20; 21; 22; 23; 24; 25; 26; 27]. These
methods rely on pre-generated databases of input-output pairs obtained through
analytical or numerical solutions. Real-time predictions for varying input
data can then be made using these databases. However, generating a reliable
training dataset requires extensive numerical computations, which can be
computationally inefficient and time-consuming. Obtaining a sufficiently large
dataset is also challenging, especially for systems with complex geometries or
nonlinear equations. Additionally, using neural networks as the solution space
poses challenges in imposing boundary values and can lead to less accurate
solutions. Researchers are actively working on developing methods to enhance
the efficiency and accuracy of these data-driven approaches [28].

Figure 1: Domain triangulation for three different complex geometries.

To address the limitations mentioned above, in the present paper, a novel
approach is proposed for solving diverse parametric PDEs. This approach
eliminates the need for a training dataset of solutions by introducing a novel
method of learning the coefficients. Specifically, we introduce the Finite
Element Operator Network (FEONet) which leverages the finite element
approximation to learn the solution operator in a novel way, eliminating the
need for a training dataset of solutions to solve various parametric PDEs. In
the FEM framework, the numerical solution is approximated by the linear
combination of nodal coefficients, $\alpha_{k}$, and the nodal basis,
$\phi_{k}({\mathbf{x}})$. The nodal basis in the FEM consists of piecewise
polynomials defined by the finite set of nodes of a mesh, which exhibit both
near orthogonality and local support such that

${u_{h}}({\mathbf{x}})=\sum\alpha_{k}\phi_{k}({{\mathbf{x}}}),\quad\mathbf{x}\in\mathbb{R}^{d},$
(1)

where ${u_{h}}$ is the FEM solution for the given PDEs. Motivated by (1), the
FEONet is able to predict numerical solutions to the PDEs when given the
initial conditions, external forcing functions, or PDE coefficients as inputs.
As the neural network only predicts the nodal coefficients, this approach
eliminates the need for extensive numerical computations to generate a large
dataset and also allows for real-time predictions for varying input data.
Therefore, the FEONet can learn multiple instances of the solutions of PDEs,
making it a versatile approach for solving different types of PDEs in complex
domains; see e.g. Figure 1. The loss function of the FEONet is designed based
on the residual quantity of the finite element approximation, similar to the
FEM [29; 30]. This allows the FEONet to accurately approximate the solutions
of PDEs, while also ensuring that the boundary conditions are exactly
satisfied. To construct an approximation for the solution of PDEs, the FEONet
borrows a framework from the FEM. This approach involves inferring
coefficients, denoted as $\widehat{\alpha}_{k}$, which are then used to
construct the linear combination $\sum\widehat{\alpha}_{k}\phi_{k}$ to
approximate the solution of PDEs. Since each basis function satisfies the
exact boundary condition, the predicted solution will also satisfy the exact
boundary condition. Additionally, since the high-order FEM is available to
generate our framework in (1), the predicted numerical solution is expected to
yield relatively small errors compared to other machine learning-based
approaches.

Another key contribution of this study is to introduce a novel learning
architecture designed to accurately solve convection-dominated singularly
perturbed problems involving strong boundary layer phenomena. These types of
problems pose significant challenges for traditional numerical methods due to
sharp transitions inside thin layers caused by a small diffusive parameter. By
adapting theory-guided methods [31], the network effectively captures the
behavior of boundary layers, an area that has been underexplored in recent
machine learning approaches.

The main contributions of the paper are summarized as follows:

  * •

Our proposed method for learning operators, FEONet (Figure 3), is novel, as it
enables the learning of multiple solutions for parametric PDEs without the
requirement of a input-output training dataset.

  * •

The proposed neural network is versatile in its application, allowing for the
inclusion of various domains, (Figure 4), as well as diverse equations with
boundary conditions, including both Dirichlet and Neumann boundary conditions
(Table 1).

  * •

The FEONet includes rigorous convergence analysis via finite element
approximation on P1 and P2 elements (Theorem 3.1), and it is expected to have
low generalization errors.

  * •

Lastly, the proposed neural network is a specialized learning architecture
tailored to accurately solve convection-dominated singular perturbation
problems that feature robust boundary layers (Figure 5).

##  2 Related works

#### Physics-informed neural network (PINN).

The PINNs, developed by Raissi et al. [19], utilize a neural network to
efficiently solve PDEs with little or no data. The key feature of PINN is its
ability to minimize PDE residual loss by enforcing physical constraints, with
output fields automatically differentiated with respect to input coordinates.
PINNs have been successfully applied in various fields such as material
science and biophysics [32; 33]. However, one of the key limitations of the
PINN method is that it is trained on a single set of input data, which
includes initial conditions, boundary conditions, and external force terms.

#### Operator learning.

Operator learning utilizes a dataset of input-output pairs obtained from an
existing solver to train a model to learn the solution operator [34; 35; 36;
37]. Operator learning aims to improve the efficiency and accuracy of solving
PDEs by learning the underlying patterns and relationships in the input-output
pairs, allowing for real-time predictions for varying input data [38].
Kovachki et al. [39] proposed an iterative neural operator structure to learn
the solution operator of PDEs. Based on the universal approximation theorem
for the operator, Lu et al. [20] developed a Deep operator network (DeepONet)
architecture. In recent works [23; 24; 40; 41; 42], a message-passing neural
PDE solver capable of handling diverse problem properties has been developed.
However, this approach faces challenges such as the requirement for large
amounts of training data, limited generalization ability, and extrapolation to
unseen conditions. Researchers are currently addressing these challenges by
incorporating physics or PDE constraints and designing more expressive models.

#### Physics-informed operator learning.

The integration of PINNs and operator learning has led to the development of
innovative methods [28; 43] such as the Physics-informed neural operator
(PINO) [44] and Physics-informed DeepONet (PIDeepONet) [22]. These methods aim
to capitalize on the strengths of both PINNs and operator learning by
embedding physical equations in the loss function of the neural operator.
However, there are certain limitations to these methods, such as the inability
to effectively address the boundary layer issue, which plays a crucial role in
real-world problems. Additionally, PIDeepONet has been observed to have
relatively high generalization errors when no training data is available.
Furthermore, PINO has mainly been applied to problems with a square-shaped
domain.

##  3 Method

In this section, we present our proposed method, the FEONet, for solving
parametric PDEs. We begin by providing a concise overview of the FEMs, which
serve as the foundation of our approach. Next, we introduce the FEONet by
demonstrating how neural network techniques can be seamlessly integrated into
FEMs to solve the parametric PDEs. Throughout our method description and
experiments, we focus on the PDEs of the following form: for uniformly
elliptic coefficients $\boldsymbol{a}(\mathbf{x})$ and $\varepsilon>0$,

$-\varepsilon\,{\rm{div}}\,(\boldsymbol{a}(\mathbf{x})\nabla
u)+\mathcal{F}(u)=f\quad{\rm{in}}\,\,D.$ (2)

Here $\mathcal{F}$ can be either linear or nonlinear, and both cases will be
covered in the experiments performed in Section 4. As will be described in
more detail later, we shall propose an operator-learning-type method which can
provide real-time solution predictions whenever the input data of the PDE
changes.

###  3.1 Finite element method (FEM)

The FEM is a general class of techniques for the numerical solution of PDEs.
FEM is based on the variational formulation of the given PDE (2): find $u\in
V$ satisfying

$B[u,v]:=\varepsilon\int_{D}\boldsymbol{a}(\mathbf{x})\nabla u\cdot\nabla
v\,\mathrm{d}\mathbf{x}+\int_{D}\mathcal{F}(u)v\,\mathrm{d}\mathbf{x}=\int_{D}fv\,\mathrm{d}\mathbf{x}=:\ell(v)\quad{\text{for
all}}\,\,v\in V,$ (3)

where the solution and test function space $V$ is usually chosen to be an
infinite-dimensional function space. In the theory of FEM, the first step is
to define a triangulation of the given domain $D\subset\mathbb{R}^{d}$. If
$d=1$ and $D=[a,b]$ for some $a$, $b\in\mathbb{R}$, we set
$a=x_{0}<x_{1}<x_{2}<\cdots<x_{K}=b$ and each interval ($1$-simplex)
$[x_{i-1},x_{i}]$ defines an element.

Figure 2: Shape of nodal basis functions and approximation of $y=\sin x$ for
(a) P1-element and (b) P2-element.

For $d=2$, the triangulation is referred to as a tessellation of $D$ into a
finite number of closed triangles ($2$-simplex) $T_{i}$, $i=1,\cdots,K$, whose
interiors are pairwise disjoint, and for $i\neq j$ with $T_{i}\cap T_{j}$ is
nonempty, $T_{i}\cap T_{j}$ is either a common vertex or a common edge of
$T_{i}$ and $T_{j}$. In this case, each triangle defines an element and the
vertices $\\{\mathbf{x}_{i}\\}$ of triangles are called nodes (Figure 1). The
triangulation in higher dimension $d\geq 3$ can be defined in a similar manner
using $d$-simplex as an element. Letting $h_{T}$ be the largest edge of a
triangle $T$, we define a finite element parameter $h>0$ as the longest among
the $h_{T}$. In addition, let us denote by $S_{h}$, the space of all
continuous functions $v_{h}$ defined on $D$ such that the restriction of
$v_{h}$ to an arbitrary triangle is a polynomial. Then we define our finite-
dimensional ansatz space as $V_{h}=S_{h}\cap V$. We shall denote the set of
all vertices in the triangulation by $\\{\mathbf{x}_{i}\\}$, and define the
so-called nodal basis $\\{\phi_{j}\\}$ for $V_{h}$, defined by
$\phi_{j}(\mathbf{x}_{i})=\delta_{ij}$. If we use piecewise linear nodal basis
functions, we call it the P$1$-element method, while the adoption of piecewise
quadratic polynomials is called as the P$2$-element method (Figure 2). Note
that the dimension of $V_{h}$ depends on the triangulation of $D$, and hence,
depends also on the finite element parameter $h>0$.

The idea of the FEM is to approximate $V$ with a finite-dimensional subspace
$V_{h}$ above, generated by the basis functions
$\\{\phi_{1},\phi_{2},\cdots,\phi_{N(h)}\\}$, so that the given problem
becomes numerically feasible. More precisely, we aim to compute the
approximate solution $u_{h}\in V_{h}$ based on the, so called, Galerkin
approximation

$B[u_{h},v_{h}]=\ell(v_{h})\quad{\text{for all}}\,\,v_{h}\in V_{h}.$ (4)

If we write the finite element solution

$u_{h}=\sum^{N(h)}_{k=1}\alpha_{k}\phi_{k},\quad\alpha_{i}\in\mathbb{R},$ (5)

the Galerkin approximation (4) is transformed in to the linear algebraic
system

$A\alpha=F\quad{\text{with}}\,\,A_{ik}:=B[\phi_{k},\phi_{i}]\,\,{\rm{and}}\,\,F_{i}:=\ell(\phi_{i}),$
(6)

where the matrix $A\in\mathbb{R}^{N(h)\times N(h)}$ is invertible provided
that the given PDEs have some suitable structure. Therefore, we can determine
the coefficients $\\{\alpha_{k}\\}^{N(h)}_{k=1}$ by solving the system of
linear equations (6), and consequently, produce the approximate solution
$u_{h}\in V_{h}$ as defined in (5).

Once we compute the approximate solution $u_{h}$ of $u$, we can measure how
close $u_{h}\in V_{h}$ is to $u\in V$. If the given PDE has some good
structure, the weak solution $u$ is sufficiently smooth, and the triangulation
is reasonable, we can derive the following error estimates: for some constants
$C_{1}$, $C_{2}>0$ independent of $h>0$, we have

(P1-element method) $\displaystyle\|u-u_{h}\|_{L^{2}(D)}\leq C_{1}h^{2},$ (7)
(P2-element method) $\displaystyle\|u-u_{h}\|_{L^{2}(D)}\leq C_{2}h^{3}.$

For a detailed discussion of the error analysis, see, for example, Theorem
5.4.8 of [30]. From the estimates (7), we can confirm that the approximate
solution $u_{h}$ computed by the FEM converges to the true solution $u$ as
$h\rightarrow 0$ and the convergence rate for the P$1$ method is $2$, while
the P$2$-method achieves the higher convergence rate $3$.

Figure 3: Schematic diagram of the Finite Element Operator Network (FEONet).

###  3.2 Finite Element Operator Network (FEONet)

The FEONet is based on the FEM described in the previous section. In our
scheme, an input for the neural network is the data for the given PDEs. In
this paper, we choose external force as a prototype input feature, but the
same scheme can be developed in a straightforward manner also with other types
of data, e.g., boundary conditions, diffusion coefficients or initial
conditions for time-dependent problems. For this purpose, the external forcing
terms are parametrized by the random parameter $\boldsymbol{\omega}$ contained
in the (possibly high-dimensional) parameter space $\Omega$. For each
realization $f(\mathbf{x};\boldsymbol{\omega})$ (and hence for each load
vector $F(\boldsymbol{\omega})$ defined in (6)), instead of computing the
coefficients from the linear algebraic system (6), we approximate the
coefficients $\alpha$ by a deep neural network. More precisely, these input
features representing external forces pass through the deep neural network and
the coefficients $\\{\widehat{\alpha}_{k}\\}$ is generated as an output of the
neural network. We then reconstruct the solution

$\widehat{u}_{h}(\mathbf{x};\boldsymbol{\omega})=\sum_{k=1}^{N(h)}\widehat{\alpha}_{k}(\boldsymbol{\omega})\phi_{k}(\mathbf{x}).$
(8)

For the training, we use a residual of the variational formulation (4) and
define the loss function by summing over all basis functions: for randomly
chosen parameters
$\boldsymbol{\omega}_{1},\cdots,\boldsymbol{\omega}_{M}\in\Omega$,

$\mathcal{L}^{M}(\widehat{u}_{h})=\frac{1}{M}\sum^{M}_{m=1}\sum^{N(h)}_{i=1}|B[\widehat{u}_{h}(\mathbf{x};\boldsymbol{\omega}_{m}),\phi_{i}(\mathbf{x})]-\ell(\phi_{i}(\mathbf{x});\boldsymbol{\omega}_{m})|^{2}.$
(9)

For each training epoch, once the neural network parameters are updated in the
direction of minimizing the loss function (9), the external force
$f(\mathbf{x};\boldsymbol{\omega})$ passes through this updated neural network
to generate more refined coefficients, and this procedure is repeated until
the sufficiently small loss is achieved. A schematic diagram of the FEONet
algorithm is depicted in Figure 3.

###  3.3 Convergence of FEONet

We address a convergence result of the FEONet which provides theoretical
justification for the proposed method. As outlined in the previous section,
the proposed method is based on the FEM and the finite element approximation
$u_{h}$ in (5) can be considered as an intermediate solution between the true
solution $u$ of (3) and the approximate solution $\widehat{u}_{h}$ predicted
by the FEONet. Consequently, in order for an error analysis, we simply split
the error $u-\widehat{u}_{h}$ into two parts; namely,

$u-\widehat{u}_{h}=(u-u_{h})+(u_{h}-\widehat{u}_{h})=:({\rm{I}})+({\rm{II}}).$
(10)

The first error (I) can be controlled via the classical theory of FEM, e.g.
through the error estimate (7). The second error (II) is a core part which
needs a novel theoretical investigation. The parameters of
$\widehat{u}_{h}=:\widehat{u}_{h,n,M}$ concerning the convergence is the
parameter representing neural network architecture and the number of sampling
points $\boldsymbol{\omega}\in\Omega$ which we denote by $n\in\mathbb{N}$ and
$M\in\mathbb{N}$ respectively. Here, larger $n\in\mathbb{N}$ implies more
approximation power of neural networks. We will prove the following theorem
concerning the second error (II), which guarantees that our method is reliable
and the theorem forms the theoretical foundation of the proposed numerical
scheme. The precise mathematical statement and proof can be found in Theorem
B.7.

######  Theorem 3.1.

If we let $u_{h}$ be the finite element approximation of the true solution $u$
and $\widehat{u}_{h,n,M}$ be the approximate solution computed by the FEONet,
then we have

$\mathbb{E}_{\boldsymbol{\omega}}\left[\|u_{h}-\widehat{u}_{h,n,M}\|_{L^{2}(D)}^{2}\right]\rightarrow
0\quad{\rm{as}}\,\,n,\,M\rightarrow\infty.$ (11)

Figure 4: Solution profiles for complex geometries.

It is noteworthy that $\widehat{u}_{h}$ shares some interesting properties
with $u_{h}$. As it will be made clear in the next section, we may conclude
from experiments that $\|{(\rm{I})}\|\gg\|(\rm{II})\|$, and the total error
$\|u-\widehat{u}_{h}\|$ is sufficiently close to the finite element error
$\|(\rm{I})\|$. Therefore, as we can see from Figure 6, for large enough $n$,
$M\in\mathbb{N}$, $\|u-\widehat{u}_{h}\|$ has the convergence rates close to
$2$ and $3$ (for P$1$ and P$2$ respectively) with respect to $h>0$, which
agrees well with the theoretical result (7) for $\|(\rm{I})\|$.

##  4 Numerical experiments

In this section, we demonstrate the experimental results of the FEONet across
various settings for PDE problems. We randomly generate 3000 input samples for
external forces and train the FEONet in an unsupervised manner without the
usage of pre-computed input-output $(f,u)$ pairs. We then evaluate the
performance of our model with another randomly generated test set consisting
of 3000 pairs of $f$ and the corresponding solutions $u$. For this test data,
we employed the FEM with a sufficiently fine mesh discretization to ensure
that the numerical solutions can be considered as true solutions. The forcing
term $f$ is randomly generated. Although there are other operator-learning-
based methods, such as PINO [44] and PIDeepONet [22], which can be trained in
an unsupervised manner, these models have certain limitations when dealing
with complex domain shapes (refer to Appendix C). As a result, our main focus
is primarily on comparing our model to DeepONet [20], considering the varying
numbers of training data. The DeepONet is known to provide accurate results
for complex domains. In our experiments, we prepare a total of 30/300/3000
pairs of $(f,u)$ by applying the FEM on sufficiently fine meshes. The
triangulation of domains and the generation of nodal basis for the FEM are
performed using FEniCS [45; 46]. A detailed description of the model’s
hyperparameters and experimental setup can be found in Appendix D.

###  4.1 Comparison in various PDE problems

To evaluate the computational flexibility of the FEONet, we conducted a
sequence of experiments covering diverse domains, boundary conditions, and
equations.

#### Various domains.

We consider the 2D convection-diffusion equation as

$\displaystyle-\varepsilon\Delta u+\textbf{v}\cdot\nabla u$
$\displaystyle=f(x,y),$ $\displaystyle(x,y)\in D,$ (12) $\displaystyle u(x,y)$
$\displaystyle=0,$ $\displaystyle(x,y)\in\partial D.$

For our numerical experiments, we fixed the values of $\varepsilon=0.1$ and
$\textbf{v}=(-1,0)$ across various domains $D$, which included a circle, a
square with a hole, and a polygon (see Figure 1). The second to fourth columns
of Table 1 display the mean errors of the test data for the FEONet
(unsupervised) and DeepONet (30/300/3000 training data). When the DeepONet is
trained with 30 or 300 data pairs, the FEONet consistently outperforms it in
terms of error rates. Moreover, even when trained with 3000 data pairs, the
FEONet achieves comparable or even lower errors compared to the DeepONet. It
is worth noting that for the relatively more complex Domain II and Domain III,
DeepONet struggles to make accurate predictions even when trained with 3000
data pairs. Figure 4 displays the solution profile, clearly demonstrating that
FEONet achieves higher accuracy in predicting the solution compared to
DeepONet, especially for complex domains.

Table 1: Mean Rel. $L^{2}$ test errors $\mathbf{(\times 10^{-2})}$ with
standard deviations $\mathbf{(\times 10^{-2})}$ for the 3000 test set on
diverse PDE problems. Five training trials are performed independently. Domain
I (circle), Domain II (square with a hole), Domain III (polygon) with Eq.
(12), BC I (Dirichlet), BC II (Neumann) with Eq. (13), Eq I (second-order
linear equation (14)) and Eq II (Burgers’ equation (15)).

Model (#Train data) | Domain I | Domain II | Domain III | BC I | BC II | Eq I | Eq II
---|---|---|---|---|---|---|---
Ours (unsupervised) | 1.24$\pm$0.00 |  1.76$\pm$0.03 |  0.51$\pm$0.00 |  0.13$\pm$0.01 |  0.32$\pm$0.03 |  0.13$\pm$0.01 | 0.54$\pm$0.07
DON (w/30) | 27.15$\pm$1.16 | 51.21$\pm$3.58 | 53.92$\pm$4.59 | 21.75$\pm$1.19 | 22.75$\pm$1.05 | 24.38$\pm$1.37 | 10.26$\pm$0.14
DON (w/300) | 2.10$\pm$0.75 | 5.62$\pm$0.37 | 6.22$\pm$0.96 | 0.68$\pm$0.11 | 0.96$\pm$0.06 | 0.76$\pm$0.10 |  0.20$\pm$0.09
DON (w/3000) |  0.69$\pm$0.17 | 4.75$\pm$0.75 | 6.20$\pm$1.00 | 0.53$\pm$0.36 | 0.33$\pm$0.09 | 0.33$\pm$0.27 | 0.24$\pm$0.13

#### Various boundary conditions.

We also tested the performance of each model, when the given PDEs are equipped
with either Dirichlet or Neumann boundary conditions:

$\displaystyle-\varepsilon u_{xx}+bu_{x}$ $\displaystyle=f(x),$ $\displaystyle
x\in D,$ (13) $\displaystyle u(x)=0\quad\text{or}\quad u_{x}(x)$
$\displaystyle=0,$ $\displaystyle x\in\partial D,$

where $\varepsilon=0.1$ and $b=-1$. For the Neumann boundary condition, we add
$u$ to the left-hand side of (13) to guarantee the uniqueness of solution up
to constant. The FEONet has an additional significant advantage to make the
predicted solutions satisfy the exact boundary condition, by selecting the
appropriate coefficients for the basis functions. Using the FEONet without any
input-output training data, we are able to obtain similar or slightly higher
accuracy compared to the DeepONet with 3000 training data, for both
homogeneous Dirichlet or Neumann boundary conditions as shown in Table 1.

#### Various equations.

We performed some additional experiments applying the FEONet to the various
equations: (i) general second-order linear PDE with variable coefficients
defined as

$\displaystyle-\varepsilon u_{xx}+b(x)u_{x}+c(x)u$ $\displaystyle=f(x),$
$\displaystyle x\in D,$ (14) $\displaystyle u(x)$ $\displaystyle=0,$
$\displaystyle x\in\partial D,$

where $\varepsilon=0.1$, $b(x)=x^{2}+1$ and $c(x)=x$; and (ii) the Burgers’
equation defined as

$\displaystyle-u_{xx}+uu_{x}$ $\displaystyle=f(x),$ $\displaystyle x\in D,$
(15) $\displaystyle u(x)$ $\displaystyle=0,$ $\displaystyle x\in\partial D.$

Since the Burgers’ equation is nonlinear, an iterative method is required for
classical numerical schemes including the FEM, which usually causes some
computational costs. Furthermore, it is more difficult than linear equations
to obtain training data for nonlinear equations. One important advantage of
the FEONet is its ability to effectively learn the nonlinear structure without
requiring any training data. Once the training process is completed, the model
can provide accurate real-time solution predictions without relying on
iterative schemes. See Appendix D.6 for the detailed explanation of the loss
function constructed for the nonlinear Burgers’ equation. The 7th and 8th
columns of Table 1 highlight the effectiveness of the FEONet, predicting the
solutions with reasonably low errors for various equations. Although the 8th
column of Table 1 shows that the DeepONet with 3000 training data has a lower
error compared to our model, it demonstrates the strength of the FEONet that
achieves a similar level of accuracy even for the nonlinear PDE, that can be
trained in an unsupervised manner.

###  4.2 Singular perturbation problem

Table 2: Comparisons of various deep-learning-based models for solving the singular perturbation problem. Each model is trained five times independently. Model | Ours | PINN | DeepONet w/30 | PIDeepONet
---|---|---|---|---
Unsupervised | Yes | Yes | NO | Yes
Multiple instance | Yes | No | Yes | Yes
Mean Rel. $L^{2}$ test error$\pm$std |  0.0132$\pm$0.0091 | 1.3827$\pm$0.7580 | 0.2303$\pm$0.0074 | 0.5713$\pm$0.0007
Figure 5: Input function $f$ and the corresponding approximate solution $u$
obtained by various deep-learning-based methods for the singular perturbation
problem with $\varepsilon=10^{-5}$.

One of the notable advantages of the FEONet, which predicts coefficients
through well-defined basis functions, is its applicability to problems
involving singularly perturbed PDEs. This is accomplished by incorporating
additional basis functions guided by theoretical considerations. For this we
consider

$\displaystyle-\varepsilon u_{xx}+bu_{x}$ $\displaystyle=f(x),$ $\displaystyle
x\in D,$ (16) $\displaystyle u(x)$ $\displaystyle=0,$ $\displaystyle
x\in\partial D,$

where $\varepsilon\ll 1$. For the implementation, we assign the values of
$b=-1$ and $\varepsilon=10^{-5}$. For instance, in the right panel of Figure
5, the red dotted line (ground-truth) shows the solution profile which
produces a sharp transition near $x=-1$. To capture the boundary layer, we
utilize an additional basis function known as the corrector function in
mathematical analysis, defined as

$\phi_{\text{cor}}(x):=e^{-(1+x)/\varepsilon}-(1-(1-e^{-2/\varepsilon})(x+1)/2);$
(17)

see Appendix D.7 for more detailed information regarding the derivation of the
additional basis function. By integrating the boundary layer element into the
finite element space, we establish the proposed enriched Galerkin space for
utilization in the FEONet. As the corrector basis is incorporated alongside
the conventional nodal basis functions in the FEM, we also predict the
additional coefficient originating from the corrector basis. Note that the
computational overhead remains minimal as the enriched basis exclusively
encompasses boundary elements. Table 2 summarizes the characteristics of
various operator learning models and the corresponding errors when solving the
singularly perturbed PDE described in (16). The performance of the FEONet in
accurately solving boundary layer problems surpasses that of other models,
even in an unsupervised manner, across multiple instances. The incorporation
of a very small $\varepsilon$ in the residual loss of the PDE makes operator
learning challenging or unstable when using PINN and PIDeepONet. In fact,
neural networks inherently possess a smooth prior, which can lead to
difficulties in handling boundary layer problems. In contrast, FEONet utilizes
theory-guided basis functions, enabling the predicted solution to accurately
capture the sharp transition near the boundary layer. The PIDeepONet faces
difficulties in effectively learning both the residual loss and boundary
conditions, while PINNs struggle to predict the solution (refer to Figure 11).
Figure 5 exhibits the results of operator learning for the singular
perturbation problem. As depicted in the zoomed-in graph of Figure 5, the
FEONet stands out as the only model capable of accurately capturing sharp
transitions, whereas other models exhibit noticeable errors.

###  4.3 Further experiments

Figure 6: Rel. $L^{2}$ error of the FEONet with respect to the number of
elements.

#### Convergence rate.

One of the intriguing aspects is that the theoretical results on the
convergence error rate of the FEONet (as presented in (7)) are also observable
in the experimental results. Figure 6 shows the relationship between the test
error and the number of elements, utilizing both P1 and P2 basis functions. In
1D problems, the convergence rate is approximately -2 for P1 basis functions
and -2.9 for P2 basis functions. Meanwhile, in 2D problems, the convergence
rate is approximately -2 for P1 basis functions and -2.5 for P2 basis
functions. These remarkable trends confirm the theoretical convergence rates
observed in the experimental results.

#### Generalization errors.

Figure 7 displays a comparison of the training and test errors between FEONet
(unsupervised) and DeepONet (with 30 training data pairs) across different
settings as performed in Section 4 (Table 1). The disparity between the
training error and the test error serves as a reliable indicator of how
effectively the models generalize the underlying phenomenon. As depicted in
Figure 7, the generalization error of FEONet is notably smaller in comparison
to that of DeepONet when trained with 30 data pairs.

##  5 Conclusion

In this paper, we introduce the FEONet, a novel deep learning framework that
approximates nonlinear operators in infinite-dimensional Banach spaces using
finite element approximations. Our operator learning method for solving PDEs
is both simple and remarkably effective, resulting in significant improvements
in predictive accuracy, domain flexibility, and data efficiency compared to
state-of-the-art techniques. Furthermore, we demonstrate that the FEONet is
capable of learning the solution operator of parametric PDEs, even in the
absence of paired input-output training data, and accurately predicting
solutions that exhibit singular behavior in thin boundary layers. Although the
external force was considered as the input function in the experiment, the
FEONet can be easily extended to operator learning for other input functions,
such as boundary conditions, variable coefficients, or initial conditions of
time-dependent PDEs (see Appendix D.4). We are confident that the FEONet can
provide a new and promising approach to simulate and model intricate,
nonlinear, and multiscale physical systems, with a wide range of potential
applications in science and engineering.

Figure 7: Rel. $L^{2}$ train and test error comparison between the FEONet and
the DeepONet with 30 training data. Results are averaged over five independent
training trials, and black bars indicate the standard deviation.

## References

  * Courant and Hilbert [1953] Richard Courant and David Hilbert.  _Methods of Mathematical Physics_.  Interscience, New York, 2nd edition, 1953.  ISBN 978-0-471-50447-4.  Translated from German: Methoden der mathematischen Physik I, 2nd ed, 1931 [15].
  * Bender [2000] Edward A Bender.  _An introduction to mathematical modeling_.  Courier Corporation, 2000.
  * Gershenfeld and Gershenfeld [1999] Neil A Gershenfeld and Neil Gershenfeld.  _The nature of mathematical modeling_.  Cambridge university press, 1999.
  * Lin and Segel [1988] Chia-Chiao Lin and Lee A Segel.  _Mathematics applied to deterministic problems in the natural sciences_.  SIAM, 1988.
  * Simon et al. [1994] Carl P Simon, Lawrence Blume, et al.  _Mathematics for economists_ , volume 7.  Norton New York, 1994.
  * Mortenson [1999] Michael E Mortenson.  _Mathematics for computer graphics applications_.  Industrial Press Inc., 1999.
  * Burden et al. [2015] Richard L Burden, J Douglas Faires, and Annette M Burden.  _Numerical analysis_.  Cengage learning, 2015.
  * Atkinson [1991] Kendall Atkinson.  _An introduction to numerical analysis_.  John wiley & sons, 1991.
  * Larson and Bengzon [2013] Mats G Larson and Fredrik Bengzon.  _The finite element method: theory, implementation, and applications_ , volume 10.  Springer Science & Business Media, 2013.
  * Strikwerda [2004] John C Strikwerda.  _Finite difference schemes and partial differential equations_.  SIAM, 2004.
  * LeVeque [2002] Randall J. LeVeque.  _Finite-Volume Methods for Hyperbolic Problems_.  Cambridge University Press, 2002.
  * Zienkiewicz and Taylor [2000] Olek C Zienkiewicz and Robert L Taylor.  _The Finite Element Method: Its Basis and Fundamentals_.  Butterworth-Heinemann, 2000.
  * Hughes [2000] Thomas J. R. Hughes.  _The Finite Element Method: Linear Static and Dynamic Finite Element Analysis_.  Dover Publications, 2000.
  * Lagaris et al. [1998] 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.
  * Lu et al. [2021a] Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis.  DeepXDE: A deep learning library for solving differential equations.  _SIAM Review_ , 63(1):208–228, 2021a.
  * Karniadakis et al. [2021] George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang.  Physics-informed machine learning.  _Nature Reviews Physics_ , 3(6):422–440, 2021\.
  * Yang et al. [2021] Liu Yang, Xuhui Meng, and George Em Karniadakis.  B-PINNs: Bayesian physics-informed neural networks for forward and inverse PDE problems with noisy data.  _Journal of Computational Physics_ , 425:109913, 2021.
  * Meng et al. [2020] Xuhui Meng, Zhen Li, Dongkun Zhang, and George Em Karniadakis.  PPINN: Parareal physics-informed neural network for time-dependent PDEs.  _Computer Methods in Applied Mechanics and Engineering_ , 370:113250, 2020.
  * Raissi et al. [2019] 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\.
  * Lu et al. [2021b] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis.  Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators.  _Nature Machine Intelligence_ , 3(3):218–229, 2021b.
  * Li et al. [2021a] Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Fourier Neural Operator for Parametric Partial Differential Equations.  In _International Conference on Learning Representations_ , 2021a.  URL https://openreview.net/forum?id=c8P9NQVtmnO.
  * Wang et al. [2021] Sifan Wang, Hanwen Wang, and Paris Perdikaris.  Learning the solution operator of parametric partial differential equations with physics-informed DeepONets.  _Science Advances_ , 7(40):eabi8605, 2021.
  * Brandstetter et al. [2022] Johannes Brandstetter, Daniel E. Worrall, and Max Welling.  Message Passing Neural PDE Solvers.  In _International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=vSix3HPYKSU.
  * Lienen and Günnemann [2022] Marten Lienen and Stephan Günnemann.  Learning the Dynamics of Physical Systems from Sparse Observations with Finite Element Networks.  In _International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=HFmAukZ-k-2.
  * Lu et al. [2022] Lu Lu, Xuhui Meng, Shengze Cai, Zhiping Mao, Somdatta Goswami, Zhongqiang Zhang, and George Em Karniadakis.  A comprehensive and fair comparison of two neural operators (with practical extensions) based on fair data.  _Computer Methods in Applied Mechanics and Engineering_ , 393:114778, 2022.
  * Fanaskov and Oseledets [2022] Vladimir Fanaskov and Ivan Oseledets.  Spectral Neural Operators.  _arXiv preprint arXiv:2205.10573_ , 2022.
  * Lee et al. [2023] Jae Yong Lee, SungWoong CHO, and Hyung Ju Hwang.  HyperDeepONet: learning operator with complex target function space using the limited resources via hypernetwork.  In _The Eleventh International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=OAw6V3ZAhSd.
  * Goswami et al. [2022a] Somdatta Goswami, Aniruddha Bora, Yue Yu, and George Em Karniadakis.  Physics-informed deep neural operators networks.  _arXiv preprint arXiv:2207.05748_ , 2022a.
  * Ciarlet [2002] Philippe G. Ciarlet.  _The Finite Element Method for Elliptic Problems_.  Society for Industrial and Applied Mathematics, 2002.  doi: 10.1137/1.9780898719208.  URL https://epubs.siam.org/doi/abs/10.1137/1.9780898719208.
  * Brenner and Scott [2008] Susanne C Brenner and Ridgway Scott.  _The Mathematical Theory of Finite Element Methods_.  Springer Science & Business Media, New York, 2008.
  * Chekroun et al. [2020] Mickaël D. Chekroun, Youngjoon Hong, and Roger Temam.  Enriched numerical scheme for singularly perturbed barotropic quasi-geostrophic equations.  _Journal of Computational Physics_ , 416:109493, 2020.
  * Goswami et al. [2022b] Somdatta Goswami, Minglang Yin, Yue Yu, and George Em Karniadakis.  A physics-informed variational DeepONet for predicting crack path in quasi-brittle materials.  _Computer Methods in Applied Mechanics and Engineering_ , 391:114587, 2022b.  ISSN 0045-7825.  doi: https://doi.org/10.1016/j.cma.2022.114587.  URL https://www.sciencedirect.com/science/article/pii/S004578252200010X.
  * Alber et al. [2019] Mark Alber, Adrian Buganza Tepole, William R Cannon, Suvranu De, Salvador Dura-Bernal, Krishna Garikipati, George Karniadakis, William W Lytton, Paris Perdikaris, Linda Petzold, et al.  Integrating machine learning and multiscale modeling—perspectives, challenges, and opportunities in the biological, biomedical, and behavioral sciences.  _NPJ digital medicine_ , 2(1):115, 2019.
  * Bhatnagar et al. [2019] Saakaar Bhatnagar, Yaser Afshar, Shaowu Pan, Karthik Duraisamy, and Shailendra Kaushik.  Prediction of aerodynamic flow fields using convolutional neural networks.  _Computational Mechanics_ , 64(2):525–545, 2019\.
  * Guo et al. [2016] Xiaoxiao Guo, Wei Li, and Francesco Iorio.  Convolutional neural networks for steady flow approximation.  In _Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining_ , pages 481–490, 2016.
  * Khoo et al. [2017] Yuehaw Khoo, Jianfeng Lu, and Lexing Ying.  Solving parametric PDE problems with artificial neural networks.  _arXiv preprint arXiv:1707.03351_ , 2017.
  * Zhu and Zabaras [2018] Yinhao Zhu and Nicholas Zabaras.  Bayesian deep convolutional encoder–decoder networks for surrogate modeling and uncertainty quantification.  _Journal of Computational Physics_ , 366:415–447, 2018\.
  * Hwang et al. [2022] Rakhoon Hwang, Jae Yong Lee, Jin Young Shin, and Hyung Ju Hwang.  Solving pde-constrained control problems using operator learning.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 36(4):4504–4512, Jun. 2022.  doi: 10.1609/aaai.v36i4.20373.  URL https://ojs.aaai.org/index.php/AAAI/article/view/20373.
  * Kovachki et al. [2021] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Neural operator: Learning maps between function spaces.  _arXiv preprint arXiv:2108.08481_ , 2021.
  * Pfaff et al. [2021] Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez, and Peter Battaglia.  Learning mesh-based simulation with graph networks.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=roNqYL0_XP.
  * Li et al. [2020] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar.  Multipole graph neural operator for parametric partial differential equations.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 6755–6766. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper_files/paper/2020/file/4b21cf96d4cf612f239a6c322b10c8fe-Paper.pdf.
  * Boussif et al. [2022] Oussama Boussif, Yoshua Bengio, Loubna Benabbou, and Dan Assouline.  MAgnet: Mesh Agnostic Neural PDE Solver.  In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, _Advances in Neural Information Processing Systems_ , 2022.  URL https://openreview.net/forum?id=bx2roi8hca8.
  * Zhu et al. [2019] Yinhao Zhu, Nicholas Zabaras, Phaedon-Stelios Koutsourelakis, and Paris Perdikaris.  Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data.  _Journal of Computational Physics_ , 394:56–81, 2019.
  * Li et al. [2021b] Zongyi Li, Hongkai Zheng, Nikola B. Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar.  Physics-informed neural operator for learning partial differential equations.  _CoRR_ , abs/2111.03794, 2021b.  URL https://arxiv.org/abs/2111.03794.
  * Alnaes et al. [2015] M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells.  The FEniCS project version 1.5.  _Archive of Numerical Software_ , 3, 2015.  doi: 10.11588/ans.2015.100.20553.
  * Logg et al. [2012] A. Logg, K.-A. Mardal, G. N. Wells, et al.  _Automated Solution of Differential Equations by the Finite Element Method_.  Springer, 2012.  doi: 10.1007/978-3-642-23099-8.
  * Cybenko [1989] G. Cybenko.  Approximation by superpositions of a sigmoidal function.  _Math. Control Signals Systems_ , 2(4):303–314, 1989.  ISSN 0932-4194.  doi: 10.1007/BF02551274.  URL https://doi.org/10.1007/BF02551274.
  * Hornik [1991] Kurt Hornik.  Approximation capabilities of multilayer feedforward networks.  _Neural Networks_ , 4(2):251–257, 1991.  ISSN 0893-6080.  doi: https://doi.org/10.1016/0893-6080(91)90009-T.  URL https://www.sciencedirect.com/science/article/pii/089360809190009T.
  * Pinkus [1999] Allan Pinkus.  Approximation theory of the MLP model in neural networks.  In _Acta numerica, 1999_ , volume 8 of _Acta Numer._ , pages 143–195. Cambridge Univ. Press, Cambridge, 1999.  doi: 10.1017/S0962492900002919.  URL https://doi.org/10.1017/S0962492900002919.
  * Kidger and Lyons [2020] Patrick Kidger and Terry Lyons.  Universal Approximation with Deep Narrow Networks.  In Jacob Abernethy and Shivani Agarwal, editors, _Proceedings of Thirty Third Conference on Learning Theory_ , volume 125 of _Proceedings of Machine Learning Research_ , pages 2306–2327. PMLR, 09–12 Jul 2020.  URL https://proceedings.mlr.press/v125/kidger20a.html.
  * Ko et al. [2022] Seungchan Ko, Seok-Bae Yun, and Youngjoon Hong.  Convergence analysis of unsupervised Legendre-Galerkin neural networks for linear second-order elliptic PDEs, 2022.
  * Bartlett and Mendelson [2002] Peter L. Bartlett and Shahar Mendelson.  Rademacher and Gaussian complexities: risk bounds and structural results.  _J. Mach. Learn. Res._ , 3(Spec. Issue Comput. Learn. Theory):463–482, 2002.  ISSN 1532-4435.  doi: 10.1162/153244303321897690.  URL https://doi.org/10.1162/153244303321897690.
  * Gnecco and Sanguineti [2008] Giorgio Gnecco and Marcello Sanguineti.  Approximation error bounds via Rademacher’s complexity.  _Appl. Math. Sci. (Ruse)_ , 2(1-4):153–176, 2008\.  ISSN 1312-885X.
  * Wainwright [2019] Martin J. Wainwright.  _High-dimensional statistics_ , volume 48 of _Cambridge Series in Statistical and Probabilistic Mathematics_.  Cambridge University Press, Cambridge, 2019.  ISBN 978-1-108-49802-9.  doi: 10.1017/9781108627771.  URL https://doi.org/10.1017/9781108627771.  A non-asymptotic viewpoint.
  * Shalev-Shwartz and Ben-David [2014] Shai Shalev-Shwartz and Shai Ben-David.  _Understanding Machine Learning - From Theory to Algorithms._ Cambridge University Press, 2014.  ISBN 978-1-10-705713-5.
  * Sirignano and Spiliopoulos [2018] Justin Sirignano and Konstantinos Spiliopoulos.  DGM: A deep learning algorithm for solving partial differential equations.  _Journal of computational physics_ , 375:1339–1364, 2018\.
  * Yu et al. [2018] Bing Yu et al.  The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems.  _Communications in Mathematics and Statistics_ , 6(1):1–12, 2018.
  * Bar-Sinai et al. [2019] Yohai Bar-Sinai, Stephan Hoyer, Jason Hickey, and Michael P Brenner.  Learning data-driven discretizations for partial differential equations.  _Proceedings of the National Academy of Sciences_ , 116(31):15344–15349, 2019.
  * Hong and Jung [2018] Youngjoon Hong and Chang-Yeol Jung.  Enriched spectral method for stiff convection-dominated equations.  _Journal of Scientific Computing_ , 74(3):1325–1346, 2018.

##  Appendix: Finite Element Operator Network for Solving Parametric PDEs

##  Appendix A Notations

The list of notations used throughout the paper is provided in Table 3.

Table 3: Notations Notation | Meaning
---|---
$u$ | ground-truth solution of the given PDE (target function)
$u_{h}$ | finite element solution of the given PDE
$\widehat{u}_{h}=\widehat{u}_{h,n,M}$ | solution predicted by the FEONet
$h$ | discretization parameter for the FEM
$n$ | number of learnable parameters in the neural network
$M$ | number of input samples
$\alpha$ | true coefficients for the FEM basis
$\widehat{\alpha}$ | predicted coefficient by the FEONet
$\phi_{\rm{cor}}$ | corrector basis for the singular perturbation problem
$B[\cdot,\cdot]$ | bilinear form of the variational form of the given PDE
$\ell(\cdot)$ | linear functional of the variational form of the given PDE
$\delta_{ij}$ | Kronecker delta
$\mathbb{N}$ | set of natural numbers
$\mathbb{R}$ | set of real numbers
$D$ | physical domain
$\Omega$ | parametric domain
$\mathbf{x}$ | (vector-valued) physical variable
$\boldsymbol{\omega}$ | (vector-valued) parametric random variable
${\rm{div}}$ | divergence operator
$\boldsymbol{a}(\mathbf{x})$, $c(\mathbf{x})$ | coefficients of PDEs
$f(\mathbf{x},\boldsymbol{\omega})$ | external force on $\mathbf{x}\in D$ parametrized by $\boldsymbol{\omega}\in\Omega$
$\mathcal{L}$ | population loss of the FEONet
$\mathcal{L}^{M}$ | empirical loss of the FEONet
$\mathcal{R}_{M}(\mathcal{G})$ | Rademacher complexity of the function class $\mathcal{G}$

##  Appendix B Convergence analysis of FEONet

In this section, we shall provide the convergence analysis of the FEONet to
establish the theoretical justification of the proposed numerical scheme. In
order to convey the idea of the proof clearly, we will focus solely on self-
adjoint equations with homogeneous Dirichlet boundary conditions:

$\displaystyle-{\rm{div}}\,(\boldsymbol{a}(\mathbf{x})\nabla
u)+c(\mathbf{x})u$ $\displaystyle=f(\mathbf{x})$
$\displaystyle{\rm{in}}\,\,D,$ (18) $\displaystyle u$ $\displaystyle=0$
$\displaystyle{\rm{on}}\,\,\partial D,$

where the coefficient $\boldsymbol{a}(\mathbf{x})$ is uniformly elliptic and
$c(\mathbf{x})\geq 0$. Note, however, that our analysis can be extended to
more general cases in a straightforward manner.

###  B.1 Mathematical formulation

As described before, the input of neural networks is an external forcing term
$f$, which is parametrized by the random parameter $\boldsymbol{\omega}$ in
the probability space $(\Omega,\mathcal{G},\mathbb{P})$. In this section, we
will regard $f(\mathbf{x};\boldsymbol{\omega})$ as a bivariate function
defined on $D\times\Omega$, and assume that

$f(\mathbf{x};\boldsymbol{\omega})\in
C(\Omega;L^{1}(D)):=\left\\{f:\Omega\rightarrow
L^{1}(D):\sup_{\boldsymbol{\omega}\in\Omega}\int_{D}|f(\mathbf{x};\boldsymbol{\omega})|\,\mathrm{d}\mathbf{x}<\infty\right\\}.$
(19)

For each $\boldsymbol{\omega}\in\Omega$, the external force
$f(\mathbf{x};\boldsymbol{\omega})$ is determined, and the corresponding weak
solution is denoted by $u(\mathbf{x};\boldsymbol{\omega})$ which satisfies the
following weak formulation:

$B[u,v]:={\int_{D}[\boldsymbol{a}(\mathbf{x})\nabla u\cdot\nabla
v+c(\mathbf{x})uv]\,\mathrm{d}\mathbf{x}}=\int_{D}f(\mathbf{x})v\,\mathrm{d}\mathbf{x}=:\ell(v)\quad\forall
v\in H^{1}_{0}(D).$ (20)

For given $h>0$, let $V_{h}\subset H^{1}_{0}(D)$ be a conforming finite
element space generated by the basis functions $\\{\phi_{k}\\}_{k=1}^{N(h)}$
and $u_{h}\in V_{h}$ be a finite element approximation of $u$ satisfying the
Galerkin approximation

$B[u_{h},v_{h}]=\ell(v_{h})\quad\forall v_{h}\in V_{h}.$ (21)

We shall write

$u_{h}(\mathbf{x},\boldsymbol{\omega})=\sum^{N(h)}_{k=1}\alpha^{*}_{k}(\boldsymbol{\omega})\phi_{k}(\mathbf{x}),$
(22)

where $\alpha^{*}$ is the solution of the linear algebraic equations

$A\alpha^{*}=F,$ (23)

with

$A_{ik}=B[\phi_{k},\phi_{i}]\quad{\rm{and}}\quad F_{i}=\ell(\phi_{i}).$ (24)

In other words, $u_{h}$ in (22) is the finite element approximation of the
solution to the equation (18), as described in Section 3.1. It is noteworthy
that instead of the equations (23), $\alpha^{*}$ can also be characterized in
the following way:

$\alpha^{*}=\operatorname*{arg\,min}_{\alpha\in
C(\Omega,\mathbb{R}^{N(h)})}\mathcal{L}(\alpha),$ (25)

where $\mathcal{L}$ is the population loss

$\mathcal{L}(\alpha)=\mathbb{E}_{\boldsymbol{\omega}\sim\mathbb{P}_{\Omega}}\bigg{[}\sum^{N(h)}_{i=1}|B[\widehat{u}(\boldsymbol{\omega}),\phi_{i}]-\ell(\phi_{i};(\boldsymbol{\omega}))|^{2}\bigg{]}=\|A\alpha(\boldsymbol{\omega})-F(\boldsymbol{\omega})\|^{2}_{L^{2}(\Omega)}.$
(26)

Next, let us denote the class of feed-forward neural networks by
$\mathcal{N}_{n}$. Here, the subscript $n$ represents the architecture of
neural networks and we assume that $\mathcal{N}_{n_{2}}$ has more expressive
power than $\mathcal{N}_{n_{1}}$ provided that $n_{1}\leq n_{2}$. For example,
$n$ can be the number of layers with bounded width or the number of neurons
when the number of layers is fixed. The class of neural networks is known to
be a suitable ansatz space for the nonlinear approximation, as supported by
the universal approximation theorem (see, for example, Cybenko [1989], Hornik
[1991], Pinkus [1999], Kidger and Lyons [2020]). For the later analysis, in
this section, we shall assume that for all the neural networks under
consideration, the last layer is activated by a bounded activation function
(e.g. sigmoid, tanh, etc.), so that the resulting neural networks are
uniformly bounded. Then by a simple scaling argument, we can show that the
universal approximation theorem still holds for this new class of neural
networks: see, for example, Theorem 2.2 in Ko et al. [2022].

Now as a neural network approximation of $\alpha^{*}$, we seek for
$\widehat{\alpha}(n):\Omega\rightarrow\mathbb{R}^{N(h)}$ solving the residual
minimization problem

$\widehat{\alpha}(n)=\operatorname*{arg\,min}_{\alpha\in\mathcal{N}_{n}}\mathcal{L}(\alpha),$
(27)

where the loss function is minimized over $\mathcal{N}_{n}$, and we shall
write the corresponding solution by

$u_{h,n}(\mathbf{x};\boldsymbol{\omega})=\sum^{N(h)}_{k=1}\widehat{\alpha}(n)_{k}(\boldsymbol{\omega})\phi_{k}(\mathbf{x}).$
(28)

Note that for the neural network under consideration
$\alpha\in\mathcal{N}_{n}$, the input vector is $\boldsymbol{\omega}\in\Omega$
which determines the external force $f(\mathbf{x};\boldsymbol{\omega})$ and
the output is the coefficient vector in $\mathbb{R}^{N(h)}$.

As a final step, let us define the solution to the following discrete residual
minimization problem:

$\widehat{\alpha}(n,M)=\operatorname*{arg\,min}_{\alpha\in\mathcal{N}_{n}}\mathcal{L}^{M}(\alpha),$
(29)

where $\mathcal{L}^{M}$ is the empirical loss defined by the Monte–Carlo
integration of $\mathcal{L}(\alpha)$:

$\mathcal{L}^{M}(\alpha)=\frac{|\Omega|}{M}\sum_{m=1}^{M}\sum^{N(h)}_{i=1}|B[\widehat{u}(\boldsymbol{\omega}_{m}),\phi_{i}]-\ell(\phi_{i};(\boldsymbol{\omega}_{m}))|^{2}=\frac{|\Omega|}{M}\sum^{M}_{m=1}|A\alpha(\boldsymbol{\omega}_{m})-F(\boldsymbol{\omega}_{m})|^{2},$
(30)

with $\\{\boldsymbol{\omega}_{n}\\}^{M}_{m=1}$ is an i.i.d. sequence of random
variables distributed according to $\mathbb{P}_{\Omega}$. Then we write the
corresponding solution as

$u_{h,n,M}(\mathbf{x};\boldsymbol{\omega})=\sum^{N(h)}_{k=1}\widehat{\alpha}(n,M)_{k}(\boldsymbol{\omega})\phi_{k}(\mathbf{x}),$
(31)

which is the numerical solution actually computed by the scheme proposed in
the present paper.

In this paper, we assume that it is possible to find the exact minimizers for
the problems (27) and (29), and the error that occurred by the optimization
can be ignored.

In order to provide appropriate theoretical backgrounds for the proposed
method, it would be reasonable to show that our solution is sufficiently close
to the approximate solution computed by the proposed scheme for various
external forces, as the index $n\in\mathbb{N}$ for a neural-network
architecture and the number of input samples $M\in\mathbb{N}$ goes to
infinity. Mathematically, it can be formulated as

$\|u-u_{h,n,M}\|_{L^{2}(\Omega;L^{2}(D))}\rightarrow
0\quad{\rm{as}}\,\,h\rightarrow 0\,\,{\rm{and}}\,\,n,M\rightarrow\infty.$ (32)

The main error can be split into three parts:

$u-u_{h,n,M}=(u-u_{h})+(u_{h}-u_{h,n})+(u_{h,n}-u_{h,n,M}).$ (33)

The first error is the one caused by the finite element approximation, which
is assumed to be negligible for a proper choice of $h$. In fact, according to
the estimate (7), we can reduce the error as much as we want if we choose
proper $h>0$. Henceforth, we shall assume that we have chosen suitable $h>0$
which guarantees a sufficiently small error for the finite element
approximation. The second error is referred to as the approximation error, as
it appears when we approximate the target function with a class of neural
networks. The third error is often called the generalization error, which
measures how well our approximation predicts solutions for unseen data. We
will focus on proving that our approximate solution converges to the finite
element solution (which is sufficiently close to the true solution) as a
neural network architecture grows and the number of input samples goes to
infinity.

###  B.2 Approximation error

We first note from (26) and (30), that $A$ defined in (22), (23) mostly
determines the structures of the loss functions and it would be advantageous
for us to analyze the loss function if we know more about the matrix. Note
that the matrix $A$ is determined by the structure of the given differential
equations, the choice of basis functions and the boundary conditions.
Therefore, the characterization of $A$ which is useful for the analysis of the
loss function and can cover a wide range of PDE settings simultaneously is
important. The next lemma addresses this issue, which is quoted from Ko et al.
[2022].

######  Lemma B.1.

Suppose that the matrix $A\in\mathbb{R}^{N(h)\times N(h)}$ is symmetric and
invertible, and let $\rho_{\min}=\min_{i}\\{|\lambda_{i}|\\}$,
$\rho_{\max}=\max_{i}\\{|\lambda_{i}|\\}$ where $\\{\lambda_{i}\\}$ is the
family of eigenvalues of the matrix $A$. Then we have for all
$\mathbf{x}\in\mathbb{R}^{N(h)}$,

$\rho_{\min}|\mathbf{x}|\leq|A\mathbf{x}|\leq\rho_{\max}|\mathbf{x}|.$ (34)

Since the equation under consideration (18) is self-adjoint, the corresponding
bilinear form $B[\cdot,\cdot]$ defined in (20) is symmetric which
automatically guarantees that the matrix $A$ in our case is symmetric.
Furthermore, due to the facts that the coefficient $\boldsymbol{a}(\cdot)$ is
uniformly elliptic and $c(\cdot)$ is non-negative, the bilinear form
$B[\cdot,\cdot]$ is coercive, which implies that $A$ is positive-definite.
Therefore, we can apply Lemma B.1 to the matrix $A$ of our interest.

We are now ready to prove the approximation error for neural networks
converges to zero, which is encapsulated in the following theorem.

######  Theorem B.2.

Assume that (19) holds. Then we have that

$\|\alpha^{*}-\widehat{\alpha}(n)\|_{L^{2}(\Omega)}\rightarrow
0\quad{\rm{as}}\,\,n\rightarrow\infty.$ (35)

###### Proof.

Since $A$ is symmetric and invertible, by Proposition B.1, we obtain

$\displaystyle\|\alpha^{*}-\widehat{\alpha}(n)\|^{2}_{2}$
$\displaystyle\leq\frac{1}{(\rho_{\rm{min}})^{2}}\|A\alpha^{*}-A\widehat{\alpha}(n)\|^{2}_{2}\leq\frac{2}{(\rho_{\rm{min}})^{2}}\left(\|A\alpha^{*}-F\|^{2}_{2}+\|A\widehat{\alpha}(n)-F\|^{2}_{2}\right)$
$\displaystyle=\frac{2}{(\rho_{\rm{min}})^{2}}\mathcal{L}(\widehat{\alpha}(n))\leq\frac{2}{(\rho_{\rm{min}})^{2}}\inf_{\alpha\in\mathcal{N}_{n}}\mathcal{L}(\alpha)=\frac{2}{(\rho_{\rm{min}})^{2}}\inf_{\alpha\in\mathcal{N}_{n}}\|A\alpha-F\|^{2}_{2}$
$\displaystyle\leq\frac{4}{(\rho_{\rm{min}})^{2}}\inf_{\alpha\in\mathcal{N}_{n}}\left(\|A\alpha-A\alpha^{*}\|^{2}_{2}+\|A\alpha^{*}-F\|^{2}_{2}\right)$
$\displaystyle\leq
4\left(\frac{\rho_{\rm{max}}}{\rho_{\rm{min}}}\right)^{2}\inf_{\alpha\in\mathcal{N}_{n}}\|\alpha-\alpha^{*}\|^{2}_{2}.$

Note that the term $\rho_{\rm{max}}/\rho_{\rm{min}}$ may depend on $h>0$. But
as we mentioned before, we assumed that $h>0$ is fixed with the sufficiently
small finite element error. Then by the universal approximation property ,
$\inf_{\alpha\in\mathcal{N}_{n}}\|\alpha-\alpha^{*}\|^{2}_{2}\rightarrow 0$ as
$n\rightarrow\infty$, which completes the proof. ∎

###  B.3 Generalization error

We first introduce the concept so-called Rademacher complexity Bartlett and
Mendelson [2002].

######  Definition B.3.

For a collection $\\{X_{i}\\}_{i=1}^{M}$ of i.i.d. random variables, we define
the Rademacher complexity of the function class $\mathcal{G}$ by

$R_{M}(\mathcal{G})=\mathbb{E}_{\\{X_{i},\varepsilon_{i}\\}^{M}_{i=1}}\bigg{[}\sup_{f\in\mathcal{G}}\bigg{|}\frac{1}{M}\sum^{M}_{i=1}\varepsilon_{i}f(X_{i})\bigg{|}\bigg{]},$

where $\varepsilon_{i}$’s denote i.i.d. Bernoulli random variables, in other
words,
$\mathbb{P}(\varepsilon_{i}=1)=\mathbb{P}(\varepsilon_{i}=-1)=\frac{1}{2}$ for
all $i=1,\cdots,M$.

Note that the Rademacher complexity is the expectation of the maximum
correlation between the random noise
$(\varepsilon_{1},\cdots,\varepsilon_{M})$ and the vector
$(f(X_{1}),\cdots,f(X_{M}))$, where the supremum is taken over the class of
functions $\mathcal{G}$. By intuition, the Rademacher complexity of
$\mathcal{G}$ is often used to measure how the functions from $\mathcal{G}$
can fit random noise. For a more detailed discussion on the Rademacher
complexity, see Gnecco and Sanguineti [2008], Wainwright [2019], Shalev-
Shwartz and Ben-David [2014].

Next, we shall find the connection between the generalization error and the
Rademacher complexity for the uniformly bounded function class $\mathcal{G}$.
In the following theorem, we assume that the function class is $b$-uniformly
bounded, meaning that $\|f\|_{\infty}\leq b$ for arbitrary $f\in\mathcal{G}$.

######  Theorem B.4.

[Theorem 4.10 in Wainwright [2019]] Assume that the function class
$\mathcal{G}$ is $b$-uniformly bounded and let $M\in\mathbb{N}$. Then for
arbitrary small $\delta>0$, we have

$\sup_{f\in\mathcal{G}}\bigg{|}\frac{1}{M}\sum^{M}_{i=1}f(X_{i})-\mathbb{E}[f(X)]\bigg{|}\leq
2R_{M}(\mathcal{G})+\delta,$

with probability at least $1-\exp(-\frac{M\delta^{2}}{2b^{2}})$.

Now we define the function class of interest

$\mathcal{G}_{n}:=\\{|A\alpha-F|^{2}:\alpha\in\mathcal{N}_{n}\\},$ (36)

where the matrix $A$ and the vector $F$ are defined in (24). Then by Lemma
B.1, we have that

$\|A\alpha-F\|_{L^{\infty}(\Omega)}\leq\|A\alpha\|_{L^{\infty}(\Omega)}+\|F\|_{L^{\infty}(\Omega)}\leq\rho_{\max}\|\alpha\|_{L^{\infty}(\Omega)}+\|f\|_{C(\Omega;L^{1}(D))}.$

Since the class of neural networks under consideration is uniformly bounded
and (19) holds, we see that for any $n\in\mathbb{N}$, $\mathcal{G}_{n}$ is
$\tilde{b}$-uniformly bounded for some $\tilde{b}>0$. The following lemma is
the direct consequence of Theorem B.4 within our setting.

######  Lemma B.5.

Let $\\{\boldsymbol{\omega}_{m}\\}^{M}_{m=1}$ be i.i.d. random samples
selected from the distribution $\mathbb{P}_{\Omega}$. Then for arbitrary small
$\delta>0$, we obtain with probability at least
$1-2\exp(-\frac{M\delta^{2}}{32\tilde{b}^{2}})$ that

$\sup_{\alpha\in\mathcal{N}_{n}}\left|\mathcal{L}^{M}(\alpha)-\mathcal{L}(\alpha)\right|\leq
2R_{n}(\mathcal{G}_{n})+\frac{\delta}{2}.$ (37)

From Lemma B.5, now we shall establish the following convergence result for
the generalization error. Note that we assume here that for any
$n\in\mathbb{N}$, the Rademacher complexity of $\mathcal{G}_{n}$ converges to
zero, which is known to be true in many cases.

######  Theorem B.6.

Suppose that (19) holds, and we further assume that for any $n\in\mathbb{N}$,
$\lim_{M\rightarrow\infty}R_{M}(\mathcal{G}_{n})=0$. Then we have with
probability $1$ that

$\lim_{n\rightarrow\infty}\lim_{M\rightarrow\infty}\|\widehat{\alpha}(n,M)-\widehat{\alpha}(n)\|_{L^{2}(\Omega)}=0.$

###### Proof.

By Proposition B.1 and the definition (27), we obtain

$\displaystyle\|\widehat{\alpha}(n)-\widehat{\alpha}(n,M)\|^{2}_{2}$
$\displaystyle\leq\frac{1}{(\rho_{\rm{min}})^{2}}\|A\widehat{\alpha}(n)-A\widehat{\alpha}(n,M)\|^{2}_{2}$
(38)
$\displaystyle\leq\frac{1}{(\rho_{\rm{min}})^{2}}\left(\|A\widehat{\alpha}(n)-F\|^{2}_{2}+\|A\widehat{\alpha}(n,M)-F\|^{2}_{2}\right)$
$\displaystyle=\frac{1}{(\rho_{\rm{min}})^{2}}\left(\mathcal{L}(\widehat{\alpha}(n))+\mathcal{L}(\widehat{\alpha}(n,M))\right)$
$\displaystyle\leq\frac{2}{(\rho_{\rm{min}})^{2}}\mathcal{L}(\widehat{\alpha}(n,M)).$

We next apply Lemma B.5 for $\delta=2M^{-\frac{1}{2}+\varepsilon}$ with
$0<\varepsilon<\frac{1}{2}$. Then with probability at least
$1-2\exp(-\frac{M^{2\varepsilon}}{8\tilde{b}^{2}})$, we have from the
minimality of $\widehat{\alpha}(n,M)$ that,

$\mathcal{L}(\widehat{\alpha}(n,M))\leq\mathcal{L}^{M}(\widehat{\alpha}(n,M))+2R_{M}(\mathcal{G}_{n})+M^{-\frac{1}{2}+\varepsilon}\leq\mathcal{L}^{M}(\widehat{\alpha}(n))+2R_{M}(\mathcal{G}_{n})+M^{-\frac{1}{2}+\varepsilon}.$

Using Lemma B.5 again gives us that

$\mathcal{L}(\widehat{\alpha}(n,M))\leq\mathcal{L}(\widehat{\alpha}(n))+4R_{M}(\mathcal{G}_{n})+2M^{-\frac{1}{2}+\varepsilon}.$

Letting $M\rightarrow\infty$ on (38), we obtain that

$\lim_{M\rightarrow\infty}\|\widehat{\alpha}(n,M)-\widehat{\alpha}(n)\|^{2}_{2}\leq\frac{1}{(\rho_{\rm{min}})^{2}}\mathcal{L}(\widehat{\alpha}(n)).$

As we did before, we conclude that

$\displaystyle\lim_{n\rightarrow\infty}\lim_{M\rightarrow\infty}\|\widehat{\alpha}(n,M)-\widehat{\alpha}(n)\|^{2}_{2}$
$\displaystyle\leq\lim_{n\rightarrow\infty}\frac{1}{(\rho_{\rm{min}})^{2}}\mathcal{L}(\widehat{\alpha}(n))$
$\displaystyle\leq\lim_{n\rightarrow\infty}2\left(\frac{\rho_{\rm{max}}}{\rho_{\rm{min}}}\right)^{2}\inf_{\alpha\in\mathcal{N}_{n}}\|\alpha-\alpha^{*}\|^{2}_{2}=0,$

which completes the proof. ∎

###  B.4 Proof of main theorem

From Theorem B.2 and Theorem B.6 combined with the triangular inequality, we
have

$\lim_{n\rightarrow\infty}\lim_{M\rightarrow\infty}\|\alpha^{*}-\widehat{\alpha}(n,M)\|_{L^{2}(\Omega)}=0.$
(39)

Let us now prove the main theorem on the convergence.

######  Theorem B.7 (Main theorem).

Suppose (19) holds, and assume that for all $n\in\mathbb{N}$,
$R_{M}(\widetilde{{\mathcal{G}}}_{n})\rightarrow 0$ as $M\rightarrow\infty$,
with
$\widetilde{{\mathcal{G}}}_{n}:=\\{|A\alpha-F|^{2}:\alpha\in\mathcal{N}_{n}\\}$.
Then for given $h>0$, we have with probability $1$ that

$\lim_{n\rightarrow\infty}\lim_{M\rightarrow\infty}\|u_{h}-u_{h,n,M}\|_{L^{2}(\Omega;L^{2}(D))}=0.$
(40)

###### Proof.

From Theorem B.2, Theorem B.6, we note that for a fixed $h>0$,

$\displaystyle\|u_{h}-u_{h,n,M}\|^{2}_{L^{2}(\Omega;L^{2}(D))}$
$\displaystyle=\int_{\Omega}\int_{D}\bigg{|}\sum^{N(h)}_{i=1}(\alpha^{*}_{i}-\widehat{\alpha}(n,M)_{i})\phi_{i}\bigg{|}^{2}\,\mathrm{d}\mathbf{x}\,\mathrm{d}{\boldsymbol{\omega}}$
$\displaystyle=\int_{\Omega}\int_{D}\bigg{|}\sum^{N(h)}_{i,j=1}(\alpha^{*}_{i}-\widehat{\alpha}(n,M)_{i})(\alpha^{*}_{j}-\widehat{\alpha}(n,M)_{j})\phi_{i}\phi_{j}\bigg{|}^{2}\,\mathrm{d}\mathbf{x}\,\mathrm{d}{\boldsymbol{\omega}}$
$\displaystyle\leq\int_{\Omega}\int_{D}\bigg{(}\sum^{N(h)}_{i,j=1}|\alpha^{*}_{i}-\widehat{\alpha}(n,M)_{i}|^{2}|\phi_{i}|^{2}+\sum^{N(h)}_{i,j=1}|\alpha^{*}_{j}-\widehat{\alpha}(n,M)_{j}|^{2}|\phi_{j}|^{2}\bigg{)}\,\mathrm{d}\mathbf{x}\,\mathrm{d}{\boldsymbol{\omega}}$
$\displaystyle\leq\int_{\Omega}\int_{D}2N(h)\sum^{N(h)}_{k=1}|\alpha^{*}_{k}-\widehat{\alpha}(n,M)_{k}|^{2}|\phi_{k}|^{2}\,\mathrm{d}\mathbf{x}\,\mathrm{d}{\boldsymbol{\omega}}$
$\displaystyle\leq
2|D|N(h)\|\alpha^{*}-\widehat{\alpha}(n,M)\|^{2}_{L^{2}(\Omega)}\rightarrow
0,$

where we have used the fact that $|\phi_{j}|\leq 1$ for all $1\leq j\leq N(h)$
together with Young’s inequality. ∎

##  Appendix C Comparison of FEONet with other models

The proposed FEONet model is an unsupervised-learning approach that does not
require labeled training data, and it has a number of advantages in learning
the solution operator as confirmed in Section 4. Unlike methods designed for a
single instance, such as PINN Raissi et al. [2019], this novel numerical
scheme enables the handling of multiple instances of parametric partial
differential equations. On the other hand, well-established models for
operator learning, such as FNO Li et al. [2021a] and DeepONet Lu et al.
[2021b], depend on labeled training data for their training process. To
address this limitation of operator-learning methods, recent advances in
physics-informed methods, such as PIDeepONet and PINO, have emerged. These
methods utilize the PDE residual loss, enabling unsupervised training of the
models. However, these methods are not suitable for complex 2D domains or
singular perturbation problems as demonstrated in the experiments conducted in
Section 4.2.

Table 4: Comparison of various machine learning approaches for solving PDEs.

| Unsupervised | PDE instances | Prediction location |  | Applicable for complex domains
---
and singularly perturbed problem
PINN Raissi et al. [2019], DGM Sirignano and Spiliopoulos [2018], DRM Yu et al. [2018] | No | Single instance | Entire domain | $\triangle$
DON Lu et al. [2021b] | No | Multiple instances | Entire domain | $\triangle$
FNO Li et al. [2021a] | No | Multiple instances | Grid points | $\times$
PIDON Wang et al. [2021] | Yes | Multiple instances | Entire domain | $\triangle$
PINO Li et al. [2021b] | Yes | Multiple instances | Grid points | $\times$
| Graph-based neural PDE solver
---
Pfaff et al. [2021], Li et al. [2020], Brandstetter et al. [2022], Lienen and
Günnemann [2022], Boussif et al. [2022].
No | Multiple instances | Grid points | -
Ours | Yes | Multiple instances | Entire domain | $\Circle$

The comparison with the FNO and PINO models was omitted as they are not
applicable to problems on domains with complex geometry, such as those
depicted in Figure 1, due to their restriction on uniform grids for input and
output.

The FEONet effectively utilizes the weak form of the given PDE by employing it
as a loss function over the triangulation of the physical domain. In contrast,
PIDeepONet requires a complex grid sampling process in a complicated 2D
domain, resulting in computational costs associated with calculating the
residual loss using automatic differentiation. Moreover, as shown in Table 2
and Figure 5, when a small diffusion coefficient $\varepsilon>0$ is used in a
PDE residual loss term, the models struggle to accurately learn the
corresponding solutions of the PDEs. This raises a significant concern
regarding reliability, as illustrated in Figure 11, which displays the
predictions from various training trials. In this regard, we conducted a
comparative analysis of the supervised DeepONet approach for PDE problems
under various conditions, considering different amounts of available data.
Additional comparisons of various deep learning-based approaches for solving
PDEs are demonstrated in Table 4.

##  Appendix D Experiment details

###  D.1 Hyperparameters

For the problems under consideration, we used the neural network, which
consists of 6 convolutional layers with swish activation followed by a fully
connected layer flattening the output. For 1D problems, we used Conv1D, while
Conv2D was used for 2D problems. The FEONet was trained with the LBFGS
optimizer along with the following hyperparameters.

  * •

Maximal number of iterations per optimization step : $10$,

  * •

Termination tolerance on first-order optimality : $10^{-15}$,

  * •

Termination tolerance on function value/parameter changes : $10^{-15}$,

  * •

Update history size : $10$.

For both the DeepONet and the PIDeepONet, a fully connected neural network
with a depth of 3 and a width of 100 was used for the trunk and branch
networks. The two models differ only in the loss function, where the training
data and the residual of the PDE were used. The PINN also used a fully
connected neural network with a depth of 3 and a width of 100. All these
models used the tanh activation function and a learning rate of $10^{-4}$ with
the commonly used ADAM optimizer. To ensure the sufficient convergence of the
training results for the boundary layer problem, we conducted five experiments
with the PINN for $5\times 10^{5}$ iterations and five experiments with the
PIDeepONet for $10^{4}$ iterations. Note that computing the residual loss of
PDE for every input function $f$ in the PIDeepONet requires significant
computation time using the automatic differentiation. We used the Intel Xeon
Gold 6226R processor and NVIDIA RTX A6000 48GB GPU.

Table 5: Parameters for each problem

| Input function | PDE dimension | P1 or P2 | # Element |  | Iteration
---
(FEONet)
| Iteration
---
(DeepONet)
Domain I |  $m_{j}\in[1,2]$ $n_{j}\in[0,\pi]$ | 2D | P2 | 392 | $5\times 10^{4}$ | $5\times 10^{5}$
Domain II | 551
Domain III | 334
BC I |  $m_{j}\in[3,5]$ $n_{j}\in[0,2\pi]$ | 1D | P2 | 24 | $1.5\times 10^{5}$ | $5\times 10^{5}$
BC II | 32 | $2\times 10^{5}$ | $5\times 10^{5}$
Eq I |  $m_{j}\in[3,5]$ $n_{j}\in[0,2\pi]$ | 1D | P2 | 32 | $1.5\times 10^{5}$ | $5\times 10^{5}$
Eq II | P1 | 128 | $3.5\times 10^{5}$ | $5\times 10^{5}$
Singular |  | $m_{j}\in[3,5]$
---
$n_{j}\in[0,2\pi]$
1D | P1 | 32 | $5\times 10^{4}$ | $5\times 10^{5}$

###  D.2 Trade-off of FEONet

The order of the polynomial basis and the number of elements employed in
FEONet for each domain, boundary condition, and equation are presented in the
third and fourth columns of Table 5. Utilizing a higher-order basis and a
greater number of elements for domain triangulation in FEONet can enhance
prediction accuracy, but it comes at the cost of increased memory usage and
longer training time. Therefore, determining the appropriate basis and element
number for FEONet involves a trade-off between prediction accuracy and
computational cost. Further details can be found in Table 5.

###  D.3 FEniCS for ground-truth solution, domain triangulation and nodal
basis

The FEniCS [Alnaes et al., 2015, Logg et al., 2012], which stands for Finite
Element Computational Software, is a widely used open-source software suite
that is specifically designed for solving PDEs using the FEM. It offers a
robust and user-friendly interface that enables solving a broad range of PDEs
in various settings. Using the FEniCS, we obtained solutions with very small
errors for each PDE problem on sufficiently fine meshes, which were used as
the ground-truth solutions for each problem. We were able to easily obtain
nodes and nodal basis by triangulating the 2D domain as shown in Figure 1 and
Figure 8. Additionally, in our experiments, we were able to compute the loss
function (9) by using the FEniCS both for the P1 and the P2 elements.

Figure 8: FEM triangulation of Domain II (square with a hole) with varying
numbers of elements.

###  D.4 Types and random generation of input functions for FEONet

This paper specifically focuses on the learning of the solution operator of a
partial differential equation (PDE) from an external forcing term. Notably,
the FEONet can be readily extended to handle various input functions, such as
boundary conditions, coefficients, or initial conditions, in order to learn
the solution operator of a given PDE. As an example, the FEONet is capable of
learning the operator mapping from the variable coefficient $c(x)$ to the
corresponding solution $u(x)$ of a PDE (14) with variable coefficients. This
demonstrates the versatility and adaptability of the FEONet approach. Figure 9
displays the solution profile obtained when variable coefficients are
considered as input. This experiment demonstrates the feasibility of the
FEONet approach for accommodating different types of input functions.

Figure 9: Solution profiles predicted by the trained FEONet for PDE (14)
obtained from variable coefficients $c(x)$ as two distinct input functions.

In order to train the network, we generate random external forcing functions
or variable coefficient functions. Inspired by Bar-Sinai et al. [2019], we
created a random signal $f(\mathbf{x};\boldsymbol{\omega})$ as a linear
combination of sine functions and cosine functions. More precisely, we use

$f(x)=m_{0}\sin(n_{0}x)+m_{1}\cos(n_{1}x)$ (41)

for a 1D case and

$f(x,y)=m_{0}\sin(n_{0}x+n_{1}y)+m_{1}\cos(n_{2}x+n_{3}y)$ (42)

for a 2D case where $m_{i}$ for $i=1,2$ and $n_{j}$ for $j=0,1,2,3$ are drawn
independently from the uniform distributions whose ranges are shown in the
second column of Table 5. It is worth noting that even when considering
different random input functions, such as those generated by Gaussian random
fields, we consistently observe similar results. This robustness indicates the
reliability and stability of the FEONet approach across various input
scenarios.

###  D.5 Description of various domains with complex geometries

Figure 1 illustrates three complex domains considered in this paper. Firstly,
Domain I is chosen as a domain with a smooth boundary, represented by a circle
centered at $(0,0)$ with radius 1, to accommodate non-rectangular shapes.
Secondly, Domain II is designed to have a challenging geometry with holes,
excluding a circle of radius 0.5 centered at $(0,0)$ from the square domain
$[-1,1]\times[-1,1]$. Lastly, Domain III is created to have an irregular
polygon with sharp turns, connecting the vertices
$(-1.0,-1.0),(1.0,-1.0),(1.0,0.0),(0.6,0.0),(0.5,-0.5),(-0.0,-0.5),(-0.6,1.0),(-1.0,1.0)$,
and $(-1.0,-1.0)$. Although this paper focuses on these three domains in the
experiments, the proposed FEONet has the potential to solve PDEs using FEM
triangulation for more complex domains such as airfoils where the classical
FEM works.

###  D.6 Burgers’ equation

For the variational loss function for the Burgers’ equation described in (15),
we identified the $uu_{x}$ term as $(\frac{1}{2}u^{2})_{x}$. Therefore, we
induce the $B[\widehat{u}_{h}(x;\boldsymbol{\omega}_{m}),\phi_{i}(x)]$ in the
loss function (9) as

$\displaystyle B[\widehat{u}_{h}(x;\boldsymbol{\omega}_{m}),\phi_{i}(x)]$
$\displaystyle=\varepsilon\int_{D}(\widehat{u}_{h})_{x}(\phi_{i})_{x}\,\mathrm{d}x+\int_{D}\widehat{u}_{h}(\widehat{u}_{h})_{x}\phi_{i}\,\mathrm{d}x$
$\displaystyle=\varepsilon\int_{D}(\widehat{u}_{h})_{x}(\phi_{i})_{x}\,\mathrm{d}x-\int_{D}\frac{1}{2}\widehat{u}_{h}^{2}(\phi_{i})_{x}\,\mathrm{d}x$
$\displaystyle=\underbrace{\varepsilon\sum_{k=1}^{N(h)}\widehat{\alpha}_{k}\int_{D}(\phi_{k})_{x}(\phi_{i})_{x}\,\mathrm{d}x}_{\text{(I)}}-\underbrace{\int_{D}\frac{1}{2}\left(\sum_{k=1}^{N(h)}\widehat{\alpha}_{k}\phi_{k}\right)^{2}(\phi_{i})_{x}\,\mathrm{d}x}_{\text{(II)}}.$

The main difference from other PDE problems is the second term (II) on the
right-hand side. Since the term (II) cannot be expressed linearly with respect
to $\widehat{\alpha}_{k}$, an iterative method should be used for the
linearization when using the classical FEM. However, the loss function can be
directly computed for the FEONet, and once the training is completed, it can
make a real-time solution prediction without any need for iteration methods
whenever the forcing function is given. In other words, the FEONet is able to
learn the coefficients $\\{\widehat{\alpha}_{k}\\}_{k=1}^{N(h)}$ effectively
even for the nonlinear Burgers’ equation.

###  D.7 Singular perturbation problem

Figure 10: Boundary layer transition for (16) with respect to different
viscosity parameters $\varepsilon>0$ is displayed. As $\varepsilon$ decreases,
the thickness of the boundary layer gets smaller.

Solving convection-dominated singularly perturbed problems in numerical
analysis is a challenging task as the small diffusive parameter
$\varepsilon>0$ generates a sharp transition inside thin layers, requiring
special treatment. The transition of the boundary layer for (16) with respect
to different viscosity parameters $\varepsilon>0$ is provided in Figure 10. ML
approaches suffer from the boundary layer problem since neural networks have a
smooth prior (see Figure 5). When solving a single singularly perturbed PDE
using the PINN, the training using residual loss becomes unstable due to a
very small diffusion coefficient, resulting in poor learning. Figure 11 shows
different solutions predicted by five different PINN training for a single PDE
in $5\times 10^{5}$ iterations.

Our FEONet overcomes this limitation by using an additional basis
$\phi_{\text{cor}}$, as defined in (17). We obtain the corrector
$\phi_{\text{cor}}$ by deriving the asymptotic equation with the formal limit
equation assuming $\varepsilon=0$, as described in Hong and Jung [2018].
Figure 12 shows the shape of the boundary layer corrector for
$\varepsilon=10^{-5}$.

As a paradigm example to derive $\phi_{\text{cor}}$, we consider the following
singularly perturbed convection-diffusion equation:

$\left\\{\begin{array}[]{rl}-\varepsilon
u^{\varepsilon}_{xx}-u^{\varepsilon}_{x}=f(x),&0<x<1,\\\
u^{\varepsilon}=0,&x=0,1.\end{array}\right.$ (43)

Our main objective is to construct the corrector basis for the singular
perturbation problem. To find the corrector, we formally replace $\varepsilon$
by $0$ in (43)1, whose limit problem of (43) at $\varepsilon=0$ is

$\left\\{\begin{array}[]{rl}-u^{0}_{x}=f(x),&0<x<1,\\\
u^{0}=0,&x=1.\end{array}\right.$ (44)

We impose the inflow boundary condition at $x=1$ for $u^{0}$, leading to the
formal limit of $u^{0}$ in the following form:

$u^{0}(x)=\int_{0}^{1}f\,\mathrm{d}x.$ (45)

Performing matching asymptotics for equation (43), we observe the presence of
a boundary layer of size $\epsilon$ near the outflow boundary at $x=0$.
Consequently, we determine the asymptotic equation, considering small
$\epsilon$, for the corrector $\varphi$ that approximates the difference
$u^{\epsilon}-u^{0}$. This equation can be written as follows:

$\left\\{\begin{array}[]{rl}-\varepsilon\varphi_{xx}-\varphi_{x}=0,&0<x<1,\\\
\varphi=-u^{0},&x=0.\end{array}\right.$ (46)

It is well-known that the corrector $\varphi$ is given in the form

$\varphi(x)=-u^{0}(0)\,e^{-x/\varepsilon}+{\rm{e.s.t.}},$ (47)

where the e.s.t. stands for an exponentially small term with respect to the
small perturbation parameter $\varepsilon>0$. To obtain the implementable
corrector basis, we modify equation (47) by incorporating the boundary
condition.

Performing the conventional energy estimates on the difference
$u^{\varepsilon}-(u^{0}+\varphi)$, we notice that

$\begin{array}[]{l}\|u^{\varepsilon}-(u^{0}+\varphi)\|_{L^{2}((0,1))}\leq\kappa\varepsilon,\\\
\|u^{\varepsilon}-u^{0}\|_{L^{2}((0,1))}\leq\kappa\varepsilon^{\frac{1}{2}},\end{array}$
(48)

for a constant $\kappa>0$ independent of $\varepsilon$. The convergence
results indicate two important observations. The convergence results highlight
two significant observations. Firstly, the diffusive solution $u^{\epsilon}$
converges to the limit solution $u^{0}$ as $\epsilon$ approaches 0, with the
convergence rate proportional to $\epsilon$. Secondly, based on these
convergence results, we can deduce that the corrector $\varphi$ effectively
captures the singular behavior of $u^{\epsilon}$ at small diffusivity
$\epsilon$. This implies that the diffusive solution $u^{\epsilon}$ can be
decomposed into the sum of a rapidly decaying component represented by
$\varphi$ and a slowly varying component represented by $u^{0}$.

Figure 11: The solution profiles of five different experiments solving the
singular perturbation problem using the PINN. Figure 12: The profile of the
boundary layer corrector which is used as the additional basis function for
the singular perturbation problem.

###  D.8 Plots on convergence rates

Let us make some further comments on the experimental details for Figure 6.
For 1D case, we consider the convection-diffusion equation

$\displaystyle-\varepsilon u_{xx}+bu_{x}$ $\displaystyle=f(x),$ $\displaystyle
x\in D=[-1,1],$ (49) $\displaystyle u(x)$ $\displaystyle=0\quad,$
$\displaystyle x\in\partial D,$

where $\varepsilon=0.1$ and $b=-1$. For the 2D case, we consider the following
equation

$\displaystyle-\varepsilon\Delta u+\textbf{v}\cdot\nabla u$
$\displaystyle=f(x,y),$ $\displaystyle(x,y)\in D=[-1,1]^{2},$ (50)
$\displaystyle u(x,y)$ $\displaystyle=0,$ $\displaystyle(x,y)\in\partial D,$

where $\varepsilon=0.1$ and $\textbf{v}=(-1,0)$. The FEONet was used to
approximate the solution of these two equations using P1 and P2 nodal basis
functions, and the experiments were conducted with varying domain
triangulation to have different numbers of elements to observe the convergence
rate. As shown in Figure 6, the observed convergence rate of the FEONet shows
the convergence rates close to 2 and 3 for P1 and P2 approximation
respectively, which are the theoretical results for the classical FEM. Since
the precision scale of the neural network we used here is about $10^{-2}$, we
cannot observe the exact trends in rates below this level of error. We expect
to see the same trend even at lower errors if we use larger scale and more
advanced models than the CNN structure.

##  Appendix E Limitations

Despite the promising results presented in this study, there are still several
technical questions that remain unresolved and necessitate further
investigation. One such question is determining the optimal features embedding
or network architecture of the FEONet for a given parametric governing law.
Understanding the most effective approach for feature representation and
network design is crucial for achieving optimal performance and generalization
in the FEONet framework. Addressing this question would contribute to
advancing the understanding and applicability of the FEONet methodology.
Furthermore, this approach holds the potential for addressing challenging
problems, including the Navier-Stokes equations, which represent one of the
most fascinating nonlinear problems in the numerical solution of partial
differential equations. By applying the FEONet to such complex scenarios as a
future direction, we can make significant advancements in the field and
uncover new insights into the behavior and dynamics of these systems.