SINDy-BVP: Sparse identification of nonlinear dynamics for boundary value problems

We develop a data-driven model discovery and system identification technique for spatially-dependent boundary value problems (BVPs). Specifically, we leverage the sparse identification of nonlinear dynamics (SINDy) algorithm and group sparse regression techniques with a set of forcing functions and corresponding state variable measurements to yield a parsimonious model of heterogeneous material systems. The technique models forced systems governed by linear or nonlinear operators of the form $L\left[u\right(x\left)\right]=f\left(x\right)$ on a prescribed domain $x\in [a,b]$. We demonstrate the approach on a range of example systems, including Sturm-Liouville operators, beam theory (elasticity), and a class of nonlinear BVPs. The generated data-driven model is used to infer the governing operator and spatially-dependent parameters that describe the heterogenous, physical quantities of the system. Our SINDy-BVP framework enables the characterization of a broad range of systems, including for instance, the discovery of anisotropic materials with heterogeneous variability.

Figure 1

SINDy-BVP studies steady-state systems subjected to a forcing function. One simple example system is a beam clamped at both ends subjected to a static load (a). The beam deflects (b) in response to the load, and the forcing function and deflection are used for data-driven modeling via SINDy-BVP to learn the parametric coefficients (c) in the governing operator. The coefficients $p\left(x\right)$ and $q\left(x\right)$ vary spatially. The coefficients are directly related to the beam's spatially-varying mechanical properties. The grey boxes in (d) indicate that error can occur in the learned coefficients near the boundaries.

Figure 2

Overview of constructing the data sets and regression for each discrete spatial point in the data set. Part (a) shows a collection of trials, where different forcings ($f_j\left(x\right)$) are applied to a system yielding different responses ($u_j\left(x\right)$). A matrix containing a library of candidate terms $\mathit{\Theta }(\mathbf{U},\mathbf{F})$ is produced for each trial, where the rows are each a spatial point $x_k$ and each column contains a candidate function, as seen in part (b). A sliding procedure is used to select rows of a single $x_k$ from the libraries in part (b) to produce an aggregated library ${\mathit{\Theta }}^{\left(k\right)}$ for each $x_k$ in the data set. In (c), the regression is performed for each $x_k$ to produce a vector of the parametric coefficients $p\left(x\right)$ and $q\left(x\right)$ at each $x_k$. The regression in (c) is a Ridge regression algorithm. However, a sparsity constraint is applied by an iterative thresholding and grouping mechanism of the algorithm.

Figure 3

Summary of the models and operators studied with SINDy-BVP. The center column shows three example trials from the training data set ($u_1\left(x\right)$, $u_2\left(x\right)$, and $u_3\left(x\right)$). The solutions are all of $\mathcal{O}\left(1\right)$. The final column shows the parametric coefficients in the operator, as well as the inferred parameters for clean data. Coefficients $p\left(x\right)$ and $q\left(x\right)$ are plotted with an offset, with markers every 30 points.

Figure 4

SINDy-BVP succeeds at operator identification in signals with high SNR. A collection of trials with varying SNR is used for data-driven operator identification. The plots in (a) and (b), from top to bottom, are the linear Sturm-Liouville, nonlinear Sturm-Liouville, Poisson, and Euler-Bernoulli beam models. Relative loss [Eq. (10)] and a count of spurious terms in the identified model are quantified. In (a), 200 trials are used for regression and the SNR is varied. In (b), the number of trials is varied while holding a constant SNR of 100.

Figure 5

Parametric coefficient estimation in noisy data. With a known operator $L$, SINDy-BVP can estimate the parametric coefficients. Coefficient error is used to compare the effect of varying SNR (a) and the number of trials used as input data (b). In (a), the number of trials used as input data are held constant at 200 and in (b) the SNR is fixed at 100. The plots in (a) and (b), from top to bottom, are the linear Sturm-Liouville, nonlinear Sturm-Liouville, Poisson, and Euler-Bernoulli beam models.

Figure 6

Determination of correct order equation to use for building the data-driven model. This example is the Euler-Bernoulli beam theory, which should use the left hand side term $d^{4}u\left(x\right)/d{x}^4$ to build the correct model. The model that best captures the relationship $Lu\left(x\right)=f\left(x\right)$ is determined from the validation error $\parallel L\left[{\mathbf{U}}_j\right]-{\mathbf{F}}_j{\parallel }_2$ from a test data set. The fourth-order model ($d^{4}u\left(x\right)/d{x}^4$) exhibits the lowest error.