Nonlinear Differential Equations:

Consider the system  ({{{, ({{x'= f(x,y), (1)}, {y'= g(x,y), (2)}})}})  where f and g are functions of two variables x and y. Assume that f and g are not both linear.

The functions x and y are functions of time t and the derivatives are also taken with respect to the variable t.  ({{{, ({{x'(t)=f(x(t),y(t)), (1)}, {y'(t)=g(x(t),y(t)), (2)}})}}) .

A linear system  ({{{, ({{x'= ax+by, (1)}, {y'= cx+dy, (2)}})}})  whose coefficient matrix (({{a, b}, {c, d}})) is invertible has only one equilibrium point (critical point, steady state solution) at (x,y)=(0,0). Analyzing the stability of this equilibrium solutions gives complete information about the Phase Portrait of the ODE (the trajectories).

Unlike a linear system, a nonlinear system could have none, one, two, three, or any number of critical points. Like a linear system, however, the critical points are found by setting x’ = y’ = 0, and solve the resulting system f(x,y)=0 and g(x,y)=0. 
Since there might be multiple critical points present on the phase portrait, each trajectory could be influenced by more than one critical point. This results in a much more chaotic appearance of the phase portrait. Consequently, the type and stability of each critical point need to be determined locally (in a small neighborhood on the phase plane around the critical point in question) on a case-by-case basis. The trajectories will appear very differently than those of the linear systems, except very near the critical points.

General Strategy :

We approximate the nonlinear system by a linear system near the equilibrium points. 
Suppose that (a, b) is an equilibrium point, that is f(a, b) = g(a, b) = 0.

({{{, ({{x'= f(x,y) ≃ f(a,b) + (∂_x)f(a,b)(x-a) + (∂_y)f(a,b)(y-b)  = (∂_x)f(a,b)(x-a) + (∂_y)f(a,b)(y-b) , (1)}, {y'= g(x,y) ≃ g(a,b) + (∂_x)g(a,b)(x-a) + (∂_y)g(a,b)(y-b)  = (∂_x)g(a,b)(x-a) + (∂_y)g(a,b)(y-b) , (2)}})}})

The critical point can be translated to 0 via the coordinate change x := x - a and y := y - b and still retains all of its properties. The approximated system becomes a linear system:

({{{, ({{x'= (∂_x)f(a,b)x + (∂_y)f(a,b)y, (1)}, {y'=(∂_x)g(a,b)x + (∂_y)g(a,b)y , (2)}})}})   with coefficient matrix  (Subscript[J, (a,b)])=(({{(∂_x)f(a,b),  (∂_y)f(a,b)}, {(∂_x)g(a,b), (∂_y)g(a,b)}})) .

The matrix (Subscript[J, (x,y)])=(({{(∂_x)f(x,y),  (∂_y)f(x,y)}, {(∂_x)g(x,y), (∂_y)g(x,y)}}))is called the Jacobian matrix of the system. It contains the first order partial derivatives of f and g evaluated at the point (x,y).

Stability of equilibrium points:

We can use the eigenvalues of the matrix Jacobian matrix to decide the type and stability of the critical point (a,b). Let λ and μ be the two eigenvalues. Assume that λ and μ are different from 0.

({{Eigenvalues, Linear System, Nonlinear System}, {λ ,μ   Real, □, □}, {λ >μ>0, Nodal Source (Unstable), Nodal Source (Unstable)}, {λ <μ<0, Nodal Sink(Asymptotically Stable), Nodal Sink( Asymptotically Stable)}, {λ >0>μ, Saddle point (Unstable), Saddle point (Unstable)}, {λ =μ>0, Degenerate Source or Nodal Source(Unstable) 
depending of the geometric multiplicity of λ ,   Source (Degenerate, Nodal, Spiral Source)(Unstable)
depending of the geometric multiplicity of λ and the nonlinear terms }, {λ =μ<0, Degenerate Sink or Nodal Sink(Stable)
depending of the geometric multiplicity of λ,   Sink (Degenerate, Nodal, Spiral Sink)(Stable)
depending of the geometric multiplicity of λ and the nonlinear terms }, {λ ,μ   Complex Conjugate, □, □}, {Re(λ)>0, Spiral Source (Unstable), Spiral Source (Unstable)}, {Re(λ)<0, Spiral Sink (Asymptotically Stable), Spiral Sink (Asymptoticaly Stable)}, {Re(λ)=0, Center (Stable),                                  Center, Spiral Sink ,Spiral Source
(Stability cannot be determined based on the linear part alone !)
}})

Hyperbolic Equilibrium Points:

The equilibrium is said to be hyperbolic if all eigenvalues of the Jacobian matrix have non - zero real parts.

Hyperbolic equilibria are robust: Small perturbations do not change qualitatively the phase portrait near the equilibria. Moreover, local phase portrait of a hyperbolic equilibrium of a nonlinear system is equivalent to that of its linearization. This statement has a mathematically precise form known as the Hartman - Grobman Theorem.

Examples of nonlinear differential equations:

Example. Solve the system  ({{{, ({{x'= x(1-y), (1)}, {y'= y(x-1), (2)}})}})

Usually an explicit solution cannot be found, so DSolve will not work, but one can ask Mathematica to find a numeric solution using the NDSolve command.

DSolve[{x'[t]==x[t](1-y[t]), y'[t]==y[t](x[t]-1),x[0]==1,y[0]==1},{x[t],y[t]},t]

NDSolve[{x'[t]==x[t](1-y[t]), y'[t]==y[t](x[t]-1),x[0]==2,y[0]==1},{x[t],y[t]},{t,0,10}]

The most reliable information here comes from the Vector Field Plot.

StreamPlot[{x(1-y),y(x-1)},{x,-5,5},{y,-5,5},StreamScale->Automatic,StreamPoints->{{{{4,1},Red},Automatic}}]

What is the nature of the equilibrium solutions?

1. Set x' = 0 and y' = 0. In this case, it is easy to solve the equation by hand:  x(1-y)=0 and y(x-1)=0  implies (x,y)=(0,0) or (x,y)=(1,1). For more complicated cases, we will use Solve[..] or NSolve[..]

EqPoints=Solve[{x(1-y)==0, y(x-1)==0}, {x,y}]

2. Next we find the Jacobian matrix of the system. In this case, it is easy to compute it by hand: 
   (Subscript[J, (x,y)])=(({{(∂_x) (x(1-y)), (∂_y) (x(1-y))}, {(∂_x) (y(x-1)), (∂_y) (y(x-1))}}))=(({{1-y, -x}, {y, x-1}}))

In more complicated cases, we will ask Mathematica to compute the Jacobian matrix for us, as follows :

The commands D[f[x], {x}] or f'[x]  give the derivative of f with respect to x. If f is a function of several variables, D[f[x, y], {x}] or (∂_y)(f[x,y]) gives the partial derivative of f with respect to x. Try them out!

(∂_y)(x y+ y^2)
D[y^2, {y}]

The following code computes the Jacobian matrix:

Clear[f,J,x,y]
f[x_,y_]:= {x(1-y), y(x-1)}
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]]
J[f]

3. We evaluate the Jacobian matrix at the two equilibrium points and find its eigenvalues:

B=J[f] /. Flatten[EqPoints[[1]]];
Print["The Jacobian Matrix at the equilibrium point (0,0) is ", MatrixForm[B]]
Print["The eigenvalues of B are ", Eigenvalues[B]]

Clear[x,y];
StreamPlot[B.{x,y}, {x,-2,2},{y,-2,2}]

T=J[f] /. Flatten[EqPoints[[2]]];
Print["The Jacobian Matrix at the equilibrium point (1,1) is ",MatrixForm[ T]]
Print["The eigenvalues are ", Eigenvalues[T]]

Clear[x,y];
StreamPlot[T.{x,y}, {x,-2,2},{y,-2,2}]

4. Conclusions:

By analyzing the Real and Imaginary Part of the Eigenvalues, we conclude the following:
- the origin (0, 0) is a saddle fixed point (the Jacobian matrix has two real eigenvalues, one smaller than 0, one greater than 0)
- the point (1, 1) is a center (the Jacobian matrix has two complex conjugate eigenvalues with Real part equal to 0, and the Phase Portrait of the nonlinear equation depicts closed circular trajectories). 
This gives us information about the shape of trajectories in a vicinity of the equilibrium solutions.

Example 2:

Consider the system  ({{{, ({{x'= y-x(x^2+y^2), (1)}, {y'= -x-y(x^2+y^2), (2)}})}}).  Study the nature of the critical point (x,y)=(0,0).

Clear[h]
(* Define the function h *)
h[x_,y_]:={y-x(x^2+y^2),-x-y(x^2+y^2)}

(* Find the critical points of the system *)
sol = Solve[h[x,y]=={0,0},{x,y}];
Print["Critical point: ", sol]

(*Find the Jacobian matrix *)
Print["The Jacobian Matrix is ", MatrixForm[J[h]]]

(* Evaluate the Jacobian matrix at the only critical point (0,0) *)
A = J[h] /. Flatten[sol];
Print["The Jacobian Matrix at (0,0) is A=", MatrixForm[A]]
Print["The eigenvalues of A are: ", Eigenvalues[A]]

The only critical point for the nonlinear system is (0, 0). The eigenvalues of the Jacobian matrix at (0, 0) are i and - i (complex conjugate with Real part equal to 0), so the origin is a center for the linearized system. Let’s see if (0,0) is a center for the nonlinear system as well.

(* Phase Portrait for the linear system *)
StreamPlot[{y,-x }, {x,-2,2}, {y,-2, 2}, StreamPoints->{{{{1,0},Red},Automatic}}]

(* Phase Portrait for the nonlinear system *)
StreamPlot[{y-x(x^2+y^2),-x-y(x^2+y^2) }, {x,-2,2}, {y,-2, 2}, StreamPoints->{{{{1,0},Red},Automatic}}]

sol=NDSolve[{x'[t]==y[t]-x[t](x[t]^2+y[t]^2),y'[t]==-x[t]-y[t](x[t]^2+y[t]^2), x[0]==1,y[0]==0},{x[t],y[t]},{t,0,1000}]
ParametricPlot[{x[t],y[t]} /. Flatten[sol],{t,0,50}, PlotPoints->100]

ParametricPlot[{x[t],y[t]} /. Flatten[sol],{t,0,100}, PlotPoints->100]

DSolve[ ..] is not very helpful for nonlinear equations, but you can try it anyway:

DSolve[{x'[t]==y[t]-x[t](x[t]^2+y[t]^2),y'[t]==-x[t]-y[t](x[t]^2+y[t]^2), x[0]==1,y[0]==0},{x[t],y[t]},{t}]
DSolve[{x'[t]==y[t]-x[t]*(((x[t])^2)+((y[t])^2)),y'[t]==-x[t]-y[t] *(((x[t])^2)+((y[t])^2))},{x[t],y[t]},t]

Example 3

Consider the system  ({{{, ({{x'= y+x(x^2+y^2), (1)}, {y'= -x+y(x^2+y^2), (2)}})}}).  Study the nature of the critical point (x,y)=(0,0).

Clear[h]
h[x_,y_]:={y+x(x^2+y^2),-x+y(x^2+y^2)}
MatrixForm[J[h]]
Eigenvalues[J[h] /. {x->0, y->0}]

StreamPlot[{y+x(x^2+y^2),-x+y(x^2+y^2) }, {x,-2,2}, {y,-2, 2}]

Conclusions : The linearized system has a center at the origin. The nonlinear system does not have a center! The origin is a spiral source.

Exercise :

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!

Exercise (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.

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?