Antoine Levitt
email: antoine.levitt@inria.fr.

Christoph Ortner
email: c.ortner@warwick.ac.uk.

Convergence and Cycling in Walker-type Saddle Search Algorithms

Algorithms for computing local minima of smooth objective functions enjoy a mature theory as well as robust and efficient implementations. By comparison, the theory and practice of saddle search is destitute. In this paper we present results for idealized versions of the dimer and gentlest ascent (GAD) saddle search algorithms that show-case the limitations of what is theoretically achievable within the current class of saddle search algorithms: (1) we present an improved estimate on the region of attraction of saddles; and (2) we construct quasi-periodic solutions which indicate that it is impossible to obtain globally convergent variants of dimer and GAD type algorithms.

Introduction

The first step in the exploration of a molecular energy landscape is usually the determination of energy minima, using an optimization algorithm. There exists a large number of such algorithms, backed by a rich and mature theory [START_REF] Nocedal | Numerical Optimization[END_REF][START_REF] Conn | Trust-Region Methods[END_REF]. Virtually all optimization algorithms in practical use today feature a variety of rigorous global and local convergence guarantees, and well-understood asymptotic rates.

As a second step, one typically determines the saddles between minima. They represent a crude description of the transitions between minima (reactions) and can be thought of as the edges in the graph between stable states of a molecule or material system. If neighboring minima are known, then methods of NEB or string type [START_REF] Jónsson | Nudged elastic band for finding minimum energy paths of transitions[END_REF][START_REF] Ren | String method for the study of rare events[END_REF] may be employed. On the other hand when only one minimum is known, then "walker methods" of the eigenvector-following methodology such as the dimer algorithm [START_REF] Henkelman | A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives[END_REF] are required. This second class of methods is the focus of the present work; for extensive reviews of the literature we refer to [START_REF] Olsen | Comparison of methods for finding saddle points without knowledge of the final states[END_REF][START_REF] Xiang | The gentlest ascent dynamics[END_REF][START_REF] Barkema | The activation-relation technique: an efficient algorithm for sampling energy landscapes[END_REF][START_REF] Gould | A dimer-type saddle search algorithm with preconditioning and linesearch[END_REF].

Since saddles represent reactions, the determination of saddle points is of fundamental importance in determining dynamical properties of an energy landscape, yet the state of the art of algorithms is very different from that for optimization: more than 15 years after the introduction of the dimer method [START_REF] Henkelman | A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives[END_REF] (the most widely used walker-type saddle search scheme), finding saddle points remains an art rather than a science. A common practice is to detect non-convergence and restart the algorithm with a different starting point. A mathematically rigorous convergence theory has only recently begun to emerge; see [START_REF] Zhang | Shrinking dimer dynamics and its applications to saddle point search[END_REF][START_REF] Gould | A dimer-type saddle search algorithm with preconditioning and linesearch[END_REF] and references therein. To the best of our knowledge all convergence results to date are local: convergence can only be guaranteed if an initial guess is sufficiently close to a (index-1) saddle. None of the existing saddle search algorithms come with the kind of global convergence guarantees that even the most basic optimization algorithms have.

The purpose of the present work is twofold: [START_REF] Barkema | The activation-relation technique: an efficient algorithm for sampling energy landscapes[END_REF] We strengthen existing local convergence results for dimer/GAD type saddle search methods by developing an improved estimate on the region of attraction of index-1 saddle points that goes beyond the linearized regime. [START_REF] Conn | Trust-Region Methods[END_REF] We produce new examples demonstrating generic cycling in those schemes, and pathological behavior of idealized versions of these algorithms. These results illustrate how fundamentally different saddle search is from optimization. They suggest that a major new idea is required to obtain globally convergent walker-type saddle search methods, and support the idea of stringof-state methods being more robust.

Local and global convergence in optimization

We consider the steepest descent method as a prototype optimization algorithm. Given an energy landscape E ∈ C 2 (R N ), the gradient descent dynamics (or gradient flow ) is ẋ = -∇E(x).

(

) 1
This ODE enjoys the property that

d dt E(x) = ẋ, ∇E(x) = -∇E(x) 2 .
If E is bounded from below, it follows that ∇E(x) → 0 and, under mild conditions (for instance, E coercive with non-degenerate critical points), x converges to a critical point, that is generically a minimum. This property can be transferred to the discrete iterates of the steepest descent method

x n+1 = x n -α n ∇E(x n ), (2)
under conditions on the step length α n (for instance the Armijo condition). In both cases, the crucial point for convergence is that E(x(t)) or E(x n ) is an objective function (also called merit or Lyapunov function) that decreases in time.

Eigenvector-following methods: the ISD and GAD

If x * is a non-degenerate index-1 saddle, then the symmetric Hessian matrix H * = ∇ 2 E(x * ) has one negative eigenvalue, while all other eigenvalues are positive. In this case, the steepest descent dynamics (1) is repelled away from x * along the mode corresponding to the negative eigenvalue.

To obtain a dynamical system for which x * is an attractive fixed point, we reverse the flow in the direction of the unstable mode. Let v 1 (x) be a normalized eigenvector corresponding to the smallest eigenvalue of ∇ 2 E(x), then for x -x * sufficiently small, the direction

-(I -2v 1 (x) ⊗ v 1 (x))∇E(x)
points towards the saddle x * . Note that this direction does not depend on the arbitrary sign of v 1 , and therefore in the rest of the paper we will talk of "the lowest eigenvector v 1 (x)" whenever the first eigenvalue of ∇ 2 E(x) is simple. This is the essence of the eigenvector-following methodology, which has many avatars (such as the dimer method [START_REF] Henkelman | A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives[END_REF], the Gentlest Ascent Dynamics [START_REF] Xiang | The gentlest ascent dynamics[END_REF], and numerous variants). In our analysis we will consider the simplest such method, which we will call the Idealized Saddle Dynamics (ISD),

ẋ = -(I -2v 1 (x) ⊗ v 1 (x))∇E(x). (3)
Under this dynamics, a linear stability analysis shows that non-degenerate index-1 saddle points are attractive, while non-degenerate minima, maxima or saddle points of index greater than 1 are repulsive (see Lemma 1). The ISD (3) is only well-defined when v 1 (x) is determined unambiguously, that is, when the first eigenvalue of ∇ 2 E(x) is simple. The singularities of this flow where ∇ 2 E(x) has repeated first eigenvalues will play an important role in this paper.

In practice, the orientation v 1 (x) has to be computed from ∇ 2 E(x). This makes the method unattractive for many applications in which the second derivative is not available or prohibitively expensive (for instance, ab initio potential surfaces, in which E(x) and ∇E(x) are readily computed but ∇ 2 E(x) requires a costly perturbation analysis). Because of this, the orientation is often relaxed and computed in alternation with the translation (3). A mathematically simple flavor of this approach is the Gentlest Ascent Dynamics (GAD):

[4] ẋ = -(I -2v ⊗ v)∇E(x), ε 2 v = -(I -v ⊗ v)∇ 2 E(x)v. (4)
At a fixed x, the dynamics for v is a gradient flow for the Rayleigh quotient v, ∇ 2 E(x)v on the unit sphere S 1 in R N , which converges to the lowest eigenvector v 1 (x). The parameter ε > 0 controls the speed of relaxation of v towards v 1 (x) relative to that of x. The ISD is formally obtained in the limit ε → 0. The practical advantage of the GAD (4) over the ISD (3) is that, once discretized in time, it can be implemented using only the action of ∇ 2 E(x) on a vector, which can be computed efficiently via finite differences. This is the basis of the dimer algorithm [START_REF] Henkelman | A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives[END_REF]. The ε scaling is analogous to common implementations of the dimer algorithm that adapt the number of rotations per step to ensure approximate equilibration of v.

Using linearized stability analysis one can prove local convergence of the ISD, GAD or dimer algorithms [START_REF] Zhang | Shrinking dimer dynamics and its applications to saddle point search[END_REF][START_REF] Gould | A dimer-type saddle search algorithm with preconditioning and linesearch[END_REF]. However, due to the absence of a global merit function as in optimization, there is no natural Armijo-like condition to choose the stepsizes in a robust manner, or indeed to obtain global convergence guarantees (however, see [START_REF] Gould | A dimer-type saddle search algorithm with preconditioning and linesearch[END_REF][START_REF] Gao | An iterative minimization formulation for saddle point search[END_REF] for ideas on the construction of local merit functions).

In this paper, we only study the ISD and GAD dynamics: we expect that the behavior we find applies to practical variants under appropriate conditions on the parameters (for instance, the dimer algorithm with a sufficiently small finite difference step and a sufficiently high number of rotation steps per translation step).

Divergence of ISD-type methods

Even though dimer/GAD type methods converge locally under reasonable hypotheses, global convergence is out of reach. We briefly summarize two examples from [START_REF] Gould | A dimer-type saddle search algorithm with preconditioning and linesearch[END_REF][START_REF] Xiang | The gentlest ascent dynamics[END_REF] to motivate our subsequent results.

One of the simplest examples is the 1D double-well [6]

E(x) = (1 -x 2 ) 2 . ( 5
)
On this one-dimensional landscape, the ISD (3) is the gradient ascent dynamics. It converges to the saddle at x = 0 if and only if started with |x 0 | < 1. If started from |x 0 | > 1, it will diverge to ±∞. This possible divergence is usually accounted for in practice by starting the method with a random perturbation from a minimum. Here, this means that the method will converge 50% of the time.

A natural extension, studied in [START_REF] Xiang | The gentlest ascent dynamics[END_REF], is the 2D double well

E(x, y) = (1 -x 2 ) 2 + αy 2 , ( 6
) 2 1 0 1 2 2 1 0 1 2 (a) α = 2 2 1 0 1 2 2 1 0 1 2 (b) α = 6 Figure 1: E(x, y) = (1 -x 2 ) 2 + αy 2 .
This energy landscape contains a saddle at (0, 0), two minima at x = ±1, y = 0 (black dots), and singularities at x = ±r c (red line). Arrows indicate the direction of the ISD. The shaded region is the index-1 region where λ 1 < 0 < λ 2 .

where α > 0, which has a saddle at (0, 0) and minima at (±1, 0). At any (x, y) ∈ R 2 ,

∇ 2 E(x, y) = 4(3x 2 -1) 0 0 2α .
At x = ±r c , with r c = 2+α 6 , ∇ 2 E(x, y) has equal eigenvalues. As x crosses ±r c , v 1 (x) jumps: for |x| < r c , v = ±(1, 0), while for |x| > r c , v = ±(0, 1).

The lines {x = ±r c } are a singular set for the ISD while, for |x| = r c the ISD is given by

ẋ ẏ = σ(x) 4x(x 2 -1) -2αy where σ(x) = 1, |x| < r c , -1, |x| > r c .
As x approaches ±r c , ẋ approaches ±-4r c (r 2 c -1). The resulting behavior of the system depends on whether r c is greater or less than 1. For r c > 1 (α > 4), the singular line is attractive, while for r c < 1 (α < 4), the line is repulsive. When the singular line is attractive, the solution of the ISD stops existing in finite time (an instance of blowup). The resulting phase portraits is shown in Figure 1. Note that, for α < 4, every trajectory started in a neighborhood of the minima diverges. For α > 4, trajectories started from a random perturbation of a minimum converge 50% of the time.

This example shows the importance of singularities for the ISD. The GAD, due to the lag in the evolution of v, does not adapt instantaneously to the discontinuity of the first eigenvector. Instead one expects that it will oscillate back and forth near a singularity, at least for ε sufficiently small.

Neither of the two examples we discussed here is generic: in the 1D example (5) both ISD and GAD reduce to gradient ascent, while in the 2D example [START_REF] Gould | A dimer-type saddle search algorithm with preconditioning and linesearch[END_REF] the set of singularities is a line, whereas we expect point singularities; we will discuss this in detail in § 3.1.

New Results: basin of attraction

The basin of convergence of the saddle {(x, y), |x| ≤ min(1, r c )} for ( 6) is fairly large and in particular includes the index-1 region {(x, y), |x| ≤ 1 √ 3 } where the first two eigenvalues λ 1 (x) Figure 2: E(x, y) = (x 2 + y 2 ) 2 + x 2 -y 2 -x + y, a coercive energy functional with a minimum at M , an attractive singularity at S 1 and a repulsive singularity at S 2 . The shaded area is the index-1 region. The arrows indicate the direction of the ISD. The black lines are contour lines of E. The ISD remains trapped in this energy well, and trajectories converge to S 1 . This example shows, in particular, that an index-1 region is insufficient to guarantee even the existence of a saddle. and λ 2 (x) of ∇ 2 E(x) satisfy λ 1 (x) < 0 < λ 2 (x). This and other examples motivate the intuition that, when started in such an index-1 region, the ISD and GAD will converge to a saddle.

Our results in Section 2 formalizes this intuition but with an added assumption: we prove in Theorem 2 that the ISD converges to a saddle if it is started in a an index-1 region Ω that is a connected component of a sublevel set for ∇E . In Theorem 2 the same result is proven for the GAD, under the additional requirement that v(0) -v 1 (x(0)) and ε are sufficiently small.

These results give some credence to the importance of index-1 regions, but only guarantee convergence under a (strong) additional hypothesis. We show in Figure 2 an index-1 region with no saddles inside, demonstrating the importance of this additional hypothesis.

New Results: singularities and quasi-periodic orbits

Let E ∈ C 2 (R N ) and, for x ∈ R N , let λ 1 (x) ≤ λ 2 (x) denote the two first eigenvalues of ∇ 2 E(x)
and v 1 (x), v 2 (x) the associated eigenvectors. The set at which eigenvalues cross is the set of singularities

S := x ∈ R N | λ 1 (x) = λ 2 (x) .
Note that v 1 (x) is well-defined only for x ∈ R N \ S. Accordingly, the ISD is defined only away from S.

In Section 3, we study the local structure of singularities in 2D. We first show that, unlike in § 1.3, singularities are generically isolated, and stable with respect to perturbations of the energy functional. We then examine the ISD around isolated singularities, in particular classifying attractive singularities such as S 1 in Figure 2, which give rise to finite-time blow-up of the ISD.

For such attractive singularities, the GAD does not have time to adapt to the rapid fluctuations of v 1 (x) and oscillates around the singularity. For ε small, we prove in special cases that the resulting behavior for the GAD is a stable annulus of radius O(ε) and of width O(ε 2 ) around the singularity. We call such a behavior "quasi-periodic". Our main result is Theorem 8, which generalizes this to the multi-dimensional setting and proves stability with respect to arbitrary small perturbations of the energy functional E.

Notation

We call N ≥ 1 the dimension of the ambient space, and (e i ) 1≤i≤N the vectors of the canonical basis. For a matrix M , we write M op = sup x∈S 1 M x its operator norm, where S 1 denotes the unit sphere in R N . In our notation, I is the identity matrix and scalars may be interpreted as matrices. Matrix inequalities are to be understood in the sense of symmetric matrices: thus, for instance, when λ ∈ R, M ≥ λ and M -λ ≥ 0 both mean that x, M x ≥ λ x 2 for all x ∈ R N . When A is a third-order tensor and u, v, w ∈ R N , we write A[u] for the contracted matrix (A[u]) ij = N k=1 A ijk u k , and similarly A[u, v] and A[u, v, w] for the contracted vector and scalar.

E will always denote an energy functional defined on R N . We will write ∇ k E(x) for the k-th-order tensor of derivatives at x. λ i (x) and v i (x) refer to the i-th eigenvalue and eigenvector (whenever this makes sense) of ∇ 2 E(x).

A 2 × 2 matrix representing a rotation of angle ω will be denoted by R ω .

2 Region of attraction

Idealized dynamics

We first consider the ISD (3), and prove local convergence around non-degenerate index-1 saddles.

Lemma 1. (a) Let E ∈ C 3 (R N ) and λ 1 (x * ) < λ 2 (x * ) for some x * ∈ R N , then F ISD (x) := -(I -2v 1 (x) ⊗ v 1 (x))∇E(x) is C 1 in a neighborhood of x * . (b) If x * ∈ R N is an index-1 saddle, then ∇F ISD (x *
) is symmetric and negative definite. In particular, x * is exponentially stable under the ISD (3).

Proof. The proof of (a) follows from a straightforward perturbation argument for the spectral decomposition, given the spectral gap λ 1 < λ 2 . As part of this proof one obtains that

x → v 1 (x) ∈ C 1 .
To prove (b), we observe that, since ∇E(x * ) = 0,

∇F ISD (x * )[h] = -(I -2v 1 ⊗ v 1 )∇ 2 E(x * )[h] + 2∇v 1 [h] v 1 , ∇E(x * ) + 2(v 1 ⊗ ∇v 1 [h])∇E(x * ) = -(I -2v 1 ⊗ v 1 )∇ 2 E(x * )[h].
Therefore, ∇F ISD (x * ), is symmetric and negative definite, which implies the result.

Next we give an improved estimate on the ISD region of attraction of an index-1 saddle.

Theorem 2. Let E ∈ C 3 (R N ), L > 0 a level and let Ω ⊂ R N be a closed connected component of {x ∈ R N | ∇E(x)
≤ L} which is bounded (and therefore compact). Suppose, further, that

λ 1 (x) < 0 < λ 2 (x) for all x ∈ Ω.
Then, for all x 0 ∈ Ω, the ISD (3) with initial condition x(0) = x 0 admits a unique global solution x ∈ C 1 ([0, ∞); Ω). Moreover, there exist an index-1 saddle x * ∈ Ω and constants

K, c > 0 such that x(t) -x * ≤ Ke -ct .
Proof. The result is based on the observation that, if x ∈ C 1 ([0, T ]) solves the ISD (3), then for 0 < t < T ,

d dt ∇E(x) 2 = 2 d dt ∇E(x), ∇E(x) = 2 ∇ 2 E(x) ẋ, ∇E(x) = -2 ∇E(x), ∇ 2 E(x)(I -2v ⊗ v)∇E(x) ≤ -2 min(-λ 1 , λ 2 ) ∇E(x) 2 .
It follows that Ω is a stable region for the ISD. Since Ω is bounded and

S ∩ Ω = ∅, if x(0) ∈ Ω, then (3) has a global solution x ∈ C 1 ([0, ∞); Ω). Because Ω is compact, inf x∈Ω -2 min(-λ 1 , λ 2 ) > 0.
It follows that ∇E(x(t)) → 0 with an exponential rate. Again by compactness, there exists x * ∈ Ω and a subsequence

t n ↑ ∞ such that x(t n ) → x * . Since ∇E(x(t)) → 0, we deduce ∇E(x * ) = 0. Since λ 1 (x * ) < 0 < λ 2 (x * ) it follows that x * is an index-1 saddle.
Since we have now shown that, for some t > 0, x(t) will be arbitrarily close to x * , the exponential convergence rate follows from Lemma 1.

Gentlest Ascent Dynamics

The analogue of Theorem 2 for the GAD (4) requires that the relaxation of the rotation is sufficiently fast and that the initial orientation v(0) is close to optimal. Theorem 3. Assume the same prerequisites as Theorem 2.

Then, for ε, δ > 0 sufficiently small, the GAD (4) with any initial condition x(0) = x 0 ∈ int(Ω) and v 0 ∈ S 1 such that v(0) -v 1 (x 0 ) < δ admits a unique global solution (x, v) ∈ C 1 ([0, ∞); Ω × S 1 ). Moreover, there exists an index-1 saddle x * ∈ Ω, and constants K, c > 0 such that

x(t) -x * + v(t) -v 1 (x * ) ≤ Ke -ct .
The proof of this result, which is more technical but at its core follows the same idea as Theorem 2, can be found in Appendix A. The additional ingredient is to control v(t)-v 1 (x(t)) , using smallness of ε and the separation of the eigenvalues λ 1 < 0 < λ 2 in Ω.

An example of global convergence and benchmark problem

An immediate corollary of Theorem 3 is the following result.

Corollary 4. Suppose that E ∈ C 3 (R N ) has the properties λ 1 (x) < 0 < λ 2 (x) ∀x ∈ R N , ∇E(x) → ∞ as |x| → ∞.
Then, for every r > 0 there exists ε r , δ r > 0 such that the ε-GAD (4) with ε ≤ ε r and initial conditions satisfying x(0) ≤ r and v(0) -v 1 (x(0)) ≤ δ r has a unique global solution (x(t), v(t)) which converges to an index-1 saddle.

We mention this result as it establishes a simplified yet still non-trivial situation, somewhat analogous to convex objectives in optimization, in which there is a realistic chance to develop a rigorous global convergence theory for practical saddle search methods that involve adaptive step size selection and choice of rotation accuracy. Work in this direction would generate ideas that strengthen the robustness and efficiency of existing saddle search methods more generally.

Singularities and (quasi-)periodic orbits

We now classify the singularities S for the ISD (3) in 2D, exhibit finite-time blow-up of the ISD and (quasi-)periodic solutions of the GAD.

Isolated singularities and the discriminant

Recall that the set of singularities for the ISD is denoted by

S = {x ∈ R 2 | λ 1 (x) = λ 2 (x)}. The ISD is defined on R 2 \ S.
Since symmetric matrices with repeated eigenvalues are a subset of codimension 2 of the set of symmetric matrices, one can expect that S contains isolated points. This phenomenon is sometimes known as the Von Neumann-Wigner no-crossing rule [START_REF] Neumann | Uber merkwürdige diskrete eigenwerte. uber das verhalten von eigenwerten bei adiabatischen prozessen[END_REF]. This is particularly easy to see in dimension 2, because the only 2 × 2 matrices with repeated eigenvalues are multiples of the identity, and therefore are a 1-dimensional subspace of the 3dimensional space of 2 × 2 symmetric matrices. To transfer this to the set S, we first note that a point x ∈ R 2 is a singularity if and only if

e 1 , ∇ 2 E(x)e 1 = e 2 , ∇ 2 E(x)e 2 , e 2 , ∇ 2 E(x)e 1 = 0.
Writing this system of equations in the form F (x) = 0, if the Jacobian ∇F (0) is invertible, then the singularity is isolated.

For i, j, k = 1, 2 we define

E ijk = ∇ 3 E(0)[e i , e j , e k ], ∆ = (E 111 E 122 + E 112 E 222 ) -E 2 112 + E 2 122 ,
then we can compute

∇F (0) = E 111 -E 122 E 112 -E 222 E 112 E 122 and det(∇F (0)) = ∆.
If ∆ = 0 (which we expect generically) then the singularity is isolated. By the implicit function theorem, this also implies that such a singularity is stable with respect to small perturbations of the energy functional (see Lemma 9 for more details). Note that this it not the case for the example E(x, y) = (1 -x 2 ) 2 + αy 2 of Section 1.3, which has a line of singularities on which ∆ = 0. This is due to the special form of the function, where the hessian is constant along vertical lines. This behavior is not generic, and under most perturbations the singularity set S will change to a discrete set (this statement can be proven using the transversality theorem).

Formal expansion of the ISD and GAD near a singularity

We consider the ISD and GAD dynamics in the neighborhood of a singularity situated at the origin. In the following, we assume ∆ = 0, so that the singularity is isolated.

Let λ := λ 1 (0) = λ 2 (0), then expanding E about 0 yields

∇E(x) = ∇E(0) + λx + O( x 2 ), ∇ 2 E(x) = λI + ∇ 3 E(0)[x] + O( x 2 ).
Inserting these expansions into the GAD (4) yields

ẋ = -(I -2v ⊗ v)∇E(0) + O( x ), ε 2 v = -(I -v ⊗ v)∇ 3 E(0)[x, v] + O( x 2 ),
and dropping the higher-order terms we obtain the leading-order GAD

ẋ = -(I -2v ⊗ v)∇E(0), ε 2 v = -(I -v ⊗ v)∇ 3 E(0)[x, v]. (7)
Since ∆ = 0, v 1 (x) is well-defined in B r (0) \ {0} for some r > 0. To leading order, v 1 (x) is given by

v 1 (x) = w 1 (x) + O( x ),
where w 1 is the eigenvector corresponding to the first eigenvalue of

∇ 3 E(0)[x]. Inserting the expansions for ∇ 2 E and v 1 into the ISD yields ẋ = -(I -2w 1 (x) ⊗ w 1 (x))∇E(0) + O( x ),
and dropping again the O( x ) term we arrive at the leading-order ISD

ẋ = -(I -2w 1 (x) ⊗ w 1 (x))∇E(0). (8)
Next, we rewrite the leading-order GAD and ISD in a more convenient format. If v = (cos φ, sin φ) ∈ S 1 , then we define v := cos(2φ), sin(2φ) .

Furthermore, we define the matrix

A := E 111 -E 122 2 E 112 -E 222 2 E 112 E 122 , (10)
which coincides with ∇F (0), up to the scaling of the first row. In particular, det A = 1 2 ∆.

Lemma 5. Suppose that ∇E(0) = (cos α, sin α), then the leading-order GAD [START_REF] Henkelman | A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives[END_REF] and ISD (8) are, respectively, given by

ẋ = R -α v, ε 2 v = -2 R π/2 v, Ax R π/2 v ( 11
)
and ẋ = - R -α Ax Ax . (12)
Proof. For y = (a, b) ∈ R 2 , we define the matrix

Q y = a b b -a .
Geometrically, if y = 1, then Q y describes a reflection with respect to the line whose directing angle is half that of y. Accordingly, for v ∈ S 1 ,

I -2v ⊗ v = -Q v, (13)
and hence the evolution of x in the leading-order GAD equation ( 7) reduces to

ẋ = Q v R α e 1 = R -α Q v e 1 = R -α v,
which establishes the first equation in [START_REF] Neumann | Uber merkwürdige diskrete eigenwerte. uber das verhalten von eigenwerten bei adiabatischen prozessen[END_REF].

To derive the second equation in [START_REF] Neumann | Uber merkwürdige diskrete eigenwerte. uber das verhalten von eigenwerten bei adiabatischen prozessen[END_REF], subtracting the average of the diagonal entries of

∇ 3 E(0)[x] yields ∇ 3 E(0)[x] = E 111 x 1 + E 112 x 2 E 112 x 1 + E 122 x 2 E 112 x 1 + E 122 x 2 E 122 x 1 + E 222 x 2 = c(x)I + Q Ax , for some c(x) ∈ R. Together with I -v ⊗ v = R π/2 v ⊗ R π/2 v and R π/2 v, cIv = c R π/2 v, v = 0 this observation implies ε 2 v = -R π/2 v, cI + Q Ax ]v R π/2 v = -R π/2 v, Q Ax v R π/2 v.
Writing Ax = (a, b) and v = (cos φ, sin φ), we obtain

v, Q Ax v = -b cos 2 φ -sin 2 φ + 2a cos φ sin φ = -b cos(2φ) + a sin(2φ) = R π/2 v, Ax
and thus arrive at

ε 2 v = -R π/2 v, Ax R π/2 v.
Using the observations v = φR π/2 v and v = φ2R π/2 v immediately yields the second equation in [START_REF] Neumann | Uber merkwürdige diskrete eigenwerte. uber das verhalten von eigenwerten bei adiabatischen prozessen[END_REF]. Finally, to obtain [START_REF] Zhang | Shrinking dimer dynamics and its applications to saddle point search[END_REF] we first observe that the stationary points for the v equation of unit norm are w± = ±Ax/ Ax , the stable one being w-= -Ax/ Ax (corresponding to v being the lowest eigenvector of Q Ax and therefore ∇ 3 E(0)[x]). Applying (13) we obtain

ẋ = Q w-R α e 1 = R -α w-= - R -α Ax Ax .

Finite-time blow-up of the ISD near singularities

We assume, without loss of generality, that ∇E(0) = 1. Then it follows from Lemma 5 that the ISD is given by

ẋ = - R -α Ax Ax + O( x ),
where A is given by [START_REF] Olsen | Comparison of methods for finding saddle points without knowledge of the final states[END_REF]. Thus, starting sufficiently close to the origin, we can study the ISD using the tools of linear stability. Observe first that

det(R -α A) = det A = ∆ 2 .
If ∆ > 0 then R -α A either has two real eigenvalues with the same sign or it has a pair of complex conjugate eigenvalues. This results in a singularity that is either attractive, repulsive or a center (see Figure 3a, 3b and 3c respectively). If ∆ < 0, R -α A has two real eigenvalues of opposite sign and hence the origin will exhibit saddle-like behavior; cf. Figure 3d.

We are specifically interested in attractive singularities such as the one in Figure 3a. In this context, we prove the following proposition: Proposition 6. Suppose that 0 ∈ S is a singularity such that ∆ > 0 and that R -α A has two eigenvalues (counting multiplicity) with positive real part. Then, for x(0) = 0 sufficiently small the corresponding maximal solution x ∈ C 1 ([0, T * )) of (3) has blow-up time T * < ∞ and x(t) → 0 as t → T * .

(x) = cos αx 1 + sin αx 2 + 1 2 (x 2 1 + x 2 2 ) + 1 2 (sx 3 1 + x 1 x 2 2 ).
Proof. We have shown in Lemma 5 that the ISD can be written in the form

ẋ = 1 Bx -Bx + g(x) ,
where B = R -α A has eigenvalues with positive real part, g(x) ≤ C x 2 for x in a neighborhood of the origin, and C > 0 a constant. According to standard ODE theory, there is a maximal solution x(t) in an interval [0, T * ). Assuming T * = ∞, we will obtain a contradiction by showing that x(t) = 0 for some finite t. Diagonalizing B (or taking its Jordan normal form), there exists an invertible P ∈ R 2×2 such that B = P DP -1 , with D of one of the following three forms:

D = λ r -λ i λ i λ r , D = λ 1 0 0 λ 2 , or D = λ ε 0 λ ,
where ε may be chosen arbitrarily small. In all cases, by the hypothesis that B has eigenvalues

(x) = cos αx 1 + sin αx 2 + 1 2 (x 2 1 + x 2 2 ) + 1 2 (x 3 1 + x 1 x 2
2 ), α = π/4; cf. Sections 3.4 and 3.5. where

C 1 = µ R P , C 2 = C P -1 R -1 P 2 ,
for y in a neighborhood of zero. It follows that, when x(0) is sufficiently small, then y 2 is decreasing and reaches zero in finite time.

Proposition 6 demonstrates how the ISD, a seemingly ideal dynamical system to compute saddle points can be attracted into point singularities and thus gives a further example of how the global convergence of the ISD fails. Next, we examine the consequences of this result for the GAD.

The isotropic case

The GAD (11) is a nonlinear dynamical system of dimension 3 (two dimensions for x, one for v), who are known to exhibit complex (e.g. chaotic) behavior. In the setting of Proposition 6, we expect that "most" solutions of the GAD converge to a limit cycle. Numerical experiments strongly support this claim, but indicate that the limit cycles can be complex; see Figure 4.

We now seek to rigorously establish the existence of (quasi-)periodic behavior of the GAD, at least in special cases. To that end we write

A = R s DR t ,
where R s , R t ∈ SO(2) and D = diag(d 1 , d 2 ), and we recall that d 1 d 2 = det D = det A = ∆/2 > 0, that is, d 1 , d 2 have the same sign. Since R s , R t , R -α , R π/2 commute, under the substitution y = R t x, w = R -s v, the leading-order GAD becomes ẏ = R s-t-α w,

ε 2 ẇ = -2 R π/2 w, Dy R π/2 w.
From Figure 4 we observe that a complex limit cycle can occur in the anisotropic case d 1 = d 2 , while the behavior when d 1 = d 2 is much simpler. In order to get a tractable system, we restrict ourselves in the following to the isotropic case d 1 = d 2 , where we can use polar coordinates to perform a stability analysis. Note that this corresponds to imposing that A is a multiple of a rotation matrix:

A 11 = A 22 , A 12 = -A 21 .
ẋ = R -α v ε 2 v = -2 R π/2 v, x R π/2 v. ( 14
)
Thus, up to a rescaling of ε, a rotation of x and a shift in α (rotation of ∇E(0)), restricting to isotropic matrices is equivalent to restricting to A = I, which corresponds to

E 111 = 3, E 112 = 0, E 122 = 1, E 222 = 0, or E(x) = cos αx 1 + sin αx 2 + λ 2 (x 2 1 + x 2 2 ) + 1 2 (x 3 1 + x 1 x 2 2 ) + O( x 4 ). ( 15
)
We consider this case in the sequel, as well as the restriction cos α > 0, which ensures that R -α has eigenvalues with positive real part and therefore that the ISD converges to zero.

Explicit solutions of the leading-order GAD in the attractive isotropic case

We now produce an explicit solution of the leading-order isotropic GAD (14), which makes precise the intuition that delayed orientation relaxation of the GAD balances the blow-up of the ISD and thus leads to periodic orbits. On substituting polar coordinates x = r(cos θ, sin θ), v = (cos φ, sin φ) and hence v = (cos 2φ, sin 2φ) in ( 14) we obtain a set of three coupled ODEs for r, θ and φ:

ṙ = cos(2φ -α -θ) r θ = sin(2φ -α -θ) ε 2 φ = r sin(2φ -θ).
(

) 16
We now analyze the behavior of this set of equations for ε 1. This corresponds to an adiabatic limit where the evolution of v is fast enough to relax instantly to its first eigenvector, so that the dynamics mimics closely the ISD. However, this is counterbalanced by the fact that the dynamics for v becomes slow as r → 0.

The r dynamics takes place at a timescale 1, the θ dynamics at a timescale r, and the φ dynamics at a timescale ε 2 r . The adiabatic approximation of fast relaxation for v (the ISD) is valid when ε 2 r r, or r ε. In this scaling we recover the ISD [START_REF] Zhang | Shrinking dimer dynamics and its applications to saddle point search[END_REF]. One the other hand, when r ε, then θ relaxes to a stable equilibrium 2φ -α -θ = 2kπ, k ∈ Z, in which case we obtain ṙ = +1. Therefore we may expect that for r ε, r decreases, while for r ε, r increases. We now examine the intermediate scaling r ∼ ε. Rescaling r = εr we obtain

ε ṙ = cos(2φ -α -θ) εr θ = sin(2φ -α -θ) ε φ = r sin(2φ -θ).
All variables now evolve at the same characteristic timescale ε, hence we rescale t = εt . For the sake of simplicity of presentation we drop the primes to obtain the system

ṙ = cos(2φ -α -θ) r θ = sin(2φ -α -θ) φ = r sin(2φ -θ), (17)
which describes the evolution ( 16) on time and space scales of order ε.

We observe that the evolution of (17) does not depend on θ and φ individually, but only on ω = 2φ -θ. Keeping only the variables of interest, r and ω, we arrive at the 2-dimensional system ṙ = cos(ω -α)

ω = 2r sin ω - 1 r sin(ω -α). (18)
Since cos α > 0, (18) has two fixed points, with associated stability matrix J ± ,

r 0 = 1 2 cos α , ω ± 0 = α ± π 2 , J ± = 0 ∓1 ±4 cos α ∓ 2 sin α √ 2 cos α . ( 19
)
The determinant of J ± is positive. The eigenvalues are either complex conjugate or both real; in both cases their real part is of the same sign as the trace,

trJ ± = ∓ 2 sin α √ 2 cos α .
If sin α > 0, then (r 0 , ω + 0 ) is stable, whereas if sin α < 0, then (r 0 , ω - 0 ) is stable. The case sin α = 0 cannot be decided from linear stability, and so we exclude it in our analysis.

In real variables, the resulting behavior is that the system stabilizes in a periodic orbit at

r 0 = ε 1 2 cos α .
θ evolves twice at fast as φ, so that ω = 2φ -θ stays constant at ω ± = α ± π 2 . Thus we have established the following result. Lemma 7. If cos α > 0, sin α = 0, then the projection (18) of the leading order isotropic GAD admits a stable circular orbit of radius

r = ε √ 2 cos α .
In the next section, we will show that this behavior survives to a threefold generalization: the re-introduction of the neglected higher-order terms, perturbations of the energy functional, as well as dimension N > 2.

Quasi-periodic solutions of GAD

The computation of Section 3.5 suggests that the GAD for the energy functional

E(x 1 , x 2 ) = (cos αx 1 + sin αx 2 ) + λ 2 (x 2 1 + x 2 2 ) + 1 2 (x 3 1 + x 1 x 2 2 ) (20)
has nearly periodic trajectories near the origin when α ∈ (-π/2, 0) ∪ (0, π/2). Any thirdorder term of the form ax 3 1 + bx 2 1 x 2 + ax 1 x 2 2 + bx 3 2 for a, b ∈ R reduces to (20) upon a suitable change of variables. We will now rigorously prove the existence of quasi-periodic behavior in the multidimensional and perturbed case. We split an N -dimensional state space V = R N into two components V = V s ⊕ V c : a two-dimensional subspace V s (singular ) on which the dynamics is the same as in the 2D case, and an (N -2)-dimensional subspace V c (converging) on which the GAD dynamics converges to zero. Let I = {1, . . . , N }, I s = {1, 2} and I c = {3, . . . , N } be the corresponding set of indices and, for x ∈ R N , x s = (x 1 , x 2 , 0, . . . , 0) ∈ V s and x c = (0, 0, x 3 , . . . , x N ) ∈ V s .

We consider a functional E = E 0 of the form,

E 0 (x) = (cos α 0 x 1 + sin α 0 x 2 ) + λ 0 2 (x 2 1 + x 2 2 ) + 1 2 (x 3 1 + x 1 x 2 2 ) + 1 2 N i,j∈Ic H 0 ij x i x j + 1 6 i,j,k∈Ic G 0 ijk x i x j x k + O( x 4 ), (21)
where

α 0 ∈ (0, π/2)∪(π/2, π), λ 0 ∈ R, H 0 ij = ∇ 2 E 0 (0)[e i , e j ], i, j ∈ I c , and G 0 ijk = ∇ 3 E 0 (0)[e i , e j , e k ], i, j, k ∈ I.
For x c = 0, E 0 coincides with (20) to within O( x 4 ) and the condition on α 0 are consistent with Lemma 7. We assume for the remainder that H 0 > max(λ 0 , 0): the requirement H 0 > λ 0 ensures that λ 0 is indeed the lowest eigenvalue, while H 0 > 0 ensures that x c → 0 as t → ∞.

An example of a functional in this class, and the resulting GAD dynamics are shown in Figure 5. Our main result is the following theorem stating that the limit cycles at r = ε/ √ 2 cos α 0 present in the 2D leading-order GAD survive in the nonlinear, multidimensional, perturbed regime. The proof is given in Appendix B. Theorem 8. Let E 0 , ∆E ∈ C 4 (R N ) with E 0 satisfying (21) with H > max(0, λ 0 )I. For δ > 0 let E δ := E 0 +δ∆E. Then there exist constants δ 0 , ε 0 , m, M > 0 such that, for all ε < ε 0 , δ < δ 0 , the following statements hold.

There exists z

δ ∈ R N with z δ ≤ M δ such that ∇ 2 E δ (z δ ) has repeated eigenvalues λ 1 = λ 2 and ∇E δ (z δ ) ∈ span{e 1 , e 2 }
, where e i are the eigenvectors corresponding to λ i .

For all

x 0 ∈ R N such that x 0 -z δ = ε/ √ 2 cos α 0 (cf. Lemma 7) there exists v 0 ∈ S 1 such that the ε-GAD (4) for E δ with initial conditions satisfying x(0) -x 0 ≤ mε(δ + ε) and v(0) -v 0 ≤ mε(δ + ε), v(0) ∈ S 1 ,
admits a unique solution, and

x(t) -z δ - ε √ 2 cos α 0 ≤ M ε(ε + δ) for all t ≥ 0.

Conclusion

In this paper we make two novel contributions to the theory of walker-type saddle search methods: Region of attraction: In Section 2 we extended estimates on the region of attraction for an index-1 saddle beyond perturbative results. Our results give some credence to the widely held belief that dimer and GAD type saddle search methods converge if started in an index-1 region.

5: E(x) = cos αx + sin αx + 1 2 (x 2 + y 2 + 1.1z 2 ) + 1 2 (x 2 y + y 3 ) + z 3 , α = 3π/4.
But we also show through an explicit example that this is not true without some additional assumptions on the energy landscape. We also highlight the global convergence result of Corollary 4, which we believe can provide a useful benchmark problem and testing ground towards a more comprehensive convergence theory for practical saddle search methods outside a perturbative regime.

Cycling behavior: Although it is already known from [START_REF] Gould | A dimer-type saddle search algorithm with preconditioning and linesearch[END_REF][START_REF] Xiang | The gentlest ascent dynamics[END_REF] that the dimer and GAD methods cannot be expected to be globally convergent, the behavior identified in those references is non-generic. In Section 3 we classify explicitly the possible generic singularities and identify a new situation in which global convergence fails, which occurs in any dimension and is stable under arbitrary perturbations of the energy functional.

In particular, our results provide a large class of energy functionals for which there can be no merit function for which the ISD or GAD search directions are descent directions.

Our results illustrate how fundamentally different saddle search is from optimization, and strengthen the evidence that dimer or GAD type saddle search methods cannot be (easily) modified to obtain globally convergent schemes. Indeed it may even prove impossible to design a globally convergent walker-type saddle search method.

We speculate that this is related to the difficulty in proving the existence of saddle points in problems in calculus of variations: while the existence of minimizers follow in many cases from variational principles, the existence of saddle points is known to be more difficult, requiring sophisticated mathematical tools such as the mountain pass theorem. Such theorems are based on the minimization of functionals of paths, and as such are conceptually closer to string-of-state methods such as the nudged elastic band and string methods [START_REF] Jónsson | Nudged elastic band for finding minimum energy paths of transitions[END_REF][START_REF] Ren | String method for the study of rare events[END_REF]. It is to our knowledge an open question to establish the global convergence of methods from this class, but the results of the present paper suggest that it may be a more promising direction to pursue than the walker type methods.

A Proof of Theorem 3 (region of attraction for the GAD)

At a fixed configuration x, the dynamics on v is a gradient descent for v, ∇ 2 E(x)v on the sphere S 1 that ensures the local convergence of v to v 1 . Our strategy is to use an adiabatic argument to show that the full GAD dynamics keeps v close to v 1 even when x is (slowly) evolving.

Assume Ω is as in the hypotheses of this theorem. Then Ω is compact, and the minimal spectral gap g of the positive-definite continuous matrix (1

-2v 1 (x) ⊗ v 1 (x))∇ 2 E(x) satisfies g = min x∈Ω min(-λ 1 (x), λ 2 (x)) > 0.
This implies that λ 2 -λ 1 ≥ 2g > 0 on Ω. Let also

M = max x∈Ω ∇ 2 E(x) op .
For v ∈ S 1 , we write P v = v ⊗ v. When v, w ∈ S 1 , then we have the following improved Cauchy-Schwarz equality

v, w = 1 - 1 2 v -w 2
and bound on projectors

P v -P w op = (v -w) ⊗ (v + w) op ≤ 2 v -w .
Step 1: variations of ∇E . For any t such that x(t) ∈ Ω, we compute (dropping the dependence on x(t), and writing

H = ∇ 2 E) d dt 1 2 ∇E 2 = -(I -2P v )H∇E, ∇E = -(I -2P v 1 )H∇E, ∇E + (P v -P v 1 )H∇E, ∇E ≤ (-g + 2M v -v 1 ) ∇E 2 . ( 22
)
Step 2: variations of v -v 1 . Similarly, when x ∈ Ω, we compute

d dt 1 2 v -v 1 2 = v -v 1 , v -v1 = - 1 ε 2 v -v 1 , (I -P v )Hv -v -v 1 , v1 .
Our goal is (25) below, which shows that the leading term in this expression is bounded by -ε -2 g v -v 1 2 , which will pull back v to v 1 when ε is small enough. We bound both terms separately. For the first term, we note that

(I -P v )Hv = (I -P v 1 )H(v -v 1 ) + (P v 1 -P v )Hv 1 + (P v 1 -P v )H(v -v 1 ) = (I -P v 1 )H(v -v 1 ) + λ 1 (v 1 -v, v 1 v) + (P v 1 -P v )H(v -v 1 ) = (I -P v 1 )H(v -v 1 ) -λ 1 v -v 1 - 1 2 v -v 1 2 v + (P v 1 -P v )H(v -v 1 ), hence it follows that -ε -2 v -v 1 , (I -P v )Hv ≤ -ε -2 v -v 1 , ((I -P v 1 )H -λ 1 I)(v -v 1 ) + ε -2 1 2 λ 1 + 2M v -v 1 3 ≤ -ε -2 g v -v 1 2 + ε -2 1 2 λ 1 + 2M v -v 1 3 . ( 23
)
For the second term, standard eigenvector perturbation theory yields

v1 = -(H -λ 1 ) + Ḣv 1 = -(H -λ 1 ) + ∇ 3 E(x)[ ẋ, v 1 ],
where (H -λ 1 ) + is the Moore-Penrose pseudo-inverse of H -λ 1 , defined by

(H -λ 1 ) + v 1 = 0 and (H -λ 1 ) + v i = 1 λ i -λ 1 v i , for i > 1. It follows that v1 ≤ g -1 ∇ 3 E(x)[ ẋ, v 1 ] and then, from ẋ ≤ L, v -v 1 , v1 ≤ L g max x∈Ω ∇ 3 E(x)[v 1 (x)] op v -v 1 . (24)
Estimates ( 23) and ( 24) imply the existence of constants C 1 , C 2 > 0 such that, when x ∈ Ω,

d dt 1 2 v -v 1 2 ≤ 1 ε 2 -g + C 1 v -v 1 v -v 1 2 + C 2 v -v 1 . ( 25
)
Step 3: conclusion. Let

δ 0 = 1 2 min g C 1 , g 2M and ε 0 = gδ 0 4C 2 .
Then, for x ∈ Ω and v -v 1 ≤ δ 0 , (22) implies that ∇E is decreasing. If, in addition, ε < ε 0 and δ 0 /2 < v -v 1 , then, (25) implies that v -v 1 is decreasing as well. Let (x, v) be the maximal solution of the GAD equations on an interval [0, T c ) with initial conditions as in the Theorem, and let

T 0 = inf{t ∈ [0, T c ), x(t) ∈ Ω or v(t) -v 1 (x(t)) > δ 0 } > 0.
Assume T 0 = T c . Then, at T 0 , either ∇E(x(T 0 )) ≥ L, in contradiction with (22), or v(T 0 ) -v 1 (x(T 0 )) ≥ δ 0 , in contradiction with (25). We can conclude, in particular, that x(t) ∈ Ω for all time. Since Ω is bounded there cannot be blow-up in finite time, hence T c = +∞.

Since ∇E(x(t)) → 0, and x(t) is bounded, a subsequence converges to a critical point x * ∈ Ω which must be an index-1 saddle. Since index-1 saddles are locally attractive for the GAD (see [START_REF] Gould | A dimer-type saddle search algorithm with preconditioning and linesearch[END_REF][START_REF] Zhang | Shrinking dimer dynamics and its applications to saddle point search[END_REF] for proofs), the exponential convergence rate follows.

B Proof of Theorem 8 (quasi-periodic solutions) B.1 Perturbation of the energy functional

We prove part 1 of Theorem 8. Heuristically, the statement is true since imposing a zero gradient on V c imposes N -2 constraints, while imposing equal eigenvalues on V s imposes 2 constraints; cf.

§ 3.1 where we showed that singularities are generically isolated in 2D. By varying the location of the singularity (N degrees of freedom) and adapting the system of coordinates, we can put the perturbed energy functional in the same functional form as E 0 , except for a perturbation of α, λ, H and of the third-order coefficients G. The latter introduces an O(δ) coupling at third order between the subspaces V s and V c . Making this precise is the content of the following lemma, which also establishes the first assertion of Theorem 8.

For the remainder of this section let (e 0 i ) i∈I be the canonical basis vectors of R N , and z 0 = 0 ∈ R N the location of the singularity with δ = 0. Lemma 9 (Perturbation of singularity). Under the conditions of Theorem 8 there exists δ 0 > 0, C > 0 such that, for every δ < δ 0 , there exist α δ , λ δ , H δ , G δ , z δ and a new orthonormal basis (e δ i ) i∈I such that, with x = z δ + N i=1 x i e δ i ,

E δ (x) = ∇E δ (z) (cos α δ x 1 + sin α δ x 2 ) + λ δ 2 (x 2 1 + x 2 2 ) + 1 2 N i,j∈Ic H δ ij x i x j + i,j,k∈I G δ ijk x i x j x k + O( x 4 ),
and moreover,

max i,j,k∈I z δ -z 0 , |e δ i -e 0 i |, |α δ -α 0 |, |λ δ -λ 0 |, |H δ ij -H 0 ij |, |G δ ijk -G 0 ijk |, | ∇E δ (z) -1| ≤ Cδ.
Proof. We need to determine a new origin z and a new orthogonal basis (e i ) i∈I that are O(δ)close to z 0 and e 0 i , such that ∇E δ (z) ∈ span{e 1 , e 2 }, and e 1 and e 2 are eigenvectors of ∇ 2 E δ (z) associated with equal (smallest) eigenvalue.

Step 1: construction of the (e i ) i∈I . Let R be the distance between λ 0 and the next-lowest eigenvalue in the spectrum of ∇ 2 E 0 (0). Let γ be the circular contour in the complex plane centered on λ 0 and of radius R/2. For any z, δ small enough, with overlap matrix O ij = ẽi , ẽj . For z, δ sufficiently small, ẽi are well-defined and O ij is positive definite, hence we can define

P (z, δ) = - 1 2πi γ (∇ 2 E δ (z) -y) -1 dy
e i = N j=1 (O -1/2 ) ij ẽj .
One can readily check that (e i ) i=1,...,N is an orthonormal basis, of class C 2 with respect to z, and that the basis vectors satisfy e i -e 0 i ≤ C(δ + z ), provided that δ, z are sufficiently small. Moreover, since O ij = 0 for i ∈ I s , j ∈ I c , we have that e 1 , e 2 ∈ Ran(P (z, δ)) and therefore (e 1 , e 2 ) are a basis of Ran(P (z, δ)).

Differentiating e i , e j = δ ij with respect to z, we obtain e i , ∇ z e j + e j , ∇ z e i = 0 for all i, j ∈ I.

Step 2: construction of z. We seek z ∈ R N , near z 0 = 0, satisfying the N equations

e 1 , ∇ 2 E δ (z)e 2 = 0, (27)
e 1 , ∇ 2 E δ (z)e 1 -e 2 , ∇ 2 E δ (z)e 2 = 0, ( 28
) e i , ∇E δ (z) = 0 for i ∈ I c . (29)
Equation ( 27) combined with ∇ 2 E δ (z)e i ∈ Ran(P (z, δ)) ⊥ e j for i ∈ I s , j ∈ I c ensures that e 1 , e 2 are eigenvectors of ∇ 2 E δ (z), and equation (28) ensures that the two associated eigenvalues are the same. We write this set of equations as F (z, δ) = 0. F is a C 2 map from a neighborhood of the origin of R N × R to R N , with F (0, 0) = 0. From (26) we obtain that the Jacobian with respect to z of this system of N equations at (z, δ) = (0, 0), in the basis (e 0 1 , e 0 2 , . . . , e 0 N ), is

∂F ∂z (0, 0) =    G 0 111 -G 0 122 G 0 112 -G 0 222 0 G 0 112 G 0 122 0 0 0 H    ,
We therefore obtain that

det ∂F ∂z (0, 0) = ∆ det H, with ∆ = (G 0 111 G 0 122 + G 0 112 G 0 222 ) -(G 0 112 ) 2 + (G 0 122 ) 2 = 2.
Since we assumed that H is positive definite, it follows that ∂F ∂z (0, 0) is invertible. From the implicit function theorem, for any δ small enough, there exists z δ in an O(δ) neighborhood of 0 satisfying F (z, δ) = 0, and the result follows.

B.2 The GAD dynamics

We are now ready to prove the second assertion of Theorem 8.

Step 1: decoupling of the singular and converging dynamics. We use Lemma 9 to change variables

x = z δ + N i=1 x i e δ i , v = N i=1
v i e δ i and then drop the primes and dependence on δ for the sake of convenience of notation. For δ small enough, we set

E(x) = E δ z + N i=1 x i e i = ∇E(0) (cos αx 1 + sin αx 2 ) + λ 2 (x 2 1 + x 2 2 ) + 1 2 i,j∈Ic H ij x i x j + i,j,k∈I G ijk x i x j x k + O( x 4 ),
with ∇E(0) = 0, sin α > 0, H > λ and H > 0. We decompose x = x s + x c , and similarly v = v s + v c . We call P s and P c the associated projectors onto the spaces V s , V c We expand the GAD equations (4) to leading order in x,

ẋs = -(1 -2v s ⊗ v s )P s ∇E(0) + O( x ) ε 2 vs = -λ -λ v s 2 -v c , Hv c v s + P s G[x, v] -G[x, v, v]v s + O( x 2 ) ẋc = -(1 -2v c ⊗ v c )Hx c + 2 v s , ∇E(0) v c + 2λ v s , x s v c + O( x 2 ) ε 2 vc = -H -λ v s 2 -v c , Hv c v c + P c G[x, v] -G[x, v, v]v c + O( x 2 ).
From the 2D case, we guess the re-scaling x = εx , t = εt . Further, since we expect v c to be small, it is convenient to rescale it as well by v c = εv c . For convenience, we drop the primes again in the following equations, and we obtain

In these equations and in what follows, the notation O is understood for with a uniform constant, as long as x and v s remain bounded: a term f (x s , v s , x c , v s ) is O(ε n ) if for every R > 0, there is K > 0 such that, when x ≤ R, v s ≤ R, then |f (x s , v s , x c , v s )| ≤ Kε n . Because in (33) P c G[x, v s ] = O(δ), we expect that the restoring force of the -(H -λ)v c term will force v c to be O(ε + δ). In turn, this will make the v s , ∇E(0) v c term in (32) to be O(ε + δ), and the restoring force of the -Hx c term will make x c to be O(ε + δ). This will decouple the dynamics on V s from that on V c : expanding for x c small, we get We now study these two equations separately, using the computations of Section 3 in the 2D case.

Step 2 : linearization of the singular dynamics. We pass to angular coordinates as in the 2D case: x s = r(cos θ, sin θ), v s = v s (cos φ, sin φ). Noting that ∇E(0) = 1 + O(δ), the ẋs and vs equations become We choose the stable solution ω 0 ∈ {ω ± } with associated Jacobian J ∈ {J ± }, and linearize about the corresponding X 0 = (r 0 , ω 0 ). Denoting X = (r -r 0 , ω -ω 0 ), we obtain

Ẋ = JX + O(ε + δ + X 2 + x c ) φ = -r cos(2φ -θ) + O(ε + δ + x c ) ( 34
)
where J is negative definite.

Step 3 : stability.

Let Ω = (X, φ, x c , v c ) | X ≤ √ ε + δ, φ ∈ R, x c ≤ 1, v c ≤ 1 .
From (34), and from the ẋc , vc equations (32) and (33), writing out fully the remainder terms as f X , f φ , f xc and f vc , we obtain the system (for ε and δ sufficiently small) Ẋ = JX + f X (X, φ, x c , v c ) φ = -r cos(2φ -θ) + f φ (X, φ, x c , v c )

1 ε ẋc = -Hx c + f xc (X, φ, x c , v c ) ε vc = -(H -λ)v c + f vc (X, φ, x c , v c ) (35)
where f X , f φ , f xc and f vc are C 1 functions satisfying

|f X (X, φ, x c , v c )| ≤ C f 2 (ε + δ + X 2 + x c ) ≤ C f (ε + δ + x c ) |f φ (X, φ, x c , v c )| ≤ C f (ε + δ + x c ) |f xc (X, φ, x c , v c )| ≤ C f (ε + δ + v c ) |f vc (X, φ, x c , v c )| ≤ C f (ε + δ)
when (X, φ, x c , v c ) ∈ Ω for some C f > 0.

Our assumptions on the initial data entail that X(0) ≤ ε + δ, φ ∈ R, x c (0) ≤ ε + δ, and v c (0) ≤ ε + δ.

Let (X, φ, x c , v c ) be a maximal solution in [0, T c ). Let also

T Ω = sup{T ∈ [0, T c ), (X, φ, x c , v c ) ∈ Ω}.
Since H > λ, e -(H-λ)t ≤ Ce -ct for some C > 0, c > 0. Thus, using Duhamel's formula for the v c equation we obtain for all t ∈ [0, T Ω ] that

v c (t) = e -H-λ ε t v c (0) + 1 ε t 0 e -H-λ ε (t-t ) f vc (X(t ), φ(t ), x c (t ), v c (t ))dt , v c (t) ≤ Ce -c ε t v c (0) + CC f (ε + δ) ε t 0 e -c ε (t-t ) dt ≤ C(ε + δ) + CC f c (ε + δ).
This shows that v c (t) ≤ K(ε + δ) for all t ∈ [0, T Ω ). Analogously, applying Duhamel's formula to the x c equation, using H > 0I, we obtain that x c (t) ≤ K (ε + δ).

Applying Duhamel's formula a third time, to the X equation, and using J < 0I, we obtain X ≤ K (ε + δ). This shows that, for ε, δ small enough, T Ω = T c and therefore T Ω = T c = +∞.

We have therefore shown that, whenever X(0) ≤ ε + δ, φ(0) ∈ R, x c (0) ≤ ε + δ, v c (0) ≤ ε + δ, then there exists a unique global solution to (35) and that X(t) ≤ K (ε + δ), x c (t) ≤ K (ε + δ), v c (t) ≤ K(ε + δ) for all t ∈ R + .

Returning to the original variables and inverting the rescaling x = z + ε N i=1 x i e i , v c = ε N i=1 (v c ) i e i v = i=1,2 (v s ) i e i + i>2 (v c ) i e i completes the proof.

4 Figure 3 :

 43 Figure 3: ISD phase planes near a singularity at 0. The energy functional is E(x) = cos αx 1 + sin αx 2 + 1 2 (x 2 1 + x 2 2 ) + 1 2 (sx 3 1 + x 1 x 2 2 ).

  ISD and GAD for E

  Leading order GAD (11) in a non-isotropic case: α = -1, A = -0.3 0.4 -0.4 0.3 , which has singular values 0.7 and 0.1. The color corresponds to the phase of v.

Figure 4 :

 4 Figure 4: Leading order ISD and GAD in an isotropic (left) and anisotropic (right) case.

Figure

  Figure 5: E(x) = cos αx + sin αx + 1 2 (x 2 + y 2 + 1.1z 2 ) + 1 2 (x 2 y + y 3 ) + z 3 , α = 3π/4.

is a projector of rank 2 and C 2

 2 with respect to both z and δ. P (z, δ) projects onto the eigenspaces of ∇ 2 E δ (z) associated to the (at most two) eigenvalues in [λ 0 -R/2, λ 0 + R/2]. Next, we define ẽi (z, δ) = P (z, δ)e 0 i if i ∈ I s , (I -P (z, δ))e 0 i if i ∈ I c .

ẋs = -( 1 -

 1 2v s ⊗ v s )∇E(0) + O(ε) (30) vs = -(I -v s ⊗ v s )P s G(x, v s ) + O(ε) (31) ẋc = -εHx c + 2ε v s , ∇E(0) v c + O(ε 2 ) (32) ε vc = -(H -λ)v c -P c G[x, v s ] + O(ε)

  ẋs = -(1 -2v s ⊗ v s )P s ∇E(0) + O(ε) vs = -(I -v s ⊗ v s )P s G[x s , v s ] + O(ε + x c ).

  ṙ = cos(2φ -α -θ) + O(ε + δ) r θ = sin(2φ -α -θ) + O(ε + δ) φ = r sin(2φ -θ) + O(ε + δ + x c )As in the 2D case, we introduce ω = 2φ -θ,

  This is equivalent to the condition E 111 = 3E 122 , E 222 = 3E 112 , i.e. the cubic terms are of the form ax 3 1 + bx 2 1 x 2 + ax 1 x 2 2 + bx 3 2 , for any a, b ∈ R.Under this hypothesis, A = dR t for some scalars d > 0, t ∈ R.

				Under the transformations
	x	R t x, ε 2	ε 2 /d, α	α -t, the leading order GAD equations (11) become

CO was supported by ERC Starting Grant 335120.


C W Lim

X S Xu

Keywords: Symplectic Elasticity: Theory and Applications eigenvalue, eigenvector, energy, Hamiltonian, symplectic elasticity, symplecticity

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Symplecticity is a mathematical concept of geometry and it is an analog of "complex" in Greek first due to Weyl ͓1,2͔. Symplectic geometry is a branch of differential geometry and differential topology. Hence, it is also called symplectic topology although, in reality, the latter is only an important subset of the former. Symplectic geometry was first developed due to its close relation to the Hamiltonian system in classical mechanics where the phase space can be constructed using symplectic manifold. The theory on symplectic geometry can be referred to Koszul and Zou ͓3͔. The symplectic group is a mathematical group, previously known as the line complex group, and a symplectic group of dimension 2n is the group of 2n ϫ 2n matrices, which preserve a symplectic form where the latter can be constructed into a canonical form using symplectic basis ͓1,2͔. Besides theoretical and classical mechanics, symplectic group has been applied in a number of fields in physics and mathematics for many years particularly in relativity and gravitation ͓4͔, quantum mechanics ͓5,6͔, and the Yang-Mills field theory ͓7͔.

2 Symplectic Numerical Analysis, Symplectic Elasticity, and Other Cross Disciplinary Areas 2.1 Symplectic Numerical Methods. Symplectic approach in applied mechanics was first developed due to the quest in accurate solutions for some nonlinear dynamical systems and stable numerical algorithm in computational mechanics constructed using Hamiltonian systems. The application of symplectic space and development of numerical methods for such systems was pioneered by Feng and his research group ͓8-34͔. Feng first proposed symplectic algorithms based on symplectic geometry for Hamiltonian systems with finite and infinite dimensions, and on dynamical systems with Lie algebraic structures, such as contact systems and source free systems via the corresponding geometry and Lie group. These algorithms are superior to conventional algorithms in many practical applications, such as celestial mechanics and molecular dynamics. A collection of his works was published later in recognition of this contribution in applied mathematics and applied mechanics ͓32,33͔ in which Vol. II ͓33͔ are the papers related to symplectic algorithm and computation of Hamiltonian systems. The contribution of Feng in the development of symplectic numerical methods was particularly significant and important as highlighted in a memorial article dedicated to him by Lax ͓35͔.

Symplectic Elasticity in Applied Mechanics.

In elasticity and applied mechanics, the development of symplectic space for analytical, engineering mechanics, and applied elasticity was attributed mainly to the research group led by Zhong ͓36-43͔. A number of research texts have been published by this group of researchers and. In particular, Zhong et al. ͓36͔ published some early works on structural mechanics and control optimization; Zhong ͓37͔ introduced a new symplectic approach for systematically deriving analytical solution in engineering elasticity; Zhong ͓38-40͔ introduced the duality system and symplectic analytical solutions approach in applied mechanics; Zhong et al. ͓39,41͔ applied the symplectic methodology for control systems and problems; Yao et al. ͓42͔ systematically described the principle of symplectic elasticity and its application in applied engineering mechanics including Timoshenko beam, plate bending, and laminated plate analysis; and Zhong ͓43͔ recently published a very concise yet simple masterpiece of symplectic elasticity by putting forward an innovative way the seemingly mathematical and abstract concept of symplecticity and relating the common analytical engineering and/or mathematical tools such as state-space method and matrix multiplication using very simple single-degree-offreedom and multiple-degree-of-freedom spring-mass systems. Most of the research works of symplectic elasticity were originally pioneered by Zhong, Xu, Zhang, and Yao from Dalian University of Technology, PRC, and subsequently researched by Lim and Leung from City University of Hong Kong in collaboration with the former group of researchers.

In general, the symplectic elasticity methodology aims to develop a new analytical elasticity approach for deriving exact analytical and asymptotic solutions to some basic problems in solid mechanics and elasticity, as described in Yao et al. ͓42͔. Some of these problems have long been bottlenecks in the development of elasticity, e.g., bending of plates with no opposite sides simply supported and bending of corner-supported plates. This symplectic elasticity approach is based on the use of the Hamiltonian principle with Legendre's transformation and subsequently analytical solutions could be obtained by expansion of eigenfunctions. It is a rational and systematic methodology with a clearly defined, step-by-step derivation procedure in which no trial functions are necessary. The advantage of this symplectic approach with respect to the classical approach, which relies on the semiinverse method, is at least threefold. First, the symplectic approach alters the classical practice and concept of solution methodology and many basic problems previously unsolvable or too complicated to be solved can be resolved accordingly in a systematic and rational procedure. For instance, the conventional approach in plate and shell theories by Timoshenko ͓44͔ has been based on the semi-inverse method with trial 1D or 2D displacement functions, such as Navier's method ͓45͔ and Lévy's method ͓46͔ for plates. The trial functions, however, do not always exist except in some very special cases of boundary conditions such as plates with two opposite sides simply supported ͓47-50͔. Using the symplectic approach, trial functions are not required. Second, it consolidates the many seemingly scattered and unrelated solutions of rigid body movement and elastic deformation by mapping with a series of zero and nonzero eigenvalues and their associated eigenvectors. Last but not least, the Saint-Venant problems for plane elasticity and elastic cylinders can be described in a new system of equations and solved. The difficulty of satisfying end boundary conditions in conventional problems, which could only be covered using the Saint-Venant principle can also be resolved.

Symplecticity for Elasticity.

Many of the early works on symplectic elasticity were published in Chinese and as a result, the early works have been unavailable and unknown to researchers worldwide. The early works on analytic symplectic elasticity begins with Zhong and co-workers ͓51-57͔ and Yao ͓58͔ who presented the plane elasticity problems based on Hamiltonian systems ͓51,53͔, the reciprocal theorem, and adjoint symplectic orthogonality ͓48͔; the Saint-Venant problems ͓54,55,58͔; and symplectic eigensolutions for cylinders and body of revolution ͓56,57͔. Recently, Zhao and Chen ͓59͔ reported some numerical studies for elasticity problems. The symplectic method with Hamiltonian state-space approach was also applied to solve analytical solutions for elastic circular solids including a circular elastic cylinder with self-weight ͓60͔, a transversely isotropic cylindrical body ͓61͔ and a circular cantilever ͓62͔.

Symplectic Methodology for Beams.

Based on the foundations of symplectic methods for plane elasticity, duality in Hamiltonian systems, symplectic orthogonality, and the solutions to Saint-Venant problems, the approach is applied to applied mechanics problems for beams. The early works was reported by Leung et al. for

Symplectic Elasticity for Plates.

Due to the practical ap-plication and a quest for exact solutions, the symplectic approach for plate analysis in applied engineering mechanics has received very much attention ͓72-101͔. Most of the research articles were based on the thin plate model ͓72-74,77,79-82,85-93,95-101͔ with a few others, which worked on thick plate problems ͓75,76,78,83,84,94͔. Among the papers, many investigated rectangular plates ͓72-79,81-84,86-91,94,95,97-101͔, while only a few focused on circular or sector plates ͓80,85,92,93,96͔, and the majority analyzed isotropic plates ͓75,76,78-80,83-101͔, while a few worked on composite or laminated plates ͓72-74,77,81,82͔. The papers can be further classified into a few categories according to the specific plate problems that were analyzed. Among the specific topics, plate bending ͓78-81, 83-86,88,90-92,94,97-99,101͔ is the most popular subject of interest; others include vibration ͓56,76,87,89,95,100͔, buckling ͓93,96͔, wave propagation ͓87͔, and the Sain-Venant problem ͓74,77,80͔.

2.6 Symplectic Elasticity for Shells. Besides plate problems, the symplectic methods have also been applied for shell problems ͓102-108͔. Due to the analytical complexity, only closed cylindrical shells have been considered because for such shells the circumferential deflection can be approximated with sinusoidal shape functions, while the other ends could be arbitrarily supported or free. Among these papers, the dynamic buckling problems for cylindrical shells subjected to axial impact were investigated by Xu and co-workers ͓102,103,105,106͔. They also studied viscoelastic hollow circular cylinders ͓104͔, torsional buckling ͓107͔, and thermal buckling ͓108͔ of cylindrical shells.

Fracture Mechanics and Stress Singularity.

The symplectic methodology has also been applied to solve other applied engineering mechanics problems and the most important issue has been in fracture mechanics for structural elements with cracks and/or stress singularity ͓109-123͔. In these articles, the authors examined plane crack element ͓109͔, bimaterial crack singularity ͓110͔, wedge body problems ͓111-114,118,119,121͔, crack tip singularity ͓115͔, multimaterial crack singularity ͓116͔, application of the Dugdale model ͓117͔, and stress intensity factors ͓122͔. Mode III crack and stress singularity problems were also investigated ͓120,123͔.

Perturbation and Symplectic

Element. The development of numerical perturbation method using the symplectic approach was also reported mainly by Zhong and co-workers ͓124-129͔. In particular, Zhong and Sun examined three different methods based on symplectic perturbation ͓124͔ and further developed symplectic conserved perturbation method ͓125͔, finite element models ͓126͔, time-space harmony element ͓127͔, and numerical integration techniques ͓128͔. They also researched into dual variable symplectic principle ͓129͔.

2.9 Viscoelasticity. Structures with viscoelasticity have received relatively less attention. There are only three papers identified ͓104,130,131͔ where one focuses viscoelastic hollow circular cylinders ͓104͔ and the other two considered two-dimensional viscoelasticity without ͓130͔ and with thermoviscoelasticity ͓131͔.

Fluid Mechanics.

Besides solid mechanics, the symplectic methodology could also be applied to fluid mechanics problems, due mainly to Xu and co-workers ͓132-136͔ particularly for Stokes flow. In specific, the eigenvalue problem and eigensolutions for Stokes flow were modeled ͓132͔, while the analytical and numerical solutions in two-dimensional domain were reported ͓133͔. The influence of inlet radius on Stokes flow in a circular tube ͓134͔ and viscous fluid in lid-driven cavities ͓135͔ were further investigated. In nonlinear water waves driven by moving plates in shallow single-or double-layer fluid layers, dual variables associated with velocity potential and the wave elevation in a Hamiltonian system were introduced and the traveling waves solutions in elliptic cosine functions were presented ͓136͔.

Control.

It is common knowledge for researcher in applied mechanics that control models bear many similarities with respect to dynamical systems in terms of governing equation, solution methodology, and numerical techniques. The very close relations and analogy with respect to symplecticity were first established by Zhong and co-workers ͓36,39,41͔ who discussed the duality systems for optimal control problems ͓39͔ and state-space control systems ͓41͔. Recently, the time-varying optimal control problem ͓137͔ and optimal control for nonlinear dynamical systems ͓138͔ were analyzed.

2.12 Thermal Effects. The thermal effects on structures could also be solved using the symplectic method. One early work was the investigation of thermal stress for laminated plates ͓72͔. The subjected did not receive much attention subsequently. Recently, the thermal effects were again reported for thermal buckling of beams ͓70,71͔, circular plates ͓93͔, cylindrical shells ͓108͔, and for analysis of thermal stress intensity factors ͓139,140͔.

2.13 Functionally Graded Effects. Materials that exhibit functional graded properties could also be solved using the symplectic method ͓141-144͔. The topics of interest include statics and dynamics of plates ͓141͔, plane elasticity problems without ͓142͔ and with piezoelectric effects ͓143͔, and coupling of functionally graded and magneto-electro-elastic effects ͓141,144͔.

2.14 Piezoelectricity. The symplectic elasticity approach has been also extended to some cross-disciplinary research areas. One example is research in symplectic piezoelectricity ͓143,145-154͔ where the Hamiltonian systems are constructed to model the coupling of pressure and electricity for actuation and sensing purposes where driving force and/or electric power could be generated and vice versa. One early attempt was presented to develop the algorithms for analyzing piezoelectric materials ͓145͔. Subsequent research works include analysis of piezoelectric twodimensional transversely isotropic piezoelectric structures ͓146,147,151͔, cantilever beams ͓148,149͔, three-dimensional piezoelectric media ͓150,152͔, boundary layer phenomena ͓151͔, wedge body and crack singularity ͓121,123͔, composite structures such as laminated and cantilever plate ͓153,154͔, and statics and dynamics of functionally graded piezoelectric structures ͓143,152͔.

Electromagnetism and Waveguide.

Electromagnetism also bears many similarities with Hamiltonian systems and hence it can be solved using the symplectic methodology. This crossdisciplinary subject was initiated by Zhong and co-workers who applied the approach and developed semi-analytical solutions to electromagnetic wave guides ͓155-158͔. Zhong and Sun ͓159͔ later attempted to develop finite elements for electromagnetic resonant cavity.

Magneto-Electro-Elasticity.

The coupling of magnetism, electricity, and elasticity was also reported ͓141,144,160,161͔. Yao established symplectic solutions with the Saint-Venant principle for antiplane ͓160͔ and plane ͓161͔ magneto-electro-elastic solids. Lately, the statics and dynamics of functionally graded and layered magneto-electro-elastic plate/pipe ͓141͔ and plane problems for materials having coupled functionally graded and magneto-electro-elastic effects ͓144͔ were also reported.

Symplectic Space, Its Properties, and Correlation With Euclidean Space

All conservative real physical processes can be described by a suitable Hamiltonian system whose common mathematical fundamentals are the symplectic spaces. A symplectic space is different from a Euclidean space, which studies the metric properties such as length. It focuses on the study of area or the study of work and this is a mathematical structure present throughout the application of symplectic elasticity. The following presents a brief account of the fundamentals of symplectic space and subsequently discusses some of the special symplectic characteristics. Some proofs and more details could be referred to "Symplectic Elasticity" by Yao et al. ͓42͔.

For an n-dimensional linear space V defined in a real number field R, and a corresponding n-dimensional dual linear space VЈ, a 2n-dimensional phase space W in a real number field R can be defined as

W = V ϫ VЈ = ͭͯ ͩ q p ͪͯqV,pVЈͮ ͑1͒
which is constructed by two n-dimensional vectors, q and p, defined in V and VЈ, respectively. Although there is no direct relation between these vectors and they could have absolutely different dimensions in actual problems, the product of their corresponding elements always has a specific physical meaning. For instance, one vector could represent displacement while the other stress and the product become work.

Let ͗x , y͘ denote the symplectic inner product for any two component vectors x = ͕x 1 , x 2 , . . . ,x 2n ͖ T and y = ͕y 1 , y 2 , . . . , y 2n ͖ T , then the symplectic space requires that the following four rules be satisfied absolutely:

͑1͒ ͗x,y͘ = -͗y,x͘ ͑ 2a͒
͑2͒ ͗kx,y͘ = k͗x,y͘ k is an arbitrary real number ͑2b͒ ͑3͒ ͗x + ␥,y͘ = ͗x,y͘ + ͗␥,y͘ ␥ is an arbitrary vector in W ͑2c͒ ͑4͒ x = 0 if ͗x,y͘ = 0 for every vector y in W ͑2d͒

Any phase space that satisfies the four rules above for symplectic inner product constitutes a symplectic phase space or simply a symplectic space. Consider a kind of symplectic inner product ͗x , y͘ defined as ͗x,y͘ = ͑x,J 2n y͒ = ͚ i=1 n ͑x i y n+ix n+i y i ͒ = x T J 2n y ͑3͒

where

J 2n = ͫ 0 I n -I n 0 ͬ ͑4͒
is the unit symplectic matrix of dimension 2n, denoted briefly as J, and I n is the common unit matrix of dimension n. It is obvious that a unit symplectic matrix has the following properties:

͉J͉ = 1 ͑5a͒ J 2 = -I ͑5b͒ J T = J -1 = -J ͑5c͒
where ͉J͉ indicates the determinant of J. It is also obvious that the symplectic inner self-product of every vector must vanish, i.e., for every vector x, ͗x,x͘ = 0 ͑6͒

It is obvious that Eq. ͑3͒ satisfies the four rules of symplectic inner product defined in Eqs. ͑2a͒-͑2d͒, and therefore it forms a 2n-dimensional symplectic space. For example, a twodimensional real linear space with vectors, the symplectic inner product defined in Eq. ͑3͒ can be expressed as

͗x,y͘ = x 1 y 2 -x 2 y 1 ͑7͒
Equation ͑7͒ forms a two-dimensional symplectic space and the symplectic inner product represents the area of a parallelogram constructed by x, y as the adjacent sides.

It should be emphasized that the rules for a symplectic space defined in Eqs. ͑2a͒-͑2d͒ are global. There are various different definitions of symplectic inner product for a phase space and therefore there exist different symplectic spaces. The symplectic inner product defined in Eq. ͑3͒ is called the normal symplectic inner product in a 2n-dimensional real vector space R 2n and this specific, normal symplectic space is always referred to throughout for all works in symplectic elasticity.

Any two component vectors x, y defined in a symplectic space are symplectic orthogonal if their symplectic inner product vanishes as ͗x,y͘ = 0 ͑8͒

Otherwise, they are symplectic adjoint. Hence, from Eq. ͑2d͒, there exists a symplectic adjoint nonzero vector for every nonzero vector, or x and Jx must be symplectic adjoint if x 0.

Another important property is the definition of a symplectic matrix S, a 2n ϫ 2n matrix S, which satisfies

S T JS = J ͑9͒

A symplectic matrix has the following properties: ͑i͒ The inverse matrix of a symplectic matrix is a symplectic matrix, ͑ii͒ the transpose matrix of a symplectic matrix is a symplectic matrix, ͑iii͒ the determinant of a symplectic matrix is equal to either 1 or Ϫ1, and ͑iv͒ the product of two symplectic matrices is a symplectic matrix. If a 2n ϫ 2n matrix H acting on arbitrary 2n-dimensional vectors x , y satisfies ͗x,Hy͘ = ͗y,Hx͘ ͑ 10͒

then matrix H is a Hamiltonian matrix. It is obvious that a Hamiltonian matrix satisfies

͑JH͒ T = JH ͑11a͒ JHJ = H T ͑11b͒
There are many other important properties, definitions, and theorems in the symplectic space, such as the adjoint symplectic orthonormal vector set and normal adjoint symplectic orthonormal vector set. Some theorems without proofs and definitions include the facts that ͑i͒ an adjoint symplectic orthonormal vector set is a linearly independent vector set. ͑ii͒ Every adjoint symplectic orthonormal vector set in a 2n-dimensional symplectic space can be extended to an adjoint symplectic orthonormal basis. ͑iii͒ Every normal adjoint symplectic orthonormal vector set in a 2n-dimensional symplectic space can be extended to a normal adjoint symplectic orthonormal basis. ͑iv͒ The transformation matrix for normal adjoint symplectic orthonormal bases is a symplectic matrix. ͑v͒ If is an eigenvalue of a Hamiltonian matrix with multiplicity m, then -is also an eigenvalue with multiplicity m.

If zero is an eigenvalue of a Hamiltonian matrix H, then the multiplicity number is even. More details are available in Ref.

͓42͔.

One useful conclusion from the theorems and properties above is as follows. There exists an adjoint symplectic orthonormal basis composed of the basic eigenvectors and Jordan form eigenvectors of the Hamiltonian matrix H in a 2n-dimensional symplectic space. Through normalization, a normal adjoint symplectic orthonormal basis can be formed. The matrix formed by the column vectors is indeed a symplectic matrix.

The readers are referred to Table 1.1 in Yao et al. ͓42͔ for a correlation between a Euclidean space and a symplectic space and it describes a clear summary of the relevant concepts and properties.

Some Applications of the Symplectic Methodology in Applied Mechanics

A few examples using the symplectic elasticity methodology is presented in this section to illustrate the application of the method for exact solutions to some common engineering mechanics problems. Contrary to the semi-inverse approach ͓44-50͔ where bending of beam is related to solving a nonhomogeneous set of equations while vibration of beam is an eigenvalue problem, from the following examples, it is obvious that bending of beam becomes an eigenvalue problem with the eigenvalue indicating the spatial wave frequency, while the vibration or wave propagation problem becomes a multiple eigenvalue problem with eigenvalues indicating the spatial and temporal frequencies.

4.1 Timoshenko Beam Theory. The Hamiltonian system for a Timoshenko beam is established first and subsequently the analytical solutions for bending and wave propagation are presented ͓42͔.

For bending of a Timoshenko beam of length L, let the x-axis be along the beam midaxis before deformation and the xz-plane be the deflection plane. The beam deflection at a point along the x-axis is denoted as w ˜͑x͒ and ˜͑x͒ is the rotation of the cross section.

To construct an analytical system, first, the x-coordinate is modeled as the time coordinate of the Lagrangian system and the Hamiltonian system. Furthermore, denote q = ͕w,͖ T , q ˙= ͕w ˙, ˙͖T ͑12͒

where an overdot indicates differentiation with respect to the x-coordinate, i.e., ͑ ͒ ˙= d / dx. Application of the variational principle yields the Lagrange equation as

d dx ͩ ץL ץq ˙ͪ - ץL ץq = 0 ͑13͒
where L͑q , q ˙͒ is the Lagrange density function and the temporal coordinate t is replaced by the spatial coordinate x.

A state vector is defined as

v = ͭ q p ͮ ͑14͒
where the dual variable of q according to Legendre's transformation is

p = ץL ץq ˙͑15͒
Then the Hamiltonian density function can be obtained as H͑q,p͒ = p T q ˙-L͑q,q ˙͒ ͑ 16͒

The Hamiltonian dual system can be expressed as

q ˙= ץH ץp ͑17a͒ p ˙= ץH ץq ͑17b͒
which can also be expressed as

v ˙= Hv + h ͑18͒
where H is the Hamiltonian matrix and h = ͕h q , h p ͖ T is a 2n-vector and h q , h p are two n-vectors, which can be obtained from Eqs. ͑14͒-͑16͒, ͑17a͒, and ͑17b͒, and they are related to external loadings. The Hamiltonian density function, also known as the mixed energy density, can be expressed in terms of the Hamiltonian matrix ͓42͔.

The systems expressed in Eqs. ͑17a͒, ͑17b͒, and ͑18͒ are sys-tems of first-order equations in which the spatial differential of the state vector is only present on the left-hand-side. Hence, the original system of n second-order differential equations ͑13͒ can be transformed into a system of 2n first-order differential equations ͑18͒ and the transformation from Lagrange system to Hamiltonian system is completed.

The physical meaning of the dual variable p can be interpreted as follows. For a Timoshenko beam, it is easily derived that

p = ͕kGA͑w ˙-͒, EI ˙͖T = ͕F s , M ͖ T ͑19͒
where k, E, G, A, I, F s , and M are the shear correction factor, Young's modulus, shear modulus, cross-sectional area, second moment of area, shear force, and bending moment, respectively. Hence, the dual variable p can be interpreted as a the generalized internal force vector and the state vector becomes v = ͕w, ; F s , M ͖ T ͑20͒

4.1.1 Bending of Timoshenko Beam. For static bending of a Timoshenko beam ͓42͔, the Hamiltonian matrix is

H = ΄ 0 1 ͑kGA͒ -1 0 0 0 0 ͑EI͒ -1 0 0 0 0 0 0 -1 0 ΅ ͑21͒
The eigenvalue and eigenvector for H can be obtained from

H = - ͑22͒
For nontrivial solution, the eigenvalue = 0 has quadruple multiplicity. The corresponding eigenvectors are

v 0 ͑0͒ = 0 ͑0͒ = ͕1, 0, 0, 0 ͖ T ͑23a͒
which represents beam rigid body translation

v 0 ͑1͒ = 0 ͑1͒ + x 0 ͑0͒ = ͕x, 1, 0, 0 ͖ T ͑23b͒
which represents beam rigid body rotation

v 0 ͑2͒ = 0 ͑2͒ + x 0 ͑1͒ + 1 2 x 2 0 ͑0͒ = ͕ 1 2 x 2 + c, x, 0, EI ͖ T ͑23c͒
which represents pure bending of beam, where c is an arbitrary constant, and

v 0 ͑3͒ = 0 ͑3͒ + x 0 ͑2͒ + 1 2 x 2 0 ͑1͒ + 1 6 x 3 0 ͑0͒ = ͭ 1 6 x 3 + cx, 1 2 x 2 + EI kGA + c; -EI, EIx ͮ T ͑23d͒
which represents constant shear force bending of beam. The basic eigenvector pairs 0 ͑0͒ , 0 ͑1͒ and 0 ͑0͒ , 0 ͑2͒ are symplectic orthogonal, while 0 ͑0͒ , 0 ͑3͒ are symplectic adjoint, i.e., ͗ 0 ͑1͒ , 0 ͑0͒ ͘ =0,

͗ 0 ͑2͒ , 0 ͑0͒ ͘ = 0, and ͗ 0 ͑3͒ , 0 ͑0͒ ͘ = EI 0. It is also easy to prove that ͗ 0 ͑1͒ , 0 ͑2͒ ͘ = EI 0 and by assigning c =-EI / 2GA, ͗ 0 ͑2͒ , 0 ͑3͒ ͘ = 0. Hence, 0 ͑0͒ , 0 ͑1͒ , 0 ͑2͒
, and 0 ͑3͒ form a set of adjoint symplectic orthonormal basis. Consequently, all basic eigenvectors and Jordan form eigenvectors for static bending of Timoshenko beam are established.

In the presence of external loading and/or external moment, h 0 in Eq. ͑18͒ and a particular solution has to be determined. This particular solution depends cannot be explicitly derived without the type of loadings being known. This will be left as an exercise to the readers.

Wave Propagation for Timoshenko Beam.

For wave propagation in a Timoshenko, the following solutions are stated without derivation. The readers are referred to Zhong ͓37͔ and Yao et al. ͓42͔ for more details.

In this case, the Hamiltonian matrix H is given by

H = ΄ 0 1 ͑kGA͒ -1 0 0 0 0 ͑EI͒ -1 -2 A 0 0 0 0 -2 I -1 0 ΅ ͑24͒
where and are the material density and the wave frequency, respectively. The eigenvalue must satisfy that 2 be real. There are three different cases for 2 , which can be classified according to a critical frequency as

cr 2 = GA I ͑25͒
͑i͒ For 2 Ͼ cr 2 , there are two negative roots for 2 . The solution in a Timoshenko beam indicates propagation of two pairs of waves with velocities / 1 and / 2 , and traveling along -x and +x, respectively. ͑ii͒ For 2 Ͻ cr 2 , there are a positive root and a negative root for 2 . The negative root indicates propagation of a pair of waves along the positive and negative directions of x, similar to case ͑i͒ above The positive root indicates local vibration capable of creating resonance. ͑iii͒ For 2 = cr 2 , there are a negative root and a zero root for 2 . The zero root indicates that there are two zero roots for and there will be Jordan form solutions. The wave propagation solutions for the roots are given in Ref. ͓42͔.

Bending of Rectangular

Plates. Consider an isotropic rectangular Kirchhoff plate with uniform thickness h, length a, width b, Young's modulus E, and Poisson's ratio . Here only thin plates having at least a pair of opposite sides simply supported are considered. The plate, bounded within a domain -a / 2 Յ x Յ a / 2 and 0 Յ y Յ b where sides y =0,b are always simply supported, as shown in Fig. 1, is subjected to uniformly distributed load q. There are totally six cases with different boundary conditions. The analytical bending solutions are presented as follows and the details can be referred to Ref. ͓90͔. Here C, S, and F denote the clamped support, the simple support, and the free support, respectively.

In this case, the Hamiltonian matrix is where and D are the Poisson ratio and the flexural rigidity of plate, and x and y are the bending moment functions, and y and xy are the curvature and twisting curvature, respectively. For a SSSS plate, the xy-midplane deflection w͑x , y͒ solution can be expressed as ͓90͔

H = ΄ 0 ץ ץy D͑1 -2 ͒ 0 - ץ ץy 0 0 2D͑1 -͒ 0 0 0 - ץ ץy 0 - 1 D ץ 2 ץy 2 ץ ץy 0 ΅ ͑26͒
w = q 24D ͑y 4 -2by 3 + b 3 y͒ + 2q Db ͚ n=1 ϱ ͓ n x sinh͑ n x͒ -cosh͑ n x͒͑2 + ␣ n tanh ␣ n ͔͒sin͑ n y͒ n 5 cosh ␣ n ͑28͒
For a SFSF plate, simply supported at y =0,b and free at x = Ϯ a / 2, the deflection w͑x , y͒ solution is ͓90͔

w = q 24D ͑y 4 -2by 3 + b 3 y͒ + 4q ͑1 -͒bD ͚ n=1 ϱ ͕ cosh͑ n x͓͒␣ n ͑1 -͒cosh ␣ n -͑1 + ͒sinh ␣ n ͔ -n x͑1 -͒ sinh ␣ n sinh͑ n x͖͒sin͑ n y͒ n 5 ͓␣ n ͑1 -͒ -͑3 + ͒cosh ␣ n sinh ␣ n ͔ ͑29͒
For a SCSC plate, simply supported at y =0,b, and clamped at x = Ϯ a / 2, the deflection w͑x , y͒ solution is ͓90͔

w = q 24D ͑y 4 -2by 3 + b 3 y͒ + 4q bD ͚ n=1 ϱ ͕͕e 4␣ n sin͑ n y͓͒-2͑3␣ n + n x͒cosh͑␣ n -n x͒ + ͑␣ n + n x͒cosh͑3␣ n -n x͒ -5␣ n cosh͑␣ n + n x͒ + 3 n x cosh͑␣ n + n x͒ + 2␣ n cosh͑3␣ n + n x͒ -2 n x cosh͑3␣ n + n x͒ + 8 cosh͑2␣ n ͒cosh͑ n x͒sinh ␣ n -4␣ n ͑␣ n -n x͒sinh͑␣ n + n x͔͖͒/͓ n 5 ͑1 + 8␣ n e 4␣ n -e 8␣ n ͔͖͒ ͑30͒
For a SFSC plate, simply supported at y =0,b, free at x =-a / 2, and clamped at x = a / 2, the deflection w͑x , y͒ solution is ͓90͔

w = q 24D ͑y 4 -2by 3 + b 3 y͒ + 8q bD ͚ n=1 ϱ ͕͕e 4␣ n sin͑ n y͕͓͒5 + 4␣ n 2 ͑ -1͒ 2 + ͑2 + ͔͒cosh͑␣ n -n x͒ + 3 cosh͑3␣ n + n x͒ + cosh͑ n x͕͓͒1 + 4␣ n 2 ͑ -1͒ + ͔cosh ␣ n + ͑-3 + ␣ n -2͒cosh͑3␣ n ͒ + n x͓-5 sinh ␣ n -͑ -1͒sinh͑3␣ n ͔͖͒ + ͕ -n x͑1 + 2͒cosh ␣ n + n x͑ -1͒cosh͑3␣ n ͒ + ͓1 + 4␣ n 2 ͑ -1͒ + ͔sinh ␣ n + ͑␣ n -1͒sinh͑3␣ n ͖͒sinh͑ n x͒ -␣ n ͕cosh͑ n x͕͒4 n x͑ -1͒cosh ␣ n -3 sinh ␣ n + ͕cosh͑3␣ n ͒ + ͓5 -4 + 4 cosh͑2␣ n ͔͒sinh␣ n ͖͖ + ͓͑3 -7͒cosh ␣ n + 2 cosh͑3␣ n ͒ + 4 n x͑ -1͒͑2 -1͒sinh ␣ n + sinh͑3␣ n ͔͒sinh͑ n x͖͒ + n x sinh͑␣ n -n x͒ -͓␣ n ͑ -3͒ -n x͑ -1͒͑ + 3͔͒sinh͑3␣ n + n x͖͖͒/͕ n 5 ͕͑ -1͒͑ + 3͒ + e 8␣ n ͑ -1͒͑ + 3͒ -2e 4␣ n ͓5 + 8␣ n 2 ͑ -1͒ 2 + ͑2 + ͔͖͖͖͒ ͑31͒
For a SSSF plate, simply supported at y =0,b and x =-a / 2, and free at x = a / 2, the deflection w͑x , y͒ solution is ͓90͔

w = q 24D ͑y 4 -2by 3 + b 3 y͒ + 2q bD ͚ n=1 ϱ ͕͕e ␣ n -n x ͕-2e 2͑␣ n + n x͒ ͓-1 + 3␣ n ͑ -1͒ + n x͑ -1͒ -͔ + e 2 n x ͑-2 + ␣ n + n x͒͑ -1͒ ϫ͑3 + ͒ -2e 4␣ n ͓1 + 3␣ n ͑ -1͒ + n x͑ -1͒ + ͔ + 2͓-1 -␣ n + n x + ͑-1 + ␣ n -n x͔͒ + e 2␣ n ͓6 + 4␣ n 2 ͑ -1͒ 2 + 3 n x͑ -1͒ 2 -␣ n ͑5 + 4 n x͒͑ -1͒ 2 + 2 2 ͔ + 2e 6␣ n +2 n x ͓1 + ␣ n ͑ -1͒ + -n x͑ -1͔͒ + e 4␣ n +2 n x ͓-4␣ n 2 ͑ -1͒ 2 + 3 n x͑ -1͒ 2 + ␣ n ͑-5 + 4 n x͒͑ -1͒ 2 -2͑3 + ͒ 2 ͔͖sin͑ n y͖͒/͕ n 5 ͑ -1͓͒8␣ n e 4␣ n ͑ -1͒ + ͑1 -e 8␣ n ͒͑3
+ ͔͖͖͒ ͑32͒

For a SSSF plate, simply supported at y =0,b and x =-a / 2, and clamped at x = a / 2, the deflection w͑x , y͒ solution is ͓90͔ highlighted here that all geometric and natural boundary conditions are satisfied and the analytical solutions can be considered exact to a very large extent except some stages toward the end of the analysis where numerical solutions have to be sought. In addition, the twisting moment at the corners derived from static bending equilibrium are also satisfied.

w = q 24D ͑y 4 -2by 3 + b 3 y͒ + 4q bD ͚ n=1 ϱ ͕͕e 4␣ n sin͑ n y͕͒-2͑3␣ n + n x͒cosh͑␣ n -n x͒ + ͑␣ n + n x͒cosh͑3␣ n -n x͒ -5␣ n cosh͑␣ n + n x͒ + 3 n x cosh͑␣ n + n x͒ + 2␣ n cosh͑3␣ n + n x͒ -2 n x cosh͑3␣ n + n x͒ + 8 cosh͑2␣ n ͒cosh͑ n x͒sinh͑␣ n ͒ -4␣ n ͑␣ n -n x͒sinh͑␣ n + n x͖͖͒/͓͑1 + 8e 4␣ n -e 8␣ n ͒ n 5 ͔͖ ͑33͒ 4
For such a corner-supported plate, the deflection w͑x , y͒ solution can be expressed as ͓91͔

w = q 24D͑-1 + ͒͑1 + ͒ 2 ͕-6a 2 ͑1 + ͒͑-x 2 + y 2 ͒ + 2b 2 ͓-͑1 + 5͒x 2 -͑-3 + ͒͑2 + ͒y 2 ͔ + ͑1 + ͓͒-x 4 + 6x 2 y 2 + ͑-2 + ͒y 4 ͔͖ + c 1 + ͚ n=1 ϱ ͭ f n e n x ͩ ͓1 + -͑3 + ͒cos 2 ͑ n b͔͒cos͑ n y͒ n D͑1 -͒ 2 - y sin͑ n y͒ D͑1 -͒ ͪ+f ¯ne ¯nx ͩ ͓1 + -͑3 + ͒cos 2 ͑ ¯nb͔͒cos͑ ¯ny͒ ¯nD͑1 -͒ 2 - y sin͑ ¯ny͒ D͑1 -͒ ͪ + f -n e -n x ͩ ͓1 + -͑3 + ͒cos 2 ͑ -n b͔͒cos͑ -n y͒ -n D͑1 -͒ 2 - y sin͑ -n y͒ D͑1 -͒ ͪ + f ¯-n e ¯-n x ͩ ͓1 + -͑3 + ͒cos 2 ͑ ¯-n b͔͒cos͑ ¯-n y͒ ¯-n D͑1 -͒ 2 - y sin͑ ¯-n y͒ D͑1 -͒ ͪͮ ͑34͒
4.4 Vibration of a Rectangular Plates. For vibration of rectangular plate, the geometry, as illustrated in Fig. 2, is adopted except that the supports are not at the corners but they are along the boundaries with various conditions to be considered. At least one-pair of the sides at x = Ϯ a are simply supported, while the other sides y = Ϯ b could be free, simply supported, or clamped. Again only the analytical frequency solutions are presented and the details can be referred to Refs. ͓89,100͔ and Ref. ͓95͔.

From the property of symplectic adjoint orthogonality of eigenvectors and expansion of eigenvectors, the state vector for free vibration can be expanded as

v = ͚ n=1 ϱ ͓f 1 ͑n͒e n x n ͑y͒ + f 2 ͑n͒e -n x -n ͑y͒ + f 3 ͑n͒e n Јx n Ј͑y͒ + f 4 ͑n͒e -n Јx -n Ј ͑y͔͒ ͑35͒
where f i ͑n͒ , i =1,2,3,4 are unknown functions depending on the boundary conditions at x = Ϯ a and

n = Ϯ ͱ n ͱ h D + n 2 2 4b 2 , 3 = ͱ-2 n ͱ h D - n 2 2 4b 2 , ͩ1= Ϯ i n 2b ͪ ͑36a͒ n Ј = Ϯ ͱ-n ͱ h D + n 2 2 4b 2 , 1 = ͱ2 n ͱ h D - n 2 2 4b 2 , ͩ3= Ϯ i n 2b ͪ ͑36b͒
The circular frequencies n are the only unknowns in the state vector in Eq. ͑35͒. Satisfaction of the boundary conditions at x = Ϯ a will lead to the frequency equation. For a SSSS plate, combining Eq. ͑35͒ with the boundary conditions yields a frequency relation sinh n a cosh n a sinh n Јa cosh n Јa = 0 ͑37͒ which can be subsequently solved for the free vibration frequency n as

n = ͫ m 2 2 ͑2a͒ 2 + n 2 2 ͑2b͒ 2 ͬ ͱ D h , ͑m,n = 1,2, . . .͒ ͑38͒
which is identical to the well-known solution in any text on thin plate vibration ͓48,50͔.

For plates with other boundary conditions at x = Ϯ a, similar procedure as described above is followed. Combining Eqs. ͑35͒, ͑36a͒, and ͑36b͒ with the relevant boundary conditions yield the following frequency relations: where k 2 = ␣ 2 / l 2 , ␣ 2 = n ͱ h / D, and l 2 = n 2 2 / 4b 2 .

• CSCS ͑ n 2 + n Ј 2 ͒sinh͑2 n a͒sinh͑2 n Јa͒ = 2 n n Ј͓cosh͑2 n a͒cosh͑2 n Јa͒ -1͔ ͑39͒ • SSCS n cosh͑2 n a͒sinh͑2 n Јa͒ -n Ј sinh͑2 n a͒cosh͑2 n Јa͒ = 0
These frequency relations have to be solved numerical for the frequency solutions. Unlike the classical semi-inverse method ͓44-50͔, the frequency relations above are derived analytically and rigorously step-by-step in a rational manner without introducing any trial functions. Furthermore, unlike the semi-inverse trial functions, which only satisfy the geometric boundary conditions, these exact relations satisfy all natural and geometric boundary conditions at the outset. However, no closed form solutions for the transcendental frequency relation for n in Eqs. ͑39͒-͑43͒ could be obtained. The numerical solutions are available in Ref. ͓95͔.

Conclusions

This paper presents a comprehensive review of the previous works on the theory and application of the symplectic methodol-ogy in applied mechanics and engineering. It begins with a description of Feng's works on numerical symplectic algorithm in the middle 1980s, a discussion on Zhong's works on analytic symplectic elasticity since the beginning of the 1900s, and the many other works using the approach for closed form or nearly closed form solutions to elasticity problems for which exact or closed formed solutions have been impossible based on the semi-inverse method of Navier, Lévy, and Timoshenko ͓44-50͔. A unique feature of this method is that bending of plate becomes an eigenvalue problem and vibration becomes a multiple eigenvalue problem. This paper also presents an introductory background and fundamentals of symplectic space and Hamiltonian dual system. This paper finally ends with a presentation of some analytical solutions to some beam and plate problems for bending and vibration for which only numerical solution are considered possible prior to symplectic methodology.

Fig. 1

 1 Fig. 1 Plan of a thin plate with support conditions

  Fig. 2 Geometry and corner support conditions of a rectangular plate

2

  Ј͓k 4 -͑1 -͒ 2 ͔ + 2 n n Ј͓k 4 + ͑1 -͒ 2 ͔cosh͑2 n a͒cosh͑2 n Јa͒ + ͑ n 2 + n Ј 2 ͓͒͑1 -2͒k 4 -͑1 -͒ 2 ͔sinh͑2 n a͒sinh͑2 n Јa͒ = 0 ͑41͒ • SSFS n ͓k 2 -͑1 -͔͒ 2 cosh͑2 n a͒sinh͑2 n Јa͒ = n Ј͓k 2 + ͑1 -͔͒ 2 sinh͑2 n a͒cosh͑2 n Јa͒ ͓k 2 -͑1 -͔͒ 4 + n Ј 2 ͓k 2 + ͑1 -͔͒ 4 ͖sinh͑2 n a͒sinh͑2 n Јa͒ = 2 n n Ј͓k 4 -͑1 -͒ 2 ͔ 2 ͓cosh͑2 n a͒cosh͑2 n Јa͒ -1͔ ͑ 43͒

  nonlinear beam vibration and dynamics ͓63-65͔, Galerkin method ͓63͔, beam finite element ͓64͔ and frames ͓65͔, Zhong et al. ͓66͔ for curved beams, Xu et al. ͓67͔ and Ma et al. ͓68͔ for elastic structural vibration, Lü et al. ͓69͔ for bending of beams resting on elastic foundation, and for nonlinear thermal buckling of elastic beams ͓70,71͔.

Acknowledgment

This work described in this paper was supported by a grant from the City University of Hong Kong ͑Project No. 9667036͒. The authors are also grateful to Professor W.A. Yao for advice and for providing some of the references.