Lotka Volterra with an additional parameter

I have attempted to model the equations provided below. the code keeps giving me an error when I try to plot the figure using the parameter k instead of t. I am trying to model the effect of varying k on the population dynamics for both prey and predator.
note: if I use plot(k, P) it gives me a complicated and entirely different graph that is unreadable.
Script Window:
tspan = [0, 400];
P0 = [200, 50];
k = [-1, 1];
[t,P] = ode45(@(t,P)lotkavolterra(t,P),tspan,P0);
plot(k,P(:,1),k,P(:,2)) %This section gives me the error. It works, however, if I use t instead of k.
function [dPdt] = lotkawithk(t,P)
alpha = 1;
beta = 0.05;
gamma = 0.05;
delta = 0.0025;
x = P(1);
y = P(2);
dPdt = zeroes(2,1);
dPdt(1) = alpha*x - (beta*x*y)/((log(exp(1)+t))^k);
dPdt(2) = (delta*x*y)/((log(exp(1)+t))^k) - (gamma*y);
end

 Réponse acceptée

Star Strider
Star Strider le 25 Sep 2020
I cannot imagine how you got that code to work at all, considering that the function you call in ode45 is not the function you posted!
Try this:
tspan = linspace(0, 400, 100);;
P0 = [200, 50];
k = [-1, 1];
for k1 = 1:numel(k)
[t,P] = ode45(@(t,P)lotkawithk(t,P,k(k1)),tspan,P0);
Pm{k1} = P;
end
figure
hold on
for k1 = 1:numel(k)
plot(t,Pm{k1}(:,1), t,Pm{k1}(:,2), 'DisplayName',sprintf('k = %d',k(k1))) %This section gives me the error. It works, however, if I use t instead of k.
end
grid
legend
function [dPdt] = lotkawithk(t,P,k)
alpha = 1;
beta = 0.05;
gamma = 0.05;
delta = 0.0025;
x = P(1);
y = P(2);
dPdt = zeros(2,1);
dPdt(1) = alpha*x - (beta*x*y)/((log(exp(1)+t))^k);
dPdt(2) = (delta*x*y)/((log(exp(1)+t))^k) - (gamma*y);
end
Note that ‘k’ is passed as a parameter.
Make appropriate changes to get different results.

6 commentaires

Khushi Patel
Khushi Patel le 25 Sep 2020
Ohh okay. Thank you! I'm just a little confused on how to edit the legend so that it indicates which graph is predator and which is prey at k =1 and at k = -1.
My pleasure!
I’m not certain which is which, so make the appropriate changes here:
figure
hold on
for k1 = 1:numel(k)
plot(t,Pm{k1}(:,1), 'DisplayName',sprintf('Predator (k = %d)',k(k1))) %This section gives me the error. It works, however, if I use t instead of k.
plot(t,Pm{k1}(:,2), 'DisplayName',sprintf('Prey (k = %d)',k(k1)))
end
hold off
grid
legend('Location','NW')
That should do what you want. (It also positions the legend to cover as little of the plotted data as possible.)
Khushi Patel
Khushi Patel le 25 Sep 2020
Yep that works! Thank you so much.
Just one last thing, I was wondering why is it not possible to plot P against changing values of k, instead of t?
My pleasure!
It is possible. However plotting a (100x1) vector dependent variables against a scalar independent variable produces uninterpretable results, so it’s not worth the effort:
figure
hold on
for k1 = 1:numel(k)
plot(k(k1),Pm{k1}(:,1), '.', 'LineWidth',1.5, 'DisplayName',sprintf('Predator (k = %d)',k(k1))) %This section gives me the error. It works, however, if I use t instead of k.
plot(k(k1),Pm{k1}(:,2), '.', 'LineWidth',1.5, 'DisplayName',sprintf('Prey (k = %d)',k(k1)))
end
hold off
grid
xlim([-1.1 1.1])
legend('Location','N')
Note that this uses a marker to plot the points, because they are plotted as individual points, not lines. If you didn’t plot markers, that’s likely the reason that nothing showed up when you attempted to plot it.
Still another option is to plot them in 3D:
figure
hold on
for k1 = 1:numel(k)
plot3(k(k1)*ones(size(t)),t,Pm{k1}(:,1), 'LineWidth',1.5, 'DisplayName',sprintf('Predator (k = %d)',k(k1))) %This section gives me the error. It works, however, if I use t instead of k.
plot3(k(k1)*ones(size(t)),t,Pm{k1}(:,2), 'LineWidth',1.5, 'DisplayName',sprintf('Prey (k = %d)',k(k1)))
end
hold off
grid
xlim([-1.1 1.1])
legend('Location','NW')
view(30, 30)
Experiment with these to see what works best for you.
Khushi Patel
Khushi Patel le 25 Sep 2020
Oh okay. Greatly appreciate the help and the clarification about plotting k!
Star Strider
Star Strider le 25 Sep 2020
As always, my pleasure!

Connectez-vous pour commenter.

Plus de réponses (1)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by