基于Matlab对费根鲍姆常数和费根鲍姆图的一些算法分析

最新推荐文章于 2024-04-21 11:03:48 发布

程序菜鸟一只

最新推荐文章于 2024-04-21 11:03:48 发布

阅读量3.2k

点赞数 3

文章标签： matlab 矩阵算法

本文链接：https://blog.csdn.net/sbsbsb666666/article/details/122636053

版权

先说明一下，本文主要有两部分内容，一部分是逻辑斯蒂映射，这部分我主要参考了费根鲍姆常数4.669…9(上)

这一篇文章，并且使用matlab编写代码实现了里面的一些内容，另外第二块就是第一个文件的5之后的使用我自己想出来的方法去绘制费根鲍姆图的理论收敛解（有些词语比如“理论双根解”是我造的哈，别太在意），当然，在我们的课堂上，讲费根鲍姆常数的另外几个同学还用python使用迭代散点图绘制实现了费根鲍姆图的动态呈现，这个我会把拍的她们的python迭代绘制代码放在文章中（这个我确实没有试她们的python代码）

先给出费根鲍姆常数大概是 $\delta = 4.66920.....$

当然了，因为我主要把教程写在实时脚本里面了，所以看的人最好还是用我的百度网盘分享下载我的两个脚本。其中本文的公式讲解图主要来自第一个脚本:

（MATLAB:R2019a）

脚本就放在这里了，虽然这个脚本没有把费根鲍姆常数算出来，但是希望后人能通过这个继续研究下去，可以获得更多关于费根鲍姆常数的一些认识和做出更多的东西，这个脚本里面也画出了两个解的理论分叉图

链接：点击这个获取完整脚本的网盘分享
提取码：s793

我就先把另外几个同学的python迭代代码放上来了（他们做的ppt上面的代码，用手机拍的不太清，这个我具体没研究过）。具体费根鲍姆常数是什么我们后面讲

下面我们开始讲讲费根鲍姆常数的由来 (下面的基本都是脚本内容，看脚本的下面可以不看)

一、逻辑斯蒂方程

逻辑斯蒂方程是一个描述种群增长的方程，我们首先从最初逻辑斯蒂方程的数学建模模拟开始，介绍逻辑斯蒂方程的建模过程：

这就是逻辑斯蒂方程微分方程的数学模型

下面给出一个逻辑斯蒂方程的求解示例

% 同时使用微分方程求解和递推的数值解法
% --------方法一、使用微分方程----------------
syms x(t)
r =0.2
% 模拟20只，增长率为0.2，K=1200,55年的增长曲线
EQ = diff(x,t,1) == x * r*(1-x/1200);
x = dsolve(EQ,x(0)==20);
simplify(x)
figure(1);
hold on
fplot(x,[0,55]);
axis([0,55,0,1300]);
%--------方法二、使用递推法进行计算------------
t = 0:1:55;
y = zeros(1,56);
N_t=20;  % 初始（第0年）
for i = 1:1:56
    y(i) = N_t;
    N_t = N_t *(1+r*(1- N_t/1200));
end
plot(t,y,'--b','LineWidth',1);

得到我们典型的种群增长曲线：

但逻辑斯蒂方程的内容不止于此，由逻辑斯蒂方程的形式，引出逻辑斯蒂映射：

对于上式中求解不动点的问题，取一个Nt0 = 0.1并且使用数值的迭代解法时（N0应当在0,1之间）当k取得不同的常数值时，得到的解是不同的，甚至可能产生混沌现象。

先贴上两个网上搜索的费根鲍姆常数的数值解所收敛到的解的图形（由于这个本来就是“计算器按出来的常数，而且目前求解得到的常数的精度也不是很高”）

经过费根鲍姆的研究，每一段分叉之间的距离之比，逐渐地趋近于一个常数 $\delta$ 约为4.669.....（讲解见图一），（当然了也有一个叫做费根鲍姆第二常数 $\alpha$ 的东西，描述的是叉之间的距离之比）

下面是我对吴振奎老师文章中一些绘制图像的Matlab实现：

% 参考：https://max.book118.com/html/2019/0801/6000041052002052.shtm
% 参考：http://www.doczj.com/doc/294915043169a4517723a340.html
%  《费恩鲍姆常数4.669...》吴振奎
%   费根鲍姆常数有费根鲍姆第一常数和费根鲍姆第二常数
clear

syms x
figure(2);
hold on
fplot(x)
fplot(3*x*(1-x))
axis([0,1,0,1])
title('以两个方程为例求解不动点')
hold off
% 演示交点的求解
clear



%------设置可变参数--------
syms x

% k=3.5;  %
x0 = 0.1;    % 设置初始值
IterNumber = 1000;   %设置迭代次数为1000
k_try=[2,2.7,3,3.35,3.47,4];  %
for fig=1:1:length(k_try)
    %--------------
    k = k_try(fig);
    x_i = x0;
    st = 1:1:IterNumber;
    x_show = zeros(1,IterNumber);
    x_ = zeros(1,2*IterNumber);
    y_ = zeros(1,2*IterNumber);
    for i = 1:2:2*IterNumber-1
        % 演示不同的k值下的不动点求解
        x_(i) = x_i;
        y_(i) = x_i;
        x_(i+1) = x_i;
        y_(i+1) = k*x_i*(1-x_i);
        %----前面是求解每一个点处的x,y坐标
        x_show((i+1)/2) = x_i;
        %-----迭代------
        x_i = y_(i+1);
    end
    
    % 注意！！！！在图像非常多，使用数字会混乱而不能新建图像，所以使用命名方法
    figure('name',['showphoto',num2str(2+2*fig-1)]);
    hold on
    fplot(x);
    fplot(k*x*(1-x));
    plot(x_,y_);
    axis([0,1,0,1]);
    title(['\fontsize{13}k=',num2str(k),'的解'])
    hold off
    
    figure('name',['showphoto',num2str(2+2*fig)]); 
    hold on
    plot(st,x_show)
    title(['\fontsize{13}k=',num2str(k),'解的震荡图像'])
    hold off
end

当然了，就像网上搜索的图形一样，在k<3时，迭代法总是可以收敛到相应的解（这一部分研究来自我文章开头给的参考文章）

(K=3时，收敛速度极慢，收敛到单个的解极其困难)

但是当K=3.35时，会收敛到两个解（抛物线上的两个点）,对应的是k>3时的两条分叉

当K=3.9....4时，几乎是无序的，没有收敛到一个确定的值，有很多的周期，这样的“混沌算法”其实有时是用于随机数的求解的

所以想绘制费根鲍姆图也成了一个难题，由于有那么多的收敛解，也几乎不可能以一个函数画出来

所以下面就说说我对费根鲍姆图绘制的一些进展：

首先，以k = 3.35的情况为例，本来数值解法是有两个解的，并且最终呈现一个在两个解之间循环交替的趋势，所以我将迭代方程设出来，所收敛到的解x应当满足方程:

$f(f(x)) = x$

采用r = 3.35进行求解，有意思的是，它不是两个解，而是3个非0解，通过放大图像发现， $\frac{87}{134}\pm\frac{\surd(609)}{134}$ 才是最终收敛到的值，那么，为什么它没有收敛到47/67这个解呢？

我还发现，对于47/67,直接代入 $f(x)=x$ ,也是成立的，也就是说，47/67是原迭代方程的精确（可以说是唯一）解，

放大k=3.35解震荡图像的前半部分,发现它刚好错过了！没有到47/67处，然后逐渐放大收敛到另外的两个根

在上面的分析中，我们暂且称那个没有收敛到的根47/67为“单根理论解”，所以只要分别取不同的r值，都应该有一个单根理论解满足 $f(x)=x$ ，所以通过取不同的r，我绘制出了单根理论解的图像：

由此递推，对于每一个比双分叉点大的r，均有一对“理论双根解”（这个除掉0解和单根解之外）满足：

$f(f(x))=x$ 但是并不仅仅是比双分叉点大的r有3个非0解，这个方程对于任何一个r都有解，但是在双分叉点以前是虚数解，假如忽略它们的虚部进行绘图，大概张这个样子（第二张图，没有把单根解画出来）：

所以后来，我将双根解的前半部分替换成原始单根解，得到了理论双根解的图像：

可以看出，这个图像的前半部分，与费根鲍姆图是一样的

那么，由此类推，必定也有“理论四根解”，“理论八根解”.......,，即满足： $f(f(f(f(x))))=x$ ，这样我们就能算出费根鲍姆图的第三个分叉来，于是，只要取足够的精度和有足够的算力（就下面这个图大概我的电脑消耗一分多钟才算得出来，而且精度也不是很高），是不是就可以绘制出更加精确的绘制费根鲍姆图并有利于算出这个常数呢？

---------------------所以接下来的一切都是未解之谜--------------------------------------

希望后面会有人会继续揭开一些天空上飘着的乌云吧，一方面作者水平实在有限，也没有用过多的精力去看更多的内容。

下面是我计算理论单根解和理论双根解的代码附上：

% 费根鲍姆图的理论稳态解
% 注意这个可能会运行好几分钟时间
% 如果第二幅图不显示打掉58行的分号试试
clear
% ----这一项是可调节参数，用来调节步长
r_0 = 1:0.025:3.9;


syms f0(x)
f0(x) =  x*(1-x);


y_0 = zeros(1,length(r_0));
i=1;
for r = r_0
    f(x) = f0(x)*r;
    eq = f(x) == x;
    A = solve(eq,x); %解得结果是矩阵
    if A(1) ~=0
        y_0(i) = A(1);
    else
        y_0(i) = A(2);
    end
    i=i+1;
    % 当达到4时，可能会无限循环
end

figure('name','费根鲍姆图的理论单根解');
plot(r_0,y_0); title('\fontsize{13}费根鲍姆图的理论单根解')
axis([0,4,0,1]);


% 第二幅图:费根鲍姆图的理论
% 这一个数据有必要和上一组大小r0保持一致，因为我是借用的上一组的y_0对y做了初始化

y = y_0;

for j = 3:4  %取得它第3,4个解（即为多解（循环解））
    i=1;
    for r = r_0
        f(x) = f0(x)*r;
        eq = f(f(x)) == x;    
        %其实，f可以多重嵌套甚至无限嵌套，对应不同的混沌解
        A = solve(eq,x); %解得结果是矩阵
        %(可以看出,在所得的四个解中，第一个是0解，第二个是理论稳定的单根解，
        %(后边的两个一开始是虚根解，到分裂点后，分裂为两个部分，都是实根解)
        
        if imag(A(j)) ~= 0
           % 什么也不做,保持原有的单根解（这里删除了虚数解） 
        else
            y(i) = A(j);  % 将对应的项设为混沌解项
        end
        %---下面的几行是用于修正解的--------
        if A(j) == y_0(i) %同解的时候，显示同解时A的情况并且
            % disp(A)
            % disp(y_0(i))---->这两句取消注释后，输出结果在下面
            if j == 3
                y(i) = A(2);
            end
        end
        % 当达到4时，可能会无限循环
        i=i+1;
    end
    figure(30)  %集中画在一幅图里面
    hold on
    plot(r_0,y);
    axis([0,4,0,1]);
    title('费根鲍姆图的理论双根解');
end
axis([0,4,0,1]);

% 注意开启了55行的“disp后”显示的是解的修正
%  0
%  2/3
%  2/3
%  2/3
%  第一个是交点，其他的都是
%     0.6667
% 
%      0
%  36/61  % 此时第二个解对应的是j=3的情况，所以需要进行修正
%  41/61
%  45/61
% ..........

至于前面有一段代码里面所说的第二个资料，我确实也没细看（他说是给了费根鲍姆常数的求解方法，但是我也没VIP哈）

参考部分和原创部分都已经说明地非常清楚了

毕竟也是后来因为算法和模型建立上，我还算是个一窍不通的小白，这一次也算当做是一个算法练习，也算是一个留给后人待发现的内容吧