Calculating orbital period for the logistical map [closed]

Question

I have the following MatLab code

clear all
close all

format long     % print all sig figs unless instructed otherwise

%% initial stuff
itotal = 10^8;      % total number of iterates
x = zeros(itotal,1);    % fill orbit with zeros
x(1) = rand;        % random initial condition
a = 1+sqrt(6)+.1;   % parameter of logistic map

%% iterate the map for a long time
for i = 1:itotal-1
  x(i+1) = a*x(i)*(1-x(i));     % logistic map, g(x)=ax(1-x)
end

T = 35;         % total points to print
%% print the last T points of the orbit
disp(' ')
disp('  Iterate        X')
disp('---------------------------------------------')
for i = itotal-T+1:itotal                   % only showing T points on the orbit
  disp(sprintf('   %10.0f  %15.15f',i,x(i)))
end
disp(' ')

%% determine if there is a periodic orbit of period less than 2^S
S = 26;                 % length of longest periodic orbit (log_2) we can find, maybe?
tol = 10^(-50);             % tolerance within which I believe I've found a repeat
inside=0;               % trigger telling me if I've been inside the 'if' loop below
for i = 0:S             % loop looking for orbit of period 2^(i)
  difference(i+1) = norm(x(end)-x(end-2^i));        % compute difference between orbit points
  if (difference(i+1) < tol) && (inside == 0)       % am I within tolerance? have I already qualified?
    disp(sprintf('This orbit repeats every %s iterates',int2str(2^i)))
    inside=1;       % don't come back inside this loop, see what happens without this
  end
end

In this code, the current value of a has a known attractive orbit.

I am attempting to test values of a close to the accumulation point of the logistic map, so I need to use higher-precision arithmetic to work with this value, specifically to find orbits of period greater than $2^{20}$, my current record with the MatLab code.

I have translated this code to what I think is equivalent Mathematica code.

itotal = 10^8;(* number of iterates *)

x = ConstantArray[0, itotal]; (* fill orbit with zeros *)

x[[1]] = RandomReal[WorkingPrecision -> 10]; (* random initial value *)
a = 1 + Sqrt[6] + .1 (* parameter of logistic map *)

(* iterate the map for a long time*)

Do[x[[i + 1]] = a*x[[i]]*(1 - x[[i]]), {i, itotal - 1}];

(* Pretty Printing *)
Print["Iterate           X"]
Print["-------------------"]
Do[Print[x[[i]]], {i, itotal - 35, itotal, 1}]

(* determine if there is a periodic orbit of period less than 2^S *)

difference[i_+1] := Norm[x[[itotal]] - x[[itotal - 2^i]]];
S = 26; (* length of longest periodic orbit (log_ 2) we can find, maybe? *)
tol = 10^-50; (* tolerance within which I believe I've found a repeat *)
inside = 0; (* trigger telling me if I've been inside the'if' loop below *)
(* loop looking for orbit of period 2^(i) *)
Do[difference[i]
  If[(difference[i] < tol) && (inside == 0), 
    {Print["This orbit repeats every " <> IntegerString[2^i] <> " iterates."],
    inside = 1}], {i, 0, S, 1}]

This retains the known a for testing.

My question has 2 parts:

1. Is this translation accurate? I'm mostly questioning the n+k patterns and the use of lists over vectors, but there may be other things I'm missing. I'm not especially familiar with MatLab, so I did a naïve translation.

2. Is there anything I can do to make this implementation more efficient? It currently runs in about 10 minutes, but I'd like to speed it up, if possible.

Answer 1

If you want to translate Matlab code into Mathematica, my advice is - don't! As programming languages, the two are very different and an idiom that works well in one is unlikely to work well in the other.

A fundamental theorem theorem in discrete dynamics states that if there's an attractive orbit, then it must attract at least one critical point. Thus, if your objective is to find an attractive orbit of the logistic map, then it makes sense to start at the lone critical point of $0.5$, rather than a random point.

There's really no need for super high precision. The parameter $a=3.569945$ agrees with the accumulation point to the millionths place and the following code finds an orbit of period 512 after about 37000 iterates.

f[x_] = 3.569945 x (1 - x);
i = 0;
orbit = Monitor[
  NestWhileList[(i++; f[#]) &, 0.5, Unequal, {1, 1025}, 50000],
i];

Here's the period.

{period} = DeleteDuplicates[
  Flatten[Differences[Position[Chop[orbit - Last[orbit]], 0]]]]

(* Out: {512} *)

Let's check it.

Nest[f, Last[orbit], period] - Last[orbit]

(* -5.60663*10^-14 *)

We expect this to work, since we're searching for stability. Thus, again, high precision is unnecessary. Of course, for parameters slightly larger than the accumulation point, we expect no attractive orbits.