Solving a class of multi-scale elliptic PDEs by means of Fourier-based
  mixed physics informed neural networks
Deep neural networks have garnered widespread attention due to their
simplicity and flexibility in the fields of engineering and scientific
calculation. In this study, we probe into solving a class of elliptic partial
differential equations(PDEs) with multiple scales by utilizing Fourier-based
mixed physics informed neural networks(dubbed FMPINN), its solver is configured
as a multi-scale deep neural network. In contrast to the classical PINN method,
a dual (flux) variable about the rough coefficient of PDEs is introduced to
avoid the ill-condition of neural tangent kernel matrix caused by the
oscillating coefficient of multi-scale PDEs. Therefore, apart from the physical
conservation laws, the discrepancy between the auxiliary variables and the
gradients of multi-scale coefficients is incorporated into the cost function,
then obtaining a satisfactory solution of PDEs by minimizing the defined loss
through some optimization methods. Additionally, a trigonometric activation
function is introduced for FMPINN, which is suited for representing the
derivatives of complex target functions. Handling the input data by Fourier
feature mapping will effectively improve the capacity of deep neural networks
to solve high-frequency problems. Finally, to validate the efficiency and
robustness of the proposed FMPINN algorithm, we present several numerical
examples of multi-scale problems in various dimensional Euclidean spaces. These
examples cover both low-frequency and high-frequency oscillation cases,
demonstrating the effectiveness of our approach. All code and data accompanying
this manuscript will be made publicly available at
\href{https://github.com/Blue-Giant/FMPINN}{https://github.com/Blue-Giant/FMPINN}.
Xi'an Li, Jinran Wu, You-Gan Wang, Xin Tai, Jianhua Xu
Solving a class of multi-scale elliptic PDEs by Fourier-based mixed physics informed neural networks

###### Abstract

Deep neural networks have garnered widespread attention due to their simplicity and flexibility in the fields of engineering and scientific calculation. In this study, we probe into solving a class of elliptic partial differential equations(PDEs) with multiple scales by utilizing Fourier-based mixed physics informed neural networks(dubbed FMPINN), its solver is configured as a multi-scale deep neural network. In contrast to the classical PINN method, a dual (flux) variable about the rough coefficient of PDEs is introduced to avoid the ill-condition of neural tangent kernel matrix caused by the oscillating coefficient of multi-scale PDEs. Therefore, apart from the physical conservation laws, the discrepancy between the auxiliary variables and the gradients of multi-scale coefficients is incorporated into the cost function, then obtaining a satisfactory solution of PDEs by minimizing the defined loss through some optimization methods. Additionally, a trigonometric activation function is introduced for FMPINN, which is suited for representing the derivatives of complex target functions. Handling the input data by Fourier feature mapping will effectively improve the capacity of deep neural networks to solve high-frequency problems. Finally, to validate the efficiency and robustness of the proposed FMPINN algorithm, we present several numerical examples of multi-scale problems in various dimensional Euclidean spaces. These examples cover both low-frequency and high-frequency oscillation cases, demonstrating the effectiveness of our approach. All code and data accompanying this manuscript will be made publicly available at [https://github.com/Blue-Giant/FMPINN](https://github.com/Blue-Giant/FMPINN).

keywords: Multi-scale; Rough coefficient; FMPINN; Fourier feature mapping; Flux variable; Reduce order
**AMS subject classifications.** 35J25, 65N99, 68T07 +
Footnote †: journal: Journal of LaTeX Templates

## 1 Introduction

Multi-scale problems, governed by partial differential equations(PDEs) with multiple scales, are prevalent in diverse scientific and engineering fields like reservoir simulation, high-frequency scattering, and turbulence modeling. This paper focuses on solving the following type of multi-scale problem.

\[\begin{cases}-\mathbf{div}\bigg{(}A^{\varepsilon}(\mathbf{x})\nabla u^{\varepsilon }(\mathbf{x})\bigg{)}=f(\mathbf{x}),&\mathbf{x}\in\Omega,\\ \mathcal{B}u^{\varepsilon}(\mathbf{x})=g(\mathbf{x}),&\mathbf{x}\in\partial\Omega.\end{cases} \tag{1.1}\]

where \(\Omega\) is a bounded subset of \(\mathbb{R}^{d}(d=1,2,3,\ldots)\) with piecewise Lipschitz boundary and satisfies the interior cone condition, \(\varepsilon\) is a small positive parameter that signifies explicitly the multiscale nature of the rough coefficient \(A^{\varepsilon}(\mathbf{x})\). \(\mathcal{B}\) is a boundary operator in \(\partial\Omega\) that imposes the boundary condition of \(u^{\varepsilon}\), such as Dirchlete, Neumman and Robin. \(\nabla\) and \(\mathbf{div}\) are the gradient and divergence operators, respectively. \(f(\mathbf{x})\in L^{2}(\Omega)\) is a given function. In addition, \(A^{\varepsilon}(\mathbf{x})\) is symmetric and uniformly elliptic on \(\Omega\). It means that all eigenvalues of \(A^{\varepsilon}\) are uniformly bounded by two strictly positive constants
\(\lambda_{\min}(A^{\varepsilon})\) and \(\lambda_{\max}(A^{\varepsilon})\). In other word, for all \(\mathbf{x}\in\Omega\) and \(\mathbf{\xi}\in\mathbb{R}^{d}\), we have

\[\lambda_{\min}(A^{\varepsilon})|\mathbf{\xi}|^{2}\leqslant\mathbf{\xi}^{T}A^{\varepsilon }(\mathbf{x})\mathbf{\xi}\leqslant\lambda_{\max}(A^{\varepsilon})|\mathbf{\xi}|^{2}. \tag{1.2}\]

The multi-scale problem(1.1) frequently arise in the fields of physical simulations and engineering applications, including the study of flow in porous media and the analysis of mechanical properties in composite materials [1, 2, 3]. Generally, the analytical solutions of (1.1) are seldom available, then solving numerically this problem through approximation methods is necessary. Lots of numerical methods focus on efficient, accurate and stable numerical schemes have gained favorable achievement, such as heterogeneous multi-scale methods [2, 3, 4], numerical homogenization [5, 6, 7], variational multi-scale methods [8, 9], multi-scale finite element methods [10, 11, 12], flux norm homogenization [13, 14], rough polyharmonic splines (RPS)[15], generalized multi-scale finite element methods [16, 17, 18], localized orthogonal decomposition [19, 20], etc. In contrast to standard numerical methods including FEM and FDM, they alleviate substantially the computational complexity in handling all relevant scales, improve the numerical stabilities and expedite the convergence. However, they still will encounter the curse of complex domain and dimensionality in general.

Deep neural networks(DNN), an efficient meshfree method without the discretization for a given interested domain, have drawn more and more attention from researchers to solve numerically the ordinary and partial differential equations as well as the inverse problems for complex geometrical domain and high-dimensional cases [21, 22, 23, 24, 25, 26, 27], due to their extraordinary universal approximation capacity[28]. Among these methods, the physics-informed neural networks (PINN) dating back to the early 1990s again attracted widespread attention of researchers and have made remarkable achievements for approximating the solution of PDEs by embracing the physical laws with neural networks, on account of the rapid development of computer science and technology[24, 29]. This method skillfully incorporates the residual of governing equations and the discrepancy of boundary/initial constraints, then formulates a cost function can be optimized easily via the automatic differentiation in DNN. Many efforts have been made to further enhance the performance of PINN are concluded as two aspects: refining the selection of the residual term and designing the manner of initial/boundary constraints. In terms of the residual term, there are XPINN [30], cPINN [31], two-stage PINN [32] and gPINN [33], and so on. By subtly encoding the I/B constraints into DNN in a hard manner, the PINN can be easy to train with low computational complexity and obtain a high-precision solution of PDEs with complex boundary conditions[34, 35, 36]. Motivated by the reduction of order in conventional methods[12], some attempts have been made to solve the high-order PDEs by reframing them as some first-order systems, this will overcome the shortcomings of the computational burden for high-order derivatives in DNN. For example, the deep mixed residual method [27], the local deep learning method[37] and the deep FOSLS method[38, 39].

Many studies and experiments have indicated that the general DNN-based algorithms are commonly used to solve a low-frequency problem in varying dimensional space, but will encounter tremendous challenge for high-frequency problems such as multi-scale PDEs(1.1). The frequency principle[40] or spectral bias[41] of DNN shows that neural networks are typically efficient for fitting objective functions with low-frequency modes but inefficient for high-frequency functions. Then, a series of multi-scale DNN(MscaleDNN) algorithms were proposed to overcome the shortcomings of normal DNN for high-frequency problems by converting high-frequency contents into low-frequency ones via a radial scale technique [42, 43, 44, 45]. After that, some corresponding mechanisms were developed to explain this performance of DNN, such as the Neural Tangent Kernel (NTK)[46, 47]. Furthermore, many researchers attempted to utilize a Fourier feature mapping consisting of sine and cosine to improve the capacity of MscaleDNN, which will alleviate the pathology of spectral bias and let neural networks capture high frequencies component effectively[48, 49, 50, 45, 46, 45, 46].

Recently, some works[52, 53] have shown that general PINN architecture is unable to capture the multi-scale property of the solution due to the effect of rough coefficient in multi-scale PDEs. In [52], Wing Tat Leung et.al proposed a Neural homogenization-based PINN(NH-PINN) method to solve (1.1), it can well overcome the unconvergence of PINN for multi-scale problems. However, NH-PINN also will encounter the dilemma of dimensional and the burden of computation, because it will convert one low-dimensional problem into a high-dimensional case. By carefully analyzing the Neural Tangent Kernel matrix associated with the PINN, Sean P. Carney et. al[53] found that the Forbehuis norm of the NTK matrix will become unbound as the oscillation factor \(\varepsilon\) in \(A^{\varepsilon}\) tends to zero. It means that the evolution of residual loss term in PINN will become increasingly stiff as \(\varepsilon\to 0\), then lead to poor training behavior for PINN.
In this paper, a Fourier-based multi-scale mixed PINN(FMPINN) structure is proposed to solve the multi-scale problems (1.1) with rough coefficients. This method consists of the general PINN architecture and the aforementioned MscaleDNN model with subnetworks being used to capture different frequencies component. To overcome the weakness of the normal PINN that failed to capture the jumping gradient information of the oscillating coefficient when tackling the governed equation in multi-scale PDEs(1.1), a (dual)flux variable is introduced to alleviate the adverse effect of the rough coefficient. Meantime, it can also reduce the computational burden of PINN for the second-order derivatives of space variables. In addition, the Fourier feature mapping is used in our model to learn each target frequency efficiently and express the derivatives of multi-frequency functions easily, it will remarkably improve the capacity for our FMPINN model to solve multi-scale problems. In a nutshell, the primary contributions of this paper are summarized as follows:

1. We propose a novel neural networks approach by combining normal PINN and MscaleDNN with subnetworks structure to address multi-scale problems, leveraging the Fourier theorem and the F-principle of DNN.
2. Inspired by the reduced order scheme for high-order PDEs, a dual (flux) variable about the rough coefficient of multi-scale PDEs is introduced to address the gradient leakage about the rough coefficient for PINN.
3. By introducing some numerical experiments, we show that the classical PINN method with MscaleDNN solver is still insufficient in providing accurate solutions for multi-scale equations.
4. We showcase the exceptional performance of FMPINN in solving a class of multi-scale elliptic PDEs with essential boundaries in various dimensional spaces. Our method outperforms existing approaches and demonstrates its superiority in addressing these complex problems.

The remaining parts of our work are organized as follows. In Section 2, we briefly introduce the underlying conceptions and formulations for MscaleDNN and the structure of PINN. Section 3 provides a unified architecture of the FMPINN to solve the elliptic multi-scale problem (1.1) based on its equivalent reduced order scheme, and gives the option of activation function as well as the error analysis of our proppsed method. Section 4 details the FMPINN algorithm for approximating the solution of multi-scale PDEs, then provide the option of activation function and the simple error analysis for FMPINN method. In Section 5, some scenarios of multi-scale PDEs are performed to evaluate the feasibility and effectiveness of our proposed method. Finally, some conclusions of this paper are made in Section 6.

## 2 Multi-scale Physics Informed Neural Networks

### Multi-scale Deep Neural Networks with ResNet technique

The basic concept and formulation of DNN are described briefly in this section, which helps audiences to understand the DNN structure through functional terminology. Mathematically, a deep neural network defines the following mapping

\[\mathcal{F}:\mathbf{x}\in\mathbb{R}^{d}\Longrightarrow\mathbf{y}=\mathcal{F}(x)\in \mathbb{R}^{c} \tag{2.1}\]

with \(d\) and \(c\) being the dimensions of input and output, respectively. In fact, the DNN functional \(\mathcal{F}\) is a nested composition of the following single-layer neural unit:

\[\mathbf{y}=\{y_{1},y_{2},\cdots,y_{m}\}\ \ \text{and}\ \ y_{l}=\sigma\left(\sum_{n=1} ^{d}w_{ln}*x_{n}+b_{l}\right) \tag{2.2}\]

where \(w_{ln}\) and \(b_{l}\) are called weight and bias of \(l_{th}\) neuron, respectively. \(\sigma(\cdot)\) is an element-wise non-linear operator, generally referred as the activation function. Then, we have the following formulation of DNN:

\[\mathbf{y}^{[d]}=\sigma\circ(\mathbf{W}^{[\ell]}\mathbf{y}^{[\ell-1]}+\mathbf{b}^{[\ell]}),\ \ \text{for}\ \ \ell=1,2,3,\cdots\cdots,L \tag{2.3}\]

and \(\mathbf{y}^{[0]}=\mathbf{x}\), where \(\mathbf{W}^{[\ell]}\in\mathbb{R}^{n_{\ell+1}\times n_{\ell}},\mathbf{b}^{[\ell]}\in \mathbb{R}^{n_{\ell+1}}\) stand for the weight matrix and bias vector of \(\ell\)-th hidden layer, respectively, \(n_{0}=d\) and \(n_{L+1}\) is the dimension of output, and "\(\circ\) " stands for the elementary-wise operation. For convenience, the output of DNN is denoted by \(\mathbf{y}(\mathbf{x};\mathbf{\theta})\) with \(\mathbf{\theta}\) standing for its all weights and biases.
Residual neural network (ResNet) [54] as a common skillful technique by introducing skip connections between adjacent or nonadjacent hidden layers can overcome effectively the vanishing gradient of parameters in the backpropagation for DNN, then make the network much easier to train and improve well the performance of DNN. Many experiment results showed that the ResNet can also improve the performance of DNNs to approximate high-order derivatives and solutions of PDEs [21; 27]. We utilize the one-step skip connection scheme of ResNet in this work. Except for the normal data flow, the data will also flow along with the skip connection if the two consecutive layers in DNN have the same number of neurons, otherwise, the data flows directly from one to the next layer. The filtered \(\mathbf{y}^{[\ell+1]}(\mathbf{x};\mathbf{\theta})\) produced by the input \(\mathbf{y}^{[\ell]}(\mathbf{x};\mathbf{\theta})\) is expressed as

\[\mathbf{y}^{[\ell+1]}(\mathbf{x};\mathbf{\theta})=\mathbf{y}^{[\ell]}(\mathbf{x};\mathbf{\theta})+ \sigma\circ\bigg{(}\mathbf{W}^{[\ell+1]}\mathbf{y}^{[\ell]}(\mathbf{x};\mathbf{\theta})+\mathbf{b} ^{[\ell+1]}\bigg{)}.\]

As we are aware, a normal DNN model is capable of providing a satisfactory solution for general problems. However, it will encounter troublesome difficulty to solve multi-scale problems with high-frequency components. Recently, a MscaleDNN architecture has shown its remarkable performance to deal with high-frequency problems by converting original data to a low-frequency space [42; 43; 44; 46]. A schematic diagram of MscaleDNN with \(Q\) subnetworks is depicted in Fig. 1.

The detailed procedure of MscaleDNN is described in the following.

Figure 1: A schematic diagram of MscaleDNN with \(Q\) subnetworks, \(\sigma\) stands for the activation function
1. Generating a scale vector or matrix with \(Q\) parts \[\Lambda=(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3}\cdots,\mathbf{k}_{Q-1},\mathbf{k}_{Q})^{T},\] (2.4) where \(\mathbf{k}_{i}(i=1,2,\ldots,Q)\) is a scalar or matrix (trainable or untrainable).
2. Converting the input data \(\mathbf{x}\) into \(\tilde{\mathbf{x}}=\Lambda\odot\mathbf{x}\) with \(\odot\) being the Hadamard product, then feeding \(\tilde{\mathbf{x}}\) into the pipeline of MscaleDNN. It is \[\begin{cases}\hat{\mathbf{x}}=\mathbf{k}_{i}\mathbf{x}\\ \mathbf{F}_{i}(\mathbf{x})=\mathcal{FCN}_{i}\left(\hat{\mathbf{x}}\right)\end{cases}\qquad i =1,2,\ldots,Q,\] (2.5) where \(\mathcal{FCN}_{i}\) stands for the \(i_{th}\) fully connected subnetwork and \(\mathbf{F}_{i}\) is its output.
3. Obtaining the result of MscaleDNN by aggregating linearly the output of all subnetworks, each scale input goes through a subnetwork. It is \[\mathbf{NN}(\mathbf{x})=\mathbf{W}_{O}\cdot[\mathbf{F}_{1}(\mathbf{x}),\mathbf{F}_{2}(\mathbf{x}),\cdots, \mathbf{F}_{Q}(\mathbf{x})]+\mathbf{b}_{O},\] (2.6) where \(\mathbf{W}_{O}\) and \(\mathbf{b}_{O}\) stand for the weights and biases of the last linear layer, respectively.

From the perspective of Fourier transformation and decomposition, the first layer of the MscaleDNN model will be treated as a series of basis in Fourier space and its output is the combination of basis functions [42; 44; 46].

### Overview of Physics-Informed Neural Networks

In the scope of PINN, a type of PDE governed by parameters as the toy to show its implementation, it is

\[\begin{split}\mathcal{N}_{\mathbf{\lambda}}[\hat{u}(\mathbf{x})]=\hat{f} (\mathbf{x}),\quad\mathbf{x}\in\Omega\\ \mathcal{B}\hat{u}\left(\mathbf{x}\right)=\hat{g}(\mathbf{x}),\qquad\mathbf{x} \in\partial\Omega\end{split} \tag{2.7}\]

in which \(\mathcal{N}_{\mathbf{\lambda}}\) stands for the linear or nonlinear differential operator with parameters \(\mathbf{\lambda}\), \(\mathcal{B}\) is the boundary operator, such as Dirichlet, Neumann, periodic boundary conditions, or a mixed form of them. \(\Omega\) and \(\partial\Omega\) respectively illustrate the zone of interest and its border. For approximating the solution of the multi-scale PDE, a multi-scale deep neural network is used. In classical PINN, the ideal parameters of the DNN can be obtained by minimizing the following composite loss function

\[Loss=Loss_{R}+\gamma Loss_{B} \tag{2.8}\]

with

\[\begin{split} Loss_{R}&=\frac{1}{N_{R}}\sum_{i=1}^ {N_{R}}\left\|\mathcal{N}_{\mathbf{\lambda}}[u_{NN}(\mathbf{x}_{I}^{i})]-\hat{f}(\mathbf{ x}_{I}^{i})\right\|^{2}\\ Loss_{B}&=\frac{1}{N_{B}}\sum_{j=1}^{N_{B}}\left\| \mathcal{B}u_{NN}\left(\mathbf{x}_{B}^{j}\right)-\hat{g}(\mathbf{x}_{B}^{j})\right\|^ {2}\end{split} \tag{2.9}\]

where \(\gamma>0\) is used to control the contribution for the corresponding loss term. \(Loss_{R}\) and \(Loss_{B}\) depict the residual of the governing equations and the loss on the boundary condition, respectively. If some additional observed data are available inside the interested domain, then a loss term indicating the mismatch between the predictions produced by DNN and the observations can be taken into account

\[Loss_{D}=\frac{1}{N_{D}}\sum_{i=1}^{N_{D}}\left\|u_{NN}(\mathbf{x}^{i})-u_{Data}^{i }\right\|^{2}. \tag{2.10}\]

## 3 Fourier-based mixed PINN to solve multi-scale problem

In this section, the unified architecture of FMPINN is proposed to overcome the adverse effect of derivative for rough coefficient \(A^{e}\) by embracing a multi-output neural network with an equivalent reduced-order formulation of the multi-scale problem (1.1).
### Failure of classical PINN

Despite the success of various PINN models in studying ordinary and partial differential equations, it has been observed in [52] that the classical PINN approach fails to provide accurate predictions for multi-scale PDEs(1.1). Furthermore, we find that a direct application of the PINN with multi-scale DNN framework on solving (1.1) still cannot provide a satisfactory solution, because of the ill-posed NTK matrix caused by rough coefficient \(A^{\varepsilon}\). For example, let us consider the following one-dimensional elliptic equation with a homogeneous Dirichlet boundary in \(\Omega=[0,1]\):

\[\begin{cases}-\dfrac{d}{dx}\bigg{(}A^{\varepsilon}(x)\dfrac{d}{dx}u^{ \varepsilon}(x)\bigg{)}=5\cos(\pi x)\\ \qquad\qquad\qquad u^{\varepsilon}(0)=u^{\varepsilon}(1)=0\end{cases},\]

in which \(A^{\varepsilon}(x)=\dfrac{1+x^{2}}{2+\sin(2\pi x/\varepsilon)}\) with \(\varepsilon>0\) being a small constant.

We employ the classical PINN method with the MscaleDNN framework(see Fig. 1) to solve (1.1), called this method as MPINN. The scale factors \(\Lambda\) for MscaleDNN is set as \((1,2,3,4,5,10,\cdots,90,95,100)\) and the size of each subnetwork is chosen as \((1,30,40,30,30,30,1)\). The activation function of the first hidden layer for all subnetworks is set as Fourier feature mapping(see Section 3.3) and the other activation functions(except for their output layer) are chosen as \(\frac{1}{2}\sin(x)+\frac{1}{2}\cos(x)\)[55], their output layers are all linear. For \(\varepsilon=\frac{1}{32},\frac{1}{64}\) and \(\frac{1}{128}\),We train the aforementioned MPINN model for \(50000\) epochs and conduct testing every \(1000\) epochs within the training cycle. The optimizer is set as Adam with

Figure 2: Left: the rough coefficient \(A^{\varepsilon}\). Middle: the MPINN approximated solution vs the reference solution. Right: \(l^{2}\) relative error varies with the testing epoch.
an initial learning rate of \(0.01\) and the learning rate will decay by \(2.5\%\) for every \(100\) epochs. Finally, the results are demonstrated in Fig. 2.

As \(\varepsilon=1/32\), the coefficient \(A^{\varepsilon}(x)\) possesses a little multi-scale information, the MPINN performs quite well. However, the permeability \(A^{\varepsilon}(x)\) will exhibits various multi-scale properties for \(\varepsilon=1/64\), the performance of MPINN deteriorates with a low relative error and the MPINN fails to converge for \(\varepsilon=1/128\). In addition, we perform the MPINN with different setups of the hyperparameters such as the learning rate and the \(\omega_{B}\) for \(Loss_{B}\) in (2.8) as well as the network size, but we still cannot obtain a satisfactory result.

### Unified architecture of FMPINN

Based on the above observation, it is necessary to seek some extra techniques to improve the accuracy of the PINN. Inspired by the mixed finite element method [12; 56] and the mixed residual method[27], we can leverage a mixed scheme to solve (1.1) by replacing the flux term \(A^{\varepsilon}\nabla u\) in (1.1) with an auxiliary variable. This strategy not only can avoid the unfavorable effect of the oscillating coefficient \(A^{\varepsilon}\), but also can reduce the computation burden of second-order derivatives in cost function when utilizing a multi-scale deep neural network to approximate the solution of (1.1). Therefore, we introduce a flux variable \(\mathbf{\phi}(\mathbf{x})=\big{(}\phi_{1}(\mathbf{x}),\ldots,\phi_{d}(\mathbf{x})\big{)}=A^{ \varepsilon}(\mathbf{x})\nabla u^{\varepsilon}(\mathbf{x})\) and rewrite the first equation in (1.1) as the following expressions:

\[\begin{split}&-\text{\bf div}\mathbf{\phi}(\mathbf{x})=f(\mathbf{x})\\ &\mathbf{\phi}(\mathbf{x})-A^{\varepsilon}(\mathbf{x})\nabla u^{\varepsilon} (\mathbf{x})=\mathbf{0}\end{split} \tag{3.1}\]

Then we turn to search a couple of functions \((u^{\varepsilon},\mathbf{\phi})\) in admissible space, rather than approximating a unique solution of the original problem (1.1). Here and thereafter, \((u^{\varepsilon},\mathbf{\phi})\in\mathcal{A}=\mathcal{H}^{1}(\Omega)\times \mathcal{H}(\text{\bf div};\Omega)\) with \(\mathcal{H}^{1}(\Omega)=\big{\{}v\in L^{2}(\Omega):\nabla v\in L^{2}(\Omega) \big{\}}\) and \(\mathcal{H}(\text{\bf div};\Omega)=\big{\{}\mathbf{\psi}\in[L^{2}(\Omega)]^{d}: \text{\bf div}\mathbf{\psi}\in L^{2}(\Omega)\big{\}}\).

When utilizing numerical solvers to address the equation (3.1), one can obtain the optimum solution by minimizing the following least-squares formula in the domain \(\Omega\):

\[u^{*},\mathbf{\phi}^{*}=\operatorname*{arg\,min}_{(u,\mathbf{\phi})\in \mathcal{H}^{1}(\Omega)\times\mathcal{H}(\text{\bf div};\Omega)}\mathcal{L}(u,\mathbf{\phi}) \tag{3.2}\]

with

\[\mathcal{L}(u,\mathbf{\phi})=\int_{\Omega}\big{|}-\text{\bf div}\mathbf{\phi}(\mathbf{x}) -f(\mathbf{x})\big{|}^{2}d\mathbf{x}+\beta\int_{\Omega}\big{|}\mathbf{\phi}(\mathbf{x})-A^{ \varepsilon}(\mathbf{x})\nabla u^{\varepsilon}(\mathbf{x})\big{|}^{2}d\mathbf{x} \tag{3.3}\]

where \(\beta>0\) is used to adjust the approximation error of the flux variable and flux term.

Generally, two independent neural networks are necessary to approximate the flux variable \(\mathbf{\phi}\) and solution \(u\), but \(\mathbf{\phi}\) is unconstrained without any coercive boundary condition. Based on the potentiality of DNN for approximating any linear and non-linear complex functions, we take a DNN with multi outputs to model ansatzes \(\mathbf{\phi}\) and \(u\), denoted by \(\mathbf{\phi}_{NN}\) and \(u_{NN}\), respectively. Fig. 3 describes the multi-output neural network for input \(\mathbf{x}\in\mathbb{R}^{2}\).

Once the expressions of auxiliary functions \(\mathbf{\phi}\) and solution \(u\) have been determined, we can discretize

Figure 3: The multi-output neural network for approximating the state and flux variables
(3.3) by the Monte Carlo method [57], then employ the PINN conception and obtain the following form

\[\mathscr{L}_{in}(S_{I};\mathbf{\theta})=\frac{|\Omega|}{N_{in}}\sum_{i=1}^{N_{in}} \bigg{[}\big{|}-\mathbf{div}\mathbf{\phi}_{NN}(\mathbf{x}_{I}^{i};\mathbf{\theta})-f(\mathbf{x} _{I}^{i})\big{|}^{2}+\beta\big{|}\mathbf{\phi}_{NN}(\mathbf{x}_{I}^{i},\mathbf{\theta})-A^{ \varepsilon}(\mathbf{x}_{I}^{i})\nabla u_{NN}(\mathbf{x}_{I}^{i},\mathbf{\theta})\big{|}^ {2}\bigg{]}, \tag{3.4}\]

for \(\mathbf{x}_{I}^{i}\in S_{I}\), here and hereinafter \(S_{I}\) stands for the collection sampled from \(\Omega\) with prescribed probability density.

Same to the traditional numerical methods such as FDM and FEM for addressing PDEs, boundary conditions play a crucial role in DNN representation as well. They serve as important constraints that ensure the uniqueness and accuracy of the solution. Consequently, the output \(u_{NN}\) of DNN should also satisfy the boundary conditions of (1.1), which means

\[\mathscr{L}_{bd}(S_{B};\mathbf{\theta})=\frac{1}{N_{bd}}\sum_{j=1}^{N_{bd}}\bigg{[} \mathcal{B}u_{NN}\big{(}\mathbf{x}_{B}^{j};\mathbf{\theta}\big{)}-g(\mathbf{x}_{B}^{j}) \bigg{]}^{2}\to 0\ \ \text{for}\ \ \mathbf{x}_{B}^{j}\in S_{B}. \tag{3.5}\]

here and hereinafter \(S_{B}\) represents the collection sampled on \(\partial\Omega\) with prescribed probability density.

According to the above results, the weights and biases of the DNN model are updated by optimizing gradually the following cost function:

\[\mathscr{L}(S_{I},S_{B};\mathbf{\theta})=\mathscr{L}_{in}(S_{I};\mathbf{\theta})+ \gamma\mathscr{L}_{bd}(S_{B};\mathbf{\theta}) \tag{3.6}\]

where \(S_{I}=\{\mathbf{x}_{I}^{i}\}_{i=1}^{N_{in}}\) and \(S_{B}=\{\mathbf{x}_{B}^{j}\}_{j=1}^{N_{bd}}\) stand for the train data of \(\Omega\) and \(\partial\Omega\), respectively. The term of \(\mathscr{L}_{in}\) composed of the residual governed by differential equations and the discrepancy with respect to flux, minimizes the residual of the PDE, whereas the term of \(\mathscr{L}_{bd}\) pushes the DNN solver to match the given boundary conditions. In addition, a constant parameter \(\gamma>0\) is introduced to forces well the \(\mathscr{L}_{bd}(S_{B};\mathbf{\theta})\to 0\) in the loss function, it is increasing gradually with training process going on.

Based on the analysis in [39], a nonconstant continuous activation function \(\sigma\) can guarantee the mapping \(\mathbf{\theta}\mapsto(u_{NN},\mathbf{\phi}_{NN})\) is continuous, then the distance between approximation functions \(\mathbf{q}_{NN}=(u_{NN},\mathbf{\phi}_{NN})\) and exact solution \(\mathbf{q}^{*}=(u^{*},\mathbf{\phi}^{*})\) will decrease by adjusting gradually the parameters of DNN, i.e.,

\[d(\mathbf{q}^{*},\mathcal{A}_{k})=\inf_{\mathbf{q}_{NN}\in\mathcal{A}_{k}}\|\mathbf{q}^{*} -\mathbf{q}_{NN}\|\to 0\ \ \text{as}\ \ k\rightarrow\infty.\]

It means the loss function \(\mathscr{L}(S_{I},S_{B};\mathbf{\theta})\) will attain the corresponding minimum when \(d\to 0\).

Hence, Our purpose is to find an optimal set of parameter \(\mathbf{\theta}^{*}\) such that the approximations \(u_{NN}\) and \(\mathbf{\phi}_{NN}\) minimize the loss function \(\mathscr{L}(S_{I},S_{B};\mathbf{\theta})\). In order to obtain the ideal \(\mathbf{\theta}^{*}\), one can update the weights and biases of DNN through the optimization methods such as gradient descent (GD) or stochastic gradient descent (SGD) during the training process. In this context, the SGD method with a "mini-batch" of training data is given by:

\[\mathbf{\theta}^{k+1}=\mathbf{\theta}^{k}-\alpha^{k}\nabla_{\mathbf{\theta}^{k}}\mathscr{ L}(\mathbf{x};\mathbf{\theta}^{k})\ \ \text{with}\ \ \mathbf{x}\in S_{I}\ \text{or}\ \mathbf{x}\in S_{B} \tag{3.7}\]

where the "learning rate" \(\alpha^{k}\) decreases with \(k\) increasing.

### Option of activation function for FMPINN and its explanation

Choosing a suitable and effective activation function is a critical concern when aiming to enhance the performance of DNN in computer vision, natural language processing, and scientific computation. Generally, an activation function such as rectified linear unit \(\text{ReLU}(\mathbf{z})\) and hyperbolic tangent function \(\tanh(\mathbf{z})\), can obviously improve the capacity and nonlinearity of neural networkS to address various nonlinear problems, such as the solution of various PDEs and classification. Recently, the works [40; 41] manifested that the DNN often captures firstly the low-frequency component for target functions, then match the high-frequency component, they called it as the spectral bias or frequency preference of DNN. Under this phenomenon, many researchers attempt to utilize a Fourier feature mapping consisting of sine and cosine as the activation function to improve the capacity of MscaleDNN, it will mitigate the pathology of spectral bias and enable networks to learn high frequencies more effectively[41; 46; 49; 50]. It is expressed as follows:

\[\zeta(\mathbf{x})=\left[\begin{array}{c}\cos(\mathbf{\kappa}\mathbf{x})\\ \sin(\mathbf{\kappa}\mathbf{x})\end{array}\right], \tag{3.8}\]
where \(\mathbf{\kappa}\) is a user-specified vector or matrix (trainable or untrainable) which is consistent with the number of neural units in the first hidden layer for DNNs. Further, the work [45] designed a soften Fourier mapping by introducing a relaxing parameter \(s\in(0,1]\) in \(\zeta(\mathbf{x})\), numerical results show that this modification will improve the performance of \(\zeta(\mathbf{x})\). Actually, this activation function is used in the first hidden layer of DNN and maps the input data in \(\Omega\) into a range of \([-1,1]\), then enhances the ability of DNN and expedites its convergence.

Therefore, a real function \(\mathcal{P}(x)\) represented by DNN can be expressed as follows

\[\mathcal{P}(\mathbf{x})=\sum_{n=0}^{\bar{N}}\bigg{(}S\left(\cos(\mathbf{k}_{n}\mathbf{x}); \mathbf{\bar{\theta}}_{n}\right)+T\left(\sin(\mathbf{k}_{n}\mathbf{x});\mathbf{\bar{\theta}}_{n }\right)\bigg{)},\]

where \(S(\cdot,\mathbf{\bar{\theta}}),T(\cdot,\mathbf{\bar{\theta}})\) are the DNNs or the sub-modules of DNNs, respectively, \(\{\mathbf{k}_{0},\mathbf{k}_{1},\mathbf{k}_{2},\cdots\}\) are the frequencies of interest for the objective function. Obviously, the first hidden layer performed by Fourier feature mapping mimics the Fourier basis function, and the remaining blocks with different activation functions are used to learn the coefficients of these functions. After performing the Fourier mapping for input points with a given scale factor, the neural network can well capture the fine varying information for multi-scale problems.

**Remark 1**.: **(Lipschitz continuous)** If an activation function \(\sigma\) is continuous(i.e., \(\sigma\in C^{1}\)) and satisfies the following boundedness condition:

\[|\sigma(x)|<1\quad\text{ and }\quad|\sigma^{\prime}(x)|<1\]

for any \(x\in\mathbb{R}\). Then, we have

\[|\sigma(x)-\sigma(y)|\leqslant|x-y|\quad\text{ and }\quad|\sigma^{\prime}(x)- \sigma^{\prime}(y)|<|x-y|\]

for any \(x,y\in\mathbb{R}\). Obviously, the activation functions \(\tanh(x)\), sigmoid(x), Fourier feature mapping \(\zeta(a)\) and \(\frac{1}{2}sin(x)+\frac{1}{2}cos(x)\) are all satisfy the above condition and have a good regularity, they will overcome the gradient explosion of parameter in the backpropagation for DNN and improve the capacity of DNN.

### Simple error analysis for FMPINN

In recent times, there have been endeavors to rigorously analyze the convergence rate of the deep mixed residual method and compare it with the deep Galerkin method (DGM) and deep Ritz method (DRM) across different scenarios[39; 58; 59]. In this study, we investigate those results of convergence again, then provide the expression of generalization error for FMPINN and some remarks of errors.

For convenience, let \(\mathbf{q}^{*}=(u^{*},\mathbf{\phi}^{*})\) be the exact solution of equation (3.1) or the minimum of cost function (3.2) with (3.3) for coercive boundary constraints. Meantime, the \(\mathbf{q_{\theta^{*}}}=(u_{\mathbf{\theta^{*}}},\mathbf{\phi_{\theta^{*}}})\) stands for the final output of DNN optimized by SGD optimizer(such as Adam or LBFGS) that attains the local minimum of (3.6). Further, we let \(\widehat{\mathcal{L}}(u,\mathbf{\phi})\) be the cost function evaluated on \(N\) points sampled from \(\Omega\) and denote the output of DNN as \(\mathbf{q_{\theta}}=(u_{\mathbf{\theta}},\mathbf{\phi_{\theta}})\). Finally, \(\mathcal{S}_{NN}\) represents the function space sapnned by the output of DNN. Then, the total error(or generalization error) between the exact solution \(\mathbf{q^{*}}\) and the output of DNN \(\mathbf{q_{\theta}}\) can be expressed as

\[\big{\|}u_{\mathbf{\theta}}-u^{*}\big{\|}_{\mathcal{H}^{1}(\Omega)}+\big{\|}\mathbf{ \phi_{\theta}}-\mathbf{\phi^{*}}\big{\|}_{\mathcal{H}^{1}(\mathbf{div},\Omega)} \leqslant C(coe)\sqrt{\delta_{app}+\delta_{est}+\delta_{opt}} \tag{3.9}\]

Figure 4: Illustration of the total error for FMPINN.
with

\[\begin{cases}\delta_{app}=\inf_{(u,\mathbf{\phi})\in\mathcal{S}_{NN}}\|u-u^{*}\|^{2}_{ \mathcal{H}^{(1)}}(\Omega)+\|\mathbf{\phi}-\mathbf{\phi}^{*}\|^{2}_{\mathcal{H}(\mathbf{div},\Omega)}\\ \delta_{est}=\sup_{(u,\mathbf{\phi})\in\mathcal{S}_{NN}}[\mathcal{L}(u,\mathbf{\phi})- \widetilde{\mathscr{L}}(u,\mathbf{\phi})]+\sup_{(u,\mathbf{\phi})\in\mathcal{S}_{NN}}[ \widetilde{\mathscr{L}}(u,\mathbf{\phi})-\mathcal{L}(u,\mathbf{\phi})]\\ \delta_{opt}=\widetilde{\mathscr{L}}(u_{\mathbf{\theta}^{*}},\mathbf{\phi}_{\mathbf{\theta }^{*}})-\widetilde{\mathscr{L}}(u_{\mathbf{\theta}},\mathbf{\phi}_{\mathbf{\theta}})\end{cases}\]

In which, the approximated error \(\delta_{app}\) indicates the difference between \((u^{*},\mathbf{\phi}^{*})\) and its projection onto \(\mathcal{S}_{NN}\), the estimation error \(\delta_{est}\) measures the difference between the continuous cost function \(\mathcal{L}\) and discrete cost function \(\widetilde{\mathscr{L}}\), the optimization error \(\delta_{opt}\) stands for the discrepancy between the output of DNN with optimizing and the output of DNN without optimizing. In Fig. 4, we depict the diagram of error for FMPINN.

**Remark 2**: _For the approximating error, it is generally dependent on the architectural design of the neural network and the choice of the activation function. Classical radial basis network[60], the vanilla DNN and extreme learning machine(ELM)[61] are the common meshless method for approximating the solution of PDEs. To address the spatio-temporal problems, some hybrid network frameworks have been designed by combining PINN with traditional numerical methods to solve PDE, such as FDM-PINN and Runge-Kutta PINN[24; 62]. Moreover, instead of soft constraints by a hard manner for the boundary or initial conditions in those methods, the approximation will automatically meet the boundary and initial conditions of PDEs, then reduce the complexity and improve the precision of NN[35]. On the other hand, a powerful activation function, such as the hyperbolic tangent activation function and Fourier feature mapping, not only enhance the nonlinearity of DNN, but also improve its approximating capacity and accuracy. In addition, some available data are generally considered as a loss term to reduce the approximating error._

**Remark 3**: _Generally, the proposed FMPINN surrogate can provide more accurate approximations as the number of random collocation points increases. However, it will lead to heavy computational costs for lots of samplings. Then, it is worthwhile to take into account the trade-off between accuracy and computational cost when designing a DNN surrogate and determining its training mode. Alternatively, one can employ some effective low-discrepancy sampling approaches to decrease the statistical error, such as the Latin hypercube sampling method[63], quasi-random sampling[64] and multilevel Monte Carlo method[65]._

**Remark 4**: _Since the cost function generally is non-convex and has several local minima, then the gradient-based optimizer will almost certainly become caught in one of them. Therefore, choosing a good optimizer is important to reduce the optimization error and get a better minimum. In many scenarios of optimizing DNN, the Adam optimization method has shown its good performance including efficiency and accuracy, it can dynamically adjust the learning rates of each parameter by using the first and second moments estimation of the gradients[66]. BFGS is a quasi-Newton method and numerically stable, it may provide a higher-precision approximated solution[67]. In an implementation, the limited memory version of BFGS(L-BFGS) is the common choice to decrease the optimization error and accelerate convergence for cases with a little amount of training data and/or residual points. Further, by combining the merits of the above two approaches, one can optimize the cost function firstly by the Adam algorithm with a predefined stop criterion, then obtain a better result by the L-BFGS optimizer._

## 4 FMPINN algorithm

For the FMPINN method with the MscaleDNN model composed of \(Q\) subnetworks as in Fig.1 being its solver, the input data for each subnetwork will be transformed by the following operation

\[\widehat{\mathbf{x}}=a_{i}*\mathbf{x},\quad i=1,2,\ldots,Q\]

with \(a_{i}\geqslant 1\) being a positive scalar factor, it means the scale vector \(\Lambda=(a_{1},a_{2},\ldots,a_{Q})\) as in (2.4). Denoting the output of each subnetwork as \(\mathbf{F}_{i}(i=1,2,\ldots,Q)\), then the overall output of the MscaleDNN model is obtained by

\[\mathbf{y}(\mathbf{x};\mathbf{\theta})=\frac{1}{Q}\sum_{i=1}^{Q}\frac{\mathbf{F}_{i}}{a_{i}}.\]
According to the above discussions, the procedure of the FMPINN algorithm for addressing the multi-scale problem (1.1) in finite-dimensional spaces is described in the following.

1. Generating the \(k_{th}\) training set \(\mathcal{S}^{k}\) includes interior points \(S^{k}_{I}=\{\mathbf{x}^{i}_{I}\}_{i=1}^{N_{in}}\) with \(\mathbf{x}^{i}_{I}\in\mathbb{R}^{d}\) and boundary points \(S^{k}_{B}=\{\mathbf{x}^{j}_{B}\}_{j=1}^{N_{bd}}\) with \(\mathbf{x}^{j}_{B}\in\mathbb{R}^{d}\). Here, we draw the random points \(\mathbf{x}^{i}_{I}\) and \(\mathbf{x}^{j}_{B}\) from \(\mathbb{R}^{d}\) with positive probability density \(\nu\), such as uniform distribution.

2. Calculating the objective function \(\mathscr{L}(\mathcal{S}^{k};\mathbf{\theta}^{k})\) for train set \(\mathcal{S}^{k}\):

\[\mathscr{L}(\mathcal{S}^{k};\mathbf{\theta}^{k})=\mathscr{L}_{in}(S^{k}_{I};\mathbf{ \theta}^{k})+\gamma\mathscr{L}_{bd}(S^{k}_{B};\mathbf{\theta}^{k})\]

with \(\mathscr{L}_{in}(\cdot;\mathbf{\theta}^{k})\) being defined in (3.4) and \(\mathscr{L}_{bd}(\cdot;\mathbf{\theta}^{k})\) being defined in (3.5).

3. Take a descent step at the random point of \(\mathbf{x}^{k}\):

\[\mathbf{\theta}^{k+1}=\mathbf{\theta}^{k}-\alpha^{k}\nabla_{\mathbf{\theta}^{k}}\mathscr{ L}(\tilde{\mathbf{x}}^{k};\mathbf{\theta}^{k})\ \ \text{with}\ \ \tilde{\mathbf{x}}^{k}\in\mathcal{S}^{k},\]

where the "learning rate" \(\alpha^{k}\) decreases with \(k\) increasing.

4. Repeat steps 1-3 until the convergence criterion is satisfied or the objective function tends to be stable.

**Algorithm 1** FMPINN algorithm for solving multi-scale PDEs(1.1)

## 5 Numerical experiments

The goal of our experiments is to show that our Fourier-based mixed physics-informed neural networks are indeed capable of approximating the analytical solution given in (1.1). For comparison purposes, the PINN method with MscaleDNN being its solver and the local deep learning method(LDLM) with normal DNN being its solver are as the baseline to solve (1.1) in varying-dimensional spaces.

### Model and training setup

In the aforementioned FMPINN and MPINN models, a standard MscaleDNN with multi subnetworks that stretch the input data via various scale factors is configurated as their solver. The MscaleDNN consists of 25 subnetworks according to the manually defined frequencies vector \(\Lambda=(1,2,3,4,5,10,\cdots,90,95,100)\). Each subnetwork contains 5 hidden layers with proper size and the activation function of the first hidden layer for each subnetwork is set as Fourier feature mapping and the other activation functions(except for the output layer) are set as \(\frac{1}{2}\sin(x)+\frac{1}{2}\cos(x)\), its output layer is linear. The overall output is a weighted sum of the outputs of all subnetworks through the relevant scale factors. In terms of the LDLM[37], two activation functions are considered for this model: LDLM1 with \(ReQU=\max\{0,x\}^{2}\) being its activation for hidden layers and LDLM2 with \(\frac{1}{2}\sin(x)+\frac{1}{2}\cos(x)\) being its activation function for hidden layers, their output are all linear.

In our numerical experiments, all training data are sampled from the domain(including its boundaries) of interest in Euclidean space \(\mathbb{R}^{d}\), the sampling probability densities are assigned as the uniform distribution. We train all neural networks by an Adam optimizer with an initial learning rate of 0.01, and the learning rate will be decayed by 2.5% for every 100 training epochs [66]. Here, the following \(l^{2}\) relative error is used to evaluate our models:

\[REL=\sqrt{\frac{\sum_{i=1}^{N^{\prime}}|\tilde{u}(\mathbf{x}^{i})-u^{*}(\mathbf{x}^{i} )|^{2}}{\sum_{i=1}^{N^{\prime}}|u^{*}(\mathbf{x}^{i})|^{2}}}\]

where \(\tilde{u}(\mathbf{x}^{i})\) and \(u^{*}(\mathbf{x}^{i})\) are the approximate solution of deep neural network and exact solution for testing points \(\{\mathbf{x}^{i}\}(i=1,2,\cdots,N^{\prime})\), respectively, and \(N^{\prime}\) represents the number of sample points for testing. In order to visualize the training process, our model will be evaluated once for every 1000 iterations in the whole training cycle and recorded the result at the end. In our codes, the penalty
parameter \(\gamma\) is set as

\[\gamma=\begin{cases}\gamma_{0},&\text{if}\ \ i_{\text{epoch}}<M_{\text{max}}*0.1\\ 10\gamma_{0},&\text{if}\ \ M_{\text{max}}*0.1<=i_{\text{epoch}}<M_{\text{max}}*0.2\\ 50\gamma_{0},&\text{if}\ \ M_{\text{max}}*0.2<=i_{\text{epoch}}<M_{\text{max}}*0.25\\ 100\gamma_{0},&\text{if}\ \ M_{\text{max}}*0.25<=i_{\text{epoch}}<M_{\text{max}}*0.5\\ 200\gamma_{0},&\text{if}\ \ M_{\text{max}}*0.5<=i_{\text{epoch}}<M_{\text{max}}*0.75\\ 500\gamma_{0},&\text{otherwise}\end{cases} \tag{5.1}\]

where the \(\gamma_{0}=10\) in all our tests and \(M_{\text{max}}\) represents the total number of epochs. We implement and perform all neural network models by means of the package of Pytorch (version 1.14.0) on a workstation (64-GB RAM, single NVIDIA GeForce RTX 4090 24-GB).

### Performance of FMPINN for solving multi-scale elliptic PDEs

**Example 5.1**.: Firstly, we consider the one-dimensional case for (1.1) with Dirichlet boundary in interval \([0,1]\), in which \(A^{\varepsilon}(x)\) is given by

\[A^{\varepsilon}(x)=\left(2+\cos\left(2\pi\frac{x}{\varepsilon}\right)\right)^ {-1} \tag{5.2}\]

with a small parameter \(\varepsilon>0\) such that \(\varepsilon^{-1}\in\mathbb{N}^{+}\) and the force term \(f(x)=1\). Under these conditions, a unique solution is given by

\[u^{\varepsilon}(x)=x-x^{2}+\varepsilon\left(\frac{1}{4\pi}\sin\left(2\pi \frac{x}{\varepsilon}\right)-\frac{1}{2\pi}x\sin\left(2\pi\frac{x}{\varepsilon }\right)-\frac{\varepsilon}{4\pi^{2}}\cos\left(2\pi\frac{x}{\varepsilon} \right)+\frac{\varepsilon}{4\pi^{2}}\right). \tag{5.3}\]

Clearly, the analytical solution induces its boundary condition \(u(0)=u(1)=0\).

In this example, we use the FMPINN, MPINN, LDLM1, and LDLM2 models to solve (1.1) when \(\varepsilon=0.1,0.01\) and \(0.001\), respectively. The size of hidden layer for each subnetwork of FMPINN and MPINN is set as \((30,40,30,30,30)\) and the balance parameter \(\beta\) in (3.4) is set as \(10\). The hidden layer's size for LDLM is set as \((300,400,300,300,300)\). Their parameters' numbers are comparable. At each training step, we randomly sample \(3000\) points inside the \([0,1]\) and \(500\) boundary points as a training dataset. In addition, the testing dataset includes \(1000\) equidistant samples from \([0,1]\). All models are trained for \(50000\) epochs. We depict the related experiment results in Figs. 5, 6 and 7, respectively. Meantime, the final relative errors and total running time are listed in Table 1.

Based on these figures, the FMPINN model can perfectly capture the oscillation of the exact solution for \(\varepsilon=0.1,0.01\) and \(0.001\), but LDLM models are not convergent for these cases. At the same time, the performance of MPINN competes with that of FMPINN when \(\varepsilon=0.1\). However, the MPINN model fails to solve the multi-scale problem for \(\varepsilon=0.1\) and \(0.01\). Compared to \(\varepsilon=0.01\), the rough coefficient \(A^{\varepsilon}\) with \(\varepsilon=0.001\) have more oscillation in the interval \([0,1]\), but the FMPINN still can keep its remarkable performance. According to the point-wise errors in Figs. 5(d), 6(d) and 7(d) and the relative error in Figs. 5(h), 6(h) and 7(h), we can conclude that the FMPINN is able to approximate high-precisely the exact solution of (1.1) in one-dimensional space. In addition, the total time in Table 1 shows the running time of FMPINN is less than that of MPINN for \(50000\) training epochs.

_Influence of hyper-parameter \(\beta\):_ In the previous tests, the parameter \(\beta\) was initially set to \(10\). Now, we study the influence of \(\beta\) for our FMPINN model. In these tests, we take \(\varepsilon=0.001\) in (5.2), and consider values of \(\beta\) equal to \(1,5,10,15,20\) and \(25\), while keeping all other parameters fixed. All models with different \(\beta\) values are trained for \(50000\) epochs. Fig. 8 plots the results of flux loss for the training

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|} \hline  & \multicolumn{4}{c|}{REL} & \multicolumn{4}{c|}{Total time(s)} \\ \hline \(\varepsilon\) & FMPINN & MPINN & LDLM1 & LDLM2 & FMPINN & MPINN & LDLM1 & LDLM2 \\ \hline \(0.1\) & 2.92e-6 & 2.60e-7 & 0.3227 & 0.3389 & 680.734 & 865.849 & 345.791 & 373.537 \\ \hline \(0.01\) & 3.43e-5 & 0.94 & 0.3397 & 0.3406 & 689.729 & 868.199 & 351.451 & 377.089 \\ \hline \(0.001\) & 9.28e-5 & 0.99 & 0.3389 & 0.3398 & 691.458 & 875.297 & 358.435 & 388.273 \\ \hline \end{tabular}
\end{table}
Table 1: The relative error and running time of FMPINN, MPINN, LDLM1 and LDLM2 for Example 5.1
Figure 5: Rough coefficient, loss of flux term and testing results for Example 5.1 when \(\varepsilon=0.1\).
Figure 6: Rough coefficient, loss of flux term and testing results for Example 5.1 when \(\varepsilon=0.01\).
Figure 7: Rough coefficient, loss of flux term and testing results for Example 5.1 when \(\varepsilon=0.001\).
process as well as the relative error for testing. Additionally, the final relative errors obtained from the tests are listed in Table 2.

According to the above results in Fig. 8 and Table 2, it can be observed that the FMPINN model exhibits remarkable and stable performance across different values of \(\beta\). The performances of the FMPINN model for \(\beta=1\) and \(\beta=5\) are slightly weaker than that of other cases. The loss of flux term is also stable and consistent with the trendlines of REL. Therefore, for the subsequent tests, we will continue to set \(\beta=10\).

**Example 5.2**.: Let us attempt to solve the following three-scale problem with Dirichlet boundary in \(\Omega=[0,1]\). In which,

\[A^{\varepsilon}(x)=\left(2+\cos\left(2\pi\frac{x}{\varepsilon_{1}}\right) \right)\left(2+\cos\left(2\pi\frac{x}{\varepsilon_{2}}\right)\right)\]

with two small parameter \(\varepsilon_{1},\varepsilon_{2}>0\) such that \(\varepsilon_{1}^{-1},\varepsilon_{2}^{-1}\in\mathbb{N}^{+}\) and an exact solution is given by

\[u^{\varepsilon}(x)=x-x^{2}+\frac{\varepsilon_{1}}{4\pi}\sin\left(2\pi\frac{x} {\varepsilon_{1}}\right)+\frac{\varepsilon_{2}}{4\pi}\sin\left(2\pi\frac{x} {\varepsilon_{2}}\right). \tag{5.4}\]

Clearly, \(u^{\varepsilon}(0)=u^{\varepsilon}(1)=0\). One can obtain the force side after careful computation, we omit it here.

We solve the above three scale problem when \(\varepsilon_{1}=0.1\) and \(\varepsilon_{2}=0.01\) by employing the aforementioned FMPINN, MPINN, LDLM1 and LDLM2 models, respectively. Their all settings are same as the Example 5.1. The training dataset includes 3000 interior random points and 500 boundary random points, and the testing dataset includes 1000 equidistant samples. The related experiment results are listed in Table 3 and plotted in Fig. 9, respectively.

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|} \hline \(\beta\) & 1 & 5 & 10 & 15 & 20 & 25 \\ \hline REL & 3.31e-4 & 1.37e-4 & 9.27e-5 & 7.75e-5 & 6.40e-5 & 8.60e-5 \\ \hline \end{tabular}
\end{table}
Table 2: The relative error of FMPINN model VS various \(\beta\) for Example 5.1.

Figure 8: The loss of flux term VS training epoch and the relative error VS testing epoch for Example 5.1 when \(\varepsilon=0.001\).
Figure 9: Rough coefficient, loss of flux term and testing results for Example 5.2 when \(\varepsilon_{1}=0.1\) and \(\varepsilon_{2}=0.01\).
Fig. 9 shows that the FMPINN model still is well able to capture all oscillations of the exact solution for the three-scale problem, the MPINN model also captures the profile of the solution of (1.1) with \(\varepsilon_{1}=0.1\) and \(\varepsilon_{2}=0.01\). However, the LDLM1 and LDLM2 all fail to fit the solution. Figs. 9(d) - 9(g) not only show the point-wise errors of FMPINN for major points that are close to zero but also reveal the point-wise error of FMPINN is very smaller than that of the MPINN and the LDLM models are all bad. Additionally, Fig. 9(h) and Table 3 illustrate that the REL of FMPINN is superior to that of MPINN by more than two orders of magnitude, and its running time is 696.537 seconds and less than that of MPINN.

From the above results, we conclude that the FMPINN model is remarkable to address the (1.1) with rough coefficient in one-dimensional space, it generally outperforms the MPINN and LDLM models.

**Example 5.3**.: We consider the following two-dimensional problem for (1.1) with Dirichlet boundary in regular domains \(\Omega=[-1,1]\times[-1,1]\). In this example, we choose the \(f(x_{1},x_{2})=5\) and provide the following two-scales coefficient with scale separation

\[A^{\varepsilon}(x_{1},x_{2})=\frac{1.5+\sin(2\pi x_{1}/\varepsilon)}{1.5+\sin (2\pi x_{2}/\varepsilon)}+\frac{1.5+\sin(2\pi x_{2}/\varepsilon)}{1.5+\cos(2 \pi x_{1}/\varepsilon)}+\sin(4x_{1}^{2}x_{2}^{2})+1. \tag{5.5}\]

where \(\varepsilon>0\) is a small parameter such that \(\varepsilon^{-1}\in\mathbb{N}^{+}\). Since the corresponding exact solution can not be expressed explicitly in this example, then a reference solution \(u^{\varepsilon}(x_{1},x_{2})\) is set as the finite element solution computed by numerical homogenization method [15] on a square grid \([-1,1]\times[-1,1]\) with mesh-size \(h=1/128\).

We solve the above two scale problem when \(\varepsilon=0.05\) by employing the aforementioned FMPINN, MPINN, LDLM1 and LDLM2 models, respectively. The size of hidden layer for each subnetwork of FMPINN and MPINN is set as (40, 60, 40,40,40) and the hidden layers' size for LDLMs is set as \((400,250,250,200,200)\). At each training step, the training dataset includes 5000 points sampled inside the \(\Omega\) and 2000 boundary points sampled from the \(\partial\Omega\), respectively. In order to test our models, the testing dataset is the collection of all grid points in domain \([-1,1]\times[-1,1]\) with mesh-size \(h=1/128\). The related experiment results are listed in Table 4 and plotted in Fig. 10, respectively.

In this example, the \(A^{\varepsilon}(x_{1},x_{2})\) have two different frequency components and is quite oscillating(seeing Fig. 10(a)), then DNN will encounter some troubles to address multi-scale PDEs(1.1). According to the results of point-wise error (Figs.10(d) - 10(g)) and relative errors(Fig.10(h)), the performance of our FMPINN model is still superior to that of the MPINN, LDLM1 and LDLM2 models, and can obtain a favorable approximation to multi-scale problems(1.1). In addition, the test REL curve in Fig. 10(h) indicates the FMPINN model is stable in the whole training cycle and its tendency is consistent with the curve of loss for flux term in Fig. 10(c). Clearly, the running time of our FMPINN model is about half of that of the MPINN model, which means the FMPINN model is efficient in solving multi-scale PDEs(1.1) with two scales coefficient.

**Example 5.4**.: We consider the following two-dimensional problem for (1.1) with Dirichlet boundary in regular domains \(\Omega=[-1,1]\times[-1,1]\). In this example, we choose the \(f(x_{1},x_{2})=1\) and provide a multi-frequency coefficient

\[A^{\varepsilon}(x_{1},x_{2})=\Pi_{i=1}^{5}\bigg{(}1+0.5\cos\big{(}2^{i}\pi(x _{1}+x_{2})\big{)}\,\bigg{)}\bigg{(}1+0.5\sin\big{(}2^{i}\pi(x_{2}-3x_{1}) \big{)}\,\bigg{)}. \tag{5.6}\]

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline Method & FMPINN & MPINN & LDLM1 & LDLM2 \\ \hline REL & 0.0139 & 0.99 & 0.2431 & 0.2401 \\ \hline Total time(s) & 2098.258 & 3885.934 & 626.685 & 689.619 \\ \hline \end{tabular}
\end{table}
Table 4: The relative error and consumed time of FMPINN, MPINN, LDLM1 and LDLM2 for Example 5.3.

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline Method & FMPINN & MPINN & LDLM1 & LDLM2 \\ \hline REL & 3.36e-4 & 5.02e-2 & 0.5341 & 0.8372 \\ \hline Total time(s) & 696.537 & 965.076 & 377.18 & 399.575 \\ \hline \end{tabular}
\end{table}
Table 3: The relative error and consumed time of FMPINN, MPINN, LDLM1, and LDLM2 for Example 5.2.
Figure 10: Rough coefficient, reference solution, loss of flux term and testing results for Example 5.3
same as the Example 5.3, a reference solution \(u^{\varepsilon}(x_{1},x_{2})\) is set as the finite element solution computed by numerical homogenization method [15] on a square grid \([-1,1]\times[-1,1]\) with mesh-size \(h=1/128\).

Figure 11: Rough coefficient, reference solution, loss of flux term and testing results for Example 5.4
By meticulously implementing the previously mentioned FMPINN, MPINN, and LDLMs models with the specified setups, we obtain the approximated solution of (1.1) with (5.6). The setup for all models is identical to that of Example 5.3. During each training step, the training dataset comprises 5000 points randomly sampled from \(\Omega\) and 2000 boundary points sampled from the boundary \(\partial\Omega\), respectively. Meantime, the testing dataset composes of grid points on the square domain \([-1,1]\times[-1,1]\) with mesh-size \(h=1/128\). The related experiment results are listed in Table 5 and plotted in Fig. 11, respectively.

In this example, the \(A^{\varepsilon}(x_{1},x_{2})\) is obviously oscillating with six different frequency components(seeing Fig. 11(a)), it will increase the difficulty for DNN to address multi-scale PDEs(1.1). The point-wise error (Figs.11(d) - 11(g)) and the relative errors(Fig.11(h)) indicate that our FMPINN model is still favorable to capture the solution of multi-scale problems with complex multi-frequency coefficient, but the MPINN, LDLM1 and LDLM2 models all performance poorly for approximating the solution of (1.1). Additionally, the test REL curve in Fig. 11(h) and the curve of loss for flux term in Fig. 11(c) are all flat indicates the FMPINN model is stable in the whole training cycle. Moreover, the running time of our FMPINN model is less than that of the MPINN model in solving multi-scale PDEs(1.1) for coefficient (5.6).

**Example 5.5**.: We next study the performance of our FMPINN model to solve the elliptic equation (1.1) with Dirichlet boundary in a cubic domain \(\Omega=[0,1]\times[0,1]\times[0,1]\). In which, we take

\[A^{\varepsilon}(x_{1},x_{2},x_{3})=2+\sin\left(\frac{2\pi x_{1}}{\varepsilon} \right)\sin\left(\frac{2\pi x_{2}}{\varepsilon}\right)\sin\left(\frac{2\pi x _{3}}{\varepsilon}\right). \tag{5.7}\]

with a small parameter \(\varepsilon>0\) such that \(\varepsilon^{-1}\in\mathbb{N}^{+}\). Also, we let the force side \(f(x_{1},x_{2},x_{3})=20\) and the boundary function \(g(x_{1},x_{2},x_{3})=0\) on \(\partial\Omega\).

We utilize the FMPINN, MPINN, LDLM1 and LDLM2 models to approximate the solution of three-dimensional multi-scale problem (1.1) with rough coefficient (5.7) when \(\varepsilon=0.1\), the setups the four models are same as the Example 5.4. The training dataset includes 7500 interior points and 1000 boundary points randomly sampled from \(\Omega\) and \(\partial\Omega\), respectively. To facilitate the process, a reference solution \(u^{\varepsilon}(x_{1},x_{2},x_{3})\) is established as the numerical solution obtained using the finite difference method on the domain \([-1,1]\times[-1,1]\times[-1,1]\) with a mesh-size \(h=1/64\). The test dataset is formed by including all grid points within the domain \([-1,1]\times[-1,1]\) with a mesh-size \(h=1/64\), while keeping the value of \(z\) fixed at 0.3125. We list the total running time and REL in Table 6 and plot the related results in Fig. 12.

Based on the results in Fig.12, we can see that our FMPINN model still outperforms the MPINN and LDLMs model for multi-scale problems in three-dimensional space. The point-wise absolute error and the relative error of the former one are much smaller than that of the latter three, the precision of FMPINN is very good with a small absolute point-wise error. Additionally, the REL curve and the loss curve of the flux term are all flat in the later period of the training process, which means the performance of FMPINN is stable. The running time of FMPINN is 5179.601 seconds and less 3800 seconds than MPINN's.

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline Method & FMPINN & MPINN & LDLM1 & LDLM2 \\ \hline REL & 0.0071 & 0.0335 & 0.8048 & 0.5326 \\ \hline Total time(s) & 5179.601 & 9271.072 & 1065.541 & 1195.233 \\ \hline \end{tabular}
\end{table}
Table 6: The relative error and running time of FMPINN, MPINN, LDLM1 and LDLM2 for Example 5.5.

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline Method & FMPINN & MPINN & LDLM1 & LDLM2 \\ \hline REL & 0.0628 & 0.99 & 0.936 & 0.9127 \\ \hline Total time(s) & 2013.258 & 3985.934 & 606.685 & 659.619 \\ \hline \end{tabular}
\end{table}
Table 5: The relative error and consumed time of FMPINN, MPINN, LDLM1 and LDLM2 for Example 5.4.
Figure 12: Rough coefficient, exact solution, loss of flux term and testing results for Example 5.5
**Example 5.6**.: We consider the following eight-dimensional problem for (1.1) with Dirichlet boundary in regular domains \(\Omega=[0,1]^{8}\). In which, we take

\[A(x_{1},x_{2},\cdots,x_{8})=1+\frac{1}{8}\bigg{[} \cos(2\pi x_{1})+\cos(4\pi x_{2})+\cos(8\pi x_{3})+\cos(16\pi x_{4 })+\] \[\cos(16\pi x_{5})+\cos(8\pi x_{6})+\cos(4\pi x_{7})+\cos(2\pi x_{8 })\bigg{]}.\]

Meantime, an exact solution satisfied (1.1) is given by

\[u(x_{1},x_{2},\cdots,x_{8})=\prod_{j=1}^{8}\sin(\pi x_{j})\]

The functions \(f(x_{1},x_{2},\cdots,x_{8})\) in \(\Omega\) and \(g(x_{1},x_{2},\cdots,x_{8})\) on \(\partial\Omega\) are easy to obtain according to the rough coefficient and exact solution, we omit it.

In this example, we only perform the FMPINN, LDLM1 and LDLM2 model to solve the (1.1) in eight-dimensional space, because the huge computation requirement of MPINN has exceeded the limitation of memory for our station. The size of hidden layers for each subnetwork of FMPINN is set as (60, 80, 60, 60, 60) and the hidden layers' size for LDLM is set as \((400,500,300,300,300)\). At each training step, we construct the training dataset by sampling 20000 interior points inside the \(\Omega\) and 5000 boundary points from the \(\partial\Omega\). A testing dataset is given that included 1600 random points distributed in \(\Omega\). The related experiment results are plotted in Fig.13 and listed in Table 7. Additionally, the point-wise error for the FMPINN model evaluated on 1600 sample points is projected into a rectangular region with mesh size \(40\times 40\). Noting that the mapping is only aimed at visualizing, it is independent of the actual coordinates of those points.
For an eight-dimensional problem, the FMPINN still can obtain a satisfactory solution for (1.1) with small point-wise absolute error and relative error. However, the LDLM1 and LDLM2 both fail to approximate the solution of (1.1). Additionally, the loss of flux term and overall REL show that the FMPINN model is also stable during the training process. The running time of LDLMs is less thah that of FMPINN in Table 7, but their performance are obviously weaker that the latter's.

## 6 Conclusion

Physics-informed neural networks (PINN) have gained significant popularity in solving both forward and inverse problems. However, the normal PINN with a multi-scale DNN framework is unable to solve multiscale PDEs with rough coefficients. Inspired by the mixed finite element method, this work designs a Fourier-based mixed PINN(dubbed FMPINN) by combining a dual (flux) technique and Fourier decomposition to solve a class of elliptic multi-scale PDEs. By incorporating the loss of the flux term into the loss function, our model achieves improved stability and robustness. To handle multi-frequency contents, a Fourier activation function has been used to address the input data transformed radially by different frequency factors, and a sub-network is designed to match the target function, this strategy can improve clearly the accuracy and convergence rate for the FMPINN method. Compared to the previous works of PINN, this novel method skillfully casts the original problem into two first-order systems, it will overcome the shortcomings of the computational burden for high-order derivatives in DNN and the ill-condition of neural tangent kernel matrix resulting from the rough coefficient. Computational results show this novel method is feasible and efficient to solve this multi-scale equation with an inhomogeneous coefficient in various dimensional spaces. In the future, we aim to extend this novel network architecture, incorporating Fourier theory and lower-order mixed schemes, to tackle more complex multiscale problems.

## Declaration of interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Credit authorship contribution Statement

Xi'an Li: Conceptualization, Methodology, Investigation, Validation, Writing - Original Draft. Jinran Wu: Investigation, Writing - Review & Editing. You-Gan Wang: Writing - Review & Editing, Xin Tai: Writing - Review & Editing, Jianhua Xu: Writing - Review & Editing.

## Acknowledgements

The authors wish to thank Prof. Dr. Zhi-Qin John Xu for valuable suggestions which improved the quality of the paper.
