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

KeyTo9_Fans

问：我想使用Matlab绘制一张差分方程极限值的分叉图，横坐标是c，范围是3到4，纵坐标是{a(n)}稳定后的2^k个极限值，k可能是1、2、3、……、inf。{a(n)}的递推方程是a(1)=rand()*c，a(i+1)=a(i)*(c-a(i))。据说这个简单的差分方程能制造出分叉和混沌现象，我想把这个现象画出来。

AI答：要绘制差分方程 a(i+1) = a(i) * (c - a(i)) 的分叉图，并观察随着参数 c 的变化 a(n) 的极限行为，你可以按照以下步骤在 MATLAB 中实现：
初始化参数 c 的范围和步长。
对于每个 c 值，初始化随机种子（可选，以确保结果的可重复性）并迭代差分方程足够多次，直到达到稳定状态或达到某个最大迭代次数。
记录 a(n) 的稳定值（或迭代过程中的最后几个值）。
将这些稳定值与对应的 c 值一起绘制在图上。
下面是一个 MATLAB 脚本的示例，用于生成分叉图：

1. c_range = 3:0.00001:4; % c 的范围和步长  
   
2. max_iterations = 10000; % 最大迭代次数  
   
3. bifurcation_values = zeros(length(c_range)*16, 1);  % 初始化用于绘制分叉图的数据  
   
4. for i = 1:length(c_range)  % 循环遍历 c 的范围  
   
5.     c = c_range(i);  % 取出一个c
   
6.     a = rand()*c; % 初始化 a(1)  
   
7.     for j = 1:max_iterations % 迭代差分方程  
   
8.         a = a * (c - a); % 计算下一项
   
9.         if j > max_iterations-16 % 记录最后16项数据
   
10.             bifurcation_values(i*16-max_iterations+j) = a;  
    
11.         end  
    
12.     end  
    
13. end  
    
14. % 绘制分叉图  
    
15. figure;  
    
16. plot(linspace(3,4,length(bifurcation_values)), bifurcation_values, 'b.', 'Markersize', 1);  
    
17. xlabel('c');  
    
18. ylabel('Stable value of a(n)');  
    
19. title('Bifurcation Diagram of the Difference Equation');  
    
20. grid on;

复制代码

这个图画出来还是挺不错的，符合我的预期：

但是追问的回答我就不满意了，最终还得人工求解……

追问：如何求出c取什么值的时候是2叉、取什么值的时候是4叉、取什么值的时候是8叉、取什么值的时候是16叉、……、取什么值的时候是2^k叉？

答：要确定参数c取什么值时，差分方程a(i+1) = a(i) * (c - a(i))的解会展现出2^k叉的分叉现象，我们需要分析这个差分方程的稳定性和不动点。
首先，不动点是满足a = a * (c - a)的a值。解这个方程，我们得到两个不动点：
a = 0（对所有c都成立）
a = c（但只有当c在[0, 1]区间内时，这个不动点才是有效的，因为a(i)必须在[0, c]内）
对于混沌和分叉现象，我们主要关注c > 1的情况，因为当c <= 1时，差分方程只有一个有效的不动点a = 0，并且没有分叉现象。
分叉现象通常发生在不动点失去稳定性的时候。为了分析稳定性，我们需要考虑差分方程在不动点附近的线性化。但是，由于这个差分方程是非线性的，我们不能直接应用线性化的方法。相反，我们需要通过数值方法或图形分析来观察分叉现象。
然而，对于某些特定的c值，我们可以观察到2^k叉的分叉现象。这些值通常与差分方程的周期轨道有关。但是，要精确找出这些值，通常需要数值方法或符号计算工具。
不过，我们可以给出一些一般性的指导：
当c接近某个值时，如果差分方程的解开始展现出多个稳定的周期轨道，那么可能发生了分叉现象。
分叉的数量（即2^k）通常与周期轨道的数量和稳定性有关。当一个新的稳定周期轨道出现时，分叉的数量可能会翻倍。
为了找到这些特定的c值，你可以使用数值方法（如二分法、牛顿法等）来搜索使差分方程展现出特定分叉数量的c值范围。
请注意，由于这个差分方程的复杂性，找到精确的c值可能是困难的，而且可能需要大量的计算。因此，在实际操作中，你可能需要使用数值工具或编程来辅助你的分析。

KeyTo9_Fans

把c=[3.562, 3.572]这个区间放大，可以发现分叉越来越快，下一个分叉区间的长度仅为上一个分叉区间的长度的1/4，如下图所示：

据此应该能大致估算出2叉、4叉、8叉、……、2^k叉所在的区间

mathe

这个属于混沌理论，没有极限，只有混沌

mathe

https://blog.csdn.net/weixin_42207434/article/details/132527164

Jack315

非线性科学丛书《从抛物线谈起——混沌动力学引论》(作者)郝柏林 上海科技教育 1993年9月第1版.pdf
https://pan.baidu.com/share/link ... ;_at_=1716534276748

Jack315

Mathematica 代码：
Lyapunov Exponents for the Logistic Map

1. Manipulate[
   
2. Column[{
   
3.    ListPlot[
   
4.     Flatten[logisticMapOrbitCF[#, n] & /[url=home.php?mod=space&uid=6175]@[/url]
   
5.       Range[a[[1]] - .04/magH, Min[1, a[[1]] + .04/magH],
   
6.        ControlActive[.6, 0.15]/magH/pts], 1],
   
7.     PlotStyle -> {RGBColor[.25, .43, .82], AbsolutePointSize[1]},
   
8.     PlotLabel -> Style["Bifurcation Diagram", 9],
   
9.     PlotRange -> {{a[[1]] - .04/magH,
   
10.        a[[1]] + .04/magH}, {a[[2]] - .52/magV, a[[2]] + .5/magV}},
    
11.     AxesOrigin -> {a[[1]] - .04/magH, Automatic},
    
12.     ImagePadding -> {{15, 12}, {20, 6}}, AspectRatio -> .65,
    
13.     ImageSize -> 500],
    
14.    ListPlot[{#, lyapListCF[#, n]} & /@
    
15.      Range[a[[1]] - .04/magH, Min[1, a[[1]] + .04/magH],
    
16.       ControlActive[.6, 0.15]/magH/pts],
    
17.     PlotLabel -> Style["Lyapunov Exponent", 9], Joined -> True,
    
18.     PlotStyle -> RGBColor[.49, 0, 0],
    
19.     PlotRange -> {{a[[1]] - .04/magH,
    
20.        a[[1]] + Min[1, .04/magH]}, {Full, .5}},
    
21.     AxesOrigin -> {a[[1]] - .04/magH, 0}, Ticks -> {Automatic, {0}},
    
22.     ImagePadding -> {{16, 12}, {12, 0}}, AspectRatio -> .14,
    
23.     ImageSize -> 500]}, Spacings -> 0, Frame -> True],
    
24. {{a, {.85, 0.5}}, {0.7, 0}, {1, 1}, Locator,
    
25.   Appearance -> Style["\[CirclePlus]", RGBColor[.49, 0, 0], Bold, 18]},
    
26. Style["zoom", Bold],
    
27. "horizontal",
    
28. {{magH, .26667, ""}, .1, 10, .1, VerticalSlider, ImageSize -> Small},
    
29. Delimiter,
    
30. "vertical",
    
31. {{magV, 1, ""}, .1, 10, .1, VerticalSlider,
    
32.   ImageSize -> Small}, Delimiter,
    
33. {{n, 1000, "number of iterates"}, ControlType -> None},
    
34. {{pts, 2000, "number of \[Lambda] values"}, ControlType -> None},
    
35. Button["reset", a = {.85, 0.5}; magH = .26667; magV = 1; n = 1000;
    
36.   pts = 2000],
    
37. Button["detail", n = 2000; pts = 4000],
    
38. ControlPlacement -> Left,
    
39. AutorunSequencing -> {{1, 10}, {2, 5}, {3, 5}},
    
40. TrackedSymbols :> {a, magH, magV, n, pts},
    
41. Initialization :> (
    
42.    logisticMapOrbitCF =
    
43.     Compile[{{\[Lambda], _Real}, {n, _Integer}}, {\[Lambda], #} & /@
    
44.       Take[NestList[4. \[Lambda] # (1. - #) &, \[Lambda], n], -24]];
    
45.    lyapListCF =
    
46.     Compile[{{\[Lambda], _Real}, {n, _Integer}},
    
47.      Mean[(Log[Abs[4. \[Lambda] (1. - 2 #)]] & /@
    
48.         NestList[4. \[Lambda] # (1. - #) &, \[Lambda], n])]]
    
49.    )]

复制代码

Jack315

Mathematica 代码：
Hénon Map Visualization

1. Henon[\[Alpha]_, \[Beta]_, x0_, y0_, n_] := Module[{h0},
   
2.    h0 [{x_, y_}] := {1 - \[Alpha] x^2 + y, \[Beta] x};
   
3.    NestList[h0, {x0, y0}, n]
   
4.    ];

复制代码

1. Manipulate[
   
2. ListPlot[Henon[\[Alpha], \[Beta], 0.0, 0.0, 20000], Axes -> False,
   
3.   Background -> GrayLevel[0.05],
   
4.   PlotStyle -> {LightOrange, Opacity[opacity], PointSize[pointSize]},
   
5.   ImageSize -> {520, 380}],
   
6. {{\[Alpha], 0.0973}, 0.0001, 0.1, Appearance -> "Labeled"},
   
7. {{\[Beta], 0.99973}, 0.99, 0.99999, Appearance -> "Labeled"},
   
8. {{opacity, 0.45}, 0.1, 0.9, Appearance -> "Labeled"},
   
9. {{pointSize, Tiny, "point size"}, {Tiny -> "tiny", Small -> "small"}},
   
10. SaveDefinitions -> True]

复制代码

Jack315

分享：
feigenbaum.rar (169.61 KB, 下载次数: 2)

2024-5-24 16:08 上传

点击文件名下载附件

KeyTo9_Fans

当c=3时，1叉分成2叉，分叉位置是(3,2)

当c=3.44948974278317809819728407470589时，2叉分成4叉，分叉位置是(c,1.51763809020504152469779767524810)和(c,2.93185165257813657349948639945779)

这个c的精确值是1+sqrt(6)

当c=3.5440903595519228536159659866...时，4叉分成8叉，这个常数的前102位对应这个数列：

https://oeis.org/A086181

当c=3.5644072660954325977735575865...时，8叉分成16叉，这个常数的前99位对应这个数列：

https://oeis.org/A091517

而16叉→32叉、32叉→64叉……的精确分叉位置，目前我还没找到，粗糙的数据如下表所示：



表格出处：https://wiki.swarma.org/index.php/Logistic%E6%98%A0%E5%B0%84

当c=3.56994567187...时，分成∞叉，这个常数目前在OEIS里已经延长到94位了，链接如下：

https://oeis.org/A098587

Jack315

【说明】
这个是分享暗线方程，不是求分岔点的参数值。

【暗线】
按 5# 的教材选择抛物线映射函数为：
\(x_{n+1}=1-\mu x_n^2\)
参变量：\(0<\mu\le 2\)
相变量：\(-1\le x\le 1\)
下图为此映射函数相应的分岔图：


分岔图中的暗线方程在 Mathematica 里生成方法：
(1) 定义映射函数

1. f[\[Mu]_, x_] := 1 - \[Mu] x^2

复制代码

(2) 取此映射唯一的超稳定点 \(P_0(\mu)=1\) 作为迭代初值。
(3) 用映射函数进行迭代：\(P_{n+1}(\mu)=f[\mu, P_n(\mu)]\) 。
迭代代码：

1. p0 = Collect[f[\[Mu], 0], \[Mu]]
   
2. p1 = Collect[f[\[Mu], p0], \[Mu]]
   
3. p2 = Collect[f[\[Mu], p1], \[Mu]]
   
4. p3 = Collect[f[\[Mu], p2], \[Mu]]
   
5. p4 = Collect[f[\[Mu], p3], \[Mu]]
   
6. p5 = Collect[f[\[Mu], p4], \[Mu]]

复制代码

下图为相应的暗线图：

1. Plot[{p0, p1, p2, p3, p4, p5},
   
2. {\[Mu], 0, 2},
   
3. AspectRatio -> Automatic,
   
4. PlotLegends -> "Expressions",
   
5. PlotLabel -> "暗线",
   
6. GridLines -> Automatic,
   
7. Frame -> True]

复制代码

迭代得到一族多项式函数，其项数随着迭代过程迅速爆增。