How to plot stability regions in matlab?
显示 更早的评论
Hi,
I would like to plot the stability region for the following system of ODE
;\\
x.=ax(1-x)-bxy,\\
y.=-dy+cxy,\\
where a b d c are positive parameters. Can I plot the division plane (a-d)?
Thank you very much in advance.
Sanaa Moussa
采纳的回答
更多回答(1 个)
clear; clc;
% ===== Parameters =====
A = 0.1;
c = 0.01;
theta = 0.02;
N = 400; % (use N=2000 for full resolution, but slower)
% ===== Fixed point (nontrivial equilibrium) =====
q1_star = (2 - theta - A) * (1 - c) / (2*(2 - theta) - A^2 * (1 - theta));
q2_star = (1 - c) * (2 - A*(1 - theta)) / (2*(2 - theta) - A^2 * (1 - theta));
% ===== Grid =====
alpha1_vals = linspace(0, 4, N);
alpha2_vals = linspace(0, 4, N);
[Alpha1, Alpha2] = meshgrid(alpha1_vals, alpha2_vals);
% ===== Storage =====
Fold = zeros(size(Alpha1));
Flip = zeros(size(Alpha1));
Neimark = zeros(size(Alpha1));
Stable = zeros(size(Alpha1));
% ===== Compute conditions =====
for i = 1:N
for j = 1:N
a1 = Alpha1(i,j);
a2 = Alpha2(i,j);
% Jacobian entries
a11 = 1 + a1*(1 - 4*q1_star - A*q2_star - c);
a12 = -a1*A*q1_star;
a21 = -a2*A*(1 - theta)*q2_star;
a22 = 1 + a2*(1 - A*(1-theta)*q1_star - 2*(2-theta)*q2_star - c);
% Trace and determinant
trJ = a11 + a22;
detJ = a11*a22 - a12*a21;
% Jury conditions
Fold(i,j) = 1 - trJ + detJ;
Flip(i,j) = 1 + trJ + detJ;
Neimark(i,j) = 1 - detJ;
% Stability check
if (Fold(i,j) > 0) && (Flip(i,j) > 0) && (Neimark(i,j) > 0)
Stable(i,j) = 1; % Stable
else
Stable(i,j) = 0; % Unstable
end
end
end
% ===== Plot =====
figure; hold on;
% --- Shade stability regions ---
% Use pcolor to show shaded background (stable vs unstable)
pcolor(Alpha1, Alpha2, Stable);
shading flat;
colormap([1 0.8 0.8; 0.7 1 0.7]); % Red-ish for unstable, green-ish for stable
clim([0 1]);
%caxis([0 1]);
% --- Overlay the bifurcation curves ---
contour(Alpha1, Alpha2, Fold, [0 0], 'r', 'LineWidth', 2);
contour(Alpha1, Alpha2, Flip, [0 0], 'b', 'LineWidth', 2);
contour(Alpha1, Alpha2, Neimark, [0 0], 'g', 'LineWidth', 2);
% --- Labels and legend ---
xlabel('\alpha_1'); ylabel('\alpha_2');
title('Stable (green) and Unstable (red) Regions with Bifurcation Curves');
legend({'Region shading','Fold curve','Flip curve','Neimark curve'}, ...
'Location','bestoutside');
grid on;
axis tight;
hold off;
类别
在 帮助中心 和 File Exchange 中查找有关 Downloads 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
