使用提供的 Matlab 代码生成分岔图-解网

问：

尊敬的各位专家、各位资深同仁：我想使用提供的代码创建下图。然而，它似乎花费了很长时间，没有明显的结果。 enter image description here

这是我的Matlab代码，它没有提供任何结果。

clc;close all;clear;
results = [];

tmax = 300;

r1 = 0.36; a = 0.1; k1 = 0.36; k2 = 0.48; r2 = 1; deltav = 0.2; deltaz = 0.036;

x0 = 0.5; y0 = 0; v0 = 0.3; z0 = 0.1; 

for b = 0:0.1:30
    ode_system = @(t, y) [
        r1 * y(1) * (1 - y(1)) - a * y(1) * y(3) - k1 * y(1) * y(4);
        a * y(1) * y(3) - k2 * y(2) * y(4) - y(2);
        b * y(2) - a * y(1) * y(3) - deltav * y(3);
        r2 * y(2) * y(4) - deltaz * y(4)
    ];

    [t, sol] = ode45(ode_system, [0, tmax], [x0, y0, v0, z0]);

    results = [results; b * ones(length(t), 1), sol(:, 1)];

    x0 = sol(end, 1);
    y0 = sol(end, 2);
    v0 = sol(end, 3);
    z0 = sol(end, 4);
end

figure;
plot(results(:, 1), results(:, 2), '.');
xlabel('b');
ylabel('x');
title('Bifurcation Diagram');
grid on;

MATLAB MATLAB-图数值方法常微程图

这样的系统通常要么爆炸，要么收敛到一个平衡点。在后一种情况下，解几乎是一条直线，但雅可比的特征值限制了显式方法的步长。更改隐式方法以获得更大的步长。请更准确地描述预期的绘图，您似乎使用（和存储）了太多数据。选择默认公差以生成求解函数的可用图，它们对于像素显示足够精确，但仅此而已。您可能需要关闭“优化”计算。ode45

0赞 Lutz Lehmann 10/13/2023

尝试加速：设置，只存储一对。要了解会发生什么，请在计算结束后直接打印一条消息，由于有一万个图形对象，绘制可能需要很长时间。tmax = 3b, sol(end,1)

0赞 ZSN 10/13/2023

感谢您的回复。我的代码适用于 tmax=32，但不适用于 tmax>32。

0赞 Lutz Lehmann 10/13/2023

你应该探索你实际得到的数据。当您给出时间跨度时，输出由所有内部步骤和每个内部段的 3 个附加点组成。使用显式方法时，由于 ODE 的刚度，步长被限制在某个小值，因此您可以在每个积分的输出中获得许多许多点。你把它们都画出来了。图形对象预计最多绘制几万个点，它没有优化以容纳更多点，因为通常更多的点不会产生明显的差异。因此，修剪数据，将时间跨度替换为采样时间列表,...

答：

0赞 Lutz Lehmann 10/14/2023 #1

你会得到两种类型的行为，趋向平衡（有）和振荡，可能会限制周期，但不排除奇怪的吸引子。也许对于参数的不同段，甚至有不同的振荡模式。在第二种情况下，您需要在图中绘制振荡的跨度，即在中绘制极值。这些极值位于的导数为零的位置。y(1)=1by(1)y(1)

Matlab ODE求解器有一个很好的机制，可以在称为“事件”的解上定位和注册零函数。因此，按照规范构建事件函数

    function [val, term, dir] = equi1(t,y)
      val = r1*(1-y(1))-a*y(3)-k1*y(4);
      term = false;
      dir = 0;
    end

将其插入到参数列表中

    opts = odeset('RelTol',1e-10,'AbsTol',1e-8);
    opts = odeset(opts, 'Events', @(t,y) equi1(t,y));

并使用参数列表调用求解器

    [t, sol, t_e, sol_e, ind_e] = ode45(ode_system, [0:0.1:tmax], [x0, y0, v0, z0], opts);
    if length(t_e)==0 % no oscillation, monotonic to the equilibrium
      results = [results; b * ones(4, 1), sol(end-3:end, 1)];
    else % there were maxima and minima
      results = [results; b * ones(length(t_e), 1), sol_e(:, 1)];
    end%if

然后给出一个像这样的情节

我是用八度音阶完成的，那里产生了一个不那么摇摆不定的图表，一路上的步长有问题。ode23sode15s

上一个：Python 图表更改背景大小

下一个：Excel VBA：图表的动态 x 轴取决于输入

使用提供的 Matlab 代码生成分岔图

Generate a bifurcation diagram using the provided Matlab code

评论