A Nonlinear, Time-Evolving Statistical FEM

A Gaussian process (GP), $$ {\xi }_{\theta }$$, is introduced inside of the governing equations, which represent an unknown forcing process in space and time, with time-varying parameters $$ \theta $$. For a general nonlinear PDE, this is given by

Coefficients of Eq. 2 are assumed to be known, and we choose to update the numerical solution to the model, acknowledging that estimating $$ \mathrm{\Lambda }$$ [using, e.g., maximum likelihood methods (20), Markov chain Monte Carlo (21), or inversion methods in general (22)] is also of utmost interest.

Discretizing Eq. 2 using finite elements in space with an implicit or explicit Euler method in time,^† denote by $$ {\mathbf{u}}^n\in {\mathbb{R}}^M$$ the FEM coefficients at time $$ n\mathrm{\Delta }t$$, for time step $$ \mathrm{\Delta }t$$. Further analysis of this discretization will be the focus of future work. For a potentially nonlinear system of equations $$ {\mathcal{F}}_{\mathrm{\Lambda }}$$, the resultant system can be expressed as (full construction is given in SI Appendix, sections 2 and 3)

Unlike the elliptic example in the Introduction, for Eq. 2 the induced probability measure on the FEM coefficients is not available in closed form, and we present two approximations, based on the extended Kalman filter (EKF) and the ensemble Kalman filter (EnKF). The first linearizes about the current solution with the Jacobian $$ \mathbf{J}$$ of the nonlinear $$ {\mathcal{F}}_{\mathrm{\Lambda }}$$ (evaluated at the current solution) to give a Gaussian approximation of $$ p\left({\mathbf{u}}^n|{\mathit{\theta }}_{1:n},\mathrm{\Lambda }\right)\approx \mathcal{N}\left({\mathbf{m}}_n,{\mathbf{C}}_n\right)$$. The second uses an ensemble in which a perturbed system is solved, with realizations from $$ {\mathbf{e}}_n$$, at each time step. Summary statistics (e.g., mean, covariance) are then computed from this ensemble.

For the prior, the deterministic FEM solution is identically equal to the mean in the EKF approach. However we have found in numerical experiments that the EKF method, due to the use of the Jacobian, inflates the covariance at points of high gradient with reduction at points approaching near-zero gradient, when using large time steps. This does not occur with the EnKF approach.

Algorithm 1: EKF algorithm
for $$ n\le n_t$$ do
(Prediction step)
Solve $$ \mathcal{F}\left({\mathbf{m}}_{n|n-1},{\mathbf{m}}_{n-1|n-1}\right)=\mathbf{0}$$.
$$ {\widehat{\mathbf{C}}}_{n|n-1}={\left({\mathbf{J}}^n\right)}^{-1}\left({\mathbf{J}}^{n-1}{\mathbf{C}}_{n-1|n-1}{\left({\mathbf{J}}^{n-1}\right)}^{\top }\right){\left.mb{J}^n\right)}^{-\top }$$
Estimate:
$$ {\text{arg max}}_{{\mathit{\theta }}_n,{\sigma }_n}\left\{\log p\left({\mathbf{y}}_n|{\mathbf{y}}_{1:n-1},{\mathit{\theta }}_{1:n},{\sigma }_{1:n}\right)\right.$$
$$ \left.+\log p\left({\mathit{\theta }}_n\right)+\log p\left({\sigma }_n\right)\right\}$$.
$$ {\mathbf{C}}_{n|n-1}={\widehat{\mathbf{C}}}_{n|n-1}+{\mathrm{\Delta }}_t{\left({\mathbf{J}}^n\right)}^{-1}\mathbf{G}\left({\mathit{\theta }}_n\right){\left.{\mathbf{J}}^n\right)}^{-\top }$$.
$$ {\mathbf{S}}_n={\mathbf{H}}_n{\mathbf{C}}_{n|n-1}{\mathbf{H}}_n^{\top }+{\sigma }_n^2\mathbf{I}$$.
(Analysis step)
$$ {\mathbf{m}}_{n|n}={\mathbf{m}}_{n|n-1}+{\mathbf{C}}_{n|n-1}{\mathbf{H}}_n^{\top }{\mathbf{S}}_n^{-1}\left({\mathbf{y}}_n-{\mathbf{H}}_n{\mathbf{m}}_{n|n-1}\right)$$.
$$ {\mathbf{C}}_{n|n}={\mathbf{C}}_{n|n-1}-{\mathbf{C}}_{n|n-1}{\mathbf{H}}_n^{\top }{\mathbf{S}}_n^{-1}{\mathbf{H}}_n{\mathbf{C}}_{n|n-1}$$.

Conditioning on Data.

Data $$ {\mathbf{y}}_n\in {\mathbb{R}}^N$$ are observed at time $$ n\mathrm{\Delta }t$$ on the grid $$ {\mathbf{x}}_{\text{obs}}$$. These data are corrupted with noise $$ {\mathit{\eta }}_n\sim \mathcal{N}\left(\mathbf{0},{\sigma }_n^2\mathbf{I}\right)$$ independent to the model $$ {\mathbf{u}}^n$$ to give the data-generating process $$ {\mathbf{y}}_n={\mathbf{H}}_n{\mathbf{u}}^n+{\mathit{\eta }}_n$$, where the linear observation operator $$ {\mathbf{H}}_n:{\mathbb{R}}^M\to {\mathbb{R}}^N$$ maps from the computed solution grid to the observation grid, using the FEM interpolant.

The filtered distribution $$ p\left({\mathbf{u}}^n|{\mathbf{y}}_{1:n},{\mathit{\theta }}_{1:n},{\sigma }_{1:n},\mathrm{\Lambda }\right)$$, where $$ {\mathbf{y}}_{1:n}=\left({\mathbf{y}}_1,{\mathbf{y}}_2,\dots ,{\mathbf{y}}_n\right)$$, is our primary object of interest.^‡ We take a Bayesian interpretation and refer to this as the filtered posterior distribution or just the posterior, when the context is clear. However, as this is a filtering problem, non-Bayesian methods are perfectly valid. We assume that all distributions are Gaussian, so the posterior can be computed with standard methods in data assimilation (23); we use the EKF (24) and the EnKF (stochastic form) (25).

The initial conditions are known (i.e., they are given a Dirac measure). For time $$ n$$, we make a tentative prediction step according to the PDE model $$ {\mathcal{F}}_{\mathrm{\Lambda }}$$, propagating uncertainty in the previous time step [described by $$ p\left({\mathbf{u}}^{n-1}|{\mathbf{y}}_{1:n-1},{\mathit{\theta }}_{1:n-1},{\sigma }_{1:n-1}\right)$$], to give the prediction measure $$ p\left({\mathbf{u}}^n|{\mathbf{y}}_{1:n-1},{\mathit{\theta }}_{1:n-1},{\sigma }_{1:n-1}\right)$$. Parameters $$ \left({\mathit{\theta }}_n,{\sigma }_n\right)$$ are then estimated, and the full prediction step is completed to estimate $$ p\left({\mathbf{u}}^n|{\mathbf{y}}_{1:n-1},{\mathit{\theta }}_{1:n},{\sigma }_{1:n}\right)$$. Data observed at time $$ n{\mathrm{\Delta }}_t$$ are then conditioned on to give the updated filtering distribution $$ p\left({\mathbf{u}}^n|{\mathbf{y}}_{1:n},{\mathit{\theta }}_{1:n},{\sigma }_{1:n}\right)$$.

We assume the parameters are independent across time [i.e., $$ p\left({\mathit{\theta }}_n|{\mathit{\theta }}_{n-1}\right)=p\left({\mathit{\theta }}_n\right)$$]. Parameters may also be time constant, which is discussed in the maximum likelihood setting in ref. 26 and in the hierarchical Bayesian setting in ref. 27. SI Appendix also contains a possible modification of the method that accounts for time-invariant parameters. The following procedure provides an overview of the method, and Algorithms 1 and 2 give pseudocode versions.

Simulation Study.

We condition on data generated from an extended Korteweg–de Vries (eKdV) equation with a cubic nonlinear term:

$w_t+\alpha w{w}_x+\beta w_{xxx}+\epsilon w^3{w}_x=0,$

Algorithm 2: EnKF algorithm
for $$ n\le n_t$$ do
(Prediction step)
for $$ i\le N_{ens}$$ do
Solve $$ \mathcal{F}\left({\mathbf{u}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{n,\left[i\right]},{\mathbf{u}}^{n-1,\left[i\right]}\right)=0$$.
Compute $$ {\mathbf{m}}_{n|n-1}$$ and $$ {\widehat{\mathbf{C}}}_{n|n-1}$$ from $$ \left\{{\mathbf{u}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{n,\left[i\right]}\right\}$$
Estimate:
$$ {\text{arg max}}_{{\mathit{\theta }}_n,{\sigma }_n}\left\{\log p\left({\mathbf{y}}_n|{\mathbf{y}}_{1:n-1},{\mathit{\theta }}_{1:n},{\sigma }_{1:n}\right)\right.$$
$$ \left.+\log p\left({\mathit{\theta }}_n\right)+\log p\left({\sigma }_n\right)\right\}$$
for $$ i\le N_{ens}$$ do
Solve $$ \mathcal{F}\left({\mathbf{u}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{n,\left[i\right]},{\mathbf{u}}^{n-1,\left[i\right]}\right)+{\mathbf{e}}_{n-1}^{\left[i\right]}=0$$.
Compute $$ {\mathbf{m}}_{n|n-1}$$ and $$ {\widehat{\mathbf{C}}}_{n|n-1}$$ from $$ \left\{{\mathbf{u}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{n,\left[i\right]}\right\}$$
$$ {\mathbf{S}}_n={\mathbf{H}}_n{\mathbf{C}}_{n|n-1}{\mathbf{H}}_n^{\top }+{\sigma }_n^2\mathbf{I}$$.
(Analysis step)
for $$ i\le N_{ens}$$ do
$$ {\mathbf{u}}^{n,\left[i\right]}={\mathbf{u}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{n,\left[i\right]}+{\mathbf{C}}_{n|n-1}{\mathbf{H}}_n^{\top }{\mathbf{S}}_n^{-1}\left({\mathbf{y}}_n+{\mathit{\eta }}_n^{\left[i\right]}-{\mathbf{H}}_n{\mathbf{u}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{n,\left[i\right]}\right)$$

setting $$ \alpha =1$$, $$ \beta =0.01$$, and $$ \epsilon =20$$. The misspecified KdV model

$u_t+\alpha u{u}_x+\beta u_{xxx}+{\xi }_{\theta }=0$ [3]

has the same coefficient values, with initial conditions set to a wave of depression $$ u\left(x,0\right)=-0.3{\text{sech}}^2\left(x-15\right)$$ on the space–time grid $$ \left(x,t\right)\in \left[\mathrm{0,20}\right]\times \left[\mathrm{0,50}\right]$$. Boundary conditions are periodic. Gaussian random forcing $$ {\xi }_{\theta }$$ has spatial covariance $$ k_{\theta }\left(x,x^{\prime }\right)={\tau }^2\exp \left(-\parallel x-x^{\prime }{\parallel }^2/\left(2{\ell }^2\right)\right)$$ (we refer to $$ \tau $$ as the scale parameter and $$ \ell $$ as the length parameter).

KdV is discretized following ref. 29, using $$ P_1$$ trial functions and $$ P_0$$ testing functions, with a Crank–Nicolson method in time. The data-generating process is simulated using Dedalus (30) with 1,024 grid points in space. This is then down sampled to 20 grid points and jittered with synthetic Gaussian observational noise (mean 0, variance $$ 0.001^2$$) to give the simulated dataset.

We assume $$ {\mathbf{y}}_n=\mathbf{H}{\mathbf{u}}^n+{\mathit{\eta }}_n$$, where $$ {\mathbf{u}}^n$$ is the Galerkin discretized solution to Eq. 3 and $$ {\mathit{\eta }}_n\sim \mathcal{N}\left(0,{\sigma }_n^2\mathbf{I}\right)$$. Hyperparameters $$ {\mathit{\theta }}_n=\left({\tau }_n,{\ell }_n\right)$$ and noise level $$ {\sigma }_n$$ are estimated at each step by maximizing the log-marginal posterior, with the weakly informative truncated Gaussian priors $$ {\tau }_n\sim {\mathcal{N}}_{+}\left(1,1^2\right)$$, $$ {\ell }_n\sim {\mathcal{N}}_{+}\left(1,1^2\right)$$, and $$ {\sigma }_n\sim {\mathcal{N}}_{+}\left(0,1^2\right)$$. The EnKF is used in this section, with $$ n_x=400$$, $$ n_t=\mathrm{2,001}$$, and $$ N_{ens}=400$$. Results are presented in Fig. 1.

Fig. 1.

Simulation study results (using EnKF). A shows a prior solution with fixed hyperparameters $$ {\mathit{\theta }}_n$$ with the data-generating process (DGP; as labeled) for two times with accompanying $$ 95\%$$ credible intervals. Note the degree of model mismatch between the two profiles. B shows the DGP and the posterior mean and $$ 95\%$$ credible intervals for the same times as the prior. Given the data, the model is highly certain about the updated mean, for which the model mismatch has been corrected. This is now using estimated hyperparameters $$ {\mathit{\theta }}_n$$ for a fully data-driven approach. C shows the prior mean, DGP, and posterior mean over the entire simulation grid in space and time, and D shows the estimated parameters $$ {\mathit{\theta }}_n$$, $$ {\sigma }_n$$ across time, with the true value for $$ {\sigma }_n$$ shown as a dashed turquoise line.

For a fixed set of hyperparameters $$ {\tau }_n=0.0025$$, $$ {\ell }_n=1$$ for all $$ n$$, the data-generating process and estimated prior are shown in Fig. 1A, appearing visually mismatched in the mean via phase shift, increased oscillations, and increased wave interactions. Note that the stochastic forcing induces an uncertainty about the PDE solution, represented by the $$ 95\%$$ credible intervals shown. Note also that the data-generating process is approximately contained within the credible intervals.

Fig. 1B shows that the posterior mean approximates the data-generating process, and the posterior uncertainty bounds have shrunk as a result of conditioning, indicating high certainty about the posterior mean values. Model discrepancy between the data and the statFEM solution has been corrected for. The space–time view of the posterior, shown in Fig. 1C, shows that the posterior has incorporated the complex soliton interactions in the data, not present in the prior.

Parameter estimates (Fig. 1D) indicate that the length and noise parameters are both stable, with the noise being slightly overestimated (i.e., $$ {\sigma }_n\approx 0.003>0.001$$). Times at which the noise is not identified result in it being set to the lower bound. The scale parameter quantifies the accuracy of the model prediction step at each time step. In this case, model predictions vary in their accuracy and appear approximately bounded to within $$ \left(10^{-5},10^{-1}\right)$$.

Case Study: Experimental Data

We now apply the method to the experimental data collected in ref. 31. Experiments were conducted to study weakly nonlinear models for internal waves in lakes and consisted of generating internal waves in a two-layer stratified system, inside of a clear acrylic tank of dimensions $$ 6\times 0.3\times 0.29\hspace{0.17em}\text{m}$$. The tank contained an upper layer of fresh water and a lower layer of saline water, with a density gradient of $$ \mathrm{\Delta }\rho =20\hspace{0.17em}{\text{kgm}}^{-3}$$. The tank was able to rotate in order to establish the initial conditions, which were an inclined plane of angle $$ 𝜗=0.5°$$. This initial condition mimics the shear induced by strong winds in lakes. At time $$ t=0$$, the tank is rotated to restore it to the horizontal.

Data were recorded at three spatially equidistant locations in the tank using ultrasonic wave gauges (Fig. 2), taking measurements approximately every $$ 0.01\hspace{0.17em}\text{s}$$, up to $$ T=\mathrm{1,000}\hspace{0.17em}\text{s}$$; we use data up to $$ T=300\hspace{0.17em}\text{s}$$. Data are measured in voltages and are postprocessed to give pynocline displacements in meters. These data are plotted in Fig. 3, where the small measurement error is visually apparent. Transient behavior is observed before steepening, and a soliton wave train forms; three such steepening events are observed in the data we analyze. As $$ T\to \mathrm{1,000}\hspace{0.17em}\text{s}$$, dissipation results in the wave profile approaching a flat steady-state profile.

Fig. 2.

Schematic diagram of the experimental apparatus. Wave gauges (WGs) are labeled WG1, WG2, and WG3, and the initial conditions are shown as a gray line, labeled with initial angle $$ 𝜗°$$.

Fig. 3.

Experimental data results; posterior computed by the EnKF. A shows the observed data, the deterministic solution to the misspecified model, and the posterior mean at the wave gauges with 95% credible intervals up to time $$ T=300\hspace{0.17em}\text{s}$$. Model mismatch between the model and data is visually apparent and has been offset in the posterior mean. B shows the posterior mean wave profile inside of the tank, with 95% credible intervals and the data, for three times. The wave profile has been reconstructed by conditioning on the observed data. C shows the estimated hyperparameters $$ {\mathit{\theta }}_n$$, and D shows the estimated posterior mean over the space–time grid.

Our physical model is an eKdV equation with—for computational simplicity—a linear dissipation term. We acknowledge that for laminar boundaries, other methods are preferred (31, 32). Including some form of dissipation is important as otherwise, the model becomes impractically mismatched by the end of the simulation. The eKdV is given by

Incorporation of reflective boundary conditions is done by solving the eKdV equation across the extended domain $$ \left[\mathrm{0,2}L\right]$$ with periodic boundary conditions and summing solutions in the (reflected) subdomains $$ \left[0,L\right]$$, $$ \left[L,2L\right]$$:

The eKdV model does not capture the observed behavior exactly. The model waves have higher velocity than the observations, and model amplitudes are slightly larger than observed amplitudes. It is conjectured that this is due to misparameterization of dissipation, but in any case, the model is misspecified. Rather than estimating the eKdV parameters using inversion techniques, we sequentially update the model with observations to give the posterior $$ p\left({\mathbf{u}}^n|{\mathbf{y}}_{1:n},{\mathit{\theta }}_{1:n},{\sigma }_{1:n}\right)$$.

As before, we assume $$ {\mathbf{y}}_n=\mathbf{H}{\mathbf{u}}^n+{\mathit{\eta }}_n$$ with known noise $$ {\mathit{\eta }}_n\sim \mathcal{N}\left(\mathrm{0,1.3588}\times 10^{-8}\mathbf{I}\right)$$. As we solve on the extended domain $$ \left[\mathrm{0,2}L\right]$$ and sum solutions, the observation operator $$ \mathbf{H}$$ is taken to be the sum of the appropriate function values given by our FEM interpolant (a linear operation). The observation points are unchanging, as is the solution mesh, so $$ \mathbf{H}$$ is constant in time. The hyperparameters of $$ {\xi }_{\theta }$$, $$ {\mathit{\theta }}_n=\left({\tau }_n,{\ell }_n\right)$$, must be estimated at each iteration by maximizing the log-marginal posterior. Due to small data in space (three observations each time step), we use a projection method to estimate hyperparameters. This linearly projects the predicted mean $$ {\mathbf{m}}_{n|n-1}$$ forward, estimated from the data points: $$ y_i^n=a_n+b_n{\mathbf{m}}_{n|n-1}\left(x_i\right)$$. Parameters $$ a_n,b_n$$ are estimated to give the best least squares linear projection from the prediction to the data. This gives a projected dataset, $$ {\tilde{\mathbf{y}}}_n$$, using the linear shift: $$ {\tilde{\mathbf{y}}}_n=a_n+b_n{\mathbf{m}}_{n|n-1}$$. The estimated hyperparameters are then given by $$ {\text{arg max}}_{{\mathit{\theta }}_n}\left\{\log p\left({\tilde{\mathbf{y}}}_n|{\mathbf{y}}_{1:n-1},{\mathit{\theta }}_{1:n},{\sigma }_{1:n}\right)+\log p\left({\mathit{\theta }}_n\right)\right\}$$, in which the observed data $$ {\mathbf{y}}_n$$ are replaced with the projected data $$ {\tilde{\mathbf{y}}}_n$$. We project to a grid of 100 points uniformly spaced across the solution grid. Note that this is only for the parameter estimation step, and we do not use this $$ {\tilde{\mathbf{y}}}_n$$ as the data in the analysis step.

We set weakly informative priors: $$ {\tau }_n\sim {\mathcal{N}}_{+}\left(1,1^2\right)$$ and $$ {\ell }_n\sim {\mathcal{N}}_{+}\left(1,1^2\right)$$. The posterior is computed using the ensemble method with $$ n_x=200$$, $$ n_t=\mathrm{1,001}$$, and $$ N_{ens}=\mathrm{2,048}$$. Results are shown in Fig. 3. The posterior mean values of $$ p\left({\mathbf{u}}^n|{\mathbf{y}}_{1:n},{\mathit{\theta }}_{1:n},{\sigma }_{1:n}\right)$$ at the wave gauges are shown in Fig. 3A and offer a close fit to the data in comparison with the eKdV solution. The credible intervals shrink about the data (compare Fig. 3B) and are not seen on the figure.

Posterior wave profiles are shown in Fig. 3B and demonstrate that given the data, the method is able to yield a sensible estimate for the underlying wave profile and is hence able to reconstruct the wave profile given sparse observations in space (Fig. 3D). Furthermore the provided uncertainty quantification is physically sensible, with bounds contracting about the data and expanding near the boundaries.

The hyperparameters, $$ {\mathit{\theta }}_n$$, are shown in Fig. 3C. The scale parameter is seen to vary between two distinct levels, indicating that the model predictions vary in their accuracy. The amplitude of these mismatch scales shows that in this case, model mismatch is a cumulative effect that takes some time before it is obviously occurring (Fig. 3A, eKdV solution). Repeated conditioning on data helps to mitigate these long timescale effects due to continual updating. A space–time view of the posterior mean wave profile is shown in Fig. 3D, demonstrating that the general behavior of the flow (e.g., reflective boundary conditions, dissipation, wave train formation) is indeed captured.

Conclusions

We present a data-driven approach to the FEM that assimilates observations into nonlinear, time-dependent PDEs by embedding model misspecification uncertainty in the governing equations and sequentially updating the discretized equations with observations in a filtering context. Examples presented using the KdV equation (and the additional systems studied in SI Appendix) demonstrate that the method can approximate the data-generating process to give an interpretable posterior distribution, which reconstructs the observed phenomena.

This work sets the foundation for future studies of embedding data within FEM models. The use of the underlying Kalman framework permits the drawing upon of ideas from high-dimensional data assimilation, which for the systems studied here, were not needed. Techniques from Bayesian inversion can also be used to provide uncertainty quantification for physical quantities of interest, which will also allow for more accurate prediction. Finally, we believe that the development of similar methodology for alternate discretizations (e.g., spectral methods) could also be of great benefit, allowing for even broader application.