{'fraction_non_alphanumeric': 0.042784316865508314, 'fraction_numerical': 0.012571062535030828, 'mean_word_length': 4.599760872814228, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
DeepXDE: A deep learning library for solving differential equations

Deep learning has achieved remarkable success in diverse applications; however, its use in solving partial differential equations (PDEs) has emerged only recently. Here, we present an overview of physics-informed neural networks (PINNs), which embed a PDE into the loss of the neural network using automatic differentiation. The PINN algorithm is simple, and it can be applied to different types of PDEs, including integro-differential equations, fractional PDEs, and stochastic PDEs. Moreover, from the implementation point of view, PINNs solve inverse problems as easily as forward problems. We propose a new residual-based adaptive refinement (RAR) method to improve the training efficiency of PINNs. For pedagogical reasons, we compare the PINN algorithm to a standard finite element method. We also present a Python library for PINNs, DeepXDE, which is designed to serve both as an education tool to be used in the classroom as well as a research tool for solving problems in computational science and engineering. Specifically, DeepXDE can solve forward problems given initial and boundary conditions, as well as inverse problems given some extra measurements. DeepXDE supports complex-geometry domains based on the technique of constructive solid geometry, and enables the user code to be compact, resembling closely the mathematical formulation. We introduce the usage of DeepXDE and its customizability, and we also demonstrate the capability of PINNs and the user-friendliness of DeepXDE for five different examples. More broadly, DeepXDE contributes to the more rapid development of the emerging Scientific Machine Learning field.

1. Introduction. In the last 15 years, deep learning in the form of deep neural networks (NNs), has been used very effectively in diverse applications [32], such as computer vision and natural language processing. Despite the remarkable success in these and related areas, deep learning has not yet been widely used in the field of scientific computing. However, more recently, solving partial differential equations (PDEs), e.g., in the standard differential form or in the integral form, via deep learning has emerged as a potentially new sub-field under the name of Scientific Machine Learning (SciML) [4]. In particular, we can replace traditional numerical discretization methods with a neural network that approximates the solution to a PDE.
To obtain the approximate solution of a PDE via deep learning, a key step is to constrain the neural network to minimize the PDE residual, and several approaches have been proposed to accomplish this. Compared to the traditional mesh-based methods, such as the finite difference method (FDM) and the finite element method (FEM), deep learning could be a mesh-free approach by taking advantage of the automatic differentiation [47], and could break the curse of dimensionality [45,18]. Among these approaches, some can only be applied to particular types of problems, such as image-like input domain [28,33,58] or parabolic PDEs [6,19]. Some researchers adopt the variational form of PDEs and minimize the corresponding energy functional [16,20]. However, not all PDEs can be derived from a known functional, and thus Galerkin type projections have also been considered [39]. Alternatively, one could use the PDE in strong form directly [15,52,30,31,7,50,47]; in this form, automatic differentiation could be used directly to avoid truncation errors and the numerical quadrature errors of variational forms. This strong form approach was introduced in [47] coining the term physics-informed neural networks (PINNs). An attractive feature of PINNs is that it can be used to solve inverse problems with minimum change of the code for forward problems [47,48,51,21,13]. In addition, PINNs have been further extended to solve integro-differential equations (IDEs), fractional differential equations (FDEs) [42], and stochastic differential equations (SDEs) [57,55,41,56].
In this paper, we present various PINN algorithms implemented in a Python library DeepXDE 1 , which is designed to serve both as an education tool to be used in the classroom as well as a research tool for solving problems in computational science and engineering (CSE). DeepXDE can be used to solve multi-physics problems, and supports complex-geometry domains based on the technique of constructive solid geometry (CSG), hence avoiding tedious and time-consuming computational geometry tasks. By using DeepXDE, time-dependent PDEs can be solved as easily as steady states by only defining the initial conditions. In addition to the main workflow of DeepXDE, users can readily monitor and modify the solution process via callback functions, e.g., monitoring the Fourier spectrum of the neural network solution, which can reveal the learning mode of the NN Figure 2. Last but not least, DeepXDE is designed to make the user code stay compact and manageable, resembling closely the mathematical formulation.
The paper is organized as follows. In section 2, after briefly introducing deep neural networks and automatic differentiation, we present the algorithm, approximation theory, and error analysis of PINNs, and make a comparison between PINNs and FEM. We then discuss how to use PINNs to solve integro-differential equations and inverse problems. In addition, we propose the residual-based adaptive refinement (RAR) method to improve the training efficiency of PINNs. In section 3, we introduce the usage of our library, DeepXDE, and its customizability. In section 4, we demonstrate the capability of PINNs and friendly use of DeepXDE for five different examples. Finally, we conclude the paper in section 5.
2. Algorithm and theory of physics-informed neural networks. In this section, we first provide a brief overview of deep neural networks and automatic differentiation, and present the algorithm and theory of PINNs for solving PDEs. We then make a comparison between PINNs and FEM, and discuss how to use PINNs to solve integro-differential equations and inverse problems. Next we propose RAR, an efficient way to select the residual points adaptively during the training process.
2.1. Deep neural networks. Mathematically, a deep neural network is a particular choice of a compositional function. The simplest neural network is the feedforward neural network (FNN), also called multilayer perceptron (MLP), which applies linear and nonlinear transformations to the inputs recursively. Although many different types of neural networks have been developed in the past decades, such as the convolutional neural network and the recurrent neural network. In this paper we consider FNN, which is sufficient for most PDE problems, and residual neural network (ResNet), which is easier to train for deep networks. However, it is straightforward to employ other types of neural networks.
Let N L (x) : R din → R dout be a L-layer neural network, or a (L − 1)-hidden layer neural network, with N neurons in the -th layer (N 0 = d in , N L = d out ). Let us denote the weight matrix and bias vector in the -th layer by W ∈ R N ×N −1 and b ∈ R N , respectively. Given a nonlinear activation function σ, which is applied element-wisely, the FNN is recursively defined as follows: see also a visualization of a neural network in Figure 1. Commonly used activation functions include the logistic sigmoid 1/(1 + e −x ), the hyperbolic tangent (tanh), and the rectified linear unit (ReLU, max{x, 0}).

Automatic differentiation.
In PINNs, it is required to compute the derivatives of the network outputs with respect to the network inputs. There are four possible methods for computing the derivatives [5,38]: (1) hand-coded analytical derivative; (2) finite difference or other numerical approximations; (3) symbolic differentiation (used in software programs such as Mathematica, Maxima, and Maple); and (4) automatic differentiation (AD, also called algorithmic differentiation). In deep learning, the derivatives are evaluated using backpropagation [49], a specialized technique of AD.
Considering the fact that the neural network represents a compositional function, then AD applies the chain rule repeatedly to compute the derivatives. There are two steps in AD: one forward pass to compute the values of all variables, and one backward pass to compute the derivatives. To demonstrate AD, we consider a FNN of only one hidden layer with two inputs x 1 and x 2 and one output y: The forward pass and backward pass of AD for computing the partial derivatives ∂y ∂x1 and ∂y ∂x2 at (x 1 , x 2 ) = (2, 1) are shown in Table 1.
We can see that AD only requires one forward pass and one backward pass to compute all the partial derivatives, no matter what the input dimension is. In contrast, using finite differences computing each partial derivative ∂y ∂xi requires two function valuations y(x 1 , . . . , x i , . . . , x din ) and y(x 1 , . . . , x i + ∆x i , . . . , x din ) for some small number ∆x i , and thus in total d in + 1 forward passes are required to evaluate all the partial derivatives. Hence, AD is much more efficient than finite difference when the input dimension is high (see [5,38] for more details of the comparison between AD and the other three methods). To compute n th -order derivatives, AD can be applied recursively n times. However, this nested approach may lead to inefficiency and numerical instability, and hence other methods, e.g., Taylor-Mode AD, have been developed for this purpose [8,9]. Finally we note that with AD we differentiate the NN and therefore we can deal with noisy data [42].

Forward pass
Backward pass where B(u, x) could be Dirichlet, Neumann, Robin, or periodic boundary conditions. For time-dependent problems, we consider time t as a special component of x, and Ω contains the temporal domain. The initial condition can be simply treated as a special type of Dirichlet boundary condition on the spatio-temporal domain.
Schematic of a PINN for solving the diffusion equation ∂u ∂t = λ ∂ 2 u ∂x 2 with mixed boundary conditions (BC) u(x, t) = g D (x, t) on Γ D ⊂ ∂Ω and ∂u ∂n (x, t) = g R (u, x, t) on Γ R ⊂ ∂Ω. The initial condition (IC) is treated as a special type of boundary conditions. T f and T b denote the two sets of residual points for the equation and BC/IC.

Procedure 2.1
The PINN algorithm for solving differential equations.
Step 2 Specify the two training sets T f and T b for the equation and boundary/initial conditions. Step 3 Specify a loss function by summing the weighted L 2 norm of both the PDE equation and boundary condition residuals.
Step 4 Train the neural network to find the best parameters θ * by minimizing the loss function L(θ; T ).
The algorithm of PINN [31,47] is shown in Procedure 2.1, and visually in the schematic of Figure 1 solving a diffusion equation ∂u ∂t = λ ∂ 2 u ∂x 2 with mixed boundary conditions u(x, t) = g D (x, t) on Γ D ⊂ ∂Ω and ∂u ∂n (x, t) = g R (u, x, t) on Γ R ⊂ ∂Ω. We explain each step as follows. In a PINN, we first construct a neural networkû(x; θ) as a surrogate of the solution u(x), which takes the input x and outputs a vector with the same dimension as u. Here, θ = {W , b } 1≤ ≤L is the set of all weight matrices and bias vectors in the neural networkû. One advantage of PINNs by choosing neural networks as the surrogate of u is that we can take the derivatives ofû with respect to its input x by applying the chain rule for differentiating compositions of functions using the automatic differentiation (AD), which is conveniently integrated in machine learning packages, such as TensorFlow [1] and PyTorch [43].
In the next step, we need to restrict the neural networkû to satisfy the physics imposed by the PDE and boundary conditions. In practice, we restrictû on some scattered points (e.g., randomly distributed points, or clustered points in the domain [37]), i.e., the training data T = {x 1 , x 2 , . . . , x |T | } of size |T |. In addition, T is comprised of two sets T f ⊂ Ω and T b ⊂ ∂Ω, which are the points in the domain and on the boundary, respectively. We refer T f and T b as the sets of "residual points".
To measure the discrepancy between the neural networkû and the constraints, we consider the loss function defined as the weighted summation of the L 2 norm of residuals for the equation and boundary conditions: and w f and w b are the weights. The loss involves derivatives, such as the partial derivative ∂û/∂x 1 or the normal derivative at the boundary ∂û/∂n = ∇û · n, which are handled via AD. In the last step, the procedure of searching for a good θ by minimizing the loss L(θ; T ) is called "training". Considering the fact that the loss is highly nonlinear and non-convex with respect to θ [10], we usually minimize the loss function by gradient-based optimizers, such as gradient descent, Adam [29], and L-BFGS [12]. We remark that based on our experience, for smooth PDE solutions L-BFGS can find a good solution with less iterations than Adam, because L-BFGS uses secondorder derivatives of the loss function, while Adam only relies on first-order derivatives. However, for stiff solutions L-BFGS is more likely to be stuck at a bad local minimum. The required number of iterations highly depends on the problem (e.g., the smoothness of the solution), and to check whether the network converges or not, we can monitor the loss function or the PDE residual using callback functions. We also note that acceleration of training can be achieved by using adaptive activation function that may remove bad local minima, see [24,23].
Unlike traditional numerical methods, for PINNs there is no guarantee of unique solutions, because PINN solutions are obtained by solving non-convex optimization problems, which in general do not have a unique solution. In practice, to achieve a good level of accuracy, we need to tune all the hyperparameters, e.g., network size, learning rate, and the number of residual points. The required network size depends highly on the smoothness of the PDE solution. For example, a small network (e.g., a few layers and twenty neurons per layer) is sufficient for solving the 1D Poisson equation, but a deeper and wider network is required for the 1D Burgers equation to achieve a similar level of accuracy. We also note that PINNs may converge to different solutions from different network initial values [42,51,21], and thus a common strategy is that we train PINNs from random initialization for a few times (e.g., 10 independent runs) and choose the network with the smallest training loss as the final solution.
In the algorithm of PINN introduced above, we enforce soft constraints of boundary/initial conditions through the loss L b . This approach can be used for complex domains and any type of boundary conditions. On the other hand, it is possible to enforce hard constraints for simple cases [30]. For example, when the boundary condition is u(0) = u(1) = 0 with Ω = [0, 1], we can simply choose the surrogate model asû(x) = x(x − 1)N (x) to satisfy the boundary condition automatically, where N (x) is a neural network.
We note that we have great flexibility in choosing the residual points T , and here we provide three possible strategies: 1. We can specify the residual points at the beginning of training, which could be grid points on a lattice or random points, and never change them during the training process. 2. In each optimization iteration, we could select randomly different residual points. 3. We could improve the location of the residual points adaptively during training, e.g., the method proposed in subsection 2.8. When the number of residual points required is very large, e.g., in multiscale problems, it is computationally expensive to calculate the loss and gradient in every iteration. Instead of using all residual points, we can split the residual points into small batches, and in each iteration we only use one batch to calculate the loss and update model parameters; this is the so-called "mini-batch" gradient descent. The aforementioned strategy (2), i.e., re-sampling in each step, is a special case of mini-batch gradient descent by choosing T = Ω with |T | = ∞.
Recent studies show that for function approximation, neural networks learn target functions from low to high frequencies [46,54], but here we show that the learning mode of PINNs is different due to the existence of high-order derivatives. For example, when we approximate the function f (x) = 5 k=1 sin(2kx)/(2k) in [−π, π] by a NN, the function is learned from low to high frequency ( Figure 2A). However, when we employ a PINN to solve the Poisson equation −f xx = 5 k=1 2k sin(2kx) with zero boundary conditions in the same domain, all frequencies are learned almost simultaneously ( Figure 2B). Interestingly, by comparing Figure 2A and Figure 2B we can see that at least in this case solving the PDE using a PINN is faster than approximating a function using a NN. We can monitor this training process using the callback functions in our library DeepXDE as discussed later.
2.4. Approximation theory and error analysis for PINNs. One fundamental question related to PINNs is whether there exists a neural network satisfying both the PDE equation and the boundary conditions, i.e., whether there exists a neural network that can simultaneously and uniformly approximate a function and its partial derivatives. To address this question, we first introduce some notation. Let Z d + be the set of d-dimensional nonnegative integers. For m = (m 1 , . . . , m d ) ∈ Z d + , we set |m| := m 1 + · · · + m d , and is continuous} is the space of continuous functions. Then, we recall the following theorem of derivative approximation using single hidden layer neural networks due to Pinkus [44].
Assume σ ∈ C m (R) and σ is not a polynomial. Then the space of single hidden layer neural nets , for any f ∈ C m 1 ,...,m s (R d ), any compact K ⊂ R d , and any ε > 0, there exists a g ∈ M(σ) satisfying max for all k ∈ Z d + for which k ≤ m i for some i. Theorem 2.1 shows that feed-forward neural nets with enough neurons can simultaneously and uniformly approximate any function and its partial derivatives. However, neural networks in practice have limited size. Let F denote the family of all the functions that can be represented by our chosen neural network architecture. The solution u is unlikely to belong to the family F, and we define u F = arg min f ∈F f −u as the best function in F close to u (Figure 3). Because we only train the neural network on the training set T , we define u T = arg min f ∈F L(f ; T ) as the neural network whose loss is at global minimum. For simplicity, we assume that u, u F and u T are well defined and unique. Finding u T by minimizing the loss is often computationally intractable [10], and our optimizer returns an approximate solutionũ T .
We can then decompose the total error E as [11] E :

Generalization error E gen
Optimization error E opt Fig. 3. Illustration of errors of a PINN. The total error consists of the approximation error, the optimization error, and the generalization error. Here, u is the PDE solution, u F is the best function close to u in the function space F , u T is the neural network whose loss is at a global minimum, andũ T is the function obtained by training a neural network.
The approximation error E app measures how closely u F can approximate u. The generalization error E gen is determined by the number/locations of residual points in T and the capacity of the family F. Neural networks of larger size have smaller approximation errors but could lead to higher generalization errors, which is called bias-variance tradeoff. Overfitting occurs when the generalization error dominates. In addition, the optimization error E opt stems from the loss function complexity and the optimization setup, such as learning rate and number of iterations. However, currently there is no error estimation for PINNs yet, and even quantifying the three errors for supervised learning is still an open research problem [36,35,25].

Comparison between PINNs and FEM.
To further explain the ideas of PINNs and to help those with the knowledge of FEM understand PINNs more easily, we make a comparison between PINNs and FEM point by point (Table 2): • In FEM we approximate the solution u by a piecewise polynomial with point values to be determined, while in PINNs we construct a neural network as the surrogate model parameterized by weights and biases. • FEM typically requires a mesh generation, while PINN is totally mesh-free, and we can use either a grid or random points. • FEM converts a PDE to an algebraic system, using the stiffness matrix and mass matrix, while PINN embeds the PDE and boundary conditions into the loss function. • In the last step, the algebraic system in FEM is solved exactly by a linear solver, but the network in PINN is learned by a gradient-based optimizer. At a more fundamental level, PINNs provide a nonlinear approximation to the function and its derivatives whereas FEM represent a linear approximation. Not available yet Partially available [14,26] 2.6. PINNs for solving integro-differential equations. When solving integrodifferential equations (IDEs), we still employ the automatic differentiation technique to analytically derive the integer-order derivatives, while we approximate integral operators numerically using classical methods (Figure 4) [42], such as Gaussian quadra-ture. Therefore, we introduce a fourth error component, the discretization error E dis , due to the approximation of the integral by Gaussian quadrature.
For example, when solving we first use Gaussian quadrature of degree n to approximate the integral and then we use a PINN to solve the following PDE instead of the original equation PINNs can also be easily extended to solve FDEs [42] and SDEs [57,55,41,56], but we do not discuss here such cases due to the page limit.  Integral & Integro-differential Differential Fig. 4. Schematic illustrating the modificaiton of the PINN algorithm for solving integrodifferential equations. We employ the automatic differentiation technique to analytically derive the integer-order derivatives, and we approximate integral operators numerically using standard methods. (The figure is reproduced from [42].) 2.7. PINNs for solving inverse problems. In inverse problems, there are some unknown parameters λ in (2.1), but we have some extra information on some points T i ⊂ Ω besides the differential equation and boundary conditions: From the implementation point of view, PINNs solve inverse problems as easily as forward problems [47,42]. The only difference between solving forward and inverse problems is that we add an extra loss term to (2.2): We then optimize θ and λ together, and our solution is θ * , λ * = arg min θ,λ L(θ, λ; T ).

2.8.
Residual-based adaptive refinement (RAR). As we discussed in subsection 2.3, the residual points T are usually randomly distributed in the domain. This works well for most cases, but it may not be efficient for certain PDEs that exhibit solutions with steep gradients. Take the Burgers equation as an example, intuitively we should put more points near the sharp front to capture the discontinuity well. However, it is challenging, in general, to design a good distribution of residual points for problems whose solution is unknown. To overcome this challenge, we propose a residual-based adaptive refinement (RAR) method to improve the distribution of residual points during training process (Procedure 2.2), conceptually similar to FEM refinement methods [2]. The idea of RAR is that we will add more residual points in the locations where the PDE residual f x; ∂û ∂x1 , . . . , ∂û ∂x d ; ∂ 2û ∂x1∂x1 , . . . , ∂ 2û ∂x1∂x d ; . . . ; λ is large, and we repeat adding points until the mean residual is smaller than a threshold E 0 , where V is the volume of Ω.

Procedure 2.2 RAR for improving the distribution of residual points for training.
Step 1 Select the initial residual points T , and train the neural network for a limited number of iterations.
Step 2 Estimate the mean PDE residual E r in (2.3) by Monte Carlo integration, i.e., by the average of values at a set of randomly sampled locations S = {x 1 , x 2 , . . . , x |S| }: Step 3 Stop if E r < E 0 . Otherwise, add m new points with the largest residuals in S to T , and go to Step 2.

DeepXDE usage and customization.
In this section, we introduce the usage of DeepXDE and how to customize DeepXDE to meet new problem requirements.

Usage. Compared to traditional numerical methods, the code written with
DeepXDE is much shorter and more comprehensive, resembling closely the mathematical formulation. Solving differential equations in DeepXDE is no more than specifying the problem using the build-in modules, including computational domain (geometry and time), PDE equations, boundary/initial conditions, constraints, training data, neural network architecture, and training hyperparameters. The workflow is shown in Procedure 3.1 and Figure 5.
In DeepXDE, the built-in primitive geometries include interval, triangle, rectangle, polygon, disk, cuboid and sphere. Other geometries can be constructed from these primitive geometries using three boolean operations: union (|), Procedure 3.1 Usage of DeepXDE for solving differential equations.
Step 1 Specify the computational domain using the geometry module.
Step 2 Specify the PDE using the grammar of TensorFlow.
Step 3 Specify the boundary and initial conditions.
Step 4 Combine the geometry, PDE and boundary/initial conditions together into data.PDE or data.TimePDE for time-independent problems or for timedependent problems, respectively. To specify training data, we can either set the specific point locations, or only set the number of points and then DeepXDE will sample the required number of points on a grid or randomly.
Step 5 Construct a neural network using the maps module.
Step 6 Define a Model by combining the PDE problem in Step 4 and the neural net in Step 5.
Step 7 Call Model.compile to set the optimization hyperparameters, such as optimizer and learning rate. The weights in (2.2) can be set here by loss weights.
Step 8 Call Model.train to train the network from random initialization or a pretrained model using the argument model restore path. It is extremely flexible to monitor and modify the training behavior using callbacks.
Step 9 Call Model.predict to predict the PDE solution at different locations.

difference (-) and intersection (&).
This technique is called constructive solid geometry (CSG), see Figure 6 for examples. CSG supports both two-dimensional and three-dimensional geometries.
DeepXDE supports four standard boundary conditions, including Dirichlet (DirichletBC), Neumann (NeumannBC), Robin (RobinBC), and periodic (PeriodicBC), and a more general BC can be defined using OperatorBC. The initial condition can be defined using IC. There are two types of neural networks available in DeepXDE: feed-forward neural network (maps.FNN) and residual neural network (maps.ResNet). It is also convenient to choose different training hyperparameters, such as loss types, metrics, optimizers, learning rate schedules, initializations and regularizations.
In addition to solving differential equations, DeepXDE can also be used to approximate functions from multi-fidelity data [40], and learn nonlinear operators [34]. 3.2. Customizability. All the components of DeepXDE are loosely coupled, and thus DeepXDE is well-structured and highly configurable. In this subsection, we discuss how to customize DeepXDE to address new problem requirements, e.g., new geometry or network architecture.
3.2.1. Geometry. As we introduced above, DeepXDE has already supported 7 basic geometries and the CSG technique. However, it is still possible that the user needs a new geometry, which cannot be constructed in DeepXDE. In this situation, a new geometry can be defined as shown in Procedure 3.2. Currently DeepXDE does not support accurately descriptions of complex curvilinear boundaries; however, a future extension could be incorporation of the non-uniform rational basis spline (NURBS) [22] for such representations. """Compute the unit normal at x for Neumann or Robin boundary conditions.""" 8 def periodic_point(self, x, component): 9 """Compute the periodic image of x for periodic boundary condition.""" 10 def uniform_points(self, n, boundary=True): 11 """Compute the equispaced point locations in the geometry.""" 12 def random_points(self, n, random="pseudo"): 13 """Compute the random point locations in the geometry.""" 14 def uniform_boundary_points(self, n): 15 """Compute the equispaced point locations on the boundary.""" 16 def random_boundary_points(self, n, random="pseudo"): def targets(self): 10 """Return the targets of the net outputs.""" 11 def build(self): 12 """Construct the network.""" 3.2.3. Callbacks. It is usually a good strategy to monitor the training process of the neural network, and then make modifications in real time, e.g., change the learning rate. In DeepXDE, this can be implemented by adding a callback function, and here we only list a few commonly used ones already implemented in DeepXDE: • ModelCheckpoint, which saves the model after certain epochs or when a better model is found.
• OperatorPredictor, which calculates the values of the operator applying on the outputs.
• FirstDerivative, which calculates the first derivative of the outpus with respect to the inputs. This is a special case of OperatorPredictor with the operator being the first derivative.
• MovieDumper, which dumps the movie of the function during the training progress, and/or the movie of the spectrum of its Fourier transform.
It is very convenient to add other callback functions, which will be called at different stages of the training process, see Procedure 3.4.  −∆u(x, y) = 1, (x, y) ∈ Ω, u(x, y) = 0, (x, y) ∈ ∂Ω. Table 3 Hyperparameters used for the following 5 examples. "Adam, L-BFGS" represents that we first use Adam for a certain number of iterations, and then switch to L-BFGS. The optimizer L-BFGS does not require learning rate, and the neural network is trained until convergence, so the number of iterations is also ignored for L-BFGS. We choose 1200 and 120 random points drawn from a uniform distribution as T f and T b , respectively. The PINN solution is given in Figure 7B. For comparison, we also present the numerical solution obtained by using the spectral element method (SEM) [27] ( Figure 7A). The result of the absolute error is shown in Figure 7C.

RAR for Burgers equation.
We first consider the 1D Burgers equation: Let ν = 0.01/π. Initially, we randomly select 2500 points (spatio-temporal domain) as the residual points, and then 40 more residual points are added adaptively via RAR developed in subsection 2.8 with m = 1 and E 0 = 0.005. We compare the PINN solution with RAR and the PINN solution based on 2540 randomly selected training data (Figure 8 A and B), and demonstrate that PINN with RAR can capture the discontinuity much better. For a comparison, the finite difference solutions using central difference scheme for space discretization and forward Euler scheme for time discretization for the Burgers equation in the conservative form are also shown in Figure 8A. Here, we also present two examples for the residual points added via the RAR method. As shown in Figure 8 C and D, the added points (green crosses) are quite close to the sharp interface, which indicates the effectiveness of RAR.
We further solve the following two-dimensional Burgers equation using the RAR: where u and v are the velocities along the x and y directions, respectively. In addition, Re is a non-dimensional number (Reynolds number) defined as Re = U L/ν, in which U and L are respectively the characteristic velocity and length, and ν is the kinematic viscosity of fluid. The exact solution can be obtained [3] as using the Dirichlet boundary conditions on all boundaries. In the present study, Re is set to be 5000, which is quite challenging due to the fact that the high Reynolds number leads to steep gradient in the solution. Initially, 200 residual points are randomly sampled in the spatio-temporal domain, and 5000 and 1000 random points are used for each initial and boundary conditions, respectively. We only add 10 extra residual points using RAR, and the total number of the residual points is 210 after convergence. For comparison, we also test the case without the RAR using 210 randomly sampled residual points. The results are displayed in Figure 8 E and F, demonstrating the effectiveness of RAR.

4.3.
Inverse problem for the Lorenz system. Consider the parameter identification problem of the following Lorenz system with the initial condition (x(0), y(0), z(0)) = (− 8,7,27), where ρ, σ and β are the three parameters to be identified from the observations at certain times. The observations are produced by solving the above system to t = 3 using Runge-Kutta (4,5) with the underlying true parameters (ρ, σ, β) = (10, 15, 8/3). We choose 400 uniformly distributed random points and 25 equispaced points as the residual points T f and T i , respectively. The evolution trajectories of ρ, σ and β are presented in Figure 9A, with the final identified values of (ρ, σ, β) = (10.002, 14.999, 2.668).

Inverse problem for diffusion-reaction systems.
A diffusion-reaction system in porous media for the solute concentrations C A , C B and C C (A + 2B → C) is described by where D = 2 × 10 −3 is the effective diffusion coefficient, and k f = 0.1 is the effective reaction rate. Because D and k f depend on the pore structure and are difficult to measure directly, we estimate D and k f based on 40000 observations of the concentrations C A and C B in the spatio-temporal domain. The identified D (1.98 × 10 −3 ) and k f (0.0971) are displayed in Figure 9B, which agree well with their true values.  with the exact solution y(x) = e −x cosh x. We solve this IDE using the method in subsection 2.6, and approximate the integral using Gaussian-Legendre quadrature of degree 20. The L 2 relative error is 2 × 10 −3 , and the solution is shown in Figure 10. 5. Concluding Remarks. In this paper, we present the algorithm, approximation theory, and error analysis of the physics-informed neural networks (PINNs) for solving different types of partial differential equations (PDEs). Compared to the traditional numerical methods, PINNs employ automatic differentiation to handle differential operators, and thus they are mesh-free. Unlike numerical differentiation, automatic differentiation does not differentiate the data and hence it can tolerate noisy data for training. We also discuss how to extend PINNs to solve other types of differential equations, such as integro-differential equations, and also how to solve inverse problems. In addition, we propose a residual-based adaptive refinement (RAR) method to improve the distribution of residual points during the training process, and thus increase the training efficiency.
To benefit both the education and the computational science communities, we have developed the Python library DeepXDE, an implementation of PINNs. By introducing the usage of DeepXDE, we show that DeepXDE enables user codes to be compact and follow closely the mathematical formulation. We also demonstrate how to customize DeepXDE to meet new problem requirements. Our numerical examples for forward and inverse problems verify the effectiveness of PINNs and the capability of DeepXDE. Scientific machine learning is emerging as a new and potentially powerful alternative to classical scientific computing, so we hope that libraries such as DeepXDE will accelerate this development and will make it accessible to the classroom but also to other researchers who may find the need to adopt PINN-like methods in their research, which can be very effective especially for inverse problems.
Despite the aforementioned advantages, PINNs still have some limitations. For forward problems, PINNs are currently slower than finite elements but this can be alleviated via offline training [58,53]. For long time integration, one can also use time-parallel methods to simultaneously compute on multiple GPUs for shorter time domains. Another limitation is the search for effective neural network architectures, which is currently done empirically by users; however, emerging meta-learning techniques can be used to automate this search, see [59,17]. Moreover, while here we enforce the strong form of PDEs, which is easy to be implemented by automatic differentiation, alternative weak/variational forms may also be effective, although they require the use of quadrature grids. Many other extensions for multi-physics and multi-scale problems are possible across different scientific disciplines by creatively designing the loss function and introducing suitable solution spaces. For instance, in the five examples we present here, we only assume data on scattered points, however, in geophysics or biomedicine we may have mixed data in the form of images and point measurements. In this case, we can design a composite neural network consisting of one convolutional neural network and one PINN sharing the same set of parameters, and minimize the total loss which could be a weighted summation of multiple losses from each neural network.