Approximating Limit Cycles of an Autonomous Delay Differential - - PowerPoint PPT Presentation

Approximating Limit Cycles of an Autonomous Delay Differential Equation

David E. Gilsinn Mathematical and Computational Sciences Division National Institute of Standards and Technology Gaithersburg, MD 20899-8910 dgilsinn@nist.gov

Acknowledgements: Dianne O’Leary, Computer Science Dept., Univ. of MD Tim Burns, MCSD, NIST

OUTLINE

• Introduction
• Main Existence Theorem
• Principal Parameter Estimation Algorithms
• Application to a Van der Pol Equation with Delay
• Final Observations

• Measurement and control of machine tool errors has

led to delay differential equation models. (Hanna, Tobias [4])

• Machine tool chatter has been established as a

Hopf bifurcation of limit cycles from stable machining. (Gilsinn [1], Nayfeh, et al. [5])

• Chatter is self sustained limit cycles caused by the

cutting tool interacting with undulations from a previous cut.

• Wish to approximate limit cycles with an analytic form

and develop a computable error bound.

INTRODUCTION

SOME BASIC FACTS

• An autonomous delay differential equation (DDE) with a fixed delay

will be written with initial condition from the space of continuous functions

• n [-h,0] .
• If satisfies a Lipschitz condition with respect to x , independent
• f y , there exists a unique solution of (1) for on [-h,0] .
• A linear DDE with fixed delay takes the form

with initial condition on [-h,0] , a column vector. (Forward Int.)

• The formal adjoint equation to (2) takes the form

for on , some initial interval, , a row vector. (Backward Int.)

• A solution of the linear DDE (2) is called a fundamental solution

if ( )

( ) ( ( ), ( )) 1 x t X x t x t h = −

• φ

( , ) X x y

φ

( )

( ) ( ) ( ) ( ) ( ) 2 x t A t x t B t x t h = + −

• φ

n

x R ∈ ( )

( ) ( ) ( ) ( ) ( ) 3 y t y t A t y t h B t h = − − + +

• ψ

[ , ] t t h +

n

R ψ ∈ ( , ) Z t s

( , ) , ( , ) 0, . Z s s I Z t s t s = = <

OBJECTIVE

• Find periodic solution of autonomous delay differential

equation (DDE)

( )

( ) ( ( ), ( )) 4 x t X x t x t h = −

• with initial condition from the space of continuous

functions on

• Period is also unknown.
• Introduce for t to get
• Look for periodic solutions of fixed period
• Desired Result: Given an approximate -periodic

solution and frequency, , of (2), wish to show that if they satisfy a certain noncriticality condition then there exists an exact frequency and -periodic solution, , in a computable neighborhood of 2 / T π ω =

/ t ω

φ

[ ,0] h −

( )

( ) ( ( ), ( )) 5 x t X x t x t h ω ω = −

• 2π

2π

( , ( )) x t ω

2π

( *, *( )) x t ω

( , ( )) x t ω

NOTATION

1 1 2 2 1 2

( ( ), ( ) ( ; ) ( ) ( ) , ( , ) x x t x t h dX x X x X x

ω ω ω ω

ω φ φ φ φ φ φ = − = + =

• If are an approximate frequency and - periodic solution then
• The variational equation with respect to the approximate solution is

Let A(t) = , B(t) =

• is a characteristic multiplier of the linear system (2)

if there exists a non-trivial solution of (2) such that

( ) ( ) ( ), ( 2 ) ( ) (6) x t X x k t k t k t

ω

ω π = + + =

• ( )

( ; ) (7) z t dX x z

ω ω

ω =

• ( , )

x ω 2π

ρ

( 2 ) ( ) (8) x t x t π ρ + =

1(

) X xω

2(

) X xω

NONCRITCAL APPROXIMATE SOLUTION

• The pair is noncritical with respect to if the

variational equation about the approximate solution (7) has a characteristic multiplier of multiplicity one with the remaining multipliers unequal to one. If is the periodic solution

• f the adjoint corresponding to then

where (Hale [3], Stokes [7])

( , ) x ω

( ) x X xω ω =

• ρ

2

( ), 1 v t v =

2π

ρ

2

( ) ( , )( ) (13) v t J x t dt

π

ω ≠

• 2

,

( , )( ) ( ) (14)

t t

J x t x hX x x

ω ω

ω ω = +

IMPORTANT LEMMAS

LEMMA (Halanay [2]): If is noncritical, -periodic, such that then there exists a unique -periodic solution of which satisfies for some , independent of

, 2 f π

( , ) x ω

2

( ) ( ) (11) v t f t dt

π

=

• 2π

( ) ( ; ) ( ) (12) z t dX x z f t

ω ω

ω = +

• 2

z M f ≤

M >

f

LEMMA (Halanay [2]): When the linear DDE coefficients A(t), B(t) are periodic the linear and adjoint systems have the same finite number of independent solutions. LEMMA (Hale [3], Halanay [2]): If is a simple characteristic multiplier of the linear DDE (2), p(t) a nontrivial -periodic solution

• f the linear DDE (2), q(t) a nontrivial - periodic solution of the

adjoint (3) and then

( , )( ) ( ) ( ) ( ) (9) J p t p t hB t p t h ω ω ω = + − ( )

2

( ) ( , )( ) 10 q t J p t dt

π

ω ≠

• ρ

2π 2π

FUNDAMENTAL THEOREM

THEOREM (Stokes [7]): Let satisfy and let . Suppose there exist and such that for Assume is noncritical (in the delay sense) and let be the appropriate solution of the adjoint to the variational equation, . Let If M is the constant from the previous lemma, let Finally, if there is a function C of computable parameters such that then there exists an exact -periodic solution and an exact frequency so that

( , ) x ω

( ) ( ) ( ) (15) x t X x k t

ω

ω = +

• k

r ≤

1

K

K

1 2

, ψ ψ

1 1 2 1 2

( ; ) (16) ( ; ) ( ; ) dX x K dX x dX x K

ω ω ω ω ω ω

φ φ ψ φ ψ φ ψ ψ φ ≤ + − + ≤ −

( , ) x ω

v

1 v =

1 2

1 ( ) ( , )( ) (17) 2 v t J x t dt

π ω

α ω π

−

• =
• 1

1 2 1

(1 ( , ) ) (18) (1 ) M J x MK M λ α ω λ λ = +

• =

+

• 1

1 2

( , , , , , , , , ) 1 (19) rC K K h x x α λ λ ω <

2π

*

x

*

ω

* 1 *

4 2 (20) x x r r λ ω ω α − ≤ − ≤

PROOF (OUTLINE)

Goal: Find -periodic so that is and exact solution of Substituting (21) into (22) where and is a function of computable parameters and Strategy: Wish to find a fixed point of a map such that the perturbation term on the right of (23) is orthogonal to the solution of the adjoint in the noncriticality definition. The Lemma can then be used To solve for a

, ( ) 2 z t β π

, ( ) ( ) ( ) (21) x t x t z t ω ω ω β ω β = + = + +

( ) ( ( ), ( )) (22) x t X x t x t h ω ω = −

• ( )

( ; ) ( , ) ( , )( ) ( ) (23) z t dX x z R z J x t k t

ω ω

ω β β ω = + − −

• (

) ( ) (24) x X x k t

ω

ω = +

• ( ,

) R z β ( , ). z β

( ( ), ( )) g R z g g β =

( ) z t

PROOF (CONTINUED)

Define the sets Construct map as a composition Define: Given solve for unique so that solution of adjoint for noncritical. By Lemma there exists a unique such that Now define Define: Given Define:

{ }

{ }

1 1

| , ( 2 ) ( ) | P g g C g t g t P P C N g P g π δ = ∈ + = ⊂ ∩ = ∈ ≤ S

( ) S L T = β >

( , )( ) ( ) ( ) g J x t k t t β ω ν − − ⊥

• ( )

t ν ( , ) x ω

1

z P ∈

2

( ) ( ; ) ( ) ( , )( ) ( ) ( , ) z t dX x z g t J x t k t z M g J x k

ω ω

ω β ω β ω = + − − ≤ − −

• ( )

( ( ), ( )) L g g z g β =

1

: L N R P → ×

g N ∈

1

: T R P N × →

1

0, ( , ) , z R P β β > ∈ × ( , ) ( , ) T z R z β β =

: , ( ) ( ( ), ( )) S N N S g R z g g β → =

PROOF (CONTINUED)

For in the definition of N there exists a bounded function Such that , which implies is a contraction. Therefore there exists a fixed point The exact solution is then given by Finally, we can show that NOTE: This provides only O(r) estimates. These may not be optimal bounds but they are computable.

r δ =

1 1 2

( , , , , , , , ) C K K h x α λ λ ω

• 1

rC <

: S N N →

* * *

( ( ), ( )) g R z g g β =

* * * * *

( ) ( ) ( ) ( )( ) ( ) g x t x t z g t g ω ω β ω ω β = + = + +

* 1 *

4 2 x x r r λ ω ω α − ≤ − ≤

APPLICATION STEPS

• Compute the approximation pair
• Verify that the pair is noncritical
• Compute M and

( , ) x ω

α

A quote from Stokes [7] “The computational difficulties here are considerably greater than in the case of ordinary differential equations,…, but they are not insurmountable” He never produced an example. This talk describes the first application.

APPROXIMATE SOLUTION

• Develop solution as a trigonometric polynomial
• Set coefficient of one term to zero, say sin(t), in order

to estimate

• Can develop Galerkin projection equations using, e.g. MAPLE,

although the expansions are nontrivial (typical over 135 terms). Not recommended in general.

• Summer student, Chris Copeland, and I have developed a

a fast projection procedure in MATLAB based on some FFT ideas. Subject for another talk.

( )

1 2 2 1 2 2

( ) cos( ) cos( ) sin( )

n k k k

x t a a t a kt a kt

− =

= + + +

• ω

VARIATION OF CONSTANTS FORMULAS

• From now on we revert to classic notation where we set
• and are -periodic.
• The variation of constants formula for the linear system

is ( solution of linear system with ) ,

• h has been normalized to 1.
• The variation of constants formula for the adjoint

is

• These formulas are developed in Halanay [2]. Significance of

the adjoint formula is that it only requires a forward integration.

1 2

( ) ( ), ( ) ( ) A t X x B t X x

ω ω

= = ( ) A t ( ) B t 2π

( ,0) Z t (0,0) , ( ,0) 0, Z I Z t t = = <

( ) ( ,0) (0) ( , ) ( ) ( ) z t Z t Z t B d

ω

φ α ω α ω φ α α

−

= + + +

• [

,0]

• C

φ ω ∈ −

( ) ( ) ( ) ( ) ( ) x t A t x t B t X t ω = + −

• ( )

( ) ( ) ( ) ( ) y t y t A t y t B t ω ω = − − + +

• 2

2

( ) (2 ) (2 , ) ( ) ( ) ( , ) y t Z t s B s Z s t ds

π ω π

ψ π π ψ ω

+

= + −

• [2 ,2

] ψ π π ω ∈ +

TESTING NONCRITICALITY CONDITION

• Test characteristic multiplier of multiplicity one and compute
• is a characteristic multiplier of the linear system if there is a

solution such that

• Halanay [2] shows that the eigenvalues of the following operator

are the multipliers of the variational equation for the linear system

• Called the monodromy operator, defined formally by

where

• U is compact with at most a countable number of eigenvalues

with 0 the only possible limit point.

• In the present case will be

α

ρ

( ) z t

( 2 ) ( ) z t z t π ρ + =

( )( ) ( 2 ,0) (0) ( 2 , ) ( ) ( ) , U s Z s Z s B d

ω

φ π φ π α ω α ω φ α α

−

= + + + + +

• (

)( ) ( 2 ; ) U t x t φ π φ = +

0[

,0]. C φ σ ∈ −

[ ,0] s σ ∈ −

σ −

ω −

MONODROMY OPERATOR MAP ON INITIAL SPACE

ω − 2π ω − 2π

Uφ

φ

[ , 0] C ω −

Translate back

Trajectory Plot Interpreted as a Mapping on Initial Space

DISCRETIZATION OF U

• The eigenvalues of are approximated by the eigenvalues
• f where is constructed as follows:
• Discretize

with equal intervals where

• Approximate the integral operator by the trapezoidal rule (other rules

could be used) by setting where satisfies the linear variational system about the approximate periodic solution,

Uφ ρφ =

N

U φ ρφ =

N

U

[ ,0] ω −

1 2 1 N

s s s ω

+

− = < < < =

• 1

/ .

i i

h s s N ω

+

= − =

1 1

, , 2, , 2

N i

h w w w h i N

+

= = = =

• 1

1

( )( ) ( 2 ,0) (0) ( 2 , ) ( ) ( )

N N j j j j j

U t Z t w Z t s B s s φ π φ π ω ω φ

+ =

= + + + + +

• ( ,0)

Z t

(0,0) , ( ,0) 0, 0. Z I Z t t = = <

Block Matrix for

N

U

• The matrix for becomes
• It is not necessary to compute for all

. For large N, say 1000 or more, this would be somewhat impractical computationally.

• Note that and

for

• For column j compute but save intermediate

integration points as is done in dde23 in Matlab.

• Interpolate the values of up the column for
• This also applies to
• This reduces the block integrations of Z to N+1 instead of (N+1)2.

N

U

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

( 2 , ) ( ) ( 2 , ) ( ) ( 2 , ) ( 2 , ) ( ) ( 2 , ) ( ) ( 2 , ) ( ) ( 2 , ) ( 2 , ) ( ) ( 2 , ) (

j j j N j N N i j i j j i N j i N N N

w Z s s B s w Z s s B s Z s S w Z s s B s w Z s s B s w Z s s B s Z s S w Z s s B s w Z s s B s π ω ω π ω ω π π ω ω π ω ω π ω ω π π ω ω π ω

+ + + + + + +

+ + + + + + + + + + + + + + + + + + + + + + + +

• 1

1 1 1 1 1 1

) ( 2 , ) ( ) ( 2 , ) ( 2 , ) ( )

j N j j N N j N N N

w Z s s B s Z s S w Z s s B s ω π ω ω π π ω ω

+ + + + + +

• +

+ + + + + + + +

• (

2 , )

i j

Z s s π ω + +

, 1, , 1 i j N = +

• 1

2 2

N

s π π

+ +

=

2 2

i

s π π < + <

1, , . i N =

(2 , )

j

Z s π ω + ( 2 , )

i j

Z s s π ω + +

1, , . i N =

( 2 ,0).

i

Z s π +

Computing the parameter - 1 α

• Need to solve adjoint equation for on

[0,2 ], t π ψ ∈ [2 ,2 ] π π ω + ( ) ( ) ( ) ( ) y t y t A y t B t ω ω = − − + +

• Halanay [2] showed that the solution satisfies

( ) y t

2 2

( ) (2 ) (2 , ) ( ) ( ) ( , ) y t Z t s B s Z s t ds

π ω π

ψ π π ψ ω

+

= + −

• Significance: Requires only forward integration to compute

( , ) Z s t

• To find: Solution of the adjoint associated with multiplier of variational eq.
• For on , a row vector, define

φ

• [

,0] ω −

( )

( ) ( ) (2 , ) ( ) ( ) (2 , ) U s Z s B Z s d

ω

φ φ ω π ω φ η η ω π η ω η

−

= − + + + + +

• Halanay [2] defines an associated operator for

[2 ,2 ] s π π ω ∈ +

( )

2 2

( ) ( 2 ; ) (2 ) (2 , 2 ) ( ) ( ) ( , 2 ) V s y s Z s B Z s d

π ω π

ψ π ψ ψ π π π ψ η η η ω π η

+

= − = − + − −

• and showed that an eigenvalue
• f is associated with a

multiplier of the adjoint.

ρ

V

• 1/ ρ

Computing the parameter - 2 α

• Halanay [2] showed that the eigenvalues of are all the same

and that the eigenvectors of are related by

, , U U V

, U V

( ) ( 2 ), [ ,0] s s s φ ψ π ω ω = + + ∈ −

• To solve for the solution of the adjoint in row form on

we need only compute the significant eigenvalue and eigenvector of

[0,2 ] π

U

• Discretizing

the j-th column is

1 1

0, /

N

s s N ω ω

+

− = < < = ∆ =

• 1

1 1 1 1 1

(2 , ) ( ) ( 2 , ) ( ) ( 2 , ) ( )( ) [ ( ), , ( ), , ( )] ( ) ( 2 , )

j j j j i i j j i N j N N j

Z s w B s Z s s w B s Z s s U s s s s w B s Z s s π ω ω π ω ω π ω φ φ φ φ ω π ω

+ + +

+ + ∆ + + +

• ∆

+ + +

• =
• ∆

+ + +

• Then for

* 1 1 * 1 1 * 1 1

(2 , ) ( 2 ) ( 2 , ) ( ) [ ( ), , ( ), , ( )] ( 2 ) ( 2 , ) ( 2 ) ( 2 , )

j j j j i N j i i j j N N j

Z t w B s Z s t y t s s s w B s Z s t w B s Z s t π π ω π φ φ φ π ω π π ω π

+ + +

• +

∆ + + +

• =

∆ + + +

• ∆

+ + +

• *

1 1

2 , 2 /

P

t t P π π

+

= < < = ∆ =

• This is the discretized form of the variation of constants formula for

the adjoint equation

Estimating M such that

2

z M f ≤

• Solution of nonhomogeneous system,

is

[0,2 ] t π ∈

( ) ( ) ( ) ( ) ( ) z t Az t B t z t f t ω = + − +

• ( )

( ,0) (0) ( , ) ( ) ( ) ( , ) ( ) (25)

t

z t Z t Z t s B s s ds Z t s f s ds

ω

φ ω ω φ

−

= + + + +

• periodic initial condition for

2π

[ ,0] s ω ∈ −

2

( ) ( 2 ,0) (0) ( 2 , ) ( ) ( 2 , ) ( ) (26)

s

s Z s Z s d Z s f d

π ω

φ π φ π α ω φ α α π α α α

+ −

= + + + + + +

• Solve for as

(0) φ

2

(0) ( (2 ,0)) (2 , ) ( ) ( ) ( (2 ,0)) (2 , ) ( ) (27) I Z Z B d I Z Z f d

π ω

φ π π α ω α ω φ α α π π α α α

+ + −

= − + + + −

• Insert (27) into (26), use sup norm and Schwarz inequality to get

1 1 2 2 1 2 2

, (1 ) m m f m m f φ φ φ

−

≤ + ≤ −

• Insert (27) into (25) to get

2

( ) [ ( ,0)( (2 ,0)) (2 , ) ( , )] ( ) ( ) [ ( ,0)( (2 ,0)) (2 , ) ( , )] ( ) z t Z t I Z Z Z t B d Z t I Z Z Z t f d

π ω

π π α ω α ω α ω φ α α π π α α α α

+ + −

= − + + + + + − +

• Again use sup norm and Schwarz inequality to get

1 3 4 2 3 1 4 2 2 2

, [ (1 ) ] z m m f z m m m m f M f φ

−

≤ + ≤ − + =

Van der Pol Equation with Unit Delay

2

( ( 1) 1) ( 1) x x t x t x λ + − − − + =

• Introduce unknown frequency by substituting for to get

2 2

( ( ) 1) ( ) x x t x t x ω ωλ ω ω + − − − + =

• /

t ω

• In vector form with , initial condition on
• The variational equation can be written as
• For this example we take

1 2

( ) ( ), ( ) ( ) x t x t x t x t = =

1 1 1 2 2 2 2 2 1 2

( ) 1 ( ) ( / ) ( ) ( ) 1/ / x x x t d x x x t x t x t dt ω ω λ ω ω ω ω λ ω −

• =

+ +

• −

− − − −

• [

,0] ω −

t ( )

2 2 1 2 1

1 ( ) , ( ) , ( 2 ) ( ) 2( / ) ( ) ( ) ( / ) 1 ( ) 1/ A t B t B t B t x t x t x t π λ ω ω ω λ ω ω ω

• =

= + =

• −

− − − − −

• ( )

( ) ( ) ( ) z t Az t B t Z t ω = + −

• 0.1

λ =

Developing an Approximate Solution

• Selected a 7 harmonic expansion and in the delay equation

( )

7 1 2 2 1 2 2

( ) cos( ) cos( ) sin( )

k k k

x t a a t a kt a kt

− =

= + + +

• The sin(t) term was dropped in order to estimate .

ω

• 15 Galerkin projection equations developed using MAPLE.
• Solving the projection equations produced

and

0.1 λ =

2 5 6 9 10 13 14

2.0185, 2.5771 3, 2.5655 2, 1.0667 4 5.2531 4, 7.1791 6, 2.2042 6 1.0012 a a e a e a e a e a e a e ω = = − = − = − = − − = − − = − − =

• Using these and MAPLE to produce an expansion of the Van der

Pol equation and then taking the sup norm gave residual r = 6.0867e-6

1 3 4 7 8 11 12

a a a a a a a = = = = = = =

Estimating the Bound and Lipschitz constant K

1

K

• Use the fact that for matrix product with then

to show that

Bx

1

max

i i n

x x

≤ ≤

=

1 1

max

n ij i n j

B b

≤ ≤ =

=

• (

)

( )

2 2 1 2 1 2 14 7 14 2 2 2 1 2 1 2 1

1 ( ; ) (1/ ) (2 / ) ( ) ( ) ( / ) 1 ( ) max 1, 1/ (2 / ) ( / ) 1

i i i i i i i

dX x x t x t x t a a i a a a φ φ ω λ ω ω ω λ ω ω ω λ ω λ ω φ

− = = =

• ≤
• −

− − − − −

• ≤

+ + + + +

• If then

0.1, λ =

• Also, working within the domain

[ ]

{ }

0,2 : 1 D x C x x π = ∈ − ≤ [ ][ ] [ ][ ]

{

[ ] [ ] }

( )

1 2 1 11 2 21 1 12 2 22 2 2 1 11 1 12 1 2

( ; ) ( ; ) (2 / ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( / ) ( ) ( ) ( ) ( ) (6 / ) 1 dX x dX x x t t x t t x t t x t t x t t x t t x

ω ω ω ω

ψ φ ψ φ λ ω ω ψ ω ψ ω ψ ω ψ λ ω ω ψ ω ψ φ λ ω ψ ψ φ + − + ≤ − + − + − − + − + + − + − − + ≤ + −

• On

If then

[0,2 ], 2.0225 x π ≤

0.1, λ =

1

2.3776. K = 1.8113. K =

Van der Pol Equation with Unit Delay

• Using the characteristic multiplier algorithm and integrating

Z(t,0), etc. using dde23 in MATLAB get

• M is estimated as 1.7411
• is estimated as 8.0665
• Other estimates

0.99806 ρ =

α

1 2

2.0225 2.0258 2.1192 ( , 2.6082 38.3722 113.2727 x x x J x ω λ λ = = = = = =

Calculation Data Flow , , , 1.0012 x x x ω =

( , ) 2.6082 J xω ω =

• 0.99806

ρ =

1.7411 M =

1

38.3722 λ =

8.0665 α =

6.0867 6 r e = −

1

2.3776 K =

2

113.2727 λ =

1.8113 K =

1 1 2

( , , , , , , , , ) 0.3299 1 rC K K h x x α λ λ ω = <

* 1

4 9.3424 4 x x r e λ − ≤ = −

*

2 9.8197 5 r e ω ω α − ≤ = −

Final Observations

• Galerkin projection symbolic calculations are very lengthy.
• Lack of general form for required numerical integration by dde23.
• Spline interpolation reduced the number of times

had to be computed.

• Integration rule possibly led to fine discretization of intervals.
• No reasonably computable Green’s Function required numerical

estimation of M .

• Numerical Procedures led to interesting algorithmic results as also

demonstrated by Urabe and Reiter [9].

• Further study of convergence questions needed.
• Recent studies show that collocation might be more efficient for

solving for characteristic multipliers.

( , ) Z s t ( , ) Z s t

References

[1] Gilsinn, D. E., “Estimating Critical Hopf Bifurcation Parameters for a Second-Order Delay Differential Equation with Application to Machine Tool Chatter”, Nonlinear Dynamics, 30, 2002, 103-154. [2] Halanay, S., Differential Equations: Stability, Oscillations, Time Lags, Academic Press, New York, 1966. [3] Hale, J, Functional Differential Equations, Springer-Verlag, New York, 1971. [4] Hanna, N. H., Tobias, S. A., “A Theory of Nonlinear Regenerative Chatter”, ASME Journal of Engineering for Industry, 96, 1974, 247-255, [5] Nayfeh, S. H., Chin, C., Pratt, J., “Perturbation Methods in Nonlinear Dynamics – Applications to Machining Dynamics”, ASME Journal of Manufacturing Science and Engineering, 119, 1997, 485-493. [6] Stokes, A., “On the Approximation of Nonlinear Oscillations”, Journal of Differential Equations, 12, 3, 1972, 535-558. [7] Stokes, A. P., “On the Existence of Periodic Solutions of Functional Differential Equations”, Journal of Mathematical Analysis and Applications, 54, 1976, 634-652. [8] Urabe, M., “Galerkin’s Procedures for Nonlinear Periodic Systems”,

• Arch. Rational Mech. Anal., 20, 1965, 120-152.

[9] Urabe, M., Reiter, A., “Numerical Computation of Nonlinear Forced Oscillations by Galerkin’s Procedure”, Journal of Mathematical Analysis and Applications, 14, 1966, 107-140.