TY - JOUR
AU - Overton, Michael L
AB - Summary We consider the problem of optimal placement of concentrated masses along a massless elastic column that is clamped at one end and loaded by a nonconservative follower force at the free end. The goal is to find the largest possible interval such that the variation in the loading parameter within this interval preserves stability of the structure. The stability constraint is nonconvex and nonsmooth, making the optimization problem quite challenging. We give a detailed analytical treatment for the case of two masses, arguing that the optimal parameter configuration approaches the flutter and divergence boundaries of the stability region simultaneously. Furthermore, we conjecture that this property holds for any number of masses, which in turn suggests a simple formula for the maximal load interval for |$n$| masses. This conjecture is strongly supported by extensive computational results, obtained using the recently developed open-source software package granso (GRadient-based Algorithm for Non-Smooth Optimization) to maximize the load interval subject to an appropriate formulation of the nonsmooth stability constraint. We hope that our work will provide a foundation for new approaches to classical long-standing problems of stability optimization for nonconservative elastic systems arising in civil and mechanical engineering. 1. Introduction Consider an elastic Euler–Bernoulli beam clamped at one end and loaded at the tip by a follower force (1, 2). The follower force is defined as a force with the line of action that always coincides with the tangent line to the neutral axis of the deformed beam at its free end, much like a rocket thrust (3). The follower force does not depend on the velocity of the beam. However, it cannot be derived from a potential: the work done by the follower force along a closed contour is nonzero (4, 5). This structure is frequently called the Beck column (1, 2). A straight form of the Beck column is in a stable equilibrium when the follower force is absent or relatively small. Nevertheless, at some sufficiently large value the follower force excites exponentially growing oscillations of the beam that are known as flutter instability (6, 7). Flutter is critically important both for safety of engineering structures interacting with fluid flows and for efficiency of energy harvesting devices that are based on the fluid-structure interactions. Recent years have seen an increasing interest in the Beck column in the modelling of biological filaments and their artificial biomimetic analogues, that is, hair-like slender microscale structures that play an important part in such biological processes as swimming, pumping, mixing and cytoplasmic streaming by performing rhythmic, wave-like motion that usually sets in via flutter instability (8–12). Structural optimization of the Beck column against instabilities, including flutter and buckling (or divergence instability), is usually formulated as a problem on a redistribution of the material of the column of a given density under an isoperimetric constraint fixing the volume of the column in order to maximize the range of variation of the follower load corresponding to the stable structure. In the literature, many specific numerically optimized shapes of the Beck column have been reported (29–34) with the maximal critical dimensionless load reaching the values of |$p\approx 100.00$| (33), |$p\approx 139.30$| (35), |$p\approx 143.59$| (36) and |$p\approx 148.62$| (37), which significantly improve upon the critical load |$p\approx 20.05$| of the uniform column with a constant cross-section (see Appendix A for the definition of |$p$|⁠). Nevertheless, none of these designs is proven to be a global or even a local optimizer. Such a proof would be difficult to obtain because the problem of structural optimization of the critical flutter load for the elastic Beck column is both nonconvex and nonsmooth (38, 39). Indeed, the elastic Beck column is a time-reversible dynamical system in which the transition from stability to flutter instability generically happens via the reversible-Hopf bifurcation, that is, through the formation of a double imaginary eigenvalue with a Jordan block at the stability boundary and its subsequent splitting when parameters enter the instability region (40). Codimension-1 parts of the stability boundary are thus smooth hypersurfaces corresponding to double imaginary eigenvalues with a Jordan block (provided that the remaining eigenvalues are simple and imaginary) (7, 41, 42). These hypersurfaces can meet each other at sets of higher codimension such as intersections, cuspidal edges and points, conical points etc.; see (7) for a full classification of generic singularities on the stability boundary of mechanical systems with nonpotential positional forces. The unavoidable singularities linked to multiple eigenvalues is the main reason for nonsmoothness of the merit functionals in the optimization of such systems, including the Beck column, with respect to stability criteria (43, 44). Many studies report on the phenomenon of overlapping of eigenvalue curves that accompanies the process of optimization of the Beck column. The eigenfrequencies plotted as functions of the load exhibit sudden crossings during the optimization that lead to transfer of instability between modes and to a discontinuous change in the merit functional (17, 30–34, 36, 37, 45–47). The high sensitivity of the optimized design to variation of parameters is caused by the nonconvexity of the stability domain (38, 39). For this reason, the unambiguous determination of the optimal design of the Beck column by numerical procedures typically used in civil and mechanical engineering remains a challenge (33, 36, 37). All of the phenomena described above were also observed in simplified settings with the uniform Beck column carrying relocatable lumped masses (24, 47–53). Nevertheless, to the best of our knowledge, no rigorously proven local or global optimal solutions or credible numerical guesses exist in the literature even in the problems of optimal localization of point masses along elastic beams loaded by the follower force. Structures loaded by follower forces have long been questioned for their practical realization (13, 14), despite an evident example given by flexible missiles (15–17). In the 1970–1990s, Sugiyama et al. (3, 18, 19) used solid rocket motors to demonstrate flutter of cantilevers under a follower thrust on relatively short (several seconds) time intervals. A mechanism recently invented by Bigoni and Noselli produces a frictional follower force (20, 54) and enables experimental realization of fluttering cantilevered rods under follower loads on virtually infinite time intervals (21, 22). These practical realizations differ from the classical Beck column; however, by the presence of a finite-size loading unit at the tip of the cantilever and therefore are better described by the model of the Pflüger column (23, 24), which is the Beck column with a point mass at the loaded end; see the left panel of Fig. A.1 in Appendix A. In recent mechanical laboratory experiments with follower forces (21, 22), the ratio of the end mass to the mass of the column was chosen to be very large, approaching the so-called Dzhanelidze limit corresponding to a massless column (6). The instability thresholds obtained in these experiments were in a very good agreement with the theoretical predictions based on the Pflüger model. In the Dzhanelidze limit, the mathematical model is reduced to a system of ordinary differential equations (55–58). The works (6, 24, 47) considered stability of a massless Pflüger column with an additional relocatable mass. A recent work (58) corrected some of the results reported in (6) and proposed extending the model to incorporate several relocatable masses. The primary purpose of our article is to study this last variant, the Pflüger model in the Dzhanelidze massless limit with relocatable point masses, in detail. One reason is that this comparatively simple but still mechanically meaningful model allows a detailed analytical treatment of the case of two masses, providing a benchmark for numerical optimization carried out for |$n$| masses. A second advantage of studying the discrete mass model instead of the classical Beck column is that it does not require Galerkin or finite element discretization, and hence the number of optimization variables is small (only |$2n-1$|⁠). Nonetheless, the problem of maximizing the load interval subject to the stability constraint is far from trivial because of the nonconvexity and nonsmoothness (in fact, non-Lipschitzness) of the constraint, so even this simplified model provides a good test of how much insight we can obtain using nonsmooth optimization techniques. Our first contribution, presented in section 3, is to give a detailed analytical treatment for the case of two masses, arguing that the optimal parameter configuration approaches the flutter and divergence boundaries simultaneously. Furthermore, we conjecture that this property holds for any number of masses, which in turn suggests a simple formula for the optimal load interval for |$n$| masses. Our second contribution, in section 4, is to present a practical numerical formulation of the stability constraint and to maximize the load interval subject to this constraint using modern techniques for nonsmooth, nonconvex optimization, employing a recently developed open-source software package, granso (GRadient-based Algorithm for Non-Smooth Optimization) (59, 60). As well as verifying our analytical solution for two masses, these computations strongly support the formula for the conjectured optimal load interval for |$n$| masses. We hope that our techniques and results will provide a foundation and inspiration for new approaches to classical long-standing problems of stability optimization for nonconservative elastic systems arising in civil and mechanical engineering. 2. A massless elastic column with |$n$| concentrated masses It is convenient to first consider the simple model of the Pflüger column without relocatable masses, with zero mass per unit length and zero point mass at the free end of the column (see Appendix A for details). Then, the boundary value problem (A.5), (A.6) takes the form $$\begin{equation}\label{dlbz} \partial_{\xi}^4 f +\kappa^2 \partial_{\xi}^2 f=0, \end{equation}$$(2.1) $$\begin{equation}\label{dlbz1} f(0)=0,\quad \partial_{\xi}f(0)=0,\quad \partial_{\xi}^2 f(1)=0,\quad \partial_{\xi}^3 f(1)=0, \end{equation}$$(2.2) where $$\begin{equation}\label{dck} \kappa^2=p, \end{equation}$$(2.3) with |$p$| given in (A.4). Following (6, 24, 58), consider the case when a concentrated constant force |$F$| is acting in a direction perpendicular to the nondeformed column at the point |$s=\alpha l$|⁠. Introducing the dimensionless version of the force parameter, |$\phi=\frac{F l^2}{EI}$|⁠, we seek the general solution to (2.1) in the form (6, 24, 58) $$\begin{eqnarray}\label{gsa} f(\xi)=u(\xi)+\left\{\begin{array}{r} 0, \quad \xi\in[0,\alpha) \\ v(\xi), \quad \xi\in[\alpha,1] \end{array}\right., \end{eqnarray}$$(2.4) where $$u(\xi)=A\sin{\kappa\xi}+B\cos{\kappa\xi}+C\xi+D$$ and $$v(\xi)=A_1\sin{\kappa \xi}+B_1\cos{\kappa \xi}+C_1\xi+D_1.$$ To determine the coefficients |$A_1$|⁠, |$B_1$|⁠, |$C_1$| and |$D_1$|⁠, we require that $$\begin{eqnarray}\label{glu} &u(\alpha)=f(\alpha),\quad \partial_{\xi}u(\alpha)=\partial_{\xi}f(\alpha),&\nonumber\\ &\partial^2_{\xi}u(\alpha)=\partial^2_{\xi}f(\alpha),\quad \partial^3_{\xi}f(\alpha)-\partial^3_{\xi}u(\alpha)=\phi.& \end{eqnarray}$$(2.5) This yields $$\begin{equation}\label{fv} v(\xi) = \frac{(\xi-\alpha)\kappa-\sin((\xi-\alpha)\kappa)}{\kappa^3}\phi. \end{equation}$$(2.6) Taking (2.6) into account in the general solution (2.4) and then substituting |$f(\xi)$| into the boundary conditions (2.2), we find the coefficients |$A$|⁠, |$B$|⁠, |$C$| and |$D$| to obtain $$\begin{equation}\label{fu} u(\xi) = \frac{\sin(\kappa \alpha)-\xi\kappa\cos(\kappa \alpha)+\sin((\xi-\alpha)\kappa)}{\kappa^3}\phi. \end{equation}$$(2.7) Fig. 1 Open in new tabDownload slide The massless Beck column loaded by the follower force |$P$| with |$n$| concentrated masses |$M_1$|⁠, |$\ldots$|⁠, |$M_i$|⁠, |$\ldots$|⁠, |$M_n$| attached (53, 58). Fig. 1 Open in new tabDownload slide The massless Beck column loaded by the follower force |$P$| with |$n$| concentrated masses |$M_1$|⁠, |$\ldots$|⁠, |$M_i$|⁠, |$\ldots$|⁠, |$M_n$| attached (53, 58). Let us now assume that the massless cantilevered column loaded by the follower force at its free end carries |$n$| concentrated masses with the mass |$M_n>0$| fixed at the loaded end; see Fig. 1. The masses |$M_i\ge 0$|⁠, |$i=1,\ldots,n-1,$| are located at the distances |$s_i<l$| from the clamped end of the column. Let |$v_i$| be a transversal displacement of the mass |$M_i$| from the equilibrium configuration, as shown in Fig. 1. Introducing the dimensionless displacements of the masses, |$w_i$|⁠, the distances, |$\alpha_i$|⁠, and the mass ratios, |$\mu_i$|⁠, as $$\begin{equation}\label{dlp} w_i=\frac{v_i}{l}, \quad \alpha_i=\frac{s_i}{l},\quad \mu_i=\frac{M_i}{M_n}, \quad i=1,\ldots,n, \end{equation}$$(2.8) we write the equations of motion of the masses (6, 24, 58) $$\begin{eqnarray}\label{emm} w_1&=&-\gamma_{11} \mu_1 \frac{d^2 w_1}{d \tau^2}-\gamma_{12} \mu_2 \frac{d^2 w_2}{d \tau^2}-\ldots -\gamma_{1n} \mu_n \frac{d^2 w_n}{d \tau^2},\nonumber\\ w_2&=&-\gamma_{21} \mu_1 \frac{d^2 w_1}{d \tau^2}-\gamma_{22} \mu_2 \frac{d^2 w_2}{d \tau^2}-\ldots -\gamma_{2n} \mu_n \frac{d^2 w_n}{d \tau^2},\nonumber\\ &\vdots & \nonumber\\ w_i&=&-\gamma_{i1} \mu_1 \frac{d^2 w_1}{d \tau^2}-\gamma_{i2} \mu_2 \frac{d^2 w_2}{d \tau^2}-\ldots -\gamma_{in} \mu_n \frac{d^2 w_n}{d \tau^2},\nonumber\\ &\vdots & \nonumber\\ w_n&=&-\gamma_{n1} \mu_1 \frac{d^2 w_1}{d \tau^2}-\gamma_{n2} \mu_2 \frac{d^2 w_2}{d \tau^2}-\ldots -\gamma_{nn} \mu_n \frac{d^2 w_n}{d \tau^2}, \end{eqnarray}$$(2.9) where the dimensionless time |$\tau$| is defined now as $$\begin{equation}\label{dt2} \tau=t\sqrt{\frac{EI}{M_n l^3}}. \end{equation}$$(2.10) Note that |$\alpha_{1}\leqslant \alpha_{2}\leqslant \ldots \leqslant \alpha_{n}=1$| and |$\mu_{n}=1.$| The coefficient |$\gamma_{ij}$| is the displacement of the mass |$\mu_i$| as a result of application to the column of a unit force |$\phi=1$| at the point |$\alpha_j$|⁠. According to (2.4) with the functions (2.6) and (2.7), the coefficient |$\gamma_{ij}$| is given by |$\delta_{ij}/\kappa^3$|⁠, where $$\begin{eqnarray}\label{deltaij} \delta_{ij}&=&{\sin(\kappa \alpha_j)-\alpha_i\kappa\cos(\kappa \alpha_j)+\sin((\alpha_i-\alpha_j)\kappa)}\nonumber\\ & & + \left\{\begin{array}{r} 0, \quad i \leqslant j \\ (\alpha_i-\alpha_j)\kappa-\sin((\alpha_i-\alpha_j)\kappa), \quad i>j. \\ \end{array}\right. \end{eqnarray}$$(2.11) Separating time with the ansatz |$w_i=u_ie^{\sigma \kappa^{3/2}\tau}$| we arrive at the eigenvalue problem $$\begin{equation}\label{eip} ({{\bf M}}\sigma^2+{{\bf K}})u=0, \end{equation}$$(2.12) where |$u=(u_1,u_2,\ldots, u_n)$|⁠, |${\bf K}$| is the |$n\times n$| unit matrix, and $$\begin{equation}\label{mass} {{\bf M}}=\left( \begin{array}{llll} \mu_1 \delta_{11} & \mu_2\delta_{12} & \cdots & \mu_n\delta_{1n}\\ \mu_1 \delta_{21} & \mu_2\delta_{22} & \cdots & \mu_n\delta_{2n}\\ \vdots & \vdots & \ddots & \vdots \\ \mu_1 \delta_{n1} & \mu_2\delta_{n2} & \cdots & \mu_n\delta_{nn}\\ \end{array} \right)\!, \end{equation}$$(2.13) where, as already noted, |$\mu_{n}=1$|⁠. The eigenvalues |$\sigma_{k}$| are given by $$\begin{equation}\label{sigma-lambda} \sigma_{k} = \pm \sqrt{-\lambda_{k}^{{-1}}}, \end{equation}$$(2.14) where the |$\lambda_{k}$| are the eigenvalues of the matrix |${\bf M}$|⁠. The trivial equilibrium of the circulatory system (2.9) is stable if and only if the eigenvalues |$\sigma_{k}$| are imaginary and semisimple (that is, the algebraic and geometric multiplicity are equal), or equivalently, the |$\lambda_{k}$| are real, positive and semisimple. Cases with a multiple imaginary eigenvalue |$\sigma_{k}$| with a Jordan block (that is, with the algebraic multiplicity exceeding the geometric multiplicity) lie on the boundary between the stability and flutter domains. In the generic case the crossing of this stability boundary is accompanied by merging of two simple imaginary eigenvalues into a double imaginary eigenvalue with a Jordan block, indicating the onset of the reversible-Hopf bifurcation or flutter (6, 7, 41, 42). Nonoscillatory instability or divergence corresponds to one or more positive real eigenvalues |$\sigma_{k}$| and in this model it generically sets in when two conjugate simple imaginary eigenvalues meet at infinity, split and turn back towards the origin along the real axis in the complex plane (26, 27, 58). Summarizing, for a given number of masses |$n$|⁠, the eigenvalue problem (2.12) is defined by (2.11) and (2.13), which depend on the given load |$\kappa$| and the parameters |$\alpha_{i}$| and |$\mu_{i}$|⁠, |$i=1,\ldots,n-1$|⁠, defined in (2.8) (as |$\alpha_{n}=\mu_{n}=1$|⁠). It is convenient to use the parameterization $$\begin{equation}\label{tanbeta} \mu_{i}=\tan\beta_{i}, \quad \beta_{i}\in[0,\pi/2), \quad i=1,\ldots,n-1. \end{equation}$$(2.15) Given |$\alpha_{i},\beta_{i}$|⁠, |$i=1,\ldots,n-1$|⁠, let us define |$\kappa_{{\rm crit}}^{\alpha,\beta}$| as the largest value such that the eigenvalues |$\sigma_{k}$| (which depend on |$\alpha_{i}$|⁠, |$\beta_{i}$| and |$\kappa$|⁠) are imaginary for all |$\kappa\in[0,\kappa_{{\rm crit}}^{\alpha,\beta}]$|⁠. Our goal is to find the supremum of |$\kappa_{{\rm crit}}^{\alpha,\beta}$| over all parameters |$\alpha_{i} \in [0,1]$| and |$\beta_{i}\in[0,\pi/2)$|⁠, |$i=1,\ldots,{n-1}$|⁠. We begin with the case |$n=2$|⁠, where we propose an analytical solution. 3. Analytical derivation of the supremal load interval for the massless column carrying two concentrated masses When |$n=2$|⁠, the massless column carries a relocatable mass |$M_{1}$| between the clamped end and the free end of the rod with mass |$M_2$| fixed at the free end. There are two parameters, |$\alpha_{1}$| and |$\beta_{1}$|⁠. Expression (2.11) allows us to find the coefficients |$\delta_{ij}$| in the explicit form, cf. (6, 58), $$\begin{eqnarray}\label{entries2} \delta_{11}&=&\sin(\kappa \alpha_1)-\kappa \alpha_1 \cos(\kappa \alpha_1),\nonumber\\ \delta_{12}&=&\sin(\kappa)-\kappa \alpha_1\cos(\kappa)-\sin(\kappa(1-\alpha_1)),\nonumber\\ \delta_{21}&=&\sin(\kappa\alpha_1)-\kappa\cos(\kappa \alpha_1)+\kappa(1-\alpha_1),\nonumber\\ \delta_{22}&=&\sin(\kappa)-\kappa\cos(\kappa). \end{eqnarray}$$(3.1) As we will see, already in this simplest possible mechanical system, the subdivision of the parameter space into the domains of stability, flutter instability and divergence instability is highly nontrivial. However, we will be able to explore it completely and find an apparent supremum of the critical load parameter defining the longest stability interval |$[0,\kappa_{{\rm crit}}^{\alpha_{1},\beta_{1}}]$| in the space of parameters |$\alpha_{1}\in[0,1]$|⁠, |$\beta_{1}\in[0,\pi/2)$|⁠. In general, the stability map for a mechanical system with the characteristic polynomial |$p(\sigma)=\det({\bf M}\sigma^2+K)$| can be obtained with the use of the Gallina criterion (7, 61, 62) that is based on the investigation of the discriminant of the polynomial. For |$n=2$|⁠, |$p(\sigma)$| is a biquadratic function $$\begin{eqnarray}\label{charp} p(\sigma)&=& \sigma^4\tan\beta_1\left\{\kappa(\alpha_1-1)(\sin\kappa-\kappa\alpha_1\cos\kappa+\sin(\kappa\alpha_1-\kappa))\right.\nonumber\\ &&\qquad -\left.\sin(\kappa(\alpha_1-1))(\sin(\kappa\alpha_1)-\kappa\cos(\kappa\alpha_1)-\kappa(\alpha_1-1))\right\}\nonumber\\ &&+ ~ \sigma^2\left[\tan\beta_1\left(\sin(\kappa\alpha_1)-\kappa\alpha_1\cos(\kappa\alpha_1)\right)-\kappa\cos\kappa+\sin\kappa\right]+1. \end{eqnarray}$$(3.2) Notice that the coefficient at the leading power of |$\sigma$| in the polynomial (3.2) is nothing else but |$\det {\bf M}$|⁠; see (7). The system loses stability by divergence as soon as |$\det {\bf M}=0$|⁠, which yields the following equation determining the divergence boundary: $$\begin{equation}\label{divb} \frac{\sin\kappa-\kappa\alpha_1\cos\kappa+\sin(\kappa\alpha_1-\kappa)}{\sin(\kappa\alpha_1)-\kappa\cos(\kappa\alpha_1)-\kappa(\alpha_1-1)} =\frac{\sin(\kappa(\alpha_1-1))}{\kappa(\alpha_1-1)}. \end{equation}$$(3.3) Note that this equation is independent of |$\beta_{1}$|⁠. The right panel of Fig. 2 shows the divergence boundary (3.3) in the |$(\alpha_1,\beta_1,\kappa)$|-space. Fig. 2 Open in new tabDownload slide The case of |$n=2$| concentrated masses. (Left) the flutter domain is a finite solid set in the |$(\alpha_1,\beta_1,\kappa)$| space, enclosed within the singular surface defined by (3.4). (Right) The divergence domain lies above the boundary set defined by (3.3). For a given |$(\alpha_1$|⁠, |$\beta_1)$|⁠, the critical value of the load parameter, |$\kappa_{{\rm crit}}^{\alpha_{1},\beta_{1}}$|⁠, is the minimal value of |$\kappa$| that satisfies either (3.4) or (3.3), as this is the length of the longest vertical line segment rising from the point |$(\alpha_{1},\beta_{1},0)$| that does not enter either the flutter or divergence domain. Consequently, this is the largest value |$\tilde\kappa$| such that the column is stable for all |$\kappa \in [0,\tilde\kappa)$|⁠. The optimization problem to be solved is to find the supremum of |$\kappa_{{\rm crit}}^{\alpha_{1},\beta_{1}}$| over all |$\alpha_{1}\in[0,1], \beta_{1}\in [0,\pi/2)$|⁠. Fig. 2 Open in new tabDownload slide The case of |$n=2$| concentrated masses. (Left) the flutter domain is a finite solid set in the |$(\alpha_1,\beta_1,\kappa)$| space, enclosed within the singular surface defined by (3.4). (Right) The divergence domain lies above the boundary set defined by (3.3). For a given |$(\alpha_1$|⁠, |$\beta_1)$|⁠, the critical value of the load parameter, |$\kappa_{{\rm crit}}^{\alpha_{1},\beta_{1}}$|⁠, is the minimal value of |$\kappa$| that satisfies either (3.4) or (3.3), as this is the length of the longest vertical line segment rising from the point |$(\alpha_{1},\beta_{1},0)$| that does not enter either the flutter or divergence domain. Consequently, this is the largest value |$\tilde\kappa$| such that the column is stable for all |$\kappa \in [0,\tilde\kappa)$|⁠. The optimization problem to be solved is to find the supremum of |$\kappa_{{\rm crit}}^{\alpha_{1},\beta_{1}}$| over all |$\alpha_{1}\in[0,1], \beta_{1}\in [0,\pi/2)$|⁠. The roots of the characteristic polynomial (3.2) are double imaginary if the discriminant of the biquadratic function vanishes: $$\begin{align}\label{discra} &(\sin(\kappa\alpha_1)-\kappa\alpha_1\cos(\kappa\alpha_1))^2(\tan\beta_1)^2\nonumber\\ &\quad{} +~2\alpha_1\kappa^2\tan\beta_1\cos\kappa\left[\cos(\kappa\alpha_1)+2(\alpha_1-1)\right]\nonumber\\ &\quad{} +~2\tan\beta_1\sin(\kappa\alpha_1)\left[2\sin(\kappa(\alpha_1-1))+\sin\kappa\right]\nonumber\\ &\quad{} -~\kappa\tan\beta_1\left[7\sin(\kappa(\alpha_1-1))(\alpha_1-1) +\sin(\kappa(\alpha_1+1))(\alpha_1+1)\right]\nonumber\\ &\quad{} -~2\kappa\tan\beta_1\left[(2\alpha_1-3)\sin\kappa+\sin(\kappa(2\alpha_1-1))\right]\nonumber\\ &\quad{} +~(\sin\kappa-\kappa\cos\kappa)^2=0. \end{align}$$(3.4) For this reason (6, 7, 42), (3.4) determines the boundary of the flutter domain that is shown in the left panel of Fig. 2. Fig. 3 Open in new tabDownload slide Stability diagrams for (left) |$\beta_{1} = {\hat\beta =} \#1.450234089$| in the |$(\alpha_1,\kappa)$|-plane and (right) |$\alpha_1= {\hat\alpha =} 0.4947347666$| in the |$(\beta_1,\kappa)$|-plane. The solid thin (blue) curves at the top of each panel designate the divergence boundary (3.3) and the solid thick (green) curves mark the flutter boundary (3.4). The flutter boundary in the left panel has a crossing at the saddle point located at |$\alpha_{1}=\hat\alpha$| and |$\kappa = 5.591633160$|⁠. The dashed (black) curves in the left panel correspond to the flutter boundaries at (upper and lower curves) |$\beta_1=\hat\beta-0.01$| and (left and right curves) |$\beta_1=\hat\beta+0.01$|⁠. In the right panel, the divergence boundary is a horizontal line with height |$\kappa=\hat\kappa=7.113918994$|⁠. (Colors in parentheses are relevant only if the article is viewed online.) Fig. 3 Open in new tabDownload slide Stability diagrams for (left) |$\beta_{1} = {\hat\beta =} \#1.450234089$| in the |$(\alpha_1,\kappa)$|-plane and (right) |$\alpha_1= {\hat\alpha =} 0.4947347666$| in the |$(\beta_1,\kappa)$|-plane. The solid thin (blue) curves at the top of each panel designate the divergence boundary (3.3) and the solid thick (green) curves mark the flutter boundary (3.4). The flutter boundary in the left panel has a crossing at the saddle point located at |$\alpha_{1}=\hat\alpha$| and |$\kappa = 5.591633160$|⁠. The dashed (black) curves in the left panel correspond to the flutter boundaries at (upper and lower curves) |$\beta_1=\hat\beta-0.01$| and (left and right curves) |$\beta_1=\hat\beta+0.01$|⁠. In the right panel, the divergence boundary is a horizontal line with height |$\kappa=\hat\kappa=7.113918994$|⁠. (Colors in parentheses are relevant only if the article is viewed online.) For a given (⁠|$\alpha_1$|⁠, |$\beta_1$|⁠), the critical value of the load parameter is given by $$\kappa_{{\rm crit}}^{\alpha_{1},\beta_{1}} = \min \{\kappa: (\kappa,\alpha_{1},\beta_{1}) \text{ satisfies either }(3.3) \text{ or }(3.4) \},$$ as this is the length of the longest vertical line segment rising from the point |$(\alpha_{1},\beta_{1},0)$| that does not enter either the flutter or the divergence domain. Consequently, the quantity $$\begin{equation}\label{kappacritdef} \kappa^{*} = \sup\{\kappa_{{\rm crit}}^{\alpha_{1},\beta_{1}}: \alpha_{1}\in[0,1],\beta_{1}\in[0,\pi/2) \} \end{equation}$$(3.5) is the supremum of all loads associated with a stable column. Note that although the divergence boundary (3.3) is smooth, the boundary of the flutter domain (3.4) is nonsmooth. Figure 3 shows cross-sections of the flutter boundary and the divergence boundary in the |$(\alpha_1,\kappa)$|- and |$(\beta_1,\kappa)$|-planes. In the left panel, for which |$\beta_1$| is fixed to |$\hat\beta \approx 1.45$|⁠, we see that the flutter boundary has a saddle point in the |$(\alpha_1,\kappa)$|-plane at |$\alpha_1=\hat\alpha\approx0.495$|⁠, |$\kappa\approx 5.59$|⁠. On the other hand, when |$\alpha_{1}$| is fixed to |$\hat\alpha$|⁠, the flutter boundary has a vertical tangent in the |$(\beta_1,\kappa)$|-plane at |$\beta_1=\hat\beta$|⁠, as is visible in the right panel of Fig. 3. Consequently, when |$\alpha_{1}=\hat\alpha$|⁠, the maximal stable load |$\kappa_{{\rm crit}}^{\alpha_{1},\beta_{1}}$| varies smoothly for |$\beta_{1}\in(0,\hat\beta)$|⁠, but when |$\beta_{1}$| reaches |$\hat\beta$| it jumps up discontinuously from the flutter boundary to the divergence boundary. For the system under study, such jumps were first described in the work (58) that corrected the classical result of Bolotin (6), whose plot in the |$(\beta_1,\kappa)$|-plane did not contain the divergence boundary at all, but provided the correct shape for the flutter boundary. Notice that such overlapping of eigenvalue branches typically accompanies optimization of nonconservative systems and was reported in numerous studies (6, 17, 19, 27, 30–37). The general theory of this effect has been developed in (7, 38, 39). Fig. 4 Open in new tabDownload slide (Left) imaginary and (right) real roots of the characteristic polynomial (3.2) for |$\alpha_1=\hat\alpha=0.4947347666$|⁠, and |$\beta_1=\hat\beta=1.450234089$| (solid (green) curves) and |$\hat\beta\pm0.01$| (dashed (black) curves). A bubble of complex eigenvalues appears for |$\beta_1=1.450234089-0.01$| and corresponds to flutter instability. The dotted (black) vertical line at |$\kappa=\hat\kappa=7.113918994$| marks the onset of divergence instability. Increasing |$\beta_1$| from |$\hat\beta-0.01$| to |$\hat\beta+0.01$| results in the disappearance of the complex eigenvalues and hence is accompanied by the transition from the overlapping eigenvalue branches to an avoided crossing that yields a jump in the critical load parameter to the maximal value that is reached at |$\kappa=\hat\kappa$| on the divergence boundary (58); see also the right panel of Fig. 3. Fig. 4 Open in new tabDownload slide (Left) imaginary and (right) real roots of the characteristic polynomial (3.2) for |$\alpha_1=\hat\alpha=0.4947347666$|⁠, and |$\beta_1=\hat\beta=1.450234089$| (solid (green) curves) and |$\hat\beta\pm0.01$| (dashed (black) curves). A bubble of complex eigenvalues appears for |$\beta_1=1.450234089-0.01$| and corresponds to flutter instability. The dotted (black) vertical line at |$\kappa=\hat\kappa=7.113918994$| marks the onset of divergence instability. Increasing |$\beta_1$| from |$\hat\beta-0.01$| to |$\hat\beta+0.01$| results in the disappearance of the complex eigenvalues and hence is accompanied by the transition from the overlapping eigenvalue branches to an avoided crossing that yields a jump in the critical load parameter to the maximal value that is reached at |$\kappa=\hat\kappa$| on the divergence boundary (58); see also the right panel of Fig. 3. We can obtain a clearer picture of the jump discontinuity by plotting the real and imaginary parts of the eigenvalues |$\sigma$| which describe the flutter boundary, as is done in Fig. 4. For |$\alpha_{1}=\hat\alpha$|⁠, when |$\beta_{1}$| is decreased from the value |$\hat\beta$|⁠, a bubble of complex eigenvalues corresponding to flutter appears, but this vanishes for |$\beta_{1}\geqslant\hat\beta$|⁠, resulting in the transition of the critical load from the flutter boundary to the divergence boundary. Fig. 5 Open in new tabDownload slide Stability diagrams for (left) |$\beta_1=\tilde\beta=0.4342999969$| in the |$(\alpha_1,\kappa)$|-plane and (right) for |$\alpha_1=\tilde\alpha=0.5810701268$| in the |$(\beta_1,\kappa)$|-plane. The solid thin (blue) curves at the top of each panel designate the divergence boundary (3.3), and the solid thick (green) curves mark the flutter boundary (3.4). The flutter boundary in the left panel has a crossing at the saddle point located at |$\alpha_{1}=\tilde\alpha$| and |$\kappa = 6.600674669$|⁠. The dashed (black) curves in the left panel correspond to the flutter boundaries at (upper and lower curves) |$\beta_{1}=\tilde\beta-0.05$| and (left and right curves) |$\beta_{1}=\tilde\beta+0.05$|⁠. In the right panel, the divergence boundary is the horizontal line with height |$\kappa=7.607584259$|⁠. The horizontal dash-dotted (red) line in the left panel shows the value |$\kappa_{0}$| given in (3.8) which is the smallest positive root of (3.6): the flutter boundary for the case |$\beta_1=0$|⁠. The other dash-dotted (red) curve in the left panel is the solution to (3.7): the flutter boundary for the case |$\beta_1=\pi/2$|⁠. Fig. 5 Open in new tabDownload slide Stability diagrams for (left) |$\beta_1=\tilde\beta=0.4342999969$| in the |$(\alpha_1,\kappa)$|-plane and (right) for |$\alpha_1=\tilde\alpha=0.5810701268$| in the |$(\beta_1,\kappa)$|-plane. The solid thin (blue) curves at the top of each panel designate the divergence boundary (3.3), and the solid thick (green) curves mark the flutter boundary (3.4). The flutter boundary in the left panel has a crossing at the saddle point located at |$\alpha_{1}=\tilde\alpha$| and |$\kappa = 6.600674669$|⁠. The dashed (black) curves in the left panel correspond to the flutter boundaries at (upper and lower curves) |$\beta_{1}=\tilde\beta-0.05$| and (left and right curves) |$\beta_{1}=\tilde\beta+0.05$|⁠. In the right panel, the divergence boundary is the horizontal line with height |$\kappa=7.607584259$|⁠. The horizontal dash-dotted (red) line in the left panel shows the value |$\kappa_{0}$| given in (3.8) which is the smallest positive root of (3.6): the flutter boundary for the case |$\beta_1=0$|⁠. The other dash-dotted (red) curve in the left panel is the solution to (3.7): the flutter boundary for the case |$\beta_1=\pi/2$|⁠. Looking at the discriminant (3.4), we notice that it degenerates into the equation $$\begin{equation}\label{nomass} \kappa\cos(\kappa)-\sin(\kappa)=0 \end{equation}$$(3.6) for |$\beta_1=0$| (that is, when |$\mu_1=0$|⁠) and reduces to the equation $$\begin{equation}\label{allmass} \sin(\kappa\alpha_1)-\kappa\alpha_1\cos(\kappa\alpha_1)=0 \end{equation}$$(3.7) in the limit |$\beta_{1}\to\pi/2$| (that is, |$\mu_{1}\to\infty$|⁠). The sets defined by (3.6) and (3.7) are shown by the dash-dotted (red) line and curve, respectively, in the left panel of Fig. 5. The flutter boundary is tangent to the planes |$\beta_1=0$| and |$\beta_1=\pi/2$| along this line and curve. Note that the height of the line is the smallest positive root of (3.6), which we denote by |$\kappa_{0}$|⁠, with $$\begin{equation}\label{kap0def} \kappa_{0}\approx 4.493409458. \end{equation}$$(3.8) Since |$\beta_{1}=0$| is the case where the mass |$M_{1}=0$|⁠, |$\kappa_{0}$| is the square root of the critical load for the Dzhanelidze column (in view of (A.9) and (2.3)). The lines |$\kappa=\kappa_0$| at |$\alpha_1=0$| and |$\alpha_1=1$| form singularities (edges) of the flutter domain. Fig. 6 Open in new tabDownload slide Stability diagrams in the |$(\alpha_1,\kappa)$|-plane for (upper left) |$\beta_1=\pi/2-0.5$|⁠, (upper right) |$\beta_1=\pi/2-0.15$|⁠, (middle left) |$\beta_1=\pi/2-0.1$|⁠, (middle right) |$\beta_1=\pi/2-0.05$|⁠, (lower left) |$\beta_1=\pi/2-0.01$| and (lower right) |$\beta_1=\pi/2-0.001$|⁠. The regions of flutter instability shrink to the solution of (3.7), shown by the dash-dotted (red) curve, as |$\beta_{1} \to \pi/2.$| The dotted (black) lines intersect at the point with the coordinates of the optimal solution: |$\alpha_{1}^{*} \approx 0.588527598$| and |$\kappa^{*}\approx 7.635002111$|⁠. Fig. 6 Open in new tabDownload slide Stability diagrams in the |$(\alpha_1,\kappa)$|-plane for (upper left) |$\beta_1=\pi/2-0.5$|⁠, (upper right) |$\beta_1=\pi/2-0.15$|⁠, (middle left) |$\beta_1=\pi/2-0.1$|⁠, (middle right) |$\beta_1=\pi/2-0.05$|⁠, (lower left) |$\beta_1=\pi/2-0.01$| and (lower right) |$\beta_1=\pi/2-0.001$|⁠. The regions of flutter instability shrink to the solution of (3.7), shown by the dash-dotted (red) curve, as |$\beta_{1} \to \pi/2.$| The dotted (black) lines intersect at the point with the coordinates of the optimal solution: |$\alpha_{1}^{*} \approx 0.588527598$| and |$\kappa^{*}\approx 7.635002111$|⁠. As soon as |$\beta_1$| starts deviating from zero, a closed region of flutter instability appears around the horizontal line |$\kappa=\kappa_0$| in the |$(\alpha_1,\kappa)$|-plane. Furthermore, another region of flutter originates above it that touches the divergence boundary. These two regions coalesce when |$\beta_{1}$| reaches |$\tilde\beta\approx 0.434$|⁠; see the left panel of Fig. 5, which shows another resulting saddle point on the flutter boundary defined by (3.4). With further growth in |$\beta_1$| the flutter region in the |$(\alpha_1,\kappa)$|-plane is simply connected, as shown in the two upper panels of Fig. 6 corresponding to |$\beta_1=\pi/2-0.5$| and |$\beta_1=\pi/2-0.15$|⁠, respectively, until this parameter passes the value |$\beta_1 \approx 1.45$|⁠, after which the flutter domain bifurcates into two parts; see the middle and the lower panels in Fig. 6. As |$\beta_1$| approaches |$\pi/2$|⁠, the upper portion of the flutter region concentrates around the solution of (3.7), shown by the dash-dotted (red) curve in Fig. 6, and coincides with this curve exactly at |$\beta_{1}=\pi/2$|⁠. At this very limit, the critical load |$\kappa$| reaches its supremal value |$\kappa^{*}$|⁠, defined in (3.5), which can be obtained by finding the intersection point of the solution of (3.7), shown by the dash-dotted (red) curve, and the divergence boundary (3.3), shown by the solid thin (blue) curve. Solving (3.3) and (3.7) simultaneously, we find $$\begin{equation}\label{supremum-kappa-alpha} \kappa^{*}\approx 7.635002112,\quad \alpha_{1}^{*} \approx 0.5885275986, \end{equation}$$(3.9) and we write $$\begin{equation}\label{supremum-beta} \beta_{1}^{*}=\frac{\pi}{2} \end{equation}$$(3.10) to indicate that the supremum occurs in the limit |$\beta_{1}\to\pi/2$|⁠. Fig. 7 Open in new tabDownload slide Stability diagrams for (left) |$\alpha_1=\alpha_{1}^{*}-0.1$|⁠, (center) |$\alpha_1=\alpha_{1}^{*}\approx 0.5885275986$| and (right) |$\alpha_1=\alpha_{1}^{*}+0.1$|⁠. The thick (green) and thin (blue) curves respectively show the flutter and divergence boundaries. In the left and center panels, the critical load reaches the divergence boundary, but this is higher in the center panel, and there it is reached only if |$\beta_{1}=\pi/2$|⁠. In the right panel, the flutter boundary prevents the critical load from reaching the divergence boundary. The dotted (black) line shows the value |$\kappa=\kappa^{*}.$| Fig. 7 Open in new tabDownload slide Stability diagrams for (left) |$\alpha_1=\alpha_{1}^{*}-0.1$|⁠, (center) |$\alpha_1=\alpha_{1}^{*}\approx 0.5885275986$| and (right) |$\alpha_1=\alpha_{1}^{*}+0.1$|⁠. The thick (green) and thin (blue) curves respectively show the flutter and divergence boundaries. In the left and center panels, the critical load reaches the divergence boundary, but this is higher in the center panel, and there it is reached only if |$\beta_{1}=\pi/2$|⁠. In the right panel, the flutter boundary prevents the critical load from reaching the divergence boundary. The dotted (black) line shows the value |$\kappa=\kappa^{*}.$| Stability diagrams in Fig. 7 presented in the |$(\beta_1,\kappa)$|-plane show the decrease in the critical load |$\kappa$| when |$\alpha_1$| deviates from the value |$\alpha_{1}^{*}$|⁠, indicating that the value |$\kappa^{*}$| is a local supremum in the parameter space |$\alpha_{1} \in [0,1]$|⁠, |$\beta_{1} \in [0, \pi/2)$|⁠. Experiments reported in the next section strongly indicate that |$\kappa^{*}$| is actually the global supremum. However, note that |${\bf M}$| is not defined at |$\beta_{1}^{*}=\pi/2$|⁠, since then the mass ratio |$\mu_{1}=M_{1}/M_{2}$| is infinite, so the supremum is not attained. Furthermore, as |$(\kappa,\alpha_{1},\beta_{1})\rightarrow (\kappa^{*},\alpha_{1}^{*},\pi/2)$|⁠, the matrix element |${\bf M}_{21}$| diverges to |$\infty$| and |${\bf M}_{12}$| converges to |$0$| (see (3.3)), but |${\bf M}_{11}$| is the product of two quantities, one diverging to |$\infty$| and the other converging to |$0$| (see (3.7)). For this reason, it is difficult to rigorously state limiting properties of the eigenvalues |$\lambda_{k}$| of |${\bf M}$| as the supremum is approached, though based on both our symbolic and numerical calculations, it seems that, under the appropriate assumptions, the eigenvalues converge to a double zero eigenvalue with a Jordan block, indicating that the parameters are on the boundary of both the flutter and divergence domains, and that the limiting eigenvalues |$\sigma_{k}$| of (2.12) coalesce into a quadruple eigenvalue at |$\infty$|⁠. A key point in the derivation above is that the supremal value of |$\kappa_{{\rm crit}}^{\alpha,\beta}$| occurs when the divergence boundary meets the flutter boundary in the limit |$\beta_{1}\to\pi/2$|⁠. We conjecture that this property holds for all |$n\geqslant 2$|⁠, not just for |$n=2$|⁠. If we substitute (3.7), which is the equation for the flutter boundary in the limit |$\beta_{1}\to\pi/2$|⁠, into the divergence boundary equation (3.3), the latter can be simplified and reduced to $$\begin{equation}\label{simp} \kappa(\alpha_1-1)\sin(\kappa(\alpha_1-1))(\cos(\kappa\alpha_1)-1)^2=0. \end{equation}$$(3.11) Writing |$\sin(\kappa(\alpha_1-1))=0$| yields |$\kappa\alpha_1-\kappa + k \pi =0$|⁠, |$k\in \mathbb{Z}$|⁠. On the other hand, the relation (3.7) can be written as |$\kappa\alpha_1=\tan(\kappa\alpha_1)$|⁠, yielding |$\kappa\alpha_1=\kappa_0$|⁠, where |$\kappa_{0}$|⁠, given by (3.8), is the smallest positive root of the equation |$\tan \kappa=\kappa$|⁠. Combining the results, we obtain |$\kappa=\kappa_{0} + \pi k,$| with |$k\in \mathbb{Z}$|⁠. For |$k=1$|⁠, we obtain |$\kappa=\kappa_0+\pi$|⁠, which agrees with the supremum in (3.9) just obtained for the optimal load for two concentrated masses |$M_{1}$| and |$M_{2}$|⁠. Let us therefore set |$k=n-1$|⁠, giving $$\begin{equation}\label{kapm} \kappa=\kappa_{0} + (n-1)\pi, \end{equation}$$(3.12) and hence, using |$\kappa\alpha_1=\kappa_0$|⁠, $$\begin{equation}\label{kapm1} \alpha_1=\frac{\kappa_0}{\kappa_0+(n-1)\pi}. \end{equation}$$(3.13) For |$n=2$|⁠, we have |$\alpha_1=\kappa_0(\kappa_0+\pi)^{-1}$|⁠, which again agrees with the optimal value |$\alpha_{1}^{*}$| given in (3.9). This suggests a conjecture that (3.12) and (3.13) are respectively the supremal value |$\kappa^{*}$| and the corresponding limiting value |$\alpha_{1}^{*}$| for all |$n \geqslant 2$|⁠, with the corresponding limiting value |$\beta_{1}^{*}$| equal to |$\pi/2$|⁠. Figure 8 shows the values (3.12) and (3.13) as defined by the intersections of equations (3.3) and (3.7), the divergence boundary equation and the flutter boundary equation in the limit |$\beta_{1}=\pi/2$|⁠, respectively. (It is perhaps worth noting that, for all |$n$|⁠, we have |$\tan(\kappa_{0}+(n-1)\pi)=\tan(\kappa_{0})=\kappa_{0}$|⁠.) Remarkably, the numerical computations reported in the next section for |$n$| concentrated masses, with |$n=2,3,4,5$|⁠, strongly indicate that the supremal load |$\kappa^{*}$| and the corresponding limiting value |$\alpha_{1}^{*}$| are precisely the values given in (3.12) and (3.13) and illustrated in Fig. 8, with the corresponding limiting value |$\beta_{1}^{*}$| equal to |$\pi/2$|⁠. While we do not have conjectured formulas for the limiting values |$\alpha_{i}^{*}$| for |$i>1$| and |$n>2$|⁠, we conjecture that the corresponding limiting values |$\beta_{i}^{*}$| are all |$\pi/2$|⁠. Indeed, the property |$\beta_{1}^{*}=\pi/2$| implies that the mass ratio |$\mu_{1}=M_{1}/M_{n}\to\infty$| as |$\kappa\to\kappa^{*}$|⁠, which implies, assuming that |$M_{1}$| is bounded above, that |$M_{n}\to 0$|⁠. Consequently, if the other masses are nonzero in the limit, all mass ratios |$\mu_{i}=M_{i}/M_{n}$| must diverge to infinity as |$\kappa\to\kappa^{*}$|⁠. Fig. 8 Open in new tabDownload slide Graphs of the solutions of (3.7) defining the flutter boundary in the limit |$\beta_{1}\to\pi/2$| as a function of |$\alpha_{1}$|⁠, shown by the dash-dotted (red) curve, and (3.3) defining the divergence boundary as a function of |$\alpha_{1}$|⁠, shown by the solid (blue) curves. The intersection points are given by the expressions (3.12) and (3.13). Fig. 8 Open in new tabDownload slide Graphs of the solutions of (3.7) defining the flutter boundary in the limit |$\beta_{1}\to\pi/2$| as a function of |$\alpha_{1}$|⁠, shown by the dash-dotted (red) curve, and (3.3) defining the divergence boundary as a function of |$\alpha_{1}$|⁠, shown by the solid (blue) curves. The intersection points are given by the expressions (3.12) and (3.13). 4. Numerical derivation of the optimal load for the massless column carrying multiple relocatable masses Recall that, as discussed in section 2, for a given number of masses |$n$|⁠, our stability constraint is defined by the eigenvalue problem |$({\bf M}\sigma^{2} + {\bf K})u=0$| (see (2.12)). Here, |${\bf K}$| is the unit matrix while |${\bf M}$| is defined by (2.11) and (2.13), which depend on the dimensionless parameters |$\alpha_{i}$| and |$\mu_{i} = \tan\beta_{i}$|⁠, |$i=1,\ldots,n-1$|⁠, defined in (2.8), as well as a given load |$\kappa$|⁠. Let us write |${\bf M}(\alpha,\beta,\kappa)$| for the matrix |${\bf M}$| defined by |$\alpha=[\alpha_{1},\ldots,\alpha_{n-1}]^{T}$|⁠, |$\beta=[\beta_{1},\ldots, \beta_{n-1}]^{T}$| and |$\kappa$|⁠. As noted in (2.14), the eigenvalues |$\sigma_{k}$| of |$({\bf M}(\alpha,\beta,\kappa)\sigma^{2} + {\bf K})u=0$| are related to |$\lambda_{k}$|⁠, the eigenvalues of the matrix |${\bf M}(\alpha,\beta,\kappa)$|⁠, by |$\sigma_{k} = \pm (-\lambda_{k}^{{-1}})^{{1/2}}$|⁠. The stability constraint requires that, for given |$(\alpha,\beta,\kappa)$|⁠, all eigenvalues |$\sigma_{k}$| should be imaginary, or equivalently, that all eigenvalue reciprocals |$\lambda_{k}^{-1}$| are real and nonnegative. Clearly, another equivalent condition is that all eigenvalues |$\lambda_{k}$| are real and nonnegative, interpreting |$1/0$| as |$+\infty$|⁠. Consequently, we define a stability violation function |$\tilde v: \mathbb R^{2n-1}\to\mathbb R_{+}$| by $$\begin{equation}\label{stabviolinit} \big (\alpha, \beta,\kappa \big ) \mapsto \max \left ( {\rm Re} \sqrt{-\lambda_{k} } \right ), \end{equation}$$(4.1) where the maximum is taken over all eigenvalues of |${\bf M}(\alpha,\beta,\kappa)$|⁠, using the principal square root, hence implying that |$\tilde v$| cannot take negative values. Besides avoiding the nonlinearity in the reciprocal, the stability violation function |$\tilde v$| has the virtue that it is continuous, though not Lipschitz continuous, at points in parameter space where a positive eigenvalue |$\lambda_{k}$| passes through the origin to the negative real axis, and hence |$\tilde v$| changes continuously from the value zero to a positive value that grows like the square root function at zero. In this case, the parameters cross the divergence boundary, since a conjugate pair of imaginary eigenvalues |$\sigma_{k}$| coalesce at |$\infty$| and split along the real axis. The function |$\tilde v$| is also continuous, though not Lipschitz continuous, at points in parameter space where two positive real eigenvalues |$\lambda_{k},\lambda_{\ell}$| coalesce and split into a complex conjugate pair, and hence again |$\tilde v$| increases from zero to a positive quantity that, generically, increases with the square root of the perturbation. In this case, the parameters cross the flutter boundary, because two simple imaginary eigenvalues |$\sigma_{k},\sigma_{\ell}$| (and also their conjugates) coalesce on the imaginary axis and split into a complex pair. We argued in section 3 that, in the case |$n=2$|⁠, the optimal parameter configuration is simultaneously at both the flutter boundary and the divergence boundary, likely with a double eigenvalue |$\lambda$| at zero (equivalently, a quadruple eigenvalue |$\sigma$| at |$\infty$|⁠) and, if this is the case, generically, the stability violation function |$\tilde v$| would grow at nearby parameter configurations with the fourth root of the perturbation. To compensate for this non-Lipschitz behavior of |$\tilde v$|⁠, we define a modified stability violation function |$v: \mathbb R^{2n-1}\to\mathbb R$| by $$\begin{equation}\label{stabviol} v(\alpha,\beta,\kappa) = \left \{ \begin{array}{r} \tilde v(\alpha,\beta,\kappa)^{\rho}, \quad ~~~ \tilde v(\alpha,\beta,\kappa) \in [0,1]\\ \rho \tilde v(\alpha,\beta,\kappa)- (\rho - 1), \quad ~ \tilde v(\alpha,\beta,\kappa) \in [1,\infty] \end{array} \right., \end{equation}$$(4.2) where |$\rho$| is a positive integer. In the situations just discussed, the choice |$\rho=2$| is sufficient to make |$v$| generically Lipschitz continuous at points where the parameters cross either the divergence or the flutter boundary separately, and |$\rho=4$| is sufficient to make |$v$| Lipschitz continuous even when the parameters cross the divergence and flutter boundaries simultaneously, at least at the proposed optimal configuration (3.9), (3.10) for |$n=2$|⁠. In our computations, we experimented with choices of |$\rho$| from 1 to 5 and we found that |$\rho=4$| gave significantly better results than |$\rho<4$|⁠, but that setting |$\rho=5$| made no further improvement. Consequently, we chose to use |$\rho=4$|⁠. Note that the specific form of |$v$| is chosen so that it does not cause blow-up when |$\tilde v(\alpha,\beta,\kappa)$| is large, and so that it is continuously differentiable where |$\tilde v(\alpha,\beta,\kappa)=1$|⁠. However, what makes this problem particularly difficult is that as any |$\beta_{i}\rightarrow\pi/2$|⁠, the coefficient |$\mu_{i}\rightarrow \infty $| in (2.12). Consider the case |$n=2$|⁠. We already mentioned in section 3 that as |$\alpha_{1}\to\alpha_{1}^{*}$|⁠, |$\beta_{1}\to \pi/2$| and |$\kappa\to\kappa^{*}$|⁠, we have |${\bf M}_{21}\to\infty$| and |${\bf M}_{12}\to 0$|⁠, while |${\bf M}_{11}$| is a product of |$\tan(\beta_{1})$| with a second factor that converges to zero. If this second factor converges to zero more slowly than |$(\tan(\beta_{1}))^{-1}$| does, so that |$|{\bf M}_{11}|\to \infty$|⁠, a change in its sign causes an eigenvalue |$\lambda_{k}$| to discontinuously pass through |$\infty$| from the positive real to the negative real axis, implying infinitely large growth in the stability violation |$v$| as the parameters cross the divergence boundary. This presents a serious difficulty as we shall see. In order to solve our optimization problem, we need to impose the stability constraint not only at a given point |$(\alpha,\beta,\kappa)$| but also at all points |$(\alpha,\beta,\nu)$| with |$\nu\in [0,\kappa]$|⁠. Although we could construct an approximation to |$v(\alpha,\beta,\cdot)$| on the interval |$[0,\kappa]$| using approximation software such as Chebfun (63), this is computationally expensive. In our optimization computations, we found that a more effective approach is to impose the stability constraint on a coarse grid of |$\tilde q$| logarithmically spaced points on |$(0,\kappa]$|⁠, defining $$\begin{equation}\label{coarsegrid} c(\alpha,\beta,\kappa) = \max_{0\leqslant j\leqslant \tilde q} ( v(\alpha,\beta,\nu_{j}): \nu_{0}=\kappa, \nu_{j} = (1 - 2^{-j})\kappa, j=1,\ldots,\tilde q ) \end{equation}$$(4.3) and imposing the constraint |$c(\alpha,\beta,\kappa) \leqslant 0$|⁠, or equivalently, |$c(\alpha,\beta,\kappa) = 0$|⁠. Then, after a potential solution is obtained by optimization, we check its stability on a much finer grid of |$q\gg \tilde q$| uniformly spaced points on |$(0,\kappa)$|⁠, rejecting it if this test is not passed. We found that using a coarse grid with |$\tilde q=10$| points and a fine grid with |$q=10,000$| points worked well, typically with the majority of the solutions obtained by optimization that are feasible for the coarse grid also passing the fine grid test. We then pose our optimization problem as $$\begin{eqnarray}\label{optprob} \sup_{\alpha\in\mathbb R^{n-1},\beta\in\mathbb R^{n-1},\kappa\in\mathbb R} & \kappa \\ {\rm subject~to~} & c(\alpha,\beta,\kappa) \leqslant 0, \nonumber \\ &0 \leqslant \alpha_{1} \leqslant \ldots \leqslant \alpha_{n-1} \leqslant 1, \nonumber\\ &0 \leqslant \beta_{i} \leqslant \pi/2,~ i=1,\ldots,n-1, \nonumber \end{eqnarray}$$(4.4) This is not an easy problem to solve, since the stability constraint is nonconvex and nonsmooth, as well as discontinuous as |$\beta_{i}\to\pi/2$|⁠. We tackled it using granso, a recently developed open-source software package for nonsmooth constrained optimization (59, 60). As its name suggests, the algorithm implemented in granso is based on employing user-supplied gradients. This might seem contradictory since it is intended for nonsmooth optimization problems, but although the constraints are not differentiable everywhere, they are differentiable almost everywhere. Specifically, the stability violation function |$v$| is differentiable at |$(\alpha,\beta,\kappa)$| if the following conditions hold: the maximum in (4.3) is attained only at one index |$j \in (0,\ldots,\tilde q)$| the maximum in (4.1) is attained only at one eigenvalue |$\lambda_{k}$| of |${\bf M}(\alpha,\beta,\nu_j)$| this eigenvalue |$\lambda_{k}$| is simple and nonzero. Thus, evaluating the gradient of |$v$| makes sense almost everywhere in parameter space. Of course, the gradient does not vary continuously, but granso is designed to exploit gradient difference information, even near points where the gradient varies discontinuously, building a model of the constraint function on the parameter space using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) quasi-Newton updating method. For more details, see (59), and for application of BFGS in other stability optimization problems, see (64) and the papers cited there. To derive the gradient of |$v$|⁠, we need to differentiate an eigenvalue |$\lambda_{k}$| with respect to changes in the matrix |${\bf M}$|⁠. Let us write |${\bf M}(t)={\bf M} + t({\bf \Delta}{\bf M})$| and let |$\lambda(t)$| denote the eigenvalues of |${\bf M}(t)$|⁠. It is well known (65) that, if |$\lambda_{k} = \lambda(0)$| is a simple eigenvalue of |${\bf M}={\bf M}(0)$| satisfying the right and left eigenvector equations |${\bf M} u=\lambda u$| and |$w^{*} {\bf M}^{*} = \lambda w^{*}$|⁠, where the asterisk denotes complex conjugate transpose, then $$\left. \frac{d}{dt} \lambda(t)\right |_{t=0} = \frac{w^{*}( {\bf \Delta}{\bf M} )u} {w^{*} u}.$$ With this in mind, deriving the gradient of |$v$| with respect to the |$2n-1$| parameters given by |$(\alpha,\beta,\kappa)$| is straightforward, employing the chain rule to incorporate the variation in the power function in (4.2), the square root in (4.1), and the formulas (2.13), (2.11) and (2.15). We now describe our experiments using granso (version 1.6.4), running in matlab (release R2020a) on a MacBook Air laptop, to solve (4.4). We used the default choice of parameters with the following exceptions: we set maxit, the limit on the iteration count, to 500, and we set the tolerances opt_tol and feas_tol to zero, to obtain the highest possible accuracy. We added bound constraints on the load variable formulated as |$0 \leqslant \kappa \leqslant \kappa^{\max}$| with |$\kappa^{\max}=1.1\times (\kappa_0+(n-1)\pi)$|⁠, that is, with a lower bound of zero and an upper bound set to 10 per cent higher than the proposed optimal value of |$\kappa$| given in (3.12). Since granso may generate iterates violating these bounds or the other bound constraints in (4.4), we defined |$v$| to be zero if |$\kappa\leqslant 0$| and replaced |$\beta_{i}$| in (2.15) by pi/2, the 16-digit rounded value of |$\pi/2$|⁠, if |$\beta_{i}$| exceeds pi/2, to avoid the discontinuity in the tangent function at |$\pi/2$| (note that tan(pi/2)|$\approx 1.6\times 10^{16}$| has the desired positive sign). Because of the difficulty of the problem, we ran the code from many randomly generated starting points, with the initial values for |$\kappa$|⁠, |$\alpha_{i}$| and |$\beta_{i}$| generated from the uniform distribution on |$[0,\kappa^{\max}]$|⁠, |$[0,1]$| and |$[0,\pi/2],$| respectively, with the |$\alpha_{i}$| then sorted into increasing order. Fig. 9 Open in new tabDownload slide Summary of results for solving (4.4) with |$n=2$|⁠, running granso from 1000 randomly generated starting points. Of the 1000 candidate solutions obtained, 734 satisfied the bound and coarse grid stability constraints imposed by granso, and of these, 691 also passed the fine grid stability test. The top panel in the figure shows the computed optimal loads |$\kappa$| for the best 500 of these feasible solutions, sorted into decreasing order; the top 100 final values all agree with |$\kappa_{0}+\pi$| to 4 digits, while the top two final values agree with |$\kappa_{0} + \pi$| to 10 digits. The second and third panels show the associated final values of |$\alpha_{1}$| and |$\beta_{1}$| computed by these same 500 runs. The fourth panel shows the eigenvalues of the final associated matrix |${\bf M}(\alpha_{1},\beta_{1},\kappa)$|⁠. The computed solutions clearly separate into four flavors; see the text for details. Fig. 9 Open in new tabDownload slide Summary of results for solving (4.4) with |$n=2$|⁠, running granso from 1000 randomly generated starting points. Of the 1000 candidate solutions obtained, 734 satisfied the bound and coarse grid stability constraints imposed by granso, and of these, 691 also passed the fine grid stability test. The top panel in the figure shows the computed optimal loads |$\kappa$| for the best 500 of these feasible solutions, sorted into decreasing order; the top 100 final values all agree with |$\kappa_{0}+\pi$| to 4 digits, while the top two final values agree with |$\kappa_{0} + \pi$| to 10 digits. The second and third panels show the associated final values of |$\alpha_{1}$| and |$\beta_{1}$| computed by these same 500 runs. The fourth panel shows the eigenvalues of the final associated matrix |${\bf M}(\alpha_{1},\beta_{1},\kappa)$|⁠. The computed solutions clearly separate into four flavors; see the text for details. 4.1 Results for |$n=2$| Our analytical discussion of the case |$n=2$| was given in section 3; the results here strongly support our claim that the optimal configuration is given by (3.9), (3.10). Figure 9 shows the results obtained by running granso from 1000 randomly generated starting points. Of the 1000 candidate solutions generated by granso, 734 satisfied the bound and coarse grid stability constraints imposed by granso, and of these, 691 also passed the fine grid stability test described above. The top panel in the figure shows the computed optimal loads |$\kappa$| for the best 500 of these feasible solutions, sorted into decreasing order, while the second and third panels show the associated final values of |$\alpha_{1}$| and |$\beta_{1}$| computed by these same 500 runs. The fourth panel shows the eigenvalues of the final associated matrix |${\bf M}(\alpha_{1},\beta_{1},\kappa)$|⁠. The highest two final values of |$\kappa$| agree with each other, and with the value |$\kappa_0 + \pi$| given in (3.9) and (3.12), to 10 digits. The final values for |$\alpha_{1}$| and |$\beta_{1}$| for these same two best results agree with the value |$\kappa_0[\kappa_0+\pi]^{-1}$| (given in (3.9) and (3.13)) and |$\pi/2$|⁠, to 10 and 12 digits, respectively. It is also worth noting that the top 100 final values for the computed optimal load agree with |$\kappa_{0}+\pi$| to 4 digits. Looking at all four panels of Fig. 9, we see that the top 500 results come in several clearly distinct flavors. The first flavor is exhibited by the best 180 or so runs which all give good approximations to |$\kappa_{0}+\pi$|⁠. However, starting with the 284th result, we find a very different second flavor: many runs find that the computed optimal load is about |$\kappa=4.493$|⁠, which agrees with |$\kappa_{0}$|⁠, the square root of the critical load for the Dzhanelidze column, to four digits. Clearly, this is a locally maximal value for (4.4); otherwise, it would not be found so frequently. If we look at the associated computed |$\alpha_{1}$| and |$\beta_{1}$| values, usually |$\alpha_{1}$| is close to zero, but if not, then |$\beta_{1}$| is close to zero. It is easily checked that, regardless of the value of |$\beta_{1}$|⁠, if |$\alpha_{1}=0$| then |${\bf M}(\alpha_{1},\beta,\kappa_{0})$| is the zero matrix, with a double semisimple zero eigenvalue, so this locally optimal parameter configuration, like the apparent globally optimal configuration (3.9), (3.10), is on both the flutter and divergence boundaries. Physically, this corresponds to the mass |$M_{1}$| being fixed at the clamped end of the column. On the other hand, regardless of the value of |$\alpha_{1}$|⁠, if |$\beta_{1}=0$|⁠, then |${\bf M}(\alpha_{1},\beta_{1},\kappa_{0})$| has all zero entries except for |${\bf M}_{12}$|⁠, and hence has a double zero eigenvalue with a Jordan block. Again, this parameter configuration is on both the flutter and divergence boundaries, and physically, it corresponds to the mass |$M_{1}$| being zero. Note that the computed eigenvalues for this second flavor of solutions are relatively small. The third flavor of results is exhibited by the results numbered approximately 180 to 280. In these cases, granso terminated prematurely, without approximating a globally or locally maximal value, and we can see also that, on average, the larger |$\kappa$| is, the closer |$\alpha_{1}$| is to |$\kappa_0[\kappa_0+\pi]^{-1}$|⁠. Investigation of these cases shows that termination occurs because of the discontinuity in the stability constraint that we described above. This is also supported by the enormous associated eigenvalues of |${\bf M}$| shown in the fourth panel. Note also that as |$\kappa$| increases towards its optimal value, these eigenvalues decrease, but they neither converge to specific values, nor do they become very small. In fact, the matrix |${\bf M}$| associated with the best computed optimal |$\kappa$| is $$\begin{bmatrix} 5.4382\times 10^{1} & -4.0893\times 10^{-10}\\ 1.3020\times 10^{12} & -6.8246 \times 10^{-1} \end{bmatrix}\!.$$ Although its eigenvalues are not close to each other or to zero, they are small relative to the norm of the matrix, and their associated right eigenvectors are almost identical, indicating the nearby presence of a double eigenvalue. Furthermore, the diagonal and upper triangular elements are very small compared to the norm of the matrix, implying that a relatively small perturbation removing them yields a Jordan block with a double zero eigenvalue. Finally there is a fourth flavor of results: those that did not even reach a good approximation to the locally optimal value |$\kappa_{0}$|⁠. A final comment on Fig. 9: the granso termination codes are plotted at the bottom of the first panel. The value 1 means that granso terminated because the limit of 500 iterations was reached, while the value 2 means that it terminated because it could not find a higher feasible value for the load. Observe that the latter termination always occurred for the runs which approximated the apparent globally optimal load |$\kappa_{0}+\pi$| well (the first flavor) and the runs that obtained loads higher than |$\kappa_{0}$| but terminated without reaching a good approximation to |$\kappa_{0}+\pi$| (the third flavor). Thus, increasing the iteration limit would not have improved any of these values. On the other hand, the runs that provided a good approximation to the locally maximal value |$\kappa_{0}$| (the second flavor) or terminated before reaching that value (the fourth flavor) sometimes, but not always, terminated by exceeding the maximum iteration limit. The physical interpretation of the proposed supremum (3.9), (3.10) is that the mass |$M_{2}$| mounted on the free end of the column is zero in the limit |$\kappa\to\kappa^{*}$| (assuming that |$M_{1}$| is bounded above). It is interesting to consider what happens if we disallow this case, putting an upper limit on |$\mu_{1}=M_{1}/M_{2}$|⁠. Figure 10 shows the results when we introduce the constraint |$\mu_{1}\leqslant 100$| (left) or |$\mu_{1}\leqslant 10$| (right) by limiting |$\beta_{1}$| to |$\arctan(100)$| or |$\arctan(10)$| respectively. For |$\mu_{1}\leqslant 100$|⁠, we now find an optimal load of |$7.6287$|⁠, and for |$\mu_{1}\leqslant 10$|⁠, we find the optimal load |$7.4666$|⁠, which are respectively just 0.1 per cent and 2 per cent lower than the apparent unconstrained supremum |$\kappa_{0}+\pi$|⁠. The biggest difference we observe from comparing Fig. 10 with Fig. 9 is that the eigenvalues of the final computed |${\bf M}$| are now dramatically reduced, from more than |$10^{16}$| to less than 100 and 15, respectively. Thus, we obtain an only slightly reduced optimal load while introducing a much more physically reasonable model with much better numerical properties. Fig. 10 Open in new tabDownload slide Solving (4.4) for |$n=2$| with the additional constraint (left) |$\beta_{1}\leqslant \arctan(100)$| and (right) |$\beta_{1}\leqslant \arctan(10)$|⁠. The best computed optimal loads are respectively just 0.1 per cent and 2 per cent lower than the apparent unconstrained supremum |$\kappa_{0}+\pi$|⁠. Note the dramatic reduction in the size of the eigenvalues of the final computed |${\bf M}$| compared to Fig. 9. See the caption of Fig. 9 and the accompanying text for more details. Fig. 10 Open in new tabDownload slide Solving (4.4) for |$n=2$| with the additional constraint (left) |$\beta_{1}\leqslant \arctan(100)$| and (right) |$\beta_{1}\leqslant \arctan(10)$|⁠. The best computed optimal loads are respectively just 0.1 per cent and 2 per cent lower than the apparent unconstrained supremum |$\kappa_{0}+\pi$|⁠. Note the dramatic reduction in the size of the eigenvalues of the final computed |${\bf M}$| compared to Fig. 9. See the caption of Fig. 9 and the accompanying text for more details. Fig. 11 Open in new tabDownload slide Solving (4.4) for |$n=3$| with (left) no additional constraints and (right) with the constraints |$\alpha_{1}=\alpha_{1}^{*}=\kappa_0[\kappa_0+2\pi]^{-1}$|⁠, |$\beta_{1}=$|pi/2. See the accompanying text for more details. Fig. 11 Open in new tabDownload slide Solving (4.4) for |$n=3$| with (left) no additional constraints and (right) with the constraints |$\alpha_{1}=\alpha_{1}^{*}=\kappa_0[\kappa_0+2\pi]^{-1}$|⁠, |$\beta_{1}=$|pi/2. See the accompanying text for more details. Fig. 12 Open in new tabDownload slide Solving (4.4) for |$n=3$| with the additional constraints (left) |$\beta_{i}\leqslant \arctan(100)$|⁠, |$i=1,2$| and (right) |$\beta_{i}\leqslant \arctan(10)$|⁠, |$i=1,2$|⁠. Fig. 12 Open in new tabDownload slide Solving (4.4) for |$n=3$| with the additional constraints (left) |$\beta_{i}\leqslant \arctan(100)$|⁠, |$i=1,2$| and (right) |$\beta_{i}\leqslant \arctan(10)$|⁠, |$i=1,2$|⁠. 4.2 Results for |$n=3$| The left panel in Fig. 11 shows the results for solving (4.4) for |$n=3$|⁠. There are five variables: |$\alpha_{1},\alpha_{2},\beta_{1},\beta_{2}$| and |$\kappa$|⁠. Of the 1000 candidate solutions generated by granso, 517 satisfied the bound and coarse grid stability constraints imposed by granso, and of these, 416 passed the fine grid test. We see immediately that the problem for |$n=3$| is significantly harder than for |$n=2$|⁠, with not many runs approximating the conjectured optimal value well. Nonetheless, the two best runs generate |$\kappa\approx 1.0776$|⁠, which agrees with the conjectured optimal value |$\kappa_{0}+2\pi$| to five digits. These two runs also generate |$\alpha_{1}\approx 0.4169$| and |$\beta_{1}\approx 1.570796$| which agree with the conjectured optimal values |$\alpha_{1}^{*}=\kappa_0[\kappa_0+2\pi]^{-1}$| and |$\pi/2$| to 4 and 7 digits, respectively. The right panel in the same figure shows the results when we fix |$\alpha_{1}= \alpha_{1}^{*}$| and |$\beta_{1}=$|pi/2 (the 16 digit rounded value of |$\pi/2$|⁠) and optimize over the remaining three variables |$\alpha_{2}$|⁠, |$\beta_{2}$| and |$\kappa$|⁠. Then the best two runs generate |$\kappa$| agreeing with |$\kappa_{0}+2\pi$| to 12 digits, and the best 100 runs agree with this to 10 digits. Together, the results reported in the left and right panels of Fig. 11 make a convincing argument that the values shown in (3.12) and (3.13) are indeed the supremal value |$\kappa^{*}$| and the corresponding limiting value |$\alpha_{1}^{*}$| when |$n=3$|⁠, and that the corresponding limiting value |$\beta_{1}^{*}$| is again |$\pi/2$|⁠. Although we do not have a conjectured formula for |$\alpha_{2}^{*}$|⁠, its computed optimal value is 0.7085. Furthermore, the limiting value |$\beta_{2}^{*}$| is again apparently |$\pi/2$|⁠, meaning the mass ratio |$\mu_{2}=M_{2}/M_{3} \to\infty$| as |$\kappa\to\kappa^{*}$|⁠, which indeed must be the case assuming the limiting value of |$M_{2}$| is nonzero and |$M_{1}$| is bounded above, since then |$\beta_{1}\to\pi/2$| implies that |$M_{3}\to 0$|⁠. Figure 12 shows the results when we introduce the mass ratio constraint |$\mu_{i}\leqslant 100$| (left) or |$\mu_{i}\leqslant 10$| (right) by limiting |$\beta_{i}\leqslant \arctan(100)$|⁠, |$i=1,2$| or |$\beta_{i} \leqslant \arctan(10)$|⁠, |$i=1,2,$| respectively. For |$\mu_{i}\leqslant 100$|⁠, we now find an optimal load of 10.589, which is only 0.5 per cent lower than |$\kappa_{0}+2\pi$|⁠. However, when we constrain |$\mu_{i}\leqslant 10$|⁠, the best optimal load found is only 7.59, which is a 30 per cent reduction from |$\kappa_{0}+2\pi$|⁠. When we repeat these runs with 10,000 starting points instead of 1000, these numbers are only slightly improved. Fig. 13 Open in new tabDownload slide (Left) solving (4.4) for |$n=4$| with the constraints |$\alpha_{1}=\alpha_{1}^{*}=\kappa_0[\kappa_0+3\pi]^{-1}$|⁠, |$\beta_{1}=$|pi/2 and (right) solving (4.4) for |$n=5$| with the constraints |$\alpha_{1}=\alpha_{1}^{*}=\kappa_0[\kappa_0+4\pi]^{-1}$|⁠, |$\beta_{1}=$|pi/2. Fig. 13 Open in new tabDownload slide (Left) solving (4.4) for |$n=4$| with the constraints |$\alpha_{1}=\alpha_{1}^{*}=\kappa_0[\kappa_0+3\pi]^{-1}$|⁠, |$\beta_{1}=$|pi/2 and (right) solving (4.4) for |$n=5$| with the constraints |$\alpha_{1}=\alpha_{1}^{*}=\kappa_0[\kappa_0+4\pi]^{-1}$|⁠, |$\beta_{1}=$|pi/2. 4.3 Results for |$n=4$| and |$n=5$| The optimization problem is so much harder for |$n=4$| and |$n=5$| that we needed 10,000 starting points to get good results, even when we set |$\alpha_{1}$| to its conjectured optimal value in (3.13) and |$\beta_{1}$| to pi/2, optimizing over the remaining 5 and 7 variables, respectively. The results are shown in the left and right panels of Fig. 13. For |$n=4$|⁠, the best 5 results agree with our conjectured optimal load |$\kappa_{0} + 3\pi$| to eight digits, while for |$n=5$|⁠, the best 25 results agree with |$\kappa_{0}+4\pi$| to seven digits. These results strongly support our conjecture regarding the supremal load |$\kappa$| for |$n$| masses given in (3.12). 5. Concluding remarks We believe we have made a convincing case that the supremal load for the strongest stable massless column with a follower load and |$n$| relocatable masses is, in the dimensionless model defined in section 2, |$\kappa_{0}+(n-1)\pi$|⁠, where |$\kappa_{0}$| is the smallest positive root of |$\tan(\kappa)=\kappa$|⁠. This conjecture has not previously appeared in the literature as far as we know, except in the case |$n=1$| where it has been known to be true for decades (6). We have given a detailed analytical derivation of this result for |$n=2$| and presented extensive computational results that support it for |$n=2,3,4,5$|⁠, using numerical nonsmooth optimization. With this model problem effectively solved, we believe it would be interesting to apply our nonsmooth optimization techniques to more realistic columns with follower loads (29), such as the Beck, Pflüger and Leipholz columns with a single free end as well as to free–free beams both with distributed and concentrated masses to get new insights about the nature of the optimal solution to these long-standing optimization problems. We believe it is also important to consider extending traditional stability constraints to more robust stability constraints based on pseudospectra (66), a topic that is beyond the scope of this article. Acknowledgements The authors thank Tim Mitchell, the author of granso, for many helpful discussions and suggestions regarding the formulation of the stability constraint. They also thank the London Mathematical Society for supporting the second author’s visit to Northumbria through the Scheme 4 Research in Pairs grant No 41820. The second author was supported in part by the U.S. National Science Foundation Grant DMS-2012250. A. Pflüger’s column To model the Pflüger column, we consider an elastic beam of length |$l$|⁠, with Young’s modulus |$E$| and mass per unit length |$m$|⁠, clamped at one end and loaded by a tangential follower force |$P$| at the other end, where a point mass |$M$| is mounted. The moment of inertia of a cross-section of the column is denoted by |$I$|⁠. Small lateral vibrations of the Pflüger column near the undeformed equilibrium are described by the linear partial differential equation (23, 25) $$\begin{equation}\label{pce} EI\frac{\partial^4 y}{\partial s^4}+P\frac{\partial^2 y}{\partial s^2}+m\frac{\partial^2 y}{\partial t^2}=0, \end{equation}$$(A.1) where |$y(s,t)$|⁠, is the amplitude of the vibrations and |$s \in [0, l]$| is a coordinate along the column. At the clamped end |$(s = 0),$| (A.1) satisfies the boundary conditions $$\begin{equation}\label{pbc1} y=0, \quad \frac{\partial y}{\partial s}=0,\quad s=0, \end{equation}$$(A.2) while at the loaded end |$(s = l)$|⁠, the boundary conditions are $$\begin{equation}\label{pbc2} EI\frac{\partial^2 y}{\partial s^2}=0,\quad EI\frac{\partial^3 y}{\partial s^3}=M\frac{\partial^2 y}{\partial t^2}, \quad s=l. \end{equation}$$(A.3) Introducing the dimensionless quantities $$\begin{equation}\label{dlq} \xi=\frac{s}{l},\quad \tau=\frac{t}{l^2}\sqrt{\frac{EI}{m}},\quad p=\frac{P l^2}{EI},\quad \mu=\frac{M}{ml}, \end{equation}$$(A.4) and separating the time variable through |$y(\xi, \tau) = lf(\xi)\exp(\lambda \tau)$|⁠, we obtain the dimensionless boundary eigenvalue problem $$\begin{equation}\label{dlbp} \partial_{\xi}^4 f +p \partial_{\xi}^2 f +\lambda^2 f=0, \end{equation}$$(A.5) $$\begin{eqnarray}\label{dlbc} &\partial_{\xi}^2 f(1)=0,\quad \partial_{\xi}^3 f(1)=\mu \lambda^2 f(1),&\nonumber \\ &f(0)=0,\quad \partial_{\xi}f(0)=0& \end{eqnarray}$$(A.6) defined on the interval |$\xi \in [0, 1]$|⁠. A solution to (A.5) with boundary conditions (A.6) is (23, 25) $$\begin{equation}\label{sbp} f(\xi)=A(\cosh(g_2 \xi)-\cos(g_1 \xi))+B(g_1\sinh(g_2 \xi)-g_2\sin(g_1 \xi)) \end{equation}$$(A.7) with $$ g_{1,2}=\sqrt{\frac{\sqrt{p^2-4\lambda^2}\pm p}{2}},$$ where the subscripts |$1$| and |$2$| correspond to the signs |$+$| and |$-$|⁠, respectively. Imposing the boundary conditions (A.6) on the solution (A.7) yields the characteristic equation |$\Delta(\lambda) = 0$| for the determination of the eigenvalues |$\lambda$|⁠, where $$\Delta(\lambda)=\Delta_1-\Delta_2\mu \lambda^2$$ and $$\begin{eqnarray} \Delta_1&=&g_1g_2(g_1^4+g_2^4+2g_1^2g_2^2\cosh g_2 \cos g_1 + g_1 g_2 (g_1^2-g_2^2)\sinh g_2 \sin g_1),\nonumber\\ \Delta_2&=&(g_1^2+g_2^2)(g_1 \sinh g_2 \cos g_1 - g_2 \cosh g_2 \sin g_1). \end{eqnarray}$$(A.8) Parameterizing the mass ratio in (A.4) by |$\mu = \tan \beta$| with |$\beta \in [0, \pi/2]$| enables the exploration of all possible ratios |$\mu=M/(ml)$| of the end mass to the mass of the column from zero (⁠|$\beta = 0$|⁠) to infinity (⁠|$\beta = \pi/2$|⁠). The former case, without the end mass, corresponds to the Beck column, whereas the latter corresponds to a massless rod with an end mass, which is known as the Dzhanelidze column (6). It is well-established that the uniform Beck column is stable against flutter if the dimensionless follower force, |$p$|⁠, is such that |$0\le p \lesssim 20.05$|⁠, (1, 2, 6, 7). In contrast, the Dzhanelidze column becomes unstable at |$p \approx 20.19$|⁠, which is the smallest positive root of (6) $$\begin{equation}\label{dcp} \tan \sqrt{p} = \sqrt{p}. \end{equation}$$(A.9) These values, representing two extreme situations, are connected by a marginal stability curve in the |$(\beta, p)$|-plane (6, 23, 25–28); see the right panel of Fig. A.1. Fig. A.1 Open in new tabDownload slide The Pflüger column and its stability diagram. The ratio of the end mass to the mass of the column, |$\mu=M/(ml)$|⁠, is parameterized by |$\mu=\tan\beta$|⁠. The Beck column corresponds to the vanishing end mass (⁠|$M=0$|⁠, so |$\beta=0$|⁠) and the massless Pflüger column (or Dzhanelidze’s column (6)) to the vanishing mass of the rod (⁠|$m=0$|⁠, so |$\beta=\pi/2$|⁠). The vertical axis of the stability diagram shows the dimensionless load |$p=\frac{Pl^2}{EI}$|⁠, where |$E$| is Young’s modulus, |$I$| is the moment of inertia of a cross-section of the column and |$l$| is the length of the column. Fig. A.1 Open in new tabDownload slide The Pflüger column and its stability diagram. The ratio of the end mass to the mass of the column, |$\mu=M/(ml)$|⁠, is parameterized by |$\mu=\tan\beta$|⁠. The Beck column corresponds to the vanishing end mass (⁠|$M=0$|⁠, so |$\beta=0$|⁠) and the massless Pflüger column (or Dzhanelidze’s column (6)) to the vanishing mass of the rod (⁠|$m=0$|⁠, so |$\beta=\pi/2$|⁠). The vertical axis of the stability diagram shows the dimensionless load |$p=\frac{Pl^2}{EI}$|⁠, where |$E$| is Young’s modulus, |$I$| is the moment of inertia of a cross-section of the column and |$l$| is the length of the column. For every fixed value |$\beta \in [0, \pi/2)$|⁠, the Pflüger column loses stability via flutter when an increase in |$p$| causes the imaginary eigenvalues of two different modes to approach each other and merge into a double imaginary eigenvalue with a Jordan block (that is, with algebraic multiplicity two and geometric multiplicity one). When |$p$| crosses the threshold, the double eigenvalue splits into two complex eigenvalues, one with positive real part, which determines a flutter-unstable mode. At |$\beta = \pi/2$| the stability boundary of the Pflüger column has a vertical tangent and the type of instability changes from flutter to divergence, that is, nonoscillatory growth of a mode corresponding to a positive real eigenvalue, for |$p \gtrsim 20.19$|⁠; see (6, 26, 27). References 1. Beck M. , Die Knicklast des einseitig eingespannten, tangential gedrückten Stäbes , Z. Angew. Math. Mech. 3 ( 1952 ) 225 – 228 . Google Scholar OpenURL Placeholder Text WorldCat 2. Carr J. and Malhardeen M. Z. M., Beck’s problem , SIAM J. Appl. Math. 37 ( 1979 ) 261 – 262 . Google Scholar Crossref Search ADS WorldCat 3. Sugiyama Y. , Langthjem M. and Katayama K., Dynamic Stability of Columns under Nonconservative Forces: Theory and Experiment. Solid Mechanics and Its Applications, Vol. 255 ( Springer , Berlin 2019 ). Google Scholar Crossref Search ADS Google Preview WorldCat COPAC 4. Ziegler H. , Linear elastic stability. A critical analysis of methods, First part , ZAMP Z. Angew. Math. Phys. 4 ( 1953 ) 89 – 121 . Google Scholar Crossref Search ADS WorldCat 5. Ziegler H. , Linear elastic stability. A critical analysis of methods, Second part , ZAMP Z. Angew. Math. Phys. 4 ( 1953 ) 167 – 185 . Google Scholar Crossref Search ADS WorldCat 6. Bolotin V. V. , Nonconservative Problems of the Theory of Elastic Stability ( Pergamon Press , Oxford 1963 ). Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC 7. Kirillov O. N. , Nonconservative Stability Problems of Modern Physics . De Gruyter Series in Mathematical Physics , Vol. 14 ( De Gruyter , Berlin, Boston 2013 ). Google Scholar Crossref Search ADS Google Preview WorldCat COPAC 8. Bayly P. V. and Dutcher S. K., Steady dynein forces induce flutter instability and propagating waves in mathematical models of flagella , J. R. Soc. Interface 13 ( 2016 ) 20160523 . Google Scholar Crossref Search ADS PubMed WorldCat 9. De Canio G. , Lauga E. and Goldstein R. E., Spontaneous oscillations of elastic filaments induced by molecular motors , J. R. Soc. Interface 14 ( 2017 ) 20170491 . Google Scholar Crossref Search ADS PubMed WorldCat 10. Fatehiboroujeni S. , Gopinath A. and Goyal S., Three-dimensional nonlinear dynamics of prestressed active filaments: flapping, swirling, and flipping , Phys. Rev. E 103 ( 2021 ) 013005 . Google Scholar Crossref Search ADS PubMed WorldCat 11. Zhu L. and Stone H. A., Propulsion driven by self-oscillation via an electrohydrodynamic instability , Phys. Rev. Fluids 4 ( 2019 ) 061701 . Google Scholar Crossref Search ADS WorldCat 12. Zhu L. and Stone H. A., Harnessing elasticity to generate self-oscillation via an electrohydrodynamic instability , J. Fluid Mech. 888 ( 2020 ) A31 . Google Scholar Crossref Search ADS WorldCat 13. Sugiyama Y. , Langthjem M. and Ryu B. J., Realistic follower forces , J. Sound Vibr. 225 ( 1999 ) 779 – 782 . Google Scholar Crossref Search ADS WorldCat 14. Sugiyama Y. , Langthjem M. and Ryu B.-J., Beck’s column as the ugly duckling , J. Sound Vibr. 254 ( 2002 ) 407 – 410 . Google Scholar Crossref Search ADS WorldCat 15. Sundararajan C. , Optimization of a nonconservative elastic system with stability constraint , J. Opt. Theory Appl. 16 ( 1975 ) 355 – 378 . Google Scholar Crossref Search ADS WorldCat 16. Park Y. P. and Mote C. D., The maximum controlled follower force on a free-free beam carrying a concentrated mass , J. Sound. Vibr. 98 ( 1985 ) 247 – 256 . Google Scholar Crossref Search ADS WorldCat 17. Kirillov O. N. and Seyranian A. P., Optimization of stability of a flexible missile under follower thrust , AIAA Paper 98-4969 ( 1998 ) 2063 – 2073 . Google Scholar OpenURL Placeholder Text WorldCat 18. Sugiyama Y. , Matsuike J., Ryu B.-T., Katayama K., Kinoi S. and Enomoto N., Effect of concentrated mass on stability of cantilevers under rocket thrust , AIAA J. 33 ( 1995 ) 499 – 503 . Google Scholar Crossref Search ADS WorldCat 19. Sugiyama Y. , Langthjem M. A., Iwama T., Kobayashi M., Katayama K. and Yutani H., Shape optimization of cantilevered columns subjected to a rocket-based follower force and its experimental verification , Struct. Multidisc. Opt. 46 ( 2012 ) 829 – 838 . Google Scholar Crossref Search ADS WorldCat 20. Bigoni D. and Noselli G., Experimental evidence of flutter and divergence instabilities induced by dry friction , J. Mech. Phys. Sol. 59 ( 2011 ) 2208 – 2226 . Google Scholar Crossref Search ADS WorldCat 21. Bigoni D. , Kirillov O. N., Misseroni D., Noselli G. and Tommasini M., Flutter and divergence instability in the Pflüger column: experimental evidence of the Ziegler destabilization paradox , J. Mech. Phys. Sol. 116 ( 2018 ) 99 – 116 . Google Scholar Crossref Search ADS WorldCat 22. Bigoni D. , Misseroni D., Tommasini M., Kirillov O. N. and Noselli G., Detecting singular weak-dissipation limit for flutter onset in reversible systems , Phys. Rev. E 97 ( 2018 ) 023003 . Google Scholar Crossref Search ADS PubMed WorldCat 23. Pflüger A. , Zur Stabilität des tangential gedrückten Stäbes , Z. Angew. Math. Mech. 35 ( 1955 ) 191 . Google Scholar Crossref Search ADS WorldCat 24. Deineko K. S. and Leonov M. I., A dynamic method for the investigation of the stability of a compressed bar , Prikl. Mat. Mekh. 19 ( 1955 ) 738 – 744 . Google Scholar OpenURL Placeholder Text WorldCat 25. Tommasini M. , Kirillov O. N., Misseroni D. and Bigoni D., The destabilizing effect of external damping: singular flutter boundary for the Pflüger column with vanishing external dissipation , J. Mech. Phys. Sol. 91 ( 2016 ) 204 – 215 . Google Scholar Crossref Search ADS WorldCat 26. Oran C. , On the significance of a type of divergence , J. Appl. Mech. 39 ( 1972 ) 263 – 265 . Google Scholar Crossref Search ADS WorldCat 27. Sugiyama Y. , Kashima K. and Kawagoe H., On an unduly simplified model in the non-conservative problems of elastic stability , J. Sound. Vib. 45 ( 1976 ), 237 – 247 . Google Scholar Crossref Search ADS WorldCat 28. Chen L. W. and Ku D. W., Eigenvalue sensitivity in the stability analysis of Beck’s column with a concentrated mass at the free end , J. Sound Vib. 153 ( 1992 ) 403 – 411 . Google Scholar Crossref Search ADS WorldCat 29. Gajewski A. and Zyczkowski M., Optimal structural design under stability constraints , ( Kluwer , Dordrecht 1988 ). Google Scholar Crossref Search ADS Google Preview WorldCat COPAC 30. Claudon J. L. , Characteristic curves and optimum design of two structures subjected to circulatory loads , J. de Mecanique 14 ( 1975 ) 531 – 543 . Google Scholar OpenURL Placeholder Text WorldCat 31. Hanaoka M. and Washizu K., Optimum design of Beck’s column , Comp. Struct. 11 ( 1980 ) 473 – 480 . Google Scholar Crossref Search ADS WorldCat 32. Kounadis A. N. and Katsikadelis J. T., On the discontinuity of the flutter load for various types of cantilevers , Intern. J. Solids Struct. 16 ( 1980 ) 375 – 383 . Google Scholar Crossref Search ADS WorldCat 33. Bogacz R. and Frischmuth K., On optimality of column geometry , Arch. Appl. Mech. 88 ( 2018 ) 317 – 327 . Google Scholar Crossref Search ADS WorldCat 34. Mahrenholtz O. and Bogacz R., On the shape of characteristic curves for optimal structures under non-conservative loads , Arch. Appl. Mech. 50 ( 1981 ) 141 – 148 . Google Scholar OpenURL Placeholder Text WorldCat 35. Langthjem M. A. and Sugiyama Y., Optimum design of cantilevered columns under the combined action of conservative and nonconservative loads Part I: The undamped case , Comp. Struct. 74 ( 2000 ) 385 – 398 . Google Scholar Crossref Search ADS WorldCat 36. Ringertz U. T. , On the design of Beck’s column , Struct. Opt. 8 ( 1994 ) 120 – 124 . Google Scholar Crossref Search ADS WorldCat 37. Temis Yu M. and Fedorov I. M., Shape optimization of nonconservatively loaded beams with a stability criterion , Probl. Strength Plast. 69 ( 2007 ) 24 – 37 . Google Scholar Crossref Search ADS WorldCat 38. Kirillov O. N. and Seyranian A. P., Metamorphoses of characteristic curves in circulatory systems , J. Appl. Math. Mech. 66 ( 2002 ), 371 – 385 . Google Scholar Crossref Search ADS WorldCat 39. Kirillov O. N. and Seyranian A. P., A non-smooth optimization problem , Moscow Univ. Mech. Bull. 57 ( 2002 ) 1 – 6 . Google Scholar OpenURL Placeholder Text WorldCat 40. O’Reilly O. M. , N. K. Malhotra and N. S. Namachchivaya, Some aspects of destabilization in reversible dynamical systems with application to follower forces , Nonlinear Dyn. 10 ( 1996 ) 63 – 87 . Google Scholar Crossref Search ADS WorldCat 41. Seyranian A. P. and Kirillov O. N., Bifurcation diagrams and stability boundaries of circulatory systems , Theor. Appl. Mech. 26 ( 2001 ) 135 – 168 . Google Scholar OpenURL Placeholder Text WorldCat 42. Kirillov O. N. and Seyranian A. P., Collapse of the Keldysh chains and stability of continuous non-conservative systems , SIAM J. Appl. Math. 64 ( 2004 ) 1383 – 1407 . Google Scholar OpenURL Placeholder Text WorldCat 43. Kirillov O. N. and Overton M. L., Robust stability at the swallowtail singularity , Front. Phys. 1 ( 2013 ) 24 . Google Scholar Crossref Search ADS WorldCat 44. Kirillov O. N. , Singularities in structural optimization of the Ziegler pendulum , Acta Polytechn. 51 ( 2011 ) 32 – 43 . Google Scholar Crossref Search ADS WorldCat 45. Langthjem M. A. and Sugiyama Y., Optimum shape design against flutter of a cantilevered column with an end-mass of finite size subjected to a non-conservative load , J. Sound Vibr. 226 ( 1999 ) 1 – 23 . Google Scholar Crossref Search ADS WorldCat 46. Katsikadelis J. T. and Tsiatas G. C., Optimum design of structures subjected to follower forces , Intern. J. Mech. Sci. 49 ( 2007 ) 1204 – 1212 . Google Scholar Crossref Search ADS WorldCat 47. Kordas Z. and Zyczkowski M., On the loss of stability of a rod under a super-tangential force , Arch. Mech. Stos. 1 ( 1963 ) 7 – 31 . Google Scholar OpenURL Placeholder Text WorldCat 48. Lee H. P. , Flutter of a cantilever rod with a relocatable lumped mass , Comput. Methods Appl. Mech. Eng. 144 ( 1997 ) 23 – 31 . Google Scholar Crossref Search ADS WorldCat 49. Hauger W. , Bemerkungen zu dem Einfluss der Massenverteilung bei nichtkonservativen Stabitätsproblemen elastischer Stäbe , Ing.-Arch. 35 ( 1967 ) 283 – 291 . Google Scholar Crossref Search ADS WorldCat 50. Leipholz H. and Lindner G., Über den Einflüss der Massenverteilung auf das nichtkonservative Knicken von Stäben , Ing.-Arch. 39 ( 1970 ) 187 – 194 . Google Scholar Crossref Search ADS WorldCat 51. Kapoor R. N. and Leipholz H., Stability analysis of a damped polygenic systems with relocatable mass along its length , Ing-Arch 43 ( 1974 ) 233 – 239 . Google Scholar Crossref Search ADS WorldCat 52. Kapoor R. N. and Leipholz H., On mass distribution, rotary inertia and external damping of a viscoelastic polygenic system , Z. Angew. Math. Mech. 54 ( 1974 ) 205 – 208 . Google Scholar Crossref Search ADS WorldCat 53. Kounadis A. N. , Stability of elastically restrained Timoshenko cantilevers with attached masses subjected to a follower force , ASME J. Appl. Mech. 44 ( 1977 ) 731 – 736 . Google Scholar Crossref Search ADS WorldCat 54. Cazzolli A. , Dal Corso F. and Bigoni D., Non-holonomic constraints inducing flutter instability in structures under conservative loadings , J. Mech. Phys. Solids 138 ( 2020 ) 103919 . Google Scholar Crossref Search ADS WorldCat 55. Bazant Z. P. and Cedolin L., Stability of Structures: Elastic, Inelastic, Fracture and Damage Theories ( World Scientific , Singapore 2010 ). Google Scholar Crossref Search ADS Google Preview WorldCat COPAC 56. Kagan-Rosenzweig L. M. , Quasi-static approach to non-conservative problems of the elastic stability theory , Int. J. Solids Struct . 38 ( 2001 ) 1341 – 1353 . OpenURL Placeholder Text WorldCat 57. Ingerle K. , Stability of massless non-conservative elastic systems , J. Sound Vibr. 332 ( 2013 ) 4529 – 4540 . OpenURL Placeholder Text WorldCat 58. Kagan-Rosenzweig L. M. , Topics in Nonconservative Stability Theory ( St.-Petersburg University Press , St. Petersburg, (In Russian) 2014 ). Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC 59. Curtis F. E. , Mitchell T. and Overton M. L., A BFGS-SQP method for nonsmooth, nonconvex, constrained optimization and its evaluation using relative minimization profiles , Optim. Methods Softw. 32 ( 2017 ) 148 – 181 . Google Scholar Crossref Search ADS WorldCat 60. Mitchell T. , GRANSO: GRadient-based Algorithm for Non-Smooth Optimization . http://www.timmitchell.com/software/GRANSO/ ( 2020 ). 61. Gallina P. , About the stability of non-conservative undamped systems , J. Sound Vibr. 262 ( 2003 ) 977 – 988 . Google Scholar Crossref Search ADS WorldCat 62. Bulatovic R. M. , A sufficient condition for instability of equilibrium of non-conservative undamped systems , Phys. Lett. A 375 ( 2011 ) 3826 – 3828 . Google Scholar Crossref Search ADS WorldCat 63. Driscoll T. A. , Hale N. and Trefethen L. N., Chebfun Guide ( Pafnuty Publications , Oxford 2014 ). Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC 64. Overton M. L. , Stability optimization for polynomials and matrices, Nonlinear Physical Systems: Spectral Analysis, Stability and Bifurcations (eds Kirillov O. & Pelinovsky D.; John Wiley & Sons Inc. , New York 2014 ) 351 – 375 . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC 65. Greenbaum A. , Li R.-C. and Overton M. L., First-order perturbation theory for eigenvalues and eigenvectors , SIAM Rev. 62 ( 2020 ) 463 – 482 . OpenURL Placeholder Text WorldCat 66. Trefethen L. N. and Embree M., Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators ( Princeton University Press , Princeton, NJ 2006 ). Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC © The Author, 2021. Published by Oxford University Press; all rights reserved. For Permissions, please email: journals.permissions@oup.com. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Finding the strongest stable massless column with a follower load and relocatable concentrated masses
JF - The Quarterly Journal of Mechanics and Applied Mathematics
DO - 10.1093/qjmam/hbab005
DA - 2021-05-10
UR - https://www.deepdyve.com/lp/oxford-university-press/finding-the-strongest-stable-massless-column-with-a-follower-load-and-8WCV88I1fh
SP - 1
EP - 1
VL - Advance Article
IS - 
DP - DeepDyve
ER -