A Novel Fourier Feature Network for Solving Partial Differential Equations

# A Novel Fourier Feature Network for Solving Partial Differential Equations
Source: [https://arxiv.org/html/2608.14733](https://arxiv.org/html/2608.14733)
\[style=chinese, orcid=0000\-0002\-8398\-7212\]

\[style=chinese\]

\[style=chinese\]

\[style=chinese\]

Qihong Yangyangqh0808@163\.comAddress:organization=School of Mathematics, Sichuan University, city=Chengdu, postcode=610065, country=ChinaYangtao Dengytdeng1998@foxmail\.comQiaolin Heqlhejenny@scu\.edu\.cn

###### Abstract

Building on the foundation of single\-hidden\-layer neural networks, Fourier Feature Networks \(FENs\) are proposed, which incorporate Fourier features usingcos\\cos,sin\\sin, or a combination of both\. Similar to Extreme Learning Machines \(ELMs\), FENs employ a single\-hidden\-layer architecture to generate a set of basis functions\. The target function is then approximated as a linear combination of these basis functions, with the coefficients determined using the least squares method\. However, unlike ELMs, which often rely on affine transformations to improve representational power, FENs can achieve high\-precision solutions without requiring such transformations on the input variables\. To evaluate the representational capacity of these networks, we search for an optimal scaling factor within a predefined range for the randomly initialized and fixed weights and biases\. By adjusting this scaling factor, we ensure a fair comparison between FENs and ELMs using various activation functions, such assigmoid,tanh\\tanh, andswish\. Our numerical experiments demonstrate that FENs consistently achieve higher accuracy than ELMs\.

###### keywords

Neural networks ,Function approximation ,Fourier features ,Least squares method ,Partial differential equations

††corresponding:Corresponding author††credit:Conceptualization of this study, Methodology, Writing \- Original draft preparation## 1Introduction

In recent years, the application of neural networks in scientific computing has become increasingly widespread\. Notably, the development of methods such as the Deep Ritz Method \(DRM\)\[[38](https://arxiv.org/html/2608.14733#bib.bib18)\], the Deep Galerkin Method \(DGM\)\[[33](https://arxiv.org/html/2608.14733#bib.bib17)\], and Physics\-Informed Neural Networks \(PINNs\)\[[28](https://arxiv.org/html/2608.14733#bib.bib16)\]has garnered significant attention\. A substantial body of work\[[39](https://arxiv.org/html/2608.14733#bib.bib19),[21](https://arxiv.org/html/2608.14733#bib.bib20),[14](https://arxiv.org/html/2608.14733#bib.bib21),[15](https://arxiv.org/html/2608.14733#bib.bib22),[1](https://arxiv.org/html/2608.14733#bib.bib23)\]has demonstrated the tremendous potential of neural networks in this domain, particularly in the numerical solution of partial differential equations \(PDEs\)\.

Current neural network\-based methods for solving PDEs can be broadly divided into two main categories\. The first category is training\-based, where the goal is to minimize the residual of the PDEs to fit the target function\. This is typically achieved using gradient\-based optimization algorithms such as Adam\[[18](https://arxiv.org/html/2608.14733#bib.bib24)\]or L\-BFGS\[[23](https://arxiv.org/html/2608.14733#bib.bib25)\]\. Methods like the DRM\[[38](https://arxiv.org/html/2608.14733#bib.bib18)\], the DGM\[[33](https://arxiv.org/html/2608.14733#bib.bib17)\], and the PINNs\[[28](https://arxiv.org/html/2608.14733#bib.bib16)\]all fall under this category\. A closely related approach is operator learning\[[19](https://arxiv.org/html/2608.14733#bib.bib26),[20](https://arxiv.org/html/2608.14733#bib.bib27),[25](https://arxiv.org/html/2608.14733#bib.bib28)\], which aims to approximate mappings between infinite\-dimensional Banach spaces using data\. Although operator learning differs in formulation and application scope, it still fundamentally relies on training\-based optimization\. The second category comprises randomized neural network approaches, which do not involve iterative training\. Instead, they generate a set of basis functions using a randomly initialized single\-hidden\-layer neural network\. The target function is represented as a linear combination of these basis functions, and the PDE is discretized into a system of linear equations in terms of the coefficients of this combination\. These coefficients are then obtained using the least squares method\. Once the network is initialized, the weights and biases remain fixed; only the linear coefficients are optimized\. This approach can be viewed either as a least\-squares\-based algorithm or as a randomized neural network method\. Representative techniques include the Random Feature Method \(RFM\)\[[4](https://arxiv.org/html/2608.14733#bib.bib33),[5](https://arxiv.org/html/2608.14733#bib.bib34),[3](https://arxiv.org/html/2608.14733#bib.bib35)\], the Randomized Neural Network with Petrov–Galerkin methods \(RNN\-PG\)\[[29](https://arxiv.org/html/2608.14733#bib.bib29),[31](https://arxiv.org/html/2608.14733#bib.bib30),[30](https://arxiv.org/html/2608.14733#bib.bib31),[36](https://arxiv.org/html/2608.14733#bib.bib32)\], and the Hidden\-Layer Concatenated Extreme Learning Machine \(HLConcELM\)\[[27](https://arxiv.org/html/2608.14733#bib.bib36)\]\.

Despite extensive research on training\-based neural networks for solving PDEs, these methods often ‌exhibit limited accuracy‌ and typically require ‌considerable computational time\[[7](https://arxiv.org/html/2608.14733#bib.bib37)\]\. In contrast, randomized neural networks have gained significant attention in scientific computing due to their ability to solve PDEs both efficiently and accurately\. However, this class of methods also presents several notable implementation challenges\.\\chaddedExtreme Learning Machine \(ELM\) is a single\-layer feedforward neural network proposed in\[[16](https://arxiv.org/html/2608.14733#bib.bib1)\], which may not yield satisfactory results likely due to their reliance on Xavier or Kaiming initialization methods\. Dong and Yang\[[10](https://arxiv.org/html/2608.14733#bib.bib40)\]presented a method for computing the optimal or near\-optimal value ofRmR\_\{m\}based on the differential evolution algorithm in ELM\. The RFM closely resembles Local Extreme Learning Machines \(locELMs\)\[[9](https://arxiv.org/html/2608.14733#bib.bib15)\], as both employ a domain decomposition strategy known as Partition of Unity \(PoU\)\. This technique enables the use of ELMs to solve PDEs through the strong form of the equations\. While effective, domain decomposition significantly increases the complexity of the algorithm’s implementation\.\\chaddedThe numerical experiments can not display good results in some complex cases\. In the RNN\-PG, the solution is formulated through the weak form of the PDEs\. As a result, mesh generation and numerical integration are essential components of the algorithm, which further increase its implementation complexity\. The HLConcELM differs from the above methods by modifying the architecture of ELMs, specifically, by adding an additional hidden layer and concatenating the outputs of the hidden layers\. This architectural enhancement improves accuracy but also introduces additional computational overhead\. In particular, when derivatives are computed using automatic differentiation\[[2](https://arxiv.org/html/2608.14733#bib.bib38)\], complex network architectures can lead to significantly longer runtimes, especially when a large number of basis functions are involved\.

Although these algorithms can achieve significantly higher accuracy than training\-based neural networks, they still fall short of achieving machine precision\. To address this while maintaining low computational cost and avoiding increased network complexity, we focus on using ELMs for solving PDEs, deliberately steering clear of complexity\-increasing operations such as domain decomposition and numerical integration\. To further enhance the accuracy of ELMs without sacrificing simplicity, we propose the integration of Fourier features into the network architecture\. Specifically, we introduce Fourier Feature Networks \(FENs\), which employ activation functions based oncos\\cos,sin\\sin, or a combination of both\. In this work, we present and evaluate three types of FENs, each corresponding to one of these activation strategies, and conduct a detailed comparison of their accuracy against traditional ELMs activated bysigmoid,tanh\\tanh, andswishfunctions\.\\chaddedThe proposed FENs do not apply affine transformations to the input variables\. Moreover, we use a uniform initialization method to initialize the weights and biases, and employ a method to search for the optimal scaling factor within a given range, which enables FENs to achieve optimal performance\.

This article is organized as follows\. The network architecture and algorithmic details of ELMs, including the activation functions ofsigmoid,tanh\\tanh, andswishare introduced in Section[2](https://arxiv.org/html/2608.14733#S2)\. We propose Fourier Feature Networks, along with the three distinct activation modes of the FENs in Section[3](https://arxiv.org/html/2608.14733#S3)\. A series of numerical experiments designed to validate the efficacy of our methods are presented in Section[4](https://arxiv.org/html/2608.14733#S4)\. The article concludes with Section[5](https://arxiv.org/html/2608.14733#S5), where we summarize our results, discuss the implications of our work, and suggest directions for future research\.

## 2Preliminaries

### 2\.1The neural feature space

Neural networks are widely recognized as nonlinear mapping functions that transformdd\-dimensional inputs into either lower or higher dimensional spaces\. Architecturally, a neural network typically comprises three components: an input layer, multiple hidden layers, and an output layer\. In this context, the neural feature space refers to the functional space represented by the outputs of the last hidden layer\.

Assume that the input vector is𝒙∈ℝd\\boldsymbol\{x\}\\in\\mathbb\{R\}^\{d\}\. The nonlinear mapping represented by the hidden layers of the neural network is given byϕi​\(𝒙\)\\phi\_\{i\}\(\\boldsymbol\{x\}\), where1≤i≤M1\\leq i\\leq M, andMMis the number of outputs of the last hidden layer\.

Consequently, the neural feature space, denoted as𝒫N​N\\mathcal\{P\}\_\{NN\}, is defined as a space spanned by the basis functionsϕi\{\\phi\_\{i\}\}, i\.e\.,

𝒫N​N=s​p​a​n​\{ϕ1,ϕ2,⋯,ϕM\}\.\\mathcal\{P\}\_\{NN\}=span\\\{\\phi\_\{1\},\\phi\_\{2\},\\cdots,\\phi\_\{M\}\\\}\.\(1\)
In neural networks, the output layer typically computes a linear combination of the outputs from the last hidden layer\. Here, we assume that the output layer has a single output\. Therefore, the function represented by the neural network can be expressed as

uM​\(𝒙\)=∑i=1Mwi​ϕi\.u\_\{M\}\(\\boldsymbol\{x\}\)=\\sum\_\{i=1\}^\{M\}w\_\{i\}\\phi\_\{i\}\.\(2\)

### 2\.2Extreme learning machine

Assume that the network hasddneurons in the input layer, corresponding to the input vector𝒙∈ℝd\\boldsymbol\{x\}\\in\\mathbb\{R\}^\{d\}\. The main idea of the ELM is to randomly initialize the weights𝑾∈ℝM×d\\boldsymbol\{W\}\\in\\mathbb\{R\}^\{M\\times d\}and the biases𝒃∈ℝM\\boldsymbol\{b\}\\in\\mathbb\{R\}^\{M\}between the input and hidden layers\. Then, the weights𝒘∈ℝM\\boldsymbol\{w\}\\in\\mathbb\{R\}^\{M\}from the hidden layer to the output layer are directly computed\. This can be mathematically formulated as follows:

uM​\(𝒙\)=𝒘T​σ​\(𝑾​𝒙\+𝒃\),u\_\{M\}\(\\boldsymbol\{x\}\)=\\boldsymbol\{w\}^\{T\}\\sigma\(\\boldsymbol\{W\}\\boldsymbol\{x\}\+\\boldsymbol\{b\}\),\(3\)whereσ\\sigmais the activation function, which acts elementwise on the vector𝑾​𝒙\+𝒃\\boldsymbol\{W\}\\boldsymbol\{x\}\+\\boldsymbol\{b\}, anduM​\(𝒙\)u\_\{M\}\(\\boldsymbol\{x\}\)is the function represented by the ELM\. We only need to calculate𝒘\\boldsymbol\{w\}\.

### 2\.3Affine transformation

As we know, normalization plays a vital role in enhancing model performance and generalization capabilities in neural network training\[[17](https://arxiv.org/html/2608.14733#bib.bib2)\]\. Normalization accelerates the training process, enhances model generalization, prevents overfitting, and improves robustness to variations in weight initialization methods\. In scientific computing, when neural networks are used for function approximation or solving PDEs, normalization serves as an essential preprocessing step to ensure numerical stability\. In this context, affine transformations are commonly applied\.

Let the input variables be𝒙∈Ω⊂ℝd\\boldsymbol\{x\}\\in\\Omega\\subset\\mathbb\{R\}^\{d\}, whereΩ\\Omegais a closed set\. The affine transformation maps the input variables𝒙\\boldsymbol\{x\}to a new vector𝒙~∈\[−1,1\]d⊂ℝd\\tilde\{\\boldsymbol\{x\}\}\\in\[\-1,1\]^\{d\}\\subset\\mathbb\{R\}^\{d\}, which is given by

Ti​\(xi\)=2​xi−xi\(l\)xi\(u\)−xi\(l\)−1,T\_\{i\}\(x\_\{i\}\)=2\\frac\{x\_\{i\}\-x\_\{i\}^\{\(l\)\}\}\{x\_\{i\}^\{\(u\)\}\-x\_\{i\}^\{\(l\)\}\}\-1,\(4\)wherexix\_\{i\}is theii\-th component of𝒙\\boldsymbol\{x\}andxi\(l\)x\_\{i\}^\{\(l\)\}andxi\(u\)x\_\{i\}^\{\(u\)\}are the lower and upper bounds ofxix\_\{i\}, respectively\.

When solving PDEs with ELM, applying affine transformations is critical to achieving superior performance\. However, studies have demonstrated that omitting these transformations can lead to significantly degraded results\. An ELM incorporating affine transformations can be expressed as follows

uM​\(𝒙\)=𝒘T​σ​\(𝑾​𝒙~\+𝒃\)=𝒘T​σ​\(𝑾​𝑻​\(𝒙\)\+𝒃\),u\_\{M\}\(\\boldsymbol\{x\}\)=\\boldsymbol\{w\}^\{T\}\\sigma\(\\boldsymbol\{W\}\\tilde\{\\boldsymbol\{x\}\}\+\\boldsymbol\{b\}\)=\\boldsymbol\{w\}^\{T\}\\sigma\\left\(\\boldsymbol\{W\}\\boldsymbol\{T\}\(\\boldsymbol\{x\}\)\+\\boldsymbol\{b\}\\right\),\(5\)where𝑻=\[T1​\(x1\),T2​\(x2\),⋯,Td​\(xd\)\]T\\boldsymbol\{T\}=\\left\[T\_\{1\}\(x\_\{1\}\),T\_\{2\}\(x\_\{2\}\),\\cdots,T\_\{d\}\(x\_\{d\}\)\\right\]^\{T\}\.

## 3Methodologies

### 3\.1Fourier Feature Networks

It is observed that neural networks employingtanh\\tanhactivation function, while capable of generating acceptable results, fundamentally fail to attain machine\-level precision\. To address this limitation in approximation capacity and elevate solution accuracy, we introduce a novel architectural enhancement through Fourier feature embedding\.

#### 3\.1\.1Fourier Feature Network with acos\\cosactivation

Gallant and White\[[11](https://arxiv.org/html/2608.14733#bib.bib4)\]pioneered the integration of Fourier features into neural network architectures through their development of the cosine squasher activation function\. Building upon this foundation, Silvescu\[[32](https://arxiv.org/html/2608.14733#bib.bib5)\]implemented Fourier feature embedding via cosine\-based activation operators\. More recently, Ngom and Marin\[[26](https://arxiv.org/html/2608.14733#bib.bib3)\]introduced a Fourier Neural Network \(FNN\) with a single hidden layer activated by the cosine function\. Although their experiments demonstrated promising results, the accuracy of the solutions still falls short of machine precision\. Furthermore, their models omit the bias term, which can limit the expressive power of the neural network\.

Inspired by these works, we propose a Fourier feature network with acos\\cosactivation\. In this context, we also focus on a single\-hidden\-layer neural network\. The target function can be represented as follows

uM​\(𝒙\)=𝒘T​cos⁡\(𝑾​𝒙\+𝒃\)\.u\_\{M\}\(\\boldsymbol\{x\}\)=\\boldsymbol\{w\}^\{T\}\\cos\(\\boldsymbol\{W\}\\boldsymbol\{x\}\+\\boldsymbol\{b\}\)\.\(6\)It is noteworthy that, unlike networks activated by thesigmoid,tanh\\tanh, orswishfunctions \(as shown in Equation \([5](https://arxiv.org/html/2608.14733#S2.E5)\)\), the input to the Fourier feature network with a cos activation does not require affine transformations\.

#### 3\.1\.2Fourier Feature Network with asin\\sinactivation

Sitzmann et al\.\[[34](https://arxiv.org/html/2608.14733#bib.bib6)\]introduced Sinusoidal Representation Networks \(SIRENs\), which utilize‌ the sine function as a periodic activation function\. This design enables neural networks to accurately represent signals and their derivatives\. Building on this work, Li et al\. proposed an enhanced model called Spatially Collaged Coordinate Networks \(SCONE\), retaining‌ the sine activation function as a ‌foundational‌ component of the architecture\. Belbute\-Peres et al\.\[[8](https://arxiv.org/html/2608.14733#bib.bib7)\]also proposed an alternative improvement to SIRENs by incorporating a learnable scaling factor, enabling‌ automatic adjustment to inputs with ‌diverse‌ frequency characteristics\.

Motivated by these studies, and following the natural progression from cosine to sine activation, we propose a Fourier feature network with a sin activation\. As before, we focus on a single\-hidden\-layer neural network\. The function to be learned is represented as follows

uM​\(𝒙\)=𝒘T​sin⁡\(𝑾​𝒙\+𝒃\)\.u\_\{M\}\(\\boldsymbol\{x\}\)=\\boldsymbol\{w\}^\{T\}\\sin\(\\boldsymbol\{W\}\\boldsymbol\{x\}\+\\boldsymbol\{b\}\)\.\(7\)Similar to the Fourier feature network with acos\\cosactivation, the input to the network with asin\\sinactivation does not require affine transformations\.

#### 3\.1\.3Fourier Feature Network withcos\\cosandsin\\sinactivations

In 2013, Liu\[[24](https://arxiv.org/html/2608.14733#bib.bib8)\]proposed a Fourier neural network activated by both cosine and sine functions\. However, these networks were primarily applied to regression and classification tasks, and their accuracy was not satisfactory\. More recent studies have incorporated Fourier features into neural networks by transforming input variables using cosine and sine functions, thereby embedding these features directly into the network’s input layer\[[35](https://arxiv.org/html/2608.14733#bib.bib9),[37](https://arxiv.org/html/2608.14733#bib.bib10),[22](https://arxiv.org/html/2608.14733#bib.bib11)\]\. In line with Liu’s work, Fourier PINNs\[[6](https://arxiv.org/html/2608.14733#bib.bib12)\]also employ both cosine and sine activations as part of their basis functions\. The crucial difference lies in how these basis functions are combined within the network to represent the learned function\.

Inspired by these approaches, we propose a Fourier feature network utilizing bothcos\\cosandsin\\sinactivations\. The output function is defined as

uM​\(𝒙\)=\(𝒘\(1\)\)T​cos⁡\(𝑾\(1\)​𝒙\+𝒃\(1\)\)\+\(𝒘\(2\)\)T​sin⁡\(𝑾\(2\)​𝒙\+𝒃\(2\)\),u\_\{M\}\(\\boldsymbol\{x\}\)=\(\\boldsymbol\{w\}^\{\(1\)\}\)^\{T\}\\cos\(\\boldsymbol\{W\}^\{\(1\)\}\\boldsymbol\{x\}\+\\boldsymbol\{b\}^\{\(1\)\}\)\+\(\\boldsymbol\{w\}^\{\(2\)\}\)^\{T\}\\sin\(\\boldsymbol\{W\}^\{\(2\)\}\\boldsymbol\{x\}\+\\boldsymbol\{b\}^\{\(2\)\}\),\(8\)where𝑾\(1\)∈ℝp×d\\boldsymbol\{W\}^\{\(1\)\}\\in\\mathbb\{R\}^\{p\\times d\},𝑾\(2\)∈ℝp×d\\boldsymbol\{W\}^\{\(2\)\}\\in\\mathbb\{R\}^\{p\\times d\},𝒃\(1\)∈ℝp\\boldsymbol\{b\}^\{\(1\)\}\\in\\mathbb\{R\}^\{p\},𝒃\(2\)∈ℝp\\boldsymbol\{b\}^\{\(2\)\}\\in\\mathbb\{R\}^\{p\},𝒘\(1\)∈ℝp\\boldsymbol\{w\}^\{\(1\)\}\\in\\mathbb\{R\}^\{p\},𝒘\(2\)∈ℝp\\boldsymbol\{w\}^\{\(2\)\}\\in\\mathbb\{R\}^\{p\}andM=2​pM=2pis the number of basis functions\.

### 3\.2Function approximation

Let us assume that a real\-valued matrix𝐀\\mathbf\{A\}of sizeN×MN\\times Mis constructed from the basis functionsΦ=\[ϕ1,ϕ2,⋯,ϕM\]\\Phi=\[\\phi\_\{1\},\\phi\_\{2\},\\cdots,\\phi\_\{M\}\], evaluated at discrete points within the datasetSS\. This matrix is defined by

𝐀=\[ϕ1​\(𝒙1\)ϕ2​\(𝒙1\)⋯ϕM​\(𝒙1\)ϕ1​\(𝒙2\)ϕ2​\(𝒙2\)⋯ϕM​\(𝒙2\)⋯ϕ1​\(𝒙N\)ϕ2​\(𝒙N\)⋯ϕM​\(𝒙N\)\],\\mathbf\{A\}=\\begin\{bmatrix\}\\phi\_\{1\}\(\\boldsymbol\{x\}\_\{1\}\)&\\phi\_\{2\}\(\\boldsymbol\{x\}\_\{1\}\)&\\cdots&\\phi\_\{M\}\(\\boldsymbol\{x\}\_\{1\}\)\\\\ \\phi\_\{1\}\(\\boldsymbol\{x\}\_\{2\}\)&\\phi\_\{2\}\(\\boldsymbol\{x\}\_\{2\}\)&\\cdots&\\phi\_\{M\}\(\\boldsymbol\{x\}\_\{2\}\)\\\\ \\vdots&\\vdots&\\cdots&\\vdots\\\\ \\phi\_\{1\}\(\\boldsymbol\{x\}\_\{N\}\)&\\phi\_\{2\}\(\\boldsymbol\{x\}\_\{N\}\)&\\cdots&\\phi\_\{M\}\(\\boldsymbol\{x\}\_\{N\}\)\\\\ \\end\{bmatrix\},\(9\)where𝒙i∈S\\boldsymbol\{x\}\_\{i\}\\in Sand1≤i≤N1\\leq i\\leq N\.

Based on Equation \([2](https://arxiv.org/html/2608.14733#S2.E2)\), the coefficient vector𝒘=\[w1,w2,…,wM\]T\\boldsymbol\{w\}=\[w\_\{1\},w\_\{2\},\.\.\.,w\_\{M\}\]^\{T\}can be determined by solving the linear system𝐀​𝒘=𝐅\\mathbf\{A\}\\boldsymbol\{w\}=\\mathbf\{F\}using the least\-squares method\. Here,𝐀=Φ⁡\(S\)\\mathbf\{A\}=\\Phi\(S\)denotes the basis function matrix formed by evaluating all basis functions at discrete data points in the datasetSS, and𝐅=\[f⁡\(𝒙1\),f⁡\(𝒙2\),…,f⁡\(𝒙N\)\]T\\mathbf\{F\}=\[f\(\\boldsymbol\{x\}\_\{1\}\),f\(\\boldsymbol\{x\}\_\{2\}\),\\ldots,f\(\\boldsymbol\{x\}\_\{N\}\)\]^\{T\}corresponds to the vector of function values at these sampled points\.

### 3\.3Searching for optimal scaling factor

The most widely used initialization methods in neural networks are Xavier initialization \(also known as Glorot initialization\)\[[12](https://arxiv.org/html/2608.14733#bib.bib13)\]and Kaiming initialization\[[13](https://arxiv.org/html/2608.14733#bib.bib14)\]\. These techniques are designed to initialize neural network weights according to specific distributions with carefully controlled variances\. In conventional practice, both weight and bias parameters are typically assumed to originate from identical distributions\. Consequently, biases are generally initialized using the same methodology as weights\.

In this work, we employ a uniform distribution with unit variance to initialize both weight𝑾\\boldsymbol\{W\}and bias𝒃\\boldsymbol\{b\}parameters\. This choice is motivated by practical observations: neither Xavier nor Kaiming initialization consistently provides optimal results for our use cases\. Crucially, the primary difference between these two methods lies in the scaling factor applied to weights and biases\. To address this limitation, we implement a systematic search protocol across a constrained scaling factor domain to identify initialization\-sensitive optimal configurations\.

Let the scaling factor be denoted byρ\\rho\. The functional mapping described by the neural network architecture in Equation \([3](https://arxiv.org/html/2608.14733#S2.E3)\) can then be reparameterized in terms ofρ\\rhoas

uM​\(𝒙\)=𝒘T​σ​\(ρ⁡\(𝑾​𝒙\+𝒃\)\)\.u\_\{M\}\(\\boldsymbol\{x\}\)=\\boldsymbol\{w\}^\{T\}\\sigma\(\\rho\(\\boldsymbol\{W\}\\boldsymbol\{x\}\+\\boldsymbol\{b\}\)\)\.\(10\)The protocol for locating the optimal scaling factorρo​p​t\\rho\_\{opt\}is formally specified in Algorithm[1](https://arxiv.org/html/2608.14733#algorithm1)\. It is important to note that although the affine transformation is not explicitly shown in Equation \([10](https://arxiv.org/html/2608.14733#S3.E10)\), it must still be applied when usingsigmoid,tanh\\tanh, andswishas activation functions\.

\\chadded

It is important to note that the optimal scaling factor search method proposed in\[[10](https://arxiv.org/html/2608.14733#bib.bib40)\]used a differential evolution algorithm, whereas in our work, we enumerate all candidate scaling factors within a given range using a step\-size approach\. The scaling factor that minimizes the error is then selected from these candidates\. This method is simple and easy to implement\. Moreover, the work\[[10](https://arxiv.org/html/2608.14733#bib.bib40)\]focused exclusively on Gaussian activation functions, and the resulting numerical accuracy was not optimal in some cases\. We argue that incorporating Fourier features is more essential ‌for achieving better performance\.

Give the range of

ρ\\rhoas

\(ρm​i​n,ρm​a​x\]\(\\rho\_\{min\},\\rho\_\{max\}\]and the step size

ρs\\rho\_\{s\}\. Let

MMbe the number of basis functions and

ρ=ρm​i​n\+ρs\\rho=\\rho\_\{min\}\+\\rho\_\{s\}\. Randomly initialize the weights of the networks\. while*ρ<ρm​a​x\\rho<\\rho\_\{max\}*do

Obtain the basis functions

Φ\\Phi\. Compute the matrix

𝐀\\mathbf\{A\}\. Construct the right\-hand vector

𝐅\\mathbf\{F\}\. Express the current approximation of the solution as

u⁡\(𝒙\)=Φ⋅𝒘u\(\\boldsymbol\{x\}\)=\\Phi\\cdot\\boldsymbol\{w\}\. Solve the linear system

𝐀⋅𝒘=𝐅\\mathbf\{A\}\\cdot\\boldsymbol\{w\}=\\mathbf\{F\}using least squares method to obtain the vector of coefficients

𝒘\\boldsymbol\{w\}\. Let

L​o​s​s=∥𝐀⋅𝒘−𝐅∥Loss=\\lVert\\mathbf\{A\}\\cdot\\boldsymbol\{w\}\-\\mathbf\{F\}\\rVertand record the

L​o​s​sLossand the current

ρ\\rhovalue\. Let

ρ=ρm​i​n\+ρs\\rho=\\rho\_\{min\}\+\\rho\_\{s\}\. end while

The optimal scale factor

ρo​p​t\\rho\_\{opt\}is the one with the smallest

L​o​s​sLoss\.

Algorithm 1Searching for optimal scaling factor###### Theorem 3\.1\.

\(Theorem 2\.2\[[29](https://arxiv.org/html/2608.14733#bib.bib29)\]\) Givenp≥1,s,k,d∈N\+p\\geq 1,s,k,d\\in N^\{\+\},s≥k\+1s\\geq k\+1\. Letσ\\sigmabe the logistic function ortanh\\tanhfunction\. Denoteℱs,p,d:=\{u∈𝒲s,p​\(\[0,1\]d\):∥u∥𝒲s,p​\(\[0,1\]d\)≤1\}\\mathscr\{F\}\_\{s,p,d\}:=\\\{u\\in\\mathcal\{W\}^\{s,p\}\(\[0,1\]^\{d\}\):\\lVert u\\rVert\_\{\\mathcal\{W\}^\{s,p\}\(\[0,1\]^\{d\}\)\}\\leq 1\\\}, and𝒩σ\(MD,BD\):=\{uM\(x\)=𝐰Tσ\(𝐖𝐱\+𝐛\)withND≤MD,and\|wi​j\|≤BD,\|bi\|≤BD\}\\mathscr\{N\}\_\{\\sigma\}\(M\_\{D\},B\_\{D\}\):=\\\{u\_\{M\}\(x\)=\\boldsymbol\{w\}^\{T\}\\sigma\(\\boldsymbol\{W\}\\boldsymbol\{x\}\+\\boldsymbol\{b\}\)\\text\{ with \}N\_\{D\}\\leq M\_\{D\},\\text\{ and \}\\lvert w\_\{ij\}\\rvert\\leq B\_\{D\},\\lvert b\_\{i\}\\rvert\\leq B\_\{D\}\\\}, whereNDN\_\{D\}is the number of non\-zero parameter elements in the hidden layer of the neural network\. For anyϵ\\epsilon¿ 0 and u∈ℱs,p,d\\in\\mathscr\{F\}\_\{s,p,d\}, there exists a neural networkuMu\_\{M\}∈𝒩σ​\(MD,BD\)\\in\\mathscr\{N\}\_\{\\sigma\}\(M\_\{D\},B\_\{D\}\)withMD≤C⋅ϵ−d/s−k−μkM\_\{D\}\\leq C\\cdot\\epsilon^\{\-d/s\-k\-\\mu k\}, andBD≤C⋅ϵ−θB\_\{D\}\\leq C\\cdot\\epsilon^\{\-\\theta\}such that

∥u−uM∥𝒲s,p​\(\[0,1\]d\)≤ϵ,\\parallel u\-u\_\{M\}\\parallel\_\{\\mathcal\{W\}^\{s,p\}\(\[0,1\]^\{d\}\)\}\\leq\\epsilon,\(11\)whereC,θC,\\thetaare constants depending ond,s,p,kd,s,p,k;μ\\muis an arbitrarily small positive number\.

In Theorem[3\.1](https://arxiv.org/html/2608.14733#S3.Thmtheorem1), it is established that for any desired approximation accuracyϵ\\epsilon, there exist suitable neural network parameters that allow the constructed network to approximate the target function within this tolerance\. However, the challenge lies in the fact that these optimal parameters are generally unknown in practice\. In the RNN\-PG framework\[[29](https://arxiv.org/html/2608.14733#bib.bib29)\], these parameters are typically obtained through repeated random initialization, which are both computationally inefficient and potentially suboptimal\.

Theorem[3\.2](https://arxiv.org/html/2608.14733#S3.Thmtheorem2)demonstrates that, for any given initialization, it is possible to improve the applicability of the network parameters by appropriately adjusting the scaling factorρ\\rho\.

###### Theorem 3\.2\.

\(Existence of the optimal scaling factor\)Assume thatu⁡\(𝐱\)u\(\\boldsymbol\{x\}\)is a continuous function andσ\\sigmais a continuous activation function\. LetuM​\(𝐱\)=𝐰T​\(ρ\)​σ​\(ρ⁡\(𝐖​𝐱\+𝐛\)\)u\_\{M\}\(\\boldsymbol\{x\}\)=\\boldsymbol\{w\}^\{T\}\(\\rho\)\\sigma\(\\rho\(\\boldsymbol\{W\}\\boldsymbol\{x\}\+\\boldsymbol\{b\}\)\)anduMu\_\{M\}is the projection ofuuonto the linear space spanned by the basis functionsσ⁡\(ρ⁡\(𝐖​𝐱\+𝐛\)\)\\sigma\(\\rho\(\\boldsymbol\{W\}\\boldsymbol\{x\}\+\\boldsymbol\{b\}\)\), where𝐰⁡\(ρ\)\\boldsymbol\{w\}\(\\rho\)is a vector consisting of a set of functions related toρ\\rho\. Then,∀ρm​a​x\>0\\forall\\rho\_\{max\}\>0,∃ρo​p​t∈\[0,ρm​a​x\]\\exists\\rho\_\{opt\}\\in\[0,\\rho\_\{max\}\]such that

∥u−uM​\(𝒙,ρo​p​t\)∥=min0≤ρ≤ρm​a​x⁡∥u⁡\(𝒙\)−uM​\(𝒙,ρ\)∥\.\\lVert u\-u\_\{M\}\(\\boldsymbol\{x\};\\rho\_\{opt\}\)\\rVert=\\min\_\{0\\leq\\rho\\leq\\rho\_\{max\}\}\\lVert u\(\\boldsymbol\{x\}\)\-u\_\{M\}\(\\boldsymbol\{x\};\\rho\)\\rVert\.

###### Proof\.

LetΦ\\Phibe the set of basis functionsσ⁡\(ρ⁡\(𝑾​𝒙\+𝒃\)\)\\sigma\(\\rho\(\\boldsymbol\{W\}\\boldsymbol\{x\}\+\\boldsymbol\{b\}\)\), andS=\[𝒙1,𝒙2,⋯,𝒙N\]TS=\[\\boldsymbol\{x\}\_\{1\},\\boldsymbol\{x\}\_\{2\},\\cdots,\\boldsymbol\{x\}\_\{N\}\]^\{T\}be the set of collocation points in the domainΩ\\Omega\. According to Equation \([9](https://arxiv.org/html/2608.14733#S3.E9)\), we substitute the set of collocation pointsSSinto it to obtain

𝐀⁡\(ρ\)​𝒘​\(ρ\)=𝐅,\\mathbf\{A\}\(\\rho\)\\boldsymbol\{w\}\(\\rho\)=\\mathbf\{F\},where𝐅=\[u⁡\(𝒙1\),u⁡\(𝒙2\),⋯,u⁡\(𝒙N\)\]T\\mathbf\{F\}=\[u\(\\boldsymbol\{x\}\_\{1\}\),u\(\\boldsymbol\{x\}\_\{2\}\),\\cdots,u\(\\boldsymbol\{x\}\_\{N\}\)\]^\{T\}\.

The continuity of the activation functionσ\\sigmaimplies continuity of𝐀⁡\(ρ\)\\mathbf\{A\}\(\\rho\)inρ\\rho\. Given𝐅\\mathbf\{F\}is likewise continuous inρ\\rho, it follows that both𝒘⁡\(ρ\)\\boldsymbol\{w\}\(\\rho\)and ultimatelyuM​\(𝒙,ρ\)u\_\{M\}\(\\boldsymbol\{x\};\\rho\)are continuous inρ\\rho\. Lete⁡\(𝒙,ρ\)=u⁡\(𝒙\)−uM​\(𝒙,ρ\)e\(\\boldsymbol\{x\};\\rho\)=u\(\\boldsymbol\{x\}\)\-u\_\{M\}\(\\boldsymbol\{x\};\\rho\), wheree⁡\(𝒙,ρ\)e\(\\boldsymbol\{x\};\\rho\)is continuous inρ\\rho\. By the continuity of the norm,∥e⁡\(𝒙,ρ\)∥\\lVert e\(\\boldsymbol\{x\};\\rho\)\\rVertis likewise continuous inρ\\rho\. Consequently,∀ρm​a​x\>0\\forall\\rho\_\{max\}\>0,∥e⁡\(𝒙,ρ\)∥\\lVert e\(\\boldsymbol\{x\};\\rho\)\\rVertattains its maximum and minimum values on the closed interval\[0,ρm​a​x\]\[0,\\rho\_\{max\}\]\. ∎

Figure[1](https://arxiv.org/html/2608.14733#S3.F1)shows how theL∞L\_\{\\infty\}error evolves with the scaling factor across different basis function counts\. The ELM with thetanh\\tanhactivation function is employed to solve the two\-dimensional Equation \([25](https://arxiv.org/html/2608.14733#S4.E25)\)\. The scaling factorρ\\rhois varied within the range\(0,10\]\(0,10\]with a step size of0\.010\.01\. As shown in the figure, increasing the number of basis functions enhances the representational capacity of the ELM, leading to a significant reduction in error\. Furthermore, while the error fluctuates as the scaling factorρ\\rhochanges, reveals the existence of an optimal scaling factor \(denoted asρo​p​t\\rho\_\{opt\}\) that minimizes the error\. This observation numerically validates Theorem[3\.2](https://arxiv.org/html/2608.14733#S3.Thmtheorem2)\.

\\chadded

It should be noted that although Figure[1](https://arxiv.org/html/2608.14733#S3.F1)displays only the error curves of the ELM withtanh\\tanhactivation function under varying scaling factorρ\\rho, a similar trend occurs in FENs\. Moreover, FENs typically operate across a broaderρ\\rhosearch range and achieve lower errors\.

![Refer to caption](https://arxiv.org/html/2608.14733v1/figures/2D_loss_infty_scale.png)

Figure 1:Error variation with the scaling factorρ\\rhofor different numbers of basis functions\.
### 3\.4Solving linear PDEs

Consider the following generic linear PDE

ℒ​u=f,in​Ω,\\displaystyle\\mathcal\{L\}u=f,\\ \\mbox\{in\}\\kern 5\.0pt\\Omega,\(12\)ℬ​u=g,on​∂Ω,\\displaystyle\\mathcal\{B\}u=g,\\ \\mbox\{on\}\\kern 5\.0pt\\partial\\Omega,\(13\)whereuuis the scalar field function to be approximated,ℒ\\mathcal\{L\}is a linear PDE operator defined within the domainΩ\\Omegaandℬ\\mathcal\{B\}is a linear boundary operator acting on the boundary∂Ω\\partial\\Omega\. The functionsffandggrepresent the source term and the boundary condition, respectively\.

LetΦ\\Phidenote the basis functions generated by the outputs of the FENs\. Applying the operatorsℒ\\mathcal\{L\}andℬ\\mathcal\{B\}to these basis functions yields

ℒ​Φ=\(ℒ​ϕ1,ℒ​ϕ2,⋯,ℒ​ϕM\),\\displaystyle\\mathcal\{L\}\\Phi=\(\\mathcal\{L\}\\phi\_\{1\},\\mathcal\{L\}\\phi\_\{2\},\\cdots,\\mathcal\{L\}\\phi\_\{M\}\),\(14\)ℬ​Φ=\(ℬ​ϕ1,ℬ​ϕ2,⋯,ℬ​ϕM\)\.\\displaystyle\\mathcal\{B\}\\Phi=\(\\mathcal\{B\}\\phi\_\{1\},\\mathcal\{B\}\\phi\_\{2\},\\cdots,\\mathcal\{B\}\\phi\_\{M\}\)\.\(15\)Assume that the collocation datasetSSconsists ofN=Nr\+NbN=N\_\{r\}\+N\_\{b\}points, withNrN\_\{r\}points distributed within the domainΩ\\OmegaandNbN\_\{b\}points allocated on the boundary∂Ω\\partial\\Omega\. These points are partitioned into two subsets:SrS\_\{r\}\(interior points\) andSbS\_\{b\}\(boundary points\)\. Under this configuration, the matrix𝐀\\mathbf\{A\}, which corresponds to the left\-hand side of the PDE system, is constructed as follows

𝐀\\displaystyle\\mathbf\{A\}=\[ℒ​Φ​\(Sr\)ℬ​Φ​\(Sb\)\]\\displaystyle=\\begin\{bmatrix\}\\mathcal\{L\}\\Phi\(S\_\{r\}\)\\\\ \\mathcal\{B\}\\Phi\(S\_\{b\}\)\\end\{bmatrix\}\(16\)=\[ℒ​ϕ1​\(𝒙1\)ℒ​ϕ2​\(𝒙1\)⋯ℒ​ϕM​\(𝒙1\)⋯ℒ​ϕ1​\(𝒙Nr\)ℒ​ϕ2​\(𝒙Nr\)⋯ℒ​ϕM​\(𝒙Nr\)ℬ​ϕ1​\(𝒙Nr\+1\)ℬ​ϕ2​\(𝒙Nr\+1\)⋯ℬ​ϕM​\(𝒙Nr\+1\)⋯ℬ​ϕ1​\(𝒙Nr\+Nb\)ℬ​ϕ2​\(𝒙Nr\+Nb\)⋯ℬ​ϕM​\(𝒙Nr\+Nb\)\]\.\\displaystyle=\\begin\{bmatrix\}\\mathcal\{L\}\\phi\_\{1\}\(\\boldsymbol\{x\}\_\{1\}\)&\\mathcal\{L\}\\phi\_\{2\}\(\\boldsymbol\{x\}\_\{1\}\)&\\cdots&\\mathcal\{L\}\\phi\_\{M\}\(\\boldsymbol\{x\}\_\{1\}\)\\\\ \\vdots&\\vdots&\\cdots&\\vdots\\\\ \\mathcal\{L\}\\phi\_\{1\}\(\\boldsymbol\{x\}\_\{N\_\{r\}\}\)&\\mathcal\{L\}\\phi\_\{2\}\(\\boldsymbol\{x\}\_\{N\_\{r\}\}\)&\\cdots&\\mathcal\{L\}\\phi\_\{M\}\(\\boldsymbol\{x\}\_\{N\_\{r\}\}\)\\\\ \\mathcal\{B\}\\phi\_\{1\}\(\\boldsymbol\{x\}\_\{N\_\{r\}\+1\}\)&\\mathcal\{B\}\\phi\_\{2\}\(\\boldsymbol\{x\}\_\{N\_\{r\}\+1\}\)&\\cdots&\\mathcal\{B\}\\phi\_\{M\}\(\\boldsymbol\{x\}\_\{N\_\{r\}\+1\}\)\\\\ \\vdots&\\vdots&\\cdots&\\vdots\\\\ \\mathcal\{B\}\\phi\_\{1\}\(\\boldsymbol\{x\}\_\{N\_\{r\}\+N\_\{b\}\}\)&\\mathcal\{B\}\\phi\_\{2\}\(\\boldsymbol\{x\}\_\{N\_\{r\}\+N\_\{b\}\}\)&\\cdots&\\mathcal\{B\}\\phi\_\{M\}\(\\boldsymbol\{x\}\_\{N\_\{r\}\+N\_\{b\}\}\)\\\\ \\end\{bmatrix\}\.
As described in the previous section, the coefficient vector𝒘\\boldsymbol\{w\}is obtained by solving the linear system𝐀​𝒘=𝐅\\mathbf\{A\}\\boldsymbol\{w\}=\\mathbf\{F\}, where the right\-hand side vector𝐅\\mathbf\{F\}is given by

𝑭=\[f⁡\(𝒙1\),⋯,f⁡\(𝒙Nr\),g⁡\(𝒙Nr\+1\),⋯,g⁡\(𝒙Nr\+Nb\)\]T\.\\boldsymbol\{F\}=\[f\(\\boldsymbol\{x\}\_\{1\}\),\\cdots,f\(\\boldsymbol\{x\}\_\{N\_\{r\}\}\),g\(\\boldsymbol\{x\}\_\{N\_\{r\}\+1\}\),\\cdots,g\(\\boldsymbol\{x\}\_\{N\_\{r\}\+N\_\{b\}\}\)\]^\{T\}\.\(17\)Here,f⁡\(𝒙i\)f\(\\boldsymbol\{x\}\_\{i\}\)for1≤i≤Nr1\\leq i\\leq N\_\{r\}corresponds to the source term evaluated at the ‌interior collocation points, whereasg⁡\(𝒙i\)g\(\\boldsymbol\{x\}\_\{i\}\)forNr\+1≤i≤Nr\+NbN\_\{r\}\+1\\leq i\\leq N\_\{r\}\+N\_\{b\}represents the boundary conditions imposed at the ‌boundary collocation points\.\\chaddedIt is worth noting that, for time\-dependent problems, we adopt a unified treatment of both initial and boundary conditions by handling them as boundary constraints\. This approach simplifies the overall formulation\. Additionally, we do not explicitly assign different weights to the governing equations and the initial or boundary conditions; instead, all constraints are treated with equal importance, with an implicit weight of 1\. The training and test sets are chosen to be identical, as our objective is to compute solution values on the given dataset\. This setup aligns with the traditional goal in computational mathematics, which is to approximate the solution at prescribed or unknown locations\.

### 3\.5Approxmation theory of FENs

Let us consider functions defined onℝd\\mathbb\{R\}^\{d\}that admit the following Fourier integral representation

u⁡\(𝒙\)=∫ℝde𝕚​𝑾r​𝒙​F~​\(d​𝑾r\),u\(\\boldsymbol\{x\}\)=\\int\_\{\\mathbb\{R\}^\{d\}\}e^\{\\mathbb\{i\}\\boldsymbol\{W\}\_\{r\}\\boldsymbol\{x\}\}\\tilde\{F\}\(\\,\\mathrm\{d\}\{\\boldsymbol\{W\}\_\{r\}\}\),\(18\)whereF~​\(d​𝑾r\)=ei​θ​\(𝑾r\)​F​\(d​𝑾r\)\\tilde\{F\}\(\\,\\mathrm\{d\}\{\\boldsymbol\{W\}\_\{r\}\}\)=e^\{i\\theta\(\\boldsymbol\{W\}\_\{r\}\)\}F\(\\,\\mathrm\{d\}\{\\boldsymbol\{W\}\_\{r\}\}\)is a unique complex\-valued measure \(referred to as the Fourier distribution\), withF⁡\(d​𝑾r\)F\(\\,\\mathrm\{d\}\{\\boldsymbol\{W\}\_\{r\}\}\)representing the magnitude distribution,θ⁡\(𝑾r\)\\theta\(\\boldsymbol\{W\}\_\{r\}\)denoting the phase, and𝑾r∈ℝ1×d\\boldsymbol\{W\}\_\{r\}\\in\\mathbb\{R\}^\{1\\times d\}\. The class𝔽\\mathbb\{F\}of functions that admit such a representation with a finite total magnitude measure is defined as

𝔽=\{u:ℝd→ℝ\|∫ℝdF\(d𝑾r\)<∞\}\.\\mathbb\{F\}=\\left\\\{u:\\mathbb\{R\}^\{d\}\\to\\mathbb\{R\}\\big\|\\int\_\{\\mathbb\{R\}^\{d\}\}F\(\\,\\mathrm\{d\}\{\\boldsymbol\{W\}\_\{r\}\}\)<\\infty\\right\\\}\.\(19\)For anyC\>0C\>0, the class𝔽C\\mathbb\{F\}\_\{C\}is defined as

𝔽C=\{u:ℝd→ℝ\|Cu=∫ℝdF\(d𝑾r\)<C\}\.\\mathbb\{F\}\_\{C\}=\\left\\\{u:\\mathbb\{R\}^\{d\}\\to\\mathbb\{R\}\\big\|C\_\{u\}=\\int\_\{\\mathbb\{R\}^\{d\}\}F\(\\,\\mathrm\{d\}\{\\boldsymbol\{W\}\_\{r\}\}\)<C\\right\\\}\.\(20\)We define the class𝔽Ω\\mathbb\{F\}\_\{\\Omega\}as the set of functions defined on a bounded domainΩ\\Omegathat admit a Fourier representation of the form \([18](https://arxiv.org/html/2608.14733#S3.E18)\) \(valid for𝒙∈Ω\\boldsymbol\{x\}\\in\\Omega\) with finite magnitude measure

𝔽Ω=\{u:Ω→ℝ\|∫ℝdF\(d𝑾r\)<∞\}\.\\mathbb\{F\}\_\{\\Omega\}=\\left\\\{u:\\Omega\\to\\mathbb\{R\}\\big\|\\int\_\{\\mathbb\{R\}^\{d\}\}F\(\\,\\mathrm\{d\}\{\\boldsymbol\{W\}\_\{r\}\}\)<\\infty\\right\\\}\.\(21\)For eachC\>0C\>0, we define

𝔽C,Ω=\{u:Ω→ℝ\|Cu,Ω=∫ℝdF\(d𝑾r\)≤C\}\.\\mathbb\{F\}\_\{C,\\Omega\}=\\left\\\{u:\\Omega\\to\\mathbb\{R\}\\big\|C\_\{u,\\Omega\}=\\int\_\{\\mathbb\{R\}^\{d\}\}F\(\\,\\mathrm\{d\}\{\\boldsymbol\{W\}\_\{r\}\}\)\\leq C\\right\\\}\.\(22\)
###### Theorem 3\.3\.

\(Approximation error of FEN with acos\\cosactivation\[[40](https://arxiv.org/html/2608.14733#bib.bib39)\]\)LetuM=∑i=1Mwi​cos⁡\(𝐖i​𝐱\+bi\)u\_\{M\}=\\sum\_\{i=1\}^\{M\}w\_\{i\}\\cos\(\\boldsymbol\{W\}\_\{i\}\\boldsymbol\{x\}\+b\_\{i\}\)for𝐖i∈ℝ1×d\\boldsymbol\{W\}\_\{i\}\\in\\mathbb\{R\}^\{1\\times d\},bi∈ℝb\_\{i\}\\in\\mathbb\{R\},\|wi\|≤CM\\lvert w\_\{i\}\\rvert\\leq\\frac\{C\}\{M\}, and𝐱∈ℝd\\boldsymbol\{x\}\\in\\mathbb\{R\}^\{d\}\. For each functionu∈𝔽C,Ωu\\in\\mathbb\{F\}\_\{C,\\Omega\}and any probability measureμ\\mu, there existsuMu\_\{M\}\(M≥1M\\geq 1\), such that

∫Ω\|u⁡\(𝒙\)−uM​\(𝒙\)\|2​μ​\(𝑑𝒙\)≤C2M\.\\int\_\{\\Omega\}\\lvert u\(\\boldsymbol\{x\}\)\-u\_\{M\}\(\\boldsymbol\{x\}\)\\rvert^\{2\}\\mu\(\\,\\mathrm\{d\}\\boldsymbol\{x\}\)\\leq\\frac\{C^\{2\}\}\{M\}\.

###### Theorem 3\.4\.

\(Approximation error of FEN with asin\\sinactivation\)LetuM=∑i=1Mwi​sin⁡\(𝐖i​𝐱\+bi\)u\_\{M\}=\\sum\_\{i=1\}^\{M\}w\_\{i\}\\sin\(\\boldsymbol\{W\}\_\{i\}\\boldsymbol\{x\}\+b\_\{i\}\)for𝐖i∈ℝ1×d\\boldsymbol\{W\}\_\{i\}\\in\\mathbb\{R\}^\{1\\times d\},bi∈ℝb\_\{i\}\\in\\mathbb\{R\},\|wi\|≤CM\\lvert w\_\{i\}\\rvert\\leq\\frac\{C\}\{M\}and𝐱∈ℝd\\boldsymbol\{x\}\\in\\mathbb\{R\}^\{d\}\. For each functionu∈𝔽C,Ωu\\in\\mathbb\{F\}\_\{C,\\Omega\}and any probability measureμ\\mu, there existsuMu\_\{M\}\(M≥1M\\geq 1\), such that

∫Ω\|u⁡\(𝒙\)−uM​\(𝒙\)\|2​μ​\(𝑑𝒙\)≤C2M\.\\int\_\{\\Omega\}\\lvert u\(\\boldsymbol\{x\}\)\-u\_\{M\}\(\\boldsymbol\{x\}\)\\rvert^\{2\}\\mu\(\\,\\mathrm\{d\}\\boldsymbol\{x\}\)\\leq\\frac\{C^\{2\}\}\{M\}\.

According to Theorem[3\.3](https://arxiv.org/html/2608.14733#S3.Thmtheorem3), the proof of Theorem[3\.4](https://arxiv.org/html/2608.14733#S3.Thmtheorem4)is straightforward\.

###### Theorem 3\.5\.

\(Approxmation error of FEN withcos\\cosandsin\\sinactivations\)Let

uM=∑i=1pwi\(1\)​cos⁡\(𝑾i\(1\)​𝒙\+bi\(1\)\)\+∑i=1pwi\(2\)​sin⁡\(𝑾i\(2\)​𝒙\+bi\(2\)\)u\_\{M\}=\\sum\_\{i=1\}^\{p\}w\_\{i\}^\{\(1\)\}\\cos\(\\boldsymbol\{W\}\_\{i\}^\{\(1\)\}\\boldsymbol\{x\}\+b\_\{i\}^\{\(1\)\}\)\+\\sum\_\{i=1\}^\{p\}w\_\{i\}^\{\(2\)\}\\sin\(\\boldsymbol\{W\}\_\{i\}^\{\(2\)\}\\boldsymbol\{x\}\+b\_\{i\}^\{\(2\)\}\)for𝐖i\(1\)∈ℝ1×d\\boldsymbol\{W\}\_\{i\}^\{\(1\)\}\\in\\mathbb\{R\}^\{1\\times d\},𝐖i\(2\)∈ℝ1×d\\boldsymbol\{W\}\_\{i\}^\{\(2\)\}\\in\\mathbb\{R\}^\{1\\times d\},bi\(1\)∈ℝb\_\{i\}^\{\(1\)\}\\in\\mathbb\{R\},bi\(2\)∈ℝb\_\{i\}^\{\(2\)\}\\in\\mathbb\{R\},\|wi\(1\)\|≤Cp\\lvert w\_\{i\}^\{\(1\)\}\\rvert\\leq\\frac\{C\}\{p\},\|wi\(2\)\|≤Cp\\lvert w\_\{i\}^\{\(2\)\}\\rvert\\leq\\frac\{C\}\{p\}, and𝐱∈ℝd\\boldsymbol\{x\}\\in\\mathbb\{R\}^\{d\}\. For each functionu∈𝔽C,Ωu\\in\\mathbb\{F\}\_\{C,\\Omega\}and any probability measureμ\\mu, there existsuMu\_\{M\}\(M=2​p≥2M=2p\\geq 2\), such that

∫Ω\|u⁡\(𝒙\)−uM​\(𝒙\)\|2​μ​\(𝑑𝒙\)≤4​C2p\.\\int\_\{\\Omega\}\\lvert u\(\\boldsymbol\{x\}\)\-u\_\{M\}\(\\boldsymbol\{x\}\)\\rvert^\{2\}\\mu\(\\,\\mathrm\{d\}\\boldsymbol\{x\}\)\\leq\\frac\{4C^\{2\}\}\{p\}\.

###### Proof\.

From Theorem[3\.3](https://arxiv.org/html/2608.14733#S3.Thmtheorem3), for any given functionu∈𝔽C,Ωu\\in\\mathbb\{F\}\_\{C,\\Omega\}, there is a sumup\(1\)=∑i=1pwi\(1\)​cos⁡\(𝑾i\(1\)​𝒙\+bi\(1\)\)u\_\{p\}^\{\(1\)\}=\\sum\_\{i=1\}^\{p\}w\_\{i\}^\{\(1\)\}\\cos\(\\boldsymbol\{W\}\_\{i\}^\{\(1\)\}\\boldsymbol\{x\}\+b\_\{i\}^\{\(1\)\}\), such that

∫Ω\|12​u​\(𝒙\)−up\(1\)​\(𝒙\)\|2​μ​\(𝑑𝒙\)≤C2p\.\\int\_\{\\Omega\}\\lvert\\frac\{1\}\{2\}u\(\\boldsymbol\{x\}\)\-u\_\{p\}^\{\(1\)\}\(\\boldsymbol\{x\}\)\\rvert^\{2\}\\mu\(\\,\\mathrm\{d\}\\boldsymbol\{x\}\)\\leq\\frac\{C^\{2\}\}\{p\}\.From Theorem[3\.4](https://arxiv.org/html/2608.14733#S3.Thmtheorem4), for any given functionu∈𝔽C,Ωu\\in\\mathbb\{F\}\_\{C,\\Omega\}, there is a sumup\(2\)=∑i=1pwi\(2\)​cos⁡\(𝑾i\(2\)​𝒙\+bi\(2\)\)u\_\{p\}^\{\(2\)\}=\\sum\_\{i=1\}^\{p\}w\_\{i\}^\{\(2\)\}\\cos\(\\boldsymbol\{W\}\_\{i\}^\{\(2\)\}\\boldsymbol\{x\}\+b\_\{i\}^\{\(2\)\}\), such that

∫Ω\|12​u​\(𝒙\)−up\(2\)​\(𝒙\)\|2​μ​\(𝑑𝒙\)≤C2p\.\\int\_\{\\Omega\}\\lvert\\frac\{1\}\{2\}u\(\\boldsymbol\{x\}\)\-u\_\{p\}^\{\(2\)\}\(\\boldsymbol\{x\}\)\\rvert^\{2\}\\mu\(\\,\\mathrm\{d\}\\boldsymbol\{x\}\)\\leq\\frac\{C^\{2\}\}\{p\}\.
SinceuM=up\(1\)\+up\(2\)u\_\{M\}=u\_\{p\}^\{\(1\)\}\+u\_\{p\}^\{\(2\)\}, then

∫Ω\|u⁡\(𝒙\)−uM​\(𝒙\)\|2​μ​\(d𝒙\)\\displaystyle\\int\_\{\\Omega\}\\lvert u\(\\boldsymbol\{x\}\)\-u\_\{M\}\(\\boldsymbol\{x\}\)\\rvert^\{2\}\\mu\(\\,\\mathrm\{d\}\\boldsymbol\{x\}\)=∫Ω\|12​u​\(𝒙\)−up\(1\)​\(𝒙\)\+12​u​\(𝒙\)−up\(2\)​\(𝒙\)\|2​μ​\(d𝒙\)\\displaystyle=\\int\_\{\\Omega\}\\lvert\\frac\{1\}\{2\}u\(\\boldsymbol\{x\}\)\-u\_\{p\}^\{\(1\)\}\(\\boldsymbol\{x\}\)\+\\frac\{1\}\{2\}u\(\\boldsymbol\{x\}\)\-u\_\{p\}^\{\(2\)\}\(\\boldsymbol\{x\}\)\\rvert^\{2\}\\mu\(\\,\\mathrm\{d\}\\boldsymbol\{x\}\)≤∫Ω2​\|12​u​\(𝒙\)−up\(1\)​\(𝒙\)\|2​μ​\(d𝒙\)\+∫Ω2​\|12​u​\(𝒙\)−up\(2\)​\(𝒙\)\|2​μ​\(d𝒙\)\\displaystyle\\leq\\int\_\{\\Omega\}2\\lvert\\frac\{1\}\{2\}u\(\\boldsymbol\{x\}\)\-u\_\{p\}^\{\(1\)\}\(\\boldsymbol\{x\}\)\\rvert^\{2\}\\mu\(\\,\\mathrm\{d\}\\boldsymbol\{x\}\)\+\\int\_\{\\Omega\}2\\lvert\\frac\{1\}\{2\}u\(\\boldsymbol\{x\}\)\-u\_\{p\}^\{\(2\)\}\(\\boldsymbol\{x\}\)\\rvert^\{2\}\\mu\(\\,\\mathrm\{d\}\\boldsymbol\{x\}\)≤2​C2p\+2​C2p=4​C2p\.\\displaystyle\\leq\\frac\{2C^\{2\}\}\{p\}\+\\frac\{2C^\{2\}\}\{p\}=\\frac\{4C^\{2\}\}\{p\}\.∎

## 4Numerical Experiments

In this section, we present numerical experiments to demonstrate the applicability and accuracy of the proposed FENs\. The experiments include two key tasks: function approximation and solving PDEs\. All experiments are conducted on a high\-performance computing server running Debian 12\. The server is equipped with an Intel Xeon Platinum 8358 CPU operating at 2\.60 GHz, and an NVIDIA A100 GPU with 80 GB of memory\. These robust hardware specifications provide the computational capacity required for the intensive operations involved in training and evaluating FENs\.

Moreover, to quantitatively evaluate the approximation capabilities of neural networks in the numerical experiments, the maximum absolute error \(L∞L\_\{\\infty\}error\) and the relativeL2L\_\{2\}error are defined as follows:

eL∞=max1≤i≤N⁡\|uM​\(𝒙i\)−ue​x​a​c​t​\(𝒙i\)\|,\\displaystyle e\_\{L\_\{\\infty\}\}=\\max\_\{1\\leq i\\leq N\}\|u\_\{M\}\(\\boldsymbol\{x\}\_\{i\}\)\-u\_\{exact\}\(\\boldsymbol\{x\}\_\{i\}\)\|,\(23\)eL2=∑i=1N\(uM​\(𝒙i\)−ue​x​a​c​t​\(𝒙i\)\)2∑i=1N\(ue​x​a​c​t​\(𝒙i\)\)2,\\displaystyle e\_\{L\_\{2\}\}=\\sqrt\{\\frac\{\\sum\_\{i=1\}^\{N\}\(u\_\{M\}\(\\boldsymbol\{x\}\_\{i\}\)\-u\_\{exact\}\(\\boldsymbol\{x\}\_\{i\}\)\)^\{2\}\}\{\\sum\_\{i=1\}^\{N\}\(u\_\{exact\}\(\\boldsymbol\{x\}\_\{i\}\)\)^\{2\}\}\},\(24\)whereuMu\_\{M\}anduexactu\_\{\\text\{exact\}\}represent the approximate and exact solutions, respectively, and𝒙i\\boldsymbol\{x\}\_\{i\}\(1≤i≤N1\\leq i\\leq N\) are the collocation points for error evaluation\.

### 4\.1Two\-dimensional function

To evaluate the capability of the FENs in function approximation, we implement them on a two\-dimensional function defined in Equation \([25](https://arxiv.org/html/2608.14733#S4.E25)\)

u⁡\(x,y\)=sin⁡\(π​x\)​sin⁡\(4​π​y\)\.u\(x,y\)=\\sin\(\\pi x\)\\sin\(4\\pi y\)\.\(25\)
In this experiment, we utilize a uniformly distributed training grid with dimensionsNx×Ny=101×101N\_\{x\}\\times N\_\{y\}=101\\times 101, whereNxN\_\{x\}andNyN\_\{y\}denote the number of points along thexx\-axis andyy\-axis, respectively\. The number of basis functionsMMis chosen to be400400,900900,16001600, and25002500, respectively\. For comparative analysis, we implement ELMs with three distinct activation functionssigmoid,tanh\\tanh, andswishto address the same problem\.

Table[1](https://arxiv.org/html/2608.14733#S4.T1)summarizes the parameters used for the approximation\. The optimal scaling factor for each activation function is determined by searching within a specified range with a small step size\. In Table[2](https://arxiv.org/html/2608.14733#S4.T2), we present the approximation errors for FENs and ELMs \(usingsigmoid,tanh\\tanh, andswishactivations\) across different basis function configurations \(M=400,900,1600,2500M=400,900,1600,2500\)\. As illustrated in Figure[2](https://arxiv.org/html/2608.14733#S4.F2), theL∞L\_\{\\infty\}error analysis reveals that FEN solutions exhibit significantly superior accuracy compared to ELM implementations\. These results demonstrate that the Fourier feature incorporation in FENs provides enhanced approximation capability for the target two\-dimensional function, as evidenced by the systematically lowerL∞L\_\{\\infty\}error metrics\.

Quantitatively, when the number of basis functions is sufficient, the smallestL∞L\_\{\\infty\}andL2L\_\{2\}errors achieved by FENs are6\.4756×10−156\.4756\\times 10^\{\-15\}and1\.4677×10−151\.4677\\times 10^\{\-15\}, respectively\. In contrast, the smallestL∞L\_\{\\infty\}andL2L\_\{2\}errors for ELMs are1\.7599×10−101\.7599\\times 10^\{\-10\}and1\.8013×10−111\.8013\\times 10^\{\-11\}, respectively\. These results highlight the significantly higher representational power of FENs compared to ELMs in this two\-dimensional approximation task\.

Table 1:Function approximation: Parameters used in approximating the two\-dimensional function defined by Equation \([25](https://arxiv.org/html/2608.14733#S4.E25)\)\.Activations\(ρm​i​n,ρm​a​x\]\(\\rho\_\{min\},\\rho\_\{max\}\]ρs\\rho\_\{s\}ρo​p​t\\rho\_\{opt\}M=400M=400M=900M=900M=1600M=1600M=2500M=2500sigmoid\(0,10\]\(0,10\]0\.011\.522\.192\.943\.81tanh\\tanh\(0,10\]\(0,10\]0\.010\.761\.471\.551\.92swish\(0,10\]\(0,10\]0\.011\.492\.613\.063\.76cos\\cos\(0,50\]\(0,50\]0\.17\.49\.59\.722\.8sin\\sin\(0,50\]\(0,50\]0\.17\.79\.110\.322\.7cos\\cos&\\&sin\\sin\(0,50\]\(0,50\]0\.19\.09\.69\.612\.6![Refer to caption](https://arxiv.org/html/2608.14733v1/figures/error_Func2D.png)

Figure 2:Function approximation:L∞L\_\{\\infty\}errors of neural networks in approximating the two\-dimensional function defined by Equation \([25](https://arxiv.org/html/2608.14733#S4.E25)\)\.Table 2:Function approximation: Performance comparison of FENs and ELMs withsigmoid,tanh\\tanh, andswishactivations in approximating the two\-dimensional function defined by Equation \([25](https://arxiv.org/html/2608.14733#S4.E25)\)\. TheL∞L\_\{\\infty\}andL2L\_\{2\}errors for each model configuration are presented\.ActivationsM=400M=900M=1600M=2500eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}sigmoid3\.2425E\-059\.2209E\-066\.3796E\-086\.2611E\-094\.0454E\-092\.7714E\-101\.7599E\-101\.8013E\-11tanh\\tanh2\.7418E\-056\.1370E\-061\.2666E\-071\.0780E\-083\.1869E\-092\.4560E\-101\.9782E\-101\.8775E\-11swish3\.3855E\-051\.0940E\-053\.5856E\-084\.5318E\-099\.6770E\-101\.6892E\-103\.7835E\-102\.0718E\-11cos\\cos1\.8362E\-077\.4465E\-081\.2727E\-142\.3325E\-159\.6481E\-151\.6508E\-153\.3154E\-143\.5549E\-15sin\\sin2\.2316E\-079\.5567E\-081\.0184E\-143\.9475E\-151\.5099E\-142\.0114E\-151\.7819E\-143\.4593E\-15cos\\cos&\\&sin\\sin2\.1903E\-077\.8576E\-081\.9601E\-144\.6683E\-156\.4756E\-151\.4677E\-152\.6069E\-142\.3770E\-15
### 4\.2Helmholtz equation

The Helmholtz equation in two dimensions is given by Equation \([26](https://arxiv.org/html/2608.14733#S4.E26)\)

\{Δ​u\+k2​u=q,in​Ω,u=h,on​∂Ω,\\begin\{array\}\[\]\{r@\{\}l\}\\left\\\{\\begin\{aligned\} \\Delta u\+k^\{2\}u&=q,&&\\mbox\{in\}\\kern 5\.0pt\\Omega,\\\\ u&=h,&&\\mbox\{on\}\\kern 5\.0pt\\partial\\Omega,\\end\{aligned\}\\right\.\\end\{array\}\(26\)whereΩ=\(0,1\)2\\Omega=\(0,1\)^\{2\}\. The exact solution is defined as

u⁡\(x,y\)=sin⁡\(a1​π​x\)​sin⁡\(a2​π​y\),u\(x,y\)=\\sin\(a\_\{1\}\\pi x\)\\sin\(a\_\{2\}\\pi y\),\(27\)with the corresponding source term given by

q⁡\(x,y\)=\(k2−\(a1​π\)2−\(a2​π\)2\)​sin⁡\(a1​π​x\)​sin⁡\(a2​π​y\),q\(x,y\)=\\left\(k^\{2\}\-\(a\_\{1\}\\pi\)^\{2\}\-\(a\_\{2\}\\pi\)^\{2\}\\right\)\\sin\(a\_\{1\}\\pi x\)\\sin\(a\_\{2\}\\pi y\),\(28\)where the parameters are set asa1=1a\_\{1\}=1,a2=4a\_\{2\}=4, andk=1k=1\.

For training both FENs and ELMs, we employ a uniform grid ofNx×Ny=101×101N\_\{x\}\\times N\_\{y\}=101\\times 101collocation points\. This training dataset comprises both interior points withinΩ\\Omegaand boundary points on∂Ω\\partial\\Omega\. Consequently, the datasetSSis divided into two distinct subsets:SrS\_\{r\}\(containing interior points\), andSbS\_\{b\}\(containing boundary points\)\.

Table[3](https://arxiv.org/html/2608.14733#S4.T3)summarizes the optimal scaling factors and search parameters used in solving the Helmholtz equation\. Table[4](https://arxiv.org/html/2608.14733#S4.T4)presents the approximation errors of FENs and ELMs with varying numbers of basis functions, all utilizing the optimized scaling configuration\. Figure[3](https://arxiv.org/html/2608.14733#S4.F3)illustrates theL∞L\_\{\\infty\}error trajectories for all models, demonstrating conclusively that FENs achieve superior approximation accuracy compared to ELMs across all parametric configurations\.

The data presented in Table[4](https://arxiv.org/html/2608.14733#S4.T4)reveal that with a sufficient number of basis functions, FENs attain minimumL∞L\_\{\\infty\}andL2L\_\{2\}errors of5\.3300×10−145\.3300\\times 10^\{\-14\}and2\.7200×10−142\.7200\\times 10^\{\-14\}, respectively\. In contrast, the best performance from ELMs yieldsL∞L\_\{\\infty\}andL2L\_\{2\}errors of4\.9477×10−104\.9477\\times 10^\{\-10\}and2\.7595×10−102\.7595\\times 10^\{\-10\}, respectively\. This result clearly demonstrates that FENs possess superior representational capacity and can achieve significantly higher accuracy than ELMs for the Helmholtz boundary value problem\.

Table 3:Helmholtz equation: Parameters when solving the Helmholtz equation \([26](https://arxiv.org/html/2608.14733#S4.E26)\) with solution in Equation \([27](https://arxiv.org/html/2608.14733#S4.E27)\)\.Activations\(ρm​i​n,ρm​a​x\]\(\\rho\_\{min\},\\rho\_\{max\}\]ρs\\rho\_\{s\}ρo​p​t\\rho\_\{opt\}M=400M=400M=900M=900M=1600M=1600M=2500M=2500sigmoid\(0,10\]\(0,10\]0\.011\.381\.662\.482\.98tanh\\tanh\(0,10\]\(0,10\]0\.010\.60\.941\.11\.5swish\(0,50\]\(0,50\]0\.011\.42\.02\.22\.4cos\\cos\(0,100\]\(0,100\]0\.15\.69\.39\.013\.1sin\\sin\(0,100\]\(0,100\]0\.15\.89\.99\.912\.0cos\\cos&\\&sin\\sin\(0,100\]\(0,100\]0\.15\.510\.510\.18\.8![Refer to caption](https://arxiv.org/html/2608.14733v1/figures/error_Helmholtz2D.png)

Figure 3:Helmholtz equation:L∞L\_\{\\infty\}errors of neural networks when solving Helmholtz equation \([26](https://arxiv.org/html/2608.14733#S4.E26)\) with solution in Equation \([27](https://arxiv.org/html/2608.14733#S4.E27)\)\.Table 4:Helmholtz equation: Performance comparison of FENs and ELMs activated bysigmoid,tanh\\tanhandswishwhen solving the Helmholtz equation \([26](https://arxiv.org/html/2608.14733#S4.E26)\) with solution in Equation \([27](https://arxiv.org/html/2608.14733#S4.E27)\)\. TheL∞L\_\{\\infty\}errors andL2L\_\{2\}errors for each model configuration are presented\.ActivationsM=400M=900M=1600M=2500eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}sigmoid3\.3875E\-038\.3152E\-046\.5565E\-075\.7295E\-078\.3819E\-092\.8206E\-096\.0754E\-103\.9546E\-10tanh\\tanh2\.7466E\-036\.7760E\-044\.4703E\-072\.0869E\-071\.1059E\-086\.0855E\-094\.9477E\-102\.7595E\-10swish6\.3477E\-031\.5769E\-037\.8976E\-072\.2362E\-076\.6357E\-093\.0735E\-091\.0987E\-093\.9585E\-10cos\\cos6\.3797E\-051\.7419E\-055\.8978E\-131\.5736E\-131\.7333E\-137\.0401E\-141\.7977E\-137\.1273E\-14sin\\sin1\.2699E\-041\.8981E\-059\.0587E\-132\.7940E\-131\.8736E\-137\.4318E\-145\.3300E\-142\.7200E\-14cos\\cos&\\&sin\\sin1\.2552E\-042\.2074E\-051\.5099E\-122\.4098E\-131\.4385E\-136\.6964E\-143\.7487E\-131\.3203E\-13\\chadded

For comparison, we also employ the finite element method \(FEM\) to solve the Helmholtz equation using a triangulation‌ with different mesh sizeshh\. The corresponding numerical results are presented in Table[5](https://arxiv.org/html/2608.14733#S4.T5)\. It can be observed that even with a large number of degree of freedom \(DoF\), the accuracy achieved by FEM is still much lower than that of FENs, which is shown in Table[4](https://arxiv.org/html/2608.14733#S4.T4)\.

\\chadded

Table 5:Helmholtz equation: Results of FEM when solving the Helmholtz equation \([26](https://arxiv.org/html/2608.14733#S4.E26)\) with solution in Equation \([27](https://arxiv.org/html/2608.14733#S4.E27)\)\. TheL∞L\_\{\\infty\}errors andL2L\_\{2\}errors for each model configuration are presented\.MethodhhDoFeL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}FEM0\.119691\.5408E\-027\.4130E\-030\.0575771\.5028E\-038\.9865E\-040\.02469051\.4795E\-044\.1033E\-050\.011859131\.8508E\-055\.0898E\-060\.0057412091\.3809E\-069\.4240E\-07\\chadded

To further assess the approximation capabilities of FENs and to compare with cos/sin\-based ELMs, we consider the following solution of Helmholtz equation \([26](https://arxiv.org/html/2608.14733#S4.E26)\),

u⁡\(x,y\)=tanh⁡\(x​y\),u\(x,y\)=\\tanh\(xy\),\(29\)with the associated source term defined by

q⁡\(x,y\)=\(k2−2​x2−2​y2\)​tanh⁡\(x​y\)\+\(2​x2\+2​y2\)​tanh3⁡\(x​y\)\.q\(x,y\)=\(k^\{2\}\-2x^\{2\}\-2y^\{2\}\)\\tanh\(xy\)\+\(2x^\{2\}\+2y^\{2\}\)\\tanh^\{3\}\(xy\)\.\(30\)\\chadded

Unlike the solution in Equation \([27](https://arxiv.org/html/2608.14733#S4.E27)\), this exact function is closely related to thetanh\\tanhactivation function\. Despite this, as shown in Table[7](https://arxiv.org/html/2608.14733#S4.T7), FENs still outperform ELMs activated by non\-trigonometric function, particularly when the number of basis functions is small\. The minimumL∞L\_\{\\infty\}andL2L\_\{2\}errors achieved by FENs are3\.2613×10−153\.2613\\times 10^\{\-15\}and2\.4689×10−152\.4689\\times 10^\{\-15\}, respectively, compared to4\.1078×10−154\.1078\\times 10^\{\-15\}and3\.8454×10−153\.8454\\times 10^\{\-15\}achieved by ELMs with trigonometric activation function\. They both achieve precise results, which demonstrate that activation function is very important during the computation\.

Table 6:Helmholtz equation: Parameters when solving the Helmholtz equation \([26](https://arxiv.org/html/2608.14733#S4.E26)\) with solution in Equation \([29](https://arxiv.org/html/2608.14733#S4.E29)\)\.Activations\(ρm​i​n,ρm​a​x\]\(\\rho\_\{min\},\\rho\_\{max\}\]ρs\\rho\_\{s\}ρo​p​t\\rho\_\{opt\}M=400M=400M=900M=900M=1600M=1600M=2500M=2500sigmoid\(0,10\]\(0,10\]0\.011\.512\.272\.473\.31tanh\\tanh\(0,10\]\(0,10\]0\.010\.731\.061\.531\.59swish\(0,50\]\(0,50\]0\.011\.62\.12\.84\.0cos\\cos\(0,100\]\(0,100\]0\.17\.69\.410\.212\.1sin\\sin\(0,100\]\(0,100\]0\.18\.39\.210\.315\.9cos\\cos&\\&sin\\sin\(0,100\]\(0,100\]0\.17\.59\.711\.223\.6Table 7:Helmholtz equation: Performance comparison of FENs and ELMs activated by different functions, when solving the Helmholtz equation \([26](https://arxiv.org/html/2608.14733#S4.E26)\) with solution \([29](https://arxiv.org/html/2608.14733#S4.E29)\)\. TheL∞L\_\{\\infty\}errors andL2L\_\{2\}errors for each model configuration are presented\.MethodsM=400M=900M=1600M=2500eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}ELMs \(OPENsigmoid\)\\text\{sigmoid\}\)2\.9184E\-081\.5916E\-082\.5498E\-111\.5287E\-112\.4958E\-131\.5588E\-132\.0040E\-142\.1436E\-14ELMs \(OPENtanh\)\\text\{tanh\}\)2\.3647E\-081\.2001E\-081\.6025E\-111\.0172E\-111\.8674E\-131\.5324E\-131\.1768E\-149\.1806E\-15ELMs \(OPENswish\)\\text\{swish\}\)2\.5343E\-081\.0115E\-081\.7768E\-111\.8566E\-113\.4794E\-133\.9641E\-133\.8497E\-145\.8373E\-14ELMs \(sin\\sin\)1\.0436E\-094\.3592E\-101\.6043E\-148\.7617E\-154\.1078E\-153\.8454E\-154\.2188E\-153\.9039E\-15FENs \(cos\\cos&\\&sin\\sin\)4\.0931E\-101\.6929E\-101\.2101E\-149\.6805E\-153\.2613E\-152\.4689E\-153\.7192E\-152\.2473E\-15
### 4\.3Diffusion equation

In this subsection, we investigate the behavior of a two\-dimensional diffusion equation that involves both spatial and temporal variables\.The initial\-boundary value problem is governed by

\{∂u∂t−ν​∂2u∂x2=f⁡\(x,t\),\(x,t\)∈\(a1,b1\)×\(0,tf\],u⁡\(a1,t\)=g1​\(t\),t∈\(0,tf\],u⁡\(b1,t\)=g2​\(t\),t∈\(0,tf\],u⁡\(x,0\)=h⁡\(x\),x∈\[a1,b1\],\\begin\{array\}\[\]\{r@\{\}l\}\\left\\\{\\begin\{aligned\} \\frac\{\\partial u\}\{\\partial t\}\-\\nu\\frac\{\\partial^\{2\}u\}\{\\partial x^\{2\}\}&=f\(x,t\),&&\(x,t\)\\in\\ \(a\_\{1\},b\_\{1\}\)\\times\(0,t\_\{f\}\],\\\\ u\(a\_\{1\},t\)&=g\_\{1\}\(t\),&&t\\in\(0,t\_\{f\}\],\\\\ u\(b\_\{1\},t\)&=g\_\{2\}\(t\),&&t\\in\(0,t\_\{f\}\],\\\\ u\(x,0\)&=h\(x\),&&x\\in\[a\_\{1\},b\_\{1\}\],\\\\ \\end\{aligned\}\\right\.\\end\{array\}\(31\)wheref⁡\(x,t\)f\(x,t\)denotes the source term,ν\>0\\nu\>0denotes the constant diffusion coefficient,g1​\(t\)g\_\{1\}\(t\)andg2​\(t\)g\_\{2\}\(t\)prescribe the time\-dependent Dirichlet boundary conditions, andh⁡\(x\)h\(x\)defines the initial condition\. The parameters are chosen asa1=0a\_\{1\}=0,b1=5b\_\{1\}=5andν=0\.01\\nu=0\.01\. The final timetft\_\{f\}considered in simulations is set to be11\.

We choose the suitable functionsf⁡\(x,t\)f\(x,t\),g1​\(t\)g\_\{1\}\(t\),g2​\(t\)g\_\{2\}\(t\), and the initial conditionh⁡\(x\)h\(x\)so that the exact solution is given by

u⁡\(x,t\)=\[2​cos⁡\(π​x\+π5\)\+32​cos⁡\(2​π​x−3​π5\)\]​\[2​cos⁡\(π​t\+π5\)\+32​cos⁡\(2​π​t−3​π5\)\]\.u\(x,t\)=\\left\[2\\cos\\left\(\\pi x\+\\frac\{\\pi\}\{5\}\\right\)\+\\frac\{3\}\{2\}\\cos\\left\(2\\pi x\-\\frac\{3\\pi\}\{5\}\\right\)\\right\]\\left\[2\\cos\\left\(\\pi t\+\\frac\{\\pi\}\{5\}\\right\)\+\\frac\{3\}\{2\}\\cos\\left\(2\\pi t\-\\frac\{3\\pi\}\{5\}\\right\)\\right\]\.\(32\)By treating time as a spatial dimension, we convert the 1D time\-dependent diffusion equation into a 2D problem\. Both FENs and ELMs are trained on aNx×Ny=101×101N\_\{x\}\\times N\_\{y\}=101\\times 101collocation points\. Table[8](https://arxiv.org/html/2608.14733#S4.T8)lists the relevant parameters and the optimal scaling factor used in this setup\. Table[9](https://arxiv.org/html/2608.14733#S4.T9)compares the approximate errors of the solutions obtained by FENs and ELMs for different numbers of basis functions\. Figure[4](https://arxiv.org/html/2608.14733#S4.F4)illustrates theL∞L\_\{\\infty\}error curves, showing significantly lower errors for FENs than ELMs\.

Furthermore, by examining Table[9](https://arxiv.org/html/2608.14733#S4.T9), we observe that the minimalL∞L\_\{\\infty\}andL2L\_\{2\}errors achieved by FENs are9\.2371×10−149\.2371\\times 10^\{\-14\}and5\.9186×10−155\.9186\\times 10^\{\-15\}, respectively\. In contrast, the minimalL∞L\_\{\\infty\}andL2L\_\{2\}errors obtained by ELMs are5\.7246×10−95\.7246\\times 10^\{\-9\}and3\.8524×10−103\.8524\\times 10^\{\-10\}, respectively\. These results demonstrate that FENs are more suitable for solving this problem compared to ELMs\.

Table 8:Diffusion equation: Parameters when solving the diffusion equation \([31](https://arxiv.org/html/2608.14733#S4.E31)\) withtf=1t\_\{f\}=1\.Activations\(ρm​i​n,ρm​a​x\]\(\\rho\_\{min\},\\rho\_\{max\}\]ρs\\rho\_\{s\}ρo​p​t\\rho\_\{opt\}M=400M=400M=900M=900M=1600M=1600M=2500M=2500sigmoid\(0,10\]\(0,10\]0\.011\.551\.862\.613\.35tanh\\tanh\(0,10\]\(0,10\]0\.010\.91\.011\.382\.03swish\(0,50\]\(0,50\]0\.011\.552\.52\.533\.96cos\\cos\(0,100\]\(0,100\]0\.15\.55\.05\.46\.4sin\\sin\(0,100\]\(0,100\]0\.14\.05\.55\.36\.6cos\\cos&\\&sin\\sin\(0,100\]\(0,100\]0\.15\.65\.25\.45\.3![Refer to caption](https://arxiv.org/html/2608.14733v1/figures/error_Diffusion.png)

Figure 4:Diffusion equation: theL∞L\_\{\\infty\}errors of neural networks when solving diffusion equation \([31](https://arxiv.org/html/2608.14733#S4.E31)\) withtf=1t\_\{f\}=1\.Table 9:Diffusion equation: Performance comparison of FENs and ELMs activated bysigmoid,tanh\\tanhandswishwhen solving the diffusion equation \([31](https://arxiv.org/html/2608.14733#S4.E31)\) withtf=1t\_\{f\}=1\. TheL∞L\_\{\\infty\}errors andL2L\_\{2\}errors for each model configuration are presented\.ActivationsM=400M=900M=1600M=2500eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}sigmoid1\.0276E\-017\.3530E\-032\.6306E\-052\.4873E\-062\.4815E\-071\.8386E\-085\.7246E\-093\.8524E\-10tanh\\tanh1\.4518E\-017\.2468E\-031\.9460E\-051\.1581E\-062\.4414E\-071\.7353E\-081\.0454E\-084\.3533E\-10swish6\.3047E\-024\.5548E\-031\.8239E\-051\.5684E\-061\.1645E\-076\.2690E\-091\.1023E\-089\.1794E\-10cos\\cos7\.4218E\-095\.5655E\-101\.3012E\-135\.9433E\-151\.1724E\-139\.1437E\-151\.8852E\-139\.1456E\-15sin\\sin1\.4405E\-088\.5286E\-109\.2371E\-145\.9186E\-151\.0303E\-138\.3162E\-152\.1139E\-136\.9327E\-15cos\\cos&\\&sin\\sin9\.6760E\-091\.3414E\-102\.1139E\-131\.1447E\-141\.1546E\-131\.3558E\-149\.7700E\-148\.3418E\-15
### 4\.4Heat equation

We consider the heat equation within the spatial\-temporal domainΩ×\(0,tf\]=\(0,1\)2×\(0,1\]\\Omega\\times\(0,t\_\{f\}\]=\(0,1\)^\{2\}\\times\(0,1\], governed by the following system of PDEs

\{ut​\(x,y,t\)−Δ​u​\(x,y,t\)=f⁡\(x,y,t\),\(x,y,t\)∈Ω×\(0,1\],u⁡\(x,y,t\)=g⁡\(x,y,t\),\(x,y,t\)∈∂Ω×\(0,1\],u⁡\(x,y,0\)=h⁡\(x,y\),\(x,y\)∈Ω\.\\begin\{array\}\[\]\{r@\{\}l\}\\left\\\{\\begin\{aligned\} u\_\{t\}\(x,y,t\)\-\\Delta u\(x,y,t\)&=f\(x,y,t\),&&\(x,y,t\)\\in\\Omega\\times\(0,1\],\\\\ u\(x,y,t\)&=g\(x,y,t\),&&\(x,y,t\)\\in\\partial\\Omega\\times\(0,1\],\\\\ u\(x,y,0\)&=h\(x,y\),&&\(x,y\)\\in\\Omega\.\\\\ \\end\{aligned\}\\right\.\\end\{array\}\(33\)We choose the suitable functionsf⁡\(x,y,t\)f\(x,y,t\),g⁡\(x,y,t\)g\(x,y,t\), and the initial conditionh⁡\(x,y,t\)h\(x,y,t\)so that the exact solution is given by

u⁡\(x,y,t\)=2​e−t​sin⁡\(π2​x\)​sin⁡\(π2​y\)\.u\(x,y,t\)=2e^\{\-t\}\\sin\\left\(\\frac\{\\pi\}\{2\}x\\right\)\\sin\\left\(\\frac\{\\pi\}\{2\}y\\right\)\.
For the training phase, a uniform grid ofNx×Ny×Nt=51×51×51N\_\{x\}\\times N\_\{y\}\\times N\_\{t\}=51\\times 51\\times 51collocation points is employed\. The parameters and results from the optimal scale search are summarized in Table[10](https://arxiv.org/html/2608.14733#S4.T10)\. Table[11](https://arxiv.org/html/2608.14733#S4.T11)presents theL∞L\_\{\\infty\}andL2L\_\{2\}errors for both FENs and ELMs across varying numbers of basis functions\. TheL∞L\_\{\\infty\}error curves depicted in Figure[5](https://arxiv.org/html/2608.14733#S4.F5)clearly demonstrate that FENs consistently yield lower approximation errors compared to ELMs, highlighting their superior performance in solving the heat equation\.

In terms of specific approximation error values, when the number of basis functions is sufficiently large, FENs can achieve the smallestL∞L\_\{\\infty\}andL2L\_\{2\}errors of2\.6090×10−142\.6090\\times 10^\{\-14\}and1\.7562×10−141\.7562\\times 10^\{\-14\}, respectively\. In contrast, ELMs reach their smallestL∞L\_\{\\infty\}andL2L\_\{2\}errors of2\.4443×10−122\.4443\\times 10^\{\-12\}and1\.4761×10−121\.4761\\times 10^\{\-12\}, respectively\. This significant difference in error magnitudes clearly demonstrates that FENs can provide solutions with considerably higher precision than ELMs, making them a more suitable method for solving this problem\.

Table 10:Heat equation: Parameters when solving the heat equation \([33](https://arxiv.org/html/2608.14733#S4.E33)\)\.Activations\(ρm​i​n,ρm​a​x\]\(\\rho\_\{min\},\\rho\_\{max\}\]ρs\\rho\_\{s\}ρo​p​t\\rho\_\{opt\}M=400M=400M=900M=900M=1600M=1600M=2500M=2500sigmoid\(0,5\]\(0,5\]0\.010\.150\.30\.460\.56tanh\\tanh\(0,5\]\(0,5\]0\.010\.070\.140\.190\.22swish\(0,5\]\(0,5\]0\.010\.120\.210\.310\.48cos\\cos\(0,50\]\(0,50\]0\.11\.11\.42\.21\.7sin\\sin\(0,50\]\(0,50\]0\.11\.21\.31\.62\.4cos\\cos&\\&sin\\sin\(0,50\]\(0,50\]0\.11\.11\.31\.62\.0![Refer to caption](https://arxiv.org/html/2608.14733v1/figures/error_Heat.png)

Figure 5:Heat equation:L∞L\_\{\\infty\}errors of neural networks when solving the heat equation \([33](https://arxiv.org/html/2608.14733#S4.E33)\)\.Table 11:Heat equation: Performance comparison of FENs and ELMs activated bysigmoid,tanh\\tanhandswishwhen solving the heat equation \([33](https://arxiv.org/html/2608.14733#S4.E33)\)\. TheL∞L\_\{\\infty\}errors andL2L\_\{2\}errors for each model configuration are presented\.ActivationsM=400M=900M=1600M=2500eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}sigmoid1\.5832E\-077\.3817E\-082\.4738E\-102\.3643E\-103\.0013E\-112\.2705E\-115\.9936E\-123\.5015E\-12tanh\\tanh1\.2910E\-075\.3362E\-081\.1386E\-105\.7319E\-111\.5007E\-117\.7614E\-124\.1442E\-123\.0037E\-12swish1\.4388E\-077\.5731E\-081\.4778E\-107\.2960E\-112\.0669E\-118\.1287E\-122\.4443E\-121\.4761E\-12cos\\cos1\.5015E\-074\.9863E\-082\.4385E\-121\.2761E\-122\.9932E\-131\.2883E\-134\.5741E\-142\.8907E\-14sin\\sin6\.6993E\-083\.8847E\-082\.8706E\-121\.1786E\-121\.2257E\-139\.7663E\-142\.6090E\-141\.7562E\-14cos\\cos&\\&sin\\sin6\.4772E\-082\.3952E\-082\.1253E\-129\.1364E\-132\.1938E\-131\.3437E\-133\.1822E\-143\.0063E\-14
### 4\.5Wave equation

We consider the following wave equation defined over the spatial domainΩ=\(0,1\)2\\Omega=\(0,1\)^\{2\}and the temporal domain\(0,1\]\(0,1\],

\{∂2u∂t2−Δ​u​\(x,y,t\)=f⁡\(x,y,t\),\(x,y,t\)∈Ω×\(0,1\],u⁡\(x,y,t\)=g⁡\(x,y,t\),\(x,y,t\)∈∂Ω×\(0,1\],u⁡\(x,y,0\)=h⁡\(x,y\),\(x,y\)∈Ω,∂u∂t​\(x,y,0\)=w⁡\(x,y\),\(x,y\)∈Ω\.\\begin\{array\}\[\]\{r@\{\}l\}\\left\\\{\\begin\{aligned\} \\frac\{\\partial^\{2\}u\}\{\\partial t^\{2\}\}\-\\Delta u\(x,y,t\)&=f\(x,y,t\),&&\(x,y,t\)\\in\\Omega\\times\(0,1\],\\\\ u\(x,y,t\)&=g\(x,y,t\),&&\(x,y,t\)\\in\\partial\\Omega\\times\(0,1\],\\\\ u\(x,y,0\)&=h\(x,y\),&&\(x,y\)\\in\\Omega,\\\\ \\frac\{\\partial u\}\{\\partial t\}\(x,y,0\)&=w\(x,y\),&&\(x,y\)\\in\\Omega\.\\\\ \\end\{aligned\}\\right\.\\end\{array\}\(34\)The exact solution is given byu⁡\(x,y,t\)=sin⁡\(π2​x\)​sin⁡\(π2​y\)​sin⁡\(π2​t\)u\(x,y,t\)=\\sin\\left\(\\frac\{\\pi\}\{2\}x\\right\)\\sin\\left\(\\frac\{\\pi\}\{2\}y\\right\)\\sin\\left\(\\frac\{\\pi\}\{2\}t\\right\), with suitable boundary conditiong⁡\(x,y,t\)g\(x,y,t\), initial conditionsh⁡\(x,y,t\)h\(x,y,t\)andw⁡\(x,y\)w\(x,y\), as well as the source termf⁡\(x,y,t\)f\(x,y,t\)\.We employ a uniform grid ofNx×Ny×Nt=51×51×51N\_\{x\}\\times N\_\{y\}\\times N\_\{t\}=51\\times 51\\times 51collocation points\.

Following the same procedure as in previous sections, we conduct an optimal scale search across varying numbers of basis functions\. The corresponding parameters are listed in Table[12](https://arxiv.org/html/2608.14733#S4.T12)\. The resulting errors obtained using FENs and ELMs with their respective optimal scaling factors are presented in Table[13](https://arxiv.org/html/2608.14733#S4.T13)\. Figure[6](https://arxiv.org/html/2608.14733#S4.F6)illustrates theL∞L\_\{\\infty\}andL2L\_\{2\}error curves for both computational approaches, clearly showing that FENs maintain superior error\-reduction capabilities throughout all test configurations when benchmarked against ELMs\.

Furthermore, as quantitatively demonstrated in Table[13](https://arxiv.org/html/2608.14733#S4.T13), the minimumL∞L\_\{\\infty\}andL2L\_\{2\}errors achieved by FENs are2\.4425×10−152\.4425\\times 10^\{\-15\}and1\.7659×10−151\.7659\\times 10^\{\-15\}, respectively\. In contrast, minimumL∞L\_\{\\infty\}andL2L\_\{2\}errors achieved by ELMs are1\.4412×10−121\.4412\\times 10^\{\-12\}and5\.1755×10−135\.1755\\times 10^\{\-13\}, respectively\. This systematic comparison conclusively establishes FENs’ superior efficacy in wave equation solutions, with demonstrably higher computational precision compared to ELMs\.

Table 12:Wave equation: Parameters when solving the wave equation \([34](https://arxiv.org/html/2608.14733#S4.E34)\)\.Activations\(ρm​i​n,ρm​a​x\]\(\\rho\_\{min\},\\rho\_\{max\}\]ρs\\rho\_\{s\}ρo​p​t\\rho\_\{opt\}M=400M=400M=900M=900M=1600M=1600M=2500M=2500sigmoid\(0,5\]\(0,5\]0\.010\.170\.270\.50\.54tanh\\tanh\(0,5\]\(0,5\]0\.010\.060\.130\.180\.24swish\(0,5\]\(0,5\]0\.010\.130\.240\.350\.49cos\\cos\(0,10\]\(0,10\]0\.11\.281\.341\.461\.7sin\\sin\(0,10\]\(0,10\]0\.11\.31\.381\.471\.39cos\\cos&\\&sin\\sin\(0,10\]\(0,10\]0\.11\.161\.281\.621\.74![Refer to caption](https://arxiv.org/html/2608.14733v1/figures/error_Wave.png)

Figure 6:Wave equation:L∞L\_\{\\infty\}errors of neural networks when solving the wave equation \([34](https://arxiv.org/html/2608.14733#S4.E34)\)\.Table 13:Wave equation: Performance comparison of FENs and ELMs activated bysigmoid,tanh\\tanhandswishwhen solving the wave equation \([34](https://arxiv.org/html/2608.14733#S4.E34)\)\. TheL∞L\_\{\\infty\}errors andL2L\_\{2\}errors for each model configuration are presented\.ActivationsM=400M=900M=1600M=2500eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}sigmoid2\.4436E\-076\.4600E\-081\.9004E\-101\.3386E\-102\.0691E\-117\.1025E\-122\.1174E\-121\.0359E\-12tanh\\tanh1\.7975E\-073\.5759E\-088\.8139E\-114\.3830E\-118\.2991E\-123\.5998E\-121\.6200E\-125\.1755E\-13swish1\.0241E\-073\.8084E\-088\.5493E\-113\.3717E\-119\.0949E\-124\.1869E\-121\.4412E\-126\.2451E\-13cos\\cos3\.7392E\-081\.6164E\-082\.1672E\-131\.2717E\-135\.2736E\-152\.3601E\-152\.4425E\-151\.7659E\-15sin\\sin6\.0666E\-082\.4112E\-084\.3354E\-132\.1115E\-136\.9944E\-153\.3327E\-155\.6621E\-153\.9232E\-15cos\\cos&\\&sin\\sin3\.5241E\-081\.0682E\-083\.5172E\-131\.5920E\-135\.3291E\-152\.0016E\-152\.8588E\-152\.6000E\-15
### 4\.6Nonlinear Helmholtz equation

For the nonlinear example, we evaluate the performance of FENs and ELMs on a boundary value problem governed by the one\-dimensional nonlinear Helmholtz equation\. The formulation is defined as follows

\{∂2u∂x2−λ​u\+β​sin⁡\(u\)=f⁡\(x\),x∈\(a,b\),u⁡\(a\)=h1,u⁡\(b\)=h2\.\\begin\{array\}\[\]\{r@\{\}l\}\\left\\\{\\begin\{aligned\} \\frac\{\\partial^\{2\}u\}\{\\partial x^\{2\}\}\-\\lambda u\+\\beta\\sin\(u\)&=f\(x\),&&x\\in\(a,b\),\\\\ u\(a\)&=h\_\{1\},\\\\ u\(b\)&=h\_\{2\}\.\\\\ \\end\{aligned\}\\right\.\\end\{array\}\(35\)The exact solution is given byu⁡\(x\)=sin⁡\(3​π​x\+3​π20\)​cos⁡\(4​π​x−2​π5\)\+32\+x10u\(x\)=\\sin\\left\(3\\pi x\+\\frac\{3\\pi\}\{20\}\\right\)\\cos\\left\(4\\pi x\-\\frac\{2\\pi\}\{5\}\\right\)\+\\frac\{3\}\{2\}\+\\frac\{x\}\{10\}, with suitable boundary conditions and source termf⁡\(x\)f\(x\)\. The constant parameters in the equation are set asa=0a=0,b=8b=8,λ=50\\lambda=50andβ=10\\beta=10\.

In this nonlinear case, we useNx=3000N\_\{x\}=3000uniform collocation points for training\. To handle the nonlinearity of the problem, we adopt the Picard iteration method, performing a total of100100iterations\. Initially, a vector of coefficients𝒘0\\boldsymbol\{w\}\_\{0\}is randomly initialized, generating‌ the approximate solutionu0=Φ⋅𝒘0u\_\{0\}=\\Phi\\cdot\\boldsymbol\{w\}\_\{0\}\. The nonlinear termβ​sin⁡\(u0\)\\beta\\sin\(u\_\{0\}\)is then computed and moved to the right\-hand side of the equation, thereby linearizing the problem ‌for the current iteration step‌\. Solving the resulting linear system yields an updated coefficient vector𝒘1\\boldsymbol\{w\}\_\{1\}\. This iterative process continues, updating the coefficients at each step until convergence is achieved\. The final approximate solution is expressed asu=Φ⋅𝒘u=\\Phi\\cdot\\boldsymbol\{w\}, where𝒘\\boldsymbol\{w\}is the ‌converged‌ coefficient vector\.

For the nonlinear Helmholtz problem, the parameters ‌along with‌ the corresponding optimal scaling factor ‌identified via scale search‌ are ‌listed‌ in Table[14](https://arxiv.org/html/2608.14733#S4.T14)\. The resultingL∞L\_\{\\infty\}andL2L\_\{2\}errors from FENs and ELMs using the optimal scaling factor‌ are ‌presented‌ in Table[15](https://arxiv.org/html/2608.14733#S4.T15)\. Figure[7](https://arxiv.org/html/2608.14733#S4.F7)displays theL∞L\_\{\\infty\}error curves for both methods\. As clearly shown in the figure‌, the error curves ‌for‌ FENs are consistently and significantly lower than ‌those for‌ ELMs, ‌demonstrating‌ the superior accuracy and robustness of FENs in solving this nonlinear problem\.

For the nonlinear Helmholtz problem ‌studied‌, the ELM with thetanh\\tanhactivation function ‌achieves‌ clearly ‌superior‌ performance compared to ‌networks using‌sigmoidandswish, ‌delivering‌ notably higher accuracy ‌while still lagging behind FENs‌\. This ‌demonstrates‌ that the choice of activation function ‌plays a critical role‌ in the efficacy of ELMs, withtanh\\tanh‌providing‌ a distinct advantage ‌in this specific context‌\. ‌However‌, FENs ‌maintain‌ superior representational ‌capacity‌ and precision, ‌highlighting‌ their robustness in addressing ‌this class of‌ nonlinear problems\.

Upon examining the approximation error values in Table[15](https://arxiv.org/html/2608.14733#S4.T15), we observe that the smallestL∞L\_\{\\infty\}andL2L\_\{2\}errors achieved by FENs are2\.2427×10−132\.2427\\times 10^\{\-13\}and1\.3705×10−151\.3705\\times 10^\{\-15\}, respectively\. In contrast, the lowestL∞L\_\{\\infty\}andL2L\_\{2\}errors obtained by ELMs are3\.7986×10−113\.7986\\times 10^\{\-11\}and4\.4912×10−124\.4912\\times 10^\{\-12\}, respectively\. These results clearly demonstrate that FENs are more suitable for solving this nonlinear problem and are capable of achieving significantly higher solution accuracy compared to ELMs\.

Table 14:Nonlinear Helmholtz equation: Parameters when solving the nonlinear Helmholtz equation \([35](https://arxiv.org/html/2608.14733#S4.E35)\)\.Activations\(ρm​i​n,ρm​a​x\]\(\\rho\_\{min\},\\rho\_\{max\}\]ρs\\rho\_\{s\}ρo​p​t\\rho\_\{opt\}M=400M=400M=900M=900M=1600M=1600M=2500M=2500sigmoid\(0,10\]\(0,10\]0\.019\.659\.999\.959\.72tanh\\tanh\(0,10\]\(0,10\]0\.019\.659\.719\.99\.98swish\(0,10\]\(0,10\]0\.0110\.09\.949\.949\.84cos\\cos\(0,100\]\(0,100\]116333464sin\\sin\(0,100\]\(0,100\]120342474cos\\cos&\\&sin\\sin\(0,100\]\(0,100\]122253371![Refer to caption](https://arxiv.org/html/2608.14733v1/figures/error_Helmholtz1D.png)

Figure 7:Nonlinear Helmholtz equation:L∞L\_\{\\infty\}errors of neural networks when solving the nonlinear Helmholtz equation \([35](https://arxiv.org/html/2608.14733#S4.E35)\)\.Table 15:Nonlinear Helmholtz equation: Performance comparison of FENs and ELMs activated bysigmoid,tanh\\tanhandswishwhen solving the nonlinear Helmholtz equation \([35](https://arxiv.org/html/2608.14733#S4.E35)\)\. TheL∞L\_\{\\infty\}errors andL2L\_\{2\}errors for each model configuration are presented\.ActivationsM=400M=900M=1600M=2500eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}sigmoid3\.4384E\-041\.6943E\-055\.2436E\-062\.4732E\-071\.5165E\-067\.1510E\-081\.4590E\-067\.7030E\-08tanh\\tanh2\.0442E\-061\.0492E\-071\.0126E\-084\.7025E\-103\.7986E\-114\.4912E\-125\.6331E\-103\.1179E\-11swish3\.1724E\-041\.9411E\-059\.2887E\-065\.5102E\-071\.4125E\-066\.6785E\-081\.4330E\-068\.1166E\-08cos\\cos2\.9874E\-121\.5546E\-132\.2427E\-131\.3705E\-149\.2371E\-134\.8289E\-147\.0832E\-133\.5773E\-14sin\\sin8\.6642E\-134\.2264E\-148\.4621E\-134\.6449E\-141\.4064E\-126\.4684E\-147\.4474E\-133\.7450E\-14cos\\cos&\\&sin\\sin6\.0973E\-132\.8716E\-145\.8753E\-133\.4123E\-141\.3309E\-126\.4362E\-146\.6214E\-133\.1905E\-14
### 4\.7Poisson equation with an oscillating solution

We consider the following one\-dimensional Poisson equation:

\{−Δ​u​\(x\)=f⁡\(x\),in​Ω,u⁡\(x\)=g⁡\(x\),on​∂Ω,\\begin\{array\}\[\]\{r@\{\}l\}\\left\\\{\\begin\{aligned\} \-\\Delta u\(x\)&=f\(x\),&&\\mbox\{in\}\\kern 5\.0pt\\Omega,\\\\ u\(x\)&=g\(x\),&&\\mbox\{on\}\\kern 5\.0pt\\partial\\Omega,\\end\{aligned\}\\right\.\\end\{array\}\(36\)whereΩ=\(0,1\)\\Omega=\(0,1\)\. We choose the suitablef⁡\(x\)f\(x\)and boundary conditiong⁡\(x\)g\(x\)so that the exact solution is given by

u⁡\(x\)=16​∑i=16sin⁡\(2i​π​x\),u\(x\)=\\frac\{1\}\{6\}\\sum\_\{i=1\}^\{6\}\\sin\(2^\{i\}\\pi x\),\(37\)which represents a superposition of sine functions with exponentially increasing frequencies, producing a highly oscillatory ‌behavior\.

We train the networks usingNx=3000N\_\{x\}=3000uniformly distributed collocation points\. The parameters associated with different neural architectures‌ during the optimal scaling factor search across different numbers of basis functions, are summarized in Table[17](https://arxiv.org/html/2608.14733#S4.T17)\. Notably, the admissible range for the activation functiontanh\\tanhis narrower than those ofsigmoidandswish\. This constraint ‌originates from‌ numerical stability requirements: ‌applying‌ excessively large scaling factors totanh\\tanh‌generates‌ ‌severely‌ ill\-conditioned matrices \(‌approaching singularity‌\), ‌which destabilizes the solving process\. ‌To mitigate this, thetanh\\tanhscaling range is ‌strategically restricted‌ to maintain ‌trainable system ‌conditions\.

Table[16](https://arxiv.org/html/2608.14733#S4.T16)compiles the optimal scaling factors ‌and their associated search parameters\. Furthermore, Figure[8](https://arxiv.org/html/2608.14733#S4.F8)shows a comparison between the exact solution and the numerical solution obtained using the FEN with thecos\\cosactivation function, where the number of basis functions is set toM=900M=900and the optimal scaling factor isρopt=130\\rho\_\{\\text\{opt\}\}=130\. It is observed that the two curves exhibit an almost perfect overlap, demonstrating that the neural network has successfully captured the exact solution with high accuracy\.

Figure[9](https://arxiv.org/html/2608.14733#S4.F9)quantifies theL∞L\_\{\\infty\}convergence of neural networks ‌equipped with‌ optimal scaling factors for solving the ‌highly oscillatory‌ Poisson equation \([37](https://arxiv.org/html/2608.14733#S4.E37)\)\. It is clearly observed that the errors of ELMs withsigmoid,swish, andtanh\\tanhactivation functions are significantly higher than those of FENs\. In Table[17](https://arxiv.org/html/2608.14733#S4.T17), we report theL∞L\_\{\\infty\}andL2L\_\{2\}errors corresponding to various numbers of basis functions\. The minimumL∞L\_\{\\infty\}andL2L\_\{2\}errors achieved by ELMs are6\.1572×10−86\.1572\\times 10^\{\-8\}and8\.0396×10−88\.0396\\times 10^\{\-8\}, respectively, whereas FENs achieve notably smaller errors of1\.3878×10−111\.3878\\times 10^\{\-11\}and3\.0395×10−113\.0395\\times 10^\{\-11\}\. These results highlight the superior representational power of FENs and demonstrate their greater suitability for solving this highly oscillatory Poisson problem\.

Table 16:Poisson equation with an oscillating solution: Parameters when solving the Poisson equation \([36](https://arxiv.org/html/2608.14733#S4.E36)\)\.Activations\(ρm​i​n,ρm​a​x\]\(\\rho\_\{min\},\\rho\_\{max\}\]ρs\\rho\_\{s\}ρo​p​t\\rho\_\{opt\}M=400M=400M=900M=900M=1600M=1600M=2500M=2500sigmoid\(0,20\]\(0,20\]0\.0118\.7918\.1218\.4619\.19tanh\\tanh\(0,10\]\(0,10\]0\.019\.159\.839\.939\.70swish\(0,20\]\(0,20\]0\.0119\.9019\.3319\.3918\.95cos\\cos\(0,1000\]\(0,1000\]1176130205564sin\\sin\(0,1000\]\(0,1000\]1147150234890cos\\cos&\\&sin\\sin\(0,1000\]\(0,1000\]1164146325903![Refer to caption](https://arxiv.org/html/2608.14733v1/figures/Poisson1D.png)

Figure 8:Poisson equation with an oscillating solution: comparison between the exact solution and the approximate solution obtained using FEN with acos\\cosactivation, whereM=900M=900andρo​p​t=130\\rho\_\{opt\}=130\.![Refer to caption](https://arxiv.org/html/2608.14733v1/figures/error_Poisson1D.png)

Figure 9:Poisson equation with an oscillating solution:L∞L\_\{\\infty\}errors of neural networks when solving the Poisson equation \([36](https://arxiv.org/html/2608.14733#S4.E36)\)\.Table 17:Poisson equation with an oscillating solution: Performance comparison of FENs and ELMs activated bysigmoid,tanh\\tanhandswishwhen solving the Poisson equation \([36](https://arxiv.org/html/2608.14733#S4.E36)\)\. TheL∞L\_\{\\infty\}errors andL2L\_\{2\}errors for each model configuration are presented\.ActivationsM=400M=900M=1600M=2500eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}sigmoid3\.6288E\-053\.4848E\-052\.9000E\-076\.0888E\-071\.6587E\-071\.8036E\-071\.0117E\-061\.0275E\-06tanh\\tanh3\.0062E\-053\.8163E\-072\.6801E\-072\.8331E\-076\.1572E\-088\.0396E\-082\.0924E\-071\.7937E\-07swish1\.7748E\-042\.4055E\-041\.4265E\-061\.5708E\-065\.9915E\-078\.1065E\-071\.9759E\-062\.1461E\-06cos\\cos3\.0724E\-119\.7238E\-111\.3878E\-113\.0395E\-117\.2102E\-119\.4554E\-111\.7749E\-113\.2850E\-11sin\\sin3\.2458E\-115\.7302E\-112\.9600E\-114\.9456E\-112\.7607E\-104\.0441E\-104\.6233E\-119\.1841E\-11cos\\cos&\\&sin\\sin1\.9459E\-113\.0421E\-116\.1339E\-111\.1078E\-101\.7186E\-112\.8037E\-111\.9275E\-113\.1528E\-11
### 4\.8Nonlinear Burgers’ equation

Consider the following 1D Burgers’ equation

\{ut\+u​ux−ϵ​ux​x=f,in​Ω×\(0,tf\],u=g,on​∂Ω×\(0,tf\],u=u0,in​Ω,\\begin\{array\}\[\]\{r@\{\}l\}\\left\\\{\\begin\{aligned\} u\_\{t\}\+uu\_\{x\}\-\\epsilon u\_\{xx\}&=f,&&\\mbox\{in\}\\kern 5\.0pt\\Omega\\times\(0,t\_\{f\}\],\\\\ u&=g,&&\\mbox\{on\}\\kern 5\.0pt\\partial\\Omega\\times\(0,t\_\{f\}\],\\\\ u&=u\_\{0\},&&\\mbox\{in\}\\kern 5\.0pt\\Omega,\\end\{aligned\}\\right\.\\end\{array\}\(38\)whereΩ=\(0,1\)\\Omega=\(0,1\)andtf=1t\_\{f\}=1\. We choose the suitablef⁡\(x,t\)f\(x,t\),g⁡\(x,t\)g\(x,t\)andu0​\(x\)u\_\{0\}\(x\)so that the exact solution is given by

u⁡\(t,x\)=11\+ex−t2​ϵ,u\(t,x\)=\\frac\{1\}\{1\+e^\{\\frac\{x\-t\}\{2\\epsilon\}\}\},\(39\)with a small value ofϵ=0\.01\\epsilon=0\.01\.

We useNx×Nt=200×200N\_\{x\}\\times N\_\{t\}=200\\times 200uniform collocation points to train the neural network\. Similarly to the one\-dimensional nonlinear Helmholtz equation, we employ Picard iteration to solve this problem, implementing100100iterations\. In Table[18](https://arxiv.org/html/2608.14733#S4.T18), we show the optimal scaling factors and the relevant parameters used during the search process\. In Figure[10](https://arxiv.org/html/2608.14733#S4.F10), we present heat maps of the exact solution, the approximate solution, and the absolute error between them\. The approximate solution is obtained using FEN with asin\\sinactivation, where the number of basis functions isM=5000M=5000and the optimal scaling factor isρo​p​t=112\\rho\_\{opt\}=112\. It can be observed that, while the approximate solution captures the details well, the error is noticeably larger compared to other problems due to the solution has large gradient in a local area\.

In Figure[11](https://arxiv.org/html/2608.14733#S4.F11), we present theL∞L\_\{\\infty\}error curves\. It can be observed that when the number of basis functionsM≤2500M\\leq 2500, none of the networks approximate the exact solution well\. However, whenM=5000M=5000, FENs outperform ELMs in approximating the exact solution\. Although the representational ability of ELMs has improved, the error achieved by FENs is still much lower than that of ELMs, indicating that ELMs are not suitable for solving this problem\. TheL∞L\_\{\\infty\}andL2L\_\{2\}errors listed in Table[19](https://arxiv.org/html/2608.14733#S4.T19)further support this conclusion\.\\chaddedIn addition, it can be observed that both FENs and ELMs exhibit a noticeable decrease in accuracy when solving this nonlinear Burgers’ equation, which is attributed to the presence of sharp variations in the solution\.

Table 18:Nonlinear Burgers’ equation: Parameters when solving the nonlinear Burgers’ equation \([38](https://arxiv.org/html/2608.14733#S4.E38)\)\.Activations\(ρm​i​n,ρm​a​x\]\(\\rho\_\{min\},\\rho\_\{max\}\]ρs\\rho\_\{s\}ρo​p​t\\rho\_\{opt\}M=900M=900M=1600M=1600M=2500M=2500M=5000M=5000sigmoid\(0,10\]\(0,10\]0\.16\.78\.15\.09\.7tanh\\tanh\(0,10\]\(0,10\]0\.16\.35\.73\.55\.9swish\(0,10\]\(0,10\]0\.19\.95\.79\.09\.6cos\\cos\(0,150\]\(0,150\]1456279111sin\\sin\(0,150\]\(0,150\]1444978112cos\\cos&\\&sin\\sin\(0,150\]\(0,150\]1506280114![Refer to caption](https://arxiv.org/html/2608.14733v1/figures/Burgers_heatmap3.png)

Figure 10:Nonlinear Burgers’ equation: the heat maps of the exact solution and the approximate solution of FEN with asin\\sinactivation, whereM=5000M=5000andρo​p​t=112\\rho\_\{opt\}=112\.![Refer to caption](https://arxiv.org/html/2608.14733v1/figures/error_Burgers_sharp.png)

Figure 11:Nonlinear Burgers’ equation:L∞L\_\{\\infty\}errors of neural networks when solving the nonlinear Burgers’ equation \([38](https://arxiv.org/html/2608.14733#S4.E38)\)\.Table 19:Nonlinear Burgers’ equation: Performance comparison of FENs and ELMs activated bysigmoid,tanh\\tanhandswishwhen solving the nonlinear Burgers’ equation \([38](https://arxiv.org/html/2608.14733#S4.E38)\)\. TheL∞L\_\{\\infty\}errors andL2L\_\{2\}errors for each model configuration are presented\.ActivationsM=900M=1600M=2500M=5000eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}sigmoid8\.8672E\-012\.0842E\-013\.6381E\-011\.3008E\-011\.1966E\-011\.6866E\-021\.1641E\-021\.1815E\-03tanh\\tanh9\.1862E\-012\.0184E\-013\.3211E\-011\.1212E\-018\.8108E\-022\.0741E\-021\.4296E\-021\.9656E\-03swish1\.1182E\+006\.4755E\-014\.7318E\-019\.8497E\-022\.0387E\-012\.9728E\-021\.2937E\-022\.9260E\-03cos\\cos5\.8175E\-011\.4085E\-012\.3469E\-015\.9826E\-024\.3138E\-025\.1853E\-039\.9221E\-044\.8742E\-05sin\\sin1\.2913E\+003\.7190E\-014\.6918E\-011\.3188E\-014\.1016E\-025\.6931E\-036\.2835E\-045\.4839E\-05cos\\cos&\\&sin\\sin6\.2089E\-013\.0694E\-012\.4240E\-016\.6556E\-024\.4349E\-025\.9443E\-038\.7493E\-046\.0509E\-05
### 4\.9High\-dimensional Poisson equation

We consider the high\-dimensional Poisson equation given by Equation \([40](https://arxiv.org/html/2608.14733#S4.E40)\)

\{−Δ​u=f⁡\(𝒙\),in​Ω,u=h⁡\(𝒙\),on​∂Ω,\\begin\{array\}\[\]\{r@\{\}l\}\\left\\\{\\begin\{aligned\} \-\\Delta u&=f\(\\boldsymbol\{x\}\),&&\\mbox\{in\}\\kern 5\.0pt\\Omega,\\\\ u&=h\(\\boldsymbol\{x\}\),&&\\mbox\{on\}\\kern 5\.0pt\\partial\\Omega,\\end\{aligned\}\\right\.\\end\{array\}\(40\)whereΩ=\(−1,1\)d\\Omega=\(\-1,1\)^\{d\}represents the spatial domain\.The exact solution to this equation is provided by

u⁡\(𝒙\)=\(1d​∑i=1dxi\)2\+sin⁡\(1d​∑i=1dxi\),u\(\\boldsymbol\{x\}\)=\\left\(\\frac\{1\}\{d\}\\sum\_\{i=1\}^\{d\}x\_\{i\}\\right\)^\{2\}\+\\sin\\left\(\\frac\{1\}\{d\}\\sum\_\{i=1\}^\{d\}x\_\{i\}\\right\),\(41\)with suitablef⁡\(𝒙\)f\(\\boldsymbol\{x\}\)and boundary conditionh⁡\(𝒙\)h\(\\boldsymbol\{x\}\)\.

For the training of both FENs and ELMs, we randomly select50,00050\{,\}000collocation points in the domainΩ\\Omegaand1000​d1000dcollocation points on the boundary∂Ω\\partial\\Omega\. Unlike in previous examples, in this problem we fix the number of basis functions toM=10,000M=10\{,\}000\. This is because, for high\-dimensional problems, a sufficient number of basis functions is essential to capture the complexity of the solution space\. To evaluate the capability of our methods in solving high\-dimensional problems, we solve the Poisson equation in dimensions55,77,1010, and1515\.

In Table[20](https://arxiv.org/html/2608.14733#S4.T20), we provide the parameters used for the optimal scale search, along with the optimal scaling factors obtained for problems of various dimensions\. Table[21](https://arxiv.org/html/2608.14733#S4.T21)summarizes theL∞L\_\{\\infty\}andL2L\_\{2\}errors achieved by the FENs and ELMs when using the optimal scaling factors\. Figure[12](https://arxiv.org/html/2608.14733#S4.F12)displays theL∞L\_\{\\infty\}error curves\.

The error distributions ‌reveal critical comparative information about‌ the dimensional scalability of FENs and ELMs\. Atd=5d=5, the errors of FENs are significantly lower than those of ELMs\. However, this trend changes as the dimensionality increases\. In higher dimensions, the ELM with thesigmoidactivation and the FEN with thesin\\sinactivation exhibit comparable precisions\. The ELM with thetanh\\tanhactivation shows the highest error, indicating its limited ability to achieve high precision\. Although the ELM with theswishactivation performs better than those withsigmoidandtanh\\tanhactivations atd=5d=5, its performance deteriorates in higher dimensions\. In these cases, its errors fall between the two and are comparable to those of FENs with thecos\\cosactivation, albeit still slightly higher than those of FENs with combinedcos\\cosandsin\\sinactivations\.

From the perspective of specific approximation error values, atd=5d=5, FENs achieve the smallestL∞L\_\{\\infty\}andL2L\_\{2\}errors of1\.2396×10−131\.2396\\times 10^\{\-13\}and3\.9641×10−143\.9641\\times 10^\{\-14\}, respectively\. In comparison, ELMs reach minimumL∞L\_\{\\infty\}andL2L\_\{2\}errors of1\.5994×10−121\.5994\\times 10^\{\-12\}and4\.5378×10−134\.5378\\times 10^\{\-13\}, respectively\. At higher dimensions \(d=7d=7,1010, and1515\), both FENs and ELMs achieveL∞L\_\{\\infty\}andL2L\_\{2\}errors of similar magnitudes: approximately10−910^\{\-9\}and10−1010^\{\-10\}ford=7d=7,10−610^\{\-6\}and10−710^\{\-7\}ford=10d=10, and10−510^\{\-5\}for both metrics atd=15d=15\. These numerical results highlight the performance of FENs and ELMs across different dimensions\. While FENs clearly outperform ELMs at lower dimensions in terms of approximation accuracy, the performance gap narrows as dimensionality increases, with both methods achieving comparable levels of precision\.

Table 20:High\-dimensional Poisson equation: Parameters when solving the high dimensional Poisson equation \([40](https://arxiv.org/html/2608.14733#S4.E40)\)\.Activations\(ρm​i​n,ρm​a​x\]\(\\rho\_\{min\},\\rho\_\{max\}\]ρs\\rho\_\{s\}ρo​p​t\\rho\_\{opt\}d=5d=5d=7d=7d=10d=10d=15d=15sigmoid\(0,1\]\(0,1\]0\.0010\.1430\.0490\.0270\.016tanh\\tanh\(0,1\]\(0,1\]0\.0010\.0460\.0150\.0140\.008swish\(0,1\]\(0,1\]0\.0010\.120\.0330\.0350\.024cos\\cos\(0,1\]\(0,1\]0\.0010\.250\.1110\.0910\.043sin\\sin\(0,1\]\(0,1\]0\.0010\.2620\.1020\.0520\.026cos\\cos&\\&sin\\sin\(0,1\]\(0,1\]0\.0010\.2650\.1050\.0690\.05![Refer to caption](https://arxiv.org/html/2608.14733v1/figures/error_PoissonHD.png)

Figure 12:High\-dimensional Poisson equation:L∞L\_\{\\infty\}errors of neural networks when solving the high dimensional Poisson equation \([40](https://arxiv.org/html/2608.14733#S4.E40)\)\.Table 21:High\-dimensional Poisson equation: Performance comparison of FENs and ELMs activated bysigmoid,tanh\\tanhandswishwhen solving the high dimensional Poisson equation \([35](https://arxiv.org/html/2608.14733#S4.E35)\)\. TheL∞L\_\{\\infty\}errors andL2L\_\{2\}errors for each model configuration are presented\.Activationsd=5d=7d=10d=15eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}eL∞e\_\{L\_\{\\infty\}\}eL2e\_\{L\_\{2\}\}sigmoid1\.8921E\-118\.7116E\-122\.7337E\-093\.8450E\-101\.0511E\-063\.1487E\-076\.2374E\-054\.0279E\-05tanh\\tanh9\.5997E\-101\.0016E\-101\.5291E\-081\.9115E\-092\.1064E\-058\.1063E\-064\.9823E\-043\.1662E\-04swish1\.5994E\-124\.5378E\-135\.2239E\-091\.0120E\-095\.5396E\-064\.2761E\-073\.3137E\-045\.9395E\-05cos\\cos1\.2396E\-133\.9641E\-145\.6010E\-097\.7473E\-105\.7360E\-063\.5996E\-073\.3140E\-045\.9459E\-05sin\\sin3\.6643E\-137\.9493E\-141\.5148E\-093\.1943E\-101\.1275E\-064\.4159E\-077\.9356E\-055\.2070E\-05cos\\cos&\\&sin\\sin2\.4564E\-135\.0408E\-141\.5341E\-093\.2616E\-102\.9906E\-064\.9836E\-072\.5601E\-043\.3998E\-05\\chadded

In Table 22, we report the computational time for solving the high\-dimensional Poisson equation using FENs and ELMs withtanh\\tanh,sigmoid, andswishactivation functions\. It can be observed that all models exhibit relatively low computational time, primarily because we do not rely on automatic differentiation but instead derive the derivatives of the basis functions analytically\. Notably, FENs achieve significantly lower runtime compared to ELMs\. This is because the derivatives of the basis functions in FENs withcos\\cos/sin\\sinactivations are analytically simple, whereas the derivatives of the basis functions in ELMs withtanh\\tanh,sigmoid, andswishactivations introduce considerable computational complexity\. Finally, it should be noted that in the searching for optimal scaling factor, if the range ofρ\\rhois too large or the step sizeρs\\rho\_\{s\}is too small, the overall computational cost can increase significantly\.

Table 22:High\-dimensional Poisson equation: Computational time for solving the high dimensional Poisson equation \([40](https://arxiv.org/html/2608.14733#S4.E40)\)\.Activationsd=5d=5d=7d=7d=10d=10d=15d=15sigmoid2\.38472\.62473\.13733\.7696tanh\\tanh2\.23542\.46112\.88023\.5063swish2\.68653\.11983\.80614\.9203cos\\cos2\.07132\.17672\.52193\.0101sin\\sin2\.01172\.21242\.65373\.1934cos\\cos&\\&sin\\sin1\.57861\.73452\.01902\.4217

## 5Conclusions

In this work, we propose Fourier Feature Networks \(FENs\) to study function approximation and the solution of linear and nonlinear PDEs\. These networks employ a single\-hidden\-layer neural network to represent the target function, where the outputs of the hidden layer correspond to a set of basis functions\. The linear combination of this set of basis functions yields the representation of the target function\. This concept is similar to that of ELMs, which also utilize a single\-hidden\-layer neural network\. However, the key difference is that, in ELMs, affine transformations are indispensable\. Without them, high\-precision solutions cannot be obtained for certain problems\. In contrast, FENs do not require affine transformations\. Additionally, while ELMs typically use activation functions such assigmoid,tanh\\tanh, andswish, FENs utilizecos\\cos,sin\\sin, or a combination of both, which naturally introduces Fourier features\.

We initialize the input\-to\-hidden weights and biases by random sampling from a uniform distribution \(variance=1=1\), which remain fixed during training\. Only the output layer’s linear combination coefficients require optimization\. To fully demonstrate the capability of neural networks in solving problems, the scaling factors for weights and biases are searched within a specified range\. The optimal scaling factors, which minimize the error in the algebraic equation, are identified\. Problems solved using these optimal scaling factors yield higher precision solutions\. In our numerical experiments, both for function approximation and solving linear or nonlinear PDEs, we observe that the solutions obtained by FENs are significantly more accurate than those obtained by ELMs usingsigmoid,tanh\\tanh, orswishactivation functions\. In high\-dimensional Poisson problems, FENs achieve higher precision solutions even in55dimensions, and for problems with even higher dimensions, FENs can perform comparably to ELMs\.

Although the proposed neural networks have shown promising results in numerical experiments, there are still areas that require further research\. The neural networks are highly dependent on the choice of scaling factors; poor selection can make it difficult to obtain a high\-precision solution\. Thus, researching more effective methods for determining the optimal scaling factor is necessary\. Both FENs and ELMs struggle with high\-dimensional problems, so the development of an algorithm capable of handling such problems with greater precision is essential\.\\chaddedFinally, we consider that analyzing the frequency\-domain superiority of trigonometric activation functions over tanh, sigmoid, and swish will constitute an important direction for future work\.

## Acknowledgment

This research is partially supported by the National Key R & D Program of China \(No\.2022YFE03040002\) and the National Natural Science Foundation of China \( No\.12371434\)\.

## Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request\.

## References

- \[1\]H\. Bararnia and M\. Esmaeilpour\(2022\)On the application of physics informed neural networks \(PINN\) to solve boundary layer thermal\-fluid problems\.International Communications in Heat and Mass Transfer132,pp\. 105890\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p1.1)\.
- \[2\]A\. G\. Baydin, B\. A\. Pearlmutter, A\. A\. Radul, and J\. M\. Siskind\(2018\)Automatic differentiation in machine learning: a survey\.Journal of Marchine Learning Research18,pp\. 1–43\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p3.1)\.
- \[3\]J\. Chen, X\. Chi, and Z\. Yang\(2022\)Bridging traditional and machine learning\-based algorithms for solving PDEs: the random feature method\.Journal of Machine Learning1\(3\),pp\. 268–298\.External Links:ISSN 2790\-2048,[Document](https://dx.doi.org/https%3A//doi.org/10.4208/jml.220726),[Link](http://global-sci.org/intro/article_detail/jml/21029.html)Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[4\]J\. Chen, W\. E, and Y\. Sun\(2024\)Optimization of random feature method in the high\-precision regime\.Communications on Applied Mathematics and Computation6\(2\),pp\. 1490–1517\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[5\]J\. Chen Y\. Luoet al\.\(2023\)The random feature method for time\-dependent problems\.arXiv preprint arXiv:2304\.06913\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[6\]M\. Cooley, V\. Shankar, R\. M\. Kirby, and S\. Zhe\(2024\)Fourier PINNs: from strong boundary conditions to adaptive Fourier bases\.arXiv preprint arXiv:2410\.03496\.Cited by:[§3\.1\.3](https://arxiv.org/html/2608.14733#S3.SS1.SSS3.p1.1)\.
- \[7\]S\. Cuomo, V\. S\. Di Cola, F\. Giampaolo, G\. Rozza, M\. Raissi, and F\. Piccialli\(2022\)Scientific machine learning through physics–informed neural networks: where we are and what’s next\.Journal of Scientific Computing92\(3\),pp\. 88\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p3.1)\.
- \[8\]F\. de Avila Belbute\-Peres and J\. Z\. KolterSimple initialization and parametrization of sinusoidal networks via their kernel bandwidth\.InThe Eleventh International Conference on Learning Representations,Cited by:[§3\.1\.2](https://arxiv.org/html/2608.14733#S3.SS1.SSS2.p1.1)\.
- \[9\]S\. Dong and Z\. Li\(2021\)Local extreme learning machines and domain decomposition for solving linear and nonlinear partial differential equations\.Computer Methods in Applied Mechanics and Engineering387,pp\. 114129\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p3.1)\.
- \[10\]S\. Dong and J\. Yang\(2022\)On computing the hyperparameter of extreme learning machines: algorithm and application to computational PDEs, and comparison with classical and high\-order finite elements\.Journal of Computational Physics463,pp\. 111290\.External Links:ISSN 0021\-9991,[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.jcp.2022.111290),[Link](https://www.sciencedirect.com/science/article/pii/S0021999122003527)Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p3.1),[§3\.3](https://arxiv.org/html/2608.14733#S3.SS3.p4.2)\.
- \[11\]Gallant and White\(1988\)There exists a neural network that does not make avoidable mistakes\.InIEEE 1988 International Conference on Neural Networks,Vol\.,pp\. 657–664 vol\.1\.External Links:[Document](https://dx.doi.org/10.1109/ICNN.1988.23903)Cited by:[§3\.1\.1](https://arxiv.org/html/2608.14733#S3.SS1.SSS1.p1.1)\.
- \[12\]X\. Glorot and Y\. Bengio\(2010\)Understanding the difficulty of training deep feedforward neural networks\.InProceedings of the thirteenth international conference on artificial intelligence and statistics,pp\. 249–256\.Cited by:[§3\.3](https://arxiv.org/html/2608.14733#S3.SS3.p1.1)\.
- \[13\]K\. He, X\. Zhang, S\. Ren, and J\. Sun\(2015\)Delving deep into rectifiers: surpassing human\-level performance on ImageNet classification\.InProceedings of the IEEE International Conference on Computer Vision \(ICCV\),Cited by:[§3\.3](https://arxiv.org/html/2608.14733#S3.SS3.p1.1)\.
- \[14\]H\. Hu, L\. Qi, and X\. Chao\(2024\)Physics\-informed neural networks \(PINN\) for computational solid mechanics: numerical frameworks and applications\.Thin\-Walled Structures,pp\. 112495\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p1.1)\.
- \[15\]B\. Huang and J\. Wang\(2022\)Applications of physics\-informed neural networks in power systems\-a review\.IEEE Transactions on Power Systems38\(1\),pp\. 572–588\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p1.1)\.
- \[16\]G\. Huang, Q\. Zhu, and C\. Siew\(2006\)Extreme learning machine: theory and applications\.Neurocomputing70\(1\-3\),pp\. 489–501\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p3.1)\.
- \[17\]L\. Huang, J\. Qin, Y\. Zhou, F\. Zhu, L\. Liu, and L\. Shao\(2023\)Normalization techniques in training dnns: methodology, analysis and application\.IEEE transactions on pattern analysis and machine intelligence45\(8\),pp\. 10173–10196\.Cited by:[§2\.3](https://arxiv.org/html/2608.14733#S2.SS3.p1.1)\.
- \[18\]D\. P\. Kingma\(2014\)Adam: a method for stochastic optimization\.arXiv preprint arXiv:1412\.6980\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[19\]N\. B\. Kovachki, S\. Lanthaler, and A\. M\. Stuart\(2024\)Operator learning: algorithms and analysis\.arXiv preprint arXiv:2402\.15715\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[20\]N\. Kovachki, Z\. Li, B\. Liu, K\. Azizzadenesheli, K\. Bhattacharya, A\. Stuart, and A\. Anandkumar\(2023\)Neural operator: learning maps between function spaces with applications to PDEs\.Journal of Machine Learning Research24\(89\),pp\. 1–97\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[21\]Z\. K\. Lawal, H\. Yassin, D\. T\. C\. Lai, and A\. Che Idris\(2022\)Physics\-informed neural network \(PINN\) evolution and beyond: a systematic literature review and bibliometric analysis\.Big Data and Cognitive Computing6\(4\),pp\. 140\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p1.1)\.
- \[22\]S\. Li, Y\. Xia, Y\. Liu, and Q\. Liao\(2023\)A deep domain decomposition method based on Fourier features\.Journal of Computational and Applied Mathematics423,pp\. 114963\.Cited by:[§3\.1\.3](https://arxiv.org/html/2608.14733#S3.SS1.SSS3.p1.1)\.
- \[23\]D\. C\. Liu and J\. Nocedal\(1989\)On the limited memory bfgs method for large scale optimization\.Mathematical programming45\(1\),pp\. 503–528\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[24\]S\. Liu\(2013\)Fourier neural network for machine learning\.In2013 international conference on machine learning and cybernetics,Vol\.1,pp\. 285–290\.Cited by:[§3\.1\.3](https://arxiv.org/html/2608.14733#S3.SS1.SSS3.p1.1)\.
- \[25\]L\. Lu, P\. Jin, G\. Pang, Z\. Zhang, and G\. E\. Karniadakis\(2021\)Learning nonlinear operators via deeponet based on the universal approximation theorem of operators\.Nature machine intelligence3\(3\),pp\. 218–229\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[26\]M\. Ngom and O\. Marin\(2021\)Fourier neural networks as function approximators and differential equation solvers\.Statistical Analysis and Data Mining: The ASA Data Science Journal14\(6\),pp\. 647–661\.Cited by:[§3\.1\.1](https://arxiv.org/html/2608.14733#S3.SS1.SSS1.p1.1)\.
- \[27\]N\. Ni and S\. Dong\(2023\)Numerical computation of partial differential equations by hidden\-layer concatenated extreme learning machine\.Journal of Scientific Computing95\(2\),pp\. 35\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[28\]M\. Raissi, P\. Perdikaris, and G\. E\. Karniadakis\(2019\)Physics\-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations\.Journal of Computational physics378,pp\. 686–707\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p1.1),[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[29\]Y\. Shang, F\. Wang, and J\. Sun\(2023\)Randomized neural network with Petrov–Galerkin methods for solving linear and nonlinear partial differential equations\.Communications in Nonlinear Science and Numerical Simulation127,pp\. 107518\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p2.1),[§3\.3](https://arxiv.org/html/2608.14733#S3.SS3.p5.1),[Theorem 3\.1](https://arxiv.org/html/2608.14733#S3.Thmtheorem1.p1.1.1)\.
- \[30\]Y\. Shang and F\. Wang\(2023\)Randomized neural networks with Petrov\-Galerkin methods for solving linear elasticity problems\.arXiv preprint arXiv:2308\.03088\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[31\]Y\. Shang and F\. Wang\(2024\)Randomized neural networks with Petrov–Galerkin methods for solving linear elasticity and navier–stokes equations\.Journal of Engineering Mechanics150\(4\),pp\. 04024010\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[32\]A\. Silvescu\(1999\)Fourier neural networks\.InIJCNN’99\. International Joint Conference on Neural Networks\. Proceedings \(Cat\. No\. 99CH36339\),Vol\.1,pp\. 488–491\.Cited by:[§3\.1\.1](https://arxiv.org/html/2608.14733#S3.SS1.SSS1.p1.1)\.
- \[33\]J\. Sirignano and K\. Spiliopoulos\(2018\)DGM: a deep learning algorithm for solving partial differential equations\.Journal of computational physics375,pp\. 1339–1364\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p1.1),[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[34\]V\. Sitzmann, J\. Martel, A\. Bergman, D\. Lindell, and G\. Wetzstein\(2020\)Implicit neural representations with periodic activation functions\.Advances in neural information processing systems33,pp\. 7462–7473\.Cited by:[§3\.1\.2](https://arxiv.org/html/2608.14733#S3.SS1.SSS2.p1.1)\.
- \[35\]M\. Tancik, P\. Srinivasan, B\. Mildenhall, S\. Fridovich\-Keil, N\. Raghavan, U\. Singhal, R\. Ramamoorthi, J\. Barron, and R\. Ng\(2020\)Fourier features let networks learn high frequency functions in low dimensional domains\.Advances in neural information processing systems33,pp\. 7537–7547\.Cited by:[§3\.1\.3](https://arxiv.org/html/2608.14733#S3.SS1.SSS3.p1.1)\.
- \[36\]F\. Wang and H\. Dang\(2024\)Randomized neural network methods for solving obstacle problems\.Banach Center Publications127,pp\. 261–276\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[37\]S\. Wang, H\. Wang, and P\. Perdikaris\(2021\)On the eigenvector bias of Fourier feature networks: from regression to solving multi\-scale PDEs with physics\-informed neural networks\.Computer Methods in Applied Mechanics and Engineering384,pp\. 113938\.Cited by:[§3\.1\.3](https://arxiv.org/html/2608.14733#S3.SS1.SSS3.p1.1)\.
- \[38\]B\. Yu and W\. E\(2018\)The Deep Ritz Method: a deep learning\-based numerical algorithm for solving variational problems\.Communications in Mathematics and Statistics6\(1\),pp\. 1–12\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p1.1),[§1](https://arxiv.org/html/2608.14733#S1.p2.1)\.
- \[39\]L\. Yuan, Y\. Ni, X\. Deng, and S\. Hao\(2022\)A\-PINN: auxiliary physics informed neural networks for forward and inverse problems of nonlinear integro\-differential equations\.Journal of Computational Physics462,pp\. 111260\.Cited by:[§1](https://arxiv.org/html/2608.14733#S1.p1.1)\.
- \[40\]A\. Zhumekenov, R\. Takhanov, A\. J\. Castro, and Z\. Assylbekov\(2021\)Approximation error of Fourier neural networks\.Statistical Analysis and Data Mining: The ASA Data Science Journal14\(3\),pp\. 258–270\.Cited by:[Theorem 3\.3](https://arxiv.org/html/2608.14733#S3.Thmtheorem3.p1.1.1)\.