NONLINEAR DYNAMICS
OF A
KITE IN FLIGHT
by
Monica Geist .
B.S., Colorado State University, 1984
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Applied Mathematics
1996
This thesis for the Master of Science
degree by
Monica Geist
has been approved
by
William L. Briggs
^1
Geist, Monica (M.S., Applied Mathematics)
Nonlinear Dynamics of a Kite in Flight
Thesis directed by Professor William L. Briggs
ABSTRACT
In this paper we mathematically model a kite in flight. There are two versions of
the model, the basic model and the effective model. We develop the basic model
by studying the forces on the kite. We then describe the method of linearization
to study the behavior near fixed points. Then we apply the linear theory to the
basic model. Numerical methods are used to generate phase portraits of the
basic model. The effective model is derived by adding the effective wind speed
and the effective angle to the basic model. Again we use the linear theory,
plot the phase portrait and examine the behavior near the fixed points. We close
by comparing the two models.
This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed
William L. Briggs
DEDICATION
To Marc
CONTENTS
Acknowledgements xi
Chapter
1. Introduction.......................................................... 1
2. Basic Model of the Kite........................................... 3
2.1 Force Equations..................................................... 3
2.2 Proposed Bifurcation Diagram........................................ 6
2.3 Dynamics of the Kite................................................ 8
3. Analysis of the Basic Model ...................................... 11
3.1 Fixed Points....................................................... 12
3.2 Theory of Linearization............................................ 13
3.3 Classification of Linearized Systems............................... 15
3.4 Linearization ..................................................... 21
3.5 Numerical Solution of the Basic Model.............................. 23
3.6 Phase Portrait of the Basic Model.................................. 25
4. Effective Model of the Kite........................................... 28
4.1 Case of the Rising Kite............................................ 28
4.2 Case of the Falling Kite .......................................... 31
4.3 Putting It All Together............................................ 32
5. Analysis of the Effective Model....................................... 34
5.1 Fixed Points....................................................... 34
5.2 Linearization .................................................... 34
5.3 Actual Bifurcation Diagram for the Effective Model................. 36
5.4 Phase Portrait of the Effective Model.............................. 40
6. Conclusion........................................................... 46
v
Appendix
A. Calculated Values and C++ Code................................. 47
References ....................................................... 58
I
VI
FIGURES
Figure
2.1 The forces on the kite consist of wind, lift, drag, string tension and
gravity............................................................... 4
2.2 This is the plot of sin#s(l 2acos6s) = /?, where 6S is the angle
at which the kite is at equilibrium................................... 7
2.3 This is the graph of Figure 2.2 with the axes inverted. Drawn in
this fashion, it becomes the proposed bifurcation diagram of the basic
model of the kite......................................................... 8
2.4 This graph shows the components of force in the tangential direction.
The positive component of force is FL cos 6. The negative components
are Fq cos 6 and FD sin6.............................................. 10
3.1 A trajectory is a curve in the (8,0) phase space......................... 11
3.2 This graph shows the trajectories in the phase space near a saddle point. 16
3.3 Trajectories in the phase space are repelled from unstable nodes. ... 17
3.4 Trajectories in the phase space are attracted to stable nodes............ 17
3.5 This graph shows a trajectory in the phase space near an unstable spiral. 18
3.6 This graph shows the trajectories in the phase space near a stable spiral. 19
3.7 This graph shows two pure imaginary eigenvalues in the complex plane.
Nonlinear terms could have the effect of shifting them off the imaginary
axis..................................................................... 20
3.8 This graph shows a neutrally stable center of a linear system. Lin
earization of a nonlinear system cannot detect centers. However, if
the system of ODEs were linear to begin with, then two imaginary
eigenvalues indicate a center............................................ 20
Vll
3.9 This is the bifurcation diagram of the basic model according to lin
earization. We have determined that one of the fixed points is an
unstable saddle point and the other has unknown stability....... 23
3.10 Fourthorder RungeKutta method. In each step the right hand side of
the ODEs is evaluated four times: once at the initial point, 9n, twice
at midpoints, and once at the endpoint, 0n+i From these four values
the final function value (shown as a filled dot) is calculated........ 24
3.11 The RungeKutta method was used to calculate the trajectories for the
phase portrait of the basic model. The fixed point 9 = 0.0121 is a
saddlepoint and the other 9 = 1.1775 a center................... 26
4.1 When the kite is rising, the wind appears to be going down from the
perspective of the kite. The effective wind speed as viewed by the kite
can be expressed as A2 = (A (vsin#))2 + (0 vcos 9)2.............. 29
4.2 This graph shows the geometric relationships used to determine the
effective angle 9e for the case of the rising kite.................... 30
4.3 When the kite is falling, the wind appears to be going up from the
perspective of the kite. The effective wind speed as viewed by the kite
can be expressed as Ag = (A usm#)2 + (0 (vcos9))2............... 31
4.4 This graph shows the geometric relationships used to determine the
effective angle 9e for the case of the falling kite................... 32
5.1 This graph shows that B(a) > 0 for a > 0. Notice that B(\) = 0
and B(a) > 0 for a > 0................................................ 38
5.2 This bifurcation diagram summarizes the signs of 4> and Â£............. 39
5.3 This bifurcation diagram summarizes the type and stability of fixed
points, based on the determinant, A, and the trace, r of the Jacobian
matrix.......................................................... 39
vm
5.4 This phase portrait of the effective model shows trajectories that start
at angles between 0 and , which is the valid region for our kite model.
The fixed points occur at 0.44733 radians and 0.81 radians with /? =
0.607432............................................................ 40
5.5 The middle graph is a zoomed in view of the stable node. The upper
left picture is an angle versus time plot of the trajectory starting at
the upper left corner of the phase portrait (6,6) = (0.805,0.001) The
bottom middle graph shows what happens when the kite is placed pre
cisely on the stable node. The kite does not move from that position.
The fixed point occurs at 0.81 radians with /? = 0.607432........... 41
5.6 The phase portrait drawn here include the saddle point. Trajectories
that start out with a positive velocity are drawn toward the stable
node on the right. Trajectories starting with negative velocity are
drawn toward a stable node not in the valid physical boundaries of the
system. The three angle versus time plots show what happens when
the kite is placed at the three indicated points. The fixed point occurs
at 0.44733 radians with (3 = 0.607432............................... 43
5.7 This drawing shows the saddle point on the left and the stable node
on the right. Trajectories leaving the saddle point along the growing
eigenvector soon are drawn into the decaying eigenvector of the stable
node................................................................ 44
5.8 This phase portrait "zooms" in on the unstable node. A trajectory that
starts in the lower right corner goes down to zero, which physically
means that the kite hits the ground. Trajectories that start at the
fixed point as well as trajectories that start out in the upper right
corner also go down to zero, just not as fast. The fixed point occurs
at 0.65 radians with /? = 0.679557.................................. 45
IX
TABLES
Table
3.1 The determinant A, and the trace r, of the Jacobian matrix are used
to determine the stability of the system local to the fixed point. When
A = 0, t = 0, or t2 4A = 0 the linearization in inconclusive. ... 21
A.l This table summarizes the specific fixed used in the phase portraits. . 48
A.2 This table summarizes the eigenvalues and eigenvectors for each fixed
point calculated from the Jacobian matrices......................... 49
x
ACKNOWLEDGEMENTS
I would like to express my deep gratitude to Dr. Briggs for his inspiration
and guidance, his helpful suggestions and encouragement, for his time, energy
and many proofreading sessions. I would like to thank John Starrett for the
programming ideas. Finally, I would like to thank Marc, Amanda, Marissa,
Bettijean, Dodie and Dennis for their moral support and many pep talks.
XI
1. Introduction
Kites have been around for at least 2,000 years. We know that ancient
Chinese as well as ancient Egyptians flew kites. Through the ages, kites have
been used for religious celebrations, competitions, pleasure and experimentation.
Kites have also had some very interesting practical applications.
Chinese legend has it that the emperor of the Han dynasty attached
bamboo pipes to kites and flew them over enemy camps on moonless nights. The
wind blowing through the pipes made an eerie sound. The emperor sent
infiltrators into the camps to spread the word that the noise was that of ghosts
proclaiming the defeat of their invasion. The invaders were so scared they left in
a hurry.
As a boy, Benjamin Franklin used a kite to pull himself across a pond
while floating on his back. Later in his life he used kites to investigate
atmospheric electricity. Kite studies were also made by physicist and inventor
Alexander Graham Bell.
In 1825, English schoolteacher George Pocock invented a lightweight,
kitedrawn vehicle, which he named a charvolant This carriage drawn by two
kites could carry five people and reach speeds of 20 to 25 miles per horn.
Kites were used during the construction of the suspension bridge just
below Niagara Falls in New York state. Kites were flown across the river with
small ropes. When enough small ropes were flown across, they were used
together to pull cables across the river. This made it possible to begin
construction.
Samuel Franklin Cody invented a system of military observation which
the British Army used in 1906. A kite raised an observer in a basket with a
telescope, camera and gun. Notes were sent to and from the ground to relay
1
information. Enemy movements could be tracked from the air using this system.
The use of kites for military observation became obsolete with the
invention of the airplane. However, kites were still used in the early twentieth
century for military purposes. Kites were flown around military establishments
to keep out spy planes. The strings and wires attached to the kites would catch
airplane wings or entangle propellers.
The kite has been used to solve many other problems, from pulling boats
across bodies of water, to raising meteorological instruments for the study of
weather, from facilitating aerial photography to being used as signals in airsea
rescue operations.
In this paper we look at the problem of flying the kite. We will
mathematically model a kite in flight. There are two versions of the model. In
section 2 we develop the basic model by studying the forces on the kite. In
section 3 we describe the method of linearization to study the behavior near
fixed points. Then we apply the theory to our basic model. After exhausting
qualitative methods, we use a computer to solve the system of equations
numerically. This gives us a big picture of the basic systems behavior. In
section 4 we refine the basic model by adding effective quantities to produce the
effective model. Then in section 5 we do both a qualitative study and a
numerical study, of the effective model. We will close by comparing the two
models.
2
2. Basic Model of the Kite
A kite is a windsupported flying device consisting of a wooden or similar
framework covered with paper, cloth or synthetic material. A kite may be
thought of as a rudimentary airfoil. In this section of the paper we will develop a
system of equations that model the basic forces on a kite. Figure 2.1 shows the
forces that we will be looking at in this first model. There is the wind, the
tension of the string, the force due to lift, the force due to drag, and the force
due to gravity. We will define the following quantities,
T = the tension in the string,
Fd = the force due to drag,
Fl = the force due to lift,
Fg = the force due to gravity,
A = the wind velocity,
0 = the angle between the ground and the kite string.
2.1 Force Equations
The movement of the kite in this model will be measured by the angle
between the string and the horizontal. The kite is assumed to be fixed to the
string and does not rotate around the string. Other ingredients needed in the
equations of motion are as follows;
3
Fl
Figure 2.1: The forces on the kite consist of wind, lift, drag, string tension and gravity.
p = the density of air,
A the surface area of the kite,
m = the mass of the kite,
CD = the coefficient of drag dependent on the angle 0,
Cl = the coefficient of lift dependent on the angle 0.
The vertical forces on the kite are the ycomponent of the tension in the
string, Ty, the force due to gravity on the kite, FG, and the force due to lift, FL.
These forces must balance, which means that
Ty + FG = Fl (1)
The horizontal forces are the ^component of the tension of the string, Tx, and
the force due to drag, F&. Again, these forces must balance, giving
Tx = Fd. (2)
4
The ^component of the tension in the string is T sin 9 and the ^component is
T cos 6. Equations (1) and (2) become
Tsm9 + FG = FL (3)
and
T cos 9 = Fd (4)
The wind that keeps the kite afloat has a speed of A. The force due to lift is
given by
Fl = \pCL\2A. (5)
Likewise the force due to drag is given by
Fd = \PCD\2A (6)
[Adomaitis, ODell, Marchaj]. The lift coefficient, Cd, is usually determined
experimentally in a wind tunnel. For an angle less than  the lift coefficient can
be modeled by
Cl = b sin 29,
where b is a constant which is determined experimentally [Adomaitis, ODell].
The drag coefficient, Cd, can be expressed as
Cp = a cos 9,
where a is also a constant, also determined experimentally [Adomaitis, ODell].
The force due to gravity on the kite, FG, is the mass of the kite, m, times the
acceleration of gravity at the earths surface, g. Substituting FG, Fl, Fd, Cl
and Cd into equations (3) and (4) gives
Tsin 9 + mg = sin 29 (7)
z
and
T cos 9 ^pa)?Acos9. (8)
Z
5
Solving equation (8) for T .and substituting into equation (7) results in
paX2A) sin 9 + mg = pb\2A sin 26.
z z
Further simplification gives
/, b o/. ~2m9
sm6sin 26 = rr.
a apXzA
We use the trigonometric identity sin 26 = 2 sin 9 cos 6 to get
an 26 a\ ~2m9
smcqlcos 6) = tttt
a ap\2A
If we now let
6 2rag
a = and p
a ap\2A
and realize that all the forces are balanced, equation (9) becomes
(9)
sin0s(l 2acos9s) = /?,
(10)
where 6S is the angle at which the kite is at equilibrium.
2.2 Proposed Bifurcation Diagram
Plotting equation (10) helps determine what information the equation
has to offer. To determine the maximum value of /? as a function of 9, the
derivative is set equal to zero:
sin 6(2a sin 6) + (1 2a cos 6)( cos 6) = 0.
dd
This condition reduces to
2a(cos2 6 sin2 6) cos 9 0.
Using the trigonometric identity cos 29 = cos2 9 sin2 9 gives
2a cos 29 = cos 9. (11)
6
Figure 2.2: This is the plot of sin0s(l 2a:cos0s) = [3, where 6S is the angle at
which the kite is at equilibrium.
Let 0* satisfy equation (11) corresponding to the maximum value of (3 which we
will call (3*.
Given certain parameters in the system, there are positions in which the
kite is not moving. These points are called fixed points or equilibrium points.
Perturbing the system slightly from these fixed points causes one of two things to
happen. If the system drifts back to this fixed point, the point is a stable point.
If the system drifts away from the fixed point, the point is an unstable point.
The system is attracted to stable points and repelled from unstable points.
Fixed points can be created, destroyed, or have their stability changed
when the parameters are changed. In dynamics, these qualitative changes are
called bifurcations The parameter values at which these changes occur are
called bifurcation points The graph in Figure 2.2 can be made into a
bifurcation diagram by inverting the axes in order to view (3 as the independent
variable as shown in Figure 2.3. We see that for 0 < (3 < /?*, there are two
possible equilibrium values (fixed points) of 6. The point (/?*, &*) in the graph
of Figure 2.3 is called the turning point.
Observe that when the wind has a very high velocity,
(3 = sin0s(l 2a cos 9S) > 0. In this limit, two conditions are possible. Either
7
Figure 2.3: This is the graph of Figure 2.2 with the axes inverted. Drawn in this
fashion, it becomes the proposed bifurcation diagram of the basic model of the kite.
9 = 0 or 1 2a cos 0 = 0 for some 9 Â£ (0,7t/2). The latter case occurs if and
only if a > . Therefore, in order to achieve multiple equilibrium states,
b 1
a = > 9
a 2j
On the other hand, when the wind has a very low velocity,
/? = sin 0S(1 2a cos 9S) oo, there are no solutions to equation (10). In fact,
this condition occurs when f3 > f3*.
2.3 Dynamics of the Kite
Newtons second law of motion, F = ma, states that the sum of all
external forces on an object is the mass of the object times its acceleration. The
tangential velocity, v, of a mass moving on a circular path is the radius of the
circle, r, times the rate of change of the angle. Letting the rate of change of the
angle with respect to time be
d9
8
the tangential velocity is
v = rr7,
and the tangential force is
rnvrj.
Computing the components of force in the tangential direction, there are
three components that make up the force of the nonstationary kite. Referring
to Figure 2.4, the positive component of force is Fl cos 8. The negative
components are Fa cos 6 and Fq sin#. The tension in the string, T, is
perpendicular to the tangent, therefore there are no components of T in the
tangential direction. Combining these forces gives,
mrr) mg cos 8 + FL cos 8 FD sin 9.
9
;
Figure 2.4: This graph shows the components of force in the tangential direction.
The positive component of force is FL cos 9. The negative components are Fc cos 9
and Fd sin0.
Using equations (5) and. (6), we see that
1 1
7777.77 = mg cos 9 + (pAX2b sin 29) cos 9 (pAX2a cos 9) sin 9.
Z Z
Recall that 77 = 9. We can introduce F(9,rj) and G{9,rj) and write the system as
e = F(e,n) = ri (12)
and
77 = G(9,77) = sin 29 cos 9 cos 9 sin 9 cos 9. (13)
2rm 2rm r
Equations (12) and (13) are the ordinary differential equations that describe the
tangential motion of the kite under the assumptions of this model.
10
3. Analysis of the Basic Model
Now that the governing system of ODEs has been defined, we observe
that analytical solutions are impossible to obtain. We will have to resort to
other methods for understanding the behavior of the system.
Suppose we set the kite in motion at a certain angle at time = 0, 0(0),
with a certain angular velocity, 0(0). We would like to be able to predict what
the kite will be doing at any future time. At a time t, the position and angular
velocity of the kite will be represented by the functions 0(Â£) and 0(Â£). By
constructing an abstract space, as in Figure 3.1, with coordinates (0,0), the
solution (0(t),0(t)) corresponds to a point moving along the curve. This curve is
called a trajectory and the space is called the phase space. Since each point in
the phase space can serve as an initial condition, the phase space is completely
filled with trajectories. [Strogatz]
Figure 3.1: A trajectory is a curve in the (0,0) phase space.
11
Our goal is to build the phase space for our system of ODEs. The first
step is to determine the location of the fixed points. Next, we will qualitatively
determine the behavior near the fixed points. Then we can give the system
several different initial conditions and plot the trajectories.
3.1 Fixed Points
A fixed point, or equilibrium point, occurs when there is no movement in
the system, that is, when 9 rj = 0. Imagine placing the kite in the air at a
fixed point and then perturbing it just a little. If the kite oscillates a few times
and then returns to the fixed position, then the fixed point is stable. If the kite
moves away from the fixed point, it is unstable.
So, how do we qualitatively determine the behavior near fixed points? We
will use the method of linearization. But first, some explanation is required. Let
(0S, rjs) be a fixed point and let (6p,T)p) be a new location relative to the fixed
point after slightly perturbing the kite. We have
0p{t) = 6(t) 6S and r)p(t) = r](t) rjs. (14)
To see whether the perturbation grows or decays, and to determine the rate at
which it grows or decays, we substitute and note 0a = 0. Recall that 6 = F(9,r]),
so
6 = Op + 9s = 9p = F(0S + 0P, T)s + Tjp).
Now we use a Taylor series expansion around the fixed point (9s,rjs),
qp rlF
. 0P = F(9s,r}s) + 9p(0s,rjs)+riP{0s,ris) + 0(0^, 0p?7p), (15)
where 0{9^,r}^,9prip) denotes quadratic terms in 9P and r]p. Since our
perturbations are assumed to be very small, these quadratic terms are extremely
12
small. Note also that F(9p,r)p) = 0 since (0s,r)s) is a fixed point. Thus, equation
(15) becomes
dF
dF
0P = ep^(0s, Vs) + Vp^r + O(0p, rjp, 9pr)p).
d9 v &q
We can differentiate r]p similarly to get
dG
dG
Vp 0p (9s,r}g) F Vp n (Os,Vs) T O(0p,Vpi OpVp)
d9 v IP
In matrix form, these equations appear as
^ ^(Os,Vs) %(0s,Vs)
\
dG i
U '
p
\Vp )
+ quadratic terms
(16)
\ 80(^5 Vs) Qv(0S,Vs) )
Dropping the quadratic terms in equation (16), since they are extremely small,
linearizes the system. The matrix
f dF dF \
de dr)
dG dG
\ de dr) )
is called the Jacobian matrix at the fixed point (9s,r)s).
3.2 Theory of Linearization
At this point it behooves us to discuss the theory of linearization after
which we will return to the kite problem.
Suppose we are given a system
M
U/
which can be written in vector form as
z = Az,
(17)
13
where
z =
U
\y
and A
f a b ^
c d
and z =
\J
The matrix A can be thought of as the Jacobian matrix of the kite system. The
solution to equation (17) has the form
z = VeAt,
where A and t are scalars and
\V2)
Calculating the derivative of z and substituting into equation (17) results in
z = VXext = AVXe
At
Dividing both sides by ext / 0 we get
AV = XV,
(18)
which is an algebraic eigenvalue problem, where A is the eigenvalue and V is the
eigenvector. We now write equation (18) as
(A A/)V = 0,
where I is the 2x2 identity matrix. Notice that
(aA b ,
(A AI)V =  V = 0.
\ c dA
In order for a nontrivial solution to exist, the determinant of (A XI) must be
zero,
(a X)(d X) be = 0.
14
This condition results in the characteristic equation
A2 (a + d) A + (ad be) = 0.
Notice that a + d is the trace, r, and ad be is the determinant, A, of A.
Equation (19) can be written as
A2 r A + A = 0.
Using the quadratic formula to solve for the eigenvalues gives
? t vV2 4A
(19)
(20)
Recall that the solution to the system will be of the form, z = ~Vext. We will
refer to this form of solution as an eigensolution. Observe that the eigenvector V
gives us the direction of motion of the trajectory, and ext indicates that the
motion is exponential, with Re (A) as the growth or decay rate.
As a result, the behavior of the trajectories near a fixed point depends on
the eigenvalues. However, it is not necessary to compute the eigenvalues. We
notice from equation (20) that we can tell everything we need to know from the
trace and the determinant of the matrix A. It is important to note, though, that
any trajectory behavior knowledge gained in this fashion, is only local to the
fixed point. We will discuss the reasons for this later in the paper.
3.3 Classification of Linearized Systems
Before linearizing the system of ODEs that model our kite, let us discuss
the different possibilities that can occur with the trajectories. Recall that,
A =
r \/r2 4A
15
First, we will discuss the case in which A < 0 which results in one positive real
and one negative real eigenvalue. Hence, the first eigensolution
z = V+eX+t
grows exponentially and the second,
z = V_e*t
decays exponentially. The fixed point in this situation is known as a saddle
point. Observe in Figure 3.2 that unless the initial condition is on the
eigenvector associated with the negative eigenvalue, the resulting trajectory is
repelled from the fixed point. Therefore, a saddle point is unstable.
Figure 3.2: This graph shows the trajectories in the phase space near a saddle point.
The second case is A > 0. We then need to look at the sign of r and
r2 4A. If r2 4A > 0 and r > 0, then both eigenvalues are positive real
numbers. In this situation, the eigensolutions grow and all of the trajectories are
repelled away from the fixed point (Figure 3.3). The fixed point is known as an
unstable node.
16

X
Figure 3.3: Trajectories in the phase space are repelled from unstable nodes.
Another situation occurs when A > 0 and r2 4A > 0, but r < 0. Both
eigenvalues are negative real numbers in this case, the eigensolutions decay and
the trajectories approach the fixed point as shown in Figure 3.4. The fixed point
in this case is known as a stable node.

x
Figure 3.4: Trajectories in the phase space are attracted to stable nodes.
The saddle point, the unstable node and the stable node all have distinct
real eigenvalues. Now, lets look at what happens when the eigenvalues are
17
Figure 3.5: This graph shows a trajectory in the phase space near an unstable spiral.
complex which occurs when the A > 0 and r2 4A < 0. Recall that complex
eigenvalues come in complex conjugate pairs as given by
Assuming that \/4A t2 ^ 0 and that r > 0, we notice the real part of A is
positive, indicating a growing solution. These eigensolutions grow by oscillation,
as in Figure 3.5. The imaginary part gives the frequency of oscillation. The fixed
point is called an unstable spiral.
On the other hand, assuming that >/4A r2 ^ 0 and r < 0, the real part
of these eigenvalues are negative which indicates decay. Again the imaginary
part gives the frequency of oscillation. The eigensolution in this case is attracted
to the fixed point, as shown in Figure 3.6. This fixed point is called a stable
spiral.
Since our system is nonlinear, we will have to drop what we assume are
small terms from the equation when linearizing the equation. So, how does this
linearized analysis apply to our nonlinear system? In the cases of the
saddlepoint, stable and unstable nodes, stable and unstable spirals, the type and
18
Figure 3.6: This graph shows the trajectories in the phase space near a stable spiral.
stability of fixed points of a nonlinear system is given by the linear analysis.
However, there are some cases in which the linear analysis can not determine the
stability of fixed points of a nonlinear system.
One such scenario occurs when both eigenvalues are pure imaginary
numbers. Perturbing the fixed point has the effect of shifting them off the
imaginary axis, as shown in Figure 3.7. Therefore, if the eigenvalues of a
linearized nonlinear system are pure imaginary, no conclusions can be drawn
about the stability [Strogatz].
However, if the system is truly linear and both eigenvalues are pure
imaginary, the fixed point is known as a center. The fixed point in this
particular case is surrounded by a family of closed orbits. The stability of a
center is called neutral since nearby trajectories are neither attracted to nor
repelled from the fixed point as shown in Figure 3.8.
The linearized analysis of a nonlinear system is summarized in Table
(3.1).
19
Im
Figure 3.7: This graph shows two pure imaginary eigenvalues in the complex plane.
Nonlinear terms could have the effect of shifting them off the imaginary axis.
Figure 3.8: This graph shows a neutrally stable center of a linear system. Linearization
of a nonlinear system cannot detect centers. However, if the system of ODEs were
linear to begin with, then two imaginary eigenvalues indicate a center.
20
Determinant A and the Trace r Eigenvalues of the Jacobian Matrix Stability of Fixed Point of Nonlinear System
A <0 Real with opposite signs Unstable Saddle Point
A > 0 r2 4A > 0 r < 0 Real, Both Negative Stable Node
r = 0 Both Pure Imaginary Unknown
T > 0 Real, Both Positive Unstable Node
r2 4A < 0 r < 0 Complex Conjugates with Negative Real Parts Stable Spiral
T = 0 Both Pure Imaginary Unknown
r > 0 Complex Conjugates with Positive Real Parts Unstable Spiral
r2 4A = 0 Both Equal Unknown
A = 0 At Least One is Zero Unknown
Table 3.1: The determinant A, and the trace r, of the Jacobian matrix are used to
determine the stability of the system local to the fixed point. When A = 0, r = 0, or
r2 4A = 0 the linearization in inconclusive.
3.4 Linearization
Now that we are equipped, to linearize a system and study the behavior
around its fixed points, lets return to equation (16). Observe that
^ = ^ = 0 andthat ^ = 1.
86 dr] dr]
The Jacobian matrix then becomes,
/ dF OF \ f 0 1 \
ae dr)
d G dG dG 0
V ae dr) ) \ d0 /
Notice that the determinant is A = ^ and the trace is r = 0. In order to
determine the sign of A we must compute
dG d (pAX2b 0/) pAX2a n n 9 a\
= sin 29 cos 9cos 9 sin 9cos 9 .
d9 d9 y 2 rm 2rm r J
(21)
21
Setting
apAX2 b 1
7 =  and a = >
2 rm a 2
equation (21) becomes,
dG
7TT = 7 (a sin 29 sin 9 + a cos 9 cos 29 cos2 9 + sin2 9\ + sin 9. (22)
o9 v 'r
Using trigonometric identities, we find that
dG~
= 7 (2a sin2 9 cos 0 + 2a: cos 9 cos 29 2 cos2 0 + l) + sin 9
r)H V / r
Factoring out cos 0 and substituting
P =
2 mg
ap\2A
yields
dG
= 7 cos 9 (2a sin2 9 + 2a cos 29 2 cos 0 \^ ^
d9 V cos 9 cos 0 y
From equation (10) we know that at equilibrium (3 = sin 0S(1 2acos0a).
Substituting this value of (3 gives us
dG(93,r)3)
d9
= 7 cos 9s(2a cos 29s cos 9S).
We need to determine the sign of Observe that the constant 7 is positive.
, <  sc
determined by the sign of
Also, assuming that 0 < 9a <  so that cos 9S > 0, the sign of ^ can be
2a cos 29 q cos 9S.
Notice also that at equilibrium
(3 = sin 9S( 1 2a cos 9S).
If we take the derivative of (3 with respect to 9S we get
= cos 9S + 2a cos 29 s.
d9a
(23)
(24)
22
Figure 3.9: This is the bifurcation diagram of the basic model according to lineariza
tion. We have determined that one of the fixed points is an unstable saddle point and
the other has unknown stability.
Notice that (23) and (24) involve the same expression! So the sign of j^ is the
same as the sign of which we recall from Figure 2.3 as,
+1 for0s<0*
1 for 6a > 6*
However, since A = ^, the sign of A is
( +1 for 6S > 6*
sgn A = <
y 1 for 6S < 6*
Notice from Table (3.1) that if r = 0 we cannot determine the stability.
So the only valid piece of information we can acquire from linearizing this
system is that when sgn A = 1, the fixed point is an unstable saddle point. So
now our bifurcation diagram, Figure 3.9, looks a little bit more interesting.
3.5 Numerical Solution of the Basic Model
We have found the fixed points and have qualitatively determined the
behavior near the fixed points. Now we can give the system several different
sgn
dfi
dd.q
=
23
1
2
0,
n
e,
n+1
3
4
Figure 3.10: Fourthorder RungeKutta method. In each step the right hand side of
the ODEs is evaluated four times: once at the initial point, 9n, twice at midpoints, and
once at the endpoint, 0n+1. From these four values the final function value (shown as
a filled dot) is calculated.
initial conditions and plot the trajectories in the phase space. We will use a
fourthorder RungeKutta method. In each step, the functions F and G, (the
right hand sides of (12) and (13)) are evaluated four times: once at the initial
point, twice at midpoints, and once at an endpoint (Figure 3.10). From these
calculations, the solution is updated [Press]. Recall the definitions of F and G,
e = F(6,ri)=n
(25)
and
ft = G(9, ri) = sin 2 6 cos 6 cos 9 sin 9 cos 9. (26)
2 rm 2 rm r
The method requires computing the following quantities at each time
step:
and
Vn+l =Vn+c(w 1 + 2w2 + 2 S3 + W4)
0
24
where h is the step size and
v1 = hF(6n,r)n),
wi = hG(6n,r]n),
v2 = hF{6n + \vUT]n + \wl),
w2 = hG(On + \vi,r)n + \wi),
v3 = hF(dn + \v2,'qn + \w2),
w3 = hG{Bn + lv2,rjn + \w2),
v4 = hF(6n + v3,r}n + w 3),
w4 = hG(6n + V3 ,TJn + w3).
The computer code was written in C++ and can be found in Appendix A.
3.6 Phase Portrait of the Basic Model
We use the fourthorder RungeKutta method to fill the phase space with
several trajectories using fixed values of a and (3 (Figure 3.11).
25
Saddle Point Center
Figure 3.11: The RungeKutta method was used to calculate the trajectories for the
phase portrait of the basic model. The fixed point 9 = 0.0121 is a saddlepoint and
the other 6 = 1.1775 a center.
The linearization indicated that one fixed point is a saddle point, which
can be seen near the origin. The linearization was not conclusive on the other
fixed point. If the system of ODEs would have been linear to begin with, this
fixed point would have shown up as a center as discussed in section 3.3. Notice
though, that the phase portrait does show a center around the other fixed point.
What is the physical interpretation of the phase portrait? When the
initial angle of the kite is less than 6*, the equilibrium is a saddle point. Some
trajectories will be attracted to the equilibrium point then turn away for it
before actually arriving at it. Some trajectories will be repelled. The only
trajectory that will be attracted to and stay attracted to the equilibrium point
26
are those that start on the decaying eigenvector. In that case, the kite will find
and stay at the equilibrium point. Otherwise, the kite will not remain stable.
On the other hand, when the initial angle of the kite is greater than 6*,
the fixed point appears to be a center. If the kite is placed precisely on the fixed
point, it will remain there. If the kite is placed slightly off of the fixed point, in
theory, it will oscillate up and down with constant period and frequency.
27
4. Effective Model of the Kite
The basic model of the kite lacks additional effective quantities
introduced by movements of the kite. As the kite moves along a fixed arc, the
effective wind speed and effective angle of the wind changes. These
quantities introduced by the dynamics of the kite can greatly effect the stability
of the system. We will now enhance our basic model with these effective
quantities.
To understand these effective quantities, we consider a simple example.
Suppose you are standing facing a wind which is blowing at 30 miles per hour.
The wind that you feel is 30 miles per hour. Now suppose you run into the wind
at 3 miles per hour. You now feel 33 miles per hour of wind. This relative wind
speed you feel can be considered the effective wind speed In vector language,
the relative velocity between two things is the difference in their vectors.
In this section we will develop a slightly different model of the kite, one
that includes these effective quantities [Adomaitis].
4.1 Case of the Rising Kite
Consider the kite increasing in altitude along a circular arc. The climbing
kite will have a positive value for 6. The wind speed, A, can be expressed in the
(x,y) coordinate system as (A, 0). The kite speed, v, is the product of the
angular velocity and the radius of the arc, v = r6. The x and y components of
the kites velocity can be expressed as (v sin 6, v cos 6) as shown in Figure 4.1.
The effective (or relative) wind speed as viewed by the kite, Ae, can be
28
vsin 9
Figure 4.1: When the kite is rising, the wind appears to be going down from the
perspective of the kite. The effective wind speed as viewed by the kite can be expressed
as A2 = (A (vsin9))2 + (0 ucos0)2.
expressed, as
Ag = (A (vsin.9))2 + (0 vcosQ)2 (27)
which simplifies to
A2 = A2+u2 + 2uAsin0. (28)
The direction of the wind as seen by the kite, or effective angle, is more
subtle to resolve. First, consider Figure 4.2 which depicts the relative wind
speed, Ae, as viewed from the kites perspective. When the kite is rising, the
wind appears to be going down from the perspective of the kite. The effective
angle, 6e, is the sum of 0 and the angle 0 where
6 = tan 1
V cos 9
A + v sin 9'
The resultant sum is
9e = 9 + tan 1
V cos 0
A + v sin 9
(29)
29
Figure 4.2: This graph shows the geometric relationships used to determine the
effective angle 9e for the case of the rising kite.
Therefore, 9e is the effective angle between the effective wind speed and the
kite string.
Observe another relationship that contains 9e that we will be using later
in the analysis of the effective model,
Ae cos 9e = Ae cos (9 + 9) Ae(cos 9 cos 9 sin 9 sin 9). (30)
Notice from Figure 4.2 that
 A + v sin 9  v cos 9
cos 9 : and sin 9 =
Ae
Substituting into equation (30) results in
A + usin 9'
Ae
Ag cos 9g Ag
which reduces to
COS0
sin0
rv cos 9'
Ag .
Ae cos 9P = A cos 9.
(31)
30
vsin 0
Figure 4.3: When the kite is falling, the wind appears to be going up from the
perspective of the kite. The effective wind speed as viewed by the kite can be expressed
as A2 = (A usin#)2 + (0 (vcosd))2.
4.2 Case of the Falling Kite
For completeness, we will look at the case where the kite is decreasing in
altitude along a circular arc. The falling kite will have a negative value for 9.
When the kite is falling, the wind appears to be going up from the perspective of
the kite. The wind speed is still (A, 0) but now the x and y components of the
kites speed can be expressed as ('usin#, v cos 6), as shown in Figure 4.3. The
effective wind speed as viewed by the kite, Ae, is now
A2 = (A v sin 0)2 + (0 (v cos 6))2.
which simplifies to
A2 = A2 2uA sin 9 + v2.
The effective angle in this case is similar to the rising kite. Figure 4.4
depicts the relative wind speed, Ae, as viewed from the kites perspective. When
the kite is falling the wind appears to be going up from the kites perspective.
31
y
X vsin 0 vsin 0
Figure 4.4: This graph shows the geometric relationships used to determine the
effective angle 0e for the case of the falling kite.
The effective angle 6e is now the difference of 6 and the angle 6, where
v cos 6
6 = tan 1
A v sin 6
The resultant difference is
0e 6 tan 1
v cos 6
A u sin 6
As with the rising kite,
Ae cos 0e = A cos 0,
can also be shown analytically.
4.3 Putting It All Together
We have been using v to represent the kite speed. In order to use rrj for
the kite speed, notice that the sign of 7/ reflects the kite direction changes. The
equation for the effective wind speed in both the case of the rising kite and
the case of the falling kite is represented by
Ag = A2 + 2rrj\sin 0 + r2rf. (32)
32
The equation for the effective angle in both case is
9e = 9 + tan 1
rr) cos 9
X + rr] sin 9
(33)
Substituting these effective quantities into all of the force terms of equation (13)
gives
pAX2eb . . pAX2ea g
77 = sin 29e cos 9ecos 9e sin 9ecos 9,
2 rm Arm r
where Xe and 9e are defined by equations (32) and (33). Substituting Ae and
rearranging gives the second equation in our effective model,
9 = F{9,rj)=r/
(34)
and
r> = G(9,ri) = (A2 + r2rt2 + 2rXn sin 9) (o; sin 29e cos 9e cos 9e sin 9e) cos 9.
v 2 rm r
(35)
33
5. Analysis of the Effective Model
Just as in the basic model, this system is impossible to solve analytically.
We will follow the same steps as in the basic model. We will find the fixed
points, then linearize the system around those fixed points. We will numerically
plot the trajectories in the phase space and close by comparing the basic model
to the effective model.
5.1 Fixed Points
Recall that a fixed point occurs when there is no movement in the
system. Since 9 and rj both describe movement, fixed points occur when
9 = rj = 0. The angle of the fixed point in this effective model satisfies
equation (10), so the coordinates of the fixed point are (9S, 0).
5.2 Linearization
The Jacobian matrix for the system (34) and (35) when 9 = 9a and
rj =r)s 0 becomes
( 1
where
Observe that
* dG , . dG
$ = and f = .
o9 orj
Ae(9,0) = A and cos9e(9,0) = cos 9.
34
Before calculating $, we need to compute two ingredients. First note that
9A
9e = cos 1 cos 9
Using the fact that Ae(9,0) = A we find that
d6e
= 1.
ae
(36)
After some computation, we have
$ = 7 a sin 29e sin 9e + 2a cos 9e cos 29e cos2 0e + sin2 9^j + sin 9.
We can eliminate the effective angles by recalling that cos 9e(0,0) = cos 9 which
gives us
$ = 7 a sin 29 sin 9 + 2a cos 9 cos 29 cos2 9 + sin2 9^j + sin 9,
which is exactly equation (22) of the basic model. Following the same steps we
used on equation (22), we arrive at
$(8s,V s) = 7 cos 9s(2a cos 29 s cos 9S).
We now turn to the computation of Â£(0S, 0). Two derivatives that we
need are
d\
d9P
enm and W((,0)
Recalling that Ae = (A2 + r2rf + 2r\r) sin 9) 2, it is straigtforward to find that
d\e
Using equation (31),
7 = r sin 9.
drj
9e = cos 1 f cos 9 J ,
35
the calculations are again straight forward and lead to
We can summarize these Jacobian calculations as follows:
$(0,0) = 7 cos 9(2a cos 29 cos 9) and Â£(9,0) = cXr cos 9 (2a cos3 9 l)
(37)
where
apA
2 rm
5.3 Actual Bifurcation Diagram for the Effective Model
Now that $ and 77 are computed, we return to the Jacobian matrix,
and notice that the determinant of this matrix is $ and the trace is Â£. In order
to determine the signs of the determinant and trace, we need to determine the
signs of $ and Â£. We will first consider
$ = 7 cos 9s(2a cos 29 s cos 9S),
and notice that 7 cos 9S > 0 assuming 0 < 9S < So the sign of $ is the sign of
J(@Sl Vs)
2a cos 29s cos 9S.
(38)
At equilibrium, f3 = sin #s(l 2a: cos 9S) and
36
From equation (38), we see that the sign of $ = the sign of Jp. Recall the graph
of (3 from Figure 2.3. When Qs < 9*, [3 is increasing and the derivative is
positive. When 9a > 9*, (3 is decreasing and the derivative is negative. Therefore,
'
sgn $ = <
+1
1
for 9S < 9*
for 9S> 9*.
However, since the determinant, A = ,
sgn A = <
+1
1
for 9S > 9*
for 9S <9*.
The sign of the trace of the Jacobian, Â£, on the other hand is more
complicated. Recall
Â£ = cXr cos 9 (2a cos3 9 1),
where
apA b
c and a = .
Arm a
Since c > 0, a > , and O<0<, Â£ = 0 implies that
2a: cos3 0 1 = 0.
This equation has roots 9b that satisfy
cos 9 b =
1
(2a:) 5
(39)
Lets set equation (39) aside for a moment and notice the turning point
in the bifurcation diagram. We have already let 9* be the angle at which the
turning point occurs; 9* satisfies
2a cos 29* cos 9* = 0.
37
To solve for 9* we use a trigonometric identity and arrive at the quadratic
equation,
4a cos2 9* cos 9* 2a = 0,
to which we can apply the quadratic formula to find that
cos 9* = +
8 a
7 1 \2 1
V8a) +2
We keep the positive root because 0 < 9* < We need to consider whether 9b is
greater than or less than 9*. Observe that if 9b > 9*, then cos 9b < cos 9*. We
then define a new function B{a) by
B(a) = cos 9* cos 9b = 77 +
8a
()2 + 
\8a/ 2
1
(2a)
Observe that B(^) = 0. Also notice that lim^oo B(a) Notice the plot of
Figure 5.1: This graph shows that B(a) > 0 for a > 0. Notice that _B() = 0 and
B(a) > 0 for a > 0.
B(a) in Figure 5.1. When a > 0, B(a) > 0. When a = then B(a) = 0. Since
a > \ then 0 < B{a) < ^= and cos 6* > cos 6b which leads us to the relationship
9b > 9*. (40)
38
Lets turn our attention back to Â£. We evaluate the sign of Â£ for 6 > 9b
and 6 < 6b,
( +1 for 9 < 6b
sgn f = <
^ 1 for 6 > 6b.
The graph in Figure 5.2 summarizes the signs of $ = ^ and Â£ =
Figure 5.2: This bifurcation diagram summarizes the signs of $ and Â£.
The signs of $ and Â£ allow us to determine the signs of the determinant, A, and
the trace, r. Figure 5.3 summarizes the stability characteristics.
Figure 5.3: This bifurcation diagram summarizes the type and stability of fixed points,
based on the determinant, A, and the trace, r of the Jacobian matrix.
39
5.4 Phase Portrait of the Effective Model
We will look at four interesting aspects of the bifurcation diagram. First,
the linearization study told us that there is a 9 > which is a fixed point and
that fixed point is either a stable node or a stable spiral. We chose a (3 that
corresponds to such a 9 (see Figure 5.3). Figure 5.4 shows the phase portrait we
plotted by using our computer program. (The computer code and specific
numerical examples are in Appendix A). It is easier to see the stable point if
Figure 5.4: This phase portrait of the effective model shows trajectories that start at
angles between 0 and , which is the valid region for our kite model. The fixed points
occur at 0.44733 radians and 0.81 radians with (3 = 0.607432.
we zoom in on that point as shown in Figure 5.5. Notice that the trajectories
come together into a fine that is attracted to the fixed point. The fixed point,
9 = 0.81 in this particular example, is a stable node. Figure 5.5 has five
40
satellite graphs that show the behavior of the angle 6 with respect to time,
from 0 to 90 seconds. In other words, if the kite were placed at
(6,0) = (0.805,0.001), it would move up to 0.81 radians and stay there. If the
kite were placed at (6,6) = (0.815,0.001) radians, it would rise slightly, then fall
down to 0.81 radians then stay there. Placing the kite exactly on the stable
node, (6,6) = (0.81,0) means that the kite does not move at all.
I 0.81
<9
Time
Figure 5.5: The middle graph is a zoomed in view of the stable node. The upper left
picture is an angle versus time plot of the trajectory starting at the upper left corner
of the phase portrait (6,6) = (0.805,0.001) The bottom middle graph shows what
happens when the kite is placed precisely on the stable node. The kite does not move
from that position. The fixed point occurs at 0.81 radians with /? = 0.607432.
41
A calculation of the eigenvalues and eigenvectors of the Jacobian matrix
evaluated at the fixed point, verify that it is indeed a stable node and not a
stable spiral. Both eigenvalues are real and negative indicating that all the
solution trajectories decay into the fixed point. If the fixed point would have
been a stable spiral, the eigenvalues would have had imaginary parts. (Appendix
A contains the eigenvalues and eigenvectors for all the fixed points.)
The second aspect of the bifurcation diagram in Figure 5.3 is when
6 < 0*. The linearization indicates that the fixed point (9s,rjs) is a saddle point.
Figure 5.6 shows the phase portrait close to the saddle point.
This saddle point is very interesting. Trajectories that start out with a
positive velocity are drawn toward the stable node on the right. Trajectories
starting with negative velocity are drawn toward a stable node not in the valid
physical boundaries of the system. A saddle point however, has trajectories that
are drawn into the fixed point before turning away. The computer generated
phase portrait in Figure 5.6 does not show any trajectories being drawn toward
the fixed point. By finding the eigenvectors of the Jacobian matrix, we discover
that the eigenvector corresponding to the decaying eigensolution is (1,0.0094),
which is very close to the saxis. A clearer picture of the saddle point is drawn
in Figure 5.7.
Recall that linearization indicates behavior that is only local to fixed
points. We can see that trajectories leaving the location of the saddle point in
the positive 6 direction, do so along the growing eigenvector. Then they are
drawn toward the stable node on the right. If the kite were place on one of these
trajectories, it would eventually find the angle that is characterized by the stable
node. The trajectories that leave the location of the saddle point in the negative
6 direction, also do so along the growing eigenvector. These trajectories are
drawn toward a stable node, that has a negative angle, which is out of the valid
42
physical region of our system. If the kite were placed on one of these
trajectories, it would hit the ground.
Figure 5.6: The phase portrait drawn here include the saddle point. Trajectories
that start out with a positive velocity are drawn toward the stable node on the right.
Trajectories starting with negative velocity are drawn toward a stable node not in the
valid physical boundaries of the system. The three angle versus time plots show what
happens when the kite is placed at the three indicated points. The fixed point occurs
at 0.44733 radians with j3 = 0.607432.
43
0
Figure 5.7: This drawing shows the saddle point on the left and the stable node on
the right. Trajectories leaving the saddle point along the growing eigenvector soon
are drawn into the decaying eigenvector of the stable node.
The third, region of interest in our bifurcation diagram is when
f3b < P < ft*. The linearization indicates that the fixed point (9S, rjs) is an
unstable node or unstable spiral. Figure 5.8 seems to indicate that the point is
an unstable node. Again, this can be verified by calculating the eigenvalues.
They are both real and positive which indicates the solutions grow and hence are
repelled from the fixed point.
The fourth and last item of interest in the bifurcation diagram is the
point (/3b,6b). The fixed point changes from a stable node to an unstable node as
/? increases. The point at which this occurs is the bifurcation point.
The Kites and Bifurcation Theory paper, [Adomaitis], suggested that
this bifurcation point was a Hopf bifurcation. A Hopf bifurcation occurs when
the real parts of two complex conjugate eigenvalues change from negative to
positive or positive to negative. All of the eigenvalues in this model were real
numbers. So this point {f3b,0b) is not a Hopf bifurcation.
44
Figure 5.8: This phase portrait "zooms" in on the unstable node. A trajectory that
starts in the lower right corner goes down to zero, which physically means that the
kite hits the ground. Trajectories that start at the fixed point as well as trajectories
that start out in the upper right corner also go down to zero, just not as fast. The
fixed point occurs at 0.65 radians with (3 = 0.679557.
45
6. Conclusion
The basic model and the effective model have very different qualitative
behavior. The linearization of the basic model indicated a saddle point, which
we found in the phase portrait. If the system had been linear, the linearization
would have conclusively determined the existence of a center. However, since the
system was nonlinear, we could not use the results from the linearization on that
fixed point because of the sensitive nature of both eigenvalues on the imaginary
axis. Much to our surprise though, the phase portrait did have a center. The
basic model was interesting to study, however, it is clearly not representative of
a kite in flight.
The effective model, on the other hand, seems to model the kite more
closely. The system of ODEs that represent the effective model is much more
difficult to analyze. The lower leg of the bifurcation diagram indicated a saddle
point, just as in the basic model. The upper leg however, was split into two
parts: one corresponding to a stable point and the other an unstable point. We
then calculated the eigenvalues which verified the stable point was a stable node
and the unstable point was an unstable node. If they have both been spirals
instead of nodes, the bifurcation point would have been a Hopf bifurcation.
46
Appendix A. Calculated Values and C++ Code
This appendix contains the specific numbers as well as the C++ code
used to plot the trajectories. The constants used were
4
a = 3, b 4, and a = ,
o
p = 1.29 kg/m3,
g 9.81 m/sec2,
A = 2.5 m2,
r = 100 m,
and
m ~ 0.5 kg.
The bifurcation point of the effective model is
(/3b,6b) = (0.63945,0.76537).
The turning point of the effective model is
(/T,0*) = (0.6803,0.6316).
47
Phase Portrait Figure Number P A e
Basic Model 3.11 0.02069 7.0 0.0121 and 1.1775
Effective Model Big Picture 5.4 0.607432 1.291993 0.44733 and 0.81
Effective Model Stable Node 5.5 0.607432 1.291993 0.81
Effective Model Saddle Point 5.6 0.607432 1.291993 0.44733
Effective Model Unstable Node 5.8 0.679557 1.2215 0.65
Table A.l: This table summarizes the specific fixed used in the phase portraits.
48
Model Fixed Points Eigenvalues Corresponding Eigenvectors
Basic Saddle Point 2.8102 (0.3353,0.9421)
2.8102 (0.3353,0.9421)
Center 0 + 2.0288* (0.4421,0 + 0.8970*)
0 2.0288* (0.4221,0 0.8970*)
Effective Stable Node 0.0604 (0.9982, 0.0603)
1.5131 (0.5514,0.8343)
Saddle Point 0.0094 (1.0000,0.0094)
11.9396 (0.0835, 0.9965)
Unstable Node 0.0033 (1.0000,0.0033)
4.3140 (0.2258, 0.9742)
Table A.2: This table summarizes the eigenvalues and eigenvectors for each fixed
point calculated from the Jacobian matrices.
This is the C++ code that implements the RungeKutta method for the
effective model.
///////////////////////////////////////////////////////////
//
// EFFECTRK.C
//
// This utility draws phase portraits for the effective
// model described in this paper. It uses a fourth order
// RungeKutta method.
//
// This utility requires a PC with video capable of
// 640 x 480 at two colors. When it runs it sets the
// video mode for graphing and resets video mode back
// to normal when finished.
//
// This utility was developed using Microsoft Visual
// C++ Ver 1.52.
//
49
#include
#inelude
#inelude
#include
#inelude
ctime.h>
#define GRAPH //set/unset for debugging purposes
#ifdef GRAPH
#include #endif
#define ITERATIONS 10000L
#define METHOD "RungeKutta"
enum boolean { FALSE, TRUE };
// Declare global variables
long il, iterations;
double rho, bigA, a, b, lambdaE, r, m, thetaE, g, lambda,
theta, eta, myx, serratio, pi, wleft, wright, wtop,
wbottom, tb, lr, tbadj, lradj, dx, dy, xdiv, ydiv,
i, j, h, xx, yy, tmpO, tmpl, tmp2, tmp3, tmp4, If,
bt, rt, tp, vl, v2, v3, v4, wl, w2, w3, w4;
// place grid on working area
static void grid(void)
{
#ifdef GRAPH
for(i=wleft; i<=wright; i += lr / xdiv)
for(j=wbottom; j<=wtop; j += tb / ydiv)
_setpixel_w(i, j);
>
}
#endif
>
// helper function for drawing lines
50
static void line(double xl, double yl, double x2, double y2)
{
#ifdef GRAPH
_moveto_w(xl, yl);
_lineto_w(x2, y2);
#endif
}
// draw working area graph boundaries
static void liner(void)
{
#ifdef GRAPH
struct tm *newtime;
time_t aclock;
char buf[100];
time( feaclock ); // Get time in seconds
newtime = localtime( ftaclock ); // Convert to struct tm
_settextposition(2,2);
sprintfCbuf, "Left = %lf, Right = %lf",
wleft, wright);
_outtext(buf);
_settextposition(3,2);
sprintf(buf, "Top = /0lf, Bottom = %lf",
wtop, wbottom);
_outtext(buf) ;
_settextposition(4,2);
sprintf(buf, "Iterations = %ld, Method = %s, \
Scrratio = 70lf", iterations, METHOD, scrratio);
_outtext(buf);
_settextposition(5,2);
sprintf(buf, "Lambda = /0lf, Time = %s",
lambda, asctime( newtime ));
_outtext(buf);
line(wleft, wbottom, wleft, wtop);
line(wleft, wbottom, wright, wbottom);
line(wright, wtop, wleft, wtop);
line(wright, wtop, wright, wbottom);
51
line(0, wbottom, 0, wtop);
line(wleft, 0, wright, 0);
#endif
}
void main(int argc, char **argv)
{
// Declare local variables
char *penv;
short stopflag=0;
//Define good ol pi
pi = 3.14159265;
if(argc > 1) // If params were passed on command line
if(argc != 7)
{
printf(M\nUsage: effectrk \
\n");
exit(1);
>
else
{
// Define x and y working area ranges
wleft = atof (argv[l]);
wright = atof(argv[2]);
wtop = atof(argv[3]);
wbottom = atof(argv[4]);
iterations = atol(argv[5]);
lambda = atof(argv[6]);
>
else // Take defaults if params were not passed on cmd line
{
// Define x and y working area ranges
wleft = pi;
52
wright = pi;
wtop = pi;
wbottom = pi;
iterations = ITERATIONS;
lambda = 1.22;
}
#ifdef GRAPH
// Set screen video mode to 640 x 480 white on black
_setvideomode(_VRES2C0L0R);
_setcolor(7);
#endif
//Define other constants
rho = 1.29;
g = 9.81;
bigA
r
m
b
a
2.5
100
0.5
4.0
3.0
0;
// The default working area is .5 of full screen size. Give user
// ability to change it via an environment variable,
ifC(penv=getenv("SCRRATI0")) != NULL) scrratio=atof(penv);
else scrratio = .5; // working area/viewing area ratio
// Calculate working area screen size
tb = fabs(wtopwbottom);
lr = fabs(wrightwleft);
// Calculate actual screen size so that displayed image will be
// scrratio of full screen.
tbadj = fabs(wtopwbottom);
tbadj = tbadj/scrratio;
tbadj = (tbadj fabs(wtopwbottom))/2.0;
lradj = fabs(wrightwleft);
lradj = lradj/scrratio;
lradj = (lradj fabs(wrightwleft))/2.0;
53
// Specify # of horz and vert grid increments
xdiv = 10;
ydiv = 10;
// Add (subtract) scrratio of image/screen size difference\
// to working screen size area.
If = wleft lradj;
bt = wbottom tbadj;
rt = wright + lradj;
tp = wtop + tbadj;
// Set the viewing area screen size (its bigger then actual
// working area)
#ifdef GRAPH
// Set size of total viewing area. The actual area the
// graph occupies is scrratio of this size and will be
// displayed exactly in the middle of the screen.
_setwindow(TRUE, If, tp, rt, bt);
#endif
// Set the incremental step size
h = .001;
// Draw grid and graph boundaries
gridO ;
liner();
// Main processing loop
for(xx=wleft; xx<=(wright+(lr / xdiv)/2.0); xx += lr / xdiv)
{
for(yy=wbottom; yy<=(wtop+(tb / ydiv)/2.0); yy += tb / ydiv)
{
theta = xx;
eta = yy;
for(il =1; il <= iterations; ++il)
54
#ifdef GRAPH
// Check as efficiently as possible if pixel can be
// output to the working screen area.
if( theta > wright) {}
else if( theta < wleft) {>
else if( eta < wbottom) Q
else if( eta > wtop) Q
else _setpixel_w(theta, eta);
#endif
// The RungeKutta Method
// Calculate vl and wl
vl = h eta;
lambdaE = sqrt((lambda lambda) + ((r r) *
(eta eta)) + (2 r eta lambda *
sin(theta)));
tmpO = (rho bigA (lambdaE*lambdaE)) / (2 r m);
tmpl = tmpO a;
tmp2 = tmpO b;
myx = (lambda / lambdaE) cos(theta);
tmp3 = (sqrt(l myx myx)) / myx;
thetaE = atan(tmp3);
wl = (tmp2*sin(2*thetaE)*cos(thetaE)tmpl*cos(thetaE)*
sin(thetaE)(g/r)*cos(theta));
wl = h*wl;
// Calculate v2 and w2
v2 = h (eta + 0.5*wl);
lambdaE = sqrt((lambda lambda) + ((r r) *
((eta + 0.5*wl) (eta + 0.5*wl))) +
(2 r (eta + 0.5*wl) lambda *
sin(theta+0.5*vl)));
tmpO = (rho bigA (lambdaE*lambdaE)) / (2 r m);
tmpl = tmpO a;
tmp2 = tmpO b;
myx = (lambda / lambdaE) cos(theta+0.5*vl);
tmp3 = (sqrt(l myx myx)) / myx;
55
thetaE = atan(tmp3);
w2 = (tmp2*sin(2*thetaE)*cos (thetaE)tmpl*cos (thetaE)*
sin(thetaE)(g/r)*cos(theta+O.5*vl));
w2 = h*w2;
// Calculate v3 and w3
v3 = h (eta + 0.5*w2);
lambdaE = sqrt((lambda lambda) + ((r r) *
((eta + 0.5*w2) (eta + 0.5*w2))) +
(2 r (eta + 0.5*w2) lambda *
sin(theta+0.5*v2)));
tmpO = (rho bigA (lambdaE*lambdaE)) / (2 r m) ;
tmpl = tmpO a;
tmp2 = tmpO b;
myx = (lambda / lambdaE) cos(theta+O.5*v2);
tmp3 = (sqrt(l myx myx)) / myx;
thetaE = atan(tmp3);
w3 = (tmp2*sin(2*thetaE)*cos(thetaE)tmpl*cos(thetaE)*
sin(thetaE)(g/r)*cos(theta+O.5*v2));
w3 = h*w3;
// Calculate v4 and w4
v4 = h (eta + w3);
lambdaE = sqrt((lambda lambda) + ((r r)*
((eta + w3) (eta + w3))) + (2 r (eta + w3)*
lambda sin(theta+v3)));
tmpO = (rho bigA (lambdaE*lambdaE)) / (2 r m) ;
tmpl = tmpO a;
tmp2 = tmpO b;
myx = (lambda / lambdaE) cos(theta+v3);
tmp3 = (sqrt(l myx myx)) / myx;
thetaE = atan(tmp3);
w4 = (tmp2*sin(2*thetaE)*cos(thetaE)tmpl*cos(thetaE*
sin(thetaE)(g/r)*cos(theta+v3));
w4 = h*w4;
// Calculate new approximations
theta = theta + (1.0/6.0)*(vl + 2.0*v2 + 2.0*v3 + v4);
56
eta = eta + (1.0/6.0)*(wl + 2.0*w2 + 2.0*w3 + w4);
>
>
#ifdef GRAPH // allow user to interrupt run by capturing
// keyboard input
if(_kbhit())
C
_getchO ;
stopflag=l;
break;
}
#endif
>
#ifdef GRAPH
// Print diagnostic if keyboard interruption occurs
_settextposi.tion(33,1) ;
if(stopflag) _outtext("Interrupted");
else _outtext("Done");
while(!_kbhit());
// Reset video mode back to normal
_setvideomode( _DEFAULTM0DE );
#endif
>
57
REFERENCES
[Adomaitis] Adomaitis, Raymond A. Kites and Bifurcation Theory.
SIAM Review Volume 31, Number 3 (September 1989):478483.
[Belsky] Belsky, Nancy Ann. Building Kites. Flying High with Math. Dale
Seymour Publications, Palo Alto, California, 1995.
[Briggs] Briggs, William. Class notes, Continuous Modeling. University of
Colorado at Denver, Fall 1994.
[Geist] Geist, Marc A. Electrical Engineer, AT&T Bell Laboratories, 1996.
[Marchaj] Marchaj, C.A. AeroHydrodynamics of Sailing. International Marine
Publishing, Camden, Maine, 1988.
[Meriam] Meriam, J.L. Engineering Mechanics Volume 2 Dynamics. John Wiley
and Sons, New York, 1978.
[ODell] ODell, C.R. The Physics of Aerobatic Flight. Physics Today
(November 1987):2430.
[Press] Press, William H. Flannery, Brian P., Teukolsky, Saul A., Vetterling,
William T. Numerical Recipes in C. Cambridge University Press,
Cambridge, Massachusetts, 1988.
[Strogatz] Strogatz, Steven. Nonlinear Dynamics and Chaos. AddisonWesley
Publishing, Reading, Massachusetts, 1994.
[Tipler] Tipler, Paul A. Physics. Worth Publisher, Inc., New York, New York,
1976.
58
