Physics Informed Extreme Learning Machine (PIELM) -- A rapid method for the numerical solution of partial differential equations

TL;DR

PIELM is a fast, physics-informed machine learning method that efficiently solves linear partial differential equations, matching or surpassing PINNs in accuracy and offering a scalable distributed version for large domains.

Contribution

This paper introduces PIELM, a rapid physics-informed neural network alternative for PDEs, and proposes DPIELM, a distributed extension for large-scale problems.

Findings

PIELM achieves comparable or better accuracy than PINNs.

DPIELM provides effective solutions for large domain PDEs.

Neural network approaches can be competitive with traditional methods.

Abstract

There has been rapid progress recently on the application of deep networks to the solution of partial differential equations, collectively labelled as Physics Informed Neural Networks (PINNs). In this paper, we develop Physics Informed Extreme Learning Machine (PIELM), a rapid version of PINNs which can be applied to stationary and time dependent linear partial differential equations. We demonstrate that PIELM matches or exceeds the accuracy of PINNs on a range of problems. We also discuss the limitations of neural network based approaches, including our PIELM, in the solution of PDEs on large domains and suggest an extension, a distributed version of our algorithm -{}- DPIELM. We show that DPIELM produces excellent results comparable to conventional numerical techniques in the solution of time-dependent problems. Collectively, this work contributes towards making the use of neural networks in the solution of partial differential equations in complex domains as a competitive alternative to conventional discretization techniques.

Tables

Table 1. Table 1: List of test cases for PIELM.

Steady/ Unsteady	1D/2D	Test case ID	Description
Steady	1D	TC-1	Linear advection
		TC-2	Linear diffusion
		TC-3	Linear advection-diffusion
	2D	TC-4	Linear advection in a star shaped computational domain
		TC-5	Linear diffusion in a star shaped computational domain
		TC-6	Linear diffusion in a complex computational domain
Unsteady	1D	TC-7	Linear advection
		TC-8	Quasi-linear advection
		TC-9	Linear advection-diffusion
	2D	TC-10	Linear advection-diffusion

Table 2. Table 2: Summary of experiments for 1D steady test cases.

TC	$[N_f,N_{bc},N^{\ast }]$	$𝒪\left(Error\right)$
TC-1	$[40,2,42]$	${10}^{-4}$
TC-2	$[40,2,42]$	${10}^{-4}$
TC-3	$[20,2,22]$	${10}^{-6}$

Table 3. Table 3: Summary of experiments for 2D steady test cases.

TC	$[N_f,N_{BC},N^{\ast }]$	$𝒪\left(Error\right)$
TC-4	$[921,240,2000]$	${10}^{-6}$
TC-5	$[921,240,2000]$	${10}^{-4}$
TC-6	$[1489,881,5000]$	${10}^{-7}$

Table 4. Table 4: List of numerical experiments that show limitations of PIELM and PINN

	Test case ID	Description
Representation of functions with PIELM	TC-11	Representation of sharp and discontinuous functions in 1D
	TC-12	Representation of a sharp peaked 2D Gaussian
Solution of PDEs with PIELM	TC-13	TC-7 with $F=e^{-100x^2}$
	TC-14	TC-3 with $\nu =0.02$
Solution of PDEs with PINN	TC-15	TC-7 with $F=e^{-5x^2}sin\left(10\pi x\right)$

Table 5. Table 5: Details of DPIELM architecture for the test cases. For 1D steady problems, architecture is given by [ N ​ B x , n ​ b x , N c ​ e ​ l ​ l ∗ ] 𝑁 subscript 𝐵 𝑥 𝑛 subscript 𝑏 𝑥 superscript subscript 𝑁 𝑐 𝑒 𝑙 𝑙 [NB_{x},nb_{x},N_{cell}^{*}] , where N ​ B x 𝑁 subscript 𝐵 𝑥 NB_{x} , n ​ b x 𝑛 subscript 𝑏 𝑥 nb_{x} and N c ​ e ​ l ​ l ∗ superscript subscript 𝑁 𝑐 𝑒 𝑙 𝑙 N_{cell}^{*} refer to number of cells, number of points in the cell and size of hidden layer of the PIELM. Similarly, for 1D and 2D unsteady problems, it is given by [ N ​ B x , N ​ B t , n ​ b x , n ​ b t , N c ​ e ​ l ​ l ∗ ] 𝑁 subscript 𝐵 𝑥 𝑁 subscript 𝐵 𝑡 𝑛 subscript 𝑏 𝑥 𝑛 subscript 𝑏 𝑡 superscript subscript 𝑁 𝑐 𝑒 𝑙 𝑙 [NB_{x},NB_{t},nb_{x},nb_{t},N_{cell}^{*}] and [ N ​ B x , N ​ B y , N ​ B t , n ​ b x , n ​ b y , n ​ b t , N c ​ e ​ l ​ l ∗ ] 𝑁 subscript 𝐵 𝑥 𝑁 subscript 𝐵 𝑦 𝑁 subscript 𝐵 𝑡 𝑛 subscript 𝑏 𝑥 𝑛 subscript 𝑏 𝑦 𝑛 subscript 𝑏 𝑡 superscript subscript 𝑁 𝑐 𝑒 𝑙 𝑙 [NB_{x},NB_{y},NB_{t},nb_{x},nb_{y},nb_{t},N_{cell}^{*}] respectively.

Test Case ID	Description	Architecture
TC-9	1D unsteady linear advection-diffusion	$[10,10,5,5,30]$
TC-10	2D unsteady linear advection-diffusion equation	$[20,20,50,3,3,3,30]$
TC-11	Representation of 1D sharp and discontinuous functions	$[50,5,5]$
TC-12	Representation of a sharp peaked 2D Gaussian with DPIELM	$[15,15,5,5,15]$
TC-13	Advection of a sharp peaked Gaussian	$[15,10,5,5,30]$
TC-14	1D steady advection-diffusion	$[20,5,20]$
TC-15	Advection of a high frequency wave packet	$[15,10,5,5,30]$

Equations

$${\overrightarrow{y}}_{ELM}^{(i)}={\mathbf{H}}^{(i)}(\overrightarrow{x})\overrightarrow{c}$$

$${\mathbf{H}}^{(i)}=[{\overrightarrow{h}}^{(i)}({\overrightarrow{x}}_1),{\overrightarrow{h}}^{(i)}({\overrightarrow{x}}_2),...,{\overrightarrow{h}}^{(i)}({\overrightarrow{x}}_N){]}^T,$$

$${\overrightarrow{h}}^{(i)}({\overrightarrow{x}}_k)=[\phi ({\overrightarrow{x}}_k;{\overrightarrow{a}}_1^{(i)},b_1^{(i)}),\phi (\overrightarrow{x};{\overrightarrow{a}}_2^{(i)},b_2^{(i)}),...,\phi (\overrightarrow{x};{\overrightarrow{a}}_{N^{\ast }}^{(i)},b_{N^{\ast }}^{(i)})],$$

$${\mathbf{H}}^{(i)}\overrightarrow{c}={\overrightarrow{y}}^{(i)}-{\overrightarrow{\xi }}^{(i)},i=1,2,...m,$$

$$J=\frac{1}{2}||\overrightarrow{c}|{|}^2+\frac{1}{2N}\lambda \sum _{i=1}^m{\overrightarrow{\xi }}^{(i{)}^T}{\overrightarrow{\xi }}^{(i)},$$

$$\frac{\partial J}{\partial c_k}=0,k=1,2,...,N^{\ast }$$

$$\frac{\partial }{\partial t}u(\overrightarrow{x},t)+\mathcal{N}u(\overrightarrow{x},t)=R(\overrightarrow{x},t),(\overrightarrow{x},t)ϵ\Omega \text{x}[0,T],$$

$$u(\overrightarrow{x},t)=B(\overrightarrow{x},t),(\overrightarrow{x},t)ϵ\partial \Omega \text{x}[0,T],$$

$$u(\overrightarrow{x},0)=F(\overrightarrow{x}),\overrightarrow{x}ϵ\Omega ,$$

$${\overrightarrow{\xi }}_f=\frac{\partial \overrightarrow{f}}{\partial t}+\mathcal{N}\overrightarrow{f}-\overrightarrow{R},(\overrightarrow{x},t)ϵ\Omega \text{x}[0,T],$$

$${\overrightarrow{\xi }}_{bc}=\overrightarrow{f}-\overrightarrow{B,}(\overrightarrow{x},t)ϵ\partial \Omega \text{x}[0,T],$$

$${\overrightarrow{\xi }}_{ic}=\overrightarrow{f}(.,0)-\overrightarrow{F},\overrightarrow{x}ϵ\Omega .$$

$$J=\frac{{\overrightarrow{\xi }}_f^T{\overrightarrow{\xi }}_f}{2N_f}+\frac{{\overrightarrow{\xi }}_{bc}^T{\overrightarrow{\xi }}_{bc}}{2N_{bc}}+\frac{{\overrightarrow{\xi }}_{ic}^T{\overrightarrow{\xi }}_{ic}}{2N_{ic}},$$

$$\frac{\partial }{\partial t}u(\overrightarrow{x},t)+\mathcal{L}u(\overrightarrow{x},t)=R(\overrightarrow{x},t),(\overrightarrow{x},t)ϵ\Omega \text{x}[0,T],$$

$$u(\overrightarrow{x},t)=B(\overrightarrow{x},t),(\overrightarrow{x},t)ϵ\partial \Omega \text{x}[0,T],$$

$$u(\overrightarrow{x},0)=F(\overrightarrow{x}),\overrightarrow{x}ϵ\Omega ,$$

$$h_k=\phi (z_k),$$

$$f(\overrightarrow{\chi })=\overrightarrow{h}\overrightarrow{c}.$$

$$\frac{{\partial }^p{f}_k}{\partial x^p}=m_k^p\frac{{\partial }^p\phi }{\partial z^p},$$

$$\frac{\partial f_k}{\partial t}=n_k\frac{\partial \phi }{\partial z}.$$

$${\overrightarrow{\xi }}_f=\frac{\partial \overrightarrow{f}}{\partial t}+\mathcal{L}\overrightarrow{f}-\overrightarrow{R},(\overrightarrow{x},t)ϵ\Omega \text{x}[0,T],$$

$${\overrightarrow{\xi }}_{bc}=\overrightarrow{f}-\overrightarrow{B,}(\overrightarrow{x},t)ϵ\partial \Omega \text{x}[0,T],$$

$${\overrightarrow{\xi }}_{ic}=\overrightarrow{f}(.,0)-\overrightarrow{F},\overrightarrow{x}ϵ\Omega .$$

$${\overrightarrow{\xi }}_f=\overrightarrow{0},$$

$${\overrightarrow{\xi }}_{bc}=\overrightarrow{0},$$

$${\overrightarrow{\xi }}_{ic}=\overrightarrow{0}.$$

$$\mathbf{H}\overrightarrow{c}=\overrightarrow{K}.$$

$$J=\frac{1}{2}||\overrightarrow{c}|{|}^2+\frac{1}{2}\lambda \left(\frac{{\overrightarrow{\xi }}_f^T{\overrightarrow{\xi }}_f}{2N_f}+\frac{{\overrightarrow{\xi }}_{bc}^T{\overrightarrow{\xi }}_{bc}}{2N_{bc}}+\frac{{\overrightarrow{\xi }}_{ic}^T{\overrightarrow{\xi }}_{ic}}{2N_{ic}}\right),$$

$$\frac{\partial J}{\partial c_k}=0,k=1,2,...,N^{\ast }$$

$$u_x=R,0<x\le 1,$$

Full text

K****eywords Partial differential equations  $\cdot$ Physics informed neural networks  $\cdot$ Extreme learning machine  $\cdot$ Advection-diffusion equation

1 Introduction

Partial differential equations (PDEs) are extensively used in the mathematical modelling of various problems in physics, engineering and finance. In practical situations, these equations typically lack analytical solutions and are solved numerically. In current practice, most numerical approaches to solve PDEs like finite element method (FEM), finite difference method (FDM) and finite volume method (FVM) are mesh based. A typical implementation of a mesh based approach involves three steps: (1) Grid generation, (2) Discretization of governing equation and (3) Solution of the discretized equations with some iterative method.

However, there are limitations to these approaches. Some of the limitations of these methods are as follows:

They cannot be used to solve PDEs in complex computational domains because grid generation (step 1) itself becomes infeasible. 2. 2.

The process of discretization (step 2) creates a discrepancy between the mathematical nature of actual PDE and its approximate difference equation [1]. Sometimes this can lead to quite serious problems [2].

One of the options to fix these issues is to use neural networks. In this approach, the data set consists of some randomly selected points in the domain and on the boundary. The governing equations and the boundary conditions are fitted using neural network. There are two main motivations for this approach. First, being universal approximators, neural networks can potentially represent any PDE. So this avoids the discretization step and thus discretization based physics errors too. Second, it is meshfree and therefore complex geometries can be easily handled [3]. Initial work in this direction can be credited to Lagaris et al. [4, 5]. Firstly, they solved the initial boundary value problem using neural networks and later they extended their work to handle irregular boundaries. Since then, a lot of work has been done in this field [6, 7, 8, 9, 10, 11, 12, 13]. In particular, we refer to the physics-informed neural networks (PINN) approach by Raissi and Karniadakis [11] and Raissi et. al [12, 13]. This approach has produced promising results for a series of benchmark nonlinear problems.

Recently, Berg et al. [3] have developed a PINN based method to solve PDEs on complex domains and produced several results. However, in spite of various advantages of using deep networks for solving PDEs, PINNs have several problems [13]. Firstly, there is no theoretical basis to know the size of neural network architecture and the amount of data needed. Then, there is no guarantee that the algorithm will not hit upon a local minima. Finally, their learning time is slower than the traditional numerical methods making them very expensive for practical problems.

We show that some of the problems mentioned above can be easily handled by using an alternative network called the extreme learning machine (ELM). The basic ELM was proposed by Huang et al. [14] for a single hidden layer feed forward networks (SLFNs) and later it was extended to generalized SLFNs. The essence of ELM is that the weights of the hidden layer of SLFNs need not to be learnt. It is much faster than the traditional gradient based optimization methods alleviating the learning time problem. Previously, ELMs have been used in approximating functions [15] and solving ordinary differential equations (ODEs) and stationary PDEs [16, 17] using Legendre and Bernstein polynomial basis functions respectively.

In this paper, we propose a new machine learning algorithm to solve stationary and time dependent PDEs in complex geometries. We have named it physics informed extreme learning machine (PIELM) because it is a combination of two algorithms namely ELM and PINN. Theoretically, there is no question over the employment of ELM as a PDE solver because it is a universal approximator [18] and therefore it can approximate any PDE. We have made our ELM “physics informed” by incorporating the information about the physics of PDE as the cost function. In addition to this, we have also proposed an extension to original PIELM called distributed PIELM that enhances the representation power of PIELM without adding any extra hidden layers. We demonstrate that both PIELM and DPIELM exibit superior performance on a range of stationary and time-dependent problems in comparison to existing methods.

This paper proceeds as follows. We give a brief review of PINN and ELM in Section 2. The proposed PIELM is described in Section 3.In Section 4, we evaluate the performance of PIELM in solving various stationary and time-dependent PDEs. To our knowledge, this is the first application of an ELM based algorithm to solve a 2D unsteady PDE. In Section 5, we discuss the limitation of PIELM to represent discontinuous functions and the functions with sharp gradient. We also illustrate a test case where even a deep PINN fails to represent a complicated function. In Section 6, DPIELM, the distributed version of PIELM for enhanced representation is described. In Section 7, the results of implementation of DPIELM algorithm in test cases involving representation of functions with sharp gradients have been presented. We have also shown that DPIELM outperforms the deep PINN in representing the complicated function described earlier. Finally, conclusion and future work are given in Section 8.

2 Brief review of ELM and PINN

The PIELM is combination of two learning algorithms: ELM and PINN. In this section, we review these two algorithms in brief.

2.1 Extreme learning machine

Traditional gradient-based learning algorithms [19] have prohibitively slow learning speed and they suffer from various problems like improper learning rate, local minima etc. Huang et al. [14] originally proposed a novel learning algorithm called ELM to fix these issues. ELM is extremely fast and mostly it shows better generalization performance than gradient-based learning approaches like back-propagation. A typical implementation of the algorithm involves the following steps:

Select a shallow neural network. 2. 2.

Fix the hidden layer weights and biases randomly. These parameters will not be learned and therefore no iterative optimization is required for them. 3. 3.

Apply a nonlinear transformation to the input data set. This gives the input to the final layer. 4. 4.

Take the linear combination of all the inputs of the final layer. This is the output of ELM. 5. 5.

Learn the output layer weights using the least squares method.

Mathematical formulation

Consider the basic ELM shown in Fig (1). It is a single layer feed forward neural network with $N^{*}$ neurons in the hidden layer. Input is a vector of size $n$ and output is the $i^{th}$ component of the output vector of size $m$. We denote the non-linear activation by $\varphi$ and the weights and biases of $j^{th}$ node of hidden layer by $\overrightarrow{a}_{j}^{(i)}$ and $b_{j}^{(i)}$respectively. The output of the hidden node $j$ is $\varphi(\overrightarrow{x};\overrightarrow{a}_{j}^{(i)},b_{j}^{(i)})$. For a given data set $\{(x_{k},y_{k})\}_{k=1}^{N}\subset\Re^{n}\text{x}\Re^{m}$ with $N$ distinct samples, the ELM output is given by

[TABLE]

where,

[TABLE]

[TABLE]

and $\overrightarrow{c}=[c_{1},c_{2},...,c_{N^{*}}]^{T}$ is vector of output layer weights.

On writing the Eq. (1) for all the $m$ components, the resulting ELM is given by

[TABLE]

where $\overrightarrow{\xi}$ is the training error vector. The ELM tends to reach the smallest training error together with the smallest norm of the output weights. Mathematically saying, the loss function to be minimized for the ELM is given by

[TABLE]

where $\lambda$ is the regularization parameter. The correct weights that minimize $J$ can be calculated by solving the normal equations as given below.

[TABLE]

2.2 Physics informed neural network

Raissi et al. [13] proposed a data efficient PINN for approximating solutions to general non-linear PDEs and validated it with a series of benchmark test cases. The main feature of the PINN is the inclusion of the prior knowledge of physics in the learning algorithm as cost function. As a result, the algorithm imposes penalty for any non-physical solution and quickly directs it towards the correct solution. This physics informed approach enhances the information content of the data. As a result, the algorithm has good generalization property even in the small data set regime.

Mathematical formulation

Consider a PDE of the following form

[TABLE]

[TABLE]

[TABLE]

where $\mathcal{\mathscr{N}}$ may be a linear or nonlinear differential operator and $\partial\varOmega$ is the boundary of computational domain $\Omega$. We approximate $u(\overrightarrow{x},t)$ with the output $f(\overrightarrow{x},t)$ of PINN. The network architecture may be shallow or deep depending upon the non-linearity $\mathcal{\mathscr{N}}$. The essence of PINN lies in the definition of its loss function. In order to make the neural network “physics informed”, the loss function is defined such that a penalty is imposed whenever the network output doesn’t respect the physics of the problem. If we denote the training errors in approximating the PDE, BCs and IC by $\overrightarrow{\xi}_{f}$, $\overrightarrow{\xi}_{bc}$ and $\overrightarrow{\xi}_{ic}$ respectively. Then, the expressions for these errors are as follows:

[TABLE]

[TABLE]

[TABLE]

For shallow networks, $\frac{\partial\overrightarrow{f}}{\partial t}$ and $\mathcal{\mathcal{\mathscr{N}}}\overrightarrow{f}$ can be determined using hand calculations. However, for deep networks, we have to use automatic differentiation [20]. The loss function $J$ to be minimized for a PINN is given by

[TABLE]

where $N_{f}$, $N_{bc}$ and $N_{ic}$ refer to number of collocation points, boundary condition points and initial condition points respectively. Finally, any gradient based optimization routine may be used to minimize $J$.

This completes the mathematical formulation of PINN. The key steps in its implementation are as follows:

Identify the PDE to be solved along with the initial and boundary conditions. 2. 2.

Decide the architecture of PINN. 3. 3.

Approximate the correct solution with PINN. 4. 4.

Find expressions for the PDE, BCs and IC in terms of PINN and its derivatives. 5. 5.

Define a loss function which penalizes for error in PDE, BCs and IC. 6. 6.

Minimize the loss with gradient based algorithms.

3 Proposed PIELM

Consider the following unsteady linear PDE

[TABLE]

[TABLE]

[TABLE]

where $\mathcal{L}$ is a linear differential operator and $\partial\varOmega$ is the boundary of computational domain $\Omega$. We approximate $u(\overrightarrow{x},t)$ with the output $f(\overrightarrow{x},t)$ of PIELM. For simplicity, we consider the 1D unsteady version of Eqns (12, 13, 14). The extension to higher dimensional problems is straightforward. The PIELM for 1D unsteady problem is schematically shown in Fig (2). The number of neurons in the hidden layers is $N^{*}$. If we define $\overrightarrow{\chi}=[x,t,1]^{T}$, $\overrightarrow{m}=[m_{1},m_{2},...,m_{N^{*}}]^{T}$, $\overrightarrow{n}=[n_{1},n_{2},...,n_{N^{*}}]^{T}$, $\overrightarrow{b}=[b_{1},b_{2},...,b_{N^{*}}]^{T}$ and $\overrightarrow{c}=[c_{1},c_{2},...,c_{N^{*}}]^{T}$ then, the output of the $k^{th}$ hidden neuron is

[TABLE]

where $z_{k}=[m_{k},n_{k},b_{k}]\overrightarrow{\chi}$ and $\varphi=tanh$ is the nonlinear activation function. The PIELM output is given by

[TABLE]

Similarly, the formulae for $\frac{\partial^{p}f}{\partial x^{p}}$ and $\frac{\partial f}{\partial t}$ are given by

[TABLE]

[TABLE]

We denote the training errors in approximating the PDE, BCs and IC by $\overrightarrow{\xi}_{f}$, $\overrightarrow{\xi}_{bc}$ and $\overrightarrow{\xi}_{ic}$ respectively. The expressions for these errors are as follows:

[TABLE]

[TABLE]

[TABLE]

Next, we put a hard constraint on $\overrightarrow{c}$ to solve the PDE exactly with zero error by setting

[TABLE]

[TABLE]

[TABLE]

Eqns (21 to 23) lead to the a system of linear equations which can be represented as

[TABLE]

The form of $\boldsymbol{H}$ and $\overrightarrow{K}$ depends on the $\mathcal{L}$, $B$ and $F$ i.e. on the type of PDE, boundary condition and initial condition. In order to find $\overrightarrow{c}$, Moore–Penrose generalized inverse [21] (also called pseudo-inverse) should be used as it works well for singular and non square $\boldsymbol{H}$ too. An additional advantage with this formulation is that we have a basis to guess the scale of the PIELM architecture. When we are solving Eqns (21 to 23) simultaneously, we know that a unique solution will exist when number of unknowns are equal to number of equations which means that $N^{*}=N_{f}+N_{bc}+N_{ic}$. This gives us an idea of the size of the hidden layer. However, in practice we get the correct solution even with lesser number of neurons. For example, if we supply a large number of points to approximate a linear function, the PIELM would not require the same number of neurons for learning.

This completes the mathematical formulation of PIELM. The key steps in its implementation are as follows:

Assign the input layer weights randomly. 2. 2.

Depending on the PDE and the initial and boundary conditions, find the expressions for $\overrightarrow{\xi}_{f}$, $\overrightarrow{\xi}_{bc}$ and $\overrightarrow{\xi}_{ic}$. 3. 3.

Assemble the three sets of equations in the form of $\boldsymbol{H}\overrightarrow{c}=\overrightarrow{K}.$ 4. 4.

Output layer weight vector is given by $pinv(\boldsymbol{H})\overrightarrow{K},$ where $pinv$ refers to pseudo-inverse.

It is to be noted that unlike conventional ELMs , we are not solving an optimization problem. The loss function $J$ to be minimized for a conventional physics informed ELM would be given by

[TABLE]

where $\lambda$ is a regularization parameter and $N_{f}$, $N_{bc}$ and $N_{ic}$ refer to number of collocation points, boundary condition points and initial condition points respectively. The correct ELM weights that minimize $J$ can be calculated by solving the normal equations which is given below.

[TABLE]

Although a PIELM can be made with this minimization approach, we have opted for the direct approach due to the following reasons:

The direct approach is straightforward to formulate and code. It saves the effort of calculating loss function and setting the derivatives equal to zero. 2. 2.

The learning of the minimization approach is comparatively less “physics informed” because physics is not being imposed in an exact sense.

4 Performance evaluation of PIELM

To evaluate the performance of PIELM, we rigorously test it on various linear and quasi-linear PDEs described in Table(1). TC-1, TC-2, TC-4, TC-5, TC-6 are taken from Berg et al. [3]. TC-8 is taken from Kopriva et al. [22]. TC-9 and TC-10 are taken from Borker et al. [23]. All the experiments are conducted in Matlab 2017b environment running in an Intel Core i5 2.20GHz CPU and 8GB RAM Dell laptop. The error is defined as the difference between the PIELM prediction and the exact solution.

4.1 1D steady cases [ TC-1, TC-2, TC-3 ]

The 1D stationary advection, diffusion and advection-diffusion equations are given by

[TABLE]

[TABLE]

[TABLE]

respectively. The expressions for $R$ and the Dirichlet boundary conditions for these cases are calculated by assuming the following exact solutions

[TABLE]

[TABLE]

[TABLE]

respectively. In order to solve these equations in PIELM framework, we have to solve Eqn(21 to 23). The expression for $\overrightarrow{\xi}_{f}$ depends on the linear differential operator $\mathcal{L}$. The definitions of $\mathcal{L}$ in these cases are $\frac{\partial}{\partial x}$, $\frac{\partial^{2}}{\partial x^{2}}$ and $\frac{\partial}{\partial x}-\nu\frac{\partial^{2}}{\partial x^{2}}$ respectively. When $\mathcal{L}$ acts on $u$, the corresponding expressions for $\overrightarrow{\xi}_{f}=\overrightarrow{0}$ may be written as follows:

[TABLE]

[TABLE]

[TABLE]

where $\overrightarrow{x}_{f}$ is collocation points vector, $\overrightarrow{I}$ is bias vector, $\boldsymbol{X_{f}}=[\overrightarrow{x}_{f},\overrightarrow{I}]$ and ’$\odot$’ refers to Hadamardt product. Referring to Fig (3a), $\boldsymbol{W}=[\overrightarrow{m},\overrightarrow{b}]$. Similarly, expression for $\overrightarrow{\xi}_{bc}=\overrightarrow{0}$ is given by

[TABLE]

where $\overrightarrow{x}_{bc}$ is boundary points vector, $\boldsymbol{X_{bc}}=[\overrightarrow{x}_{bc},\overrightarrow{I}]$ and $B$ is the boundary condition.

The results for these test cases are given in Fig(4), Fig(5) and Fig (6) respectively and the summary of the experiments is given in Tab(2).

Remark

It should be noted that the unified deep ANN algorithm [3] took 100 points to achieve an order of accuracy of $10^{-5}$ and $10^{-3}$ in TC-1 and TC-2 respectively. In comparison, PIELM took less than half of the points and still achieved an order of accuracy of $10^{-4}$ in both the cases.

4.2 2D steady cases [ TC-4, TC-5, TC-6 ]

The stationary 2D advection and diffusion equations for the three cases are given by

[TABLE]

[TABLE]

[TABLE]

where $a=1$ and $b=\frac{1}{2}$ are advection coefficients. The computational domains $\varOmega_{1}$ and $\varOmega_{2}$ are shown in respectively. The expressions for $R$ and the $B$ for these cases are constructed by choosing the following exact solutions

[TABLE]

[TABLE]

[TABLE]

respectively. PIELM equations for these problems are as follows:

$\overrightarrow{\xi}_{f}$=$\overrightarrow{0}$

• •

Case 1: Advection

[TABLE]

• •

Case 2: Diffusion

[TABLE]

where $\boldsymbol{X_{f}}=[\overrightarrow{x}_{f},\overrightarrow{y}_{f},\overrightarrow{I}]$ and $\boldsymbol{W}=[\overrightarrow{m},\overrightarrow{n},\overrightarrow{b}]$ (refer Fig (2)). 2. 2.

$\overrightarrow{\xi}_{bc}=\overrightarrow{0}$

[TABLE]

where $\boldsymbol{X_{bc}}=[\overrightarrow{x}_{bc},\overrightarrow{y}_{bc},\overrightarrow{I}]$.

The results for the advection and diffusion cases on $\varOmega_{1}$ are shown in Fig (8) and Fig (9) respectively. The result for diffusion case on $\varOmega_{2}$ is given in Fig (10). Summary of the experiments is given in

Remarks

The unified deep ANN algorithm [3] took $5500$ points to achieve an order of accuracy of $10^{-3}$ in TC-4 and TC-5. In comparison, PIELM solved these cases with merely $1161$ points and still achieved an order of accuracy of $10^{-6}$ and $10^{-4}$ respectively. 2. 2.

False diffusion is an error which gives diffusion like appearance to solution of pure advection equation when it is solved using upwind discretization. In Fig (8) these errors can be seen flowing along streamlines. However, PIELM reduces the order of these errors to an insignificant level. 3. 3.

Computational domain $\Omega_{2}$ is the map of state of Illinois in USA. This is an example of a complicated polygon which has very short line segments and fine grained details in various regions. Conventional mesh based methods are not feasible for these kind of geometries. The latitudes and longitudes of the boundary are available in MATLAB’s in-built function “usamap”. We have re-scaled the data in the range $0-1$. PIELM solved this case with just $2370$ points to an accuracy level of $10^{-7}$ in just $10$ seconds.

4.3 1D unsteady advection cases [ TC-7, TC-8 ]

The unsteady 1D advection equation is given by

[TABLE]

[TABLE]

where $a(x,t)$ is the advection coefficient. In this problem, we consider two cases: (1) constant coefficient case with $a(x,t)=1$ and (2) variable coefficient case [22] with $a(x,t)=1+x$. The two cases have periodic and inflow boundary conditions respectively. The value of $F$ is $sin(\pi x)$. The exact solutions to the 1D unsteady advection problems with constant and variable coefficient are respectively given by

[TABLE]

[TABLE]

PIELM equations

$\overrightarrow{\xi}_{f}$=$\overrightarrow{0}$

• •

TC-7: Linear case

[TABLE]

• •

TC-8: Quasi-linear case

[TABLE]

where $\overrightarrow{x}_{f},\overrightarrow{t}_{f}$ are collocation point vectors, $\boldsymbol{X_{f}}=[\overrightarrow{x}_{f},\overrightarrow{t}_{f},\overrightarrow{I}]$, $\boldsymbol{W}=[\overrightarrow{m},\overrightarrow{n},\overrightarrow{b}]$ and \boldsymbol{A_{f}}=[\begin{array}[]{cccc}a\left(\overrightarrow{x_{f}},\overrightarrow{t_{f}}\right),&...&,&a\left(\overrightarrow{x_{f}},\overrightarrow{t_{f}}\right)\end{array}]_{N\text{x}N^{*}}. 2. 2.

$\overrightarrow{\xi}_{bc}=\overrightarrow{0}$

[TABLE]

where $\overrightarrow{x}_{lbc}$, $\overrightarrow{t}_{lbc}$ are left boundary points vectors, $\overrightarrow{x}_{rbc}$, $\overrightarrow{t}_{rbc}$ are right boundary points vectors, $\boldsymbol{X_{lbc}}=[\overrightarrow{x}_{lbc},\overrightarrow{t}_{lbc},\overrightarrow{I}]$ and $\boldsymbol{X_{rbc}}=[\overrightarrow{x}_{rbc},\overrightarrow{t}_{rbc},\overrightarrow{I}]$. 3. 3.

$\overrightarrow{\xi}_{ic}=\overrightarrow{0}$

[TABLE]

where$\overrightarrow{x}_{ic},\overrightarrow{t}_{ic}$ are initial condition vectors, $\boldsymbol{X_{ic}}=[\overrightarrow{x}_{ic},\overrightarrow{t}_{ic},\overrightarrow{I}]$ and $F$ is initial condition.

The results are shown in Fig (11). PIELM predicts the exact solution correctly for both linear and quasi-linear advection. In this case, we took $N_{f}=420$, $N_{bc}=21$, $N_{ic}=20$ and $N^{*}=440$. The total learning time is within 2-3 seconds.

Remark

For advection problems, the time step of the traditional numerical schemes like upwinding can not exceed the mesh size due to stability issues. However, PIELM doesn’t impose any such restriction and we can take larger time steps.

4.3.1 1D unsteady advection-diffusion [ TC-9 ]

The 1D equivalent of the unsteady 2D advection-diffusion equation solved by Borker et al. [23] is given by

[TABLE]

where $a$ is the advection coefficient. The expressions for the initial condition $F$ and the boundary condition $B$ are constructed on the basis of the following exact solution

[TABLE]

The PIELM equations are as follows:

$\overrightarrow{\xi}_{f}$=$\overrightarrow{0}$

[TABLE]

where $\overrightarrow{x}_{f},\overrightarrow{t}_{f}$ are collocation point vectors, $\boldsymbol{X_{f}}=[\overrightarrow{x}_{f},\overrightarrow{t}_{f},\overrightarrow{I}]$ and $\boldsymbol{W}=[\overrightarrow{m},\overrightarrow{n},\overrightarrow{b}]$. 2. 2.

$\overrightarrow{\xi}_{bc}=\overrightarrow{0}$

[TABLE]

where $\overrightarrow{x}_{bc}$, $\overrightarrow{t}_{bc}$ are boundary points vectors, $\boldsymbol{X_{bc}}=[\overrightarrow{x}_{bc},\overrightarrow{t}_{bc},\overrightarrow{I}]$ and $B$ is the boundary condition. 3. 3.

$\overrightarrow{\xi}_{ic}=\overrightarrow{0}$

[TABLE]

where$\overrightarrow{x}_{ic},\overrightarrow{t}_{ic}$ are initial condition vectors, $\boldsymbol{X_{ic}}=[\overrightarrow{x}_{ic},\overrightarrow{t}_{ic},\overrightarrow{I}]$ and $F$ is initial condition.

The results are shown in Fig (12). In this case, PIELM prediction clearly goes wrong in the following grounds:

Initial and boundary conditions aren’t captured correctly. 2. 2.

Solution has unphysical oscillations throughout the domain. 3. 3.

The hump of the Gaussian decays a lot slower than the correct rate.

We increased the size of hidden layer but didn’t see any improvement in results. For example, we took $250\text{x}30$ points in the computational domain and put as many as $7780$ neurons in the hidden layer. Still the PIELM predictions were poor.

4.3.2 2D unsteady advection-diffusion [ TC-10 ]

We solve the following unsteady 2D advection-diffusion equation[23]

[TABLE]

[TABLE]

[TABLE]

where $\Omega_{xy}=[0,1]\text{x}[0,1]$. $a=cos(22.5)$ and $b=sin(22.5)$ are advection coefficients in $x$ and $y$ directions and $\nu=0.005$ is the diffusion coefficient. The expressions for the initial and the boundary conditions are constructed by choosing the following exact solution:

[TABLE]

PIELM equations to be solved are as follows:

$\overrightarrow{\xi}_{f}$=$\overrightarrow{0}$

[TABLE]

where $\overrightarrow{x}_{f},\overrightarrow{y}_{f},\overrightarrow{t}_{f}$ are collocation point vectors, $\boldsymbol{X_{f}}=[\overrightarrow{x}_{f},\overrightarrow{y}_{f},\overrightarrow{t}_{f},\overrightarrow{I}]$. Referring to Fig (3b), $\boldsymbol{W}=[\overrightarrow{p},\overrightarrow{q},\overrightarrow{r},\overrightarrow{s}]$. 2. 2.

$\overrightarrow{\xi}_{bc}=\overrightarrow{0}$

[TABLE]

where $\overrightarrow{x}_{bc}$, $\overrightarrow{y}_{bc}$, $\overrightarrow{t}_{bc}$ are boundary points vectors and $\boldsymbol{X_{bc}}=[\overrightarrow{x}_{bc},\overrightarrow{y}_{bc},\overrightarrow{t}_{bc},\overrightarrow{I}]$. 3. 3.

$\overrightarrow{\xi}_{ic}=\overrightarrow{0}$

[TABLE]

where$\overrightarrow{x}_{ic},\overrightarrow{y}_{ic},\overrightarrow{t}_{ic}$ are initial condition vectors and $\boldsymbol{X_{ic}}=[\overrightarrow{x}_{ic},\overrightarrow{y}_{ic},\overrightarrow{t}_{ic},\overrightarrow{I}]$.

The results are shown in Fig (13-14). They are 2D equivalents of 1D results. The solution is diffusive and contains unphysical oscillations throughout the domain. In spite of all these issues, it should be noted that:

The advection component of the equation is captured correctly. The two humps are moving at the same speed. 2. 2.

The non-physical oscillations don’t grow with time.

For this case, we have taken a total of 125000 data points and 1000 neurons in the hidden layer. The results show signs of improvement if we keep increasing the number of points. However, that takes a lot of time which makes the option impractical. The limitations of the PIELM will be further investigated in the next section. We close this section by summarizing the advantages of the PIELM which are as follows:

It is extremely fast as well as data efficient.(TC-1 to TC-6) 2. 2.

It can be seamlessly extended to higher dimensions. (TC-10) 3. 3.

It reduces numerical artefacts like false diffusion.(TC-4) 4. 4.

It is meshfree method and can handle the complex geometries.(TC-6) 5. 5.

It reduces the arbitrariness of the number of neurons in the hidden layer.

5 Limitations of PIELM

A potential reason for the failure of PIELM in solving advection-diffusion equation [23] could be the limited representation capacity of PIELM to represent a complex function. A PDE consists of function and its derivatives. If a neural network cannot represent the function itself, then calculation of derivatives only adds to the error.

Representation of functions with sharp gradients [ TC-11, TC-12

]

To put the representation power of PIELM to test we choose the following two functions: (1) A 1D composite function that contains both discontinuous functions and functions with sharp gradients and (2) A 2D Gaussian function with a sharp peak. The expressions for 2D Gaussian function and 1D composite function are $f(x,y,t)=e^{-20\{(x-0.25)^{2}+(y-0.25)^{2}\}}$ and

[TABLE]

respectively. The results are shown in figure (15)and Fig (16) respectively. These cases clearly expose the limitation of PIELM in representing profiles with sharp gradients and corners.

PDEs with sharp solutions [ TC-13, TC-14 ]

Due to this limitation, PIELM fails to solve any PDE which admits functions with sharp gradients. For example, we have already seen the failure of our algorithm in solving 1D and 2D advection-diffusion equation. We further illustrate the impact of this limitation on two simpler equations. Firstly, we consider pure advection of a sharp peaked Gaussian. This equation has been solved in TC-7 for a smooth sine function. Next , we take steady advection-diffusion equation (TC-3) with a low value of diffusion coefficient.

The exact and PIELM solutions for these problems are shown in Fig (17) and Fig (18) respectively.

Representation of a high frequency wavelet with PINN [TC-15]

An obvious idea to improve the representation capacity of a neural network is to make it deep. Therefore, we test the performance of PINN in TC-7 with $F=e^{-5x^{2}}sin(10\pi x)$ i.e. pure advection of a high frequency wavelet. The PINN code is freely available at https://github.com/maziarraissi/PINNs. We modified the original code by replacing the Burgers equation with pure linear advection equation. Fig(19) shows the architecture of PINN. It consists of 9 hidden layers with 20 neurons each. Each layer is activated by $\tanh$ functions.

The exact and the PINN solutions are shown in Fig (20). We can see that even a deep network is unable to capture the sharp gradients of this wavelet.

List of the experiments conducted in this section are given in Tab 4. We summarize this section as follows:

The main limitation of a PIELM is its inability to represent complex functions. 2. 2.

This limitation restricts the PIELM to solve the PDEs with sharp exact solutions. 3. 3.

Adding extra layers is not the practical solution to this problem. (TC-15)

6 Distributed PIELM

In this section, we propose a distributed version of PIELM called DPIELM. This algorithm takes motivation from finite volume methods in which the whole computational domain is partitioned into multiple cells and governing equations are solved at each cell. The solutions of these individual cells are stitched together with additional convective and diffusive fluxes conditions at the cell interfaces. We adopt a similar strategy in DPIELM. As representation of a complex function is very hard for a single PIELM or PINN in the whole domain, we divide the domain into multiple cells and install a PIELM in each cell. Therefore, each PIELM uses different representations in different portions of domain while satisfying some additional constraints of continuity and differentiability.

6.1 Mathematical formulation

Consider the following 1D unsteady problem

[TABLE]

[TABLE]

[TABLE]

where $\mathcal{L}$ is a linear differential operator and $\partial\varOmega$ is the boundary of computational domain $\Omega$. In this problem, the rectangular domain $\Omega$ is given by $\Omega=[x_{L},x_{R}]\text{x}[0,T]$. On uniformly dividing $\Omega$ into $N_{c}$ non-overlapping rectangular cells, $\Omega$ may be written as

[TABLE]

The boundary of the cell $\Omega_{i}$ is denoted by $\partial\varOmega_{i}$. For rectangular cells,

[TABLE]

where $I_{m}^{(i)}$ represents the $m^{th}$ interface of $\Omega_{i}$.

Fig (21a) shows the distribution of PIELMs in a rectangular computational domain with $NB_{x}\text{x}NB_{t}$ cells (i.e. $N_{c}=NB_{x}\text{x}NB_{t}$). We denote the PIELM on the $i^{th}$ cell by $M^{(i)}$. Fig (21b) shows a PIELM with collocation points at the interior and the boundary points at the four interfaces. The weights and output corresponding to a given $M^{(i)}$ are denoted by \text{[\overrightarrow{m}}^{(i)},\text{\overrightarrow{n}}^{(i)},\text{\overrightarrow{b}}^{(i)},\text{\overrightarrow{c}}^{(i)}] and $f^{(i)}$respectively.

At each $M^{(i)}$, we enforce additional constraints of continuity (or smoothness) of solution at the cell interfaces depending on the differential operator $\mathcal{L}$. For example, continuity of solution is sufficient for advection problems. For diffusion problem, the solution should be continuously differentiable. For the computational domain shown in the Fig (21), the system of equations to be solved in the DPIELM framework are given below.

6.1.1 Regular PIELM equations

$\overrightarrow{\xi}_{f}^{(i)}$=$\overrightarrow{0}$

[TABLE] 2. 2.

$\overrightarrow{\xi}_{bc}^{(i)}=\overrightarrow{0}$

[TABLE]

[TABLE] 3. 3.

$\overrightarrow{\xi}_{ic}^{(i)}=\overrightarrow{0}$

[TABLE]

6.1.2 Additional interface equations

Constraints for $C^{0}$ solutions i.e. $\overrightarrow{\xi}_{C^{0}}^{(i)}=\overrightarrow{0}$.

• •

Continuity along $x$ direction

[TABLE]

where $\kappa=[1,(1+NB_{x}),...,1+(NB_{t}-1)NB_{x}]^{T}$.

• •

Continuity along $t$ direction

[TABLE]

where $\kappa=[1,2,...,NB_{x}]^{T}$. 2. 2.

Constraints for $C^{1}$ solutions i.e. $\overrightarrow{\xi}_{C^{1}}^{(i)}=\overrightarrow{0}$.

• •

Smooth solutions along $x$ direction

[TABLE]

where $\kappa=[1,(1+NB_{x}),...,1+(NB_{t}-1)NB_{x}]^{T}$.

Assembly of Eqns (72 to 78) leads to a system of linear equations which can be represented as

[TABLE]

where $\overrightarrow{c}=[\overrightarrow{c}^{(1)},\overrightarrow{c}^{(2)},...,\overrightarrow{c}^{(N_{c})}]^{T}.$ The form of $\boldsymbol{H}$ and $\overrightarrow{K}$ depends on the $\mathcal{L}$, $B$ and $F$. Finally $\overrightarrow{c}$ can be found using pseudo inverse. It is to be noted that although we have shown the formulation for 1D unsteady problems, no special adjustment is needed to extend the formulation to higher dimensional problems.

This completes the mathematical formulation of DPIELM. The main steps in its implementation are as follows:

Divide the computational domain into uniformly distributed non overlapping cells and install a PIELM in each cell. 2. 2.

Depending on the cell location, PDE and the initial and boundary conditions, find the expressions for $\overrightarrow{\xi}_{f}^{(i)}$, $\overrightarrow{\xi}_{bc}^{(i)}$ and $\overrightarrow{\xi}_{ic}^{(i)}$ at each cell. 3. 3.

Depending on the PDE, find the expressions for $\overrightarrow{\xi}_{C^{0}}^{(i)}$ and $\overrightarrow{\xi}_{C^{1}}^{(i)}$ at each cell interface. 4. 4.

Assemble these equations in the form of $\boldsymbol{H}\overrightarrow{c}=\overrightarrow{K},$ where $\overrightarrow{c}=[\overrightarrow{c}^{(1)},\overrightarrow{c}^{(2)},...,\overrightarrow{c}^{(N_{c})}]^{T}.$ 5. 5.

Find the value of $\overrightarrow{c}$ using pseudo inverse.

7 Performance evaluation of DPIELM

In this section, we evaluate the performance of DPIELMs by testing it on all the cases in which regular PIELM and PINN failed to perform. The details of the architecture is given in Tab(5).

Representation of functions with sharp gradient [ TC-11, TC-12]

The mathematical formulation of DPIELM equations for 1D case are as follows:

$\overrightarrow{\xi}_{f}^{(i)}$=$\overrightarrow{0}$

[TABLE] 2. 2.

$\overrightarrow{\xi}_{bc}^{(i)}=\overrightarrow{0}$

[TABLE]

where $I_{1},I_{2}$ refer to left and right cell interfaces respectively. 3. 3.

$\overrightarrow{\xi}_{C^{0}}^{(i)}=\overrightarrow{0}$

[TABLE]

The equations for the 3D case can be written in a similar fashion. The exact and DPIELM solutions for these cases are shown in Figs (22) and (23) respectively.

1D steady advection-diffusion equation with low value of diffusion

constant [TC-14]

The DPIELM equations are as follows:

$\overrightarrow{\xi}_{f}^{(i)}$=$\overrightarrow{0}$

[TABLE] 2. 2.

$\overrightarrow{\xi}_{bc}^{(i)}=\overrightarrow{0}$

[TABLE] 3. 3.

$\overrightarrow{\xi}_{C^{0}}^{(i)}=\overrightarrow{0}$

[TABLE] 4. 4.

$\overrightarrow{\xi}_{C^{1}}^{(i)}=\overrightarrow{0}$

[TABLE]

The results of the exact and DPIELM solution is given in Fig(24).

7.1 1D unsteady advection of a sharp peaked Gaussian and a high frequency

wavelet and [ TC-13, TC-15 ]

The DPIELM equation to be solved are as follows:

$\overrightarrow{\xi}_{f}^{(i)}$=$\overrightarrow{0}$

[TABLE] 2. 2.

$\overrightarrow{\xi}_{bc}^{(i)}=\overrightarrow{0}$

[TABLE]

[TABLE] 3. 3.

$\overrightarrow{\xi}_{ic}^{(i)}=\overrightarrow{0}$

[TABLE] 4. 4.

$\overrightarrow{\xi}_{C^{0}}^{(i)}=\overrightarrow{0}$

• •

Continuity along $x$ direction

[TABLE]

where $\rho=1,2...,NB_{t}$, $j(\rho,i)=\kappa(\rho)+i,$ $i=1,2,...,NB_{x}-1$ and $\kappa=[1,(1+NB_{x}),...,1+(NB_{t}-1)NB_{x}]^{T}$.

• •

Continuity along $t$ direction

[TABLE]

where $\rho=1,2...,NB_{x}$, $j(\rho,i)=\kappa(\rho)+iNB_{x},$ $i=1,2,...,NB_{t}-1$ and $\kappa=[1,2,...,NB_{x}]^{T}$.

The results for the TC-13 and TC-15 are given in Fig (25) and Fig (26) respectively.

7.2 1D and 2D unsteady advection-diffusion

equations [ TC-9, TC-10]

In this section, we present the DPIELM equation for the 1D case. The equations for the 2D case can be formulated in a similar fashion. The equations to be solved for 1D unsteady advection-diffusion equation are as follows:

$\overrightarrow{\xi}_{f}^{(i)}$=$\overrightarrow{0}$

[TABLE] 2. 2.

$\overrightarrow{\xi}_{C^{1}}^{(i)}=\overrightarrow{0}$

[TABLE]

where $\rho=1,2...,NB_{t}$, $j(\rho,i)=\kappa(\rho)+i,$ $i=1,2,...,NB_{x}-1$ and $\kappa=[1,(1+NB_{x}),...,1+(NB_{t}-1)NB_{x}]^{T}$.

Rest of the equations i.e. $\overrightarrow{\xi}_{bc}^{(i)}=\overrightarrow{0}$, $\overrightarrow{\xi}_{ic}^{(i)}=\overrightarrow{0}$ and $\overrightarrow{\xi}_{C^{0}}^{(i)}=\overrightarrow{0}$ are same as that used in 1D unsteady linear advection. The results for the 1D and 2D cases are given in Fig (27) and Fig (28 to 30) respectively.

This brings us to the end of our numerical experiments. We close this section by highlighting the key points which are as follows:

The process of partitioning of the whole computational domain into multiple cells simplifies the representation of the complicated function ( and thus PDE ) in the individual cells. As a result, local PIELMs are able to capture not just the functions with sharp gradients, but also discontinuous functions ( TC-11 ). 2. 2.

To our knowledge, this is the first demonstration of capability of ELM based algorithms to solve 2D unsteady PDEs and produce results comparable to sophisticated numerical methods. (TC-10) 3. 3.

The distributed version of PIELM exhibits better representation ability than a deep PINN (TC-15).

8 Conclusion and future work

We have presented in this paper PIELM – an efficient method to utilize physics informed neural networks to solve stationary and time dependent linear PDEs. As PIELM inherits the unique qualities of its parent algorithms (PINN and ELM), it works very well on complex geometries, respects the inherent physics of the PDEs and is extremely fast as well. This leads to several advantages over existing conventional numerical methods; PIELM reduces numerical artefacts such as false diffusion as well can handle complex geometries in a meshfree approach. PIELM also reduces, to a certain extent, the arbitrariness of the number of neurons in typical deep PINNs. Our numerical tests also confirm that, for a fixed problem, our minimal PIELM is more accurate than prior deep NN results [3] while being faster.

We have also presented in this paper the limitations in representing complex functions using a single PIELM or PINN for the whole domain. For practical problems using PINNs can lead to very deep networks with the concomitant training problems and efficiency issues. Our proposed solution is to use a distributed PIELM (DPIELM) which uses different representations in different portions of the domain while imposing some continuity and differentiability constraints. The resultant DPIELM easily captures profiles that PINNs have difficulty with. Further, on time-dependent problems, DPIELM gives results that are comparable to sophisticated conventional numerical techniques as seen in Section 7.2.

We believe that the method, as formulated, is already very powerful for linear PDEs with constant or space-varying coefficients. Two primary areas of development remain, in our opinion. Our preliminary tests show that, for linear PDEs, unlike deep PINNs [13], our method may actually be competitive with conventional techniques in terms of speed and accuracy. However, a firm conclusion on this cannot be reached until a more thorough and fair study is done. Theoretical and numerical evidence for this efficacy is the first area which deserves attention, as the number of practical applications (such as heat conduction, etc) where numerical methods for linear equations are a staple is large. Having an efficient neural network framework for such problems can be tremendously beneficial to practitioners.

Even more importantly, PIELM’s efficacy is thanks to the linear nature of the final problem. We have, therefore, limited our present study to linear problems. Extending this method to nonlinear equations is the obvious next frontier. This may be approached in one of two ways. The first approach would be to use the PIELM structure as is and then solve the resultant, non-convex optimization problem. Another approach would be to linearize the equation around the current time step to predict a future time step. This would be the equivalent of a linearized, explicit time-stepping method in conventional time marching techniques. We are currently investigating these approaches and will report progress in future publications.