Nonlinear Differential Equations, Electron. J. Diff. Eqns., Conf. 05, 2000, pp. 183–199 http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu or ejde.math.unt.edu (login: ftp)

Multiple periodic solutions to a suspension bridge ordinary differential equation ^{∗}

P. J. McKenna & K. S. Moore

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.

1 Introduction

In [8], 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

∗Mathematics Subject Classifications: 34C25, 34A47.

Key words: Torsional oscillations, suspension bridge.

\begin{matrix}c \\ ◯ \end{matrix}2000 Southwest Texas State University.

Published October 25, 2000.

Figure 1: A simple model of the center span

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.

2 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 1). 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 θ(t) denote the angle of the rod from horizontal at time t (see figure 2).

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 2 we see that the extension in the right hand cable is (y − l sin θ), and hence the force exerted by the right hand cable is

−K(y − l sin θ)^{+} = \{ \begin{matrix}−K(y − l sin θ) & if y − l sin θ ≥ 0 \\ 0 & if y − l sin θ < 0 \end{matrix}

where u^{+} = max{u, 0}. Similarly, the extension in the left hand cable is (y + l sin θ) and the force exerted by the left hand cable is −K(y + l sin θ)^{+}.

In [8], the author showed that the torsional and vertical motion satisfy

θ′′ = \frac{3K}{ml} cos θ[(y − l sin θ)^{+} − (y + l sin θ)^{+}] − δ_{1}θ′ + f(t)  (1)
y′′ = −\frac{K}{m} [(y − l sin θ)^{+} + (y + l sin θ)^{+}] − δ_{2}y′ + g ,

where δ_{1}, δ_{2} are damping constants, g is the force due to gravity, and f(t) is the external force at time t.

P. J. McKenna & K. S. Moore 185 Figure 2: A horizontal cross section of the center span

Assuming that the cables never lose tension, we have y ± l sin θ ≥ 0 and hence (y ± l sin θ)^{+} = (y ± l sin θ). Thus, the equations (1) become uncoupled and the torsional and vertical motion satisfy

θ′′ = −\frac{6K}{m} cos θ sin θ − δ_{1}θ′ + f(t)  (2)
y′′ = −\frac{2K}{m} y − δ_{2}y′ + g  (3)

respectively.

The equation for the vertical motion, (3), is simply the equation for a damped, forced, linear harmonic oscillator and the behavior of its solutions is well known, [2]. We will study the equation for the torsional motion (2), which is the damped, forced, pendulum equation and is known to possess chaotic solutions, [2].

To choose the physical constants K, m, δ := δ_{2} and the external forcing term f(t), we rely on [1], [8], and [12]. From [1], we choose m = 2500 and δ = .01. To determine K, from [1] 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 [12], so in (2) we choose f(t) = λ sin(μt) where λ ∈ [0, 0.06] is chosen to produce the appropriate behavior near equilibrium and the frequency μ 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 μ ∈ [1.2, 1.6].

3 A Multiplicity Theorem

In this section, we prove that undamped equations of the form (2) have multiple periodic weak solutions. To some extent, this result is inspired by the earlier results of Castro and Lazer, [3], where they obtained multiple solutions of elliptic boundary value problems once a critical value, measured by f(0), passes across key eigenvalues.

Let Ω = (−π, π), H = {u ∈ L^{2}(Ω)|u(t) = −u(−t)} and define Lu = −u′′. For u ∈ H, let ∥u∥ = ∥u∥_{L^{2}(Ω)} = (∫_{Ω} |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 (2) as

Lu = b sin 2u − εh(t)  (4)
u(−π) = u(π), u′(−π) = u′(π)

We observe that the eigenvalues and eigenfunctions of L are

λ_{n} = n^{2} , ϕ_{n}(t) = sin nt  (5)

hence L^{−1} exists, is compact, and ∥L^{−1}∥ = 1.

Definition We say that u ∈ H is a solution to (4) if

u − L^{−1}(b sin 2u − εh) = 0.  (6)

We will prove the existence of multiple solutions u ∈ H to the problem (4).

Theorem 3.1 Let h ∈ H with ∥h∥ ≤ 1 and let b ∈ (\frac{1}{2}, 2). Then there exists ε_{0} > 0 such that if |ε| < ε_{0}, (4) has at least two solutions.

Proof Let ε < ε_{0}; we will determine the value of ε_{0} later. Note that by (5) and our choice of b, 2b is not an eigenvalue of L.

Define T_{1} : H → H by

T_{1}(u) = u − L^{−1}(b sin 2u − εh)

and note that, by (6), zeros of T_{1} correspond to solutions of (4). 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

(D1): there exists R_{0} > 0 such that for R > R_{0}, deg(T_{1}, B_{R}(0), 0) = 1

(D2): there exists γ ∈ (0, R_{0}) such that deg(T_{1}, B_{γ}(0), 0) = −1

Then, since deg(T_{1}, B_{γ}(0), 0) ≠ 0, there exists a zero of T_{1} (i.e., a solution of (4)) in B_{γ}(0). Moreover, by the additivity property of degree, we have deg(T_{1}, B_{R}(0)\\overline{B_{γ}(0)}, 0) ≠ 0 and hence (4) has second solution in the annulus B_{R}(0)\\overline{B_{γ}(0)}.

To establish (D1), define

T_{β}u = u − βL^{−1}(b sin 2u − εh)

for β ∈ [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_{β}, B_{R}(0), 0) is constant provided that 0 ∉ T_{β}(∂B_{R}(0)) for all β ∈ [0, 1].

Fix β ∈ [0, 1] and suppose u ∈ H solves T_{β}u = 0. We will show that u is bounded above by some R_{0} > 0 and that this bound is independent of β.

Since T_{β}u = 0, we have

∥u∥ = β∥L^{−1}(b sin 2u − εh)∥
≤ β[ε_{0} + b∥ sin 2u∥]
≤ β[ε_{0} + bm(Ω)^{1/2}]
≤ [ε_{0} + bm(Ω)^{1/2}]
= [ε_{0} + b\sqrt{2π}] < R_{0}

if we choose R_{0} > ε_{0} + b\sqrt{2π}.

Thus, for R > R_{0} we have

deg(T_{1}, B_{R}(0), 0) = deg(T_{0}, B_{R}(0), 0) = 1  (7)

and (D1) above holds.

To establish (D2), define

T_{μ}u = u − (1 − μ)L^{−1}(2bu) − μL^{−1}(b sin 2u − εh)

for μ ∈ [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 γ > 0

deg(T_{1}, B_{γ}(0), 0) = deg(T_{0}, B_{γ}(0), 0) = −1.  (8)

We will show first that for any γ > 0,

deg(T_{0}, B_{γ}(0), 0) = deg(I − 2bL^{−1}, B_{γ}(0), 0) = −1.  (9)

Consider the finite dimensional subspace M_{N} :=span{ϕ_{n}}_{n=1}^{N} of H and recall that, by compactness, 2bL^{−1} can be approximated in operator norm by the operators B_{N} : M_{N} → M_{N} given by

B_{N}(u) = 2b ∑_{n=1}^{N} \frac{c_{n}}{λ_{n}} ϕ_{n}.

By definition, for N sufficiently large,

deg(T_{0}, B_{γ}(0), 0) = deg(I − B_{N}, B_{γ}(0) ∩ M_{N}, 0)
= ∑_{u∈(I−B_{N})^{−1}(0)} sign J_{I−B_{N}} (u)

where J_{ϕ}(u) is the Jacobian determinant of ϕ at u.

Since I −B_{N} can be identified with an N ×N diagonal matrix whose entries are 1 − \frac{2b}{λ_{n}}, we have

deg(I − B_{N}, B_{γ}(0) ∩ M_{N}, 0) = sign ∏_{n=1}^{N}(1 − \frac{2b}{λ_{n}}).

Since b ∈ (\frac{1}{2}, 2) and λ_{n} = n^{2}, the only negative value of 1− \frac{2b}{λ_{n}} occurs at λ_{1} = 1. Since the eigenvalues of L are simple, we have

deg(I − B_{N}, B_{γ}(0) ∩ M_{N}, 0) = −1

and (9) holds.

To establish (8), and hence (D2), we must show that for some γ > 0,

deg(T_{1}, B_{γ}(0), 0) = deg(T_{0}, B_{γ}(0), 0).  (10)

The homotopy property of degree ensures that deg(T_{μ}, B_{γ}(0), 0) is constant provided 0 ∉ T_{μ}(∂B_{γ}(0)) for all μ ∈ [0, 1]. Observe that u = 0 is the only zero of T_{0} since 2b is not an eigenvalue of L. Fix μ ∈ (0, 1] and suppose that u ≠ 0 solves T_{μ}u = 0. Set ∥u∥ = γ̃_{μ}; we will show that γ̃_{μ} is bounded below by some γ_{μ} > 0.

Set ψ = \frac{u}{∥u∥} = \frac{u}{γ̃_{μ}}. We claim first that there exists a compact set K with ψ ∈ K. Since u is a zero of T_{μ}, u satisfies

Lu − 2bu = μ(b sin 2u − 2bu − εh)  (11)

and therefore

Lu = μ(b sin 2u − εh) + (1 − μ)2bu.

Thus,

∥Lu∥ ≤ ∥b sin 2u − εh∥ + 2b∥u∥
≤ b∥ sin 2u∥ + 2bγ̃_{μ} + ε_{0}
≤ bm(Ω)^{\frac{1}{2}} + 2bγ̃_{μ} + ε_{0} ≤ \sqrt{2π}b + 2bγ̃_{μ} + ε_{0}.

It follows that u ∈ \overline{L^{−1}(B_{\sqrt{2π}b+2bγ̃_{μ}+ε_{0}}(0))}, which is compact since L^{−1} is compact, and thus there exists a compact set K with ψ ∈ K.

Since u satisfies (11), we have

∥u − L^{−1}2bu∥ = μ∥L^{−1}(b sin 2u − 2bu − εh)∥.  (12)

Denote the left and right hand sides of (12) by LHS and RHS respectively. Since 2b is not an eigenvalue of L, we have Lψ − 2bψ ≠ 0 and hence

inf_{ψ∈K} ∥ψ − L^{−1}2bψ∥ = α > 0

and for our u we have

LHS = ∥u − L^{−1}2bu∥ ≥ αγ̃_{μ}.  (13)

Now considering the right hand side of (12), we have

RHS ≤ μ∥b sin 2u − 2bu − εh∥
≤ μ[ε_{0} + ∥b sin 2u − 2bu∥]
= μ[ε_{0} + ∥b sin(2γ̃_{μ}ψ) − 2bγ̃_{μ}ψ∥].

We claim that if γ̃_{μ} is sufficiently small, RHS < αγ̃_{μ}, which contradicts (13). To establish this, we must first prove the following lemma.

Lemma 3.2 Let K, H, b be as above and denote γ̃_{μ} by η. Then there exists a function δ : (0, ∞) → (0, ∞) such that

(L1): ∥b sin(2ηψ) − 2bηψ∥ ≤ ηδ(η)

(L2): δ(η) → 0 as η → 0

hold for all ψ ∈ K, η > 0 Proof Define δ : (0, ∞) → (0, ∞) by

δ(η) = max_{ψ∈K} ∥\frac{b}{η} sin 2ηψ − 2bψ∥

and note that (L1) above is satisfied.

To show that (L2) holds, we must show that f_{η}(ψ) := ∥\frac{b}{η} sin 2ηψ −2bψ∥ → 0 uniformly on K as η → 0. To show that f_{η} → 0 uniformly on K, we’ll show that f_{η}(ψ) → 0 for each ψ ∈ K and that the family F := {f_{η}} is equicontinuous on K.

Choose ψ ∈ K, ψ ≠ 0 and note that

\frac{b}{η} sin 2ηψ − 2bψ = b\frac{sin 2ηψ}{2ηψ} · 2ψ − 2bψ → 0

as η → 0, hence

|\frac{b}{η} sin 2ηψ − 2bψ|^{2} → 0.

Moreover, since | sin w| ≤ |w| and since w^{2} is convex, we have

|\frac{b}{η} sin 2ηψ − 2bψ|^{2} = 4|\frac{1}{2} \frac{b}{η} sin 2ηψ + \frac{1}{2}(−2bψ)|^{2}
≤ 4[\frac{1}{2}|\frac{b}{η} sin 2ηψ|^{2} + \frac{1}{2}|2bψ|^{2}]
≤ 2[\frac{b^{2}}{η^{2}} |2ηψ|^{2} + |2bψ|^{2}] = 16b^{2}ψ^{2} ∈ L^{1}

since ψ ∈ L^{2}. By the Dominated Convergence Theorem, we conclude that f_{η}(ψ) → 0 for each ψ ∈ K, as desired.

To see that the family F = {f_{η}} is equicontinuous on K, choose ε̃ > 0 and ψ, ψ̃ ∈ K. We have that

|f_{η}(ψ) − f_{η}(ψ̃)| = |∥\frac{b}{η} sin 2ηψ − 2bψ∥ − ∥\frac{b}{η} sin 2ηψ̃ − 2bψ̃∥|
≤ ∥\frac{b}{η} (sin 2ηψ − sin 2ηψ̃) − 2b(ψ − ψ̃)∥
≤ \frac{b}{η} ∥ sin 2ηψ − 2ηψ̃∥ + 2b∥ψ − ψ̃∥
≤ \frac{b}{η} ∥2ηψ − 2ηψ̃∥ + 2b∥ψ − ψ̃∥
= 4b∥ψ − ψ̃∥ < ε̃

provided ∥ψ − ψ̃∥ < \frac{ε̃}{4b}.

Since f_{η} are equicontinuous on K and converge pointwise on K, we have that f_{η} converge uniformly on K and hence (L2) holds. ♢

Returning now to the proof of the theorem and invoking the lemma above, we have that

RHS ≤ μ[ε_{0} + γ̃_{μ}δ(γ̃_{μ})].

Assume now that ε_{0} < \frac{1}{2}αγ̃_{μ}. Since δ → 0, there exists γ_{μ} > 0 such that if γ̃_{μ} < γ_{μ}, then δ(γ̃_{μ}) < \frac{1}{2}α and hence

RHS < αγ̃_{μ}.

But this contradicts (13). Therefore γ̃_{μ} is bounded below by γ_{μ} > 0. Note also that because of the factor of μ on the right hand side of the above inequalities, we can take 0 < γ_{1} ≤ γ_{μ} for μ ∈ (0, 1], and hence we can choose γ ∈ (0, γ_{1}). ♢

4 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 longterm 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 (2). Using the physical constants from the engineers’ reports of the Tacoma Narrows failure and the forcing term described in section 1 yields the ordinary differential equation

θ′′ = −2.4 cos θ sin θ − .01θ′ + λ sin μt.  (14)

In [8], the author prescribed the initial position θ(0) and velocity θ′(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, [1].

As in [8], we compute periodic solutions to (14). Moreover, we describe a numerical continuation algorithm and we use it to plot the amplitude of periodic solutions versus the amplitude λ of the external forcing term. We demonstrate that for small λ, three periodic solutions to (14) 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 λ.

4.1 The Algorithm

We wish to generate the bifurcation curve for periodic solutions to (14) as the amplitude λ of the external forcing term varies; i.e., we wish to plot the amplitude A_{λ} of the periodic solution versus λ. Define T = 2π/μ. To compute T periodic solutions to (14), we rewrite it as a first order system as follows. Let

U = [ \begin{matrix}u \\ v \end{matrix} ] = [ \begin{matrix}θ \\ θ′ \end{matrix} ] .  (15)

Then solutions θ of (14) correspond to solutions U of the first order system

u′ = v  (16)
v′ = −2.4 cos u sin u − .01v + λ sin μt

We employ a shooting algorithm to find periodic solutions to the above system. In the shooting algorithm, we search for initial conditions η = [η_{1}, η_{2}]^{T} so that the solution to the initial value problem

u′ = v
v′ = −2.4 cos u sin u − .01v + λ sin μt  (17)
u(0) = η_{1}, v(0) = η_{2}

is T -periodic, i.e., we search for initial conditions η so that the solution U to (17) satisfies

U(T) = U(0) = η.

Define F : R^{2} → R^{2} by

F (η) = η − U(T).  (18)

Then searching for initial conditions η which yield a T -periodic solution to (17) 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 η^{0} which yield a T -periodic solution to (17). We then solve the initial value problem (17) for t ∈ [0, T] and test whether F (η^{0}) = 0. If not, we update η^{0} via Newton’s method:

η^{1} = η^{0} − DF (η^{0})^{−1}F (η^{0})

where

DF = [ \begin{matrix}\frac{∂F_{1}}{∂η_{1}} & \frac{∂F_{1}}{∂η_{2}} \\ \frac{∂F_{2}}{∂η_{1}} & \frac{∂F_{2}}{∂η_{2}} \end{matrix} ]  (19)

is the Jacobian matrix of F.

To compute this matrix, observe that by our definition of F in (18), we have

DF = [ \begin{matrix}\frac{∂F_{1}}{∂η_{1}} & \frac{∂F_{1}}{∂η_{2}} \\ \frac{∂F_{2}}{∂η_{1}} & \frac{∂F_{2}}{∂η_{2}} \end{matrix} ] = [ \begin{matrix}1 − \frac{∂u}{∂η_{1}} (T) & −\frac{∂u}{∂η_{2}} (T) \\ −\frac{∂v}{∂η_{1}} (T) & 1 − \frac{∂v}{∂η_{2}} (T) \end{matrix} ] .  (20)

To compute \frac{∂u}{∂η_{i}} , \frac{∂v}{∂η_{i}} , i = 1, 2, we recall that the solution to the initial value problem is continuously differentiable with respect to initial conditions and parameters, [5], and hence we can differentiate (17) with respect to η_{1}. This yields the first order system in \frac{∂u}{η_{1}}, \frac{∂v}{η_{1}}

( \frac{∂u}{∂η_{1}} ; \frac{∂v}{∂η_{1}} )′ = ( −2.4 cos(2u)\frac{∂u}{∂η_{1}} − .01\frac{∂v}{∂η_{1}}^{\frac{∂v}{∂η_{1}}} ); ( \frac{∂u}{∂η_{1}} ; \frac{∂v}{∂η_{1}} ) (0) = ( \begin{matrix}1 \\ 0 \end{matrix} ) (21)

Hence [\frac{∂u}{∂η_{1}} , \frac{∂v}{∂η_{1}} ]^{T} solves the initial value problem

( \begin{matrix}w \\ z \end{matrix} )′ = ( \begin{matrix}z \\ −2.4 cos(2u)w − .01z \end{matrix} ); ( \begin{matrix}w \\ z \end{matrix} ) (0) = ( \begin{matrix}1 \\ 0 \end{matrix} )  (22)

where u is the solution to the initial value problem (17).

Similarly, differentiating (17) with respect to η_{2}, we see that [\frac{∂u}{∂η_{2}} , \frac{∂v}{∂η_{2}} ]^{T} solves the initial value problem

( \begin{matrix}w \\ z \end{matrix} )′ = ( \begin{matrix}z \\ −2.4 cos(2u)w − .01z \end{matrix} ); ( \begin{matrix}w \\ z \end{matrix} ) (0) = ( \begin{matrix}0 \\ 1 \end{matrix} ) .  (23)

Therefore, to compute the entries of the Jacobian matrix (20), we must solve the initial value problems (17), (22), and (23).

The above algorithm (Shooting Algorithm with Newton Update) suggests a method by which we can plot the amplitude A_{λ} of periodic solutions to (14) versus the amplitude λ of the external forcing term. Specifically, fix λ = λ̂, use the algorithm described above to compute the initial conditions which yield a T periodic solution to (17), record the amplitude Aλ̂ of the periodic solution, increment λ, and repeat.

Unfortunately, this algorithm fails at bifurcation points because, under natural parameterization, the Jacobian matrix (20) 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 [4], [7] for details.

Let (η^{0}, λ^{0}) be a zero of F and set (η^{0}, λ^{0}) = (η(s_{0}), λ(s_{0})). We will introduce the pseudoarclength normalization N given by

N(η, λ, s) = ψ \frac{∥η(s) − η(s_{0})∥^{2}}{s − s_{0}} + (1 − ψ)\frac{|λ(s) − λ(s_{0})|^{2}}{s − s_{0}} − (s − s_{0})  (24)

where ψ ∈ (0, 1) and s is a chord-length parameter, [4]. We use an algorithm similar to the Shooting Algorithm with Newton Update described above to search for zeros of the map F given by

F(η, λ, s) = ( \begin{matrix}F (η, λ) \\ N(η, λ, s) \end{matrix} ) .  (25)

The Jacobian matrix

DF = [ \begin{matrix}\frac{∂F_{1}}{∂η_{1}} & \frac{∂F_{1}}{∂η_{2}} & \frac{∂F_{1}}{∂λ} \\ \frac{∂F_{2}}{∂η_{1}} & \frac{∂F_{2}}{∂η_{2}} & \frac{∂F_{2}}{∂λ} \\ \frac{∂N}{∂η_{1}} & \frac{∂N}{∂η_{2}} & \frac{∂N}{∂λ} \end{matrix} ]  (26)

is nonsingular even if the matrix DF is singular, [7]. Finally, we observe that the entries in the third row of DF can be computed directly from the definition (24) of N and that the computation of \frac{∂F_{i}}{∂λ} is analogous to the computation of \frac{∂F_{i}}{∂η_{j}}, i, j = 1, 2.

4.2 The Results

In this section we apply the continuation algorithm described in section 4.1 to the system (17) and examine the bifurcation properties of periodic solutions. We observe that for μ ∈ [1.0, 1.4], the path of periodic solutions is S-shaped, thus for fixed λ in some interval (\underline{λ}, \overline{λ}), three periodic solutions to (17) exist. Moreover, we see that bifurcation from single to multiple periodic solutions occurs at a small value of λ = \underline{λ}. Finally, we observe that as μ increases, \underline{λ} decreases. These results are described in Experiments 1, 2, and 3.

Consider the linearization

θ′′ + .01θ′ + 2.4θ = λ sin μt

of (14) and note that for the undamped equation, resonant solutions occur at μ^{∗} ≈ 1.55. In our experiments we find that if μ is smaller than, but close to, μ^{∗}, we cannot compute large amplitude periodic solutions; the algorithm does not converge to a periodic solution (see Experiment 4). If μ > μ^{∗}, the amplitude of periodic solutions increases with λ, but bifurcation from single to multiple periodic solutions does not occur.

Experiment 1 μ = 1.0; See Figure 3 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{λ} ≈ 0.0126. If λ < \underline{λ}, (17) has a unique periodic solution of small amplitude, but if λ ∈ (\underline{λ}, .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 [8]. 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, [1].

Experiment 2 μ = 1.2; See Figure 4 Again we see that the path of periodic solutions is S-shaped and that bifurcation from single to multiple periodic solutions occurs at \underline{λ} ≈ 0.0117.

Experiment 3 μ = 1.4; See Figure 5 Again we see that the path of periodic solutions is S-shaped and that bifurcation from single to multiple periodic solutions occurs at \underline{λ} ≈ 0.0088.

Experiment 4 μ = 1.5; See Figure 6 From the curve in Figure 6, we see that multiple solutions exist for λ ∈ [.0197, .0213].

P. J. McKenna & K. S. Moore 195 Figure 3: Experiment 1 Figure 4: Experiment 2

Figure 5: Experiment 3

However, because μ = 1.5 is close to the resonant frequency μ^{∗} ≈ 1.55, the algorithm fails as the amplitude of the solution increases.

Experiment 5 μ = 1.6, μ = 1.7, μ = 2.0, μ = 2.6; See Figure 7 From top to bottom, the curves pictured correspond to μ = 1.6, μ = 1.7, μ = 2.0 and μ = 2.6. We see here that beyond the resonant frequency μ^{∗}, the amplitude of periodic solutions grows with λ, but that the growth is slower for higher frequencies. Moreover, for larger μ, the amplitude of the periodic solution grows nearly linearly with λ.

5 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 P. J. McKenna & K. S. Moore 197

Figure 6: Experiment 4 Figure 7: Experiment 2 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 [8], 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 [8], 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 [6], [9], and [10].

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 [11].

In short, we expect this subject to be a fruitful area of study (and insight into history) for some time to come.

References

[1] O.H. Amann, T. von Kármán, and G.B. Woodruff. The Failure of the Tacoma Narrows Bridge. Federal Works Agency, 1941.

[2] Blanchard, Devaney, and Hall. Differential Equations. Brooks/Cole Publishing Company, Pacific Grove, 1998.

[3] A. Castro and A. C. Lazer. Critical point theory and the number of solutions of a nonlinear Dirichlet problem. Ann. Mat. Pura Appl. (4), 120:113–137, 1979.

[4] Y.S. Choi, K.C. Jen and P.J. McKenna. The structure of the solution set for periodic oscillations in a suspension bridge model. IMA J. Appl. Math., 47:283-306, 1991.

[5] E.A. Coddington and N. Levinson. Theory of Ordinary Differential Equations. McGraw-Hill, New York, 1955.

[6] S.H. Doole and S.J. Hogan. Non-linear dynamics of the extended Lazer-McKenna bridge oscillation model. Dyn. Stab. Syst., 15 (2000), no. 1 43–58.

[7] H. B. Keller. Lectures on Numerical Methods in Bifurcation Problems. Springer-Verlag, Berlin, 1987.

[8] P.J. McKenna. Large torsional oscillations in suspension bridges revisited: fixing an old approximation. The American Mathematical Monthly, 106:1-18, 1999.

[9] P.J. McKenna and K.S. Moore. The global structure of periodic solutions of a suspension bridge mechanical model. In preparation.

In preparation.

[10] P.J. McKenna and Cilliam O’Tuama. Large torsional oscillations in suspension bridges revisited yet again: vertical forcing creates torsional response. In preparation.

[11] K.S. Moore. Large amplitude torsional oscillations in a nonlinearly suspended beam: a theoretical and numerical investigation. Dissertation, University of Connecticut, 1999.

[12] R.H. Scanlan and J.J. Tomko. Airfoil and bridge deck flutter derivatives. Proc. Am. Soc. Civ. Eng. Eng. Mech. Division, EM6, 1717-1737, 1971.

Joe McKenna Department of Mathematics University College Cork, Ireland e-mail: mckenna@math.uconn.edu

Kristen S. Moore Department of Mathematics University of Michigan Ann Arbor, MI 48109-1109, USA e-mail: ksmoore@math.lsa.umich.edu