# How to approximate the solution of a system of quartic equations

In an attempt to evaluate the point at which the cycle 2 becomes unstable for the given map:

$$x_{n+1}=μ-x_n^4=f(x_n), \quad μ \in \mathbb{R}$$

I have so far managed to locate the two numerical values $\hat{x}_1, \hat{x}_2$ for which the period 2 trajectory appears:

Manipulate[Module[{list = NestList[μ - #^4 &, x0, 100]}, list2 = list;
Column[{ListLinePlot[list, PlotRange -> {-1, 1.5},
ImageSize -> {450, 375}],
TableForm[Transpose@{Range[86, 101], list[[-16 ;;]]},
TableHeadings -> {None, {"point", "x"}}]}]], {{μ, 0.2,
"parameter μ"}, 0, 4,
Appearance -> "Labeled"}, {{x0, 0.4,
"Initial \!$$\*SubscriptBox[\(x$$, $$0$$]\)"}, 0, 1,
Appearance -> "Labeled"}]


But now I would like to locate the value of $μ$ for which this 2-cycle ceases to be stable and jumps onto the next one, a 4-cycle. The condition which has to hold is the following: $$f'(\hat{x}_1)f'(\hat{x}_2)=-1 \Leftrightarrow 16\hat{x}_1^3 \hat{x_2}^3=-1$$ where $\hat{x}_i,i=1,2$ would be an expression of $μ$.

The problem is that Mathematica is not able to solve analytically the system: $$\begin{cases} f(\hat{x}_1)=\hat{x}_2 \\ f(\hat{x}_2)=\hat{x}_1 \end{cases}$$ in order to get the trajectory of the 2-cycle in a closed form, which will depend only on the parameter $μ$, which I can then plug it in the condition and find the value of that $μ$. What I get instead for

Reduce[{y == μ - z^4, z == μ - y^4}, {y, z}, Reals]


is the solution in Root form:

 ((μ == Root[27 + 256 #1^3 &, 1] &&
y == Root[-μ + #1 + #1^4 &, 1]) || (Root[27 + 256 #1^3 &, 1] <
μ <= Root[-125 + 256 #1^3 &,
1] && (y == Root[-μ + #1 + #1^4 &, 1] ||
y == Root[-μ + #1 + #1^4 &, 2])) || (μ >
Root[-125 + 256 #1^3 &, 1] && (y == Root[-μ + #1 + #1^4 &, 1] ||
y == Root[
1 - μ^3 - μ^2 #1 - μ #1^2 - #1^3 + 3 μ^2 #1^4 +
2 μ #1^5 + #1^6 - 3 μ #1^8 - #1^9 + #1^12 &, 1] ||
y == Root[-μ + #1 + #1^4 &, 2] ||
y == Root[
1 - μ^3 - μ^2 #1 - μ #1^2 - #1^3 + 3 μ^2 #1^4 +
2 μ #1^5 + #1^6 - 3 μ #1^8 - #1^9 + #1^12 &, 2]))) &&
z == -y^4 + μ


I understand that I am not able to solve this system analytically, but how can I then approximate this particular $μ$ value.

If I plug in the condition above the numerical values of $\hat{x}_1,\hat{x}_2$, I would have no $μ$ involved.

The point of this whole procedure is to calculate somewhow the universal $\delta$ Feigenbaum constant for this representative of the quartic equivalence class of discrete mappings.

Any help would be greatly appreciated. Thanks!

• Replace the = with == in your Reduce and also check documentation for those respective symbols. Mar 10, 2016 at 22:22
• @PlatoManiac Thank you, you were right. But still, how can I then plug in those non-analytical roots into the condition mentioned above? Sorry for my perhaps naive questions I am still new to Mathematica, trying my best to make decent questions here. Mar 10, 2016 at 22:38
• p = Reduce[{y == \[Mu] - z^4, z == \[Mu] - y^4}, {y, z, \[Mu]}, Reals, Backsubstitution -> True] // FullSimplify // N; First@p Mar 10, 2016 at 23:48
• @Dr.belisarius Thank you for your help! But I cannot understand fully what did you do there. Mathematica returned to me the value $μ=6.29$. Is that the $μ$ which I would get if I knew the exact expression of $\hat{x}_1,\hat{x}_2$ and solved for $μ$ the condition: $16 {\hat{x}_1}^3{\hat{x}_2}^3=-1$ ? Mar 10, 2016 at 23:58
• @Mitscaype p = Reduce[{y == \[Mu] - z^4, z == \[Mu] - y^4}, {y, z, \[Mu]}, Reals, Backsubstitution -> True] // FullSimplify // N; Solve[p[[1]]] Mar 11, 2016 at 0:13

Let's start with a brute-force bifurcation diagram to see what's happening, modifying some code from @bbgodfrey's answer to this question:

f[x_, μ_] := μ - x^4;

res[μ_] :=
DeleteDuplicatesBy[
NestList[f[#, μ] &, Nest[f[#, μ] &, μ, 1000], 100],
Round[#, 10^-6] &];

ListPlot[Catenate[Table[{μ, #} & /@ res[μ], {μ, 0, 1.2, 0.001}]],
PlotStyle -> PointSize[Tiny], PlotRange -> All]


Looks like the μ you're looking for is near 1.1. To find a 2-cycle for a given μ numerically, we can try

μ = 1.05;
FindRoot[f[f[x, μ], μ] == x, {x, -0.5}]

(* {x -> -0.162297} *)


The bifurcation occurs when the two-cycle goes through a period-doubling bifurcation (eigenvalue=-1). We can numerically find the the bifurcation point by simultaneously solving for the two-cycle and when it loses stability:

Clear[μ];
FindRoot[{f[f[x, μ], μ] == x, D[f[f[x,μ], μ], x] == -1}, {x, -0.5}, {μ, 1.1}]

(* {x -> -0.359832, μ -> 1.11964} *)


Edit:

To generalize to higher order cycles (e.g. here's where the 4-cycle gives way to an 8-cycle):

FindRoot[{Nest[f[#, μ] &, x, 4] == x, D[Nest[f[#, μ] &, x, 4], x] == -1},
{x, -0.64}, {μ, 1.16}]

(* {x -> -0.646515, μ -> 1.16166} *)


Of course it gets more delicate, so you'll probably want better initial guesses (which I got from looking at a zoomed-in version of the bifurcation diagram).

I don't know of any way to directly calculate the Feigenbaum constant here.

Edit 2:

I applied this approach to the first 7 period-doublings. The ratio of bifurcation points as in this wikipedia page looks like {7.90629, 7.25658, 7.3238, 7.28859, 7.28785}.

• Actually that is pretty correct. I was guessing the bifurcation values of $μ$ from a $x_{n+1}$ vs $n$ plot and I was around 1.13 for a 4-cycle. But to locate the 8-cycle is pretty hard this way. Is there also a code for the Feigenbaum $\delta$ in this case? Mar 11, 2016 at 0:43
• I edited my answer to give code to generalize to finding the bifurcation of an n-cycle. I also tried calculating that δ for you: looks like around 7.288 if I didn't mess that up. Mar 11, 2016 at 2:49
• Found this article which says it's 7.28468..., which makes me feel OK. Mar 11, 2016 at 3:04
• Well, I guess I will have to thank like, forever! Yeah, I was able to approximate that $\delta$ too, but only by the first 2 iterations. I was around 8.28, which means wayyyyy out when it comes to this kind of systems. Again, many thanks. Mar 12, 2016 at 20:23
• Not a bad paper either. I guess I had missed that somehow. Again, thanks Mar 12, 2016 at 20:25