Basic Difference Equation Plot
21 次查看(过去 30 天)
显示 更早的评论
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?
0 个评论
回答(1 个)
Paul
2021-11-19
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
2 个评论
Paul
2021-11-19
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).'
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)
size(g)
[find(g~=0) find(f~=0)]
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.
另请参阅
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!