Year: 2020

Abstract:
Approximating the numerical solutions of partial differential equations (PDEs) using neural networks is a promising application of deep learning. The smooth architecture of a fully connected neural network is appropriate for finding the solutions of PDEs; the corresponding loss function can also be intuitively designed and guarantees the convergence for various kinds of PDEs. However, the rate of convergence has been considered as a weakness of this approach. This paper introduces a novel loss function for the training of neural networks to find the solutions of PDEs, making the training substantially efficient. Inspired by the recent studies that incorporate derivative information for the training of neural networks, we develop a loss function that guides a neural network to reduce the error in the corresponding Sobolev space. Surprisingly, a simple modification of the loss function can make the training process similar to Sobolev Training although solving PDEs with neural networks is not a fully supervised learning task. We provide several theoretical justifications for such an approach for the viscous Burgers equation and the kinetic Fokker–Planck equation. We also present several simulation results, which show that compared with the traditional L loss function, the proposed loss function guides the neural network to a significantly faster convergence. Moreover, we provide the empirical evidence that shows that the proposed loss function, together with the iterative sampling techniques, performs better in solving high dimensional PDEs.

heading: 1 INTRODUCTION
text: Deep learning has achieved remarkable success in many scientific fields, including computer vision and natural language processing. In addition to engineering, deep learning has been successfully applied to the field of scientific computing. Particularly, the use of neural networks for the numerical integration of partial differential equations (PDEs) has emerged as a new important application of the deep learning.
Being a universal approximator (Cybenko, 1989; Hornik et al., 1989; Li, 1996), a neural network can approximate solutions of complex PDEs. To find the neural network solution of a PDE, a neural network is trained on a domain wherein the PDE is defined. Training a neural network comprises the following: feeding the input data through forward pass and minimizing a predefined loss function with respect to the network parameters through backward pass. In the traditional supervised learning setting, the loss function is designed to guide the neural network to generate the same output as the target data for the given input data. However, while solving PDEs using neural networks, the target values that correspond to the analytic solution are not available. One possible way to guide the neural network to produce the same output as the solution of the PDE is to penalize the neural network to satisfy the PDE itself (Sirignano & Spiliopoulos, 2018; Berg & Nyström, 2018; Raissi et al., 2019; Hwang et al., 2020).
Unlike the traditional mesh-based schemes including the finite difference method (FDM) and the finite element method (FEM), neural networks are inherently mesh-free function-approximators. Advantageously, as mesh-free function-approximators, neural networks can avoid the curse of dimensionality (Sirignano & Spiliopoulos, 2018) and approximate the solutions of PDEs on complex geometries (Berg & Nyström, 2018). Recently, Hwang et al. (2020) showed that neural networks could approximate the solutions of kinetic Fokker–Planck equations under not only various kinds
of kinetic boundary conditions but also several irregular initial conditions. Moreover, they showed that the neural networks automatically approximate the macroscopic physical quantities including the kinetic energy, the entropy, the free energy, and the asymptotic behavior of the solutions. Further issues including the inverse problem were investigated by Raissi et al. (2019); Jo et al. (2020).
Although the neural network approach can be used to solve several complex PDEs in various kinds of settings, it requires relatively high computational cost compared to the traditional mesh-based schemes in general. To resolve this issue, we propose a novel loss function using Sobolev norms in this paper. Inspired by a recent study that incorporated derivative information for the training of neural networks (Czarnecki et al., 2017), we develop a loss function that efficiently guides neural networks to find the solutions of PDEs. We prove that the H1 and H2 norms of the approximation errors converge to zero as our loss functions tend to zero for the 1-D Heat equation, the 1-D viscous Burgers equation, and the 1-D kinetic Fokker–Planck equation. Moreover, we show via several simulation results that the number of epochs to achieve a certain accuracy is significantly reduced as the order of derivatives in the loss function gets higher, provided that the solution is smooth. This study might pave the way for overcoming the issue of high computational cost when solving PDEs using neural networks.
The main contributions of this work are threefold: 1) We introduce novel loss functions that enable the Sobolev Training of neural networks for solving PDEs. 2) We prove that the proposed loss functions guarantee the convergence of neuarl networks in the corresponding Sobolev spaces although it is not a supervised learning task. 3) We empirically demonstrate the effect of Sobolev Training for several regression problems and the improved performances of our loss functions in solving several PDEs including the heat equation, Burgers’ equation, the Fokker–Planck equation, and high-dimensional Poisson equation.

heading: 2 RELATED WORKS
text: Training neural networks to approximate the solutions of PDEs has been intensively studied over the past decades. For example, Lagaris et al. (1998; 2000) used neural networks to solve Ordinary Differential Equations (ODEs) and PDEs on a predefined set of grid points. Subsequently, Sirignano & Spiliopoulos (2018) proposed a method to solve high-dimensional PDEs by approximating the solution using a neural network. They focused on the fact that the traditional finite mesh-based scheme becomes computationally intractable when the dimension becomes high. However, because neural networks are mesh-free function-approximators, they can solve high-dimensional PDEs by incorporating mini-batch sampling. Furthermore, the authors showed the convergence of the neural network to the solution of quasilinear parabolic PDEs under certain conditions.
Recently, Raissi et al. (2019) reported that one can use observed data to solve PDEs using physicsinformed neural networks (PINNs). Notably, PINNs can solve a supervised regression problem on observed data while satisfying any physical properties given by nonlinear PDEs. A significant advantage of PINNs is that the data-driven discovery of PDEs, also called the inverse problem, is possible with a small change in the code. The authors provided several numerical simulations for various types of nonlinear PDEs including the Navier–Stokes equation and Burgers’ equation. The first theoretical justification for PINNs was provided by Shin et al. (2020), who showed that a sequence of neural networks converges to the solutions of linear elliptic and parabolic PDEs in L2 sense as the number of observed data increases. There also exists a study aiming to enhance the convergence of PINNs (van der Meer et al., 2020).
Additionally, several works related deep neural networks with PDEs but not by the direct approximation of the solutions of PDEs. For instance, Long et al. (2018) attempted to discover the hidden physics model from data by learning differential operators. A fast, iterative PDE-solver was proposed by learning to modify each iteration of the existing solver (Hsieh et al., 2019). A deep backward stochastic differential equation (BSDE) solver was proposed and investigated in Weinan et al. (2017); Han et al. (2018) for solving high-dimensional parabolic PDEs by reformulating them using BSDE.
The main strategy of the present study is to leverage derivative information while solving PDEs via neural networks. The authors of Czarnecki et al. (2017) first proposed Sobolev Training that uses derivative information of the target function when training a neural network by slightly modifying
the loss function. They showed that Sobolev Training had lower sample complexity than regular training, and therefore it is highly efficient in many applicable fields, such as regression and policy distillation problems. We appropriate the concept of Sobolev Training to develop a loss function for the efficient training of a neural network for solving PDEs.

heading: 3 LOSS FUNCTION
text: We consider the following Cauchy problem of PDEs:
Pu = f, (t, x) ∈ [0, T ]× Ω, (3.1) Iu = g, (t, x) ∈ {0} × Ω, (3.2) Bu = h, (t, x) ∈ [0, T ]× ∂Ω, (3.3)
where P denotes a differential operator; I and B denote the initial and boundary operators, respectively; f , g, and h denote the inhomogeneous term, and initial and boundary data, respectively. In most studies that reported the neural network solutions of PDEs, a neural network was trained on uniformly sampled grid points {(ti, xj)}Nt,Nxi,j=1 ∈ [0, T ] × Ω, which were completely determined before training. One of the most intuitive ways to make the neural network satisfy PDEs (3.1)–(3.3) is to minimize the following loss functional:
Loss(unn; p) = ‖Punn − f‖pLp([0,T ]×Ω) + ‖Iunn − g‖ p Lp(Ω) + ‖Bunn − h‖ p Lp([0,T ]×∂Ω),
where unn denotes the neural network and p = 1 or 2, as they have been the most commonly used exponents in regression problems in previous studies. Evidently, an analytic solution u satisfies Loss(u) = 0, and thus one can conceptualize a neural network that makes Loss(unn) = 0 a possible solution of PDEs (3.1)–(3.3). This statement is in fact proved for second-order parabolic equations with the Dirichlet boundary condition in Jo et al. (2020), and for the Fokker–Planck equation with inflow and specular reflective boundary conditions in Hwang et al. (2020). Both the proofs are based on the following inequality:
‖u− unn‖L∞(0,T ;L2(Ω)) ≤ CLoss(unn; 2), for some constant C, which states that minimizing the loss functional implies minimizing the approximation error.
The main concept behind Sobolev Training is to minimize the error between the output and the target function, and that between the derivatives of the output and those of the target function. However, unlike the traditional supervised regression problem, neither the target function nor its derivative is provided while solving PDEs via neural networks. Thus, a special treatment is required to apply Sobolev Training for solving PDEs using neural networks. In this and the following sections, we propose several loss functions and prove that they guarantee the convergence of the neural network to the solution of a given PDE in the corresponding Sobolev space. Therefore, the proposed loss functions play similar roles to those in Sobolev Training.
We define the loss function that depends on the Sobolev norm W k,p as follows: LossGE(unn; k, p, l, q) = ∥∥∥‖P (unn(t, ·))− f(t, ·)‖qW l,q(Ω)∥∥∥pWk,p([0,T ]) , (3.4)
LossIC(unn; l, q) = ‖Iunn(t, x)− g(x)‖qW l,q(Ω), (3.5) LossBC(unn; k, p, l, q) = ∥∥∥‖Bunn(t, ·)− h(t, ·)‖qW l,q(∂Ω)∥∥∥pWk,p([0,T ]) . (3.6)
Remark 3.1. Here, Loss(0)TOTAL(unn) = LossGE(unn; 0, 2, 0, 2) + LossIC(unn; 0, 2) + LossBC(unn; 0, 2, 0, 2) coincides with the traditional L2 loss function employed by Sirignano & Spiliopoulos (2018); Berg & Nyström (2018); Raissi et al. (2019); Hwang et al. (2020).
When we train a neural network, the loss functions (3.4)–(3.6) are computed by Monte-Carlo approximation. Because the grid points are uniformly sampled, the loss functions are approximated as follows:
LossGE(unn; k, p, l, q) ≈ T |Ω| NtNx ∑ |β|≤k Nt∑ i=1 ∣∣∣∣∣∣ d β dtβ ∑ |α|≤l Nx∑ j=1 |DαP (unn(ti, xj))−Dαf(ti, xj)|q ∣∣∣∣∣∣ p ,
LossIC(unn; l, q) ≈ |Ω| Nx ∑ |α|≤l Nx∑ j=1 |Dαunn(0, xj)−Dαg(xj)|q,
LossBC(unn; k, p, l, q) ≈ T |∂Ω| NtNB ∑ |β|≤k Nt∑ i=1 ∣∣∣∣∣∣ d β dtβ ∑ |α|≤l ∑ xj∈∂Ω |Dαunn(ti, xj)−Dαh(ti, xj)|q ∣∣∣∣∣∣ p ,
where α and β denote the conventional multi-indexes, and D denotes the spatial derivatives.

heading: 4 THEORETICAL RESULTS
text: In this section, we theoretically validate our claim that our loss functions guarantee the convergence of the neural network to the solution of a given PDE in the corresponding Sobolev spaces, and that they play a similar role to those in Sobolev Training while solving PDEs via neural networks. Throughout this section, we will denote the strong solution of each equation by u, neural network solution by unn, and Sobolev spaces W 1,2 and W 2,2 by H1 and H2, respectively. All the proofs are provided in the Appendix.

heading: 4.1 THE HEAT EQUATION AND BURGERS’ EQUATION
text: We define the following three total loss functions for the heat equation and Burgers’ equation:
Loss (0) TOTAL(unn) = LossGE(unn; 0, 2, 0, 2) + LossIC(unn; 0, 2) + LossBC(unn; 0, 2, 0, 2),
(4.1)
Loss (1) TOTAL(unn) = LossGE(unn; 0, 2, 0, 2) + LossIC(unn; 1, 2) + LossBC(unn; 0, 2, 0, 2),
(4.2)
Loss (2) TOTAL(unn) = LossGE(unn; 1, 2, 0, 2) + LossIC(unn; 2, 2) + LossBC(unn; 0, 2, 0, 2).
(4.3)
We then obtain the following convergence theorem:
Theorem 4.1. (Proofs are provided in (A.5) for the heat equation, and (A.8) for Burgers’ equation) For the following 1-D heat and Burgers’ equations:
The heat equation Burgers’ equation ut − uxx = 0 in (0, T ]× Ω, u(0, x) = u0(x) on Ω, u(t, x) = 0 on [0, T ]× ∂Ω, ut + uux − νuxx = 0 in (0, T ]× Ω, u(0, x) = u0(x) on Ω, u(t, x) = 0 on [0, T ]× ∂Ω,
there hold, provided that unn is smooth,
max 0≤t≤T
‖u(t)− unn(t)‖L2(Ω) → 0 as Loss (0) TOTAL → 0,
ess sup 0≤t≤T
‖u(t)− unn(t)‖H10 (Ω) → 0 as Loss (1) TOTAL → 0,
ess sup 0≤t≤T
‖u(t)− unn(t)‖H2(Ω) → 0 as Loss (2) TOTAL → 0.

heading: 4.2 THE FOKKER–PLANCK EQUATION
text: For the Fokker–Planck equation, we need additional parameters for a new input variable v. We define the following two total loss functions for the Fokker–Planck equation:
Loss (0;FP) TOTAL(unn) = LossGE(unn; 0, 2, 0, 2, 0, 2) + LossIC(unn; 0, 2, 0, 2)
+ LossBC(unn; 0, 2, 0, 2, 0, 2), (4.4)
Loss (1;FP) TOTAL(unn) = LossGE(unn; 0, 2, 1, 2, 1, 2) + LossIC(unn; 1, 2, 1, 2)
+ LossBC(unn; 0, 2, 0, 2, 0, 2). (4.5)
We then have the following convergence theorem: Theorem 4.2. (Proofs are provided in (A.10) and (A.12)) For the 1-D Fokker–Planck equation with the periodic boundary condition:
ut + vux − β(vu)v − quvv = 0, for (t, x, v) ∈ [0, T ]× [0, 1]× R, u(0, x, v) = u0(x, v), for (x, v) ∈ [0, 1]× R,
∂αt,x,vu(t, 1, v)− ∂αt,x,vu(t, 0, v) = 0, for (t, v) ∈ [0, T ]× R,
there hold, under assumptions (A.50) and (A.51),
sup 0≤t≤T
‖u(t)− unn(t)‖L2(Ω×[−V,V ]) → 0 as Loss (0;FP) TOTAL → 0,
sup 0≤t≤T
‖u(t)− unn(t)‖H1(Ω;L2([−V,V ])) → 0 as Loss (1;FP) TOTAL → 0.
Remark 4.3. The theorems in this section imply that the proposed loss functions guarantee the convergence of neural networks in the corresponding Sobolev spaces, thereby coinciding with the main idea of Sobolev Training.
Remark 4.4. The theorems in this section cannot be directly generalized to the high-dimensional cases because even the 2-dimensional case starts involving the convexity of the boundary. Though it has also been shown that the Fokker-Planck operator has strong hypoellipticity and the solutions to the boundary problems are smooth even in the higher dimensional case, the proof requires long rigorous mathematical analysis. For more information, see Hwang et al. (2018; 2019).
Remark 4.5. Because we cannot access the label (which corresponds to the analytic solution) on the interior grid, solving PDEs using a neural network is not a fully supervised problem. Interestingly, by incorporating derivative information in the loss function, the proposed approach enables Sobolev Training even if neither the labels nor the derivatives of the target function are provided.

heading: 5 EXPERIMENTAL RESULTS
text: In this section, we provide experimental results for toy examples that comprise several regression problems and various kinds of differential equations, including the heat equation, Burgers’ equation, the kinetic Fokker–Planck equation, and high-dimensional Poisson’s equation. We employ a fully connected neural network, which is a natural choice for function approximation. We use the hyperbolic tangent function as a nonlinear activation function. Although ReLU(x) = max(0, x) is a frequent choice in modern machine learning, it is not appropriate for solving PDEs because the second derivatives of the neural network vanish.
In appreciation of Automatic Differentiation, we can easily compute derivatives of any order of a neural network with respect to input data despite the compositional structure; see Baydin et al. (2017) and references therein. We implemented our neural network using PyTorch, a widely used deep learning library (Paszke et al., 2019). For the numerical experiments, we used a neural network with three hidden layers each of which had d-256-256-256-1 neurons, where d denotes the input dimension. We used the ADAM optimizer (Kingma & Ba, 2014), a popular gradient-based optimizer.
To see whether our loss functions performed more efficiently than the traditional L2 loss function introduced in Remark 3.1, we made everything maintain the same except the loss function. We compared the loss functions on the basis of L2(Ω) test error for the toy examples, absolute relative test error for the high-dimensional Poisson equation, and L∞(0, T ;L2(Ω)) test error for the other PDEs. For each loss function, we recorded the number of epochs required to meet a certain error threshold and the test error. Considering the randomness due to network initialization, we repeated the training a hundred times. Conversely, we initialized a hundred different neural networks with uniform initialization and trained them in the same manner. To compute the test error, we used analytic solutions for the Heat equation, Burgers’ equation, and the high-dimensional Poisson equation, and a numerical solution from Wollman & Ozizmir (2008) for the kinetic Fokker–Planck equation.

heading: 5.1 TOY EXAMPLES
text: First, we consider two simple regression problems with target functions sin(x) and ReLU(x), respectively. For these toy examples, we define the loss functions as follows:
L2 loss = ‖unn(x)− y(x)‖22, H1 loss = ‖unn(x)− y(x)‖22 + ‖u′nn(x)− y′(x)‖22, H2 loss = ‖unn(x)− y(x)‖22 + ‖u′nn(x)− y′(x)‖22 + ‖u′′nn(x)− y′′(x)‖22,
where y(x) denotes either sin(x), or ReLU(x). We uniformly sampled a hundred grid points from [0, 2π] for training sin(x). Similarly, we uniformly sampled a hundred grid points from [−1, 1] for training ReLU(x). We expected the training to become fast using higher order derivatives as many as possible when training sin(x) and ReLU(x). Figure 1 confirms our assumption to be true. Interestingly, although ReLU(x) is not twice weakly differentiable at only one point x = 0, the H2 loss does not facilitate the training.
In order to explore the nature of Sobolev Training, we design more complicated toy examples. Consider the target functions sin(kx), for k = 1, 2, ..., 5, and ReLU(kx) = max(0, kx), for k = 1, 2, 3..., 10. As k increases, the target functions and their derivatives contain drastic changes in their values, so it is difficult to learn those functions. We hypothesize that in Sobolev Training, the training becomes faster since we give explicit label for the derivatives and it becomes easier to capture the drastic changes in the derivatives. This is empirically shown to be true in Figure 2. We train neural networks to approximate sin(kx), andReLU(kx) for different k and record the number of training epochs to achieve certain error threshold which can be regarded as a difficulty of the problem. As one can see in Figure 2, the difficulty changes little to no when we train with H1 and H2 losses while the difficulty increases with k when L2 loss is used. This implies that the difficulty of training barely changes in Sobolev Training even the target function has stiff changes. The same observations are made when solving PDEs. The improvement of our loss functions compare to L2 loss function are more dramatic for Burgers’ equation (which has stiff solution (Raissi et al., 2019)) than for the heat equation, with the initial condition of f2 (which has a higher frequency) than with the initial condition of f1 initial condition in the Fokker–Planck equation, and as k increases for the high-dimensional Poisson equation 7.

heading: 5.2 THE HEAT EQUATION & BURGERS’ EQUATION
text: We now demonstrate the results of the Sobolev Training of the neural networks for solving PDEs. We begin with the 1-D heat equation, and Burgers’ equation, which is the simplest PDE that combines both the nonlinear propagation effect and diffusive effect. Burgers’ equation often appears as a
simplification of a more complex and sophisticated model, such as the Navier–Stokes equation. The equations with the homogeneous Dirichlet boundary condition read as follows:
The Heat equation Burgers’ equation ut − uxx = 0 in (0, 10]× [0, π], u(0, x) = sin(x) on [0, π], u(t, x) = 0 on [0, 10]× {0, π}. ut + uux − 0.2uxx = 0 in (0, 0.01]× [0, 1], u(0, x) = − sin(πx) on [0, 1], u(t, x) = 0 on [0, 0.01]× {0, 1}.
The heat equation attains a unique analytic solution u(t, x) = sin(x) exp(−t); an analytic solution of Burgers’ equation is provided in Basdevant et al. (1986).
Although Sirignano & Spiliopoulos (2018) indicated that iterative random sampling reduces the computational cost, we fixed the grid points before training because we aimed to compare the efficiency of our loss function with that of the traditional one. For the heat equation and the Burgers’ equation, we uniformly sampled the grid points {ti, xj}Nt,Nxi,j=1 from (0, T ] × Ω, where Nt and Nx denote the number of samples for interior t and x, respectively. For the initial and boundary conditions, we sampled the grid points from {t = 0, xj}Nxj=1 ∈ {0}×Ω and {ti, xj} Nt,NB i,j=1 ∈ [0, T ]×∂Ω, respectively, where NB denotes the number of grid points in ∂Ω. Here, we set Nt, Nx, NB = 31. The testing data were also uniformly sampled from the domain of the PDEs.
The L2, H1, and H2 losses are the Monte-Carlo approximations of (4.1), (4.2), and (4.3), respectively, for the heat equation and Burgers’ equation. Working on achieving a smooth solution, we observed that the H2 loss performed the best, followed by the H1 loss and then the L2 loss in both accuracy, and computation time. We show the corresponding results in Figure 3.

heading: 5.3 THE FOKKER–PLANCK EQUATION
text: The kinetic Fokker–Planck equation describes the dynamics of a particle whose behavior is similar to that of the Brownian particle. The Fokker–Planck operator has a strong regularizing effect not just in the velocity variable but also in the temporal and the spatial variables by the hypoellipticity. The Fokker–Planck equation has been considered in numerous physical circumstances including the Brownian motion described by the Uhlenbeck-Ornstein processes.
We provide two simulation results for different initial conditions for the 1-D Fokker–Planck equation with the periodic boundary condition. For the Fokker–Planck equation, we adopted the idea of sampling from Hwang et al. (2020). Because it is practically difficult to consider the entire space for the v ∈ R variable, we truncated the space for v as [−5, 5]. We then uniformly sampled the grid points {ti, xj , vk}Nt,Nx,Nvi,j,k=1 from (0, T ]×Ω× [−5, 5], where Nv denotes the number of samples for v. The grid points for the initial and periodic boundary conditions were accordingly sampled. The truncated equation reads as follows:
ut + vux − β(vu)v − quvv = 0, for (t, x, v) ∈ (0, 3]× [0, 1]× [−5, 5], u(0, x, v) = f(x, v), for (x, v) ∈ [0, 1]× [−5, 5],
∂αt,x,vu(t, 1, v)− ∂αt,x,vu(t, 0, v) = 0, for (t, v) ∈ [0, 3]× [−5, 5],
where f(x, v) is either f1(x, v) = exp(−v2)∫ 5
−5 exp(−v2)dv , or f2(x, v) = (1+cos(2πx)) exp(−v2)∫ 1 0 ∫ 5 −5(1+cos(2πx)) exp(−v2)dvdx ,
and β = 0.1, q = 0.1.
A numerical solution on the test data was computed by a method shown by Wollman & Ozizmir (2008) and used for computing the test error. L2 loss and H1 loss denote the Monte-Carlo approximations of (4.4) and (4.5), respectively. The values of Nt, Nx, and Nv were set to be 31, and the grid points were uniformly sampled. Expectedly, a solution of the Fokker–Planck equation could be estimated substantially faster using our loss function in both cases. We have provided the detailed results in Figure 4.

heading: 5.4 THE HIGH-DIMENSIONAL POISSON EQUATION
text: The Poisson equation serves as an example problem in the recent literature; see Weinan & Yu (2018); Hsieh et al. (2019); Zang et al. (2020). In this section, we provide empirical results to demonstrate
that the proposed loss functions perform satisfactorily when equipped with iterative sampling for solving high-dimensional PDEs; see Sirignano & Spiliopoulos (2018) for more information. Convergence result similar to those of in section 4 for the Poisson equation is given in section A.4. We consider the following high-dimensional Poisson equation with the Dirichlet boundary condition:
−4u = π 2
4 d∑ i=1 sin( π 2 xi), for x ∈ Ω = (0, 1)d,
u = d∑ i=1 sin( π 2 xi), for x ∈ ∂Ω,
where x = (x1, x2, ..., xd) ∈ Ω. One can readily prove that u(x) = ∑d i=1 sin( π 2xi) is a strong solution. We compare the following three loss functions with each other:
Loss (0;Poisson) TOTAL (unn) = LossGE(unn; 0, 2) + LossBC(unn; 0, 2), (5.1)
Loss (1;Poisson) TOTAL (unn) = LossGE(unn; 1, 2) + LossBC(unn; 0, 2), (5.2)
Loss (2;Poisson) TOTAL (unn) = LossGE(unn; 1, 2) + LossBC(unn; 1, 2). (5.3)
Notably, the aforementioned loss functions have the variable x only. Table 1 presents the relative errors on a predefined test set for d = 10, 50, and 100. Evidently, in all cases, the proposed loss functions outperform the traditional L2 loss function. More detailed experimental results are given in section B.

heading: 6 DISCUSSION AND CONCLUSION
text: Inspired by Sobolev Training, we proposed novel loss functions, which efficiently guided the training of neural networks for solving PDEs. We theoretically justified that the proposed loss functions guaranteed the convergence of a neural network to a solution of PDEs in the corresponding Sobolev spaces. We also discussed that the proposed theorems imply that the training becomes Sobolev Training by slightly modifying the loss function, although the process of estimating neural network solutions of PDEs is not fully supervised.
In addition to the toy examples, which showed the exceptional speed of Sobolev Training, we provided empirical evidences demonstrate that our loss functions expedited the training more than the traditional L2 loss function. We believe that this can solve the problem associated with the high costs involved in estimating the neural network solutions of PDEs. Moreover, our experiments on high-dimensional problems showed that the proposed loss function performed better when equipped with iterative grid sampling. The histograms in Figure 1-4 indicate that our loss function provided more stable training in that it reduced the variance in the distribution of the number of epochs (e.g., for the Burgers’ equation, L2 loss: 3651±812, H1 loss: 995±71, and H2 loss: 331±15). Thus, the training, when governed by our loss function, became robust to the random initialization of the weights.

heading: A PROOFS FOR THE THEOREMS IN SECTION 4
text: We begin with the basic definitions of Sobolev spaces. We excerpt the definitions from Evans (2010). Definition A.1. Suppose u, v ∈ L1loc(U) and α is a multiindex. We say that v is the αth-weak derivative of u, written
Dαu = v,
provided ∫ U uDαφdx = (−1)|α| ∫ U vφdx, for all test functions φ ∈ C∞c (U).
Now we define the Sobolev space. Fix 1 ≤ p ≤ ∞ and let k be a nonnegative integer. Definition A.2. The Sobolev space W k,p(U) consists of all locally integrable functions u : U → R such that for each multiindex α with |α| ≤ k, Dα exists in the weak sense and belongs to Lp(U). Note that if p = 2, we usually write Hk(U) = W k,2(U) (k = 0, 1, 2, ...).
The Sobolev norms are defined as follow: Definition A.3. If u ∈W k,p(U), we define its norm to be
‖u‖Wk,p(U) =
{ ( ∑ |α|≤k ∫ U |Dαu|pdx)1/p (1 ≤ p ≤ ∞),∑
|α|≤k ess supU |Dαu| (p =∞).
Finally, we define a notion of convergence in Sobolev spaces. Definition A.4. Let {um}∞m=1, u ∈W k,p(U). We say um converges to u in W k,p(U) provided
lim m→∞
‖um − u‖Wk,p(U) = 0.

heading: A.1 THE HEAT EQUATION
text: We denote the strong solution of the heat equation
ut − uxx = 0 in (0, T ]× Ω, u(0, x) = u0(x) on Ω, u(t, x) = 0 on [0, T ]× ∂Ω,
by u and the neural network solution by unn. Then, v = u− unn satisfies: vt − vxx = f(t, x) in (0, T ]× Ω, v(0, x) = g(x) on Ω, (A.1) v(t, x) = 0 on [0, T ]× ∂Ω,
for some f, and g. Here, we can set the boundary to be zero by multiplying B(x), where B(x) is a smooth function satisfying B(x) = {
0, x ∈ ∂Ω 6= 0, x ∈ Ω . Then the following holds:
Theorem A.5. (Theorem 7.1.5 in Evans (2010)) If g ∈ H2(Ω), ft ∈ L2(0, T ;L2(Ω)), then, max
0≤t≤T ‖v(t)‖L2(Ω) ≤ C1(‖f‖L2(0,T ;L2(Ω)) + ‖g‖L2(Ω)), (A.2)
ess sup 0≤t≤T ‖v(t)‖H10 (Ω) ≤ C2(‖f‖L2(0,T ;L2(Ω)) + ‖g‖H10 (Ω)), (A.3)
ess sup 0≤t≤T
‖v(t)‖H2(Ω) ≤ C3(‖f‖H1(0,T ;L2(Ω)) + ‖g‖H2(Ω)), (A.4)
for some C1, C2, C3.
By applying above theorem to (A.1), we get the results of Theorem 4.1. Remark A.6. The left hand sides in (A.2) - (A.4) are the errors of neural networks in corresponding norms, and the right hand sides are the losses (4.1) - (4.3) for the heat equation, respectively. This implies that the proposed loss functions are the upper bounds of the errors in the Sobolev spaces, and by minimizing them, we can expect the effect of Sobolev Training when solving PDEs with neural networks.
In the rest of this section, we will show the similar results for Burgers’ equation and the Fokker– Planck equation.

heading: A.2 BURGERS’ EQUATION
text: We consider the strong solution u of the following Burgers equation in a bounded interval Ω = [a, b],
∂tu+ u∂xu− ∂2xu = 0 in Ω, (A.5) u = 0 on ∂Ω, (A.6)
and the corresponding neural network solution unn satisfying
∂tunn + unn∂xunn − ∂2xunn = f in Ω, (A.7) unn = 0 on ∂Ω. (A.8)
with the inital data u(0, ·) and unn(0, ·), respectively. The following proposition ensures the existence of a strong solution to the initial boundary value problem (A.5)–(A.6) (see Benia & Sadallah (2016)). Here, we multiply B(x) to unn(t, x) in order to meet the boundary condition. We use the notation . where the relation A . B stands for A ≤ CB, where C denotes a generic constant. Proposition A.7. (Theorem 1.2 in Benia & Sadallah (2016)) Let u0 ∈ H10 . Then there exists a time T ∗ = T ∗(u0) > 0 such that the problem (A.5)–(A.6) with initial data u0 has a unique solution of u satisfying
u ∈ L2(0, T ∗;H2(Ω)) ∩ C([0, T ∗);H10 (Ω)), ut ∈ L2(0, T ∗;L2(Ω)).
Furthermore, if T ∗ <∞, then ‖u‖H1(Ω) →∞ as t→ T ∗.
We will show that the following theorem holds.
Theorem A.8. Let u and unn be strong solutions of (A.5)–(A.6) and (A.7)–(A.8) respectively, on the time interval [0, T ]. For w := unn − u, following statements are valid.
(1) There exists a continuous function
F0 = F0 ( ‖w(0, ·)‖22, ∫ T 0 ‖f‖22dt, ∫ T 0 ‖∂xu‖22dt ) such that
sup 0≤t≤T
‖w‖22 + ∫ T
0
‖∂xw‖22dt ≤ F0 → 0,
as
‖w(0, ·)‖22, ∫ T
0
‖f‖22dt→ 0.
(2) There exists a continuous function
F1 = F1 ( ‖w(0, ·)‖2H1 , ∫ T 0 ‖f‖22dt, ∫ T 0 ‖∂xu‖22dt ) such that
sup 0≤t≤T
‖w‖2H1 + ∫ T
0
‖∂xw‖2H1 + ‖∂tw‖22dt ≤ F1 → 0, (A.9)
as
‖w(0, ·)‖2H1 , ∫ T
0
‖f‖22dt→ 0.
(3) There exists a continuous function
F2 = F2 ( ‖w(0, ·)‖2H2 , ∫ T 0 ‖f‖22 + ‖∂tf‖22dt, sup 0≤t≤T ‖∂xu‖22 + ∫ T 0 ‖∂tu‖22dt ) such that
sup 0≤t≤T
( ‖w‖2H2 + ‖wt‖22 ) + ∫ T 0 ‖∂xw‖2H2dt ≤ F2 → 0, (A.10)
as
‖w(0, ·)‖2H2 , ∫ T
0
‖f‖22 + ‖∂tf‖22dt→ 0.
Remark A.9. By the Morrey’s embedding theorem and the Poincare’s inequality, for f ∈ H10 (Ω), we have the following inequality,
‖f‖2∞ . ‖f‖22 + ‖fx‖22 . ‖fx‖22, (A.11)
Throughout the proof, we widely use (A.11).
Proof. Subtracting (A.5) from (A.7), we get equations of w as follows.
wt − wxx + wwx + wux + uwx = f in Ω, (A.12) w = 0 on ∂Ω, (A.13)
w(0, ·) = g in Ω. (A.14)
By multiplying w to (A.12) and integrating by parts in Ω, we have
1
2
d dt ‖w‖22 + ‖wx‖22 (A.15)
= ∫ Ω fw − ∫ Ω w2wx − ∫ Ω w2ux − ∫ Ω uwwx + ∫ ∂Ω wwx,
= 5∑ k=1 I0k .
Now we estimate the terms on the right hand side of (A.15). Applying the Young’s inequality, the Hölder’s inequality, the Sobolev inequality, and thd Poincare inequality, we have
I01 ≤ ‖f‖2‖w‖2 . ‖f‖22 + ‖w‖22 . ‖f‖22 + ‖wx‖22, (A.16)
I02 = ∫ ∂Ω 1 3 w3 = 1 3 w3(b)− 1 3 w3(a) = 0, (A.17) I03 ≤ ‖w‖∞‖w‖2‖ux‖2 . ‖wx‖2‖w‖2‖ux‖2, (A.18) . ‖ux‖22‖w‖22 + ‖wx‖22,
I04 ≤ ‖u‖∞‖w‖2‖wx‖2 . ‖ux‖22‖w‖22 + ‖wx‖22, (A.19) I05 = w(b)wx(b)− w(a)wx(a) = 0, (A.20)
for any small > 0. Applying estimates (A.16)–(A.20) to (A.15), we have the following inequality
d dt ‖w‖22 + ‖wx‖22 . ‖ux‖22‖w‖22 + ‖f‖22, (A.21)
(A.21) and the Grönwall inequality imply that
sup 0≤t≤T
‖w‖22 . e ∫ T 0 ‖ux‖22dt ( ‖g‖22 + ∫ T 0 ‖f‖22dt ) . (A.22)
Let us denote ∫ T 0 ‖ux‖22dt = β
Now we integrate (A.21) between 0 and T and drop the term ‖w(T )‖22 on the left-hand side to obtain∫ T 0 ‖wx‖2dt . ‖g‖22 + ∫ T 0 ‖ux‖22‖w‖22 + ‖f‖22dt (A.23)
. (βeβ + 1) ( ‖g‖22 + ∫ T 0 ‖f‖22dt ) .
This completes the proof of (1) of Theorem A.8.
Next, by multiplying −wxx to (A.12) and integrating by parts in Ω, we obtain
1
2
d dt ‖wx‖22 + ‖wxx‖22 (A.24)
= − ∫
Ω
fwxx + ∫ Ω wwxwxx + ∫ Ω uxwwxx + ∫ Ω uwxwxx + ∫ ∂Ω wtwx,
= 5∑ k=1 I1k .
Similarly to (A.16)–(A.20), we estimate the terms on the right hand side of (A.24).
I11 ≤ ‖f‖2‖wxx‖22 . ‖f‖22 + ‖wxx‖22, (A.25) I12 ≤ ‖w‖∞‖wx‖2‖wxx‖2 . ‖wx‖4 + ‖wxx‖22, (A.26) I13 ≤ ‖w‖∞‖ux‖2‖wxx‖2 . ‖ux‖22‖wx‖22 + ‖wxx‖22, (A.27) I14 ≤ ‖u‖∞‖wx‖2‖wxx‖2 . ‖ux‖22‖wx‖22 + ‖wxx‖22, (A.28) I15 = wt(b)wx(b)− wt(a)wx(a) = 0. (A.29)
for any small > 0. Applying estimates (A.25)–(A.29) to (A.24), we have the following inequality
d dt ‖wx‖22 + ‖wxx‖22 . (‖wx‖22 + ‖ux‖22)‖wx‖22 + ‖f‖22. (A.30)
It follows from (A.23), (A.30), and the Grönwall inequality that
sup 0≤t≤T
‖wx‖22 . e ∫ T 0 ‖wx‖22+‖ux‖ 2 2dt ( ‖gx‖22 + ∫ T 0 ‖f‖22dt ) (A.31)
. eβeF0 ( ‖gx‖22 + ∫ T 0 ‖f‖22dt ) .
In a similar way to (A.23), there exists a function F̃1 = F1(‖f‖L2(0,T ;L2), ‖g‖H1 , β) such that∫ T 0 ‖wxx‖22dt . F1. (A.32)
(2) of Theorem A.8 follows from (A.31), (A.32), and the fact that
wt = f + wxx − wwx − wux − uwx. (A.33)
Finally, we differentiate (A.12) with respect to t, then we obtain
wtt − wxxt + wtwx + wwxt + wtux + wuxt + utwx + uwxt = ft in Ω, (A.34) wt = 0 on ∂Ω, (A.35)
wt(0) = f + gxx − ggx − gu0x − u0gx in Ω. (A.36)
By multiplying wt to (A.34) and integrating by parts in Ω, we have
1
2
d dt ‖wt‖22 + ‖wxt‖22 (A.37)
= ∫ Ω ftwt − ∫ Ω w2twx − ∫ Ω wwtwxt − ∫ Ω w2t ux
− ∫
Ω
wwtuxt − ∫
Ω
utwtwx − ∫
Ω
uwtwxt + ∫ ∂Ω wtwxt,
= 8∑ k=1 I2k .
Terms on the right hand side of (A.37) are estimated by
I21 ≤ ‖ft‖2‖wt‖2 . ‖ft‖22 + ‖wxt‖22, (A.38) I22 ≤ ‖wt‖∞‖wt‖2‖wx‖2 . ‖wx‖22‖wt‖22 + ‖wxt‖22, (A.39) I23 ≤ ‖w‖∞‖wt‖2‖wxt‖2 . ‖wx‖22‖wt‖22 + ‖wxt‖22, (A.40) I24 ≤ ‖wt‖∞‖wt‖2‖ux‖2 . ‖ux‖22‖wt‖22 + ‖wxt‖22, (A.41)
I25 + I 2 6 = ∫ Ω wwxtut ≤ ‖w‖∞‖ut‖2‖wxt‖2 . ‖ut‖22‖wt‖22 + ‖wxt‖22, (A.42)
I27 ≤ ‖u‖2‖wt‖2‖wxt‖2 . ‖ux‖22‖wt‖22 + ‖wxt‖22, (A.43) I28 = 0. (A.44)
for any small > 0. Applying estimates (A.38)–(A.44) to (A.37), we have the following inequality
d dt ‖wt‖22 + ‖wxt‖22 .
( ‖wx‖22 + ‖ux‖22 + ‖ut‖22 ) ‖wt‖22 + ‖ft‖22. (A.45)
It follows from (A.36), (A.45) and the Grönwall inequality that
sup 0≤t≤T
‖wt‖22 . e ∫ T 0 ‖wx‖22+‖ux‖ 2 2+‖ut‖ 2 2dt ( ‖wt(0)‖22 + ∫ T 0 ‖ft‖22 ) (A.46)
. eγeF0 ( ‖f0‖22 + ‖gxx‖22 + ‖gx‖42 + ‖u0x‖22‖gx‖22 + ∫ T 0 ‖ft‖22dt ) .
where
γ = ∫ T 0 ‖ux‖22 + ‖ut‖22dt.
In a similar way to the proof of (2) of Theorem A.8, (3) of Theorem A.8 follows from (A.33) and (A.46). This completes the proof of the Theorem.

heading: A.3 THE FOKKER–PLANCK EQUATION

heading: A.3.1 BOUNDARY LOSS DESIGN
text: Define the loss function for the periodic boundary condition as
LossBC = ∑ |α|=1 ∫ T 0 dt ∫ 5 −5 dv ∣∣∂αt,x,vfnn(t, 1, v;m,w, b)− ∂αt,x,vfnn(t, 0, v;m,w, b)∣∣2
≈ 1 Ni,k ∑ |α|=1,i,k ∣∣∂αt,x,vfnn(ti, 1, vk;m,w, b)− ∂αt,x,vfnn(ti, 0, vk;m,w, b)∣∣2 . (A.47)

heading: A.3.2 THE FOKKER–PLANCK EQUATION IN A PERIODIC INTERVAL
text: In this section, we introduce an L2 energy method for the Fokker–Planck equation and introduce a regularity inequality for the solutions to the equation. Throughout the section, we will abuse the notation and use both notations ∂zu and uz for the same derivative of u with respect to z.
We consider the Fokker–Planck equation in a periodic interval [0, 1]:
ut + vux − β(vu)v − quvv = 0, for (t, x, v) ∈ [0, T ]× [0, 1]× R, u(0, x, v) = u0(x, v), for (x, v) ∈ [0, 1]× R, and
∂αt,x,vu(t, 1, v)− ∂αt,x,vu(t, 0, v) = 0, for (t, v) ∈ [0, T ]× R, (A.48)
for any 3-dimensional multi-index α such that |α| ≤ 1 and a given initial distribution u0 = u0(x, v). Now we consider the Fokker–Planck equation that the corresponding neural network solution unn would satisfy:
(unn)t + v(unn)x − β(vunn)v − q(unn)vv = f for (t, x, v) ∈ [0, T ]× [0, 1]× [−5, 5], unn(0, x, v) = g, for (x, v) ∈ [0, 1]× [−5, 5],∑
|α|=1
∫ T 0 dt ∫ 5 −5 dv ∣∣∣∣(∂αt,x,vunn)(t, 1, v)− (∂αt,x,vunn)(t, 0, v)∣∣∣∣2 ≤ L, (A.49) for any 3-dimensional multi-index α such that |α| ≤ 1 and given f = f(t, x, v), g = g(x, v), and a constant L > 0. Suppose that f , g and h are C1 functions. Also, we suppose that the a priori solutions u and unn are sufficiently smooth; indeed, we require them to be in C 1,1,2 t,x,v .
For the a priori solution u and unn to equation A.48 and equation A.49, assume that if |v| is sufficiently large, then we have that for some sufficiently small > 0,
sup t∈[0,T ] ∥∥∂αt,x,vu(t, ·,±5)− ∂αt,x,vunn(t, ·,±5)∥∥L2x([0,1]) ≤ , (A.50) for |α| ≤ 1 and α = (0, 0, 2). Also, suppose that
|∂αt,x,vu(t, x,±5)|, |∂αt,x,vunn(t, x,±5)| ≤ C, (A.51)
for some C < ∞ for |α| ≤ 1 and α = (0, 0, 2). Now we introduce the following theorem on the energy estimates:
Theorem A.10. Let u and unn be the classical solutions to equation A.48 and equation A.49, respectively. Then we have
sup 0≤t≤T
‖unn(t)− u(t)‖22 + 2(q − ε) ∫ T
0
‖∂vunn(s)− ∂vu(s)‖22ds
≤ ( ‖g − u0‖22 + L
2
) exp [( 1 + 25β2
2ε
) T ] + ∫ T 0 ‖f(s)‖22ds+ 2q CT,
for any ε ∈ (0, q), where L, u0, f, g, β, q,m, , and C are given in equation A.48-equation A.51.
Proof. Define w def= unn − u. Then by equation A.48 and equation A.49, w satisfies
wt + vwx − βvwv − qwvv = f for (t, x, v) ∈ [0, T ]× [0, 1]× [−5, 5], w(0, x, v) = w0, for (x, v) ∈ [0, 1]× [−5, 5],
(A.52)
where w0 def = g − u0. By multiplying w to equation A.52 and integrating with respect to dxdv, we have
1
2
d
dt ∫∫ [0,1]×[−5,5] |w|2dxdv + ∫∫ [0,1]×[−5,5] vwxwdxdv − ∫∫ [0,1]×[−5,5] qwvvwdxdv
= ∫∫ [0,1]×[−5,5] fwdxdv + ∫∫ [0,1]×[−5,5] βvwvwdxdv.
Then we take the integration by parts and obtain that
1
2
d
dt ∫∫ [0,1]×[−5,5] |w|2dxdv + 1 2 ∫ 5 −5 dv v(w(t, 1, v)2 − w(t, 0, v)2)
+ q ∫∫ [0,1]×[−5,5] |wv|2dxdv
= ∫∫ [0,1]×[−5,5] fwdxdv + ∫∫ [0,1]×[−5,5] βvwvwdxdv
+ q ∫ [0,1] wv(t, x, 5)w(t, x, 5)dx− q ∫ [0,1] wv(t, x,−5)w(t, x,−5)dx
def = I1 + I2 + I3 + I4.
We first define
A(t) def = ∣∣∣∣12 ∫ 5 −5 dv v(w(t, 1, v)2 − w(t, 0, v)2) ∣∣∣∣ . We now estimate I1-I4 on the right-hand side. By the Hölder inequality and Young’s inequality, we have
|I1| ≤ ‖f‖2‖w‖2 ≤ 1
2 ‖f‖22 +
1 2 ‖w‖22,
where we denote
‖h‖2 def = ∫∫ [0,1]×[−5,5] |h|2dxdv.
Similarly, we observe that
|I2| ≤ 5β‖wv‖2‖w‖2 ≤ ε‖wv‖22 + 25β2
4ε ‖w‖22,
for a sufficiently small ε > 0 as |v| ≤ 5. By equation A.50, we have
|I3 + I4| ≤ q‖wv(t, ·, 5)‖L2x‖w(t, ·, 5)− w(t, ·,−5)‖L2x + q‖wv(t, ·, 5)− wv(t, ·,−5)‖L2x‖w(t, ·,−5)‖L2x ≤ 2q C.
Altogether, we have
d dt ‖w‖22 + 2(q − ε)‖wv‖22 ≤ ‖f‖2 +
( 1 + 25β2
2ε
) ‖w‖22 +A(t) + 2q C.
We integrate with respect to the temporal variable on [0, t] and obtain
‖w(t)‖22 + 2(q − ε) ∫ t
0
‖wv(s)‖22ds
≤ ‖w(0)‖22 + ∫ t
0
( ‖f(s)‖2 + ( 1 + 25β2
2ε
) ‖w(s)‖22 +A(s) + 2q C ) ds.
By equation A.493, we have ∫ t 0 A(s)ds ≤ L 2 . Thus, by the Grönwall inequality, we have
‖w(t)‖22 + 2(q − ε) ∫ t
0
‖wv(s)‖22ds
≤ ( ‖w0‖22 + L
2
) exp [( 1 + 25β2
2ε
) t ] + ∫ t 0 ‖f(s)‖22ds+ 2q Ct,
where w0(x, v) = g(x, v)− u0(x, v). This completes the proof for the theorem.
Regarding the derivatives ∂tw and ∂xw we can obtain the similar estimates as follows. Corollary A.11. Let u and unn be the classical solutions to equation A.48 and equation A.49, respectively. Assume that equation A.51 holds. Then for z = t or x we have
sup 0≤t≤T
‖∂zunn(t)− ∂zu(t)‖22 + 2(q − ε) ∫ T
0
‖∂v∂zunn(s)− ∂v∂zu(s)‖22ds
≤ ( ‖∂zg − ∂zu0‖22 + L
2
) exp [( 1 + 25β2
2ε
) T ] + ∫ T 0 ‖∂zf(s)‖22ds+ 2q CT,
for any ε ∈ (0, q), where L, u0, f, g, β, q,m, , and C are given in equation A.48-equation A.51.
Proof. For both ∂z = ∂t and ∂x, we take ∂z onto equation A.52 and obtain
(∂zw)t + v(∂zw)x − βv(∂zw)v − q(∂zw)vv = ∂zf for (t, x, v) ∈ [0, T ]× [0, 1]× [−5, 5], (∂zw)(0, x, v) = (∂zw)0, for (x, v) ∈ [0, 1]× [−5, 5],
(A.53)
where (∂zw)0 def = ∂zg − ∂zu0. Then the proof is the same as the one for Theorem A.10 for ∂zw replacing the role of w. This completes the proof.
Finally, we can also obtain the regularity estimates for the derivative ∂vw as follows: Theorem A.12. Let u and unn be the classical solutions to equation A.48 and equation A.49, respectively. Assume that equation A.51 holds. Then we have
sup 0≤t≤T
‖∂vunn(t)− ∂vu(t)‖22 + 2(q − ε) ∫ T
0
‖∂vvunn(s)− ∂vvu(s)‖22ds
≤ (L+ ‖∂xg − ∂xu0‖22 + ‖∂vg − ∂vu0‖22) exp [( 2 + 25β2
2ε
) T ] +
∫ T 0 (‖∂xf(s)‖22 + ‖∂vf(s)‖22)ds+ 4q CT,
for any ε ∈ (0, q), where L, u0, f, g, β, q,m, , and C are given in equation A.48-equation A.51.
Proof. we take ∂v onto equation A.52 and obtain
(∂vw)t + wx + v(∂vw)x − βv(∂vw)v − q(∂vw)vv = ∂vf, for (t, x, v) ∈ [0, T ]× [0, 1]× [−5, 5], (A.54)
where (∂vw)0 def = ∂vg − ∂vu0. By multiplying ∂vw to equation A.54 and integrating with respect to dxdv, we have
1
2
d
dt ∫∫ [0,1]×[−5,5] |∂vw|2dxdv + ∫∫ [0,1]×[−5,5] v(∂vw)x(∂vw)dxdv
− ∫∫
[0,1]×[−5,5] q(∂vw)vv∂vwdxdv
= ∫∫ [0,1]×[−5,5] (−wx + ∂vf)∂vwdxdv + ∫∫ [0,1]×[−5,5] βv(∂vw)v∂vwdxdv.
Then we take the integration by parts and obtain that
1
2
d
dt ∫∫ [0,1]×[−5,5] |wv|2dxdv
+ 1
2 ∫ [−5,5] ( v(∂vw) 2(t, 1, v)− v(∂vw)2(t, 0, v) ) dv
+ q ∫∫ [0,1]×[−5,5] |wvv|2dxdv
= ∫∫ [0,1]×[−5,5] ∂vfwvdxdv + ∫∫ [0,1]×[−5,5] βvwvvwvdxdv
+ q ∫ [0,1] wvv(t, x, 5)wv(t, x, 5)dx− q ∫ [0,1] wvv(t, x,−5)wv(t, x,−5)dx
− ∫∫
[0,1]×[−5,5] wxwvdxdv
def = I1 + I2 + I3 + I4 + I5.
We first define
B(t) def = ∣∣∣∣12 ∫ 5 −5 dv v(∂vw(t, 1, v) 2 − ∂vw(t, 0, v)2) ∣∣∣∣ .
We now estimate I1-I4 on the right-hand side. We now estimate I1-I5 on the right-hand side. By the Hölder inequality and Young’s inequality, we have
|I1| ≤ ‖∂vf‖2‖wv‖2 ≤ 1
2 ‖∂vf‖22 +
1 2 ‖wv‖22,
where we denote ‖h‖2 def = ∫∫ [0,1]×[−5,5] |h|2dxdv.
Similarly, we observe that
|I2| ≤ 5β‖wvv‖2‖wv‖2 ≤ ε‖wvv‖22 + 25β2
4ε ‖wv‖22,
for a sufficiently small ε > 0 as |v| ≤ 5. By equation A.50 and equation A.51, we have
|I3 + I4| ≤ q‖wvv(t, ·, 5)‖L2x‖wv(t, ·, 5)− wv(t, ·,−5)‖L2x + q‖wvv(t, ·, 5)− wvv(t, ·,−5)‖L2x‖wv(t, ·,−5)‖L2x ≤ 2q C.
Finally, we have
|I5| ≤ ‖wx‖2‖wv‖2 ≤ 1
2 ‖wx‖22 +
1 2 ‖wv‖22.
Altogether, we have
d dt ‖wv‖22 + 2(q − ε)‖wvv‖22 ≤ ‖∂vf‖2 + ‖wx‖22 +
( 2 + 25β2
2ε
) ‖wv‖22 + 5L
2 + 2q C.
Then we take the integration with respect to the temporal variable on [0, t] and obtain that
‖wv(t)‖22 + ∫ t
0
ds 2(q − ε)‖wvv(s)‖22
≤ ‖wv(0)‖22 + ∫ t
0
ds ( ‖∂vf(s)‖2 + ‖wx(s)‖22 + ( 2 + 25β2
2ε
) ‖wv(s)‖22 +B(s) + 2q C ) .
By equation A.493, we have ∫ t 0 ds B(s) ≤ L 2 . Thus, by the Grönwall inequality, we have
‖wv(t)‖22 + 2(q − ε) ∫ t
0
‖wvv(s)‖22ds
≤ ( L
2 + ‖∂vw0‖22
) exp [( 2 + 25β2
2ε
) t ] + ∫ t 0 (‖wx(s)‖22 + ‖∂vf(s)‖22)ds+ 2q Ct, (A.55)
where ∂vw0(x, v) = g(x, v) − u0(x, v). Then we use Corollary A.11 for an upper-bound of ‖∂xw(s)‖22 and obtain that∫ T
0
‖wx(s)‖22ds
≤ ( L
2 + ‖∂xg − ∂xu0‖22
) exp [( 1 + 25β2
2ε
) T ] + ∫ T 0 ‖∂xf(s)‖22ds+ 2q CT.
Then by equation A.55, we obtain that
sup 0≤t≤T
‖wv(t)‖22 + 2(q − ε) ∫ T
0
‖wvv(s)‖22ds
≤ (L+ ‖∂xw0‖22 + ‖∂vw0‖22) exp [( 2 + 25β2
2ε
) T ] +
∫ T 0 (‖∂xf(s)‖22 + ‖∂vf(s)‖22)ds+ 4q CT,
where ∂vw0(x, v) = ∂vg(x, v)− ∂vu0(x, v). This completes the proof for the theorem.

heading: A.4 THE POISSON EQUATION
text: We consider the Poisson equation equation with Dirichlet boundary condition:
−4u = f in Ω, u = g on ∂Ω.
Suppose there exists g̃ ∈ H2(Ω̄) s.t. g̃|∂Ω = g (A.56)
Then, the equation can be written by:
−4v = f̃ in Ω, v = 0 on ∂Ω,
where v = u−g̃, f̃ = f−4g̃. Therefore, we assume the homogeneous Dirichlet boundary condition provided (A.56).
Now, let u be a strong solution of
−4u = f in Ω, u = 0 on ∂Ω,
(A.57)
and let unn be a neural network such that
−4unn = fnn in Ω, unn = 0 on ∂Ω.
(A.58)
Here, we can set the boundary to be zero by multiplying B(x), where B(x) is a smooth function satisfying B(x) = {
0, x ∈ ∂Ω 6= 0, x ∈ Ω . By subtracting (A.58) from (A.57), we get
−4(u− unn) = (f − fnn) in Ω, u− unn = 0 on ∂Ω.
(A.59)
Then we apply below theorem to (A.59) to get the convergence results. Theorem A.13. (Theorem 6.3.5 in Evans (2010)) Let m be a nonnegative integer. Suppose that u ∈ H10 (Ω) is a weak solution of the boundary-value problem (A.57). Assume ∂Ω is Cm+2. Then,
‖u‖Hm+2(Ω) ≤ C(‖f‖Hm(Ω) + ‖u‖L2(Ω)), (A.60)
Furthermore, if u is the unique solution of A.57, then
‖u‖Hm+2(Ω) ≤ C‖f‖Hm(Ω). (A.61)
By applying (A.61) to (A.59), we obtain
‖u− unn‖Hm+2(Ω) ≤ C‖f − fnn‖Hm(Ω).
where the right hand side corresponds to LossGE(unn;m, 2).

heading: B ADDITIONAL EXPERIMENTAL RESULTS

heading: B.1 DEPENDENCY ON LEARNING RATES
text: In this subsection, we show several experiments that shows the proposed loss functions generally performs better in different learning rates. We first show the results for Burgers’ equation. In Figure 5, we show the test errors versus training epochs plot for different learning rates. We used 10−3, 10−4, 10−5 as learning rates and we observe that H2 loss performs best followed by H1 and L2 loss functions.
We next present the similar experiments for the high-dimensional Poisson equation. We trained 30 neural networks with different initializations with different learning rates. The average errors are presented in Figure 6. As the same in the Burgers’ equation, our loss functions performs better than the traditional one in all learning rates.

heading: B.2 HIGH-DIMENSIONAL POISSON EQUATION
text: In this subsection, we consider the high-dimensional Poisson equation with different boundary condition. In subsection 5.1, we pointed out that the ”difficulty” of learning sin(kx) increases as k increases. As a generalization of the argument, we consider the following PDEs:
−4u = (kπ) 2
4 d∑ i=1 sin( kπ 2 xi), for x ∈ Ω = (0, 1)d,
u = d∑ i=1 sin( kπ 2 xi), for x ∈ ∂Ω,
for d=10, k = 1,3, and 5. As one can see in Figure 7, the improvement of Sobolev training gets bigger as k increases. This observation coincides with the one in section 5, as we expected. Moreover, we present the comparison of training time to meet a certain error value for different loss functions in Figure 7. The result shows that it is advantageous to use the proposed loss functions in time, even in high-dimensional case.

Review summaries:
The idea of using neural networks to approximate the solutions of the pdes is very interesting, specially in high-dimensional setting where classical approaches fail to scale. Although there has been many efforts in this direction, there are still open venues to explore. One of the most important aspect is the choice of loss function to guide the training of the neural network. And the paper's aim is to address this issue by proposing Sobolev norm as the loss function instead of the commonly used $L^2$-norm. The Sobolev norm includes additional term about derivatives of the error. The main claim is that with the inclusion of the additional term, the convergence of the neural network training becomes faster. This is the basic promise of the paper.  

Sobolev training of neural networks, which augments the standard loss function with terms that penalize discrepancies between the derivatives of the network and target functions, has been shown empirically to improve data-efficiency. Intuitively, one would expect that it also aids generalization in settings where the target function is sufficiently smooth. This manuscript proposes augmenting the loss functions used to represent the solutions of partial differential equations with terms penalizing the Sobolev norm of the solution, its initial condition, and the boundary condition. The motivation for this approach is clear because data -efficiency is of the utmost importance in PDE learning problems where the data could be very difficult to access.