\documentclass[twoside]{article}
\usepackage{epsf} % To input PostScript figures
\pagestyle{myheadings}
\markboth{ Multiple periodic solutions to a suspension bridge O.D.E. }
{ P. J. McKenna \& K. S. Moore }
\begin{document}
\setcounter{page}{183}
\title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent
Nonlinear Differential Equations, \newline
Electron. J. Diff. Eqns., Conf. 05, 2000, pp. 183--199\newline
http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.swt.edu or ejde.math.unt.edu (login: ftp)}
\vspace{\bigskipamount} \\
%
Multiple periodic solutions to a suspension bridge ordinary
differential equation
%
\thanks{ {\em Mathematics Subject Classifications:} 34C25, 34A47.
\hfil\break\indent
{\em Key words:} Torsional oscillations, suspension bridge.
\hfil\break\indent
\copyright 2000 Southwest Texas State University.
\hfil\break\indent Published October 25, 2000. } }
\date{}
\author{ P. J. McKenna \& K. S. Moore }
\maketitle
\begin{abstract}
We present an ordinary differential equation which models the
torsional motion of a horizontal cross section of a suspension
bridge. We use Leray-Schauder degree theory to prove that the
undamped equation has multiple periodic weak solutions. We use a
numerical continuation algorithm to demonstrate the existence of three
periodic solutions (one of small amplitude and two of large amplitude)
and to examine the bifurcation properties of the periodic solutions.
\end{abstract}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\section{Introduction} \label{intro}
In \cite{kn:MCK}, the first author considered a horizontal cross section of
the center span of a suspension bridge and proposed an ordinary
differential equation model for the torsional motion of the cross
section. Using physical constants from the engineers' reports of the
Tacoma Narrows collapse, he investigated this model
numerically. By specifying the initial position and velocity of the
cross section and using the Runge-Kutta method to solve the initial
value problem over long time, he demonstrated that under the same
small periodic forcing term, small or large amplitude periodic motion
may result; the ultimate outcome depends on the initial conditions.
The methodology in the above-mentioned paper was somewhat
primitive. Different initial conditions were prescribed randomly and the
eventual behavior of the solution of the initial value problem was
observed. Sometimes the motion converged to a large amplitude solution
and sometimes to the small near-equilibrium solution.
In this paper, we present a more systematic approach to the study of
the equation for the torsional motion of a
cross section of the center span. We use Leray-Schauder degree
theory to prove that, under certain
physical assumptions, the undamped equation has multiple periodic
solutions. We demonstrate numerically that for small
forcing, multiple periodic solutions exist and that whether large or
small amplitude motion results depends only on the initial conditions.
Finally, we use a more sophisticated approach to compute periodic
solutions of the nonlinear differential equation. Using continuation
methods, we examine the bifurcation properties of periodic solutions
as the amplitude of the forcing term varies and we demonstrate that
bifurcation from single to multiple periodic solutions occurs for
small forcing.
\section{The Equation for the Torsional Motion of the Cross Section}
We treat the center span of the bridge as a beam of length $L$ and
width $2l$ suspended by cables (see figure \ref{fig1}). To model
the motion of a horizontal cross section of the beam, we
treat it as a rod of length $2l$ and mass $m$ suspended by cables.
Let $y(t)$ denote the downward distance of the center of gravity of the rod from
the unloaded state and let $\theta (t)$ denote the angle of the
rod from horizontal at time $t$ (see figure \ref{fig2}).
\begin{figure}
\begin{center}
\epsfxsize=10cm
\epsffile{fig1.eps}
\end{center}
\caption{\label{fig1} A simple model of the center span}
\end{figure}
\begin{figure}
\begin{center}
\epsfxsize=10cm
\epsffile{fig2.eps}
\end{center}
\caption{\label{fig2} A horizontal cross section of the center span}
\end{figure}
We assume that the cables do not resist compression, but resist
elongation according to Hooke's Law with spring constant $K$; i.e.,
the force exerted by the cable is proportional to the elongation in
the cable with proportionality constant $K$. In figure \ref{fig2} we see
that the extension in the
right hand cable is $(y-l\sin{\theta})$, and hence the force exerted by the
right hand cable is
\[
-K (y-l\sin \theta)^+=\left \{
\begin{array}{ll}
-K(y-l \sin \theta) & \mbox{if $y - l\sin \theta \geq 0$}\\
0 & \mbox{if $y - l\sin \theta < 0$}
\end{array}
\right .
\]
where $u^+=\max \{u,0 \}$.
Similarly, the extension in the left hand cable is $(y+l\sin\theta)$ and
the force exerted by the left hand cable is
$-K(y+l\sin\theta)^+$.
In \cite{kn:MCK}, the author showed that the torsional and vertical
motion satisfy
\begin{eqnarray}
\label{eq:torode}
\theta''
&=&\frac{3K}{ml}\cos{\theta}[(y-l\sin{\theta})^{+}-(y+l\sin{\theta})^{+}]-\delta_1\theta'
+f(t) \\
y''&=& -\frac{K}{m}[(y-l\sin{\theta})^{+}+(y+l\sin{\theta})^{+}]-\delta_2y'+g\,,
\nonumber
\end{eqnarray}
where $\delta_1,\delta_2$ are damping constants, $g$ is the force due
to gravity, and $f(t)$ is the external force at time $t$.
Assuming that the cables never lose tension, we have $y \pm l
\sin \theta \geq 0$ and hence
$(y \pm l\sin{\theta})^{+}=(y \pm l\sin{\theta})$.
Thus, the equations $(\ref{eq:torode})$ become uncoupled and the torsional
and vertical motion satisfy
\begin{eqnarray}
\label{eq:torode1}
\theta'' &=& -\frac{6K}{m}\cos{\theta}\sin{\theta}-\delta_1 \theta'
+f(t)\\
\label{eq:vertode1}
y''&=&-\frac{2K}{m}y-\delta_2y' +g
\end{eqnarray}
respectively.
The equation for the vertical motion, $(\ref{eq:vertode1})$, is simply the equation for a damped, forced,
linear harmonic oscillator and the behavior of its solutions is well
known, \cite{kn:BDH}. We will study the equation for the torsional motion
$(\ref{eq:torode1})$, which is the
damped, forced, pendulum equation and is known to possess
chaotic solutions, \cite{kn:BDH}.
To choose the physical constants $K,m,\delta:= \delta_2$ and the
external forcing term $f(t)$, we rely on \cite{kn:FWA},
\cite{kn:MCK}, and \cite{kn:Scan}. From \cite{kn:FWA}, we choose $m=2500$ and
$\delta=.01$. To determine $K$, from
\cite{kn:FWA} we know that the main span would deflect about half
a meter when loaded with 100 kgs per unit length, so we have
$100(9.8) - 2K(.5)=0$ and we take $K=1000$.
For a cross section similar to the Tacoma
Narrows bridge, wind tunnel experiments indicate that aerodynamic forces
induce approximately sinusoidal oscillations of amplitude three
degrees \cite{kn:Scan}, so in $(\ref{eq:torode1})$
we choose $f(t) = \lambda \sin
(\mu t)$ where $\lambda \in [0,0.06]$ is chosen to produce the
appropriate behavior near equilibrium and the frequency $\mu$ is
chosen to match the frequency of the oscillations observed at Tacoma
Narrows on the day of the collapse. The frequency of the torsional
motion was approximately one cycle every 4 or 5 seconds, so we take
$\mu \in [1.2,1.6]$.
\section{A Multiplicity Theorem}
In this section, we prove that undamped equations of the form
$(\ref{eq:torode1})$ have multiple periodic weak solutions. To some
extent, this result is inspired by the earlier results of Castro and
Lazer, \cite{kn:Castro}, where they obtained multiple solutions of elliptic
boundary value problems once a critical value, measured by $f(0)$,
passes across key eigenvalues.
Let $\Omega = (- \pi, \pi)$, ${\mathcal{H}} = \{u \in
{\mathcal{L}}^2(\Omega)| u(t) = -u(-t)\}$ and define $Lu = -u''$. For
$u \in {\mathcal{H}}$, let $\|u\| = \|u\|_{{\mathcal{L}}^2(\Omega)} =
(\int_\Omega |u|^2 dt)^\frac{1}{2}$.
Using $\cos u \sin u = \frac{1}{2} \sin 2u$, removing the damping
term, and imposing periodicity
conditions, we rewrite $(\ref{eq:torode1})$ as
\begin{eqnarray}
\label{eq:undamp}
&Lu = b \sin 2u - \varepsilon h(t)&\\
&u(- \pi) = u(\pi), \quad u'(- \pi) = u'(\pi)&\nonumber
\end{eqnarray}
We observe that the eigenvalues and eigenfunctions of $L$ are
\begin{equation} \label{eq:evalues}
\lambda_n = n^2\,,\quad \phi_n(t) = \sin nt
\end{equation}
hence $L^{-1}$ exists, is compact, and $\|L^{-1}\| = 1$.
\paragraph{Definition}
We say that $u \in {\mathcal{H}}$ is a solution to $(\ref{eq:undamp})$
if
\begin{equation}
\label{eq:solndef}
u - L^{-1}(b\sin2u - \varepsilon h) =0.
\end{equation}
\medskip
We will prove the
existence of multiple solutions $u \in {\mathcal{H}}$ to the
problem $(\ref{eq:undamp})$.
\begin{theorem}
Let $h \in {\mathcal{H}}$ with $\|h\| \leq 1$ and let $b \in
(\frac{1}{2},2)$. Then there exists $\varepsilon_0 >0$ such that if
$|\varepsilon|< \varepsilon_0$, $(\ref{eq:undamp})$ has at least two solutions.
\end{theorem}
\paragraph{Proof}
Let $\varepsilon < \varepsilon_0$; we will determine the value of
$\varepsilon_0$ later. Note that by $(\ref{eq:evalues})$ and
our choice of $b$, $2b$ is {\em not} an eigenvalue of $L$.
Define $T_1:{\mathcal{H}} \rightarrow {\mathcal{H}}$ by
\begin{eqnarray*}
T_1(u)=u-L^{-1}(b \sin 2u - \varepsilon h)
\end{eqnarray*}
and note that, by $(\ref{eq:solndef})$, zeros of $T_1$ correspond to
solutions of $(\ref{eq:undamp})$. Denote the Leray-Schauder
degree of the map $T_1$ in the domain $U$ at the point $p$ by
$\deg(T_1,U,p)$. To prove the theorem, we will show
\begin{description}
\item{(D1):} there exists $R_0>0$ such that for $R > R_0$,
$\deg(T_1,B_R(0),0)=1$
\item{(D2):} there exists $\gamma \in (0,R_0)$ such that $\deg(T_1, B_\gamma
(0),0)=-1$
\end{description}
\noindent Then, since $\deg(T_1, B_\gamma(0),0) \neq 0$, there exists a zero
of $T_1$ (i.e., a solution of (\ref{eq:undamp})) in
$B_\gamma(0)$. Moreover, by the additivity property of degree, we have
$\deg(T_1, B_R(0) \backslash \overline{B_\gamma(0)},0) \neq 0$ and hence
$(\ref{eq:undamp})$ has second solution in the annulus $B_R(0) \backslash
\overline{B_\gamma(0)}$.
To establish (D1), define
$$
T_\beta u=u-\beta L^{-1}(b \sin2u - \varepsilon h)
$$
for $\beta \in [0,1]$ and note that
this definition of $T_1$ is consistent with our previous definition.
Note also that $T_0$ is simply the identity map, hence for any $R>0$
we have $\deg(T_0,B_R(0),0)=1$. The homotopy property of degree
ensures that $\deg(T_\beta,B_R(0),0)$ is constant provided that $0
\notin T_\beta (\partial B_R(0))$ for all $\beta \in [0,1]$.
Fix $\beta \in [0,1]$ and suppose $u \in {\mathcal{H}}$ solves $T_\beta u
=0$. We will show that $u$ is bounded above by some $R_0 > 0$ and that
this bound is independent of $\beta$.
Since $T_\beta u=0$, we have
\begin{eqnarray*}
\|u\|&=&\beta \| L^{-1}(b \sin 2u - \varepsilon h)\|\\
&\leq& \beta [\varepsilon_0 + b \|\sin 2u\|]\\
&\leq& \beta [\varepsilon_0 +b m(\Omega)^{1/2}]\\
&\leq& [\varepsilon_0 +b m(\Omega)^{1/2}]\\
&=& [\varepsilon_0 +b \sqrt{2 \pi}] < R_0
\end{eqnarray*}
if we choose $R_0 > \varepsilon_0 +b \sqrt{2 \pi}$.
Thus, for $R> R_0$ we have
\begin{equation}
\deg(T_1,B_R(0),0)=\deg(T_0,B_R(0),0)=1
\end{equation}
and (D1) above holds.
To establish (D2), define
\begin{eqnarray*}
T_\mu u=u - (1- \mu)L^{-1}(2bu) -\mu L^{-1}(b \sin 2u- \varepsilon h)
\end{eqnarray*}
for $\mu \in [0,1]$
and note again that this definition of $T_1$ is consistent with our
previous definitions. We will again apply the homotopy property of
degree and a standard degree calculation to show that for some $\gamma >0$
\begin{equation}
\label{eq:smallball}
\deg(T_1,B_\gamma(0),0)=\deg(T_0,B_\gamma(0),0)=-1.
\end{equation}
We will show first that for any $\gamma > 0$,
\begin{equation}
\label{eq:minus1}
\deg(T_0, B_\gamma(0),0) = deg (I - 2bL^{-1},B_\gamma(0),0) = -1.
\end{equation}
Consider the finite dimensional subspace ${\mathcal{M}}_N$
:=span$\{\phi_{n}\}_{n=1}^N$ of ${\mathcal{H}}$ and recall that, by
compactness, $2bL^{-1}$ can be approximated in operator norm by the operators
$B_N:{\mathcal{M}}_N \rightarrow {\mathcal{M}}_N$ given by
\[
B_N(u)=2b \sum_{n=1}^N \frac{c_{n}}{\lambda_{n}} \phi_{n}.
\]
By definition, for $N$ sufficiently large,
\begin{eqnarray*}
\deg(T_0,B_\gamma(0),0)&=&\deg(I-B_N,B_\gamma(0) \cap {\mathcal{M}}_N,0)\\
&=& \sum_{u \in (I-B_N)^{-1}(0)} \mathop{\rm sign} J_{I-B_N}(u)
\end{eqnarray*}
where $J_\phi(u)$ is the Jacobian determinant of $\phi$ at $u$.
Since $I-B_N$ can be identified with an $N \times N$ diagonal matrix whose entries are $1-\frac{2b}{\lambda_{n}}$, we have
\[
\deg(I-B_N,B_\gamma(0) \cap {\mathcal{M}}_N,0)=\mathop{\rm sign}
\prod_{n=1}^N (1-\frac{2b}{\lambda_{n}}).
\]
Since $b \in (\frac{1}{2},2)$ and $\lambda_n=n^2$, the only negative value of
$1-\frac{2b}{\lambda_{n}}$ occurs at $\lambda_1=1$. Since the
eigenvalues of $L$ are simple, we have
\[
\deg(I-B_N,B_\gamma(0) \cap {\mathcal{M}}_N,0)=-1
\]
and $(\ref{eq:minus1})$ holds.
To establish $(\ref{eq:smallball})$, and hence (D2), we must show that
for some $\gamma >0$,
\begin{equation}
\deg(T_1, B_\gamma(0), 0) = \deg(T_0, B_\gamma(0), 0).
\end{equation}
The homotopy property of degree ensures that $\deg(T_\mu, B_\gamma(0),
0)$ is constant provided $0 \notin T_\mu(\partial B_\gamma(0))$
for all $\mu \in [0,1]$. Observe that $u=0$ is the only zero of $T_0$
since $2b$ is not an eigenvalue of $L$. Fix $\mu \in (0,1]$ and
suppose that $u \neq 0$ solves $T_\mu u =0$. Set $\|u\| =
\tilde{\gamma}_\mu$; we will show that $\tilde{\gamma}_\mu$ is bounded
below by some $\gamma_\mu >0$.
Set $\psi = \frac{u}{\|u\|}=\frac{u}{\tilde{\gamma}_\mu}$. We claim
first that there exists a compact set $K$ with $\psi \in K$.
Since $u$ is a zero of $T_\mu$, $u$ satisfies
\begin{equation}
\label{eq:punchline1}
Lu - 2bu = \mu(b\sin 2u - 2bu -\varepsilon h)
\end{equation}
and therefore
\[
Lu = \mu(b\sin 2u -\varepsilon h) + (1 - \mu)2bu.
\]
Thus,
\begin{eqnarray*}
\|Lu\| &\leq& \|b \sin 2u - \varepsilon h \| + 2b \|u\|\\
&\leq& b \|\sin 2u\| + 2b\tilde{\gamma}_\mu +\varepsilon_0\\
&\leq& b m(\Omega)^\frac{1}{2} + 2b \tilde{\gamma}_\mu + \varepsilon_0
\leq \sqrt{2 \pi} b+2b\tilde{\gamma}_\mu +\varepsilon_0.
\end{eqnarray*}
It follows that $u \in \overline{
L^{-1}(B_{\sqrt{2\pi}b+2b\tilde{\gamma}_\mu+ \varepsilon_0}(0))}$,
which is compact
since $L^{-1}$ is compact, and thus
there exists a compact set $K$ with $\psi \in K$.
Since $u$ satisfies $(\ref{eq:punchline1})$, we have
\begin{equation}
\label{eq:punchline2}
\|u - L^{-1}2bu\| = \mu \|L^{-1}(b \sin 2u -2bu - \varepsilon h)\|.
\end{equation}
Denote the left and right hand sides of $(\ref{eq:punchline2})$ by
$LHS$ and $RHS$ respectively. Since $2b$ is not an eigenvalue of $L$,
we have $L\psi - 2b \psi \neq 0$ and hence
\[
\inf_{\psi \in K}\|\psi - L^{-1}2b \psi\|= \alpha >0
\]
and for our $u$ we have
\begin{equation}
\label{eq:punchline3}
LHS = \|u - L^{-1}2bu\| \geq \alpha \tilde{\gamma}_\mu.
\end{equation}
Now considering the right hand side of $(\ref{eq:punchline2})$, we
have
\begin{eqnarray*}
RHS &\leq& \mu \|b \sin 2u - 2bu - \varepsilon h\|\\
&\leq& \mu [\varepsilon_0 + \|b\sin2u - 2bu\|]\\
&=&\mu[\varepsilon_0 + \|b\sin (2\tilde{\gamma}_\mu \psi) - 2b
\tilde{\gamma}_\mu \psi\|].
\end{eqnarray*}
We claim that if $\tilde{\gamma}_\mu$ is sufficiently small, $RHS <
\alpha \tilde{\gamma}_\mu$, which contradicts
$(\ref{eq:punchline3})$. To establish this, we must first prove the
following lemma.
\begin{lemma} Let $K, {\mathcal{H}}, b$ be as above and denote
$\tilde{\gamma}_\mu$ by $\eta$.
Then there exists a function $\delta:(0,\infty) \rightarrow (0,\infty)$
such that
\begin{description}
\item{(L1):} $\|b\sin (2 \eta \psi) - 2b \eta \psi\| \leq \eta \delta (\eta)$
\item{(L2):} $\delta(\eta) \rightarrow 0$ as $\eta \rightarrow 0$ \\
\end{description}
\noindent
hold for all $\psi \in K, \eta > 0$
\end{lemma}
\paragraph{Proof}
Define $\delta:(0,\infty) \rightarrow (0,\infty)$ by
\[
\delta(\eta) = \max_{\psi \in K} \|\frac{b}{\eta} \sin 2\eta \psi -
2b\psi\|
\]
and note that (L1) above is satisfied.
To show that (L2) holds, we must show that $f_\eta(\psi):=
\|\frac{b}{\eta} \sin 2
\eta \psi - 2b\psi\| \rightarrow 0$ uniformly on $K$ as $\eta
\rightarrow 0$.
To show that $f_\eta \rightarrow 0$ uniformly on $K$, we'll show that
$f_\eta(\psi) \rightarrow 0$ for each $\psi \in K$ and that the family
${\mathcal{F}}:=\{f_\eta\}$ is equicontinuous on $K$.
Choose $\psi \in K$, $\psi \neq 0$ and note that
\[\frac{b}{\eta}\sin 2 \eta \psi - 2b \psi = b \frac{\sin 2 \eta
\psi}{2 \eta \psi}\cdot 2 \psi - 2b \psi \rightarrow 0
\]
as $\eta \rightarrow 0$, hence
\[
|\frac{b}{\eta} \sin 2 \eta \psi - 2b \psi|^2 \rightarrow 0.
\]
Moreover, since $|\sin w| \leq |w|$ and since $w^2$ is convex, we have
\begin{eqnarray*}
|\frac{b}{\eta}\sin 2 \eta \psi - 2b \psi|^2 &=&
4|\frac{1}{2}\frac{b}{\eta}\sin 2 \eta \psi +
\frac{1}{2}(-2b\psi)|^2\\
&\leq&4[\frac{1}{2}|\frac{b}{\eta}\sin 2 \eta
\psi|^2+\frac{1}{2}|2b\psi|^2]\\
&\leq& 2[\frac{b^2}{\eta^2}|2\eta \psi|^2 +|2b\psi|^2] = 16 b^2 \psi ^2
\in {\mathcal{L}}^1
\end{eqnarray*}
since $\psi \in {\mathcal{L}}^2$. By the Dominated Convergence
Theorem, we conclude that $f_\eta(\psi) \rightarrow 0$ for each $\psi
\in K$, as desired.
To see that the family ${\mathcal{F}} = \{f_\eta\}$ is
equicontinuous on $K$, choose $\tilde{\varepsilon} > 0$ and $\psi, \tilde{\psi}
\in K$. We have that
\begin{eqnarray*}
|f_\eta(\psi)-f_\eta(\tilde{\psi})| &=& | \|\frac{b}{\eta}\sin 2 \eta
\psi -2b \psi\| - \|\frac{b}{\eta} \sin 2 \eta \tilde{\psi} - 2b
\tilde{\psi}\| |\\
&\leq& \|\frac{b}{\eta}(\sin 2 \eta \psi - \sin 2 \eta \tilde{\psi}) -
2b(\psi - \tilde{\psi})\|\\
&\leq& \frac{b}{\eta}\|\sin 2 \eta \psi - 2 \eta \tilde{\psi}\| + 2b
\|\psi - \tilde{\psi}\|\\
&\leq&\frac{b}{\eta}\|2 \eta \psi - 2 \eta \tilde{\psi}\| + 2b\|\psi -
\tilde{\psi}\|\\
&=&4b\|\psi - \tilde{\psi}\| < \tilde{\varepsilon}
\end{eqnarray*}
provided $\|\psi - \tilde{\psi}\| < \frac{\tilde{\varepsilon}}{4b}$.
Since ${f_\eta}$ are equicontinuous on $K$ and converge pointwise on
$K$, we have that $f_\eta$ converge uniformly on $K$ and hence (L2) holds.
\hfill$\diamondsuit$ \smallskip
Returning now to the proof of the theorem and invoking the lemma
above, we have that
\[
RHS \leq \mu[\varepsilon_0 + \tilde{\gamma}_\mu
\delta(\tilde{\gamma}_\mu)].
\]
Assume now that $\varepsilon_0 < \frac{1}{2}\alpha
\tilde{\gamma}_\mu$. Since $\delta \rightarrow 0$, there exists
$\gamma_\mu >0$ such that if $\tilde{\gamma}_\mu < \gamma_\mu$, then
$\delta(\tilde{\gamma}_\mu) < \frac{1}{2} \alpha$ and hence
\[
RHS < \alpha \tilde{\gamma}_\mu.
\]
But this contradicts $(\ref{eq:punchline3})$. Therefore
$\tilde{\gamma}_\mu$ is bounded below by $\gamma_\mu>0$. Note also
that because of the factor of $\mu$ on the right hand side of the
above inequalities, we can take $0 < \gamma_1 \leq \gamma_\mu$ for
$\mu \in (0,1]$, and
hence we can choose $\gamma \in (0, \gamma_1)$.
\hfill$\diamondsuit$
\section{The Bifurcation Curve of Periodic Solutions}
As promised in the introduction, we now begin a more systematic
numerical study of the structure of the solution set of periodic
solutions of the torsional equation in the range of parameter values
where the large
amplitude periodic solutions were observed at Tacoma Narrows. Instead
of randomly varying the initial conditions and hoping that
multiple periodic solutions show up as long-term behavior of the
solution, we will study the structure of the solutions as we
continuously vary the amplitude of the forcing term.
Recall that the equation for the torsional motion of a horizontal
cross section of the center span of a suspension bridge is given by
$(\ref{eq:torode1})$. Using the physical constants from the
engineers' reports of the Tacoma Narrows failure and the forcing term
described in section \ref{intro} yields the ordinary differential equation
\begin{equation}
\label{eq:torode2}
\theta''=-2.4 \cos \theta \sin \theta - .01 \theta' + \lambda \sin \mu
t.
\end{equation}
In \cite{kn:MCK}, the author prescribed the
initial position $\theta(0)$ and velocity $\theta'(0)$ of the cross
section and employed the Runge-Kutta method to solve the initial value problem
over large time. He demonstrated numerically that small or large
amplitude periodic motion may result, depending only on the initial
conditions. Moreover, the amplitude and frequency of the large
amplitude solutions matched the behavior observed at Tacoma Narrows on
the day of its collapse, \cite{kn:FWA}.
As in \cite{kn:MCK}, we compute periodic solutions to
$(\ref{eq:torode2})$. Moreover, we describe a numerical
continuation algorithm and we use it to plot the
amplitude of periodic solutions versus the amplitude $\lambda$ of the
external forcing term. We demonstrate that for small $\lambda$,
three periodic solutions to $(\ref{eq:torode2})$ exist and that
whether large or small amplitude motion results depends on the initial
position and velocity of the cross section. Moreover, we demonstrate
that bifurcation from single to multiple periodic solutions occurs for
small $\lambda$.
\subsection{The Algorithm} \label{contalg}
We wish to generate the bifurcation curve for periodic solutions to
$(\ref{eq:torode2})$ as the amplitude $\lambda$ of the external
forcing term varies; i.e., we wish to plot the amplitude $A_\lambda$
of the periodic solution versus $\lambda$. Define $T=2\pi/\mu$.
To compute $T$ periodic solutions to
$(\ref{eq:torode2})$, we rewrite it as a first order system as
follows. Let
\begin{equation}
U= \left [ \begin{array}{c}
u\\
v
\end{array} \right ] = \left [
\begin{array}{c}
\theta\\
\theta' \end{array} \right ].
\end{equation}
Then solutions $\theta$ of $(\ref{eq:torode2})$ correspond to
solutions $U$ of the first order system
\begin{eqnarray}
\label{eq:system4}
&u' = v&\\
&v' = -2.4 \cos u \sin u -.01 v + \lambda \sin \mu t &\nonumber
\end{eqnarray}
We employ a shooting algorithm to
find periodic solutions to the above system. In the shooting
algorithm, we search for initial conditions
$\eta = [\eta_1,\eta_2]^T$ so that the solution to the initial value
problem
\begin{eqnarray}
\label{eq:sysivp4}
& u' = v &\nonumber\\
& v' = -2.4 \cos u \sin u -.01 v + \lambda \sin \mu t &\\
& u(0)= \eta_1, \quad v(0)= \eta_2 \nonumber
\end{eqnarray}
is $T$-periodic, i.e., we search for initial conditions
$\eta$
so that the solution $U$ to $(\ref{eq:sysivp4})$ satisfies
\[
U(T)=U(0)=\eta.
\]
Define $F:{\mathbf{R}}^2 \rightarrow {\mathbf{R}}^2$ by
\begin{equation}
\label{eq:fdef}
F(\eta)=\eta-U(T).
\end{equation}
Then searching for initial conditions $\eta$ which yield a
$T$-periodic solution to $(\ref{eq:sysivp4})$ is equivalent to
searching for zeros of $F$; we employ Newton's method
to search for these zeros. Specifically, we take a ``good guess'' at
an initial condition $\eta^0$ which yield a
$T$-periodic solution to $(\ref{eq:sysivp4})$. We then solve the
initial value problem $(\ref{eq:sysivp4})$ for $t \in [0,T]$ and test whether
$F(\eta^0) = 0$. If not, we update
$\eta^0$ via Newton's method:
\[
\eta^1=\eta^0 - DF(\eta^0)^{-1}F(\eta^0)
\]
where
\begin{equation}
\label{eq:jacdef4}
DF=
\left[ \begin{array}{cc}
\frac{\partial F_1}{\partial \eta_1} & \frac{\partial F_1}{\partial \eta_2}\\
\frac{\partial F_2}{\partial \eta_1} & \frac{\partial F_2}{\partial
\eta_2}
\end{array} \right ]
\end{equation}
is the Jacobian matrix of $F$.
To compute this matrix, observe that by our definition of $F$ in $(\ref{eq:fdef})$, we have
\begin{equation}
\label{eq:jacobian4}
DF=
\left[ \begin{array}{cc}
\frac{\partial F_1}{\partial \eta_1} & \frac{\partial F_1}{\partial \eta_2}\\
\frac{\partial F_2}{\partial \eta_1} & \frac{\partial F_2}{\partial
\eta_2}
\end{array} \right ] =
\left [ \begin{array}{cc}
1-\frac{\partial u}{\partial \eta_1}(T) &
-\frac{\partial u}{\partial \eta_2}(T) \\
-\frac{\partial v}{\partial \eta_1}(T) &
1-\frac{\partial v}{\partial \eta_2}(T)
\end{array} \right ].
\end{equation}
To compute $\frac{\partial u}{\partial \eta_i},\frac{\partial
v}{\partial \eta_i}, i=1,2$, we recall that the solution to the
initial value problem is continuously differentiable with respect to initial
conditions and parameters, \cite{kn:Codd}, and hence we can
differentiate $(\ref{eq:sysivp4})$ with respect to $\eta_1$.
This yields the first order system in $\frac{\partial u}{\eta_1}, \frac{\partial v}{\eta_1}$
\begin{equation}
\left (
\begin{array}{c}
\frac{\partial u}{\partial \eta_1}\\
\frac{\partial v}{\partial \eta_1}
\end{array} \right )' = \left ( \begin{array}{c}
\frac{\partial v}{\partial \eta_1}\\
-2.4 \cos (2u) \frac{\partial u}{\partial \eta_1} - .01 \frac{\partial
v}{\partial \eta_1}
\end{array} \right );\\
\left ( \begin{array}{l}
\frac{\partial u}{\partial \eta_1}\\
\frac{\partial v}{\partial \eta_1}
\end{array} \right ) (0) =
\left ( \begin{array}{l}
1\\ 0
\end{array} \right )
\end{equation}
Hence $[\frac{\partial u}{\partial \eta_1},\frac{\partial v}{\partial
\eta_1}]^T$ solves the initial value problem
\begin{equation}
\label{eq:jacsys1}
\left ( \begin{array}{c}
w\\
z \end{array} \right ) ' =
\left ( \begin{array}{c}
z\\
-2.4 \cos (2u) w - .01 z
\end{array} \right );\\
\left ( \begin{array}{l}
w\\
z \end{array} \right ) (0) =
\left ( \begin{array}{l}
1\\ 0
\end{array} \right )
\end{equation}
where $u$ is the solution to the initial value problem
$(\ref{eq:sysivp4})$.
Similarly, differentiating $(\ref{eq:sysivp4})$ with respect to
$\eta_2$, we see that
$[\frac{\partial u}{\partial \eta_2},\frac{\partial v}{\partial
\eta_2}]^T$ solves the initial value problem
\begin{equation}
\label{eq:jacsys2}
\left ( \begin{array}{c}
w\\
z \end{array} \right ) ' =
\left ( \begin{array}{c}
z\\
-2.4 \cos (2u) w - .01 z
\end{array} \right );\\
\left ( \begin{array}{l}
w\\
z \end{array} \right ) (0) =
\left ( \begin{array}{l}
0\\ 1
\end{array} \right ).
\end{equation}
Therefore, to compute the entries of the Jacobian matrix
$(\ref{eq:jacobian4})$, we must solve the initial value problems
$(\ref{eq:sysivp4}), (\ref{eq:jacsys1})$, and $(\ref{eq:jacsys2})$.
The above algorithm (Shooting Algorithm with Newton Update)
suggests a method by which we can plot the
amplitude $A_\lambda$ of periodic solutions to $(\ref{eq:torode2})$
versus the amplitude $\lambda$ of the external forcing term.
Specifically, fix $\lambda = \hat{\lambda}$, use the algorithm
described above to compute the initial conditions which yield a $T$
periodic solution to $(\ref{eq:sysivp4})$, record the amplitude
$A_{\hat{\lambda}}$ of the periodic solution, increment $\lambda$,
and repeat.
Unfortunately, this algorithm fails at bifurcation points because,
under natural parameterization, the Jacobian matrix
$(\ref{eq:jacobian4})$ is singular at such points. However, we can
remedy this difficulty via pseudoarclength parameterization. We
describe the parameterization briefly below and refer the reader to
\cite{kn:Jen}, \cite{kn:Keller} for details.
Let $(\eta^0,\lambda^0)$ be a zero of $F$ and set
$(\eta^0,\lambda^0)=(\eta(s_0),\lambda(s_0))$. We will introduce the
pseudoarclength normalization $N$ given by
\begin{equation}
\label{eq:ndef}
N(\eta,\lambda,s)=\psi\frac{\|\eta(s)-\eta(s_0)\|^2}{s-s_0}+(1-\psi)\frac{|\lambda(s)-\lambda(s_0)|^2}{s-s_0}
- (s-s_0)
\end{equation}
where $\psi \in (0,1)$ and $s$ is a chord-length parameter, \cite{kn:Jen}.
We use an algorithm similar to the Shooting Algorithm with Newton
Update described above to search for zeros
of the map ${\mathcal{F}}$ given by
\begin{equation}
\label{eq:scriptfdef}
{\mathcal{F}}(\eta,\lambda,s) =
\left( \begin{array}{cc}
F(\eta,\lambda)\\
N(\eta,\lambda,s) \end{array}
\right).
\end{equation}
The Jacobian matrix
\begin{equation}
\label{eq:jacn4}
D{\mathcal{F}} = \left [
\begin{array}{ccc}
\frac{\partial F_1}{\partial \eta_1} & \frac{\partial F_1}{\partial
\eta_2} & \frac{\partial F_1}{\partial \lambda}\\
\frac{\partial F_2}{\partial \eta_1} & \frac{\partial F_2}{\partial
\eta_2} & \frac{\partial F_2}{\partial \lambda}\\
\frac{\partial N}{\partial \eta_1} & \frac{\partial N}{\partial
\eta_2} & \frac{\partial N}{\partial \lambda}
\end{array} \right ]
\end{equation}
is nonsingular even if the matrix $DF$ is singular, \cite{kn:Keller}.
Finally, we observe that the entries in the third row of
$D{\mathcal{F}}$ can be computed directly from the definition
$(\ref{eq:ndef})$ of $N$
and that the computation of $\frac{\partial F_i}{\partial \lambda}$
is analogous to the computation of $\frac{\partial F_i}{\partial
\eta_j}$, $i,j=1,2$.
\subsection{The Results}
In this section we apply the continuation algorithm described in
section \ref{contalg} to the system $(\ref{eq:sysivp4})$
and examine the bifurcation properties of
periodic solutions. We observe that for $\mu \in [1.0,1.4]$, the path
of periodic solutions is S-shaped, thus for fixed $\lambda$ in some
interval $(\underline{\lambda},\overline{\lambda})$, three periodic
solutions to $(\ref{eq:sysivp4})$ exist. Moreover, we see that
bifurcation from single to
multiple periodic solutions occurs at a small value of $\lambda =
\underline{\lambda}$. Finally, we observe that as $\mu$ increases,
$\underline{\lambda}$ decreases. These results are described in
Experiments 1, 2, and 3.
Consider the linearization
\[
\theta''+.01\theta'+2.4\theta=\lambda \sin \mu t
\]
of $(\ref{eq:torode2})$ and note that
for the undamped equation, resonant solutions occur at $\mu^* \approx
1.55$. In our experiments we find that if $\mu$ is smaller than, but
close to, $\mu^*$, we cannot compute large amplitude periodic
solutions; the algorithm does not converge to
a periodic solution (see Experiment 4). If $\mu > \mu^*$, the
amplitude of periodic solutions increases with $\lambda$, but
bifurcation from single to multiple periodic solutions does not
occur.
\paragraph{Experiment 1} $\mu = 1.0$; See Figure \ref{fig41}\\
In this experiment we see that the path of periodic solutions is
S-shaped and that bifurcation from single to multiple periodic solutions occurs
at $\underline{\lambda} \approx 0.0126$. If $\lambda <
\underline{\lambda}$, $(\ref{eq:sysivp4})$ has a unique
periodic solution of small amplitude, but if $\lambda \in
(\underline{\lambda},.623)$, multiple periodic solutions exist. Whether the
cross-section oscillates with small or large amplitude depends only on
the initial conditions. This is consistent with the results in
\cite{kn:MCK}. Moreover, we observe that the amplitude of the large
oscillations is close to one radian, as was observed at Tacoma
Narrows on the day of its collapse, \cite{kn:FWA}.
\begin{figure}
\begin{center}
\epsfxsize=10cm
\epsffile{fig41.eps}
\end{center}
\caption{\label{fig41} Experiment 1}
\end{figure}
\paragraph{Experiment 2} $\mu = 1.2$; See Figure \ref{fig42}\\
Again we see that the path of periodic solutions is S-shaped and that
bifurcation from single to multiple periodic solutions occurs at
$\underline{\lambda} \approx 0.0117$.
\begin{figure}
\begin{center}
\epsfxsize=10cm
\epsffile{fig42.eps}
\end{center}
\caption{\label{fig42} Experiment 2}
\end{figure}
\paragraph{Experiment 3} $\mu = 1.4$; See Figure \ref{fig43}\\
Again we see that the path of periodic solutions is S-shaped and that
bifurcation from single to multiple periodic solutions occurs at
$\underline{\lambda} \approx 0.0088$.
\begin{figure}
\begin{center}
\epsfxsize=10cm
\epsffile{fig43.eps}
\end{center}
\caption{\label{fig43} Experiment 3}
\end{figure}
\paragraph{Experiment 4} $\mu = 1.5$; See Figure \ref{fig44}\\
From the curve in Figure \ref{fig44}, we see that multiple solutions exist
for $\lambda \in [.0197,.0213]$. However, because $\mu=1.5$ is close
to the resonant frequency $\mu^* \approx 1.55$, the algorithm fails
as the amplitude of the solution increases.
\begin{figure}
\begin{center}
\epsfxsize=10cm
\epsffile{fig44.eps}
\end{center}
\caption{\label{fig44} Experiment 4}
\end{figure}
\paragraph{Experiment 5} $\mu = 1.6, \mu = 1.7, \mu=2.0, \mu=2.6$; See
Figure \ref{fig45}\\
From top to bottom, the curves pictured correspond to $\mu = 1.6, \mu
= 1.7, \mu=2.0$ and $\mu=2.6$.
We see here that beyond the resonant frequency $\mu^*$, the amplitude
of periodic solutions grows with $\lambda$, but that the growth is
slower for higher frequencies. Moreover, for larger
$\mu$, the amplitude of the periodic solution grows nearly linearly with
$\lambda$.
\begin{figure}
\begin{center}
\epsfxsize=10cm
\epsffile{fig45.eps}
\end{center}
\caption{\label{fig45} Experiment 2}
\end{figure}
\section{Conclusion}
In this paper, we have presented substantial progress on the
understanding of the structure of the periodic solutions for the
forced torsional equation in the parameter ranges where the
large-amplitude oscillation was observed at Tacoma Narrows.
The multiple solutions exist
over roughly the right interval of amplitude and period, consistent
with the historical evidence.
Nonetheless, some intriguing questions remain.
First, we have not presented any proof that the large amplitude
solutions persist if there is small damping. This is certainly
supported by the numerical evidence but remains open as a mathematical
question. Obviously, there must be some relationship between the small
forcing term and the small damping, since in the absence of any
forcing, all motion must eventually decay.
Second, the original paper \cite{kn:MCK}, contains two other areas
which have not been studied here. One is the observation of
subharmonic solutions at integer multiples of the
frequency of forcing studied here. More sophisticated continuation
methods may be needed to see exactly where these families of solutions
lie on the big bifurcation picture.
Third, the relationship between the vertical
and torsional motions has not been addressed. Presumably, there are
multiplicity theorems for the coupled ordinary differential equation
system. But more importantly, the role of the loosening of the cables
in large vertical motion and the instantaneous transition to
torsional motion, so dramatically portrayed in \cite{kn:MCK}, has not
been addressed here. More light will be shed on this question when
the global bifurcation picture of the coupled system is available.
The relationship between the instability of the vertical
large-amplitude solutions and the torsional solutions should prove to be
particularly interesting. Progress in this direction has been made in
\cite{dh}, \cite{mckm}, and \cite{ot}.
Fourth, the question of the spatial dependence and nodal structure of
the solutions of a nonlinear beam equation will give more insight into
the observed large-amplitude motions captured on the famous film of
the Tacoma Narrows. This question is studied in \cite{kristenthesis}.
In short, we expect this
subject to be a fruitful area of study (and insight into history) for
some time to come.
\begin{thebibliography}{99}
\bibitem{kn:FWA}
O.H. Amann, T. von K\'{a}rm\'{a}n, and G.B. Woodruff.
\newblock {\em The Failure of the Tacoma Narrows Bridge}.
\newblock Federal Works Agency, 1941.
\bibitem{kn:BDH}
Blanchard, Devaney, and Hall.
\newblock {\em Differential Equations}.
\newblock Brooks/Cole Publishing Company, Pacific Grove, 1998.
\bibitem{kn:Castro}
A. Castro and A. C. Lazer.
\newblock {\em Critical point theory and the number of solutions of a
nonlinear Dirichlet problem.}
\newblock Ann. Mat. Pura Appl. (4), 120:113--137, 1979.
\bibitem{kn:Jen}
Y.S. Choi, K.C. Jen and P.J. McKenna.
\newblock {\em The structure of the solution set for periodic
oscillations in a suspension bridge model}.
\newblock IMA J. Appl. Math., 47:283-306, 1991.
\bibitem{kn:Codd}
E.A. Coddington and N. Levinson.
\newblock {\em Theory of Ordinary Differential Equations}.
\newblock McGraw-Hill, New York, 1955.
\bibitem{dh}
S.H. Doole and S.J. Hogan.
\newblock{\em Non-linear dynamics of the extended Lazer-McKenna bridge
oscillation model.} \newblock Dyn. Stab. Syst., {\bf 15} (2000), no. 1
43--58.
\bibitem{kn:Keller}
H. B. Keller.
\newblock{\it Lectures on Numerical Methods in Bifurcation Problems}.
\newblock Springer-Verlag, Berlin, 1987.
\bibitem{kn:MCK}
P.J. McKenna.
\newblock {\em Large torsional oscillations in suspension bridges
revisited: fixing an old approximation}.
\newblock The American Mathematical Monthly, 106:1-18, 1999.
\bibitem{mckm}
P.J. McKenna and K.S. Moore.
\newblock {\em The global structure of periodic solutions of a
suspension bridge mechanical model.}
\newblock In preparation.
\newblock In preparation.
\bibitem{ot}
P.J. McKenna and Cilliam O'Tuama.
\newblock {\em Large torsional oscillations in suspension bridges
revisited yet again: vertical forcing creates torsional response.}
\newblock In preparation.
\bibitem{kristenthesis}
K.S. Moore.
\newblock {\em Large amplitude torsional oscillations in a nonlinearly
suspended beam: a theoretical and numerical investigation}.
\newblock Dissertation, University of Connecticut, 1999.
\bibitem{kn:Scan}
R.H. Scanlan and J.J. Tomko.
\newblock {\em Airfoil and bridge deck flutter derivatives}.
\newblock Proc. Am. Soc. Civ. Eng. Eng. Mech. Division, EM6, 1717-1737,
1971.
\end {thebibliography} \medskip
\noindent{\sc Joe McKenna}\\
Department of Mathematics\\
University College\\
Cork, Ireland \\
e-mail: mckenna@math.uconn.edu \medskip
\noindent{Kristen S. Moore}\\
Department of Mathematics\\
University of Michigan\\
Ann Arbor, MI 48109-1109, USA \\
e-mail: ksmoore@math.lsa.umich.edu
\end {document}