I want to plot phase diagram.I have written a code of nullcline and want plot phase diagram

clc;clear;close all;

r=0.1;

k=50;

a=0.01;

e=0.5;

m=0.05;

s=0.1;

w=0.1;

b=0.001;

q=0.03;

F=1.613;

M=50;

X=0:1:50;J=[0 50 0 20];

Pi_0=(w/(w+b*M))*(F/s);Pi_1=((w+b*M)/w)*(F/s);

D1=Pi_0+0*X;D2=Pi_1+0*X;

plot(D1,X,'k--')

axis(J)

hold on

plot(D2,X,'k--')

hold on

X1=(Pi_0/99);X2=(Pi_1-Pi_0)/99;X3=(50-Pi_1)/99;Y2=20/99;

P1=0:X1:Pi_0;H1=0:Y2:20;P2=Pi_0:X2:Pi_1;H2=0:Y2:20;P3=Pi_1:X3:50;H3=0:Y2:20;

u_1=0;u_2=((s.*P2)/(s.*P2+F)+(w*(s.*P2-F))/(b*M*(s.*P2+F)));u_3=1;

f = @(P1,H1) (r.*P1).*(1-P1/k)-(a.*P1.*H1)./(1+q*u_1*M);

fimplicit(f,[0 Pi_0 0 20],'g')

hold on

f1 = @(P1,H1) (e*a.*P1.*H1)./(1+q*u_1*M)- m.*H1;

fimplicit(f1,[0 Pi_0 0 20],'g')

hold on

f = @(P2,H2) (r.*P2).*(1-P2/k)-(a.*P2.*H2)./(1+q*((s.*P2)/(s.*P2+F)+(w*(s.*P2-F))/(b*M*(s.*P2+F)))*M);

fimplicit(f,[Pi_0 Pi_1 0 20],'r')

hold on

f1 = @(P2,H2) (e*a.*P2.*H2)./(1+q*((s.*P2)/(s.*P2+F)+(w*(s.*P2-F))/(b*M*(s.*P2+F)))*M)- m.*H2;

fimplicit(f1,[Pi_0 Pi_1 0 20],'r')

hold on

f = @(P3,H3) (r.*P3).*(1-P3/k)-(a.*P3.*H3)./(1+q*u_3*M);

fimplicit(f,[Pi_1 50 0 20],'m')

hold on

f1 = @(P3,H3) (e*a.*P3.*H3)./(1+q*u_3*M)- m.*H3;

fimplicit(f1,[Pi_1 50 0 20],'m')

Answer 1

[Abhinaya Kennedy]

This version of your code includes plotting the nullclines and the phase diagram. It uses the "quiver" function (https://www.mathworks.com/help/matlab/ref/quiver.html) to plot the phase portrait.

clc; clear; close all;

% Parameters
r = 0.1;
k = 50;
a = 0.01;
e = 0.5;
m = 0.05;
s = 0.1;
w = 0.1;
b = 0.001;
q = 0.03;
F = 1.613;
M = 50;

% Phase space
X = 0:1:50;
J = [0 50 0 20];

% Nullclines
Pi_0 = (w / (w + b * M)) * (F / s);
Pi_1 = ((w + b * M) / w) * (F / s);
D1 = Pi_0 + 0 * X;
D2 = Pi_1 + 0 * X;

figure;
plot(D1, X, 'k--')
axis(J)
hold on
plot(D2, X, 'k--')
hold on

% Define the nullclines
f1 = @(P, H, u) (r .* P) .* (1 - P / k) - (a .* P .* H) ./ (1 + q * u * M);
f2 = @(P, H, u) (e * a .* P .* H) ./ (1 + q * u * M) - m .* H;

% Plot nullclines
fimplicit(@(P, H) f1(P, H, 0), [0 Pi_0 0 20], 'g')
hold on
fimplicit(@(P, H) f2(P, H, 0), [0 Pi_0 0 20], 'g')
hold on
fimplicit(@(P, H) f1(P, H, (s .* P) ./ (s .* P + F) + (w * (s .* P - F)) ./ (b * M * (s .* P + F))), [Pi_0 Pi_1 0 20], 'r')
hold on
fimplicit(@(P, H) f2(P, H, (s .* P) ./ (s .* P + F) + (w * (s .* P - F)) ./ (b * M * (s .* P + F))), [Pi_0 Pi_1 0 20], 'r')
hold on
fimplicit(@(P, H) f1(P, H, 1), [Pi_1 50 0 20], 'm')
hold on
fimplicit(@(P, H) f2(P, H, 1), [Pi_1 50 0 20], 'm')
hold on

% Create a grid for the phase portrait
[P, H] = meshgrid(0:2:50, 0:2:20);

% Define the system of ODEs
dP = @(P, H) (r .* P) .* (1 - P / k) - (a .* P .* H) ./ (1 + q * ((s .* P) ./ (s .* P + F) + (w * (s .* P - F)) ./ (b * M * (s .* P + F))) * M);
dH = @(P, H) (e * a .* P .* H) ./ (1 + q * ((s .* P) ./ (s .* P + F) + (w * (s .* P - F)) ./ (b * M * (s .* P + F))) * M) - m .* H;

% Compute the derivatives
dP_dt = dP(P, H);
dH_dt = dH(P, H);

% Normalize the arrows for better visualization
norm_factor = sqrt(dP_dt.^2 + dH_dt.^2);
dP_dt = dP_dt ./ norm_factor;
dH_dt = dH_dt ./ norm_factor;

% Plot the phase portrait
quiver(P, H, dP_dt, dH_dt, 'b')
xlabel('P')
ylabel('H')
title('Phase Diagram with Nullclines')
legend('Nullcline 1', 'Nullcline 2', 'Nullcline 3', 'Nullcline 4', 'Nullcline 5', 'Nullcline 6', 'Phase Portrait')
hold off