\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}