- The first plot, with r,g,b,k traces, shows that N goes asymptotically to about N=32.5 when N0=2, 25, and 40. The phase plane plot shows what I interpret as a stable equilibrium at N~=27 (the value of N where dN/dtau crosses zero, with negative slope). Why is the equilibrium value for N not the same in the first plot and the phase plane plot?
- The phase plane plot shows dN/dtau<0 when N=2. By tracing the expected path on the phase plane plot, starting at N=2 , I expect N to decrease asymptotically to zero. But the first plot shows dN/dt>0 when N=N0=2. Why the mismatch?
Population Growth Model Analysis
9 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Athanasios Paraskevopoulos
am 15 Mai 2024
Kommentiert: William Rose
am 16 Mai 2024
Below I tried to solve Population Growth problem, I provide you the mathematical solution. Also, I create MATLAB code, however i feel like that i repeat this code to each part of this code as you can see. How could I improve this code?
Assume a continuous growth model where the function that gives the rate of change R is defined as
. Provide the function n, which represents the population size.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696096/image.png)
If the growth of a population follows the logistic law but there is also a term for population reduction (e.g., from a predator while the members of the population are the prey), the equation that expresses the population size is of the form
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696101/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696106/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696111/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696116/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696121/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696126/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696131/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696136/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696141/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696146/image.png)
So we have
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696151/image.png)
- We are going to find the function
:
To find
, we integrate
:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696156/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696166/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696171/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696176/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696181/image.png)
where C is the integration constant, which can be determined if an initial condition is given, such as
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696186/image.png)
- Now we are going to normalize the equation:
Let
and
. Then:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696191/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696196/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696201/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696206/image.png)
We substitute into the original equation:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696211/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696216/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696221/image.png)
- Then we are going to calculate the Stationary Points:
Set
:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696226/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696231/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696236/image.png)
For
:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696241/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696246/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1696251/image.png)
% Define the normalized logistic equation with reduction term
% K/a = 35, ar/b = 0.4
Ka = 35;
ar_b = 0.4;
% Define the function for the normalized logistic equation
f = @(tau, N) ar_b * N * (1 - N/Ka) - N / (1 + N);
% Time span for the simulation
tspan = [0 50]; % Adjust if needed to observe long-term behavior
% Initial conditions
N0 = [40, 25, 2, 0.05];
% Prepare the figure for plotting
figure(1);
hold on;
% Colors for different trajectories
colors = ['r', 'g', 'b', 'k']; % Red, Green, Blue, Black
% Solve the equation for each initial condition and plot
for i = 1:length(N0)
[T, N] = ode45(f, tspan, N0(i));
plot(T, N, 'Color', colors(i), 'LineWidth', 2);
legendInfo{i} = ['N_0 = ' num2str(N0(i))]; % Create legend entry
end
% Add graph details
title('Numerical solutions of the normalized logistic equation');
xlabel('Normalized time τ');
ylabel('Normalized population size N');
legend(legendInfo, 'Location', 'northeastoutside');
grid on;
hold off;
% Equilibrium analysis
syms N
eqn = ar_b * N * (1 - N / Ka) - N / (1 + N) == 0;
equilibria = double(solve(eqn, N));
% Stability analysis
f = @(N) ar_b * N * (1 - N / Ka) - N / (1 + N);
J = diff(f(N), N);
eigenvalues = double(subs(J, N, equilibria));
disp('Equilibrium Points and their Stability:');
for i = 1:length(equilibria)
if real(eigenvalues(i)) < 0
stability = 'Stable';
else
stability = 'Unstable';
end
fprintf('N = %f: Eigenvalue = %f, %s\n', equilibria(i), eigenvalues(i), stability);
end
% Bifurcation analysis setup
ar_b_values = linspace(0.1, 1, 100);
bifurcation_results = zeros(length(ar_b_values), 10); % Adjust the columns based on expected number of roots
% Compute equilibrium points for varying ar/b
for i = 1:length(ar_b_values)
eqn = ar_b_values(i) * N * (1 - N / Ka) - N / (1 + N) == 0;
solutions = double(solve(eqn, N));
bifurcation_results(i, 1:length(solutions)) = solutions';
end
% Plot bifurcation diagram
figure;
plot(ar_b_values, bifurcation_results, 'b*');
title('Bifurcation Diagram');
xlabel('ar/b');
ylabel('Equilibrium Points');
% Set up the ODE function for phase plane analysis
f = @(N) ar_b * 0.4 * N .* (1 - N / Ka) - N ./ (1 + N);
% Create a grid of N values and compute their time derivatives
N_values = linspace(0, 50, 400);
dN_dtau = f(N_values);
% Plot phase plane
figure;
plot(N_values, dN_dtau);
title('Phase Plane Analysis');
xlabel('Population Size N');
ylabel('dN/dτ');
grid on;
hold on;
% Mark zero crossings (equilibrium points)
zero_crossings = N_values(abs(dN_dtau) < 1e-3);
for i = 1:length(zero_crossings)
plot(zero_crossings(i), 0, 'ro');
end
hold off;
0 Kommentare
Akzeptierte Antwort
William Rose
am 16 Mai 2024
I see that you have
f = @(tau, N) ar_b * N * (1 - N/Ka) - N / (1 + N);
and
syms N
eqn = ar_b * N * (1 - N / Ka) - N / (1 + N) == 0;
...
f = @(N) ar_b * N * (1 - N / Ka) - N / (1 + N);
and
eqn = ar_b_values(i) * N * (1 - N / Ka) - N / (1 + N) == 0; % inside a for loop
and
f = @(N) ar_b * 0.4 * N .* (1 - N / Ka) - N ./ (1 + N);
I wouldn;t worry about the redundancy. The code is clear, which makes it easier to understand, for others, and for you, if you sxccome back to it in a couple of years.
I don't understand why there is "0.4" in the last equation.
I have two other questions - which reveal my lack of understanding of the model, but I'm curious:
3 Kommentare
William Rose
am 16 Mai 2024
@Sam Chak, thanks. Once again, I wish I could give a thumbs up to your explanation.
@Athanasios Paraskevopoulos, good luck with your work.
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Calculus finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!