A set-oriented path following method for the approximation of parameter dependent attractors

TL;DR

This paper introduces a set-oriented path following technique for efficiently approximating parameter-dependent attractors in dynamical systems, demonstrated on shear flow models transitioning to turbulence.

Contribution

It presents a novel method that tracks relative global attractors across parameters without repeated subdivision, improving computational efficiency.

Findings

Successfully applied to shear flow models transitioning to turbulence

Enables efficient tracking of attractors as parameters vary

Reduces computational effort compared to traditional methods

Abstract

In this work we present a set-oriented path following method for the computation of relative global attractors of parameter-dependent dynamical systems. We start with an initial approximation of the relative global attractor for a fixed parameter $\lambda_0$ computed by a set-oriented subdivision method. By using previously obtained approximations of the parameter-dependent relative global attractor we can track it with respect to a one-dimensional parameter $\lambda > \lambda_0$ without restarting the whole subdivision procedure. We illustrate the feasibility of the set-oriented path following method by exploring the dynamics in models for shear flows during the transition to turbulence.

Figures

Equations

$$\dot{x}=F(x,\lambda ),$$

$$F(x,\lambda )=0.$$

$$x_{j+1}=f(x_j,\lambda ),j=0,1,\dots ,$$

$$A_Q=\bigcap _{j\ge 0}{g}^j(Q).$$

$$\mathrm{diam}({\mathcal{B}}_k)=\underset{B\in {\mathcal{B}}_k}{\max }\mathrm{diam}(B)$$

$$\bigcup _{B\in {\hat{\mathcal{B}}}_k}B=\bigcup _{B\in {\mathcal{B}}_{k-1}}B$$

$$\mathrm{diam}({\hat{\mathcal{B}}}_k)=\theta \mathrm{diam}({\mathcal{B}}_{k-1}),$$

$${\mathcal{B}}_k=\left\{B\in {\hat{\mathcal{B}}}_k:\exists \hat{B}\in {\hat{\mathcal{B}}}_k\text{mbox}suchthat{g}^{-1}(B)\cap \hat{B}\ne \varnothing \right\}.$$

$$\operatorname{diam}({\mathcal{B}}_{k})\leq\theta^{k}\operatorname{diam}({\mathcal{B}}_{0})\rightarrow 0\quad\mbox{for $k\rightarrow\infty$.}$$

$${\mathcal{B}}_k=\left\{B\in {\hat{\mathcal{B}}}_k:\exists \hat{B}\in {\hat{\mathcal{B}}}_k\text{mbox}suchthat{g}^{-1}(B)\cap \hat{B}\ne \varnothing \right\}.$$

$$g^{-1}(B_i)\cap {\hat{B}}_j\ne \varnothing .$$

$$g(x)\notin B_i\text{mbox}foralltestpointsx\in {\hat{B}}_j.$$

$$Q_k=\bigcup _{B\in {\mathcal{B}}_k}{B}_{.}$$

$$Q_{\infty }=\bigcap _{k=0}^{\infty }{Q}_k$$

$$A_Q=Q_{\infty }.$$

$$A_{Q^1}=A=A_{Q^2},$$

$$\underset{k\to \infty }{\lim }h(Q_k^1,Q_k^2)=0,$$

$$\text{dist}\left(B,C\right)=\underset{x\in B}{\sup }\underset{y\in C}{\inf }\Vert x-y\Vert .$$

$$\underset{\lambda \to {\lambda }_0}{\lim }\text{dist}\left(A^{\lambda },A^{{\lambda }_0}\right)=0.$$

$$\underset{\lambda \to {\lambda }_0}{\lim }\text{dist}\left(A^{{\lambda }_0},A^{\lambda }\right)=0\forall \lambda \in \Lambda .$$

$$\text{dist}\left(A^{\lambda },A^{{\lambda }_0}\right)<\epsilon \forall \lambda \in U_{\delta }({\lambda }_0).$$

$$A^{\lambda }\subset U_{\epsilon }(A^{{\lambda }_0})\subset Q_k^{{\lambda }_0}\forall \lambda \in U_{\delta }({\lambda }_0).$$

$$\dot{x}$$

$$\dot{y}$$

$$\dot{z}$$

$$\text{dist}\left(A^{{\lambda }_1},A^{{\lambda }_0}\right)\le C\delta$$

$$A^{{\lambda }_1}\subset U_{C\delta }(A^{{\lambda }_0})\subset Q_K^{{\lambda }_0}$$

$$\frac{\text{mbox}d}{\text{mbox}dt}\left(\begin{array}{c}u\\ v\\ w\\ m\end{array}\right)=\frac{1}{R}\left(\begin{array}{c}0\\ 0\\ 0\\ \sigma \end{array}\right)-\frac{1}{R}\left(\begin{array}{c}\lambda u\\ \mu v\\ \nu w\\ \sigma m\end{array}\right)+\left(\begin{array}{cccc}0& 0& -\gamma w& v\\ 0& 0& \delta w& 0\\ \gamma w& -\delta w& 0& 0\\ -v& 0& 0& 0\end{array}\right)\left(\begin{array}{c}u\\ v\\ w\\ m\end{array}\right),$$

Full text

1 Introduction

Over the last two decades so-called set-oriented numerical methods have been developed in the context of the numerical treatment of finite dimensional dynamical systems (e.g., [12, 13, 14, 27, 29]). Here, the basic idea is to cover the objects of interest such as attractors, invariant manifolds or almost invariant sets by outer approximations which are created via subdivision and continuation techniques. The numerical effort depends essentially on the dimension of the global attractor, i.e., it is easier to compute a one-dimensional attractor in a ten-dimensional space than to compute a three-dimensional attractor in a four-dimensional space [13]. The set-oriented techniques have been used successfully in several different application areas such as molecular dynamics ([19, 46, 10]), astrodynamics ([17, 16, 15]) and ocean dynamics ([28]). Moreover, a set-oriented numerical framework has recently been developed which allows to perform uncertainty quantification for dynamical systems from a global point of view [18] and the computation of finite dimensional attractors or unstable manifolds of infinite dimensional dynamical systems [11, 52, 51]. In this work we connect these set-oriented algorithms with ideas from bifurcation analysis and path following methods in order to efficiently treat parameter dependencies in the underlying dynamical system.

Systems of physical interest typically depend on parameters appearing in the defining systems of equations, i.e., one considers dynamical systems of the form

[TABLE]

where $x\in\mathbb{R}^{n}$ are the variables, $\lambda\in\mathbb{R}^{p}$ are the system parameters, and $F:\mathbb{R}^{n}\times\mathbb{R}^{p}\to\mathbb{R}^{n}$ is assumed to be sufficiently smooth. By varying $\lambda$, the qualitative structure of the solutions might change significantly, e.g., stable equilibria of the dynamical system might become unstable. Such phenomena are known as bifurcations and the parameter values are called bifurcation parameters (see [34, 36, 32, 49] for a detailed overview).

Bifurcations typically come in two kinds. On the one hand, one can look at local bifurcations [34, 32], where the analysis of such bifurcations is generally performed by studying the corresponding vector fields near equilibrium points. That is, one can usually reduce the problem to the analysis of an equation of the form

[TABLE]

Starting with an initial $\lambda_{0}$ the aim is to find numerically an approximation of a solution curve depending on $\lambda$ which is implicitly defined by (2). The basic idea of path following methods which tackle this kind of problem for parameter-dependent fixed points and steady states can be found, e.g., in [1] and [2]. For the numerical treatment of bifurcation problems there exist a wide range of software tools, e.g., the well-known software package AUTO [21]. A first-version of AUTO has already been developed in 1981. More recent contributions to bifurcation analysis software are the MATLAB packages MATCONT [20] and COCO [7].

On the other hand, there are global bifurcations [31, 49] for which it is not sufficient to reduce the study to a neighborhood of an equilibrium or a fixed point. For example, these bifurcations occur when invariant sets, such as periodic orbits, collide with equilibria. Therefore, typical examples are formations of homoclinic and heteroclinic orbits which can also be found numerically with the software package AUTO [21] and COCO by formulating and solving an appropriate boundary value problem (e.g., [22]). In addition, MATCONT is able to continue homoclinic and heteroclinic orbits by a homotopy method [8].

In this work we extend the subdivision method for the approximation of the so–called relative global attractor to parameter-dependent dynamical systems, i.e., we design a path following method from a global point of view, which indirectly allows the numerical analysis of global bifurcations. Note that the relative global attractor contains all backward invariant sets and, in particular, the homoclinic or heteroclinic connections. Using previously obtained approximations of the relative global attractor we can track it with respect to a one-dimensional parameter without restarting the whole subdivision procedure. Since one computation of an attracting set might be very expensive by using the set-oriented method, the algorithm developed in this article can reduce the overall computational cost significantly. The set-oriented path following algorithm has been implemented in the dynamical systems software package GAIO (Global Analysis of Invariant Objects) [9], which is available for MATLAB (see https://github.com/gaioguy/GAIO).

We will use the set-oriented path following method to explore the dynamics of shear flows during the transition to turbulence. In contrast to convection or Taylor-Couette flows, where the transition can be captured by sequences of bifurcations that branch off the laminar state [4, 5], shear flows such as pipe flow or plane Couette flow do not show a linear instability of the laminar profile [45, 33]. The transition is then mediated by the presence of 3d states that form in saddle-node bifurcations [44, 24, 23, 25]. Very close to the bifurcation, the states of the upper branches can be stable, but they soon undergo secondary bifurcations and period doubling cascades [39, 3, 50]. Crisis bifurcations destroy the attractors and give rise to the open saddles that underlie the observed transience of the turbulent state [38, 40].

Exploring the dynamics in the full phase space is very complicated because of the high-dimensionality of the relevant subspaces, though notable progress has been made in some cases [35, 30]. In order to prepare for the application to such computationally expensive examples, we begin here with an analysis of simple models that capture much of the phenomenology of shear flows: they have a laminar state that is linearly stable and a saddle-node bifurcation that represents the transition to turbulence. The models even contain the transition from an attractor to a transient saddle.

The simplest model where the modes have a physical interpretation is the four-dimensional model proposed by Waleffe [47]: the four components represent the transversal velocity components, the vortices, the streak, and the mean velocity. With the methods described here, we are able to trace the appearance of the attractor and its growth as the Reynolds number increases. We will then extend the analysis to the case of the nine-dimensional model discussed in [42, 43].

A detailed outline of the paper is as follows: in Section 2.1 we briefly summarize the subdivision method introduced in [13]. In Section 2.2 we then develop the set-oriented path following method which allows us to compute attracting sets of parameter-dependent dynamical systems without restarting the subdivision procedure. In Section 3 we illustrate the efficiency of the path following method for two different reduced order models from fluid dynamics, namely, a four-dimensional model of self-sustained flows and a nine-dimensional model of turbulent shear flows. Conclusions are given in Section 4.

2 A set-oriented path following method

We consider parameter dependent dynamical systems of the form

[TABLE]

where $x_{j}\in\mathbb{R}^{n}$ and $f:\mathbb{R}^{n}\times\Lambda\rightarrow\mathbb{R}^{n}$ is a homeomorphism for each parameter $\lambda$ in an open subset $\Lambda\subset\mathbb{R}$ and uniformly continuous in $\lambda$ on bounded sets, e.g., $f$ could be a time-$T$-map for the system (1). In what follows, for fixed $\lambda\in\Lambda$ we use the abbreviation $f_{\lambda}(x):=f(x,\lambda)$. Moreover, we assume that $f_{\lambda}$ has a compact global attractor $A^{\lambda}\subset\mathbb{R}^{n}$ for each $\lambda\in\Lambda$, that is, the compact set $A^{\lambda}\subset\mathbb{R}^{n}$ attracts any bounded set $B\subset\mathbb{R}^{n}$. Later on, we will additionally assume that $A^{\lambda}$ is upper semi-continuous in $\lambda$.

The aim of this work is to develop a set-oriented path following method for the approximation of $A^{\lambda}$ depending on the parameter $\lambda$. Since this method is based on the framework developed in [13] we start with a short review of the related material.

2.1 Brief review of the subdivision method

Let $Q\subset\mathbb{R}^{k}$ be a compact set. For a homeomorphism $g:\mathbb{R}^{n}\to\mathbb{R}^{n}$ we define the global attractor relative to $Q$ by

[TABLE]

Observe that $A_{Q}=A$ if $A\subset Q$ is the compact global attractor of $g$.

As a first step we approximate this set with a subdivision procedure introduced in [13] for a fixed $\lambda_{0}\in\Lambda$, that is, $g=f_{\lambda_{0}}$. Roughly speaking, the idea of the algorithm is to start with a finite family ${\mathcal{B}}_{0}$ of compact subsets of $\mathbb{R}^{n}$ that are a partition of $Q$. Then we refine this partition and remove those subsets that do not contain parts of the relative global attractor. By continuing this process we generate a sequence ${\mathcal{B}}_{0},{\mathcal{B}}_{1},\ldots$ of finite collections of compact subsets of $\mathbb{R}^{n}$ such that the diameter

[TABLE]

converges to zero and the union $\bigcup_{B\in B_{k}}B$ approaches the relative global attractor $A_{Q}$ for $k\rightarrow\infty$.

More precisely, let ${\mathcal{B}}_{0}$ be the initial collection, e.g., ${\mathcal{B}}_{0}=\{Q\}$. Then we inductively obtain ${\mathcal{B}}_{k}$ from ${\mathcal{B}}_{k-1}$ for $k=1,2,\ldots$ in two steps.

Subdivision: Construct a new collection $\widehat{\mathcal{B}}_{k}$ such that

[TABLE]

and

[TABLE]

where $0<\theta<1$. 2. 2.

Selection: Define the new collection ${\mathcal{B}}_{k}$ by

[TABLE]

The first step guarantees that the collections ${\mathcal{B}}_{k}$ consist of successively finer sets for increasing $k$. In fact, by construction

[TABLE]

In the second step we remove each subset whose preimage does neither intersect itself nor any other subset in $\widehat{\mathcal{B}}_{k}$.

Remark 2.1**.**

In the application of the subdivision scheme for the computation of the relative global attractor we have to perform the selection step

[TABLE]

Thus, we have to decide whether or not the preimage of a given set $B_{i}\in{\mathcal{B}}_{k}$ has a nonzero intersection with another set $B_{j}\in\widehat{{\mathcal{B}}_{k}}$, i.e.,

[TABLE]

Numerically this is realized as follows: we discretize each box $\widehat{B}_{j}$ by a finite set of test points $x\in\widehat{B}_{j}$ and replace the condition (5) by

[TABLE]

The selection step is responsible for the fact that the unions $\bigcup_{B\in{\mathcal{B}}_{k}}B$ approach the relative global attractor. Denote by $Q_{k}$ the collection of compact sets obtained after $k$ subdivision steps, that is

[TABLE]

Note that these sets form a nested sequence $Q_{k+1}\subset Q_{k}$. Therefore the limit

[TABLE]

exists and the following result has been proven in [13]:

Proposition 2.2**.**

Let $A_{Q}$ be the global attractor relative to $Q$, and let ${\mathcal{B}}_{0}$ be a finite collection of closed subset with $Q_{0}=\bigcup_{B\in{\mathcal{B}}_{0}}=Q$. Then the sets $Q_{k}$ obtained by the subdivision algorithm contain the relative global attractor, $A_{Q}\subset Q_{k}$ for all $k=0,1,\ldots$, and moreover

[TABLE]

2.2 Path following method for the approximation of relative global attractors

In this section we develop a path following algorithm that allows to compute the relative global attractor for various parameter values $\lambda\in\Lambda$ of (3) by reusing previously obtained coverings. The idea of this method is to first approximate the relative global attractor $A_{Q}^{\lambda_{0}}$ for $\lambda_{0}\in\Lambda$ and then to use a covering of this set, denoted by $Q_{k}^{\lambda_{0}}$, as a starting point to compute the relative global attractor $A_{Q}^{\lambda_{1}}$ for $\lambda_{1}\in\Lambda$ sufficiently close to $\lambda_{0}$.

We start with the following observation.

Lemma 2.3**.**

Let $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ be a homeomorphism and assume that $f$ has a compact global attractor $A\subset\mathbb{R}^{n}$. Furthermore, let $Q^{1},Q^{2}$ be compact sets such that $A\subset Q^{1}\subset Q^{2}$. Then,

[TABLE]

where $A_{Q^{i}}$ is the global attractor relative to $Q^{i}$ for $i=1,2$.

Proof.

This statement follows directly by the definition of the global attractor relative to a compact set $Q$ (cf. (4)) and the fact that $A\subset Q^{1}\subset Q^{2}$. ∎

As a direct consequence we obtain the following result.

Proposition 2.4**.**

Let $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ be a homeomorphism and assume that $f$ has a compact global attractor $A\subset\mathbb{R}^{n}$. Let $Q^{1},Q^{2}$ be compact sets such that $A\subset Q^{1}\subset Q^{2}$ and denote by $Q_{k}^{i}$, $i=1,2$, the corresponding sets obtained after $k$ subdivision and selection steps (see (6)). Then

[TABLE]

where $h(B,C)$ denotes the usual Hausdorff distance between two compact subsets $B,C\subset\mathbb{R}^{n}$.

Proof.

This result is a direct consequence of Proposition 2.2 and Lemma 2.3. ∎

Proposition 2.4 states that by performing sufficiently many subdivision and selection steps we obtain a good approximation of the attracting set $A$ no matter how we choose the initial covering $Q$ as long as $A\subset Q$.

In what follows, let us assume that $Q\subset\mathbb{R}^{n}$ is a compact set such that $A^{\lambda}\subset Q$ for all $\lambda\in\Lambda$. We denote by $\mathcal{Q}(\lambda_{0})$ a finite family of compact sets $Q_{k}^{\lambda_{0}}$, $k=0,1,\ldots,m$, which all contain the global attractor $A_{Q}^{\lambda_{0}}$ relative to $Q$ of the map $f_{\lambda_{0}}$. Next, either we fix $\lambda_{1}\in\Lambda$ and choose $K\leq m$ sufficiently small such that $A^{\lambda_{1}}\subset Q_{K}^{\lambda_{0}}$ or we fix $K\leq m$ and choose $\lambda_{1}\in\Lambda$ sufficiently close to $\lambda_{0}$ such that $A^{\lambda_{1}}\subset Q_{K}^{\lambda_{0}}$. In any case Lemma 2.3 and Proposition 2.4 allows us to approximate the global attractor $A_{Q}^{\lambda_{1}}$ of the map $f_{\lambda_{1}}$ using the initial covering $Q_{K}^{\lambda_{0}}$ with the subdivision algorithm described in the previous section. Observe that we can always choose $K=0$ for every $\lambda_{1}$, since $Q_{0}^{\lambda_{0}}=Q$ and $A_{Q}^{\lambda_{1}}\subset Q$ for $\lambda_{1}\in\Lambda$ by assumption.

To discuss the choice of $\lambda_{1}$ for fixed $K$ we have to analyze the change of the attractor $A^{\lambda}$ when the parameter $\lambda$ is varied. To this end, we define the distance between two subsets $B,C$ of $\mathbb{R}^{n}$ as

[TABLE]

This distance allows us to define the continuity of attractors as follows:

Definition 2.5**.**

Let $\lambda_{0}\in\Lambda$. The family of attractors $A^{\lambda}$ is upper semi-continuous at $\lambda_{0}$ if

[TABLE]

The attractor $A^{\lambda}$ is called upper semi-continuous if $A^{\lambda}$ is upper semi-continuous for each $\lambda\in\Lambda$.

Note that upper semi-continuity at $\lambda_{0}$ implies that for every $\varepsilon>0$ there exists a neighborhood $U_{\delta}(\lambda_{0})\subset\Lambda$ of $\lambda_{0}$ such that $A^{\lambda}\subset U_{\varepsilon}(A^{\lambda_{0}})$ for all $\lambda\in U_{\delta}(\lambda_{0})$, where $U_{\varepsilon}(A^{\lambda_{0}})\subset\mathbb{R}^{n}$ denotes the $\varepsilon$-neighborhood of $A^{\lambda_{0}}$. Thus, the attractor cannot become larger instantaneously by varying $\lambda$ which is a naturally needed property for our proposed path following scheme. In fact, if the attracting set suddenly “explodes” we cannot assure that any previously computed covering $Q_{k}^{\lambda}$ besides $k=0$ still covers $A^{\lambda_{1}}$. For the sake of completeness we note that lower semi-continuity of $A^{\lambda}$, i.e.,

[TABLE]

prevents a sudden shrinking of $A^{\lambda}$ which holds for instance for gradient systems with hyperbolic fixed points [37]. However, this would be, in general, no issue for our proposed algorithm.

We are now in the position to prove that for each $K\in\{0,\ldots,m\}$ there is a range of parameters $\lambda\in\Lambda$ such that the attractor $A^{\lambda}$ can be approximated by the subdivision algorithm with the initial compact set $Q_{K}^{\lambda_{0}}$.

Proposition 2.6**.**

Let $\mathcal{Q}(\lambda_{0})=\{Q_{0}^{\lambda_{0}},\ldots,Q_{m}^{\lambda_{0}}\}$ be the family of sets generated by the subdivision method such that $A_{Q}^{\lambda_{0}}=A^{\lambda_{0}}\subset Q_{k}^{\lambda_{0}}$ for all $k=0,\ldots,m$. Suppose that there is $\varepsilon>0$ such that $U_{\varepsilon}(A^{\lambda_{0}})\subset Q_{k}^{\lambda_{0}}$ for all $k=0,\ldots,m$ and $A^{\lambda}$ is upper semi-continuous at $\lambda_{0}$. Then for every $k=0,\ldots,m$ there exists $\delta=\delta(k)>0$ such that $A^{\lambda}\subset Q_{k}^{\lambda_{0}}$ for all $\lambda\in U_{\delta}(\lambda_{0})$. In particular, $A^{\lambda}$ can be approximated by using the initial compact set $Q_{k}^{\lambda_{0}}$.

Proof.

Let $k\in\{0,\ldots,m\}$ be fixed. Due to the upper semi-continuity of $A^{\lambda}$ in $\lambda_{0}$ there is $\delta>0$ such that

[TABLE]

Thus, we conclude by assumption that

[TABLE]

∎

In the following we will call such $\lambda\in U_{\delta}(\lambda_{0})$ feasible and Proposition 2.6 guarantees the existence of a feasible $\lambda\in\Lambda$. Throughout the remainder of this article we now suppose that $A^{\lambda}$ is upper semi-continuous and the additional assumptions of Proposition 2.6 hold. Then the numerical realization of the set-oriented path following algorithm can be described as follows:

In Figure 1, we illustrate two steps of Algorithm LABEL:alg:pathFollowing for the Lorenz system [41] given by

[TABLE]

where we use $\beta$ as our parameter of interest. We note that the derivation of the Lorenz model shows that $\beta$ is related to the aspect ratio of the convection cells.

Remark 2.7**.**

1. (a)

Instead of following the path along $\lambda$ in positive direction one can also take $\lambda_{j+1}<\lambda_{j}$ and follow the path in the negative direction. 2. (b)

Intuitively, choosing a larger $K_{j+1}$ decreases the range of feasible parameters $\lambda_{j+1}$. However, one has to perform less subdivision and selection steps to reach the same level, i.e., the same fineness of the final box covering of $A^{\lambda_{j+1}}$. With more knowledge on the upper semi-continuity property (7) we can make this more precise in the following Proposition. 3. (c)

If $\lambda_{j+1}$ is not feasible, i.e., $\lambda_{j+1}\notin U_{\delta}(\lambda_{j})$, $Q_{K_{j}}^{\lambda_{j}}$ might not contain all of $A^{\lambda_{j+1}}$ and thus, the subdivision algorithm does not approximate the whole attractor $A^{\lambda_{j+1}}$. However, to make the algorithm more robust, we reintroduce sets if an image $f_{\lambda_{j+1}}(x)$ of a test point $x\in\mathbb{R}^{n}$ is not contained in any set of the current collection ${\mathcal{B}}_{k}(\lambda_{j})$.

After discussing the choice of $\lambda_{1}\in\Lambda$ for fixed $K\in\{0,\ldots,m\}$ the following Proposition tells us, how to choose $K_{j+1}$ for a fixed step size in the parameter space $\Lambda$, i.e., how to choose $K\in\{0,\ldots,m\}$ when $\lambda_{1}\in\Lambda$ is given.

Proposition 2.8**.**

Let $\delta>0$ and $\lambda_{1}\in\Lambda$ with $\left\lvert\lambda_{1}-\lambda_{0}\right\rvert\leq\delta$. Suppose there is a constant $C>0$ such that

[TABLE]

and $U_{C\delta}(A^{\lambda_{0}})\subset Q_{K}^{\lambda_{0}}$ for one $K\in\{0,\ldots,m\}$. Then $A^{\lambda_{1}}\subset Q_{K}^{\lambda_{0}}$ and, in particular, $\lambda_{1}$ is feasible for $K$.

Proof.

According to (9) the $C\delta$-neighborhood $U_{C\delta}(A^{\lambda_{0}})$ of $A^{\lambda_{0}}$ contains $A^{\lambda_{1}}$ and by assumption we immediately obtain

[TABLE]

as claimed. ∎

3 Numerical examples

In this section we present results of computations carried out for a four- and a nine-dimensional reduced order model from fluid dynamics.

3.1 A four-dimensional model of self-sustained flows

As a first example, we consider a four-dimensional nonlinear model of self-sustained flows introduced in [47] and further discussed in [48]. The four variables of the model represent the spanwise velocity components ($u$), the vortices ($v$), the streak ($w$), and the mean profile ($m$) (cf. [48] for more explanations of the model). Their dynamics is given by

[TABLE]

where $R>0$ is the Reynolds number, and $\lambda$, $\mu$, $\nu$, $\sigma,\gamma$, and $\delta$ are positive parameters. For the parameters, we take values $\lambda=\mu=\sigma=10$, $\nu=15$, $\delta=1$ and $\gamma=0.5$; a saddle-node bifurcation then appears for $R_{c}=98.6325$.

Our goal is to numerically analyze how the attracting sets change for the Reynolds numbers $R\in[98,400]$. We define $x=(u,v,w,m)^{\top}$ and denote by $f(x,R)$ the time-$T$-map of (10) which depends on the Reynolds number $R$ and the integration time $T$. For the following analysis we fix the integration time of $T=20$. Furthermore, we choose the initial box $Q=[-0.9,1.1]\times[-0.8,1.2]\times[-1,1]\times[-0.8,1.2]$ in which we want to approximate the relative global attractor $A_{Q}^{R}$. According to Algorithm LABEL:alg:pathFollowing we fix $m=36$ and $K=32$. Following of Remark 2.7 (a) we start the path following algorithm for the Reynolds number $R_{0}=400$ and define $R_{j+1}=R_{j}-1$ for $j=0,\ldots,302$ such that the interval $\Lambda=[98,400]$ is discretized equidistantly. A detailed analysis of (10) for $R\in\Lambda$ can be found in [48]. In what follows, we will study the attractor from a global point of view. To this end, Figure 2 and Figure 4 show three-dimensional projections of the relative global attractor obtained by Algorithm LABEL:alg:pathFollowing for different Reynolds numbers.

For a parameter value just below the saddle-node bifurcation, i.e., for $R=98$, the attracting set does not contain any invariant structures besides the laminar profile $(u,v,w,m)=(0,0,0,1)$. For $R$ just above the saddle node bifurcation, i.e., for $R=98.6325$, Figure 2 (a) shows that the attractor now contains an upper branch steady solution, as well as the lower branch saddle state. This situation persists up to $R=100.0232$, after which the upper branch undergoes a supercritical Hopf bifurcation.

The limit cycle is stable (cf. Figure 2 (b)) until $R=101.0311$ (cf. Figure 2 (c)), where it disappears in a homoclinic bifurcation. Note that the attracting set does not change significantly, although the dynamics does: in the homoclinic bifurcation the attractor rips open and becomes a transient saddle. In order to demonstrate this change in the dynamics, we compare in Figure 3 (a) and (b) the average lifetimes for each box in the attracting set for $R=101$ and $R=102$, respectively.

For $R=101$ we see a strict separation in the attracting set by the edge state: all boxes between the edge state and the laminar state have finite lifetimes, as they are attracted towards the laminar profile. All boxes inside the lobe around the upper branch state have infinite lifetimes because they are attracted to the limit cycle and cannot return to the laminar profile. Here, test points in each box were only integrated for times up to $t=300$: the choice of time is not essential, since for a closed attractor around the upper branch state no trajectory can return to the laminar flow. For $R=102$ the limit cycle is not stable anymore, and all points except for the fixed points and the limit cycle return to the laminar profile. Observe that the obtained box coverings do not only cover the laminar profile but also the unstable manifolds which is due to the fact that the relative global attractor contains all backward invariant sets.

For $102<R\leq 356$ the attractor develops a more complicated structure with folds and other features (Figure 4 (a–d)), that arise from the projection from the four-dimensional state space to the three-dimensional image plane.

Finally, for $356<R<435$ another stable limit cycle appears (cf. Figure 4 (e) and (f)).

3.2 A nine-dimensional model for turbulent shear flows

As a second example we consider a Fourier based nine-dimensional model for a turbulent shear flow with free slip boundaries and a sinusoidal base profile. A first version with 8 degrees of freedom was described by Waleffe [47, 48], and extended to 9 modes in [42, 43]. The model is based on a Fourier expansion of the velocity field, and every component can be assigned to a specific flow field. In addition to the Reynolds number, the model has two spatial parameters describing the domain size in spanwise and downstream direction. The choice we consider here is the NBC domain (Nagata, Busse, Clever [44, 6]), i.e., a box with $L_{z}=2\pi$ in spanwise and $L_{x}=4\pi$ in downstream direction. For these parameters, the first state that appears is a periodic orbit, at $R=89.76$. The periodic orbit is unstable from the onset, so that the dynamics is transient without the need for the occurrence of a crisis bifurcation. However, the lifetime increases very rapidly with Reynolds number [42].

Again, we choose the Reynolds number $R$ as our parameter of interest and compute the relative global attractor $A_{Q}^{R_{j}}$ for Reynolds numbers starting with an initial $R_{0}=200$ and then decreasing in steps of one, $R_{j+1}=R_{j}-1$, ${j=0,\ldots,100}$. Due to the complexity of the dynamics in the nine-dimensional phase space it is hard to visualize the attracting set. An impression of the set may be obtained by three-dimensional projections onto the modes $(a_{1},a_{4},a_{3})$ that correspond to $(m,u,v)$ in the four-dimensional model (10) and to the modes $(a_{1},a_{6},a_{2})$ which represent the principal components of the set.

In Figure 5 we see that these projections of the attractor look like dense balls with no discernable interior structure. A measure of the increasing complexity of the dynamics is the dimension of the attracting set, here quantified by the box-counting dimension (cf. [26]) of the box coverings obtained by Algorithm LABEL:alg:pathFollowing. Figure 6 shows that this dimension increases with Reynolds number. For $R=100$ the dimension of the attracting set is about $d=2.6$ and increases to about $3.1$ for $R=200$ and most likely even higher values for $R>200$.

4 Conclusions

The set-oriented path following method is a powerful tool for the numerical analysis of global bifurcations for attracting sets. With the method developed in this paper the actual covering of an attracting set for a parameter value $\lambda$ can be used as an initial covering for a parameter $\bar{\lambda}$ which is sufficiently close to $\lambda$. Thus, we do not need to restart the subdivision method, thereby significantly reducing the overall computational time. However, this path following method is only defined for finite dimensional dynamical systems. Therefore, it should be desirable to extend this work to the infinite dimensional setting, i.e., for the parameter dependent computation of attractors of infinite dimensional dynamical systems. In this context the results of [11, 52, 51] would have to be adapted.

The results in this article show that set-oriented methods can also be used in moderate-dimensional settings, where they can provide insight into the dynamics of shear flows during the transition to turbulence in simple models. The similarities in the behaviour of the low-dimensional models and fully resolved numerical simulations suggests that the global organization of phase space will be similar, even though such a study is computationally out of reach for high-dimensional attractors.

Acknowledgments

This work is supported by the Priority Programme SPP 1881 Turbulent Superstructures of the Deutsche Forschungsgemeinschaft.