Polynomial Stochastic Dynamical Indicators

A popular technique to study the dependency of a dynamical system on a set of uncertain quantities is Polynomial Chaos Expansions. The idea is to represent the state vector $\mathbf{z}$ as a truncated expansion in the orthogonal polynomials $\Psi_{i}(\mathbf{p})$ of the uncertain quantities $\mathbf{p}$:

where $\mathbf{c}_{i}(t)$ are time dependent coefficients. The $\Psi_{i}$ terms define a set of orthogonal polynomials up to degree $m$ (orthpoli, ). The orthogonality condition is formalised as follows:

where $\langle\cdot,\cdot\rangle$ is a shorthand of the inner product. As mentioned before when the $w$ is a distribution (6) defines the expectation operator associated to $w$. Because of the polynomial nature of the terms appearing in (6), it is straightforward to compute the non-zero terms. Then, given a particular weight function $w(\mathbf{p})$, one can use the following three terms recursion relation given in ttr to create stabilised univariate orthogonal polynomials:

In the case in which more than one source of uncertainty is present, it is still possible to construct orthogonal multivariate polynomials via tensor product rules (FEINBERG201546, ). Note that while the method proposed in this paper is applicable to any orthogonal polynomial constructed with (7) in all the examples in this paper Chebyshev basis functions of the second kind are used together with the associated weight function $w(\mathbf{p})$.

By substituting the approximation given by (5) in (1), one gets:

and by making use of the intrusive Galerkin method, one obtains the following:

from which the time variation of the coefficients can be derived:

The integrals at numerator of the right hand side of (10) need to be computed numerically, in the general case, while the integrals at denominator can be pre-computed analytically. Gauss quadrature rules (FEINBERG201546, ) can be used to approximate the integrals at numerator, as follows:

where $W_{j_{i}}$ and $\mathbf{p}_{j_{i}}$ are respectively the $N$ quadrature weights and abscissa points along each dimension $i$. Sparse quadrature schemes smolyak can be used to reduce the computational complexity of the numerical integrals with the increase in the number of dimensions.

The initial value of the coefficients $\mathbf{c}_{k}(t=0)$ is found by projecting the initial conditions $\mathbf{z}_{0}$:

which greatly simplifies in the case in which the initial state is deterministic (i.e., none of the components of $\mathbf{z_{0}}$ are components of $\mathbf{p}$): the only non-zero coefficient is $\mathbf{c_{0}}$, the one associated to the degree-zero polynomial of the orthogonal basis, whose value is the one of the deterministic initial condition.

Up to this point PCEs are simply a way to represent the state of the system $\mathbf{z}$ with a polynomial expansion of the parameters $\mathbf{p}$ and propagate this expansion forward in time. Thus regardless of whether $\mathbf{p}$ is an uncertain quantity with an associated probability distribution $w$ or a simple parameter defined on a parameter space $\Omega$, (41) provides a way to propagate the polynomial forward in time.

Furthermore, (12) can be applied at any time $t$ to calculate a polynomial expansion of the state variables with respect to the uncertainty variables. In this case (12) reads:

In both (12) and (13) the integral at denominator can be computed analytically, one time before the calculation of the coefficients. The integral at numerator of (13) can be solved numerically as in (11):

The polynomial expansion computed with (13) is called non-intrusive because one needs only samples of the state vector $\mathbf{z}(t,\mathbf{p})$ at time $t$ for different realisations of $\mathbf{p}$. These samples can be obtained from the direct forward integration of the equations of motion.

The use of a non-intrusive computation of the coefficients of the polynomial expansion is advantageous when the dynamical model is not directly accessible, the state vector is available through observations or, as it will be explained in section 5, if the integration of system (17) becomes problematic due to the presence of singularities or discontinuities in the uncertainty space. In this case restart mechanisms like the ones proposed in greco:2020 ; manziaas2020 and ozen can be effectively used to improve the propagation of the polynomial expansion. In this paper, however, we will not consider these restart mechanisms and we will show the use of (13) instead of (10) to compute two of the indicators.

Since the interest is to exploit the evolution of the coefficients of a polynomial expansion and not to exactly propagate a particular probability distribution, the weight $w$ and basis functions $\Psi$ can be arbitrarily chosen to make the numerical integration of (11) efficient. In the following we will consider the components of $\mathbf{p}$ to be independent and $\Omega$ to be a orthotope. Furthermore, integral (11) is performed after the change of coordinates:

with $p_{i}\in[a_{i},b_{i}]$ and $\mathbf{\xi}\in[-1,1]^{n}$ so that:

In this section we derived the expansion, intrusive or non-intrusive, of $\mathbf{z}$ in orthogonal polynomials of $\mathbf{p}$. Other forms of uncertainty quantification in the literature, like Taylor series expansions, for example, do not use orthogonal polynomials. However, in the definition of the stochastic dynamical indicators we will exploit the orthogonality of the polynomials. Thus, while, in principle, any polynomial representation of $\mathbf{z}$ is applicable, before computing the stochastic indicators one would need to transform the polynomial expansion into orthogonal basis as suggested in Iosto2022 .

Note also that the use of Taylor expansions to derive dynamical indicators was already proposed in Palau2015 . However, the approach introduced in this paper differs from the one in (Palau2015, ) in two important ways: i) in this paper we use the evolution of the coefficients of the polynomials to directly define the indicators and ii) the indicators proposed in this paper quantify the effect of uncertainty in the parameters defining the dynamic model. This later point is of particular importance because, as it will be explained in the remainder of the paper, the primary utility of the indicators proposed in this work is to study the effect of the uncertainty in the dynamic model.

Following milani Section 2.3 we now briefly recall the definition of Finite-Time Lyapunov Exponents. We start from the definition of the variational equations in the deterministic settings:

The FTLE emerges from the spectral analysis of the Cauchy–Green (CG) Strain Tensor:

where $\Phi$ is the state transition matrix of the system. From it, the definition of Finite-Time Lyapunov Exponent (shadden2005, ) is given by:

where $t_{f}$ is the time interval associated to the propagation, starting at $t_{0}$, and $\lambda_{max}$ is the maximum eigenvalue of the Cauchy–Green Strain Tensor.

A random-walk is a stochastic process that defines a path made of random steps. Steps can have random direction, random length and be taken at random times. One of the best known random-walks is Brownian motion. Brownian motion can be well described by a Weiner process $W_{t}$ with independent steps and each step taken from a normal distribution $\mathcal{N}(0,t)$ with zero mean and variance $Dt$.

where $D$ is the diffusion coefficient. In normal diffusion the exponent of the time $t$ is one, however, some stochastic processes can diffuse faster or slower (e.g. fractionated Brownian motion or Levy processes) (ALVES2016392, ). Thus in the general case one can write:

where $K$ is a constant and $\alpha$ is the diffusion exponent. In the next section we will make use of (22) to derive an indicator that relates the coefficients of the polynomial expansion to the diffusion exponent.

In this section we will develop two types of Stochastic Finite-Time Lyapunov Exponents. The first type replaces the FLTE with the statistical moments quantifying the uncertainty in the FTLE. If the dynamics depends on some uncertain quantities, the Strain Tensor in (18) is a random matrix with entries that are a function of the realisations of the uncertain quantities. Thus one could study the ensemble of matrices and derive a statistics over the realisations of the eigenvalues. An approach to derive the statistical moments of the FTLE can be found in Schomerus2002 . In Schomerus2002 the authors considered the case of a one-dimensional dynamical system driven by a random potential and built the statistical moments of the FTLE by computing the moments of the components of the matrix $\partial\mathbf{z}(t_{f})(t)/\partial\mathbf{z}_{0}$. In what follows, instead, we will use a polynomial chaos expansion of the FTLE with respect to the uncertain vector $\mathbf{p}$. By sampling the uncertain space $\Omega$, one can directly construct the PCE expansion of the FTLE $\sigma$ defined in (19):

where the coefficients $\sigma_{k}(t_{f})$ are computed by projection:

For the second type, as in the deterministic settings, we start from the hypervolume $d\mathbf{z}^{T}d\mathbf{z}$ and compute the time evolution of its expectation $\mathbb{E}(d\mathbf{z}^{T}d\mathbf{z})$.

We now derive an equivalent definition of variational equations (17) but in the coefficients of the PCE expansion of $d\mathbf{z}$.

By using the PCE decomposition, the second term in Eq. (35) leads to:

while the first term of Eq. (35) leads to:

By putting Eqs. (36) and (37) back into Eq. (35) one gets:

and, by making use of the orthogonality condition, one arrives at the following result:

As discussed for the deterministic formulation in gaw2007 , in order to compute the variation of the coefficients $\mathbf{c}_{i}$, it is possible to propagate a regularly spaced grid of tracers with the same dimension as the phase space. In fact, the spectral harmonics of the generalised State Transition Matrix appearing in Eq. (17) consist of partial derivatives which can be computed via central differencing of neighboring tracers, making use of the following second-order approximation:

with $\Delta\mathbf{z}_{j}=[0,\dots,0,\Delta z_{j},0,\dots,0]$. This methodology greatly reduces the computational cost associated to the generalisation of the variational equations, as it is for the deterministic case. While the accuracy of the computation of the CG tensor degrades with this approach, the authors in shadden2005 points out that: “finite differencing may unveil Lagrange Coherent Structures more reliably than obtaining derivatives of the flow analytically”.

From (28), one can now introduce the Cauchy-Green Tensor $\Delta^{c}_{ii}$ of the coefficients $\mathbf{c}_{i}$:

Indicators in Definition 1 will be called SFTLE1 in the remainder of this paper while indicators in Definition 2 will be called SFTLE2. Indicators SFTLE1 give the probability distribution of the FTLE in (19) as a function of the distribution of the uncertain parameter vector $\mathbf{p}$. Indicators SFTLE2, instead, give a measure of the divergence of the coefficients of the polynomial model of the distribution of the solution $\mathbf{z}(t,\mathbf{p})$ as a function of a variation of the initial condition $\mathbf{z}_{0}$. It should be noted how the eigenvectors associated to the parameter-dependent Cauchy–Green strain tensor are also characterised by a probability distribution. This implies that the direction of maximum strain is not deterministic, and there may be configurations in which there is an abrupt change of the maximum strain direction for different realisations of the uncertain parameter.

In order to derive the third indicator we start from the idea, introduced in section 3.3, that in a generic random-walk process, the expected value of the square of the displacement is proportional to $Kt^{\alpha}$. In the univariate case, by using Eq. (26) and exploiting the orthogonality of the basis functions, one can write the expected value of the square displacement as:

with $s_{i}=\langle\psi_{i},\psi_{i}\rangle$. One can now equate $\kappa_{2}$ to $Kt^{\alpha}$ to obtain:

The left hand side is the variance of $z$ at time $t$, which, for $\alpha=1$, is consistent with the fact that for a one-dimensional Brownian motion the second statistical moment of the Mean Square Displacement (MSD) is $2Dt+z_{0}^{2}$, with $2D=K$ the diffusion coefficient, and the MSD is equal to the second cumulant of the Gaussian distribution characterising the Brownian motion. This suggests that by looking at the variation of the coefficients of the polynomial, one can study the dynamical character of a system. Since the coefficients are subject to the same dynamic equations, see (41), they reflect the same evolution of the state. The evolution of the coefficients can be derived in other ways, for example via an algebra on the space of the polynomials greco:2020 ; Palau2015 . As long as the state can be expressed as an expansion in orthogonal polynomials one can derive Eq. (43).

Note that both intrusive and non-intrusive propagation methods can be used to compute the coefficients of the polynomials at time $t$. However, in all the examples in this paper the pseudo-diffusion exponent will be computed with a non-intrusive computation of the coefficients.

The motion of a periodically-perturbed pendulum can be written as in Palau2015 :

or as an equivalent system of first order differential equations:

with $p=a$ an uncertain parameter. One can then write the Jacobian of system (52):

The uncertain parameter $a$ is defined over the interval $a\in[2.5-0.25,\;2.5+0.25]$, with known or unknown distribution, dynamics (51) becomes uncertain and its evolution depends on the realisations of $a$. Thus we expanded the state variables in Chebyshev polynomials of parameter $a$, up to degree 4, and used the definition of the three indicators SFTLE1, SFTLE2 and $\tilde{\alpha}$ to study the evolution of the system.

All differential equations, were propagated forward in time for $t_{f}=10$, with an explicit adaptive Runge-Kutta method of order 4/5 with absolute tolerance and relative tolerance respectively $10^{-10}$ and $10^{-9}$. The three indicators were computed over a uniform grid of $200\times 200$ initial conditions over the domain $x\in[-3,3]$, $v_{x}\in[-3,3]$. The finite increment for he calculation of both the FTLE and SFTLE is $\Delta z_{j}=1\cdot 10^{-7}$.

Figure 1(a) shows the deterministic FTLE for $a=2.5$ while Figure 1(b) shows $\alpha_{1}^{1}$ for $a$ uncertain. Although the magnitude of the two indicators is slightly different, they present the same structures, as to be expected given that $\alpha_{1}^{1}$ is an average value over the realisations of $a$. Figure 1(c) represents the variance of the FTLE due to the uncertainty in $a$ and Figure 1(d) the skewness. Because the $\sin()$ is an odd function, the mapping $(x,v_{x})\mapsto(-x,-v_{x})$ is a symmetry of (52) and, because of this, the results given in Figure 1 are characterised by a central symmetry with respect to the origin. Note however that Figure 1(d) clearly shows that the realisations of the state vector at time $t_{f}$ are positively or negatively skewed depending on the initial conditions. Thus SFTLE1 provides different information on the distribution of the FTLE depending on the order of the indicator.

Figures 2 show the SFTLE2 from order 1 to 3. In this case all three indicators show the same structures but with very different ranges. To be noted that as the order increases the regions where the indicators are negative become more negative. This implies that the higher the coefficient $c$ the more two expansions starting from neighbouring initial conditions tend to behave similarly.

Figures 3 show the pseudo diffusion exponent field together with the probability of the trajectories in the ensemble to remain within a distance $\epsilon=0.1$ from the mean at time $t_{f}$ (see Eq. (3)) and the skewness of the ensemble of trajectories induced by multiple realisations of the uncertain parameter $a$. The skewness is computed only for the state component $x$. For multivariate problems one would need to compute the skewness vector KOLLO20082328 and then reduce it to a scalar indicator. This computation will be addressed in future work. Figure 3b shows the $\log 10$ of Figure 3a. Also this indicator identifies the same structures as SFTLE1 and 2 and the associated skewness is consistent with Figure 1d. Figure 3c provides some additional information. First it is interesting to note that it is the negative image of Figure 3a which is consistent with the idea that $\tilde{\alpha}$ provides a measure of the diffusion of the trajectories. Then 3c highlights how some sets of initial conditions, yellow regions, are weakly sensitive to the uncertainty in $a$.

Finally Figure 4 shows two notable trajectories, one for initial conditions $\mathbf{z}_{0}=[0.889447,\,-0.19598]$ and the other for $\mathbf{z}_{0}=[1.67337,1.19095]$, which correspond respectively to low and high values of $\tilde{\alpha}$. In this case 10 trajectories were propagated for 10 random realisations of $a$.

The double gyre model consists of a pair of counter-rotating gyres, with a time-periodic perturbation. The system is modelled as a first order system of differential equations, given by:

The functions and the coefficients appearing in the dynamics are given by:

The Jacobian of the velocity field is given by:

We generalise results from previous works (e.g., faraz2012 ), by considering the uncertain parameter $p=\eta$ to be an uncertain parameter defined over the interval $\eta\in[0.1-0.01,0.1+0.01]$.

As in the previous example we expand the state variables in Chebyshev polynomials of degree 4. Differential equations are propagated with the same adaptive Runge-Kutta integrator with the same absolute and relative tolerances. The propagation is performed for a fixed integration time $t_{f}=20$. A uniform grid of initial conditions has a size of $200\times 200$, in the domains $x\in[0,2]$, $y\in[0,1]$. The finite increment for the calculation of both the FTLE and SFTLE is $\Delta z_{j}=1\cdot 10^{-7}$.

Figures 5a and b compare the deterministic FTLE with SFTLE1. Also in this case $\alpha_{1}^{1}$ shows the same structures as the FTLE. It is interesting to note in Figure 5c, how the location of the ridges of $\alpha_{1}^{2}$, are located near the ridges of $\alpha_{1}^{1}$. This implies that for chaotic initial conditions the set of trajectories behaves qualitatively differently with different realisations of the uncertain parameter. This emerges also from 5d where the skewness of FLTE is positive or negative depending on the initial conditions. Similar consideration can be derived from Figure 6 where the SFTLE2 are represented and from Figure 8 where $\tilde{\alpha}$ is represented together with the expectation for a threshold of $\epsilon=0.25$ and the skewness of the x component of the ensemble of trajectories.Ten trajectories corresponding to ten realisation of $\eta$ are represented in Figure 7 for two initial conditions $x_{0}=1.37688$, $y_{0}=0.73869$ and $x_{0}=1.45729$, $y_{0}=0.44221$ corresponding respectively to high and low values of $\tilde{\alpha}$. Note in Figure 7a the bifurcation of the ensemble into two different groups of trajectories.

The Circular Restricted Three-Body Problem (CR3BP) is arguably one of the most studied problems in celestial mechanics. In this section we will consider the planar case with an uncertain mass parameter. The planar circular restricted three-body problem szeb is governed by:

where $J(x,y)$ is given by:

and $\mu$, the mass parameter of the system, is a function of the masses of the primaries. With this formulation, the reference frame is uniformly rotating and the primaries’ position, in such frame, is constant in time. We can again re-write the system as a first order system of differential equations:

with $v_{x}=\dot{x}$ and $v_{y}=\dot{y}$ and uncertain parameter $p=\mu$. The partial derivatives of $J$, appearing in (59), are given by:

The Jacobian of the velocity field, associated to the first-order formulation of the dynamics, is:

in which the second order derivatives of $J$ are not expanded for brevity. The energy is then defined as,

and is a constant of motion for the CR3BP. In order to reduce the dimensionality of the problem from four to two, the initial conditions are defined as:

where the value of $v_{y}$ is computed, from a given condition $(x_{i},v_{xi})$, making use of the conservation of energy:

For the results in this paper, the energy level has been set to $E_{0}=E(L_{1})+0.03715$, where $E(L_{1})=E(L_{1}^{x},0,0,0)$ is the potential energy at the first Lagrange point, $L_{1}^{x}$ being given wakker by:

with

and

We consider two cases. In case 1, the integration is performed for two full revolutions of the primaries, or $t_{f}=2$, using an explicit adaptive Runge-Kutta 4/5 method with absolute tolerance $10^{-10}$ and relative tolerance of $10^{-8}$. As in the previous two examples, we use Chebyshev orthogonal polynomials of type 2, thus the integration abscissas and weights are optimised for these polynomials. The initial conditions grid is $200\times 200$, in the domains $x\in[-0.85,-0.125]$, $v_{x}\in[-2.0,2.0]$. The value of the mass parameter is assumed to be uncertain and within the interval reported in Table 1 case 1. The finite increment for the calculation of both the FTLE and SFTLE is $\Delta z_{j}=1\cdot 10^{-7}$. Figure 9 reports the FTLE and SFTLE1 for case 1. Polynomials are expanded to order 4 and the figure represents the first three SFTLE1.