Augmented physics informed extreme learning machine to solve the biharmonic equations via Fourier expansions

This section provides a detailed introduction to the PIELM method for solving the biharmonic equation (1) with boundaries (2) or (3). To obtain an ideal solution to the biharmonic equation, one can seek the solution in trial function space spanned by the ELM model with a high-order smooth activation function. Then, an ansatz expressed by ELM is designed to approximate the solutions of (1), it is:

Substituting $u(\bm{x};\bm{W},\bm{b})$ into the equation of (1) for $u$, we can obtain the following equation

with $V_{i}(\bm{x})=\bm{x}^{T}\cdot\bm{W}^{[i]}+b_{i}$ being the linear output for $i_{th}$ hidden unit of ELM and $N$ is the number of hidden unit. Additionally, $u(\bm{x};\bm{W},\bm{b})$ still holds on the prescribed boundaries (2) or (3), we then have the following formulation

or

For solving the biharmonic equations by the aforementioned PIELM method, a collection of collocation points is sampled from both the interior and boundaries of the domain. At these collocation points, the governed function (5) and the (6) or (7) exhibit coerciveness. Let $Q$ denote the number of random collocation points in the interior of $\Omega$ and denote by $\mathcal{S}_{I}=\{\bm{x}_{I}^{q}\}_{q=1}^{Q}$, and $P$ denote the number of random collocation points on hyperface of $\partial\Omega$ and denote by $\mathcal{S}_{B}=\{\bm{x}_{B}^{p}\}_{p=1}^{P}$. Feeding the above collocation points into the pipeline of PIELM, we have the following residual terms for Dirichlet boundaries:

and the last term in (8) will be replaced by

for Navier boundaries.

To approximate the functions $u$ well, we need to adjust the output weights of the PIELM model to make each residual of (8) to near-zero or minimize the error of the following system of linear equations:

Here $\displaystyle\mathbf{H}_{f}=\mathbf{H}_{B}+\sum_{i=1}^{d}\sum_{j=1,j\neq i}^{d}\mathbf{H}_{ij}$ with

and

Additionally, the matrices $\mathbf{H}_{g}$, $\mathbf{H}_{h}$ and $\mathbf{H}_{k}$ are given by

$\displaystyle\mathbf{H}_{h}=\mathbf{H}_{1}\cdot n_{1}+\mathbf{H}_{2}\cdot n_{2}+\cdots+\mathbf{H}_{d}\cdot n_{d}$ with

and

In the above expressions, $\sigma^{\prime\prime}$ and $\sigma^{\prime\prime\prime\prime}$ are the second and fourth-order derivatives of activation function $\sigma$, respectively, and $\bigodot$ stands for the Hadamard product operator. In addition, $\mathbf{S}_{f}=\big{(}f(\bm{x}_{I}^{1}),f(\bm{x}_{I}^{2}),\cdots,f(\bm{x}_{I}^{Q})\big{)}^{T}$, $\mathbf{S}_{g}=\big{(}g(\bm{x}_{B}^{1}),g(\bm{x}_{B}^{2}),\cdots,g(\bm{x}_{B}^{P})\big{)}^{T}$, $\mathbf{S}_{h}=\big{(}h(\bm{x}_{B}^{1}),h(\bm{x}_{B}^{2}),\cdots,h(\bm{x}_{B}^{P})\big{)}^{T}$ and $\mathbf{S}_{k}=\big{(}k(\bm{x}_{B}^{1}),k(\bm{x}_{B}^{2}),\cdots,k(\bm{x}_{B}^{P})\big{)}^{T}$.

When the weights and biases of the hidden layer for the PIELM model are randomly initialized and fixed, the matrix $\mathbf{H}$ is determined for given input data, then the output weights $\bm{\beta}$ of PIELM model are obtained by solving the linear systems (10). The system equation $\mathbf{H}\bm{\beta}=\mathbf{S}$ is solvable and meets the following several cases:

• 1.

  If the coefficient matrix $\mathbf{H}$ is square and invertible, then

  $$\bm{\beta}=\mathbf{H}^{-1}\mathbf{S}.$$

  (11)

• 2.

  If the coefficient matrix $\mathbf{H}$ is rectangular and full rank, then

  $$\bm{\beta}=\mathbf{H}^{\dagger}\mathbf{S},$$

  (12)

  where $\mathbf{H}^{\dagger}$ is the pseudo inverse of $\mathbf{H}$ defned by $\mathbf{H}^{\dagger}=(\mathbf{H}^{T}\mathbf{H})^{-1}\mathbf{H}^{T}$.

• 3.

  If the coefficient matrix $\mathbf{H}$ is singular, then the Tikhonov regularization[51] can be utilized to to overcome the ill-posedness of (10), and the $\bm{\beta}$ is given by

  $$\bm{\beta}=(\mathbf{H}^{T}\mathbf{H}+\lambda\mathbf{I})^{-1}\mathbf{H}\mathbf{S},$$

  (13)

  where $\mathbf{I}$ is the identity matrix, and $\lambda>0$ is the ridge parameter.

According to the above discussions, the detailed procedures involved in the PIELM algorithm to address the biharmonic equation (1) within finite-dimensional spaces are described in the following.

In terms of function approximation, the hidden nodes governed by a given activation function in the ELM model can be considered a set of basis functions, and then the ELM model produces its output through a linear aggregation of these fundamental functions. Similar to other artificial neural network (ANN) architectures, the choice of activation function is critical for the ELM model’s performance. Since the biharmonic equations (1) studied here is a high-order PDE, activation functions with strong regularity are preferred. Previous studies, including Fabiani et al. [40], Calabrò et al. [41], and Ma et al. [52], have investigated the influence of sigmoidal functions and radial basis functions in the ELM model for solving both linear and nonlinear PDEs. These studies indicate that ELM models with sigmoid, Gaussian, and tangent activation functions achieve commendable numerical accuracy.

The sigmoid, gaussian, and tangent functions all fulfill the Lipschitz continuity and exhibit robust regularity, they ensure the boundedness of the output for hidden units in ELM. In terms of the high-order PDEs, the high-order derivate of activation is required. Then, the first, second, and fourth-order derivatives for the sigmoidal function $s(x)=1/(1+e^{-x})$, Gaussian function $\sigma(x)=e^{-x^{2}}$ and hyperbolic tangent function $tanh(x)=(e^{x}-e^{-x})/(e^{x}+e^{-x})$ are given by

and

respectively. Figs. 2 – 4 illustrate the curves of the sigmoid, Gaussian, and tangent functions and their corresponding first-order, second-order and fourth-order derivatives, respectively.

This visualization shows that the derivatives of the sigmoid, Gaussian, and tanh functions are bounded and tend to saturate for narrow-range inputs. However, their performance may degrade with large-range inputs. Inspired by the Fourier-based spectral method, a global basis function can generate responses, along with its derivatives, for both narrow- and large-range inputs. Furthermore, approximations based on such basis functions can handle data across various scales. The sine function serves as an ideal global Fourier basis function (with sine and cosine being fundamentally equivalent). In this study, we adopt the sine function as the activation function in the PIELM method, referring to this modified approach as FPIELM.

In this section, we apply the proposed FPIELM method to solve the two- and three-dimensional biharmonic equation (1) in both unitized and non-unitized domains. Additionally, PIELM models using Sigmoid, Gaussian, and Tanh as activation functions are established as baselines for solving (1) and are referred to as SPIELM, GPIELM, and TPIELM, respectively. Following the results in Dong and Yang [46], the weight matrices and bias vectors for the hidden layers of each model are initialized uniformly within the interval $[-\delta,\delta]$.

A criterion is provided to evaluate the performance of the four aforementioned methods, given by

where $\tilde{u}(\bm{x}^{i})$ and $u^{*}(\bm{x}^{i})$ denote the approximate solution obtained by numerical methods and the exact solution, respectively, for the collection set $\{\bm{x}^{i}\}_{i=1}^{N^{\prime}}$, with $N^{\prime}$ being the number of collection points.

All experiments are implemented in Numpy (version 1.21.0) on a workstation equipped with 256 GB RAM and a single NVIDIA GeForce RTX 4090 Ti GPU with 24 GB of memory.

First, we approximate the solution of (1) using the four methods on 16,384 equidistant grid points distributed across regular domains: $[-1,1]\times[-1,1]$, $[0,5]\times[0,5]$, $[-4,6]\times[-3,7]$, and $[5,15]\times[0,10]$. In addition to the grid points, the training set includes 10,000 random points within $\Omega$ and 4,000 random points on $\partial\Omega$. The number of hidden nodes for each of the four PIELM methods is set to 1,000.

Figure 5 shows the point-wise absolute error of the four PIELM models on the domain $[-1,1]\times[-1,1]$. Table 1 presents the optimal value of $\delta$, the REL, and the runtime for each PIELM method across the different domains.

Next, we examine the impact of the scale factor $\delta$ and the number of hidden nodes on the FPIELM model in the domain $[-1,1]\times[-1,1]$, with results illustrated in Fig. 6. Finally, comparison curves for REL as a function of the scale factor $\delta$ for the four PIELM methods on the domains $[-1,1]\times[-1,1]$ and $[0,5]\times[0,5]$ are shown in Fig. 7.

Based on the point-wise absolute error distribution shown in Figures 5(b) – 5(e) and the data in Table 1, it can be concluded that the FPIELM method outperforms other types of PIELM methods. As observed in the curves in Figure 6(a), the relative error does not decrease with an increasing number of hidden units, indicating that the performance of the FPIELM model stabilizes as the number of hidden neurons grows. The time-neuron curve in Figure 6(b) demonstrates a superlinear increase in runtime with the addition of hidden neurons, though this increase is gradual. Furthermore, Figure 6(c) illustrates that the FPIELM model remains stable and robust across variations in the scale factor $\sigma$ on the domain $[-1,1]\times[-1,1]$. The REL curves in Figure 7 confirm that FPIELM maintains stability and allows for a broader range of scale factor $\delta$ values across both unitized and non-unitized domains. In summary, our FPIELM method not only achieves high accuracy but also demonstrates efficiency and robustness. Additionally, it is noteworthy that the coefficient matrix in the final linear system is full rank, given that the number of samples exceeds the number of hidden nodes.

For comparison, we employ the finite difference method (FDM) with a five-point central difference scheme to solve the coupled scheme of (1) on the domain $[0,5]\times[0,5]$. The FDM is implemented in MATLAB (version 12.0.0.341360, R2019a). For the FPIELM model, the number of hidden nodes is set to 1,000, and the scale factor $\delta$ for initializing weights and biases is set to 5.0. The numerical results in Table 2 indicate that the proposed FPIELM method maintains stability across various mesh grid sizes and significantly outperforms FDM in terms of accuracy, particularly on sparse grids.

In this example, the setup for the four PIELM methods is identical to that in Example 4.1. An approximate solution to (1) is derived by applying the four methods on 19,214 collection points located within the hexagram domains $\Omega$. Numerical results are presented in Table 3 for the two scenarios, and the point-wise absolute error for the four PIELM methods is shown in Figure 8 over the hexagram domain derived from $[-\pi,\pi]\times[-\pi,\pi]$.

Additionally, Figure 9 illustrates the variation in relative error and runtime concerning the number of hidden nodes, as well as the variation in relative error with the scale factor $\delta$ for FPIELM on the hexagram domain derived from $[-\pi,\pi]\times[-\pi,\pi]$. Finally, Figure 10 plots the REL comparison curves for the four PIELM methods, showing how REL varies with the scale factor $\delta$ across the hexagram domains derived from $[-\pi,\pi]\times[-\pi,\pi]$ and $[0,3\pi]\times[\pi,2\pi]$.

Based on the data in Table 3 and the point-wise absolute errors depicted in Figs. 8(b) – 8(e), we can infer that the FPIELM method outperforms the PIELM methods with sigmoid, tanh, and Gaussian activation functions. As observed in the curves in Fig. 16(a), the performance of the FPIELM method stabilizes as the number of hidden nodes increases. The time-neuron curve in Fig. 16(b) shows that the runtime exhibits a superlinear increase with the growing number of hidden nodes, although the increase remains gradual. Lastly, Figs. 9(c) and 10 illustrate that the FPIELM method remains stable and robust across variations in the scale factor $\sigma$. In summary, our FPIELM method demonstrates not only high accuracy but also efficiency and robustness.

In this example, the number of hidden units for each of the four PIELM methods is set to 2,000. The testing set consists of 40,000 random points within the cut planes, excluding the holes, and includes planes parallel to the $xoy$ and $yoz$ coordinate planes. Additionally, 6,000 points randomly sampled from the six boundary surfaces are integrated into the testing set. The results are presented in Figs. 12 – 14.

Based on the point-wise error shown in Figs. 12(b) – 12(e), our proposed FPIELM method outperforms the other three methods in solving the biharmonic equation (1) in three-dimensional space. The relative errors presented in Fig. 14 indicate that the FPIELM method remains stable across various scale factors $\delta$ and is superior to the other methods. Furthermore, the curves of relative error and runtime in Fig. 13 demonstrate the favorable stability and robustness of the FPIELM method as both the number of hidden nodes and the scale factor $\delta$ increase. It is also noteworthy that the coefficient matrix in the final linear system is full rank for this example, eliminating the need for adjustment through Tikhonov regularization.

In this example, the setup for the four PIELM methods is identical to that in Example 4.1. An approximate solution of (1) is obtained by applying the four methods on 16,384 grid points within the rectangular domains $[0,1]\times[0,1]$ and $[0,4]\times[0,4]$. The numerical results are listed in Table 4 for these two scenarios, and the point-wise absolute error for the four PIELM methods is shown in Fig. 15 on the domain $[0,1]\times[0,1]$.

Additionally, Fig. 16 illustrates the variation in relative error and runtime concerning the number of hidden nodes, as well as the variation in relative error with the scale factor $\delta$ for FPIELM on the domain $[0,1]\times[0,1]$. Finally, Fig. 17 plots the comparison curves for REL as a function of the scale factor $\delta$ for the four PIELM methods on the rectangular domains $[0,1]\times[0,1]$ and $[0,4]\times[0,4]$, respectively.