Finite Element Operator Network for Solving Parametric PDEs

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

for some constant $C>0$. For each $\bm{\omega}\in\Omega$, the external force $f(\bm{x};\bm{\omega})$ is determined, and the corresponding weak solution is denoted by $u(\bm{x};\bm{\omega})$ which satisfies the following weak formulation:

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

We shall write

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

with

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

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

Next, we shall define a class of neural networks, that plays the role of an approximator for coefficients. For each $L\in\mathbb{N}$, we write an $L$-layer neural network by a function $f^{L}(x):\mathbb{R}^{n_{0}}\rightarrow\mathbb{R}^{n_{L}}$, defined by

where $W^{\ell}\in\mathbb{R}^{n_{\ell}\times n_{\ell-1}}$ and $b^{\ell}\in\mathbb{R}^{n_{\ell}}$ are the weight and the bias respectively for the $\ell$-th layer, and $\sigma:\mathbb{R}\rightarrow\mathbb{R}$ is a nonlinear activation function. We denote the architecture by the vector ${\vec{\bm{n}}}=(n_{0},\cdots,n_{L})$, the class of neural network parameters by $\theta:=\theta_{\vec{\bm{n}}}=\{(W^{1},b^{1}),\cdots,(W^{L},b^{L})\}$, and its realization by $\mathcal{R}[\theta](x)$. For a given neural network architecture $\vec{\bm{n}}$, the family of all possible parameters is denoted by

For two neural networks $\theta_{i}$, $i=1,2$, with architectures $\vec{\bm{n}}_{i}=\left(n^{(i)}_{0},\cdots,n^{(i)}_{L_{i}}\right)$, we write $\vec{\bm{n}}_{1}\subset\vec{\bm{n}}_{2}$ if for any $\theta_{1}\in\Theta_{\vec{\bm{n}}_{1}}$, there exists $\theta_{2}\in\Theta_{\vec{\bm{n}}_{2}}$ satisfying $\mathcal{R}[\theta_{1}](x)=\mathcal{R}[\theta_{2}](x)$ for all $x\in\mathbb{R}^{n_{0}}$.

Now let us assume that there exists a sequence of neural network architectures $\{\vec{\bm{n}}_{n}\}_{n\geq 1}$ satisfying $\vec{\bm{n}}_{n}\subset\vec{\bm{n}}_{n+1}$ for all $n\in\mathbb{N}$. We define the associated collection of neural networks by

It is straightforward to verify that $\mathcal{N}^{*}_{n}\subset\mathcal{N}^{*}_{n+1}$ for any $n\in\mathbb{N}$. We start with the following theorem. This is known as the universal approximation theorem, and known to hold in various scenarios [13, 20, 42].

In the present paper, the target function to be approximated by neural networks is $\alpha^{*}(\bm{\omega})$. Due to the assumption (10), we can verify that $\alpha^{*}(\bm{\omega})$ is continuous and bounded. This motivates us to deal with the family of bounded neural networks as an ansatz space. For this, let us write $\sup_{\bm{\omega}\in\Omega}|\alpha^{*}(\bm{\omega})|<C_{\alpha}$, and assume that $\sigma^{*}$ is a bounded activation function with continuous inverse (e.g., sigmoid or tanh activation functions) with $\sigma_{-}:=\inf_{x\in\mathbb{R}}\sigma^{*}(x)$ and $\sigma_{+}:=\sup_{x\in\mathbb{R}}\sigma^{*}(x)$. We then define a function $h:[\sigma_{-},\sigma_{+}]\rightarrow[-C_{\alpha},C_{\alpha}]$ satisfying $h(\sigma_{-})=-C_{\alpha}$, $h(\sigma_{+})=C_{\alpha}$. By the assumptions, $h$ also has a continuous inverse. With the above notations, we consider the following family of neural networks:

where we can see that

In the analysis, we will use the following version of the universal approximation theorem for bounded neural networks, which is quoted from [23].

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

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

Note that for the neural network under consideration $\alpha\in\mathcal{N}_{n}$, the input vector is $\bm{\omega}\in\Omega$ which determines the external force $f(\bm{x};\bm{\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:

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

with $\{\bm{\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

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 (25) and (27), 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

The main error can be split into three parts:

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 classical theory of the FEM (e.g. [8, 11]) if $h>0$ is sufficiently small, then we can expect to obtain a sufficiently small error. Therefore, in this paper, we shall assume that we have chosen suitable $h>0$ which guarantees a sufficiently small error for the finite element approximation and only focus on the convergence of $u_{h}-u_{h,n,M}$ with a fixed $h>0$.

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, i.e., we will show that

We first note from (17) and (28), that $A$ defined in (13), (14) 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 [23].

Since the equation under consideration (9) is self-adjoint, the corresponding bilinear form $B[\cdot,\cdot]$ defined in (11) is symmetric which automatically guarantees that the matrix $A$ in our case is symmetric. Furthermore, due to the facts that the coefficient $\bm{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 2.4 to the matrix $A$ of our interest.

We are now ready to present our first convergence result regarding the approximation error. The proof is based on Lemma 2.4 and can be found in Appendix A.1.

In order to handle the generalization error, we first introduce the concept so-called Rademacher complexity [5].

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 [18, 48, 45].

Now we define the function class of interest

where the matrix $A$ and the vector $F$ are defined in (15). In the following theorem, note that we assume 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 [45, 14, 37]. For the proof of the theorem below, see Appendix A.2.

From Theorem 2.5 and Theorem 2.7 combined with the triangular inequality, we have probability 1 over i.i.d. samples that

Based on the convergence of the approximate coefficients (35), let us now state the convergence theorem for the predicted solution, whose proof is presented in Appendix A.3.

We consider the 2D convection-diffusion equation as

For our numerical experiments, we fixed the values of $\varepsilon=0.1$ and $\bm{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 relative $L^{2}$ errors of the test set for FEONet (w/o labeled data), PIDeepONet (w/o labeled data), and DeepONet (with 30/300/3000 input-output data pairs) across various domains. When the DeepONet is trained with either 30 or 300 data pairs, FEONet consistently surpasses them in terms of numerical errors. Notably, even when DeepONet is trained with 3000 data pairs, FEONet achieves errors that are either comparable to or lower than those of DeepONet. It’s noteworthy that in more intricate domains like Domain II and Domain III, both the DeepONet (even with 3000 data pairs) and PIDeepONet find it challenging to produce accurate predictions. Figure 3 demonstrates the solution profile, emphasizing that FEONet is more precise in predicting the solution than both DeepONet and PIDeepONet, particularly in complex domains.

We tested the performance of each model when the given PDEs were equipped with either Dirichlet or Neumann boundary conditions. We consider the convection-diffusion equation as

where $\varepsilon=0.1$ and $b=-1$. For the Neumann boundary condition, we add $u$ to the left-hand side of (38) to address the uniqueness of the solution. 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.

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

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

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 any training data. Once the training process is completed, the model can provide accurate real-time solution predictions without relying on iterative schemes. 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.

The FEONet is capable of learning the operator mapping $\mathcal{G}:c(x)\mapsto u(x)$ from the variable coefficient $c(x)$ to the corresponding solution $u(x)$ of a PDE (39). Figure 4 displays the solution profile obtained when variable coefficients are considered as input. This highlights the flexibility of FEONet in adapting to a wide range of input function types, thereby underscoring its versatility.

We address the case when the input is a boundary condition. We consider the 2D Poisson equation as

where $f(x,y)=2$ and the domain $D$ is a square with a hole (see (b) in Figure 1). The FEONet is utilized to learn the operator $\mathcal{G}:g(x,y)\mapsto u(x,y)$. Figure 5 demonstrates that even in the case of an operator that takes the boundary condition $g(x,y)$ as input, we can predict the solution with high accuracy. This demonstrates the versatility and adaptability of the FEONet approach for accommodating different types of input functions.

The FEONet can be extended to handle the initial conditions as input. We consider the time-dependent problem with varying initial conditions. We consider the 1D convection-diffusion equation

where $\nu=0.1$ and $b=-1$. We aim to learn the operator $\mathcal{G}:u_{0}(x)\mapsto u(t,x)$. We use the implicit Euler method to discretize the time with $\Delta t=0.01$. We employed the marching-in-time scheme proposed by [24] to sequentially learn and predict solutions over the time domain $t=[0,1]$, divided into 10 intervals. For training FEONet, we utilized only the input function $u_{0}(x)$ and equation (42) without any additional data. Figure 6 illustrates FEONet’s effective prediction of true solutions using an initial condition from the test dataset.

We consider a more intricate yet practical fluid dynamics issue: the two-dimensional steady-state Stokes equations, which can be described as follows:

where $\bm{g}(x,y)=(3+1.7\sin(2\pi x),0)$. We aim to learn the operator $\mathcal{G}:\bm{f}=(f_{1},f_{2})\mapsto[\bm{u},p]$ using the FEONet. For this, we employ the P2 Lagrange element for velocity and the P1 Lagrange element for pressure. Together, this combination is known as the Taylor-Hood element [15]. Figure 7 demonstrates the results of the operator learning for the 2D Stokes equation using FEONet. This illustrates that FEONet is applicable to the learning of operators for the system of PDEs with vectorial inputs and outputs.

The FEONet is a model capable of learning the operator without the need for paired input-output training data. What we need for the training is the generation of random input samples, which only imposes negligible computational cost. Therefore, the number of input function samples for the FEONet training can be chosen arbitrarily. In the experiments performed in Section 3, we generated 3000 random input function samples to train the FEONet. On the other hand, the DeepONet requires paired input-output training data for supervised learning which causes significant computational costs to prepare. In this paper, we have focused our experiments on considering the DeepONet trained with 30, 300, and 3000 input-output data pairs. This allowed us to highlight the FEONet’s ability to achieve a certain level of accuracy without relying on data pairs.

We present the results comparing generalization capability in Figure 9, when the FEONet is trained only with 30 random input samples (green and yellow bar) (which should be distinguished from 30 input-out data pair for the supervised learning). Figure 9 depicts a comparison of train and test errors between the FEONet, trained through an unsupervised approach with 30 and 3000 input samples, and the DeepONet, trained through supervised learning with 30 input-output data pairs. The disparity between the training error and the test error serves as a reliable indicator of how effectively the models generalize the underlying phenomenon. Notably, the FEONet shows superior generalization compared to the DeepONet, even without requiring input-output paired data. Remarkably, utilizing only 30 input samples for training, the FEONet achieves significantly better generalization and exhibits lower test errors than that of the DeepONet.

One of the intriguing aspects is that the theoretical results on the convergence error rate of the FEONet are also observable in the experimental results. Figure 10 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.