# Rheo-SINDy: Finding a Constitutive Model from Rheological Data for Complex
Fluids Using Sparse Identification for Nonlinear Dynamics

Takeshi Sato1,∗,†, Souta Miyamoto2,†, and Shota Kato3 1Institute for Chemical
Research, Kyoto University, Uji 611-0011, Japan  2Department of Chemical
Engineering, Graduate School of Engineering, Kyoto University, Kyoto 615-8510,
Japan  3Graduate School of Informatics, Kyoto University, Kyoto 606-8501,
Japan<EMAIL_ADDRESS>
###### Abstract

Rheology plays a pivotal role in understanding and predicting material
behavior by discovering governing equations that relate deformation and
stress, known as constitutive equations. Despite the critical importance of
constitutive equations in predicting dynamics of complex fluids, a systematic
methodology for deriving these equations from available data has remained a
significant challenge in the field. To overcome the problem, we propose a
novel method named Rheo-SINDy, which employs the sparse identification of
nonlinear dynamics (SINDy) for discovering constitutive models from
rheological data. Rheo-SINDy was applied to five distinct scenarios, including
four cases with well-established constitutive equations and one without
predefined equations. Our results demonstrate that Rheo-SINDy successfully
identifies accurate models for the known constitutive equations and derives
physically plausible approximate models for the scenario with the unknown one.
These findings validate the robustness of Rheo-SINDy in handling real-world
data complexities and underscore its efficacy as a powerful tool for advancing
the development of data-driven approaches in rheology.

  * ∗Author to whom any correspondence should be addressed.

  * †These authors contributed equally to this work.

  * Keywords: sparse identification, rheology, constitutive model

##  1 Introduction

Mathematical models grounded in physical laws are indispensable across science
and engineering, offering profound insights into the behavior of complex
systems. These models clarify the underlying mechanisms governing system
dynamics and empower predictions and innovations in technology and natural
science. Traditionally, model derivation has leaned heavily on theoretical and
empirical knowledge, often requiring expert knowledge and intuition. Data-
driven methods have become capable of assisting in developing mathematical
models and constructing models that provide advanced predictions [1]. These
data-driven methods involve the sparse identification methods [2, 3, 4, 5],
symbolic regression methods [6, 7, 8, 9, 10, 11], and physics-informed machine
learning methods [12, 13, 14, 15]. These methods have emerged as powerful
tools for deriving governing equations directly from data, overcoming the
limitations of conventional expert-dependent approaches.

Rheology is one of the scientific fields that address the properties of
flowing matter, which plays a crucial role in many industries, such as designs
of chemical processes, by providing insights into flow behaviors of complex
fluids. One of the roles of rheology is to discover or derive governing
equations that relate deformation and stress, referred to as constitutive
equations [16]. From an engineering perspective, accurate constitutive
equations are necessary to predict the flows of complex fluids under complex
boundary conditions. Nevertheless, it is generally difficult to theoretically
obtain constitutive equations for complex fluids. Instead, mesoscopic coarse-
grained models, which are based on molecular theories, have been explored in
the field of rheology. For example, for polymeric liquids, standard molecular
theories have been proposed [17, 18], on which refined mesoscopic models have
been constructed [19, 20, 21]. In these models, the motion of individual
(coarse-grained) molecules is numerically tracked. Although these models
require significantly more computational time compared to constitutive
equations, they can reproduce (nearly) accurate rheological data. Despite
these advancements, a clear methodology for obtaining constitutive equations
from available data remains elusive.

Data-driven methods are powerful approaches for addressing the aforementioned
challenges in rheology. Indeed, such methods have enhanced rheological studies
such as constitutive modeling, flow predictions of complex fluids, and model
selection [22, 23]. Some applications have successfully identified
constitutive relations of complex fluids or governing equations to predict the
dynamics of fluids with knowledge of rheology. These studies have employed
neural networks (NN), including deep NN [24], graph NN [25], recurrent NN
[26], physics-informed NN [27, 28, 29], multi-fidelity NN [30], and tensor
basis NN [31]. Gaussian processes (GP) have also been employed, for example,
for strain-rate dependent viscosity [32] or for viscoelastic properties [33,
34, 35, 36].

Despite the success of NNs and GPs, their black-box nature often obscures the
underlying physics, making symbolic regression techniques more appealing for
transparency and interoperability. These methods, such as the sparse
identification of nonlinear dynamics (SINDy) [2], have been frequently
utilized to track (reduced order) dynamics in the field of fluid mechanics
[37]. Inspired by these successes, symbolic regression methods have recently
started to be used in the field of rheology as well. For example, Mohammadamin
and coworkers [38] relied on SINDy for flexibly identifying the constitutive
equations of an elasto-visco-plastic fluid. Nevertheless, although there are
several attempts along this line, a comprehensive study to test SINDy for
rheological data has not yet been conducted.

In this study, we employ SINDy to find constitutive models from rheological
data, which we call as Rheo-SINDy. After verifying the performance of Rheo-
SINDy when the constitutive equations are known, we apply Rheo-SINDy to
problems where the constitutive equations are unknown. The details are shown
below.

##  2 Methods

Figure 1: Schematic illustration of Rheo-SINDy.

We use a data-driven method known as a sparse identification of nonlinear
dynamics (SINDy), which was originally developed by Brunton and coworkers [2].
In this study, we attempt to obtain constitutive equations of complex fluids
using SINDy. Here, we briefly explain the basic concepts of SINDy.

We consider dynamical systems generally expressed by the following
differential equation:

$\frac{{\rm d}{\bm{x}}(t)}{{\rm d}t}={\bm{f}}[\bm{x}(t)],$ (1)

where the vector ${\bm{x}}(t)$ represents the state of a system at time $t$
and the function $\bm{f}[\bm{x}(t)]$ determines the dynamics of the state
${\bm{x}}(t)$. The basic idea of SINDy is to find dominant terms for
describing the dynamics out of numerous candidates using a sparse
identification method. One can determine the (sparse) representation of
$\bm{f}$ by a dataset including a collection of $\bm{x}(t)$ and
$\dot{\bm{x}}(t)$ (the time derivative of $\bm{x}(t)$). The regression to
points of $\bm{x}(t)$ and $\dot{\bm{x}}(t)$ is computed with sparsity-
promoting techniques, such as $\ell_{1}$-regularization.

In the rheological community, it is of great importance to determine a
relation between stress and strain rate. This relation is a so-called
constitutive model or constitutive equation. Most constitutive equations are
differential equations that depend on the (extra) stress tensor $\bm{\tau}$
and velocity gradient tensor $\bm{\kappa}$. In this study, we prefer to use
the so-called extra stress tensor $\bm{\tau}$ as the stress tensor because
this stress tensor is $\bm{\tau}=\bm{0}$ at equilibrium, which is convenient
for SINDy regression. The total stress tensor $\bm{\sigma}$ can be obtained by
the relation $\bm{\sigma}=\bm{\tau}+G\bm{I}$, where $G$ is the modulus and
$\bm{I}$ is the unit tensor. A general form for such constitutive equations
can be written as

$\frac{{\rm d}{\bm{\tau}}(t)}{{\rm
d}t}={\dot{\bm{\tau}}}(t)={\bm{f}}[\bm{\tau}(t),\bm{\kappa}(t)].$ (2)

Here, $\bm{\kappa}(t)$ is a control variable during rheological measurements.
We use SINDy algorithm to find constitutive equations for complex fluids and
we refer to this technique as Rheo-SINDy.

The training data needed to Rheo-SINDy are transient stress data $\bm{T}$ and
those time derivatives $\dot{\bm{T}}$, which can be summarized as the
following two matrices:

${\bm{T}}=\left[\begin{array}[]{cccc}{\bm{t}}_{xx}&{\bm{t}}_{yy}&\cdots&{\bm{t}}_{zx}\\\
\end{array}\right]=\left[\begin{array}[]{cccc}\tau_{xx}(t_{1})&\tau_{yy}(t_{1})&\cdots&\tau_{zx}(t_{1})\\\
\tau_{xx}(t_{2})&\tau_{yy}(t_{2})&\cdots&\tau_{zx}(t_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
\tau_{xx}(t_{n})&\tau_{yy}(t_{n})&\cdots&\tau_{zx}(t_{n})\end{array}\right]$
(3)

and

${\dot{\bm{T}}}=\left[\begin{array}[]{cccc}{\dot{\bm{t}}}_{xx}&{\dot{\bm{t}}}_{yy}&\cdots&{\dot{\bm{t}}}_{zx}\\\
\end{array}\right]=\left[\begin{array}[]{cccc}{\dot{\tau}}_{xx}(t_{1})&{\dot{\tau}}_{yy}(t_{1})&\cdots&{\dot{\tau}}_{zx}(t_{1})\\\
{\dot{\tau}}_{xx}(t_{2})&{\dot{\tau}}_{yy}(t_{2})&\cdots&{\dot{\tau}}_{zx}(t_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
{\dot{\tau}}_{xx}(t_{n})&{\dot{\tau}}_{yy}(t_{n})&\cdots&{\dot{\tau}}_{zx}(t_{n})\end{array}\right],$
(4)

where ${\bm{t}}_{\mu\nu}$ ($\mu\nu\in\\{xx,yy,zz,xy,yz,zx\\}$) is the column
of ${\bm{T}}$, and we take the stress data for $n$ sequential times. The time
derivatives of the stress data $\dot{\bm{T}}$ are computed by a numerical
differentiation method. In this study, we apply the velocity gradient
$\bm{\kappa}$ to systems of prescribed constitutive equations or mesoscopic
models for viscoelastic fluids to take the stress data determined by the
states of the systems. The data of the velocity gradient tensor $\bm{K}$ are
summarized as

${\bm{K}}=\left[\begin{array}[]{cccc}{\bm{k}}_{xx}&{\bm{k}}_{yy}&\cdots&{\bm{k}}_{zx}\\\
\end{array}\right]=\left[\begin{array}[]{cccc}\kappa_{xx}(t_{1})&\kappa_{yy}(t_{1})&\cdots&\kappa_{zx}(t_{1})\\\
\kappa_{xx}(t_{2})&\kappa_{yy}(t_{2})&\cdots&\kappa_{zx}(t_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
\kappa_{xx}(t_{n})&\kappa_{yy}(t_{n})&\cdots&\kappa_{zx}(t_{n})\end{array}\right],$
(5)

where ${\bm{k}}_{\mu\nu}$ ($\mu,\nu\in\\{x,y,z\\}$) is the column of
${\bm{K}}$.

In Rheo-SINDy, we construct a library matrix of functions, denoted as
$\bm{\Theta}$, which can include various nonlinear functions. $\bm{\Theta}$ is
expressed as

${\bm{\Theta}}=\left[\begin{array}[]{ccccccc}{\bm{1}}&{\bm{T}}&{\bm{K}}&({\bm{T}}\otimes{\bm{T}})&({\bm{T}}\otimes{\bm{K}})&({\bm{K}}\otimes{\bm{K}})&\cdots\end{array}\right],$
(6)

where ${\bm{T}}\otimes{\bm{K}}$, for example, denotes all possible
combinations of the products of the row components in ${\bm{T}}$ and
${\bm{K}}$ for each time $t_{i}$ ($1\leq i\leq n$). We note that
${\bm{\Theta}}$ can incorporate not only polynomials but also other functions,
such as sinusoidal functions. Using these expressions, we can substitute Eq.
(2) as

$\dot{\bm{T}}={\bm{\Theta}}\bm{\Xi},$ (7)

where $\bm{\Xi}$ is the coefficient matrix written as

$\bm{\Xi}=\left[\begin{array}[]{cccc}{\bm{\xi}}_{xx}&{\bm{\xi}}_{yy}&\cdots&{\bm{\xi}}_{zx}\\\
\end{array}\right]=\left[\begin{array}[]{cccc}\xi_{xx,1}&\xi_{yy,1}&\cdots&\xi_{zx,1}\\\
\xi_{xx,2}&\xi_{yy,2}&\cdots&\xi_{zx,2}\\\ \vdots&\vdots&\ddots&\vdots\\\
\xi_{xx,N_{\bm{\Theta}}}&\xi_{yy,N_{\bm{\Theta}}}&\cdots&\xi_{zx,N_{\bm{\Theta}}}\end{array}\right].$
(8)

where $N_{\bm{\Theta}}$ is the total number of library functions.

To determine the coefficient matrix $\bm{\Xi}$, we solve the following
optimization problem:

$\hat{\bm{\xi}}_{\mu\nu}={\mathop{\mathrm{argmin}}\limits_{{\bm{\xi}}_{\mu\nu}}}\|\dot{\bm{t}}_{\mu\nu}-{\bm{\Theta}}{\bm{\xi}}_{\mu\nu}\|_{2}^{2}+R({\bm{\xi}}_{\mu\nu}),$
(9)

where $\hat{\bm{\xi}}_{\mu\nu}$ is the optimized sparse vector,
$||\cdots||_{2}$ is the $\ell_{2}$ norm defined as

$||\bm{x}||_{2}=\left(\sum_{i}x_{i}^{2}\right)^{1/2},$ (10)

and $R({\bm{\xi}}_{\mu\nu})$ is the regularization term. To obtain a sparse
solution of $\bm{\Xi}$, we apply the following five methods [37]: (i) the
sequentially thresholded least square algorithm (STLSQ), (ii) sequentially
thresholded Ridge regression (STRidge), (iii) least absolute shrinkage and
selection operator (Lasso), (iv) Elastic-Net (E-Net), and (v) adaptive-Lasso
(a-Lasso). These methods employ different regularization terms to obtain
sparse solutions (see Sec. S1 in the supporting information for detail). Each
method has a hyperparameter $\alpha$ to penalize the solution complexity,
which is to be tuned for obtaining good predictive yet parsimonious
representations. For this purpose, we test various $\alpha$ values and pick an
appropriate value of $\alpha$ that gives a small loss value and the (nearly)
correct number of terms for known constitutive equations.

In this study, we limit ourselves to shear rheological measurements that give
fundamental rheological properties. Under shear flow, among the components of
$\bm{\kappa}$, only $\kappa_{xy}$ has non-zero values. Here, $x$ is the
velocity direction, and $y$ is the velocity gradient direction. Since the
major stress components are $\tau_{xx}$, $\tau_{yy}$, $\tau_{zz}$, and
$\tau_{xy}$ under shear flow, we only use these components to conduct Rheo-
SINDy.

##  3 Case Studies

For case studies, we first test whether Rheo-SINDy can find appropriate
constitutive equations from training data generated by phenomenological
constitutive equations, namely the Upper Convected Maxwell (UCM) model and the
Giesekus model (the details are summarized in Sec. S2 in the supporting
information). Subsequently, we apply Rheo-SINDy to constitutive models for the
Dumbbell model, which is the most basic mesoscale model of viscoelastic
fluids. This section provides a brief overview of the models used in this
study and the conditions for creating the datasets.

###  3.1 Fundamental Equations of Dumbbell Models

Figure 2: Schematic illustration of the dumbbell model.

The dumbbell-based models have been widely utilized in numerous previous
studies for the computation of viscoelastic fluids and are considered a
standard mesoscopic model for viscoelastic fluids [39]. As illustrated in Fig.
2, a dumbbell consists of two beads (indexed as $1$ or $2$) and a spring that
connects them. The Langevin equations for the positions of the two beads
${\bm{r}}_{1/2}(t)$ can be written as

$\zeta\left[\frac{{\rm d}{\bm{r}}_{i}(t)}{{\rm
d}t}-{\bm{\kappa}}\cdot{\bm{r}}_{i}(t)\right]=-h(t)\left\\{{\bm{r}}_{i}(t)-{\bm{r}}_{j}(t)\right\\}+{\bm{F}}_{{\rm
B}i}(t),$ (11)

with $(i,j)=(1,2)$ or $(2,1)$. Here, $\zeta$ is the friction coefficient,
$h(t)$ is the spring strength, and ${\bm{F}}_{{\rm B}i}(t)$ is the Brownian
force acting on the bead $i$. The time evolution equation for the end-to-end
vector $\bm{R}(t)$ ($={\bm{r}}_{2}(t)-{\bm{r}}_{1}(t)$) of the beads is thus
obtained as

$\zeta\left[\frac{{\rm d}{\bm{R}}(t)}{{\rm
d}t}-{\bm{\kappa}}\cdot{\bm{R}}(t)\right]=-2h(t){\bm{R}}(t)+\left\\{{\bm{F}}_{{\rm
B}2}(t)-{\bm{F}}_{{\rm B}1}(t)\right\\}.$ (12)

The Brownian force is characterized by the first and second-moment averages as

$\langle{\bm{F}}_{{\rm B}i}(t)\rangle={\bm{0}},$ (13)

and

$\langle{\bm{F}}_{{\rm B}i}(t){\bm{F}}_{{\rm B}j}(t^{\prime})\rangle=2\zeta
k_{\rm B}T\delta_{ij}\delta(t-t^{\prime}){\bm{I}},$ (14)

where $k_{\rm B}$ is the Boltzmann constant and $T$ is the temperature. From
the end-to-end vector $\bm{R}(t)$, the stress tensor can be expressed as

${\bm{\tau}}(t)=\rho\langle h(t){\bm{R}}(t){\bm{R}}(t)\rangle-\rho k_{\rm
B}T{\bm{I}},$ (15)

where $\rho$ is the density of dumbbells.

There are several expressions for the spring strength $h(t)$. The most basic
one is the Hookean spring, defined as

$h(t)=h_{\rm eq}=\frac{3k_{\rm B}T}{n_{\rm K}b_{\rm K}^{2}},$ (16)

where $n_{\rm K}$ is the number of Kuhn segments per spring and $b_{\rm K}$ is
the Kuhn length. To reproduce several properties of polymers, such as shear
thinning under shear flow, it is essential to address finite extensible
nonlinear elastic (FENE) effects. Although the exact expression for FENE
springs is given by the inverse Langevin function, the following empirical
expression is widely used [39]:

$h(t)=h_{\rm eq}\frac{1-\langle{R}_{\rm eq}^{2}\rangle/R_{\rm
max}^{2}}{1-{\bm{R}}^{2}(t)/R_{\rm max}^{2}},$ (17)

where $\langle R_{\rm eq}^{2}\rangle^{1/2}=(n_{\rm K})^{1/2}b_{\rm K}$ is the
equilibrium length of the springs, and $R_{\rm max}=n_{\rm K}b_{\rm K}$ is the
maximum length of the springs. As shown later in Sec. 3.4, a constitutive
equation cannot be analytically obtained for the FENE dumbbell model. To
address the FENE spring more analytically, the following approximate
expression of the FENE spring has been proposed [39]:

$h(t)=h_{\rm eq}\frac{1-\langle{R}_{\rm eq}^{2}\rangle/R_{\rm
max}^{2}}{1-\langle{\bm{R}}^{2}(t)\rangle/R_{\rm max}^{2}}=h_{\rm eq}f_{\rm
FENE}(t).$ (18)

This spring is referred to as the FENE-P spring. Here, “P” means Peterlin, who
proposed the approximate form of the FENE spring law. The average appearing in
Eq. (18) makes it possible to obtain the analytical constitutive equation.

We use $\lambda=\zeta/4h_{\rm eq}$ as the unit time and $G=\rho k_{\rm B}T$ as
the unit stress for the dumbbell models. To simplify the expressions, in what
follows, we omit the tilde representing dimensionless quantities.

###  3.2 Hookean dumbbell model

The most basic dumbbell model is the Hookean dumbbell model, where Hookean
springs are employed (cf. Eq. (16)). From Eqs. (12), (15), and (16), the
Hookean dumbbell model reduces to the constitutive equation of the UCM model
(cf. Eq. (S5) in the supporting information) in the limit of $N_{\rm
p}\rightarrow\infty$ with $N_{\rm p}$ being the number of dumbbells.

For the Hookean dumbbell model, we generate training data by Brownian dynamics
(BD) simulations with the finite numbers of dumbbells ($N_{\rm
p}\in\\{10^{3},10^{4},10^{5}\\}$). We apply the oscillatory shear flow,
$\gamma(t)=\gamma_{0}\sin(\omega t)$, with $\gamma_{0}=2$ and $\omega=0.5$,
over a period from $t=0$ to $t=100$. The simulations are run with $\Delta
t=1\times 10^{-3}$ for $0\leq t\leq 100$ and data are collected at the
interval of $\Delta t_{\rm train}=1\times 10^{-2}$. Each simulation is
conducted with five different random seeds, and their average data is used for
training. Due to the characteristics of the BD simulation, the training data
inherently include noise originating from the finite $N_{\rm p}$. We here test
whether Rheo-SINDy can find from the noisy data the following equations for
the UCM model under shear flow (cf. Eqs. (S6)–(S8) in the supporting
information):

${\dot{\tau}}_{xx}=-\tau_{xx}+2\tau_{xy}\kappa_{xy},$ (19)
${\dot{\tau}}_{yy/zz}={\tau}_{yy/zz}=0,$ (20)
${\dot{\tau}}_{xy}=-\tau_{xy}+\kappa_{xy}+\tau_{yy}\kappa_{xy}=-\tau_{xy}+\kappa_{xy}.$
(21)

###  3.3 FENE-P dumbbell model

We next address the so-called FENE-P dumbbell model, where Eq. (18) is
utilized as the spring strength. As shown below, the FENE-P dumbbell model has
an analytical solution and is utilized for various flow problems, such as
turbulent flows [40].

Due to the assumption shown in Eq. (18), a simple representation of the time
evolution for the conformation tensor
${\bm{C}}=\langle{\bm{R}}(t){\bm{R}}(t)\rangle$ can be obtained as

$\frac{{\rm d}{\bm{C}}}{{\rm
d}t}-{\bm{C}}\cdot{\bm{\kappa}}^{+}-{\bm{\kappa}}\cdot{\bm{C}}=-f_{\rm
FENE}(t){\bm{C}}+\frac{n_{\rm K}}{3}{\bm{I}},$ (22)

where $\bm{\kappa}^{+}$ is the transposed $\bm{\kappa}$. The stress tensor is
thus obtained by

${\bm{\tau}}(t)=\rho h(t){\bm{C}}(t)-\rho k_{\rm B}T{\bm{I}}.$ (23)

Under shear flow, Eq. (22) reduces to the following expressions:

${\dot{C}}_{xx}=-f_{\rm FENE}C_{xx}+2C_{xy}\kappa_{xy}+\frac{n_{\rm K}}{3},$
(24) ${\dot{C}}_{yy/zz}=-f_{\rm FENE}C_{yy/zz}+\frac{n_{\rm K}}{3},$ (25)
${\dot{C}}_{xy}=-f_{\rm FENE}C_{xy}+C_{yy}\kappa_{xy}.$ (26)

Using Rheo-SINDy, we test whether or not Eqs. (24)–(26) can be discovered from
the data.

Although it has not been as widely recognized due to its complexity, the
FENE-P dumbbell model can also be expressed in the form of the constitutive
equation (i.e., the stress expression) [41]. From the textbook of Bird and
coworkers [39], the constitutive equation for the FENE-P model is

$\frac{{\rm d}{\bm{\tau}}}{{\rm
d}t}-{\bm{\tau}}\cdot{\bm{\kappa}}^{+}-{\bm{\kappa}}\cdot{\bm{\tau}}=-f_{\rm
FENE}(t){\bm{\tau}}+2{\bm{D}}+\frac{{\rm D}\ln Z}{{\rm
D}t}({\bm{\tau}}+{\bm{I}}),$ (27)

where ${\rm D}(\cdots)/{\rm D}t$ is the substantial derivative and $Z$ is the
function expressed as

$Z=\frac{1}{1-\langle{\bm{R}}^{2}(t)/R_{\rm
max}^{2}\rangle}=1+\frac{1}{3n_{\rm K}Z^{-1}_{\rm eq}}({\tr}{\bm{\tau}}+3).$
(28)

Here, $Z_{\rm eq}$ indicates $Z$ at equilibrium. From Eq. (28), we can see
that ${\tr}{\bm{\tau}}$ is tightly related to the (squared) length of
dumbbells. Since we do not address the spatial gradient in rheological
calculations, ${\rm D}(\cdots)/{\rm D}t$ simply reduces to ${\rm
d}(\cdots)/{\rm d}t$. Using Eqs. (22), (27), and (28), the constitutive
equations for the FENE-P dumbbell model under shear flow can be expressed as

$\displaystyle\dot{\tau}_{xx}$ $\displaystyle=-\left\\{1+\frac{1}{3(n_{\rm
K}-1)}\right\\}\tau_{xx}-\frac{1}{3(n_{\rm K}-1)}(\tau_{yy}+\tau_{zz})$ (29)
$\displaystyle\quad-\frac{1}{9n_{\rm K}(n_{\rm
K}-1)}(\tr\bm{\tau})^{2}-\frac{1}{3n_{\rm K}}\left(2+\frac{1}{n_{\rm
K}-1}\right)\tr\bm{\tau}\tau_{xx}$
$\displaystyle\quad+2\left\\{1+\frac{1}{3(n_{\rm
K}-1)}\right\\}\tau_{xy}\kappa_{xy}-\frac{1}{9n_{\rm K}(n_{\rm
K}-1)}(\tr\bm{\tau})^{2}\tau_{xx}$ $\displaystyle\quad+\frac{2}{3(n_{\rm
K}-1)}\tau_{xx}\tau_{xy}\kappa_{xy},$ $\displaystyle\dot{\tau}_{yy/zz}$
$\displaystyle=-\left\\{1+\frac{1}{3(n_{\rm
K}-1)}\right\\}\tau_{yy/zz}-\frac{1}{3(n_{\rm K}-1)}(\tau_{xx}+\tau_{zz/yy})$
(30) $\displaystyle\quad-\frac{1}{9n_{\rm K}(n_{\rm
K}-1)}(\tr\bm{\tau})^{2}-\frac{1}{3n_{\rm K}}\left(2+\frac{1}{n_{\rm
K}-1}\right)\tr\bm{\tau}\tau_{yy/zz}$ $\displaystyle\quad+\frac{2}{3(n_{\rm
K}-1)}\tau_{xy}\kappa_{xy}-\frac{1}{9n_{\rm K}(n_{\rm
K}-1)}(\tr\bm{\tau})^{2}\tau_{yy/zz}$ $\displaystyle\quad+\frac{2}{3(n_{\rm
K}-1)}\tau_{yy/zz}\tau_{xy}\kappa_{xy},$ $\displaystyle\dot{\tau}_{xy}$
$\displaystyle=-\tau_{xy}+\kappa_{xy}+\tau_{yy}\kappa_{xy}-\frac{1}{3n_{\rm
K}}\left(2+\frac{1}{n_{\rm K}-1}\right)\tr\bm{\tau}\tau_{xy}$ (31)
$\displaystyle\quad-\frac{1}{9n_{\rm K}(n_{\rm
K}-1)}(\tr\bm{\tau})^{2}\tau_{xy}+\frac{2}{3(n_{\rm
K}-1)}\tau_{xy}^{2}\kappa_{xy}.$

For the derivation, please refer to Sec. S5 in the supporting information.
From Eqs. (29)–(31), the constitutive equation for the FENE-P model can be
expressed by a polynomial of up to third degree in $\bm{\tau}$ and
$\bm{\kappa}$. Here, we note that Eqs. (29)–(31) become equivalent to the UCM
model shown in Eqs. (S6)–(S8) in the supporting information in the limit of
$n_{\rm K}\rightarrow\infty$.

To generate noise-free training data, we solve Eqs. (23)–(26) with $n_{\rm
K}=10$ and $\Delta t=1\times 10^{-4}$ for $0\leq t\leq 100$. n We apply the
oscillatory shear flow with $\gamma_{0}=2$ and various $\omega$ values
($\omega\in\\{0.1,0.2,\ldots,1\\}$). From the computed stress data, we collect
data at the interval of $\Delta t_{\rm train}=1\times 10^{-2}$.

###  3.4 FENE dumbbell model

We finally address the FENE dumbbell model, where the spring strength is
represented by Eq. (17). Since the FENE dumbbell model does not use any
simplification for the spring strength (e.g., Peterlin approximation shown in
Eq. (18)), its analytical constitutive equation has not been obtained. We
apply Rheo-SINDy to this case to see if an approximate constitutive equation
can be obtained. The obtained equations are validated by comparing the data
obtained by numerically solving them with the data obtained by BD simulations.

The training data are generated by the BD simulations using Eqs. (12)–(15) and
(17) with $n_{\rm K}=10$, $N_{\rm p}=10^{4}$, and $\Delta t=1\times 10^{-4}$
for $0\leq t\leq 100$. We apply the oscillatory shear flows with the same
parameters as those in the FENE-P dumbbell model. The BD simulation results
with five different random seeds are averaged for each condition. Since we do
not use any approximation for the spring strength, the values of $h(t)$ differ
for each individual dumbbell. From the computed stress data, we collected data
at the interval of $\Delta t_{\rm train}=1\times 10^{-2}$.

##  4 Results and Discussions

In this section, we present the results of the case studies for the dumbbell
models. From the case studies on phenomenological constitutive equations shown
in Sec. S3 in the supporting information, we have made the following two
findings: (i) taking shear rheological data by an oscillatory shear test is
more appropriate than by a simple (constant) shear test, and (ii) among the
five optimization methods shown in Sec. 2, the STRidge or a-Lasso is superior
to the other three methods. Thus, in what follows, we generate the training
data using the oscillatory shear test, and employ the STRidge and a-Lasso as
optimization methods.

###  4.1 Hookean Dumbbell Model

Figure 3: (a) The number of total terms obtained by the STRidge (black) and
a-Lasso (red) for the training data generated by the Hookean dumbbell model,
and (b) the obtained constitutive equations. The exact equations are given in
Eqs. (19)–(21). The horizontal line in (a) indicates the correct number of
terms. In (a), circle, triangle, and square symbols indicate the total numbers
of terms obtained by the data for $N_{\rm p}=10^{3}$, $10^{4}$, and $10^{5}$,
respectively.

We first explain the results for the Hookean dumbbell model. We here used the
polynomial library that includes up to 2nd order terms of $\tau_{xx}$,
$\tau_{yy}$, $\tau_{xy}$, and $\kappa_{xy}$. Thus, the total number of
candidate terms is $N_{\bm{\Theta}}=15$ for each component.

Figure 3 shows the Rheo-SINDy results for the Hookean dumbbell model with the
different numbers of dumbbells. We note that the standard deviation of
$\bm{\tau}$ in the training data decreases proportionally with $N_{\rm
p}^{-1/2}$. From Fig. 3(a), as the value of $N_{\rm p}$ increases, sparser
solutions are obtained especially for Rheo-SINDy with the STRidge. Unlike the
case of the UCM model (cf. Fig. S1 in the supporting information), which can
be considered as the “noise-free” case of the Hookean dumbbell model, the
STRidge provides the correct number of terms only within a narrow range of
$\alpha$ values. Nevertheless, if we choose the appropriate $\alpha$ value,
the (nearly) correct constitutive equations can be found by the STRidge, as
shown in the upper part of Fig. 3(b). We note that the terms containing
$\tau_{yy}$ appear in the time evolution equation for $\tau_{xy}$ obtained by
the STRidge. Although these terms do not affect the predictions because
${\tau}_{yy}=0$, these terms do not appear in the correct equation. We
speculate that the appearance of these terms is due to the correlation effects
of the noise in $R_{x}$ and $R_{y}$ on the stress (cf. Eq. (15)). When
comparing the STRidge and a-Lasso, it is evident that the a-Lasso provides
stable and sparse solutions across a broader range of $\alpha$ values,
regardless of the $N_{\rm p}$ value. Furthermore, we confirm that the correct
equations can be obtained using the a-Lasso, as shown in the lower part of
Fig. 3(b). This partially suggests the effectiveness of the a-Lasso in
discovering essential terms from noisy data.

###  4.2 FENE-P dumbbell Model

We next examine whether Rheo-SINDy can find more complex differential
equations (i.e., the FENE-P dumbbell model) than the UCM model and the
Giesekus model. For the Rheo-SINDy regressions of the differential equations
for the conformation tensor $\bm{C}$ of the FENE-P dumbbell model explained in
Sec. 3.3, we prepare the following custom library:

${\bm{\Theta}}=\left[\begin{array}[]{cccc}{\bm{1}}&{\bm{\Omega}}(t_{1})&{\bm{\Omega}}^{2}(t_{1})&f_{\rm
FENE}(t_{1}){\bm{\Omega}}(t_{1})\\\
{\bm{1}}&{\bm{\Omega}}(t_{2})&{\bm{\Omega}}^{2}(t_{2})&f_{\rm
FENE}(t_{2}){\bm{\Omega}}(t_{2})\\\ \vdots&\vdots&\vdots&\vdots\\\
{\bm{1}}&{\bm{\Omega}}(t_{n})&{\bm{\Omega}}^{2}(t_{n})&f_{\rm
FENE}(t_{n}){\bm{\Omega}}(t_{n})\end{array}\right],$ (32)

where ${\bm{\Omega}}$ includes non-zero components of $\bm{C}$ under shear
flow ($C_{xx}$, $C_{yy}$, $C_{zz}$, and $C_{xy}$) and $\kappa_{xy}$, and
${\bm{\Omega}}^{2}$ is the vector composed of all the multiplied combinations
of the ${\bm{\Omega}}$ components. The total number of library functions is
thus $N_{\bm{\Theta}}=26$.

Figure 4 indicates (a) the total number of predicted terms and (b) the error
rate as a function of the hyperparameter $\alpha$ for the STRidge and the
a-Lasso.

Figure 4: (a) The total number of terms and (b) the error rate for the
conformation tensor $\bm{C}$ of the FENE-P dumbbell model obtained by Rheo-
SINDy with the STRidge (black squares) and a-Lasso (red reverse triangles).
The horizontal line in (a) indicates the correct number of terms. The training
data were generated by Eqs. (24)–(26) with $n_{\rm K}=10$.

The error rate is defined as the sum of the mean squared errors (MSEs) of
$\dot{\bm{t}}_{\mu\nu}-{\bm{\Theta}}\hat{\bm{\xi}}_{\mu\nu}$. The MSEs were
scaled so that the maximum value of each method was 1. Similar to the results
for the phenomenological constitutive equations shown in Sec. S3 in the
supporting information, the a-Lasso provides sparser solutions than the
STRidge, and the STRidge gives lower error rates than the a-Lasso. Figure 5
presents the differential equations obtained by the STRidge and a-Lasso for
two $\alpha$ values that yield the nearly correct number of terms with a small
error rate.

Figure 5: The differential equations for the conformation tensor $\bm{C}$ of
the FENE-P dumbbell model found by Rheo-SINDy with the STRidge and the
a-Lasso. The exact equations are given in Eqs. (24)–(26).

From the lower part of Fig. 5, while the a-Lasso can provide sparser
solutions, it does not guarantee that these are correct (cf. Eqs. (24)–(26)).
Specifically, in all cases for $\tau_{xx}$, $\tau_{yy}$, and $\tau_{zz}$, the
a-Lasso has failed to identify the constant term in Eqs. (24) and (25), which
is a possible source of larger errors compared to the STRidge. In the case of
the STRidge, we confirmed that by choosing the appropriate $\alpha$
($\alpha=1\times 10^{-1}$), nearly correct differential equations can be
obtained, as shown in the upper part of Fig. 5. Since the $yy$-component and
$zz$-component of the stress are equivalent, the exact equations can be
recovered by setting $C_{yy}=C_{zz}$. Thus, we found that the correct
differential equations for the FENE-P dumbbell model can be obtained if we can
prepare the proper library functions and choose the appropriate value of the
hyperparameter. Figure 6 shows the test simulation results using the
identified differential equations for $\bm{C}$ in Fig. 5 and the dimensionless
form of Eq. (23).

Figure 6: The test simulation results using the equations obtained by (a) the
STRidge and (b) a-Lasso. Here, the test simulations were conducted with
$\gamma_{0}=4$ and $\omega=1$. The black, blue, and red lines show
$\tau_{xx}$, $\tau_{yy}$, and $\tau_{xy}$. The bold, thin dotted, and thin
solid lines indicate the exact solutions, predictions with smaller $\alpha$
values ($\alpha=1\times 10^{-3}$ for the STRidge and $\alpha=1\times 10^{-7}$
for the a-Lasso), and predictions with larger $\alpha$ values ($\alpha=1\times
10^{-1}$ for the STRidge and $\alpha=1\times 10^{-4}$ for the a-Lasso).

Here, the oscillatory shear flow with $\gamma_{0}=4$ and $\omega=1$, which is
outside of the training data, was considered. From Fig. 6, the equations
obtained by the STRidge can reproduce the exact solutions even when the
equations are not exactly correct ($\alpha=1\times 10^{-3}$). In contrast, the
test simulations with the differential equations obtained by the a-Lasso show
the deviations from the test data, especially for $\tau_{xx}$. These results
emphasize the need to choose an appropriate optimization method to obtain
reasonable solutions.

We then examine whether the stress expression of the constitutive equation for
the FENE-P dumbbell model (cf. Eqs. (29)–(31)) can be found by Rheo-SINDy. For
such a purpose, we prepared the following custom library:

${\bm{\Theta}}=\left[\begin{array}[]{cccc}1&\\{\tr{\bm{\tau}}(t_{1})\\}^{p}{\bm{T}}_{\rm
s}(t_{1})&\\{\tr{\bm{\tau}}(t_{1})\\}^{2}&\\{{\bm{T}}_{\rm
s}(t_{1})\\}^{p}\kappa_{xy}(t_{1})\\\
1&\\{\tr{\bm{\tau}}(t_{2})\\}^{p}{\bm{T}}_{\rm
s}(t_{2})&\\{\tr{\bm{\tau}}(t_{2})\\}^{2}&\\{{\bm{T}}_{\rm
s}(t_{2})\\}^{p}\kappa_{xy}(t_{2})\\\ \vdots&\vdots&\vdots&\vdots\\\
1&\\{\tr{\bm{\tau}}(t_{n})\\}^{p}{\bm{T}}_{\rm
s}(t_{n})&\\{\tr{\bm{\tau}}(t_{n})\\}^{2}&\\{{\bm{T}}_{\rm
s}(t_{n})\\}^{p}\kappa_{xy}(t_{n})\\\ \end{array}\right],$ (33)

where ${\bm{T}}_{\rm s}$ includes
$\\{\tau_{xx},\tau_{yy},\tau_{zz},\tau_{xy}\\}$ and $p$ ($=0,1,2$) is the
polynomial order. Thus, the total number of library functions is
$N_{\bm{\Theta}}=29$. We prepared the library that includes at least the terms
present in Eqs. (29)–(31). Furthermore, we excluded terms that could
potentially become large, such as higher-order terms involving $\kappa_{xy}$.
When such terms are included in the solutions, the differential equations may
be unstable, and in worse cases, they may also diverge.

Figure 7: (a) The total number of terms and (b) the error rate for the
constitutive equation of the FENE-P dumbbell model obtained by the STRidge
(black squares) and the a-Lasso (red reverse triangles). The horizontal short-
dashed line in (a) indicates that the number of terms is zero. Figure 8: The
constitutive equations for the FENE-P dumbbell model obtained by the STRidge
and a-Lasso. The exact equations are given in Eqs. (29)–(31). Figure 9: Test
simulation results under the oscillatory shear flow with $\gamma_{0}=4$ and
$\omega=1$ for the constitutive equations of the FENE-P dumbbell model
obtained by (a) the STRidge and (b) a-Lasso. The black, blue, and red lines
represent the $xx$-, $yy$-, and $xy$-components of the stress tensor,
respectively. The bold lines show the exact solutions. The thin solid and
short-dashed lines indicate the results with smaller $\alpha$ values
($\alpha=1\times 10^{-2}$ for the STRidge and $\alpha=3\times 10^{-8}$ for the
a-Lasso) and with larger $\alpha$ values ($\alpha=1$ for the STRidge and
$\alpha=1\times 10^{-4}$ for the a-Lasso).

Figure 7 shows (a) the total number of terms and (b) the error rate obtained
by Rheo-SINDy with the STRidge and a-Lasso. Similar to what we noted
previously, the a-Lasso can yield sparser solutions than the STRidge. Based on
the number of terms shown in Fig. 7(a) and the error rates shown in Fig. 7(b),
we chose several $\alpha$ values with a small number of terms and a low error
rate. Figure 8 presents the equations obtained using the chosen $\alpha$. From
Fig. 8, the equations predicted by the STRidge with $\alpha=1$ and the a-Lasso
with $\alpha=1\times 10^{-4}$ are almost the same; conversely, the solutions
for small $\alpha$ values significantly differ between the two methods. For
the STRidge with $\alpha=1\times 10^{-2}$, the identified equations are close
to the correct equations (cf. Eqs. (29)–(31)). Furthermore, the coefficient
values for the correctly obtained terms are close to the correct values. For
the a-Lasso with $\alpha=3\times 10^{-8}$, several coefficients for the
correctly obtained terms, such as $\tau_{xx}$, $\tau_{xy}$, and
$\tau_{xx}\tau_{xy}\kappa_{xy}$ in the equation for $\dot{\tau}_{xx}$, are
close to the exact values, but for other several terms, such as
$\tr\bm{\tau}\tau_{xx}$ in the equation for $\dot{\tau}_{xx}$, the correct
coefficient values are not obtained. Nevertheless, from Fig. 9, which shows
the test simulation results, the equations obtained by the STRidge with
$\alpha=1\times 10^{-2}$ and the a-Lasso with $\alpha=3\times 10^{-8}$ can
well reproduce the exact solutions including the small oscillation of
$\tau_{yy}$. Although the equations obtained by the STRidge and a-Lasso
demonstrate the similar performance in the test simulations shown in Fig. 9,
the difference in predictions is quantified by their MSEs shown in Table 1.
When $\alpha$ is small, the error in $\tau_{xx}$ is of the same order for both
methods, but for predictions of $\tau_{yy}$ and $\tau_{xy}$, the STRidge
outperforms the a-Lasso. The STRidge, however, provides a sparse solution
within a narrow range of $\alpha$ values, requiring careful selection of
$\alpha$.

Table 1: The mean squared error (MSE) between predicted and exact solutions
for the FENE-P dumbbell model.

  * method | $\alpha$ | MSE ($\tau_{xx}$) | MSE ($\tau_{yy}$) | MSE ($\tau_{xy}$)
---|---|---|---|---
STRidge | $1\times 10^{-2}$ | $1.1\times 10^{-1}$ | $5.0\times 10^{-5}$ | $5.1\times 10^{-3}$
STRidge | $1$ | $2.3$ | $8.2\times 10^{-3}$ | $8.9\times 10^{-2}$
a-Lasso | $3\times 10^{-8}$ | $3.4\times 10^{-1}$ | $3.0\times 10^{-3}$ | $3.8\times 10^{-2}$
a-Lasso | $1\times 10^{-4}$ | $2.3$ | $8.2\times 10^{-3}$ | $9.0\times 10^{-2}$

###  4.3 FENE dumbbell Model

Figure 10: (a) The total number of terms and (b) the error rate for the FENE
dumbbell model predicted by the STRidge (black squares) and the a-Lasso (red
reverse triangles). The horizontal short-dashed line in (a) indicates that the
number of terms is zero. Figure 11: The predicted constitutive equations for
the FENE dumbbell model (left) and the test simulation results (right). Here,
the a-Lasso was utilized to obtain the approximate constitutive equations. For
test simulations, we solved the constitutive equations under the oscillatory
shear flows with $\gamma_{0}=3$ and $\omega=1$ (right upper panel) and
$\gamma_{0}=4$ and $\omega=1$ (right lower panel). The bold lines show the
exact solutions, and the thin solid and short-dashed lines show the results
with the smaller $\alpha$ value ($\alpha=1\times 10^{-6}$) and the larger
$\alpha$ value ($\alpha=3\times 10^{-4}$).

Finally, we address the FENE dumbbell model. As explained in Sec. 3.4, the
FENE dumbbell model does not have an analytical expression of the constitutive
equation. Thus, we here develop an approximate constitutive equation using
Rheo-SINDy.

To obtain dynamical equations by Rheo-SINDy, one first needs to design an
appropriate library $\bm{\Theta}$. To prepare $\bm{\Theta}$ for the FENE
dumbbell model, we utilize the physical insights obtained from the analytical
expression of the FENE-P dumbbell model. We here assume the constitutive
equation of the FENE-P dumbbell model is similar to that of the FENE dumbbell
model. Since the FENE-P dumbbell model is a simplified version of the FENE
dumbbell model, we believe that this is a reasonable assumption. Here, we note
that the stress expression shown in Eq. (23) is no longer applicable to the
FENE dumbbell model since the values of $h(t)$ differ for each individual
dumbbell. Thus, it is invalid to obtain stress through the conformation tensor
$\bm{C}$. Based on the above considerations, we decided to use the custom
library presented in Eq. (33), which was also used in the constitutive
equation of the FENE-P dumbbell model.

Figure 10 compares (a) the total number of terms and (b) the error rate
predicted by the STRidge and a-Lasso. Similar to the previous discussions, we
can obtain sparse solutions over a wide range of $\alpha$ values with the
a-Lasso, whereas the STRidge gives sparse solutions only within a limited
range of $\alpha$. The left table in Fig. 11 shows the equations obtained by
the a-Lasso with two $\alpha$ values chosen from the viewpoints of the
sparsity and error rate in the same way as Fig. 5. We note that the
predictions obtained by the STRidge are inferior to those obtained by the
a-Lasso shown in Fig. 11, which is discussed in Sec. S6 in the supporting
information. From the left table in Fig. 11, if $\alpha$ is appropriately
chosen, the a-Lasso can give sparse equations with coefficients of reasonable
(not excessively large) magnitudes. Comparing the equations for the FENE-P
model obtained by the a-Lasso with $\alpha=3\times 10^{-8}$ (Fig. 8) and those
for the FENE model obtained by the a-Lasso with $\alpha=1\times 10^{-6}$ (Fig.
11), the appearing terms are almost identical, which demonstrates the
similarity between these models. The difference in the coefficients thus
represents the difference between these models. The right panels in Fig. 11
show the test simulation results obtained by the equations shown in the left
table. We found that the equations obtained with $\alpha=1\times 10^{-6}$ can
reproduce well the BD simulation results outside the range of the training
data within the investigated parameters, including the oscillatory behavior of
$\tau_{yy}$. (With the large $\alpha$ ($\alpha=3\times 10^{-4}$), the
identified equation for $\tau_{yy}$ becomes ${\dot{\tau}}_{yy}=0$, which fails
to reproduce the oscillatory behavior of $\tau_{yy}$.) This success suggests
that Rheo-SINDy with the a-Lasso is effective for discovering unknown
constitutive equations. Nevertheless, we note that the equations presented in
Fig. 11 may fail to predict test data significantly outside the range of the
training data. Reproducing such highly nonlinear data would require the
nonlinear terms dropped in Fig. 11. In this sense, the constitutive equations
for the FENE dumbbell model obtained here are appropriately referred to as the
approximate constitutive equations.

Thanks to the equations obtained using Rheo-SINDy, it is possible to provide a
physical interpretation with the assistance of rheological knowledge. For
example, from the comparison of the equations obtained for the FENE-P dumbbell
model (cf. Fig. 8) and those for the FENE dumbbell model (cf. Fig. 11), the
equations for larger $\alpha$ value ($\alpha=1\times 10^{-4}$ for the FENE-P
dumbbell model and $\alpha=3\times 10^{-4}$ for the FENE dumbbell model) are
similar except for the coefficient values. Furthermore, the terms in these
equations are the same as those for the UCM model (and thus the Hookean
dumbbell model). This indicates that all of these models share the same origin
based on the dumbbell model. The linear term of stress in the constitutive
equation represents the relaxation of stress (see Eq. (S5) in the supporting
information). Since the relaxation time at equilibrium ($\lambda=\zeta/4h_{\rm
eq}$) is taken as the unit time in this study, the coefficient of this term
should be $-1$ at equilibrium (and thus for the UCM model, see Eqs.
(S6)–(S8)). From Figs. 8 and 11, the coefficient of the linear term of stress
is smaller than $-1$, which indicates $\lambda_{\rm sf}<\lambda_{\rm eq}$ with
the subscript “sf” and “eq” standing for “shear flow” and “equiliblium”,
respectively. This indicates that under shear flow, the values of spring
strength for the FENE-P and FENE dumbbell models become larger than $h_{\rm
eq}$, which implies the appearance of the FENE effects under flow. From this
discussion, it is evident that Rheo-SINDy can provide physically interpretable
constitutive equations.

##  5 Concluding Remarks

We tested that the sparse identification for nonlinear dynamics (SINDy)
modified for nonlinear rheological data, which we call Rheo-SINDy, is
effective in finding constitutive equations of complex fluids. We found that
Rheo-SINDy can successfully identify correct equations from training data
generated from known constitutive equations, as well as provide approximate
constitutive equations (or reduced order models) from training data generated
by mesoscopic models when constitutive equations are analytically unknown.

Rheo-SINDy for two phenomenological constitutive equations (i.e., the upper
convected Maxwell model and Giesekus model) revealed the following two things.
First, compared to constant shear tests, oscillatory shear tests are
appropriate for generating training data. Second, the sequentially thresholded
Ridge regression (STRidge) and adaptive Lasso (a-Lasso) are effective in
finding appropriate constitutive equations. We then examined the commonly used
mesoscopic model, namely the dumbbell model with three different
representations of spring strength: the Hookean, FENE-P, and FENE springs.
Although the Hookean and FENE-P dumbbell models have analytical constitutive
equations, for the FENE dumbbell model, there is no analytical expression of
the constitutive equation. We confirmed through the Hookean dumbbell model
that even in the presence of noise, the a-Lasso provides the correct solution
over a wide range of the hyperparameter $\alpha$. Rheo-SINDy was also
effective in discovering the complex constitutive equations of the FENE-P
dumbbell model. This case study revealed that the identification of complex
equations requires the preparation of an appropriate custom library based on
prior physical knowledge. Using physical insights obtained from the Hookean
and FENE-P dumbbell models, we attempted to find approximate constitutive
equations for the FENE dumbbell model. We found that the a-Lasso can
successfully give the approximate constitutive equations, which can be used in
predictions beyond the range of the training data.

From our investigation, Rheo-SINDy with the STRidge or a-Lasso is effective
for discovering constitutive equations from nonlinear rheological data. We
found that the STRidge is generally superior in terms of retaining correct
terms, while the a-Lasso is more robust to the selection of $\alpha$ than the
STRidge. To obtain correct constitutive equations, in addition to selecting
the appropriate optimization method, we are required to design an appropriate
library by using physical insights, namely domain knowledge. Designing such a
proper library necessitates not only including necessary terms but also
excluding unnecessary terms.

This research is expected to have an impact on fields such as rheology and
fluid mechanics. From a rheological perspective, for several systems such as
entangled polymers [42, 43] and wormlike micellar solutions [44, 45],
sophisticated mesoscopic models suitable for numerical simulations under flow
have been proposed. These mesoscopic models can generate reasonable training
data not only under shear flow but also under extensional flow. Finding new
approximate models from the data obtained by these mesoscopic simulations
would be an interesting research subject. Furthermore, it would be desirable
to conduct Rheo-SINDy for experimental data obtained by Large Amplitude
Oscillatory Shear (LAOS) experiments [46]. Since the LAOS measurements do not
provide all the major stress components under shear flow, exploring methods
for discovering the constitutive equations from experimental data would be a
future challenge. When approximate constitutive models are identified, those
models can be employed for predictions of complex flows, which would deepen
our understanding of complex fluids. We will continue our research in these
directions.

## Data availability statement

The data that support the findings of this study are available upon reasonable
request from the authors.

We thank Prof. Yoshinobu Kawahara for his valuable comments on data-driven
methods. We also thank Dr. John J. Molina for carefully reading the manuscript
and providing insightful comments. TS was partially supported by JST PRESTO
Grant Number JPMJPR22O3. SM was financially supported by the Kyoto University
Science and Technology Innovation Creation Fellowship (FSMAT), Grant Number
JPMJFS2123.

## References

## References

  * [1] Brunton S L and Kutz J M 2022 Data-Driven Science and Engineering 2nd ed (Cambridge University Press)
  * [2] Brunton S L, Proctor J L and Kutz J N 2016 Proc. Natl. Acad. Sci. USA 113 3932–3937
  * [3] de Silva B M, Higdon D M, Brunton S L and Kutz J N 2020 Front. Artif. Intell. 3 25
  * [4] Kaheman K, Brunton S L and Nathan Kutz J 2022 Mach. Learn.: Sci. Technol. 3 015031
  * [5] Fasel U, Kutz J N, Brunton B W and Brunton S L 2022 Proc. Math. Phys. Eng. Sci. 478 20210904
  * [6] Schmidt M and Lipson H 2009 Science 324 81–85
  * [7] Bongard J and Lipson H 2007 Proc. Natl. Acad. Sci. USA 104 9943–9948
  * [8] Reinbold P A K, Kageorge L M, Schatz M F and Grigoriev R O 2021 Nat. Commun. 12 3219
  * [9] Udrescu S M and Tegmark M 2020 Sci. Adv. 6 eaay2631
  * [10] Cranmer M, Sanchez Gonzalez A, Battaglia P, Xu R, Cranmer K, Spergel D and Ho S 2020 Discovering symbolic models from deep learning with inductive biases NeurIPS vol 33 ed Larochelle H, Ranzato M, Hadsell R, Balcan M and Lin H pp 17429–17442
  * [11] Lemos P, Jeffrey N, Cranmer M, Ho S and Battaglia P 2023 Mach. Learn.: Sci. Technol. 4 045002
  * [12] Karniadakis G E, Kevrekidis I G, Lu L, Perdikaris P, Wang S and Yang L 2021 Nat. Rev. Phys. 3 422–440
  * [13] Raissi M, Perdikaris P and Karniadakis G E 2019 Journal of computational physics 378 686–707 ISSN 0021-9991
  * [14] Jia X, Willard J, Karpatne A, Read J S, Zwart J A, Steinbach M and Kumar V 2021 ACM/IMS Trans. Data Sci. 2 1–26
  * [15] Rosofsky S G, Al Majed H and Huerta E A 2023 Machine Learning: Science and Technology 4 025022
  * [16] Larson R G 1988 Constitutive Equations for Polymer Melts and Solutions (Butterworths Series in Chemical Engineering)
  * [17] Rouse P E 1953 J. Chem. Phys. 21 1272–1280
  * [18] Doi M and Edwards S F 1986 The Theory of Polymer Dynamics (Oxford University Press)
  * [19] Masubuchi Y, Takimoto J I, Koyama K, Ianniruberto G, Marrucci G and Greco F 2001 J. Chem. Phys. 115 4387–4394
  * [20] Doi M and Takimoto J 2003 Phil. Trans. R. Soc. A. 361 641–652
  * [21] Likhtman A E 2005 Macromolecules 38 6128–6139
  * [22] Jamali S 2023 Rheology Bulletin 92(1) 20–24
  * [23] Miyamoto S 2024 Nihon Reoroji Gakkaishi (J. Soc. Rheol. Jpn.) 52 15–19
  * [24] Fang L, Ge P, Zhang L, E W and Lei H 2022 J. Mach. Learn. 1 114–140
  * [25] Mahmoudabadbozchelou M, Kamani K M, Rogers S A and Jamali S 2022 Proc. Natl. Acad. Sci. USA 119 e2202234119
  * [26] Jin H, Yoon S, Park F C and Ahn K H 2023 Rheol. Acta 62 569–586
  * [27] Mahmoudabadbozchelou M and Jamali S 2021 Sci. Rep. 11 12015
  * [28] Mahmoudabadbozchelou M, Karniadakis G E and Jamali S 2022 Soft Matter 18(1) 172–185
  * [29] Saadat M, Mahmoudabadbozchelou M and Jamali S 2022 Rheol. Acta 61 721–732
  * [30] Mahmoudabadbozchelou M, Caggioni M, Shahsavari S, Hartt W H, Em Karniadakis G and Jamali S 2021 J. Rheol. 65 179–198
  * [31] Lennon K R, McKinley G H and Swan J W 2023 Proc. Natl. Acad. Sci. USA 120 e2304669120
  * [32] Zhao L, Li Z, Caswell B, Ouyang J and Karniadakis G E 2018 J. Comput. Phys. 363 116–127
  * [33] Zhao L, Li Z, Wang Z, Caswell B, Ouyang J and Karniadakis G E 2021 J. Comput. Phys. 427 110069
  * [34] Seryo N, Sato T, Molina J J and Taniguchi T 2020 Phys. Rev. Res. 2 033107
  * [35] Seryo N, Molina J J and Taniguchi T 2021 Nihon Reoroji Gakkaishi (J. Soc. Rheol. Jpn.) 49 97–113
  * [36] Miyamoto S, Molina J J and Taniguchi T 2023 Phys. Fluids 35 063113
  * [37] Fukami K, Murata T, Zhang K and Fukagata K 2021 J. Fluid Mech. 926 A10
  * [38] Mahmoudabadbozchelou M, Kamani K M, Rogers S A and Jamali S 2024 Proc. Natl. Acad. Sci. USA 121 e2313658121
  * [39] Bird R B, Armstrong R C and Hassager O 1987 Dynamics of Polymeric Liquids 2nd ed vol 2 (Oxford University Press)
  * [40] Graham M D 2014 Phys. Fluids 26 101301
  * [41] Mochimaru Y 1983 J. Non-Newtonian Fluid Mech. 12 135–152
  * [42] Sato T and Taniguchi T 2019 Macromolecules 52 3951–3964
  * [43] Miyamoto S, Sato T and Taniguchi T 2023 Rheol. Acta 62 57–70
  * [44] Sato T, Moghadam S, Tan G and Larson R G 2020 J. Rheol. 64 1045–1061
  * [45] Sato T and Larson R G 2022 J. Rheol. 66 639–656
  * [46] Hyun K, Wilhelm M, Klein C O, Cho K S, Nam J G, Ahn K H, Lee S J, Ewoldt R H and McKinley G H 2011 Prog. Polym. Sci. 36 1697–1753