Population dynamics with one or two species - phase plot

Question

I am working on two equations. The logistic equation $$x'_i=r_i x_i \left(1-\frac{x_i}{K_i}\right) \tag{1}$$ and a two-species competition model $$ \begin{aligned} x'_1 &= r_1 x_1 \left(1-\frac{x_1}{K_1}\right) - \alpha_1 x_1 x_2\\ x'_2 &= r_2 x_2 \left(1-\frac{x_2}{K_2}\right) - \alpha_2 x_2 x_1 \end{aligned} \tag{2} $$ My assignment says

Write two Matlab codes: one that can solve a scalar ODE and one that can solve a system of two ODEs. Let both programs to choose several initial data points and draw corresponding integral curves both as functions of time. In the case of the system of two equations the code must draw trajectories of $(x_1(t); x_2(t))$ in the phase plane of $(x_1; x_2)$ in a separate figure.

Can I use quiver to draw a phase plot for equation $(2)$? So does quiver plot it as a function of time?

xdom = linspace(-5,5,51);
ydom = linspace(-5,5,51);
[X,Y] = meshgrid(xdom,ydom); % generate mesh of domain
U = equation here for i=1; % dx/dt 
V = equation here for i=2; % dx/dt 
quiver(X,Y,U,V)

What changes do I need to make it run for this equation? Any help is appreciated. Thank you.

Answer 1

In both cases, the question asks to compute trajectories numerically e.g. using ode45 for several initial data. Then, we are asked to plot the numerical solutions $(t; x_1(t))$, $(t; x_2(t))$ so-obtained with respect to time. For the competition model, it is possible to superimpose a quiver stream plot with such trajectories $(x_1(t); x_2(t))$ in the phase diagram. Here is a runable Matlab script:

%parameters
r1 = 0.1;
r2 = 0.2;
a1 = 0.3;
a2 = 0.1;
K1 = 0.5;
K2 = 0.5;
tfin = 100;
x10 = [0.1, 0.2, 0.3];
x20 = [0.2, 0.2, 0.2];

% logistic equation trajectories
fun = @(t,x) r1*x.*(1 - x/K1);
tsol = []; xsol = [];
for i=1:length(x10)
    [t,x] = ode45(fun,[0 tfin],x10(i));
    tsol = [tsol t];
    xsol = [xsol x];
end
figure;
plot(tsol,xsol,'b-');
xlabel('t');
ylabel('x');

% competition model trajectories
fun = @(t,x) [r1*x(1,:).*(1 - x(1,:)/K1) - a1*x(1,:).*x(2,:); ...
    r2*x(2,:).*(1 - x(2,:)/K2) - a2*x(1,:).*x(2,:)];
tsol = []; x1sol = []; x2sol = [];
for i=1:length(x10)
    [t,x] = ode45(fun,[0 tfin],[x10(i) x20(i)]);
    tsol = [tsol t];
    x1sol = [x1sol x(:,1)];
    x2sol = [x2sol x(:,2)];
end
figure;
plot(tsol,x1sol,'b-');
hold on
plot(tsol,x2sol,'r-');
xlabel('t');
ylabel('x');

% phase diagram with quiver and trajectories
x1min = min(min(x1sol));
x1max = max(max(x1sol));
x2min = min(min(x2sol));
x2max = max(max(x2sol));
figure;
[x1grid,x2grid] = meshgrid(linspace(x1min,x1max,20),linspace(x2min,x2max,20));
xp1 = r1*x1grid.*(1 - x1grid/K1) - a1*x1grid.*x2grid;
xp2 = r2*x2grid.*(1 - x2grid/K2) - a2*x2grid.*x1grid;
quiver(x1grid,x2grid,xp1,xp2);
hold on
plot(x1sol,x2sol,'k-');
xlabel('x1');
ylabel('x2');
xlim([x1min x1max]);
ylim([x2min x2max]);

Output:

Answer 2

The quiver function plots vector $(U,V)$ at gridpoint $(X,Y)$. So, in your code above, $U(X) = (\frac{dx(t)}{dt})_{x=X}$ and $V(Y) = (\frac{dy(t)}{dt})_{y=Y}$.

The quiver plot illustrates how one component of your time-dependent solution varies relative to the other component. To create a phase plot, you would need to overlay the quiver plot with trajectories. (see http://tutorial.math.lamar.edu/Classes/DE/PhasePlane.aspx ).

But, the assignment you listed appears to have nothing to do with phase plots. Rather, it seems like you just need to solve some ODEs and plot the solution(s) as functions of time. So, for a given ODE system, just select initial conditions, solve it and graph the solutions as functions of time on the same plot.