尊敬的各位专家、各位老同事:
我想使用提供的代码创建下图。然而,这似乎花费了很长的时间,但没有明显的结果。
这是我的 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;
你会得到两种类型的行为,向平衡收敛(具有
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
一路上步长存在问题。