A Method for Representing Periodic Functions and Enforcing Exactly Periodic Boundary Conditions with Deep Neural Networks
AA Method for Representing Periodic Functions and EnforcingExactly Periodic Boundary Conditions with Deep Neural Networks
Suchuan Dong ∗ , Naxian NiCenter for Computational and Applied MathematicsDepartment of MathematicsPurdue University, USA(July 14, 2020) Abstract
We present a simple and effective method for representing periodic functions and enforcing exactly theperiodic boundary conditions for solving differential equations with deep neural networks (DNN). Themethod stems from some simple properties about function compositions involving periodic functions. Itessentially composes a DNN-represented arbitrary function with a set of independent periodic functionswith adjustable (training) parameters. We distinguish two types of periodic conditions: those imposingthe periodicity requirement on the function and all its derivatives (to infinite order), and those imposingperiodicity on the function and its derivatives up to a finite order k ( k (cid:62) C ∞ periodic conditions, and the latter C k periodic conditions. We define operations that constitutea C ∞ periodic layer and a C k periodic layer (for any k (cid:62) C ∞ (or C k ) periodic layer incorporated as the second layer automatically and exactly satisfies the C ∞ (or C k )periodic conditions. We present extensive numerical experiments on ordinary and partial differentialequations with C ∞ and C k periodic boundary conditions to verify and demonstrate that the proposedmethod indeed enforces exactly, to the machine accuracy, the periodicity for the DNN solution and itsderivatives. Keywords: periodic function, periodic boundary condition, neural network, deep neural network, periodicdeep neural network, deep learning
Deep neural networks (DNN) have emerged in the past few years as a promising alternative to the classicalnumerical methods (such as finite difference and finite element) for solving ordinary and partial differentialequations (PDE). DNN-based solvers transform the PDE solution problem into an optimization problem.They typically represent the unknown field function in terms of a deep neural network, thanks to the universalapproximation property of DNNs [16, 17, 4, 21]. Then these methods compute the solution by minimizinga loss function that consist of residual norms of the governing equation and also possibly of the boundaryand initial conditions in strong or weak forms (see e.g. [19, 20, 29, 33, 11, 27, 36, 30], among others). DNNsolutions are smooth analytical functions (depending on the activation function used) and, once the networkis trained, can be evaluated for the function value and its derivatives exactly at any point inside the domain.Boundary (and initial) conditions play a critical role in the solution of PDEs and make the problem wellposed [15, 31, 8, 7, 25]. To enforce the boundary conditions (BC) in DNN solvers, an often-used approach isthe penalty method, by incorporating a penalty term consisting of the residual norm of the boundary condi-tions into the loss function. The penalty method enforces the boundary conditions only approximately, andthe penalty coefficient strongly influences the DNN training and convergence [3]. Choosing an appropriateor near-optimal penalty coefficient is largely an art, usually conducted by trial and error. ∗ Author of correspondence. Email: [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] J u l nforcing the boundary conditions exactly with deep neural networks, if feasible, would be highly desir-able. In this case, the DNN is constructed such that the boundary conditions are automatically and exactlysatisfied. The constrained optimization problem for the PDE solution will then become less constrainedor unconstrained, which will greatly facilitate the DNN training. Enforcing exactly the boundary condi-tions with deep neural networks is, however, highly non-trivial. For Dirichlet and Neumann type boundaryconditions, several researchers have investigated the problem and promising techniques are available (seee.g. [19, 20, 23, 1]). In [19] the unknown solution is decomposed into two components. One componentsatisfies the Dirichlet/Neumann boundary conditions and has no training parameters, while the other com-ponent vanishes on the boundary and is represented by a neural network in the domain. In [20] the authorsdecompose the unknown solution into a function represented by a deep neural network plus a linear com-bination of radial basis functions. The combination coefficients of the radial basis functions are determinedby solving a linear system at every iteration of the DNN evaluation to satisfy the Dirichlet or Neumannboundary conditions. This process is understandably computationally expensive [20]. In [23] the authorsrewrite the solution into two parts, similar to [19], with one part satisfying the Dirichlet/Neumann boundaryconditions and the other part vanishing on the boundary but otherwise unconstrained. In order to dealwith complex domain boundaries, the authors of [23] introduce a multiplicative length factor in front of theunconstrained part, which heuristically represents the distance of a point to the domain boundary. Theseideas are further developed by [1], where the Dirichlet boundary condition is considered. In [1] the Dirichletboundary data extension and the distance function to the boundary are both represented by low-capacitydeep neural networks and pre-trained. This simplifies the implementation of the DNN solver. However, theenforcement of the boundary condition becomes only approximate. Similar ideas for the Dirichlet/Neumanntype boundary conditions have also appeared in more recent works, see e.g. [28] for the DNN simulation ofelastodynamic problems.Periodic boundary conditions are widely encountered in computational science of various areas, especiallywhen the physical domain involved in is infinite or homogeneous along one or more directions [9, 6]. In suchcases, usually only one cell will be computed in numerical simulations, and periodic boundary conditionsare imposed on the cell boundaries. With classical numerical methods, another often-used technique for thistype of problems is to express the unknown field function in terms of Fourier expansions, leading to what isknown as the Fourier spectral or pseudo-spectral method [2, 5, 10]. While both are referred to as periodicconditions, the periodicity requirements imposed by the numerical method, when Fourier expansions areused and when they are not, are different. With the use of Fourier expansions, the method seeks a smoothperiodic function as the solution to the governing equations, which automatically satisfies the periodicityfor the solution value and all its derivatives (to infinite order) on the cell boundaries. On the other hand,when Fourier expansions are not used in the method, the periodicity needs to be imposed explicitly on thecell boundaries, and this can only be imposed for the solution value and its derivatives up to a certain finiteorder. We will distinguish these two types of periodic boundary conditions in this work. We refer to theformer as the C ∞ periodic conditions, and the latter as the C k periodic conditions, where k (cid:62) C k periodic conditions usingthe penalty for moderate or large k values. In particular, it is practically impossible to impose the C ∞ periodic conditions with the penalty method.How to enforce exactly the C ∞ and C k (for any k (cid:62)
0) periodic conditions with deep neural networksis the focus of the current work. This problem seems to have barely been investigated before. To the bestof the authors’ knowledge, the only work close in theme to the current effort is perhaps [13], in which theauthors enforce the periodicity condition for the solution value only, by adopting a similar idea to [19] andconstructing a trial function with two parts. On part enforces the periodicity for the solution value, and theother part vanishes on the boundary and is represented by a DNN [13].2n the current paper we present a method for enforcing exactly the C ∞ and C k (for any k (cid:62)
0) periodicboundary conditions with deep neural networks. The DNN resulting from the current method, by design,automatically and exactly satisfies the C ∞ or C k (for any prescribed k ) periodic conditions. This method isbased on some simple properties about function compositions involving periodic functions (Lemmas 2.1 and2.2 in Section 2), and leverages the universal approximation power of deep neural networks. It essentiallycomposes a DNN-represented arbitrary function with a set of independent known periodic functions withadjustable parameters. We consider the feed-forward neural network architecture [14] in this work, and definethe operations that constitute a C ∞ periodic layer and a C k periodic layer. To enforce the C ∞ periodicconditions, one only needs to set the second layer of the DNN (i.e. the first hidden layer) as a C ∞ periodiclayer. To enforce the C k periodic conditions, one only needs to set the second layer of the DNN as a C k periodic layer. The C ∞ periodic layer constructs a set of independent C ∞ periodic functions with a user-prescribed period, based on sinusoidal functions, affine mappings and nonlinear activation functions (such as“tanh” and “sigmoid”). The C k periodic layer constructs a set of independent C k periodic functions, basedon the generalized Hermite interpolation polynomials, affine mappings and nonlinear activation functions.The output of the overall DNN, with the C ∞ (or C k ) periodic layer incorporated therein, automatically andexactly satisfies the C ∞ (or C k ) periodic conditions.The operations involved in the C ∞ and C k C ∞ periodic boundary conditions with deep neural networks; (ii) the methodfor exactly enforcing the C k (for any k (cid:62)
0) periodic boundary conditions with deep neural networks.The rest of this paper is structured as follows. In Section 2 we define the operations that constitute the C ∞ periodic layer and the C k periodic layer in one and higher dimensions, and establish that a deep neuralnetwork with these layers incorporated as the second layer exactly satisfies the C ∞ or C k periodic boundaryconditions for a given domain. In Section 3 we present extensive numerical experiments on periodic functionapproximations, and on solving the Helmholtz equations in one and two dimensions, the diffusion equation,and the wave equation, with C ∞ and C k periodic boundary conditions. We demonstrate numerically thatthe proposed method enforces exactly the periodic boundary conditions, to the machine accuracy, for theDNN solution and its corresponding higher derivatives. Section 4 then concludes the presentation with someclosing remarks. C ∞ and C k Periodic Conditions
Consider a smooth periodic function f ( x ) with period L defined on the real axis, f ( x + L ) = f ( x ) , ∀ x ∈ ( −∞ , ∞ ) . (1)Now restrict f ( x ) to a finite interval [ a, b ], where b − a = L . Then f satisfies the following relations on theboundaries: f ( a ) = f ( b ) , f (cid:48) ( a ) = f (cid:48) ( b ) , f (cid:48)(cid:48) ( a ) = f (cid:48)(cid:48) ( b ) , . . . , f ( m ) ( a ) = f ( m ) ( b ) , . . . (2)We refer to the conditions in (2) as the C ∞ periodic conditions. Hereafter we will refer to a smooth function f ( x ) satisfying these conditions as a C ∞ periodic function on [ a, b ], or simply a periodic function.In practice, the function may not be smooth and the conditions in (2) may only be required for thederivatives up to a finite order k ( k (cid:62) f ( l ) ( a ) = f ( l ) ( b ) , (cid:54) l (cid:54) k, (3)3here f (0) ( x ) = f ( x ) by convention. We refer to the ( k + 1) conditions in (3) as the C k periodic conditions.With a slight abuse of notation, we will refer to a function f ( x ) satisfying the conditions (3) as a C k periodicfunction on [ a, b ].Our goal here is to devise a method for representing C ∞ and C k periodic functions with deep neuralnetworks such that, by design, the output of the DNN automatically and exactly satisfies the C ∞ or C k periodic conditions. Such neural networks will be referred to as C ∞ or C k periodic deep neural networks.When solving a boundary value problem or initial/boundary value problem together with the C ∞ or C k periodic boundary conditions, one can use the method developed herein to construct periodic DNNs as thetrial functions that automatically take into account the periodic boundary conditions. C ∞ Periodic Conditions with DNN
We present a method below for representing C ∞ periodic functions and enforcing exactly the C ∞ periodicconditions with DNN. The method is based on the following property about function compositions involvingperiodic functions. Lemma 2.1.
Let v ( x ) be a given smooth periodic function with period L on the real axis, i.e. v ( x + L ) = v ( x ) for all x ∈ ( −∞ , ∞ ) , and f(x) denote an arbitrary smooth function. Define u ( x ) = f ( v ) = f ( v ( x )) . Then u ( x + L ) = u ( x ) , ∀ x ∈ ( −∞ , ∞ ); (4a) u ( l ) ( a ) = u ( l ) ( b ) , l = 0 , , , . . . (4b) where a and b denote two real numbers with b − a = L . This lemma can be proven by straightforward verifications.We seek a DNN representation for an arbitrary smooth periodic function with a prescribed period L ,such that the C ∞ periodic conditions are automatically satisfied. In light of Lemma 2.1, our basic idea forthe representation is to compose an arbitrary DNN-presented function, together with a set of independentknown periodic functions with period L and adjustable (training) parameters. Let us first use a single knownperiodic function with prescribed period L for illustration. We consider the sinusoidal functions, p ( x ) = A cos( ωx + φ ) + c, with ω = 2 πL , (5)where the constants A , c and φ are scalar adjustable (training) parameters. Here ω is a fixed constant asgiven above and ensures that p ( x ) has a period L . Let σ ( · ) denote a nonlinear activation function (such as“tanh” or “sigmoid”). We define v ( x ) = σ ( p ( x )) = σ ( A cos( ωx + φ ) + c ) . (6)This step is crucial. The nonlinear function σ ( · ) will generate higher-frequency components in the output. Sowhile p ( x ) has a single frequency ω , v ( x ) contains not only the frequency ω , but components with other andhigher frequencies, all with a common period L . Finally, we consider an arbitrary function f ( x ) representedby a DNN, and define u ( x ) = f dnn ( v ( x )) , (7)where f dnn denotes the DNN-presented arbitrary function. By Lemma 2.1, this u ( x ) satisfies the C ∞ periodicconditions exactly.In practice, we would like to compose the DNN-represented arbitrary function f dnn ( · ) with a set ofindependent periodic functions v ( x ) as defined above with adjustable parameters. This leads to the idea ofa periodic layer with multiple nodes (neurons) within the layer. We consider the feed-forward deep neuralnetwork architecture [14] in the current work. Figure 1 illustrates the idea of the current DNN with a sketch.Let x denote the input layer to the network, and u ( x ) denote the output layer of the network. We use thesecond layer (behind the input x ) to implement the set of independent known periodic functions (copies of v ( x )) with period L and adjustable parameters, so that the output of the network u ( x ) satisfies exactly the C ∞ periodic conditions (4). We refer to such a layer as a C ∞ periodic layer, or simply a periodic layer.4 .... ..... ..... ∞ Figure 1: Sketch of a feed-forward deep neural network with a C ∞ or C k periodic layer incorporated as thesecond layer. .......... ∞ (a) .......... ∞ ... .......... (b)Figure 2: Sketch illustrating the internal structures of (a) 1D, and (b) 2D C ∞ or C k periodic layers.The operations within the C ∞ periodic layer are defined as follows (see Figure 2(a)). Let L p ( m, n ) denotethe C ∞ periodic layer, where n denotes the number of nodes in the output of this layer and m denotes thesize of the set of independent periodic functions v ( x ). Here both m and n are hyper-parameters of the C ∞ periodic layer L p ( m, n ). The operations within L p ( m, n ) are defined by: v i ( x ) = σ ( A i cos( ωx + φ i ) + c i ) , (cid:54) i (cid:54) m ; (8a) q j ( x ) = σ (cid:32) m (cid:88) i =1 v i ( x ) W ij + B j (cid:33) , (cid:54) j (cid:54) n. (8b)In these equations q j ( x ) (1 (cid:54) j (cid:54) n ) are the output of this layer, and the fixed constant ω is given in (5)for a prescribed period L . σ ( · ) is the nonlinear activation function, and it is used twice in this layer. Thetraining parameters of L p ( m, n ) are the constants A i , φ i , c i , W ij and B j , with 1 (cid:54) i (cid:54) m and 1 (cid:54) j (cid:54) n C ∞ periodic layer as described above. The C ∞ periodic layer is implemented asa user-defined Tensorflow/Keras layer, and it can be used in the same way as the built-in core Keras layers. Remark 2.1.
The operations of the C ∞ periodic layer defined by (8) can be extended to two, three andhigher dimensions (2D/3D) in a straightforward way. Here we outline the idea with two dimensions only(see Figure 2(b)). Let ( x , x ) denote the coordinates in two dimensions, and u ( x , x ) denote a smooth eriodic function to be approximated, with properties u ( x + L , x ) = u ( x , x ) , u ( x , x + L ) = u ( x , x ) , ∀ x , x ∈ ( −∞ , ∞ ) , (9) where L and L are the periods in the x and x directions, respectively. Equivalently, we can write theperiodicity conditions in terms of a single periodic cell [ a , b ] × [ a , b ] , ∂ α ∂x α u ( a , x ) = ∂ α ∂x α u ( b , x ) , ∀ x ∈ [ a , b ] ,∂ α ∂x α u ( x , a ) = ∂ α ∂x α u ( x , b ) , ∀ x ∈ [ a , b ] , α = 0 , , , . . . (10) where a , b , a and b are given constants satisfying b − a = L and b − a = L . In this case, we definethe 2D periodic layer, L Dp ( m, n ) , with the following operations: v i ( x ) = σ ( A i cos( ω x + φ i ) + c i ) , (cid:54) i (cid:54) m ; (11a) v i ( x ) = σ ( A i cos( ω x + φ i ) + c i ) , (cid:54) i (cid:54) m ; (11b) q j ( x , x ) = σ (cid:32) m (cid:88) i =1 v i ( x ) W (1) ij + m (cid:88) i =1 v i ( x ) W (2) ij + B j (cid:33) , (cid:54) j (cid:54) n. (11c) In these equations, m and n are hyper-parameters of the layer L Dp , q j ( x , x ) ( (cid:54) j (cid:54) n ) denote the outputof this layer, and the constants ω and ω are defined by ω = 2 πL , ω = 2 πL , (12) with prescribed periods ( L , L ) . The training parameters of L Dp ( m, n ) consist of the constants: A i , A i , φ i , φ i , c i , c i , W (1) ij , W (2) ij , B j , (cid:54) i (cid:54) m, (cid:54) j (cid:54) n. By composing an arbitrary DNN-represented function with the 2D C ∞ periodic layer defined above, we attainan overall DNN whose output automatically and exactly satisfies the 2D C ∞ periodic conditions. The C ∞ periodic layer for three and higher dimensions can be defined in a similar way. Remark 2.2.
In two or higher dimensions, if the C ∞ periodic conditions are imposed only in some (notall) directions, the C ∞ periodic layer as defined above can be modified in a simple way to accommodate thesituation. For example, consider the 2D C ∞ periodic layer defined in (11) and suppose that the C ∞ periodicconditions are imposed only in the x direction with period L , but not in the x direction. In this case, wecan retain the equations (11a) and (11c) , and replace (11b) by the following equation v i ( x ) = σ ( A i x + c i ) , (cid:54) i (cid:54) m, (13) where the constants A i and c i are the training parameters. The modified 2D periodic layer consisting ofequations (11a) , (13) and (11c) , when composed with a DNN-represented arbitrary function, will give rise toan overall DNN that automatically and exactly satisfies the C ∞ periodic conditions in the x direction. C k Periodic Conditions with DNN
We present in this subsection a method for representing C k periodic functions and enforcing exactly the C k periodic conditions (for any k (cid:62)
0) with DNN. The method is based on the following simple property aboutfunction compositions involving C k periodic functions: Lemma 2.2.
Let v ( x ) ( x ∈ [ a, b ] ) denote a given function with continuous derivatives up to the order k andsatisfying the following property, v ( l ) ( a ) = v ( l ) ( b ) , (cid:54) l (cid:54) k. (14) Let f ( x ) denote an arbitrary function defined on the real axis with continuous derivatives up to the order k .Define u ( x ) = f ( v ) = f ( v ( x )) ( x ∈ [ a, b ] ). Then u ( l ) ( a ) = u ( l ) ( b ) , (cid:54) l (cid:54) k. (15)6 roof. By induction one can show that u ( m ) ( x ) = g ( v, v (cid:48) , v (cid:48)(cid:48) , . . . , v ( m ) ) for 0 (cid:54) m (cid:54) k . In other words, u ( m ) ( x ) depends on x only through v and its derivatives. Equation (15) follows immediately from thisrelation and the conditions (14).We seek a DNN representation for an arbitrary C k periodic function on [ a, b ], such that the C k periodicconditions are automatically and exactly satisfied. In light of Lemma 2.2, our basic idea for the representationis to compose an arbitrary function represented by a DNN, together with a set of independent known C k periodic functions with adjustable (training) parameters. To construct a C k periodic function v ( x ) on [ a, b ]in Lemma 2.2, i.e. satisfying the conditions (14), we note that these conditions are reminiscent of the Hermiteinterpolation conditions. So the Hermite interpolation polynomial of degree at most (2 k + 1) can be usedto construct v ( x ). Once v ( x ) is obtained, we compose an arbitrary DNN-represented function f with v ( x ),and the resultant function satisfies the C k periodic conditions exactly.Let us now use a single C k periodic function v ( x ) ( x ∈ [ a, b ]) to illustrate the idea in some detail. Let s i (0 (cid:54) i (cid:54) k ) denote ( k +1) adjustable (training) parameters. Let h ( x ) denote the unique Hermite interpolationpolynomial of degree at most (2 k + 1) that satisfies the following (2 k + 2) interpolation conditions: h ( a ) = s , h ( b ) = s ; h (cid:48) ( a ) = s , h (cid:48) ( b ) = s ; · · · h ( k ) ( a ) = s k , h ( k ) ( b ) = s k . (16)The Newton form for h ( x ) can be computed based on the divided differences, and the explicit Lagrange formfor h ( x ) is available in e.g. [34, 35]. We then define (cid:40) p ( x ) = h ( x ) + ( r + r x )( x − a ) k +1 ( x − b ) k +1 ,v ( x ) = σ ( p ( x )) = σ (cid:0) h ( x ) + ( r + r x )( x − a ) k +1 ( x − b ) k +1 (cid:1) , (17)where r and r are two additional adjustable parameters, and σ ( · ) is a nonlinear activation function(e.g. “tanh” or “sigmoid”). It is straightforward to verify that the v ( x ) (and also p ( x )) given in (17)satisfies the conditions (14). The explicit forms for p ( x ) in (17) corresponding to k = 0, 1 and 2 are givenby: p ( x ) = s + ( r + r x )( x − a )( x − b ) , for k = 0; p ( x ) = s + s ( x − a )( b − x )( a + b − x ) + ( r + r x )( x − a ) ( x − b ) , for k = 1; p ( x ) = ξ ( x − a ) (cid:104) s + ( − s ξ + s )( x − b ) + (cid:16) s ξ − s ξ + s (cid:17) ( x − b ) (cid:105) + ξ ( b − x ) (cid:104) s + (3 s ξ + s )( x − a ) + (cid:16) s ξ + 3 s ξ + s (cid:17) ( x − a ) (cid:105) + ( r + r x )( x − a ) ( x − b ) , for k = 2 , where ξ = b − a . Let f dnn ( x ) denote an arbitrary function represented by a deep neural network. With v ( x )given by (17), we finally define u ( x ) = f dnn ( v ( x )) . (18)Then by Lemma 2.2 u ( x ) satisfies exactly the C k periodic conditions (15).In practice, we would like to compose the DNN-represented arbitrary function f dnn ( · ) with a set of inde-pendent C k periodic functions v ( x ) with adjustable parameters. This leads to the idea of a C k periodic layerwith multiple nodes (neurons) within the layer. Consider again the feed-forward neural network architecture,and let x denote the input and u ( x ) denote the output of the network. Analogous to the C ∞ periodic layerin Section 2.2, we define a C k periodic layer below, and use it as the second layer (behind the input x ) ofthe network to implement the set of independent C k periodic functions (copies of v ( x )) and enforce the C k periodic conditions. Figure 1 sketches the DNN with a C k periodic layer incorporated as the second layer.The operations within the C k periodic layer are defined as follows (see Figure 2(a)). Let L C k ( m, n )denote the C k periodic layer, where n denotes the number of nodes in the output of this layer and m denotes7he size of the set of independent C k periodic functions v ( x ). Both m and n are hyper-parameters of thislayer. Given the input x , we compute the output q j ( x ) (1 (cid:54) j (cid:54) n ) of the C k periodic layer L C k ( m, n ) by: v i ( x ) = σ (cid:0) h i ( x ) + ( r i + r i x )( x − a ) k +1 ( x − b ) k +1 (cid:1) , (cid:54) i (cid:54) m ; (19a) q j ( x ) = σ (cid:32) m (cid:88) i =1 v i ( x ) W ij + B j (cid:33) , (cid:54) j (cid:54) n. (19b)In these equations σ ( · ) is the nonlinear activation function, and h i ( x ) (1 (cid:54) i (cid:54) m ) are the Hermite interpo-lation polynomials of degree at most (2 k + 1) satisfying the conditions h ( l ) i ( a ) = s li , h ( l ) i ( b ) = s li , (cid:54) l (cid:54) k, (cid:54) i (cid:54) m. (20)The constant parameters involved in these equations, r i , r i , W ij , B j , s li , for 1 (cid:54) i (cid:54) m , 1 (cid:54) j (cid:54) n and 0 (cid:54) l (cid:54) k , are the training parameters of the C k periodic layer L C k ( m, n ). By incorporating the C k periodic layer defined by (19) as the second layer, the resultant DNN automatically and exactly satisfies the C k periodic conditions with its output. Remark 2.3.
The operations of the C k periodic layer L C k ( m, n ) defined by (19) can be extended to two,three and higher dimensions in a straightforward fashion. Here we use two dimensions only to illustrate theidea (see Figure 2(b)). Let x and x ( x ∈ [ a , b ] , x ∈ [ a , b ] ) denote the coordinates in two dimensions,and u ( x , x ) denote a 2D C k periodic function, satisfying the C k periodic conditions: ∂ α ∂x α u ( a , x ) = ∂ α ∂x α u ( b , x ) , ∀ x ∈ [ a , b ] ,∂ α ∂x α u ( x , a ) = ∂ α ∂x α u ( x , b ) , ∀ x ∈ [ a , b ] , α = 0 , , . . . , k. (21) In this case, we define the 2D C k periodic layer L DC k ( m, n ) with the following operations: v i ( x ) = σ (cid:16) h i ( x ) + (cid:16) r (1)0 i + r (1)1 i x (cid:17) ( x − a ) k +1 ( x − b ) k +1 (cid:17) , (cid:54) i (cid:54) m ; (22a) v i ( x ) = σ (cid:16) h i ( x ) + (cid:16) r (2)0 i + r (2)1 i x (cid:17) ( x − a ) k +1 ( x − b ) k +1 (cid:17) , (cid:54) i (cid:54) m ; (22b) q j ( x , x ) = σ (cid:32) m (cid:88) i =1 v i ( x ) W (1) ij + m (cid:88) i =1 v i ( x ) W (2) ij + B j (cid:33) , (cid:54) j (cid:54) n. (22c) In the above equations, m and n are the hyper-parameters of this layer, q j ( x , x ) ( (cid:54) j (cid:54) n ) denote theoutput of this layer, and h i ( x ) and h i ( x ) are the Hermite interpolation polynomials of degree at most (2 k + 1) satisfying the conditions: (cid:40) h ( l )1 i ( a ) = s (1) li , h ( l )1 i ( b ) = s (1) li , (cid:54) l (cid:54) k, (cid:54) i (cid:54) m ; h ( l )2 i ( a ) = s (2) li , h ( l )2 i ( b ) = s (2) li , (cid:54) l (cid:54) k, (cid:54) i (cid:54) m. (23) The constant parameters involved in the above equations, s (1) li , s (2) li , r (1)0 i , r (2)0 i , r (1)1 i , r (2)1 i , W (1) ij , W (2) ij , B j , (cid:54) l (cid:54) k, (cid:54) i (cid:54) m, (cid:54) j (cid:54) n, are the training parameters of the layer L DC k ( m, n ) . By using the 2D C k periodic layer L DC k ( m, n ) as thesecond layer of a DNN and with ( x , x ) as the input, the resultant DNN will automatically and exactlysatisfy the 2D C k periodic boundary conditions (21) . The C k periodic layer for three and higher dimensionscan be defined in a similar way. Remark 2.4.
In two and higher dimensions, if the C k periodic conditions are only imposed in some (notall) directions, the C k periodic layer as defined above can be modified in a simple way to accommodate thesituation. The modification is similar to what is discussed in Remark 2.2 for the modified C ∞ periodic layer. or illustration, let us consider the 2D C k periodic layer defined by (22) , and suppose that the C k periodicconditions are imposed only in the x direction, not in the x direction. In this case, we can retain theequations (22a) and (22c) , and replace equation (22b) by the following equation for v i ( x ) , v i ( x ) = σ (cid:16) r (2)0 i + r (2)1 i x (cid:17) , (24) where the constants r (2)0 i and r (2)1 i are training parameters. The modified 2D C k periodic layer consisting ofequations (22a) , (24) and (22c) , when used as the second layer of a DNN, will impose exactly the C k periodicconditions in the x direction in the output of this DNN. We present several numerical examples in what follows to demonstrate the effectiveness of the methodpresented in the previous section. We consider the approximation of periodic functions, and the solutionof the Helmholtz equation, the unsteady diffusion equation and the wave equation, together with periodicboundary conditions, using deep neural networks. We employ a variant of the deep Galerkin method [33]for solving the differential equations with DNN and also for enforcing the initial conditions with unsteadyproblems. The periodic boundary conditions (BC) are dealt with based on the method from Section 2.Note that with the current method the periodic boundary conditions are satisfied automatically by theDNN. Therefore there is no need to account for the periodic boundary conditions in the loss function. Theapplication codes for all the tests reported here are implemented using Tensorflow/Keras and Python, witheither Adam [18] or L-BFGS [26] as the optimizer. We employ the hyperbolic tangent (“tanh”) as thenonlinear activation function in all the tests.
Let us look into the DNN approximation of periodic functions using the method developed in Section 2. Weemploy three different functions to illustrate the performance characteristics of the method: a C ∞ periodicfunction, a C periodic function, and a non-periodic function. Note that in the case of the non-periodicfunction, we are essentially seeking a periodic function approximation of the non-periodic function.Consider first the function u ( x ) = sin(2 πx + 0 . π ) + cos(9 πx − . π ) − πx + 0 . π ) (25)on the domain Ω = { x | (cid:54) x (cid:54) } . This is a C ∞ periodic function on this domain. We would like toapproximate u ( x ) using a C ∞ periodic DNN, and using a C k periodic DNN with k = 0 and k = 1.To approximate u ( x ), we employ a feed-forward neural network [14] as illustrated in Figure 1. Theinput to the network is x , and the output is the approximation u ( x ). We use 3 hidden layers in between,each with a width of 30 neurons. The hyperbolic tangent (“tanh”) function is used as the activation functionfor all the hidden layers, and no activation is applied on the output layer. As discussed in Section 2, weset the second layer of the network (i.e. the first hidden layer) as a C ∞ periodic layer for the C ∞ periodicapproximation, and as a C k periodic layer for the C k periodic approximation. More specifically, we employa L p ( m, n ) with m = 11 and n = 30 for the C ∞ periodic layer, and a L C k ( m, n ) with m = 11 and n = 30for the C k -periodic layer with k = 0 and 1. In other words, a set of 11 independent periodic functions v ( x )has been used within the C ∞ periodic layer and the C k periodic layer. For the C ∞ periodic layer L p ( m, n ),the constant ω in equation (8a) is set to, according to equation (5), ω = 2 πV Ω = 2 π π, (26)where V Ω = (cid:82) Ω dx = 2 is the size of the domain Ω. 9arameter value parameter valuehidden layers depth=3, width=30 N e C ∞ periodic layer L p (11 , Q
30 ( C ∞ periodic DNN),or C k periodic layer L C k (11 ,
30) or 40 ( C k periodic DNN)activation tanh optimizer Adammaximum epochs 10000 learning rate 1 e − x ei (0 (cid:54) e (cid:54) N e −
1, 0 (cid:54) i (cid:54) Q −
1) label data u ( x ei ), or u ( x ei ), or u ( x ei ) x ei Gauss-Lobatto-Legendre quadrature ω π/V Ω Table 1: Function approximation: DNN and simulation parameters.We minimize the following loss function with this DNN,Loss = 1 V Ω (cid:90) Ω | u ( x ) − u ( x ) | dx = 1 V Ω N e − (cid:88) e =0 (cid:90) Ω e | u ( x ) − u ( x ) | dx = 1 V Ω N e − (cid:88) e =0 Q − (cid:88) i =0 | u ( x ei ) − u ( x ei ) | J e w i , (27)where N e is the number of elements (i.e. sub-intervals) we have partitioned the domain Ω into in order tocompute the integral, Q is the number of quadrature points within each element, Ω e denotes the sub-intervaloccupied by the element e (0 (cid:54) e (cid:54) N e − x ei (0 (cid:54) i (cid:54) Q −
1) are the Gauss-Lobatto-Legendre quadraturepoints within Ω e for 0 (cid:54) e (cid:54) N e − J e is the Jacobian of the element Ω e with respect to the standard element[ − , w i (0 (cid:54) i (cid:54) Q −
1) are the quadrature weights associated with the Gauss-Lobatto-Legendrequadrature points. In the numerical experiments we have employed 3 elements ( N e = 3) to partition thedomain Ω, with Ω = [0 , . = [0 . , .
4] and Ω = [1 . , Q = 30)within each element for the C ∞ periodic DNN, and 40 quadrature points ( Q = 40) within each elementfor the C k ( k = 0 ,
1) periodic DNN. The input data to the network consist of all the quadrature points x ei (0 (cid:54) i (cid:54) Q −
1, 0 (cid:54) e (cid:54) N e − u ( x ei ). The Adam optimizer has been usedto train the network for 10000 epochs for each case, with the learning rate fixed at the default value 10 − .The options of “early stopping” and “restore to best weight” have been used in Keras during the trainingprocess. The main parameters for the DNN and the simulations are summarized in Table 1.Figure 3 shows the training histories of the loss function corresponding to the C ∞ , C and C periodicDNN approximations. The loss function decreases rather slowly as the training begins. Then we observea short stage when the loss function decreases sharply. After that, the reduction in the loss function slowsdown again, resulting in a long tail in the training history curve. These characteristics seem to be commonto the DNN training for all the problems we have considered in this work. During the stage with slowreduction in the loss function (long tail), we observe that the loss value fluctuates from time to time duringthe training, resulting in a sequence of spikes in the training history curves (Figure 3). This is likely due tothe fixed learning rate employed here. Reducing the learning rate gradually as the training progresses willlikely reduce the loss fluctuations. In spite of the spikes, one can observe the trend of decreasing loss as thetraining proceeds. Since we have turned on the options of “early stopping” and “restore to best weight” inKeras, the spikes in the training curves have no effect on the simulation results we have obtained here.In Figure 4 we compare the DNN approximation results of u ( x ) (top row) obtained with C ∞ (plot (a)), C (plot (b)), and C (plot (c)) periodic conditions, together with the exact function u ( x ). The distributionsof the absolute error, | u ( x ) − u ( x ) | , corresponding to these approximations are shown in Figures 4(d,e,f),respectively. It is observed that the DNN approximations computed with all three methods agree well withthe exact function u ( x ), and the approximation function curves overlap with the exact function curve.In Table 2 we list the values of the function u ( x ) and its derivatives du dx and d u dx on the domainboundaries x = 0 and x = 2, obtained from the C ∞ , C and C periodic DNN approximations, as well asfrom the exact u ( x ) function given in (25). The function derivatives have been computed based on auto-differentiation. Once the DNN is trained, the derivatives computed in this way are exact values correspondingto the given DNN representation. We have shown 14 significant digits (double precision) for the values inthis table. It is evident that the C ∞ periodic DNN enforces exactly, to the machine accuracy, the periodicity10 poch l o ss (a) epoch l o ss (b) epoch l o ss (c)Figure 3: Approximation of the periodic function u ( x ): training loss histories corresponding to (a) C ∞ ,(b) C , and (c) C periodic DNN approximations. x S o l u t i on s DNNExact (a) x S o l u t i on s DNNExact (b) x S o l u t i on s DNNExact (c) x E rr o r (d) x E rr o r (e) x E rr o r (f)Figure 4: Approximation of periodic function u ( x ): comparison of DNN approximations (top row) andtheir errors (bottom row) with (a,d) C ∞ , (b,e) C , (c,f) C periodic conditions. The exact function is shownin the plots (a,b,c) for reference. 11 ∞ periodic DNN C periodic DNN C periodic DNN Exact value u (0) 9.2478271387350e-01 9.5998617007276e-01 9.3202741230659e-01 9.3667924347381e-01 u (2) 9.2478271387350e-01 9.5998617007276e-01 9.3202741230659e-01 9.3667924347381e-01 u (cid:48) (0) -9.2010196979120e+00 -1.1103158994625e+01 -9.3797238469930e+00 -9.2086781968974e+00 u (cid:48) (2) -9.2010196979120e+00 -7.1832211211812e+00 -9.3797238469930e+00 -9.2086781968972e+00 u (cid:48)(cid:48) (0) 6.0998633006437e+01 8.4031868694133e+01 1.0357079019194e+02 4.4301828505852e+01 u (cid:48)(cid:48) (2) 6.0998633006437e+01 8.6538176666730e+01 3.5624115949487e+01 4.4301828505854e+01Table 2: Approximations of the periodic function u ( x ): Values of the function and its first and secondderivatives at the left/right domain boundaries ( x = 0 ,
2) from the DNN approximations with C ∞ , C and C periodic conditions and from the exact function. 14 significant digits (double precision) are listed toshow that the current method enforces the periodic conditions exactly. The boxes highlight the values thatmis-match on the boundaries. x S o l u t i on DNNExact x (a) x S o l u t i on DNNExact x (b) x S o l u t i on DNNExact x (c) x E rr o r (d) x E rr o r (e) x E rr o r (f)Figure 5: Approximation of C periodic function u ( x ): Comparison of approximation results (top row)and their errors (bottom row) obtained with (a,d) C ∞ , (b,e) C , and (c,f) C periodic DNNs. The exactfunction is included for comparison. The insets are the magnified views near x = 0.for the function as well as its derivatives. On the other hand, the C periodic DNN enforces exactly theperiodicity only for the function value, and the C periodic DNN enforces exactly the periodicity only forthe function value and the first derivative. These numerical results have verified the analyses about thesemethods in Section 2.We next consider the function u ( x ) = sin πx { x | (cid:54) x (cid:54) } . This is a C periodic function on this domain, with u (0) = u (2) and u (cid:48) (0) (cid:54) = u (cid:48) (2). We would like to approximate u ( x ) with C ∞ and C k ( k = 0 and 1) periodic DNNs.We employ the same DNN and simulation parameters to approximate u ( x ) as for u ( x ); see Table 1. Theloss function is given by (27), with u ( x ) replaced by u ( x ). Figure 5 is a comparison of the approximationresults and their errors obtained with C ∞ , C and C periodic DNNs. The exact function u ( x ) has alsobeen included for comparison. The C periodic DNN produces results that are considerably more accuratethan the other two methods, as expected. The C ∞ and C periodic DNN approximations produce accurateresults in the bulk of the domain, but exhibit larger errors near/at the domain boundaries. The C ∞ and C periodic conditions appear to have the tendency of bending the function curve near the boundaries toachieve periodicity for the derivatives; see the insets of Figures 5(a,b,c).To verify the periodicity of the DNN approximations on the boundaries, we list in Table 3 the valuesof the approximated u ( x ) and its derivatives (up to order two) on the domain boundaries x = 0 and 2from different approximations and from the exact function u ( x ). Again 14 significant digits have been12 ∞ periodic DNN C periodic DNN C periodic DNN Exact value u (0) 2.0077934937638e-02 -5.6212180516792e-05 1.0675453567003e-02 0 u (2) 2.0077934937638e-02 -5.6212180516792e-05 1.0675453567003e-02 0 u (cid:48) (0) 2.0046381097937e-03 1.5738634304121e+00 -8.9688048016955e-02 1.5707963267949e+00 u (cid:48) (2) 2.0046381097902e-03 -1.5651374257149e+00 -8.9688048016955e-02 -1.5707963267949e+00 u (cid:48)(cid:48) (0) 4.8956024103701e+01 -3.0860019634346e-02 1.1444431255354e+02 0 u (cid:48)(cid:48) (2) 4.8956024103700e+01 1.7388612411815e-01 1.0788094853391e+02 0Table 3: Approximation of C periodic function u ( x ): Values of the function and its derivatives at theleft/right domain boundaries ( x = 0 and x = 2) from the C ∞ , C and C periodic DNN approximationsand from the exact function. 14 significant digits are listed to demonstrate that the current method enforcesthe periodic conditions exactly. x S o l u t i on DNNExact x (a) x S o l u t i on DNNExact x (b) x S o l u t i on DNNExact x (c) x E rr o r (d) x E rr o r (e) x E rr o r (f)Figure 6: Approximation of non-periodic function u ( x ): Approximation results (top row) and their errors(bottom row) obtained with (a,d) C ∞ , (b,e) C , (c,f) C , periodic DNNs. The exact profile of the function u ( x ) is also included. The insets are magnified views near x = 0.shown for each value. It is observed that the current methods indeed enforce exactly the periodicity forthe approximation function and its derivatives on the boundaries as expected. In the C ∞ periodic DNNapproximation, the function and its derivatives (up to order 2 considered here) have identical values on thetwo boundaries. In contrast, with the C periodic DNN approximation only the function value is identicalon the boundaries, and with the C periodic DNN approximation the function and the first derivative haveidentical values on the two boundaries.We finally consider a non-periodic function, u ( x ) = cos πx , (29)on the domain Ω = { x | (cid:54) x (cid:54) } . We would like to approximate this function using C ∞ and C k ( k = 0 , u ( x )and u ( x ) (see Table 1), except for the number of quadrature points within each element. Here for u ( x ) weemploy Q = 40 with the C ∞ periodic DNN, and Q = 50 with the C and C periodic DNNs.Figure 6 shows the approximation functions u ( x ) and their errors obtained with C ∞ , C and C periodicDNNs. In the bulk of the domain the DNN approximations appear to be in good agreement with the exact13 ∞ periodic DNN C periodic DNN C periodic DNN Exact value u (0) 1.3452550840065e-01 7.7219924514080e-02 6.4853245159018e-02 1 u (2) 1.3452550840063e-01 7.7219924514080e-02 6.4853245159018e-02 -1 u (cid:48) (0) 1.7647534129428e+02 7.0097116541954e+02 1.5259262265715e+03 0 u (cid:48) (2) 1.7647534129428e+02 1.0015322821138e+03 1.5259262265715e+03 0 u (cid:48)(cid:48) (0) -9.2430741908226e+03 8.4509087153852e+05 -4.6596618576551e+05 -2.4674011002723e+00 u (cid:48)(cid:48) (2) -9.2430741908221e+03 -9.4100497028271e+05 -4.5176214142043e+05 2.4674011002723e+00Table 4: Approximation of the non-periodic function u ( x ): Values of the function and its derivatives atthe left/right domain boundaries ( x = 0 and x = 2) from the approximations with C ∞ , C and C periodicDNNs and from the exact function. 14 significant digits are listed to demonstrate that the current methodenforces the periodic conditions exactly for the approximation function and its derivatives.function u ( x ) with all three methods. In a region near the two boundaries, the periodic DNN approximationsexhibit large errors, and one can observe fluctuations in the approximation functions (Gibbs phenomenon).Table 4 lists the values of the approximation function u ( x ) and its derivatives on the two boundaries( x = 0 ,
2) obtained with the C ∞ , C and C periodic DNNs as well as the exact function u ( x ). The resultsagain demonstrate that the current methods enforce exactly the periodicity for the approximation functionand its derivatives. In this subsection we test the performance of the proposed method with the one-dimensional (1D) Helmholtzequation, d udx − λu = f ( x ) , (30)on the domain Ω = { x | a (cid:54) x (cid:54) b } , where λ ( λ (cid:62) a and b are given constants and f ( x ) is a prescribedsource term. We impose periodic boundary conditions (BC) on the domain boundaries, x = a and b .Specifically, we consider two types of periodic boundary conditions. The first type is the C periodiccondition, u ( a ) = u ( b ) , u (cid:48) ( a ) = u (cid:48) ( b ) . (31)The second type is the C ∞ periodic condition, u ( a ) = u ( b ) , u (cid:48) ( a ) = u (cid:48) ( b ) , u (cid:48)(cid:48) ( a ) = u (cid:48)(cid:48) ( b ) , . . . , u ( m ) ( a ) = u ( m ) ( b ) , . . . (32)Note that with the C ∞ periodic condition (32), we are effectively seeking a C ∞ periodic function, with theperiod L = b − a , on the infinite domain x ∈ ( −∞ , ∞ ) that solves the equation (30). Since the Helmholtzequation is a second-order equation, imposing the C periodic condition only, i.e. u ( a ) = u ( b ), does not leadto a unique solution to the problem.For the numerical tests in this section we fix the problem parameters to the following values: λ = 10 , a = 0 , b = 4 , L = b − a = 4 . (33)We choose the source term f ( x ) such that the Helmholtz equation (30) has an analytic solution u ( x ) = sin[3 π ( x + 0 . π ( x + 0 . . (34)This is a periodic function with L = 4 as a period, and it satisfies the boundary conditions (31) and (32).To simulate this problem, we employ a feed-forward neural network (Figure 1) with 4 hidden layers, with20 nodes in each layer, apart from the input and output layers. The input to the network is the coordinate x (1 node), and the output of the network is the solution to the Helmholtz equation u (1 node). The second14 poch l o ss (a) epoch l o ss (b)Figure 7: 1D Helmholtz equation: training histories of the loss function with the (a) C ∞ and (b) C periodic boundary conditions.layer of the network (or the first hidden layer) is set to be a C ∞ periodic layer L p ( m, n ) with m = 11and n = 20, in which we set the constant ω = πL = π in equation (8a), when the C ∞ periodic boundaryconditions in (32) are imposed. When the C periodic periodic boundary conditions in (31) are imposed, weset the second layer of the network to be a C periodic layer L C ( m, n ) with m = 11 and n = 20, as detailedin Section 2.We minimize the following loss function with this DNN,Loss = 1 L (cid:90) Ω (cid:20) d udx − λu − f ( x ) (cid:21) dx = 1 L N e − (cid:88) e =0 (cid:90) Ω e (cid:20) d udx − λu − f ( x ) (cid:21) dx = 1 L N e − (cid:88) e =0 Q − (cid:88) i =0 (cid:34) d udx (cid:12)(cid:12)(cid:12)(cid:12) x ei − λu ( x ei ) − f ( x ei ) (cid:35) J e w i . (35)In this equation, N e is the number of elements (sub-intervals) we have partitioned the domain Ω into inorder to compute the integral, Ω e (0 (cid:54) e (cid:54) N e −
1) denotes the region of element e , Q is the number ofquadrature points within each element, J e is the Jacobian of Ω e with respect to the standard element [ − , x ei (0 (cid:54) i (cid:54) Q −
1) are the Gauss-Lobatto-Legendre quadrature points within element e for 0 (cid:54) e (cid:54) N e − w i (0 (cid:54) i (cid:54) Q −
1) are the weights associated with the Gauss-Lobatto-Legendre quadrature. The inputdata to the network consist of x ei (0 (cid:54) i (cid:54) Q −
1, 0 (cid:54) e (cid:54) N e − f ( x ei )(0 (cid:54) i (cid:54) Q −
1, 0 (cid:54) e (cid:54) N e − u ( x ei ) can be obtained from the output of the DNN,and d udx (cid:12)(cid:12)(cid:12) x ei can be computed by auto-differentiation.For the numerical experiments reported below, we have partitioned the domain into three elements( N e = 3), with these elements being Ω = [0 , . = [1 . , .
6] and Ω = [2 . , . Q = 40) within each element. The activation functions for all the hidden layers are thehyperbolic tangent function (“tanh”), and no activation function is applied to the output layer. The DNNis trained using the Adam optimizer for 5000 epochs with a learning rate 10 − for the C ∞ periodic BC, andfor 15000 epochs with a learning rate 5 × − for the C periodic BC. The options of “early stopping” and“restore to best weight” have been used in Tensorflow/Keras when training the DNN. The training historiesof the loss function for the C ∞ and C periodic BCs are shown in Figure 7. We observe characteristics inthe loss histories similar to those observed in Section 3.1, such as the varied loss reduction rates at differentstages and the fluctuations in the loss value with a fixed learning rate.Figure 8 shows a comparison between the exact solution and the DNN solutions obtained with the C ∞ and C periodic BCs enforced using the current method from Section 2 (left column), as well as the errors15 S o l u t i on s DNNExact (a) x A b s o l u t e e rr o r (b) x S o l u t i on s DNNExact (c) x A b s o l u t e e rr o r (d)Figure 8: 1D Helmholtz equation: (a,c), comparison between the exact and DNN solutions. (b,d), errors ofthe DNN solutions against the exact solution. The results in (a) and (b) are obtained with the C ∞ periodicBCs, an those in (c) and (d) are obtained with the C periodic BCs. Periodic BCs are enforced using themethod from Section 2. DNN C ∞ PBC (Current) DNN C PBC (Current) DNN C PBC (Penalty) Exact solution u (0) 2.4318620484799e+00 2.4323606718327e+00 2.4314502844043e+00 2.4317706231133e+00 u (4) 2.4318620484799e+00 2.4323606718327e+00 2.4317127865026e+00 2.4317706231133e+00 u (cid:48) (0) 7.1133416801319e+00 7.1449674055594e+00 7.1058987676886e+00 7.1050608901229e+00 u (cid:48) (4) 7.1133416801319e+00 7.1449674055594e+00 7.1114177007820e+00 7.1050608901229e+00 u (cid:48)(cid:48) (0) -8.8371916646277e+01 -9.6875622546165e+01 -8.6577693358885e+01 -8.8007775637807e+01 u (cid:48)(cid:48) (4) -8.8371916646277e+01 -8.3340776157218e+01 -8.3307846875545e+01 -8.8007775637806e+01 Table 5: 1D Helmholtz equation: Values of the solution and its first and second derivatives on the left/rightdomain boundaries from the exact solution and from the DNN solutions with C ∞ and C periodic BCsenforced using the current method, and with the C periodic BCs enforced using the penalty method.16 S o l u t i on s DNNExact
Figure 9: 1D Helmholtz equation: DNN solution obtained with the C periodic BC enforced using thepenalty method.of these DNN solutions against the exact solution (right column). The DNN solutions are observed to agreewith the exact solution very well. The DNN solution curves almost exactly overlap with the exact-solutioncurve, and the maximum errors in the domain is on the order of 10 − with both the C ∞ and C periodicboundary conditions.For comparison, we have also computed this problem with the C periodic BCs in another way, byenforcing the C periodic BCs based on the penalty method. Figure 9 shows the DNN solution computedusing the penalty method. Here the DNN has the same parameters (4 hidden layers, with 20 nodes in eachlayer). The C periodic BCs are enforced by including a penalty term in the loss function as follows,Loss = 1 L (cid:90) Ω (cid:20) d udx − λu − f ( x ) (cid:21) dx + θ bc (cid:32) [ u ( a ) − u ( b )] + (cid:20) dudx (cid:12)(cid:12)(cid:12)(cid:12) x = a − dudx (cid:12)(cid:12)(cid:12)(cid:12) x = b (cid:21) (cid:33) , (36)where θ bc = 10 is the penalty coefficient in front of the boundary residual terms. The DNN has been trainedwith the Adam optimizer for 15000 epochs. We observe that the DNN solution resulting from the penaltymethod also agrees well with the exact solution.In Table 5 we list the values of the DNN solution and its first and second derivatives, with 14 significantdigits shown, on the left and right domain boundaries obtained with the C ∞ and C periodic boundaryconditions enforced using the current method, together with those obtained with the C periodic BCsenforced using the penalty method. The boundary values from the exact solution (34) are also included inthe table for comparison. We observe that the current method enforces exactly, to the machine accuracy,the periodicity for the solution and its derivatives (up to order 2 shown here) with the C ∞ periodic BCs.With the C periodic BCs, the current method enforces exactly the periodicity for the solution and its firstderivative, but not for the second derivative. In contrast, the penalty method enforces the periodic conditionfor none of these quantities exactly. We next test the performance of the proposed methods using the the Helmholtz equation in two dimensions(2D), ∂ u∂x + ∂ u∂y − λu = f ( x, y ) , (37)on a rectangular domain Ω = { ( x, y ) | a (cid:54) x (cid:54) b , a (cid:54) y (cid:54) b } . Here λ , a , a , b and b are given constants, u ( x, y ) is the unknown field function to be solved for, and f ( x, y ) is a prescribed source term. We impose17eriodic boundary conditions in both the x and y directions.Specifically, we consider C and C ∞ periodic boundary conditions in 2D. The 2D C periodic BC imposesthe relations: u ( a , y ) = u ( b , y ) , ∂∂x u ( a , y ) = ∂∂x u ( b , y ) , ∀ y ∈ [ a , b ]; u ( x, a ) = u ( x, b ) , ∂∂y u ( x, a ) = ∂∂y u ( x, b ) , ∀ x ∈ [ a , b ] . (38)The 2D C ∞ periodic BC imposes the relations: u ( a , y ) = u ( b , y ) , ∂∂x u ( a , y ) = ∂∂x u ( b , y ) , ∂ ∂x u ( a , y ) = ∂ ∂x u ( b , y ) , . . . , ∀ y ∈ [ a , b ]; u ( x, a ) = u ( x, b ) , ∂∂y u ( x, a ) = ∂∂y u ( x, b ) , ∂ ∂y u ( x, a ) = ∂ ∂y u ( x, b ) , . . . , ∀ x ∈ [ a , b ] . (39)With the C ∞ periodic BC, we are effectively seeking a smooth periodic function u ( x, y ) satisfying u ( x + L , y ) = u ( x, y ) , u ( x, y + L ) = u ( x, y ) , ∀ x, y ∈ ( −∞ , ∞ ) , (40)where L = b − a and L = b − a .We specifically consider the following parameter values for the numerical tests in this section: λ = 10 , a = a = 0 , b = b = 4 , L = L = 4 . (41)We choose the source term f ( x, y ) such that the 2D Helmholtz equation (37) has the solution given by, u ( x, y ) = − [1 . πx + 0 . π ) + 2 cos(2 πx − . π )] [1 . πy + 0 . π ) + 2 cos(2 πy − . π )] . (42)This analytic solution satisfies the periodic boundary conditions (38) and (39).To simulate this problem, we employ a feed-forward DNN with 2 nodes in the input layer, one node inthe output layer, and 4 hidden layers in between. The input layer consists of the coordinates x and y , andthe output layer is the solution to the 2D Helmholtz equation u . Each of the four hidden layers contains 20nodes (neurons) in its output. For the C ∞ periodic BCs, the second layer of this DNN (i.e. the first hiddenlayer) is set to be a 2D C ∞ periodic layer L Dp ( m, n ) with m = 12 and n = 20 (see equation (11)), in whichthe constants ω and ω are set to ω = 2 πL = π , ω = 2 πL = π . (43)For the C periodic BCs, the second layer of this DNN is set to be a 2D C periodic layer L DC ( m, n ) with m = 12 and n = 20 (see equation (22)).We minimize the following loss function,Loss = 1 V Ω (cid:90) Ω (cid:20) ∂ u∂x + ∂ u∂y − λu − f ( x, y ) (cid:21) d Ω= 1 V Ω N e − (cid:88) e =0 Q − (cid:88) i,j =0 (cid:34) ∂ u∂x (cid:12)(cid:12)(cid:12)(cid:12) ( x ei ,y ej ) + ∂ u∂y (cid:12)(cid:12)(cid:12)(cid:12) ( x ei ,y ej ) − λu ( x ei , y ej ) − f ( x ei , y ej ) (cid:35) J e w ij , (44)where V Ω = (cid:82) Ω d Ω = 16 is the area of the domain, N e is the number of elements (sub-domains) we havepartitioned the domain Ω into for computing the integral, Q is the number of quadrature points in the x and y directions within each element, J e is the Jacobian of the element e (0 (cid:54) e (cid:54) N e − x ei , y ej )(0 (cid:54) i, j (cid:54) Q −
1) are the Gauss-Lobatto-Legendre quadrature points within the element e (0 (cid:54) e (cid:54) N e − w ij (0 (cid:54) i, j (cid:54) Q −
1) are the Gauss-Lobatto-Legendre quadrature weights associated with ( x ei , y ej ). Theinput data to the DNN consist of ( x ei , y ej ) (0 (cid:54) i, j (cid:54) Q −
1, 0 (cid:54) e (cid:54) N e − f ( x ei , y ej ) (0 (cid:54) i, j (cid:54) Q − (cid:54) e (cid:54) N e −
1) are passed to the DNN as the label data. In the loss expression, u ( x ei , y ej ) can be obtained fromthe output of the DNN, and the derivatives ∂ u∂x and ∂ u∂y on ( x ei , y ej ) can be computed by auto-differentiation.18 y (a) x y (b) x y (c) x y (d)Figure 10: 2D Helmholtz equation: Contours of DNN Solutions (left column) and their errors againstthe exact solution (right column), obtained with C ∞ (top row) and C (bottom row) periodic boundaryconditions.In the numerical tests below, we partition the domain into 4 elements ( N e = 4), with 2 uniform elementsin both the x and y directions. We use 30 quadrature points ( Q = 30) in each direction within each element.We employ the hyperbolic tangent (“tanh”) function as the activation function for each of the hidden layers,and no activation is applied to the output layer. The DNN is trained using the L-BFGS optimizer for 10500iterations with the C ∞ periodic BCs, and for 6500 iterations with the C periodic BCs.Figure 10 shows contours of the DNN solutions (left column) and their errors (right column) against theexact solution (42), computed with the C ∞ periodic boundary conditions (top row) and the C periodicboundary conditions. The distributions of the DNN solutions are qualitatively the same as that of the exactsolution, and no difference can be discerned visually. The maximum absolute error of the DNN solution inthe domain is less than 5 × − with the C ∞ periodic BC and less than 10 − with the C periodic BC.Figure 11 provides a quantitative comparison between the DNN solutions and the exact solution. It showsthe profiles of the DNN solutions obtained using C ∞ and C periodic BCs, as well as the exact solution,along several horizontal lines across the domain located at y = 0 .
5, 2, and 3 .
5. The error profiles of the DNNsolutions along these lines are also shown in this figure. We observe that the DNN solutions with both the C ∞ and C periodic BCs obtained using the current method agree very well with the exact solution.To examine how well the current methods enforce the periodic boundary conditions for the 2D Helmholtz19 S o l u t i on DNN, C ∞ PBCDNN, C PBCExact (a) x S o l u t i on DNN, C ∞ PBCDNN, C PBCExact (b) x S o l u t i on DNN, C ∞ PBCDNN, C PBCExact (c) x E rr o r s DNN, C ∞ PBCDNN, C PBC (d) x E rr o r s DNN, C ∞ PBCDNN, C PBC (e) x E rr o r s DNN, C ∞ PBCDNN, C PBC (f)Figure 11: 2D Helmholtz equation: comparison of profiles of the solution (top row) and its error againstthe exact solution (bottom row) along several horizontal lines located at (a,d) y = 0 .
5, (b,e) y = 2 .
0, and(c,f) y = 3 .
5, from the DNN solutions with C ∞ and C periodic boundary conditions. The profiles of theexact solution are also included for comparison.equation, we have extracted the values of the DNN solution and its partial derivatives (up to second order) onseveral corresponding points on the left and right boundaries, and on the top and bottom boundaries. Notethat these derivatives are computed by auto-differentiation, and they are the exact derivatives correspondingto the given DNN representation of the field. Table 6 lists the values of the DNN solutions and their partialderivatives on several corresponding points on the left/right boundaries and top/bottom boundaries. Thesecond-order mixed derivatives and some first derivatives are not listed in the table, such as ∂ u∂x∂y , ∂u∂y on theleft/right boundaries, and ∂u∂x on the top/bottom boundaries. These unlisted values are exactly the sameon the corresponding boundary points with both the C ∞ and C periodic BCs. The boxed values in thistable highlight the difference in the second partial derivatives on the corresponding boundary points of theDNN solution obtained with C periodic BCs. As expected, the current methods have enforced the periodicboundary conditions exactly for the solution and the corresponding higher-order derivatives. The next test problem is the unsteady diffusion equation: ∂u∂t − ν ∂ u∂x = f ( x, t ) , (45)where the constant ν > u ( x, t ) is the unknown field function to be solved for, f ( x, t ) is a prescribed source term, x is the spatial coordinate, and t is time. We consider the spatial-temporaldomain Ω = { ( x, t ) | a (cid:54) x (cid:54) b, (cid:54) t (cid:54) T } , where a , b , T are prescribed constants whose values are specifiedbelow. This equation is supplemented by the initial condition, u ( x,
0) = u in ( x ) , (46)where u in denotes the initial distribution. 20NN C ∞ PBC DNN C PBC Exact solution u (0 , .
5) 6.3375720647119e+00 6.3341020999642e+00 6.3375550504603e+00 u (4 , .
5) 6.3375720647119e+00 6.3341020999642e+00 6.3375550504603e+00 u x (0 , .
5) 8.8479796805350e+00 8.8535226659745e+00 8.8433359504634e+00 u x (4 , .
5) 8.8479796805351e+00 8.8535226659745e+00 8.8433359504634e+00 u xx (0 , .
5) -2.0815263069180e+02 -2.1045671256754e+02 -2.0841095826411e+02 u xx (4 , .
5) -2.0815263069180e+02 -2.0946606085833e+02 -2.0841095826411e+02 u (0 , .
5) 3.9958327765965e-01 3.9859455580613e-01 3.9851292703942e-01 u (4 , .
5) 3.9958327765962e-01 3.9859455580613e-01 3.9851292703942e-01 u x (0 , .
5) 5.4769599613000e-01 5.8968112421753e-01 5.5607938177296e-01 u x (4 , .
5) 5.4769599613003e-01 5.8968112421753e-01 5.5607938177296e-01 u xx (0 , .
5) -1.3286736511194e+01 -1.4294530675460e+01 -1.3105126558055e+01 u xx (4 , .
5) -1.3286736511194e+01 -1.2645789026374e+01 -1.3105126558055e+01 u (1 ,
0) -2.4020925534027e+00 -2.4038799915243e+00 -2.4031781074217e+00 u (1 ,
4) -2.40209255340274e+00 -2.4038799915243e+00 -2.40317810742171e+00 u y (1 ,
0) -3.3463550988308e+00 -3.3630705573155e+00 -3.3533612226666e+00 u y (1 ,
4) -3.3463550988308e+00 -3.3630705573155e+00 -3.3533612226666e+00 u yy (1 ,
0) 7.8621899432667e+01 7.9457755849148e+01 7.9028686655862e+01 u yy (1 ,
4) 7.8621899432667e+01 7.9025311915153e+01 7.9028686655862e+01 u (3 ,
0) -2.4018445486575e+00 -2.4022228345948e+00 -2.4031781074217e+00 u (3 ,
4) -2.4018445486575e+00 -2.4022228345948e+00 -2.4031781074217e+00 u y (3 ,
0) -3.3590252565867e+00 -3.3499280675496e+00 -3.3533612226666e+00 u y (3 ,
4) -3.3590252565867e+00 -3.3499280675496e+00 -3.3533612226666e+00 u yy (3 ,
0) 7.8580508305656e+01 7.8911124269480e+01 7.9028686655862e+01 u yy (3 ,
4) 7.85805083056566e+01 7.8910907593259e+01 7.9028686655862e+01Table 6: 2D Helmholtz equation: Values of the solution and its derivatives on selected corresponding pointsof the left/right and top/bottom boundaries, obtained from the DNN solutions with C ∞ and C periodicBCs and from the exact solution. u x = ∂u∂x , u y = ∂u∂y , u xx = ∂ u∂x , and u yy = ∂ u∂y .We impose periodic boundary conditions on the spatial boundaries x = a and b . We specifically considerthe C ∞ and C periodic BCs. The C ∞ periodic BC requires, u ( a, t ) = u ( b, t ) , ∂∂x u ( a, t ) = ∂∂x u ( b, t ) , ∂ ∂x u ( a, t ) = ∂ ∂x u ( b, t ) , . . . , ∀ t ∈ [0 , T ] . (47)The C periodic BC requires, u ( a, t ) = u ( b, t ) , ∂∂x u ( a, t ) = ∂∂x u ( b, t ) , ∀ t ∈ [0 , T ] . (48)For the numerical tests in this subsection we employ the following parameter values, ν = 0 . , a = 0 , b = 4 , T = 4 . , L = b − a = 4 . (49)We choose the source term f ( x, t ) such that the function, u ( x, t ) = (2 cos( πx + 0 . π ) + 1 . πx − . π )) (2 cos( πt + 0 . π ) + 1 . πt − . π )) , (50)is a solution to the equation (45). We choose the initial distribution u in ( x ) by using the analytic expression(50) and setting t = 0. Note that the expression (50) satisfies the boundary conditions (48) and (47). Sounder the initial condition (46) and the periodic boundary conditions, the exact solution to the diffusionequation is given by (50).We solve this initial/boundary value problem using DNN together with the method from Section 2 forenforcing the periodic BCs in x . We employ a DNN with two nodes in the input layer, which represent the21 t (a) x t (b)Figure 12: Diffusion equation: DNN solutions obtained with C ∞ (a) and C (b) periodic boundary condi-tions in the x direction.spatial coordinate x and the time t , and one node in the output layer, which represents the unknown function u ( x, t ) to be solved for. This DNN contains 3 hidden layers between the input and the output layers. Eachof the hidden layers has an output consisting of 30 nodes. Note that periodic BCs are imposed only in the x direction, not in time. For the C ∞ periodic BCs, the second layer of this DNN (or the first hidden layer)is a set to be a modified 2D C ∞ periodic layer as discussed in the Remark 2.2. For the periodic direction x , this modified 2D C ∞ periodic layer corresponds to a 1D C ∞ periodic layer L p ( m, n ) with m = 12 and n = 30 (see equation (8a)), in which the constant ω is set to ω = 2 πL = π . (51)For the C periodic BCs, the second layer of this DNN is set to be a modified 2D C periodic layer asdiscussed in Remark 2.4. For the periodic direction x , this modified 2D C periodic layer corresponds to a1D C periodic layer L C ( m, n ) with m = 12 and n = 30 (see equation (19)).We minimize the following loss function,Loss = θ eq V Ω (cid:90) Ω (cid:20) ∂u∂t − ν ∂ u∂x − f ( x, t ) (cid:21) d Ω + θ ic b − a (cid:90) ba [ u ( x, − u in ( x )] dx = θ eq V Ω N el − (cid:88) e =0 (cid:90) Ω e (cid:20) ∂u∂t − ν ∂ u∂x − f ( x, t ) (cid:21) d Ω + θ ic b − a N xel − (cid:88) e =0 (cid:90) b e a e [ u ( x, − u in ( x )] dx = θ eq ( b − a ) T N el − (cid:88) e =0 Q − (cid:88) i,j (cid:34) ∂u∂t (cid:12)(cid:12)(cid:12)(cid:12) ( x ei ,t ej ) − ν ∂ u∂x (cid:12)(cid:12)(cid:12)(cid:12) ( x ei ,t ej ) − f ( x ei , t ej ) (cid:35) J e w ij + θ ic b − a N xel − (cid:88) e =0 Q − (cid:88) i =0 [ u ( x ei , − u in ( x ei )] J ex w i , (52)where V Ω = (cid:82) Ω d Ω = ( b − a ) T is the volume of the spatial-temporal domain Ω, N el is the number ofspatial-temporal elements we have partitioned Ω into in order to compute the integral, Ω e denotes thespatial-temporal element e for 0 (cid:54) e (cid:54) N el −
1, and Q is the number of Gauss-Lobatto-Legendre quadraturepoints in both the x and t directions within each spatial-temporal element. We use N xel to denote thenumber of elements in the spatial direction, and N tel to denote the number of elements in time, and then N el = N xel N tel . The sub-interval [ a e , b e ] denotes the spatial element e for 0 (cid:54) e (cid:54) N xel −
1. The constants θ eq and θ ic are the penalty coefficients for the equation residual term and the initial condition residual term in(52). The Gauss-Lobatto-Legendre quadrature points within the spatial-temporal element Ω e are denoted22 S o l u t i on s DNN, C ∞ PBCDNN, C PBCExact (a) x S o l u t i on s DNN, C ∞ PBCDNN, C PBCExact (b) x S o l u t i on s DNN, C ∞ PBCDNN, C PBCExact (c) x E rr o r s DNN, C ∞ PBCDNN, C PBC (d) x E rr o r s DNN, C ∞ PBCDNN, C PBC (e) x E rr o r s DNN, C ∞ PBCDNN, C PBC (f)Figure 13: Diffusion equation: comparison of profiles of the solution (top row) and the absolute error(bottom row) at time instants (a,d) t = 0 .
75, (b,e) t = 2 .
25, (c,f) t = 3 .
75 from the DNN solutions with C ∞ and C periodic boundary conditions and from the exact solution.by ( x ei , t ej ), with the associated quadrature weights w ij . J e is the Jacobian associated with the element Ω e for 0 (cid:54) e (cid:54) N el − J ex is the Jacobian of the spatial element [ a e , b e ] for 0 (cid:54) e (cid:54) N xel − w i (0 (cid:54) i (cid:54) Q − x ei in the spacialdirection. The input data to the DNN consist of all the quadrature points within the domain, ( x ei , t ej ), for0 (cid:54) i, j (cid:54) Q − (cid:54) e (cid:54) N el −
1. The values of the source term on the quadrature points, f ( x ei , t ej ),are passed to the DNN as the label data. The terms ∂u∂t (cid:12)(cid:12) ( x ei ,t ej ) and ∂ u∂x (cid:12)(cid:12)(cid:12) ( x ei ,t ej ) in (52) are computed basedon auto-differentiation. It can be observed that here we are enforcing the initial condition by the penaltymethod.In the numerical tests reported below, the domain Ω is partitioned into 9 spatial-temporal elements( N el = 9), with 3 uniform elements in time ( N tel = 3) and also 3 elements in the x direction ( N xel = 3). Alongthe x direction, the two interior boundaries of the elements are located at x = 1 . .
6. We employ 20quadrature points ( Q = 20) in space and time within each spatial-temporal element. The penalty coefficientsare set to be θ eq = 0 . θ ic = 1 − θ eq = 0 .
1. The “tanh” activation function has been employed for eachhidden layer, and no activation function is applied to the output layer of the DNN. The L-BFGS optimizerhas been employed to train the DNN for 13000 iterations with the C ∞ periodic BCs, and for 14500 iterationswith the C periodic BCs.Figure 12 shows contours in the spatial-temporal ( x − t ) plane of the DNN solutions (left column) andtheir errors against the exact solution (right column), obtained using the C ∞ periodic BCs (top row) andthe C periodic BCs (bottom row). Qualitatively, no difference can be discerned between distributions of theDNN solutions and the exact solution given by (50). Figure 13 provides a quantitative comparison betweenthe DNN solutions and the exact solution. It shows profiles of the exact solution and the DNN solutionscorresponding to C ∞ and C periodic BCs at three time instants t = 0 .
75, 2 .
25 and 3 .
75. The profilesof the absolute errors of the DNN solutions are shown in this figure as well. The DNN solution profilescorresponding to both the C ∞ and C periodic BCs are observed to overlap with those of the exact solutionat different time instants. It can be observed from the error profiles that the DNN with the C ∞ periodicBCs appears to result in generally smaller errors than that with the C periodic BCs for this problem.23NN C ∞ PBC DNN C PBC Exact solution u (0 , .
6) -2.4017654769459e+00 -2.4116317389834e+00 -2.4031781074217e+00 u (4 , .
6) -2.4017654769459e+00 -2.4116317389834e+00 -2.4031781074217e+00 u x (0 , .
6) -1.0998296722163e+01 -1.0634488873369e+01 -1.0970511273439e+01 u x (4 , .
6) -1.0998296722163e+01 -1.0634488873369e+01 -1.0970511273439e+01 u xx (0 , .
6) -4.6844809230852e+00 -9.0735064404064e+00 -4.8498203327098e+00 u xx (4 , .
6) -4.6844809230855e+00 3.4743060900784e+00 -4.8498203327099e+00 u (0 , .
92) 8.8799799206032e-01 8.8941798966940e-01 8.8446899507452e-01 u (4 , .
92) 8.8799799206032e-01 8.8941798966940e-01 8.8446899507452e-01 u x (0 , .
92) 4.0457999629891e+00 3.8734589665624e+00 4.0376021450541e+00 u x (4 , .
92) 4.0457999629891e+00 3.8734589665624e+00 4.0376021450541e+00 u xx (0 , .
92) 1.1406604203965e+00 5.2000604027409e+00 1.7849345842143e+00 u xx (4 , .
92) 1.1406604203965e+00 7.0346713619307e-01 1.7849345842144e+00 u (0 , .
3) 1.7326144730460e+00 1.7391187048218e+00 1.7317627457812e+00 u (4 , .
3) 1.7326144730459e+00 1.7391187048218e+00 1.7317627457812e+00 u x (0 , .
3) 7.9202902481403e+00 7.7216199212532e+00 7.9054992498656e+00 u x (4 , .
3) 7.9202902481403e+00 7.7216199212532e+00 7.9054992498656e+00 u xx (0 , .
3) 3.2478290437743e+00 6.6397883264174e+00 3.4948463245322e+00 u xx (4 , .
3) 3.2478290437744e+00 2.3875667090972e-01 3.4948463245323e+00Table 7: Diffusion equation: values of the solution and its derivatives on the spatial boundaries ( x = 0 and4) at several time instants (14 significant digits shown), obtained from the DNN solutions with C ∞ and C periodic BCs and from the exact solution. The boxes highlight the differences in the obtained values for thesecond derivatives.To assess how well the current method enforces the periodic boundary conditions, we have extracted thevalues of the solution and its partial derivatives (up to order two) on the spatial boundaries ( x = 0 and4) at several time instants from the exact solution and the DNN solutions obtained using the C ∞ and C periodic boundary conditions. Table 7 lists these boundary values for the solution and its derivatives, with14 significant digits shown. As expected, the current method has enforced exactly the periodicity for thesolution and all its derivatives extracted here with the C ∞ periodic BCs, while with the C periodic BCs itenforces exactly the periodicity only for the solution and its first derivative. In the last numerical example, we test the proposed method using the wave equation, ∂u∂t − c ∂u∂x = 0 , (53)where the prescribed constant c represents the wave speed, u ( x, t ) is the unknown field function to be solvedfor, x is the spatial coordinate and t is the time. We consider the spatial-temporal domain Ω = { ( x, t ) | a (cid:54) x (cid:54) b, (cid:54) t (cid:54) T } for this problem, where a , b and T are prescribed constants whose values are to bespecified below. We consider the following initial condition, u ( x,
0) = u in ( x ) = 2 sech 3( x − x ) δ , ∀ x ∈ [ a, b ] , (54)where x ∈ [ a, b ] and δ are prescribed constants whose values are specified below.We impose the periodic boundary condition on the spatial boundaries of the domain, x = a and b .Specifically, we consider C ∞ , C , C and C periodic boundary conditions in this test. The C k ( k = 0 , , u ( a, t ) = u ( b, t ) , ∀ t ∈ [0 , T ]; (55) ∂∂x u ( a, t ) = ∂∂x u ( b, t ) , ∀ t ∈ [0 , T ]; (56) ∂ ∂x u ( a, t ) = ∂ ∂x u ( b, t ) , ∀ t ∈ [0 , T ] . (57)The C periodic BC imposes the condition (55). The C periodic BC imposes the conditions (55) and (56).The C periodic BC imposes the conditions (55)–(57). The C ∞ periodic BC imposes the conditions: u ( a, t ) = u ( b, t ) , ∂∂x u ( a, t ) = ∂∂x u ( b, t ) , . . . , ∂ m ∂x m u ( a, t ) = ∂ m ∂x m u ( b, t ) , . . . , ∀ t ∈ [0 , T ] . (58)This initial/boundary value problem has the solution, u ( x, t ) = (cid:26) u in ( x + ct ) = 2 sech x − x + ct ) δ , if ( x + ct ) ∈ [ a, b ] ,u ( x ± L, t ) , otherwise , ∀ ( x, t ) ∈ Ω , (59)where L = b − a . In the numerical tests reported below we have employed the following values for theparameters: c = 2 , T = 4 , a = 0 , b = 4 , L = b − a = 4 , δ = 1 , x = 2 . (60)To solve this initial/boundary value problem, we employ a feed-forward DNN together with the methodfrom Section 2 for enforcing the periodic boundary conditions. The input layer of the DNN consists of twonodes, which represent the spatial coordinate x and the time t . The output layer of the DNN consists ofone node, which represents the field function u to be solved for. We employ 3 hidden layers between theinput and the output layers. Each hidden layer has an output with 30 nodes. Since the periodic BC isonly imposed in the x direction, we employ the modified 2D periodic layers to enforce periodic boundaryconditions; see the Remarks 2.2 and 2.4. For the C ∞ periodic BCs, the second layer of this DNN (or the firsthidden layer) is set to be a modified 2D C ∞ periodic layer as discussed in Remark 2.2. For the x direction,this modified periodic layer corresponds to the 1D C ∞ periodic layer L p ( m, n ) with m = 12 and n = 30 (seeequation (8)), in which the constant ω is set to ω = πL = π . For the C k ( k = 0 , ,
2) periodic BCs, thesecond layer of this DNN is set to be a modified 2D C k periodic layer, which for the x direction correspondsto the 1D C k periodic layer L C k ( m, n ) with m = 12 and n = 30 (see equation (19)).We employ the following loss function for this DNN,Loss = θ eq V Ω (cid:90) Ω (cid:18) ∂u∂t − c ∂u∂x (cid:19) d Ω + θ ic b − a (cid:90) ba [ u ( x, − u in ( x )] dx = θ eq V Ω N el − (cid:88) e =0 (cid:90) Ω e (cid:18) ∂u∂t − c ∂u∂x (cid:19) d Ω + θ ic b − a N xel − (cid:88) e =0 (cid:90) b e a e [ u ( x, − u in ( x )] dx = θ eq ( b − a ) T N el − (cid:88) e =0 Q − (cid:88) i,j =0 (cid:20) ∂∂t u ( x ei , t ej ) − c ∂∂x u ( x ei , t ej ) (cid:21) J e w ij + θ ic b − a N xel − (cid:88) e =0 Q − (cid:88) i =0 [ u ( x ei , − u in ( x ei )] J ex w i . (61)In the above expression, V Ω = (cid:82) Ω d Ω = ( b − a ) T is the volume of the spatial-temporal domain Ω, and theconstants θ eq and θ ic are the penalty coefficients for the loss terms associated with the equation and theinitial condition, respectively. In order to compute the integrals, we have partitioned the domain Ω into N el spatial-temporal elements, with N xel elements in the spatial direction and N tel elements in time, leading to therelation N el = N xel N tel . Ω e denotes the region occupied by the spatial-temporal element e for 0 (cid:54) e (cid:54) N el − a e , b e ] denotes the region of the spatial element e for 0 (cid:54) e (cid:54) N xel − Q is the number of25 t (a) x t (b) x t (c) x t (d)Figure 14: Wave equation: Contours of the DNN solutions obtained with the C ∞ (a), C (b), C (c), C (d) periodic boundary conditions in the x direction.quadrature points in both the spatial and temporal directions within each spatial-temporal element. ( x ei , t ej )(0 (cid:54) i, j (cid:54) Q −
1) are the Gauss-Lobatto-Legendre quadrature points within the spatial-temporal elementΩ e , for 0 (cid:54) e (cid:54) N el − J e is the Jacobian of the spatial-temporal element Ω e (0 (cid:54) e (cid:54) N el − J ex is the Jacobian of the spatial element [ a e , b e ] (0 (cid:54) e (cid:54) N xel − w ij (0 (cid:54) i, j (cid:54) Q −
1) denote theweights associated with the Gauss-Lobatto-Legendre quadrature points ( x ei , t ej ). w i (0 (cid:54) i (cid:54) Q −
1) denotethe weights associated with the spatial Gauss-Lobatto-Legendre quadrature point x ei . The input data to theDNN are the quadrature points ( x ei , t ej ) for 0 (cid:54) i, j (cid:54) Q − (cid:54) e (cid:54) N el −
1. In the loss function (61), u ( x ei , t ej ) is obtained from the output of the DNN, and the terms ∂u∂t (cid:12)(cid:12) ( x ei ,t ej ) and ∂u∂x (cid:12)(cid:12) ( x ei ,t ej ) can be computedbased on auto-differentiation. It can be observed that the initial condition is enforced by the penalty method.For the numerical results reported below, we have partitioned the domain Ω into 4 spatial-temporalelements ( N el = 4), with 2 uniform elements along the spatial and temporal directions ( N xel = N tel = 2).We employ 30 quadrature points in both space and time ( Q = 30) within each spatial-temporal element.The penalty coefficients are set to be θ eq = 0 . θ ic = 1 − θ eq = 0 .
1. We use “tanh” as the activationfunctions for the hidden layers. No activation is applied to the output layer. The Adam optimizer has beenemployed to train the DNN for 60 ,
000 epochs with the C ∞ and C periodic BCs, for 90 ,
000 epochs withthe C periodic BC, and for 80 ,
000 epochs with the C periodic BCs. The options for “early stopping” and“restore to best weight” in Tensorflow/Keras are employed during the training of the DNNs.Figure 14 illustrates the distributions of the DNN solutions and their absolute errors. The plots in theleft column of this figure are contours of the DNN solutions in the spatial-temporal ( x − t ) plane, obtainedwith the C ∞ and C k ( k = 0 , ,
2) periodic BCs on the spatial boundaries. The plots in the right column arecontours of the absolute errors of these solutions against the exact solution (59). Qualitatively, no differencecan be discerned of the distributions between the DNN solutions and the exact solution. The maximumerrors in the domain are approximately on the order of magnitude 10 − .Figure 15 shows a temporal sequence of snapshots of the wave form, obtained from the DNN solutionwith the C ∞ periodic boundary conditions. One can clearly observe the propagation of the wave form in26 S o l u t i on (a) x S o l u t i on (d) x S o l u t i on (g) x S o l u t i on (b) x S o l u t i on (e) x S o l u t i on (h) x S o l u t i on (c) x S o l u t i on (f) x S o l u t i on (i)Figure 15: Wave equation: Snapshots of the wave profiles obtained with the C ∞ periodic BC at timeinstants (a) t = 0, (b) t = 0 .
5, (c) t = 1 .
0, (d) t = 1 .
5, (e) t = 2 .
0, (f) t = 2 .
5, (g) t = 3 .
0, (h) t = 3 .
5, (i) t = 4 . − x direction at a constant speed. Because of the imposed periodic conditions, as soon as the wave exitsthe left boundary ( x = 0), it re-enters the domain through the right boundary ( x = 4) in a seamless andsmooth fashion.In Figure 16 we compare the wave profiles from the exact solution (59) and from the DNN solutions withdifferent types of periodic BCs at three time instants ( t = 0 .
5, 2 and 3 . t = 0 .
5, 2 and 3 . C ∞ and C periodic BCs, the current method has enforced exactly the periodic conditions forthe solution and its first and second derivatives. With the C periodic BC, the method enforces exactly theperiodic condition only for the solution. With the C periodic BC, the method enforces exactly the periodiccondition for the solution and its first derivative, but not for its second derivative.27 S o l u t i on s DNN, C ∞ PBCDNN, C PBCDNN, C PBCDNN, C PBCExact (a) x S o l u t i on s DNN, C ∞ PBCDNN, C PBCDNN, C PBCDNN, C PBCExact (b) x S o l u t i on s DNN, C ∞ PBCDNN, C PBCDNN, C PBCDNN, C PBCExact (c) x E rr o r s DNN, C ∞ PBCDNN, C PBCDNN, C PBCDNN, C PBC (d) x E rr o r s DNN, C ∞ PBCDNN, C PBCDNN, C PBCDNN, C PBC (e) x E rr o r s DNN, C ∞ PBCDNN, C PBCDNN, C PBCDNN, C PBC (f)Figure 16: Wave equation: comparison of profiles of the solution (top row) and the absolute error (bottomrow) obtained with the C ∞ , C , C and C periodic boundary conditions and with the exact solution atseveral time instants: (a, d) t = 0 .
5, (b, e) t = 2 .
0, (c, f) t = 3 . DNN C ∞ PBC DNN C PBC DNN C PBC DNN C PBC Exact solution u (0 , .
5) 2.0142799261625e-01 2.0190756733086e-01 2.0109869263559e-01 2.0083996838890e-01 1.9865585483886e-01 u (4 , .
5) 2.0142799261625e-01 2.0190756733086e-01 2.0109869263559e-01 2.0083996838891e-01 1.9865585483886e-01 u x (0 , .
5) 5.9482715294468e-01 5.7394660776802e-01 5.9323728164392e-01 5.8746219179701e-01 5.9302035811534e-01 u x (4 , .
5) 5.9482715294468e-01 5.8248575768435e-01 5.9323728164392e-01 5.8746219179701e-01 5.9302035811534e-01 u xx (0 , .
5) 1.6980636297410e+00 1.4696979666932e+00 1.5994875423905e+00 1.6679798466138e+00 1.7526236647042e+00 u xx (4 , .
5) 1.6980636297410e+00 1.9395212811018e+00 1.5277304396975e+00 1.6679798466138e+00 1.7526236647042e+00 u (0 ,
2) 9.8284051583224e-03 9.9241674367022e-03 9.9440310013471e-03 1.1696938690667e-02 9.9149477871207e-03 u (4 ,
2) 9.8284051583222e-03 9.9241674367022e-03 9.9440310013471e-03 1.1696938690667e-02 9.9149477871207e-03 u x (0 ,
2) 1.5256678688244e-03 -4.8521014846574e-03 9.9336770382177e-04 2.0456109624086e-03 2.9744477846340e-02 u x (4 ,
2) 1.5256678688247e-03 2.3815674708337e-03 9.9336770382177e-04 2.0456109624087e-03 2.9744477846340e-02 u xx (0 ,
2) 2.4222108105432e-01 4.2985279066538e-01 1.1927388753390e-01 2.1278569005620e-01 8.9230143930769e-02 u xx (4 ,
2) 2.4222108105432e-01 2.2731343637548e-01 2.6553423293645e-01 2.1278569005620e-01 8.9230143930769e-02 u (0 , .
5) 2.0170117943758e-01 2.0200326210422e-01 2.0458115518013e-01 1.9932351260077e-01 1.9865585483886e-01 u (4 , .
5) 2.0170117943758e-01 2.0200326210422e-01 2.0458115518013e-01 1.9932351260077e-01 1.9865585483886e-01 u x (0 , .
5) -5.8538598646821e-01 -6.1113064140580e-01 -5.8981323269116e-01 -5.7592303278591e-01 -5.9302035811534e-01 u x (4 , .
5) -5.8538598646821e-01 -5.9639323869798e-01 -5.8981323269116e-01 -5.7592303278591e-01 -5.9302035811534e-01 u xx (0 , .
5) 1.6523878185422e+00 2.3592894577992e+00 1.6262664810111e+00 1.6631571855496e+00 1.7526236647042e+00 u xx (4 , .
5) 1.6523878185422e+00 1.2305332426305e+00 1.8162774394782e+00 1.6631571855496e+00 1.7526236647042e+00
Table 8: Wave equation: values of the solution and its derivatives on the left/right boundaries at severaltime instants ( t = 0 .
5, 2 and 3 . C ∞ , C , C and C periodicboundary conditions and from the exact solution. 28 Concluding Remarks
In this paper we have presented a method for enforcing exactly the C ∞ and C k (for any k (cid:62)
0) peri-odic conditions with deep neural networks. The method stems from some simple properties about functioncompositions involving periodic functions. The method essentially composes an arbitrary DNN-representedfunction with a set of independent known periodic functions with adjustable (training) parameters. Morespecifically, we have defined the operations that constitute a C ∞ periodic layer and a C k periodic layer.The DNN with a C ∞ periodic layer incorporated as the second layer of the network (behind the input)automatically and exactly satisfies the C ∞ periodic boundary conditions in its output. The DNN with a C k periodic layer incorporated as the second layer automatically and exactly satisfies the C k periodic boundaryconditions in its output. The C ∞ periodic layer comprises constructions of a set of independent C ∞ periodicfunctions with a prescribed period, based on sinusoidal functions, affine mappings, and nonlinear activationfunctions. The C k periodic layer comprises constructions of a set of independent C k periodic functions,based on the generalized Hermite interpolation polynomials, affine mappings, and nonlinear activation func-tions. We have tested the method in extensive numerical experiments with ordinary and partial differentialequations involving C ∞ and C k periodic boundary conditions. The numerical results demonstrate that theproposed method indeed enforces exactly, to the machine accuracy, the periodicity for the solution and itsderivatives.The proposed method can be implemented in a straightforward way. The C ∞ and C k periodic layersdefined herein can be implemented as user-defined Keras layers, and used in the same way as the built-incore Keras layers. All the numerical examples in the current work are implemented based on Tensorflow andKeras.Periodic functions and periodic boundary conditions have widespread applications in computational sci-ence of various disciplines. The proposed method provides an effective tool, based on deep neural networks,for representing periodic functions and enforcing exactly the periodic boundary conditions. We anticipatethat this method will be instrumental in expanding DNN-based techniques to new classes of applicationsthat are unexplored or scarcely explored so far.An outstanding question concerning the method developed herein is the following: Can an arbitraryperiodic function of a certain regularity be represented by the current periodic DNNs to arbitrary accu-racy? This is an important question and it is open at this point. Our extensive numerical experimentsseem to suggest that the answer to this question is positive. The periodic DNNs from the current methodare essentially compositions of an arbitrary DNN-represented function with a set of independent periodicfunctions with adjustable parameters. Can the universal approximation power of the original DNN carryover to the resultant periodic DNN as the set of independent periodic functions becomes sufficiently large?Can theoretical analysis establish an analogous universal approximation property for such periodic DNNswith respect to periodic functions? These are interesting questions that call for future research and shouldbe pursued by the community. Acknowledgement
This work was partially supported by NSF (DMS-1522537).
References [1] J. Berg and K. Nystrom. A unified deep artifial neural network approach to partial differential equationsin complex geometries.
Neurocomputing , 317:28–41, 2018.[2] C. Canuto, M.Y. Hussaini, A. Quarteroni, and T. Zang.
Spectral methods in fluid dynamics . Springer,New York, 1988.[3] J. Chen, R. Du, and K. Wu. A comprehensive study of boundary conditions when solving pdes by dnns. arXiv:2005.04554 , 2020. 294] N.E. Cotter. The stone-weierstrass theorem and its application to neural networks.
IEEE Transactionson Neural Networks , 4:290–295, 1990.[5] S. Dong. Direct numerical simulation of turbulent Taylor-Couette flow.
J. Fluid Mech. , 587:373–393,2007.[6] S. Dong. Turbulent flow between counter-rotating concentric cylinders: a direct numerical simulationstudy.
Journal of Fluid Mechanics , 615:371–399, 2008.[7] S. Dong. A convective-like energy-stable open boundary condition for simulations of incompressibleflows.
Journal of Computational Physics , 302:300–328, 2015.[8] S. Dong, G.E. Karniadakis, and C. Chryssostomidis. A robust and accurate outflow boundary conditionfor incompressible flow simulations on severely-truncated unbounded domains.
Journal of ComputationalPhysics , 261:83–105, 2014.[9] S. Dong, G.E. Karniadakis, A. Ekmekci, and D. Rockwell. A combined direct numerical simulation-particle image velocimetry study of the turbulent near wake.
J. Fluid Mech. , 569:185–207, 2006.[10] S. Dong and X. Zheng. Direct numerical simulation of spiral turbulence.
J. Fluid Mech. , 668:150–173,2011.[11] W. E and B. Yu. The deep Ritz method: a deep learning-based numerical algorithm for solvingvariational problems.
Communications in Mathematics and Statistics , 6:1–12, 2018.[12] A.R. Gallant and H. White. There exists a neural network that does not make avoidable mistakes.
Proceedings of the Second Annual IEEE Conference on Neural Networks , 1988.[13] F.S. Gokuzum, L.T.K. Nguyen, and M.-A. Keip. An artificial neural network based solution schemefor periodic computational homogenization of electrostatic problems.
Mathematical and ComputationalApplications , 24:40, 2019.[14] I. Goodfellow, Y. Bengio, and A. Courville.
Deep Learning . The MIT Press, 2016.[15] P.M. Gresho. Incompressible fluid dynamics: some fundamental formulation issues.
Annual Review ofFluid Mechanics , 23:413–453, 1991.[16] K. Hornik, M. Stinchcombe, and H. White. Multilayer feedforward networks are universal approxima-tors.
Neural Networks , 2:359–366, 1989.[17] K. Hornik, M. Stinchcombe, and H. White. Universal approximation of an unknown mapping and itsderivatives using multilayer feedforward networks.
Neural Networks , 3:551–560, 1990.[18] D.P. Kingma and J. Ba. Adam: a method for stochastic optimization. arXiv:1412.6980 , 2014.[19] I.E. Lagaris, A.C. Likas, and D.I. Fotiadis. Artificial neural networks for solving ordinary and partialdifferential equations.
IEEE Transactions on Neural Networks , 9:987–1000, 1998.[20] I.E. Lagaris, A.C. Likas, and D.G. Papageorgiou. Neural-network methods for boundary value problemswith irregular boundaries.
IEEE Transactions on Neural Networks , 11:1041–1049, 2000.[21] X. Li. Simultaneous approximations of mulvariate functions and their derivatives by neural networkswith one hidden layer.
Neurocomputiing , 12:327–343, 1996.[22] S. Liu. Fourier neural network for machine learning.
Proceedings of the International Conference onMachine Learning and Cybernetics (ICMLC) , 2013.[23] K.S. McFall and J.R. Mahan. Artifical neural network method for solution of boundary value problemswith exact satisfaction of arbitrary boundary conditions.
IEEE Transactions on Neural Networks ,20:1221–1233, 2009. 3024] M. Ngom and O. Marin. Approximating periodic functions and solving differential equations using anovel type of fourier neural networks. arXiv:2005.13100 , 2020.[25] N. Ni, Z. Yang, and S. Dong. Energy-stable boundary conditions based on a quadratic form: Applica-tions to outflow/open-boundary problems in incompressible flows.
Journal of Computational Physics ,391:179–215, 2019.[26] J. Nocedal and S.J. Wright.
Numerical Optimization, Second Edition . Springer, 2006.[27] M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: a deep learningframework for solving forward and inverse problems involving nonlinear partial differential equations.
Journal of Computational Physics , 378:686–707, 2019.[28] C. Rao, H. Sun, and Y. Liu. Physics informed deep learning for computational elastodynamics withoutlabeled data. arXiv:2006.0872 , 2020.[29] K. Rudd and S. Ferrari. A constrained integration (CINT) approach to solving partial differentialequations using artificial neural networks.
Neurocomputing , 155:277–285, 2015.[30] E. Samanaiego, C. Anitescu, S. Goswami, V.M. Nguyen-Thanh, H. Guo, K. Hamdia, X. Zhuang, andT. Rabczuk. An energy approach to the solution of partial differential equations in computationalmechanics via machine learning: concepts, implementation and applications.
Computer Methods inApplied Mechanics and Engineering , 362:112790, 2020.[31] R.L. Sani and P.M. Gresho. Resume and remarks on the open boundary conidtion minisymposium.
International Journal for Numerical Methods in Fluids , 18:983–1008, 1994.[32] A. Silvscu. Fourier neural networks.
Proceedings of the International Conference on Neural Networks(IJCNN99) , pages 488–491, 1999.[33] J. Sirignano and K. Spoliopoulos. DGM: A deep learning algorithm for solving partial differentialequations.
Journal of Computational Physics , 375:1339–1364, 2018.[34] A. Spitzbart. A generalization of Hermite’s interpolation formula.
The American Mathematical Monthly ,67:42–46, 1960.[35] J.F. Traub. On lagrange-hermite interpolation.
J. Soc. Indust. Appl. Math. , 12:886–891, 1964.[36] Y. Zang, G. Bao, X. Ye, and H. Zhou. Weak adversarial networks for high-dimensional partial differentialequations.
Journal of Computational Physics , 411:109409, 2020.[37] A. Zhumekenov, M. Uteuliyeva, R. Takhanov, Z. Assylbekov, A.J. Castro, and O. Kabdolov. Fourierneural networks: a comparative study. arXiv:1902.03011arXiv:1902.03011