Basic Difference Equation Plot

3 views (last 30 days)
Sidant Dev on 19 Nov 2021
Commented: Paul on 19 Nov 2021
I am trying to graph the following difference equation:
% Default plant population
P = 1800;
f(P) = 4*P;
g(P) = P/2;
% Revised Parameters
beta = 2003.4262446653865;
rI = 1.9566407977694316;
I = zeros(1, 500);
I(1) = 50;
for t = 2:100
I(t+1) = I(t) + rI * I(t) * (1 - I(t)/f(P)) - beta*(I(t)^2/((g(P))^2)+I(t)^2);
end
plot(I)
However, I am not getting a plot that updates the values past I(1). How do I fix this? Paul on 19 Nov 2021
Look carefully at these lines
I = zeros(1, 500);
I(1) = 50;
for t = 2:100
I(t+1) = I(t) + rI * I(t) * (1 - I(t)/f(P)) - beta*(I(t)^2/((g(P))^2)+I(t)^2);
end
When t = 2, the line inside the for loop is
I(3) = I(2) + rI * I(2) + (1 - I(2)/f(P)) - beta*(I(2)^2/((g(P))^2)+I(2)^2)
But I(2) is zero, so all the terms on the RHS side are zero, so I(3) stays zero. Same thing happens with t =3, and so on.
Maybe the loop should really be
for t = 1:199
Paul on 19 Nov 2021
Let's change the code as suggested and see what happens.
% Default plant population
P = 1800;
f(P) = 4*P;
g(P) = P/2;
% Revised Parameters
beta = 2003.4262446653865;
rI = 1.9566407977694316;
I = zeros(1, 500);
I(1) = 50;
for t = 1:199
I(t+1) = I(t) + rI * I(t) * (1 - I(t)/f(P)) - beta*(I(t)^2/((g(P))^2)+I(t)^2);
end
plot(I)
xlabel('N') Recall that I was initialized to zero and the loop only updates the elements of I trhough I(200). So I(201:500) = 0, which is relfected in the plot for N > 201.
What are the first few values of I?
format short e
I(1:10).'
ans = 10×1
1.0e+00 * 5.0000e+01 -5.0084e+06 -5.0255e+16 -5.0597e+36 -5.1289e+76 -5.2702e+156 -Inf -Inf -Inf -Inf
I(1) = 50 as given in the code, and then I very rapidly becomes more and more negative until it becomes too negative to represent as a double, and so becomes -Inf. Once I(7) becomes -inf, every subsequent iterate will also be -Inf. The plot shows the first six values of I rapidly decreasing and then white space wherever I == -inf, which is exactly how plot() works.
If this result is not what's expected, then either data going into the equation or the equation itself or both need to be modified.
Also, note that these lines
f(P) = 4*P;
g(P) = P/2;
cause f and g to be allocated as 1800 element arrays, with only the last element of each being non-zero
size(f)
ans = 1×2
1 1800
size(g)
ans = 1×2
1 1800
[find(g~=0) find(f~=0)]
ans = 1×2
1800 1800
Obviously the code still runs, but there's really no reason to implement it this way. Could just as easily defined scalar variables fP and gP.