APPROXIMATING PERIODIC SOLUTIONS OF AUTONOMOUS DELAY

DIFFERENTIAL EQUATIONS

DAVID E. GILSINN

Mathematical and Computational Sciences Division

National Institute of Standards and Technology

100 Bureau Drive, Stop 8910

Gaithersburg, MD 20899-8910 1

ABSTRACT. Machine tool chatter has been characterized as isolated periodic solutions or limit cycles of

delay differential equations. Determining the amplitude and frequency of the limit cycle is sometimes crucial

to understanding and controlling the stability of machining operations. In Gilsinn [9] a result was proven

that says that, given an approximate periodic solution and frequency of an autonomous delay differential

equation that satisfies a certain non-criticality condition, there is an exact periodic solution and frequency in

a computable neighborhood of the approximate solution and frequency. The proof required the estimation of

a number of parameters and the verification of three inequalities. In this paper the details of the algorithms

will be given for estimating the parameters required to verify the inequalities and to compute the final

approximation errors. An application will be given to a Van der Pol oscillator with delay in the nonlinear

terms. A MATLAB m-file implementing the algorithms discussed in the paper is given in the appendix.

AMS (MOS) Subject Classification. 34K11, 34K13, 34K28.

1. Introduction

Machine tool dynamics has been modeled using delay differential equations for a number

of years as is clear from the vast literature associated with it. For a detailed review of

machining dynamics see Tlusty [32]. For a discussion of dynamics in milling operations see

Balchandran [1] and Zhao and Balachandran [34]. For drilling operations see Stone and

Askari [29] and Stone and Campbell [30]. For an analysis of chatter occurring in turning

operations see Hanna and Tobias [16], Marsh et al. [21], and Nayfeh et al. [22]. Machine

tool chatter is undesirable self-exited periodic oscillations during machining operations. It

has been identified as a Hopf bifurcation of limit cycles from steady state solutions. For a

way of estimating the critical Hopf bifurcation parameters that lead to machine tool chatter

see Gilsinn [8].

In studying the effects of chatter it is sometimes desirable to compute the amplitude and

frequency of the limit cycle generating the chatter. This entails solving the delay differential

equations that model the machine tool dynamics. There is a large literature on numerically

solving delay differential equations. Some representative methods are described in Banks

and Kappel [2], Engelborghs et al. [5], Kemper [20], Paul [23], Shampine and Thompson

[25], and Will ́e and Baker [35]. Although these methods generate solution vectors that can

be studied by harmonic and power spectral methods to estimate the frequency of periodic

cycles, they do not directly generate a representative model of a limit cycle such as a Fourier

series representation.

It is also desirable to know whether a representation of an approximate limit cycle is

close to a true limit cycle. In other words we wish to answer the question as to whether the

approximate solution represents sufficiently well a true solution. This is answered with a set

1Contribution of the National Institute of Standards and Technology, a Federal agency, not subject to

copyright.

1
2 DAVID E. GILSINN

of test criteria by Gilsinn [9], who showed that, given a representative approximate solution

and frequency for a periodic solution to the autonomous delay differential equation

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

where x, X ∈ Rn, h > 0, X sufficiently differentiable, there are conditions, depending on

a number of parameters, for which (1) has a unique exact periodic solution and frequency

in a computable neighborhood of the approximate solution and frequency. This result was

first established in a very general manner for functional differential equations by Stokes [28]

who extended an earlier result for ordinary differrential equations in Stokes [27]. A crucial

aspect in applying the result involves verify a certain ”non-criticality” condition. However,

no computable algorithms were given in the case of functional differential equations to

estimate the various parameters or verify the ”non-criticality” condition. Only recently

have algorithms been developed to computationally verify these conditions in the fixed

delay case. A preliminary announcement of algorithms for computing these parameters and

verifying the ”non-criticality” condition was given by Gilsinn [7]. In this paper we include a

more detailed discussion of the algorithms and apply them to a Van der Pol equation with

delay in its nonlinear terms.

The result of Stokes [28] for functional differential equations depends on verifying certain

conditions that require computing various parameters. In order to apply Stokes’ result

a proof in the case of equation (1) will be given here since certain inequalities that are

developed within the proof are necessary for proving the fixed point contraction mapping

conditions and rely on the specific form of (1).

The notation used in the paper is described in Section 2. The non-criticality condition is

defined in Section 3. In Section 4 we construct an exact frequency and 2π-periodic solution

of (1) as a perturbation problem. In Section 5 we define a map that is used to prove, by a

contraction argument, the existence of an exact frequency and 2π-periodic solution of (1).

The main contraction theorem is proven in Section 6. A Galerkin algorithm to compute a 2πperiodic solution to a nonlinear autonomous delay differential equation is given in Section 7.

The Floquet Theory for delay differential equations is discussed in Section 8. An algorithm

for computing the characteristic multipliers of the variational equation of (1) with respect to

the approximate 2π-periodic solution, is outlined in Section 9. An algorithm to determe the

solution to the formal adjoint equation with respect to the variational equation of (1) with

respect to the approximate 2π-periodic solution, is outlined in Section 10. An algorithm for

estimating a critical parameter is given in Section 11. An application of these algorithms to

the Van der Pol equation with delay is given in Section 12. Conclusions are given in Section

13 and a disclaimer is given in Section 14. The derivation of the differentiation matrix (86)

is given in Appendix 1, Section 15. Certain bounds and Lipschitz conditions used in the

fixed point theorem are proven in Appendices 2 and 3 (Section 16 and Section 17). The

MATLAB code and associated support functions implementing the algorithms are given in

Appendix 4, Section 18.

2. Notation

Let Cω denote the space of continuous functions from [−ω, 0] to Cn with norm in Cω

given by |φ| = max |φ(s)| for −ω ≤ s ≤ 0, where

(2) |φ(s)| =

X3n

i=1

|φi(s)|

2

~1/2

.

Cω is a Banach space with respect to this norm. Let P be the space of continuous 2π-periodic

functions with sup norm, |·| on (−∞, ∞). Let P1 ⊂ P be the subspace of continuously differentiable 2π-periodic functions with the sup norm. Let X(x, y) be continuously differentiable

in some domain Ωn ⊂ Cn × Cn with bounded derivatives where

(3) |Xi(x, y)| ≤ B,
STABILITY 3

for i = 1, 2, (x, y) ∈ Ωn. The subscripts of X indicate derivatives with respect to the first

and second variables of X respectively. We further assume that the first partial derivatives

satisfy Lipschitz conditions given by

(4) |Xi(x1, y1) − Xi(x2, y2)| ≤ K(|x1 − x2| + |y1 − y2|),

for (x1, y1), (x2, y2) ∈ Ωn.

In order to simplify the notation for (1) we will first normalize the delay h to unity by

setting s = t/h. Then, (1) becomes

(5) dy

ds (s) = h X(y(s), y(s − 1))

where y(s) = x(s h). Therefore we will assume h = 1 in (1). We will also make one further

transformation. Since the period T = 2π/ω of a periodic solution for (1) is unknown we can

normalize the period to [0, 2π] by introducing the substitution of t/ω for t and rewriting

(1), with h = 1, in the form

(6) ω x ̇ = X(x(t), x(t − ω)).

For ψ1,ψ2 ∈ P we denote the total derivative of X(x, y) by

(7) dX(x, y;ψ1, ψ2) = X1(x, y)ψ1 + X2(x, y)ψ2.

Let A(t), B(t) be continuous 2π-periodic matrices. Then a characteristic multiplier is

defined as follows.

Definition 2.1. ρ is a characteristic multiplier of

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

if there is a nontrivial solution y(t) of (8) such that y(t + 2π) = ρ y(t). Note that if ρ = 1

then y(t) is 2π-periodic.

To simplify some of the notation we will suppress the t and write, for example, x =

x(t), xω = x(t − ω), but in other cases we will maintain the t, especially when describing

computational steps. We will also at times use the notation

(9) |x|2 =

}8 2π

0

|x(t)|

2 dt]1/2

.

3. Non-criticality Condition

Galerkin and harmonic balance methods can be used to develop 2π-periodic approximate

solutions for (6). A fast discrete Fourier series algorithm for computing an approximate

series solution and frequency, (ˆω, xˆ), has been given by Gilsinn [8]. See Section 7 below for

a brief discussion of a Galerkin method for approximating a solution. At this point, then,

we assume that we have developed a 2π-periodic approximate solution and frequency, (ˆω, xˆ)

for (6), where ˆx is 2π-periodic and

(10) ˆω  ̇xˆ = X(ˆx, xˆωˆ ) + k

where k(t) is a 2π-periodic residual bounded by

(11) |k| ≤ r.

The required size of the residual error, r, will become clear based upon estimates later in this

paper. These estimates will indicate in particular situations how accurately an approximate

solution and frequency would need to be computed.

The variational equation with respect to the approximate solution and frequency is given

by

(12) ˆωz ̇ = dX (ˆx, xˆωˆ ; z, zωˆ ).

Let Aˆ = X1 (ˆx, xˆωˆ ), Bˆ = X2 (ˆx, xˆωˆ ). The formal adjoint of (12) is given in row form by

(13) ˆωv ̇ = −vAˆ − v−ωˆB. ˆ
4 DAVID E. GILSINN

The next lemma, proven in Halanay [14], relates the number of independent 2π-periodic

solutions of (12) to those of (13).

Lemma 3.1. System (12) and (13) have the same finite number of independent 2π-periodic

solutions.

We will not give the proof of the next lemma, since it is also proven in Halanay [14]. The

result, however, will be critical to the main approximation theorem.

Lemma 3.2. The nonhomogeneous system

(14) ˆωx ̇ = Axˆ + Bxˆ ωˆ + f

has a unique 2π-periodic solution if and only if

(15) 8 2π

0

vT

0 f dt = 0

for all independent solutions v0 of period 2π of (13). Furthermore there exists an M > 0,

independent of f, such that

(16) |x| ≤ M|f |.

We will give the proof of the next lemma. Although it is stated in Hale [15] and in Halanay

[14], the proof is not generally available. The result, however, motivates the definition of a

non-critical approximate solution.

Lemma 3.3. Let ρ = 1 be a simple (i.e. multiplicity one) characteristic multiplier of (12)

Let p be a non-trivial solution of (12) associated with ρ. Define

(17) J(p,ωˆ) = p + Bpˆ ωˆ ,

then

(18) 8 2π

0

v0

T J(p, ωˆ)dt W= 0

for all independent solutions v of the adjoint (13).

Proof. Let y(t) be any 2π-periodic solution of (12) and write

(19) z = y + tp

Then, substituting (19) into (12), write

ωˆz ̇ = Azˆ + Bzˆ ωˆ + ˆω

p

p + Bpˆ ωˆ

Q

= Azˆ + Bzˆ (20) ωˆ + ˆωJ(p, ωˆ)

We will suppose

(21) 8 2π

0

vT J(p, ωˆ)dt = 0

and show a contradiction. Since ρ is a simple characteristic multiplier then p is the only

non-trivial 2π-periodic solution associated with ρ and so there is only one solution of (13)

associated with 1/ρ. Lemma 3.2 and (21) imply that there is a unique 2π-periodic solution

z of (20). With z and p both 2π-periodic let w = z − tp. w cannot be a multiple of p for

there would be a t0 such that t0p = z − tp or z = (t0 + t) p. Since z, p are 2π-periodic we

have z = (t0 + t + 2π) p. But then we would have t0 = t0 + 2π or 2π = 0, a contradiction.

QED

Lemma 3.2 will imply, in the case of a non-critical 2π-periodic approximate solution of

(6), that there is only one v0 in (22).

We can now give the definition of a non-critical approximate solution of (6).
STABILITY 5

Definition 3.4. The pair (ˆω, xˆ), where xˆ is at least twice continuously differentiable, is said

to be non-critical with respect to (6) if (1) the variational equation about the approximate

solution xˆ, given by (12), has a simple characteristic multiplier ρ0 with all of the other

characteristic multipliers not equal to one. (2) If v0, |v0|2 = 1, is the solution of (13)

corresponding to ρo, i.e. with multiplier 1/ρ0, then

(22) 8 2π

0

v0

T J( ̇

x, ˆ ωˆ)dt W= 0,

where

(23) J( ̇x, ˆ ωˆ) =  ̇xˆ + Bˆ  ̇xˆωˆ

4. A Perturbation Problem

In this paper we will look for an exact 2π-periodic solution, x, and an exact frequency,

ω, for (6) as a perturbation of the 2π-periodic approximate solution , ˆx, and approximate

frequency, ˆω, of (6). In particular, let

ω = ˆω + β

x = ˆx +

ωˆ

ω (24) z

Then, substituting (24) into (6) and using (10), we can write the equation for z and β as

(25) ˆωz ̇ = dX (ˆx, xˆωˆ ; z,zωˆ ) + R(z, β) − βJ

p

 ̇

x, ˆ ωˆ

Q

− k

where

(26) R(z, β) = }X

w

xˆ +

ωˆ

ω

z, xˆω +

ωˆ

ωzω

W

− X (ˆx, xˆωˆ )

]

− dX (ˆx, xˆωˆ ; z, zωˆ ) + βBˆ  ̇xˆωˆ .

and J

p

 ̇

x, ˆ ωˆ

Q

is given by (23).

In the next lemma we establish bounds and Lipschitz conditions for R(z, β).

Lemma 4.1. There exist functions R0(z, β) > 0, Ri(z, β, z,  ̃ β ̃) > 0, i = 1, 2, such that

R0 → 0 as (z, β) → 0 and Ri → 0 as (z, β, z,  ̃ β ̃) → 0 and

(27) |R(z, β)| ≤ R0(z, β)

e

e

eR(z, β) − R( ̃z, β ̃)

e

e

e ≤ R1(z, β, z,  ̃ β ̃)|z − z ̃| + R2(z, β, z,  ̃ β ̃)

e

e

eβ − β ̃

e

e

e

Proof: Appendix 2.

Since we will be considering |β| small, we will begin by restricting β, which could be

negative, so that

(28) ˆω + β ≥

ωˆ

2 .

We can select |β| ≤ ωˆ/2.

As a first step to establishing the existence of a 2π-periodic solution of (25) we first study

the existence of a 2π-periodic solution of

(29) ˆωz ̇ = dX (ˆx, xˆωˆ ; z, zωˆ ) + g − βJ

p

 ̇

x, ˆ ωˆ

Q

− k

where g ∈ P. For this we have the following lemma

Lemma 4.2. If (ˆω, xˆ) are non-critical with respect to (6), then (a) there exists a unique β

such that

(30) g − βJ

p

 ̇

x, ˆ ωˆ

Q

− k ⊥ v0

where v0 is the solution of (13) corresponding to the characteristic multiplier ρ0 of (12), and

(b) there exists a unique 2π-periodic solution of (29) that satisfies

(31) |z| ≤ M

e

e

eg − βJ

p

 ̇

x, ˆ ωˆ

Q

− k

e

e

e
6 DAVID E. GILSINN

for some M > 0.

Proof: Take

(32) β = α

}8 2π

0

v0

T (g − k) dt]

where

(33) α =

}8 2π

0

v0

T J

p

 ̇

x, ˆ ωˆ

Q

dt]−1

and apply Lemma 3.2.

We can now establish bounds on β, z and z ̇. For notation, designate the unique β and

z in Lemma 4.2 by β(g) and z(g) respectively, and z ̇ by z ̇(g).

Lemma 4.3. There exist three constants, designated by λi, i = 0, 1, 2, such that

|β(g)| ≤ λ0(|g| + r)

(34) |z(g)| ≤ λ1(|g| + r)

|z ̇(g)| ≤ λ2(|g| + r)

Proof: From

|g|2 ≤ √

2π |g|

|k|2 ≤ √ (35) 2π |k|

and the Cauchy-Schwarz inequality applied to (32)

(36) |β(g)| ≤ |α|

e

ev0

T e

e

2 |g − k|2 ≤ √

2π|α|(|g| + r)

from the bound |k| ≤ r.

From (31)

|z(g)| ≤ M



|g| + |k| + |β(g)|

e

e

eJ

p

 ̇

x, ˆ ωˆ

Qe

e

e

=

,

≤ M



1 + √

2π|α|

e

e

eJ

p

 ̇

x, ˆ ωˆ

Qe

e

e

=

(37) (|g| + r)

From (29)

|ωˆ| |z ̇(g)| ≤ |dX (ˆx, xˆωˆ ; z(g), z(g)ωˆ )| +

e

e

eg − βJ

p

 ̇

x, ˆ ωˆ

Q

− k

e

e

e

(38) ≤ [B (|z(g)| + |z(g)ωˆ |)]

≤ (2MB)



1 + √

2π|α|

e

e

eJ

p

 ̇

x, ˆ ωˆ

Qe

e

e

=

(39) (|g| + r)

Therefore, from (36), (37), (38),

λ0 = √

2π|α|,

λ1 = M



1 + √2π|α|

e

e

eJ

p

 ̇

x, ˆ ωˆ

Qe

e

e

=

(40)

λ2 = λ1

|ωˆ| M (1 + 2MB)

5. A Map and its Properties

In the main approximation theorem we will show that the solution of the perturbation

problem (25) is the fixed point of a particular contraction map. In this section we will define

the map and establish some properties.

We begin by defining a subset of P, designated by Nδ, as

(41) Nδ = {g ∈ P : |g| ≤ δ},
STABILITY 7

where δ > 0. Following Stokes [28] we will define a map S : Nδ → P in terms of two

mappings

L : Nδ → R × P1,

(42) T : R × P1 → P.

To define L, let g ∈ Nδ, then Lemma 4.2 assures us of the existence of a unique β(g)

satisfying (30) and a unique solution z(g) satisfying (29). Thus, define L : Nδ → R × P1 by

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

Now define T : R × P1 → P by

(44) T(β, z) = R(z, β).

Finally, define S : Nδ → P by

(45) S(g) = T(L(g)) = R(z(g), β(g)).

Lemma 5.1. For g ∈ Nδ, g ̃ ∈ Nδ there exist two functions E1(δ), E2(δ) and two positive

constants F1, F2 so that

(46) |S(g)| ≤ E1(δ),

|S(g) − S( ̃g)| ≤ E2(δ)|g − g ̃| ,

where

E1(δ) ≤ F1δ2 (47) ,

E2(δ) ≤ F2δ.

Proof: From (45) and (27) we have

(48) |S(g)| ≤ R0(z(g), β(g)),

|S(g) − S ( ̃g)| ≤ R1 (z(g), β(g), z ( ̃g), β ( ̃g))|z(g) − z ( ̃g)|

+R2 (z(g), β(g), z ( ̃g), β ( ̃g))|β(g) − β ( ̃g)| .

By Cauchy-Schwarz, the fact that e

ev0

T e

e

2 = 1, and (40)

|β(g) − β ( ̃g)| ≤ |α|

8 2π

0

e

ev0

T (g − g ̃)

e

e dt,

≤ |α|

}8 2π

0

|g − g ̃|

2 dt]1/2

(49) ,

≤ λ0 |g − g ̃| .

From Lemma 4.2 and the definition of β(g), β ( ̃g) we have

8 2π

0

v0

T p

g − β(g)J

p

 ̇

x, ˆ ωˆ

Q

− k

Q

dt = 0,

8 2π

0

v0

T p

g ̃ − β( ̃g)J

p

 ̇

x, ˆ ωˆ

Q

− k

Q

(50) dt = 0.

Then, by subtracting,

(51) 8 2π

0

v0

T 

(g − g ̃) − (β(g) − β ( ̃g)) J

p

 ̇

x, ˆ ωˆ

Q= dt = 0.

Lemma 4.2 also shows that there exists a unique  ̄z such that there exists a  ̄z such that

(52) ˆωz ̄ ̇ = dX (ˆx, xˆωˆ ;  ̄z, z ̄ωˆ ) + 

(g − g ̃) − (β(g) − β ( ̃g)) J

p

 ̇

x, ˆ ωˆ

Q= .

But from (29), z(g) − z ( ̃g) also satisfies (52), so that  ̄z = z(g) − z ( ̃g) and from (31)

|z(g) − z ( ̃g)| ≤ M

e

e

e(g − g ̃) − (β(g) − β ( ̃g)) J

p

 ̇

x, ˆ ωˆ

Qe

e

e,

(53) ≤ λ1 |g − g ̃| .
8 DAVID E. GILSINN

Then (46) follows from (48) through (53). QED

The specific forms for E1(δ) and E2(δ) are given in Appendix A.2 as equations (201) and

(204), respectively, as well as the selection of F1 and F2 as (202) and (205) respectively.

As functions of the other parameters E1(δ) and E2(δ) depend linearly on K and B, but

non-linearly on λ0, λ1, and λ2 and thus non-linearly on M.

6. Main Approximation Theorem

In the main theorem the constants F1, F2 are those from Lemma 4.2.

Theorem 6.1. If (a) (ˆω, xˆ) is non-critical with respect to (6) in the sense of Definition 3.4,

(b) δ is selected so that

(54) δ ≤ min{1/F1, 1/2F2, ωˆ/4λ0}

and (c) r ≤ δ, then there exists and exact frequency, ω∗, and solution, x∗, of (6) such that

|x∗ − xˆ| ≤ 4λ1δ,

|ω∗ (55) − ωˆ| ≤ 2λ0δ,

where λ0, λ1 are defined in (40) and δ is defined in (41).

Proof: Let β and z be defined as in (24). By substituting (24) into (6) we have (25).

Associated with (25) we consider (29). We then define the set Nδ in (41) and consider the

map S : Nδ → P defined in (45). From (200), (201), and (202) we have |S(g)| ≤ F1δ2

for g ∈ Nδ. Furthermore, from (203), (204), and (205) we have, for g, g ̃ ∈ Nδ, that

|S(g) − S( ̃g)| ≤ F2δ |g − g ̃|. Now, if we select δ as in (54) then F1δ2 ≤ δ and F2δ ≤ 1/2, S

maps Nδ to itself and is a contraction. The last inequality that δ satisfies in (54) assures

that β(g) satisfies (28) by way of Lemma 4.3, provided r satisfies r ≤ δ. Therefore, S has

a fixed point g∗ ∈ Nδ. This implies that there exists a unique (β∗, z∗), z∗ is 2π-periodic,

satisfying (25). Then, from (24), there exists a unique (ω∗, x∗), x∗ is 2π-periodic, satisfying

(6). From (24), with r ≤ δ,

|ω∗ − ωˆ| ≤ |β∗| ≤ λ0 (|g∗| + r) ≤ 2λ0δ,

|x∗ − xˆ| ≤

e

e

e

e

ωˆ

ωˆ + β∗

e

e

e

e

|z∗| ≤ 2λ1 (|g∗ (56) | + r) ≤ 4λ1δ.

QED

We need to introduce a note here on the relationship between r and δ. In practice the

process of determining them is iterative. We start by determining an approximate solution

and the residual r. We then compute all of the parameters that involve the approximate

solution and compute δ from (54). We then compare r against δ. If r ≤ δ we are finished,

otherwise we have to return and recompute another approximate solution with possibly

smaller residual r and iterate the process. The author is not familiar with any result that

guarantees that at some point r ≤ δ, although he suspects that this will eventually happen

in most practical problems.

7. Approximating a Solution and Frequency

An approximate solution and frequency for (6) can be developed by assuming a finite

trigonometric polynomial of the form

(57) ˆxm = a2 cost + 3m

n=2

[a2n cos nt + a2n−1 sin nt]

where the sin t term has been dropped so that we can estimate a1 = ˆω, the frequency.

Note that we have centered the approximate solution about the origin, since we assumed

X(0, 0) = 0. If we set  ̄a = (a1, a2, ··· , a2m), and

(58) Em(t,  ̄a) = a1  ̇

xˆm(t) − X (ˆxm(t), xˆm(t − a1))
STABILITY 9

then for a sufficiently fine mesh, specified by {ti : i = 1, 2, ··· , 2N}, in [0, 2π], where

(59) ti = 2i − 1

2N π,

the determining equations for  ̄a can be written as (see Urabe and Reiter [33])

F1( ̄a) = 1

N

3

2N

i=1

Em (ti,  ̄a) sin ti = 0

F2( ̄a) = 1

N

3

2N

i=1

Em (ti,  ̄a) costi = 0

F2n−1( ̄a) = 1

N

3

2N

i=1

(60) Em (ti,  ̄a) sin nti = 0

F2n( ̄a) = 1

N

3

2N

i=1

Em (ti,  ̄a) cos nti = 0

for n = 2, ··· , m.

These equations give 2m equations in 2m unknowns. Standard numerical solvers, using,

for example, Newton’s method, for nonlinear equations can be used to solve for  ̄a. The

number of harmonics, m, and the quadrature index, N, can be selected independently.

8. Floquet Theory for DDEs

The analysis of the stability of an approximate periodic solution for (1) usually involves

the following considerations. If x(t), x ∈ Rn is an approximate periodic solution of (1) of

period 2π, and ω an approximate frequency, then the linear variational equation about x(t)

can be written

(61) z ̇(t) = A(t)z(t) + B(t)z(t − ω),

where A(t) and B(t) were defined previously and are periodic, with period 2π.

We now define the period map U : Cω0 → Cω0 with respect to (61) by

(62) (Uφ) (s) = z(s + 2π)

where z(s) is a solution of (61) satisfying z(s) = φ(s) for s ∈ [−ω, 0]. In this paper we

assume ω < 2π. U is then a compact operator on Cω0, whose spectrum is at most countable

with 0 as the only possible limit point (Halanay [14]).

A Floquet theory for (61) has been developed by Stokes [26] . In particular, if σ (U)

represents the spectrum of U, then for each λ ∈ σ (U), Uφ = λφ. That is, the spectrum

consists of eigenvalues. Furthermore, the space Cω0 can be decomposed as the direct sum of

two invariant subspaces

(63) Cω0 = E(λ) ⊕ K(λ)

E(λ) is finite dimensional and composed of the eigenvectors with respect to λ. Furthermore,

σ (U|K) = σ (U) − {λ}. If {ψi}, i = 1, ··· , d is a basis for E(λ) and we let Ψ be the matrix

with columns ψj for j = 1, ··· , d, then there is a matrix G(λ) such that

(64) UΨ = ΨG(λ)

Thus we can think of Cω0 as being a countable direct sum of the invariant subspaces E(λi)

plus a possible remainder subspace, R. That is

(65) Cω0 = E(λ1) ⊕ E(λ2) ⊕ ··· ⊕ R

where R is a ”remainder” set in which any solution of (61) with initial condition in R decays

faster than any exponential.

For each of the E(λi) there is a basis set Ψi, and a matrix G(λi). If we define an at most

countable basis set {Ψi}, i = 1, 2, ··· , then we can think about U operating on ∞

i=1 E(λi)

as being represented by an infinite matrix G∞. This matrix is referred to as the monodromy
10 DAVID E. GILSINN

matrix. Its eigenvalues are called the Floquet or characteristic multipliers. The periodic

solution x(t) of (1) is stable if all of the eigenvalues of U are within the unit circle and

unstable if there is at least one with a positive real part. We note that if x(t) is an exact

periodic solution of (1) then one of the characteristic multipliers is exactly one.

9. Estimating Characteristic Multipliers

In this section we assume that the variational equation with respect to the approximate

solution, ˆx(t), can be written in the form

(66) z ̇(t) = Aˆ(t)z(t) + Bˆ(t)z(t − ωˆ)

where Aˆ(t) = Aˆ(t + 2π), Bˆ(t) = Bˆ(t + 2π) and we have reintroduced t to make the operator definitions more transparent. Let Z(t, s) be the solution of (66) such that Z(s, s) =

In, Z(t, s) = 0 for t<s where In is the n2 identity matrix on Cn. The solution Z(t, s)

is sometimes referred to as the “Fundamental Solution”. Using the variation of constants

formula for (66), Halanay [14] shows that the solution of (66) for the initial function φ ∈ Cωˆ

is given by

(67) z(t) = Z(t, 0)φ(0) + 8 0

−ωˆ

Z(t, α + ˆω)Bˆ(α + ˆω)φ(α) dα

Define the operator

(68) (Uφ)(s) = z(s + 2π)

where φ ∈ Cωˆ , s ∈ [−ωˆ, 0]. If there is a non-trivial solution z(t) of (66) such that z(t+ 2π) =

ρz(t) then ρ is a characteristic multiplier of (66). If we combine (67) with (68) and note

that z(α) = φ for α ∈ [−ωˆ, 0], then characteristic multipliers are the eigenvalues of

(69) (Uφ)(s) = Z(s + 2π, 0)φ(0) + 8 0

−ωˆ

Z(s + 2π, α + ˆω)Bˆ(α + ˆω)φ(α) dα

where φ ∈ C0. Halanay [14] shows that we can restrict s ∈ [−ωˆ, 0]. This operator is

sometimes referred to as the Monodromy Operator.

9.1. Approximating the Fundamental Solution by Spectral Collocation. In this

section we will use spectral methods to compute the fundamental solution of the linear

homogeneous delay differential equation (66). These methods are well known for collocating

solutions to partial differential equations and boundary value problems. See, for example,

Gottlieb [11] and Gottlieb and Turkel [12]. They are not as well known in delay differential

equations. In this section we use a spectral method suggested by Bueler [3] and Trefethen

[31]. The method has been reported earlier in Gilsinn and Potra [10].

The computation of the fundamental matrix used in the monodromy operator (69) requires the computation of a solution z(t) of (66) on some interval [a, b]. This will be done in

a stepwise manner. We first find a positive integer q such that a + qω ≥ b. Then we solve,

at the first step, t ∈ [a, a + ω],

(70) z ̇1(t) = A(t)z1(t) + B(t)z1(t − ω),

where z1(t−ω) = φ(s) for some function φ ∈ Cω0(a) and s = t− ω. Thus the initial problem

becomes an ordinary differential equation. Then, on [a + ω, a + 2ω] we solve

(71) z ̇2(t) = A(t)z2(t) + B(t)z2(t − ω),

where z2(a+ω) = z1(a+ω), z2(t−ω) = z1(s) for s ∈ [a, a+ω], s = t−ω. Again we solve (5)

as an ordinary differential equation. The process is continued so that on [a+ (i−1)ω, a+iω],

for i = 1, 2, ··· , q,

(72) z ̇i(t) = A(t)zi(t) + B(t)zi(t − ω),

with zi(a + (i− 1)ω) = zi−1(a+ (i − 1)ω). We then define z(t) on [a, b] as the concatenation

of zi(t) for t ∈ [a + (i − 1)ω, a + iω] and i = 1, 2, ··· , q.
STABILITY 11

Since we wish to use a Chebyshev collocation method, we will shift each interval [a+ (i−

1)ω, a + iω] to the interval [−1, 1]. For t ∈ [a + (i − 1)ω, a + iω], for i = 1, 2, ··· , q, we have

z ∈ [−1, 1] provided

(73) z = 2

ω t − (2a + (2i − 1)ω)ω .

For z ∈ [−1, 1] we have t ∈ [a + (i − 1)ω, a + iω] provided

(74) t = ω

2 z + (2a + (2i − 1)ω)2 .

We note that the point t ∈ [a + (i − 1)ω, a + iω] and t − ω ∈ [a + (i − 2)ω, a + (i − 1)ω] are

translated to the same z ∈ [−1, 1]. This is clear from

(75) 2

ω (t − ω) − (2a + (2i − 3)ωω = 2ω t − (2a + (2i − 1)ωω

Therefore we can shift the iterated delay problems

(76) z ̇i(t) = A(t)zi(t) + B(t)zi(t − ω),

for t ∈ [a + (i − 1)ω, a + iω] and i = 1, 2, ··· , q, into iterated ordinary differential equations

(77) uI

i(z) = ω

2

A4i(z)ui(z) + ω

2

B4i(z)ui−1(z)

where, for t ∈ [a + (i − 1)ω, a + iω] and associated z ∈ [−1, 1],

ui(−1) = ui−1(1)

ui(z) = zi(t)

A4 (78) i(z) = A(t)

B4i(z) = B(t)

ui−1(z) = zi(t − ω)

The initial function is

(79) u0(z) = z1(t − ω) = φ(t − ω)

for t − ω ∈ [a − ω, a].

We can now approximate the fundamental solution for (66) on [a, b] by first solving the

iterated differential equations (76) subject to

ui(−1) = ui−1(1)

(80) u0(z)=0, z ∈ [−1, 1]

u1(−1) = In

where In is the n × n identity matrix. We follow the spectral method given in Bueler

[3] in that the fundamental solution is solved for in n passes of the iteration process with

u1(−1) = ej , where ej = (0, ··· , 1, ··· , 0)T with 1 in the jth element, j = 1, 2, ··· , n.

To begin the solution process we take, for some positive integer N, the Chebyshev points

(81) ηk = cos wkπ

N

W

on [−1, 1], for k = 0, 1, ··· , N. The benefits of using these points has been discussed by

Salzer [24]. The Lagrange interpolation polynomials at these points are given by

(82) lj(z) = 

N

k=0

kW=j

w z − ηk

ηj − ηk

W

.

We have lj (ηk) = δjk. Then on [−1, 1] we set

(83) ui(z) = 3

N

j=0

ui (ηj)lj (z)
12 DAVID E. GILSINN

We also need to form

(84) uI

i(z) = 3

N

j=0

ui (ηj)l

I

j (z)

At the Chebyshev points we will designate

(85) Dkj = l

I

j(ηk)

The values for these derivatives are given in Gottlieb and Turkel [12] or Trefethen [31] but

we state the values for D here for completeness. The derivations are given in Section 15,

Appendix 1.

D00 = 2N2 + 1

6

DNN = −D00

Djj = −ηj

2(1 − η2

j ) (86) , j = 1, 2, ··· , N − 1

Dij = ci(−1)i+j

cj(ηi − ηj)

for i W= j, i, j = 0, ··· , N where

(87) ci =

l

2, i = 0 or N;

1, otherwise.

For notation, let

ui(z)=(ui1, ··· , uin)

T ,

A4i(z) = 

A4(i)

pq (z)

=

p,q=1,··· ,n

(88)

B4i(z) = 

B4(i)

pq (z)

=

p,q=1,··· ,n

We then write the collocation polynomial of uir, r = 1, ··· , n, as

(89) uir(z) = 3

N

k=0

w(i)

rk lk(z)

at the Chebyshev points (81) to get

uir (ηj ) = w(i)

rj

uI

ir (ηj ) = 3

N

k=0

w(i) (90) rk Djk

ui−1,r = w(i−1)

rj

The initial conditions for the iterated differential equations are

(91) uir (ηN ) = ui−1,r (η0),

or

(92) w(i)

rN = w(i−1)

r0

for r = 1, ··· , n.

The discretized differential equations are then given by

(93)

X3

N

k=0

w(i)

rk Djk~

r=1,n

= ω

2



A4(i)

rp (z)

=

r,p=1,n

p

w(i)

rj Q

r=1,n

+

ω

2



B4(i)

rp (z)

=

r,p=1,n

p

w(i−1)

rj Q

r=1,n
STABILITY 13

for j = 0, 1, ··· , N − 1. These provide nN equations but n(N − 1) unknowns. The other n

equations come from the initial conditions. We define the following vectors

wi =

p

w(i)

10 ··· w(i)

1N w(i)

20 ··· w(i)

2N ··· w(i)

n0 ··· w(i)

nN QT

(94)

wi−1 =

p

w(i−1)

10 ··· w(i−1)

1N w(i−1)

20 ··· w(i−1)

2N ··· w(i−1)

n0 ··· w(i−1)

nN QT

Then we can write the iterated differential equation as

(95) Dw4 i = ω

2

A4iwi +

ω

2

B4iwi−1

where D4 = D ⊗ In, the Kronecker product, and each D is given by

(96) D =

⎡

⎢

⎢

⎢

⎣

D00 ··· D0N

.

.

. .

.

. .

.

.

DN−1,0 ··· DN−1,N

0 ··· 1

⎤

⎥

⎥

⎥

⎦

The unit in the lower right introduces the initial condition, w(i)

rN , r = 1, ··· , n, equation.

Thus D4 is formed by n blocks of D down the diagonal.

The matrix A4i is given by

(97)

A4i =

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

A4(i)

11 (η0) 0 ··· 0 ··· 0 ··· 0 A4(i)1n(η0)00 ··· 0

0 ... 0 0 ··· 0 ··· 0 0 ... 0 ··· 0

0 0 A4(i)

11 (ηN−1) 0 ··· 0 ··· 0 0 A4(i)1n(ηN−1) 0 ··· 0

0 0 ··· 0 ··· 0 ··· 0 ··· 0 0 ··· 0

A4(i)

n1(η0) 0 ··· 0 ··· 0 ··· 0 A4(i)

nn(η0)00 ··· 0

0 ... 0 0 ··· 0 ··· 0 0 ... 0 ··· 0

0 0 A4(i)

n1(ηN−1) 0 ··· 0 ··· 0 0 A4(i)

nn(ηN−1) 0 ··· 0

0 0 ··· 0 ··· 0 ··· 0 ··· 0 0 ··· 0

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

B4i is structured in a similar manner except every (N + 1)th row includes an element 2/ω

to take care of the initial condition. Thus

(98)

B4i =

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

B4(i)

11 (η0) 0 ··· 0 ··· 0 ··· 0 B4(i)1n (η0)00 ··· 0

0 ... 0 0 ··· 0 ··· 0 0 ... 0 ··· 0

0 0 B4(i)

11 (ηN−1) 0 ··· 0 ··· 0 0 B4(i)1n (ηN−1) 0 ··· 0

2

ω0 0 ··· 0 ··· 0 ··· 0 ··· 0 0 ··· 0

B4(i)

n1 (η0) 0 ··· 0 ··· 0 ··· 0 B4(i) nn(η0)00 ··· 0

0 ... 0 0 ··· 0 ··· 0 0 ... 0 ··· 0

0 0 B4(i)

n1 (ηN−1) 0 ··· 0 ··· 0 0 B4(i) nn(ηN−1) 0 ··· 0

0 0 ··· 0 ··· 0 ··· 0 2

ω0 0 0 ··· 0

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

The linear equation (95) can be solved for wi by setting

(99) Mi =

w

D4 − ω

2

A4i

W−1

ω

2

B4i

and

(100) wi = Miwi−1

for i = 2, 3, ··· , q .
14 DAVID E. GILSINN

To solve for w1 for the fundamental solution we need to solve

(101) uI

1(z) = ω

2

A41(z)u1(z)

for z ∈ (−1, 1] and

(102) u1(−1) = In

That is, we solve n problems at each iteration, one for each of the initial conditions ei, where

ei is the standard basis vector with a unit in the ith element and zero elsewhere. For the

moment we set the initial vector as

(103) w0 = (0 ··· u010 ··· u020 ··· u0n)

T

where u0r, r = 1, ··· , n, is placed in each of the (N + 1)th elements and zero elsewhere.

Then from the previous construction of D4 and A41 we have

(104) w1 =

w

D4 − ω

2

A41

W−1

w0

Given that we have computed

(105) uir(z) = 3

N

k=0

w(i)

rk lk(z)

on [−1, 1] for r = 1, ··· , n we can compute the result for t ∈ [a + (i − 1)h, a + ih] by setting

(106) zir(t) = uir(z)

for r = 1, ··· , n, where

(107) z = 2

ω t − (2a + (2i − 1)ω)ω

or

(108) zir(t) = 3

N

k=0

w(i)

rk lk(

2

ω t − (2a + (2i − 1)ω)ω )

The initial condition is

(109) uir(ηN ) = ui−1,r(η0).

But on [a + (i − 1)ω, a + iω], zN = −1 corresponding to t = a + (i − 1)ω and on [a + (i −

2)ω, a + (i − 1)ω], z0 = 1 corresponding to t = a + (i − 1)ω, so that

(110) zir(a + (i − 1)ω) = zi−1,r(a + (i − 1)ω)

9.2. Estimating Monodromy Operator Eigenvalues. To approximate the monodromy

operator (69) we will require a quadrature rule that satisfies

(111)

P

3

+1

k=1

vkf (sk) →

8 0

−ωˆ

f(s) ds

as P → ∞ for each continuous function f ∈ Cωˆ . The rule is satisfied if

(112)

P

3

+1

k=1

|vk| ≤ Q,

for some Q > 0 and P = 1, 2, ··· . This is satisfied by, for example, Trapezoidal or Simpson

rules.

Let −ωˆ = s1 < s2 < ··· < sP +1 = 0, and define

(113) (Uφ) (s) = Z(s + 2π, 0)φ(0) +

P

3

+1

k=1

vkZ (s + 2π, sk + ˆω) B (sk + ˆω) φ (sk)

for φ ∈ Cωˆ .
STABILITY 15

Then, for each si ∈ [−ωˆ, 0],

(114) (Uφ) (si) = Z(si + 2π, 0)φ(0) +

P

3

+1

j=1

wjZ(si + 2π, sj + ˆω)B(sj + ˆω)φ(sj)

Since sP+1 = 0, (114) can be rewritten as

(Uφ) (si) = 3

P

j=1

wjZ(si + 2π, sj + ˆω)B(sj + ˆω)φ(sj)

(115) + (Z(si + 2π, 0) + wP+1Z(si + 2π, ωˆ)B(ˆω)) φ(sP +1),

where Z(s, α) is the fundamental matrix of (66). Equation (115) can be put in matrix form

(116)

⎛

⎜⎜⎜⎜⎜⎜⎝

(Uφ)(s1)

.

.

.

(Uφ)(si)

.

.

.

(Uφ)(sP+1)

⎞

⎟⎟⎟⎟⎟⎟⎠

=

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎣

U1,1 ··· U1,j ··· U1,P+1

.

.

. ··· .

.

. ··· .

.

.

Ui,1 ··· Ui,j ··· Ui,P +1

.

.

. ··· .

.

. ··· .

.

.

UP+1,1 ··· UP+1,j ··· UP+1,P +1

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎦

,

where the block elements for i = 1, ··· , P + 1, j = 1, ··· , P are Ui,j = wjZ(si + 2π, sj +

ωˆ)B(sj + ˆω). The block elements in the last column of the matrix are given by Ui,P+1 =

Z(si+ 2π, 0)+wP+1Z(si+ 2π, ωˆ)B(ˆω) for i = 1, ··· , P + 1. The relevant eigenvalue problem

becomes

(117)

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎣

U1,1 ··· U1,j ··· U1,P +1

.

.

. ··· .

.

. ··· .

.

.

Ui,1 ··· Ui,j ··· Ui,P+1

.

.

. ··· .

.

. ··· .

.

.

UP +1,1 ··· UP +1,j ··· UP+1,P+1

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎦

= λ

⎛

⎜⎜⎜⎜⎜⎜⎝

(φ)(s1)

.

.

.

(φ)(si)

.

.

.

(φ)(sP+1)

⎞

⎟⎟⎟⎟⎟⎟⎠

10. Determining Solutions of the Adjoint Equation Associated with Multipliers

of the Variational Equation

In order to estimate α in (33), let t ∈ [0, 2π] and ψ be the initial function defined on

[2π, 2π + ω]. The formal adjoint equation from (13) is given by

(118) v ̇

0(t) = −v0(t)A(t) − v0(t + ω)B(t + ω),

where v0(t) is a row vector. Ordinarily solving the adjoint equation would require a ”backward” integration. However, Halanay [14] showed that the solution of the formal adjoint

(118) on [0, 2π] is given in row vector form by

(119) v0(t) = ψ(2π)Z(2π, t) + 8 2π+ω0

2π

ψ(α)B(α)Z(α − ω, t) dα.

The significance of this representation is that only a “forward” integration is required to

solve for the fundamental solution Z. Let φ4(s) be a continuous row vector function defined

on [−ω, 0]. Halanay [14] then defined the operator

(120) pU4φ4

Q

(s) = φ4(−ω)Z(2π, s + ω) + 8 0

−ω0

φ4(α)B(α + ω)Z(2π + α, s + ω) dα.

Note the relationship to the monodromy operator (??). He also gave an associated operator

V4, defined on [2π, 2π + ω], as

(121) pV4ψ

Q

(s) = y(s−2π,ψ) = ψ(2π)Z(2π, s−2π)+8 2π+ω0

2π

ψ(α)B(α)Z(α−ω, s−2π) dα.

Again, note the relationship to (119). He further showed that an eigenvalue ρ0 of V4 is

associated with a 1/ρ0 multiplier of the formal adjoint equation (118), the eigenvalues of
16 DAVID E. GILSINN

U, U, 4 V4 are all the same, and the eigenvectors of U, 4 V4 are related by φ4(s) = ψ(s+2π+ω), s ∈

[ω, 0]. Although V4 is the operator associated with (119), the fact that the eigenvalues and

eigenvectors of V4 and U4 are the same allows algorithms developed for the characteristic

multipliers in Section 9 to be easily modified to compute the eigenvalues and eigenvectors

for (120). In particular we again partition the interval [−ω, 0] into P equal intervals of

length ∆ = ω/P by

(122) −ω = s1 < s2 < ··· < sP+1 = 0

We use a Simpson integration method to write, with the same weights as previously,

p

U4φ4

Q

(s) = φ4(s1)

p

Z(2π, s + ω) + w1B (s1 + ω)Z (2π + s1, s + ω)

Q

+

P

3

+1

j=2

φ4(sj )

p

wjB (sj + ω)Z (2π + sj , s + ω)

Q

(123)

To solve the adjoint equation in row form on [0, 2π], we need only compute the eigenvector

of U4 associated with the characteristic multiplier of U. The eigenvector is then substituted

into equation (119). From the previous section there are likely to be two complex conjugate

eigenvalues, ρ0 and ρ0, associated with complex conjugate eigenvectors v0 and v0 of (119),

so by linearity of (118) the real part forms a discretized solution of (118) and thus a single

real independent solution. We will continue to call the real part of this eigenvector φ4 so that

for t ∈ [0, 2π]

v0(t) = φ4(s1)

p

Z(2π, t) + w1B (s1 + ω)Z (s1 + 2π, t)

Q

+

P

3

+1

j=2

φ4(sj )

p

wjB (sj + ω)Z (sj + 2π, t)

Q

(124)

The j-th block column of U4φ4 is given by

(125)

p

U ̃φ ̃

Q

(sj) = 

φ ̃(s1), ··· , φ ̃(si), ··· , φ ̃(sP+1)

=

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

Z(2π, sj + ω) + w1B(s1 + ω)X(s1 + 2π, sj + ω)

.

.

.

wiB(si + ω)Z(si + 2π, sj + ω)

.

.

.

wP+1B(sN+1 + ω)Z(sN+1 + 2π, sj + ω)

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

.

The eigenvector φ ̃ of the matrix on the right associated with the multiplier of the variational

equation is computed and substituted into equation (124) to give the value of v0(t) on [0, 2π].

To compute α we need to estimate (13). Again we use a Simpson rule. We partition

[0, 2π] by equidistant intervals 0 = t1 < t2 < ··· < tO+1 = 2π and set h = 2π/O, where O

is an even integer. Again the weights will be set as u1 = uO+1 = h/3, otherwise uk = 4h/3

if k is even, uk = 2h/3 if k is odd. We then set

(126) α =

XO

3

+1

k=1

ukv0 (tk) J

p

 ̇

x,  ω

Q

(tk)

~−1

where v0 is a row vector and J is a column vector.
STABILITY 17

Using the same partition of [0, 2π] we can compute v0(t) at each tk as

(127)

v0(tk) = 

φ ̃(s1), ··· , φ ̃(si), ··· , φ ̃(sP +1)

=

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

Z(2π, tk) + w1B(s1 + 2π + ω)Z(s1 + 2π, tk)

.

.

.

wiB(si + 2π + ω)Z(si + 2π, tk)

.

.

.

wP+1B(sP+1 + 2π + ω)Z(sP+1 + 2π, tk)

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

.

We finally normalize v0 so that |v0|

2 = 1.

11. Estimating the M Parameter

From Halanay [14] the variation of constants formula for

(128) z ̇(t) = Aˆ(t)z(t) + Bˆ(t)z(t − ωˆ) + f(t),

where t ∈ [0, 2π], is given by

(129) z(t) = Z(t, 0)φ(0) + 8 0

−ωˆ

Z(t, α + ˆω)Bˆ(α + ˆω)z(α) dα +

8 t

0

Z(t, α)f(α) dα.

The 2π periodic initial function condition with s ∈ [−ωˆ, 0] is

φ(s) = Z(s + 2π, 0)φ(0)

+

8 0

−ωˆ

Z(s + 2π, α + ˆω)Bˆ(α + ˆω)φ(α) dα +

8 s+2π

0

(130) Z(s + 2π,α)f(α) dα.

The first step in computing M involves relating φ to f. Let |φ| = sup−ωˆ≤s≤0 |φ(s)| and

similarly for |f | on [0, 2π]. To eliminate φ(0) from (130), set s = 0 in (130) and solve for

φ(0) as

φ(0) = 8 0

−ωˆ

(I − Z(2π, 0))−1Z(2π, α + ˆω)Bˆ(α + ˆω)φ(α) dα

+

8 2π

0

(I − Z(2π, 0))−1 (131) Z(2π, α)f(α) dα.

Substitute (131) into (130) and combine terms as

φ(s) = 8 0

−ωˆ

J

Z(s + 2π, 0)(I − Z(2π, 0))−1Z(2π, α + ˆω)

(132) + Z(s + 2π, α + ˆω)] Bˆ(α + ˆω)φ(α) dα

+

8 2π

0

J

Z(s + 2π, 0)(I − Z(2π, 0))−1Z(2π, α) + Z(s + 2π,α)

o

f(α) dα.

where s ∈ [−ωˆ, 0].

Let −ωˆ = s1 < s2 < ··· < sP+1 = 0, ds = ωˆ

P , and 0 = t1 < t2 < ··· < tO+1 = 2π, dt = 2π

O . We can discretize (132) by setting

(133) φ (si) =

P

3

+1

j=1

H1(i, j)φ (sj ) +

O

3

+1

k=1

H2(i, j)f (tk),

where

H1(i, j) = vj



Z (si + 2π, 0) (I − Z(2π, 0))−1Z (2π, sj + ˆω)

+Z (si + 2π, sj + ˆω)

=

Bˆ (134) (sj + ˆω)

H2(i, j) = uk



Z (si + 2π, 0) (I − Z(2π, 0))−1Z (2π, tk)

+Z (si + 2π, tk)

=
18 DAVID E. GILSINN

In vector matrix form (133) can be written

(135)

⎛

⎜⎝

φ(s1)

.

.

.

φ(sP+1)

⎞

⎟⎠ = H1

⎛

⎜⎝

φ(s1)

.

.

.

φ(sP+1)

⎞

⎟⎠ + H2

⎛

⎜⎝

f(t1)

.

.

.

f (tO+1)

⎞

⎟⎠

Using a generalized inverse we can solve for the φ vector with minimum norm by

(136)

⎛

⎜⎝

φ(s1)

.

.

.

φ(sP+1)

⎞

⎟⎠ = (I − H1)

+ H2

⎛

⎜⎝

f(t1)

.

.

.

f(tO+1)

⎞

⎟⎠

In the second step the value of φ(0), given by equation (131), is substituted into equation

(129) and terms combined to give

z(t) = 8 0

−ωˆ

J

Z(t, 0)(I − Z(2π, 0))−1Z(2π, α + ˆω)

(137) + Z(t, α + ˆω)] Bˆ(α + ˆω)φ(α) dα

+

8 2π

0

J

Z(t, 0)(I − Z(2π, 0))−1Z(2π, α) + Z(t, α)

o

f (α) dα.

This can be discretized by setting

(138) z (tk) =

P

3

+1

i=1

H3(k, i)φ (si) +

O

3

+1

j=1

H4(k,j)f (tk),

where

H3(k, i) = vi



Z (tk, 0) (I − Z(2π, 0))−1Z (2π, si + ˆω)

+Z (tk, si + ˆω)

=

Bˆ (139) (si + ˆω)

H4(k, j) = uj



Z (tk, 0) (I − Z(2π, 0))−1Z (2π, tj )

+Z (tk, tj )

=

In vector matrix form (138) can be written

(140)

⎛

⎜⎝

z(t1)

.

.

.

z(tO+1)

⎞

⎟⎠ = H3

⎛

⎜⎝

φ(s1)

.

.

.

φ(sP +1)

⎞

⎟⎠ + H4

⎛

⎜⎝

f (t1)

.

.

.

f(tO+1)

⎞

⎟⎠

By substituting (136) into (140) we have

(141)

⎛

⎜⎝

z(t1)

.

.

.

z(tO+1)

⎞

⎟⎠ =



H3 (I − H1)

+ H2 + H4

=

⎛

⎜⎝

f (t1)

.

.

.

f(tO+1)

⎞

⎟⎠

Therefore

(142) |z| ≤ M|f|,

where M =

E

E

E

H3 (I − H1)

+ H2 + H4

E

E

E

∞.
STABILITY 19

0 1 2 3 4 5 6 7

−4

−3

−2

−1

0

1

2

3

4

x 10−15 Residual Error as a Function of Time over [0, 2π]

Time

Residual Error

Figure 1. Residual Error of Approximate Solution for the Van der Pol Equation.

12. Application to a Van der Pol Equation with Delay

In this section we will apply the main theorem to approximate the limit cycle of the Van

der Pol equation with unit delay, given by

(143)  ̈x + λ

D

x(t − 1)2 − 1

i

x ̇(t − 1) + x = 0.

Since the period of the limit cycle is unknown we introduce an unknown frequency by

substituting t/ω for t to obtain

(144) ω2x ̈ + ωλ D

x(t − ω)

2 − 1

i

x ̇(t − ω) + x = 0,

for t ∈ [0, 2π]. To compare with an approximation result obtained for ordinary differential

equations in Stokes [27], we take λ = 0.1.

The first step was to estimate an approximate 2π-periodic solution, frequency and residual

to (144). By using Galerkin’s method described in Section 7 the following approximate

solution was obtained

xˆ(t)=2.0185 cos(t)

+ 2.5771 × 10−3 sin(2t)+2.5655 × 10−2 cos(2t)

+ 1.0667 × 10−4 sin(3t) − 5.2531 × 10−4 (145) cos(3t)

− 7.1780 × 10−6 sin(4t) − 2.2043 × 10−6 cos(4t),

ωˆ = 1.0012.

where we have displayed only the first few harmonics. This solution was estimated based on

11 harmonics, 40,000 sampled points over [0, 2π], and 100 Chebyshev extreme points (81).

The residual was estimated by substituting (ˆω, xˆ) from equation (145) into equation (144)

and finding the maximum of the absolute values of the residuals obtained in the interval

[0, 2π]. The result was r = 3.1086 × 10−15. This residual is significantly better than the

one given in Stokes [27]. The distribution of the residuals for the current case is shown in

Figure 1. The phase plot of the approximate solution is shown in Figure 2. For t ∈ [0, 2π]

we can then immediately estimate |xˆ| ≤ 2.0436, |  ̇xˆ| ≤ 2.0279, | ̈xˆ| ≤ 2.1165.

In the second step, the values of the constants B and K were obtained in a straightforward

manner from the variational equation about the approximate frequency and solution given
20 DAVID E. GILSINN

−2 −1 0 1 2

−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Phase Plot of Approximate Solution

X

dX/dt

Figure 2. Phase Plot of Approximate Solution for the Van der Pol Equation.

by

(146) z ̇(t) = Aˆ(t)z(t) + Bˆ(t)Z(t − ωˆ),

where

z =

w z1

z2

W

, Aˆ(t) = w 0 1

−1/ωˆ2 0

W

,

Bˆ(t) = w 0 0

−2(λ/ωˆ)ˆx1(t − ωˆ)ˆx2(t − ωˆ) (λ/ωˆ)

D

1 − xˆ1(t − ωˆ)2

i

W

.

We use the fact that the natural norm of a matrix, H, associated with a vector norm

|x| = max1≤i≤n |xi| is |H| = max1≤i≤n

n

j=1 |hij |. With this definition it is not hard to

show that

|dX(ˆx; φ)| ≤

e

e

e

e

0 1

−1/ωˆ2 − 2(λ/ωˆ)ˆx1(t − ωˆ)ˆx2(t − ωˆ) (λ/ωˆ)

D

1 − xˆ1(t − ωˆ)2i

e

e

e

e

|φ|,

(147) ≤ 2.3776|φ|.

Therefore, for λ = 0.1, B = 2.3776. Working conservatively within the domain D =

{x ∈ C[0, 2π] : |x − xˆ| ≤ 1} it is not hard to show that

|dX(ˆxωˆ + ψ1; φωˆ ) − dX(ˆxωˆ + ψ2; φωˆ )|

(148) ≤ (6λ/ωˆ) (1 + |xˆ|)|ψ1 − ψ2| |φ|.

Then from (145) and (148) we can estimate K = 1.8157 and, from (23), we can estimate

e

e

eJ( ̇

x, ˆ ωˆ)

e

e

e ≤ 2.7546.

Next, we can estimate the characteristic multipliers of the variational equation relative

to the function ˆx(t). For the quadrature steps in Sections 9 and 10 P and O were taken

as 200 and 1200 respectively. These gave mesh widths of about 1/200 on both [−ωˆ, 0] and

[0, 2π]. Using the method of Section 9 we computed two simple conjugate eigenvalues with

magnitude 1.0430. All of the other eigenvalues have magnitudes near zero. These are, of

course, the eigenvalues of the monodromy operator U. The fundamental matrix Z in (69)

is computed using the collocation method of Section 9.1 (See Figure 3). The monodromy

operator is formulated as in Section 9. The eigenvalues of the monodromy operator U are

plotted in Figure 4. Note that the significant complex conjugate eigenvalues are near the
STABILITY 21

0 2 4 6 8

−1

−0.5

0

0.5

1

Z(1,1)

tt

z11

0 2 4 6 8

−1

−0.5

0

0.5

1

Z(2,1)

tt

z21

0 2 4 6 8

−2

−1

0

1

2

Z(1,2)

tt

z12

0 2 4 6 8

−2

−1

0

1

2

Z(2,2)

tt

z22

Figure 3. Fundamental Matrix for the Variational Equational relative to

the Approximate Solution for the Van der Pol Equation.

−1 −0.5 0 0.5 1

−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Eigenvalues for U

Real Part of Eigenvalue

Imaginary Part of Eigenvalue

Figure 4. Eigenvalues for the Monodromy Operator

unit circle but are not exactly on it. This is due to the fact that (145) is only an approximate

solution. The eigenvalues are complex conjugates because the left hand matrix in (117) is

real and non-symmetric since the fundamental solution Z is non-symmetric (See Figure

3). We can confirm that the eigenvalues of the operator U ̃ are the same as those of U.

Graphically this is shown in Figure 5.

In the next step we estimate the parameter α using the methods of Section 10. The

solution of the adjoint to the variational equation was computed using equation (127) and the

parameter α in (33) was estimated by simple quadrature, with ∆ = 2π/O for a sufficiently
22 DAVID E. GILSINN

−1 −0.5 0 0.5 1

−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Eigenvalues of U−tilde

Real Part of Eigenvalue

Imaginary Part of Eigenvalue

Figure 5. Eigenvalues for U ̃

large mesh, 0 = t1 < t2 < ··· < tO+1 = 2π, as

(149) α =

^

∆

e

e

e

e

e

O

3

+1

i=1

y(ti)J( ̇x, ˆ ωˆ)(ti)

e

e

e

e

e

−1

.

The absolute value of α is estimated as 3.3547.

If we now apply the methods of Section 11, using Aˆ(t) and Bˆ(t) defined in equation (146),

we can estimate M = 2.7618 × 102. These results allow us to estimate λ0, λ1 and λ2 in

Lemma 4.3 as λ0 = 8.4091, λ1 = 6.6736 × 103, and λ2 = 3.1720 × 104. Note the magnitude

of the parameters.

With the estimates above we can compute F1 = 2.5941 × 109, F2 = 1.0798 × 1010 from

(202) and (205) respectively. Then we compute δ = 4.6305 × 10−11 from (54). Then

F1δ2 = 5.5623 × 10−12 is less than δ and F2δ = 0.5. Furthermore r < δ. Therefore, the

conditions of the main theorem are satisfied and we can conclude from Theorem 6.1 that

there exists an exact solution x∗ and an exact frequency ω∗ of equation (144) such that

|x∗ − xˆ| ≤ 1.2361 × 10−6 and |ω∗ − ωˆ| ≤ 7.7877 × 10−10.

13. Conclusions

Although there seem to be a large number of parameters to be computed and inequalities

to be tested in order to produce the final error estimates the process is feasible. All of the

steps can be completed within a single code. The current code in Appendix 3, Section 18, has

also been built around the example in Section 12 and would have to be generalized for other

applications, but the code provides a template on which to proceed. From the computational

point of view the longest compute times involve the construction of the block matrices

(116) and (125). Computing the approximate solution and the fundamental solution of the

variational equation is relatively fast compared to these matrix constructions. It behooves

anyone wishing to apply the methods of this paper to spend some effort vectorizing the

matrix construction algorithms in Sections 9.1 and 10 as much as possible.

The parameter M in the Fredholm Lemma 3.2 is a significant parameter. From the

example above, it is clear that it would be desirable to obtain as small a value for M as

possible, since its magnitude affects the λi, i = 1, 2 parameters and λ1 appears in the

final error estimates. In particular, in the example above, the effect of M causes a very

fine residual r for the approximate solution (145) to produce a pessimistic error estimate
STABILITY 23

between the approximate solution and the exact solution in the end. From (40) the critical

parameter λ1 is linearly dependent on M.

14. Disclaimer

Certain commercial software products are identified in this paper in order to adequately

specify the computational procedures. Such identification does not imply recommendation

or endorsement by the National Institute of Standards and Technology nor does it imply

that the software products identified are necessarily the best available for the purpose.

15. Appendix 1

In this appendix we present the derivation of the differentiation matrix (86). The derivation is based on a discussion of pseudospectral Chebyshev methods given in Gottlieb et al.

[13], although a full derivation of the differentiation matrix is not given.

Lemma 15.1. For some positive integer N let the Chebyshev points be given by

(150) ηk = cos wkπ

N

W

on [−1, 1], for k = 0, 1, ··· , N. The Lagrange interpolation polynomials at these points are

given by

(151) lj(z) = 

N

k=0

kW=j

w z − ηk

ηj − ηk

W

.

We have lj (ηk) = δjk. At the Chebyshev points designate

(152) Dkj = l

I

j(ηk)

The values for these derivatives are then given as

D00 = 2N2 + 1

6

DNN = −D00

Djj = −ηj

2(1 − η2

j ) (153) , j = 1, 2, ··· , N − 1

Dij = ci(−1)i+j

cj(ηi − ηj)

for i W= j, i, j = 0, ··· , N where

(154) ci =

l

2, i = 0 or N;

1, otherwise.

Proof: The Chebyshev polynomial of degree N is given by

(155) TN (z) = cos D

N cos−1i

for z ∈ [−1, 1].

Define the polynomial

(156) gj (z) =

D

1 − z2

i

T I

N (z)(−1)j+1

cjN2 (z − zj )

for j = 0, ··· , N and c0 = cN = 2, cj = 1 for 1 ≤ j ≤ N − 1. Since T I

N (zj) will be shown

below to equal zero, TI

N (z)/ (z − zj) is a polynomial of degree N −2 so gj (z) is a polynomial

of degree N. Thus, if we can show that gj (zk) = δjk for k = 0, ··· , N, then by uniqueness

gj (z) = lj (z).
24 DAVID E. GILSINN

We first need to compute the following derivatives.

T I

N (z) = −N sin D

N cos−1 z

i

√

1 − z2

T II

N (z) = −N2 D

1 − z2i1/2

cos DN cos−1 z

i

− Nz sin DN cos−1 z

i

(1 − z2)

3/2

T III

N (z) = −N2



− sin DN cos−1 z

i

N D1 − z2i−1/2 D

1 − z2i−1

+ cos DN cos−1 z

i

(−1) D1 − z2i−2(−2z)

=

−N



sin DN cos−1 z

i D1 − z2i−3/2

+z cos DN cos−1 z

i

N D1 − z2i−1/2 D

1 − z2i−3/2

+z sin DN cos−1 z

i w−3

2

W D

1 − z2i−5/2

(−2z)

]

(157)

gI

j (z) = (−1)j+1

cjN2

^

(−2z)T I

N (z)

z − zj

+

D

1 − z2i

TII

N (z)

z − zj

+

D

1 − z2i

T I

N (z)

(z − zj)

2



= (−1)j+1

cjN2

^

Nz sin DN cos−1 z

i

(z − zj ) (1 − z2)

1/2 − N2 cos D

N cos−1 z

i

(z − zj )

+

N D1 − z2i1/2 sin D

N cos−1 z

i

(z − zj)

2



We will first establish that gj (z) = lj (z). Clearly, since cos−1 zk = kπ/N, TN (zk) = 0,

and therefore, for k W= j, k W= 0, N, j W= 0, N, gj (zk) = 0. For k = j, j W= 0, N, using TI

N (z)

and L’Hospital’s rule for

(158) limz→zj

sin DN cos−1 z

i

z − zj

= N(−1)j

D

1 − z2

j

i1/2 ,

we have gj (zj ) = 1. For j = 0, z0 = 1 so that

(159) g0(z) = (−1)j+2 D

1 − z2i1/2 sin D

N cos−1 z

i

2N(z − 1) .

For z = zk, k W= 0, g0 (xk) = 0. Again apply L’Hospital’s rule to show

(160) g0 (z0) = (−1)2

2N

limz→1

^

N cos DN cos−1 z

i

− z sin D

N cos−1 z

i

(1 − z2)

1/2



= 1.

For j = N, zN = −1 and gN (zk) = 0 for k = 0, 1, ··· , N − 1. For k = 0, use L’Hospital’s

rule to show

(161) gN (zN ) = (−1)N+2

2N

lim

z→−1

^

−z sin D

N cos−1 z

i

(1 − z2)

1/2 + N cos D

N cos−1 z

i



= 1.

Therefore, gj (z) = lj (z).

We now construct the entries in the differentiation matrix (177). These are given by

Djk = gI

k (zj ) for j, k = 0, 1, ··· , N. For k W= j, k W= 0, N, since sin(kπ) = 0 and

cos(kπ)=(−1)k,

(162) gI

j (zk) = ck(−1)j+1

cj (zk − zj )
STABILITY 25

where ck = 1. For j W= 0, N, k = 0, we have z0 = 1 and, by L’Hospital’s rule,

(163) gI

j (z0) = (−1)j+1

cjN2

^

N

1 − zj

limz→1

X

sin DN cos−1 z

i

(1 − z2)

1/2

~

− N2

1 − zj



= c0(−1)j

cj (1 − zj )

,

where c0 = 2. For j W= 0, N, k = N, we have zN = −1 and, by L’Hospital’s rule,

(164) gI

j (xN ) = (−1)j+1

cjN2

^

N

1 + zj

lim z→−1

X

sin DN cos−1 z

i

(1 − z2)

1/2

~

+

N2(−1)N

1 + zj



= cN (−1)j+N

cj (zN − zj )

where cN = 2. For j = 0, k W= 0, N,

(165) gI

0 (zk) = −1

coN2 [(1 + zk) TII

N (zk)] = ck(−1)k

c0 (zk − 1)

where ck = 1, c0 = 2. For j = 0, k = 0 we start with

(166) gI

0(z) = 1

2N2 [(1 + z)TI

N (z)]

so that

(167) gI

0(z) = 1

2N2 [TI

N (z) + (1 + z)TIIN (z)] .

Since gI

0 (z0) = limz→1 gI0(z) we need to find T I

N (1) and T IIN (1). From the construction of

T I

N (z) and L’Hospital’s rule,

(168) TI

N (1) = −N lim

z→1

X

sin DN cos−1 z

i

(1 − z2)

1/2

~

= N2.

Also

T II

N (1) = −N limz→1

^

N D1 − z2

i1/2

cos DN cos−1 z

i

+ z sin DN cos−1 z

i

(1 − z2)

3/2



= N D

1 − N2i

3

limz→1

X

sin DN cos−1 z

i

(1 − z2)

1/2

~

= N4 − N2

3 (169) .

Therefore

(170) gI

0 (z0) = gI0(1) = 2N2 + 1

6 .

For j W= 0, N, we use

TII

N (zj) = (−1)j+1N2

1 − z2

j

,

TIII

N (zj) = 3(−1)j+1N2zj

D

1 − z2

j

i2 (171) ,

cj = 1, and L’Hospital’s rule to show

gI

j (zj) = (−1)j+1

N2 limz→zj

^

−2zT I

N (z)

(z − zj) +

D

1 − z2i

T II

N (z)

(z − zj ) −

D

1 − z2i

T I

N (z)

(z − zj )

2



= (−1)j+1

2N2

J

−4zjT II

N (zj ) + D

1 − z2

j

i

TIII

N (zj )

o

= − zj

2 (1 − zj )

2 (172)

Finally, for j = N, k = N, cN = 2,

(173) gI

N (zN ) = (−1)N+1

2N2 lim z→−1

[−T I

N (z) + (1 − z)T IIN (z)] .
26 DAVID E. GILSINN

By L’Hospital’s rule

TI

N (−1) = −N lim

z→−1

sin DN cos−1 z

i

(1 − z2)

1/2

= −N2(−1)N (174)

Also, by L’Hospital’s rule,

TII

N (−1) = lim z→−1

^

−N2 D

1 − z2i1/2

cos DN cos−1 z

i

− Nz sin DN cos−1 z

i

(1 − z2)

3/2



= N3 − N

3

lim z→−1

X

sin DN cos−1 z

i

(1 − z2)

1/2

~

= N4 − N2

3 (−1)N (175)

Therefore

(176) gI

N (zN ) = −2N2 + 1

6 = −gI

0 (z0).

16. Appendix 2: Bounds and Lipschitz Condition for R(z, β)

In this section we give a proof of Lemma 4.1. A lengthy, but direct, calculation shows

R(z, β) = 8 1

0

}

X1

w

xˆ + s

ωˆ

ω

z, xˆωˆ + s (ˆxω − xˆωˆ ) + s

ωˆ

ωzω

W

− X1 (ˆx, xˆωˆ )

] ωˆ

ω

zds

+

8 1

0

}

X2

w

xˆ + s

ωˆ

ω

z, xˆωˆ + s (ˆxω − xˆωˆ ) + s

ωˆ

ωzω

W

− X2 (ˆx, xˆωˆ )

] ωˆ

ω zωds

+

8 1

0

}

X2

w

xˆ + s

ωˆ

ω

z, xˆωˆ + s (ˆxω − xˆωˆ ) + s

ωˆ

ωzω

W

− X2 (ˆx, xˆωˆ )

]

(177) (ˆxω − xˆωˆ ) ds

+

wωˆ

ω − 1

W

X1 (ˆx, xˆωˆ ) z +

wωˆ

ω − 1

W

X2 (ˆx, xˆωˆ ) zωˆ

+



X2 (ˆx, xˆωˆ ) (ˆxω − xˆωˆ ) + βX2 (ˆx, xˆωˆ )  ̇xˆωˆ

=

+

ωˆ

ω

X2 (ˆx, xˆωˆ ) (zω − zωˆ ).

From

(178) ˆxω − xˆωˆ =

8 1

0

 ̇

xˆ (t − ωˆ − sβ) (−β)ds

we have

(179) |xˆω − xˆωˆ | ≤ |β|

e

e

e  ̇

xˆ

e

e

e.

Similarly

(180) |zω − zωˆ | ≤ |β| |z ̇| .

Also, from

(ˆxω − xˆωˆ ) + β  ̇

xˆωˆ = −β

8 1

0



 ̇

xˆωˆ (t − ωˆ − sβ) −  ̇

xˆωˆ (t − ωˆ)

=

ds

= β2

8 1

0

8 1

0

 ̈ (181) xˆ (t − ωˆ − usβ) s du ds

we have

(182)

e

e

e(ˆxω − xˆωˆ ) + β  ̇

xˆωˆ

e

e

e ≤ β2

2

e

e

e ̈

xˆ

e

e

e
STABILITY 27

Using (177) through (182), along with (3) and (4), we have

(183) |R(z, β)| ≤ R(z, β),

where

R(z, β)=2K

e

e

e

e

ωˆ

ω

e

e

e

e

2

|z|

2 + 2|β||z|

|ω| (|ωˆ| + |β|)

+β2

2

w

K

e

e

e ̇

xˆ

e

e

e

2

+ B

e

e

e ̈

xˆ

e

e

e

W

+ B

e

e

e

e

ωˆ

ω

e

e

e

e

(184) |β| |z ̇|

To establish the Lipschitz condition we start with the inequality

|dX (ˆx + a11, xˆωˆ + a12; b11, b12) − dX (ˆx + a21, xˆωˆ + a22; b21, b22)|

(185) ≤ K (|b11| + |b12|) (|a11 − a21| + |a12 − a22|)

+K (|a12| + |a22|) (|b11 − b21| + |b12 − b22|)

+B (|b11 − b21| + |b12 − b22|)

We need to define some functions that will help simplify the relations somewhat. Let

γ = sβ + (1 − s)β ̃

(186) q = sz + (1 − s) ̃z

q ̇ = sz ̇ + (1 − s) ̃ ̇z

for 0 ≤ s ≤ 1, and define

ψ1(q, γ)=ˆx +

ωˆ

ωˆ + γ

q

ψ2(q, γ)=ˆxωˆ+γ +

ωˆ

ωˆ + γ

qωˆ+γ

φ1(q, γ) = − ωˆ

(ˆω + γ)

2 (187) q

φ2(q, γ) = −

}

 ̇

xˆωˆ+γ +

ωˆ

ωˆ + γ

qωˆ+γ

]

− ωˆ

(ˆω + γ)

2 q ̇ωˆ+γ

Since we have earlier chosen β, β ̃ so that

ωˆ + β ≥

ωˆ

2

ωˆ + β ̃ ≥

ωˆ

2 (188)

it is easy to see that

(189)

e

e

e

e

ωˆ

ωˆ + γ

e

e

e

e

≤ 2
28 DAVID E. GILSINN

From (186) we have the following integrals

8 1

0

|q|

2ds ≤

1

3 (|z| + |z ̃|)

2

8 1

0

|γ|ds ≤

1

2

p

|β| +

e

e

eβ ̃

e

e

e

Q

8 1

0

|q|ds ≤

1

2 (|z| + |z ̃|)

8 1

0

|γ||q|ds ≤

1

3

p

|β| +

e

e

eβ ̃

e

e

e

Q

(190) (|z| + |z ̃|)

8 1

0

|q ̇||q|ds ≤

1

3

D

|z ̇| + e

ez ̃ ̇

e

e

i

(|z| + |z ̃|)

8 1

0

|q ̇||γ|ds ≤

1

3

D

|z ̇| + e

e ̃ ̇

z

e

e

i p

|β| +

e

e

eβ ̃

e

e

e

Q

8 1

0

|q ̇|ds ≤

1

2

D

|z ̇| + e

e ̃ ̇

z

e

e

i

Define the function

(191) F(z, β) = X

w

xˆ +

ωˆ

ωˆ + β z, xˆωˆ+β +

ωˆ

ωˆ + β zωˆ+β

W

.

Taking partial derivatives of (191),

d1F(z, β; y) = ωˆ

ωˆ + β dX w

xˆ +

ωˆ

ωˆ + β z, xˆωˆ+β +

ωˆ

ωˆ + β zωˆ+β; y,yωˆ+β

W

d2F(z, β; η) = ηdX wxˆ +

ωˆ

ωˆ + β z, xˆωˆ+β +

ωˆ

ωˆ + β zωˆ+β;

− ωˆ

(ˆω + β)

2 z,−

}

 ̇

xˆωˆ+β +

ωˆ

ωˆ + β z ̇ωˆ+β

]

− ωˆ

(ˆω + β)

2 zωˆ+β

~

(192)

d1F(0, 0; y) = dX (ˆx, xˆωˆ ; y, yωˆ )

d2F(0, 0; η) = ηdX p

x, ˆ xˆωˆ ; 0, −  ̇xˆωˆ

Q

From the definition of R(z, β) and (191) we have

(193) R(x, β) − R

p

z,  ̃ β ̃

Q

= F(z, β) − F

p

z,  ̃ β ̃

Q

− d2F

p

0, 0; β − β ̃

Q

− d1F (0, 0; z − z ̃)

From the definition of γ and q in (186) we define the derivative with respect to s as

(194) dsF(q, γ; ds) = 

d1F (q, γ; z − z ̃) + d2F

p

q, γ; β − β ̃

Q= ds.

By the Fundamental Theorem of Calculus

(195) 8 1

0

dsF(q, γ; ds) = F(z, β) − F

p

z,  ̃ β ̃

Q
STABILITY 29

We can the write, using (187) and (192)

R(z, β) − R

p

z,  ̃ β ̃

Q

=

8 1

0

[d1F (q, γ; z − z ̃) − d1F (0, 0; z − z ̃)] ds

+

8 1

0



d2F

p

q, γ; β − β ̃

Q

− d2F

p

0, 0; β − β ̃

Q= ds

=

8 1

0

(196) [dX (ψ1(q, γ),ψ2(q, γ); ψ1(z, γ) − ψ1( ̃z, γ),ψ2(z, γ) − ψ2( ̃z, γ))

−dX (ψ1(0, 0), ψ2(0, 0); ψ1(z, 0) − ψ1( ̃z, 0), ψ2(z, 0) − ψ2( ̃z, 0))] ds

+

8 1

0



dX pψ1(q, γ), ψ2(q, γ); p

β − β ̃

Q

φ1(q, γ),

p

β − β ̃

Q

φ2(q, γ)

Q

−dX pψ1(0, 0), ψ2(0, 0); pβ − β ̃

Q

φ1(0, 0),

p

β − β ̃

Q

φ2(0, 0)Q= ds

From (187) we note that ψ1(0, 0) = ˆx and ψ2(0, 0) = ˆxωˆ .

Then, using (185) through (196) it is possible to show with some effort that

R1

p

z, β, z,  ̃ β ̃

Q

= 8K (|z| + |z ̃|) + w2K

e

e

e  ̇

xˆ

e

e

e +

B

|ωˆ|

W p

|β| +

e

e

eβ ̃

e

e

e

Q

R2

p

z, β, z,  ̃ β ̃

Q

= K

3

w 16

|ω ̃|

+ 4W(|z| + |z ̃|)

2

+

w

2K

e

e

e  ̇

xˆ

e

e

e + B

w

1 +

2

|ωˆ|

WW (|z| + |z ̃|)

+

2B

|ωˆ|

D

|z ̇| + e

ez ̃ ̇

e

e

i (197)

+

16K

3 |ωˆ|

(|z| + |z ̃|)

D

|z ̇| + e

e ̃ ̇

z

e

e

i

+

B

e

e

e ̈

xˆ

e

e

e

2

p

|β| +

e

e

eβ ̃

e

e

e

Q

17. Appendix A.3: Bounds and Lipschitz Conditions for S(g)

Let g ∈ N and let r = δ. Then from Lemma 5.1 and the selection of β so that ˆω +β(g) ≥ ωˆ

2 , we have

(198)

e

e

e

e

ωˆ

ωˆ + β(g)

e

e

e

e

≤ 2

and

|S(g)| = |R0(z(g), β(g))|

≤ 2K

e

e

e

e

ωˆ

ωˆ + β(g)

e

e

e

e

2

|z(g)|

2

+2|β(g)||z(g)|

|ωˆ + β(g)| (|ωˆ| + |β(g)|)

+β(g)2

2

w

K

e

e

e  ̇

xˆ

e

e

e

2

+ B

e

e

e ̈

xˆ

e

e

e

W

(199)

+B

e

e

e

e

ωˆ

ωˆ + β(g)

e

e

e

e

|β(g)||z ̇(g)|
30 DAVID E. GILSINN

If we combine (40), (198), and (199) we have

|S(g)| ≤

l

32Kλ1

2 +

16λ0λ1

|ωˆ| (|ωˆ| + 2λ0δ)

+2λ0

2

w

K

e

e

e  ̇

xˆ

e

e

e

2

+ B

e

e

e ̈

xˆ

e

e

e

W

+ 8Bλ0λ2

M

δ2 (200) .

Set

E1(δ) = l

32Kλ1

2 +

16λ0λ1

|ωˆ| (|ωˆ| + 2λ0δ)

+2λ0

2

w

K

e

e

e  ̇

xˆ

e

e

e

2

+ B

e

e

e ̈

xˆ

e

e

e

W

+ 8Bλ0λ2

M

δ2 (201) .

and let F1 be a positive constant such that

F1 ≥ 32Kλ1

2 +

16λ0λ1

|ωˆ| (202) (|ωˆ| + 2λ0δ)

+2λ0

2

w

K

e

e

e  ̇

xˆ

e

e

e

2

+ B

e

e

e ̈

xˆ

e

e

e

W

+ 8Bλ0λ2.

Now let g, g ̃ ∈ N and again set r = δ. Then, from (40), (46), and (197) and ,choosing

|g| ≤ δ, we have, with some algebra,

|S(g) − S ( ̃g)| ≤

^

λ1

l

8K (|z(g)| + |z ( ̃g)|)

+

w

2K

e

e

e  ̇

xˆ

e

e

e +

B

|ωˆ|

W

(|β(g)| + |β ( ̃g)|)

M

+λ0

l

K

3

w 16

|ωˆ|

+ 4W(|z(g)| + |z ( ̃g)|)

2

+

w

2K

e

e

e  ̇

xˆ

e

e

e + B

w

1 +

2

|ωˆ|

WW (|z(g)| + |z ( ̃g)|)

+

2B

|ωˆ|

(|z ̇(g)| + |z ̇ ( ̃g)|)

+

16K

3 |ωˆ|

(|z(g)| + |z ( ̃g)|) (|z ̇(g)| + |z ̇ ( ̃g)|)

+

B

e

e

ex ̈ˆ

e

e

e

2 (|β(g)| + |β ( ̃g)|)

M |g − g ̃|

≤

^

λ1

l

32Kλ1 + 4λ0

w

2K

e

e

e  ̇

xˆ

e

e

e +

B

|ωˆ|

W M

(203)

+λ0

l

16λ0

2

3

w 16

|ωˆ|

+ 4Wδ + 4λ1

w

2K

e

e

e  ̇

xˆ

e

e

e + B

w

1 +

2

|ωˆ|

WW

+

8Bλ2

|ωˆ| +

256Kλ1λ2

3 |ωˆ| δ

+2B

e

e

e ̈

xˆ

e

e

e

λ0

Mδ |g − g ̃|
STABILITY 31

Finally, we set

E2(δ) = ^

λ1

l

32Kλ1 + 4λ0

w

2K

e

e

e  ̇

xˆ

e

e

e +

B

|ωˆ|

W M

+λ0

l

16λ0

2

3

w 16

|ωˆ|

+ 4Wδ + 4λ1

w

2K

e

e

e  ̇

xˆ

e

e

e + B

w

1 +

2

|ωˆ|

WW (204)

+

8Bλ2

|ωˆ| +

256Kλ1λ2

3 |ωˆ| δ

+2B

e

e

e ̈

xˆ

e

e

e

λ0

Mδ

and let F2 be a positive constant such that

F2 ≥ λ1

l

32Kλ1 + 4λ0

w

2K

e

e

e  ̇

xˆ

e

e

e +

B

|ωˆ|

W M

+λ0

l

16λ0

2

3

w 16

|ωˆ|

+ 4Wδ + 4λ1

w

2K

e

e

e  ̇

xˆ

e

e

e + B

w

1 +

2

|ωˆ|

WW (205)

+

8Bλ2

|ωˆ| +

256Kλ1λ2

3 |ωˆ| δ

+2B

e

e

e ̈

xˆ

e

e

e

λ0

M

18. Appendix 3: Main Matlab Script

This section includes the main script and supporting functions, except for “cheb.m”,

which is available in Trefethen [31]. These scripts are included as is. They are not necessarily

the most efficient and are specifically oriented towards the Van der Pol equation example in

Section 12. A user will have to modify the scripts for their particular problem.

global m N CS M0 M2 V0 V1 V2 T lambda DM

global a_bar

global startt endt

global D_hat

global A_hat

global NC

global zj

global m

global piinvN

%****************************************************\\

%User input\\

%****************************************************\\

cal_B = input(’Bound on derivatives for right hand side of DDE. cal_B = ’);

cal_K = input(’Lipschitz condition on right hand side derivatives of DDE. cal_K = ’);

lambda = input(’Van der Pol Equation parameter lambda. lambda = ’);

m = input(’Approximate Solution Harmonics. m = ’);

N = input(’1/2 number of integration points for init. cond. function. N = ’);

NC = input(’Enter NC for NC + 1 Chebyshev points for collocation. NC = ’);

P = input(’Enter number of points for Trapezoidal integration on [-omega,0]. P = ’);

O = input(’Enter number integration points on [0,2pi], P/O = 1/6, e.g. 250/1500. O = ’);

%******************************************************\\

%Computing the initial approximation function\\

%******************************************************\\
32 DAVID E. GILSINN

%initialize arrays

T = zeros(2*N,1);

i = zeros(2*N,1);

V0 = zeros(2*m-1,1);

V1 = zeros(2*m,1);

V2 = zeros(2*m,1);

M0 = zeros(2*N,2*m-1);

M2 = zeros(2*N,2*m-1);

A = zeros(2*m,1);

A0 = zeros(2*m,1);

DM = zeros(2*N,2*m);

CS = zeros(2*m,2*N);

%Open I/O file

fid = fopen(’Est_Periodic_Sol_err_Output.txt’,’w’);

%Projection integration steps

i = (1:2*N)’;

T = (pi/(2*N))*(2*i-1);

%Set up fixed arrays

for n = 1:m

CS(2*n-1,:) = cos((2*n-1)*T)’;

CS(2*n,:) = sin((2*n-1)*T)’;

end

M0(:,1) = cos(T);

M2(:,1) = -cos(T);

for n = 1:m-1

M0(:,2*n) = cos((2*n+1)*T);

M0(:,2*n+1) = sin((2*n+1)*T);

M2(:,2*n) = (-(2*n+1)\^2)*cos((2*n+1)*T);

M2(:,2*n+1) = (-(2*n+1)\^2)*sin((2*n+1)*T);

end

for n = 1:m

DM(:,2*n-1) = cos((2*n-1)*T);

DM(:,2*n) = sin((2*n-1)*T);

end

%Initialize A0

A(1) = 1.0;

A(2) = 2.0;

a_bar = fsolve(’Numerical_Galerkin’,A)

fprintf(fid,’Frequency and Approximate Solution Coefficients \n$’);

nabar = length(a_bar);

for io = 1:nabar

fprintf(fid,’%15.8e \n’,a_bar(io));

end

%calculate the residual

r = Van_der_Pol(a_bar);

norm\_r = max(abs(r))

fprintf(fid,’\n\nResidual error r = %15.8e \n’,norm_r);
STABILITY 33

Lgt\_T = length(T);

xt = zeros(Lgt_T,1);

xdt = zeros(2,Lgt_T);

xdt\_temp = zeros(2,1);

for i = 1:2*N

xt(i,1) = Galerkin_series(T(i),a_bar,m);

xdt_temp = Derivative_series(T(i),a_bar,m);

xdt(1,i) = xdt_temp(1,1);

xdt(2,i) = xdt_temp(2,1);

end

norm_abs_xt = max(abs(xt(:,1)))

norm_abs_dxt = max(abs(xdt(1,:)))

norm_abs_ddxt = max(abs(xdt(2,:)))

fprintf(fid,’\n\nMax absolute value x(t) on [0,2pi] = %15.8e\n’,norm_abs_xt);

fprintf(fid,’\n\nMax absolute value dx/dt(t) on [0,2pi] = %15.8e\n’,norm_abs_dxt);

fprintf(fid,’\n\nMax absolute value d2xdt2(t) on [0,2pi] = %15.8e\n’,norm_abs_ddxt);

%*************************************************************

%Plot residual error

%*************************************************************

figure;

plot(T,r);

title(’Residual Error as a Function of Time over [0, 2\pi]’);

xlabel(’Time’);

ylabel(’Residual Error’);

disp(’press enter to continue’)

pause;

%*************************************************************

%Phase plot

%*************************************************************

nt = 10000;

t = zeros(19999,1);

ii = (1:nt)’;

t = (ii-1)*2*pi/nt;

x_hat = Phase_series(t,a_bar,m);

figure;

plot(x_hat(:,1),x_hat(:,2));

title(’Phase Plot of Approximate Solution’);

xlabel(’X’);

ylabel(’dX/dt’);

axis equal;

disp(’press enter to continue’)

pause;

%***************************************************************

% End of appoximation section

%

% Output from this section a_bar, norm_r

%***************************************************************

%***************************************************************

% Begin collocation section

%***************************************************************

%***************************************************************

% Get NC Chebyshev points, zj, and Create Chebyshev differentiation

% matrix, D. See Trefethen book.
34 DAVID E. GILSINN

%***************************************************************

[D,zj]=cheb(NC);

%***************************************************************

% Reset last row of D to account for initial condition row

%***************************************************************

D(NC+1,1:NC) = 0.0;

D(NC+1,NC+1) = 1.0;

D_hat = zeros(2*(NC+1),2*(NC+1));

D_hat(1:NC+1,1:NC+1) = D(1:NC+1,1:NC+1);

D_hat(NC+2:2*(NC+1),NC+2:2*(NC+1)) = D(1:NC+1,1:NC+1);

%*****************************************************************

% Set up fixed arrays for the Van der Pol problem

% Create the A_hat matrix in the Van der Pol case. It’s constant.

% Initialize the B array and a temporary array

%*****************************************************************

A_hat = zeros(2*(NC+1),2*(NC+1));

A = [0 1; -1/(a_bar(1,1)^2) 0];

for i = 1:NC

A_hat(i,i) = A(1,1);

A_hat(i,NC+1+i) = A(1,2);

A_hat(NC+1+i,i) = A(2,1);

A_hat(NC+1+i,NC+1+i) = A(2,2);

end

B = zeros(2,2); % coefficients of linear delay term for Van der Pol

Temp = zeros(2,2);

*****************************************************************

%Compute the fundamental matrix at several points from 0 to 2pi

%for plotting only

%get Lagrange weights from 0 to 2*pi for Fundamental Solution

%*****************************************************************

disp(’Plotting Fundamental matrix from 0 to 2pi’);

%Weights = colloc(0,2*pi);

a = 0;

b = 2*pi;

%debug ************************************************

Q = fix((b-a)/a_bar(1,1))+1;

B = zeros(2,2); tj = zeros(NC+1,1); g = zeros(2,1);

W0 = zeros(2*(NC+1),1); Weights = zeros(2*(NC+1),2,Q);

B_hat = zeros(2*(NC+1),2*(NC+1),Q);

%solve for weights on first step interval

omega = a_bar(1,1);

M = ((D_hat - (omega/2)*A_hat)^(-1));

W0(NC+1,1) = 1; W0(2*(NC+1),1) = 0;

Weights(:,1,1) = M*W0;

W0(NC+1,1) = 0; W0(2*(NC+1),1) = 1;

Weights(:,2,1) = M*W0;

%other intervals

for i = 2:Q

B_hat(NC+1,1,i) = 2/omega; B_hat(2*(NC+1),NC+2,i) = 2/omega;

tj = (omega/2)*zj + (1/2)*(2*a + (2*i-1)*omega);

for k = 1:length(tj)-1

g(1,1,1) = a_bar(2,1,1)*cos(tj(k));
STABILITY 35

g(2,1,1) = -a_bar(2,1,1)*sin(tj(k));

for n = 2:m

g(1,1,1) = g(1,1,1) + a_bar(n*2,1)*cos(n*tj(k)) + a_bar(n*2-1,1)*sin(n*tj(k));

g(2,1,1) = g(2,1,1) - n*a_bar(n*2,1)*sin(n*tj(k)) + n*a_bar(n*2-1,1)*cos(n*tj(k));

end

B(2,1,1) = (-2*lambda/a_bar(1,1,1))*g(1,1,1).*g(2,1,1);

B(2,2,1) = (lambda/a_bar(1,1))*(1-g(1,1,1).^2);

%fill up the B_hat_i matrix

B_hat(k,k,i) = B(1,1,1);

B_hat(k,(NC+1)+k,i) = B(1,2,1);

B_hat((NC+1)+k,k,i) = B(2,1,1);

B_hat((NC+1)+k,(NC+1)+k,i) = B(2,2,1);

end

Weights(:,1,i) = M*(omega/2)*B_hat(1:2*(NC+1),1:2*(NC+1),i)*Weights(:,1,i-1);

Weights(:,2,i) = M*(omega/2)*B_hat(1:2*(NC+1),1:2*(NC+1),i)*Weights(:,2,i-1);

end

%end debug ********************************************

[max_row, max_col, max_plane] = size(Weights);

%interpolate the fundamental solution

pp = 1000;

tt = zeros(pp+1,1);

zz = zeros(pp+1,1);

yy = zeros(NC+1,1,1);

del = 2*pi/pp;

for ii = 1:pp+1

%get a time value betwee 0 and 2*pi

tt(ii) = (ii-1)*del;

Z = vdp_interp(0,tt(ii),Weights);

yz11(ii,1) = Z(1,1);

yz21(ii,1) = Z(2,1);

yz12(ii,1) = Z(1,2);

yz22(ii,1) = Z(2,2);

end

figure;

subplot(2,2,1);

plot(tt(1:pp+1),yz11(1:pp+1,1));

title(’Z(1,1)’);

xlabel(’tt’);

ylabel(’z11’);

subplot(2,2,2);

plot(tt(1:pp+1),yz21(1:pp+1,1));

title(’Z(2,1)’);

xlabel(’tt’);

ylabel(’z21’);

subplot(2,2,3);

plot(tt(1:pp+1),yz12(1:pp+1,1));

title(’Z(1,2)’);

xlabel(’tt’);

ylabel(’z12’);

subplot(2,2,4);

plot(tt(1:pp+1),yz22(1:pp+1,1));

title(’Z(2,2)’);

xlabel(’tt’);

ylabel(’z22’);
36 DAVID E. GILSINN

disp(’press enter to continue’)

pause;

% clear out the arrays not needed

clear pp tt yy zz Z Weights

%********************************************************************

% Computing the monodromy matrix. This is a 2*(P+1) x 2*(P+1) matrix

% Eigenvalues by Collocation. Integration by Trapezoidal rule.

%********************************************************************

omega = a_bar(1,1);

delta = omega/P; % P is user selected number of Trapezoidal points

iii = 1:P+1;

s = -omega + (iii-1)*delta; %integration from -omega to 0 for monodromy matrix

U = zeros(2*(P+1),2*(P+1)); % the 2 is for the 2 x 2 blocks for VdP equat.

Z = zeros(2,2); %redefine Z again

endt = 2*pi;

B(1,1) = 0.0;

B(1,2) = 0.0;

delta1 = delta;

% Load by column

for jj = 1:P+1

disp(sprintf(’Creating Block column %d for U\n’,jj))

startt = s(jj) + omega;

%get the weights for the number of step intervals between

%startt = s(jj)+omega and endt = 2*pi

Weights_eig = colloc(startt,endt);

% Set up B matrix for integration step using trapezoidal rule

x_hat = Vdp_series(s(jj),a_bar,m);

B(2,1) = (-2*lambda/omega)*x_hat(1)*x_hat(2);

B(2,2) = (lambda/omega)*(1.0 - x_hat(1)^2);

% form row blocks for column j

for ii = 1:P+1

tt(ii) = s(ii) + 2*pi;

Z = vdp_interp(startt,tt(ii),Weights_eig);

if ((ii == 1)|(ii == P+1))

delta1 = delta/2;

end

Temp = delta1*Z*B;

% Fill in the i row block of the current column

U(2*ii-1,2*jj-1) = Temp(1,1);

U(2*ii-1,2*jj) = Temp(1,2);

U(2*ii,2*jj-1) = Temp(2,1);

U(2*ii,2*jj) = Temp(2,2);

end

clear weights_eig %need to clear since startt changes

end

%add the initial condition blocks to the last column

disp(’Adding the initial condition blocks to the last column’)

jj = P+1;

startt = 0;

Weights_eig = colloc(startt,2*pi);

for ii = 1:P+1

tt(ii) = s(ii) + 2*pi;
STABILITY 37

Z = vdp_interp(startt,tt(ii),Weights_eig);

U(2*ii-1,2*(P+1)-1) = U(2*ii-1,2*(P+1)-1) + Z(1,1);

U(2*ii,2*(P+1)-1) = U(2*ii,2*(P+1)-1) + Z(2,1);

U(2*ii-1,2*(P+1)) = U(2*ii-1,2*(P+1)) + Z(1,2);

U(2*ii,2*(P+1)) = U(2*ii,2*(P+1)) + Z(2,2);

end

eig(U)

disp(’Operator matrix U is filled, now computing eigenvalues’)

[V,Diag] = eig(U);

[row,col] = size(Diag);

abs_diag = zeros(row,1);

disp(’Eigenvalues of U by colloc and max absolute value’)

fprintf(fid,’\n\nEigenvalues of U\n’);

for i = 1:row

eigen(i,1) = Diag(i,i);

re = real(eigen(i,1));

im = imag(eigen(i,1));

fprintf(fid,’%15.8e + i %15.8e\n’,re,im);

end

for i = 1:row

abs_diag(i,1) = abs(Diag(i,i));

end

max_abs_diag = max(abs_diag)

fprintf(fid,’\n\nMaximum Absolute Value of Eigenvalue for U\n’)

fprintf(fid,’%15.8e\n’,max_abs_diag );

%plot the first 20 eigenvalues of U

figure

hold on;

ang = 0:pi/100:2*pi;

plot(sin(ang),cos(ang),’b-’)

plot(eigen,’rx’)

hold off;

title(’Eigenvalues for U’);

xlabel(’Real Part of Eigenvalue’);

ylabel(’Imaginary Part of Eigenvalue’);

axis equal;

%disp(’press enter to continue’)

%pause;

%*********************************************************************

%Computing solution of the adjoint associated with characteristic

%multiplier near the unit circle

%*********************************************************************

%First fill the matrix U_tilde. According to Halanay this should have

%the same eigenvalues as U. Trapezoidal integration from -omega to 0.

%*********************************************************************

U_tilde = zeros(2*(P+1),2*(P+1));

delta = omega/P;

iii = 1:P+1;

s = -omega + (iii-1)*delta;

delta1 = delta;

B(1,1) = 0.0;

B(1,2) = 0.0;

%load by columns
38 DAVID E. GILSINN

for jj = 1:P+1

disp(sprintf(’Creating Block column %d for U_tilde\n’,jj))

startt = s(jj) + omega;

Weights_adj_eig = colloc(startt, 2*pi);

%Load row blocks for colum jj

for ii = 1:P+1

%B Matrix for Van der Pol

x_hat = Vdp_series(s(ii),a_bar,m);

B(2,1) = (-2*lambda/omega)*x_hat(1)*x_hat(2);

B(2,2) = (lambda/omega)*(1.0-x_hat(1)^2);

%Interpolate at row ii Z(s(ii)+2pi,s(jj)+omega)

Z = vdp_interp(startt,s(ii)+2*pi,Weights_adj_eig);

if ((ii == 1)|(ii == P+1))

delta1 = delta/2;

end

Temp = delta1*B*Z;

%Fill in the row blocks ii of the current column jj

U_tilde(2*ii-1,2*jj-1) = Temp(1,1);

U_tilde(2*ii-1,2*jj) = Temp(1,2);

U_tilde(2*ii,2*jj-1) = Temp(2,1);

U_tilde(2*ii,2*jj) = Temp(2,2);

end

clear Weights_adj_eig

end

%add initial condition block to first row

disp(’Adding the initial condition blocks to the first row’)

endt = 2*pi;

for jj = 1:P+1

startt = s(jj) + omega;

Weights_adj_eig = colloc(startt,endt);

%Interpolate Z(2pi,s(jj)+omega)

Z = vdp_interp(startt,endt,Weights_adj_eig);

U_tilde(1,2*jj-1) = U_tilde(1,2*jj-1) + Z(1,1);

U_tilde(2,2*jj-1) = U_tilde(2,2*jj-1) + Z(2,1);

U_tilde(1,2*jj) = U_tilde(1,2*jj) + Z(1,2);

U_tilde(2,2*jj) = U_tilde(2,2*jj) + Z(2,2);

clear Weights_adj_eig

end

eig(U_tilde)

disp(’Operator matrix U_tilde is now filled,’);

disp(’computing eigenvalues and eigenvectors of transpose’);

disp(’Transposing U_tilde is necessary in order for Matlab’);

disp(’to get column eigenvectors’);

[PhiT,Dadj] = eig(U_tilde’);

[row,col] = size(Dadj);

abs_diag = zeros(row,1);

disp(’Eigenvalues by Colloc of U-tilde transpose and max abs val.’);

fprintf(fid,’\n\nEigenvalues of U-tilde\n’);

for i = 1:row

eadj(i,1) = Dadj(i,i);

re = real(eadj(i,1));

im = imag(eadj(i,1));

fprintf(fid,’%15.8e + i %15.8e\n’,re,im);

end
STABILITY 39

for i = 1:row

abs_diag(i,1) = abs(Dadj(i,i));

end

max_abs_diag = max(abs_diag);

fprintf(fid,’\n\nMaximum Absolute Value of Eigenvalue for U-tilde\n’)

fprintf(fid,’%15.8e\n’,max_abs_diag );

I_max_eig = find(abs_diag == max_abs_diag)

figure;

hold on;

ang = 0:pi/100:2*pi;

plot(sin(ang),cos(ang),’b-’);

plot(eadj,’rx’);

title(’Eigenvalues of U-tilde’);

xlabel(’Real Part of Eigenvalue’);

ylabel(’Imaginary Part of Eigenvalue’);

axis equal;

%disp(’press enter to continue’)

%pause;

%*****************************************************************

%The first eigenvector in PhiT is associated with the first eigenvector of

%largest magnitude. Transpose it to make it a row vector.

%*****************************************************************

pp(:,1) = PhiT(:,I_max_eig(1));

lgth = length(pp);

phi = zeros(1,lgth);

phi(1,1:lgth) = pp(:,1)’;

U_adj = zeros(2*(P+1),2*(O+1));

%N+1 blocks

iii = 1:P+1;

s = -omega + (iii-1)*delta;

sp2pi = s + 2*pi;

endt = 2*pi;

kk = 1:O+1;

t = (kk-1)*2*pi/O;

for k = 1:O+1

%Compute adjoint solution at the t(k)

%This will be the k-th column

disp(sprintf(’Computing adjoint at y(t(%d))\n’,k));

%Load up the column k

startt = t(k);

if (k == 1)

startt = 0.0;

end

%Get weights for Z(2pi,t(k))

Weights_adj = colloc(startt,endt);

dela1 = delta;

B(1,1) = 0.0;

B(1,2) = 0.0;

% Load row block for column

for ii = 1:P+1

%disp(sprintf(’Creating Block column %d for adjoint\n’,ii))

x_hat = Vdp_series(s(ii),a_bar,m);

B(2,1) = (-2*lambda/omega)*x_hat(1)*x_hat(2);

B(2,2) = (lambda/omega)*(1.0-x_hat(1)^2);
40 DAVID E. GILSINN

%Check relation between start and end time of fundamental solution

if (startt > sp2pi(ii))

Z = zeros(2,2);

elseif (startt == sp2pi(ii))

Z = eye(2);

else

%Interpolate Z(s(ii)+2pi,t(k))

Z = vdp_interp(startt,sp2pi(ii),Weights_adj);

end

if ((ii == 1) | (ii == P+1))

delta1 = delta/2; %Trapezoidal rule end weights

end

Temp = delta1*B*Z;

% Fill in the row blocks of the current column

U_adj(2*ii-1,2*k-1) = Temp(1,1);

U_adj(2*ii-1,2*k) = Temp(1,2);

U_adj(2*ii,2*k-1) = Temp(2,1);

U_adj(2*ii,2*k) = Temp(2,2);

end

clear Weights_adj

end

endt = 2*pi;

%add initial condition block to first block row of column k

for k = 1:O+1

startt = t(k);

%Interpolate Z(2pi,t(k))

if (k == 1)

startt = 0.0;

end

Weights_adj = colloc(startt,endt);

Z = vdp_interp(startt,endt,Weights_adj);

U_adj(1,2*k-1) = U_adj(1,2*k-1) + Z(1,1);

U_adj(2,2*k-1) = U_adj(2,2*k-1) + Z(2,1);

U_adj(1,2*k) = U_adj(1,2*k) + Z(1,2);

U_adj(2,2*k) = U_adj(2,2*k) + Z(2,2);

clear Weights_adj

end

%Apply row eigenvector to U_adj

%v_0 is 1 row by 2*(O+1) columns since phi is 1 row by 2*(P+1) columns

%and U_adj is P+1 rows by O+1 columns

%v_0 might be complex

v_0 = phi*U_adj;

%normalizing v_0

norm_v_0 = norm(v_0,2);

v_0 = v_0/norm_v_0;

%*********************************************************************

%Computing the alpha parameter

%*********************************************************************

% Now compute the alpha parameter in Stokes theorem

disp(’Computing alpha’);

% First compute the J(omega,x_dot) function at the mesh points in [0,2*pi].

% compute alpha. First comput J(i)

B(1,1) = 0.0;
STABILITY 41

B(1,2) = 0.0;

alpha = 0.0;

J1 = zeros(2,O+1);

J1_temp = zeros(2,1);

w_trap = zeros(O+1,1);

for k = 1:O+1

w_trap(k,1) = 2*pi/O;

if ((k==1)|(k==O+1))

w_trap(k,1) = pi/O;

end

%first set up B(t(k))

x_hat = Vdp_series(t(k)-omega,a_bar,m);

B(2,1) = (-2*lambda/omega)*x_hat(1)*x_hat(2);

B(2,2) = (lambda/omega)*(1.0-x_hat(1)^2);

%Multiply by the derivative at t(k) - omega

x_hat_dot_delay = Derivative_series(t(k)-omega,a_bar,m);

%add derivative at t(k)

x_hat_dot_t = Derivative_series(t(k),a_bar,m);

%compute J as 2 rows, 1 column

J = x_hat_dot_t + B*x_hat_dot_delay;

J_n(k,1) = norm(J,2);

v0(1,1) = v_0(1,2*k-1);

v0(1,2) = v_0(1,2*k);

%form alpha

alpha = alpha + w_trap(k,1)*v0*J;

end

alpha = 1/alpha

norm_J = max(J_n)

fprintf(fid,’\n\nPrinting estimate of norm of J on [0,2pi] = %15.8e\n’,norm_J);

fprintf(fid,’\n\nPrinting alpha\n’);

re = real(alpha);

im = imag(alpha);

fprintf(fid,’%15.8e + i %15.8e\n’,re,im);

%This alpha might be complex but we really only use abs(alpha)

disp(’Absolute value of alpha’);

abs_alpha = abs(alpha)

fprintf(fid,’\n\nAbsolute Value of alpha\n’);

fprintf(fid,’%15.8e\n’,abs_alpha);

%disp(’Press Enter to Continue’);

%pause;

%***********************************************************************

%Now compute the M bound

%***********************************************************************

% get integral discretization steps

ds = omega/P;

dt = 2*pi/O;

% Identify points

iP = 1:P+1;

sP = -omega + (iP -1)*ds;

iO = 1:O+1;
42 DAVID E. GILSINN

tO = (iO -1)*dt;

%Trapezoidal weights for integration

w = zeros(P+1);

u = zeros(O+1);

w(1:P+1) = ds;

u(1:O+1) = dt;

w(1) = ds/2;

w(P+1) = ds/2;

u(1) = dt/2;

u(O+1) = dt/2;

% Set up some fixed matrices

% Get collocation weights between 0 and 2pi

Weights1 = colloc(0,2*pi);

% Interpolate to form fundamental solution at 2pi, i.e. Z(2*pi,0)

Z = vdp_interp(0,2*pi,Weights1);

% Form (i - Z(2*pi,0))^(-1)

I2 = eye(2);

ImZ_inv = (I2 - Z)^(-1);

% Form (i - Z(2*pi,0))^(-1)*Z(2*pi,0)

ImZ_inv_Z = ImZ_inv * Z;

% Set up some cell arrays

disp(’Setting up cell arrays’);

B(1,1) = 0;

B(1,2) = 0;

for iPP = 1:P+1

Z = vdp_interp(0,sP(iPP)+2*pi,Weights1);

Z_cell(iPP,1) = {Z*ImZ_inv};

x_hat = Vdp_series(sP(iPP),a_bar,m);

B(2,1) = (-2*lambda/omega)*x_hat(1)*x_hat(2);

B(2,2) = (lambda/omega)*(1.0 - x_hat(1)^2);

B_cell(iPP,1) = {B};

end

clear Weights1;

% Initialize the H arrays

IP = eye(2*(P+1));

H1 = zeros(2*(P+1), 2*(P+1));

H2 = zeros(2*(P+1), 2*(O+1));

H3 = zeros(2*(O+1), 2*(P+1));

H4 = zeros(2*(O+1), 2*(O+1));

H5 = zeros(2*(O+1), 2*(O+1));

% Compute H1

disp(’Computing H1’);
STABILITY 43

for j1 = 1:P+1

disp(sprintf(’Column block %d of %d columns of H1\n’,j1,P+1));

Weights2 = colloc(sP(j1)+omega,2*pi);

Z1 = vdp_interp(sP(j1)+omega,2*pi,Weights2);

for i1 = 1:P+1

Temp1 = Z_cell{i1,1}*Z1;

Z2 = vdp_interp(sP(j1)+omega, sP(i1)+2*pi, Weights2);

Temp2 = (Temp1 +Z2)*B_cell{j1,1};

Temp2 = w(j1)*Temp2; %Trapezoidal rule weights

H1(2*i1-1,2*j1-1) = Temp2(1,1);

H1(2*i1-1,2*j1) = Temp2(1,2);

H1(2*i1,2*j1-1) = Temp2(2,1);

H1(2*i1,2*j1) = Temp2(2,2);

end

clear Weights2;

end

% Compute H2

disp(’Computing H2’);

for k2 = 1:O+1

disp(sprintf(’Column block %d of %d columns of H2\n’,k2,O+1));

Weights3 = colloc(tO(k2), 2*pi);

Z1 = vdp_interp(tO(k2),2*pi,Weights3);

for i2 = 1:P+1

if (tO(k2) > sP(i2)+2*pi)

Z2 = zeros(2,2);

elseif (tO(k2) == sP(i2)+2*pi)

Z2 = eye(2);

else

Z2 = vdp_interp(tO(k2),sP(i2)+2*pi,Weights3);

end

Temp3 = Z_cell{i2,1}*Z1 + Z2;

Temp3 = u(k2)*Temp3;

H2(2*i2-1,2*k2-1) = Temp3(1,1);

H2(2*i2-1,2*k2) = Temp3(1,2);

H2(2*i2,2*k2-1) = Temp3(2,1);

H2(2*i2,2*k2) = Temp3(2,2);

end

clear Weights3;

end

Weights1 = colloc(0,2*pi);

% Compute H3

disp(’Computing H3’);

for k3 = 1:O+1

disp(sprintf(’Row block %d of %d rows of H3\n’,k3,O+1));

Temp1 = vdp_interp(0,tO(k3),Weights1);

for j3 = 1:P+1

Weights2 = colloc(sP(j3)+omega,2*pi);

if (j3 == 1)

Temp4 = (Temp1*ImZ_inv_Z + Temp1)*B_cell{1,1};

Temp4 = w(1)*Temp4;

H3(2*k3-1,1) = Temp4(1,1);

H3(2*k3-1,2) = Temp4(1,2);

H3(2*k3,1) = Temp4(2,1);
44 DAVID E. GILSINN

H3(2*k3,2) = Temp4(2,2);

else

if (tO(k3) < sP(j3) + omega)

Temp2 = zeros(2,2);

elseif (tO(k3) == sP(j3) + omega)

Temp2 = eye(2);

else

Temp2 = vdp_interp(sP(j3)+omega,tO(k3),Weights2);

Temp3 = vdp_interp(sP(j3)+omega,2*pi,Weights2);

Temp4 = (Temp1*ImZ_inv*Temp3 + Temp2)*B_cell{j3,1};

Temp4 = w(j3)*Temp4;

H3(2*k3-1,2*j3-1) = Temp4(1,1);

H3(2*k3-1,2*j3) = Temp4(1,2);

H3(2*k3,2*j3-1) = Temp4(2,1);

H3(2*k3,2*j3) = Temp4(2,2);

end

end

clear Weights2;

end

end

%Computing H4

disp(’Computing H4’);

for l4 = 1:O+1

disp(sprintf(’Row block %d of %d rows of H4\n’,l4,O+1));

Weights3 = colloc(tO(l4),2*pi);

if (l4 == 1)

Temp1 = vdp_interp(0,2*pi,Weights1);

for k4 = 1:O+1

Temp2 = vdp_interp(0,tO(k4),Weights1);

Temp3 = Temp2*ImZ_inv*Temp1 + Temp2;

Temp3 = u(1)*Temp3;

H4(2*k4-1,1) = Temp3(1,1);

H3(2*k4-1,2) = Temp3(1,2);

H3(2*k4,1) = Temp3(2,1);

H3(2*k4,2) = Temp3(2,2);

end

else

Temp1 = vdp_interp(tO(l4),2*pi,Weights3);

for k4 = 1:O+1

if (tO(l4) > tO(k4))

Temp2 = zeros(2,2);

elseif (tO(l4) == tO(k4))

Temp2 = eye(2);

else

Temp2 = vdp_interp(tO(l4),tO(k4),Weights3);

Temp3 = vdp_interp(tO(l4),2*pi,Weights3);

Temp4 = vdp_interp(0,tO(k4),Weights1);

Temp5 = Temp4*ImZ_inv*Temp3 + Temp2;

Temp5 = u(l4)*Temp5;

H4(2*k4-1,2*l4-1) = Temp5(1,1);

H4(2*k4-1,2*l4) = Temp5(1,2);

H4(2*k4,2*l4-1) = Temp5(2,1);

H4(2*k4,2*l4) = Temp5(2,2);
STABILITY 45

end

end

end

clear Weights3;

end

%Computing final bound

H5 = H3*pinv(IP - H1)*H2 + H4;

M = norm(H5,inf);

fprintf(fid,’\n\nPrinting M bound\n’);

fprintf(fid,’%15.8e\n’,M);

disp(’Final M bound’)

M

%disp(’Press Enter to Continue’);

%pause;

%**********************************************************

%Final error estimates

%**********************************************************

abs_omega = abs(a_bar(1));

disp(’lambda_0’);

lambda_0 = sqrt(2*pi)*abs_alpha

fprintf(fid,’\n\n lambda_0 = ’);

fprintf(fid,’%15.8e\n’,lambda_0);

disp(’om_min_4L0r must be >= 0’);

om_min_4L0r = omega-4*lambda_0*norm_r

fprintf(fid,’\n\n om_min_4L0r = ’);

fprintf(fid,’%15.8e\n’,om_min_4L0r);

lambda_1 = M*(1+sqrt(2*pi)*abs_alpha*norm_J)

fprintf(fid,’\n\n lambda_1 = ’);

fprintf(fid,’%15.8e\n’,lambda_1);

lambda_2 = (lambda_1/(abs_omega*M))*(1+2*M*cal_B)

fprintf(fid,’\n\n lambda_2 = ’);

fprintf(fid,’%15.8e\n’,lambda_2);

E11 = 32*cal_K*lambda_1^2;

E12 = (16*lambda_0*lambda_1/abs_omega)*(abs_omega+2*lambda_0*norm_r);

E13 = 2*lambda_0^2 *(cal_K*norm_abs_dxt^2 + cal_B*norm_abs_ddxt);

E14 = 8*cal_B*lambda_0*lambda_2;

disp(’E1 must be <= 1’)

E1 = (E11 + E12 + E13 + E14)*norm_r

fprintf(fid,’\n\n E1 = ’);

fprintf(fid,’%15.8e\n’,E1);

E21 = 32*cal_K*lambda_1^2;

E22 = 4*lambda_0*lambda_1*(2*cal_K*norm_abs_dxt + cal_B/abs_omega);

E23 = (16*lambda_1^2*lambda_0/3)*((16/abs_omega)+4)*norm_r;

E24 = 4*lambda_1*lambda_0*(2*cal_K*norm_abs_dxt + cal_B*(1 +2/abs_omega));

E25 = 8*cal_B*lambda_2*lambda_0/abs_omega;

E26 = 256*cal_K*lambda_1*lambda_2*lambda_0*norm_r/(3*abs_omega);

E27 = 2*cal_B*norm_abs_ddxt*lambda_0^2;

disp(’E2 must be strictly < 1’)

E2 = (E21 + E22 + E23 + E24 + E25 + E26 + E27)*norm_r

fprintf(fid,’\n\n E2 = ’);

fprintf(fid,’%15.8e\n’,E2);

disp(’Final errors’)

Final_error_x_star = 4*lambda_1*norm_r

Final_error_omega_star = 2*lambda_0*norm_r
46 DAVID E. GILSINN

fprintf(fid,’\n\n Final_error_x_star = ’);

fprintf(fid,’%15.8e\n’,Final_error_x_star);

fprintf(fid,’\n\n Final_error_omega_star = ’);

fprintf(fid,’%15.8e\n’,Final_error_omega_star);

%**********************************************************

fclose(fid);

18.1. colloc.m.

function Weights = colloc(a,b)

global D_hat

global A_hat

global a_bar

global NC

global zj

global lambda

global m

Q = fix((b-a)/a_bar(1,1))+1;

B = zeros(2,2,1); tj = zeros(NC+1,1); g = zeros(2,1);

W0 = zeros(2*(NC+1),1); Weights = zeros(2*(NC+1),2,Q);

B_hat = zeros(2*(NC+1),2*(NC+1),Q);

%solve for weights on first step interval

omega = a_bar(1,1);

M = ((D_hat - (omega/2)*A_hat)^(-1));

W0(NC+1,1) = 1; W0(2*(NC+1),1) = 0;

Weights(:,1,1) = M*W0;

W0(NC+1,1) = 0; W0(2*(NC+1),1) = 1;

Weights(:,2,1) = M*W0;

%other intervals

for i = 2:Q

B_hat(NC+1,1,i) = 2/omega; B_hat(2*(NC+1),NC+2,i) = 2/omega;

tj = (omega/2)*zj + (1/2)*(2*a + (2*i-1)*omega);

for k = 1:length(tj)-1

g(1,1,1) = a_bar(2,1,1)*cos(tj(k));

g(2,1,1) = -a_bar(2,1,1)*sin(tj(k));

for n = 2:m

g(1,1,1) = g(1,1,1) + a_bar(n*2,1)*cos(n*tj(k)) + a_bar(n*2-1,1)*sin(n*tj(k));

g(2,1,1) = g(2,1,1) - n*a_bar(n*2,1)*sin(n*tj(k)) + n*a_bar(n*2-1,1)*cos(n*tj(k));

end

B(2,1,1) = (-2*lambda/a_bar(1,1,1))*g(1,1,1).*g(2,1,1);

B(2,2,1) = (lambda/a_bar(1,1))*(1-g(1,1,1).^2);

%fill up the B_hat_i matrix

B_hat(k,k,i) = B(1,1,1);

B_hat(k,(NC+1)+k,i) = B(1,2,1);

B_hat((NC+1)+k,k,i) = B(2,1,1);

B_hat((NC+1)+k,(NC+1)+k,i) = B(2,2,1);

end

Weights(:,1,i) = M*(omega/2)*B_hat(1:2*(NC+1),1:2*(NC+1),i)*Weights(:,1,i-1);

Weights(:,2,i) = M*(omega/2)*B_hat(1:2*(NC+1),1:2*(NC+1),i)*Weights(:,2,i-1);

end

18.2. Derivative series.m.
STABILITY 47

function g=Derivative_series(t,a_bar,m)

g(1,:) = -a_bar(2,1)*sin(t);

g(2,:) = -a_bar(2,1)*cos(t);

for n = 2:m

g(1,:) = g(1,:) - n*a_bar(n*2,1)*sin(n*t) + n*a_bar(n*2-1,1)*cos(n*t);

g(2,:) = g(2,:) - (n^2)*a_bar(n*2,1)*cos(n*t) - (n^2)*a_bar(n*2-1,1)*sin(n*t);

end

18.3. Galerkin series.m.

function g=Galerkin_series(t,a_bar,m)

g = a_bar(2,1)*cos(t);

for n = 2:m

g = g + (a_bar(n*2,1)*cos(n*t) + a_bar(n*2-1,1)*sin(n*t));

end

18.4. lagrint.m.

function yi = lagrint(x,y,xi)

% lagrint Interpolation with Lagrange polynomials of arbitrary degree

%

% Synopsis: yi = lagrint(x,y,xi)

%

% Input: x,y = tabulated data

% xi = point where interpolation is to be evaluated

%

% Output: yi = value of y at x = xi obtained via interpolation with

% polynomial of degree n-1, where length(y) = length(x) = n

dxi = xi - x; % vector of xi - x(1), xi - x(2), ... values

n = length(x); % degree of polynomial is n-1

L = zeros(size(y)); % preallocate L for speed

% Refer to section 10.2.2 in text for explanation of vectorized code

% used to compute Lagrange basis functions, L(j)

L(1) = prod(dxi(2:n))/prod(x(1)-x(2:n)); % j = 1

L(n) = prod(dxi(1:n-1))/prod(x(n)-x(1:n-1)); % j = n

for j=2:n-1

num = prod(dxi(1:j-1))*prod(dxi(j+1:n));

den = prod(x(j)-x(1:j-1))*prod(x(j)-x(j+1:n));

L(j) = num/den;

end

yi = sum(y.*L); % Evaluate Polynomial: sum of y(j)*L(j), j=1..n

18.5. Phase series.m.

function x_hat=Phase_series(t,a_bar,m)

x_hat(:,1)=a_bar(2,1)*cos(t);

x_hat(:,2)=-a_bar(2,1)*sin(t);

for n = 2:m

x_hat(:,1) = x_hat(:,1) + a_bar(n*2,1)*cos(n*t) + a_bar(n*2-1,1)*sin(n*t);

x_hat(:,2) = x_hat(:,2) - n*a_bar(n*2,1)*sin(n*t) + n*a_bar(n*2-1,1)*cos(n*t);

end
48 DAVID E. GILSINN

18.6. Van der Pol.m.

function vdp = Van_der_Pol(A)

global m N CS M0 M2 V0 V1 V2 T lambda DM

vdp = (A(1)^2)*xdd(A) + A(1)*lambda*(delay_x(A).^2 - 1).*delay_xd(A) + x(A);

18.7. vdp interp.m.

function Z = vdp_interp(startt,eval_pt,Weights)

%VDP_INTERP Lagrange interpolation using Chebyshev Points to generate a 2 x

%2 matrix Z of the fundamental solution for Van der Pol equation. Frequency in a_bar(1).

global a_bar

global NC

global zj

omega = a_bar(1,1);

kk = fix((eval_pt - startt)/omega )+1;

zz = (2/omega)*eval_pt - (2*startt + (2*kk-1)*omega)/omega ;

yy(1:NC+1,1,1) = Weights(1:NC+1,1,kk);

Z(1,1) = lagrint(zj,yy,zz);

yy(1:NC+1,1,1) = Weights(NC+2:2*(NC+1),1,kk);

Z(2,1) = lagrint(zj,yy,zz);

yy(1:NC+1,1,1) = Weights(1:NC+1,2,kk);

Z(1,2) = lagrint(zj,yy,zz);

yy(1:NC+1,1,1) = Weights(NC+2:2*(NC+1),2,kk);

Z(2,2) = lagrint(zj,yy,zz);

18.8. Vdp series.m.

function x_hat=Vdp_series(t,a_bar,m)

%set up for dde23

x_hat(1,:)=a_bar(2,1)*cos(t);

x_hat(2,:)=-a_bar(2,1)*sin(t);

for n = 2:m

x_hat(1,:) = x_hat(1,:) + a_bar(n*2,1)*cos(n*t) + a_bar(n*2-1,1)*sin(n*t);

x_hat(2,:) = x_hat(2,:) - n*a_bar(n*2,1)*sin(n*t) + n*a_bar(n*2-1,1)*cos(n*t);

end

References

[1] B. Balachandran, Nonlinear dynamics of milling processes’, Philosophical Transactions of the Royal

Society of London, Series A, 359 (2001) 793-819.

[2] H. T. Banks and F. Kappel, Spline Approximations for Functional Differential Equations, Journal of

Differential Equations, 34 (1979) 496—522.

[3] E. Bueler, Chebyshev Collocation for Linear, Periodic Ordinary and Delay Differential Equations: A

Posteriori Estimates, http://arxiv.org/abs/math.NA/0409464

[4] E. A. Butcher, H. Ma, E. Bueler, V. Averina, and Z. Szabo, Stability of linear time-periodic delaydifferential equations via Chebyshev polynomials, Int. J. Numer. Meth. Engng, 59 (2004), 895—922.

[5] K. Engelborghs, T. Luzyanina, K. J. In ’T Hout, and D. Roose, Collocation Methods for the Computation

of Periodic Solutions of Delay Differential Equations, SIAM J. Sci. Comp., 22 (2000) 1593—1609.

[6] D. E. Gilsinn, Discrete Fourier Series Approximation to Periodic Solutions of Autonomous Delay Differential Equations, DETC 2005-84038, Proceedings of IDETC/CIE 2005: ASME 2005 International

Design Engineering Technical Conferences & Computers and Information in Engineering Conference,

September 24-28, Long Beach, CA.

[7] D. E. Gilsinn, Approximating Limit Cycles of a Van der Pol Equation with Delay, Proceedings of Dynamic

Systems and Applications 4, Ed. G. S. Ladde, N. G. Medhin, M. Sambandham, Dynamic Publishers,

Inc., 2004, 270-276

[8] D. E. Gilsinn, Estimating Critical Hopf Bifurcation Parameters for a Second-Oder Delay Differential

Equation with Application to Machine Tool Chatter, Nonlinear Dynamics 30 (2002) 103-154.
STABILITY 49

[9] D. E. Gilsinn, Computable Error Bounds for Approximate Periodic Solutions of Autonomous Delay

Differential Equations, Nonlinear Dynamics, accepted for publication

[10] D. E. Gilsinn and F. A. Potra, Integral Operators and Delay Differential Equations, Journal of Integral

Equations and Applications, 18, 3 (2006) 297—336.

[11] D. Gottlieb, The Stability of Pseudospectral-Chebyshev Methods, Math. Comp., 36 (1981) 107—118.

[12] D. Gottlieb and E. Turkel, Topics in Spectral Methods, Lecture Notes in Mathematics: 1127 A. Dold

and B. Eckmann (Ed.), Springer-Verlag, Berlin, 1985.

[13] D. Gottlieb, M. Y. Hussaini, and S. A. Orszag, Theory and Applications of Spectral Methods, in

Spectral Methods for Partial Differential Equations, R. G. Voigt, D. Gottlieb, and M. Y. Hussaini (Ed.),

SIAM, Philadelphia, (1984) 1—54.

[14] S. Halanay, Differential Equations: Stability, Oscillations, Time Lags, Academic Press, New York,

1966.

[15] J. Hale, Functional Differential Equations, Springer-Verlag, New York, 1971.

[16] N. H. Hanna and S. A. Tobias, A theory of nonlinear regenerative chatter, ASME Journal of Engineering

for Industry, 96 (1974) 247—255.

[17] T. Insperger and G. St ́ep ́an, Stability of the Damped Mathieu Equation with Time Delay, Journal of

Dynamic Systems, Measurement, and Control, 125 (2003) 166—171.

[18] T. Luzyanina and K. Engelborghs, Computing Floquet Multipliers for Functional Differential Equations,

International Journal of Bifurcation and Chaos, 22 (2002 2977—2989.

[19] T. Luzyanina and D. Roose, Numerical stability analysis and computation of Hopf bifurcation points

for delay differential equations, Journal of Computational and Applied Mathematics 72 (1996) 379—392.

[20] G. A. Kemper, Spline Function Approximation for Solutions of Functional Differential Equations, SIAM

J. Num. Anal. 12 (1975) 73—88.

[21] E. R. Marsh, D. S. Yantek, M. A. Davies, and D. E. Gilsinn, Simulation and measurement of chatter

in diamond turning, Journal of Manufacturing Science and Engineering, 120 (1995) 230—235.

[22] A. H. Nayfeh, C. Chin,and J. Pratt, Perturbation Methods in Nonlinear Dynamics - Applications to

Machining Dynamics, Journal of Manufacturing Science and Engineering, 119 (1997) 485—493.

[23] C. A. H. Paul, Developing a delay differential equation solver, Applied Numerical Mathematics, 9 (1992)

403—414.

[24] H. E. Salzer, Lagrangian interpolation at the Chebyshev points xn,ν ≡ cos (νπ/n) , ν =)(1)n; some

unnoted advantages, The Computer Journal, 15 (1972) 156-159.

[25] L.F. Shampine and S. Thompson, Solving DDE’s in MATLAB, Applied Numerical Mathematics, 37

(2001) 441—458.

[26] A. Stokes, A Floquet Theory for Functional Differential Equations, Proc. Nat. Acad, Sci., 48 (1962)

1330—1334.

[27] A. Stokes, On the Approximation of Nonlinear Oscillations, Journal of Differential Equations, 12 (1972)

535—558.

[28] A.P. Stokes, On the Existence of Periodic Solutions of Functional Differential Equations, Journal of

Mathematical Analysis and Applications, 54 (1976) 634—652.

[29] E. Stone and A. Askari, Nonlinear models of chatter in drilling processes. Dynamical Systems 17 No. 1

(2002) 65—85.

[30] E. Stone and Sue Ann Campbell, Stability and bifurcation analysis of a nonlinear DDE model of drilling,

Journal of Nonlinear Science 14 No. 1 (2004) 27—57.

[31] L. N. Trefethen, Spectral Methods in Matlab, Society for Industrial and Applied Mathematics, 2000.

[32] J. Tlusty, Machine dynamics, in Handbook of High-speed Machine Technology, Robert I. King (ed.),

Chapman and Hall, New York, (1985) 49—153.

[33] M. Urabe and A. Reiter, Numerical Computation of Nonlinear Forced Oscillations by Galerkin’s Procedure, Journal of Mathematical Analysis and Applications, 14 (1966) 107—140.

[34] M. X. Zhao and B. Balachandran, Dynamics and stability of milling process, International Journal of

Solids and Structure, 38 (2001) 2233—2248.

[35] D. R. Will ́e and C. T. H. Baker, DELSOL - a numerical code for the solution of systems of delaydifferential equations, Applied Numerical Mathematics, 9 (1992) 223—234.