# Optimize this computation of Liapunov exponent

Introduction

In dynamical systems there is an important exponent called Liapunov exponent. I'm trying to create a visual representation of Liapunov exponent for one dimension finite-difference equations. For example, the most common equation is the logistic map: $$(1)\quad x_{n + 1} = rx_{n}(1 - x_{n})$$, where $$x_{n}$$ is the $$n$$-th iteration of (1), we define $$x_0 = x$$, $$n$$ goes from 0 to $$\infty$$, and $$r$$ goes from 0 to 4. Equation (1) is discrete, but we can consider the continuous function $$(2)\quad f(x) = rx(1 - x)$$ in order to do all the computations (I will not prove nothing here), then $$(3)\quad f'(x) = r(1 - 2x)$$.

Finally, Liapunov exponent is

$$\begin{equation} (4) \hspace{1em} \lambda(r) = \lim_{n \to \infty} \left\{\frac{1}{n} \sum_{i = 0}^{n - 1} \ln \left|f'(x_i)\right|\right\} \end{equation}$$

and we can define any other $$f(x)$$ depending the finite-difference equations. For (1) the graph is: The code

I want a good approximations of this exponent and the following code plots is that way, but I want a smooth graph and it takes a lot of time.

ClearAll["Global*"]

ITERATIONS = 100; (* Fidelity *)
F[x_, r_] := r*x*(1 - x); (* Logistic map *)
f[x_, r_] := r*(1 - 2*x); (* Derivate *)
λ[r_, u_] := (1/ITERATIONS)*
Sum[Log[Abs[f[Nest[F[#, r] &, u, i], r]]], {i, 0,
ITERATIONS - 1}]; (* Liapunov exponent approximation *)

Plot[
Legended[
Style[λ[r, RandomReal[]], Black, Thickness[0.001]],
Style["Fidelity = 100", FontSize -> 20]
],
{r, 3 , 4},
PlotTheme -> "Scientific",
FrameStyle -> Directive[Thickness[0.003], FontSize -> 20, Black],
ImageSize -> Large,
AspectRatio -> 1/2,
PlotRange -> {{3, 4}, {-1, 1}},
Axes -> True,
AxesStyle -> Directive[Thickness[0.002], Blue]
]   I know this problem can be solved with a "manual Nest" because Nest itself takes a lot of time with large number of iterations. How can I improve this code? How can I use a function like a counter, but, a Nest counter?

## 1 Answer

You are using symbolic Nest which is extremely slow. Anyway, it is more efficient to use RecurrenceTable. The plot below generates in less than 2 seconds with 20000 iterations.

ClearAll["Global*"]

ITERATIONS = 20000; (*Fidelity*)
F[x_, r_] := r*x*(1. - x); (*Logistic map*)
f[x_, r_] := r*(1 - 2*x);
xlist[x0_?NumericQ, r_] := RecurrenceTable[{x[n + 1] == r*x[n]*(1 - x[n]), x == x0},
x, {n, 0, ITERATIONS - 1}]
lambda[r_, xx_?NumericQ] := Mean[xlist[xx, r] // f[#, r] & // Abs // Log]

Plot[lambda[r, RandomReal[]], {r, 3, 4}] You can superimpose with the original plot (see 18767):

img = ImageCrop@Import["https://i.stack.imgur.com/dBrkG.png"]
img = ImageTrim[img, {{60, 65}, {475, 420}}]
Plot[lambda[r, RandomReal[]], {r, 3, 4}, AspectRatio -> 1/2,
ImageSize -> 1000, PlotRange -> {{3, 4}, {-1, 1}},
PlotStyle -> Thickness[0.005],
Prolog -> {Texture[img],
Polygon[{Scaled[{0, 0}], Scaled[{1, 0}], Scaled[{1, 1}],
Scaled[{0, 1}]},
VertexTextureCoordinates -> {{0, 0}, {1, 0}, {1, 1}, {0, 1}}]}] • I had never seen before the notation ?NumericQ, can you explain it to me? I found this Reference – tajiri_numero_1 Oct 30 '20 at 17:57
• It's a Pattern, it means that the function is only defined if NumericQ returns True for the argument (= the argument is numerical value). Try for instance f[x_NumericQ] := x^2 and f then f[abc]. The latter returns unevaluated. – anderstood Oct 30 '20 at 18:11
• @tajiri_numero_1 Is it correct that the system is stable iff $\lambda <0$? $r$ would be a parameter of a dynamical system (such as mass or length)? – anderstood Oct 30 '20 at 18:15
• Yes, if $\lambda < 0$ represents a stable behavior, $\lambda > 0$ represents a chaotic behavior and $\lambda = 0$ is a marginal value. $r$ in Logistic map could represent a lot of different things and it is given by the problem. In biology, for example, $r$ could represent a parameter of birth and dead of a population. – tajiri_numero_1 Oct 30 '20 at 18:28