The Tallest Column — A Dynamical System Approach Using a Symmetry Solution
A classic problem, the design of the tallest column, is solved again using a different method. By the use of a similarity solution the equations are transformed and the difficult singularity at the endpoint is peeled away. The resulting autonomous system has a critical point and the solution must be on its stable manifold. The solution is found by starting near the critical point in the direction of the stable manifold, and solving backwards numerically. This removes the need for an iterative integration method that was previously used. The method is shown to work for clamped or hinged boundary condition and can also be used for other problems involving singularities at the endpoints.
Keywords: Tallest column, eigenvalue optimization, singular ODE, similarity solution, stable manifold.
We reexamine a classic eigenvalue maximization problem: The design of the tallest column. This problem has a lineage that originates from Euler’s analysis of column buckling. In  Keller determines the tapering of a thin column with fixed volume which maximizes its buckling load. Later, in , Keller and Niordson solved the problem of designing the tallest free standing column. The equations have a nasty singularity at the top of the column. In  this singularity is neutralized by formulating an equivalent integral equation. The latter is solved by numerical iteration.
More recently, in  and , Cox and McCarthy showed that the operator involved can have a continuous spectrum and they find the tapering of the tallest column using other methods. In this paper, no claims are made as to the validity of the equations involved or the derivation of them, see , ,  and  for a discussion regarding the existence of the optimal design and the control requirements of the design. Nevertheless, Keller and Niordson’s solution is highly plausible since the problem is infinitesimally close to a problem with a discrete spectrum.
First, in Section 2, we derive the boundary value problem (BVP) which is an ordinary differential equation (ODE) system for the deflection of the column and the cross-sectional area, together with boundary conditions (BC). The ODE’s contains the buckling load eigenvalue, to be determined as part of the solution. This derivation follows the derivation in  closely and is included for completeness.
In Section 3 we solve the BVP. First we show that the ODE’s have a scaling symmetry and an exact similarity solution which has the desired asymptotic behavior at the tip of the column. This is where the singularity is. The similarity solution suggests a transformation of variables that “peels away the singularity” at the tip. This transformation yields ODE’s for the new variables which form an autonomous system (AS). In the new variables, the similarity solution is represented by a critical point of the AS. Also, due to another symmetry, the eigenvalue disappears from the equations and now appears only in the BC at the base of the column. The transformed BC dictate that the solution must start on a certain surface in the phase space of the AS and approach the critical point.
Since the critical point of the AS has both stable and unstable manifolds, solving the AS numerically from the base to the critical point is difficult—the solution tends to deviate to a growing solution. Instead, we examine the behavior of the AS near the critical point in order to solve the system backwards, from the tip (the critical point) to the base. By linearizing the AS around the critical point, we identify stable and unstable directions. Starting near the critical point, on the stable manifold, the AS is solved numerically, until the solution intersects the surface that the BC defines. The point of intersection with the surface determines the eigenvalue and the problem is solved. The numerical results agree with those of Keller and Niordson.
In Section 4 the problem is solved for other boundary conditions. This is inexpensive once all the groundwork has been done.
This new solution method applies to other eigenvalue maximization problems whose equations have certain scale-invariant structure. In a companion paper this method determines the tapering of a javelin so the its first mode has the highest frequency of vibration (subject to length and volume constraint).
2 Derivation of the Boundary Value Problem
Consider columns, all of the same volume, clamped at the base and free at the top. We want to find the shape of the tallest column that will not buckle under its own weight. As in the Keller and Niordson papers, we solve the problem for a specific class of permissible designs. We assume that the column is thin, i.e. the characteristic width is much less than the height of the column. In addition, we only allow columns with geometrically similar, equally oriented and convex cross-sections. See figure (1).
We parameterize the column by arclength , measured from the tip of the column, along its center axis. The design information is contained in a single function , the cross-sectional area at . Lastly, we concern ourselves with the bending of the column in a specified plane only. This allows us to specify the configuration of the column using a single function —the angle that the center axis of the column makes with the vertical, measured at point . See Figure (2).
To find the BVP satisfied by for a given cross-section , we write the total energy of the system as a functional of , and compute its variational equations.
2.2 Energy Minimization
The total energy of the column, elastic and gravitational, is a functional of :
Here, is the length of the column, is the gravitational acceleration, is the density of the material, and is the vertical elevation at . The bending modulus, is proportional to ,
Here, is a dimensionless constant determined by the cross-sectional shape and orientation, and is Young’s modulus. For a discussion of the choice of cross-sectional designs and derivation of the bending modulus, see . Since is the vertical elevation of a particle at ,
We adopt dimensionless variables in which , and are measured in the units given by the scaling table:
Here, is the total volume of the column. is related to by
Where is the dimensionless load per unit volume.
2.3 Buckling Load
Recall that is the dimensionless load per unit volume. For a column tapering given by , there is a critical value such that for the only solution of the BVP is . For there are non-trivial solutions as well. This is the buckling phenomenon. For such , a small perturbation of the column from the undeflected state can have a smaller energy than the undeflected one, therefore the zero solution is no longer an energy minimum. The energy minimum is one of the non-zero solutions. For slightly greater than the energy minimizing solution will be small. In this case we can study the linearization of (9) about . The linearization of the BVP (9,10) is a Sturm-Liouville eigenvalue problem:
The following analysis assumes, as in , that the Sturm-Liouville problem determines a discrete sequence of eigenvalues , and corresponding eigenfunctions which are determined up to a multiplicative constant.
2.4 Eigenvalue Maximization
What does the tallest column look like? That is, what is the distribution of cross-sectional area which maximizes the smallest eigenvalue of the linearized BVP (11) subject to the volume constraint (7)?
Since we have no other but we will drop the subscript and use to mean the lowest eigenvalue of the linearized BVP. is a functional of . To maximize subject to the volume constraint, we find the variational derivative and solve for the function so that
Here, is a Lagrange multiplier associated with the volume constraint (7).
To find , we introduce a small variation in the cross-sectional area, . As a result, and will have corresponding small variations and respectively. The linearization of equation (11) in and about a solution , results in equation (14). To reduce clutter, the notation is dropped whenever the variable in the parenthesis is , so mean respectively. Whenever other variables are used they are written explicitly.
The BC for are
This is a linear inhomogeneous BVP for . Its solvability condition determines the functional derivative . To find the solvability condition, we multiply both sides of (14) by and integrate over . Integration by parts, use of BC and rearrangement, leads to
The LHS is zero due to (11), and we get an equation relating and :
Here the integration variable in the last term on the right was changed from to for added clarity. In functional derivative form, equation (15) becomes
Substituting this into (13) results in an integro-differential equation for the maximizing cross-sectional area function and the critical load :
The integro-differential system (11, 12, 16) is still not in a form that we can easily solve. It is convenient to convert it into a system of ODE with BC. Differentiating (16) with respect to gives the ODE
In order to transform equation (11) into a ODE and get rid of the integral, the variable is introduced:
Differentiating (17) gives a BVP for :
Physically, is the amount of volume above the point .
3 Solution of the BVP
The system that we want to solve is
with boundary conditions
Heuristically, the optimal column is expected to taper to a point as approaches 0, so is expected to go to zero as . In addition, by the BC (18). Since and both multiply equation (19), one expects severe difficulties in the numerical solution as approaches 0, and this is indeed the case.
This motivates an analysis of the asymptotic structure of the solution to (19–23) as approaches 0. We need to peel away the singularity before attempting any numerical solution. In the following sections we use a similarity solution to remove this singularity, transforming the equations so that they can be solved numerically.
3.1 Similarity Solution
Equations (19–21) have a similarity solution. (For information on similarity solutions see, for example,  or .) To find it we study the scaling relations between the variables: Specifically, let and be “units” of and respectively. A balance of units in ODE’s (19–22) gives the relationships
|cancels out, and we are left with three equations for , and . In fact, one equation is redundant, and we get a simple relation between and :|
Hence, it is expected that the ODE’s are invariant under a scaling transformation, so that if are solutions, then so are for any . This suggests a similarity solution of the ODE with
|We also look for a similar behavior for :|
|Since the equations are linear in , we can choose , therefore|
|where . Setting the two expressions for equal gives a polynomial equation for :|
The solutions to this equation are Each solution for implies a solution for . By looking at equation (16) we see that must be non-negative. Indeed, when , the LHS is non-negative, and implies . Since is positive, . The values for imply the values for , respectively. Therefore, the only admissible values for are . implies , and this is physically non-interesting: When , , hence the column has no weight. Either there is no gravity, or the density of the material is zero. We want to study the first positive eigenvalue . The solution gives hence, . The similarity solution is found to be
3.2 Peeling Away the Singularity
We use the similarity solution to find a solution of the ODE that satisfies the BC at both ends. Any solution can be written as a product of the similarity solutions and new dependent variables :
Substituting these expressions for and in the differential equations (19–21) results in ODE’s for and . These ODE’s are homogeneous in , hence we use as the independent variable. The ODE for are the autonomous system (AS):
|At the tip, as , the BC (23) transform into|
The ODE’s above, although written in an implicit form, can be rewritten as explicit expressions for as functions of . Note that due to the specific choice of dependent variables, is no longer a parameter in the AS. It only appears in the BC at the base. It is easy to see that is a critical point of the AS (34–36). Setting in (33–32) recovers the similarity solution (28–30).
We seek a solution that satisfies the BC at and approaches the critical point as . This solution satisfy all the BC. A solution that approaches a critical point as its limit at is a solution on the stable manifold of the critical point.
Since the critical point has both a stable an unstable manifold, it is difficult to start on the stable manifold far from the critical point and numerically follow the solution into the critical point. The smallest numerical error will cause the solution to deviate on a growing solution that does not approach the critical point. Instead, we identify the stable manifold, and follow the solution backwards in from the critical point until the BC at are satisfied. This way, numerical errors will not be amplified. In fact, numerical errors will either decay or remain constant when solving in this direction. The following section shows that the stable manifold is one dimensional. Thus, there is only one direction in which to start the solution near the critical point. Solving the AS backwards from the critical point will follow the stable manifold. Since the BC at the base define a plane of co-dimension 1, the stable manifold is expected to intersect this plane. Once the intersection point is found, the solution is fully determined.
3.3 Stable Manifold
To identify the stable and unstable manifolds of the critical point, we linearize the AS around the critical point, resulting in the linear ODE
Here, are the deviations from . (39) is a linear system of ODE’s, we therefore look for solutions of the form
Substituting this form of solution into equation (39) yields a generalized linear eigenvalue problem:
Equation (40) has non-trivial solutions for 4 different values of , and . The solutions are spanned by the vectors given in the table below.
As we are looking for the stable manifold, the only solution of interest is the one with a negative , i.e. solution Readers noticing the conspicuous form of solution are referred to Appendix A.1 for brief discussion of it.
To find the numerical solution of (34–36), we start near the stationary point, on the line tangent to the stable manifold, and solve the AS backwards in . The BC (37) at the base define a surface in the phase space:
3.4 Numerical results
Using Matlab ’s ODE solver
ode45, which uses the fourth-order Runge-Kutta method, the AS was solved.
The solver was used with relative error tolerance of and
an absolute error tolerance of .
The initial value for was
The values for and are taken from . Since the stable manifold is one dimensional there are two distinct options for starting near the critical point. In terms of (41) the two options are choosing to be positive or negative. The solution with positive fails to satisfy the BC at the base. Geometrically, the trajectory is going in the wrong direction in the phase space and does not intersect the plane . With the BC at the base are satisfied (that is, vanishes) at . See figure 3.
Since at the base, , the value for is given from the value of . By this formula . Higher accuracy can be achieved by starting closer to the critical point (smaller ), and setting the solver to a smaller tolerance.
4 Other Boundary Conditions
The fact that appears only in the BC at the base has an interesting implication for this problem: If we want to solve the problem with a different BC at the base, we will have the same analysis at the tip and the same similarity solution. Since the stable manifold is one dimensional, we will have to start in the same way while solving backwards. The only difference will be when to stop.
For instance, consider the problem of the tallest hinged column. A hinged column is a column that instead of being constrained to be vertical at the base, is forced to have zero torque at the base.
Note: The lowest eigenvalue of the hinged column is zero. This is because at no gravity the column can be set at any angle. To find the first non-trivial solution it would appear that we need to add an additional orthogonality condition to make sure that we get the next eigenfunction. A short calculation shows that the boundary conditions and the ODE already guarantee the orthogonality of solutions to the trivial, constant-angle, solution.
Formally, the difference is in the BC. In (10), is replaced by the zero torque condition
|The translation of this BC to is|
The solution to this problem would be the design of the tallest column that is hinged at the base and free at the top. Using the same numerical method and the same initial conditions as before, a numerical solution was found. After , the BC, , is satisfied.
Again, since , we have the value of . For these BC, the value is . The cross-section and “peeled” torque are shown in Figure 4.
The results predict a larger eigenvalue for the building with hinged base than with a clamped base. This might be non-intuitive until one notices that the hinged column has within it a shorter clamped column. This column is also buckling and since it is shorter, its buckling load will be larger than the full length clamped beam. This is also the case with simple loaded beams. The hinged loaded beam will be able to support a larger load than the clamped one. The reason for this non-intuitive result is the assumption that the column actually remains vertical.
The derivation of the equations above relies on the use of a variational derivative in order to find the cross-section that produces a stationary eigenvalue . For this procedure to be valid, the operator, BVP (9, 10) should have discrete eigenvalues. As Cox and McCarthy have already shown, this is not always true. Nevertheless, the Keller and Niordson solution is very plausible. As shown in , adding a small weight to the top of the column and letting the weight tend to zero produces a sequence of optimization problems. These problems converge to the self-weight problem, and the solutions (which in this case must exist) converge to the described solution.
We confined ourselves to demonstrating the use of the self similar solution in order to simplify the numerical solution of the given equations. Regardless of their derivation, the equations themselves require a solution and the method of peeling away the singularity using the similarity solution does away with the need of a iteration scheme. This allows for any standard ODE solver to tackle the remaining problem. This method can be applied to other problems that have singularities. In a companion paper we employ a similar analysis to find the tapering of a javelin so that its lowest vibration mode has the highest possible frequency.
Appendix A Appendix
a.1 Specific Solutions in the Linearized System
We turn our attention to the conspicuous value of the eigenvalue and respective eigenvector which matches exactly with the exponents of the similarity solution. Here we show that this solution of the linearized ODE is to be expected.
|be an ODE. We further assume that|
is a solution to the ODE. Any other solution can be written as
Here, is a diagonal matrix with . In other words, must satisfy
Note that is a solution of this system. Therefore, is a critical point of the system for . The stable and unstable manifolds of this system are important for the type of solution method described in this paper. The following proposition shows the existence of a specific unstable direction about the critical point .
The linearization of (46) around has a solution
In the scaled variable, , this corresponds to an eigenvalue of 1 and an eigenvector .
The linearization of (46) is
|Here, is a diagonal matrix (the perturbation from the solution ) and is the Jacobian matrix of evaluated at . To find out more about , we differentiate the equality|
|with respect to :|
|Equation (49) is derived by substituting (44) into (43). To show that solves equation (48), we substitute (47) into (48):|
|From equation (50) we substitute for the RHS|
This shows that indeed is a solution of the linearized ODE. ∎
a.2 Calculation of the Lagrange multiplier
To calculate the value of in (16), we multiply the equation by and integrate from to :
|Since the RHS of (16) is a constant and , the RHS stays unchanged. In the LHS, integrating the first term by parts, changing the order of integration of the second term and using the BC gives|
|Lastly, using (11) we can see that|
a.3 Torque Balance
Changing the order of integration in the RHS of (51) we see that
The inner integral is the horizontal displacement of the point from the point . The sections of the column above exert a gravity generated torque on the point . This is the RHS of (52). The LHS expresses the elastic response to this torque: The curvature of the beam at is proportional to the imposed torque. See figure (5).
-  G. I. Barenblatt, Scaling, self-similarity, and intermediate asymptotics, Cambridge University Press, New York, 1996.
-  G. W. Bluman and S. C. Anco, Symmetry and integration methods for differential equations, Applied mathematical sciences, vol. 154, Springer-Verlag, 2002.
-  Steven J. Cox and C Maeve McCarthy, The shape of the tallest column, SIAM Journal of Applied Math Anals 29 (1998), no. 3, 547–554.
-  J. B. Keller, The shape of the strongest column, Arch. Rational Mech. Anal 5 (1960), 275–285.
-  J. B. Keller and F. I. Niordson, The tallest column, J. Math. Mech. 16 (1966), 433–446.
-  Philip G. Kirmser and Kuo-Kuang Hu, The shape of the ideal column reconsidered, The Mathematical Intelligencer 15 (1993), no. 3, 62–68.
-  C. Maeve McCarthy, The tallest column — optimality revisited, Journal of computational and applied mathematics (1999), no. 101, 27–37.
-  Charles A. Stuart, Buckling of a heavy tapered rod, Journal de Mathématiques Pures et Appliquées 80 (2001), no. 3, 281–337.