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

Question

尊敬的各位专家、各位老同事： 我想使用提供的代码创建下图。然而，这似乎花费了很长的时间，但没有明显的结果。

这是我的 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;

Answer 1

你会得到两种类型的行为，向平衡收敛（具有

y(1)=1

）和振荡，可能是极限循环，但不排除奇怪的吸引子。也许参数

b

的不同部分甚至有不同的振荡模式。在第二种情况下，您希望在图中显示振荡的跨度，即在

y(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

然后给出如下图

我在八度音程中进行了此操作，

ode23s

生成了一个不太摆动的图表，

ode15s

一路上步长存在问题。