ON THE OPTIMAL DESIGN OF COLUMNS AGAINST BUCKLING

Steven J. Cox† and Michael L. Overton‡

Final Version for SIAM J. Math. Anal.

† Mathematical Sciences Department, Rice University, Houston, TX 77251. This author’s research was supported by Defense Advanced Research Projects Agency Contract F49620-87-C-0065 and Air Force Office of Scientific Research Grant 89-0363.

‡ Courant Institute of Mathematical Sciences, New York University, New York, NY 10012. This author’s research was supported by National Science Foundation Grant CCR-88-02408.

Abstract – We establish existence, derive necessary conditions, and construct and test an algorithm for the maximization of a column’s Euler buckling load under a variety of boundary conditions over a general class of admissible designs. We prove that symmetric clamped–clamped columns possess a positive first eigenfunction and introduce a symmetric rearrangement that does not decrease the column’s buckling load. Our necessary conditions, expressed in the language of Clarke’s generalized gradient [10], subsume those proposed by Olhoff and Rasmussen [25], Masur [22], and Seiranian [34]. The work of [25], [22], and [34] sought to correct the necessary conditions of Tadjbakhsh and Keller [37] who had not foreseen the presence of a multiple least eigenvalue. This remedy has been hampered by Tadjbakhsh and Keller’s miscalculation of the buckling loads of their clamped-clamped and clamped–hinged columns. We resolve this issue in the appendix.

In our numerical treatment of the associated finite dimensional optimization problem we build on the work of Overton [26] in devising an efficient means of extracting an ascent direction from the column’s least eigenvalue. Owing to its possible multiplicity this is indeed a nonsmooth problem and again the ideas of Clarke [10] are exploited.

1. Introduction

We recall Todhunter’s formulation [38, p. 66] of the following problem of Lagrange,

“To find the curve which by its revolution about an axis in its plane determines the column of greatest efficiency.”

For columns of unit length and volume, efficiency here denotes the structure’s resistance to buckling under axial compression. When λ is the magnitude of the axial load and u the resulting transverse displacement we postulate the potential energy

∫_{0}^{1} EI|u′′|^{2} dx − λ ∫_{0}^{1} |u′|^{2} dx

with the two terms measuring bending and elongation respectively. Here I is the second moment of area of the column’s cross section and E is its Young’s modulus. For sufficiently small λ, the minimum of this potential energy, over all admissible displacements, is zero.

The (Euler) buckling load of the column is the greatest λ, call it λ_{1}, for which this minimum is zero. That is,

λ_{1} = inf_{u∈V} \frac{∫_{0}^{1} EI|u′′|^{2} dx}{∫_{0}^{1} |u′|^{2} dx}  (1.1)

where V is a closed subspace of H^{2}, the space of L^{2} functions on the interval (0, 1) with first and second distributional derivatives in L^{2}, distinguished by the choice of boundary conditions. The choice that has generated the greatest interest is the clamped-clamped condition u(0) = u′(0) = u(1) = u′(1) = 0. With the corresponding V denoted by H_{0}^{2} it is not difficult to show that the infimum in (1.1) is attained at some u_{1} ∈ H_{0}^{2}. First order necessary conditions then require that u_{1} satisfy

∫_{0}^{1} EIu′′_{1}v′′ dx = λ_{1} ∫_{0}^{1} u′_{1}v′ dx,  ∀ v ∈ H_{0}^{2}.  (1.2)

When I and E are smooth it follows from (1.2) that

−(EIu′′_{1})′′ = λ_{1}u′′_{1}, u_{1}(0) = u′_{1}(0) = u_{1}(1) = u′_{1}(1) = 0.  (1.3)

With this we recognize (1.1) as Rayleigh’s principle for the least eigenvalue of (1.3) and u_{1} as an associated first eigenfunction. For the problem of Lagrange the Young’s modulus is assumed constant and, as the column is a solid of revolution, each cross section’s second moment of area is simply a constant multiple of the square of its area, A, i.e., I(x) = cA^{2}(x). Fixing our attention on columns of unit volume, we require

∫_{0}^{1} A dx = 1.  (1.4)

We have reduced the problem of Lagrange to the search for that A which, subject to (1.4), maximizes the λ_{1} of (1.1). This problem, with clamped-clamped boundary conditions, was first attacked in 1962 by Tadjbakhsh and Keller [37] in the continuation of work Keller [20] had begun at the suggestion of Clifford Truesdell. The work of [37] hinges on the necessary condition that the best A, and its corresponding eigenfunction u, satisfy

A^{4}|u′′|^{2} = A^{3}  (1.5)

along the entire column. This was obtained on formally differentiating a second order analog (see equation (2.5)) of (1.3) with respect to A subject to the integral constraint. Upon reconciling (1.3) and (1.5), Tadjbakhsh and Keller arrived at the representation

A(x) = \frac{4}{3} sin^{2} θ(x), −π/2 ≤ θ ≤ 3π/2, (1.6)

θ(x) − \frac{1}{2} sin 2θ(x) + π/2 = 2πx, 0 ≤ x ≤ 1.

The most striking aspect of this claim is that it requires the cross sectional area to vanish at 1/4 and 3/4. This result should however come as no surprise, for implicit in (1.5) is the assumption that the optimal buckling load is simple, i.e., that the corresponding space of buckled configurations is one dimensional. This requires the optimal column to buckle in much the same way as the uniform column (A ≡ 1), the first eigenfunction of which is U(x) = 1 − cos(2πx). The fact that A vanishes at the inflection points of U agrees then with the heuristic (suggested by (1.1)) that the optimal column need be thick only in regions where it bends, i.e., where the magnitude of the linearized curvature |u′′| is large.

Tadjbakhsh and Keller claimed 16π^{2}/3 as the buckling load of the resulting column. It was not until 1977 that Olhoff and Rasmussen [25], observing that (1.3) does not exclude multiple eigenvalues, noted that as the least eigenvalue does not vary smoothly with A at points where its multiplicity exceeds one, the formal differentiation in [37] would be hard to justify. As evidence that Tadjbakhsh and Keller had indeed taken the wrong course, Olhoff and Rasmussen claimed, on the basis of numerical work, 30.51 for the buckling load of the column constructed according to (1.6). Unfortunately, they neglected to describe the means by which this value was arrived at. Indeed the fact that A vanishes at 1/4 and 3/4 introduces computational difficulties. Although they did go on to suggest how 16π^{2}/3 was incorrectly obtained, a number of workers have remained unconvinced, e.g., Myers and Spillers [24] and Barnes [4]. Upon fleshing out the relevant remarks of Olhoff and Rasmussen we shall see, in work relegated to an appendix, that the buckling load for the column proposed by Tadjbakhsh and Keller does not exceed π^{2}/3. These same arguments will serve to demonstrate that Tadjbakhsh and Keller’s best clamped-hinged column also has a much lower buckling load than thought previously.

Having concluded that differentiating (1.3) would lead to less than optimal columns, Olhoff and Rasmussen presented a ‘bimodal formulation’ of the problem of Lagrange, i.e., one which would accommodate double eigenvalues. Their subsequent necessary condition paired the best A with two corresponding linearly independent eigenfunctions u, v and a scalar t ∈ [0, 1] so that

A (t|u′′|^{2} + (1 − t)|v′′|^{2}) = 1,  (1.7)

along the entire column. On implementing an algorithm that enforced this optimality condition, Olhoff and Rasmussen arrived at a column whose cross sectional area was positive throughout and which could withstand loads up to 52.3563. Their methods were however no more rigorous than those of Tadjbakhsh and Keller, and moreover, solely on the basis of claims, the latter still had the stronger column, for 52.3563 < 16π^{2}/3. Those persuaded by Olhoff and Rasmussen’s criticism of the work of Tadjbakhsh and Keller then set out to put (1.7) on a solid foundation. Actually, they joined the discussion of the more general problem: What conditions are necessary for a multiple eigenvalue to attain its extremum? The greatest advances on this question have come in finite dimensions and lie in the apparently little known work of Bratus and Seiranian [6]. These conditions, later discovered independently in a more general form by Overton [26], will be discussed in detail in §5. For now, we note that Bratus and Seiranian, upon applying their finite dimensional arguments to the problem of Lagrange arrived at the conclusion that the best A, must, with two corresponding orthogonal eigenfunctions u, v, satisfy

A(δ_{1}|u′′|^{2} + δ_{2}|v′′|^{2} + δ_{3}u′′v′′) = 1, where δ_{1}δ_{2} ≥ δ_{3}^{2}/4.  (1.8)

This condition was also proposed by Masur [22] who, like Seiranian [34], went on to represent the best A via a system of transcendental equations. Their approximate solutions to these systems are in good agreement, with respective buckling loads of 52.3564 and 52.3565, with that proposed by Olhoff and Rasmussen [25]. Note that (1.7) and (1.8), with the introduction of a second buckling mode, possess mechanisms which, at least in principle, rule out the possibility of columns with vanishing cross sectional area.

Our main contribution to the problem of Lagrange is essentially twofold. We employ the generalized gradient of Clarke in (i) a rigorous derivation of the necessary conditions (1.8) and (ii) the construction of an efficient algorithm to solve the associated finite dimensional optimization problem. Our initial focus on the clamped-clamped case will be extended in §5 to each of the boundary conditions considered by Tadjbakhsh and Keller.

In our discussion of the various optimality criteria something has been conspicuously lacking: the literature contains no proof of the existence of a best A for the problem of Lagrange. Before filling this gap we establish a number of preliminary results and look to a more general problem formulation.

2. The Optimal Design Problem

The moment I is more precisely the second moment of area of the cross section about a line through its centroid normal to the plane of buckling. That is, denoting the cross section by Ω(x) with centroid at the origin, if η is a unit normal to the plane of buckling then

I(x) = ∫_{Ω(x)} |η^{T} y|^{2} dy.  (2.1)

When Ω is a circle, in fact when Ω is a regular polygon, this integral does not depend on η, and one finds that I varies as the square of the cross-sectional area, A. On considering so–called thin-walled columns we shall now see that I varies as an affine function of A. On the lateral surface of a cylinder with circular cross section of constant radius R we add a layer of variable thickness ρ(x) with ρ(x) ≤ εR, ε ≪ 1. Neglecting powers of ρ greater than one we find I(x) = πR^{3}ρ(x) + πR^{4} and A(x) = 2πRρ(x) + πR^{2}. Taking Ã(x) = A(x)−πR^{2}/2 for our design variable we find I(x) = (R^{2}/2)Ã(x). The effect of this choice on the integral constraint is trivial. Of greater interest is that Ã, by construction, must satisfy the pointwise bounds

πR^{2}/2 ≤ Ã(x) ≤ πR^{2}/2 + 2επR^{2}.  (2.2)

It is not difficult to continue this line of reasoning and collect a number of examples where I varies as some power of A. We proceed then to consider the case where EI = σ^{p} for some p > 0. Compelled by our examination of the previous special cases we admit those σ in

ad = {σ ∈ L^{∞} : 0 < α ≤ σ(x) ≤ β, ∫_{0}^{1} σ(x) dx = 1}.

The weak formulation of the buckled column equation for σ ∈ ad is

∫_{0}^{1} σ^{p}u′′v′′ dx = λ ∫_{0}^{1} u′v′ dx,  ∀ v ∈ H_{0}^{2}.  (2.3)

As σ ∈ L^{∞} and α > 0, (2.3) possesses the sequence of eigenvalues

0 < λ_{1}(σ) ≤ λ_{2}(σ) ≤ ⋯ ↑ ∞,

repeated according to their finite multiplicities and a corresponding sequence of eigenfunctions {u_{k}(σ)}_{k=1}^{∞} ⊂ H_{0}^{2}, orthonormal in terms of the bilinear form associated with the right side of (2.3). As H_{0}^{2}(0, 1) ⊂ C^{1}([0, 1]) we find u_{k} ∈ C^{1}([0, 1]). Upon integrating by parts on the right side of (2.3) we find that σ^{p}u′′_{k} differs from −λ_{k}(σ)u_{k} by an affine function of x. Hence, σ^{p}u′′_{k} ∈ C^{1}([0, 1]), and, in fact

(σ^{p}u′′_{k})(x) = (σ^{p}u′′_{k})′(0)x + (σ^{p}u′′_{k})(0) − λ_{k}(σ)u_{k}(x),  (2.4)

We collect those eigenfunctions corresponding to λ_{1}(σ) in

E(σ) = span {u_{k}(σ); λ_{k}(σ) = λ_{1}(σ)},

a subspace of H_{0}^{2} with dimension equal to the multiplicity of λ_{1}(σ). Implicit in Olhoff and Rasmussen’s bimodal formulation is the assumption that this multiplicity is at most two. Seiranian [34], has confirmed this through Kamke’s analysis of the second order problem with nonseparated boundary conditions

w′′ + λσ^{−p}w = 0, w(1) = w(0) + w′(0), w′(1) = w′(0).  (2.5)

This is the strong version of (2.3) with w = σ^{p}u′′ and was first considered in our context by Tadjbakhsh and Keller. Kamke, in [19, §4], proves that the multiplicity of each eigenvalue of (2.5) is no greater than two. Equation (2.4) however, suggests an approach that applies directly to the weak formulation.

If corresponding to λ_{k}(σ) there existed three linearly independent eigenfunctions u_{1}, u_{2}, u_{3} then one could choose scalars a, b, c not all zero such that v = au_{1} + bu_{2} + cu_{3} satisfies, in addition to v(0) = v′(0) = v(1) = v′(1) = 0, the two conditions (σ^{p}v′′)′(0) = 0 and (σ^{p}v′′)(0) = 0. From (2.4) we conclude that v satisfies the homogeneous linear second order equation with zero initial conditions

σ^{p}(x)v′′(x) + λ_{k}(σ)v(x) = 0, v(0) = v′(0) = 0.

As the only solution to this equation is the identically zero function we have established

Lemma 2.1. If σ ∈ ad then the multiplicity of λ_{k}(σ) is at most two.

As the least eigenvalue of the uniform column is 4π^{2} we find, as a consequence of the monotonicity of the Rayleigh quotient, that

4π^{2}α^{p} ≤ λ_{1}(σ) ≤ 4π^{2}β^{p},  ∀ σ ∈ ad.  (2.6)

Corresponding to the least eigenvalue λ_{1}(σ) one expects a positive eigenfunction. Indeed, this is the only type that Tadjbakhsh and Keller expected. To our knowledge however, there is no proof that a positive first eigenfunction need exist. We remark that on this point the oscillation theory of Kamke is insufficient, for it concludes only that eigenfunctions corresponding to the least nonzero eigenvalue of (2.5) possess either three or two zeros. This translates into either one or no zero(s) for eigenfunctions corresponding to λ_{1}(σ). We now improve on this situation in the case where σ is even (about 1/2), i.e., σ(x) = σ(1 − x).

Theorem 2.2. If σ ∈ L^{∞} is even and admits a positive lower bound then there exists a positive even eigenfunction corresponding to λ_{1}(σ).

Proof: We exploit the essential idea in inverse iteration, a popular technique for computing the least eigenvalue and eigenvector of a symmetric matrix. In our context this idea amounts to approximating the least eigenfunction by the solution of a related nonhomogeneous boundary value problem. Given v_{0} ∈ H_{0}^{2} we consider its expansion in the complete set of eigenfunctions {u_{k}(σ)},

v_{0}(x) = \overline{v}(x) + ∑_{k=m+1}^{∞} a_{k}u_{k}(x),

where m is the least integer for which λ_{m}(σ) < λ_{m+1}(σ) and \overline{v} is an eigenfunction corresponding to λ_{1}(σ). From v_{0} we construct the sequence {v_{j}} ⊂ H_{0}^{2} according to

∫_{0}^{1} σ^{p}v′′_{j}ϕ′′ dx = λ_{1}(σ) ∫_{0}^{1} v′_{j−1}ϕ′ dx, ∀ ϕ ∈ H_{0}^{2}.

On expanding v_{j} in {u_{k}(σ)} one finds

v_{j}(x) = \overline{v}(x) + ∑_{k=m+1}^{∞} a_{k} (\frac{λ_{1}(σ)}{λ_{k}(σ)})^{j} u_{k}(x).

As λ_{1}(σ) < λ_{k}(σ) for all k > m we find that v_{j} converges pointwise to \overline{v} as j → ∞. It remains then to produce a v_{0} whose corresponding \overline{v} is even and positive.

Our choice for v_{0} is the first eigenfunction of the uniform column, i.e., 1 − cos(2πx), a positive even function with exactly two inflection points.

Lemma 2.3. Let f be an even member of L^{∞} with a positive lower bound and v be a positive, even member of H_{0}^{2} with precisely two inflection points. If u ∈ H_{0}^{2} satisfies

∫_{0}^{1} fu′′ϕ′′ dx = ∫_{0}^{1} v′ϕ′ dx,  ∀ ϕ ∈ H_{0}^{2}  (2.7)

then u is positive, even, and possesses precisely two inflection points.

Proof of Lemma 2.3: Upon integrating by parts on the right of (2.7) we find that fu′′ differs from v by an affine function. Dividing by f and integrating twice gives

u(x) = ∫_{0}^{x}(x − y)(ay + b − v(y))g(y) dy  (2.8)

where g = 1/f and a and b are determined by u(1) = 0 and u′(1) = 0, i.e., by

a ∫_{0}^{1} xg(x) dx + b ∫_{0}^{1} g(x) dx = ∫_{0}^{1} v(x)g(x) dx  (2.9)
a ∫_{0}^{1} x^{2}g(x) dx + b ∫_{0}^{1} xg(x) dx = ∫_{0}^{1} xv(x)g(x) dx  (2.10)

That these equations uniquely determine a and b follows from Hölder’s inequality

(∫_{0}^{1} xg(x) dx)^{2} < ∫_{0}^{1} g(x) dx ∫_{0}^{1} x^{2}g(x) dx.

Our hypotheses in fact allow us to conclude that

a = 0, b = \frac{∫_{0}^{1} v(x)g(x) dx}{∫_{0}^{1} g(x) dx}.

This obviously satisfies (2.9). Regarding (2.10), recall that every even function satisfies ∫_{0}^{1} ϕ(x) dx = 2 ∫_{0}^{1} xϕ(x) dx. Consequently,

b = \frac{∫_{0}^{1} v(x)g(x) dx}{∫_{0}^{1} g(x) dx} = \frac{2 ∫_{0}^{1} xv(x)g(x) dx}{2 ∫_{0}^{1} xg(x) dx} = \frac{∫_{0}^{1} xv(x)g(x) dx}{∫_{0}^{1} xg(x) dx}

satisfies (2.10) as well. Labeling s(x) = (b − v(x))g(x), equation (2.8), u(1) = 0, and u′(1) = 0 take the form

u(x) = ∫_{0}^{x}(x − y)s(y) dy, ∫_{0}^{1} s(y) dy = 0, ∫_{0}^{1} ys(y) = 0.

With this and the fact that s is even we find

u(1 − x) = ∫_{0}^{1−x}(1 − x − y)s(y) dy
= ∫_{0}^{1}(1 − x − y)s(y) dy − ∫_{1−x}^{1}(1 − x − y)s(y) dy
= ∫_{0}^{x}(x − y)s(y) dy = u(x).

Regarding the convexity/concavity of u we observe that f(x)u′′(x) = b − v(x). That b − v(x) has at least two zeros follows from b > 0, v(0) = v(1) = 0, and 0 < b < ∥v∥_{∞}. For b − v(x) to possess more than two zeros v must admit a local minimum, a condition that requires of v no less then four inflection points. These zeros, say x_{0} and 1 − x_{0}, are the inflection points of u. As u vanishes at 0 and is convex on (0, x_{0}) it must be positive there, and, by symmetry, positive on (1 − x_{0}, 1) as well. As u is positive at x_{0} and 1 − x_{0} while concave between these points it must be positive on this interval as well.

It now follows that {v_{j}} is a sequence of positive even functions. The convergence of v_{j} to \overline{v} being pointwise we conclude that \overline{v} is itself a positive even function.

When σ is even and λ_{1}(σ) is simple we now have, up to a scalar multiple, a unique positive even first eigenfunction, call it u_{1}. When λ_{1}(σ) is double, in addition to u_{1}, there exists a first eigenfunction u_{2} for which ∫_{0}^{1} u′_{1}u′_{2} dx = 0. Consequently, u_{2} is not even, and as u_{1} and u_{2} span E(σ) we may conclude that when σ is even, there exists, up to a scalar multiple, a unique positive even first eigenfunction.

Though Theorem 2.2 applies only to even functions we shall see in the next result that this suffices for our purposes. Note that Lemma 2.3 states that the operator (d^{2}/dx^{2}(f d^{2}/dx^{2}))^{−1}(−d^{2}/dx^{2}) leaves a subcone of the positive H_{0}^{2} functions invariant when f is even. This cone is however too ‘thin’ to allow one to invoke Krein-Rutman arguments. Regarding possible improvements of Lemma 2.3 we note that even the constant coefficient operator (d^{4}/dx^{4})^{−1}(−d^{2}/dx^{2}) does not leave the positive H_{0}^{2} functions invariant. To see this we solve for b in (2.9-10) with g = 1,

b = 4 ∫_{0}^{1} v dx − 6 ∫_{0}^{1} xv dx.

Taking for v any smooth positive function supported in (2/3, 1) produces b < 0. As u(0) = u′(0) = 0 and u′′(0) = b we conclude that u is not positive.

Theorem 2.4. Given σ ∈ ad there exists an even σ_{∗} ∈ ad for which λ_{1}(σ) ≤ λ_{1}(σ_{∗}).

Proof: There is a very simple argument when 0 < p ≤ 1. Given a function ϕ on (0, 1) we denote its even part by ϕ_{s}(x) = \frac{1}{2}(ϕ(x) + ϕ(1 − x)). Consider the even function σ̃ ≡ ((σ^{p})_{s})^{1/p} and its corresponding even first eigenfunction ũ. With the normalization ∥ũ′∥ = 1 we find

λ_{1}(σ) ≤ ∫_{0}^{1} σ^{p}|ũ′′|^{2} dx = ∫_{0}^{1}(σ^{p})_{s}|ũ′′|^{2} dx = ∫_{0}^{1} σ̃^{p}|ũ′′|^{2} dx = λ_{1}(σ̃).  (2.11)

As t ↦ t^{1/p} is convex we observe that

σ̃(x) = (\frac{σ^{p}(x)}{2} + \frac{σ^{p}(1 − x)}{2})^{1/p} ≤ \frac{σ(x)}{2} + \frac{σ(1 − x)}{2} = σ_{s}(x).  (2.12)

Now σ_{s} ∈ ad and (2.11-12) imply that λ_{1}(σ) ≤ λ_{1}(σ_{s}). Our attempts to argue in a similar fashion for p > 1 with σ̃ ≡ ((σ^{−p})_{s})^{−1/p} and (2.5) have been thwarted by the fact that λ_{1}(σ) corresponds to the third eigenvalue of (2.5). What is needed is a rearrangement of σ that echoes the curvature of its corresponding first eigenfunction. To make this precise we first need the following extension of Lemma 2.3.

Recall that a function ϕ is said to be odd about the point (x_{0}, ϕ(x_{0})) on some interval containing x_{0} when

ϕ(x_{0}) − ϕ(x_{0} − x) = ϕ(x_{0} − x) − ϕ(x_{0})

for each x on the given interval. If, in addition to the original hypotheses of Lemma 2.3, we assume that f and v when restricted to (0, 1/2) are even about 1/4 and odd about (1/4, v(1/4)) respectively, we conclude that u, when restricted to (0, 1/2), is odd about (1/4, u(1/4)).

To see this we recall that ∫_{0}^{1} u′′ dx = 0 and u′′ is even about 1/2, hence ∫_{0}^{1/2} u′′ dx = 0. For the remainder of this paragraph all functions will be restricted to (0, 1/2). Recall as well that u′′ = (b − v)/f, the quotient of a function odd about (1/4, b − v(1/4)) and a function even about 1/4. Hence u′′ is odd about (1/4, b − v(1/4)). The condition that ∫_{0}^{1/2} u′′ dx = 0 now forces b = v(1/4). As u′′ is now odd about (1/4, 0) and u(0) = u′(0) = 0 we easily conclude that u is indeed odd about (1/4, u(1/4)).

If σ_{∗} is now even about 1/2 and even about 1/4 when restricted to (0, 1/2) then beginning the iteration of Theorem 2.2 with a v_{0} that is even about 1/2 and odd about (1/4, v_{0}(1/4)), e.g., 1 − cos(2πx), will produce u_{∗}, a positive eigenfunction corresponding to λ_{1}(σ_{∗}) that is even about 1/2 and odd about (1/4, u_{∗}(1/4)) on (0, 1/2). We immediately note that σ_{∗}^{p} and |u′′_{∗}|^{2} are similarly ordered, i.e.,

(σ_{∗}^{p}(x_{1}) − σ_{∗}^{p}(x_{2}))(|u′′_{∗}(x_{1})|^{2} − |u′′_{∗}(x_{2})|^{2}) ≥ 0  ∀ x_{1}, x_{2} ∈ (0, 1).  (2.13)

Given σ ∈ ad we now define its rearrangement σ_{∗}.

ℓ_{c} = {x ∈ (0, 1) : σ(x) ≥ c}
ℓ_{c}^{∗} = \{ \begin{matrix}\{x ∈ R : |x − 1/2| ≤ 1/4|ℓ_{c}|\} & if ℓ_{c} ≠ ∅ \\ ∅ & otherwise \end{matrix}

σ_{∗}(x) = \{ \begin{matrix}σ_{∗}(1/2 − x), & if 0 ≤ x ≤ 1/4 \\ sup \{c ∈ R : x ∈ ℓ_{c}^{∗}\}, & if 1/4 ≤ x ≤ 3/4 \\ σ_{∗}(1 − x), & if 3/4 ≤ x ≤ 1. \end{matrix}

In essence, this distributes half of σ’s mass in a symmetrically decreasing fashion about 1/2 on (1/4, 3/4), completing the rest via symmetry. By construction these two functions are equimeasurable, i.e.,

|{x ∈ (0, 1) : σ(x) ≥ c}| = |{x ∈ (0, 1) : σ_{∗}(x) ≥ c}|, ∀ c ∈ R,

and consequently, σ_{∗} ∈ ad. We are now in position to apply the following result of Hardy, Littlewood and Pólya, see Pólya and Szegö [31, p. 153].

If f and f_{1} are equimeasurable, g and g_{1} are equimeasurable, f ∈ L^{q}, g ∈ L^{q′}, and f_{1} and g_{1} are similarly ordered, then

∫_{0}^{1} fg dx ≤ ∫_{0}^{1} f_{1}g_{1} dx.  (2.14)

Given σ ∈ ad we now rearrange it as above into σ_{∗} and consider its corresponding u_{∗} ∈ E(σ_{∗}). Upon normalizing ∥u′_{∗}∥ = 1 we find

λ_{1}(σ) ≤ ∫_{0}^{1} σ^{p}|u′′_{∗}|^{2} dx ≤ ∫_{0}^{1} σ_{∗}^{p}|u′′_{∗}|^{2} dx = λ_{1}(σ_{∗}).

The first inequality is a consequence of Rayleigh’s principle, the second, of (2.13-14).

The stage now set, we address, in the next two sections, existence and necessary conditions for the generalized problem of Lagrange

sup_{σ∈ad} λ_{1}(σ).  (2.15)

3. Existence

We adopt the direct method of the calculus of variations and neglect to relabel subsequences. Denote by λ̂_{1} the value of (2.15) and by {σ_{n}} ⊂ ad an associated maximizing sequence, i.e., λ_{1}(σ_{n}) ↑ λ̂_{1}. Without loss we may assume that each σ_{n} is even about 1/2. We abbreviate λ_{1}(σ_{n}) to λ_{1,n} and denote by u_{1,n} a corresponding positive eigenfunction for which ∥u′_{1,n}∥ = 1 and ∫_{0}^{1} σ_{n}^{p}|u′′_{1,n}|^{2} dx = λ_{1,n} where ∥ · ∥ denotes the L^{2} norm. These normalizations, in light of (2.6), impose a uniform H^{2} bound on the sequence {u_{1,n}}. As a result, there exists a subsequence with weak H^{2} limit \overline{u} ∈ H_{0}^{2}. The imbedding of H^{2} in H^{1} being compact, we find ∥\overline{u}′∥ = 1, and so \overline{u} is not identically zero. The natural question is whether λ̂_{1} and \overline{u} are indeed an eigenvalue and eigenfunction for some column with corresponding \overline{σ} ∈ ad. If so, then \overline{σ} is necessarily the desired optimal design. This question was first addressed by Senatorov [35] in the context of a second order problem. He discovered that one must consider weak convergence of the reciprocals of the coefficients of the highest order term. This observation continues to hold for fourth order problems, the details of which we now sketch.

Consider the weak formulation

∫_{0}^{1} σ_{n}^{p}u′′_{1,n}v′′ dx = λ_{1,n} ∫_{0}^{1} u′_{1,n}v′ dx  ∀ v ∈ H_{0}^{2}.  (3.1)

Our previous remarks reveal that the right hand side converges to λ̂_{1} ∫_{0}^{1} \overline{u}′v′ dx for each such v. Regarding the left side we define ξ_{n} = σ_{n}^{p}u′′_{1,n}, and, as in (2.4), deduce from (3.1)

ξ_{n}(x) = (σ_{n}^{p}u′′_{1,n})(0) − λ_{1,n}u_{1,n}(x),

As the sequences {ξ_{n}} and {u_{1,n}} are uniformly bounded in L^{2} so to must be {(σ_{n}^{p}u′′_{1,n})(0)}. Consequently, ξ_{n} converges strongly in L^{2} to some function \overline{ξ}. The left side of (3.1) therefore converges to ∫_{0}^{1} \overline{ξ}v′′ dx. It remains to characterize this \overline{ξ}. Recalling the pointwise bounds on the σ_{n} we may assume that σ_{n}^{−p} converges in the weak* topology of L^{∞} to some function μ. Thus, ξ_{n}σ_{n}^{−p} converges weakly in L^{2} to μ\overline{ξ}. But ξ_{n}σ_{n}^{−p} = u′′_{1,n}, whose weak L^{2} limit is \overline{u}′′. Hence, \overline{ξ} = \overline{u}′′μ^{−1}. Defining \overline{σ} = μ^{−1/p} we may pass to the limit in (3.1) and obtain

∫_{0}^{1} \overline{σ}^{p}\overline{u}′′v′′ dx = λ̂_{1} ∫_{0}^{1} \overline{u}′v′ dx  ∀ v ∈ H_{0}^{2}.  (3.2)

As symmetry is preserved under weak ∗ convergence we find \overline{σ} to be even. In addition, the pointwise convergence of u_{1,n} to \overline{u} leaves \overline{u} positive and even. Now (3.2) implies that λ̂_{1} = λ_{j}(\overline{σ}) for some j. That j = 1 follows from the existence of a positive even first eigenfunction for \overline{σ} and the fact that \overline{u} is itself positive and even. We need only determine whether \overline{σ} ∈ ad. One may verify the pointwise bounds without trouble. With respect to the integral constraint we consider the convex function f : R → R, f(t) = t^{−1/p}. The integral functional φ ↦ ∫_{0}^{1} f(φ(x)) dx is then weak* lower semicontinuous on L^{∞}, see e.g., Dacorogna [12, Theorem 1.1]. As 1/σ_{n}^{p} converges weak* to 1/\overline{σ}^{p}, this allows us to conclude that

∫_{0}^{1} \overline{σ} dx = ∫_{0}^{1} f(1/\overline{σ}^{p}) dx ≤ lim inf ∫_{0}^{1} f(1/σ_{n}^{p}) dx = lim ∫_{0}^{1} σ_{n} dx = 1.  (3.3)

If indeed equality does not hold in (3.3) then there exists an even σ̂ ∈ ad such that σ̂(x) ≥ \overline{σ}(x) for almost every x ∈ (0, 1). From Rayleigh’s principle we then easily deduce λ_{1}(σ̂) ≥ λ_{1}(\overline{σ}). We have now proven

Theorem 3.1. There exists an even σ̂ ∈ ad for which λ_{1}(σ) ≤ λ_{1}(σ̂) for every σ ∈ ad.

Our choice of ad was motivated by our interest in the “shape” of the strongest column. Theorem 3.1 however, may also be applied in the search for the “composition” of the strongest column. For example, consider the design problem where one must combine two materials in fixed proportion so as to maximize the buckling load of the resulting column. The set of admissible designs is then

ad_{E} = {αχ(x) + β(1 − χ(x)) : χ is the characteristic

function of a subset of (0,1) with measure γ},

where α and β are the Youngs moduli of the respective materials with γ the volume fraction of the first. In this context, Theorem 3.1 states that λ_{1} attains its maximum on the weak* closure of ad_{E}, i.e., on

ad_{E}^{∗} = {αθ(x) + β(1 − θ(x)) : 0 ≤ θ(x) ≤ 1, ∫_{0}^{1} θ dx = γ}.

4. Necessary Conditions

We search now for a characterization of our optimal design, σ̂. Typical of many such problems, two distinct approaches are possible. Taking advantage of the variational structure, the so–called direct approach attempts to swap the order of the limits in

λ̂_{1} = λ_{1}(σ̂) = sup_{σ∈ad} inf_{u∈H_{0}^{2}} \frac{∫_{0}^{1} σ^{p}|u′′|^{2} dx}{∫_{0}^{1} |u′|^{2} dx},

inferring necessary conditions from the resulting saddle point. The indirect approach strives to determine σ̂ through knowledge of the tangents to the graph of σ ↦ λ_{1}(σ) and the normals to ad. Our implementation of these two approaches intersect in their reliance on (i) recent work of Auchmuty [2] on dual variational principles and (ii) a lopsided minimax principle.

Proposition 4.1. λ_{1}^{−1}(σ) = sup_{u∈H_{0}^{2}} A(σ, u), A(σ, u) = \sqrt{2}∥u′∥ − \frac{1}{2} ∫_{0}^{1} σ^{p}|u′′|^{2} dx.

u ↦ A(σ, u) attains its maximum only on those u ∈ E(σ) for which ∥u′∥ = \sqrt{2}λ_{1}^{−1}(σ).

Proof: In addition to being bounded above by 2λ_{1}^{−1}(σ), the map u ↦ A(σ, u) is coercive and weakly upper semicontinuous on H_{0}^{2} and therefore attains its maximum at some \overline{u} ∈ H_{0}^{2}. Necessarily, D_{2}A(σ, \overline{u}), the Gâteaux derivative of u ↦ A(σ, u) at \overline{u}, must vanish. That is,

∫_{0}^{1} σ^{p}\overline{u}′′v′′ dx = \sqrt{2}∥\overline{u}′∥^{−1} ∫_{0}^{1} \overline{u}′v′ dx ∀ v ∈ H_{0}^{2}.

As a result, \overline{u} is an eigenfunction corresponding to the eigenvalue \sqrt{2}∥\overline{u}′∥^{−1}. As \overline{u} maximizes u ↦ A(σ, u), this must be the least eigenvalue, λ_{1}(σ).

Proposition 4.2. Consider F : X × Y → R where X and Y are topological vector spaces and assume that x ↦ F (x, y) is concave and upper semicontinuous while y ↦ F (x, y) is convex and lower semicontinuous. If there exists a y_{o} ∈ Y and c_{o} ∈ R such that {x ∈ X; F (x, y_{o}) ≥ c_{o}} is compact and c_{o} < inf_{y∈Y} sup_{x∈X} F (x, y) then

sup_{x∈X} inf_{y∈Y} F (x, y) = inf_{y∈Y} sup_{x∈X} F (x, y).

Proof: This is a weakening of Theorem 3.7, Chapter 2, in Barbu and Precupanu [3].

It is with the indirect approach that we shall meet with the greatest success. For prior attempts in this context see Haug and Rousselet [18] and Choi and Haug [7]. Our principal tool is the generalized gradient of Clarke [10].

For a real valued Lipschitz function F on a Banach space X we consider the generalized directional derivative of F at x in the direction v,

F^{o}(x; v) ≡ lim sup ; {y → x \atop t ↓ 0} \frac{F (y + tv) − F (y)}{t}.

Denoting the dual of X by X^{∗} and x^{∗}(x) by ⟨x^{∗}, x⟩ when x^{∗} ∈ X^{∗} and x ∈ X, Clarke’s generalized gradient of F at x is the nonempty, convex, weak* compact set

∂F (x) ≡ {ξ ∈ X^{∗}; F^{o}(x; v) ≥ ⟨ξ, v⟩ ∀v ∈ X}.

We demonstrate that σ ↦ λ_{1}^{−1}(σ) is Lipschitz on Σ = {σ ∈ L^{∞}; ∥σ̂ − σ∥_{∞} < α/2}. Choose σ_{1}, σ_{2} ∈ Σ such that λ_{1}^{−1}(σ_{1}) > λ_{1}^{−1}(σ_{2}) and note that for u_{1} ∈ Argmax A(σ_{1}, ·), the set on which u ↦ A(σ_{1}, u) attains its maximum,

|λ_{1}^{−1}(σ_{1}) − λ_{1}^{−1}(σ_{2})| ≤ |A(σ_{1}, u_{1}) − A(σ_{2}, u_{1})|
≤ ∫_{0}^{1} |σ_{1}^{p} − σ_{2}^{p}||u′′_{1}|^{2} dx
≤ ∥u′′_{1}∥^{2}∥σ_{1}^{p} − σ_{2}^{p}∥_{∞} ≤ α^{−2p}p|2β|^{p−1}∥σ_{1} − σ_{2}∥_{∞}.

Without loss we assume that λ_{1}(σ̂) is a double eigenvalue. Then E(σ̂) is two dimensional and Argmax A(σ̂, ·) is the intersection of E(σ̂) with the sphere ∥u′∥ = \sqrt{2}λ̂_{1}^{−1}. It will be convenient to choose a basis {û_{1}, û_{2}} for E(σ̂) for which ∫_{0}^{1} û′_{i}û′_{j} dx = 2δ_{ij}λ̂_{1}^{−2}. For then,

Argmax A(σ̂, ·) = {aû_{1} + bû_{2}; a^{2} + b^{2} = 1}.  (4.1)

Regarding the Gâteaux derivative of σ ↦ A(σ, u) at σ̂ in the direction η we have

⟨D_{1}A(σ̂, u), η⟩ = −\frac{p}{2} ∫_{0}^{1} ησ̂^{p−1}|u′′|^{2} dx.  (4.2)

Denoting convex hull by ‘co’, the sense in which the gradient of a maximum is the maximum of the gradients is Theorem 4.3. ∂λ_{1}^{−1}(σ̂) = co \{−\frac{p}{2}σ̂^{p−1}(aû′′_{1} + bû′′_{2})^{2}; a^{2} + b^{2} = 1\}.

Proof: From (4.1) and (4.2) this set is precisely

co {D_{1}A(σ̂, u); u ∈ Argmax A(σ̂, ·)}.  (4.3)

Our claim does not fit neatly into Clarke’s result [10, Theorem 2.8.2] on the generalized gradient of a pointwise maximum. The contortions involved in fitting our problem to Clarke’s hypotheses are no less difficult, and far less instructive, than an independent proof.

Let us denote the set in (4.3) by Ξ. We show that Ξ ⊂ ∂λ_{1}^{−1}(σ̂). For ξ ∈ Ξ and η ∈ L^{∞}

⟨ξ, η⟩ = ∑_{i=1}^{n} μ_{i}⟨D_{1}A(σ̂, u_{i}), η⟩

= ∑_{i=1}^{n} μ_{i} lim_{t↓0} \frac{A(σ̂ + tη, u_{i}) − A(σ̂, u_{i})}{t}

≤ ∑_{i=1}^{n} μ_{i} lim sup_{t↓0} \frac{λ_{1}^{−1}(σ̂ + tη) − λ_{1}^{−1}(σ̂)}{t}

≤ (λ_{1}^{−1})^{o}(σ̂; η),

hence ξ ∈ ∂λ_{1}^{−1}(σ̂).

Regarding the reverse inclusion we define

g(σ; η) = max_{ξ∈Ξ} ⟨ξ, η⟩, σ ∈ Σ, η ∈ L^{∞}

and prove

(λ_{1}^{−1})^{o}(σ; η) ≤ g(σ; η).

Select σ_{n} → σ in L^{∞} and t_{n} ↓ 0 in R such that

q_{n} ≡ \frac{λ_{1}^{−1}(σ_{n} + t_{n}η) − λ_{1}^{−1}(σ_{n})}{t_{n}}

converges to (λ_{1}^{−1})^{o}(σ; η). Select u_{n} ∈ Argmax A(σ_{n} + t_{n}η, ·) and note that

q_{n} ≤ \frac{A(σ_{n} + t_{n}η, u_{n}) − A(σ_{n}, u_{n})}{t_{n}}

with the right side equal to ⟨D_{1}A(σ_{n} + \overline{t}_{n}η, u_{n}), η⟩, for some \overline{t}_{n} ∈ (0, t_{n}), by the Mean Value Theorem. As σ_{n} + t_{n}η → σ in L^{∞} and u_{n} ∈ Argmax A(σ_{n} + t_{n}η, ·) we recall from our work in Theorem 3.1 that u_{n} ⇀ u ∈ Argmax A(σ, ·) in H^{2} and (σ_{n} + t_{n}η)^{p}u′′_{n} → σ^{p}u′′ in L^{2}, and hence u′′_{n} → u′′ in L^{2}, i.e., u_{n} → u in H^{2}. Recalling (4.2) this establishes

⟨D_{1}A(σ_{n} + \overline{t}_{n}η, u_{n}), η⟩ → ⟨D_{1}A(σ, u), η⟩

with u ∈ Argmax A(σ, ·). As a result,

(λ_{1}^{−1})^{o}(σ; η) = lim_{n→∞} q_{n} ≤ ⟨D_{1}A(σ, u), η⟩ ≤ g(σ; η).

If ζ is now an element of ∂λ_{1}^{−1}(σ̂) then g(σ̂; η) ≥ ⟨ζ, η⟩ for each η ∈ L^{∞}. Consequently,

0 = min_{η∈L^{∞}} max_{ξ∈Ξ} ⟨ξ − ζ, η⟩.

Noting that Ξ is closed and bounded in L^{1} and finite dimensional (it lies in the span of {|û′′_{1}|^{2}, û′′_{1}û′′_{2}, |û′′_{2}|^{2}}) we find it compact in (L^{∞})^{∗}. Invoking Proposition 4.2 yields a ξ̂ ∈ Ξ for which

⟨ξ̂− ζ, η⟩ = 0 ∀ η ∈ L^{∞}.

It follows that ζ = ξ̂ and so ∂λ_{1}^{−1}(σ̂) ⊂ Ξ.

This proof, though identical in outline to Clarke [10, Theorem 2.8.2], has exploited additional properties of λ_{1} and A to make up for the missing hypotheses. Observe that when λ_{1}(σ̂) is simple the generalized gradient reduces to the singleton

∂λ_{1}^{−1}(σ̂) = {D_{1}A(σ̂, û_{1})}.

As zero is not a tangent direction to λ_{1}^{−1} at σ̂, i.e., 0 ≠ ∂λ_{1}^{−1}(σ̂), we are compelled to investigate the constraint set ad. Separating the equality from the inequality constraints brings

C ≡ {σ ∈ L^{∞}; α ≤ σ(x) ≤ β} and V (σ) ≡ ∫_{0}^{1} σ dx.

As σ̂ minimizes σ ↦ λ_{1}^{−1}(σ) subject to σ ∈ C and V (σ) = 1 we deduce from the Lagrange Multiplier Rule, [10, Theorem 6.1.1], that a nontrivial linear combination of elements in ∂λ_{1}^{−1}(σ̂) and ∂V (σ̂) is normal to C at σ̂. In particular,

(ν_{1}∂λ_{1}^{−1}(σ̂) + ν_{2}∂V (σ̂)) ∩ N_{C}(σ̂) ≠ ∅,

where ν_{1} ≤ 0, ν_{1}^{2} + ν_{2}^{2} > 0, and

N_{C}(σ̂) = {ζ ∈ (L^{∞})^{∗}; ∫_{0}^{1}(σ̂ − σ) dζ ≥ 0, ∀σ ∈ C}

is the cone of normals to C at σ̂. In light of our previous calculations, and the fact that ∂V (σ̂) = 1 ∈ L^{1}, there exists a ξ̂ ∈ ∂λ_{1}^{−1}(σ̂) for which

∫_{0}^{1}(σ̂ − σ)(ν_{1}ξ̂+ ν_{2}) dx ≥ 0,  ∀ σ ∈ C.  (4.4)

Observing that ν_{1}ξ̂ ≥ 0 we find that ν_{2} ≥ 0 requires, through (4.4), that σ̂ ≡ β, an impossibility. Likewise, should ν_{1} = 0, (4.4) would require σ̂ ≡ α (since ν_{2} < 0). Taking ℓ^{2} = ν_{2}/ν_{1} we arrive at

∫_{0}^{1}(σ̂ − σ)(ξ̂+ ℓ^{2}) dx ≤ 0, ∀ σ ∈ C.

The subsequent reduction to pointwise optimality conditions follows a well known course, see e.g., Cea and Malanowski [8]. In particular,

σ̂(x) = α ⇒ −ξ̂(x) ≤ ℓ^{2}  (4.5)
α < σ̂(x) < β ⇒ −ξ̂(x) = ℓ^{2}  (4.6)
σ̂(x) = β ⇒ −ξ̂(x) ≥ ℓ^{2}  (4.7)

for almost every x ∈ (0, 1). To appreciate this result we must recall that ξ̂ ∈ ∂λ_{1}^{−1}(σ̂) means

−ξ̂(x) = \frac{p}{2} ∑_{i=1}^{n} t_{i}σ̂^{p−1}(x)(a_{i}û′′_{1}(x) + b_{i}û′′_{2}(x))^{2}, where

t_{i} ≥ 0, ∑_{i=1}^{n} t_{i} = 1, and a_{i}^{2} + b_{i}^{2} = 1.

On expanding this sum of squares, (4.6) becomes

σ̂^{p−1} (δ_{1}|û′′_{1}|^{2} + δ_{2}|û′′_{2}|^{2} + δ_{3}û′′_{1}û′′_{2}) = 1, where (4.8)

δ_{1} = \frac{p}{2} ∑_{i=1}^{n} t_{i}a_{i}^{2}/ℓ^{2}, δ_{2} = \frac{p}{2} ∑_{i=1}^{n} t_{i}b_{i}^{2}/ℓ^{2}, δ_{3} = p ∑_{i=1}^{n} t_{i}a_{i}b_{i}/ℓ^{2}.

Observing that δ_{1}δ_{2} indeed dominates δ_{3}^{2}/4, we have recovered (1.8), the necessary condition of Bratus and Seiranian [6] and Masur [22]. If in fact δ_{1}δ_{2} = δ_{3}^{2}/4 then for û ≡ \sqrt{δ_{1}}û_{1} + \sqrt{δ_{2}}û_{2}, equation (4.8) yields σ̂^{p−1}|û′′|^{2} = 1, the optimality condition of Tadjbakhsh and Keller. As û is an eigenfunction and therefore admits at least two inflection points the pointwise bounds must become active, i.e., σ̂^{p−1}|û′′|^{2} = 1 can not hold on the entire interval. Ignoring any bound constraints, Masur [22] and Seiranian [34] found a σ and two orthogonal elements of E(σ) for which (4.8) holds with p = 2. This appears to be the design obtained by Olhoff and Rasmussen [25] and, by all indications, the one preferred by our algorithm as well (see §7 Fig. 1). It appears likely that in this case the bound constraints are inactive due to the fact that where σ̂ is less than one, σ̂^{2} is much less than one. As σ̂^{2} is the quantity that appears in the Rayleigh quotient we expect it to be as large as possible. This suggests that σ̂ is bounded away from zero, independent of α. This lower bound with the integral constraint together support the conjecture that σ̂ is in fact bounded above as well. Hence, when α and β are respectively chosen below and above these ‘natural’ bounds, condition (4.8) is free to stand on its own. Clearly these natural bounds must depend on p. In fact, we shall provide numerical evidence in §7 in favor of the argument that the natural lower (upper) bound is an increasing (decreasing) function of p for p > 1.

Unfortunately, it is not known whether (4.8) is a sufficient condition for optimality. The proof of sufficiency offered by Tadjbakhsh and Keller [37] is incorrect. They proceed as if λ_{1}(σ) corresponds to the least eigenvalue of (2.5) and accordingly admit all functions that satisfy the boundary conditions as test functions in a Rayleigh principle argument. In fact, (2.5) possesses a double zero eigenvalue, hence only those functions that are orthogonal to the first two eigenfunctions can be admitted. We remark that Ramm’s claim [32], that Tadjbakhsh and Keller mistakenly applied Hölders inequality in their sufficiency proof, is incorrect, though [37, §6 (25)] is only valid for n < 0.

Though (4.8) need not hold over the entire length of the column, we now show that where it does hold it requires that σ̂ be smooth.

Theorem 4.4. If α < σ̂(x) < β for each x ∈ (a, b) ⊂ (0, 1) then σ̂ ∈ C^{∞}(a, b).

Proof: We observed in (2.4) that

σ̂^{p}û′′_{i} = l_{i} − λ̂_{1}û_{i},  (4.9)

where l_{i} is an affine function of x. Now multiply (4.8) by σ̂^{p+1},

δ_{1}(σ̂^{p}û′′_{1})^{2} + δ_{2}(σ̂^{p}û′′_{2})^{2} + δ_{3}(σ̂^{p}û′′_{1})(σ̂^{p}û′′_{2}) = σ̂^{p+1}.  (4.10)

¿From (4.9) we find, on recalling H_{0}^{2} ⊂ C^{1}, that each term on the left of (4.10) is C^{1}, and hence that σ̂ ∈ C^{1}. Writing (4.9) in the form

û′′_{i} = \frac{l_{i} − λ̂_{1}û_{i}}{σ̂^{p}} ,

we conclude û′′_{i} ∈ C^{1}, that is û_{i} ∈ C^{3}. Repeating this exact argument leads to σ̂ ∈ C^{3} and û_{i} ∈ C^{5}. The result then follows from continued repetitions.

Having succeeded in pursuing the indirect approach we now look to the possibility (and the implications) of exchanging the limits in the characterization

λ̂_{1}^{−1} = A(σ̂, û) = inf_{σ∈ad} sup_{u∈H_{0}^{2}} A(σ, u).

Recalling Proposition 4.2, this will require convexity and lower semicontinuity of σ ↦ A(σ, u), and concavity and upper semicontinuity of u ↦ A(σ, u) as well as compactness of one of its upper level sets.

Remark 4.5. We noted the weak H^{2} upper semicontinuity of u ↦ A(σ, u) in Proposition 4.1. As {u ∈ H_{0}^{2}; A(σ, u) ≥ c} is bounded it is also weakly compact (independent of c ∈ R and σ ∈ ad). Convexity of σ ↦ A(σ, u) follows on restricting p ≤ 1.

The two remaining properties will require more work. Note that u ∈ Argmax A(σ, ·) implies A(σ, u) = A(σ, −u) = λ_{1}^{−1}(σ) while A(σ, 0) = 0. Hence, u ↦ A(σ, u) is not concave on any set that contains Argmax A(σ, ·). This suggests that we examine the half-spaces exterior to \{u ∈ H_{0}^{2}; ∥u′∥ ≤ \sqrt{2}λ_{1}^{−1}(σ)\}. Unfortunately this ball, and hence its support planes, depend on σ. Consequently, if we expect these half-spaces to vary continuously with σ we must be careful in our choosing. This choice is greatly facilitated by the assumption that σ lies in ad_{s}, those functions in ad that are even about 1/2. For in this case one may speak unambiguously of u_{1}, the positive even eigenfunction corresponding to λ_{1}(σ). We normalize ∥u′_{1}∥ = \sqrt{2}λ_{1}^{−1}(σ), and consider the associated half-space

Π_{σ} = {v ∈ H_{0}^{2}; λ_{1}^{2}(σ) ∫_{0}^{1} u′_{1}v′ dx > 2}.

Proposition 4.6. For σ ∈ ad_{s}, u ↦ A(σ, u) is concave on Π_{σ}.

Proof: The quadratic form associated with the second Gâteaux derivative of u ↦ A(σ, u) at \overline{u} ∈ Π_{σ} satisfies

⟨D_{2}^{2}A(σ, \overline{u})v, v⟩ = \sqrt{2}∥\overline{u}′∥^{−1} ∫_{0}^{1} |v′|^{2} dx − ∫_{0}^{1} σ^{p}|v′′|^{2} dx − \sqrt{2}^{−1}∥\overline{u}′∥^{−3} (∫_{0}^{1} \overline{u}′v′ dx)^{2}
≤ (\sqrt{2}∥\overline{u}′∥^{−1} − λ_{1}(σ)) ∫_{0}^{1} |v′|^{2} dx
≤ 0, ∀ v ∈ H_{0}^{2}.

This suggests that we penalize A with the indicator function of Π_{σ},

π(σ, u) = \{ \begin{matrix}0, & if u ∈ Π_{σ}; \\ ∞, & otherwise. \end{matrix}

This not only guarantees concavity but also respects lower semicontinuity.

Proposition 4.7. σ ↦ A(σ, u) − π(σ, u) is lower semicontinuous for the strong L^{∞} topology on ad_{s}.

Proof: Now σ_{n} → σ in L^{∞} clearly implies A(σ_{n}, u) → A(σ, u) for each u ∈ H_{0}^{2}. Regarding lim sup π(σ_{n}, u) ≤ π(σ, u), it suffices to show that

π(σ_{n}, u) → 0, ∀ u ∈ Π_{σ}.

¿From the proof of Theorem 3.1 it is clear that λ_{1}(σ_{n}) → λ_{1}(σ) and u_{1}(σ_{n}) ⇀ u_{1}(σ) in H_{0}^{2}. Hence

λ_{1}^{2}(σ_{n}) ∫_{0}^{1} u′_{1}(σ_{n})u′ dx → λ_{1}^{2}(σ) ∫_{0}^{1} u_{1}(σ)u′ dx > 2.

¿From this we conclude that u is eventually in each Π_{σ_{n}}, i.e., π(σ_{n}, u) = 0.

One may now modify Proposition 4.2 (see Cox and McLaughlin [11, §7]), and conclude Theorem 4.8. If p ≤ 1 then (σ̂, u_{1}(σ̂)) is a saddle point for A over ad_{s} × H_{0}^{2}. That is, denoting u_{1}(σ̂) by û,

A(σ̂, u) ≤ A(σ̂, û) ≤ A(σ, û) ∀ (σ, u) ∈ ad_{s} × H_{0}^{2}.

The latter inequality yields the following maximum principle

∫_{0}^{1} σ^{p}|û′′|^{2} dx ≤ ∫_{0}^{1} σ̂^{p}|û′′|^{2} dx ∀ σ ∈ ad_{s}.

The subsequent pointwise conditions call for an ℓ^{2} > 0 such that

σ̂(x) = α ⇒ σ̂^{p−1}(x)|û′′(x)|^{2} ≤ ℓ^{2}  (4.11)
α < σ̂(x) < β ⇒ σ̂^{p−1}(x)|û′′(x)|^{2} = ℓ^{2}  (4.12)
σ̂(x) = β ⇒ σ̂^{p−1}(x)|û′′(x)|^{2} ≥ ℓ^{2}  (4.13)

for almost every x ∈ (0, 1). As σ ↦ λ_{1}^{−1}(σ) is convex when p ≤ 1 these conditions are also sufficient. Hence, we see that where the direct method applies it gives more information. In particular, the necessary conditions (4.11-13) involve only a single buckling mode. Comparing these to the more general conditions in (4.5-7) suggests that λ_{1}(σ̂) is indeed a simple eigenvalue when p ≤ 1. We shall see numerical evidence of this in §7. The critical case, p = 1, where the optimal buckling load changes multiplicity, has received considerable attention. In this case, the right side of (4.11-13) is independent of σ̂. In particular, a number of workers have claimed that

|û′′(x)| = ℓ.  (4.14)

We remark however that in the absence of a second buckling mode the bound constraints must become active near the inflection points of û, making (4.11) and (4.13) indeed necessary. Nonetheless, Seiranian [34], who deduced (4.14) from (1.8), proceeded to solve (4.14) in conjunction with (2.3), yielding

σ̂(x) = \{ \begin{matrix}3/2(1 − 16x^{2}), & if 0 ≤ x ≤ 1/4 \\ 3/2(16x − 16x^{2} − 3), & if 1/4 ≤ x ≤ 3/4 \\ 3/2(32x − 16x^{2} − 15), & if 3/4 ≤ x ≤ 1. \end{matrix}  (4.15)

On evaluating the Rayleigh quotient with this σ̂ and a specific C^{1} test function Seiranian arrived at a buckling load of 48. This design, like that of Tadjbakhsh and Keller for p = 2, vanishes at 1/4 and 3/4. Unlike the design of Tadjbakhsh and Keller however, we are not able to show it to be suboptimal. We can only stress that lacking an existence proof for α = 0, p = 1, there is no reason to believe that (4.14) is a necessary condition for optimality.

5. Other Boundary Conditions

Intent on a clean exposition we have to this point concentrated solely on the clampedclamped boundary conditions u(0) = u′(0) = u(1) = u′(1) = 0. We now apply the work of the previous sections to the other standard sets of boundary conditions, in particular, hinged and free. A column is said to be free at a point when no conditions are prescribed, while it is hinged, or simply supported, when one requires its displacement to vanish there. As a matter of notation, the weak formulation of the buckled column equation will read

∫_{0}^{1} σ^{p}u′′v′′ dx = μ ∫_{0}^{1} u′v′ dx,  ∀ v ∈ V_{i,j}  (5.1)

where V_{i,j} is a subspace of H^{2}, with i and j chosen from {0, 1, 2} according to whether the respective end is either free, hinged, or clamped. For example,

V_{1,2} = {u ∈ H^{2}; u(0) = 0, u(1) = u′(1) = 0}

specifies the hinged-clamped column. We denote the least eigenvalue of (5.1) by μ_{i,j}(σ), and the corresponding space of eigenfunctions by E_{i,j}(σ). As before, u ∈ E_{i,j}(σ) implies that both u and σ^{p}u′′ are elements of C^{1}([0, 1]). In addition, such functions satisfy so called natural boundary conditions. In particular, if i = 1 then, in addition to u(0) = 0 we find

σ^{p}u′′(0) = 0,  (5.2)

while if i = 0 we have, in addition to (5.2),

(σ^{p}u′′)′(0) + μ_{0,j}(σ)u′(0) = 0.  (5.3)

We shall consider only those μ_{i,j}(σ) for which i + j ≥ 2, as otherwise μ_{i,j}(σ) = 0. For comparison purposes we record these eigenvalues in the case of the uniform column.

μ_{0,2}(1) = π^{2}/4, μ_{1,1}(1) = π^{2}, μ_{1,2}(1) ≈ 2.046π^{2}, μ_{2,2}(1) = 4π^{2}.  (5.4)

Clearly, μ_{i,j}(1) = μ_{j,i}(1). Analogous to (2.5), for i + j ≥ 2, (5.4) gives the uniform bounds

π^{2}α^{p}/4 ≤ μ_{i,j}(σ) ≤ 4π^{2}β^{p},  ∀ σ ∈ ad.  (5.5)

As in §2 we address the multiplicity of μ_{i,j}(σ) and the presence of positive eigenfunctions.

Lemma 5.1. For σ ∈ ad,

(a) If 2 ≤ i + j < 4, then μ_{i,j}(σ) is simple and there exists a corresponding positive eigenfunction.

(b) μ_{0,2}(σ) < μ_{1,2}(σ) and μ_{1,1}(σ) < μ_{1,2}(σ) < μ_{2,2}(σ).

Proof: (a) Seiranian noted for these boundary conditions that (5.1) is equivalent, except for the presence of a simple zero eigenvalue when the product ij equals 2, to a second order problem with separated boundary conditions. It now follows from the oscillation theory of Stürm, see, e.g., Atkinson [1], that each μ_{i,j}(σ) is simple and that for ij ≠ 2, σ^{p}u′′ is of one sign for each u ∈ E_{i,j}(σ). In case ij equals 0 or 1 this yields respectively a positive convex or concave element of E_{i,j}(σ). When ij = 2 we find that σ^{p}u′′ vanishes exactly once on (0, 1) for each u ∈ E_{i,j}(σ). Here we find an eigenfunction that is convex on (0, x_{0}) and concave on (x_{0}, 1) for some x_{0}. As this function must vanish at 0 and 1 we conclude that it must be positive on (0, 1).

(b) As V_{i+1,j} ⊂ V_{i,j} we find μ_{i,j}(σ) ≤ μ_{i+1,j}. Should equality hold, we conclude E_{i+1,j}(σ) ⊂ E_{i,j}(σ). As in (2.4), for u ∈ E_{i,j}(σ) we deduce from (5.1) that

(σ^{p}u′′)(x) = ((σ^{p}u′′)′(0) + μ_{i,j}(σ)u′(0)) x + σ^{p}u′′(0) + μ_{i,j}(σ)u(0) − μ_{i,j}(σ)u(x). (5.6)

If μ_{0,2}(σ) = μ_{1,2}(σ) then for each u ∈ E_{1,2}(σ) ⊂ E_{0,2}(σ) equation (5.6), in view of (5.3), reads

(σ^{p}u′′)(x) = −μ_{i,j}u(x).  (5.7)

On recalling that u(1) = u′(1) = 0 we see that u satisfies a linear homogeneous equation with zero terminal data, and hence, u ≡ 0.

If μ_{1,1}(σ) = μ_{1,2}(σ) then for each u ∈ E_{1,2}(σ) ⊂ E_{1,1}(σ) (5.6), in view of (5.2), reads

(σ^{p}u′′)(x) = ((σ^{p}u′′)′(0) + μ_{i,j}(σ)u′(0)) x − μ_{i,j}(σ)u(x).  (5.8)

So (σ^{p}u′′)(1) = (σ^{p}u′′)′(0) + μ_{i,j}(σ)u′(0). But (σ^{p}u′′)(1) = 0 so (5.8) reduces to (5.7) and again the clamped conditions at 1 imply that u ≡ 0.

If μ_{1,2}(σ) = μ_{2,2}(σ) then for each u ∈ E_{2,2}(σ) ⊂ E_{1,2}(σ) equation (5.6), in view of (5.2), reads

(σ^{p}u′′)(x) = (σ^{p}u′′)′(0)x − μ_{i,j}(σ)u(x).

Hence (σ^{p}u′′)(1) = (σ^{p}u′′)′(0), from which we conclude that u is either identically zero or not of one sign. This excludes the positive element of E_{1,2}(σ) established in part (a).

Thanks to the presence of positive first eigenfunctions, the existence theory of §3 applies directly to the problem of Lagrange

sup_{σ∈ad} μ_{i,j}(σ),  2 ≤ i + j ≤ 4.  (5.9)

We note that only for symmetric boundary conditions, i.e., i = j, should one expect an even optimal design. As μ_{i,j}(σ) is simple when i + j < 4 we deduce from Theorem 4.3 and conditions (4.5-7) that

σ̂_{i,j}(x) = α ⇒ σ̂_{i,j}^{p−1}(x)|û′′(x)|^{2} ≤ ℓ^{2}  (5.10)
α < σ̂_{i,j}(x) < β ⇒ σ̂_{i,j}^{p−1}(x)|û′′(x)|^{2} = ℓ^{2}  (5.11)
σ̂_{i,j}(x) = β ⇒ σ̂_{i,j}^{p−1}(x)|û′′(x)|^{2} ≥ ℓ^{2},  (5.12)

for almost every x ∈ (0, 1), where û ∈ E_{i,j}(σ̂_{i,j}). As before, σ̂_{i,j} is smooth where (5.11) holds.

The right side of (5.11) is the sole necessary condition offered by Keller [20] and Tadjbakhsh and Keller [37]. We now investigate the extent to which their claim is valid.

Recall that their analysis of the clamped–clamped column erred in neglecting (a) double eigenvalues, and (b) bounds on σ. As the previous lemma precludes the former phenomenon we need only consider the latter. The observation to be made is that (5.10) and (5.12) are only needed near the zeros of û′′. As noted above, members of E_{1,1}(σ) and E_{2,0}(σ) have second derivatives of one sign. As such, in these cases, (5.11) stands on its own (with the minor adjustment that σ be allowed to vanish at 0 and/or 1). In addition, as the related second order problems are fully equivalent, i.e., there are no spurious eigenvalues, Tadjbakhsh and Keller’s sufficiency proof is correct. In summary, Keller [20] has the correct necessary condition for the hinged–hinged column, Tadjbakhsh and Keller [37] have the correct necessary condition for the clamped–free column, and the proof of sufficiency in [37] holds for both. We now recall their analytical solutions to these problems.

Keller, in [20], with p = 2 and i = j = 1 reconciled (5.11) and (5.1) and found

σ̂_{1,1}(x) = \frac{4}{3} sin^{2} θ(x), 0 ≤ θ ≤ π, (5.13)

θ(x) − \frac{1}{2} sin 2θ(x) = πx, 0 ≤ x ≤ 1.

We have observed that this is a shortened cycloid with parametrization

\begin{aligned}x(t) & = \frac{3}{4π} (\frac{2}{3}(t − sin t)) \\ y(t) & = \frac{2}{3}(1 − cos t) \end{aligned} 0 ≤ t ≤ 2π.

This column buckles under an axial load of 4π^{2}/3. In [37], Tadjbakhsh and Keller with p = 2 and i = 2, j = 0 reconciled (5.11) and (5.1) and found

σ̂_{2,0}(x) = \frac{4}{3} sin^{2} θ(x), −π/2 ≤ θ ≤ 0, (5.14)

θ(x) − \frac{1}{2} sin 2θ(x) + π/2 = πx/2, 0 ≤ x ≤ 1,

our parametrization being,

\begin{aligned}x(t) & = \frac{3}{2π} (\frac{2}{3}(t − sin t)) + 1 \\ y(t) & = \frac{2}{3}(1 − cos t) \end{aligned} − π ≤ t ≤ 0.

This column buckles under an axial load of π^{2}/3. Having argued in favor of the existing solutions to the clamped–free and hinged–hinged problems we now turn to the clamped– hinged problem.

We saw in Lemma 5.1 that the second derivative of each function in E_{2,1}(σ) must change sign. The effect of this is that (5.11) forces σ̂_{2,1} to vanish at an interior point. In particular, when Tadjbakhsh and Keller reconciled (5.1) and (5.11) they found

σ̂_{2,1}(x) = \frac{4 sin^{2} θ(x)}{3 sin^{2} θ(0)}, θ(0) ≤ θ ≤ π, (5.15)

θ(x) − \frac{1}{2} sin 2θ(x) + \frac{1}{2} sin 2θ(0) − θ(0) = x(π + \frac{1}{2} sin 2θ(0) − θ(0)), 0 ≤ x ≤ 1,
\frac{1}{2} sin 2θ(0) − θ(0) = −\frac{2}{3} sin^{3} θ(0) cos^{−1} θ(0) − π.

Taking a = \frac{1}{2} sin 2θ(0) − θ(0) note that this σ̂_{2,1}(x) vanishes at a/(π + a). Tadjbakhsh and Keller assert that the column built according to (5.15) will not buckle under loads less than approximately 27.22 in magnitude. We show in the appendix that this column can not withstand loads exceeding π^{2}/3 – and so, in fact, is much weaker than the uniform column. In addition, as μ_{2,1}(σ) corresponds to the second eigenvalue of its associated second order problem, the sufficiency proof of [37] is invalid. Hence, (5.15) is not an optimal design. In summary, (5.11) can not stand alone in the clamped–hinged case, σ(x) ≥ α is is indeed an active constraint, so rendering (5.10) absolutely necessary. We suspect that there exists no solution to (5.9) when ij = 2 and α = 0.

6. The Finite Dimensional Problem

We discretize the interval [0, h, 2h,…, (N −1)h = 1] and approximate V_{i,j} by the finite dimensional space V_{i,j}^{h}, the subspace of V_{i,j} whose elements, when restricted to [kh, (k+1)h], are cubic polynomials (see Strang and Fix [36]). As each member of V_{i,j}^{h} is completely determined by the value of it and its derivative at each of the N mesh points, we identify V_{i,j}^{h} with R^{2N−i−j}. We next approximate ad with the class of piecewise constant functions

ad^{h} ≡ {σ ∈ R^{N−1} : α ≤ σ_{k} ≤ β, ∑_{k=1}^{N−1} σ_{k} = N − 1}.

We have refrained from labeling elements of ad^{h} by σ^{h} to avoid confusion with powers of σ. In this context, the infinite dimensional eigenvalue problem of (5.1) is now approximated by

B_{h}(σ)q_{h} = μK_{h}q_{h}, σ ∈ ad^{h}, q_{h} ∈ R^{2N−i−j}.  (6.1)

B_{h}(σ) and K_{h}, the so–called bending and stiffness matrices, are each real, (2N − i − j) × (2N − i − j), symmetric, positive definite, and banded with half bandwidth of four. Our interest is, of course, in μ_{i,j}^{h}(σ), the least eigenvalue of (6.1). For, as h → 0, one finds e.g., in [36], that μ_{i,j}^{h}(σ) → μ_{i,j}(σ). The connection between the finite and infinite dimensional problems now understood, we concentrate solely on (6.1). It should cause no confusion if, in our presentation of the finite dimensional optimization problem, we suppress most dependence on h, i, and j. With this, (6.1) becomes

B(σ)q = λKq, σ ∈ ad^{h}, q ∈ R^{n},  (6.2)

and we denote its least eigenvalue by λ_{1}(σ). Our finite dimensional problem of Lagrange is now

max_{σ∈ad^{h}} λ_{1}(σ).  (6.3)

The care that was taken in differentiating σ ↦ λ_{1}(σ) in §4 must also be exercised here. The occurrence of multiple eigenvalues is still possible. Clarke [10, Proposition 2.8.8] specifies the generalized gradient of the largest eigenvalue of a symmetric matrix in terms of a convex hull; see also [28] and [29]. Though such a characterization may suffice for an analytical description, as in §4, for computational purposes we have found it more useful to specify first order conditions in terms of, less well known, “dual matrices”. We state the result in general terms. We shall need S^{n}, the class of n × n real symmetric matrices, and the Frobenius matrix inner product, ⟨A, B⟩ = tr A^{T} B.

Theorem 6.1. Let B : R^{N−1} → S^{n} be (Fréchet) differentiable with B_{k}(σ) = ∂B(σ)/∂σ_{k} and let K be a fixed symmetric positive semidefinite matrix of the same order n. Assume σ ∈ ad^{h} is such that λ_{1}(σ) has multiplicity t, with corresponding eigenvectors given by the columns of a matrix Q_{1} ∈ R^{n×t}, normalized so that Q_{1}^{T} KQ_{1} = I. Then a necessary condition for σ to solve (6.3) is that there exist a symmetric positive semidefinite matrix U of order t, with trace equal to one, and Lagrange multipliers ν and γ_{k}, k = 1,…, N −1, such that

⟨U, Q_{1}^{T} B_{k}(σ)Q_{1}⟩ = ν + γ_{k}, and  (6.4)

σ_{k} = α ⇒ γ_{k} ≤ 0  (6.5)
α < σ_{k} < β ⇒ γ_{k} = 0  (6.6)
σ_{k} = β ⇒ γ_{k} ≥ 0  (6.7)

for each k. Furthermore, this condition is also sufficient for optimality in the case that σ ↦ B(σ) is affine.

Proof: In the following we use the notation U ≥ 0 to mean that a symmetric matrix U is positive semidefinite. Regarding λ_{1} : S^{n} → R, we invoke Rayleigh’s principle in

λ_{1} = min {⟨q, Bq⟩; q ∈ R^{n}, ⟨q, Kq⟩ = 1}
= min {⟨qq^{T} , B⟩; q ∈ R^{n}, ⟨q, Kq⟩ = 1}

Let Q ∈ R^{n×n} be any matrix satisfying

Q^{T} KQ = I.  (6.8)

It is easily shown that

co {qq^{T}; q ∈ R^{n}, ⟨q, Kq⟩ = 1} = {QÛQ^{T}; Û ∈ S^{n}, tr Û = 1, Û ≥ 0},

by using the spectral decomposition of Û, which by assumption has nonnegative eigenvalues adding to one, to obtain the requisite convex combination showing that the second set is contained in the first. It follows that

λ_{1} = min {⟨QÛQ^{T} , B⟩; Û ∈ S^{n}, tr Û = 1, Û ≥ 0}.  (6.9)

Now take Q to be a matrix whose columns are eigenvectors of (6.2), normalized so that (6.8) holds. The first t columns of Q are the columns of Q_{1} and

Q^{T} BQ = Diag(λ_{i}),

where λ_{1} ≤ λ_{2} ≤ … are the eigenvalues of (6.2), repeated according to multiplicity. Therefore, the matrices achieving the minimization in (6.9) are those defined by

Û = Q ( \begin{matrix}U & 0 \\ 0 & 0 \end{matrix} ) Q^{T} = Q_{1}UQ_{1}^{T}

where U ∈ S^{t}, with tr U = 1 and U ≥ 0. Consequently, the generalized gradient of B ↦ −λ_{1}(B) is the set of such matrices Û (see Rockafellar [33, pages 29 and 35] or Clarke [10, §2.8]), no convex hull operation being required since the set of such Û is convex.

With λ_{1}(σ) = (λ_{1} ∘ B)(σ), the desired necessary conditions now follow from (i) the chain rule for generalized gradients [10, Theorem 2.3.10], (ii) the standard Lagrange multiplier rule [10, Theorem 6.1.1], and (iii) properties of the inner product. In particular,

∂(−λ_{1}(σ)) = {v ∈ R^{N−1}; v_{k} = ⟨U, Q_{1}^{T} B_{k}(σ)Q_{1}⟩, U ∈ S^{t}, U ≥ 0, tr U = 1}. (6.10)

These necessary conditions are also sufficient in the case that σ ↦ B(σ) is affine because the composition of a concave function with an affine function is concave.

Our attention to ∂(−λ_{1}) in (6.10) and ∂λ_{1}^{−1} in Theorem 4.3 rather than simply ∂λ_{1} is merely an artifact of Clarke’s concern with functions defined as pointwise maxima rather minima. Here, it was convenient to characterize −λ_{1} as the maximum of a Rayleigh quotient where, in §4, we found it more advantageous to maximize a functional of Auchmuty, and hence to consider λ_{1}^{−1}.

A different proof of the unconstrained version of this theorem was given (when K = I) by Overton [26], following work of Fletcher [14]. The n × n matrix Û is known as a “dual matrix” by analogy with “dual variables” (Lagrange multipliers) familiar from mathematical programming. The t × t matrix U may be called a “reduced dual matrix”, but since it is the one we shall need as a computational tool we shall also refer to it as simply the dual matrix. The distinction between Û and U is analogous to the notational question of whether inactive constraints in a nonlinear program should be assigned zero Lagrange multipliers.

In the case that t = 1 and the bound constraints are inactive, the necessary condition reduces to the requirement that the gradient of λ_{1}(σ), whose elements are q_{1}^{T} B_{k}(σ)q_{1}, has the constant value ν. (Here q_{1} is the only column in Q_{1}, and U is the scalar one.) In the case that t = 2, let

U = ( \begin{matrix}δ_{1} & δ_{3}/2 \\ δ_{3}/2 & δ_{2} \end{matrix} ) ,  (6.11)

let the two columns of Q_{1} be q_{1} and q_{2} and again assume that all bound constraints are inactive. The necessary condition then becomes

δ_{1}q_{1}^{T} B_{k}(σ)q_{1} + δ_{2}q_{1}^{T} B_{k}(σ)q_{2} + δ_{3}q_{1}^{T} B_{k}(σ)q_{2} = ν

together with the trace and positive semidefinite constraints on U. Without loss of generality, δ_{1} and δ_{2} may be taken to have nonnegative sign, and the normalizing trace condition may be replaced by the assumption that ν = 1. The positive semidefinite constraint is then simply

δ_{1}δ_{2} ≥ δ_{3}^{2}/4.

This is the same necessary condition given by Bratus and Seiranian [6] and Masur [22]. We note that the derivation given here not only applies for t > 2, but is much simpler than that given by [6] and [22] for the case t = 2. In a footnote, Masur conjectured that the positive semidefinite condition on U would also be the correct necessary condition for t > 2. Bratus [7] gave a discussion of necessary and sufficient conditions for general multiplicity t, but the given necessary condition concerns the necessary sign of the directional derivative of λ_{1} for all feasible directions; the positive semidefinite condition on U was apparently not obtained.

Before discussing the algorithm that springs from Theorem 6.1 we investigate the extent to which it suggests a new tack on the infinite dimensional problem. Regarding the variational principle of (6.9) we consider K^{+}(X), the space of positive compact linear operators on a real separable Hilbert space X. Each T ∈ K^{+}(X) possesses a countable sequence of eigenvalues λ_{1}(T) ≥ λ_{2}(T) ≥ ⋯ ↓ 0 repeated according to multiplicity and a (possibly infinite) trace tr T = ∑_{i=1}^{∞} λ_{i}. In this context it is not difficult to show for symmetric T ∈ K^{+}(X) that

λ_{1}(T) = max {tr T U; U ∈ K^{+}(X), tr U = 1}.

Recall that u ↦ (σ^{p}u′′)′′ and u ↦ −u′′ are positive symmetric isomorphisms of H_{0}^{2} onto H^{−2} and H_{0}^{1} onto H^{−1} respectively. We denote these maps by A and B, remark that B^{1/2} is a positive isomorphism of H_{0}^{1} onto L^{2}, and denote by I the compact imbedding of H_{0}^{2} in H_{0}^{1}. With ϕ = B^{1/2}Ju, and ^{∗} denoting adjoint, the buckled column equation receives the formulation

(1/λ)ϕ = B^{1/2}IA^{−1}(B^{1/2}I)^{∗}ϕ.

By construction, B^{1/2}IA^{−1}(B^{1/2}I)^{∗} is a symmetric member of K^{+}(L^{2}). Although we may now proceed to compute ∂λ_{1}^{−1} as in the previous theorem, this representation suffers from its dependence on the unknown A^{−1} and B^{1/2} in contrast to Theorem 4.3 that works directly with A and B.

We now turn to the question of how to solve the finite dimensional optimization problem. Although there is a substantial literature on the generalization of gradient methods to nonsmooth problems (see the survey [28] by Polak), little attention has been given to applying nonsmooth optimization techniques to (6.3). One exception is [29], which describes an algorithm for maximizing the least eigenvalue of a variety of important structures, accounting for the presence of multiple eigenvalues. It focuses however on the clamped vibrating column, (σ^{p}u′′)′′ = λu, u ∈ H_{0}^{2}, a problem which has long been known to admit only simple eigenvalues (see e.g., Leighton and Nehari [21, Lemma 4.1]). Our algorithm differs from that of [29] in our attention to the added structure of the generalized gradient of λ_{1} as revealed in the theorem above. We use an algorithm specifically designed to exploit this structure, which is based on Overton [26], but modified to be far more efficient for moderate to large mesh size N.

Given σ ∈ ad^{h} with λ_{1}(σ) and λ_{2}(σ) the two least eigenvalues of (6.2) we normalize the corresponding eigenvectors q_{1} and q_{2} so that Q_{1} = [q_{1} q_{2}] satisfies Q_{1}^{T} KQ_{1} = I, the 2 × 2 identity matrix. These eigenvalues and eigenvectors are computed by subspace iteration with a block size of two, with the necessary linear systems solved directly using the Cholesky factors of B(σ) (see Bathe and Wilson [5] for details).

Define the approximate multiplicity t of λ_{1} by t = 2 if

λ_{2}(σ) − λ_{1}(σ) ≤ τλ_{1}(σ)

and t = 1 otherwise. Here τ is a positive tolerance which may be adjusted during the optimization process. A multiplicity higher than two is excluded (for sufficiently small h) by Lemmas 2.1 and 5.1. Now consider the following linear program (LP):

max_{d∈R^{N}} d_{N}  (LP0)

subject to

Ed = e  (LP1)
F d ≤ f  (LP2)
d_{k} = 0, k ∈ J  (LP3)
α − σ_{k} ≤ d_{k} ≤ β − σ_{k}, k = 1,… N − 1  (LP4)
|d_{k}| ≤ ρ, k = 1,…, N  (LP5)

where

(i) The first N − 1 components of d represent proposed changes to σ, while the last component approximates the corresponding change in λ_{1}(σ). Let us write d = [η^{T} ω]^{T}, with η ∈ R^{N−1}, ω ∈ R; (ii) ρ is a positive scalar, whose purpose is to ensure ∥d∥ is not too large; (iii) J is an index set, which effectively removes the corresponding variables from the linear program; (iv) The first row of the matrix E is [1,…, 1, 0], and the first element of the right-hand side vector e is 0. This ensures that the changes to σ respect the integral constraint;

(v) The second row of E is

[−q_{1}^{T} B_{1}(σ)q_{1},…, −q_{1}^{T} B_{N−1}(σ)q_{1}, 1]

and the corresponding element of e is 0. If t = 2, then E contains an additional two rows,

[−q_{2}^{T} B_{1}(σ)q_{2},…, −q_{2}^{T} B_{N−1}(σ)q_{2}, 1] and [−q_{1}^{T} B_{1}(σ)q_{2},…, −q_{1}^{T} B_{N−1}(σ)q_{2}, 0]

with corresponding elements of e set to λ_{2}(σ) − λ_{1}(σ) and 0 respectively; (vi) If t = 1, F contains the single row

[−q_{2}^{T} B_{1}(σ)q_{2},…, −q_{2}^{T} B_{N−1}(σ)q_{2}, 1]

and f is the scalar λ_{2}(σ) − λ_{1}(σ). If t = 2, F and f are empty, i.e. (LP2) may be removed.

We note that the given rows may be computed very efficiently, since the derivative matrices B_{k}(σ) are extremely sparse.

The justification for (v) and (vi) is as follows. If t = 1, the second row of E imposes a linearization of the nonlinear equation λ_{1}(σ + η) = λ_{1}(σ) + ω, while the only row in F imposes a linearization of the inequality λ_{2}(σ+η) ≥ λ_{1}(σ)+ω. Since ω is being maximized, the solution to the linear program yields the steepest ascent direction for λ_{1}(σ), projected to satisfy the integral and bound constraints, with steplength required to be short enough that the linearized value for λ_{2}(σ + η) does not drop below that for λ_{1}(σ + η), and that the various bounds are satisfied. If t = 2, the second through fourth rows of E give a linearization of the appropriate set of three nonlinear equations imposing the coalescence of λ_{1}(σ +η) and λ_{2}(σ +η), see [16]. The common linearized value, λ_{1}(σ)+ω, is maximized, subject to the given constraints.

Theorem 6.2. Suppose that τ = 0 so that the multiplicity estimate t is exact, and suppose that ρ > 0 and J is the empty set. Then d = 0 is a (nonunique) solution to the linear program given above if and only if (6.4) holds for some U ∈ S^{t} with tr U = 1, some ν ∈ R and some γ ∈ R^{N−1} satisfying (6.5-7).

Proof: By the usual Lagrange multiplier rule, the linear program admits the solution d = 0 if and only if there exist multipliers ν ∈ R, δ ∈ R^{t(t+1)/2} and γ ∈ R^{N−1} satisfying

E^{T} ( \begin{matrix}ν \\ δ \end{matrix} ) + ( \begin{matrix}γ \\ 0 \end{matrix} ) = ( \begin{matrix}0 \\ 1 \end{matrix} )

with γ subject to the standard sign condition. Setting U = δ_{1}(= 1) if t = 1 and defining U by (6.11) if t = 2, we have (6.4-7) with tr U = 1. The same argument holds in the reverse direction.

Note two points: there is no positive semidefinite condition obtained on U, and the solution d = 0 is generally not at a vertex of the feasible region, so is not unique.

Our algorithm for solving (6.3) generates a sequence in ad^{h}. Each successive approximation is obtained from the previous one by consideration of a linear program (LP) of the form given above. We first define a simple version of the algorithm, but one which is too costly for practical use. In this version, we obtain σ + η, a candidate replacement for σ, by solving the LP for d = [η^{T} ω]^{T}. If λ_{1}(σ + η) > λ_{1}(σ), σ is replaced by σ + η and the process repeated. Otherwise, σ remains unchanged, the trust region radius ρ is decreased by a factor of two, and the revised LP is considered. This kind of trust region approach can be made very effective by modifying the size of ρ according to how well the actual increase in λ_{1}(σ) agrees with the linear prediction ω, as discussed in Fletcher [14] in the context of general nonlinear programming. As recommended by Fletcher, we double ρ if the ratio of actual to predicted increase exceeds 0.75 and halve ρ if the ratio is less than 0.25. The process is terminated when ∥d∥ ≤ ϵ, a convergence tolerance.

However, the expense of obtaining the optimal solution of each linear program is not justified. Although the “limit” LP defined by σ equal to a solution of (6.2-4) has an optimal solution which is not a vertex, generically, any LP solved during the successive approximation process can be expected to have a unique solution which must be at a vertex. Since only a few constraints involve all the variables, most of the constraints defining a vertex are simple bounds, and most of these may be trust radius bounds of the form (LP5). Since the only purpose of the trust radius bounds is to avoid taking steps too large for the linearizations to be accurate, there is little to be gained by finding the exact set of active bounds. We therefore partially solve the LP as follows. We first attempt to obtain a feasible point for the LP by setting d to the least norm solution of (LP1,LP3), contracting this step if necessary to satisfy the various inequalities. This contraction effectively scales the right-hand side of the only possible inhomogeneous equality constraint in the LP, that corresponding to the third row of E in the case t = 2. The rationale here is that if the least norm step satisfying the equality constraints is not feasible, the underlying approximations are probably not good enough to justify the solution of the unmodified LP. We then start the LP solution process as in a projected gradient method, augmenting d by projected gradient steps with steplengths determined by the inequality constraints and bounds. The gradient being projected is that of the LP objective, i.e. the vector [0,…, 0, 1], while the constraints determining the projection are the equality constraints and any inequality constraint and bounds encountered during the process. This continues until either (a) a trust radius bound of the form (LP5) is encountered, or (b) the norm of the current projected gradient increment drops below the tolerance ϵ (unlikely to happen first). At this point the process of partially solving the LP is terminated. Anywhere from zero to many active bounds of the form (LP4) may be encountered by this process, as well as, possibly, the inequality (LP2) (in the case t = 1). By adding any active bounds encountered to the set J, we avoid having to process these bounds again during the (partial) solution of subsequent LP’s. However, the signs of the associated bound multipliers must be checked after the partial LP solution and bounds with the wrong sign removed from J if necessary. The entire process is very efficient, requiring QR factorizations of matrices with only two to four columns, with rows removed corresponding to active bounds. For complete details of the process, see [27].

In the case t = 2, the LP partial solution process generates four multipliers corresponding to the rows of E, namely ν, δ_{1}, δ_{2}, and δ_{3}. If the corresponding dual matrix U, defined by (6.11), is not positive semidefinite, this is a clear indication that the multiplicity estimate t is incorrect, and so τ is reduced by a factor of ten. In principle, it might be necessary to use a more sophisticated technique to recover from a multiplicity estimate which is too large. For example, if the algorithm was started at a point where all the optimality conditions except U ≥ 0 were satisfied, it would be necessary to split the multiple eigenvalue to obtain an ascent direction; this is explained further in [26] and [27]. However, this technique has not been required in our computational experiments for the Lagrange problem.

This completes our outline of the algorithm used to generate the numerical results given in the next section. For more algorithmic details, see [27]. We do not have any proof that the given algorithm will converge to a solution of (6.3), but given any approximate solution we may verify the required signs of γ and the eigenvalues of the dual matrix U, and compute the residual of the approximate equation (6.4). We have found the algorithm to be very effective in practice as the numerical results attest.

7. Computational Results

The algorithm outlined in the previous section has been implemented in Fortran 77 and tested extensively. Subroutines from the Linpack [13] and Eispack [17] libraries were used to (a) perform the QR factorizations required during the partial LP solution process (for matrices with at most four columns), (b) obtain the Cholesky factorizations of B(σ) needed for subspace iteration (these matrices have only seven nonzero diagonals) and (c) solve the reduced generalized eigenvalue problems required for subspace iteration (these have order two). Parameters were set as follows: τ, the relative multiplicity tolerance, was initialized to 0.1; ρ, the trust radius, was initialized to 5.0; ϵ, the convergence tolerance, was set to 10^{−3}. The initial σ was set to the constant one, corresponding to the uniform column. Runs were made for various values of N, the number of mesh points; p, the power of σ in the differential equation; α and β, the lower and upper bounds on σ, and the various homogeneous boundary conditions: clamped–clamped, clamped–hinged, clamped–free and hinged–hinged.

The algorithm was found to be very efficient, typically invoking subspace iteration, in the computation of the two least eigenvalues of (6.2), about 50 times prior to reaching the convergence tolerance. At the final iterate the residual of the approximate equation (6.4) was typically found to have norm about 10^{−3}. Smaller choices of ϵ required significantly more computation but did not produce a qualitatively improved solution. Other initial choices of ρ affected only the initial efficiency of the algorithm. The results were relatively insensitive to the initial choice of the multiplicity tolerance τ, although smaller choices delayed identification of the final multiplicity. There was usually no difficulty in determining the correct final multiplicity t, with corresponding positive semidefinite dual matrix U. In the cases where the final multiplicity t was two the gap between λ_{1} and λ_{2} was typically reduced to 10^{−6}. The subspace iteration was itself very efficient, requiring only one or two steps on all but the first few steps of the optimization, reflecting the good separation of λ_{2} from λ_{3} and the availability of an excellent initial two-dimensional subspace, namely the span of the eigenvectors q_{1} and q_{2} computed at the previous optimization step. (The first subspace iteration was initialized using the first two columns of the identity matrix.) The initial σ in each case was that of the uniform column, σ ≡ 1. Symmetry was not imposed on the algorithm’s subsequent choices of σ. A typical run for N = 513 required 1.5 hours on a Sparc station.

We begin our summary of the results with p = 2. Under the assumption that each transverse cross–section of the column is circular we recall that σ is proportional to the square of the cross–section’s radius. Plotting both ±\sqrt{σ} then gives a lengthwise crosssection of the associated column. With this representation one may then view the corresponding buckling mode(s) simultaneously. Our figures, generated by Matlab [23], portray the column in the piecewise fashion produced by the algorithm of §6 while using dashed curves(s) to indicate the buckling mode(s). We remark that for those optimal designs with double buckling loads the corresponding buckling modes depend on our initial choice of subspace, in subspace iteration, and σ.

Figure 1 gives our strongest clamped–clamped column and its first two buckling modes. Here p = 2, α = 0, β = 10, with a double buckling load of 52.3533. This value agrees to four figures with that obtained by Olhoff and Rasmussen [25], Masur [22], and Seiranian

On increasing α or decreasing β these bound(s) will eventually become active. Figure 2 gives our strongest clamped–clamped column and its first two buckling modes when p = 2, α = 0.25, and β = 10. The least eigenvalue is still double, though now reduced to 52.3467.

As the uniform column has a simple first eigenvalue one expects that a sufficient increase in α would produce an optimal design with a simple first eigenvalue. Figure 3 gives our best strongest clamped–clamped column and its first buckling mode when p = 2, α = 0.5, and β = 10. In this case the first two eigenvalues are 51.07086 and 62.3479.

Clearly there must exist (at least) one value of α between 1/4 and 1/2 at which the optimal buckling load switches multiplicity. Olhoff and Rasmussen [25] declare 0.28 to be the only such value. Our algorithm also indicates the presence of such a critical α in the vicinity of 0.28. We note that in addition to being able to approach 0.28 from above – proceeding until the gap between the least two eigenvalues vanishes, we have also approached from below, in this case proceeding until the least eigenvalue of the corresponding 2 by 2 dual matrix vanishes.

Figure 4 gives our strongest clamped-hinged column and its first buckling mode when p = 2, α = 0.25, and β = 10. The buckling load of 27.0762 is, as expected, simple. Although decreasing α increases the buckling load, our designs converge, as α → 0, to the Tadjbakhsh and Keller solution, (5.13). As shown in the appendix, this column buckles at π^{2}/3, and so can not possibly be optimal. The convergence of our algorithm to (5.13) only strengthens our belief that the problem, as stated by Tadjbakhsh and Keller, is without a solution. That is, σ ↦ μ_{2,1}(σ) with p = 2, does not attain its maximum on ad when α = 0.

Our numerical results also agree with Tadjbakhsh and Keller in the cases for which we have argued that they are correct. In particular, Figure 5 gives our strongest clamped–free column and first buckling mode when p = 2, α = 0, β = 10. The buckling load, again simple, is 3.2897. Figure 6 gives our strongest hinged–hinged column and first buckling mode when p = 2, α = 0, β = 10. Its simple buckling load is 13.1579.

We return to the clamped–clamped case and consider its dependence on p. Our analysis of (4.5-7) led us to believe that, for p > 1, the minimum (maximum) of the optimal design increases (decreases) with p. This is reinforced by Figure 7, whose lower (upper) curve traces the minimum (maximum) of the optimal design as a function of p. As the buckling load is double for each of these designs there must exist a curve, between the lower one and the curve that is constantly one, across which the optimal buckling load changes multiplicity. We have seen that (2, 0.28) lies near such a curve. With respect to the range of p considered in Figure 7 we have found that both the optimal buckling load and the least eigenvalue of its corresponding dual matrix increase with p. Regarding the behavior as p tends to 1 from above we have found that the minimum of the optimal design tends to zero, and, though the optimal buckling load remains double, the least eigenvalue of the corresponding dual matrix tends to zero. Below p = 1 we found optimal designs with simple buckling loads regardless of our choice of α. Figure 8 gives our strongest clamped–clamped column when p = 1, α = 0, and β = 10. The buckling load of 47.9898 and the design itself are very close to the analytical result of (4.15). Refining the mesh in neighborhoods of 1/4 and 3/4, and perhaps using piecewise linear elements for σ, would presumably bring us even closer to (4.15). We have not pursued this for two reasons. First, we argued in §4 that in the absence of an existence proof one can not fully trust (4.15), and second, in both of the physical contexts that p = 1 has arisen there is an a priori strictly positive lower bound on the admissible σ. Regarding the latter, we present in Figure 9 our strongest clamped–clamped column when p = 1, α = 0.8, and β = 1.2. Its simple buckling load is 43.4921.

Having addressed dependence on p and α at a particular level of discretization we now fix p = 2, α = 0, β = 10 and with clamped–clamped boundary conditions demonstrate the convergence of several relevant parameters as N, the number of mesh points, becomes large. In particular, Table 1 lists μ̂_{2,2} (the optimal buckling load), the least eigenvalue of the associated dual matrix U, and ∥σ̂_{N} − σ̂_{1025}∥_{∞} (the greatest difference between the optimal design on a mesh of N points and the optimal design on a mesh of 1025 points) for values of N from 65 to 1025.

\begin{matrix}N & μ̂_{2,2} & min ev(U) & ∥σ̂_{N} − σ̂_{1025}∥_{∞} \\ 65 & 52.14944 & 0.023859 & 0.1066 \\ 129 & 52.31027 & 0.043317 & 0.0415 \\ 257 & 52.33615 & 0.047034 & 0.0424 \\ 513 & 52.35332 & 0.046435 & 0.0059 \\ 1025 & 52.35548 & 0.046607 & 0.0000 \end{matrix}

We close our study with a glance at the numerical range of the buckling load over ad. To this point we have concentrated on its maximization, and, though we may compare this value to that of the associated uniform column, it would be of interest to weigh it against the minimum buckling load. Clearly, α must now be strictly positive, for one could produce a buckling load of zero by prescribing that σ vanish on some interval. Regarding the existence of a minimizer for σ ↦ μ_{i,j}(σ) over ad we note that Theorem 3.1 is insufficient. Recall in (3.3) that the limit of the maximizing sequence integrated to less than one. This was not an obstacle, for adding mass could only increase the buckling load. As our goal now is to minimize this load it appears that one must either relax the cost functional or construct the so–called G–closure of ad to obtain a minimizing design. Instead of embarking on this we modified the strongest column algorithm to minimize instead of maximize σ ↦ μ_{i,j}(σ).

The modification to the algorithm is very simple, namely changing the sign of (LP0) and requiring a decrease instead of an increase in the smallest eigenvalue. The modified algorithm generated plausible weakest designs in ad, and, though we lack a proof of optimality, we shall content ourselves with a discussion of these numerical results. In all cases (independent of p and α) the minimum buckling load was simple; this is to be expected since the minimization should tend to separate the least eigenvalue from the remainder of the spectrum. In addition, we find that the generated designs have their mass concentrated near the inflection points of their associated positive buckling mode. This too is to be expected, making the opposite argument to that made in §2.

Figure 10 gives our weakest clamped–clamped column with p = 2, α = 0.25, and β = 10. This column buckles under a simple axial load of 2.5658. The buckling load of our strongest column in this class (see Fig. 2) is 52.3467.

Figure 11 gives our weakest clamped–clamped column with p = 1, α = 0.8, and β = 1.2. This column buckles under the simple axial load of 33.5631. The buckling load of our strongest column in this class (see Fig. 9) is 43.4921.

8. Concluding Remarks

Having thoroughly discussed the work of Tadjbakhsh and Keller [37] we now take up two related issues, that, though discussed by Keller in [20], have received little rigorous attention since.

The first involves the optimal design of cylindrical columns. Here, given again a fixed amount of material to be distributed over a column of fixed length, we seek the shape of the cross–section that when used to generate a cylinder produces a column with the greatest buckling load. One is not allowed to “taper” the column as we have in past sections. Keller quickly reduced this problem to the search for that planar domain of fixed area with the greatest least second moment of area. Recall that the second moment of area of the domain Ω with centroid at the origin in the direction η (with |η| = 1) is

I(Ω, η) ≡ ∫_{Ω} |η^{T} y|^{2} dy.  (8.1)

Denoting the unit circle by S, Keller’s problem takes the form

sup_{|Ω|=A} inf_{η∈S} I(Ω, η).  (8.2)

Keller noted the existence of Ω for which this value is infinite. To exclude such Ω he restricted himself to convex domains. Within this smaller class he then argued, without proof, that the equilateral triangle possesses the greatest least second moment of area. Citing Pólya, Truesdell later observed that Keller’s conjecture was indeed true. In particular, Pólya in [30] found that among convex sets the maximum of 4I_{1}I_{2}|Ω|^{−4} occurs for triangles. Here I_{1} and I_{2} are the respective principal second moments of area, i.e., the min and max of (8.1). Keller’s result follows on noting that only sets for which I_{1} = I_{2} need be considered, for if I_{1} < I_{2} one can simply redistribute the mass in such a way that I_{1} increases while I_{2} decreases.

This reduction to convex domains is too severe. It is not the lack of convexity that allows (8.2) to grow without bound but the possibility that Ω itself may be unbounded, though of finite area. To exclude this behavior one may simply bound |∂Ω|, the length of Ω’s boundary. That this does indeed bound I(Ω, η) follows from the isoperimetric inequality

2I_{p}(Ω)/π ≤ (|∂Ω|/2π)^{4}

of Pólya and Szegö [31, §1.5], where I_{p}(Ω) is the polar moment of inertia. As the second moment of area will be independent of η for the best Ω, its value will be one half that of its polar moment of inertia. We must now consider,

sup_{{|Ω|=A \atop |∂Ω|≤L}} inf_{η∈S} I(Ω, η).  (8.3)

For fixed A the value of (8.3) will grow as one increases L, suggesting that for sufficiently large L one may produce a domain whose second moment of area exceeds that of the equilateral triangle of the same area. We shall accomplish this with A = 1 and L = 16. To produce large values in (8.3) one need only consider domains that are symmetric about the coordinate axes as well as the two diagonals and that stretch out towards infinity. The symmetry will render η ↦ I(Ω, η) constant, while the latter condition will ensure us that I(Ω, η) is large. The domain of Fig. 12 has an area of one, a boundary whose length does not exceed 16, and a second moment of area of \frac{16}{45} + \frac{71}{8^{3}90}. This value is more than three times greater than \frac{1}{6\sqrt{3}}, the second moment of area of the equilateral triangle of the same area. Though we have not solved (8.2), this example demonstrates that (8.2) produces, through the designer’s choice of L, columns with arbitrarily large buckling load.

The other issue we wish to resurrect from Keller [20] is that of post-buckling. Analysis of the buckled column must proceed from the nonlinear model. For example, in the hingedhinged case one considers

(σ^{p}ϕ′)′ + λ sin ϕ = 0, ∫_{0}^{1} sin ϕ ds = 0, ϕ′(0) = ϕ′(1) = 0,  (8.3)

where ϕ measures the angle between the column and a fixed axis in its plane of buckling. Equation (5.1) arises from linearizing sin ϕ to ϕ, identifying u′ = ϕ, and differentiating the differential equation in (8.3). The least eigenvalue of (5.1), μ_{1,1}, is indeed a bifurcation point for (8.3). Moreover, Keller showed that the direction of the solution branch emanating from μ_{1,1} was indeed rightward, i.e., that (8.3) admits only the trivial solution for λ < μ_{1,1}. It remains to determine the nature of the nontrivial solution branch(s) for the other sets of boundary conditions. Here one would also like to understand the role of imperfections in the model and/or design. In particular, the splitting of the double eigenvalue in the optimal clamped–clamped column via an unfolding of the ideal bifurcation problem in a parameter that introduces material asymmetry would be of interest.

Appendix

We have argued throughout that the necessary and sufficient conditions proposed by Tadjbakhsh and Keller are incorrect. This does not in itself however, invalidate their designs. Indeed, we argued that their solutions to the hinged–hinged and clamped–free problems are correct. This appendix serves to demonstrate that these are their only correct designs.

In particular, we show that Tadjbakhsh and Keller incorrectly calculated the buckling loads of their proposed solutions to the clamped-clamped and clamped–hinged problems. Recall their solution of the former,

A(x) = \frac{4}{3} sin^{2} θ(x), −π/2 ≤ θ ≤ 3π/2, (A.1)

θ(x) − \frac{1}{2} sin 2θ(x) + π/2 = 2πx, 0 ≤ x ≤ 1.

As Olhoff and Rasmussen [25] observed in their numerical work, the column constructed according to (A.1) tends to deform (under axial compression) on (0, 1/4) and (3/4, 1), with the center of the column experiencing only a rigid motion. To make this precise we recall from [37] that this A, when restricted to (0, 1/4), is actually the optimal design of the clamped–free problem there, with a corresponding buckling load of π^{2}/3 (as the volume of this piece is also 1/4) and first eigenfunction u_{2,0}. We normalize u_{2,0} so that

∫_{0}^{1/4} |u′_{2,0}|^{2} dx = 1,

and use it to construct a sequence {φ_{n}} ⊂ H_{0}^{2} on which

\frac{∫_{0}^{1} A^{2}|φ′′_{n}|^{2} dx}{∫_{0}^{1} |φ′_{n}|^{2} dx} → \frac{π^{2}}{3} as n → ∞.  (A.2)

First we build the continuous displacement

ũ(x) = \{ \begin{matrix}u_{2,0}(x), & if 0 ≤ x ≤ 1/4 \\ u_{2,0}(1/4), & if 1/4 ≤ x ≤ 3/4 \\ u_{2,0}(1 − x), & if 3/4 ≤ x ≤ 1. \end{matrix}

This function is not a member of H_{0}^{2} but we shall see that it suffices to introduce the cubic perturbation

p_{n}(x) = a(x − 2nx^{2} + n^{2}x^{3}), where a = u′_{2,0}(1/4),

near its singularities. As p_{n}(0) = p_{n}(1/n) = p′_{n}(1/n) = 0, the function

φ_{n}(1 − x) = φ_{n}(x) = \{ \begin{matrix}ũ(x), & if 0 ≤ x ≤ 1/4 \\ ũ(1/4) + p_{n}(x − 1/4), & if 1/4 ≤ x ≤ 1/4 + 1/n \\ ũ(1/4), & if 1/4 + 1/n ≤ x ≤ 1/2 \end{matrix}

possesses a continuous derivative. It remains to show that φ′′_{n} ∈ L^{2}. The only possible obstacle is the behavior of u′′_{2,0} near 1/4. Returning to Tadjbakhsh and Keller [37,§3 (25)] we find

u′′_{2,0}(x) = \frac{\sqrt{3}}{2 sin θ(x)}.

And, as (A.1) implies that θ(x) = O(|x − 1/4|^{1/3}) as x → 1/4 we find u′′_{2,0}(x) = 0(|x − 1/4|^{−1/3}). This singularity is clearly square integrable hence φ_{n} ∈ H_{0}^{2} and one can consider the Rayleigh quotient

\frac{∫_{0}^{1} A^{2}|φ′′_{n}|^{2} dx}{∫_{0}^{1} |φ′_{n}|^{2} dx} = \frac{2 ∫_{0}^{1/4} A^{2}|u′′_{2,0}|^{2} dx + 2 ∫_{0}^{1/n} A^{2}(x − 1/4)|p′′_{n}|^{2} dx}{2 ∫_{0}^{1/4} |u′_{2,0}|^{2} dx + 2 ∫_{0}^{1/n} |p′_{n}|^{2} dx}.  (A.3)

By construction

∫_{0}^{1/4} A^{2}|u′′_{2,0}|^{2} dx = π^{2}/3 and ∫_{0}^{1/4} |u′_{2,0}|^{2} dx = 1.

It remains to show that the remaining terms in (A.3) tend to zero with increasing n. Beginning with the numerator, from A(x) = O(|x − 1/4|^{2/3}) comes

∫_{0}^{1/n} A^{2}(x − 1/4)|p′′_{n}(x)|^{2} dx = ∫_{0}^{1/n} O(x^{4/3})|2an(3nx − 2)|^{2} dx = O(1/n^{1/3}),

while in the denominator

∫_{0}^{1/n} |p′_{n}(x)|^{2} dx = ∫_{0}^{1/n} |a(1 − 4nx + 3n^{2}x^{2})|^{2} dx = O(1/n).

As we have constructed a sequence of admissible displacements whose Rayleigh quotients tend to π^{2}/3 we conclude that the clamped–clamped column built according to (A.1) buckles at a load not exceeding π^{2}/3.

Now recall Tadjbakhsh and Keller’s solution to the clamped–hinged problem,

σ̂_{2,1}(x) = \frac{4 sin^{2} θ(x)}{3 sin^{2} θ(0)}, θ(0) ≤ θ ≤ π, (A.4)

θ(x) − \frac{1}{2} sin 2θ(x) + \frac{1}{2} sin 2θ(0) − θ(0) = x(π + \frac{1}{2} sin 2θ(0) − θ(0)), 0 ≤ x ≤ 1,
\frac{1}{2} sin 2θ(0) − θ(0) = −\frac{2}{3} sin^{3} θ(0) cos^{−1} θ(0) − π,

and the fact that it vanishes at x_{0} = y/(π + y) where y = \frac{1}{2} sin 2θ(0) − θ(0). Analogous to the above, this design is optimal for the clamped–free column on (0, x_{0}). The volume of this piece being x_{0} as well, we find that it buckles at π^{2}/3. Denoting by u_{2,0} the corresponding clamped–free eigenfunction on (0, x_{0}) whose derivative has L^{2} norm 1, we define,

ũ(x) = \{ \begin{matrix}u_{2,0}(x), & if 0 ≤ x ≤ x_{0} \\ \frac{u_{2,0}(x_{0})}{1−x_{0}}(1 − x), & if x_{0} ≤ x ≤ 1. \end{matrix}

Here we introduce the perturbation

p_{n}(x) = bx − 2n(b + c)x^{2} + n^{2}(b + c)x^{3}, where b = u_{2,0}(x_{0}), c = u′_{2,0}(x_{0}),

and the corresponding regularization

φ_{n}(x) = \{ \begin{matrix}ũ(x), & if 0 ≤ x ≤ x_{0} \\ ũ(x_{0}) + p_{n}(x − x_{0}), & if x_{0} ≤ x ≤ x_{0} + 1/n \\ ũ(x), & if x_{0} + 1/n ≤ x ≤ 1. \end{matrix}

By construction, φ_{n} ∈ C^{1}, while as in the clamped–clamped case u′′_{2,0} behaves like |x − x_{0}|^{−1/3} near x_{0}, and so φ_{n} ∈ H_{0}^{2}. Moreover, σ̂_{2,1}(x) = O(|x − x_{0}|^{2/3}) implies as above that

∫_{0}^{1/n} σ̂_{2,1}^{2}(x − x_{0})|p′′_{n}(x)|^{2} dx = O(1/n^{1/3}), and ∫_{0}^{1/n} |p′_{n}|^{2} dx = O(1/n).

Consequently,

\frac{∫_{0}^{1} σ̂_{2,1}^{2}|φ′′_{n}|^{2} dx}{∫_{0}^{1} |φ′_{n}|^{2} dx} → \frac{π^{2}/3}{1 + u_{2,0}^{2}(x_{0})/(1 − x_{0})} as n → ∞,

i.e., the clamped–hinged column built according to (A.4) is even weaker than the clamped– clamped column of (A.1).

Acknowledgements

This work began during the first author’s visit to the Courant Institute. The efforts of Bob Kohn and Cathleen Morawetz are here greatly appreciated. We thank Andre Cherkaev for introducing us to the work of Bratus and Seiranian.

References

[1] F. V. Atkinson, Discrete and Continuous Boundary Problems, Academic Press, 1964.

[2] G. Auchmuty, Dual variational principles for eigenvalue problems, in Nonlinear Functional Analysis and its Applications, F. Browder, ed., AMS, 1986, pp. 55–71.

[3] V. Barbu and T. Precupanu, Convexity and Optimization in Banach Spaces, Reidel, Boston, 1986.

[4] D. Barnes, The shape of the strongest column is arbitrarily close to the shape of the weakest column, Q. Applied Math., XLVI(1988), pp. 605–609.

[5] K. Bathe and E. Wilson, Numerical Methods in Finite Element Analysis, Prentice– Hall, Englewood Cliffs, NJ, 1976.

[6] A. Bratus and A. Seiranian, Bimodal solutions in eigenvalue optimization problems, Prikladnaia Matematika i Mekhanika, 47(1983), pp. 546–554. Applied Math. Mech., 47(1983), pp. 451–457.

[7] A. Bratus, Multiple eigenvalues in problems of optimizing the spectral properties of systems with a finite number of degrees of freedom, Zh. vychisl. Mat. mat. Fiz., 26(1986), pp. 645–654. USSR Comp. Math. and Math. Phys., 26(1986), pp. 1–7.

[8] J. Cea and K. Malanowski, An example of a max-min problem in partial differential equations, SIAM J. Control, 8(1970), pp. 305–316.

[9] K. K. Choi and E. J. Haug, E.J., Optimization of structures with repeated eigenvalues, in Optimal Design of Distributed Parameter Structures, Haug, E.J. and Cea, J., eds., Sijthoff–Noordhoff, Leyden, 1981, pp. 219–277.

[10] F. Clarke, Optimization and Nonsmooth Analysis, Centre de recherches mathematiques, Montreal, 1989.

[11] S. J. Cox and J. R. McLaughlin, Extremal eigenvalue problems for composite membranes, II, Applied Math. & Opt., 22(1990), pp. 169–187.

[12] B. Dacorogna, Weak Continuity and Weak Lower Semicontinuity of Non-Linear Functionals, Springer, Lecture Notes in Mathematics 922, New York, 1982.

[13] J. J. Dongarra et al., LINPACK Users Guide, SIAM, Philadelphia, 1978.

[14] R. Fletcher, Semidefinite constraints in optimization, SIAM J. Control and Optimization, 23(1985), 493-513.

[15] R. Fletcher, Practical Methods of Optimization, 2nd ed., John Wiley, New York, 1987.

[16] S. Friedland, J. Nocedal, and M. L. Overton, The formulation and analysis of numerical methods for inverse eigenvalue problems, SIAM J. Numer. Anal., 24(1987), pp. 634– 667.

[17] B. S. Garbow et al., Matrix Eigensystem Routines: EISPACK Guide Extension, Lecture Notes in Computer Science 51, Springer-Verlag, New York, 1977.

[18] E. J. Haug and B. Rousselet, Design sensitivity analysis in structural mechanics. II. Eigenvalue variations, J. Struct. Mech., 8(1980), pp. 161–186.

[19] E. Kamke, Neue Herleitung der Oszillationssätze für die linearen selbstadjungierten Randwetaufgaben zweiter Ordnung, Math. Zeit., 44(1938), pp. 635–658.

[20] J. Keller, The shape of the strongest column, Arch. Rat. Mech. Anal., 5(1960), pp.275– 285.

[21] W. Leighton and Z. Nehari, On the oscillation of solutions of self–adjoint linear differential equations of the fourth order, Trans. of AMS, 89(1958), pp. 325–378.

[22] E. Masur, Optimal structural design under multiple eigenvalue constraints, Int. J. Solids Struct., 20(1984), pp. 211–231.

[23] C. Moler et al., Pro–Matlab Users Guide, MathWorks, Sherborn, MA, 1987.

[24] M. Myers and W. Spillers, A note on the strongest fixed-fixed column, Q. Appl. Math., XLIV(1986), pp. 583–588.

[25] N. Olhoff and S. Rasmussen, On single and bimodal optimum buckling loads of clamped columns, Int. J. Solids Struct., 13(1977), pp. 605–614.

[26] M. L. Overton, On minimizing the maximum eigenvalue of a symmetric matrix, SIAM J. Matrix Anal. Appl., 9(1988), pp. 256–268.

[27] M. L. Overton, Large-Scale Optimization of Eigenvalues, NYU Computer Science Department technical report, no. 505, 1990.

[28] E. Polak, On the mathematical foundation of nondifferentiable optimization in engineering design, SIAM Review 29(1987), pp. 21–89.

[29] E. Polak and Y. Wardi, A nondifferentiable optimization algorithm for structural problems with eigenvalue inequality constraints, J. of Struct. Mech. 11(1983), pp. 561–577.

[30] G. Pólya, More isoperimetric inequalities conjectured and proved, Comm. Math. Helv., 29(1955), pp. 112–119.

[31] G. Pólya and G. Szegö, Isoperimetric Inequalities in Mathematical Physics, Annals of Math. Studies, no. 27, Princeton, 1951.

[32] A. Ramm, Queries, Notices AMS, 29(1982), pp. 326–329.

[33] R. T. Rockafellar, The Theory of Subgradients and its Applications to Problems of Optimization: Convex and Nonconvex Functions, Series on Research and Education in Mathematics 1, Heldermann Verlag, Berlin, 1981.

[34] A. Seiranian, On a problem of Lagrange, Inzhenernyi Zh., Mekhanika Tverdogo Tela, 19(1984), pp. 101–111. 19(1984). Mechanics of Solids, 19(1984), pp. 100–111.

[35] P. Senatorov, The stability of the eigenvalues and eigenfunctions of a Sturm-Liouville problem, Differentsial’nye Uravneniia, 7(1971), pp. 1667–1671. Differential Equations, 7(1971), pp. 1266–1269.

[36] G. Strang and G.J. Fix, An Analysis of the Finite Element Method, Prentice–Hall, Englewood Cliffs, NJ, 1973.

[37] I. Tadjbakhsh and J. Keller, Strongest columns and isoperimetric inequalities for eigenvalues, J. Appl. Mech., 29(1962), pp. 159–164.

[38] I. Todhunter and K. Pearson, A History of the Theory of Elasticity and of the Strength of Materials, I, Cambridge, 1886.