Correcting model misspecification in physics-informed neural networks
  (PINNs)
Data-driven discovery of governing equations in computational science has
emerged as a new paradigm for obtaining accurate physical models and as a
possible alternative to theoretical derivations. The recently developed
physics-informed neural networks (PINNs) have also been employed to learn
governing equations given data across diverse scientific disciplines. Despite
the effectiveness of PINNs for discovering governing equations, the physical
models encoded in PINNs may be misspecified in complex systems as some of the
physical processes may not be fully understood, leading to the poor accuracy of
PINN predictions. In this work, we present a general approach to correct the
misspecified physical models in PINNs for discovering governing equations,
given some sparse and/or noisy data. Specifically, we first encode the assumed
physical models, which may be misspecified, then employ other deep neural
networks (DNNs) to model the discrepancy between the imperfect models and the
observational data. Due to the expressivity of DNNs, the proposed method is
capable of reducing the computational errors caused by the model
misspecification and thus enables the applications of PINNs in complex systems
where the physical processes are not exactly known. Furthermore, we utilize the
Bayesian PINNs (B-PINNs) and/or ensemble PINNs to quantify uncertainties
arising from noisy and/or gappy data in the discovered governing equations. A
series of numerical examples including non-Newtonian channel and cavity flows
demonstrate that the added DNNs are capable of correcting the model
misspecification in PINNs and thus reduce the discrepancy between the physical
models and the observational data. We envision that the proposed approach will
extend the applications of PINNs for discovering governing equations in
problems where the physico-chemical or biological processes are not well
understood.
Zongren Zou, Xuhui Meng, George Em Karniadakis
# Correcting model misspecification in physics-informed neural networks (PINNs)

###### Abstract

Data-driven discovery of governing equations in computational science has emerged as a new paradigm for obtaining accurate physical models and as a possible alternative to theoretical derivations. The recently developed physics-informed neural networks (PINNs) have also been employed to learn governing equations given data across diverse scientific disciplines, e.g., in biology and fluid dynamics. Despite the effectiveness of PINNs for discovering governing equations, the physical models encoded in PINNs may be misspecified in complex systems as some of the physical processes may not be fully understood, leading to the poor accuracy of PINN predictions. In this work, we present a general approach to correct the misspecified physical models in PINNs for discovering governing equations, given some sparse and/or noisy data. Specifically, we first encode the assumed physical models, which may be misspecified in PINNs, and then employ other deep neural networks (DNNs) to model the discrepancy between the imperfect models and the observational data. Due to the expressivity of DNNs, the proposed method is capable of reducing the computational errors caused by the model misspecification and thus enables the applications of PINNs in complex systems where the physical processes are not exactly known. Furthermore, we utilize the Bayesian physics-informed neural networks (B-PINNs) and/or ensemble PINNs to quantify uncertainties arising from noisy and/or gappy data in the discovered governing equations. A series of numerical examples including reaction-diffusion systems and non-Newtonian channel and cavity flows demonstrate that the added DNNs are capable of correcting the model misspecification in PINNs and thus reduce the discrepancy between the physical models encoded in PINNs and the observational data. In addition, the B-PINNs and ensemble PINNs can provide reasonable uncertainty bounds in the discovered physical models, which makes the predictions more reliable. We also demonstrate that we can seamlessly combine the present approach with the symbolic regression to obtain the explicit governing equations upon the training of PINNs. We envision that the proposed approach will extend the applications of PINNs for discovering governing equations in problems where the physico-chemical or biological processes are not well understood.

keywords: model uncertainty, physics-informed neural networks, model misspecification, uncertainty quantification, non-Newtonian flows, symbolic regression +
Footnote †: journal:
## 1 Introduction

Accurate governing equations or physical models are essential for understanding and quantifying physical processes across diverse scientific disciplines, such as the Navier-Stokes equations in fluid dynamics and the transport equations in geosciences. Generally, the governing equations are derived based on established theories such as the conservation of mass and momentum, thermodynamic laws, and so on. However, there are also complex systems in which the governing equations are difficult to obtain theoretically due to the lack of understanding for certain physical processes, e.g., the constitutive relation in highly nonequilibrium flows [1], the forcing caused by solar and volcanic variability in climate problems [2], non-equilibrium reactions, etc. How to obtain analytical expressions that represent physical phenomena in nature and engineering remains an open question.

With the rapid growth of available data and computing power, data-driven discovery of governing equations based on machine learning algorithms has emerged as a new paradigm for obtaining accurate physical models as the alternative to theoretical derivations [3; 4; 5; 6; 7; 8; 9; 10; 11]. To name a few examples, Schmidt et al. presented a symbolic regression approach to discover equations from experimental data [3]; Brunton et al. proposed the sparse identification of nonlinear dynamical systems (SINDy) to learn both ordinary and partial differential equations from data [4; 5]. In addition, Chen et al. employed deep neural networks (DNNs), more specifically the physics-informed neural networks (PINNs), which encode the physical laws into DNNs via automatic differentiation [12], to discover governing equations from scarce data [6]. In most of the aforementioned approaches, a library of candidate terms that describe all the possible physical processes is required as prior knowledge [3; 4; 5; 6]. However, mathematical expressions or models for certain physical processes are quite challenging to obtain in some complex real-world applications, e.g., chemical reactions in combustion, the constitutive relation in highly partially ionized plasmas, etc. It is probable that we miss some details of the physical process, hence misspecifying the physical models in such scenarios [13]. An incomplete or non-comprehensive library leads to the discrepancy between the observational data and the assumed governing equations. In other words, the discovered governing equations may not be accurate enough to describe the physical processes represented by the observational data.

Several recent approaches have been developed to reduce the discrepancy between the observational data and the assumed physical models. For instance, Ebers et al. [14] proposed to learn the _missing physics_ or the discrepancy between the imperfect models and the measurements using a plurality of methods, e.g. Gaussian process regression (GPR), DNNs, etc.; Chen et al. also employed DNNs to recover the _missing physics_ in prior models [15]; in [16; 17], combinations of PINNs and symbolic regression were employed to identify the missing physics in biological problems, e.g. Alzheimer's disease; Saurabh et al. [18] and Zhu et al. [19] discovered source terms of dynamic systems in a "gray-box" fashion with NN-based numerical integrators. However, most of the existing work focuses on resolving the missing physics. Investigations on scenarios that we _misspecify the physical models_ are rare to the best of our knowledge. Further, the observational data from sensors in real-world applications are generally noisy, and may also be incomplete due to the limitation of data collection
techniques. Such noisy and incomplete or gappy data lead to uncertainties in the predictions of machine learning algorithms or more specifically the scientific machine learning algorithms [20; 21]. In addition, as widely discussed in the climate change modeling [13; 22; 23], the incomplete understanding of physical processes will also result in uncertainties in simulations. We point out that quantifying uncertainties arsing from noisy/gappy data as well as the incomplete understanding of the physical processes is an important issue but has been largely ignored in the existing literature on data-driven discovery of governing equations. In the current study, we refer to the uncertainties in the discovered physical models as _model uncertainty_.

In the present work, we aim to address the issue of misspecifying physical models as well as quantifying model uncertainty in the data-driven discovery of governing equations. Specifically, PINNs are employed as the backbone to discover the governing equations or physical models from data due to their effectiveness as well as easy implementations. In particular, we first encode the assumed physical models that may be misspecified in PINNs, and add another DNN as the correction to the misspecified model to alleviate the discrepancy between the imperfect models and the observational data. We note that the library of candidate terms used in [3; 4; 5; 6] is not required here since the added DNN can serve as a universal approximator to model the discrepancy. The Bayesian physics-informed neural networks (B-PINNs) [24] and/or ensemble PINNs [20; 21] are used to quantify the uncertainties in the discovered governing equations or physical models arising from the noisy and/or gappy data. Further, the symbolic regression is utilized to obtain the analytical expression for the unknown physics based on the trained DNNs that are used to correct the model misspecifications in PINNs.

The rest of this paper is organized as follows. In Sec. 2, we introduce the model misspecification in PINNs as well as our approach to correct it, and also the methods for quantifying uncertainties in PINN predictions. In Sec. 3, we conduct four numerical experiments to demonstrate the effectiveness of the proposed approach in addressing the model misspecification issue in PINNs. We summarize the present work in Sec. 4.

## 2 Model misspecification in PINNs for learning governing equations

### Learning governing equations from data using PINNs

Consider the following nonlinear ODE/PDE describing a physical system:

\[\mathcal{F}_{\lambda}[u](x) =f(x),x\in\Omega, \tag{1a}\] \[\mathcal{B}_{\lambda}[u](x) =b(x),x\in\partial\Omega, \tag{1b}\]

where \(\Omega\) is the domain, \(\lambda\) is the model parameter, \(u\) is the sought solution, \(f\) is the source term, \(b\) is the boundary term, \(\mathcal{F}_{\lambda}\) and \(\mathcal{B}_{\lambda}\) are operators parameterized by \(\lambda\), defining the equation and the boundary condition. The PINN method is capable of solving ODEs/PDEs given initial/boundary conditions as well as learning governing equations and identifying model parameters from data [24; 21; 25; 26; 27; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42]. A schematic view of the PINN method is presented in Fig. 1, in which a DNN with \(\sigma\) as the activation function and \(\theta\) as the parameter is built to approximate \(u\) with \(u_{\theta}\) and \(f,b\) with \(f_{\theta},b_{\theta}\), respectively, via automatic differentiation [12].
In the present study, we focus on using PINNs to learn governing equations from data. The dataset for learning governing equations with PINNs is denoted by \(\mathcal{D}\), which can be expressed as \(\mathcal{D}=\mathcal{D}_{u}\cup\mathcal{D}_{f}\cup\mathcal{D}_{b}\), where \(\mathcal{D}_{u}=\{x_{i}^{u},u_{i}\}_{i=1}^{N_{u}}\) is the set of data for \(u\), and \(\mathcal{D}_{f}=\{x_{i}^{f},f_{i}\}_{i=1}^{N_{f}}\) and \(\mathcal{D}_{b}=\{x_{i}^{b},b_{i}\}_{i=1}^{N_{b}}\) are sets of data for the physics. The model parameter \(\lambda\), which defines the equation, is then obtained by optimizing the following PINN loss function with the gradient-descent method and its variants [43]:

\[\mathcal{L}(\theta)=w_{u}\mathcal{L}_{\mathcal{D}_{u}}(\theta)+\mathcal{L}_{PDE }(\theta), \tag{2}\]

where

\[\mathcal{L}_{\mathcal{D}_{u}}(\theta) =\frac{1}{N_{u}}\sum_{i=1}^{N_{u}}||u_{\theta}(x_{i}^{u})-u_{i}|| _{2}^{2},\] \[\mathcal{L}_{PDE}(\theta) =\frac{w_{f}}{N_{f}}\sum_{i=1}^{N_{u}}||\mathcal{F}[u_{\theta}]( x_{i}^{f})-f_{i}||_{2}^{2}+\frac{w_{b}}{N_{b}}\sum_{i=1}^{N_{b}}||\mathcal{B}[u_{ \theta}](x_{i}^{b})-b_{i}||_{2}^{2},\]

\(w_{u},w_{f},w_{b}\) are belief weights for balancing different terms, and \(||\cdot||_{2}\) is the \(\ell^{2}\)-norm for finite-dimensional vector. We note that in certain cases we may not have data corresponding to the boundary/initial conditions, i.e., the last term in \(\mathcal{L}_{PDE}(\theta)\), and we can directly drop this term in the computations.

### Correcting the model misspecification in PINNs

The computational accuracy of PINNs in learning governing equations strongly depends on the encoded physical models, i.e. \(\mathcal{F}_{\lambda}\) and/or \(\mathcal{B}_{\lambda}\) in Eq. (1). However, the correct specification of the physics is not always guaranteed, leading to significant discrepancy between

Figure 1: A schematic view of PINNs for learning governing equations from data. On the left, a general framework of PINNs is displayed where \(\sigma\) denotes the activation function, \(\theta\) denotes the DNN parameters, \(u_{\theta}\) denotes the output of the DNN that is used to approximate \(u\), and \(f_{\theta}/b_{\theta}\) are computed via automatic differentiation [12]. An example on learning governing equations with PINNs is illustrated on the right panel. Here a 1D Burgers equation is learned (\(\lambda_{1}=0,\lambda_{2}=1,\lambda_{3}=-\nu\) where \(\nu\) denotes the viscosity) given data on \(u\).
the observational data and the assumed governing equations. An example is demonstrated in Fig. 2 where the flow is non-Newtonian but the physical model is misspecified as the one for Newtonian flows (see Sec. 3.3.1 for details). As a result, the physics and the data cannot be fitted well at the same time: the training does not land on the region "IV", where both the losses for the data as well as the residual of the equations are close to zero, regardless of the choice of belief weights in the loss function Eq. (2), i.e. \(w_{f}\) for the physics and \(w_{u}\) for the data. Three results from PINNs are sampled and displayed in Fig. 2: when \(u_{\theta}\) fits the data (which come from a non-Newtonian flow), it does not satisfy the governing equations for Newtonian flows; when \(u_{\theta}\) is forced to follow the Newtonian flow by setting large \(w_{PDE}\)

Figure 2: In (a), a breakdown of total uncertainty in PINNs is presented. In (b), results from PINNs encountering model misspecifications are displayed, where a physical model of Newtonian fluid is assumed but the real physics, which generates the data, is non-Newtonian. The PINN method is employed to identify the homogeneous viscosity with different belief weights in the loss function, i.e. pairs of \(w_{PDE}\) and \(w_{\mathcal{D}_{u}}\) in (2), resulting in different pairs of PDE loss \(\mathcal{L}_{PDE}\) and data loss \(\mathcal{L}_{\mathcal{D}_{u}}\), shown in the left of (b). Due to the model misspecification, the PDE loss \(\mathcal{L}_{PDE}\) and the data loss \(\mathcal{L}_{\mathcal{D}_{u}}\) cannot be minimized simultaneously. Three results are presented in the right of (b), corresponding to three cases in the left, + (data are fitted well but the physics is not satisfied), * (the physics is satisfied but data are not fitted), and x (neither data nor physics are fitted). We note that the boundary condition is hard-encoded in the modeling. Details can be found in Sec. 3.3.1.
it presents Newtonian behavior and thus cannot fit the data from the non-Newtonian flow well.

The discrepancy between observational data and the assumed governing equations is caused by the physics/model misspecification. In this regard, we propose a general approach developed upon the original PINN framework to quantify the discrepancy as well as to correct the model misspecification. In this work, we focus on the misspecification of the differential operator \(\mathcal{F}_{\lambda}\), but the proposed approach can be easily generalized to the boundary/initial operator. Here, we denote the misspecified differential operator by \(\tilde{\mathcal{F}}_{\lambda}\) and then the discrepancy can be quantified as follows:

\[s(x)=f(x)-\tilde{\mathcal{F}}_{\lambda}[u](x). \tag{3}\]

Inspired by the PINN method, in which the sought solution \(u\) is modeled with a DNN parameterized by \(\theta\), we employ another DNN parameterized by \(\psi\) to model \(s(x)\), denoted as \(s_{\psi}(x)\). The source term \(f(x)\) is approximated by \(\tilde{\mathcal{F}}_{\lambda}[u_{\theta}](x)+s_{\psi}(x)\), and hence the term penalizing the violation of the physics in the loss function (2) is reformulated as:

\[\mathcal{L}_{PDE}(\theta,\psi)=\frac{w_{f}}{N_{f}}\sum_{i=1}^{N_{u}}||\tilde{ \mathcal{F}}_{\lambda}[u_{\theta}](x_{i}^{f})+s_{\psi}(x_{i}^{f})-f_{i}||_{2 }^{2}+\frac{w_{b}}{N_{b}}\sum_{i=1}^{N_{b}}||\mathcal{B}_{\lambda}[u_{\theta} ](x_{i}^{b})-b_{i}||_{2}^{2}. \tag{4}\]

We note that \(s_{\psi}\) plays two important roles in our approach handling model misspecification: (1) it approximates the function quantifying the model discrepancy, i.e. \(s\), and (2) it allows both the data term \(\mathcal{L}_{\mathcal{D}_{u}}(\theta)\) and the PDE term \(\mathcal{L}_{PDE}(\theta,\psi)\) to be minimized to very small

Figure 3: An overview of the proposed approach for correcting the misspecified physics. Here we assume that the nonlinear differential operator \(\mathcal{F}\) is misspecified as \(\tilde{\mathcal{F}}_{\lambda}\) such that \(\tilde{\mathcal{F}}_{\lambda}[u]\neq f\), which leads to poor performance of PINNs in learning governing equations. Compared to the conventional PINN method shown in Fig. 1, an additional DNN, parameterized by \(\psi\), is used to model the discrepancy. We note that \(s_{\psi}\) serves two roles: (1) correcting the equation such that data of \(u\) can be explained, and (2) modeling the discrepancy of the physics. With posterior samples of \(\psi\), we obtain the prediction as well as the model uncertainty of the discrepancy. We also note that DNNs in this figure, namely \(u_{\theta}\) and \(s_{\psi}\), are deterministic NNs in ensemble PINNs and Bayesian NNs (BNNs) in B-PINNs.
values simultaneously in the training of PINNs such that data are fitted by \(u_{\theta}\) and the (corrected) physics are satisfied.

### Uncertainty quantification

As mentioned in Sec. 2.2, a key step in the present work is to learn the correction term, i.e., \(s_{\psi}\), that is used to fix the issue of model misspecifications in PINNs. Generally, the measurements in real-world applications are noisy and can also be incomplete, leading to the uncertainties in predictions from PINNs [20; 21]. We refer to the uncertainties in \(s_{\psi}\) as _model uncertainty_. Here, two typical approaches are employed for quantifying the uncertainties in \(s_{\psi}\), i.e., ensemble PINNs [20; 21] and B-PINNs [24]. The former is computationally efficient but has difficulties preventing overfitting when the data are noisy, while the latter is able to handle the noisy data well but is computationally more expensive. We therefore utilize the ensemble PINNs and B-PINNs for cases with clean and noisy data, respectively.

We now briefly review the ensemble PINNs and B-PINNs. The ensemble PINN method follows the deep ensemble method, proposed in [44], and employs multiple standard training of PINNs independently. That is, PINNs with the PDE loss term in Eq. (4) are trained multiple times with random initialization of the DNN parameter. Each random initialization corresponds to one initial guess for minimizing the PINN loss function (2) and each training is finding one global minimum. As a result, ensemble PINNs can be considered as identifying different maximum a posteriori (MAP) estimates of the PINN loss function [20]. On the other hand, the B-PINN method focuses on estimating the posteriors for \(\theta\) and \(\psi\) based on the Bayes' rule as follows:

\[p(\theta,\psi|\mathcal{D})\propto p(\mathcal{D}|\theta)p(\theta). \tag{5}\]

Here \(p(\theta|\mathcal{D})\) is the density function for the posterior distribution, \(p(\mathcal{D}|\theta)\) for the likelihood, and \(p(\theta)\) for the prior. Given independent and identically distributed (i.i.d.) data, the density function for the likelihood distribution can be written as follows:

\[p(\mathcal{D}|\theta)=p(\mathcal{D}_{u}|\theta)p(\mathcal{D}_{f}|\theta)p( \mathcal{D}_{b}|\theta)=\prod_{i=1}^{N_{u}}p(u_{i}|x_{i}^{u},\theta)\prod_{i= 1}^{N_{f}}p(f_{i}|x_{i}^{f},\theta)\prod_{i=1}^{N_{b}}p(b_{i}|x_{i}^{b},\theta). \tag{6}\]

If we further assume that the data for \(u\), \(f\) and \(b\) follow Gaussian distribution with mean zero and standard deviation, \(\sigma_{u}\), \(\sigma_{f}\) and \(\sigma_{b}\), respectively, then correcting it follows that:

\[p(u_{i}|x_{i}^{u},\theta) =\frac{1}{\sqrt{2\pi}\sigma_{u}}\exp(-\frac{||u_{\theta}(x_{i})- u_{i}||_{2}^{2}}{2\sigma_{u}^{2}}), \tag{7a}\] \[p(f_{i}|x_{i}^{f},\theta) =\frac{1}{\sqrt{2\pi}\sigma_{f}}\exp(-\frac{||\mathcal{F}_{ \lambda}[u_{\theta}](x_{i}^{f})+s_{\psi}(x_{i}^{f})-f_{i}||_{2}^{2}}{2\sigma_{ f}^{2}}),\] (7b) \[p(b_{i}|x_{i}^{b},\theta) =\frac{1}{\sqrt{2\pi}\sigma_{b}}\exp(-\frac{||\mathcal{B}_{ \lambda}[u_{\theta}](x_{i}^{b})-b_{i}||_{2}^{2}}{2\sigma_{b}^{2}}),. \tag{7c}\]

We note that \(\sigma_{u}\), \(\sigma_{f}\) and \(\sigma_{b}\) are often treated as scales of the additive Gaussian noise in measuring data of \(u\), \(f\)and \(b\), respectively. Various methods, e.g. Markov Chain Monte Carlo (MCMC) and variational inference, have been proposed to tackle the posterior distribution (5). In this work, we choose to use Hamiltonian Monte Carlo (HMC) method to estimate the
posterior distribution by sampling \(\theta\) and \(\psi\) from it. We also note that here we only consider the homogeneous noise, while more complicated noise models, e.g. the heteroscedastic noise [20], are also compatible with the framework.

In both the ensemble PINN and B-PINN methods, we obtain samples of the NN parameters, from which statistics of the prediction are estimated. In general, the mean and the standard deviation are often chosen to represent the prediction and the uncertainty [20], respectively. Specifically in this work, we obtain samples of \(\theta\) and \(\psi\) from ensemble PINNs and/or B-PINNs, denoted by \(\{\theta_{i}\}_{i=1}^{M}\) and \(\{\psi_{i}\}_{i=1}^{M}\), respectively, where \(M\) denotes the number of samples. Then, the predicted means of the sought solution \(u\) and the discrepancy \(s\) are estimated as follows:

\[\hat{u}(x)\approx\frac{1}{M}\sum_{i=1}^{M}u_{\theta_{i}}(x),\ \hat{s}(x)\approx \frac{1}{M}\sum_{i=1}^{M}s_{\psi_{i}}(x), \tag{8}\]

where \(\hat{u}\) and \(\hat{s}\) denote the mean of \(u\) and \(s\), respectively. Similarly, the predicted uncertainties are:

\[Var[u](x)\approx\frac{1}{M}\sum_{i=1}^{M}(u_{\theta_{i}}(x)-\hat{u}(x))^{2},\ Var[s](x)\approx\frac{1}{M}\sum_{i=1}^{M}(s_{ \psi_{i}}(x)-\hat{s}(x))^{2}, \tag{9}\]

where \(Var[u]\) and \(Var[s]\) denote the variance of \(u\) and \(s\), respectively. Interested readers are directed to [20] for a review and [21] for a comprehensive Python library termed NeuralUQ.

We note that upon the computation of \(s_{\psi}(x)\), it can be used in various way: (1) its value represents the model discrepancy and its uncertainty represents the uncertainty in the physical model which we refer to as _model uncertainty_; (2) it corrects the misspecified physical models and identifies a system which explains the data; and (3) it could also be used to obtain mathematical expressions for the missing details of certain physical process via the symbolic regression [16; 17; 45]. Specifically in the present work, we employ the symbolic regression method developed in [46] and the associated Python library PySR to regress analytic mathematical expressions for the unknown physical processes (See Sec. 3.1 for a demonstration example).

## 3 Results and discussion

In this section, we conduct four numerical experiments, i.e., (1) an ODE system, (2) a one-dimensional (1D) reaction-diffusion equation, (3) a two-dimensional (2D) channel flow, and (4) a 2D cavity flow. In the first two examples, reaction models are possibly misspecified, and in the last two examples, we assume that we misspecify the constitutive relation in non-Newtonian flows as Newtonian ones. In all the test cases, we employ DNNs to correct the misspecified models as discussed in Sec. 2. Details for the computations (e.g., training steps, NN architectures, etc.) can be found in A.

### Pedagogical example: An ODE system

We first consider the following one-dimensional ODE:

\[\begin{split}\frac{du}{dt}&=f(t)+\lambda u(1-u),\ t\in[0,1],\\ u(t&=0)=u_{0},\end{split} \tag{10}\]
where \(\lambda\) is a non-negative constant, \(u_{0}\) is the initial condition, and \(f\) is the source term. In certain real-world applications, the reaction term, i.e., \(\lambda u(1-u)\) in the right hand side of Eq. (10), may not be exactly known and herein it is misspecified as an empirical model. In this regard, we study the effect of different prior knowledge about the reaction model in the discovery of governing equation. In particular, the following cases are considered:

1. Case (A): The reaction term is correctly specified as \(\lambda u(1-u)\), where \(\lambda\geq 0\) is unknown.
2. Case (B): The reaction term is misspecified as \(\lambda\cos(u)\), where \(\lambda\geq 0\) is unknown.
3. Case (C): The reaction term is misspecified as \(\lambda\cos(u)\) (with known \(\lambda=0.2\)) but we add \(s(t)\) to correct it, i.e., \(\lambda\cos(u)+s(t)\).

We denote the reaction term as \(\phi\), and the target here is to identify \(\phi\) given data on \(u/f\) as well as differently specified models. In Case (A) and (B), the problem degenerates to identifying the model parameter \(\lambda\), and \(\phi\) in these two cases can then be computed by \(\tilde{\lambda}u_{\theta}(1-u_{\theta})\) and \(\tilde{\lambda}\cos(u_{\theta})\), respectively, where \(\tilde{\lambda}\) is the PINN estimate of \(\lambda\). In Case(C), \(\phi\) is computed by \(\tilde{\lambda}\cos(u_{\theta})+s_{\psi}(t)\) where \(s_{\psi}\) is the estimation of the discrepancy induced by the model misspecification. We note that in Case (C) we assume \(\lambda\) is known and set \(\lambda=0.2\) for simplicity. This is because the discrepancy is already modeled by \(s(t)\) and hence treating \(\lambda\) as unknown does not make any difference in correcting the misspecified model.

To reduce the effect of the noisy/gappy data on the predicted accuracy, we first test the case with clean and sufficient data. The training data are generated by solving Eq. (10) with \(f(t)=\sin(2\pi t)\) and \(u_{0}=0\) using _MATLAB ode45_[47]. We assume that we have 101 data for \(u\) and \(f\), uniformly sampled from \(t\in[0,1]\), and we employ ensemble PINNs [20, 21] to learn the governing equation and also quantify the uncertainties in predictions. In all the test cases here, we choose a DNN with two hidden layers, each of which has 50 neurons and is equipped with hyperbolic tangent as the activation function, to approximate \(u\). The equations are then encoded in PINNs via the automatic differentiation [12]. Further, we utilize an additional DNN to correct the model misspecification in Case (C), which has the same architecture as the one used for approximating \(u\). The Adam optimizer [48] is applied to train the ensemble PINNs.

As shown in Fig. 4(a), PINNs with correct physical model (Case (A)) enable accurate inference of \(\phi\) and the surrogate model \(u_{\theta}\) fits the data and satisfies the ODE simultaneously. In contrast, PINNs with misspecified physical model (Case (B)) yield completely wrong estimate of \(\phi\) as displayed in Fig. 4(b). In Case (C), we add a DNN to correct the misspecified reaction model used in Case (B). As illustrated in Fig. 4(c) as well as Table 1, the predicted

\begin{table}
\begin{tabular}{c|c|c|c|c|c} \hline \hline  & \(\tilde{\lambda}\) & Error of \(\phi\) & Error of \(u\) & Error of \(f\) & Error of \(\tilde{u}\) \\ \hline Case (A): Known model & \(1.5000\pm 0.0001\) & 0.00\% & 0.00\% & 0.06\% & 0.01\% \\ \hline Case (B): Misspecified model & \(0.2254\pm 0.0001\) & 41.46\% & 1.73\% & 6.87\% & 4.58\% \\ \hline Case (C): Corrected model & 0.2 & 0.55\% & 0.01\% & 0.10\% & 0.02\% \\ \hline \hline \end{tabular}
\end{table}
Table 1: ODE system: Predictions for \(\phi\), \(u\) and \(f\) and the reconstruction \(\tilde{u}\). The metric is relative \(L_{2}\) error, and errors are computed between the reference and the mean of 20 PINNs trained independently. The exact value of \(\lambda\) is 1.5 and \(\tilde{\lambda}\) denotes the PINN estimate of \(\lambda\) (\(\tilde{\lambda}=0.2\) is fixed in Case (C)). We note that \(\tilde{u}\) is reconstructed by solving the identified ODE system with the identified \(\phi\), the inferred \(f\) and the initial condition \(u_{0}\). Cubic spline data interpolation is employed to approximate the function \(f\) from the prediction of \(f\) on a uniform grid.
accuracy for \(\phi\) improves significantly as compared to the result in Case (B). In addition, the computational accuracy for \(u\), \(f\) and \(\phi\) in Case (C) is quite similar as in the Case (A), where the physical model is correctly specified. The above results demonstrate the effectiveness of the present approach. We note that the predicted uncertainties in all cases are quite small, which is reasonable since we use sufficient data for both \(u\) and \(f\) to determine the unknown reaction term here. Furthermore, we present further results as well as the related discussion for the case with clean but gappy data in B.1.

We now keep the same setup as in Case (A) - (C) but assume that the training data are both noisy and gappy. In these cases, 50 data of \(f\) and \(u\) are randomly sampled from \(t\in[0,1]\) and corrupted by independent additive Gaussian noises with known noise scales (0.05 for \(f\) and 0.01 for \(u\)). As discussed in Sec. 2, we will employ the B-PINNs to quantify

Figure 4: ODE system: Predictions for \(u\), \(f\) and \(\phi\) from ensemble PINNs with clean and sufficient data. Ensemble PINNs with (a) the correct physical model, (b) a misspecified physical model, and (c) the misspecified physical model plus a correction. Blue circles: Representative measurements for \(u\) and \(f\); Red dashed line: Predicted mean; Black solid line: Reference solution.
the uncertainties arising from the noise as well as the incompleteness of the training data. As shown in Fig. 5, (1) B-PINNs with correct physical model can obtain accurate predictions for both \(u\), \(f\), and \(\phi\), and the computational errors for \(\phi\) are generally bounded by the predicted uncertainties (Case (A), Fig. 5(a)); (2) B-PINNs with the misspecified physical model fails to provide accurate inference for \(u\), \(f\), or \(\phi\) (Case (B), Fig. 5(b)); and (3) B-PINNs with a DNN as correction to the misspecified model is able to improve the predicted accuracy for \(u\), \(f\), and \(\phi\), as compared to the results in Case (B) (Case (C), Fig. 5(c)). Similar as in Case (A), the computational errors for \(\phi\) are bounded by the predicted uncertainties here.

Figure 5: ODE system: Predictions for \(u\), \(f\) and \(\phi\) from B-PINNs with noisy and gappy data. B-PINNs with (a) correct physical model, (b) misspecified physical model, and (c) misspecified physical model plus a correction. Blue circles: Noisy measurements for \(u\) and \(f\); Red dashed line: Predicted mean from B-PINNs; Black solid line: Reference solution. Our approach alleviates the model misspecification issue significantly, providing more accurate predictions. Compared to correctly specified model, the increase in the predicted uncertainty is caused by the lack of knowledge in the physical model and could be considered as model uncertainty.
In addition to the relation between the computational errors and the predicted uncertainties for \(\phi\) discussed above, we would like to further discuss that: (1) it is interesting to observe in Case (A) and Case (C) that the predicted uncertainty near the right end is larger than that near the left end. We note that the results here are reasonable since we do not have data for \(f\) on the right end; (2) the predicted uncertainties of \(\phi\) for \(t\in[0,1]\) in Case (C) are larger than those in Case (A) although we use the same training data as well as prior for the BNN that is used to approximate \(u\). We remark that the difference between the uncertainties in these two cases is caused by the lack of knowledge on the physical model and reflects the uncertainty in physical models.

We further verify the learned governing equation from the proposed approach. Specifically, we solve the equation \(du/dt=f(t)+\phi(t)\) with the inferred \(f\), \(\phi\), and the initial condition \(u_{0}=0\), using MATLAB ode45 to see if the solution for \(u\) agrees with the reference solution. We denote the solution to the identified system by \(\tilde{u}\). In particular, 1000 posterior samples of \(f\) and \(\phi\) from the B-PINNs are used to compute \(\tilde{u}\), as displayed in Fig. 6. We further illustrate the relative \(L_{2}\) errors between the mean for \(\tilde{u}\) and the reference solution in Table 2. As we can see, PINNs with correct physical model leads to the lowest error in \(\tilde{u}\). Adding a DNN to correct the misspecified model can significantly improve the accuracy in \(\tilde{u}\) as compared to the case with misspecified model, and the computational errors are slightly larger than the ones from the correct physical model. We also present the computational errors for \(\tilde{u}\) from the ensemble PINNs in Table 1, which will not be discussed again since the results are similar as the above ones. Finally, we would like to point out that (1) the

\begin{table}
\begin{tabular}{c|c|c} \hline \hline  & \(\tilde{\lambda}\) & Error of \(\tilde{u}\) \\ \hline Case (A): Known model & \(1.4867\pm 0.0488\) & 1.56\% \\ \hline Case (B): Misspecified model & \(0.0000\pm 0.0000\) & 12.73\% \\ \hline Case (C): Corrected model & 0.2 & 2.61\% \\ \hline \hline \end{tabular}
\end{table}
Table 2: ODE system: Relative \(L_{2}\) error for the reconstructed \(u\) from B-PINNs. Here the computational errors are calculated using the predicted mean from B-PINNs.

Figure 6: ODE system: Solving the identified system with the inferred \(f\) and \(\phi\) from B-PINNs with the initial condition \(u_{0}=0\). The inferred \(f\) and \(\phi\) come from (a) correct physical model, (b) misspecified physical model, and (c) our approach. Red dashed line: Predicted mean; Black solid line: Reference solution. It shows that the misspecified physical model leads to wrong reconstruction \(\tilde{u}\) while our approach is able to correct it. Due to the lack of knowledge of the model, our approach yields larger uncertainty than the correct physical model.
uncertainties in the reconstructed \(u\) (i.e., \(\tilde{u}\)) are caused by the predicted uncertainties for \(f\) and \(\phi\) from B-PINNs, and (2) the computational errors for the mean of \(\tilde{u}\) in Case (A) and (C) are of the same order as the noise in the measurements, suggesting the reasonableness of the current results.

As discussed in Sec. 2, the discovered physical models or governing equations up to now are expressed using the trained DNNs (i.e., \(s_{\psi}\)), which are not explicitly presented. We can then combine the proposed method with the symbolic regression as the further step to obtain the explicit governing equations. Here we take the results in Case (C) from both ensemble PINNs and B-PINNs to demonstrate how to get the analytical expressions for the trained NNs, i.e., \(s_{\psi}\). Recall that with the present methodology the reaction term is now expressed as \(\phi(t)=\tilde{\lambda}\cos(u_{\theta}(t))+s_{\psi}(t)\). Assume that we now have the knowledge that the reaction term is autonomous with respect to the solution, i.e. \(\phi=s(u)\). We can then adopt the symbolic regression developed in [46] to find the analytical expression of \(s(u)\) from \(u_{\theta}\) and \(s_{\psi}\). The data for the symbolic regression are generated as \(\{u_{\theta}(t_{i}),s_{\psi}(t_{i})+\tilde{\lambda}\cos(u_{\theta}(t_{i}))\}_ {i=1}^{N}\), where \(t_{i},i=1,...,N\) are 101 uniform grids points on \([0,1]\), and \(u_{\theta}\) and \(s_{\psi}\) are the predicted means from the trained ensemble PINNs/B-PINNs. As displayed in Table 3, the recovered reaction model via symbolic regression agrees well with the reference solution as we have sufficient clean data. However, for the case with noisy and gappy data, the expression for the reaction model from symbolic regression does not match the reference solution. The result is expected since we can see in Fig. 5(c) that the predicted mean for \(\phi\) does not agree well with the reference solution well and the predicted uncertainties are generally large. Generally, increasing the number of training data is able to reduce the predicted uncertainty in B-PINNs [20; 21]. Here we also show that the learned reaction model agrees better with the reference solution as we increase the training data from 50 to 101 in B-PINNs. A recommendation here is to use the predicted uncertainty as guidance, i.e., we use the symbolic regression for discovering equation when the predicted uncertainty is small to achieve better accuracy. More results for the symbolic regression can be found in B.2.

### Reaction-diffusion equation

We now consider the following time-dependent PDE that describes a reaction-diffusion system:

\[\begin{split}&\frac{\partial u}{\partial t}=D\frac{\partial^{2}u }{\partial x^{2}}+\lambda g(u)+f(x,t),\ x\in[0,1],\ t\in[0,1],\\ & u(x,0)=0.5\sin^{2}(\pi x),\ x\in[0,1],\\ & u(0,t)=u(1,t)=0,\ t\in[0,1],\end{split} \tag{11}\]

where \(D=0.01\) denotes the diffusion coefficient, \(g(u)=u(1-u)\) is the reaction model, \(\lambda\) is the reaction rate, and \(f(x,t)=0\) is the source term. Similar as in Sec. 3.1, we test the following three cases:

\begin{table}
\begin{tabular}{c|c|c|c} \hline \hline Reference & Clean and sufficient data & Noisy and gappy data & Noisy and sufficient data \\ \hline \(1.5u(1-u)\) & \(1.4980u-1.4929u^{2}\) & \(u/(u+0.6213)\) & \(1.3540u-u^{2}\) \\ \hline \hline \end{tabular}
\end{table}
Table 3: ODE system: Discovered equations from \(\phi(u)\) using a symbolic regression method [46] (see Fig B.14 for plots of the discovered models). In all three cases, the symbolic regression is performed on the predicted mean from ensemble PINNs and B-PINNs.
1. Case (A): The reaction model is correctly specified and the reaction rate \(\lambda\geq 0\) is unknown.
2. Case (B): The reaction model is misspecified as \(u^{2}\) and the reaction rate \(\lambda\geq 0\) is unknown.
3. Case (C): The reaction model is misspecified as \(u^{2}\) and we correct it with a DNN, i.e., \(\lambda u^{2}+s\) (with known \(\lambda=1\)).

We denote the reaction term as \(\phi\) and we would like to identify \(\phi\) from data of \(u\), \(f\) and different physical models. In Case (A)/(B)/(C), the estimate of \(\phi\) is \(\tilde{\lambda}u_{\theta}(1-u_{\theta})/\tilde{\lambda}u_{\theta}^{2}/\tilde {\lambda}u_{\theta}^{2}+s_{\psi}\), where \(\tilde{\lambda}\) is the estimate of \(\lambda\). We note that in Case (C) we assume \(\lambda=1\) is known because the model misspecification induced by wrong \(\lambda\) can be fully covered by \(s\).

Here we only consider a general scenario, i.e., the training data are noisy and gappy. In this case, we employ B-PINNs to identify \(\phi\) with uncertainties. We assume that we have 121 measurements on \(u\) and 195 measurements of \(f\), which are evenly distributed in the temporal-spatial domain \((x,t)\in[0,1]\times[0,1]\). In addition, the measurement errors for both \(u\) and \(f\) are assumed to be Gaussian, and the noise scales are 0.02 and 0.05 for \(u\) and \(f\), respectively. Further, the reference solution is obtained by solving Eq. (11) with \(g(u)=u(1-u)\), \(\lambda=2\), and \(f(x,t)=0\). Similar as in Sec. 3.1, we use one DNN to model the sought solution \(u\) and another DNN to approximate the discrepancy \(s\) in Case (C). Details of the B-PINN method, e.g. hyperparameters of NNs and HMC, can be found in A.

As show in Fig. 7, we observe that: (1) the predicted means for \(u\) and \(\phi\) agree the best with the reference solutions as the reaction model is correctly specified (Case (A), second column in Fig. 7); (2) the predictions for \(u\) and \(\phi\) show significant discrepancy as compared to the reference solution when the reaction model is misspecified (Case (B), third column in

Figure 7: Time-dependent reaction-diffusion system: Predicted mean for \(u\) and \(\phi\) from B-PINNs with noisy data of \(u\) and \(f\) and differently specified reaction models. Misspecified model (\(u^{2}\)) impairs the performance of B-PINNs significantly in terms of the accuracy in predicting \(u\) and \(\phi\). and the proposed approach is able to alleviate the issue caused by model misspecification; see Fig. 8 for UQ.
Fig. 7); (3) the predicted mean for \(u\) in Case (C) is similar as that in Case (A), in which the reaction model is correctly specified; and (4) the predicted mean for \(\phi\) in Case (C) is not as good as the one in Case (A), but it is much better than the result in Case (B) where the reaction model is misspecified. All the results demonstrate that the added DNNs are able to correct the misspecified model which in turn enhance the computational accuracy in PINNs.

We further illustrate the predicted uncertainties for both \(u\) and \(\phi\) in Fig. 8. In particular, we only depict the uncertainties at one representative time for simplicity. As we can see, (1) the computational errors for both \(u\) and \(\phi\) are generally bounded by the predicted uncertainties in Case (A) (first column in Fig. 8) and Case (C) (third column in Fig. 8), while they cannot be bounded by the predicted uncertainties in the case with misspecified reaction model i.e. Case (B) (second column in Fig. 8), and (2) the predicted uncertainties for \(u\) and \(\phi\) are generally larger in Case (C) (third column in Fig. 8) than those in Case (A) (first column in Fig. 8), which are similar as the results in Sec. 3.1. Again, we attribute the difference between the predicted uncertainties for Case (A) and (C) to the model uncertainty.

### 2D non-Newtonian flows

In this section, we test the capability of the present approach in non-Newtonian flows where we misspecify the constitutive relation as the Newtonian one, and then we employ DNNs to correct the misspecified models. In particular, we consider the steady flows in a

Figure 8: Time-dependent reaction-diffusion system in slices (\(t=0.5,0.75\)): Predicted uncertainty for \(u\) and \(\phi\) from B-PINNs with noisy data of \(u\) and \(f\) and differently specified reaction models. Representative predicted uncertainties for (a) \(u\), and (b) \(\phi\). Red dashed line: Predicted mean from B-PINNs; Black solid line: Reference solution. Our approach is able to improve the the predicted mean and the predicted uncertainty significantly. The decrease in the accuracy of predicted mean and increase in the predicted uncertainty of our approach are caused by the lack of knowledge over the reaction model.
2D channel and a 2D cavity, and the employed non-Newtonian fluids are power-law ones. The constitutive relation in a power-law fluid is expressed as \(\mu=\mu_{0}|\mathbf{S}|^{n-1}\), where \(\mu\) is the dynamic viscosity of a fluid, \(\mu_{0}\) is a constant, \(\mathbf{S}=\nabla\mathbf{u}+(\nabla\mathbf{u})^{T}\) denotes the shear rate (\(|\mathbf{S}|=\sqrt{(\mathbf{S}:\mathbf{S})/2}\)), \(\mathbf{u}\) is the fluid velocify, and \(n\) is the power law index. Without loss of generality, the fluids are shear-thinning (pseudoplastic, \(n<1\)) and shear-thickening (dilatant, \(n>1\)) in the channel and cavity flows, respectively.

#### 3.3.1 Channel flow

The following 2D non-Newtonian channel flow is considered here and the corresponding governing equation is expressed as [49]:

\[\begin{split}&\frac{\partial}{\partial y}(\mu(y)\frac{\partial u}{ \partial y})-\frac{\partial P}{\partial x}=f(y),\ y\in[-H/2,H/2],\\ & u(-H/2)=u(H/2)=0,\end{split} \tag{12}\]

where \(H=1\), \(\frac{\partial P}{\partial x}=c\) is a constant, and \(f(y)=0\). The analytical solution reads as

\[u(y)=\frac{n}{n+1}(-\frac{1}{\mu_{0}}\frac{\partial P}{\partial x})^{1/n}[( \frac{H}{2})^{1+1/n}-|y|^{1+1/n}]. \tag{13}\]

For this specific case, we set \(n=0.25,c=-1,H=1,\mu_{0}=0.5\) to generate the reference solution as well as training data.

In the first test case, we assume that we have 30 clean measurements of \(u\) and 51 clean measurements of \(f\), randomly and uniformly sampled from \([-H/2,H/2]\), respectively. We assume the fluid is a Newtonian one and the viscosity is an unknown constant. We then employ the PINN method to identify the unknown viscosity \(\mu_{1}\) in the following equation:

\[\frac{\partial}{\partial y}(\mu_{1}\frac{\partial u}{\partial y})-\frac{ \partial P}{\partial x}=f(y), \tag{14}\]

which describes Newtonian flows. Specifically, the boundary condition is hard-encoded in the equation for \(u\) so that we only have two terms in the PINN loss function (2), i.e., \(\mathcal{L}_{\mathcal{D}_{u}}\) for the data of \(u\) and \(\mathcal{L}_{PDE}\) for the equation. As shown in Fig. 2(b), we cannot minimize \(\mathcal{L}_{\mathcal{D}_{u}}\) and \(\mathcal{L}_{PDE}\) simultaneously in PINNs, regardless of the choice of belief weights for those two terms, indicating that the NN surrogate model \(u_{\theta}\) is not able to fit the non-Newtonian data and satisfy the Newtonian physical law at the same time.

Similar as in previous two examples in this section, we can add a DNN as the correction in Eq. (14) to address the model misspecification issue in PINNs. The equation with the correction term can now be expressed as follows:

\[\frac{\partial}{\partial y}(\mu_{1}\frac{\partial u}{\partial y})-\frac{ \partial P}{\partial x}+s(y)=f(y), \tag{15}\]

where \(s\) denotes the correction term. Here we treat \(\mu_{1}\) as known and set \(\mu_{1}=0.1\) since we have a DNN to model the discrepancy. We present the results from the PINNs with a correction in Fig. 9(a), and we can see that the current predictions for \(u\) and \(f\) agree with the reference solutions much better than in the case where the model is misspecified.
We note that in this specific case we can also employ other approaches to address the issue that we cannot fit the data and satisfy the equation simultaneously in PINNs if the physical model is misspecified. A natural way is to treat the viscosity in Eq. (14) as an unknown function and employ a DNN to directly learn the constitutive relation. We present the corresponding results in Fig. 9(b), in which we employ the same setup, e.g., training data, NN architecture in PINNs, as used in Fig. 2 as well as Fig. 9(a). As shown, PINNs with a DNN for learning the viscosity now can fit \(u\) and \(f\) much better than the results in Fig. 2. We then compare the these results with the ones in Fig. 9(a). As shown, the predictions for \(u\) are similar, but the PINN with a DNN as correction is able to provide more accurate predictions for \(f\), i.e., less fluctuations are observed in Fig. 9(a) than in Fig. 9(b). Furthermore, we illustrate the loss histories for training PINNs in these two approaches in Fig. 9(c). As we can see, the training loss decreases faster and converges to smaller values in the case where we add a DNN as correction to the misspecified model, as compared to the approach in which we model the viscosity as an unknown function. We conjecture that the superior performance of the proposed approach can be attributed to the following reason: modeling the viscosity as an unknown function via a DNN requires the accurate estimate of the derivatives in Eq. (14), which is quite challenging since \(n<1\) results in very sharp gradients for the viscosity around \(y=0\). However, the correction term, i.e., \(s\) in Eq. (14), is a quite smooth function, as shown in Fig. 10. Hence, we can expect the easier training as well as better accuracy in the proposed method.

We now move to the noisy-data case and employ B-PINNs to quantify the model uncertainty. The same 30 noisy measurements of \(u\) and 51 noisy measurements of \(f\) are sampled

Figure 9: 2D non-Newtonian channel flow: Correcting the misspecified physics with two different approaches, i.e. (a) modeling the discrepancy in the equation \(s(y)\) with a DNN, and (b) modeling the unknown space-dependent viscosity \(\mu(y)\) with a DNN. In (c), the PINN training losses are presented and they converge to smaller values for both \(\mathcal{L}_{\mathcal{D}_{u}}\) and \(\mathcal{L}_{PDE}\) in the case where we add a DNN as correction as compared to the one in which we use a DNN to approximate the viscosity.
but now they are corrupted by additive Gaussian noise with \(0.005\) and \(0.05\) noise scales, respectively. In B-PINNs, we test two specific cases, i.e., (1) we assume the fluid is Newtonian and encode Eq. (14) in B-PINNs. The objective is then to infer the unknown viscosity given data, similar as the case in Fig. 2; and (2) we are not sure if the fluid is Newtonian or not, we then add a DNN as the correction to the possible misspecified model, i.e., Eq. (15). Also, we note that in this specific case \(s\) can be computed analytically by plugging the analytic solution of the power-law fluid into Eq. (15): \(s(y)=4\mu_{1}c^{4}|y|^{3}/\mu_{0}^{4}+c\), which serves as the reference in what follows.

The results from B-PINNs are illustrated in Fig. 10. As observed, the NN surrogate \(u_{\theta}\) is not able to fit the data and satisfy the misspecified physical model at the same time (Fig. 10(a)). In addition, with the additional DNN for correcting the misspecified physical model in the previous case, \(u_{\theta}\) can now fit the data and satiafy the corrected physical model simultaneously, as shown in Fig. 10(b). Also, we can see that the computational errors between the predicted mean and the reference solution for \(s\) are bounded by the predicted uncertainties in the entire computational domain.

#### 3.3.2 Cavity flow

Next, we consider steady non-Newtonian flow in a two-dimensional lid-driven square cavity (i.e., \(x,y\in[0,1]\)), which can be described by the incompressible Navier-Stokes equations

Figure 10: 2D non-Newtonian channel flow: Predictions for \(u\), \(f\) and the discrepancy \(s\) from B-PINNs with noisy data. B-PINNs with (a) misspecified model, and (b) misspecified model plus a correction.
as

\[\nabla\cdot\mathbf{u} =0, \tag{16a}\] \[\mathbf{u}\cdot\nabla\mathbf{u} =-\nabla P+\nabla\cdot[\nu(\nabla\mathbf{u}+(\nabla\mathbf{u})^{T})]+\mathbf{f}, \tag{16b}\]

where \(\mathbf{u}=(u,v)\) denotes the velocity in \(x-\) and \(y-\) directions, respectively; \(P\) is the pressure; and \(\nu\) is the kinematic viscosity. In particular, the boundary condition for the upper wall is expressed as

\[u=U\left(1-\frac{\cosh[r(x-\frac{L}{2})]}{\cosh(\frac{rL}{2})}\right),\;v=0, \tag{17}\]

where \(U\), \(r\) and \(L\) are constant. Specifically, \(r=10\), \(L=1\) is the length of the cavity. In addition, the remaining walls are stationary. Here, the viscosity model is the same as the one used in channel flow, i.e., \(\mu=\rho\nu=\mu_{0}|\mathbf{S}|^{n-1}\). The objective is to learn the governing equations, which explain the data given measurements on \(\mathbf{u}\), \(P\), and \(\mathbf{f}\). In particular, the training data as well as the reference solution in this case are generated by solving Eq. (16a) using the lattice Boltzmann equation model (LBE) in [49] with \(n=1.5\). Details for the computations in LBE can be found in C.

Similar as in Sec. 3.3.1, suppose we misspecify the governing equations as the ones for Newtonian flows, and we then employ the present approach to correct the model misspecification as follows:

\[\nabla\cdot\mathbf{u} =0, \tag{18a}\] \[\mathbf{u}\cdot\nabla\mathbf{u} =-\nabla P+\nu_{0}\nabla^{2}\mathbf{u}+\mathbf{f}+\mathbf{s}, \tag{18b}\]

where \(\mathbf{s}=(s_{x},s_{y})\) are the two components correcting the momentum equations in \(x-\) and \(y-\) directions, respectively, and \(f_{x}(x,y)=f_{y}(x,y)=0\) for \(x,y\in[0,1]\). The continuity equation Eq. (18a) is correctly specified and hence does not need correction.

We assume that we have 400 random measurements of \(u\), 400 random measurements of \(v\), and \(81\times 81\) uniformly distributed measurements of \(P\) from the reference solution. Here we only consider the case with clean data and the ensemble PINNs are employed for quantifying the model uncertainty. We would like to mention that we do not use B-PINNs here since we need to utilize several hundred training data to achieve meaningful results, which is computationally expensive as the HMC is used for posterior estimation. Similarly, we test two specific cases, i.e., (1) the viscosity in Eq. (18b) is an unknown constant and we do not use the correction \(\mathbf{s}\) in PINNs, and (2) we have an initial guess for the viscosity in Eq. (18b) and we add \(\mathbf{s}\) in PINNs for the correction for the possible misspecification. Specifically, we use one DNN to model the velocity \(\mathbf{u}\), one DNN to approximate the pressure \(p\), and one DNN to model the discrepancy in Eq. (18b), i.e. \(\mathbf{s}\). In particular, we set \(\nu_{0}=0.0001\) as the initial guess in Case (2) which will not affect the computational accuracy in PINNs since the model discrepancy is already leveraged by \(\mathbf{s}\).

We present the predictions for \(\mathbf{u}\) in Fig. 11. As we can see, PINNs with misspecified model lead to inaccurate inference for both \(u\) and \(v\) (Fig. 11(a)). The proposed approach, however, achieves much better results by adding a DNN to correct the model misspecification. We further illustrate the computational errors for both cases in Table 4, which clearly shows the effectiveness as well as superiority of the present method. In addition, we show
some representative 1D slices for the predicted \(u,v,f_{x}\), and \(f_{y}\) with uncertainties in Fig. 12. Similarly, the PINNs with misspecified model cannot fit the data and satisfy the physical model at the same time (Fig. 12(a)), while this issue can be handled well as we add a DNN to correct the misspecified model, as shown in Fig. 12(b).

\begin{table}
\begin{tabular}{c|c|c|c|c} \hline \hline  & Error of \(u\) & Error of \(v\) & Error of \(f_{x}\) & Error of \(f_{y}\) \\ \hline Misspecified model & 7.79\% & 11.28\% & \(3.9153\times 10^{-4}\) & \(3.1466\times 10^{-4}\) \\ \hline Misspecified model with correction & 1.33\% & 2.29\% & \(7.8389\times 10^{-7}\) & \(8.2283\times 10^{-7}\) \\ \hline \hline \end{tabular}
\end{table}
Table 4: 2D non-Newtonian Cavity flow: Errors of the predicted means from ensemble PINNs. The error of \(u/v\) is defined as the relative \(L_{2}\) error, while the error of \(f_{x}/f_{y}\) the mean squared error since the reference solutions are \(f_{x}=0\) and \(f_{y}=0\).

Figure 11: Non-Newtonian cavity flow: Predictions from ensemble PINNs for \(u\) and \(v\) in contours. Red dashed line: Predicted mean; Black solid line: Reference solution.

Figure 12: Non-Newtonian cavity flow: Predictions for \(u,v,f_{x},f_{y}\) in 1D slices. Red dashed line: predicted mean; black solid line: reference solution; uncertainty bound: predicted uncertainty.
## 4 Summary

We presented a general approach to improve the accuracy of physics-informed neural networks (PINNs) for discovering governing equations from data when the physical model is misspecified. In particular, we first assume the mathematical models for the specific problem we considered and encode them in PINNs via automatic differentiation. We note that the assumed models may be misspecified in certain complex systems since not all the physical processes are fully understood here. We then utilize another deep neural networks (DNN) to model the discrepancy between the imperfect models and the observational data in order to correct the misspecified model. In addition, to quantify uncertainties arising from noisy and/or gappy data, we employ two typical uncertainty-induced PINNs methods, i.e., the Bayesian physics-informed neural networks (B-PINNs) and ensemble PINNs. A series of numerical examples are conducted to show the effectiveness of the proposed method. Specifically, we considered (1) a reaction-diffusion system where we misspecified the reaction model, and (2) non-Newtonian flows in two-dimensional channel and cavity but misspecified the fluid viscosity. The results show that the present approach can achieve good accuracy since the added DNN is able to correct the misspecified physical models encoded in PINNs. Also, the B-PINNs and/or ensemble PINNs can provide reasonable uncertainty bounds for the discovered physical models, which reflects the uncertainties in physical models, i.e., _model uncertainty_. Furthermore, we showcase that we can seamlessly combine the proposed method with the symbolic regression to obtain explicit governing equations based on the trained NNs. Finally, we expect the proposed approach to be a promising tool in a wider class of applications for solving problems where the physical processes are not fully understood.

## Acknowledgement

This work was supported by the MURI/AFOSR project (FA9550-20-1-0358), the DOE-MMICS SEA-CROGS project (DE-SC0023191), the NIH grant Neural Operator Learning to Predict Aneurysmal Growth and Outcomes (R01HL168473), and the NIH grant CR-CNS: Waste-clearance flows in the brain inferred using physics-informed neural networks (R01AT012312). Z.Z. thanks Dr. Khemraj Shukla from Brown University for helpful discussion.
