Hopf Bifurcation diagram for 3D system

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;

[Umar]

Hi kdvO,

After analyzing your code, the issue causing the plotted lines to appear straight instead of showing the expected Hopf bifurcation behavior lies in how the detection of Hopf bifurcation is implemented within the loop. The current logic for detecting the bifurcation relies on comparing the difference between the last value and the current value of `x_result` to a threshold of 0.01. However, this approach may not accurately capture the oscillatory behavior characteristic of a Hopf bifurcation.

To fix this, you need to modify the condition for detecting the Hopf bifurcation by considering both the real and imaginary parts of the solution. Specifically, you can check for changes in the sign of the imaginary part to identify oscillatory behavior. Here's how you can update the code:

% Check for Hopf bifurcation (change in stability)

% Look for oscillatory behavior in the solution

if i > 1

    if sign(imag(Y(end, 1))) ~= sign(imag(x_result(i-1)))
        fprintf('Hopf bifurcation detected at eta = %f\n', eta);
    end
end

By comparing the signs of the imaginary parts of `Y(end, 1)` and `x_result(i-1)`, you can detect when there is a change in stability, indicative of a Hopf bifurcation. Additionally, to enhance visualization and better observe the bifurcation behavior, you might consider increasing the number of points in `eta_values` or adjusting other parameters related to your system dynamics.

Once you make these changes, rerun your code to generate an accurate Hopf bifurcation diagram that showcases the expected behavior with oscillations and dynamic changes in stability. This adjustment should address the issue of getting straight lines in your plot and provide a more informative representation of the system's behavior. Good luck!

[kdv0]

Hi I made the change you suggested and tried it with multiple parameter sets and intial conditions. Still getting straight lines in all of them.

[Umar]

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]

Hi kdv0,

Also, forgot to mention, please feel free to customize my code to resolve your problem.

Answer 1

[Francisco J. Triveno Vargas]

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.