A hybrid FEM-PINN method for time-dependent partial differential equations
Xiaodong Feng Haojiong Shangguan Tao Tang Xiaoliang Wan Tao Zhou


	
	UIC]Xiaodong Feng
	xiaodongfeng@uic.edu.cn
	UIC]Haojiong Shangguan
	shangguanhaojiong@uic.edu.cn
	GZ,UIC]Tao Tang
	ttang@uic.edu.cn

	lu]Xiaoliang Wan
	xlwan@lsu.edu
	LSECL]Tao Zhou
	tzhou@lsec.cc.ac.cn
	[UIC]Division of Science and Technology, BNU-HKBU United International College, Zhuhai 519087, China
	

	[GZ]School of Electrical and Computer Engineering, Guangzhou Nanfang College, Guangzhou 510970, China
	

	[lu]
		Department of Mathematics and Center for Computation and Technology, Louisiana State University, Baton Rouge 70803, USA
	

	[LSECL]Institute of Computational Mathematics and Scientific/Engineering
		Computing, Academy of Mathematics and Systems Science, Chinese Academy
		of Sciences, Beijing, China
	
	§ ABSTRACT 
		In this work, we present a hybrid numerical method for solving evolution partial differential equations (PDEs) by merging the time finite element method with deep neural networks. In contrast to the conventional deep learning-based formulation where the neural network is defined on a spatiotemporal domain, our methodology utilizes finite element basis functions in the time direction where the space-dependent coefficients are defined as the output of a neural network.  We then apply the Galerkin or collocation projection in the time direction to obtain a system of PDEs for the space-dependent coefficients which is approximated in the framework of PINN. The advantages of such a hybrid formulation are twofold: statistical errors are avoided for the integral in the time direction, and the neural network's output can be regarded as a set of reduced spatial basis functions.
		To further alleviate the difficulties from high dimensionality and low regularity, we have developed an adaptive sampling strategy that refines the training set.
		More specifically, we use an explicit density model to approximate the distribution induced by the PDE residual and then augment the training set with new time-dependent random samples given by the learned density model.
		The effectiveness and efficiency of our proposed method have been demonstrated through a series of numerical experiments.
	
	
		Evolution equation Finite element methodDeep learning Adaptive sampling method
	



§ INTRODUCTION

Evolution equations, including both time-dependent ordinary and partial differential equations (ODEs / PDEs), are used
Many numerical approaches have been developed for such problems, e.g. the finite difference method, the spectral method, and the finite element method.
Recently solving PDEs with deep learning methods has been receiving increasing attention <cit.>. Typical techniques include physics-informed neural networks (PINNs) <cit.>, the deep Ritz methods <cit.>, the weak adversarial networks <cit.>, etc.
Although deep learning-based approaches have shown a lot of potential in solving high-dimensional PDEs, there still exist many numerical issues in adapting the neural network approximation to the problem studied.


 §.§ Related work

In this work, we pay particular attention to the error of PINNs for evolution equations which may grow too fast and limit the application of PINNs in long-term integration. Many efforts have been made to address this issue. We now briefly review the relevant works.

	Improved PINNs: PINNs represent the approximate PDE solution as a single neural network, which takes a space-time tuple as input and is trained by minimizing
the PDE residual on random collocation points in the space-time domain. To improve the performance of PINNs on long-term integration, many approaches have been developed which mainly focus on seeking a more effective training strategy. In <cit.>, a marching-in-time strategy is proposed by splitting the time domain into many small segments, where the training is done segment by segment and the approximate solution at the end of one segment is used as the initial condition for the next segment.
In <cit.>, backward-compatibility PINNs (bc-PINNs) are proposed, where the obtained solution in previous time segments is used as a constraint for the model training in the current time segment.
In <cit.>, Causal PINNs are developed to incorporate causality into the training process by introducing causality weights. In <cit.>, a unified scalable framework for causal sweeping strategies is developed.

	Evolution deep neural networks (EDNNs): EDNNs <cit.> are formulated with the Dirac-Frenkel variational principle to train networks by minimizing the residual sequentially in time, where the model parameters are time-dependent rather than global in  the whole space-time domain.
Traditional time-marching methods can be used to update the model parameters. Compared to PINNs, the training of EDNNs is more expensive while it is more flexible to adapt EDNNs to constraints such as Hamiltonian conservation <cit.>. The efficiency of EDNNs is improved in <cit.> by making a part of model parameters time-dependent and in <cit.> by updating randomized sparse subsets of model parameters at each time step.

	Operator learning: The main idea is to learn an operator that maps the solution from the current time step to the next time step.
For example, physics-informed DeepONet <cit.> can be used to learn a solution operator over a short time interval t∈[0,Δ t]. Starting with n=2, the model's prediction at nΔ t can be obtained from the trained model
using the approximate solution at (n-1)Δ t as the input. Other examples include auto-regressive networks <cit.>, transformer <cit.>, etc.

	Hybrid strategies: These approaches try to hybridize classical numerical methods with deep learning techniques by either adapting neural networks to augment classical PDE solvers <cit.> or adapting classical numerical approaches to improve the performance of PINNs <cit.>. For example, in <cit.>, a coupled automatic and numerical differentiation approach is proposed to take advantage of the regularization induced by numerical discretization.
In <cit.>, a deep adaptive basis Galerkin approach is proposed where the orthogonal polynomial expansion is employed in time direction and the expansion coefficients are modeled as the output of a deep neural network.



 §.§ Our contribution

The main contributions of this work are summarized as follows:

	
  * We have developed a hybrid numerical method by merging the time finite element method with deep neural networks. The approximate solution is a linear combination of the time finite element basis functions, where the coefficients are given by the output of a neural network. We subsequently apply Galerkin or collocation projection to eliminate the time and use PINN to approximate the system of PDEs for the coefficients. This strategy has some advantages: First, the numerical difficulties induced by random sampling and causality are avoided in the time direction since the projection can be done accurately.
	      Second, all the coefficients define a set of reduced basis functions on the computation domain, which are learned through the neural network. The approximate solution can also be regarded as a time-dependent linear combination of these reduced basis functions, which shares similarities with the low-rank representation in the study of high-dimensional problems.
	
  * We have proposed a deep adaptive sampling strategy to enhance the numerical efficiency. Mesh refinement in the time direction is straightforward. Particular attention needs to be paid to the random sampling in the physical space, especially for high-dimensional and low-regularity problems. Using a spatially conditional bounded KRnet and a discrete distribution in the time direction, we have constructed a joint density model to learn the distribution induced by the PDE residual, based on which new time-dependent samples are generated to refine the training set.

The remainder of the paper is organized as follows. In Section <ref>, we provide an overview of PINN. Then we briefly introduce the existing work in solving evolution PDEs and review classic finite element methods for the first-order ordinary differential equations, including Galerkin and collocation frameworks. In Section <ref>, we introduce a hybrid FEM-PINN method
and also develop some adaptive sampling strategies. Several numerical examples are demonstrated in Section <ref> to test the performance of the proposed method. The paper is concluded in Section <ref>.



§ PRELIMINARIES




 §.§ Physics-informed neural networks (PINNs)

We begin with a brief overview of physics-informed neural networks (PINNs). 
We consider a general time-dependent PDE 

    u_t(x,t) + 𝒩[u](x,t)=f(x,t),  x∈Ω⊂ℝ^d, t∈[0,T],

subject to the initial and boundary conditions

    u(x,0) =g(x),  x∈Ω,                                       
       ℬ[u](x,t) = b(x, t),  x∈∂Ω,  t∈[0,T],

where 𝒩[·] is a linear or nonlinear differential operator, and ℬ[·] is a boundary operator corresponding to Dirichlet, Neumann, Robin or periodic boundary conditions.

Following the original work of Raissi et al. <cit.>, we represent the unknown solution u(x, t) with a deep neural network u_θ(x, t), where θ denotes all tunable parameters (e.g. weights and biases). Then, a physics-informed model can be trained by minimizing the following composite loss function

    ℒ(θ) = λ_icℒ_ic(θ) + λ_bcℒ_bc(θ) + λ_rℒ_r(θ),

where

    ℒ_ic(θ) = 1/N_ic∑_i=1^N_ic|u_θ(x_ic^i, 0) - g(x_ic^i)|^2,                                                                          
       ℒ_bc(θ) = 1/N_bc∑_i=1^N_bc|ℬ[u_θ]( x_bc^i, t_bc^i) - b(x_bc^i, t_bc^i)|^2,                                       
       ℒ_r(θ)=1/N_r∑_i=1^N_r|∂ u_θ/∂ t(x_r^i, t_r^i) + 𝒩[u_θ](x_r^i, t_r^i) - f(x_r^i, t_r^i)|^2.

Here {x_ic^i}_i=1^N_ic, {t_bc^i,x_bc^i}_i=1^N_bc and {t_r^i, x_r^i}_i=1^N_r can be the vertices of a fixed mesh or points that are randomly sampled
at each iteration of a gradient descent algorithm. The gradients with respect to both the input variables (t,x) and the model parameters θ can be efficiently computed via automatic differentiation <cit.>. Moreover, the hyper-parameters {λ_ic, λ_bc, λ_r} allow the flexibility of assigning a varying learning rate to each loss term to balance their interplay during the training process, which may be user-specified or tuned automatically.




 §.§ Continuous time finite element method

Evolution PDEs are often approximated by spatial finite elements together with a sequential time-marching scheme. Another choice is to construct a finite element approximation space on the space-time domain. We briefly review the time finite element method for first-order ordinary differential equations. Consider the following model problem:

    u'(t)  + 𝒩[u](t) = f(t),    t∈ [0,T], 
        u(0) = 0,

where 𝒩 is a linear or nonlinear operator and f(t)∈ L_2(I) with I=[0,T].


  §.§.§ Galerkin projection

We let X {u∈ H^1(I):u(0)=0} be the trial space
and Y  L_2(I) the test space.
The primal variational formulation of (<ref>) is as follows.

    {     u ∈ X  
        (u', v) + (𝒩[u], v) = (f, v), ∀ v ∈ Y,
    	.

where (·,·) indicates the inner product of two functions. For approximation, we consider the Galerkin projection with X_N⊂ X and Y_N⊂ Y, i.e.,

    {    u_N ∈ Y_N
        (u_N' , v_N)  + (𝒩[u_N], v_N) = (f, v_N) ,  ∀ v_N ∈ Y_N.
    	.

Define

    X_N = span{ϕ_j(t)|0≤ j≤ N},    Y_N =  span{ψ_j(t)|0≤ j≤ N},

where ϕ_j(t) and ψ_j(t) are finite element basis functions. Let

    u_N(t) = ∑_i=0^Nũ_iϕ_i(t) ∈ X_N

be the approximate solution with undetermined coefficients ũ_i. Taking v_N=ψ_j for j=0,⋯,N in (<ref>) leads to the following system

    ∑_i=0^N(∂_t ϕ_i(t),ψ_j(t)) ũ_i+ ( N[∑_i=0^Nϕ_i(t)ũ_i],ψ_j(t) )=(f(t),ψ_j(t)),   ∀ j=0,1, ⋯, N.

The inner products in (<ref>) can be accurately evaluated with the Gaussian quadrature formulas, where the degree of exactness is determined by the nonlinearilty of 𝒩.


  §.§.§ Collocation projection

Collocation projection provides a flexible strategy especially when 𝒩 is nonlinear <cit.>.Let {s_k}_k=1^K be Gaussian-type quadrature points on the reference interval [0,1]

    0≤ s_1<s_2<⋯<s_K≤ 1.

Consider a partition of [0,T] with

    0=t_0<t_1<⋯ <t_M̂=T,    h_i = t_i - t_i-1,   i=1,⋯,M̂.

Define

    s_m,k = t_m-1 + h_ms_k,   1≤ k ≤ K,   1≤ m≤M̂.

We seek the approximate solution by enforcing the equation on the collocation points, i.e.,

    {    u∈ X_N∩ C^1(I)  
       ∂ _t u(s) + 𝒩[u](s) = f(s),    s∈∪_m=1^M̂{s_m,k}_k=1^K,
    	.

where X_N∩ C^1(I) defines a finite element approximation space with C^1 elements and N+1 is equal to the total number of collocation points. It is shown in <cit.> by selecting the collocation points carefully collocation projection yields the same order of accuracy as Galerkin projection.  Typical piecewise polynomials with at least C^1 regularity include piecewise cubic Hermite polynomials and cubic spline functions.




§ METHODOLOGY
 
Now we are ready to present our approach for evolution equations (<ref>).
We aim to seek an approximate solution of the following form

    u_N(x,t;θ) = ∑_i=0^Nω_i(x;θ)ϕ_i(t),

where {ϕ_i(t)}_i is a pre-specified set of time finite element basis functions, ω_i: ℝ^d→ℝ are modeled by the output of a neural network ω(x,θ):ℝ^d→ℝ^N+1, and θ includes all tunable model parameters.
More precisely, ω(x,θ) is a fully-connected neural network defined as

    ω(x;θ)  a^T h_L-1∘ h_L-2∘⋯∘ h_1(x)  for x∈ℝ^d,

where L∈ℕ^+, a∈ℝ^M_L-1× (N+1),   h_ℓ(x_ℓ)σ(W_ℓx_ℓ+ b_ℓ) with W_ℓ∈ℝ^M_ℓ× M_ℓ-1  (M_0 d) and b_ℓ∈ℝ^M_ℓ for ℓ=1,2,⋯,L-1. Then θ{a, W_ℓ, b_ℓ:1≤ℓ≤ L-1 }.
σ(x) is an activation function which acts on x componentwisely to return a vector of the same size as x. We let M_ℓ=M be a fixed number for all ℓ and ℱ_L,M the set consisting of all ω with depth L and width M.


 §.§ A hybrid FEM-PINN method

We consider the following hypothesis space

    𝒰_N{
    	u_N(x,t;θ) = ∑_i=0^Nω_i(x;θ)ϕ_i(t),  ω=(ω_0,⋯,ω_N)∈ℱ_L,M}.

The Galerkin projection along the time direction yields that

    {[                                                                                   Find  u_N∈𝒰_N; (∂_t u_N(x,·), v_N) + (𝒩[u_N](x,·), v_N) = (f(x,·), v_N)   ∀ v_N∈span{ψ_j(t)|0≤ j≤ N},   ∀ x∈Ω. ].


where (·,·) indicates the inner product of two functions with respect to time. More specifically, we have

    ∑_i=0^N(∂_t ϕ_i(t),ψ_j(t)) ω_i(x;θ) + ( N[∑_i=0^Nϕ_i(t)ω_i(x;θ)],ψ_j(t) )=(f(x,t),ψ_j(t)),   ∀ j=0,1, ⋯, N,  ∀x∈Ω.

If 𝒩 is linear with respect to ω and time-independent, the above system can be further simplified as follows:

    ∑_i=0^N(∂_t ϕ_i(t),ψ_j(t)) ω_i(x;θ)+
    	∑_i=0^N(ϕ_i(t),ψ_j(t))𝒩(ω_i(x;θ))
    	=(f(x, t),ψ_j(t)),   ∀ j=0,1, ⋯, N.


Now let us turn to the collocation projection along the time direction.
Let S_t=∪_m=1^M̂{s_m,k}_k=1^K, where s_m,k is defined in equation (<ref>). The collocation formulation can be written as

    {    u_N ∈𝒰_N 
       ∂_t u_N(x,s;θ) + 𝒩[u_N](x,s;θ) = f(x, s),∀ s∈ S_t, ∀x∈Ω.
    	.

More specifically,

    ∑_i=0^N∂_t ϕ_i(s) ω_i(x;θ) + 𝒩[∑_i=0^Nϕ_i(s)ω_i(x;θ)] = f(x,s),   ∀ s∈ S_t, ∀x∈Ω.



	The Galerkin and collocation projections yield respectively two systems of PDEs for ω_i(x;θ). Due to the hybrid form of u_N(x,t;θ), all integrals for the Galerkin projection in the time direction can be done accurately by Gaussian quadrature formulas. Since the temporal basis functions are polynomials, collocation projection is also effective <cit.>. We then mainly focus on the integration in the physical space.


We subsequently approximate the PDE systems (<ref>) and (<ref>) in the framework of PINNs.
More specifically, we consider the following minimization problem

    min_θℒ(θ) = ℒ_r(θ) + γ_1 ℒ_ic(θ) + γ_2 ℒ_bc(θ),

where

    ℒ_ic(θ) = ‖ u(x,0;θ) - g(x)‖ _L_2(Ω)^2,  ℒ_bc(θ) = ‖ℬ[u](x,t;θ) - b(x,t)‖ _L_2(∂Ω× [0,T])^2,

with 0<γ_1,γ_2<∞ being penalty parameters. For system (<ref>), we define

    ℒ_r(θ) = ∑_j=0^Nℒ^g_r,j(θ),                                                                                                                                                                                                                                                        
       ℒ^g_r,j(θ) = ‖∑_i=0^N(∂_t ϕ_i(t),ψ_j(t)) ω_i(x;θ)+ ( N[∑_i=0^Nϕ_i(t)ω_i(x;θ)],ψ_j(t) )-(f(x, t),ψ_j(t))‖^2_L_2(Ω)=r_j^g(x;θ)^2_L_2(Ω).

For system (<ref>), we define

    ℒ_r(θ) = ∑_j=1^|S_t|ℒ^c_r,j(θ),                                                                                                                                                                                                                
       ℒ^c_r,j(θ) = ‖∑_i=0^N∂_t ϕ_i(s_j)ω_i(x;θ) + 𝒩[∑_i=0^Nω_i(x;θ)ϕ_i(s_j)] - f(x, s_j)‖^2_L_2(Ω)=r_j^c(x;θ)^2_L_2(Ω),

where we order all collocation points in S_t as s_j with j=1,…,|S_t|.

We note that |S_t|≥ (N+1) in general since there are N+1 time finite element basis functions. For simplicity, we let |S_t|=N+1 and consider the following form

    ℒ_r(θ) = ∑_j=0^Nℒ_r,j(θ),  ℒ_r,j(θ) = r_j(x;θ)^2_L_2(Ω),

which are shared by both the Galerkin and collocation projections, i.e., r_j=r_j^g or r_j^c.
The loss functional (<ref>) is usually discretized numerically before the optimization with respect to θ is addressed. In practice, one often chooses uniformly distributed collocation points S_r = {S_r,j}_j=0^N={{x_r,j^(i)}_i=1^N_r,j}_j=0^N, S_ic={x_ic^(i)}_i=1^N_ic on Ω and S_bc = {(x_bc^(i), t_bc^(i))}_i=1^N_bc on ∂Ω×[0,T] for the discretization of the three terms in the objective functional (<ref>), leading to the following empirical loss

    ℒ(θ) = ∑_j=0^N‖ r_j(x;θ)‖^2_N_r,j + γ̂_1 ‖ℬ[u](x,t;θ) - b(x,t)‖ _N_ic ^2 + γ̂_2 ‖ u(x,0;θ) - g(x)‖^2_N_bc,

where 0<γ̂_1, γ̂_2<∞, and

    ‖ u(x) ‖_N_r,j = (1/N_r,j∑_i=1^N_r,ju^2(x_r,j^(i)))^1/2, ‖ u(x,0)‖ _N_ic = (1/N_ic∑_i=1^N_icu^2(x_ic^(i),0))^1/2, ‖ u(x,t)‖_N_bc = (1/N_bc∑_i=1^N_bcu^2(x_bc^(i), t_bc^(i)))^1/2.

We then seek an estimator θ̂ by minimizing the empirical loss (<ref>) via stochastic gradient descent methods, i.e.,

    θ̂=_θℒ(θ).

As suggested by <cit.>, we can also employ a time-marching strategy to reduce optimization difficulties. Specifically, we partition the temporal domain [0,T] into sub-domains [0,Δ t], [Δ t, 2Δ t],⋯, [T-Δ t, T].
Neural networks are then trained on each sub-domain successively with the initial conditions given by the same model trained on previous sub-domains.
The schematic of the proposed approach is shown in Figure <ref>, and the corresponding algorithm is summarized as follows.





  §.§.§ Some remarks on the hybrid form


PINN is formulated as a least-square method in terms of the hypothesis space given by the neural network. The error of u_N(x,t;θ̂) satisfies

    𝔼u_exact(x,t)-u_N(x,t;θ̂)≤u_exact(x,t)-u_N(x,t;θ^*)+𝔼u_N(x,t;θ^*)-u_N(x,t;θ̂)

for a proper norm, where θ^* is the minimizer of ℒ(θ), θ̂ is the minimizer of ℒ(θ) and the expectation 𝔼[·] is with respect to random samples. On the right-hand side, the first term is the approximation error determined by the hypothesis space and the second term is the statistical error introduced by the random samples.

It is well known that PINN may fail to predict convection when the frequency is large although the hypothesis space is capable of yielding a good approximate solution <cit.>. According to the inequality (<ref>), the reason is twofold: the non-convex structure of the loss landscape and the statistical error. First, the best approximation may not be obtained due to non-convexity of the loss function even the statistical error is zero. The change of the loss landscape can be achieved by adding a regularization term. For example, bc-PINNs have a penalty term to force the model to remember what was learned before <cit.>. Second, the available strategies that improve the performance of PINNs for time integration can be understood through the reduction of statistical error. Assume that N_t random samples are used in the time direction. The most straightforward strategy is to divide the time interval as [0,T]=∪_i=0^n-1[iΔ T,(i+1)Δ T] with Δ T=T/n and train the model sequentially in each time segment. After such a decomposition, the variation in time is reduced implying that the Monte Carlo approximation of the loss given by the random samples is more accurate due to variance reduction. The better the loss is discretized by random samples, the smaller the statistical error is. Another strategy, called causal training, is proposed in <cit.>. A weighted residual loss function is defined as

    ℒ_r(θ)=1/N_t∑_i=1^N_tλ_iℒ_r(t_i,θ),

where ℒ_r(t_i,θ) is the residual loss at t=t_i, and

    λ_i=exp(-ϵ∑_j=1^i-1ℒ_r(t_j,θ)),   i=2,3,…,N_t

with 0<ϵ<∞. The intuition is that the model will not be trained until the model is well trained for small t_i, which is consistent with the causality induced by evolution. Note that

    1/N_t∑_i=1^N_tλ_iℒ_r(t_i,θ)≈1/T∫_t λ(t)ℒ_r(t,θ)dt

is the Monte Carlo approximation of a weighted loss with N_t uniform random samples. If λ(t)>0 and the exact solution is included in the hypothesis space, the same θ^* will be reached.  λ(t) is a decreasing function by definition while ℒ_r(t,θ) is in general an increasing function due to the accumulation of errors with time. If λ(t) and ℒ_r(t,θ) are well balanced, their product varies much less in time, corresponding to a small variance in terms of the uniform distribution. Such a balance is mainly achieved by the selection of the so-called causality parameter ϵ. If ϵ fails to introduce a variance reduction for λ(t)ℒ_r(t,θ), the statistical error will not be reduced, implying that the training results may get worse. This explains the sensitivity of the training strategy on ϵ.

Based on the above observations, we intend to use the hybrid form  (<ref>) to alleviate the difficulties induced by the statistical errors in the time direction. We also note that the coefficients for the time finite element basis functions are given by the outputs of a neural network, which corresponds to learning a set of reduced basis functions in the physical space since the output layer of the neural network is a linear combination of these basis functions.




 §.§ Deep adaptive sampling method

Random samples are used for the integration in the physical space. To reduce the statistical errors, we consider the adaptive sampling method <cit.>. For simplicity, we only consider the interior residual ℒ_r(θ), and the time interval is [0,1]. As suggested in <cit.>, we relax the objective function ℒ_r(θ) as:

    ℒ_r(θ) = ∑_i=0^Nλ_iℒ_r,i(θ) = ∑_i=0^Nλ_i∫_Ω r_i^2(x;θ)p_i(x)dx≈1/N_r∑_i=0^N∑_j=1^N_rλ_ir_i^2(x^(i)_j;θ),

where λ_i>0, ∑_i=0^Nλ_i=1, the set {x^(i)_j}_j=1^N_r is generated with respect to the probability density function p_i(x)>0 instead of a uniform distribution. We associate ℒ_r,i(θ) with a weight λ_i. The minimizer of ℒ_r(θ) is also the solution to the problem if the exact solution is included in the hypothesis space. To reduce the error induced by the Monte Carlo approximation, we may adjust p_i(x) to make 
the residuals r_i^2(x;θ) nearly uniform.

To do this, we refine the training set gradually by adding new samples according to the distribution induced by r_i^2(x;θ^(k)), where k indicates the adaptivity iteration and θ^(k) is the optimal model parameter given by the previous training set. Once the training set is updated, the model will be retrained, based on which the training set will be refined again. In a nutshell, the model and the training set are updated alternately.

The main problem is that we need to obtain new samples from N+1 distributions induced by r_i^2(x;θ^(k)).

To handle this issue, we adopt a time-dependent density estimation strategy. Specifically, we augment the spatial variable in different terms (ℒ_r,i) with an extra dimension s,  and consider a set {s_i}_i=0^N of grid points on the time interval.
For the Galerkin approach, s_i can be interpreted as pre-defined nodes of finite element mesh; for the collocation approach, s_i can be viewed as a reordering of pre-defined Gaussian nodes s_m,k (see Section <ref>). We define
a weighted empirical measure

    δ_λ(A)=∑_i=0^Nλ_iδ_s_i(A)

for any A⊂[0,1] with δ_s_i being the Dirac measure
and let r(x,s;θ) be an interpolation function satisfying

    r(x,s;θ) = r_i(x;θ)   if  s = s_i,    i=0,⋯,N.

Let p_X,S(x,s)=p_X|S(x|s)p_S(s) be a joint PDF. Choosing p_S(s)ds=δ_λ(ds).

We have

    ∫∫ _Ω r^2(x,s;θ)p_X,S(x,s)dxds=∑_i=0^N∫_Ωr_i^2(x;θ)p_X|S(x|s_i)λ_idx,

which is consistent with equation (<ref>) if p_X|S(x|s_i)=p_i(x). Using p_X,S(x,s), the objective functional is discretized as

    ℒ_r(θ)≈1/N_r∑_i=1^N_rr^2(x^(i),s^(i);θ),

where {(x^(i),s^(i))}_i=1^N_r are sampled from p_X,S(x,s). We will use a density model with the form p_X|S(x|s)p_S(s) to approximate the distribution induced by r^2(x,s;θ). New samples from the trained density model will be added to the training set for refinement.





  §.§.§ Model p_S(s)

Without loss of generality, we assume that s∈[0,1]. We aim to find a invertible transformation z=f(s) such that

    p_S(s) = p_Z(f(s))|∇_s f|,    Z∼𝒰[0, 1],

where 𝒰 denotes the uniform distribution. We use the bounded polynomial spline layer f_poly <cit.> to parameterize f. Specifically, let 0=l_0<l_1<⋯<l_m-1<l_m=1 be a given partition of the unit interval and {k_j}_j=0^m be the corresponding weights satisfying ∑_j k_j=1. A piecewise linear polynomial can be defined as follows:

    p(s) =k_j+1-k_j/l_j+1-l_j(s-l_j) + k_j,  ∀ s∈[l_j,l_j+1].

Then the corresponding cumulative probability function f_poly admits the following formulation:

    f_poly(s) = k_j+1-k_j/2(l_j+1-l_j)(s-l_j)^2+k_j(s-l_j)+∑_i=0^j-1k_i+k_i+1/2(l_i+1-l_i),  ∀ s∈[l_j,l_j+1].

To satisfy ∫_0^1p(s) ds=1, we can model {k_j}_j=0^m as

    k_j = exp(k̃_j)/C,  ∀ j = 0, …, m,

where θ_f,1 = {k̃_j}_j=0^m are trainable parameters and C is a normalization constant:

    C = ∑_i=0^m-1(exp(k̃_i)+exp(k̃_i+1))(l_i+1-l_i)/2.


Notice that the polynomial spline layer f_poly(·;θ_f,1) (<ref>)-(<ref>) yields explicit monotonous expressions, and its inverse can be readily computed. Then an explicit PDF model p_poly(s;θ_f,1) can be obtained by letting f=f_poly, i.e.,

    p_poly(s;θ_f,1) = p_Z(f_poly(s))|∇ _sf_poly|.



  §.§.§ Model p_X|S(x|s)

For x∈ℝ^d, we seek a invertible transformation z=f(x, s)∈ℝ^d for any given s such that

    p_X|S(x|s) = p_Z|S(z|s)| ∂ f(x,s)/∂x|,   Z|S∼𝒰[-1, 1]^d,  ∀ s.

Here we employ conditional bounded KR-net f_B-KRnet(·, s) <cit.> to parameterize f(·, s). The basic idea of conditional bounded KRnet is to define the structure of f(x,s) in terms of the Knothe-Rosenblatt rearrangement. The transformation f(·, s) inspired by the Knothe-Rosenblatt (K-R) rearrangement <cit.> has a low-triangular structure

    z = f(x,s) = [[       f_1(x_1,s);   f_2(x_1,x_2,s);                ⋮; f_d(x_1,⋯,x_d,s) ]].

The sub-transformations f_1,⋯,f_d consist of polynomial spline layers and coupling layers <cit.>. More details can be found in <cit.>. Let f_B-KRnet(·,s;θ_f,2) indicate the conditional invertible transport map induced by  bounded KR-net, where θ_f,2 includes the model parameters. Then an explicit PDF model p_B-KRnet(x,s;θ_f,2) can be obtained by letting f=f_B-KRnet in equation (<ref>)

    p_B-KRnet(x|s;θ_f,2) = p_Z(f_B-KRnet(x,s)) |∇_x f_B-KRnet|.






  §.§.§ Adaptive sampling approach

Now we model a continuous joint density distribution p_θ_f(x,t)

    p_θ_f(x,t) = p_poly(t;θ_f,1) p_B-KRnet(x|t;θ_f,2),

where θ_f = {θ_f,1, θ_f,2}. To seek the "optimal" parameter θ_f, we can minimize the following objective

    D_KL(r̂_θ(x,t) || p_θ_f(x,t))
    		     = D_KL(r̂_θ(x,t)|| p_poly
    		(t;θ_f,1)p_B-KRnet(x|t;θ_f,2))                                              
        =∬r̂_θ(x,t) log(r̂_θ(x,t)) dxdt
    		- ∬r̂_θ(x,t) log( p_poly(t;θ_f,1)p_B-KRnet(x|t;θ_f,2)) dxdt,

where D_KL indicates the Kullback-Leibler (KL) divergence and r̂_θ(x,t)∝ r^2(x,t;θ) is the induced measure by continuous residual squared r^2(x,t;θ)

    r(x,t;θ) = ∂_t u_N(x,t;θ) - 𝒩[u_N](x,t;θ) - f(x,t).



The first term on the right-hand side in (<ref>) corresponds to the differential entropy of
r̂_θ(x,t), which does not affect the optimization with respect to
θ_f. So minimizing the KL divergence is equivalent to minimizing the cross
entropy between r̂_θ(x,t) and p_θ_f(x,t)

    H(r̂_θ(x,t), p_θ_f(x,t))    
    		= - ∬r̂_θ(x,t) log( p_poly(t;θ_f,1)
    		p_B-KRnet(x|t;θ_f,2)) dxdt,

Since the samples from r̂_θ(x,t) are not available, we approximate the cross entropy using the importance sampling technique:

    H(r̂_θ(x,t),p_θ_f(x,t))≈ -1/N_r∑_i=1^N_rr̂_θ(x_i, t_i)/p_poly(t_i;θ̂_f,1)
    		p_B-KRnet(x_i|t_i;θ̂_f,2)(log p_poly(t_i;θ_f,1) + log p_B-KRnet
    	(x_i|t_i;θ_f,2)),

where

    t_i∼ p_poly(·;θ̂_f,1),  x_i∼ p_B-KRnet(·|t_i;θ̂_f,2).

The choice of θ̂_f={θ̂_f,1,θ̂_f,2} will be specified as follows.

Once obtaining the well-trained parameters θ_f^* = {θ_f,1^*, θ_f,2^*}, one can refine initial training set
via adding new samples and update the weights in (<ref>). Specifically, note that the
residuals r_i(x;θ) in (<ref>) are only needed at time s_i,
we need to construct the corresponding discrete distributions.
For system (<ref>) in Galerkin projection, we use the following discrete distribution

    p_dis(s) ={[ ∫_Ωr_i^2(x;θ)dx/∑_i=0^N∫_Ωr_i^2(x;θ)dx,                           s=s_i,0≤ i≤ N;                                      0,                              otherwise. ].


For system (<ref>) in collocation projection, we simply use the following discrete distribution

    p_dis(s;θ_f,1^*) = {[               ∫_0^(s_0+s_1)/2p_poly(s;θ_f,1^*)ds,                                            s=s_0,; ∫_(s_i-1+s_i)/2^(s_i+s_i+1)/2p_poly(s;θ_f,1^*)ds,                                     s = s_i,0<i<N;             ∫_(s_N-1+s_N)/2^1p_poly(s;θ_f,1^*)ds,                                            s=s_N,;                                                0,                                        otherwise. ].
Then the weights in (<ref>) can be determined via the discrete distribution p_dis(·;θ_f,1^*).
For collocation projection, it is straightforward to refine the
training set. We first generate N_new samples {t_j}_j=1^N_new via the
well-trained model
p_dis(·;θ^*_f,1), then we generate the corresponding spatial samples
x_j∼ p_B-KRnet (·|t_j;θ^*_f,2) for each t_j.


After that we reorder newly generated data into the original spatial space

    x_j ∈ S^(i)_r,new,   if  t_j = s_i, ∀ i=0,⋯,N,  j=1,⋯,N_new.

For Galerkin projection, recall that we need to use the same training set for the basis function
with same support (see (<ref>)). Hence for each s_i in (<ref>), one
firstly identies its corresponding basis function and support, here it is denoted as Ĩ_i. Then we generate
spatial samples x_i,k∼ p_B-KRnet(·|t_i,k;θ_f,2^*), where t_i,k are Gaussian quadrature points in the interval
Ĩ_i. The total number of generated spatial points for each s_i should be proportional to the discrete weights p_dis(s_i). Again, we recoder newly
generated data into the original spatial space

    x_i,k∈ S^(i)_r,new,   ∀ i=0,⋯,N.


We are now ready to present our algorithms.
Let {S_r,0^(i)}_i=0^N, S_ic and S_bc be three sets of collocation points that are
uniformly sampled from Ω^N+1, Ω×{0} and ∂Ω× [0,T], respectively.
Using {S_r,0^(i)}_i=0^N, S_ic and S_bc, we minimize the empirical
loss (<ref>) to obtain u(x,t;θ^*, (1)). With
u(x,t;θ^*,(1)), we minimize the cross entropy (<ref>)
to get p_θ^*,(1)_f(x,t), where we use uniform samples for importance
sampling. To refine the training set, a new set
S_r,new={S^(i)_r,new}_i=0^N is generated according to the model p_θ_f^*,(1)(x,t), and new training set can be updated as
S_r,1 = S_r,0∪ S_r,new. Then we continue to update the approximate solution u(x,t;θ^*,(1)) using S_r,1 as the training set, which yields a refined model u(x,t;θ^*,(2)). Starting from k=2, we seek p_θ^*,k(x,t) using the previous approach. We repeat the procedure to obtain an adaptive algorithm (see Algorithm <ref>).




§ NUMERICAL EXPERIMENTS

In this section, we conduct some numerical experiments to demonstrate the effectiveness of the proposed method, including one convection equation, one Allen-Cahn equation, one two-dimensional and low regularity test problems, and two high-dimensional linear or nonlinear problems. Throughout all benchmarks, we will employ the fully-connected neural network equipped with hyperbolic tangent activation functions (Tanh) and initialized using the Glorot normal scheme <cit.>. All neural networks are trained via stochastic gradient descent using the Adam optimizer with default settings <cit.> and an exponential learning rate decay with a decay-rate of 0.9 every 1,000 training iterations.  All experiments are implemented by JAX <cit.>.

In order to test the validity of the method, we use the following relative L_2 error:

    err_L_2 = √(∑_i=1^num|u(x_i,t_i;θ) - u(x_i,t_i)|^2)/√(∑_i=1^num|u(x_i,t_i)|^2),

where num represents the total number of test points chosen randomly in the domain, and u(x_i,t_i;θ) and u(x_i, t_i) represent the predicted and the exact solution values, respectively.


 §.§ Convection equation

We start with a one-dimensional linear convection equation of the form

    ∂ u/∂ t + β∂ u/∂ x = 0,    (x,t)∈ [0, 2π] × [0,1],

subject to periodic boundary conditions and an initial condition u(0,x)=sin(x).
The precise solutions for varying values of β are depicted in Figure <ref>. One can observe that as β increases, the solution exhibits increasingly pronounced temporal variations.


We initially evaluate the performance with  linear finite element basis functions in Galerkin projection. The latent coefficients, denoted as ω_i(x;θ), are represented using a fully-connected neural network with tanh activation function, 4 hidden layers and 128 neurons per hidden layer. To simplify the training object, we strictly impose the periodic boundary conditions by embedding the input coordinates into Fourier expansion (see <ref>). Notice that the above equation is linear, we use the linear form of loss function (<ref>) and set N_segment=1. We create a uniform mesh of size 400 in the spatial computational domain [0,2π], yielding 400 initial points and 400 collocation points for enforcing the PDE residual. We proceed by training the resulting model via full-batch gradient descent using the Adam optimizer for 40,000 iterations. As shown in Figure <ref>, for fixed N, the relative L_2 gradually increases as β increases; for fixed β=50, the relative L_2 exhibits linear convergence concerning the number of basis functions (N). Particularly, when β=30, the proposed method attains a remarkable relative L_2 error of 2.85e-3.


Furthermore, we investigate the impact of linear and quadratic finite element basis functions on the performance of the proposed model. Specifically, we set β to 50 and vary the number of mesh elements N, ranging from 32 to 128. We then train the proposed model under same hyperparameter configurations. Figure <ref> and Table <ref> present a summary of the relative L_2 errors observed in the trained models. It is not surprising that the error diminishes when we replace linear basis functions with quadratic ones,   in accordance with the classical results of finite element theory.


To investigate the effects of increasing the number of basis functions, we compare the relative L_2 errors given by various settings. It is seen in Table <ref> that relative L_2 error decreases as N increases for all cases.

In Table <ref> we compare the performance of the proposed method with the DABG <cit.>
at different β. One can observe that for these two methods, when β is relatively small, as N increases, the error first decreases and then increases; when β is relatively large, the error consistently decreases as N increases. We conjecture that reason for this phenomenon is that when β is relatively small (i.e., the solution is smoother), small N can achieve good accuracy, while larger N may lead to greater optimization challenges. Another message from Table <ref> is that the proposed approach exhibits reduced sensitivity to temporal frequency variations in the solution compared to DABG.


We also compare the performance and computational time of the proposed method with causal PINNs. We set the same hyper parameters to compare these two methods, including the architecture of neural network, the initial learning rate, decay rate, the number of training points in space, and iterations. The ϵ in causal PINNs is set to {1e-2, 1e-1, 1e0, 1e1, 1e2}. All runtime statistics were computed on the same hardware, a Nvidia Tesla V100 w/ 32 GB memory. It is shown in Table <ref>  that as β increases, the proposed method demonstrates superior accuracy compared to causal PINNs. The proposed method (Galerkin projection) runs approximately 100 times faster than causal PINNs. Such a significant difference in efficiency is due to the distinct computational graphs adopted by these two methods. For the proposed method, because the residual originates from the linear system defined in (<ref>), we can pre-calculate the finite element sparse matrix along the time direction. We consider the computational cost of automatic differentiation in calculating the residual. For any given spatial point x, causal PINNs need n back-propagation computations for n temporal points, while our method only requires a single back-propagation computation and one matrix-vector multiplication of the above-mentioned finite element sparse matrix and the neural network output.




 §.§ Allen-Cahn equation

The next example aims to illustrate the effectiveness of time marching technology in our proposed method. Consider the following Allen-Cahn equation

    [   u_t -c_1^2 ∇ ^2 u +f(u)                         = 0,    x∈[-1, 1],t∈ [0,1],;                      f(u)                         =            c_2 (u^3 - u),;                    u(x,0)                         =             x^2 cos(π x),;                    u(1,t)                         =                 u(-1, t),;                 u_x(1, t)                         =               u_x(-1, t), ]

where c_1^2=0.0001 and c_2=5.
We take the number of linear mesh elements as 100 and represent the latent coefficients {ω_i(x;θ)} by a fully-connected neural network with tanh activation function, 4 hidden layers and 128 neurons per hidden layer. Similarly, we strictly impose the periodic boundary conditions by embedding the input coordinates into Fourier expansion (see <ref>). We create a uniform mesh of size 1,000 in the spatial computational domain [-1, 1], yielding 1,000 initial points and 1,000 collocation points for enforcing the PDE residual. We proceed by training the resulting model via full-batch gradient descent using the Adam optimizer for
the total 100,000 iterations. The resulting L_2 errors for different numbers of time segments are shown in Table <ref>. We observe the relative L_2 error significantly decreases as the number of time segments increases. One can see that the predicted solution achieves an excellent agreement with the ground truth, yielding an approximation error of 7.67e-3 measured in the relative L_2 norm.




In addition, we employ quadratic basis functions and set N to 30. All other parameters remain consistent with the previous configuration. In Table <ref> similar results are observed compared to the previous scenario.  Moreover, we also present a representative predicted solution in Figure <ref>, one can see that the predicted solution is in good agreement with the reference solution.




We summarize some relative L_2 errors for N_segment=2 in
Table <ref>. We observe that the collocation framework with piecewise spline functions has achieved the best accuracy.   Moreover, we compare the Galerkin projection with causal PINNs in Table <ref>. Again, we set the same hyper parameters to compare these two methods. The ϵ in causal PINNs is set to {1e-2, 1e-1, 1e0, 1e1, 1e2}. It is shown that the accuracy of these two methods is comparable. Note that due to the presence of non-linear term f(u), we must employ a Gaussian quadrature formula of a sufficient degree of exactness for the Galerkin projection in the time direction. Even so, the proposed method runs approximately 40 times faster than causal PINNs.




 §.§ High-dimensional linear problem


To demonstrate the effectiveness of the proposed approaches in tackling high-dimensional PDEs, we consider the following parabolic equation with non-constant coefficients

    u_t - ∇· (a(x)∇ u)=f,   Ω× (0,1], 
        u(x,0)=0,  Ω,                                          
        u = 0,  ∂Ω× [0,1],

with a(x)=1+ |x|^2/2. The domain is set to be the 20-D unit ball Ω = {x∈ℝ^20| |x|<1}, and the true solution is set to be

    u(x,t) = sin(sin (2πω t)(|x|^2 - 1)).

We set ω to be 3 and N_segment to be 1. For simplicity, here we only test the performance of the proposed approach with quadratic finite element basis functions. We represent the latent coefficients {ω_i(x;θ)} by a fully-connected neural network with activation function tanh, 5 hidden layers, and 128 neurons per hidden layer. We impose exactly the Dirichlet boundary conditions by transforming the output into the following form:

    ω_i(x;θ) = ω_i(x;θ) (|x| -1).

In addition, we set ũ_0(x;θ) to be zero such that the initial conditions are exactly satisfied. To obtain a set of training data for evaluating PDE residual, we randomly sample 100,000 collocation points in Ω. Since the problem is linear, equation (<ref>) is used to define the loss. We set the size of mini-batch to 5,000 and train the model via mini-batch stochastic descent with the Adam optimizer for 40,000 iterations. The corresponding results for different numbers of mesh elements are summarized in Table <ref>. We observe that the resulting relative L_2
error is 9.67e-4, which is more accurate than the ones in recent work <cit.> (3.07e-3 for the deep adaptive basis Galerkin approach and 4.54e-2 for the PINNs).





 §.§ High-dimensional nonlinear problem

In this case, we solve the Allen-Cahn equation

    u_t - Δ u + u^3 - u = f_AC,   Ω× [0,1], 
        u(x,0) = 0,   Ω,                               
        u = 0,   ∂Ω× [0,1].

The domain is set as Ω = {x∈ℝ^20: |x|<1}, and the true solution is set to be

    u(x,t) = sin(sin(2πω t)(|x|^2 -1)).

Here we set ω to be 3 and N_segment to be 1. For simplicity, We only test the performance of the proposed approach with quadratic finite element basis functions. We represent the latent coefficients {ũ_i(x;θ)} by a fully-connected neural network with tanh activation function, 5 hidden layers and 128 neurons per hidden layer. The initial condition and Dirichlet boundary condition can be exactly embedded into the network structure (as described in Section <ref>). We set K in equation (<ref>) to be 30 and randomly sample 100,000 collocation points for evaluating PDE residual. We set mini-batch to be 5,000 and proceed by training the resulting model via stochastic gradient descent using the Adam optimizer for 40,000 iterations. The obtained errors for different number of finite element N are shown in Table <ref>. The resulting relative L_2 error is 1.16e-3, which is more accurate than the recent ones in <cit.> (2.07e-3 for the deep adaptive basis Galerkin approach and 7.79e-2 for the PINNs).




 §.§ Low regularity test problems

Our last example aims to demonstrate the effectiveness of the proposed adaptivity strategy. We consider the following low-regularity equation

    u_t - Δ u =f  Ω× (0,T]                       
        u(x_1, x_2;0) = g(x_1, x_2)  Ω×{0},           
        u(x_1, x_2;t) = h(x_1,x_2;t)  ∂Ω× (0,T].

where Ω= [0, 1]^2, T=0.5 and the true solution is chosen as

    u(x_1,x_2,t)=exp(-β[(x_1 - t - 1/4)^2 + (x_2 - t - 1/4)^2]).

This solution has one peak at the point (t+1/4, t+1/4) and decreases rapidly away from (t+1/4,t+1/4), see Figure <ref>.


We first consider the case with β=200. We represent the latent coefficients
{ω̃_i(x;θ)} by a fully-connected neural network with
activation function tanh, 5 hidden layers and 128 neurons per hidden layer.
To simplify the training objective, we precisely embed the Dirichlet
boundary conditions into the neural network architecture through the
following transformation

    ω̃_i(x;θ) =  x_1x_2(1-x_1)(1-x_2)ω̃_i(x;θ).

For B-KRnet, we take 8 CDF coupling layers, and two fully connected layers with 32
neurons for each CDF coupling layer. The corresponding activation function is tanh. To assess the effectiveness of our adaptive sampling
approach, we generate a uniform meshgrid with size 256× 256×
	100 in [0,1]^2× [0,0.5]  and compute relative L_2 error on these grid points.

For collocation projection, we set the number of piecewise cubic spline basis functions to 20, and randomly sample 1,000 collocation points in Ω as our initial training set. We let N_adaptive=5 and
N_new=500. The number of epochs for training
u(x,t;θ) and p_θ_f(x,t) is set to 40,000 and
5,000, respectively.



In Figure <ref>, we plot the approximation errors with respect to the adaptivity iterations. It is clearly seen that the error rapidly decreases as the adaptivity iteration step k increases.
Figure <ref> shows the evolution of temporal distribution p_dis(t)
for k=1,3,5. It is seen that the largest density at the first adaptivity iteration is around terminal time. After the training set is refined,
the error profile becomes more flat at k=3 and 5. Similar results are observed in Figure <ref>, which indicates that the
distribution of residual becomes more uniform as adaptivity iteration increases.




For Galerkin projection, we set the number of piecewise quadratic basis functions
to 20, the remaining hyper-parameter settings remain the same. The corresponding numerical results, included in <ref>, are very similar to the results for collocation projection.

We also investigate the case when β=1000. Here we take the Galerkin projection as an example.

In Figure <ref>, we compare the performance of uniform sampling and adaptive sampling using the relative L_2 error, which demonstrates the necessity of adaptive sampling for this case.
Figure <ref> shows the evolution of temporal weights λ_i for k=1,3,5. The observed temporal weights exhibit an alternating high-low pattern, which can be attributed to the properties of the test basis function in time dimension. One also finds that the error profiles become more flat as the adaptivity iterations increase.
Figure <ref> shows the
generated samples of bounded KRnet with respect to adaptivity iterations k=1, 5.
It is seen that the largest density of generated samples at the first update is
around the peak of the reference solution, which is consistent with the residual-induced
distribution. Moreover, as k increases, we expect the tail of the residual-induced
distribution becomes heavier since the adaptivity tries to make the residual-induced
distribution more uniform. It is illustrated by the last update in the Figure
<ref>, which is consistent with previous findings
reported in <cit.>.









§ CONCLUSION

In this paper we have developed a hybrid numerical method for solving evolution partial differential equations
(PDEs) by merging the time finite element method with deep neural networks. The key idea of the proposed method is to represent the
solution as a tensor product comprising a series of
predetermined local finite element basis functions in
time and a sequence of unspecified neural networks in
space. Subsequently, we apply the Galerkin or collocation
formulation to the original evolution equation, eliminating
the temporal variable and resulting in a differential system
exclusively involving unknown neural networks with respect to
the spatial variable. Furthermore, to address the evolution
problems characterized by high dimensionality and low regularity, we have developed an adaptive
sampling strategy where the training set and the model are updated alternately to improve both efficiency and accuracy.
Numerical experiments have demonstrated the effectiveness of our proposed approaches.
Several issues deserve further investigation. First, a rigorous error analysis is still missing especially for the collocation projection. Second, while herein we only consider linear and quadratic finite element methods over a uniform mesh, the strategy can also be generalized to hp-finite element method in the temporal domain. Finally, the proposed methods need to be further investigated for convection-dominated problems.


§ ACKNOWLEDGEMENTS

T. Tang was partially supported by the Guangdong Provincial Key Laboratory of Interdisciplinary Research and Application for Data Science under UIC 2022B1212010006 and National Natural Science Foundation of China (Grants Nos. 11731006 and K20911001).


§ EXACT PERIODIC BOUNDARY CONDITIONS
 
Following the work from <cit.>, we can exactly enforce C^∞ periodic boundary conditions by applying a Fourier feature encoding to the spatial input of the network. The spatial encoding is

    v(x) = { 1,cos(ω x),sin(ω x) ,⋯,cos(Mω x),sin(Mω x)}

where ω =2π/L, L = x_max - x_min, and M is a non-negative integer representing the sinusoidal frequency of the input. In this work, we take M=10.



§ SUPPLEMENTARY FIGURES











