Solving Partial Differential Equations with Point Source Based on Physics-Informed Neural Networks

Here we introduce two well-known time-(in)dependent PDEs with point source forcing terms. The Poisson’s equation with a point source has the following form:

	$\displaystyle-\Delta u$	$\displaystyle=\delta(\mathbf{x}-\mathbf{x}_{0}),\quad\mathbf{x}\in\Omega,$		(2)
	$\displaystyle u$	$\displaystyle=0,\qquad\qquad\,\,\,\,\mathbf{x}\in\partial\Omega$

where $\delta$ denotes the Dirac delta function. Although it looks simple, Poisson’s equation is a fundamental model for describing diffusion process, such as fluid flow, heat transfer, and chemical transport. In the field of electromagnetic simulations, the time domain Maxwell’s equations with excitation can be expressed as

	$\displaystyle\nabla\times E$	$\displaystyle=-\mu\dfrac{\partial H}{\partial t}-J(\mathbf{x},t),$		(3)
	$\displaystyle\nabla\times H$	$\displaystyle=\epsilon\dfrac{\partial E}{\partial t}.$

where $\epsilon$ is the permittivity, $\mu$ is the permeability. The forcing term $J(\mathbf{x},t)=\delta(\mathbf{x}-\mathbf{x}_{0})g(t)$ excites a wave propagation at $\mathbf{x}_{0}$ with a particular type of pulse $g(t)$. Maxwell’s equations survived the revolutionary changes in physics in the 20th century and they underpin a huge range of natural phenomena.

In order to solve PDE problems with deep learning technology the PINNs method is proposed to approximate the solution $u(\mathbf{x},t)$ with a deep neural network $u(\mathbf{x},t;\theta)$, where $\theta$ represents the trainable parameters of the deep neural network. The loss function given by such approximation is defined as

	$\displaystyle L_{total}(\theta)$	$\displaystyle=L_{r}(\theta)+\lambda_{ic}L_{ic}(\theta)+\lambda_{bc}L_{bc}(\theta),$		(4a)
	$\displaystyle L_{r}(\theta)$	$\displaystyle=\dfrac{1}{N_{r}}\sum_{i=1}^{N_{r}}||\mathcal{N}(u(\mathbf{x}_{i},t_{i};\theta))-f(\mathbf{x}_{i},t_{i})||_{2}^{2},$		(4b)
	$\displaystyle L_{ic}(\theta)$	$\displaystyle=\dfrac{1}{N_{ic}}\sum_{i=1}^{N_{ic}}||u(\mathbf{x}_{i},0;\theta)-h(\mathbf{x}_{i})||_{2}^{2},$		(4c)
	$\displaystyle L_{bc}(\theta)$	$\displaystyle=\dfrac{1}{N_{bc}}\sum_{i=1}^{N_{bc}}||u(\mathbf{x}_{i},t_{i};\theta)-g(\mathbf{x}_{i},t_{i})||_{2}^{2}.$		(4d)

where $L_{r}$, $L_{ic}$, $L_{bc}$ correspond to the loss term of PDE residual, initial conditions, and boundary conditions, respectively. The hyperparameters $\lambda_{ic}$ and $\lambda_{bc}$ aim to balance the interplay of different terms in the loss function. $N_{r}$, $N_{ic}$, and $N_{bc}$ represent the number of samples in the entire region, initial conditions, and boundary conditions, respectively.

When $f(\cdot),g(\cdot),h(\cdot)$ are continuous, the network weights $\theta$ can be optimized by minimizing the total training loss $L_{total}(\theta)$ via standard gradient descent procedures used in deep learning. Unfortunately, due to the singularity brought by $\delta(\mathbf{x})$, the PINNs method cannot be directly used to solve the PDEs with a point source.

One can use the Deep Ritz method [12] to solve Eq.(2). Through mathematical derivation, Eq.(2) can be converted to the minimum value problem:

$\displaystyle\min\int_{\Omega}{\frac{1}{2}\|\Delta u(\mathbf{x})\|^{2}d\mathbf{x}-u(\mathbf{x}_{0})}+\int_{\partial\Omega}{\|u(\mathbf{x})\|^{2}}d\mathbf{x}$

(5)

where the Dirac function term $\delta(\mathbf{x}-\mathbf{x}_{0})$ in Eq.(2) will generate the constant term $u(\mathbf{x}_{0})$ in Eq.(5) after integration. However, application of the Deep Ritz method is restricted because it’s unable to handle PDEs like Eq.(3) which cannot be derived from any variational problems. The NVIDIA SimNet [19] proposes to solve PDEs with a point source using the variational formulation. However, the performance heavily depends on the selection of the test functions and the computational time increases linearly with the number of test functions.

Zhang et al. [16] and Moseley et al. [18] use the classical numerical method to simulate the early period ($0\leq t\leq T_{0}<T$) of wave equations to avoid the point source singularity, and then applied the PINNs method to solve the PDEs with given initial conditions at $t=T_{0}$. Bekele et al. [20] attempted to leverage the PINNs loss by adding a new loss term with a significant amount of known labels when solving the Barry and Mercer’s source problem with time-dependent fluid injection. However, this method relies on classical computational methods to generate the labeled data which is expensive or even infeasible in many situations.

In this work, we propose a universal solution based on the PINNs method which does not rely on any labeled data or variational forms. Furthermore, our method greatly improves the accuracy and convergency speed compared with current deep learning based approaches.

The Dirac delta function $\delta(x)$ is mathematically defined by:

$$\delta(x)=\begin{cases}+\infty,&x=0\\ 0,&x\neq 0\end{cases}$$

(6)

and which is also constrained to satisfy the identity $\int_{-\infty}^{+\infty}\delta(x)dx=1$. Taking $\eta(x)$ to be any symmetric unimodal continuous probability density function centered at 0, $\eta(x)=\alpha^{-1}\eta(\dfrac{x}{\alpha})$ can be a natural approximation to $\delta(x)$ as the kernel width $\alpha\rightarrow 0$. For $d$-dimensions space, $\eta_{\alpha}(\mathbf{x})=\alpha^{-d}\eta(\dfrac{\mathbf{x}}{\alpha})$ can be used instead. In our experiments, we use Gaussian distribution, Cauchy distribution, or Laplacian distribution with sufficiently small kernel width $\alpha$ to approximate $\delta(\mathbf{x})$.

Our approximation to $\delta(x)$ eliminates function discontinuity at the origin, but it also brings a difficulty for the PINNs method to optimize its training loss. As the probability density function is highly concentrated at the origin as $\alpha\rightarrow 0$, the residual samples close to the origin tend to dominate the whole residual loss, making it difficult for the network to learn the governing equations and therefore the standard optimizer may not converge to the solution. To address this problem, we split the feasible region into two subdomains such that:

$$\Omega=\Omega_{0}\cup\Omega_{1}\quad\Omega_{0}\cap\Omega_{1}=\emptyset$$

where $\Omega_{0}$ represents the subdomain that contains the origin and $\Omega_{1}$ covers the complement region. By taking this decomposition, the residual loss $L_{r}$ in Eq.(4) is divided into two parts:

$$L_{r}(\theta)=\lambda_{r,0}L_{r,0}(\theta)+\lambda_{r,1}L_{r,1}(\theta)$$

(7)

where $L_{r,0}(\theta)$ and $L_{r,1}(\theta)$ correspond to the residual losses in $\Omega_{0}$ and $\Omega_{1}$, respectively. In our experiment, for the point source excited at $\mathbf{x}_{0}\in\Omega$, we split $\Omega$ as

$$\Omega_{0}=\{\mathbf{x}_{0}+\mathbf{x}\in\Omega\,|\,\|\mathbf{x}\|\leq 3\alpha\},\quad\Omega_{1}=\Omega\backslash\Omega_{0}.$$

Together with the IC&BC conditions, the total PINNs loss in Eq.(4a) changes to be:

$\displaystyle L_{total}(\theta)=\lambda_{r,0}L_{r,0}(\theta)+\lambda_{r,1}L_{r,1}(\theta)+\lambda_{ic}L_{ic}(\theta)+\lambda_{bc}L_{bc}(\theta).$

(8)

Properly setting the weight vector $\lambda=(\lambda_{r,0},\lambda_{r,1},\lambda_{ic},\lambda_{bc})$ is critical to enhance the trainability of constrained neural networks [22, 23], but searching the optimal weight vector through manual tuning is infeasible. We adapt the uncertainty weighing method in [24] to solve the problem. For a total loss function including $m$ loss terms $L(\theta)=\sum_{i=1}^{m}\lambda_{i}L_{i}(\theta)$, we add a nonnegative $\epsilon^{2}$ to $w^{2}$, and get the total PINNs loss function as follows:

$$L_{total}(\theta;\mathbf{w})=\sum_{i=1}^{m}\dfrac{1}{2(\epsilon^{2}+w_{i}^{2})}L_{i}(\theta)+\log(\epsilon^{2}+w_{i}^{2})$$

(9)

The uncertainty of the $i$-th loss term is approximated by $\epsilon^{2}+w_{i}^{2}$, whose lower bound is constrainted by $\epsilon^{2}$ when decreasing $w_{i}$ to 0. We argue that such reformulation is nontrivial since it provides a lower bound for the uncertainty estimation of each loss term and thus increases the stability of the self-adaptive weighting algorithm. In our experiment, taking $\epsilon=0.01$ leads to robust model outputs and outperforms the original formulation ($\epsilon=0$).

Fig.1 shows the $L_{2}\ error$s when different $\epsilon$ is used in solving Maxwell’s equations with a point source (i.e., Eq.(3)). For detailed experimental settings, please refer to Sec.4.1. When $\epsilon=0.01$ and $\epsilon=0.001$, the lower bound constrained uncertainty weighting method outperforms the original uncertainty weighting method ($\epsilon=0$) with not only lower final $L_{2}\ error$s, but also faster convergence speed.

When $\Omega=(0,\pi)^{2}$ and the point source is located at $(x_{0},y_{0})$, the analytical solution of the 2-D Poisson’s equation with a point source (i.e., Eq.(2)) is listed as follows:

$$u(x,y,x_{0},y_{0})=\frac{4}{\pi^{2}}\sum_{m,n=1}^{\infty}\frac{\sin{mx}\sin{mx_{0}}\sin{ny}\sin{ny_{0}}}{m^{2}+n^{2}}$$

(10)

The analytical solution is formed by superposition of multiple frequency components. This multi-frequency phenomenon is common in solutions of many PDEs. Recent development in deep learning theories [25, 26, 21] reveals that the neural network fits the hidden physical law from low to high frequency components, which differentiates from the classical iterative numerical algorithms (e.g., Gauss-Seidel method). Liu et al. [27] proposes multi-scale deep neural networks (MscaleDNNs) based on this theory which exhibit faster convergence in higher frequency components.

Periodic nonlinearities have been investigated repeatedly over the past decades. Recently, the SIREN architecture [28] proposes to use Sine function as a periodic activation function and outperforms alternative activation functions when dealing with implicit neural representations and their derivatives. Compared with other activation functions, e.g., Tanh, ReLU, the derivative of a Sine function is Cosine, a phase-shifted Sine. Therefore, the derivatives of a SIREN inherit the properties of SIREN, thus enabling us to supervise any derivative of SIREN with complex natural signals from PDEs. However, the SIREN architecture may yield poor accuracy when applied to PDEs with a point source whose solution are formed by superposition of multiple frequency components.

Borrowing ideas from both SIREN and MscaleDNNs, we propose Multi-Scale SIREN (MS-SIREN) that combines a multi-scale DNN with Sine activation function, as shown in Fig.2. Skip connections are added between consecutive hidden layers to accelerate training and improve accuracy, which make the neural network much easier to train since they help to avoid the vanishing gradient problem [29]. Fig.3 shows the effectiveness of MS-SIREN compared with other alternative architectures.

Maxwell’s equations are a set of coupled partial differential equations related to the fundamental physical modeling of electromagnetism. We use 2-D Maxwell’s equations to validate the robustness of our method. The governing equations of TE wave are given as:

	$\displaystyle\frac{\partial E_{x}}{\partial t}$	$\displaystyle=\frac{1}{\epsilon}\frac{\partial H_{z}}{\partial y},$		(12a)
	$\displaystyle\frac{\partial E_{y}}{\partial t}$	$\displaystyle=-\frac{1}{\epsilon}\frac{\partial H_{z}}{\partial x},$		(12b)
	$\displaystyle\frac{\partial H_{z}}{\partial t}$	$\displaystyle=-\frac{1}{\mu}(\frac{\partial E_{y}}{\partial x}-\frac{\partial E_{x}}{\partial y}+J).$		(12c)

where $E_{x}$ and $E_{y}$ are the electric field in the Cartesian coordinate, and $H_{z}$ is the magnetic field perpendicular to the horizontal plane. The constants $\mu\approx 4\pi\times 10^{-7}H/m$ and $\epsilon\approx 8.854\times 10^{-12}F/m$ are known as the permeability and permittivity of the free space. $J$ is a known source function that represents the energy passing through the source node. Without loss of generality, a Gaussian pulse containing multiple frequency components is used as the source function that is expressed as:

$\displaystyle J(x,y,t)=e^{-(\frac{t-d}{\tau})^{2}}\delta(x-x_{0})\delta(y-y_{0}).$

(13)

where $d$ is the temporal delay and $\tau$ is a pulse-width parameter. Currently, we take $\tau=3.65\times\sqrt{2.3}/(\pi f)$ with the characteristic frequency $f$ set to $1GHz$. Suppose the initial electromagnetic field is static, which corresponds to the initial conditions as follows:

	$\displaystyle E_{x}(x,y,0)$	$\displaystyle=E_{y}(x,y,0)=0,$		(14a)
	$\displaystyle H_{z}(x,y,0)$	$\displaystyle=0.$		(14b)

To solve the equations uniquely over a finite vacuum domain, the standard Mur’s second-order absorbing boundary condition is utilized. For a rectangular truncated domain, boundary conditions for the left, right, upper and lower sides, can be expressed as follows:

	$\displaystyle\frac{\partial H_{z}}{\partial x}-\frac{1}{c}\frac{\partial H_{z}}{\partial t}+\frac{c\epsilon}{2}\frac{\partial E_{x}}{\partial y}$	$\displaystyle=0,$		(15a)
	$\displaystyle\frac{\partial H_{z}}{\partial x}+\frac{1}{c}\frac{\partial H_{z}}{\partial t}-\frac{c\epsilon}{2}\frac{\partial E_{x}}{\partial y}$	$\displaystyle=0,$		(15b)
	$\displaystyle\frac{\partial H_{z}}{\partial y}-\frac{1}{c}\frac{\partial H_{z}}{\partial t}+\frac{c\epsilon}{2}\frac{\partial E_{y}}{\partial x}$	$\displaystyle=0,$		(15c)
	$\displaystyle\frac{\partial H_{z}}{\partial y}+\frac{1}{c}\frac{\partial H_{z}}{\partial t}-\frac{c\epsilon}{2}\frac{\partial E_{y}}{\partial x}$	$\displaystyle=0.$		(15d)

where $c=1/\sqrt{\mu\epsilon}$ is the light speed.

We solve the problem in a rectangular domain $\Omega=[-1,1]^{2}$, and the point source is located at the origin $(x,y)=(0,0)$. The total simulation time is set to $4ns$ so that the pulse signals can propagate throughout the whole space of truncated domain. The reference solution is obtained through the finite-difference time-domain (FDTD) method [30] at a spatial resolution of $\Delta=0.005$. For model training, the spatial points are sampled from uniform distribution so that the average spacing between two adjacent random points is comparable with that in the reference FDTD mesh, and the average temporal step is chosen to satisfy the $CFL$ convergence condition [30]. The spatial point source is approximated by the Gaussian distribution with kernel width $\alpha=2\Delta=0.01$. The model is trained using Adam optimizer with an initial learning rate of 0.001, and the learning rate attenuates 10 times when the training process reaches 40%, 60%, and 80%, respectively. The total number of iterations is $200$k, and the batch size is $N_{r,0}=N_{r,1}=N_{ic}=N_{bc}=8192$ for each iteration. The instantaneous electromagnetic fields at the time of $2.4ns$ are showed in Fig.4. We can see that the predicted electromagnetic fields by our method are almost indistinguishable from the reference solution.

The MS-SIREN architecture uses skip connections between consecutive hidden layers. To validate its effectiveness, we compare MS-SIREN architectures with and without skip connections, and show the results in Fig.6. It is clear that skip connections make $L_{2}\ error$ decrease faster and reach a lower value.

Activation function is another key factor affecting the prediction accuracy. To validate the effectiveness of Sine periodic activation function, we compare it with ReLU and Tanh, and the results are shown in Fig.7. Sine produces the best results in terms of convergence speed and the final mean $L_{2}\ error$. In comparison, Tanh performs slightly worse than Sine, and ReLU fails to let the network learn any useful information.

We apply our approach to solve Eq.(2) in the region of $\Omega=[0,\pi]\times[0,\pi]$ with the point source located at $\mathbf{x}_{0}=(\frac{\pi}{2},\frac{\pi}{2})$. The point source is approximated using Gaussian distribution with fixed $\alpha=0.01$ and the lower bound hyperparameter $\epsilon$ is set to 0.01. The model is trained for $50k$ iterations using Adam optimizer with an initial learning rate of 0.001, and the learning rate attenuates 10 times when the training process reaches 40%, 60%, and 80%, respectively. In each iteration, the batch size is set to $(N_{r_{,}0},N_{r,1},N_{bc})=(8192,8192,8192)$. Fig.8 shows that the absolute error between our solution and the analytical solution is small.

For comparisons, we also apply the Deep Ritz method and the variational method in SimNet to solve the same problem. In this set of experiments, we implement our method and Deep Ritz method in MindSpore framework, and implement the variational method using the public code in SimNet with a default selection of 15 test functions. All the experiments are run on an NVIDIA Tesla V100 GPU card. Fig.9(a) shows that the convergence speed of our method is close to that of Deep Ritz method, but faster than that of SimNet. Fig.9(b) shows the final $L_{2}\ error$s obtained by different methods, and obviously our method achieves the highest accuracy.

The third problem we solve is the nondimensionalized poroelastic Barry and Mercer’s source problem with time-dependent fluid injection. The governing equations are given by

	$\displaystyle\frac{\partial^{2}u}{\partial t\partial x}+\frac{\partial^{2}v}{\partial t\partial z}-\frac{\partial^{2}p}{\partial x^{2}}-\frac{\partial^{2}p}{\partial z^{2}}-\beta Q$	$\displaystyle=0,$		(16a)
	$\displaystyle(\eta+1)\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial z^{2}}+\eta\frac{\partial^{2}v}{\partial x\partial z}-(\eta+1)\frac{\partial p}{\partial x}$	$\displaystyle=0,$		(16b)
	$\displaystyle\frac{\partial^{2}v}{\partial x^{2}}+(\eta+1)\frac{\partial^{2}v}{\partial z^{2}}+\eta\frac{\partial^{2}u}{\partial x\partial z}-(\eta+1)\frac{\partial p}{\partial z}$	$\displaystyle=0.$		(16c)

where $u$ and $v$ are the deformations of porous medium along the x and z directions, $p$ is the pore fluid pressure, $Q$ is a fluid source or sink term, $\eta$ and $\beta$ are nondimensional parameters which are functions of the real material parameters.

An illustration of the problem is shown in Fig.10. An oscillating point source $Q$ located at $(x_{0},z_{0})$ in the rectangular domain $[0,a]\times[0,b]$ is given by:

$$Q(x,z,t)=\delta(x-x_{0})\delta(z-z_{0})sin(\omega t)$$

(17)

where $\omega$ is the oscillation frequency. The parameters are set as follows: $a=b=1$, $t\in[0,2\pi]$, $\beta=2$, $\eta=1.5$, $(x_{0},z_{0})=(0.25,0.25)$ and $\omega=1$. The model is trained for $250k$ iterations using Adam optimizer with an initial learning rate of 0.001, and the learning rate attenuates 10 times when the training process reaches 40%, 60%, and 80%, respectively. In each iteration, the batch size is set to $(N_{r_{,}0},N_{r,1},N_{bc},N_{ic})=(8192,8192,8192,8192)$.

Fig.11 shows the solution at $t=\pi/5$ obtained by our method (middle) and the analytical solution (top), and the absolute error between them (bottom). The time averaged mean $L_{2}\ error$ of the predicted solution is 0.017, demonstrating the high accuracy of our method. In addition to Gaussian distribution, we also use Cauchy distribution and Laplacian distribution to approximate the Dirac delta function, and obtain similar accuracy with $L_{2}\ error$ being 0.051 and 0.025, respectively.

Fig.13 compares the convergence speed of the mean $L_{2}\ error$s when different activation functions are used. The periodic activation function Sine exhibits superior performance over ReLU and Tanh similar to the case of Maxwell’s equations.

In summary, approximating the Dirac function with the symmetric unimodal continuous probability density function is necessary for all experiments in order for PINNs to solve the singularity problem posed by the point source. For the effectiveness of the MS-SIREN architecture, the multi-scale network exhibits better accuracy and faster convergence speed than the single-scale network in both Maxwell’s equations and Poisson’s equation experiments. However, the muti-scale network does not bring significant gain in Barry and Mercer’s source problem due to the fact that the PDE solutions do not contain multiple frequency components. Moreover, the Sine activation function exhibits better accuracy than ReLU and Tanh in all three sets of experiments. For the Maxwell’s equations and the Barry and Mercer’s source problem, the lower bound constrained uncertainty weighting algorithm can obtain a lower $L_{2}\ error$ relative to the original uncertainty weighing method. It is worth noting that all three of our proposed techniques need to be utilized to obtain the desired accuracy when solving the Maxwell’s equations.