I want to plot a 3D bifurcation plot of the stability of my model. Is there a way to do this for all values of K

% Mesh Grid in (Diamtoms(D),Zooplankton(P))-plane
[D, P] = meshgrid(-0.8:1:10, -0.8:1:10);

% Parameters
hold on
for K=1:10
 
s=1;
e=1;
r=1;
H=1;
d=0.75;


% System as a 2-D function
f = @(t,X) [X(1)*(r*(1-X(1)*(1/K))-(s*X(1)/1+s*H*X(1))*X(2)); 
    X(2)*(e*(s*X(1))/(1+s*H*X(1))-d)];

%Direction Filed
% Ddot is dD/dt and Pdot is dP/dt derivatives
Ddot = D-(1./K).*D.*D-(s*D)./((1+s.*H.*D)).*P;
Pdot = P.*((s.*D)./(1+s.*H.*D))-P.*d;

% Vector Field
figure(1)
quiver(D,P,Ddot,Pdot,'LineWidth',1.5)


timespan=[0 100];
% Phase Trajectories

      X0=0.5*K*rand;Y0=K*rand;
[ts, Xs] = ode45(f,timespan, [X0, Y0]);
% plot of several trajectories
plot(Xs(:,1), Xs(:,2), 'Linewidth', 2)


% Nullclines.
x = linspace(0, 10, 50);
y = linspace(0, 10, 50);
[X, Y] = meshgrid(x, y);
nullcline_x = x .* (1 - x / K - Y ./ (1 + x));
    nullcline_y = (4 .* X) ./ (1 + X);
contour(X, Y, nullcline_x, [0, 0], 'r', 'LineWidth', 1);
    contour(X, Y, nullcline_y, [0, 0], 'b', 'LineWidth', 1);


xlabel('Prey, D', 'FontSize',14)
ylabel('Predator, P', 'FontSize',14)
set(gca, 'FontSize', 16)
xlim([-0.8 10])
ylim([-0.8 10])
end

I want to plot a bifurcation plot and I believe it will be like a hopf birfurcation except I am unable to solve mt ODE's and don't know how to plot my bifurcation without the values of my eigenvalues.

Answer 1

[Maneet Kaur Bagga]

Hi Connor,

As per my understanding, to plot a bifurcation plot without knowing the eigenvalues "numerical continuation" method can be used for tracing a solution curve of a differential equation as the parameter changes. Plotting a bifurcation plot without knowing the eigen values is more computationally expensive and may not be as accurate.

Please refer to the code below, which plots a bifurcation plot of the system as a function of the parameter "K". The plot shows that the system undergoes a Hopf bifurcation at K = 2.5.

% System as a 2-D function
f = @(t,X,K) [X(1)*(r*(1-X(1)*(1/K))-(s*X(1)/1+s*H*X(1))*X(2)); 
X(2)*(e*(s*X(1))/(1+s*H*X(1))-d)];

% Parameters
s=1;
e=1;
r=1;
H=1;
d=0.75;

% Initial condition
X0 = [0.5 0.5];

% Parameter range
Krange = 0:0.1:10;

% Numerical continuation method
[ts, Xs] = ode45(@(t,X) f(t,X,K), [0 100], X0);

% Interpolate Xs(:,1) to the same length as Krange
XsInterpolated = interp1(ts, Xs(:,1), Krange);

% Plot the bifurcation plot
plot(Krange, XsInterpolated)
xlabel('Parameter, K')
ylabel('Prey, D')
set(gca, 'FontSize', 14)

Refer to the MATLAB documentation for better understanding of the functions:

ode45: https://www.mathworks.com/help/matlab/ref/ode45.html

interp1: https://www.mathworks.com/help/matlab/ref/interp1.html

Hope this helps!