Basic Difference Equation Plot

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?

Answer 1

[Paul]

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  

[Sidant Dev]

This doesn't seem to solve the issue of the plot making little/no sense. It appears as though Matlab isn't updating the values at all.

[Paul]

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.