## Abstract

The conjugate point theory of the calculus of variations is extended to apply to the buckling of an elastic rod in an external field. We show that the operator approach presented by Manning (Manning *et al*. 1998 *Proc. R. Soc. A* **454**, 3047–3074) can be used when the second-variation operator is an integrodifferential operator, rather than a differential operator as in the classical case. The external field is chosen to model two parallel ‘soft’ walls. We consider the examples of two-dimensional buckling under both pinned–pinned and clamped–clamped boundary conditions, as well as the three-dimensional clamped–clamped problem, where we consider the importance of the rod cross-section shape as it ranges from circular to extreme elliptical. For each of these problems, we find that in the appropriate limit, the soft-wall solutions approach a ‘hard-wall’ limit, and thus we make conjectures about these hard-wall contact equilibria and their stability. In the two-dimensional pinned–pinned case, this allows us to assign stability to the configurations reported by Holmes (Holmes *et al*. 1999 *Comput. Methods Appl. Mech. Eng.* **170**, 175–207) and reconsider the experimental results discussed therein.

## 1. Introduction

In the classical theory of the calculus of variations, the stability of an equilibrium of a functional is investigated by searching for conjugate points using the Euler–Jacobi equation (Gelfand & Fomin 1963). In a refinement owing to Morse (1951), the number of conjugate points is seen to be a stability index giving the dimension of the space of allowed variations that lower the energy. Here, we follow the functional analytic approach of Manning *et al*. (1998) to this classical problem in order to extend the stability theory to the case of an elastic rod buckling in an external field. The standard Euler–Jacobi theory does not apply to this problem, as the second-variation operator is an integrodifferential operator.

The standard problem of interest is to find critical points (equilibria) of a functionalsubject to the conditions * q*(0)=

**q**_{i},

*(1)=*

**q**

**q**_{f}. Equilibria satisfy the Euler–Lagrange equationStability is determined by analysing the second variation of

*E*, which is conveniently written in the formwhere is a second-order differential operator, and

*is an allowed variation in*

**ζ***, which must satisfy*

**q***(0)=*

**ζ***(1)=*

**ζ****0**. A condition for stability is that the equilibrium have no conjugate point, where a conjugate point is defined to be a value

*σ*<1 for which the Euler–Jacobi equation

*=*

**ζ****0**has a solution

*(not identically zero) with*

**ζ***(0)=*

**ζ***(*

**ζ***σ*)=

**0**. If there are conjugate points, then their number is called the ‘index’, which gives the dimension of the space of allowed variations on which the second variation

*δ*

^{2}

*E*is negative.

In Manning *et al*. (1998), this theory is extended to give analogous conjugate point results for the case of an isoperimetrically constrained problem, with the only change that the second-variation operator is replaced by the operator , where is a projection onto the orthogonal complement of the linearized constraints. This theory is thus couched in functional analytic language, but the computation of conjugate points is easily translated into an initial-value problem (IVP) for a system of ordinary differential equations (ODEs). The unknowns in this IVP are assembled into a ‘stability matrix’, and a conjugate point occurs when the determinant of this matrix vanishes.

Here, we show that this functional analytic approach can be used to extend the stability theory to problems not covered by the standard calculus of variations approach. In particular, we look at the example of a rod in an external field, where the second-variation operator contains an integral operator term in addition to the second-order differential operator . After a few basic spectral properties are proven for this operator, the theory from Manning *et al*. (1998) immediately applies, and provides the definition of ‘conjugate point’ appropriate to this problem. Furthermore, as in the isoperimetrically constrained example, the required computations to find conjugate points can be reduced to an IVP for a system of ODEs, despite the fact that the second variation contains an integral operator.

We illustrate this stability theory on several variants of the Euler buckling problem, in which one end of an intrinsically straight elastic rod is subjected to a compressive load *λ*. We consider both clamped–clamped and pinned–pinned boundary conditions, and we add the constraint that the only translational freedom of motion of the two ends is a compression towards each other. To this classic problem, we add an external field that models two parallel ‘soft’ walls surrounding the rod. Each wall exerts a repulsive force on each infinitesimal segment of the rod, and this force increases to positive infinity as the segment's lateral displacement approaches a wall spacing parameter *a*. For sufficiently large *λ*, the rod can buckle into non-trivial equilibrium configurations, which are dependent on the strength of, and spacing between, the walls, as well as on the shape of the rod cross-section in three-dimensional rods (we consider both circular and elliptical cross-sections).

Holmes *et al*. (2000) have considered a similar rod-in-a-field problem, for a two-dimensional pinned–pinned rod in a quadratic potential, combining theoretical analysis of the equilibria via a dynamical systems perspective with extensive numerical simulations of the global bifurcation diagram using the parallel simplex algorithm. In contrast, our focus is on determining stability via the theory summarized above, which we illustrate on a small set of branches of equilibria computed using the parameter continuation package (Doedel *et al*. 1991).

In addition, we may consider the hard-wall limit of our soft-wall potential as a means to impute stability properties to hard-wall equilibria, such as those computed in the two-dimensional pinned–pinned case by Holmes *et al*. (1999). As we show in §5*a*, this stability information explains, among other things, the experimentally observed hysteresis reported in Holmes *et al*. (1999). A similar approach could, in principle, be used to determine the stability of equilibria in other ‘hard constraint’ settings, such as the asymmetric wall problem considered by Roman & Pocheau (1999) or even more complicated geometries such as buckling on a cylinder, as considered recently by van der Heijden (2001). In these hard-constraint situations, a direct conjugate point determination is beyond the scope of the present theory; among other technical obstacles, the space of allowed variations is no longer a vector space, but only a half-space, as certain variations in the rod configuration would violate the hard constraint.

In addition to providing a useful tool for understanding hard constraints, one can also view the soft-wall problem considered here as an important step towards the development of a theory of stability for rod *self*-contact configurations, since a Debye–Huckel potential is a widely used model for the self-repulsion of DNA in supercoiled configurations.

## 2. Equilibrium equations for an elastic rod

Here, we summarize the standard theory of inextensible and unshearable elastic rods (Kirchhoff 1859; Love 1927; Dill 1992; Antman 1995), and present the calculus of variations problem relevant for the buckling of such a rod in the presence of an external field with potential energy *V*. We begin with the full three-dimensional theory, and then turn at the end of the section to the simpler two-dimensional case.

The configuration of a rod is described by a centre-line * r*(

*s*)=(

*x*(

*s*),

*y*(

*s*),

*z*(

*s*))—written as a function of arclength

*s*—and a set of directors {

**d**_{1}(

*s*),

**d**_{2}(

*s*),

**d**_{3}(

*s*)} that form an orthonormal frame giving the orientation of the cross-section of the rod at location

*s*. We choose

**d**_{1}to lie along one axis of the elliptical cross-section,

**d**_{2}to lie along the other axis and

**d**_{3}to be orthogonal to the cross-section. For convenience, we choose a length-scale so that 0≤

*s*≤1. Inextensibility and unshearability of the rod are imposed by requiring that

**d**_{3}equals the tangent vector to the centre-line, or, letting ′ denote differentiation with respect to

*s*,

The bending and twisting energy of the rod is defined in terms of the *strains*The directors and strains can be written as rational functions of Euler parameters where

### (a) Energy

We make the common assumption that *E* depends quadratically on the strains, and add a potential energy term *V*,Here, *K*_{i} are the bending (*i*=1, 2) and twisting (*i*=3) stiffnesses of the rod, *λ* is the imposed load, and *d*_{33} is the *z*-component of the director **d**_{3}. We letto represent a pair of soft walls parallel to the *y*–*z*-plane and located at *x*=*a* and *x*=−*a*.

For convenience, we define

### (b) Buckling constraints

We impose the following clamped–clamped constraints on the rod(2.1)The constraints on * r* require that the ends of the rod lie on the

*z*-axis (with the

*s*=0 end at the origin). The constraints on

*require that the directors at the ends of the rod are aligned with the standard axes so that, in particular, the tangent vectors at the ends are vertical. Figure 1 illustrates two sample three-dimensional configurations consistent with these constraints.*

**q**We note that since we begin our computations with a straight untwisted rod and never introduce relative twist between the ends, our computations will miss many highly twisted solutions consistent with (2.1), for example, those we would find by starting with a straight rod with an integer number of full twists before buckling begins. We note also that with these constraints, lateral forces in the *x*- and *y*-directions may be present, and that while the *z*-force is known to equal *λ*, these lateral forces are unknowns of the problem (and, indeed, the *x*-component depends on *s* due to the wall potential).

The stiffnesses *K*_{1} and *K*_{2} are proportional to the moments of inertia of the elliptical cross-section about **d**_{1} and **d**_{2}, respectively, and the boundary conditions on * q* mandate that in the unbuckled configuration of the rod,

**d**_{1}is perpendicular to the walls. Hence, when

*K*

_{2}<

*K*

_{1}, the unbuckled rod has its major axis parallel to the wall; when

*K*

_{2}>

*K*

_{1}, its major axis is perpendicular to the wall; and when

*K*

_{2}=

*K*

_{1}, the cross-section is circular.

### (c) The variational problem useful for computation

It will be helpful to view this variational problem in two slightly different ways. First, we can consider *E* to be a functional of both * r* and

*. Equilibria of the rod are then critical points of this functional*

**q***E*[

*,*

**r***] subject to the pointwise constraint and the boundary conditions (2.1). These equilibria can be found by solving the standard Euler–Lagrange equations, which can be written in the following first-order Hamiltonian form (Dichmann*

**q***et al*. 1996), with the variable conjugate to

*(and representing the force in the rod) and the variable conjugate to*

**r***:Here, and throughout, subscripting by*

**q***or*

**r***denotes partial differentiation, for example, the term (*

**q**

**d**_{3})

_{q}represents the 3×4 matrix ∂

**d**_{3}/∂

*. We solve this system of first-order ODEs subject to the 13 boundary conditions from equation (2.1), plus the boundary condition*

**q***n*

_{3}(1)=−

*λ*representing the imposed load. Solutions of this boundary-value problem (BVP) were obtained using Auto, beginning with the unbuckled configuration

*(*

**r***s*)=(0, 0,

*s*),

*(*

**q***s*)=(0, 0, 0, 1), which is an equilibrium for any value of

*λ*. Parameter continuation in

*λ*was then performed to track families of equilibria, including the detection of bifurcating branches.1

### (d) The variational problem useful for stability theory

A second perspective that we will use for determining stability is to view *E* as a functional of * q* alone. This perspective is possible because

*is completely determined from*

**r***via the constraint combined with the initial condition*

**q***(0)=*

**r****0**,In this view, we rewrite the boundary conditions

*x*(1)=

*y*(1)=0 as isoperimetric constraintsThus, we seek critical points of the functionalsubject to

### (e) The two-dimensional case

In the two-dimensional case, the rod can be described by a single unknown—the angle *θ*(*s*) that the rod's tangent vector makes with the upward axis at arclength position *s*. For consistency with the three-dimensional coordinates, we denote the centre-line * r*(

*s*) by (

*x*(

*s*),

*z*(

*s*)), and the inextensibility–unshearability constraint simply takes the form . In place of the three strains, we have the quantity

*θ*′(

*s*), which represents the rate of bending at position

*s*.

The energy takes the form(a stiffness parameter *K*_{1} that multiplies the bending energy term *θ*′(*s*)^{2}/2 was removed by choice of units). We consider two sets of boundary conditions: the clamped–clamped case(2.2)and the pinned–pinned case(2.3)

The equations used to determine equilibria aresubject to either the five boundary conditions in equations (2.2) or those in (2.3).

For the variational problem relevant for stability, we consider the energy as a functional of *θ* only, which takes the formsubject to either the clamped–clamped conditionsor the pinned–pinned conditions

## 3. Computing the index for the soft wall

We now turn to the stability theory that is the focus of this article. As before, we begin with a detailed presentation of the full three-dimensional problem, and turn then to a summary of the simpler two-dimensional analogue.

Assume that we have found an equilibrium **q**_{0} of *E*[* q*]. To determine its stability, we consider

*allowed variations δ*of

**q***; i.e. functions*

**q***δ*for which

**q**

**q**_{0}+

*ϵδ*obeys the boundary conditions and constraints on

**q***to first order in*

**q***ϵ*(3.1)where a superscript 0 denotes evaluation at the equilibrium

**q**_{0}.

We define the index to be the dimension of a maximal (deleted) subspace of allowed variations *δ q* on which for

*ϵ*>0 sufficiently small. If the index is zero (i.e. for all allowed variations

*δ*, for

**q***ϵ*sufficiently small), then we will say that

**q**_{0}is a

*stable*equilibrium.

We note first a property of the functional *E* peculiar to the example at hand, namely that it is invariant to a scaling of * q*, i.e. for any function

*α*(

*s*),This degeneracy arises from our use of to represent the three-dimensional space

*SO*(3) of directors. If we write

*δ*(

**q***s*) as

*β*(

*s*)

**q**_{0}(

*s*)+

*(*

**v***s*), where at each

*s*,

*(*

**v***s*)⊥

**q**_{0}(

*s*), thenTherefore, we can say that for

*ϵ*is sufficiently small if, and only if, for

*ϵ*is sufficiently small. Thus, it is sufficient to consider only those allowed variations

*δ*that are orthogonal to

**q**

**q**_{0}at each

*s*. Because is a basis for the orthogonal complement of

**q**_{0}(

*s*) in , we define at each

*s*a projectionand letfor some . Because projected constraints will appear throughout §3, we defineand for

*j*=1, 2, 3, we letthe

*j*th column of

*. By direct computation, we may show that the third row of*

**M***vanishes.*

**M**Because of the conditions (3.1) imposed on *δ q*, we have the following conditions imposed on

*:(3.2)*

**ζ**### (a) The second variation

Next, we consider the second variation of *E*, which is obtained by substituting * q*=

**q**_{0}+

*ϵδ*into

**q***E*[

*] and computing the*

**q***ϵ*

^{2}term in the Taylor expansion(3.3)where . Inserting

*δ*=

**q***into equation (3.3), and integrating by parts terms starting with*

**Πζ***′(*

**ζ***s*)

^{T}, we findwhere is the Sturm–Liouville operatorwith coefficient matrices

Note that the integral operator that defines *δ r*,is a linear functional of

*, which we denote by Thus, the second variation can be rewritten as*

**ζ**### (b) An alternative form of the second variation

The adjoint operator ^{T} has an inconvenient form that does not readily yield a computational implementation via an IVP, so we now transform to a slightly different second variation. First, we note that since the function *V*(* r*) depends only on

*x*, the only non-zero entry in is the upper-left entry

*V*

_{xx}

^{0}. We define a new linear operator byand observe that, for those

*that satisfy equation (3.2),Therefore, since the third column of is zeroThenFinally, note thatTherefore,which means that we may rewrite the second variation asThus, we have a problem in the form of a classic isoperimetrically constrained calculus of variations problem, but with the integrodifferential second-variation operatorThe theory presented in Manning*

**ζ***et al*. (1998) is an operator-based derivation of the stability index for isoperimetrically constrained calculus of variations problems, and we now show that this theory may readily be extended to .

### (c) Properties of the operator

We begin by defining some notation. For any *σ*≤1, we let (*σ*) denote the Sobolev space of -functions with square-integrable weak second derivatives on the interval (0,*σ*), andLet *L*^{2}(0,*σ*) denote the space of all square-integrable -functions on (0,*σ*). The operator has a natural restriction to *L*^{2}(0,*σ*), namelyFirst, we determine ^{T} on *L*^{2}(0,*σ*). For arbitrary functions * f*,

*∈*

**g***L*

^{2}(0,

*σ*), we computeTherefore, the adjoint operator isand it has a domain of definition

*L*

^{2}(0,

*σ*).

*The operator* *is self-adjoint*.

As is well known, the Sturm–Liouville operator is self-adjoint on _{d}(*σ*). Furthermore, because multiplication by is a self-adjoint operator on *L*^{2}(0,*σ*), we can see that is also a self-adjoint operator *L*^{2}(0,*σ*), which implies that is self-adjoint on _{d}(*σ*). ▪

*The spectrum of* *on* _{d}(*σ*) *consists of isolated eigenvalues ρ*_{1}(*σ*)≤*ρ*_{2}(*σ*)≤… *each with finite* (*geometric*) *multiplicity*.

(Follows closely the proof of property 1 in Manning *et al*. (1998)): Equivalently, we show that has no essential spectrum.2 Assume *ρ* is in the essential spectrum of . Then, by definition (Birman & Solomajak 1987, p. 206), there exists a sequence with

inf

_{j}‖*ζ*_{j}‖>0,(weak convergence of

**ζ**_{j}to**0**),(strong convergence to

**0**in*L*^{2}(0,*σ*)).

The operator is compact (Naylor & Sell 1982, p. 385), and, therefore, so is ^{T} (Kato 1980, p. 159). Thus, is a compact operator (Kato 1980, p. 158), and therefore transforms the weakly convergent sequence **ζ**_{j} into a strongly convergent sequence (Birman & Solomajak 1987, p. 39), that is, . Applying this to (iii), we see that , which, combined with (i) and (ii), shows that has *ρ* as part of its essential spectrum. However, under the assumption that * P* is positive definite (which is true for our choice of

*L*), a Sturm–Liouville operator has no essential spectrum (e.g. Renardy & Rogers 1993). Thus, we arrive at a contradiction and conclude that the essential spectrum of is empty. ▪

We denote the *L*^{2}(0,*σ*)-orthonormal eigenfunctions of by , that isNext, for *m*≥1, define the following spaces:

*Each eigenvalue ρ*_{j}(*σ*) *is a monotonically decreasing function of σ*.

(Closely follows the proof of property 2 in Manning *et al*. 1998). The operator is self-adjoint, bounded below and satisfies lemma 3.2. Further, _{d}(*σ*) is dense in *L*^{2}(0,*σ*). Therefore, the eigenvalue *ρ*_{j}(*σ*) are given by a varitional principle (Weinstein & Stenger 1972, p. 6)(3.4)Consider *σ*_{1}<*σ*_{2}. Let the functions *ψ*_{i}(*s*), *i*=1, …, *m* be defined on 0≤*s*≤*σ*_{2} byNow, determine *a*_{k} so that and such that is orthogonal to every element of span , and hence an element of _{m}(σ_{2}); because we have *m*−1 orthogonality conditions to satisfy and *m* coefficients *a*_{k} at our disposal, this can always be done by solving an (*m*−1)-by-*m* linear system and rescaling. Inserting * ζ*(

*s*) into equation (3.4), we haveFinally, we exploit the fact that {

**ζ**_{i}(

*s*;

*σ*

_{1})} is a set of orthonormal eigenfunctions of to find ▪

*For σ sufficiently close to* 0, *ρ*_{j}(*σ*)>0, ∀*j*.

(Closely follows the proof of property 3 in Manning *et al*. 1998). Standard unconstrained quadratic form theory for Sturm–Liouville operators (e.g. Schulman 1981) shows that for *σ* sufficiently smallThen, for *σ* sufficiently small, since is positive semi-definite,For any *j*, we know that *ρ*_{j} satisfies for . Premultiplying this expression by **ζ**_{j}^{T} and integrating from 0 to *σ*, we see that for sufficiently small *σ*, *ρ*_{j}(*σ*)>0. ▪

### (d) Handling the isoperimetric constraints

We have shown thatand that the space of allowed variations is , whereFurther, we have proven several facts about the spectrum of in lemmas 3.1–3.4.

A definition of conjugate points appropriate to this isoperimetrically constrained setting is given in Manning *et al*. (1998). Because the space of allowed variations _{d}^{cons}(1) is not dense in the range space *L*^{2}(0,1) of , we introduce a new operator , where is the (self-adjoint) orthogonal projection from *L*^{2}(0,1) to its subspace ^{⊥}. The range space of is ^{⊥} and _{d}^{cons}(1) is dense in ^{⊥}. Then, in analogy with the notation of §3*c*, we define the spacesandIn Manning *et al*. (1998), it is then shown that satisfies the corresponding versions of lemmas 3.1–3.4. These lemmas then show that the index equals the number of solutions *σ*<1 toso that this equation is the appropriate definition of conjugate point for this setting. Written out explicitly, without the projection operator, the definition of conjugate point is(3.5)for some constants , .

### (e) Computing conjugate points

We now show that the equation may be written as a first-order system of ODEs. We observe first thatand thusgiven that *M*_{31}=0. Therefore, may be rewritten asor, defining and , as the first-order system(3.6)Thus, analogously to the presentation in Manning *et al*. (1998), the general solution to , * ζ*(0)=

**0**may be written as , where {

**ζ**_{1},

**ζ**_{2},

**ζ**_{3}} are a basis of solutions to

*=*

**ζ****0**,

*(0)=*

**ζ****0**and for

*j*=1, 2, is a solution to

*=*

**ζ**

**T**_{j},

*(0)=*

**ζ****0**. Each of these functions may be computed as a solution to an IVP. For

**ζ**_{j}, we solve equation (3.6) with and the initial conditionswhile for we solve equation (3.6) with and the initial conditionsConjugate points occur when, at some

*σ*<1, a non-trivial member of this general solution vanishes and obeys the linearized integral constraints. This combination of conditions may be written in compact form aswhereThe integrals in

*and may similarly be computed via an IVP, for example,*

**F**### (f) The two-dimensional case

Recall that, for the two-dimensional case, there is one unknown *θ*(*s*). The allowed variations *ζ* of *θ* obey either the linearized clamped–clamped conditionsor the linearized pinned–pinned conditionsThe projection required with the Euler parameters * q* is no longer required.

The second variation again takes the formbut now with andIn the clamped–clamped case, the analysis of the second-variation operator proceeds exactly as in the three-dimensional case, as does the incorporation of the isoperimetric constraints, leading to the conjugate point equationThe conjugate points are determined by an analogous (but simpler) solution of IVP and determinant condition.

For the pinned–pinned case, a technical complication arises, because the simple extension idea used to construct *ψ* introduces a discontinuity that makes *ψ* an illegal test function. Therefore, eigenvalues *ρ*_{m}(*σ*) can potentially be increasing when they cross zero at a conjugate point, thus eliminating the simplistic notion that the number of conjugate points equals the number of negative eigenvalues. This technical computation is overcome by numerically determining d*ρ*_{m}(*σ*)/d*σ* at each conjugate point, which allows an accounting of whether an eigenvalue is becoming negative or positive at each conjugate point. For an example of how this computation can be implemented, and further discussion of this topic, see Manning (submitted).

## 4. Verification of the index theory

We used Auto to compute bifurcation diagrams of equilibrium rod configurations for the three-dimensional problem, and we determined the stability index for each equilibrium by the procedure described in §3. The bifurcation diagrams plot the load *λ* on the vertical axis and the compression 1−*z*(1) on the horizontal axis. The line styles represent the index: solid black is stable (index 0), followed by solid grey (index 1), solid light grey (index 2), dashed black (index 3), dashed grey (index 4) and dashed light grey (index 5). One such bifurcation diagram is shown in figure 2. Before analysing this diagram and how its structure depends on the rod and wall parameters, we first observe that the general properties of the index on this diagram serve to verify the theory from §3. Specifically, according to general bifurcation theory, changes in stability should generically occur exactly at bifurcation points (BPs) and folds, and this is the case in figure 2. A fold occurs where a branch is horizontal; that is, *λ* changes from increasing to decreasing, or vice versa. A bifurcation point is a point of intersection of two branches on the diagram. Owing to the symmetry of the circular rod, all of the branches in figure 2 (except for the unbuckled branch) are double-covered, so all of the labelled BPs are standard pitchfork bifurcations. Usually, if we follow a branch as it bypasses a bifurcation point and continue to increase the compression, then the rod will lose stability (cf. the two BPs on branch *y* in figure 2), although in some instances, we did observe a bifurcation point with the opposite behaviour, including one case where a rod regained stability at a bifurcation point (see figure 12). Note that in figure 12, as in all our future bifurcation diagrams, we do not show all bifurcating branches, but their existence can be inferred from index changes that occur away from folds.

Note also that no stable branches are shown above a certain maximal value of *λ*. Our model does not include any self-contact forces for the rod itself, and for sufficiently high forces, the stable configurations involve rod self-contact. Our focus is on configurations at lower forces before self-contact is relevant, particularly the effect of the wall in stabilizing certain configurations. In the two-dimensional problem (with the walls sufficiently close together), wall contact occurs before self-contact can take place, and we will see a wider range of *λ* values with stable solutions, given that in the two-dimensional setting, the rod can not bend or slide parallel to the walls to reach a self-contact configuration.

## 5. Results

We first present results for the two-dimensional pinned–pinned rod, then for the two-dimensional clamped–clamped rod and, finally, for the three-dimensional clamped–clamped rod (with the circular and elliptical cross-section cases presented separately).

### (a) Two-dimensional pinned–pinned rods

In figure 3, we show three bifurcation diagrams for the two-dimensional pinned–pinned problem, two for a soft wall with *b*=10^{−4} and 10^{−5} and one for a hard wall (the hard-wall data was traced from Holmes *et al*. (1999) using the image-processing package Sigmascan).

As *b*→0, the soft-wall diagrams appear to approach the hard-wall diagram (ignoring the extra branch in the soft-wall diagram not computed in the hard-wall case). Therefore, it seems plausible to use the stability results for the soft-wall case (with *b*≈0) to impute stability to the hard-wall case. For example, in figure 4, we show the soft-wall diagram with *b*=10^{−5} with stability indicated by line style. For *λ*<*π*^{2}, the straight configuration (zero compression) is stable, as expected from standard rod theory. At *λ*=*π*^{2}, this configuration becomes unstable in favour of a branch of stable configurations bent towards either wall. As *λ* grows, the configuration becomes more mashed up against the wall (in the hard-wall limit, a line of contact emerges along this branch (Holmes *et al*. 1999), but in the soft-wall case, this behaviour is smoothed somewhat). At *λ*≈15*π*^{2} (around the point at which the portion of the rod near the wall begins to buckle away from the wall), this branch folds and becomes unstable. The second bifurcating branch off the zero-compression branch contains configurations with two turning points. These are initially unstable, but contact with the wall stabilizes them when the compression reaches approximately 0.023.

These results confirm the experimental observations reported by Holmes *et al*. (1999). As shown in the dashed branch in figure 4, as a polypropylene fibre was compressed, it initially buckled towards one of the walls and then switched to a configuration with two turning points. Then, as shown in the dotted branch, upon relaxation, this two-turning-point solution was retained past the point at which it originally emerged, rejoining the single-turning-point solution via the roughly horizontal curve. The stability results in figure 4 are consistent with this sort of hysteresis loop; once *λ* exceeds 15*π*^{2}, there is no more one-turning-point stable solution, so we would expect the rod to ‘snap’ to the two-turning-point branch. Similarly, upon relaxation, we would expect to remain on the two-turning-point branch until *λ* drops below 4*π*^{2}, when it would snap to the one-turning-point branch.

Clearly, there are limits to using this idealized two-dimensional rod as a model for the polypropylene fibre. Upon compression, instead of the hysteretic horizontal jump that the theory would predict, we have a smoother transition via a sloping branch that furthermore reaches higher compression than the theoretical two-turning-point branch. Similarly, upon relaxation, instead of a horizontal jump from two turning points to one, Holmes *et al*. report that the polypropylene fibre passes through an asymmetric two-turning-point state. In fact, there is a branch of such asymmetric states not shown in figure 4, almost parallel to the branch of symmetric two-turning-point states around *λ*=4*π*^{2} (it begins at the bifurcation point suggested by the change in stability at compression 0.023). Neither the branch of symmetric nor asymmetric two-turning-point configurations is stable below compression 0.023 in the theoretical model, so it seems there is some perturbation in the physical system that favours the asymmetric states.

### (b) Two-dimensional clamped–clamped rods

In figure 5, we show three bifurcation diagrams for the two-dimensional clamped–clamped problem, for a soft wall with *b*=10^{−3}, 10^{−4} and 10^{−5}. Once again, as *b*→0, the soft-wall bifurcation diagrams appear to approach a limit, presumably the hard-wall diagram. In figure 6, we show the soft-wall diagram with *b*=10^{−5} with stability indicated by line style. We see a pattern reminiscent of that seen in the pinned–pinned case, but different in detail. The one-turning-point branch is initially stable but loses stability at a bifurcation point, after which solutions have three turning points. These solutions are initially unstable, but become stabilized by contact with the wall around compression 0.06. The two-turning-point branch is initially unstable, but becomes stable after a bifurcation point, only to become unstable later on (but off the edge of the diagram, it becomes a branch with stable four-turning-point configurations).

### (c) Three-dimensional clamped–clamped rods

We conclude by presenting results for three-dimensional buckling subject to clamped–clamped boundary conditions. In §5*c*(i) and (ii), respectively, we consider rods with circular and elliptical cross-sections, with spacing parameter *a*=0.1 and strength *b*=10^{−4}. Some of the most common configurations the rod assumes are portrayed in figure 7. In the bifurcation diagrams throughout §5, branches have labels *a*1, *a*2, *b*1, *b*2, *c*, *d* and *e* corresponding to the configurations in figure 7.

We note that the simple configurations shown in figure 7 all appear to display ‘flip symmetry’ about *s*=1/2, as defined by Domokos & Healey (2001), who showed that all non-contact ring equilibria have flip symmetry. However, the presence of wall contact in our problem appears to break this symmetry, as there are configurations that do not have flip symmetry. These asymmetric configurations lie on secondary or tertiary branches that for simplicity have been omitted from our bifurcation diagrams. For example, in figure 10, asymmetric configurations lie on a branch that bifurcates from the *x* branch where the index changes from 1 to 2, and also on a branch that bifurcates from the *xy* branch where the index changes from 0 to 1.

#### (i) Rods with circular cross-sections

A rod with a circular cross-section initially buckles in the plane parallel to the soft walls. Buckling perpendicular to the walls requires the rod to overcome the wall strength, however slight, while buckling in the *y*–*z*-plane is free of resistance from the wall. This is illustrated in figure 8, where branch *y* (containing configurations parallel to the walls) occurs just before branch *x* (containing configurations perpendicular to the walls).

On branch *y*, the rod initially forms a stable one-turning-point bend (index 0) parallel to the walls, but then becomes unstable at a secondary bifurcation point. Table 1 displays the value of the compression at this secondary bifurcation point for several values of *a* and *b*. The data suggest that for a wall of any strength, this loss of stability occurs for compression 0.86 or higher, at which point the rod has come in contact with itself (see configuration B2 in figure 7). By this time, a true rod would be subject to self-contact effects not accounted for here.

On branch *x*, the rod initially forms a one-turning-point bend perpendicular to the wall. These initial configurations are unstable (index 1), and after a significant increase in *λ*, there are secondary BPs, after which the rod becomes more unstable (index 2, 3, 4). Where the dashed grey segment of branch *x* turns sharply right, the rod forms a three-turning-point configuration perpendicular to the wall. On the first secondary bifurcating branch (branch *xy*), we have two-turning-point configurations with bends that slide along the wall in opposite directions as seen in configuration *E* of figure 7. This branch initially has index 1, but becomes increasingly unstable. Therefore, it seems that the only stable configurations for a rod with circular cross-section in the soft-wall model lie on the plane parallel to the walls.

If the wall force at *x*=0 is sufficiently large (e.g. for a close pair of walls *a*=0.007 5, *b*=10^{−4}), then the behaviour on branch *x* can be slightly different, in that the rod may exhibit multiple turning points when initially buckling towards the wall (data not shown). However, this can only make the perpendicular configurations more unstable.

Finally, note that figure 8 includes a portion of the bifurcation diagram in the absence of the wall, to allow comparisons of this and later diagrams with this standard. As one would expect, the no-wall diagram includes a branch of uncompressed solutions containing an infinite family of BPs. Owing to the circular symmetry of the no-wall problem, buckling can occur in two dimensions, and hence the index increases by two, at each bifurcation point. The stability of the no-wall diagram is not shown in figure 8, but by the above argument, the no-wall zero-compression branch has the same 0–2–4–… index pattern seen for the wall problem in figure 8 (although recall that, for example, the first two BPs in the wall diagram actually have a tiny separation between them, in which the index is one). The first bifurcating branch in the no-wall diagram is identical to the *y* branch in the wall diagram (with the same index pattern), not surprisingly, because the two wall forces cancel out on the *y* branch. The no-wall diagram has a secondary bifurcating branch near the wall branch *yx*, but with a different shape, which is to be expected given that these nonplanar solutions do experience the wall when it is present.

#### (ii) Rods with elliptical cross-sections

We now consider rods with elliptical cross-sections whose major axes are parallel to the walls. If the walls are relatively weak, then the first primary bifurcating branch will consist of one-turning-point bends perpendicular to the wall (branch *x* in figures 9–11) that are initially stable. For higher loads, there will also be a bifurcating branch consisting of bends parallel to the wall (branch *y* in figures 9–11), but these configurations will be unstable.

For a rod with a slightly elliptical cross-section (e.g. *K*_{1}=0.1, *K*_{2}=0.09 as in figure 9), the primary bifurcating branches *x* and *y* are in close proximity. Branch *x* almost immediately has a secondary bifurcating branch, branch *xy*, which consists of configurations that bend at an angle towards the wall. These configurations look like configuration D in figure 7, one-turning-point stable bends that exist liminally between the wall and the plane parallel to the wall. The rod seems to find a balance between the force of the wall and its slightly elliptical shape. Intriguingly, these stable bends are nearly identical in shape to the unstable bends in branch *y*. Thus, the bifurcation diagram in figure 9 displays a stable and an unstable branch that share nearly identical paths.

Looking at a more elliptical rod (e.g. *K*_{1}=0.1, *K*_{2}=0.05 as in figure 10), the secondary bifurcation point on branch *x* occurs when the rod is in near contact with the wall. This secondary bifurcating branch, branch *xy*, contains stable one-turning-point bends that are ‘folded’ along the wall as in configuration C of figure 7. After the secondary bifurcation point, branch *x* continues with unstable three-turning-point configurations.

Finally, for an even more extreme ellipse (e.g. *K*_{1}=0.1, *K*_{2}=0.05 as in figure 11), the formation of three-turning-point configurations occurs closer to emergence of the *xy* branch of folded bends. Indeed, for the appropriate stiffness parameters, the rod can form a stable three-turning-point configuration perpendicular to the wall before this *xy* branch (e.g. the black segment on branch *x* in figure 12).

#### (iii) Conjectures for the hard-wall problem

The bifurcation diagrams for walls with *b*=10^{−4} and 10^{−5} are nearly identical when *a* is relatively large and the rod cross-section sufficiently elliptical (e.g. for our diagrams with *a*=0.1 and *K*_{2}/*K*_{1}≤0.9). This indicates that the equilibrium configurations and stability analysis have probably already reached their *b*→0 limits. This provides inductive evidence that the soft-wall problem might be used as an indicator of stability in the hard-wall problem, in which a rod buckles into a pair of walls that exert no force for |*x*|<*a* and disallow configurations with |*x*|>*a*. Therefore, for figures 9–12, the bifurcation diagrams for the hard wall should closely resemble those we show for the soft wall.

In cases when the walls are very close, or when the elliptical rod cross-section is nearly circular, the asymptotic approach of the soft-wall results to the hard-wall results will generally be slower. For example, for a rod with elliptical cross-section with major axis parallel to the hard walls, buckling perpendicular to the walls will occur before buckling parallel to the walls. However, if the soft walls are close to the rod, then the wall strength will tend to reverse this ordering, so that we only obtain the correct ordering for very small *b*. Similarly, when the cross-section is nearly circular, the perpendicular branch occurs only slightly before the parallel branch in the hard-wall problem, and so it takes only a very small soft-wall force to reverse the ordering. Nevertheless, even in these cases, one should still be able to extract a good approximation to the hard-wall diagram by taking *b* sufficiently small in the soft-wall problem.

In contrast, the case of a circular rod is somewhat special. For example, when bordered by hard walls, a rod with a circular cross-section is equally likely to buckle in every direction, while when bordered by soft walls, it can only buckle perpendicular to or parallel to the wall, no matter how small *b* is. Thus, the circular-rod hard-wall solution set is harder to extract as a simple limit of soft-wall results.

A few conclusions are clear. As in the soft-wall model, a rod that bends parallel to a hard wall will form a stable one-turning-point configuration until after it comes in contact with itself, as discussed in §5*c*(i). If the rod buckles towards the wall, then it would initially form a one-turning-point configuration, but unlike similar configurations in the soft-wall model, these configurations would initially be stable (by the circular symmetry of the problem before the wall contact begins). Just as in the soft-wall case, we then expect a secondary bifurcating branch like our soft-wall *xy* branches (though it is less clear whether configurations will be more like C or D—configurations such as D can certainly arise if the initial rod buckling is diagonal to the wall, but their appearance in figure 9 appears to be a result of a balance between the elliptic cross-section and the soft-wall force, an effect that would be absent in the circular-rod hard-wall case). In addition, if we continue past this secondary bifurcation point, then the rod should form unstable configurations with three turning points like configuration A2 (although, if the walls are sufficiently close, these would presumably emerge as stable solutions before the secondary bifurcation, as we have seen in our soft-wall simulations).

## Acknowledgments

R.S.M. was supported by NSF grants DMS-9973258 and DMS-0384739.

## Footnotes

As this paper exceeds the maximum length normally permitted, the authors have agreed to contribute to production costs.

↵In the numerical computations, the boundary condition

*q*_{4}(0)=1 was replaced by the condition*μ*_{4}(0)=0. This change removes a gauge symmetry of the problem that results from the representation of elements (**d**_{1},**d**_{2},**d**_{3}) of the three-dimensional group*SO*(3) by∈R**q**^{4}. Without this change, there would be an entire line of solutions to the BVP at each value of*λ*, which is undesirable for numerical computations. The condition*q*_{4}(0)=1 can safely be dropped, since the quantity || is conserved by the differential equations, so the remaining boundary conditions on**q**will guarantee that**q***q*_{4}(0)=1 without it being explicitly enforced in the computations.↵There is some variation in how the essential spectrum is defined. Here, we take the definition set forth in Birman & Solomajak (1987).

- Received August 20, 2004.
- Accepted January 18, 2005.

- © 2005 The Royal Society