Paper deep dive

A Novel Fourier Feature Network for Solving Partial Differential Equations

Qihong Yang, Zhijie Su, Yangtao Deng, Qiaolin He

Year: 2026Venue: arXiv preprintArea: cs.LGType: PreprintEmbeddings: 95

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%

Last extracted: 8/18/2026, 4:29:34 AM

Summary

The paper proposes Fourier Feature Networks (FENs) as a novel single-hidden-layer neural network architecture for solving Partial Differential Equations (PDEs). FENs utilize Fourier features via cosine, sine, or combined activations to generate basis functions, approximating the target function through a linear combination solved via the least squares method. Unlike Extreme Learning Machines (ELMs), FENs do not require affine transformations on input variables. The study introduces a method to search for an optimal scaling factor for randomly initialized weights and biases to maximize accuracy. Numerical experiments demonstrate that FENs consistently achieve higher accuracy than traditional ELMs using sigmoid, tanh, and swish activation functions.

Entities (10)

Partial Differential Equations · problem-domain · 99%Fourier Feature Network · model · 98%Extreme Learning Machine · model · 95%Least Squares Method · algorithm · 92%Scaling Factor · hyperparameter · 90%Cosine Activation · activation-function · 88%Sine Activation · activation-function · 88%Sigmoid · activation-function · 85%Swish · activation-function · 85%Tanh · activation-function · 85%

Relation Signals (9)

Fourier Feature Network → solves → Partial Differential Equations

confidence 95% · A Novel Fourier Feature Network for Solving Partial Differential Equations

Fourier Feature Network → uses → Least Squares Method

confidence 92% · The target function is then approximated as a linear combination of these basis functions, with the coefficients determined using the least squares method.

Fourier Feature Network → optimizes → Scaling Factor

confidence 90% · we search for an optimal scaling factor within a predefined range for the randomly initialized and fixed weights and biases.

Fourier Feature Network → outperforms → Extreme Learning Machine

confidence 90% · Our numerical experiments demonstrate that FENs consistently achieve higher accuracy than ELMs.

Fourier Feature Network → usesactivation → Sine Activation

confidence 88% · Fourier Feature Network with a sin activation

Fourier Feature Network → usesactivation → Cosine Activation

confidence 88% · Fourier Feature Network with a cos activation

Extreme Learning Machine → usesactivation → Sigmoid

confidence 85% · using various activation functions, such as sigmoid, tanh, and swish.

Extreme Learning Machine → usesactivation →

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Building on the foundation of single-hidden-layer neural networks, Fourier Feature Networks (FENs) are proposed, which incorporate Fourier features using $\cos$, $\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 as $\text{sigmoid}$, $\tanh$, and $\text{swish}$. Our numerical experiments demonstrate that FENs consistently achieve higher accuracy than ELMs.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

92,607 characters extracted from source content.

Expand or collapse full text

[style=chinese, orcid=0000-0002-8398-7212] [style=chinese] [style=chinese] [style=chinese] A Novel Fourier Feature Network for Solving Partial Differential Equations Qihong Yang yangqh0808@163.com Address: organization=School of Mathematics, Sichuan University, city=Chengdu, postcode=610065, country=China Zhijie Su zhijiesu164@163.com Yangtao Deng ytdeng1998@foxmail.com Qiaolin He qlhejenny@scu.edu.cn Abstract Building on the foundation of single-hidden-layer neural networks, Fourier Feature Networks (FENs) are proposed, which incorporate Fourier features using cos , 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 as sigmoid, tanh , and swish. 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 1 Introduction 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], the Deep Galerkin Method (DGM) [33], and Physics-Informed Neural Networks (PINNs) [28] has garnered significant attention. A substantial body of work [39, 21, 14, 15, 1] 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] or L-BFGS [23]. Methods like the DRM [38], the DGM [33], and the PINNs [28] all fall under this category. A closely related approach is operator learning [19, 20, 25], 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, 5, 3], the Randomized Neural Network with Petrov–Galerkin methods (RNN-PG) [29, 31, 30, 36], and the Hidden-Layer Concatenated Extreme Learning Machine (HLConcELM) [27]. Despite extensive research on training-based neural networks for solving PDEs, these methods often ‌exhibit limited accuracy‌ and typically require ‌considerable computational time [7]. 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. Learning Machine (ELM) is a single-layer feedforward neural network proposed in [16], which may not yield satisfactory results likely due to their reliance on Xavier or Kaiming initialization methods. Dong and Yang [10] presented a method for computing the optimal or near-optimal value of RmR_m based on the differential evolution algorithm in ELM. The RFM closely resembles Local Extreme Learning Machines (locELMs) [9], 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. 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], 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 on cos , 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 by sigmoid, tanh , and swish functions. 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 of sigmoid, tanh , and swish are introduced in Section 2. We propose Fourier Feature Networks, along with the three distinct activation modes of the FENs in Section 3. A series of numerical experiments designed to validate the efficacy of our methods are presented in Section 4. The article concludes with Section 5, where we summarize our results, discuss the implications of our work, and suggest directions for future research. 2 Preliminaries 2.1 The neural feature space Neural networks are widely recognized as nonlinear mapping functions that transform d-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 x ^d. The nonlinear mapping represented by the hidden layers of the neural network is given by ϕi() _i( x), where 1≤i≤M1≤ i≤ M, and M is the number of outputs of the last hidden layer. Consequently, the neural feature space, denoted as NNP_N, is defined as a space spanned by the basis functions ϕi _i, i.e., NN=spanϕ1,ϕ2,⋯,ϕM.P_N=span\ _1, _2,·s, _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( x)= _i=1^Mw_i _i. (2) 2.2 Extreme learning machine Assume that the network has d neurons in the input layer, corresponding to the input vector ∈ℝd x ^d. The main idea of the ELM is to randomly initialize the weights ∈ℝM×d W ^M× d and the biases ∈ℝM b ^M between the input and hidden layers. Then, the weights ∈ℝM w ^M from the hidden layer to the output layer are directly computed. This can be mathematically formulated as follows: uM()=Tσ(+),u_M( x)= w^Tσ( W x+ b), (3) where σ is the activation function, which acts elementwise on the vector + W x+ b, and uM()u_M( x) is the function represented by the ELM. We only need to calculate w. 2.3 Affine transformation As we know, normalization plays a vital role in enhancing model performance and generalization capabilities in neural network training [17]. 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 x∈ ^d, where Ω is a closed set. The affine transformation maps the input variables x to a new vector ~∈[−1,1]d⊂ℝd x∈[-1,1]^d ^d, which is given by Ti(xi)=2xi−xi(l)xi(u)−xi(l)−1,T_i(x_i)=2 x_i-x_i^(l)x_i^(u)-x_i^(l)-1, (4) where xix_i is the i-th component of x and xi(l)x_i^(l) and xi(u)x_i^(u) are the lower and upper bounds of xix_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( x)= w^Tσ( W x+ b)= w^Tσ ( W T( x)+ b ), (5) where =[T1(x1),T2(x2),⋯,Td(xd)]T T= [T_1(x_1),T_2(x_2),·s,T_d(x_d) ]^T. 3 Methodologies 3.1 Fourier Feature Networks It is observed that neural networks employing tanh activation 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.1 Fourier Feature Network with a cos activation Gallant and White [11] 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] implemented Fourier feature embedding via cosine-based activation operators. More recently, Ngom and Marin [26] 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 a cos activation. In this context, we also focus on a single-hidden-layer neural network. The target function can be represented as follows uM()=Tcos⁡(+).u_M( x)= w^T ( W x+ b). (6) It is noteworthy that, unlike networks activated by the sigmoid, tanh , or swish functions (as shown in Equation (5)), the input to the Fourier feature network with a cos activation does not require affine transformations. 3.1.2 Fourier Feature Network with a sin activation Sitzmann et al. [34] 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] 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()=Tsin⁡(+).u_M( x)= w^T ( W x+ b). (7) Similar to the Fourier feature network with a cos activation, the input to the network with a sin activation does not require affine transformations. 3.1.3 Fourier Feature Network with cos and sin activations In 2013, Liu [24] 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, 37, 22]. In line with Liu’s work, Fourier PINNs [6] 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 both cos and sin activations. The output function is defined as uM()=((1))Tcos⁡((1)+(1))+((2))Tsin⁡((2)+(2)),u_M( x)=( w^(1))^T ( W^(1) x+ b^(1))+( w^(2))^T ( W^(2) x+ b^(2)), (8) where (1)∈ℝp×d W^(1) ^p× d, (2)∈ℝp×d W^(2) ^p× d, (1)∈ℝp b^(1) ^p, (2)∈ℝp b^(2) ^p, (1)∈ℝp w^(1) ^p, (2)∈ℝp w^(2) ^p and M=2pM=2p is the number of basis functions. 3.2 Function approximation Let us assume that a real-valued matrix A of size N×MN× M is constructed from the basis functions Φ=[ϕ1,ϕ2,⋯,ϕM] =[ _1, _2,·s, _M], evaluated at discrete points within the dataset S. This matrix is defined by =[ϕ1(1)ϕ2(1)⋯ϕM(1)ϕ1(2)ϕ2(2)⋯ϕM(2)⋯ϕ1(N)ϕ2(N)⋯ϕM(N)],A= bmatrix _1( x_1)& _2( x_1)&·s& _M( x_1)\\ _1( x_2)& _2( x_2)&·s& _M( x_2)\\ & &·s& \\ _1( x_N)& _2( x_N)&·s& _M( x_N)\\ bmatrix, (9) where i∈S x_i∈ S and 1≤i≤N1≤ i≤ N. Based on Equation (2), the coefficient vector =[w1,w2,…,wM]T w=[w_1,w_2,...,w_M]^T can be determined by solving the linear system =A w=F using the least-squares method. Here, =Φ⁡(S)A= (S) denotes the basis function matrix formed by evaluating all basis functions at discrete data points in the dataset S, and =[f⁡(1),f⁡(2),…,f⁡(N)]TF=[f( x_1),f( x_2),…,f( x_N)]^T corresponds to the vector of function values at these sampled points. 3.3 Searching for optimal scaling factor The most widely used initialization methods in neural networks are Xavier initialization (also known as Glorot initialization) [12] and Kaiming initialization [13]. 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 W and bias 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 ρ. The functional mapping described by the neural network architecture in Equation (3) can then be reparameterized in terms of ρ as uM()=Tσ(ρ⁡(+)).u_M( x)= w^Tσ(ρ( W x+ b)). (10) The protocol for locating the optimal scaling factor ρopt _opt is formally specified in Algorithm 1. It is important to note that although the affine transformation is not explicitly shown in Equation (10), it must still be applied when using sigmoid, tanh , and swish as activation functions. It is important to note that the optimal scaling factor search method proposed in [10] 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] 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 ρ as (ρmin,ρmax]( _min, _max] and the step size ρs _s. Let M be the number of basis functions and ρ=ρmin+ρsρ= _min+ _s. Randomly initialize the weights of the networks. while ρ<ρmaxρ< _max do Obtain the basis functions Φ . Compute the matrix A. Construct the right-hand vector F. Express the current approximation of the solution as u⁡()=Φ⋅u( x)= · w. Solve the linear system ⋅=A· w=F using least squares method to obtain the vector of coefficients w. Let Loss=∥⋅−∥Loss= · w-F and record the LossLoss and the current ρ value. Let ρ=ρmin+ρsρ= _min+ _s. end while The optimal scale factor ρopt _opt is the one with the smallest LossLoss. Algorithm 1 Searching for optimal scaling factor Theorem 3.1. (Theorem 2.2 [29]) Given p≥1,s,k,d∈N+p≥ 1,s,k,d∈ N^+, s≥k+1s≥ k+1. Let σ be the logistic function or tanh function. Denote ℱs,p,d:=u∈s,p([0,1]d):∥u∥s,p([0,1]d)≤1 F_s,p,d:=\u ^s,p([0,1]^d): u _W^s,p([0,1]^d)≤ 1\, and σ(MD,BD):=uM(x)=Tσ(+) with ND≤MD, and |wij|≤BD,|bi|≤BD N_σ(M_D,B_D):=\u_M(x)= w^Tσ( W x+ b) with N_D≤ M_D, and w_ij ≤ B_D, b_i ≤ B_D\, where NDN_D is the number of non-zero parameter elements in the hidden layer of the neural network. For any ϵε ¿ 0 and u ∈ℱs,p,d∈ F_s,p,d, there exists a neural network uMu_M ∈σ(MD,BD)∈ N_σ(M_D,B_D) with MD≤C⋅ϵ−d/s−k−μkM_D≤ C·ε^-d/s-k-μ k, and BD≤C⋅ϵ−θB_D≤ C·ε^-θ such that ∥u−uM∥s,p([0,1]d)≤ϵ, u-u_M _W^s,p([0,1]^d)≤ε, (11) where C,θC,θ are constants depending on d,s,p,kd,s,p,k; μ is an arbitrarily small positive number. In Theorem 3.1, it is established that for any desired approximation accuracy ϵε, 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], these parameters are typically obtained through repeated random initialization, which are both computationally inefficient and potentially suboptimal. Theorem 3.2 demonstrates that, for any given initialization, it is possible to improve the applicability of the network parameters by appropriately adjusting the scaling factor ρ. Theorem 3.2. (Existence of the optimal scaling factor) Assume that u⁡()u( x) is a continuous function and σ is a continuous activation function. Let uM()=T(ρ)σ(ρ⁡(+))u_M( x)= w^T(ρ)σ(ρ( W x+ b)) and uMu_M is the projection of u onto the linear space spanned by the basis functions σ⁡(ρ⁡(+))σ(ρ( W x+ b)), where ⁡(ρ) w(ρ) is a vector consisting of a set of functions related to ρ. Then, ∀ρmax>0∀ _max>0, ∃ρopt∈[0,ρmax]∃ _opt∈[0, _max] such that ∥u−uM(,ρopt)∥=min0≤ρ≤ρmax⁡∥u⁡()−uM(,ρ)∥. u-u_M( x; _opt) = _0≤ρ≤ _max u( x)-u_M( x;ρ) . Proof. Let Φ be the set of basis functions σ⁡(ρ⁡(+))σ(ρ( W x+ b)), and S=[1,2,⋯,N]TS=[ x_1, x_2,·s, x_N]^T be the set of collocation points in the domain Ω . According to Equation (9), we substitute the set of collocation points S into it to obtain ⁡(ρ)(ρ)=,A(ρ) w(ρ)=F, where =[u⁡(1),u⁡(2),⋯,u⁡(N)]TF=[u( x_1),u( x_2),·s,u( x_N)]^T. The continuity of the activation function σ implies continuity of ⁡(ρ)A(ρ) in ρ. Given F is likewise continuous in ρ, it follows that both ⁡(ρ) w(ρ) and ultimately uM(,ρ)u_M( x;ρ) are continuous in ρ. Let e⁡(,ρ)=u⁡()−uM(,ρ)e( x;ρ)=u( x)-u_M( x;ρ), where e⁡(,ρ)e( x;ρ) is continuous in ρ. By the continuity of the norm, ∥e⁡(,ρ)∥ e( x;ρ) is likewise continuous in ρ. Consequently, ∀ρmax>0∀ _max>0, ∥e⁡(,ρ)∥ e( x;ρ) attains its maximum and minimum values on the closed interval [0,ρmax][0, _max]. ∎ Remark 1. In practice, the search for the optimal scaling factor typically excludes zero as an initial candidate, since empirical evidence indicates that ρ=0ρ=0 is rarely an optimal choice. Figure 1 shows how the L∞L_∞ error evolves with the scaling factor across different basis function counts. The ELM with the tanh activation function is employed to solve the two-dimensional Equation (25). The scaling factor ρ is varied within the range (0,10](0,10] with a step size of 0.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 ρ changes, reveals the existence of an optimal scaling factor (denoted as ρopt _opt) that minimizes the error. This observation numerically validates Theorem 3.2. It should be noted that although Figure 1 displays only the error curves of the ELM with tanh activation function under varying scaling factor ρ, a similar trend occurs in FENs. Moreover, FENs typically operate across a broader ρ search range and achieve lower errors. Figure 1: Error variation with the scaling factor ρ for different numbers of basis functions. 3.4 Solving linear PDEs Consider the following generic linear PDE ℒu=f,inΩ, =f,\ in 5.0pt , (12) ℬu=g,on∂Ω, =g,\ on 5.0pt∂ , (13) where u is the scalar field function to be approximated, ℒL is a linear PDE operator defined within the domain Ω and ℬB is a linear boundary operator acting on the boundary ∂Ω∂ . The functions f and g represent the source term and the boundary condition, respectively. Let Φ denote the basis functions generated by the outputs of the FENs. Applying the operators ℒL and ℬB to these basis functions yields ℒΦ=(ℒϕ1,ℒϕ2,⋯,ℒϕM), =(L _1,L _2,·s,L _M), (14) ℬΦ=(ℬϕ1,ℬϕ2,⋯,ℬϕM). =(B _1,B _2,·s,B _M). (15) Assume that the collocation dataset S consists of N=Nr+NbN=N_r+N_b points, with NrN_r points distributed within the domain Ω and NbN_b points allocated on the boundary ∂Ω∂ . These points are partitioned into two subsets: SrS_r (interior points) and SbS_b (boundary points). Under this configuration, the matrix A, which corresponds to the left-hand side of the PDE system, is constructed as follows =[ℒΦ(Sr)ℬΦ(Sb)] = bmatrixL (S_r)\\ B (S_b) 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)]. = bmatrixL _1( x_1)&L _2( x_1)&·s&L _M( x_1)\\ & &·s& \\ L _1( x_N_r)&L _2( x_N_r)&·s&L _M( x_N_r)\\ B _1( x_N_r+1)&B _2( x_N_r+1)&·s&B _M( x_N_r+1)\\ & &·s& \\ B _1( x_N_r+N_b)&B _2( x_N_r+N_b)&·s&B _M( x_N_r+N_b)\\ bmatrix. As described in the previous section, the coefficient vector w is obtained by solving the linear system =A w=F, where the right-hand side vector F is given by =[f⁡(1),⋯,f⁡(Nr),g⁡(Nr+1),⋯,g⁡(Nr+Nb)]T. F=[f( x_1),·s,f( x_N_r),g( x_N_r+1),·s,g( x_N_r+N_b)]^T. (17) Here, f⁡(i)f( x_i) for 1≤i≤Nr1≤ i≤ N_r corresponds to the source term evaluated at the ‌interior collocation points, whereas g⁡(i)g( x_i) for Nr+1≤i≤Nr+NbN_r+1≤ i≤ N_r+N_b represents the boundary conditions imposed at the ‌boundary collocation points. 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.5 Approxmation theory of FENs Let us consider functions defined on ℝdR^d that admit the following Fourier integral representation u⁡()=∫ℝderF~(dr),u( x)= _R^de^i W_r x F(\,d W_r), (18) where F~(dr)=eiθ(r)F(dr) F(\,d W_r)=e^iθ( W_r)F(\,d W_r) is a unique complex-valued measure (referred to as the Fourier distribution), with F⁡(dr)F(\,d W_r) representing the magnitude distribution, θ⁡(r)θ( W_r) denoting the phase, and r∈ℝ1×d W_r ^1× d. The class F of functions that admit such a representation with a finite total magnitude measure is defined as =u:ℝd→ℝ|∫ℝdF(dr)<∞.F= \u:R^d | _R^dF(\,d W_r)<∞ \. (19) For any C>0C>0, the class CF_C is defined as C=u:ℝd→ℝ|Cu=∫ℝdF(dr)<C.F_C= \u:R^d |C_u= _R^dF(\,d W_r)<C \. (20) We define the class ΩF_ as the set of functions defined on a bounded domain Ω that admit a Fourier representation of the form (18) (valid for ∈Ω x∈ ) with finite magnitude measure Ω=u:Ω→ℝ|∫ℝdF(dr)<∞.F_ = \u: | _R^dF(\,d W_r)<∞ \. (21) For each C>0C>0, we define C,Ω=u:Ω→ℝ|Cu,Ω=∫ℝdF(dr)≤C.F_C, = \u: |C_u, = _R^dF(\,d W_r)≤ C \. (22) Theorem 3.3. (Approximation error of FEN with a cos activation [40]) Let uM=∑i=1Mwicos⁡(i+bi)u_M= _i=1^Mw_i ( W_i x+b_i) for i∈ℝ1×d W_i ^1× d, bi∈ℝb_i , |wi|≤CM w_i ≤ CM, and ∈ℝd x ^d. For each function u∈C,Ωu _C, and any probability measure μ, there exists uMu_M (M≥1M≥ 1), such that ∫Ω|u⁡()−uM()|2μ()≤C2M. _ u( x)-u_M( x) ^2μ(\,d x)≤ C^2M. Theorem 3.4. (Approximation error of FEN with a sin activation) Let uM=∑i=1Mwisin⁡(i+bi)u_M= _i=1^Mw_i ( W_i x+b_i) for i∈ℝ1×d W_i ^1× d, bi∈ℝb_i , |wi|≤CM w_i ≤ CM and ∈ℝd x ^d. For each function u∈C,Ωu _C, and any probability measure μ, there exists uMu_M (M≥1M≥ 1), such that ∫Ω|u⁡()−uM()|2μ()≤C2M. _ u( x)-u_M( x) ^2μ(\,d x)≤ C^2M. According to Theorem 3.3, the proof of Theorem 3.4 is straightforward. Theorem 3.5. (Approxmation error of FEN with cos and sin activations) Let uM=∑i=1pwi(1)cos⁡(i(1)+bi(1))+∑i=1pwi(2)sin⁡(i(2)+bi(2))u_M= _i=1^pw_i^(1) ( W_i^(1) x+b_i^(1))+ _i=1^pw_i^(2) ( W_i^(2) x+b_i^(2)) for i(1)∈ℝ1×d W_i^(1) ^1× d, i(2)∈ℝ1×d W_i^(2) ^1× d, bi(1)∈ℝb_i^(1) , bi(2)∈ℝb_i^(2) , |wi(1)|≤Cp w_i^(1) ≤ Cp, |wi(2)|≤Cp w_i^(2) ≤ Cp, and ∈ℝd x ^d. For each function u∈C,Ωu _C, and any probability measure μ, there exists uMu_M (M=2p≥2M=2p≥ 2), such that ∫Ω|u⁡()−uM()|2μ()≤4C2p. _ u( x)-u_M( x) ^2μ(\,d x)≤ 4C^2p. Proof. From Theorem 3.3, for any given function u∈C,Ωu _C, , there is a sum up(1)=∑i=1pwi(1)cos⁡(i(1)+bi(1))u_p^(1)= _i=1^pw_i^(1) ( W_i^(1) x+b_i^(1)), such that ∫Ω|12u()−up(1)()|2μ()≤C2p. _ 12u( x)-u_p^(1)( x) ^2μ(\,d x)≤ C^2p. From Theorem 3.4, for any given function u∈C,Ωu _C, , there is a sum up(2)=∑i=1pwi(2)cos⁡(i(2)+bi(2))u_p^(2)= _i=1^pw_i^(2) ( W_i^(2) x+b_i^(2)), such that ∫Ω|12u()−up(2)()|2μ()≤C2p. _ 12u( x)-u_p^(2)( x) ^2μ(\,d x)≤ C^2p. Since uM=up(1)+up(2)u_M=u_p^(1)+u_p^(2), then ∫Ω|u⁡()−uM()|2μ(d) _ u( x)-u_M( x) ^2μ(\,d x) =∫Ω|12u()−up(1)()+12u()−up(2)()|2μ(d) = _ 12u( x)-u_p^(1)( x)+ 12u( x)-u_p^(2)( x) ^2μ(\,d x) ≤∫Ω2|12u()−up(1)()|2μ(d)+∫Ω2|12u()−up(2)()|2μ(d) ≤ _ 2 12u( x)-u_p^(1)( x) ^2μ(\,d x)+ _ 2 12u( x)-u_p^(2)( x) ^2μ(\,d x) ≤2C2p+2C2p=4C2p. ≤ 2C^2p+ 2C^2p= 4C^2p. ∎ 4 Numerical 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_∞ error) and the relative L2L_2 error are defined as follows: eL∞=max1≤i≤N⁡|uM(i)−uexact(i)|, e_L_∞= _1≤ i≤ N|u_M( x_i)-u_exact( x_i)|, (23) eL2=∑i=1N(uM(i)−uexact(i))2∑i=1N(uexact(i))2, e_L_2= _i=1^N(u_M( x_i)-u_exact( x_i))^2 _i=1^N(u_exact( x_i))^2, (24) where uMu_M and uexactu_exact represent the approximate and exact solutions, respectively, and i x_i (1≤i≤N1≤ i≤ N) are the collocation points for error evaluation. 4.1 Two-dimensional function To evaluate the capability of the FENs in function approximation, we implement them on a two-dimensional function defined in Equation (25) u⁡(x,y)=sin⁡(πx)sin⁡(4πy).u(x,y)= (π x) (4π y). (25) In this experiment, we utilize a uniformly distributed training grid with dimensions Nx×Ny=101×101N_x× N_y=101× 101, where NxN_x and NyN_y denote the number of points along the x-axis and y-axis, respectively. The number of basis functions M is chosen to be 400400, 900900, 16001600, and 25002500, respectively. For comparative analysis, we implement ELMs with three distinct activation functions sigmoid, tanh , and swish to address the same problem. Table 1 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, we present the approximation errors for FENs and ELMs (using sigmoid, tanh , and swish activations) across different basis function configurations (M=400,900,1600,2500M=400,900,1600,2500). As illustrated in Figure 2, the L∞L_∞ 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 lower L∞L_∞ error metrics. Quantitatively, when the number of basis functions is sufficient, the smallest L∞L_∞ and L2L_2 errors achieved by FENs are 6.4756×10−156.4756× 10^-15 and 1.4677×10−151.4677× 10^-15, respectively. In contrast, the smallest L∞L_∞ and L2L_2 errors for ELMs are 1.7599×10−101.7599× 10^-10 and 1.8013×10−111.8013× 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). Activations (ρmin,ρmax]( _min, _max] ρs _s ρopt _opt M=400M=400 M=900M=900 M=1600M=1600 M=2500M=2500 sigmoid (0,10](0,10] 0.01 1.52 2.19 2.94 3.81 tanh (0,10](0,10] 0.01 0.76 1.47 1.55 1.92 swish (0,10](0,10] 0.01 1.49 2.61 3.06 3.76 cos (0,50](0,50] 0.1 7.4 9.5 9.7 22.8 sin (0,50](0,50] 0.1 7.7 9.1 10.3 22.7 cos &\& sin (0,50](0,50] 0.1 9.0 9.6 9.6 12.6 Figure 2: Function approximation: L∞L_∞ errors of neural networks in approximating the two-dimensional function defined by Equation (25). Table 2: Function approximation: Performance comparison of FENs and ELMs with sigmoid, tanh , and swish activations in approximating the two-dimensional function defined by Equation (25). The L∞L_∞ and L2L_2 errors for each model configuration are presented. Activations M=400 M=900 M=1600 M=2500 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 sigmoid 3.2425E-05 9.2209E-06 6.3796E-08 6.2611E-09 4.0454E-09 2.7714E-10 1.7599E-10 1.8013E-11 tanh 2.7418E-05 6.1370E-06 1.2666E-07 1.0780E-08 3.1869E-09 2.4560E-10 1.9782E-10 1.8775E-11 swish 3.3855E-05 1.0940E-05 3.5856E-08 4.5318E-09 9.6770E-10 1.6892E-10 3.7835E-10 2.0718E-11 cos 1.8362E-07 7.4465E-08 1.2727E-14 2.3325E-15 9.6481E-15 1.6508E-15 3.3154E-14 3.5549E-15 sin 2.2316E-07 9.5567E-08 1.0184E-14 3.9475E-15 1.5099E-14 2.0114E-15 1.7819E-14 3.4593E-15 cos &\& sin 2.1903E-07 7.8576E-08 1.9601E-14 4.6683E-15 6.4756E-15 1.4677E-15 2.6069E-14 2.3770E-15 4.2 Helmholtz equation The Helmholtz equation in two dimensions is given by Equation (26) Δu+k2u=q,inΩ,u=h,on∂Ω, array[]r@l \ aligned u+k^2u&=q,&& in 5.0pt ,\\ u&=h,&& on 5.0pt∂ , aligned . array (26) where Ω=(0,1)2 =(0,1)^2. The exact solution is defined as u⁡(x,y)=sin⁡(a1πx)sin⁡(a2πy),u(x,y)= (a_1π x) (a_2π 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)= (k^2-(a_1π)^2-(a_2π)^2 ) (a_1π x) (a_2π y), (28) where the parameters are set as a1=1a_1=1, a2=4a_2=4, and k=1k=1. For training both FENs and ELMs, we employ a uniform grid of Nx×Ny=101×101N_x× N_y=101× 101 collocation points. This training dataset comprises both interior points within Ω and boundary points on ∂Ω∂ . Consequently, the dataset S is divided into two distinct subsets: SrS_r (containing interior points), and SbS_b (containing boundary points). Table 3 summarizes the optimal scaling factors and search parameters used in solving the Helmholtz equation. Table 4 presents the approximation errors of FENs and ELMs with varying numbers of basis functions, all utilizing the optimized scaling configuration. Figure 3 illustrates the L∞L_∞ 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 reveal that with a sufficient number of basis functions, FENs attain minimum L∞L_∞ and L2L_2 errors of 5.3300×10−145.3300× 10^-14 and 2.7200×10−142.7200× 10^-14, respectively. In contrast, the best performance from ELMs yields L∞L_∞ and L2L_2 errors of 4.9477×10−104.9477× 10^-10 and 2.7595×10−102.7595× 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) with solution in Equation (27). Activations (ρmin,ρmax]( _min, _max] ρs _s ρopt _opt M=400M=400 M=900M=900 M=1600M=1600 M=2500M=2500 sigmoid (0,10](0,10] 0.01 1.38 1.66 2.48 2.98 tanh (0,10](0,10] 0.01 0.6 0.94 1.1 1.5 swish (0,50](0,50] 0.01 1.4 2.0 2.2 2.4 cos (0,100](0,100] 0.1 5.6 9.3 9.0 13.1 sin (0,100](0,100] 0.1 5.8 9.9 9.9 12.0 cos &\& sin (0,100](0,100] 0.1 5.5 10.5 10.1 8.8 Figure 3: Helmholtz equation: L∞L_∞ errors of neural networks when solving Helmholtz equation (26) with solution in Equation (27). Table 4: Helmholtz equation: Performance comparison of FENs and ELMs activated by sigmoid, tanh and swish when solving the Helmholtz equation (26) with solution in Equation (27). The L∞L_∞ errors and L2L_2 errors for each model configuration are presented. Activations M=400 M=900 M=1600 M=2500 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 sigmoid 3.3875E-03 8.3152E-04 6.5565E-07 5.7295E-07 8.3819E-09 2.8206E-09 6.0754E-10 3.9546E-10 tanh 2.7466E-03 6.7760E-04 4.4703E-07 2.0869E-07 1.1059E-08 6.0855E-09 4.9477E-10 2.7595E-10 swish 6.3477E-03 1.5769E-03 7.8976E-07 2.2362E-07 6.6357E-09 3.0735E-09 1.0987E-09 3.9585E-10 cos 6.3797E-05 1.7419E-05 5.8978E-13 1.5736E-13 1.7333E-13 7.0401E-14 1.7977E-13 7.1273E-14 sin 1.2699E-04 1.8981E-05 9.0587E-13 2.7940E-13 1.8736E-13 7.4318E-14 5.3300E-14 2.7200E-14 cos &\& sin 1.2552E-04 2.2074E-05 1.5099E-12 2.4098E-13 1.4385E-13 6.6964E-14 3.7487E-13 1.3203E-13 For comparison, we also employ the finite element method (FEM) to solve the Helmholtz equation using a triangulation‌ with different mesh sizes h. The corresponding numerical results are presented in Table 5. 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. Remark 2. In this work, FENs employ a single-hidden-layer neural network with trigonometric activation functions as global basis functions, and solve PDEs using the collocation method. In contrast, the FEM uses local basis functions to approximate the solution piecewise, and is based on the variational principle to obtain a weak-form solution of the PDEs, whereas FENs aim to solve the strong form solution. While spectral methods also use global basis functions to approximate the strong form solution, they typically require the basis functions to be orthogonal. A key distinction between spectral methods and FENs lies in how they construct basis functions for multi-dimensional problems. For 2D or 3D PDEs, FENs can generate basis functions simply by randomly sampling the hidden-layer weights as in the 1D case, whereas spectral methods may generate basis functions using tensor products of 1D orthogonal basis functions. Table 5: Helmholtz equation: Results of FEM when solving the Helmholtz equation (26) with solution in Equation (27). The L∞L_∞ errors and L2L_2 errors for each model configuration are presented. Method h DoF eL∞e_L_∞ eL2e_L_2 FEM 0.1 1969 1.5408E-02 7.4130E-03 0.05 7577 1.5028E-03 8.9865E-04 0.02 46905 1.4795E-04 4.1033E-05 0.01 185913 1.8508E-05 5.0898E-06 0.005 741209 1.3809E-06 9.4240E-07 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), u⁡(x,y)=tanh⁡(xy),u(x,y)= (xy), (29) with the associated source term defined by q⁡(x,y)=(k2−2x2−2y2)tanh⁡(xy)+(2x2+2y2)tanh3⁡(xy).q(x,y)=(k^2-2x^2-2y^2) (xy)+(2x^2+2y^2) ^3(xy). (30) Unlike the solution in Equation (27), this exact function is closely related to the tanh activation function. Despite this, as shown in Table 7, FENs still outperform ELMs activated by non-trigonometric function, particularly when the number of basis functions is small. The minimum L∞L_∞ and L2L_2 errors achieved by FENs are 3.2613×10−153.2613× 10^-15 and 2.4689×10−152.4689× 10^-15, respectively, compared to 4.1078×10−154.1078× 10^-15 and 3.8454×10−153.8454× 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) with solution in Equation (29). Activations (ρmin,ρmax]( _min, _max] ρs _s ρopt _opt M=400M=400 M=900M=900 M=1600M=1600 M=2500M=2500 sigmoid (0,10](0,10] 0.01 1.51 2.27 2.47 3.31 tanh (0,10](0,10] 0.01 0.73 1.06 1.53 1.59 swish (0,50](0,50] 0.01 1.6 2.1 2.8 4.0 cos (0,100](0,100] 0.1 7.6 9.4 10.2 12.1 sin (0,100](0,100] 0.1 8.3 9.2 10.3 15.9 cos &\& sin (0,100](0,100] 0.1 7.5 9.7 11.2 23.6 Table 7: Helmholtz equation: Performance comparison of FENs and ELMs activated by different functions, when solving the Helmholtz equation (26) with solution (29). The L∞L_∞ errors and L2L_2 errors for each model configuration are presented. Methods M=400 M=900 M=1600 M=2500 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 ELMs (OPENsigmoid)sigmoid) 2.9184E-08 1.5916E-08 2.5498E-11 1.5287E-11 2.4958E-13 1.5588E-13 2.0040E-14 2.1436E-14 ELMs (OPENtanh)tanh) 2.3647E-08 1.2001E-08 1.6025E-11 1.0172E-11 1.8674E-13 1.5324E-13 1.1768E-14 9.1806E-15 ELMs (OPENswish)swish) 2.5343E-08 1.0115E-08 1.7768E-11 1.8566E-11 3.4794E-13 3.9641E-13 3.8497E-14 5.8373E-14 ELMs (sin ) 1.0436E-09 4.3592E-10 1.6043E-14 8.7617E-15 4.1078E-15 3.8454E-15 4.2188E-15 3.9039E-15 FENs (cos &\& sin ) 4.0931E-10 1.6929E-10 1.2101E-14 9.6805E-15 3.2613E-15 2.4689E-15 3.7192E-15 2.2473E-15 4.3 Diffusion 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], array[]r@l \ aligned ∂ u∂ t-ν ∂^2u∂ x^2&=f(x,t),&&(x,t)∈\ (a_1,b_1)×(0,t_f],\\ u(a_1,t)&=g_1(t),&&t∈(0,t_f],\\ u(b_1,t)&=g_2(t),&&t∈(0,t_f],\\ u(x,0)&=h(x),&&x∈[a_1,b_1],\\ aligned . array (31) where f⁡(x,t)f(x,t) denotes the source term, ν>0ν>0 denotes the constant diffusion coefficient, g1(t)g_1(t) and g2(t)g_2(t) prescribe the time-dependent Dirichlet boundary conditions, and h⁡(x)h(x) defines the initial condition. The parameters are chosen as a1=0a_1=0, b1=5b_1=5 and ν=0.01ν=0.01. The final time tft_f considered in simulations is set to be 11. We choose the suitable functions f⁡(x,t)f(x,t), g1(t)g_1(t), g2(t)g_2(t), and the initial condition h⁡(x)h(x) so that the exact solution is given by u⁡(x,t)=[2cos⁡(πx+π5)+32cos⁡(2πx−3π5)][2cos⁡(πt+π5)+32cos⁡(2πt−3π5)].u(x,t)= [2 (π x+ π5 )+ 32 (2π x- 3π5 ) ] [2 (π t+ π5 )+ 32 (2π t- 3π5 ) ]. (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 a Nx×Ny=101×101N_x× N_y=101× 101 collocation points. Table 8 lists the relevant parameters and the optimal scaling factor used in this setup. Table 9 compares the approximate errors of the solutions obtained by FENs and ELMs for different numbers of basis functions. Figure 4 illustrates the L∞L_∞ error curves, showing significantly lower errors for FENs than ELMs. Furthermore, by examining Table 9, we observe that the minimal L∞L_∞ and L2L_2 errors achieved by FENs are 9.2371×10−149.2371× 10^-14 and 5.9186×10−155.9186× 10^-15, respectively. In contrast, the minimal L∞L_∞ and L2L_2 errors obtained by ELMs are 5.7246×10−95.7246× 10^-9 and 3.8524×10−103.8524× 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) with tf=1t_f=1. Activations (ρmin,ρmax]( _min, _max] ρs _s ρopt _opt M=400M=400 M=900M=900 M=1600M=1600 M=2500M=2500 sigmoid (0,10](0,10] 0.01 1.55 1.86 2.61 3.35 tanh (0,10](0,10] 0.01 0.9 1.01 1.38 2.03 swish (0,50](0,50] 0.01 1.55 2.5 2.53 3.96 cos (0,100](0,100] 0.1 5.5 5.0 5.4 6.4 sin (0,100](0,100] 0.1 4.0 5.5 5.3 6.6 cos &\& sin (0,100](0,100] 0.1 5.6 5.2 5.4 5.3 Figure 4: Diffusion equation: the L∞L_∞ errors of neural networks when solving diffusion equation (31) with tf=1t_f=1. Table 9: Diffusion equation: Performance comparison of FENs and ELMs activated by sigmoid, tanh and swish when solving the diffusion equation (31) with tf=1t_f=1. The L∞L_∞ errors and L2L_2 errors for each model configuration are presented. Activations M=400 M=900 M=1600 M=2500 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 sigmoid 1.0276E-01 7.3530E-03 2.6306E-05 2.4873E-06 2.4815E-07 1.8386E-08 5.7246E-09 3.8524E-10 tanh 1.4518E-01 7.2468E-03 1.9460E-05 1.1581E-06 2.4414E-07 1.7353E-08 1.0454E-08 4.3533E-10 swish 6.3047E-02 4.5548E-03 1.8239E-05 1.5684E-06 1.1645E-07 6.2690E-09 1.1023E-08 9.1794E-10 cos 7.4218E-09 5.5655E-10 1.3012E-13 5.9433E-15 1.1724E-13 9.1437E-15 1.8852E-13 9.1456E-15 sin 1.4405E-08 8.5286E-10 9.2371E-14 5.9186E-15 1.0303E-13 8.3162E-15 2.1139E-13 6.9327E-15 cos &\& sin 9.6760E-09 1.3414E-10 2.1139E-13 1.1447E-14 1.1546E-13 1.3558E-14 9.7700E-14 8.3418E-15 4.4 Heat equation We consider the heat equation within the spatial-temporal domain Ω×(0,tf]=(0,1)2×(0,1] ×(0,t_f]=(0,1)^2×(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)∈Ω. array[]r@l \ aligned u_t(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)∈ .\\ aligned . array (33) We choose the suitable functions f⁡(x,y,t)f(x,y,t), g⁡(x,y,t)g(x,y,t), and the initial condition h⁡(x,y,t)h(x,y,t) so that the exact solution is given by u⁡(x,y,t)=2e−tsin⁡(π2x)sin⁡(π2y).u(x,y,t)=2e^-t ( π2x ) ( π2y ). For the training phase, a uniform grid of Nx×Ny×Nt=51×51×51N_x× N_y× N_t=51× 51× 51 collocation points is employed. The parameters and results from the optimal scale search are summarized in Table 10. Table 11 presents the L∞L_∞ and L2L_2 errors for both FENs and ELMs across varying numbers of basis functions. The L∞L_∞ error curves depicted in Figure 5 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 smallest L∞L_∞ and L2L_2 errors of 2.6090×10−142.6090× 10^-14 and 1.7562×10−141.7562× 10^-14, respectively. In contrast, ELMs reach their smallest L∞L_∞ and L2L_2 errors of 2.4443×10−122.4443× 10^-12 and 1.4761×10−121.4761× 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). Activations (ρmin,ρmax]( _min, _max] ρs _s ρopt _opt M=400M=400 M=900M=900 M=1600M=1600 M=2500M=2500 sigmoid (0,5](0,5] 0.01 0.15 0.3 0.46 0.56 tanh (0,5](0,5] 0.01 0.07 0.14 0.19 0.22 swish (0,5](0,5] 0.01 0.12 0.21 0.31 0.48 cos (0,50](0,50] 0.1 1.1 1.4 2.2 1.7 sin (0,50](0,50] 0.1 1.2 1.3 1.6 2.4 cos &\& sin (0,50](0,50] 0.1 1.1 1.3 1.6 2.0 Figure 5: Heat equation: L∞L_∞ errors of neural networks when solving the heat equation (33). Table 11: Heat equation: Performance comparison of FENs and ELMs activated by sigmoid, tanh and swish when solving the heat equation (33). The L∞L_∞ errors and L2L_2 errors for each model configuration are presented. Activations M=400 M=900 M=1600 M=2500 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 sigmoid 1.5832E-07 7.3817E-08 2.4738E-10 2.3643E-10 3.0013E-11 2.2705E-11 5.9936E-12 3.5015E-12 tanh 1.2910E-07 5.3362E-08 1.1386E-10 5.7319E-11 1.5007E-11 7.7614E-12 4.1442E-12 3.0037E-12 swish 1.4388E-07 7.5731E-08 1.4778E-10 7.2960E-11 2.0669E-11 8.1287E-12 2.4443E-12 1.4761E-12 cos 1.5015E-07 4.9863E-08 2.4385E-12 1.2761E-12 2.9932E-13 1.2883E-13 4.5741E-14 2.8907E-14 sin 6.6993E-08 3.8847E-08 2.8706E-12 1.1786E-12 1.2257E-13 9.7663E-14 2.6090E-14 1.7562E-14 cos &\& sin 6.4772E-08 2.3952E-08 2.1253E-12 9.1364E-13 2.1938E-13 1.3437E-13 3.1822E-14 3.0063E-14 4.5 Wave equation We consider the following wave equation defined over the spatial domain Ω=(0,1)2 =(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)∈Ω. array[]r@l \ aligned ∂^2u∂ t^2- 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)∈ .\\ aligned . array (34) The exact solution is given by u⁡(x,y,t)=sin⁡(π2x)sin⁡(π2y)sin⁡(π2t)u(x,y,t)= ( π2x ) ( π2y ) ( π2t ), with suitable boundary condition g⁡(x,y,t)g(x,y,t), initial conditions h⁡(x,y,t)h(x,y,t) and w⁡(x,y)w(x,y), as well as the source term f⁡(x,y,t)f(x,y,t).We employ a uniform grid of Nx×Ny×Nt=51×51×51N_x× N_y× N_t=51× 51× 51 collocation 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. The resulting errors obtained using FENs and ELMs with their respective optimal scaling factors are presented in Table 13. Figure 6 illustrates the L∞L_∞ and L2L_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, the minimum L∞L_∞ and L2L_2 errors achieved by FENs are 2.4425×10−152.4425× 10^-15 and 1.7659×10−151.7659× 10^-15, respectively. In contrast, minimum L∞L_∞ and L2L_2 errors achieved by ELMs are 1.4412×10−121.4412× 10^-12 and 5.1755×10−135.1755× 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). Activations (ρmin,ρmax]( _min, _max] ρs _s ρopt _opt M=400M=400 M=900M=900 M=1600M=1600 M=2500M=2500 sigmoid (0,5](0,5] 0.01 0.17 0.27 0.5 0.54 tanh (0,5](0,5] 0.01 0.06 0.13 0.18 0.24 swish (0,5](0,5] 0.01 0.13 0.24 0.35 0.49 cos (0,10](0,10] 0.1 1.28 1.34 1.46 1.7 sin (0,10](0,10] 0.1 1.3 1.38 1.47 1.39 cos &\& sin (0,10](0,10] 0.1 1.16 1.28 1.62 1.74 Figure 6: Wave equation: L∞L_∞ errors of neural networks when solving the wave equation (34). Table 13: Wave equation: Performance comparison of FENs and ELMs activated by sigmoid, tanh and swish when solving the wave equation (34). The L∞L_∞ errors and L2L_2 errors for each model configuration are presented. Activations M=400 M=900 M=1600 M=2500 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 sigmoid 2.4436E-07 6.4600E-08 1.9004E-10 1.3386E-10 2.0691E-11 7.1025E-12 2.1174E-12 1.0359E-12 tanh 1.7975E-07 3.5759E-08 8.8139E-11 4.3830E-11 8.2991E-12 3.5998E-12 1.6200E-12 5.1755E-13 swish 1.0241E-07 3.8084E-08 8.5493E-11 3.3717E-11 9.0949E-12 4.1869E-12 1.4412E-12 6.2451E-13 cos 3.7392E-08 1.6164E-08 2.1672E-13 1.2717E-13 5.2736E-15 2.3601E-15 2.4425E-15 1.7659E-15 sin 6.0666E-08 2.4112E-08 4.3354E-13 2.1115E-13 6.9944E-15 3.3327E-15 5.6621E-15 3.9232E-15 cos &\& sin 3.5241E-08 1.0682E-08 3.5172E-13 1.5920E-13 5.3291E-15 2.0016E-15 2.8588E-15 2.6000E-15 4.6 Nonlinear 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. array[]r@l \ aligned ∂^2u∂ x^2-λ u+β (u)&=f(x),&&x∈(a,b),\\ u(a)&=h_1,\\ u(b)&=h_2.\\ aligned . array (35) The exact solution is given by u⁡(x)=sin⁡(3πx+3π20)cos⁡(4πx−2π5)+32+x10u(x)= (3π x+ 3π20 ) (4π x- 2π5 )+ 32+ x10, with suitable boundary conditions and source term f⁡(x)f(x). The constant parameters in the equation are set as a=0a=0, b=8b=8, λ=50λ=50 and β=10β=10. In this nonlinear case, we use Nx=3000N_x=3000 uniform collocation points for training. To handle the nonlinearity of the problem, we adopt the Picard iteration method, performing a total of 100100 iterations. Initially, a vector of coefficients 0 w_0 is randomly initialized, generating‌ the approximate solution u0=Φ⋅0u_0= · w_0. The nonlinear term βsin⁡(u0)β (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 w_1. This iterative process continues, updating the coefficients at each step until convergence is achieved. The final approximate solution is expressed as u=Φ⋅u= · w, where 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. The resulting L∞L_∞ and L2L_2 errors from FENs and ELMs using the optimal scaling factor‌ are ‌presented‌ in Table 15. Figure 7 displays the L∞L_∞ 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 the tanh activation function ‌achieves‌ clearly ‌superior‌ performance compared to ‌networks using‌ sigmoid and swish, ‌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, with 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, we observe that the smallest L∞L_∞ and L2L_2 errors achieved by FENs are 2.2427×10−132.2427× 10^-13 and 1.3705×10−151.3705× 10^-15, respectively. In contrast, the lowest L∞L_∞ and L2L_2 errors obtained by ELMs are 3.7986×10−113.7986× 10^-11 and 4.4912×10−124.4912× 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). Activations (ρmin,ρmax]( _min, _max] ρs _s ρopt _opt M=400M=400 M=900M=900 M=1600M=1600 M=2500M=2500 sigmoid (0,10](0,10] 0.01 9.65 9.99 9.95 9.72 tanh (0,10](0,10] 0.01 9.65 9.71 9.9 9.98 swish (0,10](0,10] 0.01 10.0 9.94 9.94 9.84 cos (0,100](0,100] 1 16 33 34 64 sin (0,100](0,100] 1 20 34 24 74 cos &\& sin (0,100](0,100] 1 22 25 33 71 Figure 7: Nonlinear Helmholtz equation: L∞L_∞ errors of neural networks when solving the nonlinear Helmholtz equation (35). Table 15: Nonlinear Helmholtz equation: Performance comparison of FENs and ELMs activated by sigmoid, tanh and swish when solving the nonlinear Helmholtz equation (35). The L∞L_∞ errors and L2L_2 errors for each model configuration are presented. Activations M=400 M=900 M=1600 M=2500 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 sigmoid 3.4384E-04 1.6943E-05 5.2436E-06 2.4732E-07 1.5165E-06 7.1510E-08 1.4590E-06 7.7030E-08 tanh 2.0442E-06 1.0492E-07 1.0126E-08 4.7025E-10 3.7986E-11 4.4912E-12 5.6331E-10 3.1179E-11 swish 3.1724E-04 1.9411E-05 9.2887E-06 5.5102E-07 1.4125E-06 6.6785E-08 1.4330E-06 8.1166E-08 cos 2.9874E-12 1.5546E-13 2.2427E-13 1.3705E-14 9.2371E-13 4.8289E-14 7.0832E-13 3.5773E-14 sin 8.6642E-13 4.2264E-14 8.4621E-13 4.6449E-14 1.4064E-12 6.4684E-14 7.4474E-13 3.7450E-14 cos &\& sin 6.0973E-13 2.8716E-14 5.8753E-13 3.4123E-14 1.3309E-12 6.4362E-14 6.6214E-13 3.1905E-14 4.7 Poisson equation with an oscillating solution We consider the following one-dimensional Poisson equation: −Δu(x)=f⁡(x),inΩ,u⁡(x)=g⁡(x),on∂Ω, array[]r@l \ aligned - u(x)&=f(x),&& in 5.0pt ,\\ u(x)&=g(x),&& on 5.0pt∂ , aligned . array (36) where Ω=(0,1) =(0,1). We choose the suitable f⁡(x)f(x) and boundary condition g⁡(x)g(x) so that the exact solution is given by u⁡(x)=16∑i=16sin⁡(2iπx),u(x)= 16 _i=1^6 (2^iπ x), (37) which represents a superposition of sine functions with exponentially increasing frequencies, producing a highly oscillatory ‌behavior. We train the networks using Nx=3000N_x=3000 uniformly 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. Notably, the admissible range for the activation function tanh is narrower than those of sigmoid and swish. This constraint ‌originates from‌ numerical stability requirements: ‌applying‌ excessively large scaling factors to tanh ‌generates‌ ‌severely‌ ill-conditioned matrices (‌approaching singularity‌), ‌which destabilizes the solving process. ‌To mitigate this, the tanh scaling range is ‌strategically restricted‌ to maintain ‌trainable system ‌conditions. Table 16 compiles the optimal scaling factors ‌and their associated search parameters. Furthermore, Figure 8 shows a comparison between the exact solution and the numerical solution obtained using the FEN with the cos activation function, where the number of basis functions is set to M=900M=900 and the optimal scaling factor is ρopt=130 _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 quantifies the L∞L_∞ convergence of neural networks ‌equipped with‌ optimal scaling factors for solving the ‌highly oscillatory‌ Poisson equation (37). It is clearly observed that the errors of ELMs with sigmoid, swish, and tanh activation functions are significantly higher than those of FENs. In Table 17, we report the L∞L_∞ and L2L_2 errors corresponding to various numbers of basis functions. The minimum L∞L_∞ and L2L_2 errors achieved by ELMs are 6.1572×10−86.1572× 10^-8 and 8.0396×10−88.0396× 10^-8, respectively, whereas FENs achieve notably smaller errors of 1.3878×10−111.3878× 10^-11 and 3.0395×10−113.0395× 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). Activations (ρmin,ρmax]( _min, _max] ρs _s ρopt _opt M=400M=400 M=900M=900 M=1600M=1600 M=2500M=2500 sigmoid (0,20](0,20] 0.01 18.79 18.12 18.46 19.19 tanh (0,10](0,10] 0.01 9.15 9.83 9.93 9.70 swish (0,20](0,20] 0.01 19.90 19.33 19.39 18.95 cos (0,1000](0,1000] 1 176 130 205 564 sin (0,1000](0,1000] 1 147 150 234 890 cos &\& sin (0,1000](0,1000] 1 164 146 325 903 Figure 8: Poisson equation with an oscillating solution: comparison between the exact solution and the approximate solution obtained using FEN with a cos activation, where M=900M=900 and ρopt=130 _opt=130. Figure 9: Poisson equation with an oscillating solution: L∞L_∞ errors of neural networks when solving the Poisson equation (36). Table 17: Poisson equation with an oscillating solution: Performance comparison of FENs and ELMs activated by sigmoid, tanh and swish when solving the Poisson equation (36). The L∞L_∞ errors and L2L_2 errors for each model configuration are presented. Activations M=400 M=900 M=1600 M=2500 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 sigmoid 3.6288E-05 3.4848E-05 2.9000E-07 6.0888E-07 1.6587E-07 1.8036E-07 1.0117E-06 1.0275E-06 tanh 3.0062E-05 3.8163E-07 2.6801E-07 2.8331E-07 6.1572E-08 8.0396E-08 2.0924E-07 1.7937E-07 swish 1.7748E-04 2.4055E-04 1.4265E-06 1.5708E-06 5.9915E-07 8.1065E-07 1.9759E-06 2.1461E-06 cos 3.0724E-11 9.7238E-11 1.3878E-11 3.0395E-11 7.2102E-11 9.4554E-11 1.7749E-11 3.2850E-11 sin 3.2458E-11 5.7302E-11 2.9600E-11 4.9456E-11 2.7607E-10 4.0441E-10 4.6233E-11 9.1841E-11 cos &\& sin 1.9459E-11 3.0421E-11 6.1339E-11 1.1078E-10 1.7186E-11 2.8037E-11 1.9275E-11 3.1528E-11 4.8 Nonlinear Burgers’ equation Consider the following 1D Burgers’ equation ut+uux−ϵuxx=f,inΩ×(0,tf],u=g,on∂Ω×(0,tf],u=u0,inΩ, array[]r@l \ aligned u_t+u_x-ε u_x&=f,&& in 5.0pt ×(0,t_f],\\ u&=g,&& on 5.0pt∂ ×(0,t_f],\\ u&=u_0,&& in 5.0pt , aligned . array (38) where Ω=(0,1) =(0,1) and tf=1t_f=1. We choose the suitable f⁡(x,t)f(x,t), g⁡(x,t)g(x,t) and u0(x)u_0(x) so that the exact solution is given by u⁡(t,x)=11+ex−t2ϵ,u(t,x)= 11+e x-t2ε, (39) with a small value of ϵ=0.01ε=0.01. We use Nx×Nt=200×200N_x× N_t=200× 200 uniform collocation points to train the neural network. Similarly to the one-dimensional nonlinear Helmholtz equation, we employ Picard iteration to solve this problem, implementing 100100 iterations. In Table 18, we show the optimal scaling factors and the relevant parameters used during the search process. In Figure 10, 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 a sin activation, where the number of basis functions is M=5000M=5000 and the optimal scaling factor is ρopt=112 _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, we present the L∞L_∞ error curves. It can be observed that when the number of basis functions M≤2500M≤ 2500, none of the networks approximate the exact solution well. However, when M=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. The L∞L_∞ and L2L_2 errors listed in Table 19 further support this conclusion. 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). Activations (ρmin,ρmax]( _min, _max] ρs _s ρopt _opt M=900M=900 M=1600M=1600 M=2500M=2500 M=5000M=5000 sigmoid (0,10](0,10] 0.1 6.7 8.1 5.0 9.7 tanh (0,10](0,10] 0.1 6.3 5.7 3.5 5.9 swish (0,10](0,10] 0.1 9.9 5.7 9.0 9.6 cos (0,150](0,150] 1 45 62 79 111 sin (0,150](0,150] 1 44 49 78 112 cos &\& sin (0,150](0,150] 1 50 62 80 114 Figure 10: Nonlinear Burgers’ equation: the heat maps of the exact solution and the approximate solution of FEN with a sin activation, where M=5000M=5000 and ρopt=112 _opt=112. Figure 11: Nonlinear Burgers’ equation: L∞L_∞ errors of neural networks when solving the nonlinear Burgers’ equation (38). Table 19: Nonlinear Burgers’ equation: Performance comparison of FENs and ELMs activated by sigmoid, tanh and swish when solving the nonlinear Burgers’ equation (38). The L∞L_∞ errors and L2L_2 errors for each model configuration are presented. Activations M=900 M=1600 M=2500 M=5000 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 sigmoid 8.8672E-01 2.0842E-01 3.6381E-01 1.3008E-01 1.1966E-01 1.6866E-02 1.1641E-02 1.1815E-03 tanh 9.1862E-01 2.0184E-01 3.3211E-01 1.1212E-01 8.8108E-02 2.0741E-02 1.4296E-02 1.9656E-03 swish 1.1182E+00 6.4755E-01 4.7318E-01 9.8497E-02 2.0387E-01 2.9728E-02 1.2937E-02 2.9260E-03 cos 5.8175E-01 1.4085E-01 2.3469E-01 5.9826E-02 4.3138E-02 5.1853E-03 9.9221E-04 4.8742E-05 sin 1.2913E+00 3.7190E-01 4.6918E-01 1.3188E-01 4.1016E-02 5.6931E-03 6.2835E-04 5.4839E-05 cos &\& sin 6.2089E-01 3.0694E-01 2.4240E-01 6.6556E-02 4.4349E-02 5.9443E-03 8.7493E-04 6.0509E-05 4.9 High-dimensional Poisson equation We consider the high-dimensional Poisson equation given by Equation (40) −Δu=f⁡(),inΩ,u=h⁡(),on∂Ω, array[]r@l \ aligned - u&=f( x),&& in 5.0pt ,\\ u&=h( x),&& on 5.0pt∂ , aligned . array (40) where Ω=(−1,1)d =(-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( x)= ( 1d _i=1^dx_i )^2+ ( 1d _i=1^dx_i ), (41) with suitable f⁡()f( x) and boundary condition h⁡()h( x). For the training of both FENs and ELMs, we randomly select 50,00050,000 collocation points in the domain Ω and 1000d1000d collocation points on the boundary ∂Ω∂ . Unlike in previous examples, in this problem we fix the number of basis functions to M=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 dimensions 55, 77, 1010, and 1515. In Table 20, we provide the parameters used for the optimal scale search, along with the optimal scaling factors obtained for problems of various dimensions. Table 21 summarizes the L∞L_∞ and L2L_2 errors achieved by the FENs and ELMs when using the optimal scaling factors. Figure 12 displays the L∞L_∞ error curves. The error distributions ‌reveal critical comparative information about‌ the dimensional scalability of FENs and ELMs. At d=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 the sigmoid activation and the FEN with the sin activation exhibit comparable precisions. The ELM with the tanh activation shows the highest error, indicating its limited ability to achieve high precision. Although the ELM with the swish activation performs better than those with sigmoid and tanh activations at d=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 the cos activation, albeit still slightly higher than those of FENs with combined cos and sin activations. From the perspective of specific approximation error values, at d=5d=5, FENs achieve the smallest L∞L_∞ and L2L_2 errors of 1.2396×10−131.2396× 10^-13 and 3.9641×10−143.9641× 10^-14, respectively. In comparison, ELMs reach minimum L∞L_∞ and L2L_2 errors of 1.5994×10−121.5994× 10^-12 and 4.5378×10−134.5378× 10^-13, respectively. At higher dimensions (d=7d=7, 1010, and 1515), both FENs and ELMs achieve L∞L_∞ and L2L_2 errors of similar magnitudes: approximately 10−910^-9 and 10−1010^-10 for d=7d=7, 10−610^-6 and 10−710^-7 for d=10d=10, and 10−510^-5 for both metrics at d=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). Activations (ρmin,ρmax]( _min, _max] ρs _s ρopt _opt d=5d=5 d=7d=7 d=10d=10 d=15d=15 sigmoid (0,1](0,1] 0.001 0.143 0.049 0.027 0.016 tanh (0,1](0,1] 0.001 0.046 0.015 0.014 0.008 swish (0,1](0,1] 0.001 0.12 0.033 0.035 0.024 cos (0,1](0,1] 0.001 0.25 0.111 0.091 0.043 sin (0,1](0,1] 0.001 0.262 0.102 0.052 0.026 cos &\& sin (0,1](0,1] 0.001 0.265 0.105 0.069 0.05 Figure 12: High-dimensional Poisson equation: L∞L_∞ errors of neural networks when solving the high dimensional Poisson equation (40). Table 21: High-dimensional Poisson equation: Performance comparison of FENs and ELMs activated by sigmoid, tanh and swish when solving the high dimensional Poisson equation (35). The L∞L_∞ errors and L2L_2 errors for each model configuration are presented. Activations d=5 d=7 d=10 d=15 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 eL∞e_L_∞ eL2e_L_2 sigmoid 1.8921E-11 8.7116E-12 2.7337E-09 3.8450E-10 1.0511E-06 3.1487E-07 6.2374E-05 4.0279E-05 tanh 9.5997E-10 1.0016E-10 1.5291E-08 1.9115E-09 2.1064E-05 8.1063E-06 4.9823E-04 3.1662E-04 swish 1.5994E-12 4.5378E-13 5.2239E-09 1.0120E-09 5.5396E-06 4.2761E-07 3.3137E-04 5.9395E-05 cos 1.2396E-13 3.9641E-14 5.6010E-09 7.7473E-10 5.7360E-06 3.5996E-07 3.3140E-04 5.9459E-05 sin 3.6643E-13 7.9493E-14 1.5148E-09 3.1943E-10 1.1275E-06 4.4159E-07 7.9356E-05 5.2070E-05 cos &\& sin 2.4564E-13 5.0408E-14 1.5341E-09 3.2616E-10 2.9906E-06 4.9836E-07 2.5601E-04 3.3998E-05 In Table 22, we report the computational time for solving the high-dimensional Poisson equation using FENs and ELMs with tanh , sigmoid, and swish activation 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 with cos /sin activations are analytically simple, whereas the derivatives of the basis functions in ELMs with tanh , sigmoid, and swish activations introduce considerable computational complexity. Finally, it should be noted that in the searching for optimal scaling factor, if the range of ρ is too large or the step size ρs _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). Activations d=5d=5 d=7d=7 d=10d=10 d=15d=15 sigmoid 2.3847 2.6247 3.1373 3.7696 tanh 2.2354 2.4611 2.8802 3.5063 swish 2.6865 3.1198 3.8061 4.9203 cos 2.0713 2.1767 2.5219 3.0101 sin 2.0117 2.2124 2.6537 3.1934 cos &\& sin 1.5786 1.7345 2.0190 2.4217 5 Conclusions 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 as sigmoid, tanh , and swish, FENs utilize cos , 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 using sigmoid, tanh , or swish activation functions. In high-dimensional Poisson problems, FENs achieve higher precision solutions even in 55 dimensions, 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. , 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 Transfer 132, p. 105890. Cited by: §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 Research 18, p. 1–43. Cited by: §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 Learning 1 (3), p. 268–298. External Links: ISSN 2790-2048, Document, Link Cited by: §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 Computation 6 (2), p. 1490–1517. Cited by: §1. [5] J. Chen Y. Luo et al. (2023) The random feature method for time-dependent problems. arXiv preprint arXiv:2304.06913. Cited by: §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. [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 Computing 92 (3), p. 88. Cited by: §1. [8] F. de Avila Belbute-Peres and J. Z. Kolter Simple initialization and parametrization of sinusoidal networks via their kernel bandwidth. In The Eleventh International Conference on Learning Representations, Cited by: §3.1.2. [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 Engineering 387, p. 114129. Cited by: §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 Physics 463, p. 111290. External Links: ISSN 0021-9991, Document, Link Cited by: §1, §3.3. [11] Gallant and White (1988) There exists a neural network that does not make avoidable mistakes. In IEEE 1988 International Conference on Neural Networks, Vol. , p. 657–664 vol.1. External Links: Document Cited by: §3.1.1. [12] X. Glorot and Y. Bengio (2010) Understanding the difficulty of training deep feedforward neural networks. In Proceedings of the thirteenth international conference on artificial intelligence and statistics, p. 249–256. Cited by: §3.3. [13] K. He, X. Zhang, S. Ren, and J. Sun (2015) Delving deep into rectifiers: surpassing human-level performance on ImageNet classification. In Proceedings of the IEEE International Conference on Computer Vision (ICCV), Cited by: §3.3. [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, p. 112495. Cited by: §1. [15] B. Huang and J. Wang (2022) Applications of physics-informed neural networks in power systems-a review. IEEE Transactions on Power Systems 38 (1), p. 572–588. Cited by: §1. [16] G. Huang, Q. Zhu, and C. Siew (2006) Extreme learning machine: theory and applications. Neurocomputing 70 (1-3), p. 489–501. Cited by: §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 intelligence 45 (8), p. 10173–10196. Cited by: §2.3. [18] D. P. Kingma (2014) Adam: a method for stochastic optimization. arXiv preprint arXiv:1412.6980. Cited by: §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. [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 Research 24 (89), p. 1–97. Cited by: §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 Computing 6 (4), p. 140. Cited by: §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 Mathematics 423, p. 114963. Cited by: §3.1.3. [23] D. C. Liu and J. Nocedal (1989) On the limited memory bfgs method for large scale optimization. Mathematical programming 45 (1), p. 503–528. Cited by: §1. [24] S. Liu (2013) Fourier neural network for machine learning. In 2013 international conference on machine learning and cybernetics, Vol. 1, p. 285–290. Cited by: §3.1.3. [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 intelligence 3 (3), p. 218–229. Cited by: §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 Journal 14 (6), p. 647–661. Cited by: §3.1.1. [27] N. Ni and S. Dong (2023) Numerical computation of partial differential equations by hidden-layer concatenated extreme learning machine. Journal of Scientific Computing 95 (2), p. 35. Cited by: §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 physics 378, p. 686–707. Cited by: §1, §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 Simulation 127, p. 107518. Cited by: §1, §3.3, Theorem 3.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. [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 Mechanics 150 (4), p. 04024010. Cited by: §1. [32] A. Silvescu (1999) Fourier neural networks. In IJCNN’99. International Joint Conference on Neural Networks. Proceedings (Cat. No. 99CH36339), Vol. 1, p. 488–491. Cited by: §3.1.1. [33] J. Sirignano and K. Spiliopoulos (2018) DGM: a deep learning algorithm for solving partial differential equations. Journal of computational physics 375, p. 1339–1364. Cited by: §1, §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 systems 33, p. 7462–7473. Cited by: §3.1.2. [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 systems 33, p. 7537–7547. Cited by: §3.1.3. [36] F. Wang and H. Dang (2024) Randomized neural network methods for solving obstacle problems. Banach Center Publications 127, p. 261–276. Cited by: §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 Engineering 384, p. 113938. Cited by: §3.1.3. [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 Statistics 6 (1), p. 1–12. Cited by: §1, §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 Physics 462, p. 111260. Cited by: §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 Journal 14 (3), p. 258–270. Cited by: Theorem 3.3.

Instant research discovery