A Novel Neural Network-Based Approach Comparable to High-Precision Finite Difference Methods

Abstract

Deep learning methods using neural networks for solving partial differential equations (PDEs) have emerged as a new paradigm. However, many of these methods approximate solutions by optimizing loss functions, often encountering convergence issues and accuracy limitations. In this paper, we propose a novel deep learning approach that leverages the expressive power of neural networks to generate basis functions. These basis functions are then used to create trial solutions, which are optimized using the least-squares method to solve for coefficients in a system of linear equations. This method integrates the strengths of streaming PINNs and the traditional least-squares method, offering both flexibility and a high accuracy. We conducted numerical experiments to compare our method with the results of high-order finite difference schemes and several commonly used neural network methods (PINNs, lbPINNs, ELMs, and PIELMs). Thanks to the mesh-less feature of the neural network, it is particularly effective for complex geometries. The numerical results demonstrate that our method significantly enhances the accuracy of deep learning in solving PDEs, achieving error levels comparable to high-accuracy finite difference methods.

1. Introduction

Partial differential equations (PDEs) play a central role in the natural sciences, engineering, mathematics, and related fields, providing an effective tool for describing and analyzing complex phenomena and systems. These equations not only help us understand natural laws, but they also lay the foundation for scientific research and engineering practices [1,2]. Since many PDEs do not have analytical solutions, traditional numerical methods, such as the finite element method (FEM), finite difference method (FDM), finite volume method (FVM), etc., are widely used to compute approximate solutions, assisting researchers and engineers in obtaining meaningful results for various problems.

High-order compact FDMs are widely used for solving PDEs due to their high accuracy and effectiveness. These schemes typically achieve a fourth-order or higher spatial accuracy with smaller truncation errors compared to traditional methods. They require fewer mesh points to achieve the same accuracy, which is particularly beneficial for high-dimensional problems and helps reduce computational costs [3,4,5,6,7]. However, the construction and implementation of compact difference formats are complex, especially when dealing with intricate boundary conditions. They often assume that the solution is smooth, which may not be valid for problems with sharp gradients or discontinuities, potentially leading to unstable or inaccurate results. Additionally, many compact methods are implicit, requiring the solution of larger systems of linear equations to be solved, which increases the implementation complexity and computation time, particularly for large-scale problems.

In recent years, deep learning has gained traction in solving PDEs, leading to the development of various innovative methods, such as physically informed neural networks (PINNs) [8], the deep Ritz method (DRM) [9], deep Galerkin methods (DGMs) [10], and extreme learning machines (ELMs) [11], Self-adaptive loss balanced Physics-informed neural networks (lbPINN) [12]. These approaches demonstrate significant advantages in addressing complex PDE problems. PINNs integrate physical laws into the neural network’s loss function, enabling the model to learn from both data and the governing PDE during training. This method leverages the neural network’s strong fitting capabilities to represent the PDE solution, minimizing the error between the predicted and actual solutions. PINNs are particularly effective for problems with complex boundary conditions and multi-scale features, often yielding good results even with limited training data [13,14]. DGMs combine deep learning with Galerkin techniques to approximate PDE solutions through neural networks that satisfy the weak form of the equations. ELMs are single hidden-layer feedforward neural networks known for their simple and rapid training process, ideal for large-scale data. Unlike traditional neural networks, ELMs randomly select hidden layer weights and determine output weights using least-squares, significantly reducing the training time. They can quickly adapt to varying boundary and initial conditions when solving PDEs.

In addition, the extrapolation—driven network architecture—proposed by Wang, Y. [15] optimizes the physics—informed deep—learning model by combining extrapolation techniques, further enhancing the model’s accuracy and generalization ability. Zeng, B. [16] developed a method of integrating physical encoding into the message—the passing graph network (PhyMPGN), which is specifically designed to solve the spatiotemporal PDE system. These methods demonstrate the great potential of deep learning in the field of numerically solving PDEs, especially when traditional numerical methods encounter bottlenecks.

These deep learning methods offer fresh perspectives and tools for tackling PDEs, especially in scenarios where traditional numerical methods struggle. They effectively handle complex boundary conditions and nonlinear characteristics while minimizing mesh dependency. However, these approaches often require substantial data and computational resources, raising concerns about their feasibility and efficiency in practical applications [17,18,19,20]. Moreover, most deep learning-based methods necessitate optimizing the loss function to approximate the solutions accurately. The complexity of the networks and the number of hyperparameters can lead to optimization challenges, potentially resulting in local optima and an insufficient accuracy. To address these issues, researchers need to employ refined strategies in model design and training, such as advanced optimization algorithms, network structure adjustments, and regularization techniques, to enhance the generalization and accuracy [21,22,23,24].

Compared to the optimization loss function-based approach of PINNs, the ELM-based basis function method effectively avoids convergence issues in the optimization process and demonstrates a significant improvement in accuracy. This strategy reduces the reliance on traditional optimization algorithms by utilizing randomly generated weights and biases, resulting in faster computations and an enhanced model stability. Additionally, the ELM-based model exhibits a superior adaptability in managing complex boundary and initial conditions, making it more advantageous for addressing practical problems. Dwivedi [25] developed physics-informed extreme learning machines (PIELMs), which are fast versions of PINNs that can be applied to static and time-varying linear partial differential equations to meet or exceed the accuracy of PINNs on a range of problems. Quan [26] proposed a novel learning method based on the extreme learning machine algorithm, which utilizes the ELM algorithm to solve a system of linear equations, thereby determining the network parameters for solving PDEs. Wang et al. [27] presented an efficient extension of the extreme learning machine (ELM) method, moving from low dimensions to solving high-dimensional partial differential equations (PDEs). A neural network-based solution for linear and nonlinear partial differential equations has been proposed by integrating concepts from extreme learning machines (ELMs), domain decomposition, and local neural networks [28]. A high-precision subspace method based on neural networks has been proposed for solving partial differential equations. The core idea involves using neural network-based functions as the basis functions to span a subspace, allowing for the approximation of solutions within that subspace [29]. While neural networks hold significant promise for numerically solving PDEs, challenges like model generalization remain. Integrating neural networks with traditional numerical methods offers new avenues for achieving high-precision solutions, likely resulting in more efficient and accurate computations across various fields.

In light of the findings in the literature, we propose a novel deep learning strategy that integrates the ELM with the basis function space method. This approach leverages the advantages of the ELM to enhance the solution process and improve the accuracy of solving PDEs. Firstly, the ELM can be seen as selecting specific basis functions in the function space through randomly generated weights and biases. This not only enriches the expressiveness of the function space, but also significantly boosts the problem-solving efficiency, allowing for the better capture of the PDE solution characteristics and an improved model accuracy. Secondly, the ELM demonstrates an excellent adaptability, enabling flexible adjustments to various boundary and initial conditions. This adaptability is crucial in basis function space methods, as many real-world problems often involve complex and variable conditions. By incorporating the ELM, our model can maintain a strong performance across different physical scenarios, broadening its applicability. Finally, the ELM’s training time is substantially shorter than that of traditional neural networks, greatly enhancing the computational efficiency when addressing complex PDEs. This rapid training process reduces the computational costs and facilitates real-time or near real-time applications, particularly in engineering and scientific computing contexts that demand quick responses. In summary, our proposed deep learning strategy that combines an ELM with the basis function space method not only enhances the efficiency and accuracy of PDE solving, but also offers new insights and tools for tackling complex physical problems.

This paper is organized as follows: Section 2 introduces the model for partial differential equations and the fundamental principles of PINNs. Section 3 explores a high-precision basis function method based on neural networks. Section 4 evaluates the accuracy and effectiveness of the numerical validation method. The conclusion is presented in Section 5.

2. Mathematical Model and Physics-Informed Neural Network

2.1. Mathematical Model

Consider the following generalized form of PDE, which can be expressed as

(1)$f(x;{\displaystyle \frac{\partial u}{\partial x_1}},\cdots {\displaystyle \frac{\partial u}{\partial x_d}},{\displaystyle \frac{{\partial }^{2}u}{\partial x_1\partial x_1}},\cdots {\displaystyle \frac{{\partial }^{2}u}{\partial x_1\partial x_d}})=0,x\in \Omega$

where $x=\left\{x_1,x_2,\cdots x_d\right\}$ are the independent variables, $u(x)$ is the unknown partial function to be solved, and $u_{x_i}$ is the partial derivative of $u$ with respect to $x_i$. PDEs can be classified into different types, such as elliptic, hyperbolic, parabolic, etc., depending on their nature. The exact form and type depend on the structure and boundary conditions of the equation.

To simplify the discussion, the PDEs and their boundary conditions can be transformed into the operator form.

(2)$\left\{\begin{array}{cc}Lu(x)=f(x),\hfill & \hfill x\in \Omega \\ Bu(x)=g(x),\hfill & \hfill x\in \partial \Omega \end{array}\right.$

where $x={(x_1,x_2,\cdots ,x_d)}^T$, $\Omega$ is a bounded domain in ${ℝ}^d$, $\partial \Omega$ is the $\Omega$ boundary, $L$ is a differential operator describing the physical process, and $B$ is a boundary operator that specifies the conditions required at the boundary. $f$ and $g$ are the given functions.

2.2. Physics-Informed Neural Networks (PINNs)

PINNs integrate the laws of physics into neural network training, combining data-driven and physical models. They are particularly effective for optimizing the learning of complex systems through strong constraints and have recently emerged as a novel extensible approach to solving PDEs. The process of PINNs for solving PDEs involves three key elements, as shown in Figure 1: (1) the neural network construction trial function; (2) the design of the loss function, which includes PDE residuals, boundary conditions, and labeled data; and (3) the optimization of the training process. These elements interact to create a comprehensive framework for solving PDEs using PINNs.

First, a feedforward neural network (FNN) is defined with $N$ hidden layers to construct the trial function. The trial function $u(x,\Theta )$ is used as a hypothetical solution to approximate the solution of a partial differential equation. The accuracy of the solution is optimized by tuning the network parameters $\Theta$ (i.e., all weights and biases) through training.

(3)$u(x,\Theta ):={\mathrm{W}}_N\sigma ({\mathrm{W}}_{N-1}\sigma (\cdots {\mathrm{W}}_2\sigma ({\mathrm{W}}_{1}x+b_1)+b_2\cdots )+b_{N-1})+b_N$

where $x$ is the input data, usually containing independent variables (e.g., location and time); ${\mathrm{W}}_{\mathrm{i}}$ and ${\mathrm{b}}_i$ are the weight matrix and bias vector of i-th layer, with $i=1,2,\cdots ,N$ respectively; $\Theta =\{{\cup }_{k=1}^N({\mathrm{W}}^k,{\mathrm{b}}^k),k=1,2,\cdots ,N\}$ is the set of tunable parameters of the neural network; and $\sigma$ is the activation function.

This trial function (2) constructed from neural networks was used in the PINN framework to solve the corresponding partial differential equations, combining the loss function and the training optimization process to achieve a sufficient accuracy. To constrain the optimization direction of the trial function, PINNs define an objective loss function composed of four components: (1) PDE residual loss, which ensures that the trial function $u(x,\Theta )$ minimizes the error in satisfying the given PDE. (2) Boundary condition loss, which ensures that the trial function meets practical requirements under the boundary conditions, such as Dirichlet or Neumann conditions. (3) Initial condition loss (if applicable), which ensures that the trial function meets the specified initial conditions. (4) Data loss term, which accounts for errors relative to the provided data, if labeled data are available. Together, these components guide the optimization process effectively.

(4)$L(x,u,\Theta ,\omega )={\omega }_f{L}_{PDE}(x,u,\Theta )+{\omega }_i{L}_{IC}(x,u,\Theta )+{\omega }_b{L}_{BC}(x,u,\Theta )+{\omega }_d{L}_{Data}(x,u,\Theta )$

where $L_{PDE}$, $L_{BC}$, $L_{IC}$, and $L_{data}$ are defined as follows:

$\begin{array}{l}{L}_{PDE}(x,u,\Theta )={\displaystyle \frac{1}{\left|M_{\Omega }\right|}}{\displaystyle \sum _{i=1}^{M_{\Omega }}{\left|Lu(x_i,\Theta )-f(x_i)\right|}^2},\begin{array}{cc}& x_i\in \widehat{\Omega }\end{array}\\ L_{BC}(x,u,\Theta )={\displaystyle \frac{1}{\left|M_{\partial \Omega }\right|}}{\displaystyle \sum _{i=1}^{M_{\partial \Omega }}{\left|Bu(x_i,\Theta )-g(x_i)\right|}^2,\begin{array}{cc}& x_i\in \partial \widehat{\Omega }\end{array}}\\ L_{IC}(x,u,\Theta )={\displaystyle \frac{1}{\left|M_I\right|}}{\displaystyle \sum _{i=1}^{M_I}{\left|Iu(x_i,\Theta )-u_0(x_i)\right|}^2,\begin{array}{cc}& x_i\in \widehat{\Omega }\end{array}}\\ L_{Data}(x,u,\Theta )={\displaystyle \frac{1}{\left|M_{data}\right|}}{\displaystyle \sum _{i=1}^{M_{data}}{\left|u(x_i,\Theta )-u_d(x_i)\right|}^2,\begin{array}{cc}& x_i\in Data\end{array}}\end{array}$

and $\omega =\{{\omega }_f,{\omega }_i,{\omega }_b,{\omega }_d\}$ are weighting factors to balance the contribution of each component to the total loss. $M_{\Omega }, M_{\partial \Omega }, M_I$, and $M_{data}$ represent the number of collocations in the interior of the domain, on the boundary, at the initial condition, and the number of labeled data points, respectively.

Finally, the network weight and bias parameters are updated using an optimization algorithm (e.g., Adam, L-BFGS, a combination of Adam and L—BFGS, or stochastic gradient descent) to minimize the overall loss function (3). The training process concludes when the loss reaches a specified threshold, it shows minimal decrease, or the maximum number of iterations is reached.

(5)${\Theta }^{*}=\underset{\Theta }{\arg \min }L(x,u,\Theta ,\omega ), u(x)\approx u(x,{\Theta }^{*})$

Although PINNs offer a versatile framework for solving partial differential equations (PDEs), they face significant challenges in optimizing the loss function to accurately approximate the PDE solutions. This optimization process often requires substantial CPU time and can lead to issues such as convergence to local optima due to the non-convex nature of the loss function. These challenges not only hinder the convergence of the training process, but also compromise the accuracy of the approximations produced by the network.

3. The Neural Network-Based Basis Function Space Method

The universal approximation theorem demonstrates the remarkable expressive power of neural networks, indicating that even a simple two-layer neural network can approximate complex functions to an arbitrary accuracy. Building on this foundation, we utilized the concept of extreme learning machines (ELMs) to generate basis functions through a neural network framework. These basis functions served as the building blocks for constructing our trial solution. Once the trial solution was formulated from the basis functions, we enhanced its accuracy by incorporating traditional numerical methods, such as least-squares optimization, as shown in Figure 2. This hybrid approach allows us to leverage the strengths of both machine learning and numerical analyses, resulting in a more robust and accurate solution to the problem at hand. By effectively combining neural network-based basis functions with established numerical techniques, we can achieve significant improvements in accuracy while maintaining computational efficiency. This methodology not only streamlines the process of finding solutions to complex equations, but also expands the potential applications of neural networks in solving various mathematical and engineering problems. Through this integration, we aimed to harness the full potential of neural networks and traditional numerical methods, paving the way for the more effective and accurate approximation of complex functions.

A feedforward neural network is defined with $N$ hidden layers, where the $j$-th hidden layer contains $N_j$ neurons, the mathematics of which can be expressed by the following equation:

(6)$\left\{\begin{array}{l}{h}_0=x,\hfill \\ h_k={\sigma }_k(W_k\cdot h_{k-1}+b_k),k=1,2,\cdots N-1,\hfill \\ \Phi =h_{N-1},\hfill \\ u(x,\Theta )=\Phi \cdot W_N,\hfill \end{array}\right.$

Given an input $x\in {ℝ}^d$, the output of the network $u(x,\Theta )$ is computed through a series of weights and activation functions.$W_k$ and $b_k$ represent the weight and bias, respectively. $\sigma (\cdot )$ is the activation function, $\theta =\left\{W_1,W_2,\cdots ,W_N,b_1,b_2,\cdots ,b_{N-1}\right\}$ is a set of parameters in the neural network, and $u(x,\theta )$ is the output with parameters $\theta$. The initial set of parameters for a neural network are usually randomly generated.

From Equation (6), it is evident that the neural network’s output can be viewed as a combination function of the last layer’s weight coefficients and the output of the penultimate hidden layer’s neurons. This combination allows the network to capture complex input–output relationships and perform nonlinear mapping. Consequently, this multi-layer structure facilitates layer-wise feature extraction, enhancing the model’s expressive power and learning efficiency.

Assuming that the number of neurons in the penultimate hidden layer is M, then the output of M neurons can be viewed as a set of basis functions and denoted as $\Phi$

(7)$\Phi =\{{\phi }_1(x,\Theta ),{\phi }_2(x,\Theta ),\cdots ,{\phi }_M(x,\Theta )\}$

and the function space V can be written as

(8)$V=span\{{\phi }_1(x,\Theta ),{\phi }_2(x,\Theta ),\cdots ,{\phi }_M(x,\Theta )\}$

The weights of the last hidden layer can be expressed as

(9)$\omega =\{{\omega }_1,{\omega }_2,\cdots ,{\omega }_M\}$

Thus, the output of the neural network can be understood as a function formed by a linear combination of the basis functions, represented as

(10)$u(x,\Theta ,\omega )=\omega \bullet \Phi ={\displaystyle \sum _{k=1}^M{\omega }_k}\phi (x,\Theta )$

Given the diverse types of PDEs and their corresponding solutions, it is essential to train the basis functions appropriately to enable them to express the solutions of PDEs more effectively. This training process focuses on optimizing the parameters of the basis functions, allowing them to more accurately capture the unique characteristics of the solutions for specific PDEs. We used the PINN approach to construct the loss function for pre-training the neural network. Due to the complexity of the PINN optimization process, we established a stopping condition: training was halted when the loss function decreased by two orders of magnitude. This strategy ensured that the network achieved a higher accuracy while minimizing the risk of overfitting, thereby enhancing the model’s generalization ability. Generally, the loss function is as follows:

(11)$L_D={\displaystyle \frac{1}{N_{\Omega }}}{\displaystyle \sum _{j=1}^{N_{\Omega }}(Lu(x_j,\Theta })-f(x_j){)}^2+{\displaystyle \frac{1}{N_{\partial \Omega }}}{\displaystyle \sum _{j=1}^{N{\partial }_{\Omega }}(\mathrm{B}u(x_j,\Theta })-g(x_j){)}^2$

Sufficient training was conducted on the loss function (11) until the maximum number of training steps was reached or the following condition was met:

(12)${\displaystyle \frac{L_D}{L_{D_0}}}\le \epsilon .$

where $L_{D0}$ is the initial loss error. We adopted $\epsilon ={10}^{-3}$ and the maximum number of training steps was $N_{\max }=1000$.

(13)${\Theta }^{*}=\arg \underset{\Theta }{\min }{L}_D$

After completing the pre-training of the network (12), we obtained an initial set of optimal weights that satisfactorily met the requirements of the partial differential control equations and their boundary conditions (13). Using these weights, we constructed customized basis functions tailored for specific boundary value problems involving partial differential equations.

(14)${\Phi }^{*}=\{{\phi }_1(x,{\Theta }^{*}),{\phi }_2(x,{\Theta }^{*}),\cdots ,{\phi }_M(x,{\Theta }^{*})\}$

These basis functions are crucial for their adaptability, as they are specifically designed to effectively address particular problems and accurately capture the solution characteristics. Accordingly, the trial solution output by the neural network can be expressed as

(15)$u(x,\omega )=\omega \bullet {\Phi }^{*}={\displaystyle \sum _{k=1}^M{\omega }_k}\phi (x,{\Theta }^{*})$

A set of neural network-based basis functions was obtained by optimizing the loss function, as shown in Equation (14), and then the approximation solution could be expressed as a linear combination of this set of basis functions, as shown in Equation (15). Next, the combination coefficients were optimized to give a higher approximation accuracy.

Substituting the trial solution in Equation (15) into the governing equation, Equation (1), and making the control equation satisfied at each collocation point yielded the following system of equations:

(16)$\left\{\begin{array}{l}{\displaystyle \sum _{j=1}^M{\omega }_{j}A{\phi }_j(x_i)=f(x_i),\begin{array}{cc}& x_i\in \widehat{\Omega },\end{array}i=1,2,\cdots ,N_{\Omega },}\\ {\displaystyle \sum _{j=1}^M{\omega }_j\mathrm{B}{\phi }_j(x_i)=g(x_i),\begin{array}{cc}& x_i\in \partial \widehat{\Omega },\end{array}i=1,2,\cdots ,N_{\partial \Omega }.}\end{array}\right.$

This is a system of linear equations for the combined weight coefficients $\omega ={[{\omega }_1,{\omega }_2,\cdots ,{\omega }_M]}^{\mathrm{T}}$, which can be written in matrix form as

(17)$\left[\begin{array}{ccccc}A{\phi }_1(x_1)& A{\phi }_2(x_1)& A{\phi }_3(x_1)& \cdots & A{\phi }_M(x_1)\\ A{\phi }_1(x_2)& A{\phi }_2(x_2)& A{\phi }_3(x_2)& \cdots & A{\phi }_M(x_2)\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ A{\phi }_1(x_{N_{\Omega }})& A{\phi }_2(x_{N_{\Omega }})& A{\phi }_3(x_{N_{\Omega }})& \cdots & A{\phi }_M(x_{N_{\Omega }})\\ \mathrm{B}{\phi }_1(x_1)& \mathrm{B}{\phi }_2(x_1)& \mathrm{B}{\phi }_3(x_1)& \cdots & \mathrm{B}{\phi }_M(x_1)\\ \mathrm{B}{\phi }_1(x_2)& \mathrm{B}{\phi }_2(x_2)& \mathrm{B}{\phi }_3(x_3)& \cdots & \mathrm{B}{\phi }_M(x_2)\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ \mathrm{B}{\phi }_1(x_{N_{\partial \Omega }})& \mathrm{B}{\phi }_2(x_{N_{\partial \Omega }})& \mathrm{B}{\phi }_3(x_{N_{\partial \Omega }})& \cdots & \mathrm{B}{\phi }_M(x_{N_{\partial \Omega }})\end{array}\right]\left[\begin{array}{c}{\omega }_1\\ {\omega }_2\\ {\omega }_3\\ {\omega }_4\\ {\omega }_5\\ {\omega }_6\\ ⋮\\ {\omega }_M\end{array}\right]=\left[\begin{array}{c}f(x_1)\\ f(x_2)\\ ⋮\\ f(x_{N_{\Omega }})\\ g(x_1)\\ g(x_2)\\ ⋮\\ g(x_{N_{\partial \Omega }})\end{array}\right]$

It can be further abbreviated as

(18)$\mathrm{H}\omega =\mathrm{R}$

where ${\mathrm{H}}_{(N_{\Omega }+N_{\partial \Omega })\times M}$, $R_{(N_{\Omega }+N_{\partial \Omega })\times 1}$, and ${\omega }_{M\times 1}$.

This is a system of hyper-deterministic linear equations, and we can solve this system using the least-squares method to obtain the optimal combination of coefficients.

(19)${\omega }^{*}{=(\mathrm{H}}^{\mathrm{T}}{\mathrm{H})}^{-1}{\mathrm{H}}^{\mathrm{T}}\mathrm{R}$

The optimized trial solution will approximate the analytical solution of the PDEs with a high accuracy.

(20)$\left|u(x,{\omega }^{*})-u^{*}(x)\right|<\epsilon$

4. Numerical Examples

In this section, we demonstrate the high-precision characteristics of the present method through a series of numerical experiments.

The definition of error was as follows:

${‖e‖}_{L2}=\sqrt{{\displaystyle \frac{{\displaystyle \sum _{i=1}^n{(u_{\theta }(x_i)-{\widehat{u}}_{\theta }(x_i))}^2}}{{\displaystyle \sum _{i=1}^n{u}_{i\theta }(x_i)}}}}, {‖e‖}_{L\infty }=\underset{0\le i\le n}{\max }\left|u_{\theta }(x_i)-{\widehat{u}}_{\theta }(x_i)\right|$

The activation function was set to the sin function in the numerical examples. The parameter updates were based on the Adams method, with a standard learning rate of 0.001. This network model contained 100 neurons in each hidden layer, and the number of subspace layer basis functions was M. The stopping condition for the training basis function was $1\times {10}^{-3}$, and training points were sampled uniformly for a fixed random seed at 1.

4.1. Reaction–Diffusion Equation

Consider the reaction–diffusion equation in the spatio-temporal domain $\Omega =[a,b]\times (0,T]$,

$\left\{\begin{array}{l}{u}_t=u_{xx}+f(x,t),\hfill \\ u(a,t)=g_1(t),\hfill \\ u(b,t)=g_2(t),\hfill \\ u(x,0)=h(x).\hfill \end{array}\right.$

where $g_1(t), g_2(t)$ is the boundary condition and $h(x)$ is the initial condition. When $a=0,b=1,T=1$, the exact solution is $u(x)={\displaystyle \frac{1}{1+e^{\theta (x-\theta t)}}},\theta ={{\displaystyle \frac{1}{\sqrt{2}}}}_{.}$.

In Table 1, the depth of the four hidden layers remains fixed. By varying the number of training points and the quantity M of basis functions within the subspace layer, the optimal error value is achieved when $M=600$, with the error reaching as low as ${10}^{-14}$. Table 2 presents the error under different conditions of hidden layer depths and training points. It can be observed that the error of the present method ranged from ${10}^{-12}$ to ${10}^{-13}$, and when there were four hidden layers and 600 basis functions, the error reached as low as $1.75\times {10}^{-14}$. On the whole, regardless of changing the number of basis functions or the depth of the hidden layer, the error for the present method remained stable. Therefore, the optimal neural network structure obtained by this method in this example was $\left[2,100,100,100,100,M=600,1\right]$.

Table 3 compares the performance of different activation functions. It can be seen that the present method was able to stably achieve a high-precision solution under different activation functions. In the following tests, we used sin as the activation function.

Table 4 shows that different initial weights, biases, and so on are obtained by changing the random seeds. It was found that the present method is not sensitive to these changes and the error remains stable.

Table 5 compares the errors of PINNs, lbPINNs, ELMs, PIELMs, the traditional FDM, and the present method. The FDM had errors at $\tau =1/128$ and $h=1/128$, the errors of the other five methods with the same internal training points of $20\times 20$.

Among them, PINNs, lbPINNs, and the existing method are under a network structure with four hidden layers and a subspace layer $M=600$, while ELMs and PIELMs have only one hidden layer with a value of 600. It can be seen that this method achieves results comparable to the traditional FDM method and other neural network methods with fewer training points and iteration times, further demonstrating the advantages of this method in solving general problems.

Figure 3 shows the loss function evolution curves of different methods. It can be seen that PINNs and lbPINNs basically reached convergence and stability at 6000 steps. In contrast, the present method only needed to drop the loss function by three orders of magnitude, so iterations of tens of steps can achieve the purpose, avoiding the inefficiency problem caused by long iterations. The method in this paper spent far fewer iteration steps and much less CPU time than the PINN and lbPINN methods, but it achieved a much higher accuracy than them.

Figure 4 presents the exact solutions, predicted solutions, and error diagrams of the present method under different values of $\theta$. It can be observed that, when $\theta$ increases, the solutions will become steeper, which makes the function $u(x)$ change more drastically in the region close to $x=\theta t$. Meanwhile, the difficulty of solving will also increase. Nevertheless, our method can still maintain a good accuracy during the solving process.

4.2. Damped Wave Equation

This example aimed to demonstrate the effect of damping on the propagation of waves.

${\displaystyle \frac{{\partial }^{2}u}{\partial t^2}}+\alpha {\displaystyle \frac{\partial u}{\partial t}}={\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}+(1+\alpha +k^2{\pi }^2)e^t\sin (k\pi x),\left(x,t\right)\in \left[0,1\right]\times \left[0,T\right]$

The exact solution to this problem is $u\left(x,t\right)=e^t\sin (k\pi x)$.

Table 6 shows the error values under different numbers of internal training points and $M$, from which it can be seen that the minimum error corresponds to $M=600$, consistent with the optimal network structure obtained in Example 1. Therefore, in Table 5, Table 6, Table 7 and Table 8, only the optimal structure conditions are changed with different internal training points to obtain the optimal number of training points and the corresponding losses and errors. Table 6 and Table 7 show the ${‖e‖}_{L\infty }$ errors when $k=1,T=1$,$\alpha$ takes different values. It can be seen that the error at $k=1,\alpha =0.2,T=1$ remains at ${10}^{-6}-{10}^{-8}$, Figure 5 shows the exact solution, the predicted solution, and the error plot with $50\times 50$ internal training points. And the ${‖e‖}_{L\infty }$, ${‖e‖}_{L2}$ error at internal training points $40\times 40$ remains at ${10}^{-7}-{10}^{-9}$, when $k=1 ,\alpha =0.4, T=1$.

Table 8 shows that the error is optimal at the internal training points of $60\times 60$, and is $1.89\times {10}^{-7}$ when $k=1,\alpha =0.2,T=5$, Figure 6 shows the exact solution, predicted solution, and error plot under this parameter. Table 9 shows that the error is optimal at the internal training points of $100\times 100$, and is $8.32\times {10}^{-7}$ when $k=1,\alpha =0.4,T=5$. By changing the number of internal training points, the ${‖e‖}_{L\infty }$ error is very stable at ${10}^{-6}-{10}^{-7}$.

Table 10 depicts the error at different numbers of internal points when $k=1,T=5$, $\alpha =$ 30, 40, 50. We can see that the optimal internal points are $100\times 100$, $100\times 100$, and $150\times 150$ when $\alpha =$30, 40, 50. The number of required training points increases as the value of $\alpha$ increases.

Table 11 compares the error values of six methods when $\alpha ,T$ are taken with different values. The FDM method takes the spatiotemporal grid as $h=1/50,\tau =h$ when $\alpha =0.2$, $T=1$. When $\alpha =0.2,T=5$, the grid is $h=1/60,\tau =3.2h^2$. When $\alpha =0.4$ and $T=1$, the spatiotemporal grid is $h=1/40,\tau =h$; when $T=5$, the grid is $h=1/100,\tau =3.2h^2$. For PINNs, lbPINNs, and the present method, when the parameters are different, the internal training points are $50\times 50, 60\times 60, 40\times 40, 100\times 100$ and the network structure is $[2, 100, 100, 100, 100, M=600, 1]$; for the ELM and PIELM network structure of $[2,M=600,1]$. It can be seen that, under the same parameter conditions, the method proposed in this paper has a better error performance compared with those of other methods, and its precision is also very stable, truly reaching a level comparable to that of traditional methods.

Figure 7 shows the exact solution for different $k$ values, the numerical solution of the method in this paper, and the corresponding error. It can be seen that the frequency of vibration of the function values on the x-axis accelerated as k increased and the solution became more complex. However, the present method was still able to solve the case of complex solutions well and with a relatively high accuracy.

4.3. High-Wave-Number Helmholtz Equation

$u_{xx}+k^{2}u=f(x),x\in \left[0,1\right]$

$\mathrm{Exact} \mathrm{solution}:u(x)=\sin (m\pi x)$

Table 12 compares the errors under four hidden layers with 100 neurons per hidden layer at different training points and M. From the results, it seems that the present method was able to achieve a stable and high accuracy under different numbers of training points and different numbers of subspace neurons.

Table 13 shows the ${‖e‖}_{L\infty }$ error for different values of $k$ in the optimal network structure ($\left[1,100,100,100,100,M=600,1\right]$) when $m=2$. It can be seen from Table 13 that, under the conditions of different wave numbers, when the number of sampling points is changed, the performance of the method proposed in this study was quite stable and had a very high precision, with an ${‖e‖}_{L\infty }$ error within ${10}^{-13}$.

Table 14, Table 15 and Table 16 show the error values for different values of $m$. As m increased, the function graph was compressed in the $x$-axis direction and became steeper, and the difficulty of solving also increased. It can be seen that the method in this paper still maintained a high level of accuracy as m increased.

Table 17 compares the errors of PINNs, lbPINNs, ELM, PIELM, FDM, and the present method under the same parameters at different values of $m$. The PINNs, lbPINNs, and present method were under the optimal network structure of $\left[1,100,100,100,100,M=600,1\right]$; the ELM and PIELM contained only one hidden layer and took 600. By presenting the error values under different wave numbers, it can be seen that, regardless of how $m$ changed, the present method performed well, with a high precision. The ${‖e‖}_{L2}$ error values of the FDM at different $m$ values and the spatial grids $h=1/30,1/100,1/500,1/2000$ show that, compared with the present method, our method had a significantly better accuracy when the $m$ value was small. As the $m$ value increased, the accuracy of the two methods was basically the same. This also indicates that the present method performed better than traditional methods and other neural network methods.

Figure 8 shows the prediction results of the exact solution sums with different parameters and the corresponding error plots. The results show that the solution exhibited high-frequency characteristics when the parameter m increased; when m = 50, a significant jump occurred. The reason for this is that, due to complex geometric and physical phenomena, the neural network may fail to capture all the key features of the problem during the training process (underfitting). In this case, it can be eliminated simply by increasing the dataset. And the present method was able to perfectly match and capture the characteristics of the solution in the case of different forms of solutions and exhibited a superb accuracy.

4.4. The Poisson Equation Within the Complex Region

This example used the pigeon-shape interface, and had a very complex parametric interface. The interface equation is

$\begin{array}{l}x(\theta )=r\cos (\theta ) \\ y(\theta )=r\sin (\theta )\end{array}$

where $r = 0.6 + 0.216\sin (3\theta ) + 0.096\cos (2\theta ) + 0.24\cos (5\theta ), \theta \in \left[0, 2\pi \right]$.

The Poisson equation for this example is provided below:

$\begin{array}{l}-\Delta u=f(x,y),(x,y)\in \Omega ,\\ u(x,y)=g(x,y),(x,y)\in \partial \Omega .\end{array}$

And $f(x,y),g(x,y)$ satisfies certain conditions within the above parameterized interface, whose exact solution is

$u(x,y)=\sin \left(\pi x\right)\cdot \sin \left(\pi y\right).$

Figure 9 illustrates the non-regular computation region and the distribution of the training points used. There were 3000 sample points on the boundary and 30,000 sample points inside the region.

Table 18 shows the error and CPU time for different numbers of training points. The results under different training points achieved a very high accuracy, and the CPU time did not increase significantly with an increase in the number of training points, which shows a very high accuracy and efficiency.

Figure 10 illustrates the exact solution, the numerical solution, and the corresponding error for a complex regional scenario. It can be seen that the method in this paper demonstrated an excellent capability on the basis of complex regions and achieved an impressive accuracy. This synergy between efficiency and accuracy gives the method great potential for tackling complex regional challenges, paving the way for future widespread applications and heralding a new era of problem solving in related fields.

5. Discussion

The neural network-based basis function method proposed in this paper employs a two-step approach to achieve a high-accuracy approximation of PDE solutions. In the first step, the loss function is constructed based on the PDE control equations. The neural network is then trained, and the basis functions are generated through a small number of iterative steps. By excluding additional solution information, such as boundary conditions, the complexity and difficulty of the optimization problem are reduced, improving the efficiency of the basis function training process. In the second step, the base functions are used to construct the trial solution, which satisfies the control equations and solution conditions at the training points. This yields a linear system of equations for the combined coefficients of the trial solution. Highly accurate approximate solutions can then be obtained by optimizing the combination coefficients using the least-squares method. The proposed method combines the neural network approach and the traditional least-squares method, leveraging the advantages of both. It does not require mesh generation or extensive training, and it solves a system of algebraic equations to obtain the continuous function solution, rather than just discrete points. Compared to traditional methods, this approach achieves a high accuracy and efficiency.

6. Conclusions

Neural networks have achieved significant breakthroughs in various fields, particularly in solving complex PDEs. However, deep learning methods that optimize loss functions often face challenges like long training times, local optima, and a limited accuracy. This paper introduces a novel strategy inspired by extreme learning machines that uses neural networks to generate basis functions and optimize trial solutions. Our method demonstrates substantial improvements in both efficiency and accuracy compared to conventional approaches like finite difference methods and other neural network techniques such as PINNs, lbPINNs, ELM, and PIELM. The results indicate that our approach matches the accuracy of traditional methods while significantly reducing the number of iterative steps required. This advancement not only enhances the efficiency of solving PDEs, including complex domain problems, but it also challenges preconceived notions about the accuracy of neural networks. Overall, this breakthrough opens new avenues for research and applications in using neural networks for PDE solutions.

Though the proposed method has demonstrated advantages in terms of accuracy and efficiency in solving smooth PDE problems with single equations, further exploration is needed to address the challenges posed by a broader range of PDE problems with discontinuous solutions, as well as those exhibiting strongly nonlinear characteristics.

Author Contributions

Conceptualization, F.P., F.C. and Y.G.; methodology, F.P., F.C. and Y.G.; software, F.P. and F.C.; validation, F.P., F.C. and Y.G.; formal analysis, F.P. and F.C.; writing—original draft preparation, F.P. and F.C.; writing—review and editing, F.C. and Y.G. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The authors declare the availability of the data.

Acknowledgments

This work was partially supported by the National Natural Science Foundation of China (12161067), the National Natural Science Foundation of China Project (12261067), and the National Natural Science Foundation of Ningxia (2022AAC02023, 2023AAC03002).

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. The process framework for PINNs solving PDEs.

Figure 2. Neural network-based basis function construction process.

Figure 3. The loss decrease curves of PINNs, lbPINNs, and the present method.

Figure 4. Exact solution, predicted solution, and error plot under different [Forumla omitted. See PDF.] values.

Figure 5. The (left) predicted solution, (middle) exact solution, and (right) error at [Forumla omitted. See PDF.].

Figure 6. The predicted solution (a), exact solution (b), and error (c) at [Forumla omitted. See PDF.].

Figure 7. The exact solutions, predicted solutions, and error graphs corresponding to different values of [Forumla omitted. See PDF.] when [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 8. The exact solution, predicted solution, and error plot for different values of [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.].

Figure 8. The exact solution, predicted solution, and error plot for different values of [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.].

Figure 9. Computational domain and scatter plot.

Figure 10. Predicted solution (left), exact solution (middle), and error (right) for [Forumla omitted. See PDF.].

Table 1

Error with fixed hidden layer depth.

PointsM	400	500	600	700	800
$20\times 20$	5.94 $\times {10}^{-14}$	3.70 $\times {10}^{-14}$	1.75 $\times {10}^{-14}$	5.84 $\times {10}^{-14}$	3.13 $\times {10}^{-14}$
$25\times 25$	9.03 $\times {10}^{-14}$	4.47 $\times {10}^{-14}$	7.76 $\times {10}^{-14}$	4.71 $\times {10}^{-14}$	3.60 $\times {10}^{-14}$
$30\times 30$	1.77 $\times {10}^{-13}$	9.80 $\times {10}^{-14}$	7.00 $\times {10}^{-14}$	3.77 $\times {10}^{-14}$	3.57 $\times {10}^{-14}$
$35\times 35$	2.57 $\times {10}^{-13}$	8.69 $\times {10}^{-14}$	7.11 $\times {10}^{-14}$	6.11 $\times {10}^{-14}$	5.16 $\times {10}^{-14}$
$40\times 40$	1.35 $\times {10}^{-13}$	6.30 $\times {10}^{-14}$	9.79 $\times {10}^{-14}$	7.73 $\times {10}^{-14}$	6.94 $\times {10}^{-14}$
$60\times 60$	1.28 $\times {10}^{-13}$	8.09 $\times {10}^{-14}$	1.81 $\times {10}^{-13}$	1.05 $\times {10}^{-13}$	7.94 $\times {10}^{-14}$

Table 2

Error for different hidden layer depths.

Hidden LayerM	400	500	600	800
3	$4.10\times {10}^{-14}$	$1.83\times {10}^{-14}$	$1.80\times {10}^{-14}$	$1.77\times {10}^{-14}$
4	$5.94\times {10}^{-14}$	$3.70\times {10}^{-14}$	$\mathbf{1.75}\times {10}^{-14}$	$3.13\times {10}^{-14}$
5	$4.43\times {10}^{-13}$	$4.91\times {10}^{-13}$	$4.39\times {10}^{-13}$	$2.47\times {10}^{-13}$
6	$2.20\times {10}^{-12}$	$6.93\times {10}^{-13}$	$3.89\times {10}^{-12}$	$6.91\times {10}^{-13}$

Table 3

Errors corresponding to different activation functions.

Activation Function	Sin	Cos	Tanh	Sigmoid	Swish
${‖e‖}_{L2}$	7.69 $\times {10}^{-15}$	1.63 $\times {10}^{-11}$	3.86 $\times {10}^{-13}$	2.70 $\times {10}^{-10}$	1.86 $\times {10}^{-13}$
${‖e‖}_{L\infty }$	1.75 $\times {10}^{-14}$	5.14 $\times {10}^{-11}$	7.56 $\times {10}^{-13}$	4.35 $\times {10}^{-10}$	2.43 $\times {10}^{-13}$

Table 4

Errors corresponding to different random seeds.

Random Seed	1	2	5	7	10	20
${‖e‖}_{L2}$	7.69 $\times {10}^{-15}$	1.29 $\times {10}^{-14}$	2.25 $\times {10}^{-14}$	2.12 $\times {10}^{-14}$	4.10 $\times {10}^{-14}$	2.31 $\times {10}^{-14}$
${‖e‖}_{L\infty }$	1.75 $\times {10}^{-14}$	2.35 $\times {10}^{-14}$	3.47 $\times {10}^{-14}$	3.42 $\times {10}^{-14}$	7.20 $\times {10}^{-14}$	3.84 $\times {10}^{-14}$

Table 5

Comparison of the errors and CPU time of different methods.

Method	Network Structure	${‖\mathit{e}‖}_{\mathit{L}\mathit{\infty }}$	${‖\mathit{e}‖}_{\mathit{L}2}$	Epocsh	CPU
PINNs	$\left[2,100,100,100,100,M=600,1\right]$	3.04 $\times {10}^{-4}$	2.44 $\times {10}^{-4}$	6000	476.84s
lbPINNs	$\left[2,100,100,100,100,M=600,1\right]$	2.09 $\times {10}^{-4}$	1.98 $\times {10}^{-4}$	6000	496.44s
ELMs	$\left[2,M=600,1\right]$	7.18 $\times {10}^{-13}$	4.54 $\times {10}^{-13}$	0	9.34s
PIELMs	$\left[2,M=600,1\right]$	1.47 $\times {10}^{-13}$	2.88 $\times {10}^{-13}$	327	16.53s
FDM [30]	-	6.70 $\times {10}^{-14}$	4.10 $\times {10}^{-14}$	3203	0.16s
Present	$\left[2,100,100,100,100,M=600,1\right]$	1.75 $\times {10}^{-14}$	7.69 $\times {10}^{-15}$	81	16.36s

Table 6

Error at $k=1, \alpha =0.2, T=1$ with different M and training points.

PointsM	200	400	600	800	1000
$20\times 20$	9.15 $\times {10}^{-7}$	2.39 $\times {10}^{-7}$	7.31 $\times {10}^{-8}$	1.84 $\times {10}^{-7}$	1.38 $\times {10}^{-7}$
$30\times 30$	1.31 $\times {10}^{-6}$	1.58 $\times {10}^{-7}$	4.77 $\times {10}^{-8}$	1.88 $\times {10}^{-7}$	2.11 $\times {10}^{-7}$
$40\times 40$	1.19 $\times {10}^{-6}$	1.21 $\times {10}^{-7}$	1.05 $\times {10}^{-7}$	2.08 $\times {10}^{-7}$	1.90 $\times {10}^{-7}$
$50\times 50$	7.15 $\times {10}^{-7}$	4.73 $\times {10}^{-8}$	1.05 $\times {10}^{-8}$	1.81 $\times {10}^{-7}$	1.08 $\times {10}^{-7}$

Table 7

Error at $k=1, \alpha =0.4, T=1$ under different training points.

Points	$20\times 20$	$30\times 30$	$40\times 40$	$50\times 50$	$60\times 60$
${‖e‖}_{L2}$	5.47 $\times {10}^{-9}$	9.56 $\times {10}^{-9}$	2.09 $\times {10}^{-9}$	1.29 $\times {10}^{-8}$	5.38 $\times {10}^{-8}$
${‖e‖}_{L\infty }$	1.78 $\times {10}^{-8}$	2.71 $\times {10}^{-8}$	6.72 $\times {10}^{-9}$	3.87 $\times {10}^{-8}$	1.61 $\times {10}^{-7}$

Table 8

Error at different numbers of training points and subspace dimensions when $k=1, \alpha =0.2, T=5$.

Points	$50\times 50$	$60\times 60$	$80\times 80$	$100\times 100$	$120\times 120$
${‖e‖}_{L\infty }$	1.33 $\times {10}^{-5}$	1.43 $\times {10}^{-6}$	6.75 $\times {10}^{-6}$	8.16 $\times {10}^{-6}$	8.52 $\times {10}^{-6}$

Table 9

Error at different numbers of training points and subspace dimensions when $k=1, \alpha =0.4, T=5$.

Points	$50\times 50$	$60\times 60$	$80\times 80$	$100\times 100$	$120\times 120$
${‖e‖}_{L\infty }$	3.31 $\times {10}^{-6}$	5.44 $\times {10}^{-6}$	3.44 $\times {10}^{-6}$	8.32 $\times {10}^{-7}$	7.11 $\times {10}^{-6}$

Table 10

Error in the number of training points and subspace dimension at $k=1 ,T=5$.

$\mathit{\alpha }$Points	$80\times 80$		$100\times 100$		$120\times 120$		$150\times 150$
	${‖\mathit{e}‖}_{\mathit{L}2}$	${‖\mathit{e}‖}_{\mathit{L}\mathit{\infty }}$	${‖\mathit{e}‖}_{\mathit{L}2}$	${‖\mathit{e}‖}_{\mathit{L}\mathit{\infty }}$	${‖\mathit{e}‖}_{\mathit{L}2}$	${‖\mathit{e}‖}_{\mathit{L}\mathit{\infty }}$	${‖\mathit{e}‖}_{\mathit{L}2}$	${‖\mathit{e}‖}_{\mathit{L}\mathit{\infty }}$
30	1.07 $\times {10}^{-9}$	3.55 $\times {10}^{-7}$	5.42  $\times {10}^{-10}$	1.80  $\times {10}^{-7}$	8.25 $\times {10}^{-10}$	3.12 $\times {10}^{-7}$	1.23 $\times {10}^{-9}$	2.82 $\times {10}^{-7}$
40	1.28 $\times {10}^{-9}$	4.11 $\times {10}^{-7}$	1.23  $\times {10}^{-9}$	3.81  $\times {10}^{-7}$	1.39 $\times {10}^{-9}$	3.16 $\times {10}^{-7}$	2.33 $\times {10}^{-9}$	6.27 $\times {10}^{-7}$
50	1.11 $\times {10}^{-9}$	2.38 $\times {10}^{-7}$	1.15 $\times {10}^{-9}$	3.56 $\times {10}^{-7}$	1.39 $\times {10}^{-9}$	2.89 $\times {10}^{-7}$	9.38  $\times {10}^{-10}$	2.13  $\times {10}^{-7}$

Table 11

Comparison of the errors of different methods.

Method	Network Structure	$\mathit{\alpha }=0.2$				$\mathit{\alpha }=0.4$
		$\mathit{T}=1$		$\mathit{T}=5$		$\mathit{T}=1$		$\mathit{T}=5$
		${‖\mathit{e}‖}_{\mathit{L}2}$	${‖\mathit{e}‖}_{\mathit{L}\mathit{\infty }}$	${‖\mathit{e}‖}_{\mathit{L}2}$	${‖\mathit{e}‖}_{\mathit{L}\mathit{\infty }}$	${‖\mathit{e}‖}_{\mathit{L}2}$	${‖\mathit{e}‖}_{\mathit{L}\mathit{\infty }}$	${‖\mathit{e}‖}_{\mathit{L}2}$	${‖\mathit{e}‖}_{\mathit{L}\mathit{\infty }}$
PINNs	$\left[2,100,100,100,100,M=600,1\right]$	1.02 $\times {10}^{-2}$	2.46 $\times {10}^{-2}$	1.17 $\times {10}^{-2}$	2.30 $\times {10}^{+0}$	1.91 $\times {10}^{-2}$	4.60 $\times {10}^{-2}$	8.91 $\times {10}^{-3}$	2.31 $\times {10}^{+0}$
lbPINNs	$\left[2,100,100,100,100,M=600,1\right]$	3.68 $\times {10}^{-2}$	9.33 $\times {10}^{-2}$	3.95 $\times {10}^{-3}$	5.11 $\times {10}^{-1}$	4.78 $\times {10}^{-3}$	1.12 $\times {10}^{-2}$	7.46 $\times {10}^{-3}$	7.98 $\times {10}^{-1}$
ELM	$\left[2,M=600,1\right]$	1.31 $\times {10}^{-5}$	3.56 $\times {10}^{-5}$	8.24 $\times {10}^{-6}$	7.46 $\times {10}^{-4}$	3.98 $\times {10}^{-6}$	1.22 $\times {10}^{-5}$	1.98 $\times {10}^{-5}$	2.27 $\times {10}^{-3}$
PIELM	$\left[2,M=600,1\right]$	2.53 $\times {10}^{-6}$	6.60 $\times {10}^{-6}$	1.11 $\times {10}^{-6}$	9.06 $\times {10}^{-5}$	1.67 $\times {10}^{-6}$	4.47 $\times {10}^{-6}$	1.04 $\times {10}^{-6}$	9.69 $\times {10}^{-5}$
FDM [31]	-	1.82 $\times {10}^{-7}$	2.58 $\times {10}^{-7}$	2.94 $\times {10}^{-6}$	4.16 $\times {10}^{-6}$	4.41 $\times {10}^{-7}$	6.23 $\times {10}^{-7}$	3.75 $\times {10}^{-7}$	5.31 $\times {10}^{-7}$
Present	$\left[2,100,100,100,100,M=600,1\right]$	3.09  $\times {10}^{-9}$	1.05 $\times {10}^{-8}$	1.85  $\times {10}^{-8}$	1.43  $\times {10}^{-6}$	2.09  $\times {10}^{-9}$	6.72  $\times {10}^{-9}$	8.86  $\times {10}^{-9}$	8.32 $\times {10}^{-7}$

Table 12

Error at $k=1, m=2$ for different numbers of training points and M.

PointsM	300	400	500	600
30	5.34 $\times {10}^{-13}$	5.82 $\times {10}^{-13}$	4.44 $\times {10}^{-13}$	3.30 $\times {10}^{-13}$
40	7.53 $\times {10}^{-13}$	6.23 $\times {10}^{-13}$	7.53 $\times {10}^{-13}$	6.35 $\times {10}^{-13}$
50	4.41 $\times {10}^{-13}$	5.87 $\times {10}^{-13}$	5.55 $\times {10}^{-13}$	8.11 $\times {10}^{-13}$
60	4.03 $\times {10}^{-13}$	6.60 $\times {10}^{-13}$	6.04 $\times {10}^{-13}$	5.18 $\times {10}^{-13}$
80	4.44 $\times {10}^{-13}$	6.21 $\times {10}^{-13}$	5.54 $\times {10}^{-13}$	3.64 $\times {10}^{-13}$

Table 13

Error of $m=2$ at different values of $k$.

Points	60	80	90	100	140	180
$k=3$	1.38 $\times {10}^{-12}$	2.79 $\times {10}^{-13}$	1.76 $\times {10}^{-13}$	5.63 $\times {10}^{-13}$	4.46 $\times {10}^{-13}$	4.12 $\times {10}^{-13}$
$k=5\times {10}^2$	3.28 $\times {10}^{-13}$	2.87 $\times {10}^{-13}$	3.97 $\times {10}^{-13}$	2.59 $\times {10}^{-13}$	3.57 $\times {10}^{-13}$	3.35 $\times {10}^{-13}$
$k={10}^3$	4.31 $\times {10}^{-13}$	3.19 $\times {10}^{-13}$	3.82 $\times {10}^{-13}$	3.35 $\times {10}^{-13}$	3.15 $\times {10}^{-13}$	2.33 $\times {10}^{-13}$

Table 14

Errors under different training points when $m=10, k=1$.

Points	80	100	150	200
${‖e‖}_{L2}$	4.60 $\times {10}^{-11}$	1.07 $\times {10}^{-12}$	4.21 $\times {10}^{-12}$	5.57 $\times {10}^{-12}$
${‖e‖}_{L\infty }$	7.68 $\times {10}^{-11}$	2.34 $\times {10}^{-12}$	6.64 $\times {10}^{-12}$	6.43 $\times {10}^{-12}$

Table 15

Errors under different training points when $m=50, k=1$.

Points	300	500	800	1000
${‖e‖}_{L2}$	1.63 $\times {10}^{-4}$	2.61 $\times {10}^{-8}$	4.22 $\times {10}^{-7}$	3.34 $\times {10}^{-8}$
${‖e‖}_{L\infty }$	2.39 $\times {10}^{-4}$	3.16 $\times {10}^{-8}$	6.30 $\times {10}^{-7}$	3.59 $\times {10}^{-8}$

Table 16

Errors under different training points when $m=100, k=1$.

Points	1000	1500	2000	5000
${‖e‖}_{L2}$	2.05 $\times {10}^{-7}$	2.54 $\times {10}^{-7}$	8.14 $\times {10}^{-8}$	2.93 $\times {10}^{-7}$
${‖e‖}_{L\infty }$	2.85 $\times {10}^{-7}$	3.29 $\times {10}^{-7}$	8.52 $\times {10}^{-8}$	2.46 $\times {10}^{-7}$

Table 17

Comparison of the errors and CPU time of different methods.

Method	Network Structure	$\mathit{m}=2$		$\mathit{m}=10$		$\mathit{m}=50$		$\mathit{m}=100$
		${‖\mathit{e}‖}_{\mathit{L}\mathit{\infty }}$	${‖\mathit{e}‖}_{\mathit{L}2}$	${‖\mathit{e}‖}_{\mathit{L}\mathit{\infty }}$	${‖\mathit{e}‖}_{\mathit{L}2}$	${‖\mathit{e}‖}_{\mathit{L}\mathit{\infty }}$	${‖\mathit{e}‖}_{\mathit{L}2}$	${‖\mathit{e}‖}_{\mathit{L}\mathit{\infty }}$	${‖\mathit{e}‖}_{\mathit{L}2}$
PINNs	$\left[1,100,100,100,100,M=600,1\right]$	3.37 $\times {10}^{-4}$	3.13 $\times {10}^{-4}$	4.37 $\times {10}^{-1}$	5.65 $\times {10}^{-1}$	2.42 $\times {10}^{+1}$	2.73 $\times {10}^{+1}$	1.63 $\times {10}^{+0}$	1.27 $\times {10}^{+0}$
lbPINNs	$\left[1,100,100,100,100,M=600,1\right]$	6.89 $\times {10}^{-4}$	3.54 $\times {10}^{-4}$	2.00 $\times {10}^{-3}$	1.40 $\times {10}^{-3}$	2.64 $\times {10}^{-1}$	1.42 $\times {10}^{-1}$	1.01 $\times {10}^{+0}$	7.15 $\times {10}^{-1}$
ELM	$\left[1,M=600,1\right]$	4.83 $\times {10}^{-7}$	2.69 $\times {10}^{-7}$	6.28 $\times {10}^{+1}$	4.50 $\times {10}^{+1}$	4.80 $\times {10}^{+3}$	4.33 $\times {10}^{+3}$	2.06 $\times {10}^{+3}$	2.06 $\times {10}^{+3}$
PIELM	$\left[1,M=600,1\right]$	1.92 $\times {10}^{-9}$	8.42 $\times {10}^{-10}$	4.69 $\times {10}^{+1}$	4.16 $\times {10}^{+1}$	1.99 $\times {10}^{+2}$	1.35 $\times {10}^{+2}$	1.99 $\times {10}^{+2}$	1.35 $\times {10}^{+2}$
FDM [32]	-	8.20 $\times {10}^{-7}$	4.63 $\times {10}^{-7}$	1.22 $\times {10}^{-5}$	6.97 $\times {10}^{-6}$	1.24 $\times {10}^{-5}$	7.07 $\times {10}^{-6}$	1.23 $\times {10}^{-7}$	7.03 $\times {10}^{-8}$
Present	$\left[1,100,100,100,100,M=600,1\right]$	3.30 $\times {10}^{-13}$	1.82 $\times {10}^{-13}$	2.34 $\times {10}^{-12}$	1.07 $\times {10}^{-12}$	3.16 $\times {10}^{-8}$	2.61 $\times {10}^{-8}$	8.52 $\times {10}^{-8}$	8.14 $\times {10}^{-8}$

Table 18

Comparison of the errors and CPU time at the different collocation points.

${\mathit{N}}_{\mathit{i}\mathit{n}}/{\mathit{N}}_{\mathit{b}\mathit{c}}$	100/100	200/200	300/300	400/400	500/500	600/600
${‖e‖}_{L2}$	1.09 $\times {10}^{-8}$	7.90 $\times {10}^{-10}$	9.84 $\times {10}^{-10}$	1.09 $\times {10}^{-9}$	6.39 $\times {10}^{-10}$	1.15 $\times {10}^{-9}$
${‖e‖}_{L\infty }$	3.28 $\times {10}^{-8}$	1.63 $\times {10}^{-9}$	1.78 $\times {10}^{-9}$	2.44 $\times {10}^{-9}$	1.37 $\times {10}^{-9}$	1.92 $\times {10}^{-9}$
epochs	792	817	740	892	733	719
CPU(s)	13.22	18.02	19.18	31.97	29.44	31.76

References

1. Ames, W.F. Numerical Methods for Partial Differential Equations; Academic Press: Cambridge, MA, USA, 2014.

2. Morton, K.W.; Mayers, D.F. Numerical Solution of Partial Differential Equations; Cambridge University Press: Cambridge, UK, 1995.

3. Spotz, W.F.; Carey, G.F. High-order compact finite difference methods. Proceedings of the Preliminary Proceedings International Conference on Spectral and High Order Methods; Houston, TX, USA, 5–9 June 1995; pp. 397-408.

4. Liu, S.E.; Ge, Y.; Wang, S. Conservative higher-order finite difference scheme for the coupled nonlinear Schrödinger equations. Commun. Nonlinear Sci. Numer. Simul.; 2024; 131, 107797. [DOI: https://dx.doi.org/10.1016/j.cnsns.2023.107797]

5. Wang, S.; Ge, Y.; Liu, S.E. Numerical solutions of the nonlinear wave equations with energy-preserving sixth-order finite difference schemes. Comput. Math. Appl.; 2024; 168, pp. 100-119. [DOI: https://dx.doi.org/10.1016/j.camwa.2024.05.028]

6. Liu, S.E.; Ge, Y.; Qi, Y. A sixth-order finite difference method for the two-dimensional nonlinear advection diffusion reaction equation. J. Comput. Phys.; 2024; 517, 113361. [DOI: https://dx.doi.org/10.1016/j.jcp.2024.113361]

7. Tian, Z.F.; Dai, S.Q. High-order compact exponential finite difference methods for convection–diffusion type problems. J. Comput. Phys.; 2007; 220, pp. 952-974. [DOI: https://dx.doi.org/10.1016/j.jcp.2006.06.001]

8. Raissi, M.; Perdikaris, P.; Karniadakis, G.E. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J. Comput. Phys.; 2019; 378, pp. 686-707. [DOI: https://dx.doi.org/10.1016/j.jcp.2018.10.045]

9. Yu, B. The deep Ritz method: A deep learning-based numerical algorithm for solving variational problems. Commun. Math. Stat.; 2018; 6, pp. 1-12.

10. Sirignano, J.; Spiliopoulos, K. DGM: A deep learning algorithm for solving partial differential equations. J. Comput. Phys.; 2018; 375, pp. 1339-1364. [DOI: https://dx.doi.org/10.1016/j.jcp.2018.08.029]

11. Huang, G.B.; Zhu, Q.Y.; Siew, C.K. Extreme learning machine: A new learning scheme of feedforward neural networks. Proceedings of the 2004 IEEE International Joint Conference on Neural Networks; Budapest, Hungary, 25–29 July 2004; Volume 2, pp. 985-990.

12. Xiang, Z.; Peng, W.; Liu, X. Self-adaptive loss balanced Physics-informed neural networks. Neurocomputing; 2022; 496, pp. 11-34. [DOI: https://dx.doi.org/10.1016/j.neucom.2022.05.015]

13. Cuomo, S.; Di Cola, V.S.; Giampaolo, F.; Rozza, G.; Raissi, M.; Piccialli, F. Scientific machine learning through physics–informed neural networks: Where we are and what’s next. J. Sci. Comput.; 2022; 92, 88. [DOI: https://dx.doi.org/10.1007/s10915-022-01939-z]

14. Markidis, S. The old and the new: Can physics-informed deep-learning replace traditional linear solvers?. Front. Big Data; 2021; 4, 669097. [DOI: https://dx.doi.org/10.3389/fdata.2021.669097]

15. Wang, Y.; Yao, Y.; Gao, Z. An extrapolation-driven network architecture for physics-informed deep learning. Neural Netw.; 2025; 183, 106998. [DOI: https://dx.doi.org/10.1016/j.neunet.2024.106998] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/39657525]

16. Zeng, B.; Wang, Q.; Yan, M.; Liu, Y.; Chengze, R.; Zhang, Y.; Liu, H.; Wang, Z.; Sun, H. PhyMPGN: Physics-encoded Message Passing Graph Network for spatiotemporal PDE systems. arXiv; 2024; arXiv: 2410.01337

17. Wang, S.; Yu, X.; Perdikaris, P. When and why PINNs fail to train: A neural tangent kernel perspective. J. Comput. Phys.; 2022; 449, 110768. [DOI: https://dx.doi.org/10.1016/j.jcp.2021.110768]

18. Urbán, J.F.; Stefanou, P.; Pons, J.A. Unveiling the optimization process of Physics Informed Neural Networks: How accurate and competitive can PINNs be?. arXiv; 2024; arXiv: 2405.04230[DOI: https://dx.doi.org/10.1016/j.jcp.2024.113656]

19. Rathore, P.; Lei, W.; Frangella, Z.; Lu, L.; Udell, M. Challenges in training PINNs: A loss landscape perspective. arXiv; 2024; arXiv: 2402.01868

20. Yan, L.; Zhou, Y.; Liu, H.; Liu, L. An improved method for Physics-Informed Neural Networks that accelerates convergence. IEEE Access; 2024; 12, pp. 23943-23953. [DOI: https://dx.doi.org/10.1109/ACCESS.2024.3354058]

21. Zou, Z.; Meng, X.; Karniadakis, G.E. Correcting model misspecification in physics-informed neural networks (PINNs). J. Comput. Phys.; 2024; 505, 112918. [DOI: https://dx.doi.org/10.1016/j.jcp.2024.112918]

22. Müller, J.; Zeinhofer, M. Achieving high accuracy with PINNs via energy natural gradient descent. Proceedings of the International Conference on Machine Learning; Budapest, Hungary, 25–29 July 2004; pp. 25471-25485.

23. Hua, J.; Li, Y.; Liu, C.; Wan, P.; Liu, X. Physics-informed neural networks with weighted losses by uncertainty evaluation for accurate and stable prediction of manufacturing systems. IEEE Trans. Neural Netw. Learn. Syst.; 2024; 35, pp. 11064-11076. [DOI: https://dx.doi.org/10.1109/TNNLS.2023.3247163]

24. Mishra, S.; Molinaro, R. Estimates on the generalization error of physics-informed neural networks for approximating PDEs. IMA J. Numer. Anal.; 2023; 43, pp. 1-43. [DOI: https://dx.doi.org/10.1093/imanum/drab093]

25. Dwivedi, V.; Srinivasan, B. Physics informed extreme learning machine (pielm)–a rapid method for the numerical solution of partial differential equations. Neurocomputing; 2020; 391, pp. 96-118. [DOI: https://dx.doi.org/10.1016/j.neucom.2019.12.099]

26. Quan, H.D.; Huynh, H.T. Solving partial differential equation based on extreme learning machine. Math. Comput. Simul.; 2023; 205, pp. 697-708. [DOI: https://dx.doi.org/10.1016/j.matcom.2022.10.018]

27. Wang, Y.; Dong, S. An extreme learning machine-based method for computational PDEs in higher dimensions. Comput. Methods Appl. Mech. Eng.; 2024; 418, 116578. [DOI: https://dx.doi.org/10.1016/j.cma.2023.116578]

28. Dong, S.; Li, Z. Local extreme learning machines and domain decomposition for solving linear and nonlinear partial differential equations. Comput. Methods Appl. Mech. Eng.; 2021; 387, 114129. [DOI: https://dx.doi.org/10.1016/j.cma.2021.114129]

29. Xu, Z.; Sheng, Z. Subspace method based on neural networks for solving the partial differential equation. arXiv; 2024; arXiv: 2404.08223

30. Sheng, L.; Yong, G.; Ying, Q. High-precision finite-difference method for the nonlinear reaction-diffusion equation. Appl. Math.; 2023; 36, pp. 747-755.

31. Meng, W.; Yong, G.; Zhi, W. Explicit high-order compact difference method for solving nonlinear hyperbolic equations with three types of boundary conditions. Wave Motion; 2023; 118, 103120.

32. Zhi, W.; Yong, G. A high precision compact difference method for solving the Helmholtz equation of variable coefficient. J. Chongqing Norm. Univ.; 2019; 36, pp. 92-97.