Systems of first order differential equations:
Equilibrium solutions (Critical point, Steady State Solution)
Types of critical points
Stability of equilibrium solutions

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) and the function y(t) as representing a population of foxes that naturally decline (that is, exponential decay in the absence of prey). 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.

This is a nonlinear system, so an explicit solution is very hard to find.

Systems of first order linear differential equations

Consider a system S of linear differential equations  ({{{, ({{x'=ax+by, (1)}, {y'=cx+dy, (2)}})}})

The functions x and y are functions of time t and the derivatives are also taken with respect to the variable t. The numbers a,b,c,d are real numbers and they are called the coefficients, or the parameters of the system. Our system above can also be written as
 ({{{, ({{x'(t)=ax(t)+by(t), (1)}, {y'(t)=cx(t)+dy(t), (2)}})}})   or  ({{{, ({{(dx/dt)=ax(t)+by(t), (1)}, {(dy/dt)=cx(t)+dy(t), (2)}})}}) ,   but the dependence on t is usually not explicitly written.

Definition: The phase portrait of a systems of linear differential equations is a representative set of its solutions (trajectories, integral curves), plotted as parametric curves (with t as the parameter) on the Cartesian plane xy, tracing the path of each particular solution (x, y) = (x(t), y(t)), −∞ < t < ∞ .

Phase Portrait in Mathematica: The Phase Portrait is obtained in Mathematica using the command StreamPlot[..].

Similar to a vector field (obtained in Mathematica using the command VectorPlot), a phase portrait is a graphical tool to visualize how the solutions of a given system of differential equations would behave in the long run. We can classify the behavior of solutions and the type and stability of the equilibrium solution of a given linear system by looking at the phase portrait and analyzing the properties of the coefficient matrix A = (({{a, b}, {c, d}})).

Advantage: We can solve these systems explicitly!

Equilibrium Solution

Definition: An equilibrium solution of the system S is a point (x,y) where x’=0 and y’=0. An equilibrium solution is a constant solution of the system, and is sometimes also called a critical point, or a steady state solution.

How to find the equilibrium solutions of a system :

For our linear system of differential equations, an equilibrium solution occurs at each solution of the system (of homogeneous algebraic equations) ax+by=0 and cx+dy=0. Equivalently, we need to solve the matrix equation A(({{x}, {y}}))=(({{0}, {0}})), where A=(({{a, b}, {c, d}})) is called the coefficient matrix of the system. As we have seen, such a system has exactly one solution, located at the origin (({{0}, {0}})), if and only if the matrix A is invertible. If A is not invertible, then there are infinitely many solutions.

Recall also that the following statements are equivalent:  
A is invertible ⇔ the determinant of the matrix det (A) is different from 0 ⇔ 0 is an eigenvalue of A.

Classification and Stability of Critical Points

Assume that det(A) ≠ 0, therefore the system has a unique critical point at the origin (0,0). Let λ and μ be the two eigenvalues of the matrix A. We distinguish the following three cases:

Case 1: The eigenvalues λ and μ are real and different.

The general solution of the system of linear differential equations is 
(({{x(t)}, {y(t)}})) = (Subscript[c, 1]) (e^(λt))v + (Subscript[c, 2]) (e^(μt))w, 

where (Subscript[c, 1]) and (Subscript[c, 2]) are constant parameters, and v and w are eigenvectors corresponding to the eigenvalues λ and μ.

The stability type of the critical point (0, 0) is dictated by the sign of the eigenvalues λ and μ.

a) The eigenvalues λ and μ are both negative.
The origin is called a sink (attracting point, nodal sink, asymptotically stable node). All trajectories converge to 0 asymptotically when t→∞.

Example: Consider the system  ({{{, ({{x'=-3x+(Sqrt[2])y, (1)}, {y'=(Sqrt[2])x-2y, (2)}})}})

A={{-3,(Sqrt[2])},{(Sqrt[2]),-2}};
Print["The eigenvalues are: ", Eigenvalues[A]]
Print["The eigenvectors are: ", Eigenvectors[A]]

StreamPlot[{-3x+(Sqrt[2])y,(Sqrt[2])x-2y},{x,-5,5},{y,-5,5},StreamScale->Automatic,StreamPoints->{{{{(1/(Sqrt[2])),1},Red},{{((-1)/(Sqrt[2])),-1},Red},{{-(Sqrt[2]),1},Green}, {{(Sqrt[2]),-1},Green}, , Automatic}}]

b) The eigenvalues λ and μ are both positive.
The origin is called a source (repelling point, unstable node, nodal source). All trajectories diverge away from 0 when t→∞.

Example: Consider the system  ({{{, ({{x'=10x, (1)}, {y'=2y, (2)}})}})

B={{10,0},{0,2}};
Print["The eigenvalues are: ", Eigenvalues[B]]
Print["The eigenvectors are: ", Eigenvectors[B]]

StreamPlot[{10x,2y},{x,-2,2},{y,-2,2},StreamScale->Automatic,StreamPoints->{{{{0,1},Red},{{0,-1},Red},{{1,0},Orange}, {{-1,0},Orange}, , Automatic}}]

c) The eigenvalues λ and μ have different signs, say λ>0 and μ<0. The origin is unstable, and it is called a saddle.

The trajectories given by the eigenvectors of the negative eigenvalue μ<0 initially start at infinitedistance away, but converge to the origin exponentially fast (x(t),y(t)) = (Subscript[c, 2]) (e^(μt))w. The trajectories given by the eigenvectors of the positive eigenvalue λ>0 diverge from the origin exponentially fast (x(t),y(t)) = (Subscript[c, 1]) (e^(λt))v.  Every other trajectory moves towards the origin following the attracting direction w, but never converges to the origin, as it changes direction and diverges from the origin following the repelling direction v.

Example: Consider the system  ({{\[Piecewise], ({{x'=x+y, (1)}, {y'=4x+y, (2)}})}})

The eigenvalues of the coefficient matrix A are real and distinct (3 and -1), so the general solution of the system of differential equations is
(({{x(t)}, {y(t)}})) = (Subscript[c, 1]) (e^(3t))(({{1}, {2}})) + (Subscript[c, 2]) (e^(-t))(({{-1}, {2}})), where (Subscript[c, 1]) and (Subscript[c, 2]) are constant parameters. 
There are two special solutions when (Subscript[c, 1])=0, and respectively when (Subscript[c, 2])=0. We will color these special solutions in red, and respectively in orange.

A={{1,1},{4,1}}; (* coefficient matrix of the system *)
Eigensystem[A]
L=Eigenvectors[A];
Rc1 = {L[[1]],Orange}
Lc1={-L[[1]],Orange}
Rc2 ={L[[2]], Red}
Lc2={-L[[2]],Red}
StreamPlot[{x+y,4x+y},{x,-5,5},{y,-5,5},StreamScale->Automatic,StreamPoints->{{Rc1,Rc2,Lc1,Lc2,Automatic}}]

Case 2: The eigenvalues λ and μ are complex conjugate.

The general solution of the system of linear differential equations is 
(({{x(t)}, {y(t)}})) = (Subscript[c, 1]) (e^at)(u cos(bt) - vsin(bt)) + (Subscript[c, 2]) (e^at)(u sin(bt) + vcos(bt)), 

where (Subscript[c, 1]) and (Subscript[c, 2]) are constant parameters, 
λ = a+ib, μ = a-ib,
and wλ = u+iv and wμ = u-iv are complex eigenvectors corresponding to the complex conjugate eigenvalues λ and μ.

a) The eigenvalues λ and μ have real part Re(λ)=Re(μ)<0.
The origin is called a spiralling sink. All trajectories spiral toward 0 when t→∞.

b) The eigenvalues λ and μ have real part Re(λ)=Re(μ)>0.
The origin is called a spiralling source. All trajectories spiral away from 0 when t→∞.

Example: Consider the system  ({{\[Piecewise], ({{x'=y, (1)}, {y'=-x+y, (2)}})}})

A={{0,1},{-1,1}};
Print["The Eigenvalues are complex conjugate", N[Eigenvalues[A],1]," with positive real part"]
Print["Eigenvectors: ",e=Eigenvectors[A]]

StreamPlot[{y,-x+y},{x,-5,5},{y,-5,5},StreamScale->Automatic,StreamPoints->{{{{1,0},Red},{{-1,0},Red}, {{0,1},Green},{{1,1},Green},{{-3,2},Orange},Automatic}}]

c) The eigenvalues λ and μ have real part Re(λ)=Re(μ)=0.
The origin is called a center. All trajectories are closed circles around the origin.

Example: Consider the system  ({{\[Piecewise], ({{x'=y, (1)}, {y'=-x, (2)}})}})

A={{0,1},{-1,0}};
Print["The eigenvalues are ", Eigenvalues[A]]
Eigenvectors[A]

StreamPlot[{y,-x},{x,-5,5},{y,-5,5},StreamScale->Automatic,StreamPoints->{{{{0,1},Red},{{1,1},Green},Automatic}}]

Case 3: The eigenvalues λ and μ are real and λ = μ .

a) The eigenvalue λ has algebraic and geometric multiplicity 2 (there exists two linearly independent eigenvectors v and w of λ).
The general solution is (({{x(t)}, {y(t)}})) = (Subscript[c, 1]) (e^(λt))v + (Subscript[c, 2]) (e^(λt))w. 
All the theory from Case 1 holds true, and the origin is called a sink if λ<0 and respectively a source if λ>0. In the Phase Portrait, every nonzero solution looks like a straightline in the direction given by the vector (Subscript[c, 1])v + (Subscript[c, 2])w.

Example: Consider the system  ({{\[Piecewise], ({{x'=7x, (1)}, {y'=7y, (2)}})}})

A={{7,0},{0,7}};
Print["The command Eigenvalues returns ", Eigenvalues[A], " so 7 is the only eigenvalue and it has algebraic multiplicity 2"]
Print["The command Eigenvectors returns two eigenvectors",e=Eigenvectors[A],"so 7 has geometric multiplicity 2"]

StreamPlot[{7x,7y},{x,-5,5},{y,-5,5},StreamScale->Automatic,StreamPoints->{{{{1,0},Red},{{-1,0},Red}, {{0,1},Green},{{0,-1},Green},Automatic}}]

b) The eigenvalue λ has algebraic multiplicity 2 and geometric multiplicity 1.
The general solution is (({{x(t)}, {y(t)}})) = (Subscript[c, 1]) (e^(λt))v + (Subscript[c, 2]) ((e^(λt))w + t (e^(λt))v), where v is an eigenvector, and w is a generalized eigenvector (w is called a generalized eigenvector if it a solution of the equation (A-λI)w=v, where v is an eigenvector for the eigenvalue λ).
When λ<0, the origin is called a degenerate sink. All trajectories converge to the origin. With only one linearly independent eigenvector, the Phase Portrait looks like a combination between Case 1(a) and Case 2(a).  The origin is called a degenerate source when λ>0.

Example: Consider the system  ({{\[Piecewise], ({{x'=x+y, (1)}, {y'=y, (2)}})}})

A={{1,1},{0,1}};
Print["The command Eigenvalues returns ", Eigenvalues[A], " so 1 is the only eigenvalue and it has algebraic multiplicity 2"]
Print["The command Eigenvectors returns ",e=Eigenvectors[A], " but the only eigenvector is ", e[[1]]]

StreamPlot[{x+y,y},{x,-5,5},{y,-5,5},StreamScale->Automatic,StreamPoints->{{{{1,0},Red},{{-1,0},Red}, {{0,1},Green},{{1,1},Green},{{-3,2},Orange},Automatic}}]

Interactive examples with StreamPlot

Interactive plot of the vector field of various systems of linear differential equations.

Manipulate[
Row[{
Column[{
Text["  MATRIX  "],
MatrixForm[m],
Text["  EIGENVALUES  "], 
N[Eigenvalues[m],1],
Text["  EIGENVECTORS  "], 
Row[{
MatrixForm[N[Eigenvectors[m][[1]],1]],
MatrixForm[N[Eigenvectors[m][[2]],1]] 
}]
}, Spacings->2], 
StreamPlot[m.{x,y},{x,-1,1},{y,-1,1},StreamScale->Large,StreamColorFunction->"Rainbow",ImageSize->Large]
}],
{{m,({{1,0},{0,2}})},
{({{1,0},{0,2}})->"Nodal source",({{1,1},{0,1}})->"Degenerate source",({{0,1},{-1,1}})->"Spiral source",
({{-1,0},{0,-1}})->"Nodal sink",({{-1,0},{0,-2}})->"Nodal sink",({{-1,1},{0,-1}})->"Degenerate sink",({{0,1},{-1,-1}})->"Spiral sink",({{0,1},{-1,0}})->"Center",({{1,0},{0,-2}})->"Saddle"}}
]