DeepXDE: A deep learning library for solving differential equations

TL;DR

DeepXDE is a Python library that facilitates solving various types of differential equations using physics-informed neural networks, offering an accessible, flexible, and efficient tool for research and education in scientific machine learning.

Contribution

The paper introduces DeepXDE, a comprehensive Python library for PINNs that supports complex geometries, inverse problems, and various PDE types, enhancing usability and educational value.

Findings

DeepXDE effectively solves forward and inverse PDE problems.

PINNs with RAR improve training efficiency.

DeepXDE demonstrates versatility across multiple PDE examples.

Abstract

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.

Figures

Tables

Table 1. Table 1: Example of AD to compute the partial derivatives ∂ y ∂ x 1 𝑦 subscript 𝑥 1 \frac{\partial y}{\partial x_{1}} and ∂ y ∂ x 2 𝑦 subscript 𝑥 2 \frac{\partial y}{\partial x_{2}} at ( x 1 , x 2 ) = ( 2 , 1 ) subscript 𝑥 1 subscript 𝑥 2 2 1 (x_{1},x_{2})=(2,1) .

Forward pass	Backward pass
$x_1=2$	$\frac{\partial y}{\partial y}=1$
$x_2=1$
$v=-2x_1+3x_2+0.5=-0.5$	$\frac{\partial y}{\partial h}=\frac{\partial \left(2h-1\right)}{\partial h}=2$
$h=\tanh v\approx -0.462$	$\frac{\partial y}{\partial v}=\frac{\partial y}{\partial h}\frac{\partial h}{\partial v}=\frac{\partial y}{\partial h}{\text{sech}}^2\left(v\right)\approx 1.573$
$y=2h-1=-1.924$	$\frac{\partial y}{\partial x_1}=\frac{\partial y}{\partial v}\frac{\partial v}{\partial x_1}=\frac{\partial y}{\partial v}\times \left(-2\right)=-3.146$
	$\frac{\partial y}{\partial x_2}=\frac{\partial y}{\partial v}\frac{\partial v}{\partial x_2}=\frac{\partial y}{\partial v}\times 3=4.719$

Table 2. Table 2: Comparison between PINN and FEM.

	PINN	FEM
Basis function	Neural network (nonlinear)	Piecewise polynomial (linear)
Parameters	Weights and biases	Point values
Training points	Scattered points (mesh-free)	Mesh points
PDE embedding	Loss function	Algebraic system
Parameter solver	Gradient-based optimizer	Linear solver
Errors	${ℰ}_{\text{app}}$, ${ℰ}_{\text{gen}}$ and ${ℰ}_{\text{opt}}$ (Section 2.4)	Approximation/quadrature errors
Error bounds	Not available yet	Partially available [14, 26]

Table 3. 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.

Example	NN Depth	NN Width	Optimizer	Learning rate	# Iterations
1	4	50	Adam, L-BFGS	0.001	50000
2	3	20	Adam, L-BFGS	0.001	15000
3	3	40	Adam	0.001	60000
4	3	20	Adam	0.001	80000
5	4	20	L-BFGS	-	-

Equations

$${\mathcal{N}}^0(\text{x})=\text{x}\in {\mathbb{R}}^{d_{\text{in}}},$$

$${\mathcal{N}}^{\ell }(\text{x})=\sigma ({\mathbf{W}}^{\ell }{\mathcal{N}}^{\ell -1}(\text{x})+{\mathbf{b}}^{\ell })\in {\mathbb{R}}^{N_{\ell }},\text{for}1\le \ell \le L-1,$$

$${\mathcal{N}}^L(\text{x})={\mathbf{W}}^L{\mathcal{N}}^{L-1}(\text{x})+{\mathbf{b}}^L\in {\mathbb{R}}^{d_{\text{out}}};$$

$$v$$

$$h$$

$$y$$

$$f\left(\mathbf{x};\frac{\partial u}{\partial x_1},\dots ,\frac{\partial u}{\partial x_d};\frac{{\partial }^{2}u}{\partial x_1\partial x_1},\dots ,\frac{{\partial }^{2}u}{\partial x_1\partial x_d};\dots ;\mathbf{\lambda }\right)=0,\mathbf{x}\in \Omega ,$$

$$\mathcal{B}(u,\mathbf{x})=0\text{on}\partial \Omega ,$$

$$\mathcal{L}(\mathbf{\theta };\mathcal{T})=w_f{\mathcal{L}}_f(\mathbf{\theta };{\mathcal{T}}_f)+w_b{\mathcal{L}}_b(\mathbf{\theta };{\mathcal{T}}_b),$$

$${\mathcal{L}}_f(\mathbf{\theta };{\mathcal{T}}_f)$$

$${\mathcal{L}}_b(\mathbf{\theta };{\mathcal{T}}_b)$$

$$D^{\mathbf{m}}:=\frac{{\partial }^{|\mathbf{m}|}}{\partial x_1^{m_1}\dots \partial x_d^{m_d}}.$$

$$\mathcal{M}(\sigma ):=span\{\sigma (\mathbf{w}\cdot \mathbf{x}+b):\mathbf{w}\in {\mathbb{R}}^d,b\in \mathbb{R}\}$$

$$C^{{\mathbf{m}}^1,\dots ,{\mathbf{m}}^s}({\mathbb{R}}^d):={\cap }_{i=1}^s{C}^{{\mathbf{m}}^i}({\mathbb{R}}^d),$$

$$\underset{\mathbf{x}\in K}{\max }|D^{\mathbf{k}}f(\mathbf{x})-D^{\mathbf{k}}g(\mathbf{x})|<\epsilon ,$$

$$\mathcal{E}:=\Vert {\tilde{u}}_{\mathcal{T}}-u\Vert \le \underset{{\mathcal{E}}_{\text{opt}}}{\underbrace{\Vert {\tilde{u}}_{\mathcal{T}}-u_{\mathcal{T}}\Vert }}+\underset{{\mathcal{E}}_{\text{gen}}}{\underbrace{\Vert u_{\mathcal{T}}-u_{\mathcal{F}}\Vert }}+\underset{{\mathcal{E}}_{\text{app}}}{\underbrace{\Vert u_{\mathcal{F}}-u\Vert }}.$$

$$\frac{dy}{dx}+y(x)={\int }_0^x{e}^{t-x}y(t)dt,$$

$${\int }_0^x{e}^{t-x}y(t)dt\approx \sum _{i=1}^n{w}_i{e}^{t_i(x)-x}y(t_i(x)),$$

$$\frac{dy}{dx}+y(x)\approx \sum _{i=1}^n{w}_i{e}^{t_i(x)-x}y(t_i(x)).$$

$$\mathcal{I}(u,\mathbf{x})=0,\text{for}\mathbf{x}\in {\mathcal{T}}_i.$$

$$\mathcal{L}(\mathbf{\theta },\mathbf{\lambda };\mathcal{T})=w_f{\mathcal{L}}_f(\mathbf{\theta },\mathbf{\lambda };{\mathcal{T}}_f)+w_b{\mathcal{L}}_b(\mathbf{\theta },\mathbf{\lambda };{\mathcal{T}}_b)+w_i{\mathcal{L}}_i(\mathbf{\theta },\mathbf{\lambda };{\mathcal{T}}_i),$$

$${\mathcal{L}}_i(\mathbf{\theta },\mathbf{\lambda };{\mathcal{T}}_i)=\frac{1}{|{\mathcal{T}}_i|}\sum _{\mathbf{x}\in {\mathcal{T}}_i}\Vert \mathcal{I}(\hat{u},\mathbf{x}){\Vert }_2^2.$$

$${\mathcal{E}}_r=\frac{1}{V}{\int }_{\Omega }\Vert f\left(\mathbf{x};\frac{\partial \hat{u}}{\partial x_1},\dots ,\frac{\partial \hat{u}}{\partial x_d};\frac{{\partial }^2\hat{u}}{\partial x_1\partial x_1},\dots ,\frac{{\partial }^2\hat{u}}{\partial x_1\partial x_d};\dots ;\mathbf{\lambda }\right)\Vert d\mathbf{x}$$

$${\mathcal{E}}_r\approx \frac{1}{|\mathcal{S}|}\sum _{\mathbf{x}\in \mathcal{S}}\Vert f\left(\mathbf{x};\frac{\partial \hat{u}}{\partial x_1},\dots ,\frac{\partial \hat{u}}{\partial x_d};\frac{{\partial }^2\hat{u}}{\partial x_1\partial x_1},\dots ,\frac{{\partial }^2\hat{u}}{\partial x_1\partial x_d};\dots ;\mathbf{\lambda }\right)\Vert .$$

$$-\Delta u(x,y)$$

$$u(x,y)$$

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=\nu \frac{{\partial }^{2}u}{\partial x^2},x\in [-1,1],t\in [0,1],$$

$$u(x,0)=-\sin (\pi x),u(-1,t)=u(1,t)=0.$$

$${\partial }_{t}u+u{\partial }_{x}u+v{\partial }_{y}u=\frac{1}{\text{mbox}Re}({\partial }_x^{2}u+{\partial }_y^{2}u),$$

$${\partial }_{t}v+u{\partial }_{x}v+v{\partial }_{y}v=\frac{1}{\text{mbox}Re}({\partial }_x^{2}v+{\partial }_y^{2}v),$$

$$x,y\in [0,1],\text{mbox}andt\in [0,1],$$

$$u(x,y,t)=\frac{3}{4}-\frac{1}{4\left[1+\exp ((-4x+4y-t)\text{mbox}Re/32)\right]},$$

$$v(x,y,t)=\frac{3}{4}+\frac{1}{4\left[1+\exp ((-4x+4y-t)\text{mbox}Re/32)\right]},$$

$$\frac{dx}{dt}=\rho (y-x),\frac{dy}{dt}=x(\sigma -z)-y,\frac{dz}{dt}=xy-\beta z,$$

$$\frac{\partial C_A}{\partial t}=D\frac{{\partial }^2{C}_A}{\partial x^2}-k_f{C}_A{C}_B^2,\frac{\partial C_B}{\partial t}=D\frac{{\partial }^2{C}_B}{\partial x^2}-2k_f{C}_A{C}_B^2,x\in [0,1],t\in [0,10],$$

Full text

{AMS}

65-01, 65-04, 65L99, 65M99, 65N99

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 DeepXDE111Source code is published under the Apache License, Version 2.0 on GitHub. https://github.com/lululxvi/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 (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 Fig. 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 feed-forward 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 $\mathcal{N}^{L}(\mathbf{x}):\mathbb{R}^{d_{\text{in}}}\to\mathbb{R}^{d_{\text{out}}}$ be a $L$-layer neural network, or a $(L-1)$-hidden layer neural network, with $N_{\ell}$ neurons in the $\ell$-th layer ($N_{0}=d_{\text{in}}$, $N_{L}=d_{\text{out}}$). Let us denote the weight matrix and bias vector in the $\ell$-th layer by $\bm{W}^{\ell}\in\mathbb{R}^{N_{\ell}\times N_{\ell-1}}$ and $\mathbf{b}^{\ell}\in\mathbb{R}^{N_{\ell}}$, respectively. Given a nonlinear activation function $\sigma$, which is applied element-wisely, the FNN is recursively defined as follows:

[TABLE]

see also a visualization of a neural network in Fig. 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\}$).

2.2 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$:

[TABLE]

The forward pass and backward pass of AD for computing the partial derivatives $\frac{\partial y}{\partial x_{1}}$ and $\frac{\partial y}{\partial x_{2}}$ 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 $\frac{\partial y}{\partial x_{i}}$ requires two function valuations $y(x_{1},\dots,x_{i},\dots,x_{d_{\text{in}}})$ and $y(x_{1},\dots,x_{i}+\Delta x_{i},\dots,x_{d_{\text{in}}})$ for some small number $\Delta x_{i}$, and thus in total $d_{\text{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].

2.3 Physics-informed neural networks (PINNs) for solving PDEs

We consider the following PDE parameterized by $\bm{\lambda}$ for the solution $u(\mathbf{x})$ with $\mathbf{x}=(x_{1},\dots,x_{d})$ defined on a domain $\Omega\subset\mathbb{R}^{d}$:

[TABLE]

with suitable boundary conditions

[TABLE]

where $\mathcal{B}(u,\mathbf{x})$ could be Dirichlet, Neumann, Robin, or periodic boundary conditions. For time-dependent problems, we consider time $t$ as a special component of $\mathbf{x}$, and $\Omega$ contains the temporal domain. The initial condition can be simply treated as a special type of Dirichlet boundary condition on the spatio-temporal domain.

The algorithm of PINN [31, 47] is shown in Procedure 1, and visually in the schematic of Fig. 1 solving a diffusion equation $\frac{\partial u}{\partial t}=\lambda\frac{\partial^{2}u}{\partial x^{2}}$ with mixed boundary conditions $u(x,t)=g_{D}(x,t)$ on $\Gamma_{D}\subset\partial\Omega$ and $\frac{\partial u}{\partial\mathbf{n}}(x,t)=g_{R}(u,x,t)$ on $\Gamma_{R}\subset\partial\Omega$. We explain each step as follows. In a PINN, we first construct a neural network $\hat{u}(\mathbf{x};\bm{\theta})$ as a surrogate of the solution $u(\mathbf{x})$, which takes the input $\mathbf{x}$ and outputs a vector with the same dimension as $u$. Here, $\bm{\theta}=\{\bm{W}^{\ell},\bm{b}^{\ell}\}_{1\leq\ell\leq L}$ is the set of all weight matrices and bias vectors in the neural network $\hat{u}$. One advantage of PINNs by choosing neural networks as the surrogate of $u$ is that we can take the derivatives of $\hat{u}$ with respect to its input $\mathbf{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 $\hat{u}$ to satisfy the physics imposed by the PDE and boundary conditions. In practice, we restrict $\hat{u}$ on some scattered points (e.g., randomly distributed points, or clustered points in the domain [37]), i.e., the training data $\mathcal{T}=\{\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{|\mathcal{T}|}\}$ of size $|\mathcal{T}|$. In addition, $\mathcal{T}$ is comprised of two sets $\mathcal{T}_{f}\subset\Omega$ and $\mathcal{T}_{b}\subset\partial\Omega$, which are the points in the domain and on the boundary, respectively. We refer $\mathcal{T}_{f}$ and $\mathcal{T}_{b}$ as the sets of “residual points”.

To measure the discrepancy between the neural network $\hat{u}$ 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:

[TABLE]

where

[TABLE]

and $w_{f}$ and $w_{b}$ are the weights. The loss involves derivatives, such as the partial derivative $\partial\hat{u}/\partial x_{1}$ or the normal derivative at the boundary $\partial\hat{u}/\partial\mathbf{n}=\nabla\hat{u}\cdot\mathbf{n}$, which are handled via AD.

In the last step, the procedure of searching for a good $\bm{\theta}$ by minimizing the loss $\mathcal{L}(\bm{\theta};\mathcal{T})$ is called “training”. Considering the fact that the loss is highly nonlinear and non-convex with respect to $\bm{\theta}$ [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 second-order 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 $\mathcal{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 $\Omega=[0,1]$, we can simply choose the surrogate model as $\hat{u}(x)=x(x-1)\mathcal{N}(x)$ to satisfy the boundary condition automatically, where $\mathcal{N}(x)$ is a neural network.

We note that we have great flexibility in choosing the residual points $\mathcal{T}$, and here we provide three possible strategies:

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. 2.

In each optimization iteration, we could select randomly different residual points. 3. 3.

We could improve the location of the residual points adaptively during training, e.g., the method proposed in Section 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 $\mathcal{T}=\Omega$ with $|\mathcal{T}|=\infty$.

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)=\sum_{k=1}^{5}\sin(2kx)/(2k)$ in $[-\pi,\pi]$ by a NN, the function is learned from low to high frequency (Fig. 2A). However, when we employ a PINN to solve the Poisson equation $-f_{xx}=\sum_{k=1}^{5}2k\sin(2kx)$ with zero boundary conditions in the same domain, all frequencies are learned almost simultaneously (Fig. 2B). Interestingly, by comparing Fig. 2A and Fig. 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 $\mathbb{Z}_{+}^{d}$ be the set of $d$-dimensional nonnegative integers. For $\mathbf{m}=(m_{1},\dots,m_{d})\in\mathbb{Z}_{+}^{d}$, we set $|\mathbf{m}|\coloneqq m_{1}+\dots+m_{d}$, and

[TABLE]

We say $f\in C^{\mathbf{m}}(\mathbb{R}^{d})$ if $D^{\mathbf{k}}f\in C(\mathbb{R}^{d})$ for all $\mathbf{k}\leq\mathbf{m}$, $\mathbf{k}\in\mathbb{Z}_{+}^{d}$, where $C(\mathbb{R}^{d})=\{f:\mathbb{R}^{d}\to\mathbb{R}|f\text{ 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].

Theorem 2.1**.**

Let $\mathbf{m}^{i}\in\mathbb{Z}_{+}^{d}$, $i=1,\dots,s$, and set $m=\max_{i=1,\dots,s}|\mathbf{m}^{i}|$. Assume $\sigma\in C^{m}(\mathbb{R})$ and $\sigma$ is not a polynomial. Then the space of single hidden layer neural nets

[TABLE]

is dense in

[TABLE]

i.e., for any $f\in C^{\mathbf{m}^{1},\dots,\mathbf{m}^{s}}(\mathbb{R}^{d})$, any compact $K\subset\mathbb{R}^{d}$, and any $\varepsilon>0$, there exists a $g\in\mathcal{M}(\sigma)$ satisfying

[TABLE]

*for all $\mathbf{k}\in\mathbb{Z}_{+}^{d}$ for which $\mathbf{k}\leq\mathbf{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 $\mathcal{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 $\mathcal{F}$, and we define $u_{\mathcal{F}}=\arg\min_{f\in\mathcal{F}}\|f-u\|$ as the best function in $\mathcal{F}$ close to $u$ (Fig. 3). Because we only train the neural network on the training set $\mathcal{T}$, we define $u_{\mathcal{T}}=\arg\min_{f\in\mathcal{F}}\mathcal{L}(f;\mathcal{T})$ as the neural network whose loss is at global minimum. For simplicity, we assume that $u$, $u_{\mathcal{F}}$ and $u_{\mathcal{T}}$ are well defined and unique. Finding $u_{\mathcal{T}}$ by minimizing the loss is often computationally intractable [10], and our optimizer returns an approximate solution $\tilde{u}_{\mathcal{T}}$.

We can then decompose the total error $\mathcal{E}$ as [11]

[TABLE]

The approximation error $\mathcal{E}_{\text{app}}$ measures how closely $u_{\mathcal{F}}$ can approximate $u$. The generalization error $\mathcal{E}_{\text{gen}}$ is determined by the number/locations of residual points in $\mathcal{T}$ and the capacity of the family $\mathcal{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 $\mathcal{E}_{\text{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].

2.5 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.

2.6 PINNs for solving integro-differential equations

When solving integro-differential 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 (Fig. 4) [42], such as Gaussian quadrature. Therefore, we introduce a fourth error component, the discretization error $\mathcal{E}_{\text{dis}}$, due to the approximation of the integral by Gaussian quadrature.

For example, when solving

[TABLE]

we first use Gaussian quadrature of degree $n$ to approximate the integral

[TABLE]

and then we use a PINN to solve the following PDE instead of the original equation

[TABLE]

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.

2.7 PINNs for solving inverse problems

In inverse problems, there are some unknown parameters $\bm{\lambda}$ in Eq. 1, but we have some extra information on some points $\mathcal{T}_{i}\subset\Omega$ besides the differential equation and boundary conditions:

[TABLE]

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 Eq. 2:

[TABLE]

where

[TABLE]

We then optimize $\bm{\theta}$ and $\bm{\lambda}$ together, and our solution is $\bm{\theta}^{*},\bm{\lambda}^{*}=\arg\min_{\bm{\theta},\bm{\lambda}}\mathcal{L}(\bm{\theta},\bm{\lambda};\mathcal{T})$.

2.8 Residual-based adaptive refinement (RAR)

As we discussed in Section 2.3, the residual points $\mathcal{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), 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 $\left\|f\left(\mathbf{x};\frac{\partial\hat{u}}{\partial x_{1}},\dots,\frac{\partial\hat{u}}{\partial x_{d}};\frac{\partial^{2}\hat{u}}{\partial x_{1}\partial x_{1}},\dots,\frac{\partial^{2}\hat{u}}{\partial x_{1}\partial x_{d}};\dots;\bm{\lambda}\right)\right\|$ is large, and we repeat adding points until the mean residual

[TABLE]

is smaller than a threshold $\mathcal{E}_{0}$, where $V$ is the volume of $\Omega$.

3 DeepXDE usage and customization

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

3.1 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 and Fig. 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 (|), difference (-) and intersection (&). This technique is called constructive solid geometry (CSG), see Fig. 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 4. 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.

3.2.2 Neural networks

DeepXDE currently supports two neural networks: feed-forward neural network (maps.FNN) and residual neural network (maps.ResNet). A new network can be added as shown in Procedure 5.

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 6.

4 Demonstration examples

In this section, we use PINNs and DeepXDE to solve different problems. In all examples, we use the $\tanh$ as the activation function, and the other hyperparameters are listed in Table 3. The weights $w_{f}$, $w_{b}$ and $w_{i}$ in the loss function are set as 1. The codes of all examples are published in GitHub.

4.1 Poisson equation over an L-shaped domain

Consider the following two-dimensional Poisson equation over an L-shaped domain $\Omega=[-1,1]^{2}\setminus[0,1]^{2}$:

[TABLE]

We choose 1200 and 120 random points drawn from a uniform distribution as $\mathcal{T}_{f}$ and $\mathcal{T}_{b}$, respectively. The PINN solution is given in Fig. 7B. For comparison, we also present the numerical solution obtained by using the spectral element method (SEM) [27] (Fig. 7A). The result of the absolute error is shown in Fig. 7C.

4.2 RAR for Burgers equation

We first consider the 1D Burgers equation:

[TABLE]

Let $\nu=0.01/\pi$. 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 Section 2.8 with $m=1$ and $\mathcal{E}_{0}=0.005$. We compare the PINN solution with RAR and the PINN solution based on 2540 randomly selected training data (Fig. 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 Fig. 8A. Here, we also present two examples for the residual points added via the RAR method. As shown in Fig. 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:

[TABLE]

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 $=UL/\nu$, in which $U$ and $L$ are respectively the characteristic velocity and length, and $\nu$ is the kinematic viscosity of fluid. The exact solution can be obtained [3] as

[TABLE]

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 Fig. 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

[TABLE]

with the initial condition $(x(0),y(0),z(0))=(-8,7,27)$, where $\rho$, $\sigma$ and $\beta$ 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 $(\rho,\sigma,\beta)=(10,15,8/3)$. We choose 400 uniformly distributed random points and 25 equispaced points as the residual points $\mathcal{T}_{f}$ and $\mathcal{T}_{i}$, respectively. The evolution trajectories of $\rho$, $\sigma$ and $\beta$ are presented in Fig. 9A, with the final identified values of $(\rho,\sigma,\beta)=(10.002,14.999,2.668)$.

4.4 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\rightarrow C$) is described by

[TABLE]

where $D=2\times 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\times 10^{-3}$) and $k_{f}$ (0.0971) are displayed in Fig. 9B, which agree well with their true values.

4.5 Volterra IDE

Here, we consider the first-order integro-differential equation of the Volterra type in the domain $[0,5]$:

[TABLE]

with the exact solution $y(x)=e^{-x}\cosh x.$ We solve this IDE using the method in Section 2.6, and approximate the integral using Gaussian-Legendre quadrature of degree 20. The $L^{2}$ relative error is $2\times 10^{-3}$, and the solution is shown in Fig. 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.

Acknowledgments

This work is supported by the DOE PhILMs project (No. de-sc0019453), the AFOSR grant FA9550-17-1-0013, and the DARPA-AIRA grant HR00111990025.