Computing invariant sets of random differential equations using polynomial chaos

TL;DR

This paper introduces a polynomial chaos-based method to efficiently compute and visualize invariant sets of random differential equations, enabling uncertainty quantification and faster sampling in complex dynamical systems.

Contribution

It presents a novel approach using polynomial chaos to analyze invariant sets of RODEs, enhancing computational efficiency and visualization capabilities.

Findings

Efficient computation of invariant sets for RODEs using PC.

Successful application to predator-prey and Lorenz systems.

Demonstrated numerical efficiency with adaptive PC methods.

Abstract

Differential equations with random parameters have gained significant prominence in recent years due to their importance in mathematical modelling and data assimilation. In many cases, random ordinary differential equations (RODEs) are studied by using Monte-Carlo methods or by direct numerical simulation techniques using polynomial chaos (PC), i.e., by a series expansion of the random parameters in combination with forward integration. Here we take a dynamical systems viewpoint and focus on the invariant sets of differential equations such as steady states, stable/unstable manifolds, periodic orbits, and heteroclinic orbits. We employ PC to compute representations of all these different types of invariant sets for RODEs. This allows us to obtain fast sampling, geometric visualization of distributional properties of invariants sets, and uncertainty quantification of dynamical output such as periods or locations of orbits. We apply our techniques to a predator-prey model, where we compute steady states and stable/unstable manifolds. We also include several benchmarks to illustrate the numerical efficiency of adaptively chosen PC depending upon the random input. Then we employ the methods for the Lorenz system, obtaining computational PC representations of periodic orbits, stable/unstable manifolds and heteroclinic orbits.

Equations

$$\frac{\text{d}x}{\text{d}t}=\dot{x}=f(x,p),x=x(t)\in {\mathbb{R}}^m,$$

$$X(p)=\sum _{n\in {\mathbb{N}}^d}{X}_n{\varphi }_n(p).$$

$$X^N(p)=\sum _{n=0}^N{X}_n{\varphi }_n(p),$$

$$\underset{p\in [-1,1]}{\sup }\Vert X(p)-X^N(p)\Vert \le \sum _{n=N+1}^{\infty }\Vert X_n\Vert ,$$

$$\Vert \mathbb{E}\left(X(p)\right)-\mathbb{E}\left(X^N(p)\right)\Vert \text{or}\Vert \mathbb{V}\left(X(p)\right)-\mathbb{V}\left(X^N(p)\right)\Vert ,$$

$$\Vert \mathbb{E}\left(X(p)\right)-\mathbb{E}\left(X^N(p)\right)\Vert \le \sum _{n=N+1}^{\infty }\Vert X_n\Vert \left|\mathbb{E}\left({\varphi }_n(p)\right)\right|.$$

$${\int }_{\mathbb{R}}{s}^n\rho (s)\text{d}s<\infty \forall n\in \mathbb{N}.$$

$${⟨{\varphi }_{n_1},{\varphi }_{n_2}⟩}_{\rho }:={\int }_{\mathbb{R}}{\varphi }_{n_1}(s){\varphi }_{n_2}(s)\rho (s)\text{d}s=h_{n_1}{\delta }_{n_1,n_2},\forall n_1,n_2\in \mathbb{N}.$$

$${\int }_{\mathbb{R}}{g}^2(s)\rho (s)\text{d}s<\infty ,$$

$$g=\sum _{n\in \mathbb{N}}{g}_n{\varphi }_n,$$

$$\mathbb{E}\left(X(p)\right)=\sum _{n=0}^{\infty }{X}_n\mathbb{E}\left({\varphi }_n(p)\right)=\sum _{n=0}^{\infty }{X}_n⟨{\varphi }_n,1{⟩}_{\rho }=X_0,$$

$$\mathbb{V}\left(X(p)\right)=\mathbb{E}\left(X(p{)}^2\right)-\mathbb{E}{\left(X(p)\right)}^2=\sum _{m=0}^{\infty }\sum _{n=0}^{\infty }{X}_n{X}_m⟨{\varphi }_m,{\varphi }_n{⟩}_{\rho }-X_0^2=\sum _{n=1}^{\infty }{X}_n^2.$$

$${\varphi }_n(s):=\prod _{j=1}^d{\varphi }_{n_j}^{(j)}(s_j),$$

$$X_n=\frac{1}{h_n}⟨X,{\varphi }_n{⟩}_{\rho },$$

$$F(X(p),p)=0,$$

$${⟨F\left(\sum _{n_1=0}^{\infty }{X}_{n_1}{\varphi }_{n_1}(p),p\right),{\varphi }_{n_2}(p)⟩}_{\rho }=0\forall n_2\in \mathbb{N}.$$

$$X^N(p)=\sum _{n_1=0}^N{X}_{n_1}{\varphi }_{n_1}(p),$$

$$F(X(p_i),p_i)=0,\forall 1\le i\le I.$$

$$X_n=\sum _{i=1}^{I}X(p_i){\varphi }_n(p_i){\alpha }_i,$$

$$\sum _{i=1}^{I}g(p_i){\alpha }_i\approx \int g(s)\rho (s)\text{d}s,$$

$$X^N(p)=\sum _{n=0}^N{X}_n{\varphi }_n(p)$$

$$X(t,p)=\sum _{k\in \mathbb{Z}}{X}_k(p){\text{e}}^{\text{i}k\Omega (p)t},X_k(p)\in {\mathbb{R}}^m,\Omega (p)\in \mathbb{R},$$

$$\text{i}k\Omega (p)X_k(p)=f_k(X(t,p),p)\forall k\in \mathbb{Z},$$

$$\left(X(0,p)-u(p)\right)\cdot v(p)=0,$$

$$X_k(p)=\sum _{n\in \mathbb{Z}}{X}_{k,n}{\varphi }_n(p)\text{and}\Omega (p)=\sum _{n\in \mathbb{N}}{\Omega }_n{\varphi }_n(p),$$

$$X(t,p)=\sum _{k\in \mathbb{Z}}\sum _{n\in \mathbb{N}}{X}_{k,n}{\varphi }_n(p){\text{e}}^{\text{i}kt\Omega (p)}.$$

$$X_k^N(p)=\sum _{n=0}^N{X}_{k,n}{\varphi }_n(p)\text{and}\Omega (p)=\sum _{n=0}^N{\Omega }_n{\varphi }_n(p).$$

$$Q(\theta ,p)=\sum _{|k|\ge 0}{Q}_k(p){\theta }^k,\theta =\left(\begin{array}{c}{\theta }_1\\ ⋮\\ {\theta }_{m_s}\end{array}\right)\in {\mathbb{R}}^{m_s},Q_k(p)=\left(\begin{array}{c}{Q}_k^{(1)}(p)\\ ⋮\\ Q_k^{(m)}(p)\end{array}\right)\in {\mathbb{R}}^m,$$

$$Q_0(p)=\hat{X}(p)\text{and}{Q}_{e_i}(p)={\gamma }_i{v}^{(i)}(p),\forall 1\le i\le m_s,$$

$${\phi }_p(t,Q(\theta ,p))=Q({\text{e}}^{\Lambda (p)t}\theta ,p),\forall {\Vert \theta \Vert }_{\infty }\le 1,\forall t\ge 0,$$

Full text

**Keywords: invariant manifold, periodic orbit, heteroclinic orbit, Lorenz system, polynomial chaos, random differential equation. **

1. 1 Introduction

2. 2 A quick review on PC

3. 3 Computation of periodic orbits, invariant manifolds and connecting orbits

4. 3.1 Periodic orbits via Fourier series

5. 3.2 Local invariant manifolds via Taylor series and the parameterization method

6. 3.3 Heteroclinic orbits via Chebyshev series and projected boundaries

7. 4 First example: a Lotka-Volterra system

8. 4.1 Analytical results (for benchmarking)

9. 4.1.1 Deterministic framework

10. 4.1.2 Stochastic framework

11. 4.2 Numerical results

12. 4.2.1 Analysis of convergence on the steady state problem

13. 4.2.2 Computation of eigenvalues and eigenvectors

14. 4.2.3 Computation of local stable/unstable manifolds

15. 5 Second example: the Lorenz system

16. 5.1 Computation of periodic orbits

17. 5.2 Computation of heteroclinic orbits

18. 6 Conclusions & Outlook

19. 7 Appendix

20. 7.1 Linearization formulas for products

21. 7.2 Treatment of higher order terms

22. 7.3 Faster computations

1 Introduction

In this work, we study random nonlinear dynamical systems. More precisely, we focus on nonlinear random ordinary differential equations (RODEs) of the form

[TABLE]

where $f:\mathbb{R}^{m}\times\mathbb{R}^{d}\to\mathbb{R}^{m}$ is a smooth vector field and $p=p(\omega)\in\mathbb{R}^{d}$ denotes random parameters with given distributions. From the viewpoint of applications, it is frequently natural to assume that the parameters $p$ are only known from measurements, which naturally carry an associated probability distribution. Then the challenge is to quantify the uncertainty in the “output” of the RODE (1) based upon the random input. Yet, we are still far away to fully understand the nonlinear dynamics of such RODEs. In classical uncertainty quantification problems [12, 19], one is often interested [21, 33] in the moments as an output of the solution $\mathbb{E}[x(t)^{k}]$ for $k\geq 1$, where $\mathbb{E}$ denotes the expectation. Here we take a dynamical systems perspective focusing on the invariant sets (e.g. steady states, periodic orbits, invariant manifolds, connecting orbits, etc.) of (1), and especially in understanding their dependence on the noisy parameters $p$. In this paper, we study the invariant sets from a numerical viewpoint. However, the framework that we develop is well suited to the usage of rigorous numerics [30], and in particular of a-posteriori validation techniques [32], which could applied to obtain rigorous results about these stochastic invariant sets. This idea will be presented in a forthcoming work.

Let us start by mentioning that there exists many well developed techniques to numerically study invariant sets of deterministic ODEs. Therefore, a natural way of studying the invariant sets of (1) would be to use a Monte-Carlo type approach: consider a large sample of values $p_{i}$ taken according to the distribution of $p$, and for each $i$ study the invariant sets of the deterministic ODE $\dot{x}=f(x,p_{i})$. However, this approach is known to be very costly, because it requires a large sample to accurately represent the statistics of the invariant sets [10]. In this work we make use of a different technique, namely polynomial chaos (PC) expansions [35, 11], to accurately compute invariant sets of (1). Roughly speaking, we view each invariant set of (1) as a curve parametrized by $p$ (or as a manifold if $p$ is more than one-dimensional), and compute such parameterization explicitly via a PC expansion. This can be thought of as a parameter continuation in $p$, but in an astute way that allows us to obtain not only the geometrical object (i.e. the curve/manifold of invariant sets), but also statistical properties of this object (e.g. mean position, variance, …). Furthermore, our techniques naturally extend to numerical continuation algorithms [16, 17] if the probability distributions of $p$ contain further parameters, which is a direction we will pursue in future work.

Let $X=X(p)$ be an invariant set of (1), i.e., trajectories starting inside $X$ remain in $X$ for all $t\in\mathbb{R}$. In this work, we focus on the cases where $X(p)$ represents a steady state, a stable/unstable manifold, a periodic orbit or a heteroclinic orbit, which are among the most important objects in nonlinear dynamics [18, 13]. One key observation is that PC expansions can be used to unify the computational framework for invariant sets and that they provide a natural deterministic structure to view the dynamics of nonlinear random ODEs. Our goal is to find a series expansion of $X$ as a function of $p$:

[TABLE]

Note that the coefficients $X_{n}$ have to be found as solutions of suitable nonlinear problems to correctly represent the different invariant sets; we develop the details for each type of invariant set in this work. In practice, the choice of the expansion basis $\left(\phi_{n}\right)_{n\in\mathbb{N}^{d}}$ is of course also crucial, as it determines the quality of the approximation, or more precisely the number $N$ of coefficients required to obtain a good enough approximation. For the rest of this discussion, we assume for simplicity that $p$ is one-dimensional.

If $f$ is analytic with respect to $p$, then we can expect $X$ to also be analytic as a function of $p$, at least around values of $p$ that are not bifurcation values. Therefore, in many cases one could think about writing $X$ as a Taylor series by taking $\phi_{n}(p)=p^{n}$. However, it is well known in approximation theory that faster convergence can be achieved by considering instead Chebyshev series or Legendre series (i.e. taking $\phi_{n}(p)=T_{n}(p)$ or $\phi_{n}(p)=L_{n}(p)$, where $T_{n}$ and $L_{n}$ respectively denote Cheybshev and Legendre orthogonal polynomials), mainly because these expansions are less sensitive to potential poles of $X$ in the complex plane; see e.g. [29]. Chebyshev or Legendre expansions also have the advantage of being convergent even when $X$ is only of class ${\mathcal{C}}^{k}$, where the coefficients $X_{n}$ decay at an algebraic rate, rather than geometric in the analytic case. If $p$ is deterministic, then looking at the decay rates of the coefficients $X_{n}$ is a relevant benchmark, because it is related to the error in ${\mathcal{C}}^{0}$-norm. For instance, if $X^{N}$ is the truncated series given by

[TABLE]

then for a Taylor expansion, a Chebyshev expansion or a Legendre expansion alike one has

[TABLE]

because for each of these choices one has $\sup_{p\in[-1,1]}\left|\phi_{n}(p)\right|\leq 1$. However, when $p=p(\omega)$ is a random variable, one is more interested in controlling different quantities such as

[TABLE]

where $\mathbb{V}$ denotes the variance. To minimize the error in the moments, the choice of the expansion is critical not only because it influences the decay of the coefficients, but also because the $\phi_{n}$ themselves appear in the error term. For instance, one has

[TABLE]

Therefore, two different expansions leading to the same decay of the coefficients may not lead to an error of the same order for the expectation or of the variance. Of course, one could change these weights $\left|\mathbb{E}\left(\phi_{n}(p)\right)\right|$ by rescaling the $\phi_{n}$, but this only shifts the problem because the rescaling would then affect the decay of the coefficients.

In this work, we use the PC paradigm to minimize such quantities, by choosing an expansion basis $\phi_{n}$ that is adapted to the distribution of the noisy parameter $p$. PC expansions have become a very important tool in uncertainty quantification in the last decades, and a review of its many applications is far beyond the scope of the present work. We instead refer the interested reader to the survey [38] and the book [19].

The paper is structured as follows. First, we review the PC methodology in Section 2. In Section 3 we introduce some classical techniques to numerically study periodic orbits, invariant manifold and connecting orbits of deterministic ODEs. Then we proceed to the main contributions of our work. We explain, how to combine the numerical study of invariant sets with PC expansions to study the dynamics of nonlinear RODEs. In Section 4 we focus on a relatively simple example, where many quantities can also be computed analytically, which allows us to benchmark our numerical computations in the context of steady states. Yet, we also go beyond explicit structures and compute the distribution of stable/unstable manifolds using the “parameterization method” in combination with PC. Furthermore, we explain the implications of the random structure of the invariant sets and how to gain information from the moments of the invariants sets very efficiently. In Section 5 we then turn our attention to another example, where analytic computations are no longer available, and showcase the potential of our approach on the Lorenz system. For this system, we compute periodic orbits using a combined Fourier and PC ansatz, we compute again invariant stable/unstable manifolds, as well as heteroclinic connections.

2 A quick review on PC

Let $\rho:\mathbb{R}\to\mathbb{R}$ be probability distribution function (PDF) having finite moments, i.e.

[TABLE]

Given normalization constants $\left(h_{n}\right)_{n\in\mathbb{N}}$, $h_{n}>0$ for all $n\in\mathbb{N}$, there exists a unique family of orthogonal polynomials $\left(\phi_{n}\right)_{n\in\mathbb{N}}$ associated to the weight $\rho$, i.e. satisfying

[TABLE]

The most classical examples are:

• •

The Hermite polynomials $H_{n}$, which correspond to $\rho(s)=\frac{1}{\sqrt{2\pi}}{\textnormal{e}}^{-\frac{s^{2}}{2}}$ and $h_{n}=n!$;

• •

The Laguerre polynomials $L_{n}$, which correspond to $\rho(s)=\mathds{1}_{s\in[0,+\infty)}{\textnormal{e}}^{-s}$ and $h_{n}=1$;

• •

The Jacobi polynomials $P_{n}^{\alpha,\beta}$, $\alpha,\beta>-1$, which correspond to $\rho(s)=\mathds{1}_{s\in(-1,1)}\frac{(1-s)^{\alpha}(1+s)^{\beta}}{2^{\alpha+\beta+1}B(\alpha+1,\beta+1)}$ and $h_{n}=\frac{B(n+\alpha+1,n+\beta+1)}{(2n+\alpha+\beta+1)B(\alpha+1,\beta+1)B(n+1,n+\alpha+\beta+1)}$, where $B(x,y)=\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}$ is the Euler beta function.

Within the class of the Jacobi polynomials, we list a few remarkable cases (sometimes having different normalizations) that we make use of in this work:

• •

The Legendre polynomials $P_{n}$, which correspond to $\rho(s)=\frac{1}{2}\mathds{1}_{s\in(-1,1)}$ and $h_{n}=\frac{1}{2n+1}$;

• •

The Chebyshev polynomials of the first kind $T_{n}$, which correspond to $\rho(s)=\mathds{1}_{s\in(-1,1)}\frac{1}{\pi\sqrt{1-s^{2}}}$ and $h_{0}=1$, $h_{n}=\frac{1}{2}$, for $n\geq 1$;

• •

The Chebyshev polynomials of the second kind $U_{n}$, which correspond to $\rho(s)=\mathds{1}_{s\in(-1,1)}\frac{2}{\pi}\sqrt{1-s^{2}}$ and $h_{n}=1$;

• •

The Gegenbauer or ultraspherical polynomials $C_{n}^{\mu}$, $\mu>-\frac{1}{2}$, $\mu\neq 0$, which correspond to $\rho(s)=\mathds{1}_{s\in(-1,1)}\frac{2^{2\mu-1}\mu B(\mu,\mu)}{\pi}(1-s^{2})^{\mu-\frac{1}{2}}$ and $h_{n}=\frac{\mu}{n(n+\mu)B(n,2\mu)}$.

For a more complete description of PC choices and their relations to the Askey scheme, see [39].

Remark 2.1**.**

In this work, we only consider parameters having a PDF with bounded support. Indeed, in most applications these parameters have a physical meaning, for instance they could represent a quantity which must always be nonnegative, and having PDF with unbounded supports like Gaussian would mean that said parameters would be negative with a positive probability, which is not realistic. However, many sources of uncertainties are still expected to have a Gaussian-like behavior, in the sense that their PDF should be concentrated around a point. In that case, a good compromise is to use Beta distributions, which are the weights associated to the Gegenbauer polynomials and provide good bounded approximations of Gaussian distributions, at least for small variances (see e.g. [36], or the comparison on Figure 1). Another widely used option is to use truncated Gaussian distributions.

We recall that if $\rho$ has compact support, or decays at least exponentially fast at infinity (see e.g. [8]), then $\left(\phi_{n}\right)_{n\in\mathbb{N}}$ is a Hilbert basis of $L^{2}(\rho{\textnormal{d}}s)$, i.e. any mesurable function $g$ such that

[TABLE]

admits a unique series expansion of the form

[TABLE]

where the series converges in $L^{2}(\rho{\textnormal{d}}s)$.

Now, assume that the noisy parameter $p=p(\omega)$ (still assumed to be one dimensional for the moment) has a PDF given by $\rho$. The PC paradigm then tells us that we should use the orthogonal polynomials associated to $\rho$ as a basis for the expansion (2) of $X(p)$. Notice that with such a choice, the coefficients $X_{n}$ of the expansion directly provide us with the mean and variance of $X(p)$. Indeed, by orthogonality (and assuming $h_{n}=1$ for simplicity) we have

[TABLE]

and similarly

[TABLE]

If $p=(p^{(1)},\ldots,p^{(d)})$ consists of several independent random variables, each with respective PDF $\rho^{(j)}$, $j=1,\ldots,d$, one can consider the PC basis $\left(\phi_{n}\right)_{n\in\mathbb{N}^{d}}$ constructed as a tensor product of the univariates bases. That is, for all $n=(n_{1},\ldots,n_{d})\in\mathbb{N}^{d}$ and all $s=(s_{1},\ldots,s_{d})\in\mathbb{R}^{d}$,

[TABLE]

where $(\phi_{n_{j}}^{(j)})_{n_{j}\in\mathbb{N}}$ is a basis of orthogonal polynomials associated to the weight $\rho^{(j)}$. In practice, there are several ways to compute the coefficients $X_{n}$ of a PC expansion (2). Notice that directly using the orthogonality relations to get

[TABLE]

is usually not one of them, as $X$ is not known a-priori but is actually what we want to compute via the PC expansion. To formalize the discussion, assume that the quantity of interest $X(p)\in\mathbb{R}^{m}$ for which we want to find a PC expansion, solves a problem depending on a parameter $p\in\mathbb{R}^{d}$, of the form

[TABLE]

where $F:\mathbb{R}^{m}\times\mathbb{R}^{d}\to\mathbb{R}^{m}$. One common way to find the coefficients $X_{n}$ is to solve the system obtained by Galerkin projection:

[TABLE]

Of course, in practice one only solves for a truncated expansion

[TABLE]

by considering a finite-dimensional projection of associated dimension, i.e. $\forall~{}n_{2}\leq N$ in (3). This is the technique we make use of in this work. Another possible option is to use a collocation/interpolation approach [37]. This technique is based on first solving the deterministic problem several times for a well chosen sample of parameter values $\left(p_{i}\right)_{1\leq i\leq i}$, i.e. computing $X(p_{i})$ that solves

[TABLE]

The polynomials chaos coefficients $X_{n}$ are then constructed by interpolation:

[TABLE]

where $\alpha_{i}$ are weights associated to the interpolation points $p_{i}$ in such a way that, for any smooth function $g$

[TABLE]

where $\rho$ is the PDF of $p$. For a detailed discussion about these two approaches and their respective merits and limitations, we refer to the survey [38], the book [19] and the references therein. In practice, after an accurate PC approximation

[TABLE]

of $X$ has been computed, we can then do Monte Carlo simulations for a very large sample of values of $p$. Indeed, instead of having to solve the deterministic problem $F(X(p_{i}),p_{i})=0$ for each value $p_{i}$ of the sample, we only have to evaluate $\phi_{n}(p_{i})$ to obtain $X^{N}(p_{i})\approx X(p_{i})$.

3 Computation of periodic orbits, invariant manifolds and connecting orbits

In this section, we review some classical techniques to numerically study periodic orbits, invariant manifolds and connecting orbits of deterministic ODEs. An exhaustive review of these techniques is far beyond the scope of this work, and we only focus on one technique for each case, although alternative methods could certainly also be used; see e.g. [15] for a comparison of methods for computing stable/unstable manifolds. All the techniques presented here are based on series expansions, and this choice is motivated by two main reasons. The first and most important one is that series expansions can easily and efficiently be combined with PC expansions, once we go to the stochastic setting. This is particularly relevant for limit cycles, and is to be contrasted with more classical long-term integration approaches, for which PC expansions are known to be ill-behaved. The second one is that series expansions are particularly well suited for a-posteriori validation techniques [32], which we plan on developing in this context in a future work. In each case, we explain how an extra layer of PC expansion can be added, to keep track of the stochastic nature of the parameters. Illustrations for all these methods are presented in Section 4 and Section 5.

3.1 Periodic orbits via Fourier series

We start by considering the parameter dependent problem (1) from a deterministic point of view. Limit cycles of (1) can be efficiently studied and computed using Fourier series. That is, for a fixed $p$, we write a $T(p)$-periodic solution $t\mapsto X(t,p)$ of (1) as

[TABLE]

where $\Omega(p)=2\pi/T(p)$. To numerically find a periodic orbit, we thus solve for the coefficients $X_{k}(p)$ and $\Omega(p)$, which satisfy the system obtained by plugging (4) into (1), namely

[TABLE]

where $f_{k}(X(t,p),p)$ are the Fourier coefficients of $t\mapsto f(X(t,p),p)$.

Remark 3.1**.**

In practice, it is helpful to complement the above system, with a phase condition. For simplicity we choose to fix a section through which the solution has to pass at time [math], but other options are available [18]. More precisely, we add a scalar equation of the form

[TABLE]

where $u(p)$ is some (approximate) point on the orbit, $v(p)\approx f(u(p))$, and we use the dot product to denote the scalar product on $\mathbb{R}^{n}$ (and avoid potential confusion with the $L^{2}(\rho)$ scalar product $\langle,\rangle_{\rho}$). The phase condition allows to isolate the solution by eliminating time-shifts, which makes the system easier to solve in practice, especially using iterative methods. Indeed, if $t\mapsto X(t,p)$ is a periodic orbit of (1), then so is any function of the form $t\mapsto X(t+\tau,p)$, but there is (locally) only one $\tau$ that also satisfies the phase condition.

This approach of solving for the Fourier coefficients and the frequency has several advantages, compared to numerically integrating (1) and trying to find an (approximately) closed orbit. Indeed, the approach is not sensitive to the linear stability or instability of the limit cycle, and is therefore not susceptible of diverging if one tries to approximate an unstable periodic orbit, which is a definitive concern for time-integration based techniques. We illustrate this point in Section 5 by computing unstable limit cycles belonging to the chaotic attractor of the Lorenz system. The Fourier series approach is also particularly well adapted to continuation algorithms, especially to compute limit cycles originating from a Hopf bifurcation, where the linearized analysis can predict the first Fourier coefficient and the frequency.

If we now consider that $p$ is random and has a given PDF $\rho_{p}$, it is natural to still consider a Fourier ansatz [22]. It is straightforward to extend the Fourier series approach by expanding each Fourier coefficient, together with the frequency, with PC. Namely, write

[TABLE]

or equivalently

[TABLE]

In practice, we of course consider truncated expansions

[TABLE]

for which we solve using (5). Each $X_{k}$ now belongs to $\mathbb{R}^{mN}$ instead of $\mathbb{R}^{m}$, and $\Omega$ belongs to $\mathbb{R}^{N}$ instead of $\mathbb{R}$. An explicit example of such system together with numerical solutions is given in Section 5. In this setting, the Fourier series approach also has the notable advantage of separating the random period (or equivalently the random frequency $\Omega(p)$), from the description of the random cycle given by the coefficients $X_{k}(p)$. In particular, we avoid the usual pitfalls related to phase-drift and broadening of the spectrum, which are the main reasons why limit cycles are hard to compute using PC combined with time integration [7]. A similar idea was introduced in [20, 27], where a random time rescaling is used to compensate for the random period, which significantly improves the long time behavior of the PC expansions and allows to better capture stable limit cycles. Yet a Fourier series approach accomplishes that naturally, and can also be used to study unstable limit cycles.

3.2 Local invariant manifolds via Taylor series and the parameterization method

We again start by considering the parameter dependent problem (1) from a deterministic point of view, but now focus on studying local stable and unstable manifold attached to equilibrium points. In this section we only discuss the case of stable manifolds, but unstable manifolds can of course be studied with the same techniques. Let $\hat{X}(p)$ be an equilibrium point of (1), i.e. $f(\hat{X}(p),p)=0$, and assume that the derivative ${\textnormal{D}}_{x}f(\hat{X}(p),p)$ has exactly $m_{s}\leq m$ eigenvalues $\lambda^{(1)}(p),\ldots,\lambda^{(m_{s})}(p)$ with negative real part. For simplicity, we also assume that each of these eigenvalues is simple, and denote by $v^{(1)}(p),\ldots,v^{(m_{s})}(p)$ associated eigenvectors. Our goal is to find a parameterization $Q(p):\mathbb{R}^{m_{s}}\to\mathbb{R}^{m}$ of the local stable manifold of $\hat{X}(p)$. We look for a power series representation of $Q$:

[TABLE]

with the classical multi-indexes notations $|k|=k_{1}+\cdots+k_{m_{s}}$ and $\theta^{k}=\theta_{1}^{k_{1}}\cdots\theta_{m_{s}}^{k_{m_{s}}}$. Since we want $\theta\mapsto Q(\theta,p)$ to be a parameterization of the local stable manifold of $\hat{X}(p)$, we must have

[TABLE]

where $\gamma_{i}$ are scaling that can be adjusted. To obtain the higher order terms, we follow the idea of the parameterization method, introduced in [3, 4, 5] (see also the recent book [14]). We want to obtain a parameterization that conjugates the dynamics on the stable manifold with the stable dynamics of the linearized system. More precisely, introducing the diagonal matrix $\Lambda(p)$ with diagonal entries $\lambda^{(1)}(p),\ldots,\lambda^{({m_{s}})}(p)$, we want $p$ to satisfy (see Figure 2)

[TABLE]

where $\varphi_{p}$ is the flow generated by the vector field $f(\cdot,p)$ and $\left\|\theta\right\|_{\infty}\leq 1=\max_{1\leq i\leq{m_{s}}}|\theta_{i}|$.

Remark 3.2**.**

If some of the stable eigenvalues are complex conjugate, say $\lambda^{(1)}=\bar{\lambda}^{(2)},\ldots,\lambda^{(2r-1)}=\bar{\lambda}^{(2r)}$, it is more convenient to first look for a complex valued parameterization $Q$ with $Q_{k}\in\mathbb{C}^{m}$, and then recover a real valued parameterization via

[TABLE]

see e.g. [31].

Finding a parameterization $Q$ satisfying (9) is interesting because it provides us with not only a local stable manifold but also an explicit description of the dynamics on this manifold. However, the formulation (9) is not the most convenient one to work with in order to determine the higher order coefficients of $Q$, because it involves the flow. To get rid of it, one can take a time derivative of (9) and evaluate at $t=0$, to obtain the following invariance equation

[TABLE]

One can check that, if $Q$ solves (10) and is such that $Q(0,p)=\hat{X}(p)$, then $Q$ satisfies (9), therefore $Q(\cdot,p)$ is indeed a parameterization of the local unstable manifold of $\hat{X}(p)$. The invariance equation (10) is the one we are going to use to numerically find the coefficients $Q_{k}$. Plugging the expansion (7) in the invariance equation (10), we obtain

[TABLE]

where $f_{k}(Q(\theta,p),p)$ are the Taylor coefficients of $\theta\mapsto f(Q(\theta,p),p)$ and $\lambda(p)\cdot k=\lambda^{(1)}(p)k_{1}+\ldots+\lambda^{({m_{s}})}(p)k_{m_{s}}$, therefore the coefficients $Q_{k}$ must satisfy

[TABLE]

Notice that (8) already ensures that (12) is satisfied for $|k|\leq 1$. In order to solve (12) for $|k|\geq 2$, let us first describe, how $f_{k}(Q(\theta,p),p)$ depends on $Q_{k}(p)$. Given a Taylor series of the form (7), we define for all $K\in\mathbb{N}$ the truncated series

[TABLE]

Notice also that, for any $k\in\mathbb{N}^{m_{s}}$, $f_{k}(Q(\theta,p),p)=f_{k}(\pi_{|k|}Q(\theta,p),p)$. Besides, using a Taylor expansion of $f$ in the $x$ variable, we have for all $k\in\mathbb{N}^{m_{s}}\setminus\{0\}$

[TABLE]

and looking at the coefficient of degree $k$ on each side we get

[TABLE]

Therefore, assuming (8), having (12) for all $|k|\geq 2$ is equivalent to having

[TABLE]

Assuming the following non-resonance condition is satisfied:

[TABLE]

we see that (13) has a unique solution that can be computed recursively via

[TABLE]

since the right-hand side only depends on $Q_{l}(p)$ for $|l|<|k|$. We can therefore compute a truncated parameterization $\pi_{K}Q$ of arbitrary order, starting from (8) and then computing (14) recursively for $k$ as large as desired. In practice, the weights $\gamma_{i}$ in (8) are chosen in order to obtain a reasonnable decay of the coefficients $Q_{k}(p)$ (we refer to [2] for a detailed explanation of how this choice can be optimized). Explicit examples are presented in Sections 4 and 5.

Remark 3.3**.**

In cases where resonant eigenvalues are present, a similar approach can still be used, but the conjugacy condition (9) defining $Q$ has to be adapted [31].

If we now consider that $p$ is random and has a given PDF $\rho_{p}$, this approach based on the parameterization method can also be easily generalized by adding a layer of PC expansion. Namely, we write

[TABLE]

or equivalently

[TABLE]

In practice we consider a truncation

[TABLE]

which we compute via

[TABLE]

where, compared to (14), $Q_{k}^{N}(p)$ is now a vector of size $mN$ rather than $m$, and

[TABLE]

can be interpreted as a block matrix, with $m\times m$ blocks having each size $N\times N$. For explicit computations we again refer to Sections 4 and 5.

3.3 Heteroclinic orbits via Chebyshev series and projected boundaries

We go back to considering the parameter dependent problem (1) from a deterministic point of view, and extend the discussion of the previous subsection by looking at more global solutions, namely connecting orbits between equilibrium points.

Let $\hat{X}(p)$ and $\check{X}(p)$ be two equilibrium points of (1). Assume that $\hat{X}(p)$ has an unstable manifold of dimension $\hat{r}$ and that $\check{X}(p)$ has a stable manifold of dimension $\check{r}$, such that $\hat{r}+\check{r}=m+1$. Then, if the two manifolds intersect, we can generically expect this intersection to be transverse in the phase space $\mathbb{R}^{m}$, in which case there exists a transverse heteroclinic orbit between $\hat{X}(p)$ and $\check{X}(p)$. We recall that a heteroclinic orbit between $\hat{X}(p)$ and $\check{X}(p)$ (or homoclinc orbit if $\hat{X}(p)=\check{X}(p)$) is a solution $t\mapsto x(t,p)$ of (1) such that

[TABLE]

In this work, we compute such solution by solving a boundary value problem [1] between the unstable manifold of $\hat{X}(p)$ and the stable manifold of $\check{X}(p)$, for which we first compute local parameterization as in Section 3.2. This allows us to only solve (1) on a finite time interval, and recover the remaining parts of the orbit via the conjugacy satisfied by the parameterizations. More precisely, we want to find an orbit $t\mapsto X(t,p)$ such that

[TABLE]

where $W^{\textnormal{u}}(\hat{X}(p))$ and $W^{\textnormal{s}}(\check{X}(p))$ respectively denote the unstable manifold of $\hat{X}(p)$ and the stable manifold of $\check{X}(p)$. Note that, since we only consider autonomous vector fields in this work, only the length $\tau(p)$ of the time interval $[0,\tau(p)]$ is relevant, and the whole interval itself could of course be shifted. To numerically compute such an orbit, we use piece-wise Chebyshev series. In order to do so, we first introduce some notations.

For $k\in\mathbb{N}$, we denote by $T_{k}$ the Chebyshev polynomial of order $k$, defined for instance by $T_{k}(\cos(\theta))=\cos(k\theta)$. We introduce a partition $0=\tilde{t}^{(0)}<\tilde{t}^{(1)}<\ldots<\tilde{t}^{(J)}=1$ of $[0,1]$. We denote by $\tau(p)$ the (unknown) time the orbit spends between the two local manifolds and consider the partition of $[0,\tau(p)]$ given by $0=t^{(0)}(p)<t^{(1)}(p)<\ldots<t^{(J)}(p)=\tau(p)$, where $t^{(j)}(p)=\tau(p)\tilde{t}^{(j)}$ for all $1\leq j\leq J$. For all $k\in\mathbb{N}$ and $j=1,\ldots,J$, we also introduce the rescaled Chebyshev polynomial $T^{(j)}_{k}(p)$, defined as

[TABLE]

We can then write the orbit between the two manifolds as

[TABLE]

Finally, we assume that, following the methodology presented in Section 3.2, a truncated parameterization $\hat{Q}(p)$ of the local unstable manifold of $\hat{X}(p)$ as well as a truncated parameterization $\check{Q}(p)$ of the local stable manifold of $\check{X}(p)$ have been computed. Rewriting the differential equation in (LABEL:eq:BVP) as an integral one, plugging in the expansion (17) and using well known properties of the Chebyshev polynomials, namely

[TABLE]

we obtain

[TABLE]

where $f^{(j)}_{k}(X(t,p),p)$ are the Chebyshev coefficients of $t\mapsto f(X(t,p),p)$ on $(t^{(j-1)}(p),t^{(j)}(p))$. The first line in (18) corresponds to the differential equation $\dot{X}(t,p)=f(X(t,p),p)$ on each subinterval $(t^{(j-1)}(p),t^{(j)}(p))$, the second line insures that the solution connects continuously between two consecutive subintervals, and the last two lines correspond to the boundary conditions on each manifold.

Remark 3.4**.**

In (18), besides the unknown Chebyshev coefficients $X_{k}^{(j)}(p)$, we have $\hat{r}+\check{r}+1$ additional scalar unknowns: $\tau(p)\in\mathbb{R}$, $\hat{\theta}(p)\in\mathbb{R}^{\hat{r}}$ and $\check{\theta}(p)\in\mathbb{R}^{\check{r}}$. Still assuming $\hat{r}+\check{r}=m+1$, the system is then underdetermined and we can fix two of these $m+2$ parameters to recover a unique solution.

Similarly to the two previous cases of Fourier and Taylor series, when we consider $p$ as a random parameter having a given distribution we only need to expand every unknown in (18) using PC. We again refer to Section 5 for an explicit example.

4 First example: a Lotka-Volterra system

In this section, we compute steady states and invariant manifolds of a Lotka-Volterra system of competition type

[TABLE]

where $b=b(\omega)$ is a bounded random variable having a given distribution. This basic example allows us to easily study the quality of our numerical computations by comparing them to analytic results.

4.1 Analytical results (for benchmarking)

4.1.1 Deterministic framework

Assuming $ab\neq 1$, system (19) has a non trivial equilibrium given by

[TABLE]

The eigenvalues of the Jacobian ${\textnormal{D}}f(x_{\textnormal{eq}},y_{\textnormal{eq}})$ are

[TABLE]

and associated eigenvectors are given by

[TABLE]

In the case $a,b>1$, $\lambda^{(2)}$ is positive and the equilibrium is of saddle type. Its stable manifold is the line $y=\frac{b-1}{a-1}x$, $x>0$ and its unstable manifold is a (nonlinear) curve connecting to $(0,1)$ and $(1,0)$.

4.1.2 Stochastic framework

We now assume that the parameter $b=b(\omega)$ is a bounded random variable. We write

[TABLE]

where $\sigma\geq 0$ and $p$ is a random variable taking values in $[-1,1]$. We assume $a\bar{b}\neq 1$ and $\sigma\leq\left|\frac{a\bar{b}-1}{a}\right|$, so that $ab(\omega)-1\neq 0$ a.e.. For specific distributions of $p$, one can compute analytically the first two moments of $x_{\textnormal{eq}}=\frac{a-1}{ab-1}$:

• •

If $p$ is uniformly distributed on $[-1,1]$, that is if its PDF is given by $\rho_{p}(s)=\frac{1}{2}\mathds{1}_{s\in(-1,1)}$, one has

[TABLE]

• •

If $p$ is beta-distributed with parameters $(-\frac{1}{2},-\frac{1}{2})$, that is if its PDF is given by $\rho_{p}(s)=\mathds{1}_{s\in(-1,1)}\frac{1}{\pi\sqrt{1-s^{2}}}$, one has

[TABLE]

• •

If $p$ is beta-distributed with parameters $(\frac{1}{2},\frac{1}{2})$, that is if its PDF is given by $\rho_{p}(s)=\mathds{1}_{s\in(-1,1)}\frac{2}{\pi}\sqrt{1-s^{2}}$, one has

[TABLE]

In order to focus the amount of comparisons and illustrations in the next section, we only consider the first component $x_{\textnormal{eq}}$ of the equilibrium, but similar analytical and numerical computations could of course also be carried out for $y_{\textnormal{eq}}$.

4.2 Numerical results

In this section, we compute using PC the equilibrium $(x_{\textnormal{eq}},y_{\textnormal{eq}})$ and its stable and unstable manifold. For given $a$, $\bar{b}$ and $\sigma$, we use the techniques presented in Section 3 on

[TABLE]

where for convenience we use $(x,y)$ instead of $(x_{1},x_{2})$.

4.2.1 Analysis of convergence on the steady state problem

We first solve, for various choices of expansions, the steady state problem $f(x(p),y(p),p)=0$, and analyze how the choice of expansion together with the distribution of $p$ affect the convergence rates. More precisely, for a given expansion basis $\phi_{n}$ and truncation level $N$, we look for coefficients $\left(x_{n}\right)_{0\leq n<N}$ and $\left(y_{n}\right)_{0\leq n<N}$ such that

[TABLE]

where

[TABLE]

By taking the scalar product with respect to each $\phi_{n}$, $0\leq n<N$, we obtain the following system

[TABLE]

Here $\ast$ denotes the product associated to the basis $\phi_{n}$, i.e. $(x\ast y)_{0\leq n<2N-1}$ is the unique sequence such that

[TABLE]

If $\phi_{n}$ is chosen so that $\phi_{1}(p)=p$, which will often be the case in practice, then the expansion of $p$ only has a single non-zero coefficient: $p_{n}=\delta_{n,1}$.

Remark 4.1**.**

We emphasize that $\ast$ of course depends on the choice of the basis $\phi_{n}$. We refer to the Appendix for more details about this product structure, and discussions on how it can be computed in practice.

We solve (23) using Newton’s method, and the initial conditions given by the deterministic equilibrium, i.e., for $\sigma=0$ so that

[TABLE]

The results for different choices of basis $\phi_{n}$ are presented in Figure 3. To make the comparison fair, we consider rescaled version of the Chebyshev polynomials of the second kind and of the Gegenbauer polynomials (see the Appendix), so that all $\phi_{n}$ used satisfy $\sup_{s\in[-1,1]}|\phi_{n}(s)|=1$. The decay of the coefficients is therefore indicative of the truncation error in ${\mathcal{C}}^{0}$-norm. As expected, the Chebyshev and Legendre expansions converge faster than the Taylor expansions, with the Gegenbauer one lying somewhere in between. As mentioned previously, we focus here only on the first component $x$, but the same behavior is observed for the second component $y$.

However, as mentioned in the introduction, we are more interested in understanding the distribution of the solution, assuming $p$ has a prescribed distribution in $[-1,1]$. In Figure 4 and Figure 5 we display the relative error for the first two moments

[TABLE]

depending on the truncation level $N$, for several PDF $\rho_{p}$ of $p$ and again for several expansions.

The limiting moments $\mathbb{E}(x(p))$ and $\mathbb{E}(x(p)^{2})$ are computed analytically when possible (see Section 4.1.2) and using numerical integration otherwise. The truncated moments $\mathbb{E}(x^{N}(p))$ and $\mathbb{E}(x^{N}(p)^{2})$ are computed from the coefficients of the expansion

[TABLE]

If $\phi_{n}$ is the orthogonal basis associated to the distribution $\rho_{p}$, then as explained in Section 2

[TABLE]

Otherwise, we have

[TABLE]

and we first compute the above integrals to obtain the truncated moments. As expected from the paradigm of PC, we obtain faster convergence for these statistics when using the PC expansions associated to the distribution $\rho_{p}$ of the parameter $p$, i.e., a family of polynomial $\phi_{n}$ orthogonal with respect to $\left\langle\cdot,\cdot\right\rangle_{\rho_{p}}$. In particular, if $p$ has a uniform distribution then the Legendre expansion converges the fastest (Figure 4(a) and Figure 5(a)), if $p$ has an arcsine distribution then the Chebyshev expansion of the first kind converges the fastest (Figure 4(b) and Figure 5(b)), if $p$ has a Wigner semicircle distribution then the Chebyshev expansion of the second kind converges the fastest (Figure 4(c) and Figure 5(c)) and finally if $p$ has a beta distribution of parameter $(20,20)$ then the Gegenbauer expansion of parameter $\mu=20$ converges the fastest (Figure 4(d) and Figure 5(d)).

Remark 4.2**.**

These simple experiments confirm that it is in general a good option to choose the expansion basis associated to the distribution of the parameter, and we will systematically do so in the sequel. Nonetheless, we believe that in some cases, especially highly nonlinear ones, the cost of actually computing the nonlinear terms should also be taken into account (see the discussion in Section 7.3).

4.2.2 Computation of eigenvalues and eigenvectors

Now that we found an accurate representation of the steady state, we turn our attention to the dynamics around it. Before computing parameterizations of the local stable and unstable manifolds, which will be done in the next subsection, we first focus on the linearized system around the equilibrium. More precisely, we are interested in the usual eigenvalue problem

[TABLE]

where $(x(p),y(p))$ is the previously computed equilibrium, and we solve for $u(p)$, $v(p)$ and $\lambda(p)$, the last equation being a normalization of the eigenvector allowing us to have a (locally) unique solution. In the previous subsection we obtained truncated expansions $x^{N}$ and $y^{N}$ of the equilibrium, and we now aim at doing the same for the eigendata. Namely, we write

[TABLE]

and plug these expansions back in (24). For our explicit example (22), this yields the following system for the coefficients $u_{n}$, $v_{n}$ and $\lambda_{n}$:

[TABLE]

Again we solve this system using Newton’s method and the initial conditions given by numerically obtained eigenvalue and eigenvector for the deterministic problem (20)-(21). We obtain PC expansions for the eigenvalues $\lambda^{(1)}$ and $\lambda^{(2)}$ and associated eigenvectors $V_{\lambda^{(1)}}=\left(u^{(1)},v^{(1)}\right)$ and $V_{\lambda^{(2)}}=\left(u^{(2)},v^{(2)}\right)$.

From these PC expansions, we can then easily get statistical information about the eigenvalues and eigenvectors. For instance, as mentioned in the introduction, once the PC coefficients have been computed one can carry out Monte Carlo simulations basically for free (the only cost being the evaluation of the basis polynomials $\phi_{n}$, for $0\leq n<N$, at the sampled values of $p$). For a given distribution of $p$, this allows us to numerically compute the PDF of $\lambda$ (see Figure 6). Another approach, maybe more adapted to quantities of interest that are more than one-dimensional, is to use the mean and variance that can be obtained from the PC expansion to get a sense of where objects may lie in phase space (see Figure 7). Since the PC expansion also allows us to compute higher order moments (see Section 7.2), this geometrical description can be made more quantitative by using these moments to obtain concentration inequalities.

4.2.3 Computation of local stable/unstable manifolds

We are now ready to compute parameterizations of the stable and unstable manifolds of the equilibrium, as described in Section 3.2. That is, we look for a parameterization of the form

[TABLE]

Notice that for our example (22) the phase space is two-dimensional and both the stable and the unstable manifolds that we are interested in are one-dimensional. Therefore, $\theta$ will be one-dimensional and the coefficients $Q_{k,n}=\left(Q^{(1)}_{k,n},Q^{(2)}_{k,n}\right)$ belong to $\mathbb{R}^{2}$. It will be convenient to use the following notation

[TABLE]

According to (8) we define

[TABLE]

where $\left(x_{n}\right)_{n}$ and $\left(y_{n}\right)_{n}$ are the PC coefficients of the equilibrium computed in Section 4.2.1, and $(u^{(i)}_{n})_{n}$ and $(v^{(i)}_{n})_{n}$ are the PC coefficients of the eigenvectors computed in Section 4.2.2. Observe that $i=1$ if we want to compute a parameterization of the stable manifold (resp. $i=2$ for the unstable manifold). To compute the higher order coefficients $Q_{k}$, $2\leq k<K$ we use (15), which we now specialize to our example (22). We need to introduce some more notations. Given an expansion basis $\phi_{n}$ and $z=\left(z_{n}\right)_{0\leq n<N}\in\mathbb{R}^{N}$, we denote by $M_{z}$ the $N\times N$ * product matrix* such that

[TABLE]

This matrix can be computed explicitly from $z$ and the linearization coefficients $\alpha$ (see the Appendix). Given an expansion basis $\phi_{n}$ and double expansions of the form

[TABLE]

we denote by $\circledast$ the product associated to the basis $\phi_{n}(p)\theta^{k}$, that is $\left(q\circledast q^{\prime}\right)_{k,n}$ is the unique sequence such that

[TABLE]

We again refer to the Appendix for more details about this product structure, and discussions on how it can be computed in practice. Similarly to (25) we write

[TABLE]

We are now ready to go back to (15) and give explicit and commutable formulas for the coefficients $Q_{k}$ $k\geq 2$ of the parameterization in the case of our example (22). Assuming the coefficients $Q^{(1)}_{l}$ and $Q^{(2)}_{l}$ have been computed for $0\leq l\leq k-1$, the next coefficients are given by

[TABLE]

where $I_{N}$ is the $N\times N$ identity matrix and again $i=1$ if we are computing a parameterization of the stable manifold (resp. $i=2$ for the unstable manifold). Here $p$ is identified with the sequence $p_{k,n}=\delta_{k,0}\delta_{n,1}$ (again assuming the $\phi_{n}$ are normalized such that $\phi_{1}(p)=p$). Once the coefficients $Q^{(1)}_{k}$ and $Q^{(2)}_{k}$ have been computed up to the desired order, we have access to the statistical properties of the manifold. Indeed we can for instance compute

[TABLE]

and

[TABLE]

which describe the mean and variance of the first component of the manifold with respect to $\theta$ (see Figures 8 and 9). We can also carry out cheap sampling to get a sense of how the manifolds are distributed in phase space (see again Figures 8 and 9).

5 Second example: the Lorenz system

In this section, we compute periodic orbits, invariant manifolds and heteroclinic orbits for the Lorenz system

[TABLE]

where $\varrho=\varrho(\omega)$ is a bounded random variable having a given distribution. We renormalize it by writing

[TABLE]

where $\sigma>0$ and $p$ is a random variable taking values in $[-1,1]$.

5.1 Computation of periodic orbits

Using the framework presented in Section 3.1, we compute periodic orbits of (27) via a Fourier$\times$PC expansion. That is, we write

[TABLE]

and similarly for $y^{K,N}$ and $z^{K,N}$, where

[TABLE]

In this subsection, we use $\circledast$ to denote the product associated to the basis $\phi_{n}(p){\textnormal{e}}^{{\textnormal{i}}k\Omega(p)t}$, i.e. $(x\circledast y)_{k,n}$ is the unique sequence such that

[TABLE]

It will again be convenient to write

[TABLE]

and similarly

[TABLE]

and so on. Plugging the expansions (28) for $x$ , $y$ and $z$ in (27) and using the notations (29) and (30), we obtain the following system of $(2K-1)\times N$ equations:

[TABLE]

where again $p$ is identified with the sequence $p_{k,n}=\delta_{k,0}\delta_{n,1}$ and $\Omega$ with the sequence $\Omega_{k,n}=\delta_{k,0}\Omega_{n}$. As mentioned in Section 3.1, we complement this system with a phase condition, that we also expand using PC. That is, given

[TABLE]

we append to (31) the $N$ equations

[TABLE]

In practice, we solve (31)-(32) using Newton’s method. More precisely, we assume that we are given a deterministic periodic orbit (i.e. a solution for $\sigma=0$) and do a predictor-corrector continuation in $\sigma$ (using Newton’s method as the corrector step) until we reach the desired value. In Figure 10 we illustrate the output of this procedure.

5.2 Computation of heteroclinic orbits

In this subsection, we detail how our approach can be used to compute heteroclinic orbits for the Lorenz system, going between

[TABLE]

and the origin. We focus on the case where

[TABLE]

in which the origin has a two-dimensional stable manifold and (33) has a two-dimensional unstable manifold.

As explained in Section 3.3, we use these local manifolds to set up our boundary value problem for the heteroclinic orbits. A parameterization of these local manifolds is computed as described in Section 3.2, using a Taylor$\times$PC expansion. This procedure was already described in details in Section 4.2.3, so we omit these details here. We denote by $\hat{Q}^{\hat{K},N}$ (resp. $\check{Q}^{\check{K},N}$) a local parameterization of the unstable manifold of (33) (resp. of the stable manifold of the origin) of the form

[TABLE]

Notice that since both manifolds are two-dimensional, we have $\theta=(\theta_{1},\theta_{2})$. We recall that our goal is to find $\tau(p)$, $\hat{\theta}(p)$, $\check{\theta}(p)$ and an orbit $X(t,p)$ that solves the Lorenz system on $[0,\tau(p)]$, and satisfies the boundary conditions $X(0,p)=\hat{Q}^{\hat{K},N}(\hat{\theta}(p),p)$ and $X(\tau(p),p)=\check{Q}^{\hat{K},N}(\check{\theta}(p),p)$. As mentioned in Remark 3.4, this system is underdetermined and we can in fact fix $\hat{\theta}_{1}$ and $\check{\theta}_{1}$, and only solve for $\tau(p)$, $\hat{\theta}_{2}(p)$, $\check{\theta}_{2}(p)$ and $X(t,p)$ to recover a unique solution. In practice, we solve for PC expansions of $\tau(p)$, $\hat{\theta}_{2}(p)$ and $\check{\theta}_{2}(p)$:

[TABLE]

To compute the heteroclinic orbit (or to be more precise, the part of that orbit that connects the two local manifolds), we use piece-wise Chebyshev$\times$PC expansions, as exposed in Section 3.3. That is, we now write

[TABLE]

and similarly for $y^{K,N}$ and $z^{K,N}$, where we use the notations of Section 3.3 for partition of $[0,\tau]$ and the associated rescaled Chebyshev polynomials. In this subsection, we use $\circledast$ to denote the product associated to the basis $\phi_{n}(p)T_{|k|}^{(j)}(t,p)$, i.e. $(x\circledast y)_{k,n}$ is the unique sequence such that

[TABLE]

It will again be convenient to write

[TABLE]

and similarly

[TABLE]

[TABLE]

and so on. To write all three components at once, we also use

[TABLE]

With these notations, system (18) for the Lorenz vector field is then given by

[TABLE]

where, for $k=(k_{1},k_{2})$, $\hat{\theta}^{k}$ must be understood as

[TABLE]

Again, we solve this system for the desired noise level $\sigma$ by doing predictor corrector steps with Newton iterations. In Figure 11 we illustrate the output of this procedure.

We point out that, while our method is also successful when $\varrho$ takes values around the classical value $\varrho=28$, the computation is then more challenging due to the proximity of the Hopf bifurcation at $\varrho=\frac{\varsigma\left(\varsigma+\beta+3\right)}{\varsigma+\beta-1}$ ($\approx 24.7$ for the classical values $\varsigma=10$ and $\beta=8/3$), and therefore we are only able to handle smaller noise level (i.e. $\sigma\approx 1$, see Figure 12). For such parameter values, the computations (in both the stochastic and the deterministic framework) are intrinsically more difficult because of the strong oscillatory behavior induced by the pair of eigenvalues having a close to zero real part.

6 Conclusions & Outlook

PC expansions have proven very successful in studying evolution equations (be it ODEs or PDEs) with random coefficients, but mostly from the point of view of time integration. In this work we developed a complementary approach, also based on PC but aimed at studying invariant sets of such systems, and applied it to investigate steady states, stable and unstable manifolds, periodic orbits and connecting orbits for ODEs. This approach is driven by the paradigm of nonlinear dynamics to study the structure of invariant sets to understand the relevant effects. Once a PC expansion representation of an invariant set is computed, this expansion contains explicit information about the random invariant sets. Using fast sampling or geometric visualization of the moments, then allows us to understand, how likely different phase space structures are going to be based upon the random input.

We conclude by providing a brief discussion about some potential generalizations and further directions of research connected to this work.

Extension to PDEs:

A natural extension of this work would be to use similar techniques to study invariant sets of time-dependent PDEs with random coefficients. The steady state case has already been extensively investigated, but we believe that our new framework could allow to complement the already existing studies on periodic orbits, and to explore new problems related to invariant manifolds and connecting orbits for PDEs with random coefficients.

Bifurcations:

If the invariant state we are studying undergoes a bifurcation for some value of the random parameter, then the curve or manifold of invariant state that we are looking for may not be smooth, and the PC expansions will then converge slowly, if at all. Several multiresolution analysis schemes, based on subdivising the random space or on different bases such as wavelets, were developed to handle such situations (see e.g. [34],[19, Chapter 8] and the references therein), and could be used also in our setting to study bifurcation problems.

Rigorous computation:

A posteriori error analysis for PC expansions is of course critical, as the quantification of uncertainty provided by those expansions is only relevant if the error coming from the discretization/truncation can also be controlled. Techniques of rigorous computations have been developed to obtain certified a-posteriori error estimates about numerically computed invariant sets of deterministic systems [9], and we aim at generalizing them for random systems in a future work.

7 Appendix

We discuss here the computations of products of PC expansions, and detail some implementation aspects.

7.1 Linearization formulas for products

For any family of orthogonal polynomials $\phi_{n}$ associated to a weight $\rho$, and any $m$ and $n$ in $\mathbb{N}$, the product $\phi_{m}\phi_{n}$ can be written in the original basis:

[TABLE]

This is often called a linearization formula. By orthogonality, the linearization coefficients satisfy

[TABLE]

from which we infer that $\alpha^{m,n}_{k}=\alpha^{n,m}_{k}$ and that, if $k+n<m$ or $k+m<n$, then $\alpha^{m,n}_{k}=0$, i.e.

[TABLE]

Besides, if the weight $\rho$ is even, then any even function is orthogonal to any odd function and therefore, for all $k$ having a different parity than $m+n$, $\alpha^{m,n}_{k}=0$. In such case, it can be convenient to eliminate all the coefficients that are a priori equal to zero, and rewrite (35) as

[TABLE]

Remark 7.1**.**

In practice, it is convenient to precompute and store the linearization coefficients $\alpha^{m,n}_{k}$ (or $\tilde{\alpha}^{m,n}_{k}$). These coefficients can be obtained by numerically computing the integrals $\left\langle\phi_{m}\phi_{n},\phi_{k}\right\rangle_{\rho}$, for instance using quadrature rules. For the classical orthogonal polynomials, the linearization coefficients can also be computed in closed form (see e.g. [24]). For Legendre polynomials $P_{n}$, we have

[TABLE]

where

[TABLE]

For Chebyshev polynomials of the first kind $T_{n}$, we have

[TABLE]

For Chebyshev polynomials of the second kind $U_{n}$, we have

[TABLE]

For Gegenbauer polynomials $C^{\mu}_{n}$, we have

[TABLE]

Finally, we point out that linearization coefficients $\alpha^{m,n}_{k}$ can of course also be defined for arbitrary (i.e. non necessarily orthogonal) bases of polynomials. In particular, for the canonical basis $\phi_{n}(s)=s^{n}$ associated to Taylor expansions we have

[TABLE]

Definition 7.2**.**

Given to sequences $u=\left(u_{n}\right)_{n\in\mathbb{N}}$ and $v=\left(v_{n}\right)_{n\in\mathbb{N}}$, we define their convolution product $u\ast v$ by

[TABLE]

with the convention $\alpha^{m,n}_{k}=0$ for all $k>m+n$. Notice that this definition depends on the weight $\rho$, or equivalently on the family $\left(\phi_{n}\right)_{n\in\mathbb{N}}$, via the coefficients $\alpha^{m,n}_{k}$.

Lemma 7.3**.**

If the weight $\rho$ is even, the convolution of $u$ and $v$ can also be written as

[TABLE]

This definition is the natural one to describe the product of two functions in the basis given by $\left(\phi_{n}\right)_{n\in\mathbb{N}}$. Indeed, writing

[TABLE]

one has, at least formally,

[TABLE]

With the notations of Definition 7.2, one has

[TABLE]

where $\left\|u\right\|_{1}=\sum_{n=0}^{\infty}|u_{n}|$. Therefore, as soon as

[TABLE]

the space $\ell^{1}$ of sequences with finite $\left\|\cdot\right\|_{1}$ norm is stable under the convolution product, i.e. $\left(\ell^{1},\ast\right)$ is a Banach algebra.

Remark 7.4**.**

If the family $\phi_{n}$ is such that $\sup_{s\in[-1,1]}|\phi_{n}(s)|\leq 1$ for all $n\in\mathbb{N}$, then the $\ell^{1}$-norm of the coefficients controls the ${\mathcal{C}}^{0}$-norm of the function:

[TABLE]

In particular, if the coefficients associated to $u$ and $v$ have finite $\ell^{1}$-norm and (39) is satisfied, then the sum in (38) is guaranteed to converge, and (38) holds not only formally, but also in ${\mathcal{C}}^{0}$.

For each family of polynomials used in this work, namely the Legendre polynomials $P_{n}$, the Chebyshev polynomials of the first kind $T_{n}$, the canonical basis $X^{n}$, and (suitable renormalization of) the Chebyshev polynomials of the second kind $\tilde{U}_{n}=U_{n}/U_{n}(1)$ and the Gegenbauer polynomials $\tilde{C}^{\mu}_{n}=C^{\mu}_{n}/C^{\mu}_{n}(1)$, we have

[TABLE]

Indeed, in all these cases $\phi_{n}(1)=1$ for all $n\in\mathbb{N}$, therefore (35) yields

[TABLE]

Besides, in all these cases the linearization coefficients $\alpha^{m,n}_{k}$ are nonnegative for all $m,n,k\in\mathbb{N}$, and hence we get (40), which implies (39). Notice that in all those cases $\sup_{s\in[-1,1]}|\phi_{n}(s)|=1$ for all $n\in\mathbb{N}$, and therefore Remark 7.4 applies. For a more in depth discussion about the convolution products and Banach algebra structures that can be associated to orthogonal polynomials, we refer to the lecture notes [28] and the references therein.

In practice we need to implement the linear operator representing the (truncated) convolution with a given sequence (see e.g. Section 4.2.3). Given $u=\left(u_{n}\right)_{0\leq n<N}$, and linearization coefficients $\alpha$ associated to a basis $\left(\phi_{n}\right)$, we therefore define the $N\times N$ matrix $M_{u}$ by

[TABLE]

It follows from (36) that, for any vector $v=\left(v_{n}\right)_{0\leq n<N}$

[TABLE]

Let us now consider two basis $\left(\phi^{(1)}_{n_{1}}\right)_{n_{1}\in\mathbb{N}}$ and $\left(\phi^{(2)}_{n_{2}}\right)_{n_{2}\in\mathbb{N}}$ having linearization coefficients $\alpha^{(1)}$ and $\alpha^{(2)}$ and associated convolution products $\ast^{(1)}$ and $\ast^{(2)}$. We encountered this situation in Section 4 and Section 5, for instance when dealing with Taylor$\times$PC expansions, but what follows could also be used for multivariate PC expansions. Given bidimensional sequences $u=\left(u_{n}\right)_{n\in\mathbb{N}^{2}}$ and $v=\left(v_{n}\right)_{n\in\mathbb{N}^{2}}$ associated to expansions

[TABLE]

we introduce the bidimensional product $\circledast$ such that

[TABLE]

The coefficients $\left(u\circledast v\right)_{n}$ can of course be computed from the univariate convolution products, for instance by looking at $u$ and $v$ as functions of one variable (say $x$) having coefficients depending on the other variable (say $y$):

[TABLE]

Denoting $u^{(1)}=\left(u^{(1)}_{n_{1}}\right)_{n_{1}\in\mathbb{N}}$, where for all $n\in\mathbb{N}$ $u^{(1)}_{n_{1}}=\left(u_{n_{1},n_{2}}\right)_{n_{2}\in\mathbb{N}}$, and slightly abusing the notation $\ast^{(1)}$ we get

[TABLE]

Of course the role of each variable/basis is interchangeable and we also have

[TABLE]

7.2 Treatment of higher order terms

Let us now consider situation where product of three or more expansions have to be computed for a given family $\phi_{n}$. Of course one could consider the coefficients $\beta^{k,m,n}_{l}$ such that

[TABLE]

and precompute (and store) them using the formula:

[TABLE]

However, this approach becomes impracticable very fast even for univariate bases when the degree of the nonlinearity increases. Therefore in practice it is more efficient to stick with only the coefficients $\alpha$ for quadratic products, and compute nonlinear term of higher order recursively. That is, given expansions

[TABLE]

we get the coefficients of $u(s)v(s)w(s)$ by first computing $u\ast v$ and then $(u\ast v)\ast w$, which by associativity can simply be denoted $(u\ast v\ast w)$.

Remark 7.5**.**

It should be noted that associativity is often only approximately true in practice because of truncation errors, but this is negligible as soon as we use enough coefficients for the truncation error to be small (see e.g. [19, Section 4.5.1.2]).

We point out that, even for systems with nonlinear terms of low order, these considerations are also relevant if one wishes to compute higher order moments associated to PC expansions. Indeed the $k$-th moment of $u$ is given by

[TABLE]

where we have

[TABLE]

Finally, let us point out that handling non polynomial functions of PC expansions is also possible, although less straightforward, see e.g. [6].

7.3 Faster computations

The usual paradigm of PC is that, to minimize computational cost one should use the expansion basis associated to the PDF of the random inputs, as this minimize the number of coefficients needed to reach a given accuracy (see Section 4.2.1). However, the cost of computing the nonlinear terms for a given expansion should also be taken into account, and we briefly discuss it here. Let us consider truncated PC expansions of size $N$:

[TABLE]

and evaluate the cost of computing $(u\ast v)_{n}$ for all $0\leq n<N$, assuming the linearization coefficients $\alpha$ have been precomputed and stored. From (36) we have

[TABLE]

therefore, using directly this formula the total cost of computing $(u\ast v)_{k}$ for all $0\leq k<N$ is of order $N^{3}$. However, this cost can easily be reduced by one order for some specific expansions, for which we know a priori that most of the linearization coefficients are zero. In particular, for Taylor expansions the convolution product (usually called Cauchy product in that case) writes

[TABLE]

and thus the total cost of computing $(u\ast v)_{k}$ for all $0\leq k<N$ is of order $N^{2}$. Chebyshev expansions of the first kind also enjoy a similar property. Indeed, let us consider $\phi_{n}$ defined as $\phi_{0}=T_{0}$ and $\phi_{n}=2T_{n}$ for all $n\geq 1$. Then the associated convolution product is related to the classical discrete convolution product and writes

[TABLE]

and the total cost of computing $(u\ast v)_{k}$ for all $0\leq k<N$ is again of order $N^{2}$. Finally, when $N$ becomes large one may want to try and reduce this cost even further. A natural way to do so for Chebyshev or Taylor expansions is to compute the convolution via a Fast Fourier Transform, which then brings down the cost even lower, to an order of $N\log N$ (see e.g. [23, 26]). This idea was partially extended to other bases such as the Gegenbauer polynomials, see [25] and the references therein for more details.

Acknowledgments:

MB and CK have been supported by a Lichtenberg Professorship of the VolkswagenStiftung.