Advanced machine learning technique for solving elliptic partial differential equations using Legendre spectral neural networks

1.   Introduction

Among the most frequent outcomes of mathematical modeling for describing physical, biological, or chemical phenomena are partial differential equations (PDEs), which are used in all branches of engineering and science. In addition, a dramatic increase has been seen in recent years in the use of these PDEs in fields like chemistry, computer sciences (especially in relation to image processing and graphics), and economics (finance) ^[1]. These PDEs are subject to some appropriate initial and boundary conditions. Solving these PDEs and analyzing their solutions is one of the key factors. In most cases, the closed form solution of these PDEs is not possible or is very hard to obtain. Therefore, one must look for some approximate method to find their solution. Most of these methods are based on predictor–corrector, Runge–Kutta, finite difference, and finite element ^[2,3,4]. In order to use these numerical approaches, the domain must be discretized into a number of finite domains or locations where the functions are locally approximated. However, such methods still have many limitations, despite their strength and potential to help us approximate solutions to issues where analytical approaches are unable to. Finite difference techniques involve the creation of a grid on which the solution is to be found and can be difficult to specify over complicated geometries. Similar to the finite difference approach, the finite element method is simple to define over complicated geometries but can struggle with nonlinear issues and may become slow or unfeasible when dealing with high-dimensional problems or when an extraordinarily small error is needed. There is still a need for innovative numerical techniques to solve PDEs because of these constraints ^[5].

Elliptic PDEs are a class of PDEs that explain the global and steady-state phenomena that arise in many areas of computational sciences, such as fluid dynamics, heat transfer, electromagnetism, geophysics, biology, and other application areas. They are used by coastal engineers to approximate the motion of the sea, and are used to model electric potential. Many elliptical problems can be solved easily, but those with complicated meshes or complex geometries are more difficult to solve. Only in very rare situations can closed analytical expressions for their solutions be found, most of which are of not much theoretical or practical significance. Apart from this, many elliptical PDEs are defined in higher-dimensional domains, which cause computational complexities for traditional numerical schemes. Therefore, it is natural that mathematicians and scientists look for techniques for estimating solutions. In fact, the development of computational mathematics, which focuses primarily on the development and mathematical analysis of numerical methods for the approximate solution of elliptic PDEs, has been accelerated by the introduction of digital computers.

Neural network techniques are emerging as effective tools for resolving PDEs, especially elliptical PDEs when alternative approaches are impractical or may not work. If one can figure out a way to sample from the underlying domain, neural networks can be used to solve high-dimensional problems without suffering too much from the curse of dimensionality, as well as nonlinear partial differential equations, and even domains with complex geometries. Neural network-based approaches were employed whenever the need for highly precise solutions to PDEs was in demand. In terms of accuracy, the neural networks that were able to acquire any nonlinear function performed more effectively than the conventional techniques. More than two and a half decades ago, Lagaris et al. ^[6] applied neural networks to solve ordinary and PDEs, which marked the beginning of the use of neural networks in the solution of differential equations. Since then, there have been a number of advancements in the use of neural networks in the numerical solution of differential equations. For example, Mall and Chakraverty in ^[7], and Dufera in ^[8] introduced a deep neural network-based approach to solve ordinary differential equations (ODEs). The neural network was trained using the traditional gradient descent approach. Rivera et al. in ^[9], proposed a method to solve the PDEs using neural networks after determining the loss function using the quadrature approach. The quadrature weight was added to the loss function using this application of the quadrature rule. Tan et al. ^[10] solved two-dimensional elliptic PDEs using a machine learning technique based on the unsupervised form of wavelet neural networks. Recently, Sabir et al. ^[11] used artificial neural network procedures to solve the susceptible-exposed-infectious-recovered (SEIR) mathematical model for Zika virus transmission, together with the hybridization effectiveness of local search schemes and global swarming. An accurate, fast, and reliable neural network approach based on physics that solves differential equations was used in ^[12]. For more details about the use of neural network technique to solve integer-order differential equations, we refer the reader to ^{[13,14,15,16,17,18]}. Apart from solving the integer-order differential equations, the neural network technique has also been effectively applied to solve fractional-order differential equations ^[19,20].

The development of techniques for solving PDEs has advanced significantly in scientific computing over the past years. Because of their natural hierarchical structure, links to approximation theory, and advantageous convergence qualities, spectral methods are an important component of scientific computing's tools. For the purpose of guaranteeing that the underlying PDE is satisfied in a suitable sense, spectral methods often expand the solution of a PDE as a linear combination of basis functions and estimate the coefficients of the linear combination. Despite their potential potency, spectral approaches' are highly dependent on the basis function selection, which is not always clear for practical uses. The geometry of the domain where the answer is to be approximated can be a source of difficulties. For instance, while commonly used basis functions, such as orthogonal polynomials, are appropriate only for regular domains, fluid dynamics applications usually use complicated domains. The existence of highly localized characteristics in the solution, such as excessively steep gradients, might also cause issues. The approximation of the order function characterizing the evolving sharp phase boundary is one example of an application in phase field modeling. Unless the specifics of the application are taken into account while building the basis functions, the resolution of such localized features can reduce the effectiveness of a spectrum approach because the basis functions employed in spectral methods are global in nature ^{[21,22,23,24,25]}. For more details about spectral methods, we refer the reader to ^{[26,27,28,29,30,31]}.

Among the significant sequences of orthogonal polynomials that have been extensively investigated and are utilized in computational fluid dynamics, interpolation and approximation theory, numerical integration, and the solution of differential equations, particularly elliptic differential equations, are the Legendre polynomials. These are an essential part of numerical integration and the approximate solution of differential and integral equations, along with being an effective tool for approximating hard-to-calculate functions. Excellent error characteristics are demonstrated by the Legendre spectral methods when approximating a smooth function. Orthogonal polynomial expansion appears in many real-world scenarios and applications and is crucial to many branches of physics and mathematics ^[32]. Typical neural network structures have several disadvantages, despite the fact that they have generated outstanding outcomes when it comes to solving differential equations. For example physics-informed neural networks (PINNs) have a weaker theoretical underpinning and a lack of useful error analysis tools due to their recent inception as a numerical method. Since the provided loss function in PINNs is frequently non-convex, it might be difficult to comprehend the convergence criteria in the optimization process, where solutions may become trapped in the local minima. It is necessary to overcome problems such as weak generality and low solution precision, and, in some cases, training can be difficult, particularly when the equation's solution has high-frequency features ^[33]. Other issues with conventional neural networks, including the local minimum, poor convergence speed, and initial weight values. As a result, numerous scientists are investigating ways to overcome these obstacles by adjusting the network architecture, and training methodologies, as well as other strategies to overcome these challenges.

Traditional spectral methods, such as Fourier and Chebyshev spectral methods, rely on solving systems of linear or nonlinear equations derived from the discretization of the governing PDE. This process involves matrix assembly and factorization, which can become computationally expensive for high-dimensional problems or irregular domains. Additionally, the global nature of traditional spectral basis functions may require a dense representation of the solution domain, further increasing computational effort. In contrast, the Legendre spectral neural network (LSNN) bypasses the need for matrix assembly by directly minimizing a loss function through an optimization process. The elimination of hidden layers and the use of Legendre–Gauss–Lobatto points significantly reduce the computational overhead. Moreover, the LSNN employs automatic differentiation for evaluating derivatives, which simplifies the implementation and avoids the overhead associated with numerical differentiation. Inspired by ^[34], a novel single-layer neural network is developed in this study that is based on the Legendre spectral method, which is a natural choice due to its exponential rate of convergence combined with its inherited property to satisfy the boundary conditions without making any adjustment, and a robust strategy for sampling in the form of Legendre–Gauss–Lobatto points enhances the training speed and accuracy. Neural network methods based on Legendre polynomials have been extensively used for solving differential equations in ^{[35,36,37,38,39,40,41]}. Most recently, the PINN technique for solving the PDEs, mainly the elliptic PDEs has been successfully used in ^{[42,43,44,45,46,47,48]}. Among these methods, the LSNN method holds significant potential for real-world applications across various fields. In fluid dynamics, it can efficiently model steady-state flows and heat transfer in complex geometries. In material science, LSNNs can simulate stress-strain analysis and thermal conduction with high precision.

The remainder of the paper is structured as follows. In Section 2, we demonstrate the proposed scheme for both regular and irregular geometries in details. Section 3 provides some error analysis, followed by a number of numerical examples in Section 4. Section 5 presents the concluding remarks.

2.   Description of the Legendre spectral neural network method

Consider the general differential equation for which we seek the solution $v(\boldsymbol{x})$, defined over a domain $\mathcal{D} \subset \mathbb{R}^d$, of the form:

where $\mathcal{G}$ is the operator defining the structure of the differential equation, $v(\boldsymbol{x})$ is the unknown solution, and $\nabla, \nabla^2, \dots, \nabla^m$ represent the gradient, Hessian, and higher-order derivatives of $v(\boldsymbol{x})$. Equation (2.1) is subject to boundary conditions (BCs) given by (2.2)

where $\mathcal{H}$ represents the boundary operator, imposing conditions such as Dirichlet, Neumann, or Robin.

To solve the model equation (2.1), subject to the BCs given by Eq (2.2) using the LSNN method, we construct the approximate function space for the solution based on Legendre polynomials of the form:

where $L_k(x)$ denotes the $k$-th Legendre polynomial and $\gamma_k, \delta_k$ are coefficients determined by the boundary conditions. Legendre–Gauss–Lobatto (LGL) points are special sampling points that include the endpoints of the domain, $-1$ and $1$, ensuring that the boundary conditions are automatically satisfied. These points are the roots of the equation:

The corresponding weights for the LGL points are given by:

where $x_i$ are the LGL points, $k$ is the polynomial degree, and $w_i$ is the weight associated with $x_i$.

2.1.   LSNNs for rectangular geometry

Consider a one-dimensional boundary value problem with homogeneous boundary conditions:

If the boundary conditions are non-homogeneous, the problem is homogenized by introducing a function $W(x)$ such that:

The solution $v(x)$ is then written as:

where $\hat{v}(x)$ satisfies the homogeneous boundary conditions. In this case, the coefficients $\gamma_k, \delta_k$ given in Eq (2.1) are given by:

where the determinant $\Delta_k$ is given by:

For homogeneous Dirichlet boundary conditions, the basis functions simplify to:

The approximation space for one-dimensional problems is represented as:

For higher-dimensional domains, such as $\mathcal{D} = [a, b] \times [c, d]$, the basis functions are formed by tensor products:

The detailed structure of an LSNN is shown in Figure 1. In the context of LSNNs, the residual points for training are chosen as LGL points, which ensures that the spectral properties of Legendre polynomials are fully utilized and that the solution space achieves high accuracy.

The residual loss function for the LSNN is formulated as:

where $\Phi$ represents the trainable parameters, and $\mathcal{T}$ is the set of residual points sampled at the LGL points. By minimizing this loss, the trainable parameters $\Phi$ are updated to approximate the solution $v(\boldsymbol{x})$. Automatic differentiation (AD) is used to compute derivatives of the neural network output $\hat{v}(x; \Phi)$, allowing the evaluation of all differential terms in the governing equations $\mathcal{G}$. For a given output $\hat{v}(x; \Phi)$, AD is used to compute ^[49].

which eliminates the need for manual differentiation or finite difference approximations.

In order to minimize the loss and ensure convergence to the desired solution, an advanced optimizer based on the Adam optimizer is used, which combines momentum and adaptive learning rates for efficient training and updates the parameters by ^[50].

where $\hat{m}_t$ and $\hat{v}_t$ are the bias-corrected moment estimates of the gradients. This procedure is summarized in Algorithm 1.

2.2.   LSNNs for complex geometry

For complex geometries, the LSNN employs coordinate transformations to map irregular domains into rectangular ones. Consider a domain described in polar coordinates:

The mapping is defined as:

This transforms the original domain to:-

LGL points are used in the radial direction $z$, while uniform or Fourier sampling is applied in the angular direction $\theta$. This approach ensures spectral-level resolution in the radial direction while preserving periodicity in the angular direction. Agian, AD is used for computing mixed partial derivatives with respect to the transformed variables. For the neural network output $\hat{v}(z, \theta; \Phi)$, AD computes:

which ensures efficient handling of nonlinear terms and mixed derivatives arising from transformations. Algorithm 2, summarizes this procedure.

3.   Error analysis

Consider the model in Eq (2.1), subject to the BCs given in Eq (2.2). The total error $E$ in can be expressed as:

where $\hat{v}$ denotes the neural network solution, $v_N$ corresponds to an approximation using the Legendre polynomial basis, and $v$ is the exact solution.

The term $\| \hat{v} - v_N \|$ accounts for the error introduced by the neural network optimization process, while $\| v_N - v \|$ is the approximation error inherent in the Legendre polynomial representation.

For a one-dimensional boundary value problem of the form:

with the boundary conditions:

the solution $v_N(x)$ satisfies the variational formulation:

Let $v \in H^m(I)$ and $f \in H^k(I)$, where $H^m(I)$ is the Sobolev space. The error bounds for the Legendre spectral approximation are:

In Eq (3.5), $\| \cdot \|_{1}$ refers to the $H^1$-norm, which measures the error in both the function $v$ and its first derivative and is defined by:

The term $N^{1-m} \| \partial_x^m v \|$ indicates that the approximation error decreases with increasing smoothness $m$ of the solution $v$, while $N^{-k} \| \partial_x^k f \|$ accounts for the smoothness $k$ of the source term $f$, which demonstrates that both the solution $v$ and the source term $f$ significantly influence the overall error in the $H^1$-norm.

In Eq (3.6), $\| \cdot \|_2$ refers to the $L^2$-norm, defined as:

In this equation, $N^{-m} \| \partial_x^m v \|$ represents the contribution of the smoothness $m$ of $v$ to the error, which decays faster than the $H^1$-norm, as derivatives are not considered. The term $N^{-k} \| \partial_x^k f \|$ reflects the influence of the regularity $k$ of the source term $f$ on the approximate error. The constants in both Eqs (3.5) and (3.6) depend on the Sobolev norms of $v$ and $f$. The $H^1$-norm includes derivative contributions, which may lead to a slower convergence rate compared with the $L^2$-norm, which only accounts for the function values.

4.   Numerical examples

To confirm the exponential rate of convergence, in this section a number of numerical examples taken from ^[34]. In our computations, we use a PC with a 12th Gen Intel(R) Core(TM) i7-1255U processor, and 16 GB RAM. The accuracy of the LSNN method is evaluated using the following three standard error metrics:

Example 4.1. We consider the following second-order differential equation:

subject to the Robin boundary conditions:

The exact solution of this problem is given as:

To solve this problem using the LSNN approach, we use different values of $N = 12, 14, 16, 18, 20, 22$ to construct the networks. Each network is trained for $300 \times N$ epochs. The initial learning rate is set to $0.1$ and dynamically adjusted using a scheduler that reduces it by a factor of $0.7$ if the loss function does not decrease for 600 consecutive iterations. The computational domain is uniformly sampled with 501 test points within the interval $[-1, 1]$. The results demonstrate that as $N$ increases, the solution error decreases rapidly. This highlights the efficiency and accuracy of the LSNN method compared with traditional approaches. The comparison between the exact and LSNN solution is shown in Figure 2. The errors obtained for various $N$ values confirm the exponential convergence characteristic of the LSNN, showcasing its ability to achieve high precision even with relatively small networks, as shown in Figure 3.

Example 4.2. Consider the two-dimensional Poisson equation defined as:

where the source term $g(x, y)$ is determined by the exact solution:

The boundary conditions are Dirichlet and given as:

where these conditions satisfy the following continuity relations:

The trial solution is constructed as:

where $\tilde{v}(x, y; \theta)$ represents the output from the LSNN and $A(x, y)$ is a function that satisfies the nonhomogeneous boundary conditions:

The Adam optimizer is used for training with an initial learning rate of 0.1 and a maximum of 2000 iterations. A grid of 300 equally spaced points in both $x$ and $y$ directions is used to evaluate the solution. Figure 4 provides a comparison of the LSNN and the exact solution. The numerical experiments reveal that the LSNN approach achieves high accuracy, with the errors decreasing significantly as the network complexity $N$ increases, as shown in Table 1.

Example 4.3. Consider the following two-dimensional (2D) Helmholtz equation on an annular domain:

where the source term $g(x, y)$ and the Dirichlet boundary conditions are determined by the exact solution:

The annular domain is transformed into a rectangular parameter space using the coordinate mapping:

where $s, t \in [-1, 1]$. In the transformed domain, the Helmholtz equation becomes:

In the $s$-direction, Legendre basis functions are used, with sampling points chosen as LGL nodes. For the $t$-direction, due to periodicity, trigonometric basis functions are employed, with uniform sampling over the interval $[-1, 1]$. The loss function is defined as:

where $\mathcal{T}$ is the training set of sampled points. The initial learning rate is set to 0.01 and dynamically adjusted by reducing it by a factor of 0.6 if the loss does not decrease after 500 iterations. Using an LSNN network with $N = 20$, the training converges after 2000 iterations. The results are illustrated in Table 2, which shows the error behaviors between the exact solution and numerical solution at different collocation points. Figure 5 illustrates the LSNN and the exact solution.

To compare the efficiency and accuracy of the LSNN with traditional methods, we also applied the PINNs approach. The PINNs' setup included a network with four hidden layers, each with 128 neurons, and the same learning rate of 0.001. After 2000 iterations, the convergence of PINNs' solution with that of the LSNN is shown in Figure 6. This experiment highlights the ability of LSNNs to solve Helmholtz equations efficiently on non-standard domains, achieving high accuracy with significantly reduced computational effort.

Example 4.4. In this example, we address the Poisson equation on a circular domain:

where the source term $g(x, y)$ and the Dirichlet boundary conditions are determined by the exact solution:

To simplify the problem, the circular domain is mapped to a rectangular parameter space. Using polar coordinates, the transformation is defined as:

where $r \in [0, 1]$ and $\theta \in [0, 2\pi)$. In this transformed domain, the Poisson equation takes the following form:

The circular domain introduces an additional challenge: the value at the center of the circle ($r = 0$) is undetermined. To address this, an auxiliary single-parameter neural network $\hat{v}_0$ is introduced to approximate the value at the center. This parameter is optimized alongside the main LSNN. The loss function is defined as:

where $\mathcal{T}$ is the set of training points. Using the same parameter values as in Example 4.3, Figure 7 provides a comparison of the LSNN and true solution, while Figure 8 shows that error behavior for different collocation points.

The success of the LSNN method on this domain depends significantly on the learning rates for both the main and auxiliary networks. Improper learning rate settings can result in slow convergence or failure to reach the global minimum. To mitigate this issue, periodic corrections to the auxiliary network parameter $\hat{v}_0$ are applied, ensuring better alignment with the true solution.

Example 4.5. Consider a three-dimensional (3D) elliptic PDE with variable coefficients:

where the coefficients are given as:

The source term $g(x, y, z)$ and Dirichlet boundary conditions are determined by the exact solution:

The computational domain is a rolled-up cylindrical shell defined as:

where $r, \theta$ are the cylindrical coordinates. The domain is mapped to a rectangular region using the coordinate transformation:

where $r \in [0.5, 1]$, $\theta \in [0, 2\pi)$, and $z \in [-1, 1]$. In the transformed domain, the elliptic operator becomes:

In the $r$- and $z$-directions, Legendre basis functions are used with LGT points for sampling. In the $\theta$-direction, due to its periodicity, trigonometric basis functions are employed with uniform sampling. The loss function is defined as:

where $\mathcal{T}$ represents the training points. Again keeping the same parameter values as in the previous example, the true and computed solution is plotted in Figure 9. The error behaviors for different norms and at different collocation points are shown in Figure 10, which confirms the exponential convergence of our proposed scheme.

The LSNN method demonstrates excellent performance in solving this 3D problem, even for irregular geometries like cylindrical shells. The method achieves high accuracy and rapid convergence with modest computational resources, highlighting its robustness for variable-coefficient PDEs in complex domains.

5.   Conclusions

In this work, we demonstrate the ability of LSNNs to handle the problems of regular as well as irregular domains. We show, through a number of numerical examples that our proposed scheme is more stable and computationally efficient compared with other neural network techniques, e.g., PINNs. In particular, our scheme overcomes issues related to training instability and error oscillation and achieves an exponential rate of convergence. Our numerical examples, especially those related to complex geometries of annular and cylindrical shells, further enhance its potential to deal with real-world problems across various scientific disciplines. Moreover, the network structure is made simple by eliminating the presence of hidden layers using a Legendre-based neuron layout that improves training efficacy without compromising the solution's accuracy. The present study highlights the effectiveness of the LSNN method; however, its convergence may be impacted when using a large number of nodes due to numerical instabilities and optimization challenges. Future work could focus on adaptive strategies and improved optimization techniques to address these limitations and extend it to time-dependent and nonlinear equations.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the Deanship of Scientific Research, Vice Presidency for Graduate Studies and Scientific Research, King Faisal University, Saudi Arabia [Grant No. KFU250330].

Conflict of interest

The author declares there are no conflicts of interest.