Anali PAZU - Letnik 3, leto 2013, številka 1 
  
15 
 
The Control of Nonlinear Oscillatory Systems with Delay  
Upravljanje nelinearnih nihajočih sistemov z zakasnitvami 
Rudolf Pušenjak 
1,
*, Maks Oblak 
1 
 
1
 Faculty of Industrial Engineering Novo mesto / Šegova ulica 112, SI-8000 Novo mesto 
E-Mails: rudi.pusenjak@teleing.com; maks.oblak@fl.uni-mb.si 
* Full Prof. Dr. Rudolf Pušenjak, Kidričeva 5, 9240 Ljutomer, Slovenia; Tel.: +386-2-582-14-06 
 
Abstract: This paper treats the adaptation of the Extended Lindstedt-Poincare Method with multiple 
time scales (EL-PM) for analysis of stationary and nonstationary resonances of harmonically excited 
Duffing oscillator with time delay feedback control. The fundamental and ⅓ subharmonic resonance, 
respectively of Duffing oscillator with feedback control are analyzed in details. The comparison of 
stationary resonances with nonstationary resonances due to the slowly varying excitation frequency 
is presented by means of examples, which are computed by programming tool Mathematica
®
.  
Key words: Duffing oscillator, time delay feedback control, osnovna in subharmonična resonanca, 
EL-P method. 
Povzetek: Članek obravnava priredbo Razširjene Lindstedt-Poincarejeve Metode z več časovnimi 
skalami (RL-PM) za analizo stacionarnih in nestacionarnih resonanc harmonično vzbujanega 
Duffingovega nihala z zakasnitvijo povratno zančnega upravljanja. Analizi osnovne in ⅓ 
subharmonične resonance Duffingovega nihala z upravljanjem sta prikazani v podrobnostih. V članku 
je s pomočjo zgledov izvedena primerjava stacionarnih resonanc z nestacionarnimi resonancami, ki 
jih povzroča počasno spremenljiva vzbujevalna frekvenca, pri čemer so rezultati dobljeni s pomočjo 
programskega orodja Mathematica
®
.  
Ključne besede: Duffingovo nihalo, upravljanje s časovno zakasnitvijo, osnovna in subharmonična 
resonanca, RL-P metoda.  
 
1. Introduction 
The control of nonlinear dynamical systems, which 
can exhibit periodic, aperiodic and even chaotic 
oscillations, represents a great research challenge due 
to the theoretical achievements and practical 
applications. Realization of feedback control is 
connected with inevitable delays, which occur in data 
acquisition and signal processing. Mathematical 
modelling of such dynamical systems is described in 
the form of nonlinear delayed differential equations 
(DDE's), which cannot be handled by the standard 
methods. Inspired by the power of the Extended 
Lindstedt-Poincare Method with multiple time scales 
(EL-PM) [1], the goal of this paper is to extend its 
applicability to the problems of time delay nonlinear 
control in the steady-state of dynamical systems, 
which are harmonically excited. In the paper, the 
method is applied to the feedback control of the 
harmonically excited Duffing oscillator.  
2. Nomenclature 
w      =   deflection of the beam 
u      =   displacement 
x         =   dimensionless spatial variable 
t       =   time 
c
*
        =   coefficient of viscous damping 
        =  axial force 
K       =   beam stiffness 
ω L      =  natural angular frequency 
ω 0      =  linear angular frequency 
ω        =  excitation angular frequency 
β        =  coefficient of cubic nonlinearity 
θ        =  time delay 
ε         =   perturbation parameter 
a        =  gain of  the proportional part of the controller 
b       =  gain of the differential part of the controller 
p         =   excitation amplitude 
η       =  subharmonic factor 
σ       =   detuning 
τ 1       =   fast time scale 
τ 2        =   slow time scale 
A(τ 2)  =  amplitude of oscillations as 
                function of the slow time scale 
Ф 0(τ 2) =  phase angle of oscillations as 
                function of the slow time scale 
γ(τ 2)   =  auxiliary variable 
r        =  rate of detuning change 

Anali PAZU - Letnik 3, leto 2013, številka 1 
  
16 
 
3. Buckling of the hinged-hinged beam: An example 
of derivation of Duffing equation 
Duffing oscillator is a dynamical system with one 
degree of freedom and a cubic nonlinearity. Examples 
of Duffing oscillator can be found elsewhere in science 
and engineering. In mechanics, for example, vibrations 
of active suspension system, which takes into account 
the nonlinearity of tires, vibrations of offshore 
platforms, vibrations of hinged-hinged beams, 
clamped-hinged beams [2], etc. are modelled by means 
of Duffing oscillator. In modelling rotordynamics, 
even coupled nonlinear oscillators of Duffing type are 
used [3]. In this paper, we consider the buckling 
phenomenon of the hinged-hinged beam as an 
example. Deflection of the beam is presented on the 
Figure 1 and for the sake of simplicity, the beam is 
analyzed through the dimensionless spatial variable 
x=l/L on the interval 0 ≤ x ≤ 1, where L is the length of 
the beam and l is the variable distance, measured from 
its left end. The beam is at both ends subjected to the 
compressive axial force  and externally forced by a 
harmonic excitation P(x,t), which depends on the 
spatial variable x and on the time t. Excitation P(x,t) 
can be realized in practice, when the beam supports a 
machine with a rotating imbalance. Viscous damping 
of the beam is characterized by the coefficient c
*
 and 
the beam stiffness is denoted by K. By using Hamilton 
principle [4], we can derive the governing PDE of 
small deflections w(x,t), where longitudinal motions of 
the beam are neglected:
    
2 4 2 2
1
*
4 2 2
0
, d ,
w w w w w
K t c P x t
xt
x x t
  

      
    



   


,  (1) 
where 
    
   
22
22
0, 1,
0, 1, 0, 0
w t w t
w t w t
xx

   

 (2)
are the corresponding boundary conditions of the 
hinged-hinged beam at its supports. By simplifying the 
analysis assuming the symmetry of P(x,t) around the 
spatial coordinate x=0.5, the beam is in the first mode 
oscillation and both excitation P(x,t) and deflection 
w(x,t), respectively have a sinusoidal shape in x. 
However, excitation P(x,t) additionally varies 
cosinusoidal in dependence on the time. Taking both 
spatial and temporal variations into account, we seek 
the solution of Eq. (1) by means of ansatz:  
      , cos sin , , sin P x t p t x w x t u t x     , (3) 
By substituting Eq. (3) into Eq. (1), by separation of 
variables and integration involved and by using 
boundary conditions (2), we can derive the Duffing 
equation: 
 
* 2 2 4 3
1
2
cos u c u u K u p t           , (4) 
which governs the temporal vibration u(t) of the beam 
displacement. By using notations:
 
 
1
2 * 2 2 4
1
2
0
2 , , ,
k
k
Lk
k
c c K p p         


      , (5)
which can be used to introduce small quantities, 
depending on the small perturbation parameter ε (with 
an exception of introducing the natural angular 
frequency ω L), Eq. (4) can be rewritten on the standard 
form: 
1
2 3
0
2 cos
k
k
Lk
k
u cu u u p t      


     , (6) 
which is analyzed in the rest of paper in the general 
case.

Anali PAZU - Letnik 3, leto 2013, številka 1 
  
17 
 
 
Fig. 1. Deflection w of the hinged-hinged beam
4. Harmonically excited Duffing oscillator with 
time delay feedback control 
Despite of simple harmonically excitation, Duffing 
oscillator is able to produce various periodic, 
apperiodic or even chaotic oscillations in dependence 
on parameter values. Chaotic oscillations are in 
general unwanted and therefore give rise to the use of 
active feedback control in order to prevent them. 
Active control is advantageous over the passive control 
in many respects, however it opens numerous 
problems of technical nature. Due to the active control, 
the dynamics of the system becomes more complicated 
as the consequence of inevitable time delays in 
controllers, transducers and actuators such as analogue 
filters, hydraulic and pneumatic actuators. Time delays 
can cause instability of the system, which is controlled 
by the use of feedback loop. In the presence of time 
delays, the classical control theory has a limited use in 
linear systems, however it fails in the case of nonlinear 
systems. In order to facilitate the analysis of the 
nonlinear dynamical systems with time delay 
feedback, the Extended Lindstedt-Poincare method 
with multiple time scales (EL-PM) [1] is proposed in 
this paper in the case of Duffing oscillator. The 
analysis of Duffing oscillator with time delay control 
is performed by assumption of a small modal damping, 
weak cubic nonlinearity and small gain coefficients of 
controller. Characteristic for the use of EL-PM is, that 
all of these parameters are expressed in dependence on 
a small perturbation parameter ε. In the case of 
fundamental resonance, the external harmonic 
excitation is also small and depending on the 
perturbation parameter ε, however, this restriction is 
removed in the case of subharmonic resonance. In 
accordance with this assumptions, the governing 
equation of Duffing oscillator with harmonic 
excitation and cubic nonlinearity (6), which applies the 
time delay feedback control, is rearranged in the form:
 
   
       
2
1
23
2
0
dd d
2 2 2 cos
dd
d
k
k
Lk
k
u t u t
c u t u t au t b u t p t
tt
t
         


         , (7)
where u(t) denotes displacements of the dynamical 
system, 2εc represents modal damping of the system, 
ω L means natural (angular) frequency, β means 
coefficient of cubic nonlinearity, a and b mean the gain 
coefficients of the proportional and differential part of 
the controller, respectively, θ means the time delay, ω 
denotes the excitation (angular) frequency and where 
p k , (k=0,1) are excitation amplitudes, which are 
adapted in such a way, that p 0=0 , p 1 =p holds in the 
case of fundamental resonance and p 0=p , p 1 =0 holds 
in the case of subharmonic resonance. In the sequel, 
EL-PM with multiple time scales [1] will be applied in 
the analysis of nonstationary resonances of Duffing 
oscillator with time delay control. Nonstationary 
resonances are the phenomenon, which occur at slowly 
varying excitation frequency or excitation amplitude, 
respectively. In this paper, only slowly varying 
excitation frequency will be taken into account. 
Stationary resonances will be also studied as the 
special cases, when the excitation frequency has a role 
of an adjustable, but otherwise constant parameter. 
In order to master both fundamental as well as 
subharmonic resonances, one introduces the so called 
subharmonic factor , which has the value =1 at the 
fundamental (or even at the superharmonic) resonance, 
but it becomes an adequate natural number, such as 
=3 at the ⅓ subharmonic resonance, for example. 
Deviation of the excitation frequency ω from the linear 
frequency ω 0 can be expressed by equation:  
 
0
      ,   (8) 
where the parameter σ stands for the detuning. On the 
basis of Eq. (8), the fast time scale τ 1 and the slow time 
scale τ 2, respectively, are introduced: 
 

Anali PAZU - Letnik 3, leto 2013, številka 1 
  
18 
 
 
10
t    (9.a)   
 
2
t   , (9.b) 
which give rise to the replacement of time derivatives 
d/dt in 
22
d /dt with differential operators: 
                     
0
12
d
dt





 (10.a)
 
            
2 2 2 2
22
00
2 2 2
12
1 2
d
2
dt
   

 
  
  

 
.  (10.b)
 
Using partial derivatives (10.a,b), Eq. (7) is rewritten into nonlinear partial differential equation of the form: 
 
       
2 2 2
2 2 2 3
0 0 0
22
1 2 1 2
1 2
2
1 0 2 0 1 0 2 1 0 2 1 2
12
22
2 , 2 , 2 , cos
L
u u u u u
c u u
au u b u p
       
   
 
                   

     
     

   
  

       

 (11.a) 
in the case of fundamental resonance and into equation 
 
       
2 2 2
2 2 2 3
0 0 0
22
1 2 1 2
1 2
2
1 0 2 0 1 0 2 1 0 2 1 2
12
22
2 , 2 , 2 , cos
L
u u u u u
c u u
au b u b u p
       
   
 
                  

     
     

   
  

       

 (11.b)
in the case of subharmonic resonance. The 
approximate solution of Eqs. (11.a,b) in dependence of 
time scales τ 1 and τ 2 can be found  by applying the 
perturbation procedure, when one expresses 
displacements u(t) in the form of the power series: 
  
12
0
,
k
k
k
uu   


  . (12) 
5. Analysis of nonstationary fundamental 
resonance of Duffing oscillator with delay feedback 
control 
When power series (12) is introduced into Eq. 
(11.a) and terms of like powers of ε on both sides of 
equation are equated, we get the following set of 
partial differential equations, that are linear, if are 
successively solved:
 
                                         
2
22 0
00
2
1
0
L
u
u 




   (13) 
                   
 
   
2 2
3 22 00 1
0 1 1 2 0 0
2
1 2 1
1
0 1 0 2 0 0 1 0 2
1
cos 2
2 , 2 ,
L
uu u
u p c u
au b u
      
  

        


 
      

  



   

 (14)
 
The general solution of Eq. (13) is sought in the form: 
     
0 1 2 2 1 0 2
, cos uA        

, (15) 
where both the amplitude A(τ 2) and phase Φ 0(τ 2), 
respectively, are modulated in dependence on the slow 
time scale and are currently undetermined. They will 
be determined in the next step of the perturbation 
procedure by elimination of so called secular terms in 
the solution u 1(τ 1,τ 2) of Eq. (14). Due to the modulated 
amplitude A(τ 2) and phase Φ 0(τ 2) in Eq. (15), 
oscillations of Duffing oscillator are aperiodic. By 
substitution of solution (15) into Eq. (13) and by 
considering the fact, that frequencies ω 0, ω L can take 
only positive values, the following important relation 
is obtained: 
0 L
  (16) 
By substitution of the general solution (15) into Eq. 
(14) and by considering relation (16) one has:
 
0
: 
1
: 

Anali PAZU - Letnik 3, leto 2013, številka 1 
  
19 
 
 
       
 
       
 
       
2
2 1
1 1 0 2 2 0 2 1 0 2 2 0 2
2
1
2 0 2
2 1 0 2 2 1 0 2
22
3
2 1 0 2 1 0 2
cos cos sin sin
2 sin cos
1
3cos cos 3 2
4
L
L
u
u p p
dA d
cA A
dd
A aA
              

  
        

        


                
       



   
       
       
   
          
     
     
                 
2 1 0 2
2 1 0 2 2 1 0 2 2 1 0 2
cos cos
2 sin sin 2 sin cos 2 cos sin
L
L L L L L
aA b A b A
    
                   


           
     
  (17)
The right hand side of Eq. (17) contains secular terms, 
which cause the unlimited growth of the solution. In 
order to obtain an uniform solution, secular terms must 
be eliminated by assembling terms, which appear at 
cos[τ 1-Φ 0(τ 2)] and sin[τ 1-Φ 0(τ 2)], respectively and by 
equating the obtained expressions with zero:
 
   
 
         
02 3
2 0 2 2 2 2 2
2
3
cos 2 2 cos 2 sin 0
4
L L L L
d
p A A aA b A
d

              

      

, (18) 
 
 
 
         
2
2 0 2 2 2 2
2
sin 2 2 sin 2 cos 0
L L L L
dA
p cA aA b A
d

            


       
 

.  (19) 
By introduction of the new variable: 
    
2 2 0 2
       ,  (20) 
Eqs. (18) and (19) are rewritten on the form: 
 
   
 
         
2 3
2 2 2 2 2 2
22
d3
cos 2 2 cos 2 sin 0
d4
L L L L
d
p A A aA b A
d
 
              


       
 

,      (21) 
 
 
 
         
2
2 2 2 2
2
sin 2 2 sin 2 cos 0
L L L L
dA
p cA aA b A
d

          


      
 

.  (22)
Equations (21) and (22) represent the system of 
ordinary nonlinear differential equations. Solution of 
both equations give us the course of amplitude and 
phase, respectively of nonstationary fundamental 
resonance of Duffing oscillator with time delay 
feedback control. Solution of the system of Eqs. (21) 
and (22) can be obtained by various methods of 
numerical integration, where often the Runge-Kutta 
method is used.  
5.1. Stationary fundamental resonance of Duffing 
oscillator with delay feedback control 
Stationary oscillations of Duffing oscillator appear, 
when A(τ 2)=A=const, Φ 0(τ 2)= Φ 0=const, γ(τ 2)= γ = 
const . Consequently, it holds 
 
2
2
0
dA
d


 ,
 
2
2
0
d
d


 , 
 
2
2
d
d
0


  and Eqs. (21), (22) are simplified to the 
nonlinear algebraic equations:
    
3
3
cos 2 2 cos 2 sin 0
4
L L L L
p A A aA b A               ,  (23) 
     sin 2 2 sin 2 cos 0
L L L L
p c A aA b A             . (24)
By squaring Eqs. (23) and (24) and adding, we get the 
following equation of the amplitude-detuning response 
of Duffing oscillator with time delay feedback control:
 
 
       
2 2
2
22
3
cos sin sin cos
84
L L L L L L L L
p
A a b c a b A              

 
       

 



.    (25) 
 
 
 

Anali PAZU - Letnik 3, leto 2013, številka 1 
  
20 
 
The corresponding course of the phase of Duffing 
oscillator is obtained, when Eq. (24) is divided by Eq. 
(23):
 
   
   
2
cos sin sin
tan
3
cos
cos sin
8
L L L L
L L L L
b c a p
p
A a b
      


       


  
. (26)
When the amplitude A and phase angle γ, 
respectively are determined, the steady-state solution 
of the zeroth order for displacements u(t) can be 
computed in the form: 
       cos u t A t O       . (27) 
From this solution we can see, that oscillations of 
Duffing oscillator are nearly harmonic (a more 
accurate analysis would be show, that oscillations are 
in fact periodic due to the contribution of higher 
harmonics in the solution u 1(τ 1,τ 2)). When the control 
of Duffing oscillator is not applied, that is, when 
0 ab  , then Eqs. (25), (26) reduce on 
relationships: 
  
2 2
2 22
3
84
LL
p
A c A    


   




, (28) 
 
2
tan
3
8
L
L
c
A


  


. (29) 
The obtained result describes the stationary 
fundamental resonance of Duffing oscillator that is in 
accordance with result, obtained by Nayfeh in Mook 
[5], which use multiple scales as well as with result of 
Leung in Fung [6], which apply the incremental 
harmonic balance method. On the basis of this 
comparison we can conclude that zeroth order 
approximation of nonlinear oscillations of Duffing 
oscillator by EL-PM with two time scales offers 
results, which are accurate enough.  
6. Analysis of nonstationary subharmonic 
resonance of Duffing oscillator with time delay 
feedback control 
Analysis of subharmonic resonances is performed 
by the similar procedure as in the case of fundamental 
resonance starting from Eq. (11.b). In this paper, we 
restrict oneself to the ⅓ subharmonic resonance, where 
the variable excitation frequency ω can be expressed 
by means of detuning σ as follows: 
 
0
3     . (30) 
In the case of ⅓ subharmonic resonance, Eq. (11.b) is 
rewritten on the form:
 
 
     
2 2 2
2 2 2 3
0 0 0 1 2
22
1 2 1 2
1 2
2
1 0 2 0 1 0 2 1 0 2
12
2 2 cos 3
2 , 2 , 2 ,
L
u u u u u
c u u p
au b u b u
          
   
 
               

     
       

   
  

     

            (31)
By substituting the power series (12) into Eq. (31) 
and by equating terms of like powers of ε on both sides 
of equation, one obtains the following system of linear 
partial differential equations (PDE's): 
  
 
2
22 0
0 0 1 2
2
1
cos 3
L
u
up     


  

,             (32) 
  
   
2 2
3 22 00 1
0 1 0 0 0 1 0 2 0 0 1 0 2
2
1 2 1 1
1
2 2 , 2 ,
L
uu u
u c u au b u             
   


 
         

   


,            (33) 
 
The general solution of Eq. (32) is equal to the sum of 
the general solution of the corresponding homogenous 
PDE and the particular solution
        
0 1 2 2 1 0 2 0 2 1
0
, cos cos u A E

       


   



, (34)
 
0
: 
1
: 

Anali PAZU - Letnik 3, leto 2013, številka 1 
  
21 
 
where    
2 1 0 2
cos A     

  describes the solution 
of the homogeneous equation and  
0 2 1
0
cos E






 
is the corresponding particular solution. When the 
solution    
2 1 0 2
cos A     

  is substituted into 
homogeneous equation, relation (16) is obtained again. 
Because the particular solution satisfies Eq. (32), too, 
we can compute the unknown amplitude E 0(τ 2) by 
substituting the particular solution into Eq. (32):
   
 
 
   
02
2 2 2 2 2
86 3
L L L LL
p p p
E 
        
   
   
.      (35)
By means of Eq. (34) we ascertain that the amplitude 
of particular solution E 0(τ 2) is modulated in the case of 
nonstationary resonance (due to the variability of 
excitation frequency ω), while it is constant in the case 
of stationary resonance. Solution (34) can be rewritten 
by considering Eq. (35) on the form:
      
 
 
0 1 2 2 1 0 2 1 2
2
, cos cos 3
86
LL
p
uA         
     
    


. (36) 
By substituting Eq. (36) into Eq. (33) and by considering relation (16) we get: 
 
 
     
 
 
 
 
 
     
 
2
2 0 2 22 1
1 1 0 2 2 1 0 2
2
22
1
0
1 2 1 2 2 1 0 2
2 2
2
2
2
dd
2 sin cos
dd
12 3 3d
cos 3 sin 3 2 sin
d
86
86
3
sin 3
86
L L L
L
LL
LL
LL
A u
uA
p p
cA
p
  
         


    
          

     
     

     
 
        
    
 

 
      
  

 
 

 




     
 
 
   
 
 
   
 
 
3
1 2 2 1 0 2 1 2
2
2 1 0 2 1 2
2
2 1 0 2 1 2
2
cos cos 3
86
2 cos cos 3 3
86
3
2 sin sin 3 3
86
LL
LL
LL
L L L
LL
p
A
p
aA
p
bA
         
     
          
     
           
     
  
  
     

  


      



      
 


 
  


      




 





      (37)
Equation (37) contains secular terms on its right hand 
side. By using a similar procedure as in the case of 
fundamental resonance, secular terms are eliminated, 
when terms, which appear at cos[τ 1-Φ 0(τ 2)] and sin[τ 1-
Φ 0(τ 2)], respectively, are collected and equated by 
zero:
 
 
 
 
         
2
22
2 2 0 2 2 2
2
2
d3
2 2 sin 3 2 sin 2 cos 0
d
4 8 6
L L L L L
LL
A A p
cA aA b A
  
             

     
      




 , (38)
 
 
 
 
 
 
 
   
2
0 2 2 2
2 2 0 2
2 2
2
2
d3 33
2 cos 3
d4
4 8 6
2 8 6
2 cos 2 sin 0
L
LL
LL
L L L
pA p
A
ab
    
      

   
   
    
     


 



  
.        (39) 
By introducing a new variable 
    
2 2 0 2
3        ,    (40) 
Eqs. (38) and (39) can be rewritten on the form: 
 
 
 
 
 
         
2
22
2 2 2 2
2
2
d3
sin sin cos
d
8 8 6
LL
L
L L L
A A p a
cA A bA
  
        

      
     




,             (41) 

Anali PAZU - Letnik 3, leto 2013, številka 1 
  
22 
 
 
 
 
 
 
 
   
2
22 2
2 2 2
2 2
2
22
d9 d 9 9 3
cos cos 3 sin
d d 8
8 8 6
4 8 6
LL
LL
L L L
L L L
pA pa
Ab
      
        
   
    
    
       


 



.  (42)
Equations (41) and (42) represents the amplitude 
response and phase response equation, respectively, of 
⅓ nonstationary subharmonic resonance of Duffing 
oscillator with time delay control. The corresponding 
amplitude response and phase response equation, 
respectively, of ⅓ stationary subharmonic resonance is 
obtained, when amplitude and phase does not change 
in dependence on slow time scale τ 2, that is, when 
A(τ 2)=A=const, Φ 0(τ 2)= Φ 0=const, γ(τ 2)=γ=const. 
Consequently, 
     
2 2 2 2 2 2
d /d 0, d /d 0, d /d 0 A            
and Eqs. (41) and (42) are reduced on the form:
 
 
     
2
2
3
sin sin cos 0
8 8 6
LL
L
L L L
A p a
cA A bA

    

    
   



,       (43) 
 
 
     
2
2
2 2
2
9 9 9 3
cos cos 3 sin 0
8
8 8 6
4 8 6
LL
LL
L L L
L L L
p pA a
Ab
  
     

    
    
     

 



.   (44)
By squaring Eqs. (43) and (44) and adding, we get the 
equation of the amplitude response of stationary 
subharmonic resonance with time delay control
   
 
   
 
2
2
2
22
2
2
2 2 4
2
22
9 9 3
9 sin cos cos 3 sin
8
4 8 6
81
64 8 6
L L L L
L L L
L L L
L L L
a p a
c b A b A
pA

        
  
      

      



 
      
 
   
  





.     (45) 
Equation (45) possesses a trivial solution A=0 and a 
nontrivial solution, which is obtained by solving the 
equation:
 
 
   
 
   
2
2
2
2
2
2
2 2 2
2
22
9 9 3
cos 3 sin
8
4 8 6
81 3
3 sin 3 cos 0
64 8 6
LL
LL
L L L
LL
L
L L L
pa
Ab
p A a
cb

    

    

   

    


   


 
 

    

 


.     (46) 
When Duffing oscillator is not controlled, that is, when 
a=b=0, then amplitude response of stationary 
subharmonic resonance is reduced on the form:
 
 
   
2
2 2 2 2
22
22
2 2 2
9 9 81
90
8
4 8 6 64 8 6
L
L L L L L L
p p A
Ac
  


         


    

   
    
    
.           (47)
7. Computational aspects 
The present method for both fundamental and ⅓ 
subharmonic resonance of Duffing oscillator with time 
delay control is implemented in programming 
environment Mathematica
®
 to apply symbolic as well 
as numeric computational capabilities. For example, 
stationary resonances (Eqs. (25), (28) for fundamental 
and (46), (47) for ⅓ subharmonic resonance, 
respectively) are computed algebraically. On contrary, 
the corresponding nonstationary resonances (Eq. (21), 
(22) for fundamental and Eq. (41), (42) for ⅓ 
subharmonic resonance, respectively) are computed by 

Anali PAZU - Letnik 3, leto 2013, številka 1 
  
23 
 
numerical integration, using embedded Runge-Kutta 
method. The question is, why not use the numerical 
integration directly on Duffing equation (7), however 
there are many reasons, why we not use this approach. 
The first reason is that we are looking for the steady-
state oscillations and thus we must integrate Eq. (7) 
each time until the transient phenomenon dies. 
Because we must repeat this long lasting computation 
for each particular value of detuning (and time delay 
as well), the approach is evidently impractical. This 
inefficiency becomes even more apparent in 
computation of nonstationary resonances, where 
detuning is slowly varied and a special strategy must 
be applied to achieve the elimination of the transient 
phenomenon. The second reason is, that numerical 
integration of Eq. (7) is able to compute stable 
solutions only. Therefore, unstable branches of 
resonance curves cannot be obtained by numerical 
integration of Eq. (7), however they are easily 
computed by the present method. If we limit only on 
stationary resonances with delay control, we can 
alternatively apply the incremental harmonic balance 
(IHB) method [3], which will be presented by authors 
for a wide kinds of dynamic structures in a 
forthcoming paper. The preliminary results shows a 
good agreement between EL-P and IHB method in the 
case of Duffing oscillator with delay control. 
8. Results  
Now look on results of EL-P analysis of 
fundamental and ⅓ subharmonic resonance of Duffing 
oscillator with time delay control. Both analysis are 
performed by using programming envinronment 
Mathematica
®
, where the following values of 
parameters c=0.05, ω L=1, β=0.05 are choosen in both 
cases. The excitation amplitude in the case of 
fundamental resonance is choosen to be equal p=0.5, 
while in the case of subharmonic resonance it is p=30. 
Analysis of stationary resonances in both cases are 
computed at first, while the corresponding analysis of 
nonstationary resonances follows to bring out 
important differences. Analysis of stationary 
resonance without control, that is, with values a=0, 
b=0 of gain parameters is performed in accordance 
with Eqs. (28),(29) in the case of fundamental 
resonance and by using Eq. (47) in the case of ⅓ 
subharmonic resonance. The stationary resonance of 
Duffing oscillator with time delay control is computed 
with values a=0.05, b= -0.05 of gain parameters and 
various time delays taking values θ=0, θ=π/4 and θ=π, 
respectively, where Eq. (25) is used in the case of 
fundamental and Eq. (46), respectively in the case of 
⅓ subharmonic resonance. Courses of stationary 
fundamental resonance of Duffing oscillator without 
control and with applied control are plotted on the 
Fig.2.a and corresponding nonstationary resonances 
on the Fig. 2.b. The variable excitation frequency ω 
and corresponding detuning σ are assumed to be linear 
functions of the slow time scale: 
                   
2 0 2
r      , (48) 
where σ 0 denotes the detuning in the time τ 2=0 and r 
denotes the rate of detuning change. In the analysis of 
both fundamental and ⅓ subharmonic resonance, 
respectively, the rate of detuning change r=±0.05 is 
used (the sign + corresponds to the passage with 
increasing of detuning (or increasing of excitation 
frequency) and the sign – corresponds to the passage 
with decreasing of detuning (or decreasing of the 
excitation frequency). From both Figures 2.a,b is 
apparent, that a small time delay (such as θ=π/4, for 
example) causes the reduction of the peak amplitude 
of the fundamental resonance curve, while a great time 
delay (θ=π in the study) causes the magnification of 
the peak amplitude. Thus for the case of fundamental 
resonance we can conclude, that the desired peak 
amplitude can be controlled by the proper choice of 
delay. However, the courses of nonstationary 
fundamental resonances show a drastic difference, 
which is characterized in change of amplitude and 
appearance of many amplitude peaks. After the first 
amplitude peak, heights of subsequent amplitude peaks 
are smaller and smaller. Stationary fundamental 
resonances of Duffing oscillator are ploted in Fig. 2.b 
for comparison. Response curves obtained for 
stationary and nonstationary ⅓ subharmonic 
resonances are ploted on Figures 3.a,b. Fig. 3.a shows, 
that delays in the case of stationary ⅓ subharmonic 
resonance of Duffing oscillator cause a different 
behavior of the controlled system in respect to the 
fundamental resonance. When the feedback control of 
Duffing oscillator is applied and the delay is equal 
zero, the response curve has a slightly larger amplitude 
as in the uncontrolled case. On contrary, when the 
feedback control of Duffing oscillator is applied and 
the delay has a great value, such as θ=π, the response 
curve has a slightly smaller amplitude as in the 
uncontrolled case. Therefore, changing of delay cannot 
be so effective in reduction of amplitude as in the case 
of fundamental resonance. Nonstationary ⅓ 
subharmonic resonances of Duffing oscillator show a 
drastic change in respect to corresponding stationary 
resonances, which are ploted in Fig. 3.b for 
comparison.
  

Anali PAZU - Letnik 3, leto 2013, številka 1 
  
24 
 
 
 a)  b) 
Fig 2. Fundamental resonance of Duffing oscillator with time delay control. a) stationary resonance 
without control (a=b=0) and with applied control at different values of delay θ. b) the corresponding 
nonstationary resonances with rate of detuning change r = ±0.05. 
 
 a)  b) 
Fig 3. ⅓ subharmonic resonance of Duffing oscillator with time delay control. a) stationary 
resonance without control (a=b=0) and with applied control at different values of delay θ. b) the 
corresponding nonstationary resonances with rate of detuning change r = ±0.05.
In both Figures 2.b and 3.b, nonstationary resonances 
are drawn with continuous lines when detuning 
increases and by broken lines, when detuning 
decreases. 
9. Conclusion 
In the paper it was shown, that EL-PM can be 
successfuly applied to the control of Duffing oscillator 
with time delay and therefore can be extended to 
problems of nonlinear theory of control. Moreover, it 
comes out, that EL-PM proves as useful in solving 
delayed nonlinear differential equations. The practical 
outcome of present analysis is that delays, which are 
troublesome in general, in fundamental resonance can 
be used with benefit, when the peak amplitude of 
resonance should be maintained in the prescribed 
range.  
References 
1. R. R. Pušenjak, Extended Lindstedt-Poincare 
method for non-stationary resonances of 
dynamical systems with cubic nonlinearities, J. 
Sound and Vib. 314, 194-216, 2008. 
2. R. Pušenjak, Razvejitve pri Van der Pol-
Duffingovem nihalu, Stroj. vestnik 49, št. 7/8, 
370-384, 2003. 
3. R. R. Pušenjak, M. M. Oblak. Incremental 
harmonic balance method with multiple time 
variables for dynamical systems with cubic non-
linearities, Int. j. numer. Methods eng., vol. 59, 
iss. 2, 255-292, 2004. 
4. A. Preymount, Twelve Lectures on Structural 
Dynamics, Springer, Dordrecht, 2013. 
5. A. H. Nayfeh, D. T. Mook, Nonlinear 
Oscillations, Wiley, New York, 1979. 
6. A. Y. T. Leung, T. C. Fung, Phase increment 
analysis of damped Duffing Oscillators, Int. J. 
Numer. Meth. Eng. 28, 193—209, 1989.