Example for system of ODEs: rabbits and foxes

We already discussed this problem on February 26.

clearvars

You first have to download vectorfield.m and put this file in the folder where your mlx files are.

Population of rabbits and foxes

Let

denote the population of rabbits, let

denote the population of foxes at time t.

We have for the rates of change

Rabbits (prey)

where

Foxes (predator)

where

We obtain the nonlinear ODE system

syms x y
f = (2-y)*x;
g = (-1+x)*y;

Find all stationary points:

Find

with

[xs,ys] = solve([f;g]==[0;0],[x;y])

xs =

ys =

m = length(xs) % number of solutions we found

m = 2

We found the two stationary points (0,0) and

Find Jacobian matrix J:

J = jacobian([f;g],[x;y])

J =

For each stationary point : find matrix

Find eigenvalues and eigenvectors

For LINEARIZED PROBLEM: find type and stability of stationary point, sketch phase portrait

Result:

(0,0) is saddle (unstable, not repelling)
(1,2) is center (stable, not attracting)

for k=1:m

stationary_point = [xs(k),ys(k)]

A = subs(J,{x,y},{xs(k),ys(k)}) % evaluate Jacobian at stationary point

[eigenvectors,eigenvalues] = eig(A)

figure

phaseportrait(A);

title(sprintf('linearized problem for stationary point (%g,%g)',stationary_point))

end

stationary_point =

A =

eigenvectors =

eigenvalues =

stationary_point =

A =

eigenvectors =

eigenvalues =

For NONLINEAR PROBLEM: What can you say about type and stability of the stationary points?

stationary point (0,0): the linearized problem has a saddle point (unstable, not repelling). Therefore also the nonlinear problem has a saddle point (unstable, not repelling), and we get the same tangential directions for the trajectories given by the eigenvectors
stationary point (1,2): the linearized problem has a center (stable, not attracting). But this does not allow us to draw conclusions about the type and stability for the nonlinear ODE system.

Plot phase portrait using numerical ODE solver ode45

We use

in place of

and define the function

as follows:

f = @(t,y) [ (2-y(2))*y(1) ; (-1+y(1))*y(2) ];

For many initial points

: use ode45 and

solve ODE with initial condition for t going from 0 to 5, plot trajectory in phase plane
solve ODE with initial condition for t going from 0 to , plot trajectory in phase plane

for a1 = 0:.25:2

for a2 = 0:.25:3

[ts,ys] = ode45(f,[0,5],[a1;a2]); plot(ys(:,1),ys(:,2),'c'); hold on

[ts,ys] = ode45(f,[0,-5],[a1;a2]); plot(ys(:,1),ys(:,2),'c')

end

vectorfield(f,0:.25:2,0:.25:3);

axis([0 2 0 3]); hold off

We can see

saddle at (0,0)
center at (1,2)

We already knew that we have a saddle at (0,0).

But how can we prove that we have a center at (1,2)? See next section about "integral methods".

Note that we have closed trajectories corresponding to periodic solutions, see discussion from February 26.

Analyzing the autonomous ODE system with integral methods

The ODE system

states that the tangential direction

is orthogonal to the vector

, i.e.,

This is called orbital equation. Formally multiplying this by dt gives the orbital equation in the form

This is a first order ODE for the trajectory curves in the

plane.

If we can use integration to turn this into an equation

we say that the ODE is conservative.

This is true if

the ODE is exact, i.e.,
after multiplication by an integrating factor we obtain an exact ODE.

We know two special cases where we can obtain obtain an exact ODE by multiplying the ODE with a factor:

if the ODE is linear
if the ODE is separable

For our rabbits & foxes problem the orbital equation

This is clearly a separable ODE: dividing by xy gives

Note that we only consider

We obtain

with the function

Use Matlab to solve the orbital equation and find

syms y(x) y1 C1
f = (2-y)*x;
g = (-1+x)*y;

Solve the orbital equation

sol = simplify( dsolve(-g+f*diff(y,x)==0,'Implicit',true) )

sol =

The first entry sol(1) is the implicit solution we want. The second entry states that the horizontal line

is a solution.

H = rhs(isolate(sol(1),C1)) % rewrite as C1=...

H =

H = subs(H,y(x),y1) % write y1 in place of y(x)

H =

We can now make a contour plot for

using a step size of 0.1 for the contour values.

fcontour(H,[0 2 0 3],'LevelStep',.1) % make contour plot for x in [0 2], y1 in [0 3]

Note that these curves are the same as in the phase portrait computed with ode45.

Function for drawing phase portrait

function phaseportrait(A)
syms y(t) [2 1]
syms C1 C2
sol = dsolve(diff(y)==A*y);             % solve ODE 
ysol = subs(y(t),sol);
f = @(t,y) A*y;                         % define function f(t,y) for vectorfield
vectorfield(f,-5:5,-5:5); hold on
for c1=[-1 0 1]
  for c2=[-1 0 1]
    ys = subs(ysol,{C1,C2},{c1,c2});    % use C1=c1, C2=c2
    fplot(ys(1),ys(2),[-5 5]);          % plot trajectory for t=-5...5
  end
end
hold off; axis equal; axis([-5 5 -5 5]); xlabel('y_1'); ylabel('y_2')
end