Answers to the exercises from the last lecture:

Exercise 1:

Consider the system  ({{{, ({{x'= y, (1)}, {y'=-y-2sin(x), (2)}})}}).  Study the nature of the critical points (x, y) of the system for which -2π < x,y < 2π. Use only Mathematica commands, even when computing the critical points!

Solution:

StreamPlot[{y ,-y-2Sin[x]}, {x,-2Pi,2Pi}, {y,-2Pi, 2Pi}, StreamPoints->Fine]

The Phase Portrait indicates 5 critical points, all located on the x - axis. Two of them are located exactly on the vertical boundaries of the visualization frame, so they might be harder to see. We can redraw the Phase Portrait on a slightly larger interval for the x-coordinate, like [-3 Pi, 3 Pi].

StreamPlot[{y ,-y-2Sin[x]}, {x,-3Pi,3Pi}, {y,-2Pi, 2Pi}, StreamPoints->Fine]

We find the critical points by setting x’=y’=0 and using Solve[..] or NSolve[..] to solve the resulting system of algebraic equations.

Clear[x,y]
criticalpoints=Solve[{y==0 && -y-2*Sin[x]==0}, {x,y}]

The output of Solve should be understood as follows : all the critical points are of the form {x, y}, with y = 0 and x = 2 π C[1] or π + 2 π C[1], where C[1] is any integer number. Therefore x can be any odd or even integer of π . To analyze the stability of the critical points, it is enough to consider C[1]=1, since the Phase Portrait is periodic with period 2π with respect to the x-coordinate.

Clear[h,J,x,y]
h[x_,y_]:={y,-y-2*Sin[x]}
(* J[f_] computes the Jacobian matrix of a function f:(ℝ^2)→(ℝ^2) *)
J[f_]:=
Module[{A=Table[0,2,2]},
A[[1,1]]=D[f[x,y][[1]],{x}];
A[[1,2]]=D[f[x,y][[1]],{y}];
A[[2,1]]=D[f[x,y][[2]],{x}];
A[[2,2]]=D[f[x,y][[2]],{y}];
Return[A]]
Print["The Jacobian matrix of the function ",h[x,y], " is ", MatrixForm[J[h]]]

We can now compute the eigenvalues of the Jacobian matrix at the critical points, by using the replacement operator /. (slash - dot).

B1=J[h]/. criticalpoints[[1]] /. C[1]->1;
MatrixForm[B1]
Eigenvalues[B1]

In conclusion, at critical points of the form {x, y} = {0, 2π C[1]}, where C[1] is any integer number, the eigenvalues of the Jacobian matrix are complex conjugate with negative real parts. Therefore these critical points are spiral sinks. They are asymptotically stable.

B2=J[h]/. criticalpoints[[2]] /. C[1]->1;
MatrixForm[B2]
Eigenvalues[B2]

In conclusion, at critical points of the form {x, y} = {0, π +2π C[1]}, where C[1] is any integer number, the eigenvalues of the Jacobian matrix are real and distinct and with opposite signs. Therefore these critical points are saddle fixed points. They are unstable.

S=StreamPlot[{y ,-y-2Sin[x]}, {x,-3Pi,3Pi}, {y,-2Pi, 2Pi}, StreamPoints->{{{{1,0},Red}, {{-1,0},Red}, {{-2.3,0},Red},{{2.5,-1},Red}, {{2.5,0},Red},{{5.38,-4},Green},{{-0.9975,6},Green}, Automatic}}];
saddlepoints = {x,y} /. criticalpoints[[2]] /. {{C[1]->0}, {C[1]->-1}};
G1=Graphics[{PointSize[Large],Green, Point[saddlepoints]}];
spiralsinks = {x,y} /. criticalpoints[[1]] /. {{C[1]->-1}, {C[1]->0},{C[1]->1}};
G2=Graphics[{PointSize[Large],Red, Point[spiralsinks]}];
Show[S,G1,G2]

Based on the pictures and the analysis of the equilibrium points, we can deduce the following : almost all trajectories converge to the spiral sinks. There are some exceptional trajectories going into and coming out of the saddle fixed points (if the system were linear with constant coefficients, they would corresponding to the eigenvectors of the Jacobian matrix), which are called separatrices, because they have exceptional behavior and separate the regions of generic behavior from each other. The seaparatrices are harder to visualize on a Stream Plot (the green curve represents an approximation of one of the separatrices).

Exercise 2 (predator-prey system):

Consider the system usually described as a predator - prey model (Voltera, 1931):
  (dx/dt)=0.2x-0.7xy,  (dy/dt)=-0.4y+0.5xy.

where you can think of the function x(t) as representing a population of rabbits that naturally grow at a rate proportional to their population (that is, exponential growth in the absence of predators, at the rate (dx/dt)=0.2x) and the function y(t) as representing a population of foxes that naturally decline (that is, exponential decay in the absence of prey, at a rate (dy/dt)=-0.4y). The term xy in both equations is proportional to the number of likely encounters of the two population in a certain environment. Encounters are detrimental to the rabbits and beneficial to the foxes.

Plot the vector field using StreamPlot[...] in the region that makes sense for the given problem. On the same plot, plot the trajectory starting at the point (0, 1) in Red, the trajectory starting at the point (1, 0) in Green, the trajectory starting at (1,1) in Orange, and the trajectory starting at (5, 2) in Purple. What do you notice? How do the solutions look like?

Find the critical points of the system and describe their nature.

Assume that initially there are 5 rabbits and 2 foxes. Use NDSolve[...] or DSolve[...] to find the solution corresponding to this initial condition. Try plotting the solution that NDSolve returns. What does this say about the long term dynamics of the rabbit/fix population?

s=StreamPlot[{0.2x-0.7x y,-0.4y+0.5x y},{x,-1,1},{y,-1,1},StreamScale->Automatic,StreamPoints->{{{{0,1},Red},{{0,-1},Red},{{1,0},Green},{{-1,0},Green},{{1,1},Orange}, {{5,2},Purple},Automatic}}]

It is very easy to see that the origin is a saddle fixed point, even without coding in Mathematica. One of the (possibly multiple) solutions of the equation 
 0.2x-0.7xy=0,  -0.4y+0.5xy=0. is x=y=0. 
 The linear part of the nonlinear system at (0,0) is simply given by (dx/dt)=0.2x,  (dy/dt)=-0.4y.

The Jacobian matrix at (0, 0) is therefore (({{0.2, 0}, {0, -0.4}})) and has eigenvalues 0.2 and -0.4. Therefore the origin is a saddle fixed point. In this case, it is very easy to identify the exceptional trajectories going into and coming out of the saddle fixed points: they are located on the y-axis, and respectively on the x-axis. The positive y-axis describes the evolution of a habitat with y foxes and 0 rabbits. In the absence of rabbits, the population of foxes converges to 0. The positive x-axis describes the evolution of a habitat with 0 foxes and x rabbits. In the absence of foxes, the population of rabbits grows to infinity.

Since we are interested in the dynamics of the populations of rabbits and foxes, we don' t care much for the negative values, so we will just focus on the first quadrant.

s1=StreamPlot[{0.2x-0.7x y,-0.4y+0.5x y},{x,0,1},{y,0,1},StreamScale->Automatic,StreamPoints->{{{{0,1},Red},{{1,0},Green},{{1,1},Orange}, {{5,2},Purple},Automatic}}];
s2= StreamPlot[{0.2x-0.7x y,-0.4y+0.5x y},{x,0,9},{y,0,9},StreamScale->Automatic,StreamPoints->{{{{0,1},Red},{{1,0},Green},{{1,1},Orange}, {{5,2},Purple},Automatic}}];
GraphicsRow[{s2,s1},ImageSize->{900,400}]

critpts=NSolve[0.2x-0.7x y==0 && -0.4y+0.5x y == 0]

There are two critical points, at (0, 0) and (0.8, 0.285714). We have already studied the eigenvalues of the Jacobian matrix of the system at (0, 0), and decided that the origin is a saddle fixed point. We will now do the same thing for the second critical point.

f[x_,y_]:={0.2x-0.7x y,-0.4y+0.5x y};
Print["The Jacobian matrix of f is "]
MatrixForm[J[f]]
Print["The Jacobian matrix of f at ", {x,y}/. critpts[[1]], "is"]
M=J[f]/. critpts[[1]] ;
MatrixForm[M]
Print["and the eignevalues are ..."]
Eigenvalues[M]

These eigenvalues are computed numerically. They are complex conjugate and the real part - 1.3877787807814457`*10^-17, which is equivalent of 0 (with respect to a numerical precision of 17 digits). The real part of the eigenvalues is in fact exactly 0, so the linear system defined by the Jacobian matrix has a center at (0.8, 0.285714). The Phase Portrait of the linear system is depicted below. However, the general theory doesn’t allow us to immediately conclude that (0.8, 0.285714) is a center for the nonlinear system as well. We could have three possibilities: spiral sink, spiral source, or center, depending on the nonlinear terms.

linear=StreamPlot[M.{x,y},{x,-1,1},{y,-1,1},StreamScale->Automatic]

If we take a look again at the Phase Portrait of the nonlinear system, it suggests that we have a lot of closed trajectories in the first quadrant, which is a strong indicator of the fact that (0.8, 0.285714) is a center for the nonlinear system as well.

The computations below show the particular solution of the system of differential equations corresponding to an initial population of 5 rabbits and 2 foxes. The choice of {t,0,100} in NDSolve[..] means that we solve x and y as functions of time t, and we are able to predict the functions for 100 units of time.

Clear[x,y,x1,y1,t];
sol5=NDSolve[{x'[t]==0.2* x[t]-0.7*x[t]*y[t], y'[t]==-0.4 *y[t]+0.5 *x[t]* y[t], x[0]==5, y[0]==2},{x[t],y[t]},{t,0,100}, MaxSteps->10000]

{x1[t_],y1[t_]}= {x[t] ,y[t]}/. Flatten[sol5];
ParametricPlot[{x1[t],y1[t]},{t,0.1,100}]
Plot[{x1[t],y1[t]},{t,0,100}]

It' s very hard to realize that this will actually be a closed trajectory if we let t increase. It is however true, as we saw from the StreamPlot[..]. We need to prove this mathematically. We will turn the system of two differential equations into one differential equation, which models y as a function of x.

(dx/dt)=0.2x-0.7xy,  (dy/dt)=-0.4y+0.5xy.

We will compute

(dy/dx)=((dy/dt)/(dx/dt))=((-0.4y+0.5xy)/(0.2x-0.7xy))

The new differential equation contains less information than the initial system, as it only models the population of foxes (y) as a function of the population of rabbits (x), but does not describe how both populations evolve as functions of time (t). Nonetheless, this is a separable differential equation that we can solve explicitly. We are only interested in positive values of x and y so we may assume that x>0 and y>0. We first separate the variables and write the equation as:

((0.2-0.7y)/y)dy=((-0.4+0.5x)/x)dx

Now we can integrate the left hand side with respect to y and the right hand side with respect to x. We obtain

⇒ 0.2 ln(y)-0.7y=-0.4ln(x)+0.5x+c
 ⇒ ln ((y^0.2)(x^0.4))=0.5x+0.7y+c
 ⇒ (y^0.2)(x^0.4)(e^(-(0.5x+0.7y)))=c

where C=(e^c) is a random positive constant. Therefore the trajectories of the Voltera system in the first quadrant are the level curves of the function F(x,y)=(y^0.2)(x^0.4)(e^(-(0.5x+0.7y))). 
In conclusion, to understand the trajectories (but not how a point moves along a trajectory as a function of time t), all we need to do is understand the function F.

The trajectory corresponding to an initial population of x = 5 rabbits and y = 2 foxes is given by F(x,y)=C, where C=F(5,2).

Clear[x,y,F]
F[x_,y_]:=y^0.2*x^0.4*Exp[-(0.5 x+0.7 y)]
F[5,2]
ContourPlot[F[x,y]==F[5,2],{x,0,7},{y,0,7} ]

We can plot several other trajectories as well:

ContourPlot[{F[x,y]==0.1, F[x,y]==0.2, F[x,y]==0.3},{x,0,7},{y,0,7} ]