# Multi-element flow-driven spectral chaos (ME-FSC) method for uncertainty
quantification of dynamical systems

Hugo Esquivel<EMAIL_ADDRESS><EMAIL_ADDRESS>Arun Prakash
<EMAIL_ADDRESS>Guang Lin<EMAIL_ADDRESS>Lyles School of Civil
Engineering, Purdue University, West Lafayette, Indiana, USA  Department of
Civil and Environmental Engineering, Universidad de la Costa, Barranquilla,
Colombia  Department of Mathematics, School of Mechanical Engineering, Purdue
University, West Lafayette, Indiana, USA

###### Abstract

The flow-driven spectral chaos (FSC) is a recently developed method for
tracking and quantifying uncertainties in the long-time response of stochastic
dynamical systems using the spectral approach. The method uses a novel concept
called _enriched stochastic flow maps_ as a means to construct an evolving
finite-dimensional random function space that is both accurate and
computationally efficient in time. In this paper, we present a multi-element
version of the FSC method (the ME-FSC method for short) to tackle (mainly)
those dynamical systems that are inherently discontinuous over the probability
space. In ME-FSC, the random domain is partitioned into several elements, and
then the problem is solved separately on each random element using the FSC
method. Subsequently, results are aggregated to compute the probability
moments of interest using the law of total probability. To demonstrate the
effectiveness of the ME-FSC method in dealing with discontinuities and long-
time integration of stochastic dynamical systems, four representative
numerical examples are presented in this paper, including the Van-der-Pol
oscillator problem and the Kraichnan-Orszag three-mode problem. Results show
that the ME-FSC method is capable of solving problems that have strong
nonlinear dependencies over the probability space, both reliably and at low
computational cost.

###### keywords:

stochastic discontinuities; stochastic dynamical systems; uncertainty
quantification; long-time integration; stochastic flow map; multi-element
flow-driven spectral chaos (ME-FSC).

††journal: Journal of Computational Physics (Elsevier)

## Highlights

  * •

A multi-element flow-driven spectral chaos method.

  * •

ME-FSC handles discontinuities over the probability space adequately.

  * •

ME-FSC quantifies uncertainties in the long-time response reliably.

  * •

ME-FSC is highly accurate and computationally efficient.

##  1 Introduction

The spectral approach has gained increasing popularity in the past few decades
as a powerful tool for solving stochastic problems at low-computational cost.
This assertion is especially true for problems where the dimensionality of the
probability space is relatively low. For problems where the dimensionality of
the probability space is relatively high, this assertion may not always hold,
especially in cases where the curse of dimensionality at the random-function-
space (RFS) level cannot be alleviated noticeably. This is, for instance, the
case of gPC-based methods (e.g. [1, 2, 3, 4, 5]) where the curse-of-
dimensionality issue arises naturally whenever a suitable random function
space (aka RFS in this paper) is defined over the probability space to
represent approximately the stochastic part of the system’s response. In other
words, when the random input consists of mutually independent random
variables, the stochastic part of the solution space is customarily spanned
via a tensor product of vector spaces with each of these vector spaces
representing a single space spanned by a finite set of univariate orthogonal
polynomials. As a consequence, the higher the dimensionality of the
probability space, the higher the dimensionality of the resulting RFS. It is
worth mentioning that an RFS constructed this way can only ensure exponential
convergence to the solution in the early times of the simulation, unless the
RFS is kept updated frequently during the simulation.

A brief historical review is given next. Wan and Karniadakis [6] developed the
multi-element generalized polynomial chaos (ME-gPC) method to handle long-time
integration and discontinuities satisfactorily in problems with random input
data. In this approach, the random domain is decomposed into smaller
subdomains (aka random elements) whenever the relative error in variance
exceeds a predefined threshold value during the simulation. For computational
efficiency, each subdomain uses a relatively low-degree polynomial chaos from
the Askey scheme to ensure exponential convergence to the solution locally
[1]. The effectiveness of ME-gPC to deal with long-time integration was
studied in detail in [7] (see also [8]), and subsequently the method was
further generalized for inputs with arbitrary probability density functions
[9]. To benefit from the computational ease of sampling-based methods, Foo et
al. [10] developed the multi-element probabilistic collocation (ME-PCM)
method. This approach uses the decomposition of the random domain and a tensor
product grid in each subdomain to deal with long-time integration and
stochastic discontinuities more efficiently. In [11], the rate of convergence
was improved by introducing an analysis of variance (ANOVA) into the scheme.
Furthermore, Zheng et al. [12] developed an adaptive multi-element polynomial
chaos for discrete probability density functions, and Kawai and Oyama [13]
proposed a multi-element stochastic Galerkin method based on edge detection to
circumvent the Gibbs phenomenon in the presence of stochastic discontinuities.

Asokan and Zabaras [14] developed a variational multiscale framework to solve
elliptic problems with heterogeneous random media. The framework combines
ideas from the variational multiscale method, the multiscale finite element
method, and the generalized polynomial chaos method to derive a variationally
consistent upscaling technique to model diffusion in problems involving
heterogeneous random media. Witteveen et al. [15] formulated an adaptive
stochastic finite elements approach based on the Newton-Cotes quadrature rule.
Here, the stochastic part of the solution is approximated with piecewise
polynomial chaos by subdividing the random domain into simplex elements. Nouy
and Clément [16] combined the extended finite element method with spectral
stochastic methods to simulate heterogeneous materials with random material
interfaces. Here, to enhance the convergence rate of the scheme, an extension
of the partition of unity method was proposed in the spectral stochastic
framework.

Ganapathysubramanian and Zabaras [17] presented an application of sparse grid
collocation schemes based on the Smolyak algorithm to solve high-dimensional
stochastic convection problems more effectively. Ma and Zabaras [18] proposed
an adaptive sparse grid collocation strategy using piecewise multi-linear
hierarchical basis functions. To automatically detect stochastic
discontinuities, and thus, adaptively refine the collocation points in the
problematic subdomains, hierarchical surpluses were used as error indicators.
Bhaduri and Graham-Brady [19] improved this technique by avoiding unnecessary
function evaluations in smoother regions of the random domain by using
successive derivative estimations along all random dimensions (see also [20]).
Agarwal and Aluru [21] proposed an adaptive sparse grid collocation strategy
to investigate the performance of micro-electromechanical systems under
uncertainties.

Marzouk et al. [22] presented a reformulation of the Bayesian approach to
inverse problems to accelerate the Bayesian inference via the gPC method.
Mohan et al. [23] developed the multi-element stochastic reduced basis methods
(ME-SRBMs) for solving linear stochastic partial differential equations. In
their approach, the random domain is decomposed into several subdomains and
the solution is approximated on each subdomain using a set of basis vectors
spanning a preconditioned stochastic Krylov subspace. Sarrouy et al. [24]
proposed a piecewise polynomial chaos expansion to perform a stability
analysis of a linear brake system with uncertainties using the complex
eigenvalue analysis method and decomposition of the random domain to realize a
low degree piecewise polynomial approximation. Jakeman et al. [25] developed
the minimal multi-element stochastic method to quantify uncertainties of
discontinuous functions more efficiently than standard locally-adaptive sparse
grid methods. This approach consists of two major steps. First, a
discontinuity detector is introduced into the scheme to partition the random
domain into elements of high regularity. Then, an adaptive technique based on
the least orthogonal interpolant is employed to construct a gPC approximation
on each element.

Maître et al. [26] introduced the Wiener-Haar expansion in the context of
uncertainty propagation to investigate situations in which the solution can
vary rapidly over the probability space. By combining concepts from the
generalized polynomial chaos methods [1, 6] and stochastic Galerkin methods,
this wavelet-based technique was demonstrated to be suitable for resolving
stochastic systems with multiple stochastic discontinuities (see also [27]).
To guarantee optimal convergence with respect to the full stochastic and
spatial discretization, adaptive wavelets were used in [28] to construct
simultaneously the stochastic and spatial bases and omit the intermediate
semidiscrete approximation stage. Gittelson [29] derived an adaptive solver
for elliptic boundary value problems with random coefficients. Though this
technique uses ideas from adaptive wavelet methods, orthonormal polynomials
were used instead.

Cho et al. [30] developed an adaptive discontinuous Galerkin method for
response-excitation PDF equations using a nonconforming adaptive discontinuous
Galerkin method for the response space and a probabilistic collocation method
for the excitation space. Heuveline and Schick [3] proposed a hybrid
generalized polynomial chaos method to tackle problems featuring strong
nonlinear dependencies on the stochastic inputs (with or without stochastic
discontinuities). This was achieved by combining ideas from TD-gPC [2] and ME-
gPC [6] methods. Chen et al. [31] developed a localized polynomial chaos
expansion for high-dimensional stochastic problems to circumvent the
difficulty of having a high-dimensional probability space in the mathematical
model of the system. Li and Stinis [32] presented a unified framework for mesh
refinement in random and physical space. The benefit of the framework is that
it does not require any explicit knowledge of a reduced model to perform the
mesh refinement in both these spaces. Finally, B-splines basis functions have
also been used for uncertainty propagation to overcome the limitations of
classical approaches in the presence of discontinuities (see e.g. [33, 34]).

Lastly, the flow-driven spectral chaos (FSC) [35, 36] is a new numerical
method developed by the authors of the present paper to tackle mainly the
long-time integration issue found in the gPC method using the spectral
approach. The method uses the concept of _enriched_ stochastic flow maps to
track the evolution of a finite-dimensional RFS efficiently in time. In FSC,
the dimensionality of the random phase space is deliberately increased to
allow both the system’s state and its first few time derivatives to be pushed
forward over time. Then, the enriched state of the system is used as a germ to
construct a suitable RFS for use within the current time step of the
simulation. It is worth mentioning that for dynamical systems with an
$n$-tuple state and driven by a stochastic flow map of order $M$, the
cardinality of the random basis is bounded from above by $n+M+1$. This
boundedness from above is what makes the FSC method be curse-of-dimensionality
free at the RFS level, even when the probability space is high-dimensional.

In this paper, we present a multi-element version of the FSC method to deal
with stochastic discontinuities and long-time integration of stochastic
dynamical systems more efficiently. This new technique is called the _multi-
element flow-driven spectral chaos_ (ME-FSC) method. In ME-FSC, the random
domain is partitioned into several elements, and then the problem is solved
separately on each random element using the FSC method. Then, the results are
aggregated to compute the mean and variance of the response with the law of
total probability. This approach is similar to the multi-element gPC (ME-gPC)
method [6], with the only difference being that the gPC method is not employed
on each random element to fully annihilate the curse of dimensionality at the
RFS level. The benefit of using ME-FSC is threefold. First, the simulation can
be run simultaneously on machines with multiple CPU cores (or if needed on
separate machines) to reduce the excessive computational burden associated
with the simulation. Second, suppose the random input is discontinuous over
the probability space. In that case, the random domain can be partitioned into
several elements so that the discontinuity only appears on regions of measure
zero. Third, if an adaptive criterion is introduced within the ME-FSC scheme
(such that the number of elements gets smaller on-the-fly whenever a threshold
value is exceeded), the errors can be kept to a minimum during the simulation.

This paper has been structured in the following way. The setting and notation
used in the manuscript are formally introduced in Section 2, and then the
problem statement we are interested in is presented in Section 3. A quick
overview of the concept associated with ‘enriched stochastic flow maps’ is
provided in Section 4 to allow the construction of an optimal, low-dimensional
RFS during the simulation. Subsequently, in Section 5 an overview of the FSC
method is imparted to highlight the key idea behind the FSC method. In Section
6 the ME-FSC method is presented in detail along with an outline of the
numerical scheme used in this work. Finally, in Section 7, four numerical
examples are explored to test the performance of the ME-FSC method (against
well-established UQ techniques) and document the findings.

##  2 Setting and notation

###  2.1 Basic spaces needed in the ME-FSC method

Definitions 1 to 5 presented in two previous works by the authors [35, 36] are
also considered herein. The _temporal space_ is a topological space defined by
$\mathfrak{T}=[0,T]$, where $T\in\mathbb{R}^{+}$ symbolizes the duration of
the simulation. The _probability space_ is defined by
$(\Omega,\boldsymbol{\Omega},\lambda)$, where $\Omega$ is the sample space,
$\boldsymbol{\Omega}\subset 2^{\Omega}$ is the $\sigma$-algebra on $\Omega$,
and $\lambda:\boldsymbol{\Omega}\to[0,1]$ is the probability measure on
$\boldsymbol{\Omega}$. Because this probability space can be abstract, a
$d$-valued random variable
$\xi=(\xi^{1},\ldots,\xi^{d}):\Omega\to\mathbb{R}^{d}$ given by
$\xi=\xi(\omega)$ is defined for computational purposes. Thus, the _random
space_ is a measure space defined by $(\Xi,\boldsymbol{\Xi},\mu)$, where
$\Xi=\xi(\Omega)\subset\mathbb{R}^{d}$ is a set representing the random domain
of the system, $\boldsymbol{\Xi}=\mathcal{B}_{\mathbb{R}^{d}}\cap\Xi$ is the
$\sigma$-algebra on $\Xi$, $\mu:\boldsymbol{\Xi}\to[0,1]$ is the probability
measure on $\boldsymbol{\Xi}$ given by
$\mu=\xi_{*}(\lambda):=\lambda\circ\xi^{-1}$ (i.e. the push-forward of
$\lambda$ by $\xi$). Naturally, $d$ symbolizes the dimensionality of the
random space. The _temporal function space_ is a continuous $n$-differentiable
function space defined by $\mathscr{T}(n)=C^{n}(\mathfrak{T};\mathbb{R})$. In
other words, this is the space of all functions $f:\mathfrak{T}\to\mathbb{R}$
that have continuous first $n$ derivatives on $\mathfrak{T}$. Moreover, the
_random function space_ , defined by
$\mathscr{Z}=(L^{2}(\Xi,\boldsymbol{\Xi},\mu;\mathbb{R}),\langle\,\cdot\,,\cdot\,\rangle)$,
is a Lebesgue square-integrable space equipped with its standard inner
product. This is the space of all (equivalence classes of) measurable
functions $f:\Xi\to\mathbb{R}$ that are square $\mu$-integrable on $\Xi$. The
standard inner product
$\langle\,\cdot\,,\cdot\,\rangle:L^{2}(\Xi,\boldsymbol{\Xi},\mu;\mathbb{R})\times
L^{2}(\Xi,\boldsymbol{\Xi},\mu;\mathbb{R})\to\mathbb{R}$ is given by $\langle
f,g\rangle=\int fg\,\mathrm{d}\mu$. For computational purposes, this space is
spanned by a complete orthogonal basis,
$\\{\Psi_{j}:\Xi\to\mathbb{R}\\}_{j=0}^{\infty}$ such that $\Psi_{0}(\xi)=1$
for all $\xi\in\Xi$. Finally, the _solution space_ and the _root space_ are
defined by $\mathscr{U}=\mathscr{T}(n)\otimes\mathscr{Z}$ and
$\mathscr{V}=\mathscr{T}(0)\otimes\mathscr{Z}$, respectively, making
$\mathscr{U}$ and $\mathscr{V}$ be two tensor products of vector spaces.

Four remarks are in order. _First_ , it is well-known that $\mathscr{Z}$ forms
a Hilbert space because it is complete under the metric induced by
$\langle\,\cdot\,,\cdot\,\rangle$. _Second_ , if $f\in\mathscr{Z}$, then it
can be represented by the Fourier series $f=\sum_{j=0}^{\infty}f^{j}\Psi_{j}$,
where $f^{j}$ denotes the $j$-th coefficient of the series with the
superscript not symbolizing an exponentiation. _Third_ , the orthogonality
property of the basis in $\mathscr{Z}$ means that
$\langle\Psi_{i},\Psi_{j}\rangle:=\int\Psi_{i}\Psi_{j}\,\mathrm{d}\mu=\langle\Psi_{i},\Psi_{i}\rangle\,\delta_{ij}$,
where $\delta_{ij}$ is the Kronecker delta. _Fourth_ , the dual space of
$\mathscr{Z}$, which is denoted by $\mathscr{Z}^{\prime}$ herein, is the space
spanned by the continuous linear functionals
$\\{\Psi^{i}:\mathscr{Z}\to\mathbb{R}\\}_{i=0}^{\infty}$ defined by
$\Psi^{i}[f]=\langle\Psi_{i},f\rangle/\langle\Psi_{i},\Psi_{i}\rangle\equiv
f^{i}$.

Fig. 1 depicts the relationship between probability space and random space for
$d=2$, and a case scenario of a random space with 9 elements. As shown, the
components of $\xi=(\xi^{1},\ldots,\xi^{d})$ will always be assumed to be
mutually independent, which means that the random domain $\Xi$ is a
rectangular hyper-solid of $d$ dimensions constructed by performing a $d$-fold
Cartesian product of intervals $\bar{\Xi}_{i}:=\xi^{i}(\Omega)$.

Figure 1: Relationship between probability space and random space for $d=2$
(along with a case partition of the random space)

###  2.2 Additional spaces needed in the ME-FSC method

The following two definitions are needed in the development of the ME-FSC
method.

######  Definition 1 (Partitioned random space).

Let $\\{\Xi_{e}\\}_{e=1}^{E}$ be a partition of the random domain, where
$\Xi_{e}\neq\emptyset$ represents the $e$-th element of the partition, and
$E\in\mathbb{N}_{2}$ is the number of random elements employed in the
partition. Let $\boldsymbol{\Xi}_{e}=\boldsymbol{\Xi}\cap\Xi_{e}$ be the
$\sigma$-algebra on $\Xi_{e}$, and $\mu_{e}=\mu|_{\boldsymbol{\Xi}_{e}}$ be
the restriction of $\mu$ to $\boldsymbol{\Xi}_{e}$. Moreover, let
$\hat{\mu}_{e}:\boldsymbol{\Xi}_{e}\to[0,1]$ be the probability measure on
$\boldsymbol{\Xi}_{e}$ given by

$\hat{\mu}_{e}=\frac{\mu_{e}}{\mu(\Xi_{e})},\quad\text{or
equivalently,}\quad\mathrm{d}\hat{\mu}_{e}=\frac{\mathrm{d}\mu_{e}}{\mu(\Xi_{e})}.$

As a result, in this work the _partitioned random space_ is defined as a
finite sequence of disjoint random spaces
$\\{(\Xi_{e},\boldsymbol{\Xi}_{e},\hat{\mu}_{e})\\}_{e=1}^{E}$. A prototype
depiction of such a partition is displayed in Fig. 1.

######  Definition 2 (Partitioned random function space).

If $\\{(\Xi_{e},\boldsymbol{\Xi}_{e},\hat{\mu}_{e})\\}_{e=1}^{E}$ is a
partition of the random space, then let $\\{\mathscr{Z}_{e}\\}_{e=1}^{E}$ be
its associated _partitioned random function space_ , where $\mathscr{Z}_{e}$
represents the $e$-th RFS for the random element $\Xi_{e}$. That is,
$\mathscr{Z}_{e}=(L^{2}(\Xi_{e},\boldsymbol{\Xi}_{e},\hat{\mu}_{e};\mathbb{R});\langle\,\cdot\,,\cdot\,\rangle_{e})$,
where

$\langle\,\cdot\,,\cdot\,\rangle_{e}:L^{2}(\Xi_{e},\boldsymbol{\Xi}_{e},\hat{\mu}_{e};\mathbb{R})\times
L^{2}(\Xi_{e},\boldsymbol{\Xi}_{e},\hat{\mu}_{e};\mathbb{R})\to\mathbb{R}\qquad:\Leftrightarrow\qquad\langle
f,g\rangle_{e}=\int fg\,\mathrm{d}\hat{\mu}_{e}.$

This is the space of all (equivalence classes of) measurable functions
$f:\Xi_{e}\to\mathbb{R}$ that are square $\hat{\mu}_{e}$-integrable on
$\Xi_{e}$. As in $\mathscr{Z}$, the complete orthogonal basis in
$\mathscr{Z}_{e}$,
$\\{\Psi_{j.e}:(\Xi_{e},\boldsymbol{\Xi}_{e})\to(\mathbb{R},\mathcal{B}_{\mathbb{R}})\\}_{j=0}^{\infty}$,
is defined such that $\Psi_{0.e}(\xi)=1$ for all $\xi\in\Xi_{e}$.

Note that the four remarks made in Section 2.1 are also applicable for each of
the $\mathscr{Z}_{e}$’s mentioned above, provided that the following changes
are made:

$\mathscr{Z}\mapsto\mathscr{Z}_{e},\quad\mu\mapsto\hat{\mu}_{e},\quad\langle\,\cdot\,,\cdot\,\rangle\mapsto\langle\,\cdot\,,\cdot\,\rangle_{e},\quad
f\mapsto f_{.e},\quad f^{k}\mapsto
f{{}^{k}}_{.e},\quad\Psi_{k}\mapsto\Psi_{k.e}\quad\text{and}\quad\Psi^{k}\mapsto\Psi{{}^{k}}_{.e}$

with $k$ symbolizing $i$ or $j$.

##  3 Problem statement

For convenience, the same stochastic problem described in [35] is considered
herein. Namely:

Find the real-valued stochastic process $u:\mathfrak{T}\times\Xi\to\mathbb{R}$
in $\mathscr{U}$, such that ($\mu$-a.e.):

$\displaystyle\mathcal{L}[u]=f$ on $\mathfrak{T}\times\Xi$ (1a)
$\displaystyle\big{\\{}\mathcal{B}_{k}[u](0,\cdot\,)=b_{k}\big{\\}}_{k=1}^{n}$
$\displaystyle\qquad\text{on $\\{0\\}\times\Xi$},$ (1b)

where the partial differential operators
$\mathcal{L}:\mathscr{U}\to\mathscr{V}$ and
$\mathcal{B}_{k}:\mathscr{U}\to\mathscr{Z}$ are of order $(n,0)$ and
$(n-1,0)$, respectively. The meaning of $\mathcal{B}_{k}[u](0,\cdot\,)$ is
that once the operator $\mathcal{B}_{k}$ is applied to $u$ the resulting
function is evaluated at $t=0$. The functions
$f:\mathfrak{T}\times\Xi\to\mathbb{R}$ and $b_{k}:\Xi\to\mathbb{R}$ are
elements of $\mathscr{V}$ and $\mathscr{Z}$, respectively, and they are given
by $f=f(t,\xi)$ and $b_{k}=b_{k}(\xi)$.

However, for illustration purposes, we very often consider the much simpler
case111We highlight the fact that this work has been generalized to deal with
problems as complicated as (1) in Section 6. where the mathematical model of
the system can be represented ($\mu$-a.e.) by an undamped oscillator under the
action of an external force:

$\displaystyle\mathcal{L}[u]=f\qquad$ $\displaystyle:\Leftrightarrow\qquad
m\ddot{u}+ku=p$ (1a*)
$\displaystyle\big{\\{}\mathcal{B}_{k}[u](0,\cdot\,)=b_{k}\big{\\}}_{k=1}^{n=2}\qquad$
$\displaystyle:\Leftrightarrow\qquad\big{\\{}u(0,\cdot\,)=\mathscr{u},\,\dot{u}(0,\cdot\,)=\mathscr{v}\big{\\}},$
(1b*)

where $m:\Xi\to\mathbb{R}^{+}$ is the mass of the system,
$k:\Xi\to\mathbb{R}^{+}$ is the stiffness of the system,
$p:\mathfrak{T}\times\Xi\to\mathbb{R}$ is the external force applied to the
system, and $\mathscr{u},\mathscr{v}:\Xi\to\mathbb{R}$ are the prescribed
initial conditions of the system. Observe that $\dot{u}:=\partial_{t}u$ and
$\ddot{u}:=\partial_{t}^{2}u$ are the first and second partial derivatives of
$u$ with respect to time. It is assumed that
$m,k,\mathscr{u},\mathscr{v}\in\mathscr{Z}$, $p\in\mathscr{V}$ and of course
that $u\in\mathscr{U}$.

##  4 Ordinary and enriched stochastic flow maps

The attention is now turned to the two flow maps considered in this
manuscript: the ordinary stochastic flow map, $\varphi$, and the enriched
stochastic flow map, $\hat{\varphi}$.222Please refer to [35] for a
generalization of the two stochastic flow maps considered in this section.

###  4.1 Ordinary stochastic flow map

First of all, the stochastic system given by (1) is assumed to be sufficiently
regular over $\mathfrak{T}\times\Xi$. Thus,

$\displaystyle\partial_{t}^{2}u(t,\xi):=\mathscr{f}(t,\xi,s(t,\xi))=\bar{p}(t,\xi)-\bar{k}(\xi)\,u(t,\xi)$
on $\mathfrak{T}\times\Xi$ (2a)
$\displaystyle\big{\\{}u(0,\xi)=\mathscr{u}(\xi),\,\dot{u}(0,\xi)=\mathscr{v}(\xi)\big{\\}}$
$\displaystyle\qquad\text{on $\\{0\\}\times\Xi$},$ (2b)

where $s=(u,\dot{u})\in\prod_{j=1}^{2}\mathscr{T}(3-j)\otimes\mathscr{Z}$ is
the configuration state of the system over $\mathfrak{T}\times\Xi$,
$\mathscr{f}:\mathfrak{T}\times\Xi\times\mathbb{R}^{2}\to\mathbb{R}$ is a
noisy, non-autonomous function (which can also be regarded as a function in
$\mathscr{V}$) defining the response $\ddot{u}=\partial_{t}^{2}u$,
$\bar{p}=p/m:\mathfrak{T}\times\Xi\to\mathbb{R}$ is the external force applied
to the system per unit mass, and $\bar{k}=k/m:\Xi\to\mathbb{R}^{+}$ is the
stiffness of the system per unit mass. In the context of (1), by ‘sufficiently
regular over $\mathfrak{T}\times\Xi$’ we mean that $\bar{p}\in\mathscr{V}$ and
$\bar{k}\in\mathscr{Z}$ (and thus, $\mathscr{f}\in\mathscr{V}$).

If the solution is further assumed to be analytic on $\mathfrak{T}$ for all
$\xi\in\Xi$, then $u$ can be expanded around $t=t_{i}$ by the Taylor series:

$u(t_{i}+h,\xi)=\sum_{j=0}^{\infty}\frac{h^{j}}{j!}\partial_{t}^{j}u(t_{i},\xi)=\sum_{j=0}^{M}\frac{h^{j}}{j!}\partial_{t}^{j}u(t_{i},\xi)+O(h^{M+1})(\xi),$

where $t_{i}\in\mathfrak{T}$ is the time instant of the simulation, and
$h:=t-t_{i}$ is the associated time-step size (once $t$ is fixed).

Hence, a (local) stochastic flow map of order $M$,
$\varphi(M):\mathbb{R}\times\mathscr{Z}^{2}\to\mathscr{Z}^{2}$, can be
stipulated for (2) to be:

$\varphi(M)(h,s(t_{i},\cdot\,)):=\big{(}u(t_{i}+h,\cdot\,),\dot{u}(t_{i}+h,\cdot\,)\big{)}-O(h^{M+1}),$
(3)

where
$\varphi^{1}(M),\varphi^{2}(M):\mathbb{R}\times\mathscr{Z}^{2}\to\mathscr{Z}$
are given by:

$\displaystyle\varphi^{1}(M)(h,s(t_{i},\cdot\,)):=u(t_{i},\cdot\,)+h\,\dot{u}(t_{i},\cdot\,)+\tfrac{1}{2}h^{2}\,\partial_{t}^{2}u(t_{i},\cdot\,)+\cdots+\tfrac{1}{M!}h^{M}\,\partial_{t}^{M}u(t_{i},\cdot\,)$
$\displaystyle\varphi^{2}(M)(h,s(t_{i},\cdot\,)):=\dot{u}(t_{i},\cdot\,)+h\,\partial_{t}^{2}u(t_{i},\cdot\,)+\tfrac{1}{2}h^{2}\,\partial_{t}^{3}u(t_{i},\cdot\,)+\cdots+\tfrac{1}{M!}h^{M}\,\partial_{t}^{M+1}u(t_{i},\cdot\,).$

It is worth mentioning that the second time derivative
$\partial_{t}^{2}u(t_{i},\cdot\,)$ is computed using (2a), and that the next
time derivatives $\\{\partial_{t}^{j+2}u(t_{i},\cdot\,)\\}_{j=1}^{M-1}$ are
computed using the recursive expression:

$\partial_{t}^{j+2}u(t,\xi):=\mathrm{D}_{t}^{j}\mathscr{f}(t,\xi,s(t,\xi))=\partial_{t}^{j}\bar{p}(t,\xi)-\bar{k}(\xi)\,\partial_{t}^{j}u(t,\xi)\qquad\forall
j\in\\{1,2,\ldots,M-1\\}.$ (4)

In simple terms, the goal of $\varphi(M)$ is to push the system’s state
$s(t_{i},\cdot\,)=(u(t_{i},\cdot\,),\dot{u}(t_{i},\cdot\,))$ through time
using the local information at $t=t_{i}$. We recall that this push is
_forward_ if $h>0$, it is _backward_ if $h<0$, and it is _still_ if $h=0$.

###  4.2 Enriched stochastic flow map

The goal of the associated enriched flow map,
$\hat{\varphi}(M):\mathbb{R}\times\mathscr{Z}^{M+2}\to\mathscr{Z}^{M+2}$, is
to push the system’s state
$s(t_{i},\cdot\,)=(u(t_{i},\cdot\,),\dot{u}(t_{i},\cdot\,))$, and the first
$M-1$ time derivatives of $\mathscr{f}(t_{i},\cdot\,,s(t_{i},\cdot\,))$ at
$t=t_{i}$ (including the function itself). This is the reason why the $k$-th
component of $\hat{\varphi}(M)$ is given by:

$\hat{\varphi}^{k}(M)(h,\hat{s}(t_{i},\cdot\,))=:\hat{s}^{k}(t_{i}+h,\cdot\,)=\begin{cases}\varphi^{k}(M)(h,s(t_{i},\cdot\,))&\text{for
$k\in\\{1,2\\}$}\\\
\mathrm{D}_{t}^{k-3}\mathscr{f}(t_{i}+h,\cdot\,,s(t_{i}+h,\cdot\,))&\text{for
$k\in\\{3,4,\ldots,M+2\\}$}\end{cases}$

where
$\hat{s}=(u,\partial_{t}u,\ldots,\partial_{t}^{M+1}u)\in\prod_{j=1}^{M+2}\mathscr{T}(M-j+2)\otimes\mathscr{Z}$
is the enriched configuration state of the system over
$\mathfrak{T}\times\Xi$. Note that
$\mathrm{D}_{t}^{0}\mathscr{f}:=\mathscr{f}$ is nothing but the function given
by (2a), and that $\\{\mathrm{D}_{t}^{k-3}\mathscr{f}\\}_{k=4}^{M+2}$ is the
set of functions given by (4) with $j=k-3$.

##  5 Overview of the FSC method

The FSC method provides a new alternative for constructing an evolving finite-
dimensional RFS efficiently in time. The construction is said to be efficient
because the dimensionality of the RFS is bounded from below by $n+2$ and from
above by $n+M+1$, regardless of the dimensionality of the probability space,
and where $n$ denotes the order of the governing ODE with respect to time, and
$M$ the order of the stochastic flow map.

To illustrate the key idea behind the FSC method333For more details, please
refer to [35, 36]., consider the stochastic dynamical system defined by (1).
From (3), we learn that the state of the system around $t=t_{i}$ can be
expanded as

$u(t,\xi)=\sum_{j=0}^{M}\frac{(t-t_{i})^{j}}{j!}\partial_{t}^{j}u(t_{i},\xi)\qquad\text{and}\qquad\dot{u}(t,\xi)=\sum_{j=0}^{M}\frac{(t-t_{i})^{j}}{j!}\partial_{t}^{j+1}u(t_{i},\xi),$
(5)

where $h=t-t_{i}$, but with $t$ not necessarily fixed.

From (5) we discover that the solution actually consists of two parts: the
deterministic part $(t-t_{i})^{j}/j!\in\mathscr{T}$, and the non-deterministic
part $\partial_{t}^{j}u(t_{i},\xi)\in\mathscr{Z}$. If, for efficiency reasons,
$\\{\partial_{t}^{j}u(t_{i},\cdot\,)\\}_{j=0}^{M+1}$ is orthogonalized with
respect to the measure defined in $\mathscr{Z}$, then $M+2$ orthogonal
functions are obtained, namely: $\\{\Psi_{j}\\}_{j=1}^{M+2}$. Hence, a more
convenient way of writing (5) would be

$u(t,\xi)=\sum_{j=1}^{M+2}u^{j}(t)\,\Psi_{j}(\xi)\qquad\text{and}\qquad\dot{u}(t,\xi)=\sum_{j=1}^{M+2}\dot{u}^{j}(t)\,\Psi_{j}(\xi).$
(6)

However, since the identically-equal-to-one function may not always be an
element of $\\{\Psi_{j}\\}_{j=1}^{M+1}$, this calls to rewriting (6) in the
following manner:

$u(t,\xi)=\sum_{j=0}^{M+2}u^{j}(t)\,\Psi_{j}(\xi)\qquad\text{and}\qquad\dot{u}(t,\xi)=\sum_{j=0}^{M+2}\dot{u}^{j}(t)\,\Psi_{j}(\xi)$
(7)

with $\Psi_{0}\equiv 1$ symbolizing the identically-equal-to-one function, as
prescribed in Section 2.1.

Put differently, in FSC the orthogonalization process is always started with
$\Psi_{0}\equiv 1$ and then the germs
$\\{\partial_{t}^{j}u(t_{i},\cdot\,)\\}_{j=0}^{M+1}$ are orthogonalized with
respect to each other to produce $M+3$ (orthogonal) basis vectors. The
resulting set of basis vectors is then used to construct a _complete_ RFS for
the system’s state at $t=t_{i}$:

$\mathscr{Z}^{[M+2]}=\mathrm{span}\\{\Psi_{j}\\}_{j=0}^{M+2}.$

######  Remark 1.

For the more general case, i.e. for a dynamical system with an $n$-tuple state
and driven by a stochastic flow map of order $M$, the aforementioned RFS is
constructed as $\mathscr{Z}^{[n+M]}=\mathrm{span}\\{\Psi_{j}\\}_{j=0}^{n+M}$.

However, since an exact representation of the system’s state is not always
needed to obtain accurate results, we highlight the fact that the FSC method
is also capable of considering as few as $n+2$ basis vectors, but at the cost
of losing some precision along the way. Naturally, the longer the duration of
the simulation, the greater this loss of precision will be. In fact, this loss
of precision is noticeable in the numerical problems explored in Section 7,
because the duration of the simulations was chosen so that this loss would
show up in the numerical results. Therefore, a more effective way to write the
RFS required for the simulation is to simply define it as

$\mathscr{Z}^{[P]}=\mathrm{span}\\{\Psi_{j}\\}_{j=0}^{P}\quad\text{with}\quad
P\in\\{n+1,n+2,\ldots,M+2\\}.$

##  6 Multi-element flow-driven spectral chaos (ME-FSC) method

This section is devoted to presenting the multi-element version of the FSC
method in detail. First, we review the multi-element concept behind the ME-FSC
method, and then we present the computation of the (global) probability
moments using the local information available on each random element. The
section is concluded with an outline of the ME-FSC scheme used in this work.

###  6.1 Overview of the ME-FSC method

Because the solution is in $\mathscr{U}$ by assumption, then $u$ can be
represented with the Fourier series:

$u(t,\xi)=\sum_{j=0}^{\infty}u^{j}(t)\,\Psi_{j}(\xi)\quad\text{on
$\mathfrak{T}\times\Xi$},$ (8)

where $u^{j}$ is a temporal function in $\mathscr{T}(2)$ symbolizing the
$j$-th random mode of $u$.

Now, let $\\{(\Xi_{e},\boldsymbol{\Xi}_{e},\hat{\mu}_{e})\\}_{e=1}^{E}$ be a
partitioned random space, and let $\\{\mathscr{Z}_{e}\\}_{e=1}^{E}$ be its
associated partitioned RFS. Then, a $p$-discretization for each of these
$\mathscr{Z}_{e}$’s can be stipulated by letting
$\mathscr{Z}_{e}^{[P_{e}]}=\mathrm{span}\\{\Psi_{j.e}\\}_{j=0}^{P_{e}}$ be a
finite subspace of $\mathscr{Z}_{e}$ with $P_{e}+1\in\mathbb{N}_{1}$ denoting
the dimensionality of $\mathscr{Z}_{e}^{[P_{e}]}$.

Moreover, if $u_{.e}(t,\cdot\,)=u(t,\cdot\,)|_{\Xi_{e}}$ is defined to be the
restriction of $u(t,\cdot\,)$ to $\Xi_{e}$ for all $t\in\mathfrak{T}$, then an
approximation $u_{.e}^{[P_{e}]}(t,\cdot\,)$ of $u_{.e}(t,\cdot\,)$ can be
represented in $\mathscr{Z}_{e}^{[P_{e}]}$ by

$u_{.e}(t,\xi)\approx
u_{.e}^{[P_{e}]}(t,\xi)=\sum_{j=0}^{P_{e}}u{{}^{j}}_{\\!.e}(t)\,\Psi_{j.e}(\xi)\equiv
u{{}^{j}}_{\\!.e}(t)\,\Psi_{j.e}(\xi)\quad\text{on
$\mathfrak{T}\times\Xi_{e}$},$ (9)

where the summation sign has been omitted in the last equality for notational
simplicity, and the summation index is taken over $j\in\\{0,1,\ldots,P_{e}\\}$
unless indicated otherwise. Similarly as before, $u{{}^{j}}_{\\!.e}$
represents a temporal function in $\mathscr{T}(2)$ symbolizing the $j$-th
random mode of $u_{.e}$.

######  Remark 2.

A dot symbol is introduced in (9) to distinguish between tensor indices and
identification indices (aka ID indices). Consequently, indices $j$ and $e$ are
to be interpreted herein as tensor and identification indices, respectively.
Moreover, no hidden summation signs should ever be expected for identification
indices.

######  Remark 3.

In writing (9), $\\{\Psi_{j.e}\\}_{j=0}^{\infty}$ was assumed to be a well-
graded orthogonal basis in order to ensure that the approximation of $u_{.e}$
could be carried out that way. This assumption does not represent an issue in
this work, since all random bases utilized in ME-FSC are well-graded by
construction.

The problem now reduces to find $E$ independent solutions to (1) by looping
across the random domain from $e=1$ to $e=E$ and using the following
procedure.

Substituting (9) into (1) gives

$\displaystyle\mathcal{L}[u{{}^{j}}_{\\!.e}\Psi_{j.e}]=f$ on
$\mathfrak{T}\times\Xi_{e}$ (10a)
$\displaystyle\big{\\{}\mathcal{B}_{k}[u{{}^{j}}_{\\!.e}\Psi_{j.e}](0,\cdot\,)=b_{k}\big{\\}}_{k=1}^{n}$
$\displaystyle\qquad\text{on $\\{0\\}\times\Xi_{e}$}.$ (10b)

Projecting (10) onto $\mathscr{Z}_{e}^{[P_{e}]}$ yields a system of $P_{e}+1$
ordinary differential equations of second order in the variable $t$, where the
unknowns are the random modes
$u\indices{{}^{j}_{.e}}=u\indices{{}^{j}_{.e}}(t)$ and their first $n-1$ time
derivatives:

$\displaystyle\Psi{{}^{i}}_{.e}\big{[}\mathcal{L}[u{{}^{j}}_{\\!.e}\Psi_{j.e}]\big{]}=\Psi{{}^{i}}_{.e}[f]$
on $\mathfrak{T}$ (11a)
$\displaystyle\big{\\{}\Psi{{}^{i}}_{.e}\big{[}\mathcal{B}_{k}[u{{}^{j}}_{\\!.e}\Psi_{j.e}](0,\cdot\,)\big{]}=\Psi{{}^{i}}_{.e}[b_{k}]\big{\\}}_{k=1}^{n}$
on $\\{0\\}$ (11b)

with $i,j\in\\{0,1,\ldots,P_{e}\\}$. Notice that this projection is done here
by applying on both sides of each equation the linear functionals
$\\{\Psi{{}^{i}}_{.e}\in\mathscr{Z}_{e}^{\prime}\\}_{i=0}^{P_{e}}$ one by one.

For the specific case of a system given by (1), it is clear that (11) would
reduce to

$\displaystyle
m\indices{{}^{i}_{j.e}}\ddot{u}\indices{{}^{j}_{.e}}+k\indices{{}^{i}_{j.e}}u\indices{{}^{j}_{.e}}=p\indices{{}^{i}_{.e}}$
on $\mathfrak{T}$ (11a*)
$\displaystyle\big{\\{}u\indices{{}^{i}_{.e}}(0)=\mathscr{u}\indices{{}^{i}_{.e}},\,\dot{u}\indices{{}^{i}_{.e}}(0)=\mathscr{v}\indices{{}^{i}_{.e}}\big{\\}}$
$\displaystyle\qquad\text{on $\\{0\\}$},$ (11b*)

where $i,j\in\\{0,1,\ldots,P_{e}\\}$ (summation sign implied only over
repeated index $j$), and:

$\displaystyle
m\indices{{}^{i}_{j.e}}=\langle\Psi_{i.e},m\Psi_{j.e}\rangle_{e}/\langle\Psi_{i.e},\Psi_{i.e}\rangle_{e},\quad
k\indices{{}^{i}_{j.e}}=\langle\Psi_{i.e},k\Psi_{j.e}\rangle_{e}/\langle\Psi_{i.e},\Psi_{i.e}\rangle_{e},$
$\displaystyle
p\indices{{}^{i}_{.e}}(t)=\langle\Psi_{i.e},p(t,\cdot\,)\rangle_{e}/\langle\Psi_{i.e},\Psi_{i.e}\rangle_{e},$
$\displaystyle\mathscr{u}\indices{{}^{i}_{.e}}=\langle\Psi_{i.e},\mathscr{u}\rangle_{e}/\langle\Psi_{i.e},\Psi_{i.e}\rangle_{e}\quad\text{and}\quad\mathscr{v}\indices{{}^{i}_{.e}}=\langle\Psi_{i.e},\mathscr{v}\rangle_{e}/\langle\Psi_{i.e},\Psi_{i.e}\rangle_{e}.$

###  6.2 Computation of probability moments

In modeling notation, the system specified by (1) can be written as:

$y=\boldsymbol{\mathcal{M}}[u][x]\quad\text{subject to initial
condition}\quad\boldsymbol{\mathcal{I}}[u],\\\ $

where $\boldsymbol{\mathcal{M}}[u]:\mathscr{V}^{r}\to\mathscr{V}^{s}$ is as in
(1a) and denotes the mathematical model of the system,
$x=(x_{1},\ldots,x_{r}):\mathfrak{T}\times\Xi\to\mathbb{R}^{r}$ is the
$r$-tuple input of $\boldsymbol{\mathcal{M}}[u]$,
$y=(y_{1},\ldots,y_{s}):\mathfrak{T}\times\Xi\to\mathbb{R}^{s}$ is the
$s$-tuple output of $\boldsymbol{\mathcal{M}}[u]$, and
$\boldsymbol{\mathcal{I}}[u]$ is as in (1b) and denotes the initial condition
of the system.

Now, define $z=y_{k}$ to be the $k$-th component of output
$y=\boldsymbol{\mathcal{M}}[u][x]$. If $z\in\mathscr{V}$, then it is clear
that such a function can be represented approximately over the $e$-th random
element using the following expansion:

$z_{.e}(t,\xi)\approx
z_{.e}^{[P_{e}]}(t,\xi)=\sum_{j=0}^{P_{e}}z{{}^{j}}_{\\!.e}(t)\,\Psi_{j.e}(\xi)\equiv
z{{}^{j}}_{\\!.e}(t)\,\Psi_{j.e}(\xi)\quad\text{on
$\mathfrak{T}\times\Xi_{e}$},$

where $z_{.e}(t,\cdot\,)=z(t,\cdot\,)|_{\Xi_{e}}$ denotes the restriction of
$z(t,\cdot\,)$ to $\Xi_{e}$ for all $t\in\mathfrak{T}$.

In this section, we are interested in computing the mean and variance of $z$
over $\mathfrak{T}\times\Xi$ using the local information available on each
random element. To do so, we first make the following observation regarding
the (local) mean and (local) variance of $z$ given that $\xi\in\Xi_{e}$. As it
should be easy to verify, the (local) mean of $z$ given that $\xi\in\Xi_{e}$,
$\mathbf{E}_{e}[z]:\mathfrak{T}\to\mathbb{R}$, is nothing but

$\mathbf{E}_{e}[z](t)\equiv\mathbf{E}[z\mid\Xi_{e}](t):=\int
z(t,\cdot\,)\,\mathrm{d}\hat{\mu}_{e}=\int
z_{.e}(t,\cdot\,)\,\mathrm{d}\hat{\mu}_{e}=z{{}^{0}}_{.e}(t),$ (13)

and that the (local) variance of $z$ given that $\xi\in\Xi_{e}$,
$\mathrm{Var}_{e}[z]:\mathfrak{T}\to\mathbb{R}_{0}^{+}$, is

$\mathrm{Var}_{e}[z](t)\equiv\mathrm{Var}[z\mid\Xi_{e}](t):=\int\big{(}z(t,\cdot\,)-\mathbf{E}_{e}[z](t)\big{)}^{2}\,\mathrm{d}\hat{\mu}_{e}=\int
z^{2}(t,\cdot\,)\,\mathrm{d}\hat{\mu}_{e}-\mathbf{E}_{e}[z]^{2}(t)\\\ =\int
z{{}_{.e}}^{2}(t,\cdot\,)\,\mathrm{d}\hat{\mu}_{e}-\mathbf{E}_{e}[z]^{2}(t)=\sum_{j=1}^{P_{e}}\langle\Psi_{j.e},\Psi_{j.e}\rangle_{e}\,z{{}^{j}}_{\\!.e}(t)\,z{{}^{j}}_{\\!.e}(t).$
(14)

Therefore, the (global) mean of $z$,
$\mathbf{E}[z]:\mathfrak{T}\to\mathbb{R}$, which is given by

$\mathbf{E}[z](t):=\int
z(t,\cdot\,)\,\mathrm{d}\mu=\sum_{e=1}^{E}\int_{\Xi_{e}}z(t,\cdot\,)\,\mathrm{d}\mu=\sum_{e=1}^{E}\int
z(t,\cdot\,)\,\mathrm{d}\mu_{e},$

simplifies to:

$\mathbf{E}[z](t)=\sum_{e=1}^{E}\mu(\Xi_{e})\int
z(t,\cdot\,)\,\mathrm{d}\hat{\mu}_{e}=\sum_{e=1}^{E}\mu(\Xi_{e})\,\mathbf{E}_{e}[z](t).$
(15)

This final form of $\mathbf{E}[z]$ is also known as the _law of total
expectation_.

Moreover, the (global) variance of $z$,
$\mathrm{Var}[z]:\mathfrak{T}\to\mathbb{R}_{0}^{+}$, is given by

$\mathrm{Var}[z](t):=\int\big{(}z(t,\cdot\,)-\mathbf{E}[z](t)\big{)}^{2}\,\mathrm{d}\mu=\int
z^{2}(t,\cdot\,)\,\mathrm{d}\mu-\mathbf{E}[z]^{2}(t)=\sum_{e=1}^{E}\int_{\Xi_{e}}z^{2}(t,\cdot\,)\,\mathrm{d}\mu-\mathbf{E}[z]^{2}(t),$
(16)

The resulting expression can be further simplified if we recognize that
_first_ :

$\int_{\Xi_{e}}z^{2}(t,\cdot\,)\,\mathrm{d}\mu=\mu(\Xi_{e})\int
z^{2}(t,\cdot\,)\,\mathrm{d}\hat{\mu}_{e}=\mu(\Xi_{e})\int
z{{}_{.e}}^{2}(t,\cdot\,)\,\mathrm{d}\hat{\mu}_{e}=\mu(\Xi_{e})\,\big{(}\mathrm{Var}_{e}[z](t)+\mathbf{E}_{e}[z]^{2}(t)\big{)}$

and that _second_ :

$\mathbf{E}[z]^{2}(t)=\sum_{e_{1}=1}^{E}\sum_{e_{2}=1}^{E}\mu(\Xi_{e_{1}})\,\mu(\Xi_{e_{2}})\,\mathbf{E}_{e_{1}}[z](t)\,\mathbf{E}_{e_{2}}[z](t)\\\
=\sum_{e=1}^{E}\mu(\Xi_{e})^{2}\,\mathbf{E}_{e}[z]^{2}(t)+2\sum_{e_{1}=2}^{E}\sum_{e_{2}=1}^{e_{1}-1}\mu(\Xi_{e_{1}})\,\mu(\Xi_{e_{2}})\,\mathbf{E}_{e_{1}}[z](t)\,\mathbf{E}_{e_{2}}[z](t).$

Hence, by replacing these two expressions into (16) and simplifying, we get

$\mathrm{Var}[z](t)=\sum_{e=1}^{E}\mu(\Xi_{e})\,\mathrm{Var}_{e}[z](t)\\\
+\sum_{e=1}^{E}\mu(\Xi_{e})\,(1-\mu(\Xi_{e}))\,\mathbf{E}_{e}[z]^{2}(t)-2\sum_{e_{1}=2}^{E}\sum_{e_{2}=1}^{e_{1}-1}\mu(\Xi_{e_{1}})\,\mu(\Xi_{e_{2}})\,\mathbf{E}_{e_{1}}[z](t)\,\mathbf{E}_{e_{2}}[z](t).$
(17)

This final form of $\mathrm{Var}[z]$ is also known as the _law of total
variance_.

###  6.3 ME-FSC scheme

Consider the stochastic dynamical system given by (1). Let
$\\{\mathfrak{T}_{i}\times\Xi_{e}\\}_{i=0,e=1}^{N-1,E}$ be a partition of the
temporal-random domain $\mathfrak{T}\times\Xi$ (aka system’s domain), where
$\mathfrak{T}_{i}\neq\emptyset$ is the $i$-th element of the temporal domain,
$\Xi_{e}\neq\emptyset$ is the $e$-th element of the random domain, and
$\mathfrak{T}_{i}\times\Xi_{e}$ is the $(i,e)$-th element of the system’s
domain. It is worth mentioning that this partition gives rise to the
partitioned random space
$\\{(\Xi_{e},\boldsymbol{\Xi}_{e},\hat{\mu}_{e})\\}_{e=1}^{E}$ and associated
partitioned RFS $\\{\mathscr{Z}_{e}\\}_{e=1}^{E}$, as defined in Section 2.
Moreover, for notational convenience, we define
$s_{.ie}=s|_{\mathrm{cl}(\mathfrak{T}_{i})\times\Xi_{e}}$ to be the
restriction of $s$ to
$\mathfrak{R}_{ie}:=\mathrm{cl}(\mathfrak{T}_{i})\times\Xi_{e}$, and
$s_{.\cdot e}=s|_{\mathfrak{T}\times\Xi_{e}}$ to be the restriction of $s$ to
$\mathfrak{T}\times\Xi_{e}$. For illustration purposes, Fig. 2 depicts the
evolution of a dynamical system by applying successively a stochastic flow map
of order $M$ over a partitioned random domain with $E$ elements.

Figure 2: Evolution of a dynamical system via a stochastic flow map of order
$M$ (with $h_{i}>0$) over a partitioned random domain

In this sense, if we assume that the system is driven by a stochastic flow map
of order $M$, we can proceed as follows. (Recall that $n$ denotes the order of
the governing ODE with respect to time, as displayed in (1).)

  1. 1.

Loop across the temporal domain from $i=0$ to $i=N-1$.

    1. i.

Loop across the random domain from $e=1$ to $e=E$. Note that this loop can be
parallelized from a computational standpoint since each iteration is
independent of the others.

      1. (a)

Define a solution representation for the configuration state $s_{.ie}$ in the
following way.

         * •

Take $\Phi_{0.ie}\equiv 1$ and
$\\{\Phi_{j.ie}:=\hat{\varphi}^{j}(M)(0,\hat{s}_{.\cdot
e}(t_{i},\cdot\,))\\}_{j=1}^{P_{e}}$ to be an ordered set of linearly
independent functions in $\mathscr{Z}_{e}$ with $n+1\leq P\leq n+M$. Note that
$\hat{\varphi}(M)(0,\hat{s}_{.\cdot
e}(t_{i},\cdot\,))\equiv\hat{s}_{.ie}(t_{i},\cdot\,)=\hat{s}_{.(i-1)e}(t_{i},\cdot\,)$
for $i\geq 1$. However, if $i=0$, then $\hat{\varphi}(M)(0,\hat{s}_{.\cdot
e}(t_{0},\cdot\,))\equiv\hat{s}_{.\cdot e}(0,\cdot\,)$.

         * •

Orthogonalize the set $\\{\Phi_{j.ie}\\}_{j=0}^{P_{e}}$ using the Gram-Schmidt
process [37], so that the resulting set $\\{\Psi_{j.ie}\\}_{j=0}^{P_{e}}$ is
an orthogonal basis in $\mathscr{Z}_{e}$. Specifically, for
$j\in\\{0,1,\ldots,P_{e}\\}$:

$\Psi_{j.ie}:=\Phi_{j.ie}-\sum_{k=0}^{j-1}\frac{\langle\Phi_{j.ie},\Psi_{k.ie}\rangle_{e}}{\langle\Psi_{k.ie},\Psi_{k.ie}\rangle_{e}}\Psi_{k.ie}.$

         * •

Define
$\mathscr{Z}_{ie}^{[P_{e}]}=\mathrm{span}\\{\Psi_{j.ie}\\}_{j=0}^{P_{e}}$ to
be an $(h,p)$-discretization of $\mathscr{Z}$ over the region
$\mathfrak{R}_{ie}$. Then, because $\mathscr{Z}_{ie}^{[P_{e}]}$ is an evolving
function space over $\Xi_{e}$, expansion (9) is now to be read as:

$u_{.ie}(t,\xi)\approx
u^{[P_{e}]}_{.ie}(t,\xi)=\sum_{j=0}^{P_{e}}u{{}^{j}}_{\\!.ie}(t)\,\Psi_{j.ie}(\xi)\equiv
u{{}^{j}}_{\\!.ie}(t)\,\Psi_{j.ie}(\xi).$ (9’)

Thus, to compute the $l$-th component of the configuration state,
$s{{}^{l}}_{\\!.ie}$, we simply need to take the $(l-1)$-th time derivative of
(9’). Here $l\in\\{1,2,\ldots,n\\}$.

      2. (b)

Transfer at $t=t_{i}$ the random modes from the old definition of the system’s
configuration state

$s_{.(i-1)e}(t_{i},\cdot\,)=(u_{.(i-1)e}(t_{i},\cdot\,),\partial_{t}u_{.(i-1)e}(t_{i},\cdot\,),\ldots,\partial_{t}^{n-1}u_{.(i-1)e}(t_{i},\cdot\,))$

to the new definition of the system’s configuration state

$s_{.ie}(t_{i},\cdot\,)=(u_{.ie}(t_{i},\cdot\,),\partial_{t}u_{.ie}(t_{i},\cdot\,),\ldots,\partial_{t}^{n-1}u_{.ie}(t_{i},\cdot\,)),$

given that $i\geq 1$. To do so, we implement the FSC-2 approach presented in
[35] to obtain the random modes of each of the components of $s_{.ie}$ at
$t=t_{i}$. The reason why the FSC-2 is chosen over the FSC-1 is that the
former can transfer the probability information exactly at any given time.
Thus, by resorting to Theorem 1 of [35], one can show that the $j$-th random
mode of the $l$-th component of $s_{.ie}(t_{i},\cdot\,)$ is given by:

$(s^{l}){{}^{j}}_{\\!.ie}(t_{i})=\begin{cases}\mathbf{E}[\Phi_{l.ie}]&\text{for
$j=0$}\\\\[5.38193pt]
\displaystyle\frac{\det\triangle_{j}(l)}{\det\square_{j}}&\text{for
$0<j<l$}\\\\[5.38193pt] 1&\text{for $j=l$}\\\ 0&\text{otherwise,}\end{cases}$
(11b’)

from where we have taken $\\{\Phi_{l.ie}:=\varphi^{l}(M)(0,s_{.\cdot
e}(t_{i},\cdot\,))\equiv
s{{}^{l}}_{\\!.ie}(t_{i},\cdot\,)=s{{}^{l}}_{\\!.(i-1)e}(t_{i},\cdot\,)\\}_{l=1}^{n}$,
and

$\displaystyle\square_{j}=$
$\displaystyle\begin{bmatrix}\mathrm{Cov}[\Phi_{1.ie},\Phi_{1.ie}]&\cdots&\mathrm{Cov}[\Phi_{1.ie},\Phi_{j.ie}]\\\
\vdots&\ddots&\vdots\\\
\mathrm{Cov}[\Phi_{j.ie},\Phi_{1.ie}]&\cdots&\mathrm{Cov}[\Phi_{j.ie},\Phi_{j.ie}]\end{bmatrix}\in\mathscr{M}(j\times
j,\mathbb{R}),$ $\displaystyle\triangle_{j}(l)=$
$\displaystyle\begin{bmatrix}\mathrm{Cov}[\Phi_{1.ie},\Phi_{1.ie}]&\cdots&\mathrm{Cov}[\Phi_{1.ie},\Phi_{j.ie}]\\\
\vdots&\ddots&\vdots\\\
\mathrm{Cov}[\Phi_{(j-1).ie},\Phi_{1}]&\cdots&\mathrm{Cov}[\Phi_{(j-1).ie},\Phi_{j.ie}]\\\
\mathrm{Cov}[\Phi_{l.ie},\Phi_{1.ie}]&\cdots&\mathrm{Cov}[\Phi_{l.ie},\Phi_{j.ie}]\end{bmatrix}\in\mathscr{M}(j\times
j,\mathbb{R})$

with $\triangle_{1}(l)=\mathrm{Cov}[\Phi_{l.ie},\Phi_{1.ie}]\in\mathbb{R}$ and

$\triangle_{2}(l)=\begin{bmatrix}\mathrm{Cov}[\Phi_{1.ie},\Phi_{1.ie}]&\mathrm{Cov}[\Phi_{1.ie},\Phi_{2.ie}]\\\
\mathrm{Cov}[\Phi_{l.ie},\Phi_{1.ie}]&\mathrm{Cov}[\Phi_{l.ie},\Phi_{2.ie}]\end{bmatrix}\in\mathscr{M}(2\times
2,\mathbb{R})$

as special cases of $\triangle_{j}$ by definition.

      3. (c)

Substitute (9’) into (1) to obtain (10).

      4. (d)

Project (10a) onto $\mathscr{Z}_{ie}^{[P_{e}]}$ to obtain (11a) subject to
(11b’). Note that if $i=0$, (11a) is subject to (11b).

      5. (e)

Integrate (11) over time, provided that a suitable time integration method has
been chosen for the system of equations in hand. This step requires finding
the random modes of each of the components of the configuration state
$s_{.ie}$ at $t=t_{i+1}$; that is,
$\\{(s^{l}){{}^{j}}_{.ie}(t_{i+1})\\}_{l=1,j=0}^{n,P_{e}}$

      6. (f)

Compute both the (local) mean and the (local) variance of each of the
components of output $y=\boldsymbol{\mathcal{M}}[x][u]$ over
$\mathfrak{R}_{ie}$, by resorting to the formulas prescribed by (13) and (14).

    2. ii.

Aggregate results to compute over
$\mathfrak{R}_{i}:=\mathfrak{T}_{i}\times\Xi$ the (global) mean and (global)
variance of $y=\boldsymbol{\mathcal{M}}[x][u]$ using the formulas stipulated
by (15) and (17).

  2. 2.

Post-process results.

##  7 Numerical results

The local and global errors, $\epsilon:\mathscr{T}\to\mathscr{T}$ and
$\epsilon_{G}:\mathscr{T}\to\mathbb{R}$, are defined using the following
expressions:

$\displaystyle\epsilon[f](t)=|f(t)-f_{\mathrm{exact}}(t)|$
$\displaystyle\epsilon_{G}[f]=\frac{1}{T}\int_{\mathfrak{T}}|f(t)-f_{\mathrm{exact}}(t)|\,\mathrm{d}t\approx\frac{\Delta
t}{T}\sum_{i=0}^{N}|f(t_{i})-f_{\mathrm{exact}}(t_{i})|,$

where $t_{i}\in\mathfrak{T}$ is the time instant of the simulation, $\Delta t$
is the associated time-step size, and $N$ is the number of time steps utilized
in the simulation (such that $t_{0}=0$ and $t_{N}=N\,\Delta t=T$).

In an effort to reduce significantly the source of errors coming from the
discretization of $\mathscr{T}$, the time-step size used in the simulations is
taken as $\Delta t=0.001$ s for Problems 1 and 2 and $\Delta t=0.005$ s for
Problems 3 and 4. To integrate (11) over time, we use the RK4 method over each
random element, and in order to obtain accurate results, the random function
space is updated at every time step. For simplicity, the partition of the
random domain is such that all elements are the same size. The temporal domain
employed in Problems 1, 2 and 3 is $\mathfrak{T}=[0,150]$ s, and in Problem 4
is $\mathfrak{T}=[0,50]$ s. To prevent ill-conditioned matrices in the early
times of a simulation with deterministic initial conditions, the gPC method is
used till the first second of the simulation. This is to ensure that the
system’s state is well-developed for the analysis with ME-FSC. Specifically,
for Problems 1 and 2 we employ the gPC method with $P=7$, and for Problem 3 we
employ the gPC method with $P=9$.

The inner products are computed using a set of Legendre quadrature rules on
each random axis defined by

$\mathrm{Uniform}\sim\text{Gauss-Legendre (10 points/element)}.$

This means that for distributions other than uniform, the probability
distribution function must be included in the integrand in order to correctly
obtain the numerical value of the inner product.

All problems are run using Apple’s Foundation and Accelerate frameworks [38]
on a 2020 MacBook Air with Apple M1 chip (_8-Core_ CPU at _3.20 GHz_ ,
_8-Core_ GPU, _16-Core_ Neural Engine, and _16 GB_ unified memory) and _1 TB
Apple-Fabric SSD_ storage (APFS-formatted), running macOS Big Sur (version
11.2). The code is written entirely in the Swift 5.3 language [39].

###  7.1 Problem 1: A linear system governed by a 2nd-order stochastic ODE

We first consider the problem of an undamped single-degree-of-freedom system
under free vibration. The law of motion for this system is defined by

$m\ddot{u}+ku=0,$

where the mass of the system is $m=100$ kg, and the stiffness of the system,
$k:\Xi\to\mathbb{R}^{+}$, is stochastic and given by $k(\xi)=\xi$. The system
has an initial displacement of $u(0,\cdot\,)=\mathscr{u}\equiv 0.05$ m, and an
initial velocity of $\dot{u}(0,\cdot\,)=\mathscr{v}\equiv 0.20$ m/s. Formally,
one can express this problem in the following way:

Find the displacement of the system $u:\mathfrak{T}\times\Xi\to\mathbb{R}$ in
$\mathscr{U}$, such that ($\mu$-a.e.):

$\displaystyle m\ddot{u}+ku=0$ on $\mathfrak{T}\times\Xi$ (18a)
$\displaystyle\big{\\{}u(0,\cdot\,)=\mathscr{u},\,\dot{u}(0,\cdot\,)=\mathscr{v}\big{\\}}$
$\displaystyle\qquad\text{on $\\{0\\}\times\Xi$}.$ (18b)

This problem statement is similar to (1), with the only difference being that
$p\equiv 0$ and $m,\mathscr{u},\mathscr{v}$ are real numbers.

Three different probability distributions are considered for $\xi$. The first
distribution is a _uniform distribution_ ,
$\mathrm{Uniform}\sim\xi\in\Xi=[a,b]$, the second distribution is a _beta
distribution_ , $\mathrm{Beta}(\alpha_{1},\beta_{1})\sim\xi\in\Xi=[a,b]$, and
the third distribution is a _gamma distribution_ ,
$\mathrm{Gamma}(\alpha_{2},\beta_{2})\sim\xi\in\Xi=[a,\infty)$. Because gamma
is a distribution with non-compact support, the simulations are run with
_truncated gamma distribution_ ,
$\mathrm{Gamma}(\alpha_{2},\beta_{2},b_{2})\sim\xi\in\Xi=[a,b_{2}]$. But
still, the exact solution is obtained with the gamma distribution. The
parameters for each of the distributions mentioned above are taken as:
$(a,b)=(340,460)$ N/m, $(\alpha_{1},\beta_{1})=(2,5)$ and
$(\alpha_{2},\beta_{2},b_{2})=(10,0.1,920\text{ N/m})$.

In Fig. 3 we show the evolution of the mean and variance of the system’s
displacement. These results are obtained using 8 elements in the random domain
and 7 basis vectors per element. We see that when this particular
discretization is used, the ME-FSC method is able to reproduce the exact
response with high fidelity, explaining why the two plots look
indistinguishable from each other.

(a) Mean

(b) Variance

Figure 3: _Problem 1_ — Evolution of $\mathbf{E}[u]$ and $\mathrm{Var}[u]$ for
the case when the $(h,p)$-discretization level of RFS is $(P,E)=(6,8)$ and
$\mu\sim\mathrm{Uniform}$

Fig. 4 presents the local errors in mean and variance of the system’s
displacement using different numbers of elements and basis vectors. For
brevity, we only present the case when the probability measure is uniform,
although similar convergence trends are also achievable when the probability
measure is beta or gamma. For sake of comparison, we also include the case
when $P=2$, even though this is not allowed by the FSC method.444This is
because in FSC the lower bound for $P$ is always $n+1$, where $n$ is the order
of the governing ODE with respect to time. These plots show that the results
get better if the number of elements is increased. However, this is not always
the case whenever the number of basis vectors is increased. For instance, when
the number of basis vectors is increased from 3 to 5, we observe that the
results improve noticeably, but when the number of basis vectors is increased
from 5 to 7 they do not. This is chiefly because the quadrature rule employed
in this work is not optimal—recall that Gaussian-based quadrature rules only
ensure exponential convergence to the sought integral if the selected
quadrature points are the roots of an orthogonal polynomial that is concordant
with the measure defined on the integral’s domain (which is manifestly not the
case here). The above assertion even holds for the case when 1 element is used
in the simulations (i.e. when there is no partition in the random domain),
since distributing only 10 Legendre quadrature points in the random domain is
not sufficient enough to allow the higher basis vectors to play a role in the
global response of the system. However, if more quadrature points were to be
employed in the simulations, an improvement between the plots with $P=4$ and
$P=6$ could have been discerned for the case of using 1 element, but virtually
no improvement for the case of using more than 1 element. It is worth pointing
out that the accuracy of these results is limited by machine precision. As a
result, better results than those depicted in Fig. 4 are difficult to obtain
for other values of $P$ and $E$.

(a) Mean error for $\mathscr{Z}^{[2]}$

(b) Variance error for $\mathscr{Z}^{[2]}$

(c) Mean error for $\mathscr{Z}^{[4]}$

(d) Variance error for $\mathscr{Z}^{[4]}$

(e) Mean error for $\mathscr{Z}^{[6]}$

(f) Variance error for $\mathscr{Z}^{[6]}$

Figure 4: _Problem 1_ — Local error evolution of $\mathbf{E}[u]$ and
$\mathrm{Var}[u]$ for different $(h,p)$-discretization levels of RFS and for
$\mu\sim\mathrm{Uniform}$

Fig. 5 depicts the convergence of global errors as a function of the number of
elements and the number of basis vectors used. Included in this figure are the
cases where the probability measure is uniform, beta or gamma. Overall,
algebraic convergence can be attained if the number of elements is increased.
In particular, when the number of basis vectors is 5, the convergence is
observed to be much steeper than when it is, say, 3 or 4. However, when the
number of basis vectors is greater than 5, no improvement in the accuracy of
the results can be achieved. This same outcome occurs consistently with the
three distributions, which means that when a Legendre quadrature rule with 10
points per element and 5 basis vectors are used in the simulations, the
maximum accuracy allowed by the machine and the ME-FSC method is ultimately
reached. Plots for the convergence of global errors using ME-gPC (with no
adaptation) are also provided in the same figure. From these, it is observed
that the convergence rate is slower than ME-FSC and that for the case of
uniform and beta distributions (distributions with compact support), the
global errors are capable of reaching the smallest value achievable by the
machine. This, however, is not the case for the gamma distribution (a
distribution with non-compact support), from which it can be observed that the
variance error obtained with ME-FSC is more than 3 orders of magnitude more
accurate. This shows the great benefit of using ME-FSC versus ME-gPC to solve
stochastic problems that have strong nonlinear dependencies over the
probability space.

(a) Mean error for $\mu\sim\mathrm{Uniform}$

(b) Variance error for $\mu\sim\mathrm{Uniform}$

(c) Mean error for $\mu\sim\mathrm{Beta}$

(d) Variance error for $\mu\sim\mathrm{Beta}$

(e) Mean error for $\mu\sim\mathrm{Gamma}$

(f) Variance error for $\mu\sim\mathrm{Gamma}$

Figure 5: _Problem 1_ — Global error of $\mathbf{E}[u]$ and $\mathrm{Var}[u]$
for different $(h,p)$-discretization levels of RFS

###  7.2 Problem 2: A linear system governed by a 3rd-order stochastic ODE

We next consider the problem of a linear mechanical system governed by a
third-order stochastic ODE. The governing differential equation for this
system is defined by

$\partial_{t}^{3}u+\tfrac{1}{2}\partial_{t}^{2}u+k\,\partial_{t}u+u=0,$

where $k:\Xi\to\mathbb{R}$ is a stochastic mechanical parameter given by
$k(\xi)=\xi$, and $u:\mathfrak{T}\times\Xi\to\mathbb{R}$ denotes the
displacement of the system with
$\partial_{t}u,\partial_{t}^{2}u,\partial_{t}^{3}u$ representing the velocity,
acceleration and jerk of the system, respectively. The initial conditions of
the system are deterministic and given by: $u(0,\cdot\,)\equiv 1$ m,
$\partial_{t}u(0,\cdot\,)\equiv-1$ m/s, and
$\partial_{t}^{2}{u}(0,\cdot\,)\equiv 2$ m/s2.

As with Problem 1, three different probability distributions are explored for
$\xi$. Namely, a _uniform distribution_ defined by
$\mathrm{Uniform}\sim\xi\in\Xi=[a,b]$, a _beta distribution_ defined by
$\mathrm{Beta}(\alpha,\beta)\sim\xi\in\Xi=[a,b]$, and a _normal distribution_
defined by $\mathrm{Normal}(\mu,\sigma^{2})\sim\xi\in\Xi=\mathbb{R}$, whence
$(a,b)=(2,3)$ N/m, $(\alpha,\beta)=(2,5)$ and $(\mu,\sigma)=(2.5,0.125)$ N/m.
However, since the normal distribution has non-compact support, the
simulations are run with _truncated normal distribution_ ,
$\mathrm{Normal}(\mu,\sigma^{2},c,d)\sim\xi\in\Xi=[c,d]$, where
$(c,d)=(1.4,3.6)$ N/m. The exact solution, nonetheless, is obtained with the
normal distribution.

Fig. 6 shows the evolution of the mean and variance of the system’s
displacement for the particular case where the probability measure is uniform
and 8 elements and 8 basis vectors are used. Like before, the results obtained
with ME-FSC are again indistinguishable from the exact response.

(a) Mean

(b) Variance

Figure 6: _Problem 2_ — Evolution of $\mathbf{E}[u]$ and $\mathrm{Var}[u]$ for
the case when the $(h,p)$-discretization level of RFS is $(P,E)=(7,8)$ and
$\mu\sim\mathrm{Uniform}$

In Fig. 7 we present the local errors in mean and variance of the system’s
displacement using different numbers of elements and basis vectors. For
brevity, the results are only presented for $\mu\sim\mathrm{Uniform}$. Once
again, it is observed that the accuracy of the results improves as the number
of elements increases, but it necessarily does not as the number of basis
vectors increases. This is exemplified in Figs. 7(c) to 7(f) from where it is
deduced that the results do not improve if the number of basis vectors is
increased from 6 to 8. This once again is due to the quadrature rule used to
compute the inner products (i.e. 10 Legendre quadrature points per element)
and the limited precision of the machine.

(a) Mean error for $\mathscr{Z}^{[3]}$

(b) Variance error for $\mathscr{Z}^{[3]}$

(c) Mean error for $\mathscr{Z}^{[5]}$

(d) Variance error for $\mathscr{Z}^{[5]}$

(e) Mean error for $\mathscr{Z}^{[7]}$

(f) Variance error for $\mathscr{Z}^{[7]}$

Figure 7: _Problem 2_ — Local error evolution of $\mathbf{E}[u]$ and
$\mathrm{Var}[u]$ for different $(h,p)$-discretization levels of RFS and for
$\mu\sim\mathrm{Uniform}$

Finally, Fig. 8 depicts the global errors in mean and variance of the system’s
displacement. The results are presented as a function of the number of
elements and the number of basis vectors used for each of the probability
distributions explored for $\xi$. Again, algebraic convergence is overall
attainable for ME-FSC if the number of elements is increased. However, as
shown, using more than 4 basis vectors in the simulations does not help
improve the overall accuracy of the results, except of course when 4 elements
and 5 basis vectors are used for the mean response. This figure therefore
suggests that in some situations, it might be beneficial to refine the
partition of the random domain (instead of making the random function space
bigger) to obtain more accurate results. As for the simulations with ME-gPC,
the convergence trend shown in Problem 1 is also showcased in this problem
through the same type of plot. For those distributions with compact support
(i.e. uniform and beta distributions), the global errors can eventually reach
the smallest value attainable by the machine. In fact, this same conclusion
applies to the distribution with non-compact support (i.e. normal
distribution), which, as observed, is also capable of achieving global errors
as small as those produced by ME-FSC.

(a) Mean error for $\mu\sim\mathrm{Uniform}$

(b) Variance error for $\mu\sim\mathrm{Uniform}$

(c) Mean error for $\mu\sim\mathrm{Beta}$

(d) Variance error for $\mu\sim\mathrm{Beta}$

(e) Mean error for $\mu\sim\mathrm{Normal}$

(f) Variance error for $\mu\sim\mathrm{Normal}$

Figure 8: _Problem 2_ — Global error of $\mathbf{E}[u]$ and $\mathrm{Var}[u]$
for different $(h,p)$-discretization levels of RFS

###  7.3 Problem 3: A nonlinear system governed by a 2nd-order stochastic ODE
(the Van-der-Pol oscillator)

In this section, we investigate the nonlinear behavior of a single-degree-of-
freedom system with Van-der-Pol damping. The governing differential equation
for this system is given by

$m\ddot{u}-(1-\rho u^{2})c\dot{u}+ku=0,$

where $m=100$ kg is the mass of the system, $\rho=150$ m-2 is the contributing
factor to the nonlinearity of the system, $c:\Xi\to\mathbb{R}^{+}$ is a
coefficient representing the strength of the damping in the system, and
$k:\Xi\to\mathbb{R}^{+}$ is the stiffness of the system. The functions $c$ and
$k$ are assumed to be given by $c(\xi)=\xi^{1}$ and $k(\xi)=\xi^{2}$. The
initial conditions of the system are: $u(0,\cdot\,)\equiv 0.20$ m and
$\dot{u}(0,\cdot\,)\equiv 0.30$ m/s. Notice here that
$u:\mathfrak{T}\times\Xi\to\mathbb{R}$ denotes the displacement of the system,
and that $\dot{u}:=\partial_{t}u$ and $\ddot{u}:=\partial_{t}^{2}u$ represent,
respectively, the velocity and acceleration of the system.

For concreteness, we take $\xi^{1}$ to be _uniformly distributed_ in
$\bar{\Xi}_{1}=[150,450]$ kg/s, and $\xi^{2}$ to be _beta distributed_ with
parameters $(\alpha,\beta)=(2,5)$ in $\bar{\Xi}_{2}=[340,460]$ N/m. Hence, the
random space is two-dimensional and defined by
$\Xi=\bar{\Xi}_{1}\times\bar{\Xi}_{2}$ with
$\mu\sim\mathrm{Uniform}\otimes\mathrm{Beta}$.

Fig. 9 depicts the evolution of the mean and variance of the system’s
displacement using ME-FSC and a Monte Carlo simulation with one million
realizations. The reason why a Monte Carlo simulation is used this time as the
reference solution is that a closed-form solution for $u$ is not available.
One drawback of using Monte Carlo as the reference solution is that the degree
of accuracy of ME-FSC will not be able to be measured adequately, since the
solution obtained with ME-FSC can be far more accurate than the one obtained
with Monte Carlo. Yet, from this figure, we learn that when 64 elements and 5
basis vectors are employed to run the simulation, the ME-FSC method is capable
of reproducing the Monte Carlo solution quite well. To compare the level of
accuracy of ME-FSC with respect to Monte Carlo, Fig. 10 presents the local
errors in mean and variance of the system’s displacement. In general, good
convergence can be observed when the number of elements and the number of
basis vectors are both increased. This observation can be better verified
using Fig. 11 which plots the global errors in mean and variance of the
system’s displacement. As expected, the results only improve if the number of
basis vectors is increased up to a certain number, which in this case happens
to be 5. However, contrary to what we observed in Problems 1 and 2, increasing
the number of elements does not necessarily improve the accuracy of the
results, as can be confirmed when $P+1$ is set to 4. This suggests that due to
the nonlinear nature of the system’s ODE, higher basis vectors can play a
major role in the description of the system’s state over time.

(a) Mean

(b) Variance

Figure 9: _Problem 3_ — Evolution of $\mathbf{E}[u]$ and $\mathrm{Var}[u]$ for
the case when the $(h,p)$-discretization level of RFS is $(P,E)=(4,64)$ and
$\mu\sim\mathrm{Uniform}\otimes\mathrm{Beta}$

(a) Mean error for $\mathscr{Z}^{[2]}$

(b) Variance error for $\mathscr{Z}^{[2]}$

(c) Mean error for $\mathscr{Z}^{[3]}$

(d) Variance error for $\mathscr{Z}^{[3]}$

(e) Mean error for $\mathscr{Z}^{[4]}$

(f) Variance error for $\mathscr{Z}^{[4]}$

Figure 10: _Problem 3_ — Local error evolution of $\mathbf{E}[u]$ and
$\mathrm{Var}[u]$ for different $(h,p)$-discretization levels of RFS and for
$\mu\sim\mathrm{Uniform}\otimes\mathrm{Beta}$

(a) Mean error

(b) Variance error

Figure 11: _Problem 3_ — Global error of $\mathbf{E}[u]$ and $\mathrm{Var}[u]$
for different $(h,p)$-discretization levels of RFS and for
$\mu\sim\mathrm{Uniform}\otimes\mathrm{Beta}$

###  7.4 Problem 4: A nonlinear system governed by a system of 1st-order
stochastic ODEs (the Kraichnan-Orszag three-mode problem)

In this last problem, we explore the Kraichnan-Orszag three-mode problem [40]
to test the performance of ME-FSC more thoroughly. This problem is
particularly challenging for methods based on the spectral approach because
the solution is known to be discontinuous over the random domain. It has been
used as a benchmark problem in various works (e.g. [6, 2, 3]), and this is the
reason why we opt to study it in this work as well.

The three-mode problem considered herein is the same as that presented in [6]
(Page 635, Section 4.3.4). The system’s governing differential equation is
defined as

$\displaystyle\dot{u}_{1}$ $\displaystyle=u_{1}u_{3},$
$\displaystyle\dot{u}_{2}$ $\displaystyle=-u_{2}u_{3},$
$\displaystyle\dot{u}_{3}$ $\displaystyle=-u_{1}^{2}+u_{2}^{2},$

where $u_{1},u_{2},u_{3}:\mathfrak{T}\times\Xi\to\mathbb{R}$ represent the
three modes of the system, and
$\dot{u}_{1}:=\partial_{t}u_{1},\dot{u}_{2}:=\partial_{t}u_{2},\dot{u}_{3}:=\partial_{t}u_{3}$
are the corresponding velocities. In this problem, we take the initial
conditions of the system to be stochastic and given by:
$u_{1}(0,\xi)=\xi^{1}$, $u_{2}(0,\xi)=\xi^{2}$, and $u_{3}(0,\xi)=\xi^{3}$.

Two probability distributions are investigated for
$\xi=(\xi^{1},\xi^{2},\xi^{3})$. The first one is a _uniform distribution_
defined by $\mathrm{Uniform}^{\otimes 3}\sim\xi\in\Xi=[a,b]^{3}$, and the
second one is a _beta distribution_ defined by
$\mathrm{Beta}(\alpha,\beta)^{\otimes 3}\sim\xi\in\Xi=[a,b]^{3}$, from where
$(a,b)=(-1,1)$ and $(\alpha,\beta)=(2,5)$. The random space is thus three-
dimensional.

In Figs. 12 and 13 we depict the evolution of the mean and variance of $u_{1}$
and $u_{3}$ for the two distributions chosen for $\xi$. The results are
obtained for the case of using 512 elements and 7 basis vectors per element,
and then compared against a Monte Carlo simulation with one million
realizations. As observed, the ME-FSC solution is in good agreement with the
Monte Carlo solution for all cases considered, except for the mean of $u_{1}$
and $u_{3}$ under uniform distributions (Figs. 12(a) and 12(c)). It is worth
commenting that the exact mean of $u_{1}$, $u_{2}$, and $u_{3}$ are known to
be identically equal to zero when the probability measure is uniform. This
demonstrates that the ME-FSC solution can be far more accurate than a Monte
Carlo simulation with one million realizations.

(a) Mean for $u_{1}$

(b) Variance for $u_{1}$

(c) Mean for $u_{3}$

(d) Variance for $u_{3}$

Figure 12: _Problem 4_ — Evolution of $\mathbf{E}[u_{1}]$,
$\mathrm{Var}[u_{1}]$, $\mathbf{E}[u_{3}]$ and $\mathrm{Var}[u_{3}]$ for the
case when the $(h,p)$-discretization level of RFS is $(P,E)=(6,512)$ and
$\mu\sim\mathrm{Uniform}^{\otimes 3}$

(a) Mean for $u_{1}$

(b) Variance for $u_{1}$

(c) Mean for $u_{3}$

(d) Variance for $u_{3}$

Figure 13: _Problem 4_ — Evolution of $\mathbf{E}[u_{1}]$,
$\mathrm{Var}[u_{1}]$, $\mathbf{E}[u_{3}]$ and $\mathrm{Var}[u_{3}]$ for the
case when the $(h,p)$-discretization level of RFS is $(P,E)=(6,512)$ and
$\mu\sim\mathrm{Beta}^{\otimes 3}$

Moreover, Figs. 14 to 16 present the local and global errors in mean and
variance of $u_{1}$ and $u_{3}$ for the two distributions chosen for $\xi$.
From these plots, it is observed that increasing the number of elements helps
improve the accuracy of the ME-FSC results. However, because the Monte Carlo
solution is not exact, the errors tend to stagnate around $10^{-3}$ and
$10^{-4}$ when the ME-FSC results are compared to Monte Carlo. This is the
reason why, in the case of the mean, the accuracy of the ME-FSC results does
not improve as the number of elements increases from 8 to 512; but, in the
case of the variance, they do improve because the errors obtained with ME-FSC
are above $10^{-4}$. Therefore, to achieve comparable solution accuracy to
Monte Carlo, around 512 elements are needed for ME-FSC.

(a) Mean error for $u_{1}$

(b) Variance error for $u_{1}$

(c) Mean error for $u_{3}$

(d) Variance error for $u_{3}$

Figure 14: _Problem 4_ — Local error evolution of $\mathbf{E}[u_{1}]$,
$\mathrm{Var}[u_{1}]$, $\mathbf{E}[u_{3}]$ and $\mathrm{Var}[u_{3}]$ for
different $(h,p)$-discretization levels of RFS and for
$\mu\sim\mathrm{Uniform}^{\otimes 3}$

(a) Mean error for $u_{1}$

(b) Variance error for $u_{1}$

(c) Mean error for $u_{3}$

(d) Variance error for $u_{3}$

Figure 15: _Problem 4_ — Local error evolution of $\mathbf{E}[u_{1}]$,
$\mathrm{Var}[u_{1}]$, $\mathbf{E}[u_{3}]$ and $\mathrm{Var}[u_{3}]$ for
different $(h,p)$-discretization levels of RFS and for
$\mu\sim\mathrm{Beta}^{\otimes 3}$

(a) Mean error for $u_{1}$

(b) Variance error for $u_{1}$

(c) Mean error for $u_{3}$

(d) Variance error for $u_{3}$

Figure 16: _Problem 4_ — Global error of $\mathbf{E}[u_{1}]$,
$\mathrm{Var}[u_{1}]$, $\mathbf{E}[u_{3}]$ and $\mathrm{Var}[u_{3}]$ for
different $(h,p)$-discretization levels of RFS

For this problem, the computational cost associated with a simulation with 512
elements and 7 basis vectors per element was 82 seconds, whereas the
computational cost associated with a Monte Carlo simulation with one million
realizations was 97 seconds. This outcome indicates that simulations conducted
with ME-FSC were at least 15% faster than those conducted with Monte Carlo.
For higher-dimensional probability spaces, it would be necessary to implement
a different quadrature rule that is not sensitive to the curse-of-
dimensionality issue in order to speed up the computation of the inner
products and reduce the computational cost associated with ME-FSC.

##  8 Conclusion

In this paper we have presented an extension of the FSC method [35, 36],
called the _multi-element flow-driven spectral chaos_ (ME-FSC) method, in
order to deal with discontinuities and long-time integration of stochastic
dynamical systems more efficiently. In ME-FSC, the stochastic problem is
divided into multiple sub-problems so that every sub-problem can be solved
individually using the FSC method. In a subsequent step, the results are
aggregated and the probability moments of interest are computed using the law
of total probability. An underlying characteristic of the FSC method makes it
possible for the probability information to be transferred exactly from one
random function space to another. Furthermore, the postulated ME-FSC scheme is
capable of reducing the error propagation by several orders of magnitude.

Four representative problems were presented in this paper to show the
effectiveness of the ME-FSC approach. The first two problems dealt with
systems governed by a linear stochastic ODE. They were selected because an
exact solution was available for each of these problems. The third problem
dealt with a system governed by a nonlinear stochastic ODE (the Van-der-Pol
oscillator). This problem was utilized to study the effectiveness of ME-FSC in
the resolution of nonlinear problems. The last problem dealt with a nonlinear
system of stochastic ODEs (the Kraichnan-Orszag three-mode problem). It was
employed to test the performance of ME-FSC when stochastic discontinuities are
present in the system. Based on our findings, we conclude that the ME-FSC
method is capable of producing accurate solutions as compared to the exact
solution (when available). For problems with no closed-form solutions, the ME-
FSC method is capable of reproducing a well-resolved Monte Carlo simulation
with one million realizations at a fraction of the computational cost required
to do so.

The ME-FSC method is particularly useful when dealing with large stochastic
dynamical systems, since it allows the analyst to decompose the problem into
several smaller sub-problems and solve them separately. In the future, the ME-
FSC method can be augmented with adaptivity criteria in order to identify
those regions in the probability space where the stochasticity of the system
is evolving faster. Moreover, if each sub-problem were to be solved
simultaneously on multiple CPUs, the computation time needed to solve each
problem could be reduced considerably.

## Acknowledgements

The first author gratefully acknowledges the scholarship granted by
Colciencias and Atlántico Department (Colombia) under Call 673 to pursue a
Ph.D. degree in Structural Engineering at Purdue University, West Lafayette,
IN. The authors also gratefully acknowledge the support of the National
Science Foundation (CNS-1136075, DMS-1555072, DMS-1736364, DMS-2053746, and
DMS-2134209), Brookhaven National Laboratory Subcontract 382247, and U.S.
Department of Energy (DOE) Office of Science Advanced Scientific Computing
Research program DE-SC0021142).

## Credit author statement

H.E. conceived the mathematical models, implemented the methods, designed the
numerical experiments, interpreted the results, and wrote the paper. A.P. and
G.L. supported the study, edited, and reviewed the final manuscript. All the
authors gave their final approval for publication.

## References

  * [1] D. Xiu, G. E. Karniadakis, The Wiener-Askey polynomial chaos for stochastic differential equations, SIAM journal on scientific computing 24 (2) (2002) 619–644.
  * [2] M. Gerritsma, J.-B. Van der Steen, P. Vos, G. Karniadakis, Time-dependent generalized polynomial chaos, Journal of Computational Physics 229 (22) (2010) 8333–8363.
  * [3] V. Heuveline, M. Schick, A hybrid generalized polynomial chaos method for stochastic dynamical systems, International Journal for Uncertainty Quantification 4 (1) (2014).
  * [4] D. M. Luchtenburg, S. L. Brunton, C. W. Rowley, Long-time uncertainty propagation using generalized polynomial chaos and flow map composition, Journal of Computational Physics 274 (2014) 783–802.
  * [5] H. C. Ozen, G. Bal, Dynamical polynomial chaos expansions and long time evolution of differential equations with random forcing, SIAM/ASA Journal on Uncertainty Quantification 4 (1) (2016) 609–635.
  * [6] X. Wan, G. E. Karniadakis, An adaptive multi-element generalized polynomial chaos method for stochastic differential equations, Journal of Computational Physics 209 (2) (2005) 617–642.
  * [7] X. Wan, G. E. Karniadakis, Long-term behavior of polynomial chaos in stochastic flow simulations, Computer methods in applied mechanics and engineering 195 (41-43) (2006) 5582–5596.
  * [8] X. Wan, G. E. Karniadakis, Error control in multi-element generalized polynomial chaos method for elliptic problems with random coefficients, Communications in Computational physics 5 (2-4) (2009) 793–820.
  * [9] X. Wan, G. E. Karniadakis, Beyond Wiener-Askey expansions: Handling arbitrary PDFs, Journal of Scientific Computing 27 (1-3) (2006) 455–464.
  * [10] J. Foo, X. Wan, G. E. Karniadakis, The multi-element probabilistic collocation method (ME-PCM): Error analysis and applications, Journal of Computational Physics 227 (22) (2008) 9572–9595.
  * [11] J. Foo, G. E. Karniadakis, Multi-element probabilistic collocation method in high dimensions, Journal of Computational Physics 229 (5) (2010) 1536–1557.
  * [12] M. Zheng, X. Wan, G. E. Karniadakis, Adaptive multi-element polynomial chaos with discrete measure: Algorithms and application to spdes, Applied Numerical Mathematics 90 (2015) 91–110.
  * [13] S. Kawai, A. Oyama, Multi-element stochastic galerkin method based on edge detection for uncertainty quantification of discontinuous responses, Journal of Verification, Validation and Uncertainty Quantification 5 (4) (2020) 041004\.
  * [14] B. V. Asokan, N. Zabaras, A stochastic variational multiscale method for diffusion in heterogeneous random media, Journal of Computational Physics 218 (2) (2006) 654–676.
  * [15] J. A. Witteveen, A. Loeven, H. Bijl, An adaptive stochastic finite elements approach based on newton–cotes quadrature in simplex elements, Computers & Fluids 38 (6) (2009) 1270–1288.
  * [16] A. Nouy, A. Clement, extended stochastic finite element method for the numerical simulation of heterogeneous materials with random material interfaces, International Journal for Numerical Methods in Engineering 83 (10) (2010) 1312–1344.
  * [17] B. Ganapathysubramanian, N. Zabaras, Sparse grid collocation schemes for stochastic natural convection problems, Journal of Computational Physics 225 (1) (2007) 652–685.
  * [18] X. Ma, N. Zabaras, An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations, Journal of Computational Physics 228 (8) (2009) 3084–3113.
  * [19] A. Bhaduri, L. Graham-Brady, An efficient adaptive sparse grid collocation method through derivative estimation, Probabilistic Engineering Mechanics 51 (2018) 11–22.
  * [20] A. Bhaduri, Y. He, M. D. Shields, L. Graham-Brady, R. M. Kirby, Stochastic collocation approach with adaptive mesh refinement for parametric uncertainty analysis, Journal of Computational Physics 371 (2018) 732–750.
  * [21] N. Agarwal, N. R. Aluru, A domain adaptive stochastic collocation approach for analysis of MEMS under uncertainties, Journal of Computational Physics 228 (20) (2009) 7662–7688.
  * [22] Y. M. Marzouk, H. N. Najm, L. A. Rahn, Stochastic spectral methods for efficient bayesian solution of inverse problems, Journal of Computational Physics 224 (2) (2007) 560–586.
  * [23] P. S. Mohan, P. B. Nair, A. J. Keane, Multi-element stochastic reduced basis methods, Computer Methods in Applied Mechanics and Engineering 197 (17-18) (2008) 1495–1506.
  * [24] E. Sarrouy, O. Dessombz, J.-J. Sinou, Piecewise polynomial chaos expansion with an application to brake squeal of a linear brake system, Journal of Sound and Vibration 332 (3) (2013) 577–594.
  * [25] J. D. Jakeman, A. Narayan, D. Xiu, Minimal multi-element stochastic collocation for uncertainty quantification of discontinuous functions, Journal of Computational Physics 242 (2013) 790–808.
  * [26] O. Le Maıtre, O. Knio, H. Najm, R. Ghanem, Uncertainty propagation using wiener–haar expansions, Journal of computational Physics 197 (1) (2004) 28–57.
  * [27] O. Le Maıtre, H. N. Najm, R. Ghanem, O. Knio, Multi-resolution analysis of wiener-type uncertainty propagation schemes, Journal of Computational Physics 197 (2) (2004) 502–531.
  * [28] C. J. Gittelson, Adaptive wavelet methods for elliptic partial differential equations with random operators, Numerische Mathematik 126 (3) (2014) 471–513.
  * [29] C. Gittelson, An adaptive stochastic galerkin method for random elliptic operators, Mathematics of Computation 82 (283) (2013) 1515–1541.
  * [30] H. Cho, D. Venturi, G. E. Karniadakis, Adaptive discontinuous galerkin method for response-excitation pdf equations, SIAM Journal on Scientific Computing 35 (4) (2013) B890–B911.
  * [31] Y. Chen, J. Jakeman, C. Gittelson, D. Xiu, Local polynomial chaos expansion for linear differential equations with high dimensional random inputs, SIAM Journal on Scientific Computing 37 (1) (2015) A79–A102.
  * [32] J. Li, P. Stinis, A unified framework for mesh refinement in random and physical space, Journal of Computational Physics 323 (2016) 243–264.
  * [33] A. Abdedou, A. Soulaimani, A non-intrusive b-splines bézier elements-based method for uncertainty propagation, Computer Methods in Applied Mechanics and Engineering 345 (2019) 774–804.
  * [34] C. Eckert, M. Beer, P. D. Spanos, A polynomial chaos method for arbitrary random inputs using b-splines, Probabilistic Engineering Mechanics 60 (2020) 103051\.
  * [35] H. Esquivel, A. Prakash, G. Lin, Flow-driven spectral chaos (FSC) method for simulating long-time dynamics of arbitrary-order non-linear stochastic dynamical systems, Journal of Computational Physics (2020) 110044.
  * [36] H. Esquivel, A. Prakash, G. Lin, Flow-driven spectral chaos (fsc) method for long-time integration of second-order stochastic dynamical systems, Journal of Computational and Applied Mathematics (2021) 113674.
  * [37] W. Cheney, D. Kincaid, Linear algebra: Theory and applications, 2nd Edition, Jones & Bartlett Learning, 2010.
  * [38] Apple Inc., Apple’s frameworks through the Xcode IDE (2021).
URL https://developer.apple.com/documentation

  * [39] Apple Inc. et al., Swift programming language (2021).
URL https://swift.org

  * [40] S. A. Orszag, L. Bissonnette, Dynamical properties of truncated Wiener-Hermite expansions, The Physics of Fluids 10 (12) (1967) 2603–2613.