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

vwkv1x7d 于 7个月前发布在 Matlab

关注(0)|答案(1)|浏览(114)

尊敬的各位Maven和资深同事，我想用提供的代码创建下图。然而，它似乎需要很长的时间，没有明显的结果。

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

来源：https://stackoverflow.com/questions/77285678/generate-a-bifurcation-diagram-using-the-provided-matlab-code

1条答案

按热度按时间

b09cbbtk1#

你会得到两种类型的行为，收敛到一个平衡（有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

这就给出了一个类似于