8

I would like to calculate the $δ$ and the $α$ Feigenbaum constants for the Logistic Map:

\begin{equation} x_{n+1}=4μx_n(1-x_n)=f(x_n) \end{equation} as part of an undergraduate assigment of mine. I have managed so far to calculate explicitly the 2-cycle and the value of $\mu$ for which it becomes unstable, leading to the bifurcation of the 4-cycle. My attempt of calculating the $δ$, which is defined as \begin{equation} δ=\lim_{n \to \infty}\frac{μ_{n}-μ_{n+1}}{μ_{n+1}-μ_{n+2}} \end{equation} is given below 

f[x_Real, mu_Real, n_Real] := Block[{y = x}, Do[y = 4 mu y (1 - y), {n}]; y]
mu /. FindRoot[f[1/2, mu, 2^2] == 1/2, {mu, 0.9}]
delta[n_] := (mu[n - 1] - mu[n - 2])/(mu[n] - mu[n - 1])
mu[n_] := mu[n] == mu /. FindRoot[f[1/2, mu, 2^n] == 
1/2, {mu, {mu[n - 1] + (mu[n - 1] - mu[n - 2])/delta[n - 1], 
mu[n - 1] + (mu[n - 1] - mu[n - 2])/(delta[n - 1] + .01)}}, 
AccuracyGoal -> 10]
Table[mu[n], {n, 0, 7}]
Table[delta[n], {n, 2, 7}]

I am literally a rookie in the Mathematica software and tbh I have searched online a lot in order to write this piece of code above. Anyway, I am getting an error of the form

$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>

once I try to calculate the various values of the $μ$-cycles. (for example: mu[3]).

I would really appreciate your assistance. Thanks!

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Spyridoula
  • 83
  • 5
  • Welcome! I suggest the following:
    1. As you receive help, try to give it too, by answering questions in your area of expertise.
    2. Take the tour and check the faqs!
    3. When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign!
     –  Apr 03 '16 at 17:08
  • Seems like you are putting "x=0.5" as your starting point, right? – MathX Apr 03 '16 at 18:08
  • @MathX That is correct yes. – Spyridoula Apr 03 '16 at 18:13
  • 2
    your function for mu[n] is recursive, but you didnt define the base case mu[1] anywhere. – AccidentalFourierTransform Apr 03 '16 at 18:36
  • @AccidentalFourierTransform But how should that be done. Is it a particular command to be used here or I just define what mu[1] should be? Also, shouldn't that be a result of the function? When I start iterating it should calculate, among the rest, also mu[1] right? Sorry, but as I stated above I am all new to this. Thank you for your help! – Spyridoula Apr 03 '16 at 18:49
  • 1
    @Spy93 I cannot help you further because I dont know anything about the Feigenbaum Constants. About your code I can say: the definition of mu[n] includes mu[n-1]. This means that if you ask Mathematica to calculate mu[10], the result will depend on mu[9] which Mathematica doesnt know its value. Therefore, it will compute mu[9] to get the value, but this depends on mu[8] which Mathematica doesnt know its value - ad infinitum. If at the begining of your code you include mu[1]=1 (for example), then the process ends when the recursion hits the base case mu[1] (hope this helps) – AccidentalFourierTransform Apr 03 '16 at 18:57
  • Check out this recent question for a related problem. – Chris K Apr 03 '16 at 22:36
  • What exactly does FindRoot[f[1/2, mu, 2^2] == 1/2, {mu, 0.9}] do here, and why's this the right thing to do? (It seems that it is, but I can't see why). – Chris K Apr 04 '16 at 02:25

1 Answers1

8

Several changes are required to obtain the desired results. First, the syntax error mu[n] == mu must be replaced by mu[n] = mu. Next, initial conditions must be provided for the recurrence in n:

mu[0] := mu /. FindRoot[f[1/2, mu, 2^0] == 1/2, {mu, 0.9}]
mu[1] := mu /. FindRoot[f[1/2, mu, 2^1] == 1/2, {mu, 0.9}]
mu[2] := mu /. FindRoot[f[1/2, mu, 2^2] == 1/2, {mu, 0.9}]

However, even with these changes, FindRoot in lines 4 - 6 of the code in the question does not evaluate mu[n - 1], etc in the initial guesses for its search. This issue can be circumvented by moving the evaluation of these quantities outside FindRoot. Altogether, the revised code becomes

f[x_Real, mu_Real, n_Real] := Block[{y = x}, Do[y = 4 mu y (1 - y), n]; y]
mu[0] := mu /. FindRoot[f[1/2, mu, 2^0] == 1/2, {mu, 0.9}]
mu[1] := mu /. FindRoot[f[1/2, mu, 2^1] == 1/2, {mu, 0.9}]
mu[2] := mu /. FindRoot[f[1/2, mu, 2^2] == 1/2, {mu, 0.9}]
delta[n_] := (mu[n - 1] - mu[n - 2])/(mu[n] - mu[n - 1])
mu[n_] := mu[n] = Module[{mu1 = mu[n - 1] + (mu[n - 1] - mu[n - 2])/delta[n - 1], 
    mu2 = mu[n - 1] + (mu[n - 1] - mu[n - 2])/(delta[n - 1] + .01)}, 
    mu /. FindRoot[f[1/2, mu, 2^n] == 1/2,  {mu, mu1, mu2}, AccuracyGoal -> 10]]
Table[mu[n], {n, 0, 7}]
Table[delta[n], {n, 2, 7}]

(* {0.5, 0.809017, 0.87464, 0.88866, 0.891667, 0.892311, 0.892449, 0.892478} *)
(* {4.70894, 4.68077, 4.66296, 4.6684, 4.66895, 4.66916} *)
bbgodfrey
  • 61,439
  • 17
  • 89
  • 156
  • 2
    How does f[1/2, mu 2^n] even evaluate, since 1/2 and 2^n aren't Real? When I run f[1/2, 0.7, 2] it doesn't evaluate, yet these FindRoots work. – Chris K Apr 04 '16 at 02:27
  • 3
    @ChrisK Interesting observation. Perhaps, FindRoot passes real numbers. In any case, the Real conditions can be eliminated, if Evaluated -> False is added to the FindRoot call. – bbgodfrey Apr 04 '16 at 04:10
  • @bbgodfrey Thank you very much for your answer. It was really helpful! – Spyridoula Apr 05 '16 at 19:53