How can I plot a three-species Lotka-Volterra model?

Hi everyone!

I am a student of cultural studies and thus really out of my depth with MatLab. I have a seminar on Modeling Systems using linear and non-linear differential equations. I am now doing my exam project but I am quite confused by MatLab (this program is completely new to me).

I am trying to model a Lotka-Volterra system with three "species". I have already found my interaction parameters (using the least squares method). The equations now look like this:

dx(t)/dt = x(t)*(a0+a1*x(t)+a2*y(t)+a3*z(t))

dy(t)/dt = y(t)*(b0+b1*x(t)+b2*y(t)+b3*z(t))

dz(t)/dt = z(t)*(c0+c1*x(t)+c2*y(t)+c3*z(t))

a0, a1, a2, a3, b0, b1, b2, b3, c0, c1, c2, c3 are the parameters and constants. I have the values for them. I have tried to figure out how to write this function as a MatLab function but I have failed so far (and thus I cannot plot this system).

I tried naming x, y and z as the first three values of a vector (v), so that I have a function that looks a little like this:

function LV=lotka (t,v)

LV(1)=v(1)*(.......);
LV(2)=v(2)*(....);
LV(3)=v(3)*(.....);

end

I also used "ode45" (where I specified the time interval t) and the "figure" command to try to plot this but something along the way failed. Any help with writing this function would be appreciated.

I hope I made the problem clear. If needed I can provide the full code I have tried so far. However, as it is obviously wrong, I don't think it will help here.

Thank you!

Answer 1

[William Rose]

@Hagen Gersie, The answer from @Sam Chak is excellent . It has the virtues of conciseness and clarity.

You may be wondering why your way did not work. It is clear that you made a solid attempt. Your function lotka.m is actually a very good start. You got farther than a lot of people would. Please share your a,b,c,d constants, and the rest of your code, if you want to see how close you did or did not come to succeeding.

[William Rose]

@Hagen Gersie,

One problem I notice is that your lotka() returns a row vector. ode45() needs it to be a column vector. You can force LV to be a column vector by including the command

LV=zeros(3,1); 

to pre-allocate the vector as a column, before you compute the elements of LV.

Or you can do this. Notice the semicolons inside the [...], to force new rows.

LV = [v(1)*(a0+a1*v(1)+a2*v(2)+a3*v(3));...
      v(2)*(b0+b1*v(1)+b2*v(2)+b3*v(3));...
      v(3)*(c0+c1*v(1)+c2*v(2)+c3*v(3))];

Try it.

[Hagen Gersie]

Thanks for your answer! Already helped a lot. For the first time a plotted graph showed up (even though it doesn't make any sense whatsoever, so apparently I got still some mistakes). The problem could of course also be the underlying data (maybe I have to use a bigger dataset and see that the data fits more "harmoniously" if that makes sense).

This is my current function file (with the error report which doesn't seem to make any sense but maybe I am missing something? I thought the time span would encompass 0 to 100 with steps of 1 in between):

tspan = [0 100];
init = [12.72 141.5 5.7];
[t,v] = ode45(@lotka1,tspan,init);

Warning: Failure at t=3.879637e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.105427e-15) at time t.

plot(t,v)
title ('Interaction');
legend('share price','percentage of float shorted', 'days to cover');
xlabel('Time');
ylabel('Population');

function LV=lotka1(t,v)
a0=-0.5461; a1=8.2919e-04; a2=0.0238; a3=-0.1979;
b0=1.3843; b1=-0.0048; b2=-0.0068; b3=-0.0644;
c0=0.1534; c1=-0.0071; c2=-0.0289; c3=0.6837;

LV=zeros(3,1);

LV(1)=v(1)*(a0+a1*v(1)+a2*v(2)+a3*v(3));
LV(2)=v(2)*(b0+b1*v(1)+b2*v(2)+b3*v(3));
LV(3)=v(3)*(c0+c1*v(1)+c2*v(2)+c3*v(3));

end

[William Rose]

@Hagen Gersie, here is a script that uses your lotka() function to solve the system and plot results.

It makes the plot below. When you run the script and generate the plot, you will be able to click and drag within the 3D plot, to rotate it.

[Hagen Gersie]

Thanks! You are a treasure!

[William Rose]

@Hagen Gersie,

Your script stopped at 3.9 seconds because the equations were blowing up, and the integration routine, ode45() could not complete another time step. This is a result of the particular values for the constants a0..a3, b0..b3, c0..c3, and the initial condition.

[Hagen Gersie]

Thanks, I figured that out as well and adapted accordingly. It now runs smoothly.

Answer 2

[Sam Chak]

Hi @Hagen Gersie

Read doc ode45.

You can also follow the example and do something like this:

sigma = 10;
beta = 8/3;
rho = 28;

% Ordinary differential equations of the Lorenz system
f = @(t, x) [-sigma*x(1) + sigma*x(2); ...
	rho*x(1) - x(2) - x(1)*x(3); ...
	-beta*x(3) + x(1)*x(2)];

tspan = [0 100];

init = [1 1 1];		% initial condition

% Runge-Kutta Dormand-Prince 4/5 solver
[t, x] = ode45(f, tspan, init);

plot3(x(:,1), x(:,2), x(:,3))	% plot lines in 3D space
view(30, 30)		% adjust view angles (in deg)

[Hagen Gersie]

Thanks, Sam! Your answer helped a lot already. I have read the ode45 explainer but as it only ever mentioned the two variables (at some point there even is a example of a lotka-volterra equation but it only uses two species), I had to figure out the rest for myself and I think I made too many mistakes while doing it.

Support's super appreciated.

[Sam Chak]

You are welcome, @Hagen Gersie. Please vote 👍 if you like the support, but consider accepting ✓ @William Rose's Answer for his kind effort in correcting and writing the MATLAB code for you. Have a nice day.

[William Rose]

@Sam Chak, @Hagen Gersie, thank you.