Creating phase portraits using quiver

I am trying to use quiver() to build a phase portrait of a funtion and I am having trouble. I know our points of equilibrium and I simply want to illistrate it with the phase portrait to demonstrate stability. Luckily this is a very classic ode model so it shouldnt be difficult I have just never dealt with creating a phase portrait in matlab before. Any help is very appreciated.

global  a b r K;
a = 2;
b = 2;
r= 0.48;
K = 17;
syms    P;
steady_states = vpa(solve( r*(1-P/K)-a*P/(b^2+P^2), P))

% Pv = -1:0.1:12;
% y =  r*Pv.*(1-Pv/K)-a*Pv.^2./(b^2+Pv.^2);

P = 0:0.1:20;
y1 = r.*(1-P/K);
ya = a*P./(b^2+P.^2);
plot(P,y1,'k -',P,ya,'b--','linewidth',3);

quiver(y1,ya) %failed

[x y]=meshgrid(...)
   quiver(x,y,y1,ya)%failed 

Answer 1

[Jaswanth]

Hi,

To create a phase portrait in MATLAB using “quiver” function, you need to first set up a grid of points over which you need to calculate the vector field.

In your case, you're working with a single-variable ODE, so you will be plotting a 1D phase portrait. “quiver” function is generally used for 2D vector fields and to set up a 2D grid that represents your system, you can start by:

• Defining a mesh grid over the area of interest.
• Calculate the derivative at each point in the grid.
• Use “quiver” function to plot the vector field.

Please refer to the following modified code using “quiver” to build a phase portrait of a function you have mentioned:

global a b r K;
a = 2;
b = 2;
r = 0.48;
K = 17;

% Define the function for the derivative
dPdt = @(P) r * P .* (1 - P / K) - a * P.^2 ./ (b^2 + P.^2);

% Create a grid of points for P
P = linspace(0, 20, 20);

% Calculate the derivative at each point
dP = dPdt(P);

% Create a 2D grid for plotting
[P_grid, dP_grid] = meshgrid(P, dP);

% Use quiver to plot the direction field
figure;
quiver(P_grid, dP_grid, ones(size(P_grid)), dP_grid, 'AutoScale', 'on');
xlabel('P');
ylabel('dP/dt');
title('Phase Portrait');
grid on;

% Plot the nullclines
hold on;
y1 = r .* P .* (1 - P / K);
ya = a * P ./ (b^2 + P.^2);
plot(P, y1, 'k-', 'LineWidth', 2);
plot(P, ya, 'b--', 'LineWidth', 2);
legend('Vector Field', 'Nullcline 1', 'Nullcline 2');

I hope the solution provided above is helpful.