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

KeyTo9_Fans · 发表于 2024-5-24 10:12:36

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？欢迎注册

×

问：我想使用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 脚本的示例，用于生成分叉图：

c_range = 3:0.00001:4; % c 的范围和步长
max_iterations = 10000; % 最大迭代次数
bifurcation_values = zeros(length(c_range)*16, 1); % 初始化用于绘制分叉图的数据
for i = 1:length(c_range) % 循环遍历 c 的范围
c = c_range(i); % 取出一个c
a = rand()*c; % 初始化 a(1)
for j = 1:max_iterations % 迭代差分方程
a = a * (c - a); % 计算下一项
if j > max_iterations-16 % 记录最后16项数据
bifurcation_values(i*16-max_iterations+j) = a;
end
end
end
% 绘制分叉图
figure;
plot(linspace(3,4,length(bifurcation_values)), bifurcation_values, 'b.', 'Markersize', 1);
xlabel('c');
ylabel('Stable value of a(n)');
title('Bifurcation Diagram of the Difference Equation');
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 · 发表于 2024-5-24 10:43:48

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

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

mathe · 发表于 2024-5-24 13:14:27

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

mathe · 发表于 2024-5-24 13:28:12

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

Jack315 · 发表于 2024-5-24 15:06:21

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

Jack315 · 发表于 2024-5-24 15:35:07

Mathematica 代码：
Lyapunov Exponents for the Logistic Map

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

复制代码

Jack315 · 发表于 2024-5-24 15:41:21

Mathematica 代码：
Hénon Map Visualization

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

复制代码

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

复制代码

Jack315 · 发表于 2024-5-24 16:08:23

分享：

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

KeyTo9_Fans · 发表于 2024-5-26 13:50:38

当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

账号		自动登录	找回密码
密码			欢迎注册

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

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。