%% 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 x(t) denote the population of rabbits, let y(t) denote the population of foxes at time t.
% We have for the rates of change x'(t),\; y'(t)
% Rabbits (prey)
%      x\prime(t) = R_1\cdot x  where R_1 = 2-y
% Foxes (predator)
%       y\prime(t)=R_2\cdot y  where R_2 = -1 + x
% We obtain the nonlinear ODE system \left[ \matrix{
%        x\prime \cr  y\prime
% }\right]
% =
% \left[ \matrix{
%    f(x,y) \cr g(x,y)
% }\right]
%          \left[ \matrix{
%        x\prime \cr  y\prime
% }\right]
% =
% \left[ \matrix{
%    (2-y)x \cr (-1+x)y
% }\right]

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

%% Find all stationary points:
% Find \left[ \matrix{
%        x \cr  y
% }\right]  with \left[ \matrix{
%        f(x,y) \cr  g(x,y)
% }\right]
% =
% \left[ \matrix{
%        0 \cr  0
% }\right] :

[xs,ys] = solve([f;g]==[0;0],[x;y])
m = length(xs)          % number of solutions we found

% We found the two stationary points (0,0) and (1,2).
%% Find Jacobian matrix J:

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

%% For each stationary point (x_*,y_*): find matrix A=J(x_*,y_*)
%% 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

%
%% 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 \left[ \matrix{1 \cr 0} \right],\;
% \left[ \matrix{0 \cr 1} \right]
% 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 y_1, y_2 in place of x,y and define the function \vec f(t,\vec y) as follows:

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

% For many initial points \vec a: use ode45 and
% solve ODE with initial condition \vec y(0)=\vec a for t going from 0 to 5, plot trajectory in phase plane
% solve ODE with initial condition \vec y(0)=\vec a for t going from 0 to -5, 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
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 \left[ \matrix{
%        x\prime \cr  y\prime
% }\right]
% =
% \left[ \matrix{
%    f(x,y) \cr g(x,y)
% }\right]  states that the tangential direction \left[ \matrix{
%        x\prime \cr  y\prime
% }\right]
%  is orthogonal to the vector \left[ \matrix{
%   -g(x,y) \cr f(x,y)
% }\right] , i.e.,
%           -g(x,y) \frac{dx}{dt}+f(x,y)\frac{dy}{dt} = 0
% This is called orbital equation. Formally multiplying this by dt gives the orbital equation in the form
%%           -g(x,y) \,dx+f(x,y)\,dy = 0
% or
%%            -g(x,y)+f(x,y)\,y\prime(x)=0
% This is a first order ODE for the trajectory curves in the x,y plane.
% If we can use integration to turn this into an equation
%%             H(x,y)=C
% we say that the ODE is conservative.
% This is true if
% the ODE -g(x,y) \,dx+f(x,y)\,dy = 0 is exact, i.e., \partial_y\left[-g(x,y)\right] = \partial_x f(x,y)
% 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 -g(x,y) \,dx+f(x,y)\,dy = 0 is
%%          (1-x)y\,dx + (2-y)x\,dy = 0
% This is clearly a separable ODE: dividing by xy gives
%            (x^{-1}-1)\,dx+(2y^{-1}-1)\,dy=0  \\
% \int(x^{-1}-1)\,dx+\int(2y^{-1}-1)\,dy=C \\
% \bigl(\log x-x\bigr) + \bigl( 2\log y-y \bigr)= C
% Note that we only consider x>0, y>0.
% We obtain H(x,y)=C with the function
%%           H(x,y) = \log x-x+2\log y-y
%% Use Matlab to solve the orbital equation and find H(x,y)

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

% Solve the orbital equation -g(x,y)+f(x,y)y\prime(x)=0:

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

% The first entry sol(1) is the implicit solution we want. The second entry states that the horizontal line y(x)=0 is a solution.

H = rhs(isolate(sol(1),C1))            % rewrite as C1=... 
H = subs(H,y(x),y1)                    % write y1 in place of y(x)

% We can now make a contour plot for x\in[0,2], y\in[0,3] 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