THE SOLUTION STABILITY OF VAN DER POL'S EQUATION .

Maplesoft

Application Demonstration

wwH.map1esaft.cam

THE SOLUTION STABILITY OF
VAN DER POL'S EQUATION.

by CO. H . TRAN .

HUI- NCU HCMC

Vietnam

cothl23@math.com & cohtran @ math . com

Copyright 2004

November 06 2004

Introduction

,5 Signature

** Abstract : The Van der Pol differential quation is solved by averaging method
** Subjects: Vibration Mechanics , The Differential equations .

This worksheet demonstrates Maple's capabilities in finding the graphical solution and dealing with the stability of the steady s

tate solution of Van der Pol 's differential equation .

All rights reserved. Copying or transmitting of this material without the permission of the authors is not allowed .

We consider the Van Der Pol differential equation :_

Two topics that we will be in discussion are

- Finding the steady state solution of this equation by averaging method

- Estimating the stability of solution obtained .

x"+CO 2 x - [i (1 - he 2 ) jc' =

(0)

1 .Define the model Of problem I We examine the effect of non-linear system under external force caused by the
AC generator ( see fig 1. )

£a

E cos ut

C

R
L

i

AC generator
Capacity
Resistance
Inductor
Current intensity

(fig I- )

The differential equation of this model is given in the form

X

Lx"+(x - a)x'-\ — = E cos vt
c

(1)

2/. .2

... ,., . , . . x"+k (x -l)x'+x = ecosvt

Alter simplifying we obtain : v J

2 . Construct the algorithm .

Generally the Van Der Pol differential equation can be expressed by

(2)

APPROVED

By CO HONG TRAN at 1:38 pm, Jun 24, 2009

x"+f(x,x')x'+g(x) = F(t)

In the special case let f(x,x') = jLl(l-?UC ) and g(x)=CO X
the equation will be rewritten as : X +CO X — \\,(\ — hx )X =

By using the transform

X = *fkx;T = cor; \i x = —

CO

(4)

Then

x -* co ^" = * ( ° 2

and substitute these relations to (2) it follows

co 2 X" ! co 2 X

vx vx

Or

|Ll x CO

(

1-

XX

2^

V

raX"_

J

vx

=0

X"+X -\i x (l-X 2 )X'=0

2\„!

x"+jc-ji (1-jc z K=0

Thus we begin with :

To normalize this equation we find the solution which is expressed in the form

x = acos ( t + y )

It is advantageous to write (p = t + y then the solution will be x = acosq*

From the transform x= acosty , x' = - asinty .

We have x" = dx'/dt = -a'sinty - acoscp . (p '
By substituting to ( 6 ) , it gives

(5)
(6)

(7)

-a'sincp -acoscpxp'+acoscp -(0,(1 -a 2 cos 2 (p).(-asincp) = , „,

In the other hand dx/dt = x' = a'costy - asinty .($>' = -asinty (9)

Obviously we reach to the following system of differential equations

-a'sincp -acos(pxp'= -|0.(1-
a'coscp -<2sin(p.(p'=-asin9

a'sincp -acoscp.(p'= -\i{\-a 2 cos 2 (p).asin(p -acoscp

(10)

> A:=Matrix( [ [-sin (phi) , -a*cos (phi) ] , [cos (phi) , -a*sin(phi) ] ] ) ;
-sin(tj)) -acos((j))"

A :=

cos((|)) -asin((|))

> u:=vector ( [-mu*a* (1- (a*cos (phi) ) A 2) *sin(phi) -a*cos (phi) , -a*sin(phi) ] ) ;

u := [-\la (I - a 2 cos((|)) 2 ) sin((|)) - a cos((|)), -a sin((|))]

> S:=simplify (linsolve (A,u) ) ;

S := [a (-1 + cos((|)) 2 ) (I (-1 + a 1 cos((|)) 2 ), [i sin((|)) cos((|)) - (I a 1 sin((|)) cos((|)) 3 + 1 ]

> a(tt) :=simplify (S[l]) ;

a( tt):= a( -1 + cos( ) 2 ) u. ( -1 + a 1 cos( ) 2 )

> phi(tt) :=simplify (S [2]) ;

())(«) := [lsin((|))cos((|))-(ifl 2 sin((|)) cos((|)) 3 + 1

By using the symbols a(tt) = a'(t) , §(tt) = §'(t) (11)

Execute the averaging method for a'(t) and (j)'(^ , we have

> aO (tt) :=normal (int (mu*a* (l-a A 2*cos (phi) A 2) *sin(phi) A 2 , phi=gamma . .gamma+2*Pi) / (2*Pi) ) ;
aO(ff) := --|la (-4 + <r)

o

> phiO (tt) :=int (mu*a* (l-a A 2*cos (phi) A 2) *sin(phi) *cos (phi) ,phi=gamma. .gamma+2*Pi) / (2*Pi) ;

0O(«):=O

The expressions of a ' and (p ' are calculated in the forms :
a ' = |i.< a(ff):=a(-l + cos((|)) 2 ) ^(-1 +<r cos((|)) 2 ) >

1 ?

aO(tt) :=--(! a (-4 + a )

o

(12)

(p' = |i,< §(tt) := (lsin((|))cos((|))-(ia 2 sin((|)) cos((|)) 3 + 1 > _ <t>0(«):-0

The steady state solution occurs when aO = or aO = 2

If aO = then x = x' = this is a trivial solution ( equilibrium ) .

If aO = 2 then x = 2costy and x' = -2sinq> .

The necessary and sufficient condition for solution stability includes

a' =da/dt = *¥(a) with W(ao) = and ¥ 'faoj <

> Psi(a) :=a0 (tt) ;

1 i

¥(a) :=--(! a (-4 + a")

o

> DaohamcuaPsi (a) : =normal (dif f (Psi (a) , a) ) ;

1 3 ,
DaohamcuaPsi( a ) := - (I - - ji a'

1 o

(13)

(14)

ii ii

> print ("Dao ham cua ham a' (t) la

1 3 ,
"Dao ham cua ham a'(t) la : ',' " a"(t)= ",-\i--[ia 2

2 8

> subs (a=2 , DaohamcuaPsi (a) ) ;

-v-

> print ("Gia tri cua a' ' (t) tai a = 2 la

"Gia tri cua a"(t) tai a = 2 la : a"(2) = ',' -|i

> subs (a=0, DaohamcuaPsi (a) ) ;

1

a ' ' (t) = " , DaohamcuaPsi (a) ) ;

a' ' (2) = ", subs (a=2, DaohamcuaPsi (a) )) ;

> print("Gia tri cua a"(t) tai a = la : a''(0) = ", subs (a=0, DaohamcuaPsi (a) )) ;

"Gia tri cua a"(t) tai a = la : a"(0) = ',' ■= \i

Thus if ao = , a"(0) = t-|J. then the solution is not stable

If ao = , a"(0) = -|i then the solution is stable asymtotically.

Note : By solving the equation ( 12 ) for the vibration amplitude a(t) (" slowly varying coefficients ") then finding the solution

expression x(t) of Van Der Pol differential equation , we get

> aO(tt) ;

1 2

--\la(-4 + a )

o

> diff_eq:= dif f (a (t) , t) =-mu*a (t) * (-4+ (a (t) A 2) ) /8;

a i ,

diff-eq := ^a(f ) = - - \i a(f) (-4 + a(0 l )

> init_con:=a (0) =ao;
init_con := a(0 ) = ao

> biendo:= [dsolve ({diff_eq, init_con} , {a(t)})];

a( = 2 j a( = -2 ■

biendo :-

1 -

e Qo 2 -4)
ao 2

> biendo [ 1 ] : =simplif y (biendo [ 1 ] ) ;

biendo := a(0 = 2

1

2 l-t 1 ') 2 h (~^ f )

-ao* + e ao - 4 e

1 -

e (ao 2 -4)
aa 2

aa

( accepted )

> biendo [2 ] : =simplif y (biendo [2 ] ) ;

biendo 9 := a(?) = — 2 ■

1

2 , (-M) 2 ^ <"^ f )

-ao + e ao" - 4 e

ao

( eliminated )

> x:=biendo [1] *cos (phi) ;

x := cos((|)) a(?) = 2

cos((|))

2 H 1 ') 2 ,, <"M

-ao +e ao -4e

ao

with (p = r + y • ( 15 )

> x:=t->2*cos (t + gamma) /sqrt( (a A 2-a A 2*exp(-mu*t)+4*exp(-mu*t) ) /a A 2) ;

cos(? + y)

x := t-> 2

2 2 C-f-0 „ <"^"

a - a e +4e

The graphical solution of Van Der Pol's equation :

> a:=0.5;

a := .5

> y:=t->2*cos(t + gamma) /sqrt ( (a A 2-a A 2*exp (-0 . l*t) +4*exp (-0 . l*t) ) /a A 2) ;

cos(? + 7)

J := t -> 2 t7

2 2 (-.11) (-.If)

a* — a~ e + 4 e

> plot(y(t) ,t=0. .4*Pi) ;

> animate (2 *cos(x*t + gamma) /sqrt ( (a A 2-a A 2*exp (-0 . l*t) +4*exp (-0 . l*t) ) /a A 2) , x=-6*Pi .. 6*Pi,
t=0. .10) ;

0.B-

0.2-

I I I I I I I I I I I I I I I I JtU

-15 -10 -5

-0.2
-0.4
-0.6-

i i i i I i i i i I i i i i I i i i

5 10 15
x

Use graphical method to consider the solution stability of Van Der Pol' s differential equation :

> with(DEtools) :mu:=0.5;

> DEplot({ (D@@2) (x) (t)+x(t)-mu*(l-x(t) A 2)*D(x) (t)=0},{x(t) },t=-

10.. 50, [ [x(0)=l,D(x) (0)=1] ] , stepsize=0.05,title= N Nghiem on dinh cua pt Van Der Pol s ) ;

3 . Conclusion .

From graphical results , we reach to conclusion that the steady state solution stability of Van Der Pol diffential equation must be

precise

and estimating the property of solution obtained is very necessary .

As presented above , we might also use the normalization to ( 2 ) by determining the non-trivial solution in the form

x = Mcoskt + Nsinkt + x*

(16)
x' = -kMsinkt + kNcoskt + x *'

Wl

th k = 1

dxWdt = x" = -M'sinr - M cost + N'cost - Nsint + jc* "

Otherwise

dx/dt = x' = M'cost -Msint + N'sint + Ncost + x* '

= -Msint + Ncost + jc* '

Substitute x " to ( 6 ) after simplifying it follows :

M'cost + N'sint =

- M's'mt + N' cost = [i[l-(M cost + Nsint) 2 ].[-M s'mt + Ncost]

> A:=Matrix( [ [cost, sint] , [-sint, cost] ] ) ;

" cost sint'
A :=

-sint cost

> f :=vector( [ [0] , [mu*F] ] ) ;
/:=[[0],[|iF]]

> V:=linsolve(A, f) ;

sint [-[iF] - [0] cos? sint [0] + [-(IF] cast

V:--

sint 2 + cost 2

sint 2 + cost 2

We rewrite the expressions : M' =
N' = [ M-F] c-o^

(17)

With

F = [I- (M cost + Nsint) 2 ].[-M sint + Ncost]

Use averaging method for ( 17 ) we get :

M' = \H<-Fsint> = --M(4-M 2 + 2N 2 -3 N)

N' = \L< Fcost> = --N (-4-5 M 2 +N)

(18)

The steady state solution exists when M = and N = ( trivial solution

APPROVED

By CO HONG TRAN at 1:38 pm, Jun 24, 2009

Or M= 2 or M= -2 and N=

Thus x = 2cost , x = - 2cost

If M = and N = or N = 4 Then x = 4sint

The steady state solutions can be formed generally :

x = 2cost + 4sinf
x = -2cosf + 4siru

But these forms are equivalent to x= 2costy [see (14)]

Next we begin with Krylov - Bogoliubov approximate method for the Van der Pol's differential equation

x"+(D 2 x-\i(l-hc 2 )x'=0

with •* \ x > x ) = — (I ~ "JC ) x , and (J, is a small constant .

The solution will be estimated by x = r(t) costy(t)

By taking the first order approximate terms ( neglecting the second and third order errors of constant \i ) we can find the
amplitude r(t) and the global phase function (p(t) from the following system

APPROVED

By CO HONG TRAN at 1:38 pm, Jun 24, 2009

2%

ll r r

r\t) = -^-— f (r cos u-r(d sin u) sin udu = — aAr)

2710) J 2

(19>

(p'(0 = C0H / (r cos w,-rco sin u) cos udu = Ja 2 (r)

27irco J

If the initial condition r(0) = ro is satisfied then the solution of ( 15 ) will be equivalent to the solution of second order linear

x"(t) + a 1 (r o )x\t) + a 2 (r o )x(t) =

with the error of estimation based on the order of

(20)

1^

The first order approximate is noticeable in the case of periodic vibration , because the equivalent linear differential equation gives us
the accumulation and dissipation of energy based on vibration period which we might obtain from the given non-linear differential equation
Therefore it is useful to apply the equivalent second linear differential equation to observe the non-linear resonant phenomenon .

REFERENCES

[ 1 ] Menh . C. Nguyen , Dao dong phi tuyen , Manuscript , Vien Co hoc , Hanoi 2002 .

[2] Korn . G , Korn . T ., So tay toan hoc , Trans Vietnamese , NXB DH-THCN , Hanoi 1978

Legal Notice: @ Maplesoft, a division of Waterloo Maple Inc. 2009. Maplesoft and Maple are trademarks of Waterloo Maple Inc.
Neither Maplesoft nor the authors are responsible for any errors contained within and are not liable for any damages resulting
from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the authors for
p ermission if vnu wi sh to use this application in for-profit activities.

Signature

Digitally signed by COHONGTRAN
DN: cn=COHONGTRAN, c=VN, o=MMI,
ou=NCU HUI, email=cohtran@math.com
Reason: I am the author of this document
Location: HCMC
Date: 2009.06.24 13:35:35 +07'00'