Population Growth Model Analysis

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.

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

where the term is given here to be and it represents a term for population reduction. The constants are positive. Normalize the equation by setting and .Find the steady solutions of the problem (real and positive) and then implement a program in MATLAB to numerically solve the dimensionless form of the equation. Make the graphical representation of the numerical solution for , with initial conditions , and respectively. What do you observe?

So we have

• We are going to find the function :

To find , we integrate :

where C is the integration constant, which can be determined if an initial condition is given, such as .

• Now we are going to normalize the equation:

Let and . Then:

We substitute into the original equation:

• Then we are going to calculate the Stationary Points:

Set :

\]

For :

% 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:');

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

N = 0.000000: Eigenvalue = -0.600000, Stable
N = 1.621444: Eigenvalue = 0.217420, Unstable
N = 32.378556: Eigenvalue = -0.340979, Stable

 % 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*');

Warning: Imaginary parts of complex X and/or Y arguments ignored.

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;

Answer 1

[William Rose]

@Athanasios Paraskevopoulos,

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:

• 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?

[Sam Chak]

Hi @William Rose, you're right. That duplicate '0.4' causes the equilibrium to shift. Removing that will fix the issue.

Ka   = 35;
ar_b = 0.4;

% Set up the ODE function for phase plane analysis
f = @(N) ar_b * N .* (1 - N / Ka) - N ./ (1 + N);

% Mark zero crossings (equilibrium points)
% zero_crossings = N_values(abs(dN_dtau) < 1e-3)
zero_crossings = equilibria

zero_crossings = 3x1

0 1.6214 32.3786

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

for i = 1:length(zero_crossings)
    plot(zero_crossings(i), 0, 'ro');
end
xline(2, '--', 'initial value, N = 2')
hold off;

[Athanasios Paraskevopoulos]

Thank you both of you this is my revised code

% 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(@(tau, N) f(tau, N), 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 \tau');
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_eq = @(N) ar_b * N * (1 - N / Ka) - N / (1 + N);
J = diff(f_eq(N), N);
eigenvalues = double(subs(J, N, equilibria));

disp('Equilibrium Points and their Stability:');

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

N = 0.000000: Eigenvalue = -0.600000, Stable
N = 1.621444: Eigenvalue = 0.217420, Unstable
N = 32.378556: Eigenvalue = -0.340979, Stable

% 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*');

Warning: Imaginary parts of complex X and/or Y arguments ignored.

title('Bifurcation Diagram');
xlabel('ar/b');
ylabel('Equilibrium Points');

% Set up the ODE function for phase plane analysis
f_phase = @(N) ar_b * 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_phase(N_values);

% Plot phase plane
figure;
plot(N_values, dN_dtau);
title('Phase Plane Analysis');
xlabel('Population Size N');
ylabel('dN/d\tau');
grid on;
hold on;

% Mark zero crossings (equilibrium points)
zero_crossings = equilibria; % Use calculated equilibria
for i = 1:length(zero_crossings)
    plot(zero_crossings(i), 0, 'ro');
end

% Mark initial value N = 2
xline(2, '--', 'initial value, N = 2');
hold off;

[William Rose]

@Sam Chak, thanks. Once again, I wish I could give a thumbs up to your explanation.

@Athanasios Paraskevopoulos, good luck with your work.