- 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 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Athanasios Paraskevopoulos
il 15 Mag 2024
Commentato: William Rose
il 16 Mag 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 Commenti
Risposta accettata
William Rose
il 16 Mag 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 Commenti
William Rose
il 16 Mag 2024
@Sam Chak, thanks. Once again, I wish I could give a thumbs up to your explanation.
@Athanasios Paraskevopoulos, good luck with your work.
Più risposte (0)
Vedere anche
Categorie
Scopri di più su Calculus in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!