Periodic orbit theory of chaos

About the lectures

These are notes for four lectures to be given in the school of mathematics of the university of Bristol in May 2015. I am grateful to the University of Bristol and especially Jonathan Robbins for the warm hospitality I received during my stay in Bristol, and to the Leverhulme Trust for generous financial support.

It is well known that many dynamical systems exhibit chaotic behaior, which renders their long-time dynamics inherently unpredictable. For chaotic systems it is therefore reasonable to abandon the idea of describing their dynamics exactly, and focus instead on probabilistic descrption, even when the dynamical laws are deterministic.

However, deriving the statistical distribution from the dynamical laws is a difficult problem that has been solved only in special cases. The goal of these lectures is to introduce periodic orbit theory, the strongest method that has been found so far to approach this problem.

The limited scope of the lectures allows me to introduce only the fundamental ideas of the method, along with the necessary basic ideas from dynamical systems theory (and only those). The level of rigour is dictated by my scientific background, and is similar to what one find in theoretical physics literature.

The main source for the material presented here is the free webbook “Chaos: Classical and Quantum” by P. Cvitanovic et al., available at  chaosbook.org . Other sources will be published together with these notes and other course material on the course website TBD.

Course plan

The Lorenz system: Introduction to chaos

The Lorenz system consists of 3 ordinary differential equations for the variables x,y,z. It was put forward by Lorenz in 1963 as a truncation of climate model equations exhibiting aperiodic dynamics.

x'(t)==σ (y(t)-x(t))
y'(t)==r x(t)-x(t) z(t)-y(t)
z'(t)==x(t) y(t)-b z(t)

σ,r, and b are real constant parameters. The values assumed by Lorenz σ=10, b=(8/3), and r=28, lead to strongly chaotic dynamics, and are also the most studied. We will use these parameter values exclusively here.

A typical orbit, that is, a solution (x(t),y(t),z(t)) of the Lorenz equations looks like this:

Module[{σ=10,b=(8/3),r=28,x0=5,y0=-5,z0=3,td=100,lx,ly,lz,gr,nseg},{lx,ly,lz}=NDSolveValue[{x'[t]==σ (-x[t]+y[t]),y'[t]==r x[t]-y[t]-x[t] z[t],z'[t]==x[t] y[t]-b z[t], x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,td}];gr=ParametricPlot3D[{lx[t],ly[t],lz[t]},{t,0,td},PlotRange->{{-30,30},{-30,30},{0,50}},PlotPoints->200,PlotStyle->Thin]; gr/.Line[l_]:>(nseg=Length[l]);Animate[Show[gr/.Line[l_]:>Line[Take[l,n]]],{n,1,nseg,1},AnimationRate->300,AnimationRepetitions->Infinity]]

The dynamics can be roughly described as follows: After an initial transient, the orbit is attracted to a butterfly shaped region in the three-dimensional state space. Each wing of the butterfly looks like a slightly curved membrane, but is has a zero-measure transverse component—the attractor is a fractal with dimension slight larger than 2.

The dynamics in each wing follows an ourward spiral—clockwise in the left wing and counter-clockwise in the right, when viewed from above—and leaves the wing when the spiral radius increases beyond a certain radius, landing at some point in the other wing. The process then repeats itself. The wings are related by the exact symmetry (x,y,z)→(-x,-y,z) of the Lorenz equations.

Although the orbits converge in directions transverse to the attractor, they diverge tangent to it, a phenomenon called sensitive dependence on initial conditions, or ‘the butterfly effect.’

Module[{σ=10,b=(8/3),r=28,x0=5,y0=-5,z0=3,td=40,lx,ly,lz,gr,lxa,lya,lza,dz=0.01},{lx,ly,lz}=NDSolveValue[{x'[t]==σ (-x[t]+y[t]),y'[t]==r x[t]-y[t]-x[t] z[t],z'[t]==x[t] y[t]-b z[t], x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,td}];
{lxa,lya,lza}=NDSolveValue[{x'[t]==σ (-x[t]+y[t]),y'[t]==r x[t]-y[t]-x[t] z[t],z'[t]==x[t] y[t]-b z[t], x[0]==x0,y[0]==y0,z[0]==z0+dz},{x,y,z},{t,0,td}];gr=ParametricPlot3D[{{lx[t],ly[t],lz[t]},{lxa[t],lya[t],lza[t]}},{t,0,td},PlotRange->{{-30,30},{-30,30},{0,50}},PlotPoints->100,PlotStyle->Thin]; gr/.Line[l_]:>(nseg=Length[l]);Animate[Show[gr/.Line[l_]:>Line[Take[l,n]]],{n,1,nseg,1},AnimationRate->300,AnimationRepetitions->1]]

This means that two orbits with small initial separation δ(0) (in the example above δ(0)=0.01in the z direction), grows exponentially δ(t)=(e^(λt))δ(0) as long as δ remains small, with a Lyapunov exponent λ>0. At some point, depending logarithmically on the smallness of the initial perturbation, the orbit separation becomes large enough that the orbits split into different wings, and then any remaining correlation between them is lost.

The simultaneous stretching and compression mixes the states on the attractor very efficiently, so that a small piece of it is quickly spread over the entire butterfly, although its volume shrinks very fast, the dynamics being dissipative. The stretching and folding, together with sensitive dependence on initial values are the fundamental properties of chaotic dynamics.

The chaotic properties of the dynamics have an important consequence. In a coarse-grained description, the Lorenz dynamics consists of transition between two states, left, and right which occur after a sequence of switching times (Subscript[t, 1]),(Subscript[t, 2]),… The state of the system can be labelled for example by the sign of x+y.

Module[{σ=10,b=(8/3),r=28,x0=5,y0=-5,z0=5,td=100,lx,ly,lz},{lx,ly,lz}=NDSolveValue[{x'[t]==σ (-x[t]+y[t]),y'[t]==r x[t]-y[t]-x[t] z[t],z'[t]==x[t] y[t]-b z[t],x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,td}];Plot[lx[t],{t,0,td},PlotRange->All]]

A typical sequence of switching times looks like this

Module[{σ=10,b=(8/3),r=28,x0=5,y0=-5,z0=5,td=1000,lx,ly,lz,tintervals,tx},{lx,ly,lz}=NDSolveValue[{x'[t]==σ (-x[t]+y[t]),y'[t]==r x[t]-y[t]-x[t] z[t],z'[t]==x[t] y[t]-b z[t],x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,td}];tintervals=Select[Range[(1/2),1000,(1/2)],Negative[lx[#-(1/2)]lx[#]]&];tx=t/.(First@FindRoot[lx[t],{t,#-(1/2),#}])&/@tintervals;ListPlot[Drop[tx-RotateRight[tx],1]]]

Since the underlying dynamics is deterministic, one could have guessed that the switching time sequence obeys well-defined rules, but this is not the case. Sensitive dependence on initial values implies that there is no way to predict the future (Subscript[t, n]) given the past beyond a few steps because the initial value is only known to a finite precision. The inevitable initial ignorance makes the dynamics effectively random. The purpose of these lectures is to understand this random process.

Unimodal maps

. Discrete dynamical systems, the logistic map

The essence of the Lorenz dynamics is captured by the return map, which describes the discrete steps of climbing on the outward spirals and wing switching. Since the Lorenz attractor is almost confined to a two-dimensional manifold, the discrete dynamics takes place on a one-dimensioal manifold. The stretch and fold action implies that the mapping has at least one turning point.

Rather than attempt a direct description of the Lorenz dynamics, we will look at the simplest type of dynamical systems with these properties: unimodal maps, that is, piecewise smooth maps f(x) of the unit interval to itself with a single extremum, such as

Plot[Sin[π x](1-(x/2)),{x,0,1},PlotRange->{0,1}]

The dynamics of unimodal maps will be mostly illustrated by the logistic map family (Subscript[l, r])(x)=r x(1-x), parameterized by r. For 0<r ⩽4 it maps the unit interval into itself.

. Fixed points, bifurcations, periodic orbits

A simple way to visualize the dynamics of one-dimensional maps is to draw the graph of the map overlayed with the graph of the identity function. An orbit is then obtained by vertical and horizontal lines bouncing between the two graphs. For example, a typical orbit of (Subscript[l, 0.8]) looks like

Show[plotOrbit[lm[0.8,#]&,0.6,15],PlotRange->{0,0.4}]

Evidently, the orbit tend to zero, (lim_(n->∞))(Subscript[l, r]^n)(x)=0, (Subscript[l, r]^n) being the nth iterate of (Subscript[l, r]), for any initial x. In other words, x=0 is a globally attracting fixed point. It follows that it is a (linearly) stable fixed point, |(Subscript[l, r])'(0)|<1. Fixed point are visible in orbit plots as intersections between the map and identity graphs, and when f is increasing they are linearly stable iff the map slope is smaller than the unit slope.

x=0 is the only fixed point, and is globally attracting for all 0<r<1, but when r=1 there is a (transcritical) bifurcation, where it becomes unstable, and a new fixed point is created at (Subscript[x, 1])=1-(1/r). The new fixed point is stable and globally attracting for 1<r<3. A typical orbit for r=1.5 is

Show[plotOrbit[lm[1.5,#]&,0.9,15],PlotRange->{0,0.6}]

When r=(Subscript[r, 1])=3 another bifurcation occurs when the (Subscript[x, 1]) fixed point becomes unstable, but this time instead of a new fixed point, a period-2 periodic orbit  is created (Subscript[x, 2]^((1,2)))=((r+1+(Sqrt[(r-3) (r+1)]))/(2 r)). The orbit is stable if for small δ

|(Subscript[l, r]) ((Subscript[l, r]) (δ+(Subscript[x, 2])))-(Subscript[x, 2]) |<δ ⟺ |(Subscript[l, r])'((Subscript[x, 2]^((1))))(Subscript[l, r])'((Subscript[x, 2]^((2))))|<1

that is, for 3<r<1+(Sqrt[6])≈3.45. A typical orbit for r=3.1 is

Show[plotOrbit[lm[3.1,#]&,0.9,30],PlotRange->{0,0.8}]

The period doubling bifurcation can be visualised by plotting the orbit of (Subscript[l, r]^2), the second iterate of (Subscript[l, r])

Show[plotOrbit[Composition[lm[3.2,#]&,lm[3.2,#]&],0.9,15],PlotRange->{0,0.9}]

The fixed points of (Subscript[l, r]^2) consist of the unstable fixed points 0,(Subscript[x, 1]) of (Subscript[l, r]) and the points of the period-2 orbit (Subscript[x, 2]^((1,2))), which are stable as fixed points of (Subscript[l, r]^2).

When r=(Subscript[r, 2])=1+(Sqrt[6]) a second period doubling bifurcation takes place: The period 2 orbit becomes unstable, and a stable period 4 orbit is created.

Show[plotOrbit[lm[3.5,#]&,0.83,50],PlotRange->{0,1}]

Show[plotOrbit[Composition[lm[3.53,#]&,lm[3.53,#]&,lm[3.5,#]&,lm[3.5,#]&],0.9,15],PlotRange->{0,1}]

The period-4 orbit then becomes unstable at (Subscript[r, 3])≈3.544 with a third period doubling bifurcation, and the process continues indefinitely with decreasing spacings between bifurcations. After each bifurcation, in which a stable (2^n)-periodic orbit is created, the previous (2^k)-periodic orbits, k<n, remain and are all unstable. It can be shown that unimodal maps have at most one stable periodic orbit.

. Chaos, strange attractors, repellers

In the limit of the sequence of bifurcation r=(Subscript[r, ∞])=(lim_n) (Subscript[r, n])≈3.5699456 the system has infinitely many periodic orbits but none of them is stable. Instead, a typical orbit is aperiodic and approaches a “strange attractor” 𝒜, a (fractal) set of zero measure on which the dynamics. The orbit of a typical attractor point (Subscript[x, 𝒜]) and the set of periodic orbits are dense in the attractor, 𝒜=(OverBar[(⋂_k)(⋃_(n=k)^∞)(Subscript[l, r]^n)((Subscript[x, 𝒜]))])=(OverBar[⋃(Subscript[x, n])]).

When r>(Subscript[r, ∞]) (for most r) the attractor has positive measure, and the dynamics on it has sensitive dependence on initial conditions. The sensitivity can be quantified by the Lyapunov exponent λ =(lim_n)(1/n)(Subscript[Σ, n=1]^∞)log((Subscript[l, r])'((Subscript[l, r]^(n-1))(x)))>0. It follows from this and the density properties that the dynamics is chaotic. λ is independent of initial x except for a zero-measure set, namely the periodic points.

plotOrbits[lm[3.9,#]&,{0.6,0.6001},10,{Black,Red}]

When r=4 the logistic map is conjugate to the piecewise linear symmetric tent map: (Subscript[l, 4])◦ g=g ◦ (Subscript[t, 2]), (Subscript[t, r])=(r/2)(1-2|x|) with the bijection g(x)=(sin^2)(((πx)/2)). For both maps the image of the map and the attractor are the entire interval, and the dynamics is maximally chaotic, being conjugate to the Bernoulli shift map (see below).

When r>4 the image of the (Subscript[l, r])([0,1])⊄[0,1]. The closure ℛ of the set of (unstable) periodic orbits is still an invariant set, but is not an attractor: Points arbitrarily close to the invariant set are mapped outside the unit interval in a finite number of iterations. Rather ℛ is a repeller, which means that it is stabilised by the inverse map (Subscript[l, r]^(-1)) (as a map on sets). The probability (Subscript[p, n]) that a randomly chosen initial point survives the dynamics after n steps decays exponentailly: (Subscript[p, n])(Underscript[∼, n->∞])(e^(-γ n)) . γ is the escape rate from the repeller.

Module[{seed=20150503,r=RandomReal[1,9]},SeedRandom@seed;plotOrbits[lm[4.5,#]&,r,15,Hue/@(0.05lifeSpan[lm[4.5,#]&,#,100]&/@r)]]

(this is also true for some (Subscript[r, ∞])<r<4, except that repelled points are attracted to a stable periodic orbit with a period which is not a power of 2 instead of leaving the interval).

. Invariant measures

In the chaotic regime the finite precision of initial conditions makes it impossible to calculate actual future orbits indefinitely with any fixed degree of precision. Instead one has to adopt a probabilistic approach. Starting from a probability measure (Subscript[μ, 0]) of initial conditions concentrated near an initial value (Subscript[x, i]) (say) we would like to know the later probability measure (Subscript[μ, n]) such that (Subscript[μ, n])(A) is the probability that (f^n)(x)∈A.

It follows from the definition that (Subscript[μ, n+1])(A)=(Subscript[μ, n])((f^(-1))(A)). If we write the μs in terms of densities (Subscript[μ, n])(A)=(∫_A)(Subscript[ρ, n])(x)dx then

(Subscript[ρ, n+1])(x)=∫(Subscript[ρ, n])(ξ)δ(x-f(ξ))dξ=(∑_(k,f((Subscript[ξ, k]))=x))(((Subscript[ρ, n])((Subscript[ξ, k])))/(|f'((Subscript[ξ, k]))|))≡(Subscript[ℒρ, n])(x)

where ℒ is the evolution operator, also called Perron-Frobenius operator, and δ(x) is the Dirac delta ‘function’ (or distribution), defined by ∫h(x)δ(x)dx=h(0) for sufficiently nice h.

For dynamics approaching an attractor one may expect the measure to tend asymptotically to a limit measure that is invariant under ℒ. Actually, there are many invariant measures. For example, any fixed point and periodic orbit carries a singular invariant measure (1/p)(∑_(k=1)^p)δ(x-(Subscript[x, p]^((k)))). We are interested in the natural measure that is the one arrived from smooth initial measures. For systems with escaping orbits the natural measure is only conditionally invariant, that is ℒρ=(e^(-γ))ρ where γ is the escape rate as before.

Cycle expansion for the escape rate

The basic idea of cycle expansions is simple: Since periodic orbits are dense in the chaotic attractor/repeller, mean quantities can be calculated by averaging them on cycles with weights given by the relative time a typical orbit spends near each periodic orbit. In this section we apply this idea in the simplest setting of calculating escape rates, proceeding first formally, and then giving a heuristic justification.

The classical trace formula for maps

The survival ratio (e^(-γ)) is the largest eigenvalue of ℒ. To get access to the spectrum of ℒ we first express tr (ℒ^n) in terms of cycles

((trℒ)^n)=∫δ(x-(f^n)(x))dx=(∑_(k,(f^n)((Subscript[ξ, k]))=x))(1/(|1-(Subscript[Λ, p])|))=(∑_(n-cycles p))((Subscript[n, p])/(|1-(Subscript[Λ, p])|))

where (Subscript[(Subscript[Λ, p])=(d/dx)(f^n)(x)|, p])=f'((Subscript[x, p]^((1))))⋯f'((Subscript[x, p]^((n)))), and (Subscript[n, p]) is the number of distinct points in the cycle p. When the cycle contains no periodic sub-orbit, the cyle is prime and (Subscript[n, p])=n; otherwise the cycle is composite and n=(Subscript[n, p])r with integer r≥2.

When the map acts on d-dimensional space this is generalised to

((trℒ)^n)=(∑_(n-cycles p))((Subscript[n, p])/(|det(1-(Subscript[M, p]))|))=(∑_(n-cycles p))((Subscript[n, p])/(|1-(Subscript[Λ, p,1])|⋯|1-(Subscript[Λ, p,d])|))

where (Subscript[(Subscript[M, p])=((∂(f^n))/(∂x))|, p])=(Subscript[(Subscript[((∂f)/(∂x))|, (Subscript[x, p]^((n)))])…((∂f)/(∂x))|, (Subscript[x, p]^((1)))]) is the monodromy matrix, and (Subscript[Λ, p,k]) are its eigenvalues.

The trace formula is a generating function for ()

(∑_(n=1)^∞)(z^n)((trℒ)^n)=tr ((zℒ)/(1-zℒ))=(∑_n)(∑_(n-cycles p))(((z^n)(Subscript[n, p]))/(|det(1-(Subscript[M, p]))|))=(∑_(prime cycles p))(∑_r)(((Subscript[n, p]) (z^(r (Subscript[n, p]))))/(|det(1-(Subscript[M, p]^r))|))

Given that our expression for ((trℒ)^n) involves sums over all n-cycles, and that the number of such cycles grows exponentially with n, the size of term in the sum must be exponentially small for the sum to converge. This will happen if the map is uniformly hyperbolic, which means that 
1. There is at least one expanding monodromy eigenvalue (Subscript[Λ, p,k])>1 for some k,  
2. The are no neutral eigenvalues \[LeftBracketingBar](Subscript[Λ, p,k])\[RightBracketingBar]≠1 for all k, and 
3. The expanding eigenvalues (Subscript[Λ, p,e]) grow exponentially and contracting eigenvalues  (Subscript[Λ, p,c]) (if any) diminish exponentially, that is, there is a positive λ such that \[LeftBracketingBar](Subscript[Λ, p,e])\[RightBracketingBar]>(e^(λ (Subscript[n, p]))), and \[LeftBracketingBar](Subscript[Λ, p,c])\[RightBracketingBar]<(e^(-λ (Subscript[n, p]))). 
This guarantees that individual terms remain small. Neutrally stable eigenvalues \[LeftBracketingBar](Subscript[Λ, p,k])\[RightBracketingBar]=1 are not good because even if (Subscript[Λ, p,k])≠1, because arbitrary positive integral powers of (Subscript[Λ, p,k]) will appear for repeats of p, and these will eventually come close to 1.

The spectral determinant for maps

Equation () expresses the spectrum of the evolution operator in terms of periodic orbits, but is not convenient to use as is, becuase the eigenvalues appear in it as poles. Instead we will use it to derive a cycle expansion for the spectral determinant, det (1-zℒ), which has zeros at inverse eigenvalues. Starting from the identity log(det (1-zℒ))=tr(log(1-zℒ)) we get

Δ(z)=det(1-zℒ)=exp(-(∑_(n=1)^∞)(((z^n)((trℒ)^n))/n))=exp(- (∑_(prime cycles p))(∑_r)((z^(r (Subscript[n, p])))/(r|det(1-(Subscript[M, p]^r))|)))

where r as before is the number of repetitions of the cycle (=1 for prime cycles). The Taylor expansion of Δ

Δ(z)=1-z tr ℒ-((z^2)/2)(tr (ℒ^2)-(tr^2)ℒ)- ⋯

is most easily calculated from

((-z(d/(dz))det (1-zℒ))/(det (1-zℒ)))=-z(d/(dz))log(det (1-zℒ))=tr ((zℒ)/(1-zℒ))

so that if (Subscript[Δ, n]) is the nth Taylor of Δ(z) then

(1+(Subscript[Δ, 1])z+(Subscript[Δ, 2])(z^2)+⋯)(z tr ℒ+(z^2) (tr^2)ℒ+⋯)=-(Subscript[Δ, 1])z-2(Subscript[Δ, 2])(z^2)-3(Subscript[Δ, 3])(z^3)-⋯

and the resulting hierarchy of equations is solved recursively by

(Subscript[Δ, n])=-(1/n)(((trℒ)^n)+(Subscript[Δ, 1])((trℒ)^(n=1))+⋯+(Subscript[Δ, n-1])trℒ)

If (Subscript[Δ, n]) decay fast enough for large n, when can use polynomial truncations of the spectral determinant to obtain approximations for γ, whose accuracy increases as a function of n. If (Δ^((n)))(z) is the nth order Taylor expansion of Δ then the nth order approximation for the escape rate is (Subscript[γ, n])=log((Subscript[z, 0]^((n)))) where (Subscript[z, 0]^((n))) is the smallest positive real zero of (Δ^((n))). The results of numerical calculations of ((trℒ)^n), (Subscript[Δ, n]), and the nth order approximation to escape rate (Subscript[γ, n]) for the logistic map with r=4.5 are

TableForm@Prepend[Transpose@{Range[7],trace[lm[4.5,#]&,ilm[4.5,##]&,#]&/@Range[7],spectralDet[lm[4.5,#]&,ilm[4.5,##]&,7],NumberForm[#,16]&/@escapeRateTable[lm[4.5,#]&,ilm[4.5,##]&,7]},{"n","trace","det","escape rate"}]

The escape rate agrees beautifully with the decay rate of a histogram of survival based on one hundred thousand orbits starting from random initial points.

Module[{r=4.5,γ},γ=escapeRate[lm[4.5,#]&,ilm[4.5,##]&,6];SeedRandom[20150505];Show[Histogram[lifeSpan[lm[4.5,#]&,#,40]&/@RandomReal[{0.,1.},100000],{1},{"Log","Count"}],LogPlot[45000(E^(-γ n)),{n,1,30}]]]

The dynamical zeta function

The magic of the cycle expansion for the escape rate lies in the exquisite precision obtained using only a small number of orbits. To try to understand this phenomenon let us consider a slightly modifed cycle expansion.  We start by separating the expanding eigenvalues (Subscript[Λ, e,k]) and the contracting eigenvalues (Subscript[Λ, c,k]) appearing in the denominator in equation ()

|1-(Subscript[Λ, p,1])|⋯|1-(Subscript[Λ, p,d])|=|(Subscript[Λ, p])|(∏_k)|1-(1/(Subscript[Λ, e,k]))|(∏_k)|1-(Subscript[Λ, c,k])|

where (Subscript[Λ, p])=(∏_k)(Subscript[Λ, e,k]). For uniformly hyperbolic systems |1-(1/(Subscript[Λ, e,k]))| and  |1-(Subscript[Λ, c,k])| are very close to 1, so it makes sense to replace ((trℒ)^n) by

(Subscript[Γ, n])=(∑_(n-cycles p))((Subscript[n, p])/(|(Subscript[Λ, p])|))

in () obtaining

(1/(ζ(z)))=exp(-(∑_(n=1)^∞)(((z^n)(Subscript[Γ, n]))/n))=exp(- (∑_(prime cycles p))(∑_r)((z^(r (Subscript[n, p])))/(r((|(Subscript[Λ, p])|)^r))))=exp(- (∑_p)log(1-((z^(Subscript[n, p]))/(|(Subscript[Λ, p])|))))=(∏_p)(1-((z^(Subscript[n, p]))/(|(Subscript[Λ, p])|)))

It is possible to express ζ as a ratio of two spectral determinants, showing that (1/ζ) has the same zeros as the spectral determinant for nice enough systems.

The dynamical zeta function ζ(z) is a generalisation of the Riemann zeta function defined for Re s>1 by

(Subscript[ζ, R])(s)=(∑_(n=1)^∞)(1/(n^(-s)))=(∑_((Subscript[r, 1]),(Subscript[r, 2]),…))(1/((Subscript[p, 1]^(-s (Subscript[r, 1])))(Subscript[p, 2]^(-s (Subscript[r, 2])))⋯))=(∏_(primes p))(1/(1-(p^(-s))))

where (Subscript[p, 1]),(Subscript[p, 2]),… are the prime integers, and (Subscript[r, 1]),(Subscript[r, 2]),… are summed from zero to infinity.

We will calculate ζ explicitly for the piecewise linear repelling tent map

t(x)=({{{, ({{(Subscript[s, 0])x, x≤((Subscript[s, 1])/((Subscript[s, 0])+(Subscript[s, 1])))}, {(Subscript[s, 1])(1-x), x≥((Subscript[s, 1])/((Subscript[s, 0])+(Subscript[s, 1])))}})}})

with positive slopes (Subscript[s, 0]),(Subscript[s, 1]) such that (Subscript[s, 0])(Subscript[s, 1])>(Subscript[s, 0])+(Subscript[s, 1]). To find the periodic orbits we calculate the inverse images of the unit interval: (t^(-1))([0,1])=[0,(1/(Subscript[s, 0]))]⋃[1-(1/(Subscript[s, 1])),1], (t^(-2))([0,1])=[0,(1/(Subscript[s, 0]^2))]⋃⋯ etc. (t^(-n))([0,1]) is the union of (2^n) disjoint interval, each of which can be labelled by a sequence of branches (Subscript[b, n])(Subscript[b, n-1])⋯ (Subscript[b, 1]) of (t^(-1)) (read from right to left).

survPlot[tent[{2,-3},#]&,itent[{2,-3},##]&,5,3]

Each interval contains precisely one point of an n-cycle, namely (lim_(r->∞))((((Subscript[t, (Subscript[b, n])]^(-1))◦(Subscript[t, (Subscript[b, n-1])]^(-1))(Subscript[◦⋯◦t, (Subscript[b, 1])]^(-1))))^r)[0,1]. The monodromy of a (Subscript[b, n])(Subscript[b, n-1])⋯ (Subscript[b, 1]) cycle is ±(∏_(k=1)^n)(Subscript[s, (Subscript[b, k])]) and each cycle (prime or composite) is represented by (n/r) points so that

(Subscript[Γ, n])=(∑_((Subscript[b, 1])(Subscript[b, 2])(Subscript[…b, n])))(1/( (∏_(k=1)^n)(Subscript[s, (Subscript[b, k])])))=((((1/(Subscript[s, 0]))+(1/(Subscript[s, 1]))))^n)

and

(1/ζ)=exp(-(∑_n)((z^n)/n)((((1/(Subscript[s, 0]))+(1/(Subscript[s, 1]))))^n))=(e^(log(1-(z/(Subscript[s, 0]))-(z/(Subscript[s, 1])))))=1-((1/(Subscript[s, 0]))+(1/(Subscript[s, 1])))z

so that γ=log((1/(Subscript[s, 0]))+(1/(Subscript[s, 1]))).

How is it that the infinite product () is equal to a linear polynomial? Since each prime cycle in the product is weighted by (z^(Subscript[n, p])) we can find the low order Taylor coefficients of (1/ζ) by looking at the short prime cycles. The prime cycles of order up to 3 are 0, 1, 01, 001, 011, so that

(1/ζ)=(1-(z/(Subscript[s, 0])))(1-(z/(Subscript[s, 1])))(1-((z^2)/((Subscript[s, 0])(Subscript[s, 1]))))(1-((z^3)/((Subscript[s, 0]^2)(Subscript[s, 1]))))(1-((z^3)/((Subscript[s, 0])(Subscript[s, 1]^2))))⋯=1-((1/(Subscript[s, 0]))+(1/(Subscript[s, 1])))z-((1/((Subscript[s, 0])(Subscript[s, 1])))-(1/(Subscript[s, 0]))(1/(Subscript[s, 1])))(z^2)-((1/((Subscript[s, 0]^2)(Subscript[s, 1])))+(1/((Subscript[s, 0])(Subscript[s, 1]^2)))-((1/(Subscript[s, 0]))+(1/(Subscript[s, 1])))(1/((Subscript[s, 0])(Subscript[s, 1]))))(z^3)+⋯=1-((1/(Subscript[s, 0]))+(1/(Subscript[s, 1])))z

The exact cancellation of high order Taylor coefficients of (1/ζ) is a special property of piecewise linear maps, but the qualitative structure we found for the tent map is valid for any unimodal map repelling map f.

survPlot[lm[4.5,#]&,ilm[4.5,##]&,5,3]

As with the tent map the nth inverse image of the interval (f^(-n))([0,1]) consists of (2^n) disjoint intervals labelled by binary sequences of length n, each containing an n cycle, so that

(1/ζ)=(1-(z/(Subscript[Λ, 0])))(1-(z/(Subscript[Λ, 1])))(1-((z^2)/(Subscript[Λ, 01])))(1-((z^3)/(Subscript[Λ, 001])))(1-((z^3)/(Subscript[Λ, 011])))⋯=1-((1/(Subscript[Λ, 0]))+(1/(Subscript[Λ, 1])))z-((1/(Subscript[Λ, 01]))-(1/(Subscript[Λ, 0]))(1/(Subscript[Λ, 1])))(z^2)-((1/(Subscript[Λ, 001]))+(1/(Subscript[Λ, 011]))-((1/(Subscript[Λ, 0]))+(1/(Subscript[Λ, 1])))(1/(Subscript[Λ, 01])))(z^3)+⋯

The coefficient of (z^n) no longer cancel exactly in this case, but if f is smooth (1/(Subscript[Λ, (Subscript[p, 1])(Subscript[p, 2])]))-(1/(Subscript[Λ, (Subscript[p, 1])]))(1/(Subscript[Λ, (Subscript[p, 2])])), where (Subscript[p, 1])(Subscript[p, 2]) is the prime cycle obtained by juxtaposing the binary representations of (Subscript[p, 1]) and (Subscript[p, 2]), decays exponentially in the length of (Subscript[p, 1])(Subscript[p, 2]), and the Taylor expansion of (1/ζ) has a radius of convergence larger than 1 (Figures and table shown for the logistic map (Subscript[l, 4.5])).

rectPlot[lm[4.5,#]&,ilm[4.5,##]&,5,3]

TableForm@Prepend[Transpose@{Range[10],pcinvζ[lm[4.5,#]&,ilm[4.5,##]&,10],spectralDet[lm[4.5,#]&,ilm[4.5,##]&,10],NumberForm[#,16]&/@escapeRateζTable[lm[4.5,#]&,ilm[4.5,##]&,10],NumberForm[#,16]&/@escapeRateTable[lm[4.5,#]&,ilm[4.5,##]&,10]},{"n","(1/ζ)","det","escape rate (ζ)","escape rate (det)"}]

Dynamical averages

Symbolic dynamics and pruning

In this section our purpose is to calculate expectation values on the natural invariant measure of chaotic unimodal maps with no escaping orbits. The spectral determinant and  1/ζ necessarily have a zero at z=1 for such systems. However, before proceeding to study more interesting quantities we have to face a new issue.

For this purpose note that each infinite orbit of the repelling maps studied above was uniquely determined by its binary symbol sequence, and that each binary sequence corresponds to an admissible orbit. This is also true for the Ulam map (Subscript[l, 4]),

Module[{r=4.},prunePlot[lm[r,#]&,ilm[r,##]&,lmx[r],4,3]]

except that here the inverse images cover the entire interval, since the invariant set, an attractor, covers the entire interval. This is the reason that the dynamics is conjugate to the Bernoulli shift. In particular, every repeating binary sequence corresponds to an actual periodic orbit, and the cancellation of terms in the expansion of the dynamical ζ function is good.

For r<4 the situation is no longer as simple. Since max (Subscript[l, r])=(r/4)<1 the points in the subinterval ](r/4),1] are transient, and so are all their inverse images. Consequently, not all symbolic itineraries are allowed

Module[{r=3.831},prunePlot[lm[r,#]&,ilm[r,##]&,lmx[r],4,3]]

Here, for example, r=3.831, and the interval labelled 100 is pruned. The 100 periodic orbit is therefore also pruned, and there is only one prime 3-cycle. The pruning spoils the cancellation of the cycle expansions, which converge much more slowly, and nonuniformly:

Module[{n=10,r1=4.0,r2=3.9},TableForm@Prepend[Transpose@{Range[n],1+Accumulate@pcinvζ[lm[r1,#]&,ilm[r1,##]&,lmx[r1],n],1+Accumulate@pcinvζ[lm[r2,#]&,ilm[r2,##]&,lmx[r2],n]},{"n","\!\(\*FractionBox["1", Cell[TextData[{
Cell[BoxData[
FormBox[
SubscriptBox["ζ", 
RowBox[{"r", "=", "4"}]], TraditionalForm]],
FormatType->"TraditionalForm"],
"(1)"
}]]]\)","\!\(\*FractionBox["1", Cell[TextData[{
Cell[BoxData[
FormBox[
SubscriptBox["ζ", 
RowBox[{"r", "=", "3.9"}]], TraditionalForm]],
FormatType->"TraditionalForm"],
"(1)"
}]]]\)"}]]

The convergence improves significantly when the pruning is described by a finite set of pruning rules. Then the dynamics is conjugate to a finite Bernoulli subshift, and cancellations becomes effective once periods longer than the longest pruning rule are taken into account. An example of this type is the logistic map with r=(Subscript[r, gm])≈3.68, where (1/2) is mapped in four step to the unstable fixed point (Subscript[x, 1])

plotOrbit[lm[rgm,#]&,0.5,3]

The symbolic dynamics is such that every itinierary without two consecutive 0 symbols is allowed. To get good convergence the (OverBar[0]) cycle that is not in the attractor has to be pruned as well.

Module[{n=16,r1=4.0,r2=3.9,r3=rgm},TableForm@Prepend[Transpose@{Range[n],1+Accumulate@pcinvζ[lm[r1,#]&,ilm[r1,##]&,lmx[r1],n],1+Accumulate@pcinvζ[lm[r2,#]&,ilm[r2,##]&,lmx[r2],n],1+Accumulate@pcinvζdz[lm[r3,#]&,ilm[r3,##]&,lmx[r3],n]},{"n","\!\(\*FractionBox["1", Cell[TextData[{
Cell[BoxData[
FormBox[
SubscriptBox["ζ", 
RowBox[{"r", "=", "4"}]], TraditionalForm]],
FormatType->"TraditionalForm"],
"(1)"
}]]]\)","\!\(\*FractionBox["1", Cell[TextData[{
Cell[BoxData[
FormBox[
SubscriptBox["ζ", 
RowBox[{"r", "=", "3.9"}]], TraditionalForm]],
FormatType->"TraditionalForm"],
"(1)"
}]]]\) (with (OverBar[0]))","\!\(\*FractionBox["1", Cell[TextData[{
Cell[BoxData[
FormBox[
SubscriptBox["ζ", 
SubscriptBox["r", "gm"]], TraditionalForm]],
FormatType->"TraditionalForm"],
"(1)"
}]]]\) (w/out (OverBar[0]))"}]]

Operator eigenvalues as a generating function

Our next goal is to use cycle expansions to calculate expectation values of the form ⟨a⟩=∫a dμ, where a is a function of the dynamical variable x and μ is the natural invariant measure of a closed discrete dynamical system. For a chaotic system, this expectation value is equal to the long-time limit of a dynamical average taken on a typical orbit ⟨a⟩=(lim_(n->∞))(1/n)(∑_(k=1)^n)(Subscript[a, n]), with (Subscript[a, n])=a((f^n)((Subscript[x, 0]))).
The idea of cycle expansions for dynamical averages is to construct a weighted evolution operator (Subscript[ℒ, β]) whose leading eigenvalue is a moment generating function. Namely, let (Subscript[ℒ, β])ρ(x)=∫(e^(β a(f(ξ))))ρ(ξ)δ(x-f(ξ))dξ, so that

(Subscript[ℒ, β]^2)ρ(x)=∫(e^(β (a(f((Subscript[ξ, 1])))+a(f((Subscript[ξ, 2]))))))δ(x-f((Subscript[ξ, 1])))ρ((Subscript[ξ, 2]))δ((Subscript[ξ, 1])-f((Subscript[ξ, 2])))d(Subscript[ξ, 1])d(Subscript[ξ, 2])=∫(e^(β (a((f^2)(ξ))+a(f(ξ)))))δ(x-(f^2)(ξ))ρ(ξ)dξ

and

∫d(Subscript[xℒ, β]^n)ρ(x)=∫(e^(β (a(f(ξ))+a((f^2)(ξ))+⋯a((f^n)(ξ)))))ρ(ξ)dξ=⟨(e^(β (Subscript[A, n])))⟩

defining (Subscript[A, n])(x)=(∑_(k=1)^n)(Subscript[a, n]) so that ⟨a⟩=(lim_(n->∞))((⟨(Subscript[A, n])⟩)/n)=(lim_(n->∞))(1/n)(d/(dβ))log(Subscript[⟨(e^(β (Subscript[A, n])))⟩, β=0]).

Now, if (e^(s(β))) is the largest eigenvalue of (Subscript[ℒ, β]) then ⟨(e^(β (Subscript[A, n])))⟩(Underscript[∼, n->∞])(e^(n s(β))) and (Subscript[⟨a⟩=((ds)/(dβ))|, β=0]). We already know how to use cycle expansion to calculate (e^s): It is the inverse of smallest positive zero of (1/(Subscript[ζ, β])), the inverse zeta function corresponding to (Subscript[ℒ, β]). Retracing the steps leading to () we obtain that (Subscript[Γ, β,n])=(∑_(n-periods p))(((Subscript[n, p])(e^(β (Subscript[A, p]))))/(|(Subscript[Λ, p])|)), where (Subscript[A, p])=(∑_k)a(f((Subscript[x, p]^((k))))) and therefore

(1/((Subscript[ζ, β])(z)))=exp(-(∑_(n=1)^∞)(((z^n)(Subscript[Γ, β,n]))/n))=(∏_(prime cycles p))(1-(((z^(Subscript[n, p]))(e^(β (Subscript[A, p]))))/(|(Subscript[Λ, p])|)))

To get ⟨a⟩ we can differentiate the implicit relation (1/((Subscript[ζ, β])((e^(-s(β))))))=0 obtaining

(Subscript[((∂/(∂β))-s'(β)z(∂/(∂z)))(1/(Subscript[ζ, β]))|, β=0,z=1])=0 ⟹ ⟨a⟩=s'(0)=(((∑_p)(∏_(q≠p))(1 - (1/(|(Subscript[Λ, q])|)))((Subscript[A, p])/(|(Subscript[Λ, p])|)))/((∑_p)(∏_(q≠p))(1 - (1/(|(Subscript[Λ, q])|)))((Subscript[n, p])/(|(Subscript[Λ, p])|))))=(([A])/([n]))

defining []  for additive functions (Subscript[B, p])=(∑_n) b((Subscript[x, p]^((n)))), n being additive for the constant function 1.

Examples: Right dwell fraction and the Lyapunov exponent

Two example dynamical averages will be shown for the logistic map. The right dwell fraction is the probability  to visit the decreasing part of the map, 𝕡(x>(r/8))=⟨Θ(x-(r/8))⟩, where Θ is the unit step function

Module[{n=20,r1=4.0,r2=3.9,r3=rgm,x0=0.6,nn=1000000},TableForm@Join[{{"𝕡(x>(r/8))","r=4","r=3.9","(Subscript[r, gm]) (w/out (OverBar[0]))"}},Transpose@{Range[n],dynAvg[N[UnitStep[#-(r1/8)]]&,lm[r1,#]&,ilm[r1,##]&,lmx[r1],n],dynAvg[N[UnitStep[#-(r2/8)]]&,lm[r2,#]&,ilm[r2,##]&,lmx[r2],n],dynAvgDZ[N[UnitStep[#-(r3/8)]]&,lm[r3,#]&,ilm[r3,##]&,lmx[r3],n]},{{"(10^6) orbit",orbitAv[N[UnitStep[#-(r1/8)]]&,lm[r1,#]&,x0,nn],orbitAv[N[UnitStep[#-(r2/8)]]&,lm[r2,#]&,x0,nn],orbitAv[N[UnitStep[#-(r3/8)]]&,lm[r3,#]&,x0,nn]},{"2×(10^6) orbit",orbitAv[N[UnitStep[#-(r1/8)]]&,lm[r1,#]&,x0,2nn],orbitAv[N[UnitStep[#-(r2/8)]]&,lm[r2,#]&,x0,2nn],orbitAv[N[UnitStep[#-(r3/8)]]&,lm[r3,#]&,x0,2nn]},{"Exact",0.5,"",""}}]]

The second example is the Lyapunov exponent λ=⟨log|f'(x)|⟩

Module[{n=18,r1=4.0,r2=3.9,r3=rgm,x0=0.6,nn=100000},TableForm@Join[{{"lyapunov","r=4","r=3.9","(Subscript[r, gm]) (w/out (OverBar[0]))"}},Transpose@{Range[n],dynAvg[Log[Abs[Derivative[0,1][lm][r1,#]]]&,lm[r1,#]&,ilm[r1,##]&,lmx[r1],n],dynAvg[Log[Abs[Derivative[0,1][lm][r2,#]]]&,lm[r2,#]&,ilm[r2,##]&,lmx[r2],n],dynAvgDZ[Log[Abs[Derivative[0,1][lm][r3,#]]]&,lm[r3,#]&,ilm[r3,##]&,lmx[r3],n]},{{"(10^6) orbit",lyapunov[lm[r1,#]&,x0,nn],lyapunov[lm[r2,#]&,x0,nn],lyapunov[lm[r3,#]&,x0,nn]},{"2×(10^6) orbit",lyapunov[lm[r1,#]&,x0,2nn],lyapunov[lm[r2,#]&,x0,2nn],lyapunov[lm[r3,#]&,x0,2nn]},{"Exact","Log[2] = "<>ToString[Log[2.]],"",""}}]]

The cycle expansions seem to work, but the performance is mediocre. Although the convergence is exponential - the number of significant digits grows linearly with the maximum cycle length, the performance is comparable to that of simple orbit averaging, because the number of periodic orbits grows exponentially with their length. In practice, it is therefore advantageous to accelerate the convergence by replacing the ζ function woth a spectral determinant, control of the pruning process to allow for better cancellations, tail resummation, or other methods. These are beyond the present scope.

The Lorenz attractor

Qualitative dynamics

The Lorenz dynamics has three fixed points, one at the origin, and two, related by symmetry, at the centres of the wings. For the standard parameter values all fixed points are unstable, that is, at least one of the eigenvalues of the Jacobian matrix ∂f/∂x has positive real part. The Jacobian eigenvectors associated with the unstable eigenvalues span unstable directions. Because the Lorenz system is strongly dissipative there are also stable directions at each fixed point, so they are all hyperbolic.

The wing fixed points have a complex conjugate pair of unstable eigenvalues which generate the outward spiralling, and a strongly stable eigenvalue which compresses the dynamics to the almost two-dimensional wings. The fixed point at the origin has one unstable direction, tangent to the intersection of the attractor with the xy plane, a stable direction along the z axis, and a strongly stable direction transverse to the attractor

Module[{σ=10,b=(8/3),r=28,x0=5,y0=-5,z0=3,td=100,lx,ly,lz,gr,sq=25},{lx,ly,lz}=NDSolveValue[{x'[t]==σ (-x[t]+y[t]),y'[t]==r x[t]-y[t]-x[t] z[t],z'[t]==x[t] y[t]-b z[t], x[0]==0.0001,y[0]==0,z[0]==r-1},{x,y,z},{t,0,td}];Show[ParametricPlot3D[{lx[t],ly[t],lz[t]},{t,0,td},PlotRange->{{-30,30},{-30,30},{0,50}},PlotPoints->200,PlotStyle->Thin],Graphics3D[{ColorData[97,2],PointSize[Large],Point[{0,0,0}],ColorData[97,3],Point[{(Sqrt[b(r-1)]),(Sqrt[b(r-1)]),r-1}],Point[{-(Sqrt[b(r-1)]),-(Sqrt[b(r-1)]),r-1}]}]]]

This information is sufficient to understand the qualitative dynamics. Since orbits leave the origin in the unstable direction, the unstable manifold of the origin fixed point delimits the outer boundaries of the wings. These orbits are outside the wing switching radius, so they continue to the other wing, where they reach the as close as possible to the wing fixed points, and start to spiral out close to the unstable manifold of the wing fixed points. The intersection of the wing with the two dimensional stable manifold of the origin fixed point is the separatrix that separates between staying and swicthing orbits

The Poincaré map

A Poincaré section of a d-dimensional continuous dynamical system is simply a set of intersection points of an orbit or orbits with a d-1-dimensional submanifold. The key observation is that the flow induces a map on the section, the Poincaré or return map, which maps each intersection point to the next one. A good section for the Lorenz system is the plane z=r-1 which passes through the wing fixed points.

Module[{σ=10,b=(8/3),r=28,x0=5,y0=-5,z0=3,td=100,lx,ly,lz,gr,sq=25},{lx,ly,lz}=NDSolveValue[{x'[t]==σ (-x[t]+y[t]),y'[t]==r x[t]-y[t]-x[t] z[t],z'[t]==x[t] y[t]-b z[t], x[0]==0.0001,y[0]==0,z[0]==r-1},{x,y,z},{t,0,td}];Show[ParametricPlot3D[{lx[t],ly[t],lz[t]},{t,0,td},PlotRange->{{-30,30},{-30,30},{0,50}},PlotPoints->200,PlotStyle->Thin],Graphics3D[Polygon[{{-sq,-sq,r-1},{-sq,sq,r-1},{sq,sq,r-1},{sq,-sq,r-1}}]]]]

The section of the orbit shown above, which starts on the section plane close to the z axis is

Module[{σ=10,b=(8/3),r=28,fpPos},fpPos={(Sqrt[b(r-1)]),(Sqrt[b(r-1)]),r-1};Show[ListPlot[{Select[losection,Abs[#[[1]]]>(Sqrt[b(r-1)])&],Select[losection,Abs[#[[1]]]<(Sqrt[b(r-1)])&]},PlotStyle->PointSize[Tiny]],Graphics[MapIndexed[{PointSize[Large],GrayLevel[0.1First[#2],0.6],Point[#1]}&,Take[losection,8]]],Graphics[{{ColorData[97,4],PointSize[Large],Point[{0,0}],ColorData[97,4],Point[-Take[fpPos,2]-{lomax,iUsmPlus[lomax]}],Point[Take[fpPos,2]+{lomax,iUsmPlus[lomax]}]},{PointSize[Large],ColorData[97,3],Point[{(Sqrt[b(r-1)]),(Sqrt[b(r-1)])}]},{PointSize[Large],ColorData[97,3],Point[{-(Sqrt[b(r-1)]),-(Sqrt[b(r-1)])}]}}],AspectRatio->1]]

Points that intersect the plane going up are marked in blue and going down in orange. The two green points are the wing fixed points, and the red points are on the stable manifold of the fixed point at the origin. The point in the center is the initial point, which is not on the attractor. The other starting points of the orbit are marked in grey: the orbit first visits the right wing, going up and then down at its outer boundary, and then hits the inner part of the left wing, where it starts spiralling outward.

It is convenient to define the return map between upward going points. Then a mapped point remains on the same wing if it is inside the origin stable manifold intersection (outer red points on the graph, with approximate displacement 4.76, 9.92 from the wing fixed point). In this inner region the distance of the mapped point from the fixed point is monotonically increasing as a function of the distance of the original point. Beyond the separatrix points are mapped to the other wing, and the distance of the mapped point from the fixed point (of the new wing) is monotonically decreasing as a function of the distance of the original point. This explains the structure of the return map, here parameterized by the x distance from the fixed point

ListPlot[Array[lomap,800,{0.01,8}],Joined->True,DataRange->{0.01,8}]

To be sure, the Poincare map arises from an invertible flow, so it is also invertible, and the two inverse images of 4, say, are not mapped to precisely the same point; rather their images are very slightly spaced in the transverse, fractal, direction of the wing. This distinction however hardly matters for the statistical properties we’re studying here.

Switching rate & the Lyapunov exponent

The Lorenz system dynamics having been effectively reduced to a unimodal map, can now be analyzed using the ideas and tools developed above. I will present two examples, the first of which is the switching rate (Subscript[γ, d]) between the two wings. Given that points on the far side of the attractor, with x>(Subscript[x, m])≈4.76 switch and points with x<(Subscript[x, m]) stay, (Subscript[γ, d])=𝕡(x>(Subscript[x, m])). This gives the discrete version of the switching rate, where ‘time’ is measured with attractor cycles. The continuous switching rate (Subscript[γ, c]) can be calculated by the flow-version of the dynamical ζ function, which was not presented here, is equal to (([Θ(x-(Subscript[x, m]))])/([T])) instead of (([Θ(x-(Subscript[x, m]))])/([n])) for (Subscript[γ, d]) where (Subscript[T, p]) is the flow time of cycle p.

Module[{n=16,x0=5.,nn=100000},TableForm@Join[{{"n","(Subscript[γ, d])","(Subscript[γ, c])"}},Transpose@{Range[n],dynAvg[N[UnitStep[#-4.76]]&,ilomap,iilomap[#1][#2]&,lomin,lomax,n],((dynAvg[N[UnitStep[#-4.76]]&,ilomap,iilomap[#1][#2]&,lomin,lomax,n])/(dynAvg[ilomapt,ilomap,iilomap[#1][#2]&,lomin,lomax,n]))},{{"(10^5)/t=(10^4) orbit",orbitAv[N[UnitStep[#-4.76]]&,ilomap,x0,nn],switchRate[(nn/10.),x0]},{"2×(10^5)/t=2×(10^4) orbit",orbitAv[N[UnitStep[#-4.76]]&,ilomap,x0,2nn],switchRate[(nn/5.),x0]}}]]

The second example is the Lyapunov exponent, again with discrete version (Subscript[λ, d]) and continuous version (Subscript[λ, c]).

Module[{n=16,x0=5.,nn=100000},TableForm@Join[{{"n","(Subscript[λ, d])","(Subscript[λ, c])"}},Transpose@{Range[n],dynAvg[Log[Abs[ilomap'[#]]]&,ilomap,iilomap[#1][#2]&,lomin,lomax,n],((dynAvg[Log[Abs[ilomap'[#]]]&,ilomap,iilomap[#1][#2]&,lomin,lomax,n])/(dynAvg[ilomapt,ilomap,iilomap[#1][#2]&,lomin,lomax,n]))},{{"(10^5)/t=(10^3) orbit",lyapunov[ilomap,x0,nn],lolyapunov[(nn/100.),x0]},{"2×(10^5)/t=2×(10^3) orbit",lyapunov[ilomap,x0,2nn],lolyapunov[(nn/50.),x0]}}]]

Auxiliary definitions

lm[r_,x_]:=r x (1-x)

ilm[r_,s_,y_]:=(1/2)(1+(2s-1)(Sqrt[1-4(y/r)]))

lmx[r_]:=(r/4)

rgm=r/.FindRoot[Last@(Differences@NestList[lm[r,#]&,0.5,4]),{r,3.68,3.67}]

tent[{ls_,rs_},x_]:=Piecewise[{{ls x,x≤((-rs)/(ls-rs))},{rs(x-1),x>((-rs)/(ls-rs))}}]

itent[{ls_,rs_},0,y_]:=(y/ls)

itent[{ls_,rs_},1,y_]:=1+(y/rs)

plotOrbit[f_,s_,n_]:=Module[{orbit},orbit=NestWhileList[f[#]&,s,#<1&,1,n];Show[Plot[{f[x],x},{x,0,1},PlotRange->{0,1},PlotStyle->Thick],Graphics[Line[Join[{{orbit[[1]],0}},Flatten[Table[{{orbit[[k-1]],orbit[[k]]},{orbit[[k]],orbit[[k]]}},{k,2,Length[orbit]}],1]]]]]]

plotOrbits[f_,s_List,n_,st_List]:=Module[{orbits},orbits=NestWhileList[f[#]&,#,#<1&,1,n]&/@s;Show[Plot[{f[x],x},{x,0,1},PlotRange->{0,1},PlotStyle->Thick],Graphics[Transpose[{st,Line[Join[{{#[[1]],0}},Flatten[Table[{{#[[k-1]],#[[k]]},{#[[k]],#[[k]]}},{k,2,Length[#]}],1]]]&/@orbits}]]]]

lifeSpan[f_,s_,cutoff_]:=Module[{n=0},NestWhileList[(n++;f[#])&,s,#<1&,1,cutoff];n]

imapInterval[if_,0,int_List,fmx_]:=if[0,#]&/@int

imapInterval[if_,1,int_List,fmx_]:=Module[{ii},ii=if[1,#]&/@int;If[ii[[2]]>fmx,Throw["pruned"]];Reverse[MapAt[Min[#,fmx]&,ii,1]]]

imapInterval[if_,0,int_List,fmn_,fmx_]:=Module[{ii},ii=if[0,#]&/@int;If[ii[[2]]<fmn,Throw["pruned"]];MapAt[Max[#,fmn]&,ii,1]]

imapInterval[if_,1,int_List,fmn_,fmx_]:=Module[{ii},ii=if[1,#]&/@int;If[ii[[2]]>fmx,Throw["pruned"]];Reverse[MapAt[Min[#,fmx]&,ii,1]]]

unimodalPO[if_,bits_,iter_]:=Fold[if[#2,#1]&,0.5,Flatten@ConstantArray[bits,Ceiling@(iter/(Length[bits]))]]

unimodalPO[if_,bits_,fmx_,iter_]:=Catch[Mean@Fold[imapInterval[if,#2,#1,fmx]&,{0,fmx},Flatten[ConstantArray[bits,Ceiling[(iter/(Length[bits]))]]]]]

unimodalPO[if_,bits_,fmn_,fmx_,iter_]:=Catch[Mean@Fold[imapInterval[if,#2,#1,fmn,fmx]&,{fmn,fmx},Flatten[ConstantArray[bits,Ceiling[(iter/(Length[bits]))]]]]]

monodromy[f_,x_,n_]:=Module[{λ=1},NestWhile[(λ*=f'[#];f[#])&,x,#<1&,1,n];λ]

orbitAv[b_,f_,x_,n_]:=Module[{av=0},NestList[(av+=b[#];f[#])&,x,n];(av/n)]

lyapunov[f_,x_,n_]:=orbitAv[Log[Abs[f'[#]]]&,f,x,n]

lyapunov[f_,x_,n_]:=Module[{λ=0},NestList[(λ+=Log[Abs[f'[#]]];f[#])&,x,n];(λ/n)]

trace[f_,if_,n_]:=With[{iter=100},Total@(1/(Abs[monodromy[f,#,n]&/@unimodalPO[if,n,iter]-1]))]

spectralDet[f_,if_,n_]:=Module[{tr,det={}},tr=trace[f,if,#]&/@Range[n];Do[AppendTo[det,-((tr[[k]]+Reverse[det].Take[tr,k-1])/k)],{k,n}];det]

escapeRate[f_,if_,n_]:=Log[x/.FindRoot[1+spectralDet[f,if,n].(x^(Range[n])),{x,0}]]

escapeRateTable[f_,if_,n_]:=Module[{sd},sd=1+spectralDet[f,if,n].(x^(Range[n]));Log[x/.FindRoot[Take[sd,#+1],{x,0}]]&/@Range[n]]

primeCycles[n_]:= (* With cyclic permutations *)primeCycles[n]=Delete[Range[(2^(n-1)),(2^n)-1],List/@Flatten[((((2^n)-1)/((2^(#))-1)))primeCycles[#]&/@Drop[Divisors[n],-1]]-(2^(n-1))+1]

primeCycles[1]={1}

removeCyclic[1]={0,1}

removeCyclic[n_]:= Cases[DeleteDuplicatesBy[Flatten[Transpose@Table[{k,RotateRight[#,k]}&/@IntegerDigits[primeCycles[n],2],{k,0,n-1}],1],Last],{0,l_List}:>FromDigits[l,2]]

discardZero[k_]:=removeCyclic[k]

discardZero[1]={1}

pcinvζ[f_,if_,n_]:=Module[{iter=100,k,ζ,tfactors},ζ=Prepend[ConstantArray[0,n],1];Do[For[tfactors=-(1/(Abs[monodromy[f,#,k]]))&/@(Fold[if[#2,#1]&,0.5,Flatten@ConstantArray[#,Ceiling@(iter/k)]]&/@IntegerDigits[removeCyclic[k],2]),tfactors≠{},tfactors=Rest[tfactors],ζ+=Take[PadLeft[ζ First[tfactors],n+1+k],n+1]],{k,n}];Rest[ζ]]

pcinvζ[f_,if_,fmx_,n_]:=Module[{iter=100,k,ζ,tfactors},ζ=Prepend[ConstantArray[0,n],1];Do[For[tfactors=-(1/(Abs[monodromy[f,#,k]]))&/@Select[(unimodalPO[if,#,fmx,iter]&/@IntegerDigits[removeCyclic[k],2]),NumberQ],tfactors≠{},tfactors=Rest[tfactors],ζ+=Take[PadLeft[ζ First[tfactors],n+1+k],n+1]],{k,n}];Rest[ζ]]

pcinvζ[f_,if_,fmn_,fmx_,n_]:=Module[{iter=100,k,ζ,tfactors},ζ=Prepend[ConstantArray[0,n],1];Do[For[tfactors=-(1/(Abs[monodromy[f,#,k]]))&/@Select[(unimodalPO[if,#,fmn,fmx,iter]&/@IntegerDigits[removeCyclic[k],2]),NumberQ],tfactors≠{},tfactors=Rest[tfactors],ζ+=Take[PadLeft[ζ First[tfactors],n+1+k],n+1]],{k,n}];Rest[ζ]]

pcinvζdz[f_,if_,fmx_,n_]:=Module[{iter=100,k,ζ,tfactors},ζ=Prepend[ConstantArray[0,n],1];Do[For[tfactors=-(1/(Abs[monodromy[f,#,k]]))&/@Select[(unimodalPO[if,#,fmx,iter]&/@IntegerDigits[discardZero[k],2]),NumberQ],tfactors≠{},tfactors=Rest[tfactors],ζ+=Take[PadLeft[ζ First[tfactors],n+1+k],n+1]],{k,n}];Rest[ζ]]

escapeRateζ[f_,if_,n_]:=Log[x/.FindRoot[1+pcinvζ[f,if,n].(x^(Range[n])),{x,0}]]

escapeRateζTable[f_,if_,n_]:=Module[{sd},sd=1+pcinvζ[f,if,n].(x^(Range[n]));Log[x/.FindRoot[Take[sd,#+1],{x,0}]]&/@Range[n]]

dynAvg[b_,f_,if_,fmx_,n_]:=Module[{iter=100,k,ζ,da,na,xp,tfactors,bterms},ζ=Prepend[ConstantArray[0,n-1],1];na=da=ConstantArray[0,n];Do[For[xp=Select[(unimodalPO[if,#,fmx,iter]&/@IntegerDigits[removeCyclic[k],2]),NumberQ];tfactors=-(1/(Abs[monodromy[f,#,k]]))&/@xp;bterms=Total[b/@NestList[f,#,k-1]]&/@xp;,tfactors≠{},tfactors=Rest[tfactors];bterms=Rest[bterms],da+=Take[PadLeft[da First[tfactors],n+k],n]-Take[PadLeft[First[bterms]First[tfactors]ζ,n+k-1],n];na+=Take[PadLeft[na First[tfactors],n+k],n]-Take[PadLeft[k First[tfactors]ζ,n+k-1],n];ζ+=Take[PadLeft[ζ First[tfactors],n+k],n]],{k,n}];((Accumulate[da])/(Accumulate[na]))]

dynAvg[b_,f_,if_,fmn_,fmx_,n_]:=Module[{iter=100,k,ζ,da,na,xp,tfactors,bterms},ζ=Prepend[ConstantArray[0,n-1],1];na=da=ConstantArray[0,n];Do[For[xp=Select[(unimodalPO[if,#,fmn,fmx,iter]&/@IntegerDigits[removeCyclic[k],2]),NumberQ];tfactors=-(1/(Abs[monodromy[f,#,k]]))&/@xp;bterms=Total[b/@NestList[f,#,k-1]]&/@xp;,tfactors≠{},tfactors=Rest[tfactors];bterms=Rest[bterms],da+=Take[PadLeft[da First[tfactors],n+k],n]-Take[PadLeft[First[bterms]First[tfactors]ζ,n+k-1],n];na+=Take[PadLeft[na First[tfactors],n+k],n]-Take[PadLeft[k First[tfactors]ζ,n+k-1],n];ζ+=Take[PadLeft[ζ First[tfactors],n+k],n]],{k,n}];((Accumulate[da])/(Accumulate[na]))]

dynAvgDZ[b_,f_,if_,fmx_,n_]:=Module[{iter=100,k,ζ,da,na,xp,tfactors,bterms},ζ=Prepend[ConstantArray[0,n-1],1];na=da=ConstantArray[0,n];Do[For[xp=Select[(unimodalPO[if,#,fmx,iter]&/@IntegerDigits[discardZero[k],2]),NumberQ];tfactors=-(1/(Abs[monodromy[f,#,k]]))&/@xp;bterms=Total[b/@NestList[f,#,k-1]]&/@xp;,tfactors≠{},tfactors=Rest[tfactors];bterms=Rest[bterms],da+=Take[PadLeft[da First[tfactors],n+k],n]-Take[PadLeft[First[bterms]First[tfactors]ζ,n+k-1],n];na+=Take[PadLeft[na First[tfactors],n+k],n]-Take[PadLeft[k First[tfactors]ζ,n+k-1],n];ζ+=Take[PadLeft[ζ First[tfactors],n+k],n]],{k,n}];((Accumulate[da])/(Accumulate[na]))]

survIntervals[if_,n_]:=Partition[#,2]&@(Sort@(Fold[if[#2,#1]&,0.,#]&/@Tuples[{0,1},n]))

itinerary[bits_List]:=Module[{w=bits},Do[If[w[[i]]==1,w[[i+1;;]]=1-w[[i+1;;]]],{i,Length[w]}];w]

survPlot[f_,if_,ng_,ns_]:=Show[Plot[f[x],{x,0,1},PlotStyle->Black],NumberLinePlot[Table[Interval@@survIntervals[if,d],{d,ng}],Spacings->{-0.2}],Table[ListPlot[Partition[NestList[f,unimodalPO[if,#,100],k],2,1]&/@IntegerDigits[removeCyclic[k],2],Filling->-1, PlotStyle->ColorData[97][k+1]],{k,1,ns}],Graphics[Table[MapThread[Text,{StringJoin/@Map[ToString,itinerary/@Tuples[{0,1},k],{2}],Outer[Sequence,Mean/@survIntervals[if,k+1],{-.2k-.15}]}],{k,ns}]],PlotRange->{-0.2(ng+1),0.9}]

rectPlot[f_,if_,ng_,ns_]:=Show[Plot[f[x],{x,0,1},PlotStyle->Directive[Thin]],NumberLinePlot[Table[Interval@@survIntervals[if,d],{d,ng}],Spacings->{-0.2}],Table[ListPlot[Partition[NestList[f,unimodalPO[if,#,100],k],2,1]&/@IntegerDigits[removeCyclic[k],2], Filling->-1,PlotStyle->ColorData[97][k+1]],{k,1,ns}],Graphics[Table[{ColorData[97][k+1],Line/@(tangent[f,#,(0.5/(2^k))]&/@Partition[NestList[f,unimodalPO[if,#,100],k],2,1]&/@IntegerDigits[removeCyclic[k],2])},{k,1,ns}]],Graphics[Table[MapThread[Text,{StringJoin/@Map[ToString,itinerary/@Tuples[{0,1},k],{2}],Outer[Sequence,Mean/@survIntervals[if,k+1],{-.2k-.15}]}],{k,ns}]],PlotRange->{-0.2(ng+1),0.9}]

pruneIntervals[if_,fmx_,n_]:=Catch[{#,Fold[imapInterval[if,#2,#1,fmx]&,{0,fmx},#]}]&/@Tuples[{0,1},n-1]

prunePlot[f_,if_,fmx_,ng_,ns_]:=Module[{symbols,intervals,po},Do[{symbols[k],intervals[k]}=Transpose[SortBy[#[[2,1]]&]@Cases[pruneIntervals[if,fmx,k],_List]];po[k-1]=Select[unimodalPO[if,#,fmx,100]&/@IntegerDigits[removeCyclic[k-1],2],NumberQ],{k,2,ng}];Show[Plot[f[x],{x,0,1},PlotStyle->Black],NumberLinePlot[Table[Interval/@intervals[k],{k,2,ng}],Spacings->{-0.2}],Table[ListPlot[Partition[NestList[f,#,k],2,1]&/@po[k],Filling->-0.2 (ng-1), PlotStyle->ColorData[97][k+1]],{k,1,ns}],Graphics[Table[MapThread[Text,{StringJoin/@Map[ToString,Reverse/@symbols[k+1],{2}],Outer[Sequence,Mean/@intervals[k+1],{-.2k+.05}]}],{k,ns}]],PlotRange->{-0.2(ng),1.0}]]

lo[x_]:={σ (-x[1]+x[2]),r x[1]-x[2]-x[1]x[3],x[1]x[2]-b x[3]}

lofp=Solve[MapThread[#1==#2&,{lo[x],{0,0,0}}],Array[x,3]]

lojac=Grad[#,Array[x,3]]&/@lo[x]

lopar={σ->10,b->(8/3),r->28};

unstable=Drop[Eigenvectors[lojac/.lofp[[2]]]/.N[lopar],1]

us1=Chop[-I (Subtract@@unstable)]

usmPlus=Module[{σ=10,b=(8/3),r=28,fpPos,nudge=0.01,us1,td=65.,usm,refine,usm1,usm2,nr=0},fpPos={(Sqrt[b(r-1)]),(Sqrt[b(r-1)]),r-1};us1=Chop[-I (Subtract@@unstable)];usm=Select[Reap[NDSolve[{x'[t]==σ (-x[t]+y[t]),y'[t]==r x[t]-y[t]-x[t] z[t],z'[t]==x[t] y[t]-b z[t], x[0]==fpPos[[1]]+nudge us1[[1]],y[0]==fpPos[[2]]+nudge us1[[2]],z[0]==fpPos[[3]],WhenEvent[z[t]==r-1,Sow[{x[t],y[t]}-Take[fpPos,2]]],WhenEvent[x[t]<0,"StopIntegration"]},{},{t,0,td}]][[-1,1]],Positive[Total[#]]&];While[nr++<450&&(refine=First@FirstPosition[Total[(#[[2]]-#[[1]])]&/@Partition[usm,2,1],_?(#>0.01&),{"None"}])≠ "None",usm1=Select[Reap[NDSolve[{x'[t]==σ (-x[t]+y[t]),y'[t]==r x[t]-y[t]-x[t] z[t],z'[t]==x[t] y[t]-b z[t], x[0]==fpPos[[1]]+((usm[[refine-nr+1,1]]+usm[[refine-nr,1]])/2),y[0]==fpPos[[2]]+((usm[[refine-nr+1,2]]+usm[[refine-nr,2]])/2),z[0]==fpPos[[3]],WhenEvent[z[t]==r-1,Sow[{x[t],y[t]}-Take[fpPos,2]]],WhenEvent[x[t]<0,"StopIntegration"]},{},{t,0,td}]][[-1,1]],Positive[Total[#]]&];usm=Sort@Join[usm,usm1]];usm];

iUsmPlus=Interpolation[usmPlus]

lomap[xi_?NumericQ]:=Module[{σ=10,b=(8/3),r=28,fpPos,td=65.,lx,ly,lz,gr,usm,refine,usm1,usm2,nr=0,yi},yi=iUsmPlus[xi];fpPos={(Sqrt[b(r-1)]),(Sqrt[b(r-1)]),r-1};Reap[NDSolve[{x'[t]==σ (-x[t]+y[t]),y'[t]==r x[t]-y[t]-x[t] z[t],z'[t]==x[t] y[t]-b z[t], x[0]==fpPos[[1]]+xi,y[0]==fpPos[[2]]+yi,z[0]==fpPos[[3]],WhenEvent[z[t]==r-1&&(x[t]>fpPos[[1]]||x[t]<-fpPos[[1]]),Sow[If[x[t]>0,x[t]-fpPos[[1]],-x[t]-fpPos[[1]]]];"StopIntegration"]},{},{t,0,td}]][[-1,1,1]]]

lomapt[xi_?NumericQ]:=Module[{σ=10,b=(8/3),r=28,fpPos,td=65.,lx,ly,lz,gr,usm,refine,usm1,usm2,nr=0,yi},yi=iUsmPlus[xi];fpPos={(Sqrt[b(r-1)]),(Sqrt[b(r-1)]),r-1};Reap[NDSolve[{x'[t]==σ (-x[t]+y[t]),y'[t]==r x[t]-y[t]-x[t] z[t],z'[t]==x[t] y[t]-b z[t], x[0]==fpPos[[1]]+xi,y[0]==fpPos[[2]]+yi,z[0]==fpPos[[3]],WhenEvent[z[t]==r-1&&(x[t]>fpPos[[1]]||x[t]<-fpPos[[1]]),Sow[t];"StopIntegration"]},{},{t,0,td}]][[-1,1,1]]]

switchRate[td_,xi_]:=Module[{σ=10,b=(8/3),r=28,fpPos,lx,ly,lz,gr,nr=0,yi},yi=iUsmPlus[xi];fpPos={(Sqrt[b(r-1)]),(Sqrt[b(r-1)]),r-1};NDSolve[{x'[t]==σ (-x[t]+y[t]),y'[t]==r x[t]-y[t]-x[t] z[t],z'[t]==x[t] y[t]-b z[t], x[0]==fpPos[[1]]+xi,y[0]==fpPos[[2]]+yi,z[0]==fpPos[[3]],WhenEvent[x[t]==-y[t],nr++]},{},{t,0,td}];(nr/td)]

lolyapunov[td_,xi_]:=Module[{σ=10,b=(8/3),r=28,fpPos,lx,ly,lz,lv,lλ,gr,yi,lor},yi=iUsmPlus[xi];fpPos={(Sqrt[b(r-1)]),(Sqrt[b(r-1)]),r-1};{lx,ly,lz,lv,lλ}=NDSolveValue[{x'[t]==σ (-x[t]+y[t]),y'[t]==r x[t]-y[t]-x[t] z[t],z'[t]==x[t] y[t]-b z[t],λ'[t]==Log[Abs[v[t].{{-σ,σ,0},{r-z[t],-1,-x[t]},{y[t],x[t],-b}}.v[t]]],v'[t]==({{-σ,σ,0},{r-z[t],-1,-x[t]},{y[t],x[t],-b}}-v[t].{{-σ,σ,0},{r-z[t],-1,-x[t]},{y[t],x[t],-b}}.v[t]IdentityMatrix[3]).v[t], x[0]==fpPos[[1]]+xi,y[0]==fpPos[[2]]+yi,z[0]==fpPos[[3]],λ[0]==0,v[0]=={1,0,0}},{x,y,z,v,λ},{t,0,td}];((lλ[td])/td)]

ilomap=ListInterpolation[Array[lomap,801,{0.0,8}],{0.0,8},InterpolationOrder->3,Method->"Spline"]

ilomapt=ListInterpolation[Array[lomapt,801,{0.0,8}],{0.0,8},InterpolationOrder->3,Method->"Spline"]

iilomap[0]=Interpolation@Array[{lomap[#],#}&,477,{0.0,4.76}]

iilomap[1]=Interpolation@Array[{lomap[#],#}&,800-475,{4.76,8.}]

lomin=1.41;lomax=7.84;