绘制出分叉和混沌图后，如何精确求解分叉位置？

Jack315

【周期 6】
迭代方程：

1. f6[r_, x_] := ((Nest[f[r, #] &, t, 6] - t) (Nest[f[r, #] &, t, 1] - t))/((Nest[f[r, #] &, t, 3] - t) (Nest[f[r, #] &, t, 2] - t)) /. t -> x // Expand // Simplify

复制代码

判别式：

1. d6 = Factor[Discriminant[f6[r, x], x]]

复制代码

判别式根：

1. Solve[d6 == 0, r]

复制代码

周期 6 分岔点参数：\(r_6=3.626553162, 3.841499008, 3.937516419, 3.977760441, 3.997582524\) 。

Jack315

周期数再提高就开始烧 CPU 了

这是计算结果汇总：


表中的符号序列参考 5# 的教材 47 页表 10.1 周期 7 以内的揉序列。
或许用“字提升法”可以计算出周期数更高的分岔点参数……努力学习中。

【一点体会】

• 西方文明在数学和自然科学的积累确实丰富；
• 中国的科学家真了不起。
  

Jack315

【字提升法：超稳定不动点处的参数值】
超稳定不动点是指 n-周期各轨道点处导数的乘积为零：\(\prod_{i=1}^n\frac{df(\mu,x)}{dx}\vert_{x=x_i}=0\)
这意味着轨道中必然含有导数为零的点 C：\(\frac{df(\mu,x)}{dx}\vert_{x=C}=0\)

具体步骤：
1. 由初始符号序列（参考 22# 表格）生成倍周期的符号序列\(\Sigma\)。
2. 求出 \(C\) 点的 \(x\) 值以及 \(f(\mu,C)\) 的值。
3. 求出映射函数 \(y=f(\mu,x)\)的两支单调逆函数 \(x_L=f_L^{-1}(\mu,y)\) 和 \(x_L=f_L^{-1}(\mu,y)\) 。
4. 根据符号序列\(\Sigma\)调用相应的逆函数进行迭代 \(f(C)=\Sigma(C)\)，得到超稳定不动点的参数值 \(\mu\) 。

映射函数 \(f(\mu,x)=x(\mu-x), \mu\in[0,4],x\in[0,4]\) （参考 15# 的分岔图）的迭代函数：

1. iterParam[seq_, pwr_, r0_, ndigits_, maxIter_] := Module[
   
2.   (*【输入参数】*)
   
3.   (* seq: 轨道 p（质数）的初始符号序列。*)
   
4.   (* pwr: 轨道 p\[Times]2^pwr 的幂次。*)
   
5.   (* r0: 映射函数中参量值的迭代初值*)
   
6.   (* ndigits: 数据精度，即小数点后的位数。*)
   
7.   (* maxIter: 允许的最大迭代次数。*)
   
8.   
   
9.   {n = ndigits + 10, ks, len, err, iters, r, rlast, status, y},
   
10.   (*【中间变量】*)
    
11.   (* n: 计算过程中的数据精度。*)
    
12.   (* ks: 生成的符号序列列表。*)
    
13.   (* len: 【】列表 ks 的长度，即轨道周期数。*)
    
14.   (* err: 【】两次迭代之间的误差*)
    
15.   (* iters: 【】迭代次数。*)
    
16.   (* r: 【】迭代得到的参量值。*)
    
17.   (* rlast: 前一次迭代得到的参量值。*)
    
18.   (* status: 【】迭代次数是否超过充许次数的状态。*)
    
19.   (* y: 映射函数逆函数的输入值。*)
    
20.   
    
21.   (*【返回参数】*)
    
22.   (* status, iters, len, r, err*)
    
23.   
    
24.   (*【生成符号序列列表】*)
    
25.   ks = Characters[Nest[StringReplace[#, {"R" -> "RL", "C" -> "RC", "L" -> "RR"}] &, seq, pwr]];
    
26.   len = Length[ks];
    
27.   
    
28.   (*【进行迭代】*)
    
29.   err = SetPrecision[4, n];
    
30.   iters = 0;
    
31.   r = SetPrecision[r0, n];
    
32.   rlast = r;
    
33.   status = "迭代成功结束";
    
34.   While[err > 10^-n,
    
35.    For[i = Length[ks], i > 0, i--,
    
36.     SetPrecision[Switch[ks[[i]],
    
37.        "R", y = 1/2 (r + Sqrt[r^2 - 4 y]),
    
38.        "L", y = 1/2 (r - Sqrt[r^2 - 4 y]),
    
39.        "C", y = r/2], n];
    
40.     ]; (*End For*)
    
41.    r = SetPrecision[2 Sqrt[y], n];
    
42.    
    
43.    (*【更新数据】*)
    
44.    err = SetPrecision[Abs[r - rlast], n];
    
45.    rlast = r;
    
46.    
    
47.    (*【检查迭代次数】*)
    
48.    iters++;
    
49.    If[iters > maxIter, status = "超出最大迭代次数"; Break[]];
    
50.    ]; (*End While*)
    
51.   
    
52.   (*【返回迭代结果】*)
    
53.   {r, err} = SetPrecision[{r, err}, ndigits];
    
54.   {status, iters, len, r, err}
    
55.   ] (*End Module*)

复制代码

函数调用示例：\(3\times2^5=96\) 周期，初始参数值 3.99，16 位精度，最多迭代 2W 次。

1. iterParam["RLC", 5, 3.99, 16, 20000]

复制代码

运行结果：
\(\begin{array}{|c|c|c|c|c|}
\hline \text{序号}(n)&\text{周期}2^n&2^n\text{超稳定不动点参数}(\mu_n)&\text{周期}3*2^n&3*2^n\text{超稳定不动点参数}(\mu_n)\\
\hline 0&1&2.000000000000000&3&3.831874055283316\\
\hdashline 1&2&3.236067977499790&6&3.627557529515523\\
\hdashline 2&4&3.498561699327702&12&3.582229835820361\\
\hdashline 3&8&3.554640862768825&24&3.572577293246184\\
\hdashline 4&16&3.566667379856269&48&3.570509238189459\\
\hdashline 5&32&3.569243531637110&96&3.570066370442875\\
\hdashline 6&64&3.569795293749945&192&3.569971521752127\\
\hdashline 7&128&3.569913465422349&384&3.569951208122691\\
\hdashline 8&256&3.569938774233305&768&3.569946857566441\\
\hdashline 9&512&3.569944194608065&1536&3.569945925810607\\
\hdashline 10&1024&3.569945355486469&3072&3.569945726257044\\
\hdashline 11&2048&3.569945604111078&6144&3.569945683518782\\
\hdashline 12&4096&3.569945657358856&12288&3.569945674365555\\
\hline \end{array}\)

两个分岔点之间夹一个超稳定点；两个超稳定点之间夹一个分岔点。
因此分岔点参数与超稳定点参数具有相同的极限 \(\mu_\infty\)
根据 A098587 的结果，\(\mu_\infty=3.5699456718709449...\) （小数点后 94 位）。
这里的计算结果当 n=12 时，有小数点后 7~8 位的精度。

\(\mu_\infty\) 是一个分水岭，在 \(\mu<\mu_\infty\) 的区间是 \(2^n\) 周期的轨道；
所有其他周期以及混沌轨道，则在 \(\mu>\mu_\infty\) 的区间。

求解 \(\mu_\infty\) 的值真不是件容易的事，再好的 CPU + 算法最多也只能提高点精度。
由于当 \(n\rightarrow\infty\) 时符号序列长度也按指数趋于无穷，
所以数值方法计算的结果无法得到很高精度的 \(\mu_\infty\) 值。

有点好奇这些牛人是如何计算到小数点后 94 位的。
只要有 CPU 烧外加足够的时间，字提升法可以算到足够高的精度。
但这未免有点笨 (Bao) 拙 (Li) ，是否还有更好的办法呢？

Jack315

来个比较对称的分岔图。
映射函数：
\(f(r, x)=rx(1-x^2)\)
其中：
\(0<r\le 3\)
\(-\sqrt{1+\frac{1}{r}}\le x \le \sqrt{1+\frac{1}{r}}\)

映射函数图：

分岔图：

周期 3 是混沌的必由之路。
这个图里周期 6 比较明显，周期 3 却不容易看到，会在哪里呢？

Jack315

【23# 相关补充和说明】
映射函数为 \(f(r, x)=x(r-x)\)。
\(\frac{df(r,x)}{dx}=r-2x\) 。令其为 0 解得 \(C:x=\frac{r}{2}\)
代入到映射函数得：\(f(r,C)=(\frac{r}{2})(r-\frac{r}{2})=\frac{r^2}{4}\)

求 \(y=x(r-x)\) 的反函数得：\(x=\frac{1}{2}(r\pm\sqrt{r^2-4y})\)

迭代函数 \(f(r, C)=\Sigma(C)\rightarrow\frac{r^2}{4}=\Sigma(\frac{r}{2})\rightarrow r=2\sqrt{\Sigma(\frac{r}{2})}\)
其中：\(\Sigma\) 为由 \(\{L,R\}\) 组成的符号序列。

说明：
代码中39行：\(y=C=\frac{r}{2}\)
代码中37、38行：使用反函数 \(x=\frac{1}{2}(r\pm\sqrt{r^2-4y})\) 计算符号序列 \(\Sigma\) 中 \(L,R\) 。
代码中41行：由字符序列迭代结果计算 \(r=2\sqrt{...}\) 。