University of Southern California Dissertations and Theses

Nonlinear dynamics and nonlinear dynamical systems

(USC Thesis Other)

Nonlinear dynamics and nonlinear dynamical systems

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

NONLINEAR DYNAMICS AND NONLINEAR DYNAMICAL SYSTEMS

by

Byungrin Han

------------------------------------------------------------------------------------------------------

A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(MECHANICAL ENGINEERING)

August 2007

Copyright 2007 Byungrin Han

ii
ACKNOWLEDGEMENTS

First, I would like to state my deepest appreciation to my advisor, Professor
Firdaus Udwadia, for his priceless guidance, indefatigable assistance, and
constant encouragement. I am also grateful to the other dissertation committee
members, Professor Marijan Dravinski, Professor Wlodek Proskurowski,
Professor Henryk Flashner, and Professor Bingen Yang for their willingness to
serve on my dissertation committee. I wish to thank my parents to whom I
dedicate this dissertation for their endless support and patient expectations. I
would like to send my deepest appreciate to my wife who has always been
besides me with her grateful encouragement. I wish to thank all my families for
their earnest considerations. A special thanks to David Yum for his grateful
advice. Also, I would like to state my warmest gratitude to my grandmother
who is in heaven. Most importantly, I give thanks to God who makes all things
possible.

iii
TABLE OF CONTENTS

ACKNOWLEDGEMENTS ii

LIST OF FIGURES v

ABSTRACT vii

CHAPTER 1: EXPLICIT EQUATIONS OF MOTION FOR CONSTRAINED
MECHANICAL SYSTEMS WITH NON-SINGULAR MASS MATRIX
1.1 Introduction 1
1.2 Derivation 4
1.3 Equivalence with the previous equation of motion
when the mass matrix symmetric and positive definite. 12
1.4 Modifications through Examples 14
1.4.1 Example 14
1.4.2 Modification 16
1.5 Conclusions 20

CHAPTER 2: FLYING OBJECT NONLINEAR TRACKING CONTROL
USING FUNDAMENTAL EQUATION OF MECHANICS
2.1 Introduction 23
2.2 Model Description 25
2.3 The Fundamental Equation 26
2.4 Equations of Motion of the Following Object (SO) 29
2.5 Numerical Example 35
2.6 Conclusions 39

CHAPTER 3: A VERIFICATION OF THE RECURSIVE FORMULAE FOR
COMPUTING THE MOORE-PENROSE M-INVERSE OF A MATRIX
3.1. Introduction 42
3.2 Verification 44
3.2.1. The Formulae of the Generalized M-Inverse of a Partitioned Matrix
by Adding a Row Vector. 44
3.2.2. The Formulae of the Generalized M-Inverse of a Partitioned Matrix
by Adding a Column Vector. 52
3.3 Conclusions 63

CHAPTER 4: A DYNAMICAL MODEL OF TERRORISM
4.1. Introduction 64
4.2. The Dynamical Model 66
4.3. Model Dynamics 72
4.3.1 On the Orbits of the Dynamical System. 73
4.3.2 Fixed Points of the Dynamical System. 75
4.3.3 Stability of Fixed Points. 78
4.3.4 Baseline Case 81
4.3.5 Effects of Nonviolent Intervention 92
4.3.6 Effects of Military/Police Intervention 96

iv
4.3.7 Effects of Combined Military/Police and Nonviolent Intervention 104
4.4 Conclusions 111

APPENDICES
APPENDIX A 118
APPENDIX B 124
APPENDIX C 131

REFERENCES 134

v
LIST OF FIGURES

Figure 1.1: A two degree-of-freedom multi-body system. 15
Figure 1.2: A two degree-of-freedom system using new coordinate system 16
Figure 2.1a: The trajectory of the first object and the second object (1) 37
Figure 2.1b: The trajectory of the first object and the second object (2) 37
Figure 2.2: The control forces acting on the first object and the second object 38
Figure 2.3a: The figure shows the distance between two objects 38
Figure 2.3b: The distance between two objects in 200 sec. 39
Figure 4.1: Schematic showing the Dynamical System 71
Figure 4.2a: Phase plot for the flow field and the asymptotic behavior (1) 83
Figure 4.2b: The time evolution of the dynamical system 84
Figure 4.3: Phase plot for the flow field and the asymptotic behavior (2) 85
Figure 4.4a: Phase plot for the flow field and the asymptotic behavior (3) 86
Figure 4.4b: Phase plot for the flow field and the asymptotic behavior (4) 87
Figure 4.5: Time history of the dynamical system 87
Figure 4.6: Phase portrait for the vector field and the periodic orbits 89
Figure 4.7: The x-axis becomes a line of a fixed points 90
Figure 4.8a: Stable limit cycle behavior when c = 0.2 < 1. 91
Figure 4.8b: Time histories of the dynamics when c = 1. 91
Figure 4.8c: Unstable behavior showing an explosion in the Terrorist population
92

Figure 4.9: The presence of nonviolent intervention 93
Figure 4.10: A dramatic drop in the equilibrium population of Terrorists 94
Figure 4.11: Phase plot showing the ineffectiveness of nonviolent intervention 95
Figure 4.12: A stable spiral is seen when 1 < c , 0 < A . 98

vi
Figure 4.13: A stable node is seen when 1 < c , 0 > A . 98
Figure 4.14: The effect of the parameter c in its influence on
the steady state values and the dynamical behavior of the system. 99

Figure 4.15: Time histories of response of the dynamical system with c = 0.7 100
Figure 4.16: Time histories of response of the dynamical system with c = 1.2 100
Figure 4.17: A comparison with Figures 4.14(a) and 4.15 104
Figure 4.18: The effect of both military and nonviolent action (1) 106
Figure 4.19: The effect of both military and nonviolent action (2) 107

vii
ABSTRACT

The concept of Dynamics exists everywhere in our lives. As an integral
component of daily reality, engineers and scientists research, analyze, and
observe the behaviors of dynamics and dynamical systems and their application
in our real world activities. As a general rule, dynamical systems are
characterized by their nonlinear quality. By their very nature, such systems
have presented scientists with enormous difficulty in finding unequivocal
solutions to systems, and even when a solution is possible, the path to finding one
is extremely complicated. Oftentimes, it has proven virtually impossible to find
these solutions. Given the inherent challenges in coming up with solutions to
systems, the research that we’ve focused on has involved four separate but related
areas.
In chapter 1 we will present the new general explicit equations of motion for
constrained discrete dynamical systems. These systems are characterized by
having holonomic and/or non-holonomic constraints, which satisfy D’Alembert’s
principle at each instant of time. These new equations provide a new and simple
method of handling of constrained mechanical systems with the mass matrix that
should not be a symmetric and/or positive definite.
In chapter 2 we present a simple and robust approach for tracking the motion of
nonlinear mechanical system by dynamically controlling it so that one flying
objects follows the motion of another flying object. Recent development in
analytical dynamics has led us to establish a set of control forces to create such
tracking between nonlinear dynamical systems. By providing an example, the
simplicity and efficiency of the method are illustrated.

viii
In chapter 3 we confirm the recursive formulae for the determination of the
generalized Moore-Penrose M-inverse of a matrix. We will provide a verification
of the formulae using four MP M-properties.
In chapter 4, as part of an application of dynamical systems in a real world setting,
we present an approach to understanding the dynamics of terrorism.
1
CHAPTER 1

EXPLICIT EQUATIONS OF MOTION FOR CONSTRAINED
MECHANICAL SYSTEMS WITH NON-SINGULAR MASS MATRIX

1.1 Introduction

Since Lagrange’s epoch (in 1787) over 200 years ago, the determination of the
equation of motion for a constrained discrete mechanical system has been one of the
core ideas in analytical dynamics. The initial formulation of discrete mechanical
system was developed by Lagrange [37] who introduced the concept that later
became known as the Lagrange multiplier. The method is often difficult to apply
when the system has a large number of degrees of freedom and/or non-integrable
constraints.

Commonly called today as Gauss’s principle, this approach is one of the most
integral principles in the field of mechanics introduced by Gauss [22]. He also
emphasized that a system of particles at each instant of time chooses one set of
accelerations among possible sets and the selected set of accelerations minimized are
often called Gaussian.

In the late 19
th
century, Gibbs [24] and Appell [6] independently developed the so-
called Gibbs-Appell equations of motion using the concept of quasi-coordinates.
2
They also have struggled with that same issue whenever faced with a large number
of degrees of freedom and/or non-integrable constraints. Using Poisson brackets,
Dirac [16] offers a recursive scheme for Hamiltonian systems with singular
Lagrangians where the constraints do not explicitly depend on time

The research set out by Udwadia and Kalaba [63] established a simple and fresh
approach to constrained mechanical systems with a symmetric and positive definite
mass matrix. Their explicit equations of motion supply a distinctive perspective on
constrained motion. By integrating the notion of generalized inverses in describing
such motion and, through the application of these ideas, they postulate a simple and
general explicit equation of motion for constrained mechanical systems that is
consequently free of the need for the notion of Lagrange multipliers. Their equations
are even capable of handling conditions with holonomic and/or non-holonomic
constraints that are not necessarily independent.

These equations are applicable to general mechanical systems and include situations
where the constraints may be nonlinear functions of the velocities, explicitly
dependent on time, and, functionally dependent. However, they deal only with
systems where the constraints are ideal and satisfy D’Alembert’s principle. The
principle indicates that at every instant of time the sum total of the work done under
virtual displacements by the forces of constraint is zero.

3
The set of Lagrange equations usually yields a mass matrix which is symmetric and
positive definite [42]. However, an asymmetric and/or non-positive definite mass
matrix can arise in Newtonian mechanics and, as not a mechanics field, Neutrino
physics [39, 38].

Very recently in 2006, more general types of explicit equations of motion for
constrained mechanical systems were developed by Udwadia and Phohomsiri [65].
Their equations of motion are even applicable to the singular case of the mass matrix,
which can arise in the modeling complex multi-body mechanical systems.

In this chapter we extend the general and explicit equations of motion for constrained
discrete dynamical system obtained by Udwadia and Kalaba [63] to the mass matrix
in a given mechanical system that should not be symmetric and/or positive definite.
The extended equations are also applicable to systems with holonomic and/or
nonholonomic constraints, where the constraint forces are ideal. Compared with the
previously obtained equations of motion [65], the extended equations of motion
obtained here developed independently and in some ways a simpler form when the
given mass matrix is non-singular.

However, these extend equations of motion have equivalence with previous results
[63] for symmetric, positive definite mass matrix. The straightforward and simple
proof of the equivalence will be provided with some useful relationship of
generalized Moore-Penrose inverse of, somehow, specifically defined matrix.
4

This newly obtained extended equation, in general, can be applied except the system
has an asymmetric and non-positive definite mass matrix which is defined by
coordinate transformation. In this case, the extended equations of motion should be
modified. We provide examples to show how systems with non-positive definite and
asymmetric mass matrices can be produced with coordinate transformation. We also
provide some modifications of the newly obtained equations of motion to explain
how to handle this specific case.

1.2. Derivation

Consider first an unconstrained, discrete dynamical system which is described by the
n generalized coordinates . ] , . . . , , [
2 1
T
n
q q q q = By ‘unconstrained’ we mean that
the components,
i
q & , of the velocity of the system can be independently defined at
any given initial time
0
t t = . Its equation of motion can be obtained, using newtonian
or lagrangian mechanics, by the relation

) , , ( ) , ( t q q Q q t q M & & & = , (1.1)

where the n by n matrix M is not-singular but it should not be a symmetric nor a
positive definite matrix. The matrix M and the generalized force n-vector (n by 1
5
matrix), ) , , ( t q q Q & , are known. The generalized acceleration of the unconstrained
system, which we denote by the n-vector a, is then given by

) , , (
1
t q q a Q M q & & & = =
−
. (1.2)

We next suppose that the system is subjected to h holonomic constraints of the form

0 ) , ( = ϕ t q
i
, i = 1, 2, . . . , h, (1.3)

and h m − nonholonomic constraints of the form

0 ) , , ( = ϕ t q q
i
& , m h h i , . . . , 2 , 1 + + = . (1.4)

The initial conditions =
0
q ) (
0
t t q = and
0
q & ) (
0
t t q = = & are assumed to satisfy these
constraints so that 0 ) , (
0 0
= ϕ t q
i
, i = 1, 2, . . . , h, and 0 ) , , (
0 0 0
= t q q
i
& ϕ , , 1 + = h i
m h , . . . , 2 + . We note that the constraints may be explicit functions of time, and the
nonholonomic constraints may be nonlinear in the velocity components
i
q & . Under
the assumption of sufficient smoothness, we can differentiate equations (1.3) twice
with respect to time and equations (1.4) once with respect to time to obtain the
consistent equation set

) , , ( ) , , ( t q q b q t q q A & & & & = (1.5)
6

where the constraint matrix, A, is a known m by n matrix and b is a known m-vector.
It is important to note that for a given set of initial conditions, equation set (1.5) is
equivalent to equations (1.3) and (1.4), which can be obtained by appropriately
integrating the set (1.5).

The presence of the constraints (1.5) imposes additional forces of constraint on the
system that alter its acceleration so that the explicit equation of motion of the
constrained system becomes

) , , ( ) , , ( t q q Q t q q Q q M
c
& & & & + = . (1.6)

The additional term,
c
Q , on the right hand side arises by virtue of the imposed
constraints prescribed by equations (1.5).

We assume that for any virtual displacement vector, ), (t v the total work done can be
described as

0 = =
c T
Q v W . (1.7)

We begin by stating our result for the constrained system described above. The
explicit equation of motion for the constrained system is:
7

) ( ) (
1
b B Ma M U b B q M
T + − + +
− + = & & (1.8)

where the m by n matrix
1 −
= AM B ,
+
B denotes the generalized Moore-Penrose
inverse [23] of the matrix B, and the n by n matrix ) ( ) )( (
1
B B I M B B I U
T + − +
− − = .
In equation (1.8) we have noted that the acceleration a, at the time t, corresponding
to the unconstrained motion does not satisfy the constraint equation (1.5).

Since the properties of the generalized Moore-Penrose inverse,
+
B , shall shortly use,
we provide here some of its properties. Given an m by n matrix B, the n by m
matrix
+
B is a unique matrix that satisfies the following four relations:

(1) B B BB =
+
; (1.9a)
(2)
+ + +
= B BB B ; (1.9b)
(3)
+
BB is symmetric ; and, (1.9c)
(4) B B
+
is symmetric. (1.9d)

The derivation of our result is as follows.

The acceleration, q & & , of the constrained system must satisfy two requirements. It
must be such that :
(1) at each instant of time t it must satisfy the constraints equation (1.5), and,
8
(2) the work done under any virtual displacement, ), (t v by the force of constraint,
c
Q , must, at each instant of time t, be as prescribed by relation (1.7).

Since inverse of the matrix M can be defined, the constraints equation (1.5) can be
written as following form:

b q BM = & & (1.10)

where again m by n matrix
1 −
= AM B . We require the acceleration of the constrained
system to satisfy the consistent set of equations (1.10), we have, from the theory of
generalized inverses,

z B B I b B q M ) (
+ +
− + = & & (1.11)

so that

) , , ( ) , , ( ) ( t q q Q t q q Q z B B I b B q M
c
& & & & + = − + =
+ +
(1.12)

where z is any arbitrary n-vector, and
+
B is the generalized Moore-Penrose inverse
of the matrix B and the properties of
+
B are described in equations (1.9a)-(1.9d).

9
To explicitly find
c
Q , we next determine the second member on the right in equation
(1.11) in such a way guarantee that the second of the above-mentioned requirements
is satisfied. A virtual displacement at time t is any displacement that satisfies the
relation

0 = Av (1.13)

at that time [63]. Similarly in equation (1.10), the equation (1.13) can be described
as

0 = μ B (1.14)

where Mv = μ . But the explicit solution of this homogeneous set of equations (1.14)
is simply

w B B I ) (
+
− = μ (1.15)

where w is any arbitrary n-vector. Using the relationship (1.9d) we obtain
) ( ) ( B B I B B I
T + +
− = − and the relation, Mv = μ , gives
T T T
M v ) (
1 −
= μ . So that we
get
T T T
M B B I w v ) )( (
1 − +
− = . From relation (1.7) we require that

0 ) )( (
1
= − = =
− + c T T c T
Q M B B I w Q v W . (1.16)
10

We observe that from relation (1.6) and (1.12) the force of constraint,
c
Q , can have
the form as

Ma z B B I b B t q q Q q M t q q Q
c
− − + = − =
+ +
) ( ) , , ( ) , , ( & & & & . (1.17)

When we substitute equation (1.17) in equation (1.16), we get

0 ] ) ( [ ) )( (
1
= − − + −
+ + − +
Ma z B B I b B M B B I w
T T
(1.18)

which, because w is arbitrary and when we denote ) ( ) )( (
1
B B I M B B I U
T + − +
− − = ,
yields

0 ) )( ( ) )( (
1 1
= − − + −
− + + − +
Ma M B B I Uz b B M B B I
T T
(1.19)

and

) ( ) )( (
1
b B Ma M B B I Uz
T + − +
− − = . (1.20)

Again, from the theory of generalized inverses, we observe

η ) ( ) ( ) )( (
1
U U I b B Ma M B B I U z
T + + − + +
− + − − = (1.21)
11

where η is arbitrary n vector. Since
+ + +
= − U B B I U ) ( (see Property 1 in Appendix
A), substituting equation (1.21) in equation (1.11), we then have

] ) ( ) ( ) ( )[ (
1
η U U I b B Ma M U B B I b B q M
T + + − + + +
− + − − + = & & . (1.22)

Finally, since 0 =
+
BU (see Property 5 in Appendix A), equation (1.22) is simply

η )] ( [ ) ( ) (
1
B B U U I b B Ma M U b B q M
T + + + − + +
+ − + − + = & & . (1.23)

When we define a (m+n) by n matrix
⎥
⎦
⎤
⎢
⎣
⎡
=
U
B
H , since 0 =
+
UB and 0 =
+
BU (see
Property 4 and 5 in Appendix A), we observe ] [
+ + +
= U B H . So that the last term
of equation (1.23) becomes H H I B B U U I
+ + +
− = + − ) ( . We should note that, in
general, because of the third term on the right-hand side of equation (1.23), the
explicit equation of motion of the constrained system with non-symmetric mass
matrix is not necessarily unique. However, when the previously defined (m+n) by n
matrix, H , has full rank (rank=n), the third member in equation (1.23) vanishes for
then
T T
H H H H
1
) (
− +
= so that I H H =
+
. So the equation of motion becomes
unique. Therefore the confirmation of the condition, the matrix H has full rank, is
required before the relation (1.23) applied as a unique equation of motion.

12
1.3. Equivalence with the previous equation of motion when the mass matrix
symmetric and positive definite.

In this chapter, we will announce that when the mass matrix, M , has symmetric and
positive definite condition, the equation of motion given by equation (1.23) is simply
the equation of motion given in [63]. When the matrix M is symmetric and positive
definite matrix the matrix H has full rank and the last term on the right in equation
(1.23) disappears since I H H =
+
. The result given by

) ( ) (
2 / 1 2 / 1
Aa b AM M Ma q M − + =
+ −
& & (1.24)

and the equation (1.24) can be rewritten as

b AM M a A AM M M q M
+ − + −
+ − = ) ( ) ) ( (
2 / 1 2 / 1 2 / 1 2 / 1
& & . (1.25)

Also, when
1 1
) (
− −
= M M
T
and I H H =
+
, the equation (1.23) can be rearranged as

b B M U I a U q M
+ − + +
− + = )] ( [
1
& & . (1.26)

But when the mass matrix, M , is symmetric and positive definite, the Moore-
Penrose generalized matrix of U can be determined as (see Property 17 in Appendix
A)
13

A AM M M U
+ − +
− = ) (
2 / 1 2 / 1
. (1.27)

When we substitute relation (1.27) into relation (1.26), it is shown that the first terms
of each equation are identical. Using the relationship of
1 −
= AM B and equation
(1.27), the second member of equation (1.26) can be described as

b AM M A AM M M I b B M U I
+ − − + − + − +
− − = − ) )]( )( ) ( ( [ )] ( [
1 1 2 / 1 2 / 1 1

b AM AM AM M
+ − − + −
= ) ( ) (
1 1 2 / 1 2 / 1
. (1.28)

In equation (1.28), when the matrix, M , is symmetric and positive definite, we can
apply the relationship as
+ − − + − −
= ) ( ) (
2 / 1 2 / 1 1 1
AM AM AM AM (see Property 8 in
Appendix A) so that equation (1.28) can be rewritten as

b AM M b B M U I
+ − + − +
= − ) ( )] ( [
2 / 1 2 / 1 1
, (1.29)

which makes the second member of equation (1.26) identical with the last term of
equation (1.25). Hence the relation (1.26) is the same as relation (1.25) alternatively
equation (1.24). One thing we should note that when the mass matrix, M , is not a
positive definite matrix,
2 / 1 −
M can not be determined, and the equation of motion in
equation (1.23) becomes inapplicable to the equation (1.24) or alternatively (1.25).

14
1.4. Modifications through Examples.

Though the equation (1.23) can be formulated without hassle when the control
problem is clearly defined, we can not guarantee it gives us correct equations of
motion for a given system made by coordinate transformation. This is because the
virtual displacement will be changed when the coordinate transformation occurs,
equation (1.7), alternatively (1.13), can not remain the same. Although this is a very
important result, it is not easy to realize. In this section, through the use of an
example we shall demonstrate how a system with non-symmetric mass matrix can
arise in classical mechanics and how it can be handled. Also we shall provide the
way in which the modification in equation (1.23) is required and how this artificially
created example can be handled.

1.4.1 Example

We first provide equations of motion with symmetric mass matrix using equation
(1.24) which is previously obtained [63]. Then, using coordinate transformation of
the same system, we made an artificial non-symmetric mass matrix and provide
equations of motion by the newly obtained equation (1.23) with modification.

Consider a system of two masses,
1
m and
2
m , connected with three springs,
1
k ,
2
k ,
and
3
k , as shown in Figure 1.1. We choose the coordinates
1
x and
2
x . The
constraint equation is given by
15

t x x ω cos
2 1
= − (1.30)

where w is frequency of the constraint of the system.

Figure 1.1: A two degree-of-freedom multi-body system.

Using Newtonian mechanics, the equation of motion for the unconstrained system in
matrix form can be described as

⎥
⎦
⎤
⎢
⎣
⎡
=
⎥
⎦
⎤
⎢
⎣
⎡
⎥
⎦
⎤
⎢
⎣
⎡
+ −
− +
+
⎥
⎦
⎤
⎢
⎣
⎡
⎥
⎦
⎤
⎢
⎣
⎡
0
0
0
0
2
1
3 2 2
2 2 1
2
1
2
1
x
x
k k k
k k k
x
x
m
m
& &
& &
. (1.31)

here we denote the symmetric and positive definite mass matrix
⎥
⎦
⎤
⎢
⎣
⎡
=
2
1
0
0
m
m
M .
Differentiating the equation (1.30) twice we get

[] t
x
x
ω ω cos 1 1
2
2
1
− =
⎥
⎦
⎤
⎢
⎣
⎡
−
& &
& &
(1.32)

1
m
2
m
1
k
2
k
3
k
1
x
2
x
16
which is the form of equation (1.5) and [ ] 1 1 − = A . Using equation (1.31), (1.32)
and (1.24), we obtain the equation of motion for the constrained system shown in
Figure 1 as

⎥
⎦
⎤
⎢
⎣
⎡
− +
+ +
+
− =
⎥
⎦
⎤
⎢
⎣
⎡
t m x k x k
t m x k x k
m m x
x
ω ω
ω ω
cos
cos 1
2
1 2 3 1 1
2
2 2 3 1 1
2 1 2
1
& &
& &
. (1.33)

1.4.2 Modification

Figure 1.2: A two degree-of-freedom system shown in Figure 1 using new coordinate system,
1
y and
2
y .

Consider a system shown in Figure 1.2 which is the same as Figure 1.1 except the
coordinate transformation from
1
x and
2
x to
1 1
x y = and
2 1 2
x x y − = . The
relationship between two coordinates in a matrix form is represented by

⎥
⎦
⎤
⎢
⎣
⎡
=
⎥
⎦
⎤
⎢
⎣
⎡
2
1
2
1
x
x
T
y
y
(1.34)
1
m
2
m
1
k
2
k
3
k
1
y
2
y
17

where T is a coordinate transformation matrix which is defined by
⎥
⎦
⎤
⎢
⎣
⎡
−
=
1 1
0 1
T . So
that the equation of motion for the unconstrained system in equation (1.31) can be
changed in matrix from

⎥
⎦
⎤
⎢
⎣
⎡
=
⎥
⎦
⎤
⎢
⎣
⎡
⎥
⎦
⎤
⎢
⎣
⎡
+ −
+
⎥
⎦
⎤
⎢
⎣
⎡
⎥
⎦
⎤
⎢
⎣
⎡
− 0
0
) (
0
2
1
3 2 3
2 1
2
1
2 2
1
y
y
k k k
k k
y
y
m m
m
& &
& &
& &
& &
. (1.35)

here we denote the non-symmetric mass matrix
⎥
⎦
⎤
⎢
⎣
⎡
−
=
2 2
1
0
~
m m
m
M (we use M
~
to
avoid confusion with M .) One thing we should note that since the mass matrix is no
longer positive definite matrix,
2
1
~
−
M is not well defined, and we can not apply the
equation (1.24).

Using this transformation matrix, equation (1.32) can be represented as

[] [ ] t
y
y
y
y
T ω ω cos 1 0 1 1
2
2
1
2
1 1
− =
⎥
⎦
⎤
⎢
⎣
⎡
=
⎥
⎦
⎤
⎢
⎣
⎡
−
−
& &
& &
& &
& &
(1.36)

so that [] 1 0
~
− = A (we use A
~
to avoid confusion with A ) and t b ω ω cos
2
− = .

18
Before we apply the equation (1.23), one thing should be noted that the previously
mentioned statement of requirement, according to this example, the work done under
any virtual displacement, ) (
2 , 1
t v
x
, by the force of constraint,
) , , (
2 , 1 2 , 1
t x x Q
c
&
, must, at each instant of time t, be as prescribed by
0 ) , , (
2 , 1 2 , 1
2 , 1
= = t x x Q v W
c T
x
& should be changed to the work done under any virtual
displacement, ) (
2 , 1
1
t v T
y
−
, by the force of constraint, ) , , (
2 , 1
1
2 , 1
1
t y T y T Q
c
&
− −
, must, at
each instant of time t, be as prescribed by

0 ) , , ( )] ( [
2 , 1
1
2 , 1
1 1
2 , 1
= =
− − −
t y T y T Q t v T W
c T
y
& . (1.37)

Using equation (1.37) , A
~
, and, M
~
, equation (1.16) can be modified as

0 ) , , ( )
~
)(
~ ~
(
2 , 1
1
2 , 1
1 1 1
= − =
− − − − +
t y T y T Q M T B B I w W
c T T
& (1.38)

where
1
~ ~ ~
−
= M A B . Using equation (1.38) and following the same procedure of
derivation which is done above in section 1.2, a new modified form of explicit
equation of motion for the constrained system is then given by

η )]
~ ~ ~ ~
( [ )
~ ~
( )
~
(
~ ~ ~
1 1
B B U U I b B a M M T U b B q M
T + + + − − + +
+ − + − + = & & (1.39)

where now
19

)
~ ~
( )
~
)(
~ ~
(
~
1 1
B B I M T B B I U
T + − − +
− − = (1.40)

which is not an identical form with U . Now we can apply equation (1.39) to get the
correct explicit equations of motion for this example.

Since [] [ ]
1 1
2 1
1 1
/ 1 / 1
/ 1 / 1
0 / 1
1 0
~ ~ ~
m m
m m
m
M A B − =
⎥
⎦
⎤
⎢
⎣
⎡
−
= =
−
, we obtain

⎥
⎦
⎤
⎢
⎣
⎡
− +
=
+
2
1
2
2
2
1
2
2
2
1
/ 1
/ 1
~
m
m
m m
m m
B . (1.41)

Using relation (1.41), T , and , M
~
, we obtain

⎥
⎦
⎤
⎢
⎣
⎡
+
+
=
2
2 1 2 1
2 1
2 2
2
2
1
2
1 2 1
) / ( /
/ 1
) (
) ( ~
m m m m
m m
m m
m m m
U (1.42)

and

⎥
⎦
⎤
⎢
⎣
⎡
+
=
+
2
1 2 1 2
1 2
2 1
2
1
) / ( /
/ 1
) (
~
m m m m
m m
m m
m
U . (1.43)

20
Note that the Moore-Penrose inverse of the matrix B
~
and U
~
can be easily obtained
in most computing environments, such as MATLAB or MAPLE.

After some algebra the relationships (1.41)-(1.43) and B
~
gives )
~ ~ ~ ~
( B B U U
+ +
+ is a 2
by 2 identity matrix so that the last term of equation (1.39) is vanished. Hence the
equation (1.39) becomes unique.

Finally using relations (1.35)-(1.36) and substituting equations (1.41) and (1.43) into
equation (1.39), the explicit equations of motion for this non-symmetric mass matrix
is obtained by

⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
− + +
+ − =
⎥
⎦
⎤
⎢
⎣
⎡
t
y k y k k t m
m m
y
y
ω ω
ω ω
cos
] ) ( cos [
1
2
2 3 1 3 1
2
2
2 1
2
1
& &
& &
. (1.44)

When we substitute
1 1
x y = and
2 1 2
x x y − = , it can be an easily confirmed equation
(1.44) and the previously obtain equation (1.33) is identical.

1.5 Conclusions

The new extended explicit equations of motion for constrained discrete dynamical
systems are hereby obtained. These equations of motion are applicable even when
21
the mass matrix in the given system is non-singular but should not be symmetric
and/or positive-definite. The main contributions of this chapter are as follows.

1. The previously obtained explicit equations of motion [63] are not applicable to
mechanical systems whose unconstrained equations of motion have asymmetric
and/or non-positive mass matrix. In this chapter we present a new, slightly more
general and explicit form of the equations of motion for constrained mechanical
systems. The equations explicitly provide the acceleration and apply to the systems
that contain holonomic and/or non-holonomic constraints. These equations of motion
have different form with previously obtained equations [63] and [65]. It may be
noted that our derivation does not require the notion of Lagrange multipliers.

2. As we stated during our derivation in section 1.2, when the previously mentioned
(m+n) by n matrix H matrix has full rank (rank=n), the equation (1.23) can have a
unique form. So that the requirement of the matrix, H , should be ensured before
applying the equation (1.23). When the condition, H has full rank, satisfying the
relation (1.23) can be used distinctively.

3. The equation (1.23), in general, is applicable except in the case of a system not
having a symmetric and positive definite mass matrix which is defined by coordinate
transformation. In this special case, the extended equations of motion should be
modified. Otherwise it can not be guaranteed that the equation of motion for a given
system is properly defined. We provide modifications and examples in section 1.4 to
22
show how systems with non-positive definite and asymmetric mass matrix can arise
with coordinate transformation in Newtonian mechanics. We also provide some
modifications of the newly obtained equations of motion to explain how to handle
this specific case. Through this modification, it can be shown that the equations of
motion using the relation (1.23) transform into acceptable form.

4. In section 1.3 we provided a simple and straightforward proof of equivalence of
our extended equations of motion with previously obtained equations of motion [63]
when the mechanical system has a symmetric and positive definite mass matrix. It
should be noted that since the other pre-obtained equations of motion [65] also have
identical result with [63] when the mass is symmetric and positive definite, as a
result, given the same circumstances our outcome and the previous findings in [65]
have an identical result.

5. During the proof of equivalence in section 1.3 we fortunately obtained the
generalized Moore-Penrose inverse of the matrix U . In a mathematical point of view,
obtaining the generalized inverse of the matrix U by itself is not a trivial work.
Combining mechanics, in our case the relation (1.24) and (1.26), we have shown a
relatively painless way to solve this issue.
23
CHAPTER 2

FLYING OBJECT NONLINEAR TRACKING CONTROL USING
FUNDAMENTAL EQUATION OF MECHANICS

2.1. Introduction

In general the tracking control has great utility in real word applications such as
aircraft, vehicles, spacecraft, gyrodynamics, Global Positioning System (GPS), and
medical science. Recently, the tracking control problem has established itself as one
of the central issues in the area of mechanics. Emphasizing air defense in the 21
st

century, the tracking of a flying object necessitates an interest in the engineering
field as a global concern. Many investigators, in particular the U.S. government,
have closely examined this tracking control problem, spending significant resources
and effort in studying this issue.

The tracking of a flying object is an important problem in nonlinear science and it is
investigated via several methods. Some of these include using fuzzy cerebellar
model arithmetic computer neural networks [23], interacting multiple model
algorithm [14], knowledge-based adaptive thresholding[30], and extended interval
kalman filter [54]. Because this nonlinear tracking control has uncountable
uncertainties, it is very difficult to exactly replicate the motion of the flying object.
24
Our study proposes a new and simple strategy that is inspired by recent advances in
analytical dynamics [63] for implementation in a nonlinear mechanical system.

In this chapter we place focus on a fresh and simple approach to track a flying object,
what we’ve labeled the first object (FO), in light of new developments in analytical
dynamics [63]. In our system, the FO displays a nonlinear behavior. This system is
approached as a tracking control problem where the other flying object, say the
second object (SO), exactly follows motion of FO after some time. The tracking
problem is then reconfigured as a constrained motion problem in which we seek to
have the nonlinear forces of SO to replicate the same motion of FO. From here, the
explicit closed form analytical control given by the fundamental equation [63] is
utilized to give us the control force that provides an exact tracking motion of the FO.
This approach to pursue the flying FO will show us several advantages. First, the
control forces found are continuous functions of time. These forces are in a closed
form resulting in an efficient determination. The minimum forces are necessary for
tracking. The result of exact tracking can be accomplished with relative ease. Finally,
the calculation of obtaining the tracking control forces can be computed in almost
real time.

This chapter is organized as follows. In section 2.2 we provide a short description of
a simple model of FO. We assume that the FO will reach a constant speed and height
exponentially. In section 2.3 we provide brief information of the fundamental
equation. This equation deals with explicit equation of motion for a general nonlinear
25
mechanical system which satisfies a set of constraints. In section 2.4 we apply the
fundamental equation to control and track the SO which can exactly follow the
motion of FO. In section 2.5 an example which represents can be a realistic problem
will be illustrated. In the last section we present our conclusions.

2.2. Model Description

Consider two objects in the same inertial frame of reference in a gravity field. The
locations of the two objects are defined as ) 3 , 2 , 1 ( , ) ( = i t x
i
, and ) 3 , 2 , 1 ( , ) ( = i t y
i
, at
time t. The first object (FO) launched from a certain location, ) 3 , 2 , 1 ( ,
0
= i x
i
, with
initial velocity, ) 3 , 2 , 1 ( , ) 0 (
0
= = = i x t x
i i
& & , after that the second object (SO) launched
from another location, ) 3 , 2 , 1 ( ,
0
= i y
i
, with initial velocity, ) 3 , 2 , 1 ( , ) 0 (
0
= = = i y t y
i i
& & ,
to follow the previously launched object, FO. We assume that after some time the
first object will travel with constant velocity,
1
x
v , in
1
x direction and
2
x
v in
2
x direction. During the flight, after some time the first object will have constant
height, h . Those three conditions can be represented by the following three
constraint equations;

0 ) (
1
1 1 1
= − +
x
v x x & & & α , (2.1)

0 ) (
2
2 2 2
= − +
x
v x x & & & α , (2.2)

26
and

0 ) (
3 3 3
= − + + h x x x γ β & & & , (2.3)

where
1
α ,
2
α , β , and γ are non-zero constants and the dots refer to differentiation
with respect to time.

Our goal is for the SO to follow the FO motion and after some time the SO will be
located within a certain area around the still flying FO. We will provide the force
acting on the SO using the fundamental equation [63] and it should be noted that the
applied force given by the equation is minimized at each instant of time.

2.3. The Fundamental Equation

As we will be utilizing the fundamental equation of mechanics to control the given
system, we will briefly examine relevant points in this section. This equation
involves the explicit equation of motion for a mechanical system when the system is
constrained to satisfy a set of consistent constraints. Consider an unconstrained
discrete mechanical system whose equation of motion is described by the equations

) , , ( ) , ( t q q f q q t M & & & = , , ) 0 (
0
q q =
0
) 0 ( q q & & = (2.4)

27
where, M is an n by n symmetric, positive definite matrix, the n-vector q represents
the generalized coordinates used to describe the configuration of the system and the
right hand side is a known function of q q & , and t. The dots refer to differentiation
with respect to time. By unconstrained we mean here that the components of the
initial velocity
0
q & can be arbitrarily specified. Equation (2.4) results from the
application of Lagrange’s equations to a mechanical system, or from Newtonian
mechanics.

Let this system be subjected to a set of s constraints of the form

0 )) ( ( = t q h , (2.5)

that are satisfied by the initial conditions so that

0 ) (
0
= q h , and 0 ) , (
0 0
= q q h &
&
. (2.6)

Here h is an s-vector. Differentiating equation (2.5) twice with respect to time we
obtain the set of matrix equations

) , ( ) , ( q q b q q q A & & & & = (2.7)

28
where A is an s by n matrix. The equation of motion of the constrained system that
satisfies these constraints exactly is then explicitly given by [63],

) , , ( ) , , ( t q q F t q q f q M
c
& & & & + = , (2.8)

where,

) ( ) ( ) , , (
1 2 / 1 2 / 1
f AM b AM M t q q F
c − + −
− = & . (2.9)

Here,
+
X denotes the Moore-Penrose (MP) inverse of the matrix X (see, Ref. [63]).
We shall denote the n components of the n-vector
c
F by
c
i
f , i = 1, 2, . . . , n.

We should recognize that the constraint (2.5) is actually implemented as equation
(2.7). In what follows, we shall suppress the arguments of the various quantities
unless needed for clarity.

When relations (2.6) are not satisfied by the initial conditions one could replace the
equation set (2.7) by any other system of constraint equations [58] whose solution
asymptotically tends to , 0 = h as ∞ → t . For example, the system of equations

0 = Σ + Δ + h h h
& & &
, (2.10)

29
where Δ and Σ are diagonal matrices with positive entries, would lead to 0 → h
exponentially, as ∞ → t , and could be used by placing it in the form given in
equation (2.7). It should be pointed out that the force
c
F given by equation (2.9)
minimizes, at each instant of time, the quantity
c T c
F M F
1
) (
−
--the weighted norm of
the active control force
c
F [63].

The general results obtained in analytical mechanics (see Ref. [63] for more details)
are far more extensive than those presented above; here, we have particularized them
to only cover the present problem of interest—the following object (see Ref. [58] for
a more extensive treatment).

2.4. Equations of Motion of the Following Object (SO)

Consider the unconstrained motion of the first object (FO) the equations of motion
can be put in the form of equation (2.4) so that

⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
−
=
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
g x
x
x
m
m
m
0
0
0 0
0 0
0 0
3
2
1
1
1
1
& &
& &
& &
, (2.11)

where
1
m is the mass of the FO and g is the gravitational constant. When we
consider the constraint equations (2.1)-(2.3) of the FO, by the equation (2.7), we
obtain
30

⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
− − −
− −
− −
=
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
) (
) (
) (
1 0 0
0 1 0
0 0 1
3 3
2 2 2
1 1 1
3
2
1
h x x
v x
v x
x
x
x
x
x
γ β
α
α
&
&
&
& &
& &
& &
, (2.12)

where again
1
α ,
2
α , β , and γ are non-zero constant. Using the equation (2.9) the
constrained equations of motion of the system can be represented by

⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
− − −
− −
− −
=
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
) (
) (
) (
3 3
2 2 2
1 1 1
3
2
1
h x x
v x
v x
x
x
x
x
x
γ β
α
α
&
&
&
& &
& &
& &
. (2.13)

Our aim is to have the motion of the second object (SO) ‘pursue’ the motion of the
FO. Hence we require,

) ( ) ( t y t x
i i
= , , 3 , 2 , 1 = i (2.14)

where ) (t x
i
and ) (t y
i
are the positions of the two different objects at time t . We
note that the equation set (2.13) constitutes a set of independent conditions. The
problem of the following object can be interpreted as one of ensuring that the
tracking conditions (2.14) are satisfied by the equations of motion of the FO.
Alternatively, we think of this problem as one in which the equations of motion of
the FO represent an unconstrained mechanical system on which the three
31
independent constraints (2.14) are required to be imposed. In fact, we can modify
this set of constraints to

0 )) ( ) ( ( ) ( = − = t y t x t h
i i i
, . 3 , 2 , 1 = i (2.15)

Enforcing these constraints would make the motion of the two objects identical.
Noting that in general the initial conditions of the FO may not satisfy the constraints
(2.14) (or alternatively (2.15)), we further modify the constraints (2.15) to

, 0 = + +
i i i
kh h h
& & &
δ . 3 , 2 , 1 = i (2.16)

where δ and k are positive constants [58]. Since the solution of the set of equations
given by (2.16) satisfies the condition that 0 →
i
h as ∞ → t , we have asymptotic
(and exponential) convergence towards the satisfaction of the constraints (2.15), and
hence obtain asymptotically identical tracks of the two different objects.

It is important to point out that by altering the parameters δ and k in equation (2.16)
one can describe different ‘paths’ taken by the system of the objects towards their
eventual harmonization. For simplicity, we have chosen the same constants δ and k
for each equation of the set (2.16). In general, we could have used different values of
δ and k for the different directions, since the values of δ and k for each of the
equations in the set (2.16) control the rate and nature of convergence of ) (t h
ij
to zero.
32

Equations (2.16) can be put in the form of equation (2.7). Since we want to remain
‘untouched’ motion for the FO, the equation (2.13) also should be included. So that

) , ( : q q b kAq q A q A & & & & = − − = δ (2.17)

the matrix A is a 6 by 6 matrix, containing 0’s, 1’s and -1’s and b is a 6 by 1 matrix.
In this three dimensional situation, the 6 by 6 matrix A takes the form

⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
−
−
−
=
1 0 0 1 0 0
0 1 0 0 1 0
0 0 1 0 0 1
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
A , (2.18)

and the 6 by 1 matrix b can be obtained as

⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
− − − −
− − − −
− − − −
− − −
− −
− −
=
) ( ) (
) ( ) (
) ( ) (
) (
) (
) (
3 3 3 3
2 2 2 2
1 1 1 1
3 3
2 2 2
1 1 1
y x k y x
y x k y x
y x k y x
h x x
v x
v x
b
x
x
& &
& &
& &
&
&
&
δ
δ
δ
γ β
α
α
. (2.19)

33
Note that the first three rows of the matrix A and b imply the motion of the FO will
remain ‘untouched’ and the last three rows inform that the set of equation (2.13) is
satisfied. When we consider that the matrix M that describes the unconstrained
motion of the mechanical system consisting of the two objects is given by

⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
=
2
2
2
1
1
1
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
m
m
m
m
m
m
M , (2.20)

where
1
m is the mass of the FO and
2
m is the mass of the SO. From equation (2.9),
the explicit generalized control force,
c
F , is required to enforce the constraint set
(2.17) and the matrix A given in (2.18), we easily determine (this can be done using
Matlab or Maple)

⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
−
−
−
=
+ −
2 2
2 2
2 2
1
1
1
2 / 1
1
0 0
1
0 0
0
1
0 0
1
0
0 0
1
0 0
1
0 0 0
1
0 0
0 0 0 0
1
0
0 0 0 0 0
1
) (
m m
m m
m m
m
m
m
AM , (2.20)
34

that when substituted in relation (2.9) will yield the explicit control forces to exactly
satisfy the constraint equations (2.17), or alternatively (2.16). One thing we should
note that, since
2 / 1 −
AM has full rank and
1 −
= A A , in this system we have
MA A M M AM M = =
− + − 1 2 / 1 2 / 1 2 / 1 2 / 1
) ( so that the equation (2.9) can be rewritten as

) ( ) , , (
1
f AM b MA t q q F
c −
− = & (2.21)

In equation (2.21), the 6 by 1 vector, f M
1 −
, can be determined as

⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
−
− − −
− −
− −
=
−
g
h x x
v x
v x
f M
x
x
0
0
) (
) (
) (
3 3
2 2 2
1 1 1
1
γ β
α
α
&
&
&
, (2.22)

so that using equations (2.18), (2.19), and (2.22), the generalized control force,
c
F ,
is given by

35
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
+ − − − − + −
− − − + −
− − − + −
=
g h x x y x k y x
v x y x k y x
v x y x k y x
M F
x
x
c
) ( ) ( ) (
) ( ) ( ) (
) ( ) ( ) (
0
0
0
3 3 3 3 3 3
2 2 2 2 2 2 2
1 1 1 1 1 1 1
γ β δ
α δ
α δ
& & &
& & &
& & &
. (2.23)

In the equation (2.23), we can easily confirm that there is no control force acting on
the FO so that the FO remains ‘untouched’ after applying the control force to pursue
the FO. Finally, using the equation (2.8), the entire equations of motion of the system
are obtained by

⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
− − − − + −
− − − + −
− − − + −
− − −
− −
− −
=
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
) ( ) ( ) (
) ( ) ( ) (
) ( ) ( ) (
) (
) (
) (
3 3 3 3 3 3
2 2 2 2 2 2 2
1 1 1 1 1 1 1
3 3
2 2 2
1 1 1
3
2
1
3
2
1
h x x y x k y x
v x y x k y x
v x y x k y x
h x x
v x
v x
y
y
y
x
x
x
x
x
x
x
γ β δ
α δ
α δ
γ β
α
α
& & &
& & &
& & &
&
&
&
& &
& &
& &
& &
& &
& &
. (2.24)

As we can see the result of the equation (2.24), the first three equations are identical
to equation (2.13).

2.5. Numerical Example

In this section we consider a numerical example. The example simulates the first
object which is launched in a certain location and the second object which is
36
following the first object, each with its own physical characteristics. The simulation
contains their trajectories, the control force witch is acting on the second object, and
the distance between the two objects. For simplicity, we assume there is no time
delay in the launching time of the first object and the second object and there is no
limitation of the magnitude of the control force.

Consider the first object launched from a certain location, A, and travels to the
calculated region which is located km 100 west and km 1000 north from the location
A with initial velocities hr km x x x / 1
0
3
0
2
0
1
= = = & & & . After some time it will have the
constant velocity,
1
x
v , in
1
x direction (from east to west),
2
x
v , in
2
x direction (from
south to north) and the height, h , in
3
x direction. We assume
1
x
v is hr km / 122 ,
2
x
v
is hr km/ 1220 and h is km 10 and the mass of this object is kg 5000 . At the same
time the second object which has the mass, kg 500 , will launch from the location
which is km 50 west and km 100 north from the location A with initial velocities
hr km y y y / 1
0
3
0
2
0
1
= = = & & & . The following simulation (see Figure 2.1a and Figure 2.1b)
will illustrate this situation. In this simulation we choose
1
α ,
2
α , β , γ , δ , and k
are 1.
37

Figure 2.1a: The solid red line represents the trajectory of the first object and the circle shape blue
line implies the trajectory of the second object within 200 sec.

Figure 2.1b: The dotted blue lines display the time history in ) 3 , 2 , 1 ( = i y
i
direction and the solid
red lines represent the response of the first flying object.

38
0 50 100 150 200
-4
-2
0
2
x 10
-12
F
c
1 on x
1
0 50 100 150 200
-10
-5
0
x 10
-13
F
c
2 on x
2
0 50 100 150 200
-5
0
5
10
x 10
-13
F
c
3 on x
3
0 50 100 150
-3
-2
-1
0
x 10
4
F
c
4 on y
1
0 50 100 150 200
-6
-4
-2
0
x 10
4
F
c
5 on y
2
0 50 100 150 200
0
2000
4000 F
c
6 on y
3

Figure 2.2: The figures on the left hand side inform us of the control forces acting on the first object
and show that the forces are close to zero. The graphs on the right side represent the control forces
exerted on the second object to follow the first one.

Figure 2.3a: The figure shows the distance between two objects in 200 sec. The red circle represents
a starting point of magnitude of distance between two objects. At the beginning the distance
exponentially increases and just after that the distance will decreases dramatically.
39
0 50 100 150 200 250
0
5
10
15
x 10
9
Distance Between x and y
Distance(m)
Time(sec)
125 126 127 128 129 130 131 132 133
0
1000
2000
3000
Distance Between x and y
Distance(m)
Time(sec)

Figure 2.3b: The top figure shows the distance between two objects in 200 sec. The bottom figure
depicts the magnitude of the top figure illustrating that the two objects will be reach within m 3 at
approximately 132 sec.

The Figure 2.2 shows that the control force applied to the first object is close to zero
allowing the first object to fly in an ‘untouched’ condition.

In Figure 2.3a and Figure 2.3b we observe the distance between the two objects
reduce almost exponentially until around 100 seconds, and after 128 seconds the
second object gradually closes to the first object. In our simulation, at about 132
seconds, the distance between the objects is less than m 3 .

2.6 Conclusions

40
Through our examination of the flying object tracking control problem we have
presented an analytical dynamics based approach. This novel approach yields
explicit generalized active control forces where the designated SO can follow and
overlap the motion of the FO. From our findings, we have determined the following
conclusions:

1. The thrust of our strategy involved a two part process: first, we present the flying
object tracking system issue as a tracking control problem, and second, we
reconfigure this tracking control problem as a problem of constrained motion of
nonlinear dynamical systems. Then the motion of the SO is constrained to exactly
follow the motion of the FO, resulting in attaining the exact control forces required
for the tracking motion of SO with the FO. Applying the newly developed general
theory of constrained motion of nonlinear mechanical systems, we obtain the control
forces needed for exact tracking explicitly and in a closed form.

2. In applying our method, control forces for tracking a flying object are found to
have several advantages. The control forces are continuous in time and they are
obtained explicitly in closed form. In addition the forces lead to exact tracking of the
flying object and they provide the minimum forces.

3. By illustrating an example, we show the efficacy of the methodology. We
illustrate a simple model as the first flying object (FO) with realistic assumptions.
The second object (SO) will exactly follow the motion of FO around 132 sec so that
41
the SO can be located within a certain area that we denote m 3 for example. For
simplicity, we assume the SO will start to fly at the same time when the FO launches
and the SO is not greatly restricted in its speed to emphasize the tracking control
forces which have a minimum quantity. In fact, this assumption easily can be
modified from a computation programming point of view. If more detailed
characteristics of the flying object are provided, the simulation will be closer to the
real word circumstances.

4. In our example the FO was designated with certain assumptions. Based on these
assumptions, we obtain the tracking control forces for the SO using fundamental
equations of motion. One of the most important advantages of this procedure is that
the calculation can be done within almost real time. Furthermore, if the location and
speed of the FO can be detected, say, using satellite or other strong observation
equipment, the assumption used in our example for the FO would no linger be
required. So when a little modification in our technique is presented, the tracking
control forces for the SO can be computed within almost real time in this more
realistic situation.
42
CHAPTER 3

A VERIFICATION OF THE RECURSIVE FORMULAE FOR COMPUTING
THE MOORE-PENROSE M-INVERSE OF A MATRIX

3.1. Introduction

The unique generalized inverse of a matrix A was defined by Moore[40] in
1920 and Penrose [43] in 1955 and it is now generally called the Moore-Penrose
(MP) inverse. Because the concept of the generalized inverse extends the idea of the
inverse of a matrix from nonsingular matrices to all matrices, both square and
rectangular, it is useful in numerous application areas, such as, in solving linear
systems of equations [29], in multivariate linear modeling [33], in statistics [56], and
in system identification [32].

There are a number of papers that deal with determining the MP inverse [10,
1, 48, 27, 41, 67, and 35] of a rectangular matrix. Among them, in 1960 Greville [28]
introduced a recursive algorithm for obtaining the generalized inverse of a matrix by
successively adding a column to a matrix, Albert and Sittler [4] derived the time-
recursive method for the least squares inverse in 1965, and Fletcher [19] mentioned a
recursive deletion algorithm in 1966. The Gerville scheme, in particular, can update
the generalized inverse of matrix whenever new information added, and that is why
it is widely extended in the field of statistical inference [46], filtering theory [31],
43
estimation theory [26], system identification [34 and 69], network theory, active
noise and control [11 and 68], linear associative memories [49], information theory
[55], optimization theory, and, more recently, analytical dynamics [63]. His scheme
is referred to in numerous specialized textbooks which concentrate on generalized
inverses [2, 3, 45, 57, and 25].

In 1997, Udwadia and Kalaba [62] provided an alternative and simple proof
of the Greville’s scheme. They also provided a constructive proof for other types of
generalized inverses [61, 60, and 64]. Recently, some improved recursive methods
for obtaining MP inverse have been obtained [69, 49, and 70].

Then, the concept of the generalized MP M-inverse of a matrix, a more general type
of MP inverse, was introduced in early 1970s [47]. This MP M-inverse has recently
found its new application in analytical dynamics [59]. It appears explicitly in the
general equations of motion for constrained mechanical systems. Because they are
slightly different from those of the regular MP inverse, all four properties of the MP
M-inverse of a matrix, namely B, are provided here

B B BB
M
=
+
, (3.1)

+ + +
=
M M M
B BB B, (3.2)

+
M
BB is symmetric, (3.3)
and
44
( ) B B M
M
+
is symmetric, (3.4)

where
+
M
B is the MP M-inverse of the matrix B. We note that the four properties of
the MP M-inverse and the regular MP inverse are different while the first three are
the same. If we consider a set of linear equations, b Bx ≈ , where B is an m by n
matrix, b is an m-vector, and x is an n-vector, the MP M-inverse of the matrix B is
the matrix such that the solution x, uniquely given by b B x
M
+
= , minimizes both the
least squares solution,
2
b Bx − , and the M-norm of the solution x,
2
M
x .

Recently in 2005, the recursive formulae of the MP M-inverse were obtained
by Udwadia and Phohomsiri [66]. They provide four formulae depending on
whether a row vector or a column vector is appended to a matrix. Their proof is
constructive. The formulae that they obtain are purposely set out to obtain the MP M-
inverse of a matrix recursively (either when a row vector or a column vector is
added). In this chapter, we provide an alternative proof through direct verification of
their formulae by using the above-mentioned four properties of MP M-inverse.

3.2 Verification

3.2.1. The Formulae of the Generalized M-Inverse of a Partitioned Matrix by
Adding a Row Vector.

45
Result 3.1.

For any given matrix

⎥
⎦
⎤
⎢
⎣
⎡
=
a
A
B
ˆ
ˆ
, (3.5)

its MP M-inverse formulae are given by

[]
+ + + +
+
+
− =
⎥
⎦
⎤
⎢
⎣
⎡
=
M M M M
M
M
c A a c A
a
A
B
ˆ
ˆ
ˆ
ˆ
ˆ
, when 0 ≠ c , (3.6)

⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
+ +
− =
+ +
+
T
T
M
T
T
M
M
uu
u A
uu
u u A
A
1
ˆ
1
ˆ
ˆ
, when 0 = c , (3.7)

where A
ˆ
is an 1 − m by n matrix, a ˆ is an additional n row-vector, )
ˆ ˆ
( ˆ A A I a c
M
+
− =
and
+
=
M
A a u
ˆ
ˆ .

Note that when 0 ≠ c , the added row vector is linearly independent of the rows of A,
and when 0 = c , the added row vector is linearly dependent of the rows of A [66].
We shall then consider each case separately: when 0 ≠ c and when 0 = c .

Proof

46
Case 1: (when 0 ≠ c )

Because
+
M
BB and B B
M
+
are repeatedly utilized to verify all four properties
of the MP M-inverse, we shall first evaluate
+
M
BB and B B
M
+
. Using equations (3.5)
and (3.6), we have

[]
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
−
−
= −
⎥
⎦
⎤
⎢
⎣
⎡
=
+
+
+ + +
+ + +
+ + + + +
M
M
M M M
M M M
M M M M M
c a
c A
A a c a A a
A a c A A A
c A a c A
a
A
BB
ˆ
ˆ
ˆ
ˆ ) ˆ (
ˆ
ˆ
ˆ
ˆ )
ˆ
(
ˆ ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
. (3.8)

Since 0
ˆ
=
+
M
c A and 1 =
+
M
ac (see Properties 2 and 3 in Appendix B), we have

⎥
⎦
⎤
⎢
⎣
⎡
=
⎥
⎦
⎤
⎢
⎣
⎡
−
=
+
+ +
+
+
1
0
0
ˆ ˆ
1
0
ˆ
ˆ
ˆ
ˆ
ˆ ˆ
M
M M
M
M
A A
A a A a
A A
BB . (3.9)

Using equations (3.5) and (3.6), we get

[][ ] a c A A a c A A
a
A
c A a c A B B
M M M M M M M M M
ˆ
ˆ ˆ
ˆ
ˆ ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
+ + + + + + + + +
+ − =
⎥
⎦
⎤
⎢
⎣
⎡
− = . (3.10)

Since )
ˆ ˆ
( ˆ A A I a c
M
+
− = , we obtain

]
ˆ ˆ
[ )]
ˆ ˆ
( ˆ
ˆ ˆ
[ c c A A I A A a c A A B B
M M M M M M
+ + + + + +
+ = − − = . (3.11)
47

We now verify the four properties of the MP M-inverse.

MP M-inverse Condition 1: B B BB
M
=
+

Using equations (3.5) and (3.11), we obtain

[]
⎥
⎦
⎤
⎢
⎣
⎡
+
+
= +
⎥
⎦
⎤
⎢
⎣
⎡
= =
+ +
+ +
+ + + +
c c a A A a
c c A A A A
c c A A
a
A
B B B B BB
M M
M M
M M M M
) ˆ (
ˆ ˆ
ˆ
)
ˆ
(
ˆ ˆ ˆ
ˆ ˆ
ˆ
ˆ
) ( . (3.12)

Using the fact that A A A A
M
ˆ ˆ ˆ ˆ
=
+
, 0
ˆ
=
+
M
c A , 1 ˆ =
+
M
c a (see Property 2 and 3 in Appendix
B), and )
ˆ ˆ
( ˆ A A I a c
M
+
− = , we have

B
a
A
A A a a A A a
A
A A I a A A a
A
c A A a
A
B BB
M M M M M
M
=
⎥
⎦
⎤
⎢
⎣
⎡
=
⎥
⎦
⎤
⎢
⎣
⎡
− +
=
⎥
⎦
⎤
⎢
⎣
⎡
− +
=
⎥
⎦
⎤
⎢
⎣
⎡
+
=
+ + + + +
+
ˆ
ˆ
ˆ
ˆ ˆ
ˆ ˆ
ˆ
ˆ
)
ˆ
( ˆ
ˆ ˆ
ˆ
ˆ
ˆ ˆ
ˆ
ˆ
. □

MP M-inverse Condition 2:
+ + +
=
M M M
B BB B

Using equations (3.6) and (3.9), we have

48

[]
[].
ˆ ˆ ˆ
ˆ
ˆ ˆ ˆ
1
0
0
ˆ ˆ
ˆ
ˆ
ˆ
) (
+ + + + + +
+
+ + + + + + + +
− =
⎥
⎦
⎤
⎢
⎣
⎡
− = =
M M M M M M
M
M M M M M M M M
c A A A a c A A A
A A
c A a c A BB B BB B
(3.13)

Because
+ + +
=
M M M
A A A A
ˆ ˆ ˆ ˆ
, we obtain

+ + + + + + +
= − =
M M M M M M M
B c A a c A BB B ]
ˆ
ˆ
ˆ
[ . □

MP M-inverse Condition 3:
+
M
BB is symmetric.

Since
+
M
A A
ˆ ˆ
is symmetric, by equation (3.9) we obtain,

T
M
T
M M
M
BB
A A A A
BB ) (
1
0
0
ˆ ˆ
1
0
0
ˆ ˆ
+
+ +
+
=
⎥
⎦
⎤
⎢
⎣
⎡
=
⎥
⎦
⎤
⎢
⎣
⎡
= . □

MP M-inverse Condition 4: ( ) B B M
M
+
is symmetric.

Since A A M
M
ˆ ˆ
+
and c Mc
M
+
are symmetric, using equation (3.11) we have

[ ][]
T
M
T
M M M M M
B B M c c M A A M c c M A A M B B M ) ( ) ( )
ˆ ˆ
( ) ( )
ˆ ˆ
( ) (
+ + + + + +
= + = + = . □

49
As shown, all four MP M-inverse conditions are satisfied. Hence, the formula (3.6) is
verified. □

Case 2 : (when 0 = c )

As before,
+
M
BB and B B
M
+
are frequently utilized to verify all four properties
of the MP M-inverse, we shall first evaluate
+
M
BB and B B
M
+
. By equations (3.5) and
(3.7) we have

⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
+
+
+
−
+
−
=
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
+ +
−
⎥
⎦
⎤
⎢
⎣
⎡
=
+
+
+
+
+
+
+ +
+ +
T
T
M
T
T
M
T
T
M
M
T
T
M
M
T
T
M
T
T
M
M M
uu
u A a
uu
u A A
uu
u u A a
A a
uu
u u A A
A A
uu
u A
uu
u u A
A
a
A
BB
1
)
ˆ
ˆ (
1
)
ˆ ˆ
(
1
)
ˆ
ˆ (
)
ˆ
ˆ (
1
)
ˆ ˆ
(
ˆ ˆ
1
ˆ
1
ˆ
ˆ
ˆ
ˆ
.
(3.14)

Since,
+
=
M
A a u
ˆ
ˆ and
T T
M
u u A A =
+
ˆ ˆ
(see Property 4 in Appendix B), we get

⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
+
+
+
+
−
=
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
+
+
+
−
+
−
=
+ +
+
T
T
T
T
T
T
T
M
T
T
T
T
T
T
T
T
M
M
uu
uu
uu
u
uu
u
uu
u u
A A
uu
uu
uu
u
uu
u uu
u
uu
u u
A A
BB
1
1
1
1
ˆ ˆ
1
1
1
1
ˆ ˆ
. (3.15)

Using relations (3.5), (3.7), and the fact that A u A A a a
M
ˆ ˆ ˆ
ˆ ˆ = =
+
(see Property 5 in
Appendix B), we obtain
50

A A
uu
a u A
uu
A u u A
A A
a
A
uu
u A
uu
u u A
A B B
M T
T
M
T
T
M
M T
T
M
T
T
M
M M
ˆ ˆ
1
ˆ
ˆ
1
)
ˆ
(
ˆ
ˆ ˆ
ˆ
ˆ
1
ˆ
1
ˆ
ˆ
+
+ +
+
+ +
+ +
=
⎥
⎦
⎤
⎢
⎣
⎡
+
+
+
− =
⎥
⎦
⎤
⎢
⎣
⎡
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
+ +
− = .(3.16)

We next verify the four properties of the MP M-inverse.

MP M-inverse Condition 1: B B BB
M
=
+

Since A A a a
M
ˆ ˆ
ˆ ˆ
+
= and A A A A
M
ˆ ˆ ˆ ˆ
=
+
, (see Property 5 in Appendix B), using equations
(3.5) and (3.16) we get

[] B
a
A
A A a
A A A
A A
a
A
B B B B BB
M
M
M M M
=
⎥
⎦
⎤
⎢
⎣
⎡
=
⎥
⎦
⎤
⎢
⎣
⎡
=
⎥
⎦
⎤
⎢
⎣
⎡
= =
+
+
+ + +
ˆ
ˆ
ˆ ˆ
ˆ
ˆ ˆ ˆ
ˆ ˆ
ˆ
ˆ
)(. □

MP M-inverse Condition 2:
+ + +
=
M M M
B BB B

Using equations (3.7) and (3.16), we obtain

[]
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
+ +
− = =
+ +
+ + + + + +
T
T
M
T
T
M
M M M M M M
uu
u A
uu
u u A
A A A B B B BB B
1
ˆ
1
ˆ
ˆ ˆ ˆ
) ( . (3.17)

51
Since
+ + +
=
M M M
A A A A
ˆ ˆ ˆ ˆ
, we have

.
1
ˆ
1
ˆ
ˆ
1
ˆ ˆ ˆ
1
ˆ ˆ ˆ
ˆ ˆ ˆ
+
+ +
+
+ + + +
+ + + +
=
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
+ +
− =
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
+ +
− =
M
T
T
M
T
T
M
M
T
T
M M
T
T
M M
M M M M
B
uu
u A
uu
u u A
A
uu
u A A A
uu
u u A A A
A A A BB B
□

MP M-inverse Condition 3:
+
M
BB is symmetric.

Since
+
M
A A
ˆ ˆ
is symmetric, using the equation (3.15) we get

()
T
M
T
T
T
T
T
T
T
T
M
T
T
T
T
T
T
T
M
M
BB
uu
uu
uu
u
uu
u
uu
u u
A A
uu
uu
uu
u
uu
u
uu
u u
A A
BB
+
+ +
+
=
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
+
+
+
+
−
=
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
+
+
+
+
−
=
1
1
1
1
ˆ ˆ
1
1
1
1
ˆ ˆ
.
□

MP M-inverse Condition 4: ( ) B B M
M
+
is symmetric.

Since A A M
M
ˆ ˆ
+
is symmetric, by the equation (3.16) we observe

() ( )
T
M
T
M M M
B MB A A M A A M B MB
+ + + +
= = =
ˆ ˆ ˆ ˆ
. □

52
As shown, all four MP M-inverse conditions are satisfied. Hence, the formula (3.7) is
verified. □

3.2.2. The Formulae of the Generalized M-Inverse of a Partitioned Matrix by
Adding a Column Vector.

Result 3.2.

For any given matrix

[ ] a A B =, (3.18)

its MP M-inverse formulae are given by

[]
⎥
⎦
⎤
⎢
⎣
⎡ − −
= =
+
+ + + +
+ +
− −
d
pd ad A A
a A B
M M
M M
, when 0 ≠ d , (3.19)

⎥
⎦
⎤
⎢
⎣
⎡ − −
=
+ +
− −
h
ph ah A A
M M
, when 0 = d , (3.20)

53
where A is an m by (n-1) matrix, a is an additional column vector of m
components, , a AA I d
M
) (
+
−
− = , MU q h
T
β
1
= , Mq q
T
= β ,
⎥
⎦
⎤
⎢
⎣
⎡
=
+
−
xm
M
A
U
1
0
,
⎥
⎦
⎤
⎢
⎣
⎡
−
+
=
1
p v
q , m M A A I p
M
~
) (
1 −
−
+
−
− = and a A v
M
+
−
= . Also, we denote

⎥
⎦
⎤
⎢
⎣
⎡
=
−
m
m
m
M
M
T
~
~
. (3.21)

where the matrix M is a positive definite n by n matrix,
−
M is the positive definite
(n-1) by (n-1) principal minor of M, m
~
is the column vector of (n-1) components, and
m is the scalar.
It should be noted that when 0 ≠ d , the added column vector is not a linear
combination of the column of A, and when 0 = d , the added column vector is a
linear combination of the columns of A [66]. As before we shall obtain the formulae
for each case: when 0 ≠ d and when 0 = d .

Proof

Case 1: (when 0 ≠ d )

54
Once more the
+
M
BB and B B
M
+
are continually used to validate all four
properties of the MP M-inverse, we shall first evaluate
+
M
BB and B B
M
+
. By equations
(3.18) and (3.19), we have

[]
+ + + + +
+
+ + + +
+
+ − − =
⎥
⎦
⎤
⎢
⎣
⎡ − −
=
− −
− −
ad d Ap ad AA AA
d
pd ad A A
a A BB
M M
M M
M
) ( .
(3.22)

Since 0 = Ap (see Property 6 in Appendix B) and a AA I d
M
) (
+
−
− = , we get

+ + + + + + + + + +
+ = − − = + − =
− − − − −
dd AA ad I AA AA ad ad AA AA BB
M M M M M M
) ( .(3.23)

By equations (3.18)and (3.19), we obtain

[] a A
d
pd ad A A
B B
M M
M ⎥
⎦
⎤
⎢
⎣
⎡ − −
=
+
+ + + +
+
− −
, (3.24)

⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
− − − −
=
+
+ + + +
+
+ + + +
+
− − − −
a d
a d p a d a A a A
A d
A d p A d a A A A
B B
M M M M
M
) ( ) ( ) ( ) (
.(3.25)

Since, 0 =
+
A d and 1 =
+
a d (see Properties 9 and 10 in Appendix B), we have

55

⎥
⎦
⎤
⎢
⎣
⎡ −
=
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
− −
=
+ + + +
+
− − − −
1 0 1 0
p A A p a A a A A A
B B
M M M M
M
. (3.26)

We now verify the four properties of the MP M-inverse.

MP M-inverse Condition 1: B B BB
M
=
+

Using equations (3.18) and (3.23), we obtain

( ) [ ]
. ] ) ( ) ( [
) (
a d d a AA A d d A AA
a A dd AA B BB B BB
M M
M M M
+ + + +
+ + + +
+ + =
+ = =
− −
−
(3.27)

Because A A AA
M
=
+
−
, 0 =
+
A d , 1 =
+
a d (see properties 9 and 10 in Appendix B) and
a AA I d
M
) (
+
−
− = , we have

[ ] B a A a AA I a AA A B BB
M M M
= = − + =
+ + +
− −
] ) ([. □

MP M-inverse Condition 2:
+ + +
=
M M M
B BB B

Using equations (3.19) and (3.23), we get

56
() []
+ +
+
+ + + +
+ + + +
+
⎥
⎦
⎤
⎢
⎣
⎡ − −
= =
−
− −
dd AA
d
pd ad A A
BB B BB B
M
M M
M M M M
, (3.28)

() ( ) ( )
()
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
+
− − − − +
=
+ + + +
+ + + + + + + + + + + + +
−
− − − − − − −
dd d A A d
dd pd A A d p ad A A A d a A d d A AA A
M
M M M M M M M
.
(3.29)

Since 0 =
+
A d , 1 =
+
a d , 0 =
+
−
d A
M
(see Properties 8 and 9 in Appendix B),
+ + +
= d dd d and,
+ + +
− − −
=
M M M
A AA A , we then have

+
+
+ + + +
+ +
=
⎥
⎦
⎤
⎢
⎣
⎡ − −
=
− −
M
M M
M M
B
d
pd ad A A
BB B . □

MP M-inverse Condition 3:
+
M
BB is symmetric.

Since
+
−
M
AA and
+
dd are symmetric, using equation (3.23), it can be shown that

( ) ( )
T
M
T
M M M
BB dd AA dd AA BB
+ + + + + +
= + = + =
− −
. □

MP M-inverse Condition 4: ( ) B B M
M
+
is symmetric

Using equations (3.20) and (3.26), we observe

57
()
⎥
⎦
⎤
⎢
⎣
⎡ −
⎥
⎦
⎤
⎢
⎣
⎡
=
+
− +
−
1 0
~
~
p A A
m
m
m
M
B B M
M
T M
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
+ −
+ −
=
−
+
+
−
−
−
m p m
m p M
A A m
A A M
T
M
T
M
~
~
~
⎥
⎦
⎤
⎢
⎣
⎡
=
2 , 2 1 , 2
2 , 1 1 , 1
E E
E E
.
(3.30)

where
1 , 1
E ,
2 , 1
E ,
1 , 2
E , and
2 , 2
E represent the elements (1,1), (1,2), (2,1), and (2,2) of
the matrix ( ) B B M
M
+
, respectively, Note that
1 , 1
E is the (n-1) by (n-1) matrix,
2 , 1
E is
the column vector of (n-1) components,
1 , 2
E is the row vector of (n-1) components
and
2 , 2
E is scalar. As we can see,
1 , 1
E = A A M
M
+
−
−
is symmetric and
2 , 2
E is scalar
which is symmetric. Thus, for ( ) B B M
M
+
to be symmetric, we need to show that
2 , 1
E
is the transpose of
1 , 2
E .
Using m M A A I p
M
~
) (
1 −
−
+
−
− = , the element
2 , 1
E can be written as

m A A m M A A M m m M I A A M m p M
T
M M M
~
) (
~
) (
~ ~
) (
~ 1 1 + −
−
+
−
−
−
+
− −
− − −
= = + − = + − ,
(3.31)

which is the transpose of the element of
1 , 2
E . Hence, ( ) B B M
M
+
is symmetric. □

As shown, all four MP M-inverse conditions are satisfied. Hence, formula (3.19) is
verified. □

Case 2 : (when 0 = d )
58

Again, because the forms of
+
M
BB and B B
M
+
are repetitively used for
verifying all four properties of the MP M-inverse, we shall first evaluate
+
M
BB and
B B
M
+
. Using equations (3.18) and (3.20), we have

[] [] ah h Ap h a AA AA
h
ph ah A A
a A BB
M M
M M
M
+ − − =
⎥
⎦
⎤
⎢
⎣
⎡ − −
=
+ +
+ +
+
− −
− −
) ( ) (.
(3.32)

Since a a AA
M
=
+
−
and 0 = Ap (see Properties 6 and 7 in Appendix B), we get

+ + +
− −
= + − =
M M M
AA ah ah AA BB ] [ . (3.33)

Again, using equations (3.18) and (3.20), we have

[]
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
− − − −
=
⎥
⎦
⎤
⎢
⎣
⎡ − −
=
+ + + + + +
+
− − − − − −
ha
pha aha A a A
hA
phA ahA A A A
a A
h
ph ah A A
B B
M M M M M M
M
.
(3.34)

We next verify the four properties of the MP M-inverse.

MP M-inverse Condition 1: B B BB
M
=
+

59
Using the fact that A A AA
M
=
+
−
and a a AA
M
=
+
−
(see Property 7 in Appendix B), by
equations (3.18) and (3.33), we obtain

() [ ] [] B a A a AA A AA a A AA B BB B BB
M M M M M
= = = = =
+ + + + +
− − −
] [ . □

MP M-inverse Condition 2:
+ + +
=
M M M
B BB B

By equation (3.20) and (3.33), we have

() []
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡ + −
=
⎥
⎦
⎤
⎢
⎣
⎡ − −
= =
+
+ + +
+
+ +
+ + + +
−
− − −
−
− −
M
M M M
M
M M
M M M M
hAA
hAA p a A A
AA
h
ph ah A A
BB B BB B
) )( (
.
(3.35)

Because
+ + +
− − −
=
M M M
A AA A and h hAA
M
=
+
−
(see Property 14 in Appendix B), we get

+
+ +
+ +
=
⎥
⎦
⎤
⎢
⎣
⎡ − −
=
− −
M
M M
M M
B
h
ph ah A A
BB B . □

MP M-inverse Condition 3:
+
M
BB is symmetric.

Because
+
−
M
AA is symmetric, by equation (3.33), we have

60
() ( )
T
M
T
M M M
BB AA AA BB
+ + + +
= = =
− −
. □

MP M-inverse Condition 4: ( ) B B M
M
+
is symmetric.

Using equations (3.20), (3.34), and a A v
M
+
−
= , we obtain

⎥
⎦
⎤
⎢
⎣
⎡ − − − −
⎥
⎦
⎤
⎢
⎣
⎡
=
+
− +
−
ha
pha vha v
hA
phA vhA A A
m
m
m
M
B B M
M
T M
~
~
) ( ,

⎥
⎦
⎤
⎢
⎣
⎡
=
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
+ − −
+ − −
+ − −
+ − −
=
−
+
+
−
−
−
2 , 2 1 , 2
2 , 1 1 , 1
) ( ) (
~
) (
~
) (
) ( ) (
~
) (
~
) (
E E
E E
ha m pha vha v m
ha m pha vha v M
hA m phA vhA A A m
hA m phA vhA A A M
T
M
T
M
.
(3.36)

where
1 , 1
E ,
2 , 1
E ,
1 , 2
E , and
2 , 2
E represent the elements (1,1), (1,2), (2,1), and (2,2) of
the matrix ( ) B B M
M
+
, respectively. Note that
1 , 1
E is the (n-1) by (n-1) square matrix,
2 , 1
E is the column vector of (n-1) components,
1 , 2
E is the row vector of (n-1)
components and
2 , 2
E is scalar. For ( ) B B M
M
+
to be symmetric, we need to show that
1 , 1
E is symmetric,
2 , 1
E is the transpose of
1 , 2
E , and
2 , 2
E is symmetric.

Let us first show
1 , 1
E is symmetric. We can rewrite
1 , 1
E as

61

1 , 1
E = ) )(
~
( hA m p M v M A A M
M
− + −
− −
+
−
−
, (3.37)

Using the fact that

)
~
(
1
A A m M v hA
M
T T +
−
−
− =
β
(3.38)

and

() m A A m p M
T
M
~ ~ +
−
−
− = − (3.39)

(see Property 13 and 16 in Appendix B), equations (3.37) can be written as

1 , 1
E = ( ) ]
~
[ ]
~
[
1
A A m M v m A A v M A A M
M
T T
T
M M
+
−
+
−
+
−
− − −
− − −
β
. (3.40)

Since ()
T
M M
A A M A A M
+
−
+
−
− −
= and

() ]
~
[ ]
~
[ A A m M v m A A v M
M
T T
T
M
+
−
+
−
− −
− −
{ }
T
M
T T T
M
A A m M v m A A v M ]
~
[ ]
~
) ( [
+
−
+
−
− −
− − = , (3.41)

the first and the second term on the right-hand side of equation (3.40) are symmetric.
62
Next, we will show that
2 , 1
E is transpose of
1 , 2
E . Let us rewrite
2 , 1
E , as

2 , 1
E = ) )(
~
( ha m p M v M v M − + −
− − −
. (3.42)

Using equation (3.39) and ( ) v m v M v ha
T T ~
1
− =
−
β
(see Property 15 in Appendix B), in
equation (3.42), we obtain

2 , 1
E = ( )
β
1
]
~
][
~
[ v m v M v m A A v M v M
T T
T
M
− − −
−
+
− −
−
. (3.43)

On the other hand,
1 , 2
E , can be written as

1 , 2
E = ) )(
~ ~
( ) (
~
hA m p m v m A A m
T T
M
T
− + −
+
−
. (3.44)

Using the property β − − = − +
−
m v v M v m p m v m
T T T T ~ ~ ~
(see Property 17 in Appendix
B) and equation (3.38) in equation (3.44), we get

1 , 2
E =
β
β
1
]
~
][
~
[ ) (
~
A A m M v m v v M v A A m
M
T T T T
M
T +
− −
+
− −
− − − − , (3.45)

which can be simplified to

63

1 , 2
E =
β
1
]
~
][
~
[ A A m M v m v v M v M v
M
T T T T T +
− − −
−
− − − . (3.46)

It can be seen from equations (3.43) and (3.46)
2 , 1
E is the transpose of
1 , 2
E . Since
we have already shown that
1 , 1
E is symmetric and since
2 , 2
E is scalar which is
symmetric, the symmetry of ( ) B B M
M
+
is verified. □

As shown, all four M M-inverse conditions are satisfied. Hence, formulae (3.20) is
verified. □

3.3 Conclusions

The recursive formulae for the determination of the Moore-Penrose (MP) M-
inverse of a given matrix were first given by [66]; however, they do not provide any
verification of their formulae. In this chapter, it is shown that their recursive
formulae satisfy all four conditions of the MP M-inverse. Hence, the validity of the
recursive formulae is confirmed.

64
CHAPTER 4

A DYNAMICAL MODEL OF TERRORISM

4.1. Introduction

The 21
st
century has been marked by a new kind of warfare—terrorism. The roots of
terrorism can be identified as lying in economic, religious, psychological,
philosophical, and political aspects of society. Incidents of terrorism can often be
sparked by the deterioration of some local conditions, in a spatial sense, as perceived
by a small segment of society. Several qualitative models of terrorism are related to
different ways of prioritization of terrorist targets, security measures available, and
the psychological impact that such activity may have on the local population in a
given area, and in the world at large. Other qualitative models argue in terms of
group behavior that is differentiated as being ideological, grievance-driven, and
understandable [44]. Several attempts to study the phenomenon of terrorism using
the case-study method have been carried out, attempting to look at the similarities
and differences between various acts of terrorism while placing them in a historical
context. Some models that are more quantitative utilize the economic approach using
utility functions to model terrorist behavior, and to assess the losses created by
terrorist actions [9, 13, and 21], while others use optimal control methods to model
governmental actions aimed at maximizing security under the constraints posed by
the optimal trajectories selected by terrorist [17]. While it is somewhat questionable
65
whether such rational behavior can be imputed to fundamentalist extremist groups,
the main difficulties posed by such models appear to be the ad hoc nature of the
utility functions and ‘costs’ assigned to terrorist groups and their pursuers [5 and 8].
Also whether extreme events such as suicide-bombings fall within the framework of
utility theory appears questionable. Still others try to develop static ‘rational-actor’
models for negotiating and bargaining with the demands of transnational terrorists [7,
52, and 53]. Strategic response to terrorist activities using game theoretic approaches
has been looked at by [15, 20, 50, and 51]. Probabilistic assessment of terrorist
activities through the analysis of historical data is another approach that has been
used in modeling terrorist activity. Work on rudimentary models of terrorist
activities using game theory by developing attack-defense strategies when multiple
targets are to be defended under resource constraints has also been initiated [12].
One of the main aims for developing such models of terrorism is the determination
of suitable measures to counteract and control it, and to be able to get a prognosis of
the environment in terms its level of security and safety. Yet there appear to be very
few models that are truly dynamical in nature and which therefore attempt to look at
the time-evolution of terrorist activity in a manner that can be usefully employed to
yield actionable information. (For a thorough overview of the growing literature in
this field, see [17].)

In this paper we present a simple dynamical model of terrorism in terms of the
dynamics of the population of individuals who engage in terrorist activities. We
imagine the population of a certain area (say the Gaza Strip, or the West Bank) as
66
divided into two categories: Terrorists (T) and Non-terrorists (NT). The Non-
terrorists are further divided into those that are susceptible to terrorist propaganda—
this segment of the population (perhaps, the Wahabis in certain areas of the world, or
those educated in madrasaas) we call the Susceptibles (S)—and those that are not
susceptible to such propaganda, which we refer to as Non-susceptibles (NS). We
stipulate a reasonable model for the dynamics which includes the effect of
military/police action against Terrorists and the effect of non-violent means to wean
away the susceptiple population from turning to terrorism. (Another stream of
literature concerns the presence of any causal connections between factors like
education and poverty and the arising and growth of terrorism. See [36].) One of the
main motivations for the development of this dynamical model is the insights it
provides to help understand the dynamical evolution of these different populations,
to understand the different regimes of dynamical behavior that arise, and to point us
in the right direction for asking the proper questions in order to predict and interdict
terrorist activity.

4.2. The Dynamical Model

Let us say that the number of Terrorists (T) in a certain geographical region (say, a
city) at time τ is x( τ ). As mentioned before, we shall think of the non-terrorist
population in the area as being made up of the population of Susceptibles (S), y( τ ),
and of nonsusceptibles (NS), z( τ ) .

67
The number of Terrorists in a given period of time can change because of several
reasons: (1) direct recruitment by the Terrorists of individuals from the Susceptible
population; the effectiveness of this is taken to be proportional to the product of the
number of Terrorists and the number of Susceptibles; (2) effect of anti-terrorist
measures that are directed directly at reducing the terrorist population, such as
military and police action/intervention which we assume increases rapidly with, and
as the square of, the number of Terrorists in the region under concern; (3) number of
Terrorists that die from natural causes, or are killed in action, and/or self-destruct (as
in the case of suicide bombers), which we assume to be proportional to the Terrorist
population itself; and (4) increase in the Terrorist population primarly through the
appeals by Terrorists (in the region under concern) to other terrorist groups, through
global propaganda using news media, and/or through the organized or voluntary
recruitment/movement of Terrorists from other regions into the region of concern,
and also through population growth in this section of the population; this brings
about an increase in the Terrorist population that we assume is proportional to the
number of Terrorists. We capture these four effects then through the following
differential equation that we posit for the evolution of the Terrorist population in the
geographical region of concern:

x c c x b xy a
d
dx
) ˆ ˆ (
ˆ
ˆ
2 1
2
− + − =
τ
, (4.1)

68
where we assume for convenience that the parameters
2 1
ˆ and , ˆ ,
ˆ
, ˆ c c b a are constant
over the time-horizon of interest, and non-negative. We denote time by the parameter
τ . The term containing
2
ˆ c refers to the death/destruction of members of the
Terrorist population either through natural causes or through suicide-bombings, and
the term containing
1
ˆ c refers to their increase either through recruitment from among
their own, or the importation of Terrorists from other geographical areas. The
effectiveness of Terrorists to attract Susceptibles to their cause is described by the
parameter a ˆ , and the effectiveness of military/police action in reducing the numbers
of Terrorists is characterized by the parameter b
ˆ
.

The change in the number of the Susceptibles (S) in a given interval of time is
likewise caused by several factors. (1) Depletion in their population caused by their
direct contact with Terrorists whose point of view they adopt. This is just the number
of Susceptibles that entered the ranks of Terrorists, given before by xy a ˆ . (2)
Depletion in the population of Susceptibles caused by: nonviolent propaganda done
by governmental and non-governmental authorities that convinces members of this
population to use peaceful methods of engagement; and/or concessions made to
disgruntled groups of Susceptibles --these ‘carrots’ offered may be of an economic,
political, or other nature--to convince them to enter the ranks of the non-susceptibe
(NS) population. We model this effect by assuming that the propaganda and/or
concessions are targeted to the Susceptible population, and that this propaganda
intensifies rapidly as the number of Terrorists in the geographical area under concern
69
increases. We assume that the change this causes is proportional to the product y x
2
.
(3) Increase in the population of Susceptibles caused by the propaganda that is
created through the notoriety and publicity of terrorist acts that are broadcast on
global information channels, like television and printed media, that causes some
members of the NS population to become Susceptibles. (4) Increase in the
Susceptible population when individuals from outside the geographical area of
concern are incited to move into the area, first as Susceptibles (S), perhaps later
going on to become Terrorists. We assume that the changes in the S population
attributable to this cause and the previous one are proportional to the number of
Terrorists in the region under concern. (5) The increase in the Susceptible population
proportional to its own size (for example, children of individuals educated in
madrasaas being educated, likewise, in madrasaas). The evolution of the Susceptible
population adduced from these effects can be expressed by the differential equation:

, ˆ )
ˆ ˆ
( ˆ ˆ
2 1
2
y g x f f y x e xy a
d
dy
+ + + − − =
τ
(4.2)

where we again assume that the parameters g f f e ˆ and
ˆ
,
ˆ
, ˆ
2 1
are each a constant and
non-negative over the time-horizons of interest. The parameter e ˆ signifies the
effectiveness of nonviolent means in weaning away Susceptibles into the NS
(pacifist) population. The effect of individuals from the NS population moving to the
S population is given by the term x f
1
ˆ
; the effect of individuals from outside of the
region of concern being attracted to the region and becoming Susceptibles is
70
indicated by the term x f
2
ˆ
. The growth rate of the Susceptible population is given by
g ˆ .

Lastly, the change in the number of nonsusceptibles (NS) in a given interval of time
is described by (1) those members of the Susceptible population that become NS by
virtue of having altered their persuasions because of the nonviolent
actions/propaganda/inducements of governmental and non-governmental authorities,
(2) those who become Susceptibles due to effects of global propaganda done by
Terrorists through news media, and the like, and (3) the increase in the NS
population, which is proportional to their population numbers. This then may be
described by the equation

z h x f y x e
d
dz
ˆ ˆ
ˆ
1
2
+ − =
τ
. (4.3)

where we assume, for simplicity again, that the parameter h
ˆ
, which is the growth
rate of the NS population, is constant and non-negative. The dynamical system is
schematically illustrated in Figure 4.1. For the purposes of our analysis we shall
assume that y x z , >> .

71
x c
1
ˆ +
2
ˆ
x b −
xy a ˆ −
x f
2
ˆ
+
y x e
2
ˆ −
x f
1
ˆ
+
x f
1
ˆ
−
x c
1
ˆ +
y g ˆ +
z h
ˆ
+
Extended
Neighborhood
Nonviolent
Intervention
Military
Intervention
x c
2
ˆ −
T
S
xy a ˆ +
y x e
2
ˆ +
NS
x c
1
ˆ +
2
ˆ
x b −
xy a ˆ −
x f
2
ˆ
+
y x e
2
ˆ −
x f
1
ˆ
+
x f
1
ˆ
−
x c
1
ˆ +
y g ˆ +
z h
ˆ
+
Extended
Neighborhood
Nonviolent
Intervention
Military
Intervention
x c
2
ˆ −
T
S
xy a ˆ +
y x e
2
ˆ +
NS

Figure 4.1: Schematic showing the Dynamical System described by equations (1), (2), and (3).

We begin by dividing equations (4.1)-(4.3) by
2
ˆ c and using the dimensionless time
τ
2
ˆ c t = . This yields the equations

x c bx axy
dt
dx
) 1 (
2
− + − = , (4.4)

,
2
gy fx y ex axy
dt
dy
+ + − − = (4.5)

and,
72

hz x f y ex
dt
dz
+ − =
1
2
(4.6)

where, all the constants are normalized so that , ˆ / ˆ
2
c a a = , ˆ /
ˆ
2
c b b = , ˆ / ˆ
2 1
c c c =
, ˆ / ˆ
2
c e e = , ˆ /
ˆ
2 1 1
c f f = , ˆ /
ˆ
2 2 2
c f f = ), (
2 1
f f f + = , ˆ / ˆ
2
c g g =

and
2
ˆ /
ˆ
c h h = . We
thus have a nonlinear system of three differential equations containing a total of 8
constant parameters all of which we shall assume, for the purposes of this analysis,
to be non-negative.

4.3. Model Dynamics

In this section our aim is to understand the unfolding of the dynamics of the system,
to study its various regimes of behavior in the phase space (x, y, z), and investigate
the manner in which the behavior changes as the values of the 8 parameters change.

We begin by noting that equation (4.6) is uncoupled from equations (4.4) and (4.5),
and hence the entire dynamics of the evolution of the various populations is
dependent only on the latter two equations. Once ) (t x and ) (t y are known, the
dynamics of the NS population, ), (t z is determined from equation (4.6). Since the
system dynamics is then only dependent on the two coupled equations (4.4) and (4.5),
we can rule out the possibility of having populations ) (t x and ) (t y that, in the strict
technical sense, chaotically change with time. Because ) (t x and ) (t y represent the
73
number of Terrorists and the number of Susceptibles, the region in phase space that
is of interest to us is limited to 0 ) ( ), ( ≥ t y t x . We shall now show that in this
quadrant of the phase space the orbits of the nonlinear dynamical system described
by equations (4.4) and (4.5) cannot be closed for 0 , > f b , and hence the nonlinear
system is devoid of any limit cycles.

4.3.1 On the Orbits of the Dynamical System.

Result 4.1: The nonlinear dynamical system described by equations (4.4) and (4.5)
above does not have any limit cycles (closed orbits) in the first quadrant of the phase
plane for 0 , > f b .
Proof: Consider the function

(4.7)

Setting 1 − = = n m , we then obtain

2
) , (
y
f
y
b
y x q − − = . (4.8)

[ ] [ ]
. ) 1 ( ) 1 ( ) 1 (
) 1 )( 1 ( ) 2 ( ) 1 (
} { } ) 1 ( { ) , (
1 1 2 1
1 1
2 2
n m n m n m n m
n m n m n m
n m n m
y x n g y fnx y x n e y x n a
y x m c y x m b y x m a
gy fx y ex axy y x
y
x c bx axy y x
x
y x q
+ + + + − + −

+ − + + − + =
+ + − −
∂
∂
+ − + −
∂
∂
=
− + + +
+ +
74
Since ) , ( y x q is negative in the first quadrant for 0 , > f b , by Dulac’s Criterion
there can be no closed orbits.

Result 4.2: All orbits that start at 0 = t in the first quadrant 0 , ≥ y x remain in that
quadrant for all time 0 > t .
Proof: We begin by noting that the origin of the phase plane is always a fixed point;
we therefore only need to concentrate on the flow on the x- and y-axes of the
positive quadrant. Consider the flow on the x-axis. By equation (4.5) we see that the
flow velocity in the y-direction is given by 0 ≥ = fx
dt
dy
. Since the flow does not
have a negative component of velocity at any point along the positive x-axis, it
cannot cross it. Similarly, along the y-axis, the x-component of the flow velocity is
zero, so the flow cannot leave the first quadrant.

Furthermore, the flow at any point of the positive x-axis is pointed in either the
positive or negative x-direction. When, in addition, 0 = = f b , the flow is pointed in
the positive x-direction for 1 > c and along the negative x-direction for 1 < c ; when
1 = c , the x-axis becomes a line of fixed points. When 0 = f , 0 > b and 1 > c the
flow field at any point on the positive x-axis is pointed in the positive x-direction for
b c x / ) 1 ( − < and in the negative x-direction for b c x / ) 1 ( − > . When 1 = c , the
fixed point at b c x / ) 1 ( − = moves to the origin, and the flow at any point on the
entire positive x-axis is along the negative x-direction. For 1 < c , 0 ≥ b the x-
75
component of the flow velocity at any point of the x-axis is always directed along
negative x-direction.

4.3.2 Fixed Points of the Dynamical System.

We begin by understanding the long-term evolutionary dynamics of the population
of Terrorists (T) and Susceptibles (S) by identifying the fixed points of this nonlinear
dynamical system described by equations (4.4) and (4.5). We observe that the point
0
0 0
= = y x is always a fixed point of the dynamical system. We now look at the
other fixed points of the nonlinear equations (4.4) and (4.5) that lie in the positive
quadrant of the phase space, and to begin with, we differentiate between four
different possible situations. These cases are provided to introduce the different
scenarios of interest in observing the dynamics of terrorism, and later on we shall
take up each of them in greater detail.

Case I: When no intervention (military or nonviolent) is undertaken, that is, when the
parameters . 0 = = e b The fixed point of equation (4.4) and (4.5) is then given by

a
c
y
f c a
c g
x
−
=
− −
−
=
1
,
) 1 (
) 1 (
0 0
. (4.9)

We note that for the fixed point to lie in the first quadrant, we require 1 ≤ c
and 1 ) ( ≤ + c f . We shall treat this case where there is no intervention as the
76
‘baseline situation,’ and in what follows in this section we shall assume 1 < c and
1 < + c f . What happen when 1 ≥ c and/or

1 ≥ + c f will be taken up later on.

Case II: When nonviolent intervention is carried out against Terrorists while
abstaining from military/police intervention, the parameter 0 = b . The fixed point of
the dynamical system is now

a
c
y
c e
c eg c f a c f a
x
−
=
−
− + − + + − +
=
1
,
) 1 ( 2
) 1 ( 4 ) 1 ( ) 1 (
0
2 2 2
0
. (4.10)

We note that the restriction 1 ) ( < + c f is no longer required for
0
x to be positive.
However, as in the previous case, for
0
y to be nonnegative, we require 1 ≤ c . We
observe that the effect of nonviolent intervention does not affect the steady state
value (
0
y ) of the population of Susceptibles when compared to that with no
intervention at all.

Case III: When military/police intervention against terrorism is carried out in the
absence of nonviolent actions, 0 = e , and the fixed point moves to

a
bx c
y
ab
c gab c f a gb c f a gb
x
0
0
2
0
1
,
2
) 1 ( 4 )] 1 ( [ ) 1 ( + −
=
− + − − − + − − −
= . (4.11)

77
Military/police intervention appears to increase the steady state value of the
Susceptible population when compared with that for no intervention at all.

However, this case is a bit more complicated, as we shall see later on, and it could
lead to two fixed points when , 0 > g 0 = f and 1 > c , one of which will be shown
to be unstable.

Case IV: When both military and nonviolent intervention against terrorism is
implemented, that is, when 0 , > e b , the fixed point is located at

a
bx c
y
0
0
1 + −
= , (4.12)

where
0
x is the real positive root of the cubic equation

0 ) 1 ( ] ) 1 ( [ ] ) 1 ( [
2 3
= − − − − − + + − + c g x gb f c a x ab c e ebx . (4.13)

We observe that the addition of nonviolent intervention to military action does not
affect the steady state population of Susceptibles. Also, when 1 < c the first two
coefficients of equation (4.13) (corresponding to the cubic term and the square term)
are always positive, the last coefficient is always negative, and the coefficient of x is
sign indefinite. Hence, from Descarte’s rule of signs there can be only one positive
root of this cubic equation. We will take up the case when 1 > c later on.
78

4.3.3 Stability of Fixed Points.

We next inquire into the stability of the fixed points related to each of the above-
mentioned cases.

Linearization of the nonlinear equations around the fixed point ) , (
0 0
y x so that
) ( ) (
0
t u x t x + = , ) ( ) (
0
t v y t y + = leads to the linearized equations

⎥
⎦
⎤
⎢
⎣
⎡
⎥
⎦
⎤
⎢
⎣
⎡
− − − −
− − −
=
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
v
u
ex ax g exy ay f
ax c bx ay
dt
dv
dt
du
y x ) , (
2
0 0
2
) 1 ( 2
=
⎥
⎦
⎤
⎢
⎣
⎡
v
u
J
y x ) , (
0 0
. (4.14)

In relation (4.14) we have denoted the Jacobian evaluated at the fixed point ) , (
0 0
y x
by
) , (
0 0
y x
J . Since
⎥
⎦
⎤
⎢
⎣
⎡ − −
=
g f
c
J
0 ) 1 (
) 0 , 0 (
, its eigenvalues are ) 1 ( c − − = λ and
g = λ , so that the fixed point (0,0) is an unstable saddle point as long as 0 > g (later
on, we shall briefly consider the case 0 = g which is nonhyperbolic). The
eigenvectors corresponding to the stable and unstable manifolds are then
T
f
g c
⎥
⎦
⎤
⎢
⎣
⎡
−
+ −
1
1
and []
T
1 0 respectively. The y-axis is thus part of the unstable
manifold of the fixed point (0,0).

79
Also, for the fixed point ) , (
0 0
y x , we have, using the equations that govern the fixed
point of the system (4.4) and (4.5),

⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
−
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
+ −
−
=
0
0
0 0
0
0
0 0
) , (
0 0
y
fx
y ex
x
gy
ax bx
J
y x
,
0
x, 0
0
> y . (4.15)

The eigenvalues, λ , of the matrix
) , (
0 0
y x
J are the roots of the characteristic equation

0
0
2
0
0 0
0
0
0 0
0
0 2
= +
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
+ +
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
− − −
y
bfx
y ex
x
gy
ax bx
y
fx
λ λ . (4.16)

Noting that for b, f > 0, the trace of the matrix in (3.15) is negative and its
determinant is positive, we find that the fixed point ) , (
0 0
y x is either a stable spiral
or a stable node.

Thus the fixed point at ) , (
0 0
y x is:

1. a stable spiral when
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
+ <
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
0 0
0
0
0
2
0
0
0
4 y ex
x
gy
ax bx
y
fx
(4.17)
and
2. a stable node when
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
+ >
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
0 0
0
0
0
2
0
0
0
4 y ex
x
gy
ax bx
y
fx
. (4.18)
80

For each of the four cases considered earlier the location of the fixed points is
provided by the relation (4.9)-( 4.13); they are functions of the six parameters that
specify the particular dynamical system being considered. As an example, for our
baseline situation with , 0 = = e b the fixed point given by relation (4.9) is a stable
node when g f c c f / ) 1 )( 1 ( 4
2 2
− − − > , and, a stable spiral when
g f c c f / ) 1 )( 1 ( 4
2 2
− − − < .

It should be observed that the system dynamics look very different when the
parameter 0 = g . For then, the entire line 0 = x constitutes a line of fixed points.
The Jacobian matrix given in relation (4.14), then becomes

) 0 ; , 0 ( = g y
J =
⎥
⎦
⎤
⎢
⎣
⎡
−
− −
0
0 ) 1 (
ay f
c ay
(4.19)

making the stability of these non-isolated fixed points difficult to ascertain. The
eigenvalues of the matrix in (4.19) are 0 = λ , and ) 1 ( c ay − − = λ . And so we may
only conjecture, since the fixed points are nonhyperbolic, that for a c y / ) 1 ( − < the
fixed points along the y-axis in the phase plane are stable (so that 0 < λ ), and that for
a c y / ) 1 ( − > they become unstable. Numerical simulation shows that this
conjecture is correct.

81
We observe in the above analysis (see equation (4.9)) that for the baseline situation,
for which b = e = 0, we assumed in addition that 1 < c and 1 ) ( < + c f . In order to
obtain a better understanding of the ‘baseline situation’ in which there is no
military/police intervention as well as no non-violent intervention, so that we can
later on compare the behavior of the dynamics when we introduce interventions, we
look at this situation here in some greater detail.

4.3.4 Baseline Case (b = e = 0).

We begin by looking at the feasible range of parameters that might be of interest in a
more-or-less realistic situation. We begin with setting the time scale (recall, we are
using τ
2
ˆ c t = to denote our dimensionless time). Imagine a population of 400,000
inhabitants in a certain region (say, the city of Ar Ramadi in Iraq). Let us say that in
a month, 5% of the Terrorists in this city resort to successful suicide bombings. Thus
2
ˆ c =0.05/month, and each unit of dimensionless time, t , then corresponds to an
actual duration, , τ of 20 months ( 05 . 0 / 1 = τ ). In this period of 20 months, let us
assume that, on average, each terrorist persuades 1 member of the S population of
this city from every 10,000 of its members to join the terrorist’s cause. Thus the
parameter . 10 000 , 10 / 1
4 −
= = a Let us say that during these 20 months, the
Susceptible population grows by 20 individuals for every 1,000 members of the S
population, so that 02 . 0 = g . Furthermore, we assume that 5 individuals, on average,
from the Non-Susceptible (NS) population are ‘softened up’ and become
82
Susceptibles by the global notoriety/propaganda (through the news media, say)
created by the collective actions/destruction wreaked by every 100 Terrorists in this
period of time (20 months), so that 05 . 0
1
= f . Also, since the region under
consideration has ‘porous boundaries’ and attracts Terrorists from its neighboring
regions (cities) we assume that the destructive activities of a 100 Terrorists within
the region cause 10 individuals to move from a neighboring region into the area and
join the Susceptibles camp, so that 1 . 0
2
= f . Lastly, we assume that during the 20
month unit of time under consideration, 120 trained Terrorists (on average) move
from outside the area into the area under consideration either voluntarily or through
an organized network of world-wide Terrorists as a consequence of the propaganda
and/or requests made by every 100 Terrorists who reside within the region; this
makes 2 . 1 = c We then see that the parameters , 10
4 −
= a 2 . 1 = c , = + =
2 1
f f f 0.15,
and 02 . 0 = g , may well be within the realm of possibilities. Were the boundaries of
the region to be assumed ‘secure,’ then we might expect only a trickle of individuals
coming in from neighboring regions/cities into our region of concern (say, the city of
Ar Ramadi) so that both the parameters
2
f and c would be much less than unity.
However, in reality, at present it appears unlikely that they would be zero.

Consider the dynamics for 0 = = e b when 1
ˆ
ˆ
2
1
> =
c
c
c , so that the rate of increase of
Terrorists exceeds their destruction/death as seen from the linear term on the right
hand side of equation (4.1). Then the dynamical system, as seen from equation (4.9),
points out that there is now no fixed point in the first quadrant of the phase space.
83
The phase portrait of the nonlinear system can now best be understood geometrically
by noting that equation (4.4) indicates that the x-component of the phase velocity is
always positive in the first quadrant. The nullcline, which describes the curve for
which the y-component of the phase velocity goes to zero is given by the relation
gy fx axy − − = 0, whose slope goes to zero when a f y / = . Above this nullcline the
y-component of the velocity is negative, while below the nullcline it is positive
causing the phase trajectories to asymptote to the line a f y / = . The positive x-
component of the velocity of the phase particle causes x to be unbounded as ∞ → t .
We illustrate this dynamical behavior in Figure 4.2.

Figure 4.2a: Phase plot showing the flow field and the asymptotic behavior of the dynamical system,
1 > c . Trajectories starting from different initial conditions are shown in color. The solid black line is
the nullcline gy fx axy − − = 0 . Notice that the terrorist population increases constantly while the S
population reaches steady state given by 1500 / = = a f y .
84

Figure 4.2b: The time evolution of the dynamical system starting from 2 ) 0 ( = x , 000 , 10 ) 0 ( = y
shows the run-away Terrorist population; the Susceptible population reduces to its asymptotic value
of 1,500.

We observe that when 1 > c , there appears a rapid terrorist expansion; note the
increase in the S population before its sudden drop in Figure 4.2b. When 0 = f , the
S population asymptotically tends to zero.

When 1 < c , but 1 ) ( > + c f , equation (4.9) again points out that there is no fixed
point in the positive first quadrant. This situation might arise when the borders of the
region are maintained relatively secure, so that the influx of Terrorists entering our
region of concern is limited, which nonetheless may have a high number of people
from the NS population (say, disgruntled pacifists who may no longer want to ‘sit on
the sidelines’) who are persuaded to become Susceptibles so that the value of
1
f may
be large. Another instance of such dynamical behavior might arise when the
boundaries are relatively secure from known Terrorists to prevent their free
movement into a region, but there are a substantial number of ‘motivated’
Susceptibles (so that the value of
2
f is now large), who are harder to keep track of,
85
and who enter the geographical region, relatively speaking, unchecked. Figure 4.3
shows the dynamical behavior of the system. We notice that the main geometrical
difference between the phase portraits in Figures 4.2 and 4.3 is the presence of the
nullcline shown in green, which now moves to the positive quadrant. The asymptotic
value of x is unbounded, while that of y again approaches a f / , as illustrated.

Figure 4.3: Phase plot showing the flow field and the asymptotic behavior of the dynamical system,
when 1 , 1 > + < c f c . Trajectories starting from different initial conditions are shown in color. The
solid black line is the nullcline gy fx axy − − = 0 ; the solid green line shows the nullcline
c ay − = 1 . The phase particle has a negative x-component of velocity below this green line,
indicating a drop in the Terrorist population, and a positive component above it. Notice that the
Terrorist population asymptotically increases constantly while the S population again reaches a steady
state given by 1000 / = = a f y .

Figures 4.2 and 4.3 make the bifurcation from unstable to stable behavior clear.
When , 1 > c the green nullcline a c y / ) 1 ( − = , which is horizontal in the phase plane
intersects the y axis at a negative ordinate. Hence there is no fixed point since the
86
other nullcline 0 = − − gy fx axy cannot intersect it in the first quadrant. When 1 = c ,
the horizontal nullcline coincides with the x-axis. As c further decreases this
nullcline moves into the positive quadrant. However, since the asymptote to the
nullcline 0 = − − gy fx axy is at a f y / = , until 1 < + c f , the two nullclines cannot
cross, and hence there is no fixed point. The phase flow indicates that the system is
unstable when 1 > + c f , with the population of Terrrorists continually increasing.
When 1 = + c f , the nullclines cross at ∞ → x , and the system remains unstable.
When 1 < c and 1 ) ( < + c f , equation (4.9) points out that we have a fixed point in
the first quadrant that may be either a stable node, or a stable spiral. The two
nullclines now intersect each other in the first quadrant. We illustrate this in Figure
4.4a, where we show the stable spiral correctly predicted, with 450
0
= x , 4500
0
= y .
Here we take 1 . 0 = c and 7 . 0 = f . Condition (4.17) is satisfied by the parameters
and we obtain a stable spiral. Figure 4.4b shows the dynamics when 7 . 0 = c and
1 . 0 = f , so that for both these simulations . 8 . 0 = + g f We note that an increase in
the value of c paradoxically reduces the asymptotic populations of Terrorists and
Susceptibles to 150
0
= x , 1500
0
= y .

87
Figure 4.4a: Phase plot showing the flow field (with b = e = 0) and the asymptotic behavior of the
dynamical system, when 1 , 1 < + < c f c , with 1 . 0 and , 8 . 0 = = + c c f . The fixed point is at
4500 , 450
0 0
= = y x . Trajectories starting from different initial conditions are shown in color.

Figure 4.4b: Phase plot showing the flow field (with b = e = 0) and the asymptotic behavior of the
dynamical system, when 1 , 1 < + < c f c , with 7 . 0 and , 8 . 0 = = + c c f . The fixed point is at
1500 , 150
0 0
= = y x . Trajectories starting from different initial conditions are shown in color.
The same initial conditions as in Figure 4.4a are used.

Increasing the parameter a by a factor of 10 yields the time histories shown in Figure
4.5. Comparing Figure 4.4 and 4.5, we notice that the time to reach an effective
steady state has reduced dramatically, somewhat paradoxically.

(a) (b)
Figure 4.5: (a) Time history of the dynamical system shown in Figure 4.4b with 002 . 0 = a . (b)
Time history of the dynamical system shown in Figure 4.4b with 02 . 0 = a .
88

Lastly, we point out that when

0 = = = f e b , the conditions in Result 4.1 of Section
4.3.1 are no longer satisfied and therefore we can no longer guarantee that there are
no limit cycles in the first quadrant of the phase space (see relation (4.8)). We then
have the simplified dynamical system given by

x c axy
dt
dx
) 1 ( − − = (4.20)

and

gy axy
dt
dy
+ − = , (4.21)

which has the invariant

) 1 (
0
) ( c g y x a
y x c e
− +
= (4.22)

where
0
c is a constant, determined from the initial conditions, ) 0 ( x and ) 0 ( y .
Relation (4.22) results in closed orbits for , 0 > g 1 < c . The fixed point given by
(4.9) is now no longer hyperbolic, and because of the closed orbits around it, it is a
nonlinear center. Figure 4.6 shows the periodic orbits for this situation, the periods
being functions of
0
c .
89

Figure 4.6: Phase portrait showing the vector field for 1 , 0 < = = = c f e b , and the periodic
orbits around the nonlinear center at 100
0
= x , 1000
0
= y , as given by equation (4.9). The orbits
are described by the closed form relation (4.22).

However, when 1 > c , the orbits are no longer closed as seen from the vector field,
since the x-component (see equation (4.4)) of the phase particle’s velocity is now
always positive. The y-component of the field is negative for a g x / > and positive
for a g x / < pointing to the fact that 0 → y , as ∞ → t . We note from equation
(4.21) that y may not vary monotonically with time if a g x / ) 0 ( < . As the limiting
case of the situation shown in Figure 4.2a, when , 0 → f . 0 ) ( = ∞ → t y

90

Figure 4.7: The x-axis becomes a line of a fixed points when c = 1, and all trajectories end up on the
x-axis. The black solid line shows the nullcline at 200 /
0
= = a g x . Initial conditions that fall to
the left of the nullcline show an increase in the S population before it eventually fades away to zero.

The phase portrait undergoes a dramatic change when c moves from less than unity
to greater than unity, and we have a bifurcation when c = 1. Thus when the rate of
death/self-destruction of Terrorists equals the rate at which new Terrorists are
imported into the area from its neighborhood, the behavior of the differential
equations shows that the fixed point
0
y moves from 1/a when c = 0, to
0
y = 0 when
c =1. When c = 1, this fixed point thus moves to the x-axis. Also, the fixed points
now lie along the x-axis, 0 ≥ x . The nullcline is given by a g x / = , and for initial
values of the Terrorist population, a g x / ) 0 ( < , the population of Susceptibles (S)
initially increases before going to zero (see Figure 4.7). The phase trajectories all end
on the x-axis, the final (steady state) population of Terrorists, ,
f
x being given by the
implicit relation (see equation (4.22))

91
) 0 ( )
) 0 (
ln( )) 0 ( ( ay
x
x
g x x a
f
f
+ = − . (4.23)

To appreciate the difference caused by the parameter c in the dynamics, we show the
time trajectories below for the two regimes of behavior when 1 0 < ≤ c , when 1 > c ,
as well as at the bifurcation point 1 = c (see Figure 4.8).

Figure 4.8a: Stable limit cycle behavior when c = 0.2 < 1. Note the soliton-like behavior in the
dynamics of the Terrorist population (T) when the initial population of Susceptibles (S) far exceeds
that of T. The parameter values (except for the value of c) are the same as those shown on the phase
plot in Figure 4.7.

Figure 4.8b: Time histories of the dynamics when c = 1 starting from different initial conditions. The
parameters (except for the value of c) are the same as those for Figure 4.8a. The initial conditions
(same as in Figure 4.8a) can be more clearly seen here.
92

Figure 4.8c: Unstable behavior showing an explosion in the Terrorist population for c =1.2 >1; the
population of Susceptibles declines to zero eventually.

Having now understood the dynamics of the baseline situation we are ready to study
the effect caused by nonzero values of the parameters b and e, which represent
intervention by different means.

4.3.5 Effects of Nonviolent Intervention (0 , 0 > = e b )

We begin by considering the locations of the fixed points when 0 = b and 0 > e , and
comparing it with the baseline situation. We note that the existence of fixed points in
the baseline situation require (see equation (4.9)) that 1 < c , and 1 < + f c . Relations
(4.9) and (4.10) point out that the presence of nonviolent intervention will have no
effect on the steady state value of the S population, which will remain a c / ) 1 ( − .
Denoting the steady state Terrorist populations
0
0
0 0
~
>
= =
e
b x x and
0
0
0 0
=
= =
e
b x x we can
show, after some algebra, that

93

) 1 ( )
~
)( 1 (
) 1 (
~
0 0
0
0
0 0
f c a x x c e
x c e
x
x x
− − + + −
−
=
−
. (4.24)

From equation (4.24) we see that if 0 > g , 1 < c , and 1 < + f c , we find that
0 0
~
x x > . Thus nonviolent intervention causes the steady state value of the Terrorist
population to decrease when compared to the situation with no intervention of any
kind.

Furthermore, unlike what happens when , 0 = e in the presence of nonviolent
intervention, when 1 < c and 1 > + f c , the dynamical behavior generates an
attracting fixed point. Figure 4.9 when contrasted with Figure 4.3 shows that the
presence of intervention causes the Terrorist population (and the S population) to be
bounded.

(a) (b)
Figure 4.9: The presence of nonviolent intervention causes the dynamical system to exhibit a fixed
point. Unlike the explosion in the Terrorist population that occurs when 0 = = e b , and 1 > + f c ,
with 1 < c , the Terrorist population in the presence of nonviolent intervention is bounded and the
94
dynamical trajectories here are spirals in phase space. The fixed point is at 350 , 16
0 0
= ≈ y x . The
parameters are the same as those used in Figure 4.3 except that 0 > e .

The reason for this is that the asymptote to the nullcline 0
2
= − − + gy fx y ex axy as
∞ → x is now the x-axis, and hence the horizontal nullcline a c y / ) 1 ( − = will
intersect it as long as 1 < c , irrespective of the value of c f + , and a stable fixed
point results in the first quadrant 0 , > y x . When 1 = c , the horizontal nullcline lies
along the x-axis, and the two nullclines intersect at ∞ → x , and the system is
unstable.

For 1 < c , and 1 < + f c , we see a drastically lower steady state population of
Terrorists when nonviolent intervention is provided, as shown in Figure 4.10. For
comparison between the situation with and without nonviolent intervention we show
the dynamics using the same parameters as in Figure 4.4a except that we now have
0 > e .

(a) (b)
95
Figure 4.10: Using the same parameters as for Figure 4.4a except that 0001 . 0 = e , we see a
dramatic drop in the equilibrium population of Terrorists through nonviolent intervention (compare
with Figure 4.4a). The fixed point is now at 4500 , 14
0 0
= ≈ y x .

Lastly, we observe from relation (4.10) that for large values of the parameter e , and
1 < c , e x / 1
0
∝ , and so the T population can be, theoretically speaking, driven
down to zero; note, however, that the S population is unaffected by the presence of
nonviolent intervention.

When 1 > c , the effect of nonviolent intervention cannot stop the explosion in the
Terrorist population (see Figure 4.11), though it asymptotically brings the population
of Susceptibles under control. The reason is the same as for the situation with no
intervention—the x-component of the velocity of the phase flow as seen in equation
(4.4) is always positive in the positive quadrant and hence there can be no fixed point
with 0
0
> x .

(a) (b)
Figure 4.11: Phase plot showing the ineffectiveness of nonviolent intervention in curbing the
Terrorist population. Parameter values are same as those for Figure 4.2 except that e = 0.0001. The
vector field has a positive x-component of velocity and hence the Terrorist population continually
increases.
96

4.3.6 Effects of Military/Police Intervention (b > 0, 0 = e )

The equations that determine the fixed points of the dynamical system are

a
bx c
y
0
0
1 + −
= , (4.25)

0 ) 1 ( ] ) 1 ( [ ) (
2
= − − − − − + = c g x gb f c a abx x r . (4.26)

The roots of equation (4.26) are given by

ab
B A
x ,
2
0
±
= (4.27)

where

A = af c a gb + − − ) 1 ( , B = abfg c f a gb 4 )] 1 ( [
2
+ − − + . (4.28)

We begin by noting that for 0 , , > g f b , 1 < c , 1 < + c f , we have, after some
algebra,
0
0
0
0
0
0
=
>
=
= >
e
b
e
b x x . We next consider 2 cases.

A. When 0 , > g f
97
We begin by noting that in the presence of military intervention ( 0 > b ) a fixed point
exists when 0 , > g f . This is because B is always positive. Three regimes of
behavior now surface, as far as the fixed points are concerned. (i) When c < 1 and A
< 0, then A B > , and we get one positive root 0
0
> x . (ii) c < 1 and A ≥ 0, we
again get one positive root 0
0
> x . (iii) When c > 1, we find that A > 0 and
A B < < 0 , and we get two real positive roots of the equation 0 ) ( = x r with the
two roots 0
) 1 (
0
) 2 (
0
> > x x . When c < 1, and 0
0
> x , we always have, by relation
(4.25) 0
0
> y . However, when c > 1, for 0
0
> y , we require by (4.25) that
b
c
x
1
0
−
> .
Furthermore, when
b
c
x
1
0
−
= , we find that 0 ) (
0
< x r (see Appendix C, part a), so
only the larger of the two roots,
) 2 (
0
x , of the equation 0 ) ( = x r will yield a
corresponding value of 0
0
> y . Thus military intervention causes the dynamics to
always have a fixed point in the first quadrant, which is stable.

We observe from equation (4.11) that when c < 1, the steady state value of the S
population is always larger than when nonviolent intervention is solely used and also
when no intervention of any kind is used. For values of b that are large, the fixed
point a g x /
0
∝ , and as b increases the steady state value of the S population is
larger by ) /(
2
ga b when compared to that obtained in the presence of no military
intervention (or with nonviolent intervention).
98

(a) (b)
Figure 4.12: A stable spiral is seen when 1 < c , 0 < A . Note the long time that the trajectory stays
in the vicinity of the y-axis. Same parameters used as in Figure 4.4b (except for value of b). The fixed
point occurs at . 1573 , 146
0 0
= = y x

(a) (b)
Figure 4.13: A stable node is seen when 1 < c , 0 > A . Compare with Figure 4.4b. The fixed point
occurs at . 3000 , 120
0 0
= = y x

Figures 4.12 and 4.13 show the dynamical behavior when 1 < c . When 0 < A , we
see a stable spiral. As the value of b increases, A becomes positive and the dynamical
behavior switches to a stable node. For comparison, we have used the same
parameter values we used in Figure 4.4b, except for the value of b.

99
The influence of the parameter c can be substantial in determining the steady state
values of the T and S populations. For example, with a = b = 0.0001, f = 0.15, and g
= 0.02, the steady state T populations for c = 0.7 and for c = 1.2 are about 361 and
3588, respectively; the steady state S populations are about 3361 and 1588,
respectively. Figure 4.14 shows the phase portraits for the two cases. Both cases lead
to stable spirals. However, the time taken to reach these steady states also differs
substantially, the system with the larger c value reaching a steady state much more
rapidly (see Figures 4.15 and 4.16 that show the time histories). We have used the
same parameters in Figure 4.14(b) as were used in Figure 4.11 except for the values
of b and e. This comparison illustrates that the presence of military action will cause
a stable fixed point to appear even when 1 > c , something that does not happen when
only nonviolent intervention is undertaken.

(a) (b)
Figure 4.14: The effect of the parameter c in its influence on the steady state values and the
dynamical behavior of the system. (a) c = 0.7. The fixed point is at 3361 , 361
0 0
≈ ≈ y x ; (b) c =
1.2. The fixed point is at 1588 , 3588
0 0
≈ ≈ y x .
100

(a) (b)
Figure 4.15: Time histories of response of the dynamical system with c = 0.7, and the parameters
shown in Figure 4.14(a). (a) Terrorist Population versus Time (b) Susceptible Populations versus
Time

(a) (b)
Figure 4.16: Time histories of response of the dynamical system with c = 1.2, and the parameters
shown in Figure 4.14(b). (a) Terrorist Population versus Time (b) Susceptible Populations versus
Time. Note the much larger steady state Terrorist population as compared to Figure 4.15.

B. When 0 , 0 = > f g
The situation when 0 = f is interesting in that we can have two fixed points in the
first quadrant. We find that

101
2
) 1 (
0
) 1 (
0
1
,
a
bg
a
c
y
a
g
x +
−
= = , and, 0 ,
1
) 2 (
0
) 2 (
0
=
−
= y
b
c
x , when, c
a
gb
> + 1
(4.29)

and

0
~
,
1
~
) 1 (
0
) 1 (
0
=
−
= y
b
c
x , when c
a
gb
< +1 . (4.30)

It can be shown (see Appendix C, part b) that the fixed point ) , (
) 1 (
0
) 1 (
0
y x is either a
stable node or a spiral. For 0 ,
~
) 2 (
0
) 1 (
0
> x x , we require 1 > c . Thus when 1 > c and
c
a
gb
> +1 , we have two fixed points, as shown in relations (4.29). The fixed point
) , (
) 2 (
0
) 2 (
0
y x is, however, unstable. The fixed point )
~
,
~
(
) 1 (
0
) 1 (
0
y x which occurs when
c
a
gb
< + < 1 1 is a stable node, as shown in Appendix C, part b.

We next address the question of comparing the steady state population of Terrorists
when we only have nonviolent intervention as compared to when we only have
military intervention, keeping the parameters c, f and g the same. We note that when
0 = b , 0 > e , we have a stable fixed point 0 ,
0 0
> y x as long as 1 < c . On the other
hand when 0 > b , 0 = e , we always have a fixed point in the first quadrant with
0
~
,
~
0 0
> y x for all values of 0 > c . In order to compare the steady state values of the
102
T and S populations in the two situations, we therefore restrict our attention to the
situation when 1 < c . We then have the following result.

Result 4.3. When 1 < c , and 0 ) 1 ( > − − ab c e , then
0 0
~
x x > , when the parameters c, f
and g are kept unchanged. This means that under these circumstances, the steady
state population of Terrorists is larger when only military/police intervention is
utilized, as compared to when only nonviolent intervention is used.
Proof: From equation (4.13) we see that 0
0
> x satisfies the relation

) 1 ( ) 1 ( ) 1 (
0
2
0
c g x f c a x c e − = − − + − (4.31)

and 0
~
0
> x satisfies the relation

) 1 (
~
] ) 1 ( [
~
0
2
0
c g x gb f c a x ab − = − − − + (4.32)

Substituting for ) 1 ( c g − from equation (4.31) in equation (4.32) gives

0
~
)
~
)( 1 ( )] 1 ( [ )
~
)(
~
(
0 0 0
2
0 0 0 0 0
= − − − − + − − + + − x gb x x f c a x c e ab x x x x ab ,
(4.33)

where we have added and subtracted the quantity
2
0
abx . This results in

103

) 1 ( )
~
(
] ) 1 ( [
~
)
~
(
0 0
2
0 0
0 0
f c a x x ab
x ab c e x gb
x x
− − + +
− − +
= − . (4.34)

Since from equation (4.31) we have
0 0
2
0
~
) 1 (
~
)] 1 ( [
~
x gb c g x f c a x ab + − = − − + ,
equation (4.34) becomes upon multiplying the numerator and denominator on both
sides by
0
~
x ,

0 0 0
2
0 0
2
0
0 0
~
) 1 ( )
~
(
~
] ) 1 ( [
~
)
~
(
x gb c g x x ab
x x ab c e x gb
x x
+ − +
− − +
= − (4.35)

The denominator is always positive, and to ensure that
0 0
~
x x > , the result now
follows. From the proof we note that the condition stated in Result 4.3 when 1 < c is
a sufficient condition for
0 0
~
x x > , but not a necessary condition.

Figure 4.17 illustrates Result 4.3, where we compare the effect of military
intervention with that produced by nonviolent intervention. Using the same
parameter values as in Figure 4.14(a), and here 0 ) 1 ( > − − ab c e . We find that the
fixed point has moved from 3361
~
, 361
~
0 0
≈ ≈ y x when using military intervention to
3000 , 14
0 0
= ≈ y x when using nonviolent intervention, though the time needed to
achieve the steady state value is much greater in the latter case. While the change in
104
the population of Susceptibles is not appreciable, that in the Terrorist population is
substantial, as seen from the two figures.

(a)

(b) (c)
Figure 4.17: A comparison with Figures 4.14(a) and 4.15 showing that nonviolent intervention brings
about lower steady state Terrorist populations than those brought about by military/police intervention,
when 0 ) 1 ( > − − ab c e . The fixed point occurs at 3000 , 14
0 0
= ≈ y x .

4.3.7 Effects of Combined Military/Police and Nonviolent Intervention (e, b > 0)

105
When both military/police intervention and nonviolent intervention are present, the
fixed point
0
x satisfies relation (4.13). When 1 < c , as pointed out in section 4.3.2,
we have a fixed point ) , (
0 0
y x , where
0
x is positive real root of

0 ) 1 ( ] ) 1 ( [ ] ) 1 ( [ ) (
2 3
= − − − − − + + − + = c g x gb f c a x ab c e ebx x r . (4.36)

This fixed point could be a stable spiral or a stable node depending on whether
relation (4.17) or (4.18) is satisfied.

However, when 1 > c , the coefficient of the
2
x term in (4.36) is sign indefinite while
the coefficient of the x term is negative. Also,

b c af
b
c
r / ) 1 ( )
1
( − − =
−
(4.37)

is negative when 0 > f and 1 > c . Hence Descarte’s rule points out that we would
have 2 positive roots,
) 1 (
0
) 2 (
0
x x > , of the equation 0 ) ( = x r , and one negative root.
Since we require
b
c
x
1
0
−
> , only the root
) 2 (
0
x will then yield a corresponding
0
y
that is positive. Hence we will have only one fixed point ) , (
0 0
y x for which
0 ,
0 0
> y x . As pointed out in Section 4.3.3, the fixed point will again be either a
stable spiral or a stable node, depending on whether relation (4.17) or (4.18) is
satisfied.
106

(a)

(b) (c)

Figure 4.18: The effect of both military and nonviolent action are illustrated here for 1 < c . The
fixed point is at 3013 , 14
0 0
≈ ≈ y x , and the trajectories appear to reach the equilibrium point
marginally faster than with only nonviolent intervention. All the parameters are the same for Figures
4.14(a), 4.15, 4.17, and 4.18, except for the values of b and e. A soliton-like behavior is observed
from the phase plot.

Figure 4.18 shows the effect of combined military/police and nonviolent intervention
when 1 < c . For 1 > c , we see that as in the case of only military intervention, we
107
have a stable node, that is reached relatively rapidly when compared with 1 < c , as
seen in Figure 4.19.

Figure 4.19: The parameters used here are the same as those for Figure 4.18, except that 2 . 1 = c .
The fixed point is a stable node with 1 , 2000
0 0
≈ ≈ y x .

A comparison of Figures 4.14(b) (and Figure 4.16) with Figure 4.19 indicates the
difference that nonviolent intervention makes when added to military/police
intervention. We note that the steady state value of the Terrorist population is
reduced substantially, as is the population of Susceptibles, while the time required to
reach steady state is increased considerably.

From our previous understanding of the various cases, we can heuristically say that
the presence of military intervention causes the Terrorist population to be
asymptotically bounded, something that cannot occur with only nonviolent
intervention; it also reduced the steady state value of the Susceptible population. The
addition of nonviolent intervention reduces the steady state value of the Terrorist
108
population as compared with what it might have been only with military/police
intervention, but is also increases the time required to reach the steady state.

Result 4.4: The addition of military/police intervention to nonviolent intervention
always reduces the steady state population of Terrorists from what it would be with
only nonviolent intervention. Hence,

0
0
0
0
0
0
=
>
>
> <
b
e
b
e x x . (4.38)

Proof: We note that a steady state is reached with only nonviolent intervention when
1 < c . Let 0
0
> x be the steady state Terrorist population when only nonviolent
intervention is used, and let 0
~
0
> x be the steady state Terrorist population when
both military/police and nonviolent interventions are employed. Thus we have, with
1 < c ,

0 ) 1 ( ) 1 ( ) 1 (
0
2
0
= − − − − + − c g x f c a x c e (4.39)

and

0 ) 1 (
~
] ) 1 ( [
~
] ) 1 ( [
~
0
2
0
3
0
= − − − − − + + − + c g x gb f c a x ab c e x eb . (4.40)

Substituting for ) 1 ( c g − from equation (4.39) into equation (4.40) we get
109

) 1 (
~
] ) 1 ( [ )
~
(
~
) 1 /(
~
~
0 0
2
0 0 0 0
0
2
0
0 0
c g x x ab c e x x x x eb
c x abfx
x x
− + + − + +
−
− = − (4.41)

from which the result follows.

Result 4.5: The addition of nonviolent intervention to military/police intervention
always reduces the steady state population of Terrorists from what it would be with
only nonviolent intervention. Hence, with , 0 , > g f we have

0
0
0
0
0
0
>
=
>
> <
b
e
b
e x x . (4.42)

Proof: We note that a steady state is reached with only military/police intervention
for any 0 > c . Let 0
0
> x be the steady state Terrorist population when only
military/police intervention is used, and let 0
~
0
> x be the steady state Terrorist
population when both military/police and nonviolent interventions are employed.
Thus we have,

0 ) 1 ( ] ) 1 ( [ ) (
0
2
0 0
= − + − − − + ≡ c g x gb f c a abx x r (4.43)

and

110
0 ) 1 (
~
] ) 1 ( [
~
] ) 1 ( [
~
0
2
0
3
0
= − + − − − + + − + c g x gb f c a x ab c e x eb . (4.44)

Substituting for ) 1 ( c g − from equation (4.44) into equation (4.43) we get

) 1 (
~
)] 1 (
~
[
~
~
0 0
0
2
0 0
0 0
c g x abx
c x b x ex
x x
− +
− +
= − , (4.45)

which is clearly positive when 1 0 ≤ ≤ c , since both
0
x and
0
~
x are positive.

As pointed out in Section 4.3.6, when 1 > c and 0 , > g f , equation (4.43) has two
positive roots of which the larger yields 0
~
0
> y . And for this to happen, from
relation (4.12) we see that
b
c
x
1
~
0
−
> , so that the numerator in relation (4.45) is
positive. To prove that the denominator in (4.45) is positive when 1 > c we note that
since 0 ) ( < − = gf
a
g
r , a g x /
0
> , and hence

) 1 (
1
~
0 0
− =
⎟
⎠
⎞
⎜
⎝
⎛ −
⎟
⎠
⎞
⎜
⎝
⎛
> c g
b
c
a
g
ab x abx . (4.46)

Thus we have shown that
0 0
~
x x − > 0 when . 0 ≥ c

111
We have thus shown that a combination of military/police action and nonviolent
intervention yields a fixed point with a lower Terrorist population that with any one
of these interventions excluded.

4.4. Conclusions

In this paper we have developed a preliminary dynamical systems approach to
understanding the dynamics of terrorism. We have conceptualized the population in
a given region in terms of Terrorists, Susceptibles, and Non-susceptibles. We
consider Terrorists to be individuals who are firmly committed to terrorism, who will
not yield to nonviolent actions such as persuasion through economic, political, social,
or other means. They are what might be called, for example, the ‘insurgent’
population in areas like Mosul and Bagdad. The segment we term as Susceptibles
constitute those that can be influenced, either into becoming Terrorists through
terrorist propaganda or into becoming pacifists that constitute what we call the Non-
Susceptible population. The population segment of Susceptibles would, for example,
include individuals who have been educated in madrasaas, or those who might be
influenced from a young age by Wahabi thinking. We assume that military/police
intervention is used for reducing the rate of expansion of the Terrorist population,
while nonviolent action is used to target Susceptibles into becoming Non-terrorists.
The flow between these three populations and the effect of military/police action as
well as nonviolent action is illustrated in Figure 4.1. We assume that military/police
action increases as the square of the Terrorist population, and that nonviolent action
112
is proportional to the product of the size of the Susceptible population and the square
of the Terrorist population.

The model we have used here is normative, and the usefulness of the results adduced
herein would naturally depend on the extent to which our model approximates reality.
In this respect, our model is somewhat different from other socio-politico-economic
models, such as those that may use game theory or optimal control where matters
like game set-ups and problem frameworks, player payoffs, and cost (utility)
functions are difficult to assess, along with the basic assumptions underlying these
approaches. The determination of the model parameters/functions in these situations
is, at best, also difficult. The model we have proposed contains parameters, which
could be estimated with greater ease, and from data that may be more widely
available. We intend to take the estimation of these parameters from field data in a
follow-on study so that further improvements in the description of the dynamical
characterization of the system can be incorporated.

The differential equations that model the system have 8 (effectively 6) parameters
and it is somewhat remarkable that, though highly nonlinear, the dynamical regimes
of behavior can still be analyzed, in general, in a relatively straightforward manner.
Our study here has basically concentrated on the fixed points of the dynamical
system, their existence, and their stability. Several issues such as the transient
behavior of the system, which may be important from a policy and strategic planning
113
viewpoint, though touched upon, have not been emphasized, nor has the possibility
of the parameters of the model being functions of time been entertained.

The dynamical behavior of our model indicates several general features, many of
which point to new strategic ways of planning to quell terrorism. While some of
these features may appear reasonable from an intuitive stance, there are many that do
not immediately suggest themselves purely on intuitive grounds. More importantly
the study points to the various conditions under which different strategies might be
meaningful and effective, and their possible effects. We summarize some these
observations below.

1. There is no possibility of the populations of the Terrorists and susceptibles
evolving chaotically in a manner that precludes prediction because of extreme
sensitivity to the small uncertainties in determining the population numbers of the
two populations at any given time. From a mathematical standpoint, dynamical chaos
is precluded in the model that has been proposed herein to describe the behavior of
the system.

2. In the presence of military/police intervention against Terrorists with or without
nonviolent intervention, the dynamical system always reaches a unique equilibrium
state as long as the Susceptible population keeps growing; that is, if 0 , > g a , and
0 ≥ e , a unique (nonzero) equilibrium population of Terrorists and susceptibles will
emerge as long as 0 > b and/or 0 > f .
114
3. As long as 0 ) 1 ( > − − ab c e , the equilibrium population of Terrorists is less when
solely nonviolent action is used instead of military/police intervention. This result is
somewhat unexpected and points out new directions in strategic planning for curbing
terrorism.

4. Were no intervention to be made so that 0 = = e b , the system becomes unstable if
1 > c , or if 1 > + c f . The first of these conditions, 1 > c , appears intuitive for it
simply states that Terrorists are increasing in number faster than they are dying/self-
destructing. The latter condition, however, is less intuitive and points out that it is
not necessary for 1 > c for the terrorist population to explode; for, the increase in
Susceptibles caused by terrorist propaganda can also lead to such an explosion, the
Susceptibles being seduced into terrorist behavior through their recruitment by
Terrorists, as per the axy term in equation (4.4). On the other hand, were either
1 > + c f and/or 1 > c , the terrorist situation would go out of hand unless some sort
of interdiction is made.

5. Perhaps the most critical parameter which emerges from this dynamical
investigation is the value of c, the ratio of the rate of increase of Terrorists either
through recruitment from among their own population or importation of Terrorists
from other geographical areas to the rate at which terrorist die/self-destruct:
(a) If this ratio is less than unity, nonviolent action will cause the populations
eventually to reach a steady state. The time to reach steady state is large and the
population of Susceptibles may fluctuate considerably, however the orbits are seen to
115
lie in a narrow region around the y-axis, pointing to the possible advantages of using
nonviolent action in curbing Terrorist populations. However, if this ratio exceeds
unity ( 1 > c ), nonviolent action will not be able to stop the upward increase in the
Terrorist population. Thus the presence of porous boundaries in an area of terrorist
activity can lead to substantially different dynamical behavior, and our model points
to the fact that sealing the boundaries of a region to prevent the flow of Terrorists
into it should be a matter of considerable strategic importance in curbing terrorist
activities.
(b) Military/police response always results in the two populations reaching an
eventual steady state. This occurs even if 1 > c .

6. The combined action of both military/police and nonviolent initiatives, in general,
brings the equilibrium population of Terrorists considerable lower than with only
military intervention, though it appears to increase the duration of time necessary to
reach the steady state.

7. The variation of the populations of Terrorists and susceptibles may show a soliton
like behavior in time, even though these populations reach steady states. This is
because the dynamical flow is slow in the vicinity of the y-axis. This means that the
Terrorist population can decrease, and can remain decreased, for extended periods of
time, before possibly increasing substantially again. Quiescent periods with low
terrorist activity, because of low Terrorist populations, seem to be an inherent
feature of the dynamics of the system. They are unrelated to explanations like the
116
increased effectiveness of military or nonviolent actions. This points to the fact that
curbing terrorism has to be a long-term strategy which must continually proceed
despite periods of time when terrorist activities may be seen to be of negligible
occurrence. This is because the dynamics itself has this nature, that is, long
‘quiescent’ periods of time could elapse between ‘explosions’ of Terrorist
populations.

8. The equilibrium populations, when they exist, are shown to be stable. Furthermore,
from a dynamical standpoint, these equilibrium points are shown to be either stable
spirals or nodes. This would mean that even in the presence of military and/or
nonviolent intervention, the population of Terrorists, together with their activities,
can increase before the equilibrium population is eventually reached. This is
controlled by the initial population numbers of the Terrorists and the Susceptibles in
relation to the location (and stability) of the equilibrium populations. Thus, the fact
that terrorist activities (especially over the short haul) may increase in spite of
military/police and/or nonviolent intervention does not indicate that such
intervention is ineffective.

9. In addition to Results 4.3-4.5 that are related to the steady state population of
Terrorists when it exists, the steady state population of susceptibles,
0
y , when it
exists, has the following property (see Appendix C, part c).

117
0
0
0
0
0
0
0
0
0
0
0
0
=
>
>
>
>
=
=
= < < =
e
b
e
b
e
b
e
b y y y y .

10. Military/police intervention appears to be necessary for allowing Terrorist
populations to reach stable values in areas where the boundaries are porous and 1 > c .
The addition of nonviolent intervention appears to reduce the steady state Terrorist
population substantially. In areas where the relative increase in the Terrorist
population is such that 1 < c and 1 < + c f , the use of nonviolent action might be
preferable to military/police intervention if the intention is to bring down the
eventual value of the Terrorist population.
118
APPENDICES

APPENDIX A

Property 1:
+ + +
= − U B B I U ) (.
Proof. Since ) ( ) )( (
1
B B I M B B I U
T + − +
− − = and
T T
U U U U
+ +
= ) ( , we have

) ( ) ( ) ( B B I U U U B B I U
T T + + + +
− = −
) )( )( )( ( ) (
1
B B I B B I M B B I U U
T + + − + +
− − − =
T T T
U U U B B I M B B I U U
+ + − + +
= − − = ) ( ) )( )( ( ) (
1

.
+
= U □

Property 2:
+ + +
= − U U B B I ) (.
Proof. Since ) ( ) )( (
1
B B I M B B I U
T + − +
− − = and
+ +
= ) (
T T
UU U U , we have

+ + + +
− = − ) ( ) ( ) (
T T
UU U B B I U B B I
+ + − + +
− − − = ) )( )( )( )( (
1 T
UU B B I M B B I B B I
+ + + − +
= − − = ) ( ) )( )( )( (
1 T T T
UU U UU B B I M B B I
.
+
= U □

Property 3: 0 = BU
119
Proof. Since ) ( ) )( (
1
B B I M B B I U
T + − +
− − = and B B BB =
+
, we obtain
) ( ) )( (
1
B B I M B B I B BU
T + − +
− − =
) ( ) )( ( ) ( ) )( (
1 1
B B I M B B B B I M B BB B
T T + − + − +
− − = − − =
0 = . □

Property 4: 0 =
+
UB
Proof. Since ) ( ) )( (
1
B B I M B B I U
T + − +
− − = and
+ + +
= B BB B , we have
+ + − + +
− − = B B B I M B B I UB
T
) ( ) )( (
1

) ( ) )( ( ) ( ) )( (
1 1 + + − + + + + − +
− − = − − = B B M B B I BB B B M B B I
T T

0 = . □

Property 5: 0 =
+
BU
Proof. Since ) ( ) )( (
1
B B I M B B I U
T + − +
− − = ,
+ +
= ) (
T T
UU U U , and B B BB =
+
, we
get

+ + − + +
− − = ) )( )( )( (
1 T
UU B B I M B B I B BU
+ + − +
− − = ) )( )( )( (
1 T
UU B B I M B BB B
+ + −
− − = ) )( ( ) )( (
1 T T
UU B B I M B B
0 = . □

Property 6: 0 =
+ +
B U
120
Proof. Since ) ( ) )( (
1
B B I M B B I U
T + − +
− − = ,
T T
U U U U
+ +
= ) (, and
+ + +
= B BB B , it
yields

+ + − + + + +
− − = B B B I M B B I U U B U
T
) )( )( ( ) (
1

) )( )( ( ) (
1 + + + − + +
− − = BB B B M B B I U U
T

) )( )( ( ) (
1 + + − + +
− − = B B M B B I U U
T

0 = . □

Property 7: When I CC =
+
,
+ +
= BB BC BC ) )((.
Proof. From the third property of Moore-Penrose generalized inverse, we have

[] [ ] [ ]
T T
T
T
T T
B C BC BC BC BC BC BC BC
+ + + +
= = = ) ( ) ( ) ( ) )( ( ) )((.

Also
T
B can be written as
+ + + +
= = = = BB B BB B B BB B B B B
T T T T T T T T
) ( ) ( ) ( ) )((.
Using this relationship, and since I CC =
+
, the above condition can be described as

[]
T T
T
B C BC BC BC
+ +
= ) ( ) )((
[ ] [ ]
+ + + + +
= = B BCC B C BC BB B C BC
T T
T
T T
T
) ( ) (
[ ] [ ]
+ + + + + +
= = B C BC BC BC B C BC BC BC
T
T
T
) ( ) )( ( ) ( ) ( )(
+ + + + + +
= = = BB B C BC B C BC BC BC ) ( ) ( ) )( (. □

121
Property 8: When M is symmetric and positive definite matrix we can define
+ +
= ) ( ) )( (
2 / 1 2 / 1 T T
BMB BMB BM BM .
Proof. Since for any m by n matrix A , whatever its rank,
+ +
= ) (
T T
AA A A , we have

+ + +
= = ) ( ) )( ( ) )( (
2 / 1 2 / 1 2 / 1 2 / 1 2 / 1 2 / 1 T T T T
BMB BMB B M BM B M BM BM BM . □

Property 9: Any m by n matrix B its transpose,
T
B , can be written as
T T T
BB B BB B B
+ +
= = .
Proof. Because
T T
B BB B ) (
+
= , we have

T T T T T
BB B B B B B BB B
+ + +
= = = ) ( )(

and

+ + +
= = = BB B BB B B BB B
T T T T T
) ( )(. □

Property 10: When M is symmetric and positive definite matrix, the relationship
+ +
BMB BMB B
T T
) ( is simply
+
B .
Proof. Using property 8, the relationship can be written as

T T T T T T T T T
B BMB BMB B B BMB BMB B BMB BMB B ] ) ( ) [( ] ) ( [ ) (
+ + + + + + +
= =
122
+ + + + + + + +
= = = = B BB B B B B B BB B B
T T T T
] ) [( ] ) [( . □

Property 11: When we define a matrix BM BMB MB M E
T T +
− = ) ( , we have
E B B I E = −
+
) (.
Proof. When we expand the relation ) ( B B I E
+
− using property 10, we get

) ]( ) ( [ ) ( B B BM BMB MB M E B EB E B B I E
T T + + + +
− − = − = −
) ( ) ( ) ( B B BM BMB MB B B M E
T T + + +
− − =
E B B M B B M E = − − =
+ +
) ( ) (. □

Property 12: When M is symmetric and positive definite matrix, B B I EU
+
− = .
Proof. Using the property 11, we get

) )( ( ) )( )( (
1 1
B B I M E B B I M B B I E EU
+ − + − +
− = − − =
) )( ]( ) ( [
1
B B I M BM BMB MB M
T T + − +
− − =
) ]( ) ( [ B B I B BMB MB I
T T + +
− − =
) ( ) ( ) ( B B B BMB MB B BMB MB B B I
T T T T + + + +
+ − − =
B B I
+
− = . □

Property 13: E E B B I = −
+
) (.
123
Proof. When we expand the relation E B B I ) (
+
− using property 8, we get

] ) ( )[ ( ) ( BM BMB MB M B B E BE B E E B B I
T T + + + +
− − = − = −
BM BMB BMB B M B B E
T T + + +
− − = ) )( ( )(
E M B B M B B E = − − =
+ +
) ( ) (. □

Property 14: When M is symmetric and positive definite matrix, B B I UE
+
− = .
Proof. Using the property 13 and 9 we get

E M B B I E B B I M B B I UE ) )( ( ) )( )( (
1 1 − + + − +
− = − − =
] ) ( )[ )( (
1
BM BMB MB M M B B I
T T + − +
− − =
] ) ( )[ ( BM BMB B I B B I
T T + +
− − = .
] ) ( )[ ( BM BMB BB B I B B I
T T + + +
− − =
] ) ( )[ ( BM BMB BB B B B I B B I
T T + + + +
− − − =
B B I
+
− = . □

Property 15: When M is symmetric and positive definite matrix, U UEU = .
Proof. Using property 14 and ) )( )( (
1
B B I M B B I U
+ − +
− − = we have

U U B B I UEU = − =
+
) (. □

124
Property 16: When M is symmetric and positive definite matrix, E EUE = .
Proof. Using property 14 and 13 we have

E E B B I EUE = − =
+
) ( . □

Property 17: When M is symmetric and positive definite matrix, the Moore-
Penrose generalized inverse of the given matrix U is given by
A AM M M U
+ − +
− = ) (
2 / 1 2 / 1
.
Proof. The matrix BM BMB MB M E
T T +
− = ) ( satisfied the following four
properties:

1) U UEU = (see Property 15);
2) E EUE = (see Property 16);
3) UE is symmetric (see Property 14); and
4) EU is symmetric (see Property 12).

Hence, when we substitute
1 −
= AM B into BM BMB MB M E
T T +
− = ) ( , the result
follows. □

APPENDIX B

125
In this section, we provide some properties that are used for verifying the
recursive formulae for determining of the Moore-Penrose M-inverse of a matrix.

Property 1:
T
T
M
c cM
c M
c
1
1
−
−
+
= (see [60]).
Proof: If
T
T
c cM
c M
1
1
−
−
is the MP M-inverse of c, it must satisfy all equations (2.1)-(2.4)
as follows:

i. c c
c cM
c M
c c cc
T
T
M
=
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
=
−
−
+
1
1

ii.
+
−
−
−
−
−
−
−
−
−
−
+ +
= =
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
= =
M
T
T
T
T
T
T
T
T
T
T
M M
c
c cM
c M
c cM
c M
c
c cM
c M
c cM
c M
c
c cM
c M
cc c
1
1
1
1
1
1
1
1
1
1
,
iii.
+
−
−
−
−
+
= =
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
=
M
T
T
T
T
T
T
M
cc
c cM
c M
c
c cM
c M
c cc
1
1
1
1
)(,
iv.
1 1
1
1
1
1
1
1
) (
− + −
−
−
−
−
−
−
+
=
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
= =
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
= cM Mc M c
c cM
c M
M
c cM
cM c
c
c cM
c M
c c
M
T
T
T
T
T
T
T
T
M
.

Thus,
T
T
c cM
c M
1
1
−
−
is the MP M-inverse of c. □

Property 2: 0
ˆ
=
+
M
c A .
126
Proof : Using
T
T
M
c cM
c M
c
1
1
−
−
+
= (see Property 1 in Appendix A), and )
ˆ ˆ
( ˆ A A I a c
M
+
− = ,
we get

T
T T
M
M
c cM
a A A I M
A c A
1
1
)
ˆ ˆ
(
ˆ ˆ
−
+ −
+
−
=
T
T
M
c cM
a M A A M I M
A
1
1 1
)
ˆ ˆ
(
ˆ
−
− + −
−
=
0
)
ˆ ˆ
(
ˆ
1
1
=
−
=
−
− +
T
T
M
c cM
a M A A I A
. □

Property 3: 1 =
+
M
ac .
Proof : Using Property 1 above and )
ˆ ˆ
( A A I a c
M
+
− = , we obtain

T
T T
M
M
c cM
a A A I M
a ac
1
1
)
ˆ ˆ
(
−
+ −
+
−
=
T
M M
T
M
a AM A M I M A A I a
a M A A M I aM
)
ˆ
( )
ˆ ˆ
(
)
ˆ ˆ
(
1 1
1 1
− + − +
− + −
− −
−
=

T
M M
T
M
a M A A I A A I a
a M A A M M aM
1
1 1
)
ˆ ˆ
)(
ˆ ˆ
(
)
ˆ ˆ
(
− + +
− + −
− −
−
= 1
)
ˆ ˆ
(
)
ˆ ˆ
(
1
1
=
−
−
=
− +
− +
T
M
T
M
a M A A I a
a M A A I a
. □

Property 4:
T T
M
u u A A =
+
ˆ ˆ
(when 0 )
ˆ ˆ
( ˆ = − =
+
A A I a c
M
).
Proof : Since
+
=
M
A a u
ˆ
, we have

() ( ) ( )
T
T
M
T
M M
T
M
T
M
u A a A A A a A A u u A A = = = =
+ + + + +
ˆ
ˆ
ˆ ˆ ˆ
ˆ
ˆ ˆ ˆ ˆ
. □

Property 5: A u A A a a
M
ˆ ˆ ˆ
ˆ ˆ = =
+
.
Proof : Since 0 )
ˆ ˆ
( ˆ = − =
+
A A I a c
M
, we get ( ) A u A A a A A a a
M M
ˆ ˆ ˆ
ˆ
ˆ ˆ
ˆ ˆ = = =
+ +
. □
127

Property 6: 0 = Ap .
Proof : Since m M A A I p
M
~
) (
1 −
−
+
−
− = , we obtain 0
~
) (
1
= − =
−
−
+
−
m M A A I A Ap
M
. □

Property 7: a AA a
M
+
−
= (when 0 = d ).
Proof : Because 0 ) ( = − =
+
−
a AA I d
M
, we obtain a AA a
M
+
−
= . □

Property 8: 0 =
+
−
d A
M
.
Proof : Since a AA I d
M
) (
+
−
− = , we obtain 0 ) ( = −
+ +
− −
a AA I A
M M
. □

Property 9: 0 =
+
A d .
Proof : Since
d d
d
d
T
T
=
+
and a AA I d
M
) (
+
−
− = , we have
0
) (
) (
) ( ) (
) (
=
−
−
=
− −
−
= =
+
+
+ +
+
+
−
−
− −
−
a AA I a
A AA I a
a AA I AA I a
A AA I a
d d
A d
A d
M
T
M
T
M
T
M
T
T
M
T
T
T
□

Property 10: 1 =
+
a d .
Proof : Since
d d
d
d
T
T
=
+
and a AA I d
M
) (
+
−
− = , we have 1
) (
) (
=
−
−
=
+
+
+
−
−
a AA I a
a AA I a
a d
M
T
M
T
.
□

128
Property 11: ()
+
−
−
− =
M
T T
A m M v h
~
1
β
.
Proof : Since MU q h
T
β
1
= , where Mq q
T
= β ,
⎥
⎦
⎤
⎢
⎣
⎡
=
+
−
xm
M
A
U
1
0
,
⎥
⎦
⎤
⎢
⎣
⎡
−
+
=
1
p v
q , a A v
M
+
−
= ,
m M A A I p
M
~
) (
1 −
−
+
−
− = , and
⎥
⎦
⎤
⎢
⎣
⎡
=
−
m
m
m
M
M
T
~
~
, we have
⎥
⎦
⎤
⎢
⎣
⎡
⎥
⎦
⎤
⎢
⎣
⎡
=
+
−
−
m
M
T T
T
A
m
m
m
M
Mq q
q
h
0
~
~

Mq q
A m A M p v
A m
A M
Mq q
q
T
M
T
M
T T
M
T
M
T
T
+ +
−
+
+
−
− −
−
−
− +
=
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
=
~
) (
~
.
Since 0 ) (
~
= − =
+ + +
− M M
T
M
T
A A A I m A M p , we obtain ( )
+
−
−
− =
M
T T
A m M v h
~
1
β
. □

Property 12:
−
+
−
=
−
M v A A M v
T
M
T
.
Proof : Since a A v
M
+
−
= and A A M
M
+
−
−
is symmetric (or A A M
M
+
−
−
=
−
+
−
M A A
T
M
) ( ),
we have
− −
+ +
−
+ + +
−
= = =
− − − − −
M v M a AA A M A A a A A A M v
T T
M M
T
M
T
M M
T
) ( ) ( ) ( □

Property 13: () A A m M v hA
M
T T +
−
−
− =
~
1
β
.
Proof : Since ()
+
−
−
− =
M
T T
A m M v h
~
1
β
and
−
+
−
=
−
M v A A M v
T
M
T
(see Property 11 and
12), we have )
~
(
1
)
~
(
1
A A m M v A A m A A M v hA
M
T T
M
T
M
T +
−
+ +
−
− − −
− = − =
β β
, □

Property 14: h hAA
M
=
+
−

129
Proof : Since ()
+
−
−
− =
M
T T
A m M v h
~
1
β
(see Property 11), we have

h A m A M v AA A m AA A M v hAA
M
T
M
T
M M
T
M M
T
M
= − = − =
+ +
−
+ + + +
−
+
− − − − − − −
)
~
(
1
)
~
(
1
β β
.
□

Property 15: )
~
(
1
v m v M v ha
T T
− =
−
β
.
Proof : Because v a A
M
=
+
and ( )
+
−
−
− =
M
T T
A m M v h
~
1
β
(see Property 11), we have

)
~
(
1
)
~
(
1
v m v M v a A m a A M v ha
T T
M
T
M
T
− = − =
−
+ +
−
β β
. □

Property 16: m A A m p M
T
M
~
) (
~ +
−
−
= − .
Proof : Since m M A A I p
M
~
) (
1 −
−
+
−
− = , we get

m m M A A M M m m M A A I M m p M
M M
~ ~
) (
~ ~
) (
~ 1 1
− − = − − = −
−
−
+
− −
−
−
+
− −
− −

( ) m A A m m AM A M m M M
T
M M
~ ~ ~ ~ 1 1 + −
−
+
−
−
− −
− −
− = − − = □

Property 17: β − − = − +
−
m v v M v m p m v m
T T T T ~ ~ ~
.
130
Proof : Since Mq q
T
= β ,
⎥
⎦
⎤
⎢
⎣
⎡
−
+
=
1
p v
q , and
⎥
⎦
⎤
⎢
⎣
⎡
=
−
m
m
m
M
M
T
~
~
, we obtain
⎥
⎦
⎤
⎢
⎣
⎡
−
+
⎥
⎦
⎤
⎢
⎣
⎡
⎥
⎦
⎤
⎢
⎣
⎡
−
+
= =
−
1
~
~
1
p v
m
m
m
M p v
Mq q
T
T
T
β
m p m p M p v m v M p v M v
T T T T T
+ − + − + =
− − −
~
2
~
22.
Above, we have used the fact that p M v v M p
T T
− −
= , m v v m
T T ~ ~
= , and m p p m
T T ~ ~
=
because they are scalars. Since m M A A I p
M
~
) (
1 −
−
+
−
− = and a A v
M
+
−
= , we have
−
−
−
+
−
−
− = M m M A A I M p
T
M
T
]
~
) [(
1
) (
~
] ) ( [
~ 1
A A I m M M A A I M M m
M
T
M
T +
− −
+
−
−
−
− −
− = − = ,
so that p m p M p
T T ~
=
−
and 0 =
−
v M p
T
. Thus, we obtain
m p m v m v M v
T T T
+ − − =
−
~ ~
2 β . Using the fact that m v v m
T T ~ ~
= , we have,

β − − = − +
−
m v v M v m p m v m
T T T T ~ ~ ~

m p m v m m p m v m v M v m v v M v
T T T T T T T
− + = + − − − − =
− −
~ ~
)
~ ~
2 (
~
. □
131
Appendix C

a. When 0 , > g f , and c > 1, for 0
0
> y we require
b
c
x
1
0
−
> . We know that the
equation 0 ) ( = x r given in equation (26) has two real positive roots. Since for
b
c
x
1 −
= , 0
) 1 (
)
1
( <
−
=
−
b
c af
b
c
r , we find that only the larger of the two roots of
0 ) ( = x r will yield a 0
0
> y . Hence, while there are two values of 0
0
> x for which
0 ) (
0
= x r , only the larger of two values of
0
x will result in a 0
0
> y .

b. When 0 > g , 0 = f , the Jacobian evaluated at the fixed point ) , (
) 1 (
0
) 1 (
0
y x is
given by

) , (
) , (
01 01
01 01
) 1 ( 2
y x
y x
ax g ay
ax c bx ay
J
⎥
⎦
⎤
⎢
⎣
⎡
− −
− − −
= =
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
⎟
⎠
⎞
⎜
⎝
⎛
+ −
−
0 1
a
bg
c
g
a
gb
.

Its eigenvalues, λ , are given by

0 1
2
= − ⎟
⎠
⎞
⎜
⎝
⎛
+ + ⎟
⎠
⎞
⎜
⎝
⎛
− − cg
a
bg
g
a
gb
λ λ .

132
Since the trace of the Jacobian is negative and its determinant is positive if
) / ( 1 a bg c + < , the fixed point ) , (
) 1 (
0
) 1 (
0
y x is either a stable node or a stable spiral.

The Jacobian at ) , (
) 2 (
0
) 2 (
0
y x , when 1 > c and 0 1 >
⎟
⎠
⎞
⎜
⎝
⎛
− + c
a
bg
, is given by

) , (
) , (
02 02
02 02
) 1 ( 2
y x
y x
ax g ay
ax c bx ay
J
⎥
⎦
⎤
⎢
⎣
⎡
− −
− − −
=
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
⎟
⎠
⎞
⎜
⎝
⎛
− +
− −
=
c
a
bg
b
a
b c a c
1 0
/ ) 1 ( 1
,

which shows that ) , (
) 2 (
0
) 2 (
0
y x is unstable saddle.

The Jacobian at )
~
,
~
(
) 1 (
0
) 1 (
0
y x is given by

) , (
) , (
0 0
0 0
) 1 ( 2
y x
y x
ax g ay
ax c bx ay
J
⎥
⎦
⎤
⎢
⎣
⎡
− −
− − −
=
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
⎟
⎠
⎞
⎜
⎝
⎛
− +
−
−
=
c
a
bg
b
a
b
c a
c
1 0
) 1 (
1
,

whose eigenvalues are both negative since c
a
gb
< + < 1 1 . Hence the fixed point is a
stable node.

133
c. Result: The steady state population of Terrorists when it exists, the steady state
population of Susceptibles,
0
y , when it exists, has the following property.

0
0
0
0
0
0
0
0
0
0
0
0
=
>
>
>
>
=
=
= < < =
e
b
e
b
e
b
e
b y y y y .

Proof: When 1 < c , from equations (3.9)-(3.12), we have
a
c
y y
e
b
e
b
−
= =
>
=
=
=
1
0
0
0
0
0
0
,
0
0
0
0
0
0
1
>
>
>
> +
−
=
e
b
e
b x
a
b
a
c
y and
0
0
0
0
0
0
1
=
>
=
> +
−
=
e
b
e
b x
a
b
a
c
y . Since we have
0
0
0
0
0
0
>
=
>
> <
b
e
b
e x x ,
from the Result 3.5, we obtain
0
0
0
0
0
0
1
>
>
>
> +
−
=
e
b
e
b x
a
b
a
c
y <
0
0
0
0
0
0
1
=
>
=
> +
−
=
e
b
e
b x
a
b
a
c
y ,
when 0 , > g f .

When 1 < c and 0 = f , from equation (11) we have

ab
c agb c a gb c a gb
x
e
b
2
) 1 ( 4 )] 1 ( [ ) 1 (
2
0
0
0
− + − − + − −
=
=
>
ab
c a gb c a gb
ab
c a gb c a gb
2
) 1 ( ) 1 (
2
)] 1 ( [ ) 1 (
2
− + + − −
=
− + + − −
=
a
g
= .
It yields, from equation (3.43),
0
0
0
0
0
0
>
=
>
> <
b
e
b
e x x . If and only if when 1 < c and
0 = = f g , we get 0
0
0
0
0
0
0
= =
>
=
>
>
b
e
b
e x x , which gives
0
0
0
0
0
0
0
0
0
0
0
0
=
>
>
>
>
=
=
= = = =
e
b
e
b
e
b
e
b y y y y .
134
REFERENCES

1. Adi, B. I., An Iterative Method for Computing the Generalized Inverse of an
Arbitrary Matrix, Mathematics of Computation, 19 (91), 452-455. Jul., 1965.

2. Adi, B. I., Greville, T. N. E., Generalized Inverses 2
nd
Edition, Springer-Verlag
New York, Inc., 2003.

3. Albert, A., Regression and the Moore-Penrose Pseudoinverse, Academic Press,
New York, 1972.

4. Albert, A. and Sittler, R. W., A Method for Computing Least Squares Estimation
using Matrix Pseudoinverse that keep up with the Data, SIAM J. Contr., [A] , 3,
384-417, 1965.

5. Anderson, C. H. and Carter, J. R., On Rational Choice Theory and the Study of
Terrorism, Defense and Peace Economics, 16 (4), 275-82, 2005.

6. Appell, P., Su rune Forme Generale des Equations de la Dynamique, C. R. Acad.
Sci., Paris, 129, 459-460, 1899.

7. Atkinson S., Sandler T., and Tschirhart T., Terrorism in a Bargaining
Framework, Journal of Law and Economics, Vol. 30 (1), 1-21, 1987.

8. Blomberg, S. B., Hess, G. D., and Orphanides, A., The Macroeconomic
Consequences of Terrorism, Journal of Monetary Economics, 51 (5), 1007-32,
2004.

9. Blomberg, S. B., Hess, G. D., and Weerapana, A., An Economic Model of
Terrorism, Conflict Management and Peace Science, 21 (1), 17-28, 2004.

10. Boot, J. C. G., The Computation of the Generalized Inverse of Singular or
Rectangular Matrices, The American Mathematical Monthly, 70 (3), 302-303,
Mar., 1963.

11. Bouchard, M., Recursive Least-Squares Algorithms with Improved Numerical
Stability and Constrained Least-Squares Algorithms for Multichannel Active
Noise Control System, Journal of the Canadian Acoustical Association, 28 (3),
66-67, 2000.

12. Carpenter, G., Advanced Techniques for Modeling Terrorism Risk, Guy
Carpenter Co., 2000.

135
13. Chen, A. H. and Siems, T. F., The Effects of Terrorism on Global Capital
Markets, European Journal of Political Economy, 20 (2), 349-66, 2004.

14. Cooperman, R. L., Tactical ballistic missile tracking using the interacting
multiple model algorithm, Information Fusion, 2002. Proceedings of the Fifth
International Conference on, 2, 824-831, 2002.

15. D’Artigues, A. and Vignolo, T., Why Global Integration May Lead to
Terrrorism: An Evolutionary Theory of Mimetic Rivalry, Economic Bulletin, 6
(11), 1-8, 2003.

16. Dirac, P. A. M., Lectures in Quantum Mechanics, Yeshiva University, New York,
1964.

17. Enders, W. and Sander, T., The Political Economy of Terrorism, Cambridge,
Cambridge University Press, 2005.

18. Faria, J. R., Terror Cycles, Studies in Nonlinear Dynamics and Econometrics, 7
(1), article 3, 2004.

19. Fletcher, R., A Technique for Orthogonalization, J. Inst. Math. Appl., 5, 162-166,
1969.

20. Frey, B. S. and Luechinger, S., How to Fight Terrorism: Alternatives to
Deterrence, Defence and Peace Economics, 14 (4), 237-49, 2003.

21. Frey, B. S., Luechinger, S., and Stutzer, A., Calculating Tragedy: Assessing the
Costs of Terrorism, CESifo Working Paper Series No. 1341; IEW Working Paper
No. 205, Nov., 2004.

22. Gauss, C. F., Uber Ein Neues Allgemeines Grundgsetz der Mechanik. J. Reine
Ang. Math., 4, 232-235, 1829.

23. Geng, Z. J., Missile Control Using Fuzzy Cerebellar Model Arithmetic Computer
Neural Networks, Journal of Guidance, Control, and Dynamics, 20 (3), 557-565,
1997.

24. Gibbs, J. W., On the fundamental formulae of dynamics, Am. J. Math., 2, 49-64,
1879.

25. Graybill, F. A., Introduction to Matrices with Applications in Statistics,
Wadsworth, Belmont, California, 1969.

26. Graybill, F., Matrices and Applications to Statistics, 2
nd
Edition, Wadsworth,
Belmont, California, 1983.
136

27. Graybill, F. A., Meyer, C. D., and Painter, R. J., Note on the Computation of the
Generalized Inverse of a Matrix, SIAM Review, 8 (4), 522-524, Oct., 1966.

28. Greville, T. N. E., Some Applications of the Pseudoinverse of a Matrix, SIAM
Review, 2, 15-22, 1960.

29. Greville, T. N. E., The Pseudoinverse of a Rectangular or Singular Matrix and its
Application to the Solution of System of Linear Equations, SIAM Review, 1, 38-
43, Jan., 1959.

30. Haker, S., Sapiro, G., Tannenbaum, A., and Washburn, D., Missile Tracking
using Knowledge-based Adaptive Thresholding, Image Processing, 2001.
Proceedings. 2001 International Conference on, 1, 786 – 789, 2001.

31. Haykin, S., Adaptive filter theory, 3
rd
Edition, Prentice Hall, Englewood Cliffs
(NJ), 1996

32. Hsia, T. C., On Least Squares Algorithms for System Parameter Identification,
IEEE Transactions on Automatic Control, 104-108, Feb., 1976.

33. John, J. A., Use of Generalized Inverse Matrices in MANOVA, Journal of the
Royal Statistical Society. Series B (Methodological), 32 (1),137-143, 1970.

34. Kalaba, R. E. and Udwadia, F. E., Associative memory approach to the
identification of structural and mechanical systems, Journal of Optimization
Theory and App., 76, 207-223, 1993.

35. Kalaba, R. E., Xu, R., and Feng, W., Solving Shortest Length Least-Squares
Problems via Dynamic Programming, Journal of Optimization Theory and
Applications, 85 (3), 613-132, June, 1995.

36. Krueger, A. B. and Maleckove, J., Education, Poverty and Terrorism: Is There a
Causal Connection?, Journal of Economic Perspectives, 17 (4), 119-44, 2003.

37. Lagrange, J. L., Mechanique Analytique, Mme Ve Courcier, Paris, 1787.

38. Leontaris, G. K. and Rizos, J., A D-brane inspired Trinification model, Journal
of Physics: Conference Series, 53, 722-730, 2006.

39. Mizue, H. Ryo, T. and Morimitsu T., Embedding the Texture of the Neutrino
Mass Matrix into the MaVaNs Scenario, arXiv:hep-ph/0510018v2 9Dec2005, 1-
16, 2005.

137
40. Moore, E. H., On the Reciprocal of the General Algebraic Matrix, Bulletin of the
American Mathematical Society, 26, 394-395, 1920.

41. Noble, B., A Method for Computing the Generalized Inverse of a Matrix, SIAM
Journal on Numerical Analysis, 3 (4), 582-584, Dec., 1966.

42. Pars, L. A., A treatise on Analytical Dynamics, Oxbow Press, Woodridge, CT,
1979.

43. Penrose, R. A., Generalized Inverse for Matrices, Proceedings of the Cambridge
Philosophical Society, 51, 406-413, 1955.

44. Pitchford, I., http://www.interdisciplines.org/terrorism/papers/1/1/3 , 2005.

45. Pringle, R. M., and Rayner, A. A., Generalized Inverse Matrices, Hafner, New
York, New York, 1971.

46. Rao, C. R., A Note on a Generalized Inverse of a Matrix with Applications to
Problems in Mathematical Statistics, Journal of the Royal Statistical Society
Series B, 24, 152-158, 1962.

47. Rao C. and Mitra S., Generalized inverse of matrices and its applications, John
Wiley and Sons, 1971.

48. Rust, B., A Simple Algorithm for Computing the Generalized Inverse of a Matrix,
Communications of the ACM, 9 (5), May, 1966.

49. Saleem, M. and Vladimer, C., On recursive Calculation of the Generalized
Inverse of a Matrix, ACM Transactions on Mathematical Software, 17 (1), 130-
147, Mar., 1991.

50. Sandler, T., Collective versus Unilateral Response to Terrorism, Public Choice,
124 (1-2), 75-93, 2005.

51. Sandler T. and Arcy, D., Terrorism and Game Theory, Simulation and Gaming,
34, 2003.

52. Sandler, T. and Enders, W., An Economic Perspective on Transnational
Terrorism, European Journal of Political Economy, 20 (2), 301-16, 2004.

53. Sandler T., Tschirhart T., and Cauley J., A Theoretical Analysis of Transnational
Terrorism, American Political Science Review, 77 (4), 1983.

138
54. Siouris, G. M., Guanrong C., and Jianrong W., Tracking an Incoming Ballistic
Missile using an Extended Interval Kalman Flter, Aerospace and Electronic
Systems, IEEE Transactions on, 33 (1), 232 – 240, 1997.

55. Smith, S. E. and Yau, S. S., Linear Sequential Pattern Classification, IEEE
Transaction on Information Theory, 673-678, Sep. 1972.

56. Speed, F. M., An Application of the Generalized Inverse to the One-Way
Classification, The American Statistician, 28 (1), 16-18, Feb., 1974.

57. Thomas, L. B. and Patrick, L. O., Generalized Inverse Matrices, Wiley-
Interscience, 1971.

58. Udwadia, F. E., A New Perspective on the Tracking Control of Nonlinear
Structural and Mechanical Systems, Proceedings of the Royal Society of London,
Series A, 459, 1783-1800, 2003.

59. Udwadia, F. E., On Constrained Motion, Applied Mathematics and Computations,
164, 313-320, 2005.

60. Udwadia, F. E. and Kalaba, R. E., A General Forms for the Recursive
Determination of Generalized Inverses: Unified Approach, Journal of
Optimization Theory and Applications, 101 (3), 509-521, June, 1999.

61. Udwadia, F. E. and Kalaba, R. E., A Unified Approach for the Recursive
Determination of Generalized Inverses, Computers and Mathematics with
Applications, 37, 125-130, 1999.

62. Udwadia, F. E. and Kalaba, R. E., An Alternative Proof of the Greville Formula,
Journal of Optimization Theory and Applications, 94 (1), 23-28, 1997.

63. Udwadia, F. E. and Kalaba, R. E., Analytical Dynamics: A New Approach,
Cambridge University Press, England, 1992.

64. Udwadia, F. E. and Kalaba, R. E., Sequential Determination of the {1,4}-Inverse
of a Matrix, Journal of Optimization Theory and Applications, 117 (1), 1-7,
2003.

65. Udwadia, F. E. and Phohomsiri, P., Explicit Equations of Motion for Constrained
Mechanical Systems with Singular Mass Matrices and Applications to Multi-
Body Dynamics, Proc. R. Soc. of London A, 462, 2097-2117, 2006.

66. Udwadia, F. E. and Phohomsiri, P., Recursive Determination of the Generalized
Moore-Penrose M-Inverse of a Matrix, Journal of Optimization Theory and
Applications, 127 (3), 639-663, 2005.
139

67. Willner, L. B., An Elimination Method for Computing the Generalized Inverse,
Mathematics of Computation, 21 (98), 227-229, Apr. 1967.

68. Yu, F. and Bouchard, M., Recursive Least-Squares Algorithms with Good
Numerical Stability for Multichannel Active Noise Control, Proceedings of IEEE
International Conference on Acoustics Speech and Signal Processing (ICASSP),
5, 3221-3224, 2001.

69. Zhang, X. and Takeda, H., An Order Recursive Generalized Least Squares
Algorithm for System Identification, IEEE Transactions on Automatic Control,
AC-30 (12), Dec., 1985.

70. Zhou, J., Zhu, Y., Rong, L. X., and You, Z., Variant of the Greville Formula with
Applications to Exact Recursive Least Squares, SIAM J. Matrix Anal. Appl., 24
(1), 150-164, 2002.

Abstract (if available)

Abstract The concept of Dynamics exists everywhere in our lives. As an integral component of daily reality, engineers and scientists research, analyze, and observe the behaviors of dynamics and dynamical systems and their application in our real world activities. As a general rule, dynamical systems are characterized by their nonlinear quality. By their very nature, such systems have presented scientists with enormous difficulty in finding unequivocal solutions to systems, and even when a solution is possible, the path to finding one is extremely complicated. Oftentimes, it has proven virtually impossible to find these solutions. Given the inherent challenges in coming up with solutions to systems, the research that we've focused on has involved four separate but related areas.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

On the synthesis of dynamics and control for complex multibody systems

PDF

An analytical dynamics approach to the control of mechanical systems

PDF

Computation of algebraic system from self-assembly

PDF

On the synthesis of controls for general nonlinear constrained mechanical systems

PDF

Nonlinear control of flexible rotating system with varying velocity

PDF

New approaches in modeling and control of dynamical systems

PDF

Periodic solutions of flexible systems under discontinuous control

PDF

Resource allocation in dynamic real-time systems

PDF

Data-driven H∞ loop-shaping controller design and stability of switched nonlinear feedback systems with average time-variation rate

PDF

Scattering of a plane harmonic wave by a completely embedded corrugated scatterer

PDF

Dynamic modeling and simulations of rigid-flexible coupled systems using quaternion dynamics

PDF

Robust control of periodically time-varying systems

PDF

Organizing complex projects around critical skills, and the mitigation of risks arising from system dynamic behavior

PDF

Designing an optimal software intensive system acquisition: a game theoretic approach

PDF

Discrete-synapse recurrent neural network for nonlinear system modeling and seismic signal classification

PDF

Modeling human reaching and grasping: cortex, rehabilitation and lateralization

PDF

Characterization of environmental variability in identified dynamic properties of a soil-foundation-structure system

PDF

Multiple degree of freedom inverted pendulum dynamics: modeling, computation, and experimentation

PDF

Dynamical system approach to heart rate variability: QT versus RR interval

PDF

Techniques for increasing number of users in dynamically reconfigurable optical code division multiple access systems and networks

Asset Metadata

Creator Han, Byungrin (author)

Core Title Nonlinear dynamics and nonlinear dynamical systems

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Mechanical Engineering

Publication Date 06/07/2007

Defense Date 04/30/2007

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag least square,nonlinear control,nonlinear dynamical systems,nonlinear dynamics,OAI-PMH Harvest,terrorism,tracking control

Language English

Advisor Udwadia, Firdaus E. (committee chair), Dravinski, Marijan (committee member), Proskurowski, Wlodek (committee member)

Creator Email byunghan@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m515

Unique identifier UC1383761

Identifier etd-Han-20070607 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-501795 (legacy record id),usctheses-m515 (legacy record id)

Legacy Identifier etd-Han-20070607.pdf

Dmrecord 501795

Document Type Dissertation

Rights Han, Byungrin

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu