demonstrated using the mathematical modeling package MLAB. A Poincar'e

return map for the chaotic case is also given.

An ideal pendulum (*i.e.* with no friction) will swing back
and forth (or loop in a full circle) forever if there is no outside force
other than gravity acting upon it. Moreover, a pendulum with friction will
come to rest if there is no other outside force besides gravity acting upon
it. A more general forced damped pendulum with a periodic driving force
pushing it shows more interesting asymptotic behavior than these two trivial
cases. The angular position in radians as a function of time **theta(t)** of
a forced damped pendulum is described by the following second order
differential equation.

**
**

theta''t + nu*theta't + sin(theta) = rho*sin(2*pi*f*t) theta(0) = theta0 theta't(0) = s

where **d^2*theta*dt^2** represents the inertia,
**nu*d(theta)/dt** represents friction at the pivot,
**sin(theta)** represents gravity, and **rho*sin(2*pi*ft)**
represents a sinusoidal frequency **f** driving torque applied at the pivot.
**theta0** is the initial angular position and **s** is the
initial angular velocity of the pendulum.

Numerical solutions show that both chaotic and periodic solutions of
the forced damped pendulum equation are possible depending on the
particular choice of system parameters **nu, rho** and **f**.

"pendulum.do --- Solve the ODE for a forced pendulum" "the differential equation for the forced damped pendulum" fct theta"t(t) = rho*sin(t) - c1*theta't - sin(theta) "rho*sin(t) is the driving force, c1*theta't is friction, sin(theta) is gravity" "Read-in parameter values" type "input the friction coefficient c1"; c1 = kread(); type "input the driven force amplitude" rho = kread(); type "input n: number of units of time to solve for"; n = kread(); "Read-in the initial conditions" type "Specify the initial conditions" type "input theta(0)"; theta0 = kread(); initial theta(0) = theta0 type "input theta'(0)"; theta1 = kread(); initial theta't(0) = theta1 "Specify the time values where we want the solution" data = 0:n:.1; "get an array from 0 to n at every .05" /* Now, integrate the differential equation. This is done implicitly inside the operator 'points', the result matrix is put in the array m. The first column of m is the data array generated above, the second column is the theta value at the time corresponding to the time value in the first column, the third column is the theta't value at the time corresponding to the time value in the first column. */ m = points(theta, theta't,data) "draw the trajectory (angular position vs. time)" draw m col 1;2, color red top title "Regular" left title "Cumulative Angular Position"; bottom title "Time" view; "view the picture"

The following pictures are generated by this **MLAB** *dofile*. (with
modified graphics). They are plots of time against the angular position. The
left picture is the trivial regular case with a periodic orbit where **c1 =
0** and **rho = 0**. The right picture shows a chaotic orbit where **c1 =
0.3** and **rho = 2.5**. In all cases, the initial point
corresponds to the initial conditions denoted by **theta0** and
**theta1** in the dofile.

Here are another pair of pictures that are generated by the same
*dofile*. The left picture shows a transient chaotic behavior which
settles down to a periodic orbit after awhile. The parameter values
are **c1 = 0.1, rho = 1.5**. The right picture shows another kind of
behavior. It is not periodic, since the angle keeps increasing, it is
not chaotic either since the the angle is increasing in a regular
pattern, thus the Lyapunov exponent is not positive. The parameter
values are **c1 = 0.2**, **rho = 2.5**.

To see a trajectory of a point in the phase-space
**(t, theta(t), theta'(t))** for the forced damped pendulum, we will
construct an object called the *Poincar'e* map which is often
used to reduce a continuous time system (or "flow") to a discrete
time map with one less dimension. The basic idea of a *Poincar'e* map is
to choose some appropriate hypersurface in the phase
space and observe the intersection of the orbit in the phase-space
with the surface. Since the solution of the system is unique with a
given initial point, when we neglect numerical round-off error, each
intersection point will uniquely determine the successive point.
Thus, a continuous flow is reduced to a discrete map with one less
dimension.

The picture above is the *Poincar'e* map of the
trajectory of the point **(0.5, 2.0)** in the **(theta(t), theta'(t))**
reduced-dimension phase space for our forced damped pendulum system
with the parameter **rho = 2.5** and **c1 = 0.3** which is the same as
the chaos picture above. We used **n = 10000**. It is the *Poincar'e*
return map with the section surface **t = 0**, which is equivalent to
**t = 2k*pi** for any integer **k** since the variable **t** only
appears in **sin(t)**, and for any integer value of **k**,
**sin(2k*pi)** will have the same value 0.

To compute the above *Poincar'e* map, we used almost the same
MLAB *dofile* with only a few minor changes. We only want
the solution at every time when **t = 2k*pi** (a multiple of **2*pi**),
Thus, we need to substitute the line for computing the time-list data
in the above MLAB *dofile* with the following line:

data = 0:(pi2*n):pi2; "get an array from 0 to pi2*n at every pi2"Also, we need to plot different columns of the solution matrix

/* The variable theta is the angle of the arm of the pendulum. The angle can range from 0 to pi2, any angle that is outside this region has a corresponding angle in this region. We define a function to map all the angles into 0:pi2 */ fct ft(theta) = mod(theta, pi2) m col 2 = ft on (m col 2); "map col 2 into the region 0 to pi2" "draw a dot at each phase plane at each point (theta,theta't)" draw m col 2:3, color red, lt none, pt dotpt top title "Forced Damped Pendulum" left title "Angular Velocity"; bottom title "Angle" view; "view the picture"

One can easily generate further pictures with other parameter values. More information about MLAB is available from:

Civilized Software Inc.

12109 Heritage Park Circle

Silver Spring, Md. 20906

Phone: (301)962-3711

FAX: (301)962-3712

Email: csi@civilized.com