A puzzle on the existence of Lyapunov functions for Limit cycle system
Xiao-Liang Gan Hao-Yu Wang Ping Ao







firstaddress,secondaryaddress,thirdaddress]Xiao-Liang Gan
ganxiaoliang21@163.com
firstaddress,secondaryaddress,thirdaddress]Hao-Yu Wang

secondaryaddress,thirdaddress,fourthaddress]Ping Aomycorrespondingauthor
[mycorrespondingauthor]Corresponding author
aoping@sjtu.edu.cn




[firstaddress]Department of Mathematics, Shanghai University, Shanghai 200444, China.

[secondaryaddress]Institute of Systems Science, Shanghai University, Shanghai 200444, China.

[thirdaddress]Shanghai Center for Quantitative Life Sciences, Shanghai University, Shanghai 200444, China.

[fourthaddress]Department of Physics, Shanghai University, Shanghai 200444, China.







Although the limit cycle have been studied for more than 100 years, the existence of its Lyapunov function is still poorly understood.
By considering a common limit cycle system, a puzzle related to the existence of Lyapunov functions for Limit cycle system is proposed and studied in this paper: The divergence is not equal to zero, while the trajectory can be infinite loop on limit cycle.
It will be discussed from three aspects.

Firstly, the definition of Lyapunov function is concerned, and a new version of definition is given which is equivalent to the usual one and better adapts to the properties of Lyapunov function.
Secondly, two criteria(divergence and dissipation power) are discussed from the perspective of dissipation, and it reaches that they are not consistent in representing the dissipation.
Thirdly, by studying the motion of charged massless particle in electromagnetic field, it obtains that the limit cycle is an isopotential line on which the charged particles can move infinitely.
Such discussions above may provide an understanding on the puzzle about the existence of Lyapunov functions for limit cycle system.






limit cycleLyapunov functiondissipation dissipation power divergence isopotential line
[2010] 34C-0737B-25 37N-05








§ INTRODUCE A PUZZLE

Consider a common limit cycle system<cit.>

    {[ ẋ_1 = f_1(x_1,x_2) =  - x_2 + x_1[1 - (x_1^2 + x_2^2)];    ẋ_2 = f_2(x_1,x_2) = x_1 + x_2[1 - (x_1^2 + x_2^2)] ]..

By the polar transformation{[ x_1 = rcosθ; x_2 = rsinθ ]., the (<ref>) can be translated into

    {[ dr/dt = r - r^3;     dθ/dt =   1 ]..

Obviously, it has a limit cycle x_1^2 + x_2^2= 1.

The divergence on the limit cycle of system(<ref>) is

    . ∇· f|_x_1^2 + x_2^2= 1   =   . [ ∂f_1/∂x_1 + ∂f_2/∂x_2]|_x_1^2 + x_2^2 = 1
       =   . [ (1 - 3x_1^2 - x_2^2) + (1 - x_1^2 - 3x_2^2)]|_x_1^2 + x_2^2 = 1
       =   . [ 2(1 - 2x_1^2 - 2x_2^2)]|_x_1^2 + x_2^2 = 1  =  - 2  0.


A puzzle arises: The divergence is not equal to zero, while the trajectory can be infinite loop on limit cycle.
For this phenomenon, researchers have two problems:


  
  Question 1: Whether the system(<ref>) has a Lyapunov function, or not?
  
  Question 2: If the Lyapunov function exists, how to use it to analysis and understand the infinite motion on limit cycle?


Some researchers have explored these two questions and gaven some opinions.
For the Question 1, such as, Strogatz<cit.> believes that for a system, if the Lyapunov function exists, there is no closed orbit.
Hirsch, Smale and Devaney<cit.> confirm that if the planar system has a strict Lyapunov function for a planar system, then there are no limit cycles.
Teschl<cit.> has the similar view with Hirsch.
In a word, the Lyapunov function they mentioned is just for fixed points<cit.>.
For the Question 2, Zhu and Cai<cit.> think that although the limit cycle is also a periodic orbital, the periodic solution and the limit cycle here have different meanings.
Hirsch, Smale and Devaney<cit.> point a property of limit cycle that when t→ +∞ the limit cycle is a limit set. Not all closed orbits have this property, such as the center of linear system.
Lefschetz<cit.> points out that a limit-cycle is a closed path which is not a member of a continuous family of closed paths, i.e. an isolated closed path.

However, they gave no specific explanation to the reason why orbitals move infinitely on limit cycle.
To find an appropriate answers to questions 1 and 2 is difficult, is there an appropriate strategy or approach to address the above problems?  The answer is positive.
Ao<cit.> divided the dynamic system ẋ = f(x) into three parts from the perspective of mechanics:

    (S + T)ẋ =  - ∇ϕ,

here, x = [ x_1, ⋯ ,x_n]^τ, τ is the transpose symbol, the friction matrix S corresponds to dissipation, the transverse matrix T corresponds to Lorentz force and the potential function ϕ is equivalent to Lyapunov function.
Based on the new construction(<ref>), Yuan et al.<cit.> have answered the first question and find a Lyapunov function for system(<ref>)

    ϕ (x_1,x_2) = 1/4(x_1^2 + x_2^2)(x_1^2 + x_2^2 - 2).

The Figure 1 and Figure 2 are the trajectory picture in vector field and Lyapunov function of (<ref>), respectively.


Furthermore, we can verify that the Lyapunov function is decreasing by following trajectory

    dϕ (x_1,x_2)/dt =  - (x_1^2 + x_2^2)(x_1^2 + x_2^2 - 1)^2≤ 0.



In the following, combining with the new construction(<ref>) and results of Yuan et al.<cit.>, we will continue to analyze and explain the puzzle from three aspects in section 2. The conclusion is in section 3.


§ ANALYZE AND EXPLAIN THE PUZZLE



 §.§ From the definition of the Lyapunov function

Firstly, the definition of Lyapunov function in most textbooks is given.

Definition 1<cit.>: Consider the following vector field

    ẋ = f(x),x ∈ℝ^n.

Let x^ * be a fixed point of (<ref>) and let V:Ω→ℝ be a C^1 function defined on some neighborhood Ω of x^ * such that

  
  (i) V(x^ * ) = 0 and V(x) > 0 if x x^ *.
  
  (ii) V̇(x) ≤ 0 in Ω - {x^ * }.

Then x^ * is stable, and the V(x) is called a Lyapunov function for x^ *.

However, there are three kinds of dynamic systems: fixed point, limit cycle and chaos.
Wolfram<cit.> classified the dynamical systems in this way, too.
The Lyapunov function V(x) just can be used for the qualitative analysis of fixed point. So the promotion of Lyapunov function has a deeper and broader practical significance.

Here, only the fixed point and limit cycle are considered, and chaos is not discussed.


Definition 2
For a smooth autonomous system ẋ = f(x), let x^ ** be a limit set which is a fixed point or a limit cycle. And ϕ(x) is a continuous differentiable function from the state space to ℝ.  If it satisfies

  
  (1) for all x in the state space, the Lyapunov function does not increase along the trajectories, that is
  
    ϕ̇(x) = dϕ/dt≤ 0;

  
  (2) for all x, ϕ (x) has an infimum.

Then, ϕ(x) is the generalized Lyapunov function.
















Remark  From the above generalized definition, it can be seen that the fixed point and limit cycle are generalized as the limit set.
The condition (1) ensures that it doesn't increase along trajectory;
The condition (2) generalizes the definition of positive definite.
This definition of Lyapunov function is logically consistent with the one in dynamical system textbooks.








 §.§ Analyze the cause of the puzzle from the perspective of system dissipation

It's easy to understand that divergence does not equal zero for the Lyapunov function of fixed point.
However, when the limit set is limit cycle, it is hard to imagine.
Combined with the decomposition (<ref>), we will analyze the cause from the perspective of system dissipation.

There are two criteria of dissipation:

Criterion 1 Dissipation is defined as dissipation power<cit.>:

    H_P   =   F_friction· ( - ẋ) 
       =    ( - Sẋ) · ( - ẋ)  
       =   ẋ^τSẋ,

here, S is a positive semidefinite friction matrix. 
This method is often used in physics.

Criterion 2 Dissipation is defined as divergence<cit.>:

    div f(x)   =   ∑∂f_i/∂x_i,

this method is often used in mathematics.

Remark  Since the equivalence between potential function and Lyapunov function has been proved by Yuan et al.<cit.>, we can further verify that the rate of change of the Lyapunov function| dϕ/dt| is equal to dissipation powerH_p. That is

    | dϕ/dt|=H_p.



In the following, just the planar dynamical system is taken into consideration, and the divergence and dissipated power will be derived out and compared to see whether they are consistent in representing the dissipation or not.









2.2.1 ẋ = f(x) is a linear system




The two dimension linear dynamical system can be expressed as

    ẋ = f(x) = Ax, A  is  a  constant  matrix, x = [ x_1,x_2]^τ.


Kwon, Ao and Thouless<cit.> have discussed the construction of Lyapunov function of linear system. And some necessary formulas are given

    A=  - [ D + Q]U =  - [ S + T]^ - 1U,


    [D + Q]^ - 1 = S + T,


    AQ + QA^τ = AD - DA^τ.

here, Q is an antisymmetric matrix, D is a semidefinite diffusion matrix, S semidefinite friction matrix,
T is an antisymmetric transverse matrix, U is a symmetric matrix. And only the Q is unknown, the S and T can be derived by D and Q.

By (<ref>)and (<ref>), it derives

    U =  - [ D + Q]^ - 1A =  - A^τ[ D - Q]^ - 1.



Set A = [ [ a_11 a_12; a_21 a_22 ]],Q = [ [      0   Q_12; - Q_12      0 ]],D = [ [ d_11 d_12; d_12 d_22 ]],
and it satisfies
d_11,d_22≥ 0,d_11d_22 - d_12^2 ≥ 0.

By(<ref>), it has

    (a_11 + a_22) Q_12=  - a_21d_11 + (a_11 - a_22)d_12 + a_12d_22,

here, only Q_12 is unknown, a_11,a_12,a_21,a_22 are known, and d_11,d_12,d_22 satisfy
d_11,d_22≥ 0,d_11d_22 - d_12^2 ≥ 0.


Then, it has

    {[ (a_11 + a_22)Q_12 =  - a_21d_11 + (a_11 - a_22)d_12 + a_12d_22;                             d_11,d_22≥ 0,d_11d_22 - d_12^2 ≥ 0 ]..




The matrix A = [ [ a_11 a_12; a_21 a_22 ]]  has the following four types of Jordan's normal form<cit.>





    [ [ λ _1    0;    0 λ _2 ]],[ [ λ _1    0;    0 λ _1 ]],[ [ λ _1    0;    1 λ _1 ]],[ [   α   β; - β   α ]],

here, λ_1, λ_2,α±β i are the different eigenvalues of A. α±β i is a pair of conjugate eigenvalues, and β 0.

In the following, let matrix A be one of the forms (<ref>), and get the corresponding divergence div f(x) and dissipation power H_P.

  
  (1) When A = [ [ λ _1    0;    0 λ _2 ]], (<ref>) can be rewritten as
  
    {[ (λ _1 + λ _2)Q_12 = (λ _1 - λ _2)d_12;    d_11,d_22≥ 0,d_11d_22 - d_12^2 ≥ 0 ].,

  
    
  (i)  If λ _1 + λ _2 0, it can get Q_12 by (<ref>)
    
    Q_12 = λ _1 - λ _2/λ _1 + λ _2d_12,

    then, it has
    
    Q    =   λ _1 - λ _2/λ _1 + λ _2d_12[ [   0   1; - 1   0 ]],  
    
        D + Q    =   [ [                  d_11 2λ _1/λ _1 + λ _2d_12; 2λ _2/λ _1 + λ _2d_12                  d_22 ]], 
    [ D + Q]^ - 1   =   (λ _1 + λ _2)^2/d_11d_22(λ _1 + λ _2)^2 - 4λ _1λ _2d_12^2[ [                    d_22 - 2λ _1/λ _1 + λ _2d_12; - 2λ _2/λ _1 + λ _2d_12                    d_11 ]] 
    
    
    
    
       S    =   [ D + Q]^ - 1 + {[ D + Q]^ - 1}^τ/2
       =   (λ _1 + λ _2)^2/d_11d_22(λ _1 + λ _2)^2 - 4λ _1λ _2d_12^2[ [   d_22 - d_12; - d_12   d_11 ]] ,


    U    =     - (λ _1 + λ _2)^2/(λ _1 + λ _2)^2d_11d_22 - 4λ _1λ _2d_12^2[ [                    λ _1d_22 - 2λ _1λ _2/λ _1 + λ _2d_12; - 2λ _1λ _2/λ _1 + λ _2d_12                    λ _2d_11 ]].

  By  ∇ϕ (x) = Ux, it can derive the Lyapunov function
    
    ϕ  =  - (λ _1 + λ _2)^2(d_22λ _1x_1^2 - 4λ _1λ _2/λ _1 + λ _2d_12x_1x_2 + d_11λ _2x_2^2)/d_11d_22(λ _1 + λ _2)^2 - 4λ _1λ _2d_12^2,

    Then, it can verify
    
    dϕ/dt   =   ∂ϕ/∂x_1ẋ_1 + ∂ϕ/∂x_2ẋ_2
       =    - (λ _1 + λ _2)^2(d_22λ _1^2x_1^2 - 2λ _1λ _2d_12x_1x_2 + d_11λ _2^2x_2^2)/(λ _1 + λ _2)^2d_11d_22 - 4λ _1λ _2d_12^2≤ 0.

   Then, the corresponding divergence div f(x) and dissipation powerH_P are derived
  
    H_P   =   ẋ^τSẋ = (λ _1 + λ _2)^2(d_22λ _1^2x_1^2 - 2d_12λ _1λ _2x_1x_2 + d_11λ _2^2x_2^2)/d_11d_22(λ _1 + λ _2)^2 - 4λ _1λ _2d_12^2 0,
    
        div f(x)     =   ∇· f = trace (A)=λ _1 + λ _2 0.

  And then, it can verify
   
    | dϕ/dt|=H_p.

     In this case, div f(x) and H_P are not zero at the same time. And they all show that the system is dissipative. So these two criteria are consistent in representing the dissipation of the system.
    
  (ii) If λ _1 + λ _2=0, it can get Q_12 by (<ref>)
    
    d_12=0, Q_12 is  an  arbitrary  real  number.

  then, it has
  
    D    =   [ [ d_11    0;    0 d_22 ]],d_11,d_22≥ 0 , 
    
       Q     =   [ [      0   Q_12; - Q_12      0 ]], 
    
       D + Q    =   [ [   d_11   Q_12; - Q_12   d_22 ]], 
    [ D + Q]^ - 1   =   1/d_11d_22 + Q_12^2[ [   d_22 - Q_12;   Q_12   d_11 ]], 
    
       S    =   1/d_11d_22 + Q_12^2[ [ d_22    0;    0 d_11 ]],

  
    U =  - [ D + Q]^ - 1A =  - λ _1/d_11d_22 + Q_12^2[ [   d_22   Q_12;   Q_12 - d_11 ]],

   By  ∇ϕ (x) = Ux, it can derive the Lyapunov function
    
    ϕ  =  - λ _1(d_22x_1^2 + 2Q_12x_1x_2 - d_11x_2^2)/d_11d_22 + Q_12^2,

   Then, it can verify
    
    dϕ/dt = ∂ϕ/∂x_1ẋ_1 + ∂ϕ/∂x_2ẋ_2 =  - λ _1^2(d_22x_1^2 + d_11x_2^2)/d_11d_22 + Q_12^2≤ 0.

  The corresponding divergence div f(x) and dissipation powerH_P are derived
  
    H_P   =   ẋ^τSẋ =λ _1^2(d_22x_1^2 + d_11x_2^2)/d_11d_22 + Q_12^2 0,  
    
       div f(x)    =   ∇· f = trace (A)=λ _1 + λ _2=0.

  
  Then, it can verify
   
    | dϕ/dt|=H_p.





  In this case, div f(q)=0 doesn't mean H_P=0. So these two criteria are not consistent in representing the dissipation of the system.
  
  (2) WhenA = [ [ λ _1    0;    0 λ _1 ]], (<ref>) can be rewritten as




   
    {[                      2λ _1Q_12 = 0; d_11,d_22≥ 0,d_11d_22 - d_12^2 ≥ 0 ].

   
     
  (i) Ifλ_1  0, it can get Q_12 by (<ref>)
   
    Q_12 = 0,

  then, it has
  
    Q     =    0,  
    
      D + Q    =   [ [ d_11 d_12; d_12 d_22 ]], 
    [ D + Q]^ - 1   =   1/d_11d_22 - d_12^2[ [   d_22 - d_12; - d_12   d_11 ]], 
    
        S    =   1/d_11d_22 - d_12^2[ [   d_22 - d_12; - d_12   d_11 ]],

   
    U =  - [ D + Q]^ - 1A =  - λ _1/d_11d_22 + d_12^2[ [   d_22 - d_12; - d_12   d_11 ]],

   By ∇ϕ (x) = Ux, it can derive the Lyapunov function
   
    ϕ  =  - λ _1/d_11d_22 + d_12^2(d_22x_1^2 - 2d_12x_1x_2 + d_11x_2^2),

  Then, it can verify
   
    dϕ/dt = ∂ϕ/∂x_1ẋ_1 + ∂ϕ/∂x_2ẋ_2 =  - λ _1^2(d_22x_1^2 - 2d_12x_1x_2 + d_11x_2^2)/d_11d_22 + d_12^2≤ 0.

  The corresponding divergence div f(x) and dissipation powerH_P are derived
  
    H_P   =   ẋ^τSẋ = λ _1^2(d_22x_1^2 - 2d_12x_1x_2 + d_11x_2^2)/d_11d_22 - d_12^2 0 . 
    
       div f(x)   =   ∇· f = trace (A)=2λ _1 0 .

  Then, it can verify
   
    | dϕ/dt|=H_p.

  In this case, div f(x) and H_P are not zero at the same time. And they all show that the system is dissipative. So these two criteria are consistent in representing the dissipation of the system.
     
  (ii) If λ_1 =0, the system is conservative. It's easy to know
  
    H_P   =   ẋ^τSẋ≡ 0,  
    
        div f(x)     =   ∇· f = trace (A)=2λ _1≡ 0 .

  Then, it can verify
   
    | dϕ/dt|=H_p.

    In this case, they are equal to zero at the same time. div f(x)≡ 0 and H_P≡0 are consistent in representing the system being conservative.
   
  
  (3)When A = [ [ λ _1    0;    1 λ _1 ]], (<ref>) can be rewritten as




   
    {[                2λ _1Q_12 =  - d_11; d_11,d_22≥ 0,d_11d_22 - d_12^2 ≥ 0 ]..

   
     
  (i) If λ_1  0,  it can get Q_12 by (<ref>)

     
    Q_12 =  - d_11/2λ _1,

  then, it has
  
    Q     =    - d_11/2λ _1[ [   0   1; - 1   0 ]], 
    
      D + Q    =   [ [              d_11 d_12 - d_11/2λ _1; d_12 + d_11/2λ _1              d_22 ]], 
    [ D + Q]^ - 1   =   4λ _1^2/4λ _1^2d_11d_22 - 4λ _1^2d_12^2 + d_11^2[ [                d_22 - d_12 + d_11/2λ _1; - d_12 - d_11/2λ _1                d_11 ]], 
    
       S     =   4λ _1^2/4λ _1^2d_11d_22 - 4λ _1^2d_12^2 + d_11^2[ [   d_22 - d_12; - d_12   d_11 ]],

   
    U    =    - 4λ _1^2/4λ _1^2d_11d_22 - 4λ _1^2d_12^2 + d_11^2[ [ d_22λ _1 - d_12 + d_11/2λ _1          - d_12λ _1 + d_11/2;          - d_12λ _1 + d_11/2                     d_11λ _1 ]],

   By ∇ϕ (x) = Ux, it can derive the Lyapunov function
     
    ϕ  =  - 4λ _1^2[ (d_22λ _1 - d_12 + d_11/2λ _1)x_1^2 + 2( - d_12λ _1 + d_11/2)x_1x_2 + d_11λ _1x_2^2]/4λ _1^2d_11d_22 - 4λ _1^2d_12^2 + d_11^2,

  Then, it gets



     
    dϕ/dt =- 4λ _1^2[ (d_22λ _1^2 - 2d_12λ _1 + d_11)x_1^2 + 2( - d_12λ _1 + d_11)λ _1x_1x_2 + d_11λ _1^2x_2^2]/4λ _1^2d_11d_22 - 4λ _1^2d_12^2 + d_11^2,

      the root discriminant of (d_22λ _1^2 - 2d_12λ _1 + d_11)x_1^2 + 2( - d_12λ _1 + d_11)λ _1x_1x_2 + d_11λ _1^2x_2^2=0 is
     
    Δ   =    4( - d_12λ _1 + d_11)^2λ _1^2 - 4(d_22λ _1^2 - 2d_12λ _1 + d_11)d_11λ _1^2 
       =     4λ _1^4[ d_12^2 - d_11d_22] ≤ 0,

     it has (d_22λ _1^2 - 2d_12λ _1 + d_11)x_1^2 + 2( - d_12λ _1 + d_11)λ _1x_1x_2 + d_11λ _1^2x_2^2≥ 0, then it obtains
     
    dϕ/dt=∂ϕ/∂x_1ẋ_1 + ∂ϕ/∂x_2ẋ_2≤ 0.

  The corresponding divergence div f(q) and dissipation powerH_P are derived
  
    H_P   =   4λ _1^2[(λ _1^2d_22^2 - 2λ _1d_12 + d_11)x_1^2 + 2( - λ _1^2d_12 + λ _1d_11)x_1x_2 +
          λ _1^2d_11x_2^2]/4λ _1^2d_11d_22 - 4λ _1^2d_12^2 + d_11^2
       ≠    0, 
    
       div f(x)     =   ∇· f = trace (A)=2λ _1 0 .

   Then, it can verify
   
    | dϕ/dt|=H_p.

  In this case, div f(x) and H_P are not zero at the same time. And they all show that the system is dissipative. So these two criteria are consistent in representing the dissipation of the system.
     
  (ii) If λ_1 = 0, it can get Q_12 by (<ref>)
     
    d_11=d_12=0, Q_12 is  an  arbitrary  real  number.

  then, it has
  
    D    =   [ [    0    0;    0 d_22 ]],  
    
       Q    =   [ [      0   Q_12; - Q_12      0 ]],   
    
       D + Q    =   [ [      0   Q_12; - Q_12   d_22 ]],  
    [ D + Q]^ - 1   =   1/Q_12^2[ [   d_22 - Q_12;   Q_12      0 ]], 
    
       S    =   1/Q_12^2[ [ d_22    0;    0    0 ]],

  
    U =  - [ D + Q]^ - 1A = 1/Q_12[ [ 1 0; 0 0 ]],

   By ∇ϕ (x) = Ux, it can derive the Lyapunov function
     
    ϕ  = 1/Q_12x_1^2,

  Then, it can verify that the Lyapunov function is decreasing by trajectory
     
    dϕ/dt = ∂ϕ/∂x_1ẋ_1 + ∂ϕ/∂x_2ẋ_2 = 2x_1ẋ_1/Q_12= 2λ _1x_1^2/Q_12 = 0,

  then, the corresponding divergence div f(x) and dissipation powerH_P are derived
  
    H_P   =   ẋ^τSẋ=  x^τA^τSAx 
       =   1/Q_12^2[ [ x_1 x_2 ]][ [ 0 1; 0 0 ]][ [ d_22    0;    0    0 ]][ [ 0 0; 1 0 ]][ [ q_1; q_2 ]] 
       ≡    0, 
    
       div f(x)     =   ∇· f = trace (A)=2λ _1≡ 0 .

      Then, it can verify
   
    | dϕ/dt|=H_p.

  In this case, they are equal to zero at the same time. div f(x)≡ 0 and H_P≡0 are consistent in representing the system being conservative.
   
  
  (4)When A=[ [   α   β; - β   α ]], (<ref>) can be rewritten as

    {[               2αQ_12 = β (d_11 + d_22); d_11,d_22≥ 0,d_11d_22 - d_12^2 ≥ 0,β 0 ]..


  
  (i) If α 0, it can get Q_12 by (<ref>)
  
    Q_12 = β (d_11 + d_22)/2α,

  then, it has
  
    Q    =   β (d_11 + d_22)/2α[ [   0   1; - 1   0 ]], 
    
       D + Q    =   [ [                      d_11 d_12 + β (d_11 + d_22)/2α; d_12 - β (d_11 + d_22)/2α                      d_22 ]],  
    [ D + Q]^ - 1   =   4α ^2/4α ^2(d_11d_22 - d_12^2) + β ^2(d_11 + d_22)^2[ [                        d_22 - d_12 - β (d_11 + d_22)/2α; - d_12 + β (d_11 + d_22)/2α                        d_11 ]], 
    
       S    =   4α ^2/4α ^2(d_11d_22 - d_12^2) + β ^2(d_11 + d_22)^2[ [   d_22 - d_12; - d_12   d_11 ]],

  
  
    U    =    - 4α ^2/4α ^2(d_11d_22 - d_12^2) + β ^2(d_11 + d_22)^2[ [ d_22α  + d_12β  + β ^2(d_11 + d_22)/2α           - d_12α  + β (d_22 - d_11)/2;           - d_12α  + β (d_22 - d_11)/2 d_11α  - d_12β  + β ^2(d_11 + d_22)/2α ]],

  
   By ∇ϕ (x) = Ux, it can derive the Lyapunov function
     
     
    ϕ  =  - 4α ^2{[d_22α  + d_12β  + β ^2(d_11 + d_22)/2α]x_1^2 + 2[ - d_12α  + β (d_22 - d_11)/2]x_1x_2 + [d_11α  - d_12β  + β ^2(d_11 + d_22)/2α]x_2^2}/4α ^2(d_11d_22 - d_12^2) + β ^2(d_11 + d_22)^2,

     
     then, it gets
     
     
    dϕ/dt   =    - 4α ^2{(d_22α ^2 + 2d_12αβ  + d_11β ^2)x_1^2 + 2[d_12(β ^2 - α ^2) + αβ (d_22 - d_11)]x_1x_2 + (d_11α ^2 + d_22β ^2 - 2d_12αβ )x_2^2}/4α ^2(d_11d_22 - d_12^2) + β ^2(d_11 + d_22)^2,

     
     the root discriminant of equation {  *  } = 0 corresponding to the formula in (<ref>) is

     
     
    Δ   =    4[d_12(β ^2 - α ^2) + αβ (d_22 - d_11)]^2 - (d_22α ^2 + 2d_12αβ  + d_11β ^2)(d_11α ^2 + d_22β ^2 - 2d_12αβ )
       =    4(α ^2 + β ^2)^2(d_12^2 - d_11d_22) ≤ 0,

     
     then, it has
     
    dϕ/dt = ∂ϕ/∂x_1ẋ_1 + ∂ϕ/∂x_2ẋ_2≤ 0,

  then, the corresponding divergence div f(x) and dissipation powerH_P are derived
     
    
    H_P   =   ẋ^τSẋ
       =   4α ^2{(α ^2d_22 + 2αβd_12 + β ^2d_11)x_1^2 + 2[αβd_22 - (α ^2 - β ^2)d_12 - αβd_11]x_1x_2 + (β ^2d_22 - 2αβd_12 + α ^2d_11)x_2^2}/4α ^2(d_11d_22 - d_12^2) + β ^2(d_11 + d_22)^2
           0  , 
    
         div f(x)    =   ∇· f = trace (A)=2α 0  .

    
    Then, it can verify
   
    | dϕ/dt|=H_p.

   In this case, div f(x) and H_P are not zero at the same time. And they all show that the system is dissipative. So these two criteria are consistent in representing the dissipation of the system.
  
  (ii) If α=0, it can get Q_12 by (<ref>)
  
    Q_12 is  an  arbitrary  real  number,d_11 = d_22 = d_12 = 0,

  then, it has
  
    Q    =   [ [      0   Q_12; - Q_12      0 ]],
    
       D     =    0, 
    
      D + Q    =   [ [      0   Q_12; - Q_12      0 ]], 
    [D + Q]^ - 1   =   1/Q_12^2[ [      0 - Q_12;   Q_12      0 ]], 
    
       S     =   [ D + Q]^ - 1 + {[ D + Q]^ - 1}^τ/2≡ 0,

   
    U =  - [ D + Q]^ - 1A =  - β/Q_12[ [ 1 0; 0 1 ]],

  By ∇ϕ (x) = Ux, it can derive the Lyapunov function
  
    ϕ  =  - β/Q_12(x_1^2 + x_2^2),

   The Lyapunov function is decreasing by following trajectory
  
    dϕ/dt   =   ∂ϕ/∂x_1ẋ_1 + ∂ϕ/∂x_2ẋ_2 =  - 2β/Q_12(x_1ẋ_1 + x_2ẋ_2)
       =    - 2β/Q_12[x_1βx_2 + x_2( - βx_1)] = 0,

  then, the corresponding divergence div f(x) and dissipation powerH_P are derived
  
    H_P   =   ẋ^τSẋ≡ 0. 
    
       div f(x)    =   ∇· f = trace (A)≡0  .

  Then, it can verify
   
    | dϕ/dt|=H_p.

  In this case, they are equal to zero at the same time. div f(x)≡ 0 and H_P≡0 are consistent in representing the system being conservative.


In a word, when the fixed point is a saddle point, the divergencediv f(x)=0 while the dissipation power H_P 0.
So these two criteria are not always consistent in representing the dissipation of the linear system.






2.2.2 ẋ = f(x) is a non-linear system



In this case, we will analyze it in combination with system (<ref>).
By the paper <cit.>, S(x) can be constructed
if the Lyapunov function ϕ(x)is known for the nonlinear system

    S(x) =  - ∇ϕ· f/f · fI,

here, I represents the identity matrix.


Combined with (<ref>) and (<ref>), it can obtain the S(x) of system (<ref>)

    S = (1 - x_1^2 - x_2^2)^2/1 + (1 - x_1^2 - x_2^2)^2[ [ 1 0; 0 1 ]].

Furthermore,

    H_P   =    (x_1^2 + x_2^2)(x_1^2 + x_2^2 - 1)^2,  
    
    div f   =   ∇· f = ∂f_1/∂x_1 + ∂f_2/∂x_2 = 2(1 - 2x_1^2 - 2x_2^2).

Then, combined with (<ref>)and (<ref>), it can verify

    | dϕ/dt|=H_p.

We can find that on limit cycle x_1^2 + x_2^2 = 1, it gets . ∇· f|_x_1^2 + x_2^2 = 1 = -2  0, H_P =0.





So these two criteria are not consistent in representing the dissipation of the nonlinear dynamical systems.

To sum up, these two criteria are not consistent in representing the dissipation.







 §.§  Analyze the movement on the limit cycle from the physical point of view

In the following, we will discuss and analyze the movement on the limit cycle from the physical point of view.


Yuan et al.<cit.> have given a physical point of view on dynamical system. For the sake of completeness, the process is given, and we will use it to analyze and explain the motion on the limit cycle.


For a charged massless particle in the electromagnetic field, its motion Newton equation is

    F_driving = mẍ = 0,

here, F_driving is usually divided into two parts: friction force and conservative force

    F_driving = F_conservative + F_friction.


The F_friction corresponds to the dissipative part in (<ref>), and it is expressed as

    F_friction =  - Sẋ,

 S is a positive semi-definite matrix which ensures the property that friction impedes the motion of an object. So it's called the friction matrix.

The conservative force is derived by lorentz force eẋ× B and electric field force - ∇ϕ (x) which is the potential gradient

    F_conservative = eẋ× B + [ - ∇ϕ (x)],

here, ∇× x = ∂ _ix_j - ∂ _jx_i.

By (<ref>), it can get

    Sẋ + eB ×ẋ= - ∇ϕ (x) .


Change the lorentz force B ×ẋ to Tẋ, T is an anti-symmetric matrix and named Lorentz-Ao matrix. Then (<ref>) can be represented as


    [S(x) + T(x)]ẋ =- ∇ϕ (x) .


There is another form of (<ref>) in paper <cit.>


    ẋ =  - [D(x) + Q(x)]∇ϕ (x).

here, D(x) is a semidefinite diffusion matrix, Q(x) an anti-symmetric matrix which denotes Poisson bracket.

















Combine with <cit.> and the Lyapunov function (<ref>), it can derive the corresponding S, T, D, Q

    S = (1 - x_1^2 - x_2^2)^2/1 + (1 - x_1^2 - x_2^2)^2[ [ 1 0; 0 1 ]],

here, the friction matrix S corresponds to dissipation.

    T = (1 - x_1^2 - x_2^2)/1 + (1 - x_1^2 - x_2^2)^2[ [   0   1; - 1   0 ]],

here, the Lorentz-Ao matrix corresponds to the action of the magnetic field.

    - ∇ϕ  =  - [ [ ∂∇ϕ/∂ x_1; ∂∇ϕ/∂ x_2 ]] = [ [ x_1(1 - x_1^2 - x_2^2); x_2(1 - x_1^2 - x_2^2) ]],

here, the potential gradient - ∇ϕ corresponds to the action of the electric field.

    D(x) =   [ [ 1 0; 0 1 ]],

here, D(x) is the diffusion matrix.

    Q(x) = 1/1 - x_1^2 - x_2^2[ [   0   1; - 1   0 ]].

here, Q(x) corresponds to the Poisson bracket.



































Here, analyze the motion of charged particles on limit cycles

  
  * The friction part:F_friction =  - Sẋ=0 , so H_P = ẋ^τSẋ = 0 ;
  
  * The conservative force part:
  
  
    
  * on the limit cycle, the electric potential is equal, so the electric field forces do no work.
    
  * Lorentz forces are everywhere perpendicular to the direction of motion of the particle and do no work.
  
  
  To sum up, the system is conservative on limit cycle, and the charged particle in electromagnetic field can move along the isopotential line.
  Here, the isopotential line is the limit cycle, so it moves in circle.

  Remark: By(<ref>), the magnetic field intensity is proportional to the electric field intensity on the limit cycle: T(x)ẋ =  - ∇ϕ (x).
   It only affects the size of the circle of motion, not changes the energy of the system.






§ CONCLUSIONS

Based on the decomposition of the dynamical system of Ao, we analyze the puzzle about the existence of Lyapunov functions for limit cycle system.
The generalized definition of Lyapunov function is logically consistent with the definition in dynamical systems textbooks, and it better adapts to the properties of Lyapunov function. By discussing the two criteria of dissipation(divergence and dissipation power), we obtain that they are not consistent in representing the dissipation. The movement on limit cycle is analyzed from the physical point of view, and it obtains that the system is conservative on the limit cycle and the limit cycle is an isopotential line, so it can be infinite motion on limit cycle.
What's more, it may provide an understanding on the puzzle about the existence of Lyapunov function for limit cycle system.



§ ACKNOWLEDGEMENTS

Thanks members in the Institute of Systems Science of Shanghai University for discussion, particularly, with Xin-Jian Xu.

                                                                                              





§ REFERENCES




