Solution Plot of the Constant Harvesting Function with Directional Field
3 views (last 30 days)
Show older comments
I'm trying to plot solutions along with the directional field as in the image given below, of the constant harvesting model given by
dP/dt = r*P(1 - (P/M)) - h , P(0) =p0 , where r = 0.8, M = 8000 and h = 0 in this particular case.
I attempted using the following code, which does not result in something as in the image given below. I would appreciate your help with making the graph looks like as in the image above.
% Constant Harvesting Model
% Parameters
r = 0.8; % Intrinsic growth rate
K = 8000; % Carrying capacity
% Define the system of differential equations
dN_dt = @(N) r * N .* (1 - N/K);
% Define the range of values for N and t
N_range = 0:2000:14000;
t_range = 0:5:20;
% Create a grid of N and t values
[N, t] = meshgrid(N_range, t_range);
% Compute the rate of change of N at each point on the grid
dN = dN_dt(N);
% Create a figure
figure;
hold on;
% Plot the directional field
quiver(t, N, ones(size(dN)), dN, 'k');
% Set plot properties
xlabel('t');
ylabel('N');
title('Constant Harvesting Model');
legend('Directional Field');
grid on;
% Invert the t-axis
set(gca, 'XDir', 'reverse');
% Display the figure
hold off;
0 Comments
Answers (1)
Sam Chak
on 4 Jun 2023
Hi @Gith
I'm unsure, but I think that it was the quiver scaling issue when the carrying capacity K is too large (on the Y-axis), compared to the time scale (on the X-axis). One workaround is to rescale on the Y-axis, and in your case, it is related to the value K.
r = 0.8;
K = 8000/1000;
[T, N] = meshgrid(0:20/14:20, 0:K/14:K);
S = r*N.*(1 - N/K);
L = sqrt(1 + S.^2);
U = 1./L;
V = S./L;
quiver(T, N, U, V, 0.3)
axis square
hold on
% differential equation
f = @(T, N) r*N.*(1 - N/K);
tspan = 0:0.01:20;
init = 10/1000;
[T, N] = ode45(f, tspan, init);
plot(T, N, 'r', 'linewidth', 1.5)
hold off
grid on
xlabel({'$t$'}, 'Interpreter', 'latex')
ylabel({'$N(t), \quad (\times 1000)$'}, 'Interpreter', 'latex')
0 Comments
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!