Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
what the is the meaning of that "difference(i+1)/difference[i+1]"? It looks like you don't need the + 1 at all. –  amr Jan 31 '13 at 4:34
    
I'd recommend browsing the classic examples at the demo site. –  Jens Jan 31 '13 at 4:35
    
For starters, replacing Do by NestList will speed up the creation of x by more than an order of magnitude. To find periodic orbits, consider a Fourier analysis of the data. (How do you know it must have a period equal to a power of $2$? That seems unlikely.) Instead, it appears this dynamical system has eight discrete limiting points. What property of this system are you really looking for? –  whuber Jan 31 '13 at 4:41
    
No offense, but I think this question is too localized, although @MarkMcClures advice is very sound. –  Yves Klett Jan 31 '13 at 20:00
add comment

closed as too localized by Yves Klett, Oleksandr R., Simon Woods, Sjoerd C. de Vries, rm -rf Feb 1 '13 at 14:22

This question is unlikely to help any future visitors; it is only relevant to a small geographic area, a specific moment in time, or an extraordinarily narrow situation that is not generally applicable to the worldwide audience of the internet. For help making this question more broadly applicable, visit the help center.If this question can be reworded to fit the rules in the help center, please edit the question.

1 Answer

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.

share|improve this answer
    
Manipulate[ListPlot[NestList[f, seed, 3000], PlotStyle->PointSize[Tiny]], {seed, 0, 1}] for one way of visualizing –  amr Jan 31 '13 at 19:10
    
There being precisely one attractive orbit was exactly the point of the code. My goal was to use values close to the accumulation point for a and I thought it would be easier to use higher precision transparently in Mathematica than in MatLab. I've ended up just working on an implementation in FORTRAN. –  Jack Henahan Feb 4 '13 at 0:25
    
@JackHenahan Given that there is a single attractive orbit, higher precision is really not particularly necessary here. Also, Mathematica's automatic error precision tracking tends to be overly pessimistic in this situation. It would make a lot more sense to use high precision numbers in a case like $a=4$. If you can modify your question to clearly indicate what you hope to accomplish using high precision numbers, I might be able to assist on that issue. (I didn't notice any mention of precision in your original post, other than the one instance of WorkingPrecision in the Mathematica code.) –  Mark McClure Feb 4 '13 at 1:24
    
@MarkMcClure The case given in the code is just a testing value, where the attractive orbit is known. I am attempting to use it for values close to the accumulation point. I will edit to reflect this. –  Jack Henahan Feb 4 '13 at 23:51
    
@JackHenahan As my latest edit shows, you don't need high precision to find reasonably long attractive orbits. –  Mark McClure Feb 5 '13 at 16:11
show 3 more comments

Not the answer you're looking for? Browse other questions tagged or ask your own question.