Hopf Bifurcation diagram for 3D system

47 vues (au cours des 30 derniers jours)
kdv0
kdv0 le 11 Juil 2024
Hello I want to plot the Hopf bifurcation diagram, I have this code but I'm getting nothing other than straight lines.
% Parameters
r=3;
a=0.2;
b=2;
gamma=0.3;
m=0.92;
c=2.5;
alpha=10;
p=100.001;
q=0.01;
% Define the range of parameter values for bifurcation diagram
eta_values = linspace(0, 3, 1000); % Range of 'eta' values
% Arrays to store results
x_result = zeros(size(eta_values));
y_result = zeros(size(eta_values));
z_result = zeros(size(eta_values));
% Initial conditions (for detecting the Hopf bifurcation)
x0 = 1000;
y0 = 200;
z0 = 600;
% Iterate over each 'eta' value
for i = 1:length(eta_values)
eta = eta_values(i);
% Solve the system using ODE45
tspan = [0 1000]; % Time span for integration
odefun = @(t, Y) [r*Y(1)*(1-Y(1)/Y(3))-a*(1-m)*Y(1)*Y(2)/(1+gamma*(1-m)*Y(1));
b*(1-m)*Y(1)*Y(2)/(1+gamma*(1-m)*Y(1))-c*Y(2);
(Y(3)-alpha)*(p-Y(3))/(q+Y(3));];
[~, Y] = ode45(odefun, tspan, [x0; y0; z0]);
% Check for Hopf bifurcation (change in stability)
% Look for oscillatory behavior in the solution
if i > 1
if abs(Y(end, 1) - x_result(i-1)) > 0.01
fprintf('Hopf bifurcation detected at eta = %f\n', eta);
end
end
% Record the final values (or a characteristic of the solution)
x_result(i) = Y(end, 1);
y_result(i) = Y(end, 2);
z_result(i) = Y(end, 3);
end
% Plot the bifurcation diagram
figure;
plot(eta_values, x_result, '.k', 'MarkerSize', 1);
hold on;
plot(eta_values, y_result, '.b', 'MarkerSize', 1);
plot(eta_values, z_result, '.r', 'MarkerSize', 1);
xlabel('\eta');
ylabel('Steady State Values');
legend('x', 'y', 'z');
title('Hopf Bifurcation Diagram');
grid on;
  4 commentaires
Umar
Umar le 11 Juil 2024

Hi kdv0,

I was experimenting with your code and was able to print Hopf Bifurcation diagram for 3D system. One potential issue in the code is the initialization of the x0, y0, and z0 variables for the initial conditions. These values are set too high, which might lead to numerical instability during the integration process. Please see attached plot along with snippet code.

Umar
Umar le 11 Juil 2024
Hi kdv0,
Also, forgot to mention, please feel free to customize my code to resolve your problem.

Connectez-vous pour commenter.

Réponses (1)

Francisco J. Triveno Vargas
Francisco J. Triveno Vargas le 15 Juil 2024
Modifié(e) : Francisco J. Triveno Vargas le 15 Juil 2024
Hi kdv0
Values of x_result, y_result and z_result are near of constants:
x_result = 31.4496
y_result = 98.6978
z_result = 100.0010
your plot take a fixed values. You need to change the last lines by:
figure;
plot3(Y(:,1),Y(:,2),Y(:,3))
grid on
because Y is the solution of your diferential equation, after that you obtain the figure below, also you can use surf or mesh. Regards.

Catégories

En savoir plus sur Nonlinear Dynamics dans Help Center et File Exchange

Produits


Version

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by