How to plot a phase potrait of a system of ODEs (in Mathematica)

Question

I am asked to use mathematica to plot the fixed-points of the following system

$$ \frac{dN}{dt} = -\gamma N \left( 1 - \left( \beta M + N\right) \right) $$ $$\frac{dM}{dt} = M \left( 1 - \left( \alpha N + M\right) \right) $$

for the case where $\alpha = 2,\ \beta =2,\ \gamma=1$.

I assume it will involve using NDSolve, but I'm not sure what to take for initial conditions or really the format of the statement.

Any help is appreciated

so far I have the following:

sol = NDSolve[{x'[t] == x[t] - \[Sigma] (y[t] x[t]) - x[t]^2, 
   y'[t] == \[Rho] (\[Beta] y[t] x[t]) - \[Rho] y[
       t] - \[Rho] y[t]^2, \[Sigma] == 2, \[Beta] == 2, \[Rho] == 1}, 
  x, y, {t, 0, 50}]

Answer 1

To plot the phase portrait in Mathematica, you can use StreamPlot:

\[Alpha] = 2;
\[Beta] = 2;
\[Gamma] = 1;

p1 = StreamPlot[{-\[Gamma]*x (1 - (\[Beta]*y + x)), 
    y (1 - (\[Alpha]*x + y))}, {x, -0.5, 2}, {y, -0.5, 1.5}];

---

If you want the nullclines as well then do,

\[Alpha] = 2;
\[Beta] = 2;
\[Gamma] = 1;

p1 = StreamPlot[{-\[Gamma]*x (1 - (\[Beta]*y + x)), 
    y (1 - (\[Alpha]*x + y))}, {x, -0.5, 2}, {y, -0.5, 1.5}];

p2 = Plot[{(1 - x)/2, 0}, {x, -0.5, 2}];

p3 = ContourPlot[{x == (1 - y)/2, x == 0}, {x, -0.5, 2}, {y, -0.5, 1.5}];

Show[p1, p2, p3]

Where the lines cross of the same color, you have a fixed point as well as at an orange and blue line but those two should be apparent.

Answer 2

For those that use Matlab,

% m file name odes.m
% change the parameters for your ode where y(1) = x, y(2) = y, and yp denotes
% derivative

function yp = odes(t, y);
a = 2; 
b = 2; 
g = 1; 

yp(1) = -g*y(1)*(1 - (b*y(2) + y(1)));
yp(2) = y(2)*(1 - (a*y(1) + y(2)));
yp = yp';

---

% m file name phase_portrait.m
% Phase Plot Program to use this function, do the following:   
% >> phase_portrait(x1, x2, y1, y2, tfinal, 'F', N); 
% In this example, I ran  phase_portrait(-0.5, 2, -0.5, 2, 5, 'odes', 8)


function [] = phase_portrait(x1, x2, y1, y2, tfinal, odes, N); 

%   x1 is the x-min value
%   x2 is the x-max value
%   y1 is the y-min value
%   y2 is the y-max value
%   tfinal is the length of time interval
%   F is the system function input as a string 'F'
%   NxN: is the number of orbits plotted

figure; 
hold on;
axis([x1 x2 y1 y2]);

Options = odeset('RelTol', 1e-6, 'AbsTol', [1e-10 1e-10]);

dx = (x2 - x1)/N; 
dy = (y2 - y1)/N;
x = x1:dx:x2; 
y = y1:dy:y2;

for i = 1:length(x)
    for j = 1:length(y)
        Y0 = [x(i); y(j)];
        grad = feval(odes, 0, Y0);
        plot_arrow(Y0, grad, x2 - x1, y2 - y1,'r');
        [T, Y] = ode45(odes, [0 tfinal], Y0, Options);
        plot(Y(:, 1), Y(:, 2));
        [T, Y] = ode45(odes, [0 -tfinal], Y0, Options);
        plot(Y(:, 1),Y(:, 2));
    end
end

function [] = plot_arrow(pnt, grd, xlngth, ylngth, c);

% put arrowhead on trajectory to indicate direction
% arrowhead is not scaled to indicate velocity
% so all arrowheads are the same size
% pnt is a point as a vector, [x y]
% grd is the gradient of the trajectory at pnt, e.g.[1 .5]
% c is the color of the arrowhead as a string

nrm = norm(grd);                          % get norm of gradient
scaler = min(abs(xlngth), abs(ylngth));   % get length of shorter axis

if nrm > 1e-6                             % don't put arrows at fixed points
    grd = .02*scaler*grd./nrm;            % scale norm of gradient to axes
    A = [0, -1; 1, 0];                    % 90 degree rotation matrix
    rgrd = A*grd;                         % a vector perp to gradient
    h = pnt + .5*grd;                     % the point of the arrow head
    tb = pnt - .5*grd + .5*rgrd;          % the top of the base of the arrow head
    bb = pnt - .5*grd - .5*rgrd;          % the bottom of the base of the arrow head
    xs = [h(1); tb(1); bb(1)];            % x values of arrow head vertices
    ys = [h(2); tb(2); bb(2)];            % y values of arrow head vertices
    patch(xs, ys, c);                     % draw the arrow head
end