The Phase Plane Analysis and Stability
Prashant K Srivastava, Ph.D., MNASc
Associate Professor,
Department of Mathematics, IIT Patna
pksri@iitp.ac.in
P.K.Srivastava (IIT Patna)
General First Order Differential Equation
Scaler equation:
Let 𝑥(𝑡) ∈ 𝑅 be an unknown function of 𝑡 and 𝑓: 𝑅 × 𝑅 → 𝑅 be a continuous function with
continuous partial derivatives w.r.t 𝑥 on some closed interval 𝐼, then
𝒅𝒙
𝒙′ = = 𝒇(𝒕, 𝒙)
𝒅𝒕
represents a first order differential equation.
If function 𝑓 is independent of 𝑡 then the equation is known as autonomous equation.
If the function 𝑓 is linear in 𝑥 then we call the equation as linear, otherwise nonlinear.
The solution of linear equation 𝑥 ′ = 𝑎𝑥, 𝑥(𝑡0 ) = 𝑥0 can be easily verified as 𝑥 (𝑡) =
𝑥0 𝑒 𝑎(𝑡−𝑡0) .
The differential equation along with the initial condition 𝑥(𝑡0 ) = 𝑥0 has a unique solution in
a neighbourhood of 𝑡0 . There is no general method to find a solution, but qualitative
methods help to understand pattern of solutions.
P.K.Srivastava (IIT Patna)
System of equations:
Consider the following system of differential equations,
𝑑𝒙
= 𝒇(𝑡, 𝒙),
𝑑𝑡
where 𝒙(𝑡) = (𝑥1 , 𝑥2 , … , 𝑥𝑛 )𝑇 ∈ 𝑅 𝑛 and 𝒇 = (𝑓1 , 𝑓2 , … , 𝑓𝑛 )𝑇 . Suppose 𝒇 =
(𝑓1 , 𝑓2 , … , 𝑓𝑛 )𝑇 satisfies existence and uniqueness conditions i.e. there exists a unique
solution for the system with some initial conditions.
As discusses earlier, if 𝒇 is independent of 𝑡, system is autonomous.
Similarly linearity of 𝒇 w.r.t. 𝒙 will give a linear system, otherwise nonlinear system.
In this talk we shall consider only autonomous equations i.e.
𝑑𝒙
= 𝒇(𝒙)
𝑑𝑡
We shall use 𝑥 for scaler case and 𝒙 for vector case (system).
P.K.Srivastava (IIT Patna)
Example:
(A) Scaler equation : 𝑥 ′ = 𝑠𝑖𝑛𝑥
(B) System of equations:
𝑥1′ = 𝑥1 + 𝑥1 𝑥2 = 𝑓1 (𝑥1 , 𝑥2 )
𝑥2′ = 𝑥12 + 𝑥2 = 𝑓2 (𝑥1 , 𝑥2 ).
or equivalently
𝒙′ = 𝒇(𝒙).
Here 𝒙 = (𝑥1 , 𝑥2 )𝑇 and 𝒇 = (𝑓1 , 𝑓2 )𝑇 .
Note: A higher order differential equation can be written as system of first order differential
equation. For example, the equation: 𝑥 ′′ + 𝑡𝑥 ′ + 𝑠𝑖𝑛𝑥 = 0, by using transformation, 𝑥 = 𝑥1 ,
𝑥 ′ = 𝑥2 , can be expressed as 𝒙′ = 𝒇(𝒙), where, 𝒙 = (𝑥1 , 𝑥2 )𝑇 and 𝒇 = (𝑥2 , −𝑠𝑖𝑛𝑥1 − 𝑡𝑥2 )𝑇 .
P.K.Srivastava (IIT Patna)
Existence and Uniqueness Theorem:
P.K.Srivastava (IIT Patna)
Important note: Although existence and uniqueness theorem only gives
information about the existence of solution of differential equation, however, the
idea of proof can be used to construct the solution. Let us consider following
example:
Example 1. 𝒙′ = 𝒕 + 𝒙, 𝒙(𝟎) = 𝟏.
From the differential equation, we get Picard’s iterates as,
𝑡
𝑥𝑛 = 𝑥 (0) + ∫ 𝑓 (𝑡, 𝑥𝑛−1 )𝑑𝑡.
0
So that we get
𝑡2 𝑡 3
𝑥0 = 1, 𝑥1 = 1 + 𝑡 + , 𝑥2 = 1 + 𝑡 + 𝑡 2 + ,
2 6
3 4
𝑡 𝑡
𝑥4 = 1 + 𝑡 + 𝑡 2 + +
3 24
P.K.Srivastava (IIT Patna)
3 𝑛
𝑡 𝑡
𝑥𝑛 = 1 + 𝑡 + 𝑡 2 + + ⋯ + .
3 𝑛!
As 𝑛 → ∞, it can be easily seen that solution becomes: 𝒙 = −𝟏 − 𝒕 + 𝟐𝒆𝒕 .
The above process helps in constructing solutions but it has a drawback. The
integrand, if it is complicated function, it may lead to tedious calculation and
even for simpler cases closed form solution might be difficult.
However, this is an important tool for constructing solution of an IVP.
The above example can also be solved simply considering integrating factor
and integration as it is a linear differential equation.
Let us now consider following example:
Example 2. 𝒙′ = 𝒕𝟐 + 𝒙𝟐 , 𝒙(𝟎) = 𝟏.
P.K.Srivastava (IIT Patna)
Geometric Method (Qualitative solution):
Consider the following differential equation 𝒙′ = 𝒇(𝒕, 𝒙)
Notice that in two dimensional (t,x) plane 𝑥′ is nothing but
slope at the point (t,x) which is given by the function value
𝑓(𝑡, 𝑥 ). Hence in general we can evaluate slope at each
point in the two dimensional plane. Choosing a starting
point, depending on the slope, we can draw a solution
curve. This curve will be solution because at any point, the
tangential vector is same as 𝑓value.
There is a systematic way to do so.
1. Find all isoclines (curves with same inclination), 𝑓(𝑡, 𝑥 ) = 𝑐,
2. Measure the slope on each of these curves, 𝑡𝑎𝑛𝜃 = 𝑐.
3. Draw a curve from any starting point and follow slope/ gradient.
P.K.Srivastava (IIT Patna)
Example 1. 𝒙′ = 𝒕 + 𝒙, 𝒙(𝟎) = 𝟏.
For solution, the isoclines are given by
𝑡+𝑥 =𝑐
These are straight lines, depending on
value of 𝑐 ∈ 𝑅 we get parallel straight lines
with slope 450 .
Mark the slope as 𝑡𝑎𝑛𝜃 = 𝑐 on each line.
Here we have chosen 0, ±1, ±2, ±3 for
plotting isoclines.
Draw solution curves from x(0)=1 as initial
point.
P.K.Srivastava (IIT Patna)
Example 2. 𝒙′ = 𝒕𝟐 + 𝒙𝟐 , 𝒙(𝟎) = 𝟏.
For solution, the isoclines are given by
𝑡2 + 𝑥2 = 𝑐2
These are circles, depending on value of c
from 0 to a large number we get concentric
circles.
Mark the slope as 𝑡𝑎𝑛𝜃 = 𝑐 2
Draw solution curves from any initial point.
Reference: S.L. Ross, Differential equations/ GF
simmons, Differential equations
P.K.Srivastava (IIT Patna)
Solution of system of Linear Ordinary Differential
Equations:
We shall analyse for system of two linear equations. Same technique can be
applied to find solution of 𝑛 differential equations. Consider as an example:
𝒅𝒙
= 𝒂𝒙 + 𝒃𝒚,
{ 𝒅𝒕
𝒅𝒚
= 𝒄𝒙 + 𝒅𝒚.
𝒅𝒕
Assume that (𝑦𝑥 ) = (𝐵𝐴)𝑒 𝜆𝑡 is a solution of above equations.
P.K.Srivastava (IIT Patna)
Hence on substitution, we obtain:
𝑎−𝜆 𝑏 𝐴 𝜆𝑡 0
( )( )𝑒 = ( )
𝑐 𝑑−𝜆 𝐵 0
Which gives following system of linear algebraic equations:
𝑎−𝜆 𝑏 𝐴 0
( )( ) = ( )
𝑐 𝑑−𝜆 𝐵 0
This system will have non-trivial solution (note that we need 𝐴 ≠ 0, 𝐵 ≠ 0,)
provided
𝑎−𝜆 𝑏
| |=0
𝑐 𝑑−𝜆
Which gives,
𝜆2 − (𝑎 + 𝑑 )𝜆 + (𝑎𝑑 − 𝑏𝑐 ) = 0.
P.K.Srivastava (IIT Patna)
This equation has possibility of
1. Two distinct real roots
2. Complex roots
3. Real and repeated roots.
We shall analyse all these cases, with help of examples.
Case-I: Two distinct real roots
Consider the system of equations:
𝒅𝒙
= −𝒙 + 𝟐𝒚,
𝒅𝒕
𝒅𝒚
= −𝟑𝒚.
𝒅𝒕
Substituting (𝑦𝑥 ) = (𝐵𝐴)𝑒 𝜆𝑡 , we obtain the following characteristic equation of
system
P.K.Srivastava (IIT Patna)
𝜆2 + 4𝜆 + 3 = 0
Corresponding eigen values are: 𝜆1,2 = −1, −3(both are real and negative).
The eigen vector corresponding to 𝜆 = −1 is given by (10)
1
eigen vector corresponding to 𝜆 = −3 is given by (−1 ).
Therefore, the general solution is given as:
𝑥 1 −𝑡 1
( ) = 𝑐1 ( ) 𝑒 + 𝑐2 ( ) 𝑒 −3𝑡 .
𝑦 0 −1
In the following we shall provide detailed graphing of these solutions.
P.K.Srivastava (IIT Patna)
The idea here is to see how the solutions obtained above are represented in
two dimensional plane. For this we shall plot specific trajectories and then try
and fill the phase plane.
Note when 𝑐1 = 0, 𝑥 = 𝑐2 𝑒 −3𝑡 and 𝑦 =
−𝑐2 𝑒 −3𝑡 giving 𝑦 = −𝑥 corresponding
to 𝑐2 = 1. Notice in this case point
corresponding to 𝑡 = 0 is 𝑥 = 1,𝑦 =
−1 which lies in 4th quadrant on line
𝑦 = −𝑥. Further as 𝑡 → ∞ both 𝑥, 𝑦 →
0 and 𝑡 → −∞ 𝑥 → ∞ and 𝑦 → −∞.
Hence the trajectory in this case is
represented as straight line 𝑦 = −𝑥 with direction towards origin. (Trajectory-
1)
P.K.Srivastava (IIT Patna)
Similarly for the case when 𝑐2 = −1 the entire process is same and the
trajectory is part of straight line 𝑦 =
−𝑥 lying in 2nd quadrant and direction
is again towards origin as both 𝑥, 𝑦 →
0 for 𝑡 → ∞. (Trajectory-2)
Again, notice that for 𝑐2 = 0 we
have𝑥 = 𝑐1 𝑒 −𝑡 and 𝑦 = 0 giving 𝑦 =
0 corresponding to 𝑐1 = 1. Notice
again in this case point corresponding
to 𝑡 = 0 is 𝑥 = 1, 𝑦 = 0 which lies on
ight half of x-axis. Further as 𝑡 → ∞ both 𝑥, 𝑦 → 0 and 𝑡 → −∞ 𝑔𝑖𝑣𝑒𝑠 𝑥 → ∞
and 𝑦 = 0. Hence the trajectory in this case is represented as straight line 𝑦 =
0 with direction towards origin. (Trajectory-3)
P.K.Srivastava (IIT Patna)
Similarly for the case when 𝑐1 = −1
the entire process is same and the
trajectory is part of left half of x-axis
and direction is again towards origin
as both 𝑥, 𝑦 → 0for 𝑡 → ∞.
(Trajectory-4)
Further, we shall plot trajectories
corresponding to both 𝑐1 , 𝑐2 ≠ 0. In
this case we notice from general solution that for large 𝑡,𝑒 −𝑡 dominates and in
case when 𝑡 is small 𝑒 −3𝑡 dominates. Also, 𝑥, 𝑦 → 0 for 𝑡 → ∞. Thus in this
case we notice that trajectories will be initiating parallel to line𝑦 = −𝑥 as
𝑡starts from −∞and approach the direction parallel to x-axis as 𝑡 becomes
P.K.Srivastava (IIT Patna)
larger before finally entering the origin as 𝑡 → ∞.(Trajectories in green in
attached fig)
P.K.Srivastava (IIT Patna)
Case-II: Complex roots
Consider the system of equations:
𝒅𝒙
= 𝟒𝒙 − 𝟑𝒚,
𝒅𝒕
𝒅𝒚
= 𝟑𝒙 + 𝟒𝒚.
𝒅𝒕
Substituting (𝑦𝑥 ) = (𝐵𝐴)𝑒 𝜆𝑡 , we obtain the following characteristic equation of
system
𝜆2 − 8𝜆 + 25 = 0
and corresponding eigen values are: 𝜆1,2 = 4 ± 3𝑖(Complex conjugate).
The eigen vector corresponding to 𝜆 = 4 − 3𝑖 is given by(1𝑖).
Therefore the general solution is given as:
P.K.Srivastava (IIT Patna)
𝑥 cos 3𝑡 4𝑡 − sin 3𝑡 4𝑡
( ) = 𝑐1 ( ) 𝑒 + 𝑐2 ( )𝑒 .
𝑦 sin 3𝑡 cos 3𝑡
The phase plane for this system is as given in Figure. Note that each solution
spirals counter clockwise and emanates from
origin.
This can be seen by a solution at two different
𝜋
points, say at 𝑡 = 0, and 𝑡 = .
6
This can be easily seen from presence of sine
and cosine terms and exponential factor. Also,
as 𝑡 → −∞ solutions enter origin and move
away as 𝑡 increases.
P.K.Srivastava (IIT Patna)
Case-III: Purely Imaginary roots
Consider the system of equations:
𝒅𝒙
= 𝒚,
𝒅𝒕
𝒅𝒚
= −𝟒𝒙.
𝒅𝒕
Substituting (𝑦𝑥 ) = (𝐵𝐴)𝑒 𝜆𝑡 , we obtain the following characteristic equation of
system
𝜆2 + 4 = 0
and corresponding eigen values are: 𝜆1,2 = ±2𝑖(purely imaginary roots).
1
The eigen vector corresponding to 𝜆 = 2𝑖 is given by (2𝑖 ).
The solution is given as: (𝑦𝑥 ) = (2𝑖
1 cos 2𝑡
)𝑒 2𝑖 = (−2sin 2𝑡
) + 𝑖( sin 2𝑡
2cos 2𝑡
).
P.K.Srivastava (IIT Patna)
Therefore, the general solution is given as: (𝑦𝑥 ) = 𝑐1 (−2sin
cos 2𝑡
2𝑡
) + 𝑐 sin 2𝑡
2 2cos 2𝑡).
(
Note that from solutions, we get,
𝑥2 𝑦2
2 2+ ( 2 2 ) = 1,
𝑐1 + 𝑐2 4 𝑐1 + 𝑐2
for all constants 𝑐1 , 𝑐2 . These are ellipses with y as major axis, i.e. closed orbits
around the origin.
For direction, we look at a solution at two different points, say at 𝑡 = 0, and
𝜋
𝑡 = and notice that solutions are moving clockwise.
4
P.K.Srivastava (IIT Patna)
Case-IV: Repeated roots
Consider the system of equations:
𝒅𝒙
= 𝒙 − 𝟑𝒚,
𝒅𝒕
𝒅𝒚
= 𝟑𝒙 + 𝟕𝒚.
𝒅𝒕
Substituting (𝑦𝑥 ) = (𝐵𝐴)𝑒 𝜆𝑡 , we obtain the following characteristic equation of
system
(𝜆 − 4)2 = 𝜆2 + 16𝜆 + 16 = 0
and corresponding eigen values are: 𝜆1,2 = 4,4(repeated roots).
1
The eigen vector corresponding to 𝜆 = 4 is given by(−1 ) and hence we have
just one eigen vector. We need another vector, and in this case it is given by
P.K.Srivastava (IIT Patna)
generalized eigen vector given by (−1/3
0
). Therefore the general solution is
given as:
𝑥 1 4𝑡
𝑡 − 1/3 4𝑡
( ) = 𝑐1 ( ) 𝑒 + 𝑐2 ( )𝑒 .
𝑦 −1 −𝑡
With 𝑐2 = 0 we note
that the curve 𝒙 =
−𝒚 represents the
solutions. Further
solution trajectories
move parallel to this
curve curve 𝒙 = −𝒚 as
𝑡 increases.
P.K.Srivastava (IIT Patna)
Note:
1. With above analysis, we can understand the dynamics of system of linear
equations in two dimensions. The same analysis can be extended to higher
dimensions.
𝑑𝒙
2. The vector linear equation = 𝐴𝒙 has solution of the form 𝒙 = 𝒄𝑒 𝜆𝑡 .
𝑑𝑡
General solution of this linearized system is given by
𝑟 𝑛−𝑟
𝒙 = ∑ 𝑨𝒊 𝑡 𝑖−1 𝑒 𝜆1 𝑡 + ∑ 𝑨𝒓+𝒊 𝑒 𝜆𝑖𝑡
𝑖=1 𝑖=2
Where 𝜆𝑖 are eigen values of the matrix A and 𝜆1 is repeated r times.
3. In general, the eigen values of linear system can tell us about stability of
system.
P.K.Srivastava (IIT Patna)
Critical Points:
A point 𝑥 = 𝑥 ∗ is said to be critical point (also known as equilibrium point
𝑑𝑥
or steady state), if | = 0.
𝑑𝑡 𝑥=𝑥 ∗
Note: The same definition translates to system of differential equations.
Example: 1) 𝑥 ′ = sin 𝑥, critical points are 𝑥 = 𝑛𝜋, 𝑛 ∈ ℤ.
2) (𝑥1′ , 𝑥2′ )′ = (𝑥2 , − sin 𝑥1 )′, has (𝑛𝜋, 0), as critical points.
Important: Isolated critical points.
P.K.Srivastava (IIT Patna)
Stability of critical points:
𝑑𝒙
A critical point 𝒙∗ of the autonomous system = 𝒇(𝒙) is said to be stable if,
𝑑𝑡
given any 𝜖 > 0, there exists a 𝛿 > 0 such that every solution 𝒙 = 𝝓(𝑡) of the
system, which at 𝑡 = 𝑡0 satisfies, ||𝝓(𝑡0 ) − 𝒙∗ || < 𝛿, exists for all 𝑡 ≥ 𝑡0 and
satisfies ||𝝓(𝑡) − 𝒙∗ || < 𝜖.
∗ 𝑑𝒙
A critical point 𝒙 of the autonomous system = 𝒇(𝒙) is said to be
𝑑𝑡
asymptotically stable if, it is stable and in addition 𝐿𝑖𝑚𝑡→∞ 𝝓(𝑡) = 𝒙∗
A critical point, that is not stable is unstable.
Note: The 𝛿 in the definition is always less than or equal to 𝜖.
P.K.Srivastava (IIT Patna)
Note: Thus, in case of stability, solutions starting near (in 𝛿 -nbd) the critical point
remain near critical point (in 𝜖 -nbd). Whereas for the asymptotic stability, solutions not
only remain close but eventually enters the critical point. An interpretation in two
dimensions is given in Figure.
Figure: Illustration of Stability in Two Dimensions
P.K.Srivastava (IIT Patna)
Physical interpretation of stability:
The mathematical definition reflects the sustainability of systems for
perturbations!
If we give a small perturbation, from critical point, and solution trajectory
either return back or remain close to critical point, then the critical point is
stable.
Take example of a ball kept on top of a hill or in bottom of a valley (see
figure). If we perturb the ball on hill, very little, it will move away and never
come back. If we perturb ball in valley, it will move in valley and eventually
reach to its original position.
For system of nonlinear equations, it is difficult to establish stability of critical
points using definition. The two important tools used are Lyapunov function
method and method of linearization.
P.K.Srivastava (IIT Patna)
Stability Analysis for scaler case (Linearization):
For Linear Stability analysis we consider the autonomous first order differential equation
𝑥 ′ = 𝑓 (𝑥 )
In order to describe the behaviour of solutions near equilibrium we introduce the
process of Linearization. If 𝑥 ∗ is an equilibrium of the differential equation 𝑥 ′ = 𝑓(𝑥 ) so
that 𝑓 (𝑥 ∗ ) = 0,we make the change of variable 𝑢(𝑡) = 𝑥(𝑡) − 𝑥 ∗ ,representing
deviation of the solution from the equilibrium value. Substitution gives
𝑢′ (𝑡) = 𝑓(𝑥 ∗ + 𝑢(𝑡))
An application of Taylor’s Theorem gives
𝑓 ′′ (𝑐 )
𝑢′ (𝑡) = 𝑓 (𝑥 ∗ ) + 𝑓 ′ (𝑥 ∗ )𝑢(𝑡) + 𝑢2 (𝑡),
2!
for some 𝑥 ∗ < 𝑐 < 𝑥 ∗ + u(t).
P.K.Srivastava (IIT Patna)
𝑓′′ (𝑐) 2
We use 𝑓 (𝑥 ∗ ) = 0 and write ℎ(𝑢) = 𝑢 (𝑡). So, we may rewrite the differential
2!
equation 𝑥 ′ = 𝑓(𝑥)in the equivalent form .
𝑢′ (𝑡) = 𝑓 ′ (𝑥 ∗ )𝑢(𝑡) + ℎ(𝑢).
ℎ(𝑢)
The function ℎ(𝑢) is small for small |𝑢| in the sense that → 0 as 𝑢 → 0.
𝑢
More precisely ,for every 𝜖 > 0 there exists 𝛿 > 0 such that |ℎ(𝑢)| < 𝜖|𝑢| whenever
|𝑢| < 𝛿.The linearization of the differential equation at the equilibrium 𝑥 ∗ is defined to
be the linear homogeneous differential equation.
𝑢′ (𝑡) = 𝑓 ′ (𝑥 ∗ )𝑢(𝑡)
obtained by neglecting the higher order term ℎ(𝑢).
The importance of the linearization lies in the fact that the behaviour of its solutions is
easy to analyse, and its behaviour also describes the behaviour of solutions of the
original equation near the equilibrium.
P.K.Srivastava (IIT Patna)
If all solutions of the linearization of 𝒙′ = 𝒇(𝒙) at an equilibrium 𝑥 ∗ tend to zero as 𝑡 →
∞ then all solutions of equation 𝒙′ = 𝒇(𝒙) with 𝑥(0)sufficiently close to 𝑥 ∗ tend to the
equilibrium 𝑥 ∗ as 𝑡 → ∞.
Theorem: An equilibrium 𝒙∗ of equation 𝒙′ = 𝒇(𝒙) with 𝒇′ (𝒙∗ ) < 𝟎 is asymptotically
stable, while an equilibrium 𝒙∗ with
𝒇′ (𝒙∗ ) > 𝟎 is unstable.
Note: For nonlinear scaler equation 𝒙′ =
𝒇(𝒙) it is very easy to judge the stability
if one observes the graph of the function
𝒇(𝒙). At the points where this function
intersects the 𝑥-axis, they are the critical
points and if the graph is below 𝑥-axis, we mark arrows to the left and if the graph is above
𝑥-axis, we mark arrows to the right. If at a critical point, both side arrows are inward, then
it is asymptotically stable, otherwise unstable. Example: 𝒙′ = 𝒓𝒙(𝒙 − 𝟏)(𝟐 − 𝒙).
P.K.Srivastava (IIT Patna)
Linear Stability for System of equations:
Theorem for Linear System: The critical point 𝒙 = 𝟎 of the linear system
𝑑𝒙
= 𝐴𝒙
𝑑𝑡
is
asymptotically stable if the eigenvalues of matrix 𝑨 have negative
real part;
stable, but not asymptotically stable, if (a zero Or a pair of pure
imaginary) and rest all with negative real part;
unstable, if at least one of them is with positive real part.
P.K.Srivastava (IIT Patna)
Routh Hurwitz criterion:
Theorem: A necessary and sufficient condition for the negativity of real
parts of the roots of a polynomial
𝐿(𝜆) = 𝜆𝑛 + 𝑎1 𝜆𝑛−1 +𝑎2 𝜆𝑛−2 + ⋯ + 𝑎𝑛−1 𝜆 + 𝑎𝑛 = 0,
with real coefficients is the positivity of all the principle diagonal minors of
the Hurwitz matrix:
𝑎1 1 0 0 0 0
𝑎3 𝑎2 𝑎1 ⋯ 0 0 0
𝐻𝑛 = 𝑎5 𝑎4 𝑎3 0 0 0
⋮ ⋱ ⋮
[ 0 0 0 ⋯ 0 0 𝑎𝑛 ]
𝑎1 1
The principle diagonal minors are given as, 𝐷1 = 𝑎1 , 𝐷2 = | | , . . 𝐷𝑛 = 𝐷𝑒𝑡(𝐻𝑛 ).
𝑎3 𝑎2
P.K.Srivastava (IIT Patna)
Example 1. 𝜆2 + 𝑎1 𝜆 + 𝑎2 = 0. RH conditions- 𝑎1 > 0, 𝑎2 > 0.
Example 2. 𝜆3 + 𝑎1 𝜆2 + 𝑎2 𝜆 + 𝑎3 = 0.
RH conditions- 𝑎1 > 0, 𝑎2 > 0, 𝑎3 > 0, 𝑎1 𝑎2 − 𝑎3 > 0.
Example 3. 𝜆3 + 5𝜆2 + 9𝜆 + 5 = 0.
RH conditions- 𝑎1 𝑎2 − 𝑎3 = 5 × 9 − 5 = 40 > 0. So all roots will be
with negative real parts.
Note: 𝑎1 > 0, 𝑎2 > 0, … 𝑎𝑛 > 0 are necessary conditions, whereas the
positivity of determinants are sufficient conditions (which also imply
necessary conditions). Thus, for n-th degree polynomial, we get n conditions,
referred as RH conditions or Routh-Hurwitz conditions.
P.K.Srivastava (IIT Patna)
Linearization of System of First Order Ordinary
Differential Equations:
Consider the following system of differential equations,
𝑑𝒙
= 𝒇(𝒙)
𝑑𝑡
Linearization around a critical point 𝒙∗ :
𝒇(𝒙 + 𝒙∗ ) = 𝑨𝒙 + 𝒉. 𝒐. 𝒕.
Here
𝜕(𝑓1 , 𝑓2 , … , 𝑓𝑛 )
𝐴=
𝜕(𝑥1 , 𝑥2 , … , 𝑥𝑛 )
is Jacobian Matrix evaluated at x*.
So that the system of equations is transformed into, taking (new 𝑥 = 𝑜𝑙𝑑 𝑥 − 𝑥 ∗ )
P.K.Srivastava (IIT Patna)
Linearization of a two-dimensional nonlinear system of
equations:
𝑑(𝑥1 +𝑥1∗ )
= 𝑓1 (𝑥1 + 𝑥1∗ , 𝑥2 + 𝑥2∗ )
Example: x old---> x new+x* {𝑑(𝑥𝑑𝑡+𝑥 ∗ )
2 2
= 𝑓2 (𝑥1 + 𝑥1∗ , 𝑥2 + 𝑥2∗ )
𝑑𝑡
𝑑𝑥1 𝜕𝑓1 𝜕𝑓1
= | 𝑥1 + | 𝑥2 + ℎ. 𝑜. 𝑡
𝑑𝑡 𝜕𝑥1 ∗ 𝜕𝑥2 ∗
𝑑𝑥2 𝜕𝑓2 𝜕𝑓2
= | 𝑥1 + | 𝑥2 + ℎ. 𝑜. 𝑡
{ 𝑑𝑡 𝜕𝑥1 ∗ 𝜕𝑥2 ∗
𝜕𝑓1 𝜕𝑓1
𝑑𝑥1
𝜕𝑥1 𝜕𝑥2
( 𝑑𝑡
𝑑𝑥2 ) = 𝐴 (𝑥𝑥1 ) where 𝐴 = 𝜕𝑓2 𝜕𝑓2
|
2
𝑑𝑡
𝜕𝑥1 𝜕𝑥2 (𝑥∗ ,𝑥∗ )
1 2
Notice that this linearized system has (0,0) as critical point.
P.K.Srivastava (IIT Patna)
Hyperbolic critical point:
The critical point is called hyperbolic critical point if none of the eigen values of
linearized matrix have zero real parts.
Hartman Grobman Theorem:
This theorem tells that near a “hyperbolic critical point” the linearized system
𝒙′ = 𝑨𝒙 has same topological structure as the non-linear system 𝒙′ = 𝒇(𝒙).
Therefore, this theorem is indispensable tool for analysing nonlinear systems. We
perform phase plain analysis for linearized system and conclude that near the
critical point (if it is hyperbolic), the structure will be same for nonlinear system.
In the following, we discuss the entire process with help of examples.
P.K.Srivastava (IIT Patna)
Qualitative solutions of system of nonlinear differential
equations
Example (A competition model):
𝒅𝒙𝟏 𝟏 𝟐
= 𝟏𝟒𝒙𝟏 − 𝒙𝟏 − 𝒙𝟏 𝒙𝟐
𝒅𝒕 𝟐
𝒅𝒙𝟐 𝟏 𝟐
= 𝟏𝟔𝒙𝟐 − 𝒙𝟐 − 𝒙𝟏 𝒙𝟐
𝒅𝒕 𝟐
For critical points, we solve
1 2
14𝑥1 − 𝑥1 − 𝑥1 𝑥2 = 0,
2
1 2
16𝑥2 − 𝑥2 − 𝑥1 𝑥2 = 0.
2
Which gives : (𝟎, 𝟎), (𝟎, 𝟑𝟐), (𝟐𝟖, 𝟎), (𝟏𝟐, 𝟖) as critical points.
P.K.Srivastava (IIT Patna)
First we find Jacobian matrix:
𝟏𝟒 − 𝒙𝟏 − 𝒙𝟐 −𝒙𝟏
𝑱=( )
−𝒙𝟐 𝟏𝟔 − 𝒙𝟏 − 𝒙𝟐
Linear Analysis for critical point (0,0):
The Jacobian matrix 𝐽 at (0,0), i.e. 𝐽(0,0) =
14 0
( )
0 16
has eigen values 𝜆1 = 14 with eigen vector
[1 0]𝑇 and 𝜆2 = 16 with eigen vector [0 1]𝑇
Thus (0,0) is unstable and is as depicted in
Figure (locally) the solutions move away from
origin.
P.K.Srivastava (IIT Patna)
Linear Analysis for critical point (0,32):
The Jacobian matrix 𝐽 at (0,32), i.e.
−18 0
𝐽(0,32) = ( )
−32 −16
has eigen values 𝜆1 = −18 with eigen
vector [1 16]𝑇 and 𝜆2 = −16 with eigen
vector [0 1]𝑇
Thus (0,32) is stable (as eigen values are
negative) and is as depicted in Figure
(locally) the solutions move to origin
(which is (0,32) for linearized system).
P.K.Srivastava (IIT Patna)
Linear Analysis for critical point (28,0):
−14 −28
The Jacobian matrix 𝐽 at (28,0), i.e. 𝐽(28,0) = ( )
0 −12
has eigen values 𝜆1 = −14 with eigen vector [1 0]𝑇 and 𝜆2 = −12 with eigen
vector [−14 1]𝑇
Thus (28,0) is stable (as eigen
values are negative) and is as
depicted in Figure (locally) the
solutions move to origin (which
is (28,0) for linearized system).
P.K.Srivastava (IIT Patna)
Linear Analysis for critical point (12,8):
−6 −12
The Jacobian matrix 𝐽 at (12,8), i.e. 𝐽(12,8) =( ) has characteristic
−8 −4
equation as 𝜆2 + 10𝜆 − 72 = 0. Thus we
have eigen values 𝜆1 = −5 − √97 < 0
𝑇
1+√97
with eigen vector [ 1] and
8
𝜆2 = −5 + √97 > 0 with eigen vector
𝑇
1−√97
[ 1]
8
Thus (12,8) is a saddle (eigen values are of
opposite signs) and unstable and is as
depicted in Figure (locally).
P.K.Srivastava (IIT Patna)
The complete global picture can be obtained by assembling all these local
pictures and put them together.
Then smoothly join the flow so as to
complete the phase plane diagram
as given in Figure. Notice that in this
case the region is divided into two
separate sets where trajectories
initiating from different regions tend
to move to respective critical points.
This shows that, for a competing
species, a species may outnumber
other species and survive.
The separating curve of two regions is known as separatrix. So for this model the
coexistence is not possible.
P.K.Srivastava (IIT Patna)
Now, consider another example:
𝒅𝒙𝟏
𝟐
𝒅𝒕 = 𝟏𝟒𝒙𝟏 − 𝟐𝒙𝟏 − 𝒙𝟏 𝒙𝟐
𝒅𝒙𝟐 = 𝟏𝟔𝒙𝟐 − 𝟐𝒙𝟐 𝟐 − 𝒙𝟏 𝒙𝟐
𝒅𝒕
Proceeding same as above, we observe that the solution trajectories may be
given as in Figure. We note that in this case the coexisting critical point (4,6) is
stable as both eigen values are negative. This shows that, in this case, the both
species peacefully coexist.
P.K.Srivastava (IIT Patna)
P.K.Srivastava (IIT Patna)
Lyapunov Stability
Suppose that the autonomous system has an isolated critical point at 𝒙∗ . If there
exists a function 𝑉such that,
(i) 𝑉 (𝒙∗ ) = 0,
𝑑𝑉
& for 𝑥 ≠ 𝑥 ∗ (ii) 𝑉 (𝑥) positive definite, and (iii) is negative definite
𝑑𝑡
(on some domain D in ℝ𝑛 containing 𝒙∗ )
A. then 𝒙∗ is an asymptotically stable critical point.
𝑑𝑉
B. If is negative semidefinite, then 𝒙∗ is a stable critical point.
𝑑𝑡
This method is useful particularly when linearization does not give enough
information. Also, Lyapunov function is not unique.
P.K.Srivastava (IIT Patna)
Example:
Consider the differential equation: 𝑥 ′ = −𝑥 3 .
Note that 𝑉 (𝑥) = 𝑥 2 will serve as Lyapunov functions as 𝑉 ′ = 2𝑥𝑥 ′ = −2𝑥 4 < 0.
Note that any function such as 𝑉 (𝑥) = 𝑋 𝑒 , where e is even integer will serve the
purpose and hence we will have infinitely many Lyapunov function.
It must be emphasized that finding Lyapunov functions is really difficult, in
general.
An example is provided to illustrate it for case when linearization cannot
determine statbility.
P.K.Srivastava (IIT Patna)
𝑑𝑥1
= −2𝑥2 + 𝑥2 𝑥3 − 𝑥13
𝑑𝑡
𝑑𝑥2
Example: = 𝑥1 − 𝑥1 𝑥3 − 𝑥23
𝑑𝑡
𝑑𝑥3 3
= 𝑥1 𝑥2 − 𝑥3
𝑑𝑡
In this case, one can find that origin is a critical point. The eigen values are
0, ±2𝑖. So we cannot say anything about stability of critical point.
Define a function 𝑉 (𝒙) = 𝑥12 + 2𝑥22 + 𝑥32 , so that 𝑉 (𝒙) > 0 and
𝑉′(𝒙) = −2(𝑥14 + 2𝑥24 + 𝑥34 ) < 0, for 𝒙 ≠ 𝟎.
Thus in this case we note that origin is asymptotically stable.
P.K.Srivastava (IIT Patna)
P.K.Srivastava (IIT Patna)
Reference:
1. S.L. Ross, Differential Equations, John Wiley and Sons, 1984.
2. G.F.Simmons, Differential equations with application and historical
notes, TMH, 2003.
3. D. S. Jones, M. J. Plank, B. D. Sleeman, Differential Equations and
Mathematical Biology, CRC Press, 2009.
P.K.Srivastava (IIT Patna)
आ नो भद्राः क्रतवो यन्तु ववश्वताः ||
Let the noble thought come from all directions!
P.K.Srivastava (IIT Patna)