Phase Portraits of Nonlinear Differential Equations

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.

The critical points (equlibrium points, steady state solutions) are found by setting x' = y' = 0, 
and solving the resulting system f(x, y) = 0 and g(x, y) = 0.

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

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

We can use the eigenvalues of the matrix Jacobian matrix to decide the type and stability of the critical point (a,b) (sink, source, saddle, etc.).

Jacobian Matrix

The following code from last class helps compute the Jacobian matrix of a function f of two variables x and y:

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

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[{y-x(x^2+y^2)==0,-x-y(x^2+y^2)==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 ", MatrixForm[A]]
Print["Eigenvalues: ", 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}}]

(* Find the trajectory that passes through the point (1,0) *)
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,100}, PlotPoints->100]

(* Poincare-Bendixon criterion *)
(∂_x)(y-x(x^2+y^2))+(∂_y)(-x-y(x^2+y^2))

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

Example 3

Consider the system  ({{{, ({{x'= x+y-x(x^2+y^2), (1)}, {y'= -x+y-y(x^2+y^2), (2)}})}}).  Find all critical points and study their nature.

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

(* Find the critical points *)
NSolve[{y-x(x^2+y^2)==0,-x-y(x^2+y^2)==0},{x,y}]

MatrixForm[J[h]];
Eigenvalues[J[h] /. {x->0, y->0}]
StreamPlot[h[x,y], {x,-2,2}, {y,-2, 2}]

Conclusions : The origin is an spiral source (unstable spiral point) for both the linear system and the nonlinear system. One may therefore think that all trajectories different from 0 should spiral out to infinity. This is clearly not the case, as far away from the origin, the arrows are pointing inward. Moreover, there is a closed trajectory (a circle) and all trajectories different from (0,0) spiral towards this circle (limit cycle). You can find the equation of this limit cycle by switching to polar coordinates x=rcosθ, y=rsinθ.

Polar Coordinates

i. We do the coordinate change:

x = r Cos[θ] , y = r Sin [θ] , 0≤ r < ∞, 0≤ θ < 2π

where r and θ  are functions of time t.

ii. We need to compute x' and y'  as functions of r, θ, r’, θ’

a) We can do it by hand, using the Chain Rule. We get x’ = r’Cos[θ] - r Sin[θ]θ’  and  y’ = r’Sin[θ] + r Cos[θ]θ’
b) We use transformations rules in Mathematica to do the change of variables:

x'[t] /.  x->(r[#]Cos[θ[#]]&)
y'[t] /.  y ->(r[#]Sin[θ[#]]&)

Changing the system from Example 3 to polar coordinates

We know that (x^2)+(y^2)=(r^2) , so can use the symmetries of the system ({{{, ({{x'= x+y-x ((x^2)+(y^2)), (1)}, {y'= -x+y-y ((x^2)+(y^2)), (2)}})}})  to create more terms of the form (x^2)+(y^2) that can be easily “translated” into  (r^2). Therefore, we can look at the following system:
({{{, ({{xx'+yy' = x(x+y-x((x^2)+(y^2)))+y(-x+y-y((x^2)+(y^2)))=(x^2)+(y^2)-((x^2)+(y^2))((x^2)+(y^2))=(r^2)(1-(r^2)), (1)}, {yx'-xy'= y(x+y-x((x^2)+(y^2)))-x (-x+y-y((x^2)+(y^2)))=(x^2)+(y^2) =(r^2), (2)}})}})

Let' s see how the left hand sides xx’ + yy’ and yx’ - xy’ look like in polar coordinates:

rule = {x->(r[#]Cos[θ[#]]&),  y->(r[#]Sin[θ[#]]&) };
Simplify[x[t] x'[t]+ y[t] y'[t] /.  rule]
Simplify[y[t] x'[t] - x[t] y'[t] /. rule]

In polar coordinates, the system becomes:

({{{, ({{r r'=(r^2)(1-(r^2)), (1)}, {-(r^2)θ'=(r^2), (2)}})}})    which can be simplified to  ({{{, ({{r'=r(1-(r^2)), (1)}, {θ'=-1, (2)}})}})

Analyze the system written in polar coordinates :

The equation θ'=-1 , solves to 
θ(t)=-t + c, where c is any real constant.

The critical points of r'=r(1-(r^2)) are: 
r=0 (the origin of the old phase portrait) and 
r=1 (the unit circle in the old phase portrait). 

Therefore, a periodic solution of the system is 
r(t)=1, θ(t)=-t + c.  
(As time t increases, a point moves clockwise on the unit circle)

From the equation r'=r(1-(r^2)), we also see that:
r’ > 0 when r<1 (thus inside the unit circle the trajectories go outwards) 
r’ < 0 when r>1 (thus outside the unit circle the trajectories are directed inwards) 

What about the other trajectories? 
When r ≠ 1 and r ≠ 0, we can use separation of variables to solve the equation: (dr/dt)=r(1-(r^2)). 
We write (dr/(r(1-(r^2))))=dt and we can integrate each term separately.

Integrate[1/(r(1-r^2)),r]

The equation becomes log(r)-(1/2)log(1-(r^2)) = t + k, where k is any real constant, with solution given by r(t)=((e^(k+t))/(Sqrt[1+(e^(2 (k+t)))]))

solr= Simplify[Solve [Log[r]-(1/2)Log[1-(r^2)] ==t+k, r]]
Print["The function r is always positive, so the solution is: ", r /. solr[[2]]]

Poincare - Bendixon Theorem

The literature on the existence and properties of limits cycles is very rich: 

Poincare - Bendixon Theorem 
Suppose f and g have continuous partial derivatives. 
a) Every closed trajectory must enclose at least one critical point. If it encloses exactly one critical point, then this cannot be a saddle point.
b) If (∂_x)f+ (∂_y)g has the same sign in a disk D, then there can be no closed trajectory of the system lying entirely inside D.
c) Suppose that U is a connected open subset that contains no critical point of the system. Let (OverBar[U]) denote the set U together with its boundary. If the vector field in (OverBar[U]) points towards the interior of U, then U contains a closed trajectory.
d) Suppose that U is a connected open subset that contains no critical point of the system. Let (OverBar[U]) denote the set U together with its boundary. If there exists some trajectory that never leaves the set U (is trapped inside U),  then U must necessarily contain a closed trajectory.
e) Every trajectory is attracted by a fixed point, or it goes to infinity, or it accumulates on a limit cycle. 

This theorem has no analogue in higher dimensions. In systems with phase space of dimension strictly greater than 2, more things can happen: trajectories may be attracted to complex geometrical objects (strange attractors), which leads to chaotical behavior.

Nonlinear Equations with several variables :

In 1963, Lorenz studied a very simple model of the atmosphere. This model exhibits chaotic behavior:

Lorenz equation

tmax=20;
sol=NDSolve[{x'[t]==-10 (x[t]-y[t]),y'[t]==-x[t] z[t]+28 x[t]-y[t],z'[t]==x[t] y[t]-(8/3) z[t],x[0]==30,y[0]==10,z[0]==40},{x,y,z},{t,0,tmax},MaxSteps->5000];

The 3D vector field can be plotte using the command VectorPlot3D

VectorPlot3D[{-10 (x-y),-x *z+28 x-y,x* y-(8/3) z}, {x, -10,10}, {y, -10,10},{z,-10,10}, VectorScale->{Small,Automatic,None}]

The solution can the be plotted using ParametricPlot3D as follows:

Clear[x,y,z,t]
ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.sol],{t,0,tmax},PlotPoints->3000,Axes->False,PlotRange->All]

Clear[x,y,z,t]
Manipulate[
ParametricPlot3D[
Evaluate[{x[t],y[t],z[t]}/.NDSolve[{x'[t]==-10 (x[t]-y[t]),y'[t]==-x[t] z[t]+28 x[t]-y[t],z'[t]==x[t] y[t]-(8/3) z[t],x[0]==30,y[0]==10,z[0]==40},{x,y,z},{t,0,T},MaxSteps->5000]],{t,0,T},PlotStyle->Red],
{{T,2},0.1, 20},SynchronousUpdating->False,
SaveDefinitions->True]