Ioannis Z Emiris

Christos Konaxis

Ilias S Kotsireas

Clément Laroche

Matrix Representations by Means of Interpolation

Keywords: Matrix representation, sparse resultant, Chow form, ray shooting, randomized algorithm, implicitization, space curve, Maple implementation

We examine implicit representations of parametric or point cloud models, based on interpolation matrices, which are not sensitive to base points. We show how interpolation matrices can be used for ray shooting of a parametric ray with a surface patch, including the case of highmultiplicity intersections. Most matrix operations are executed during pre-processing since they solely depend on the surface. For a given ray, the bottleneck is equation solving. Our Maple code handles bicubic patches in ≤ 1 sec, though numerical issues might arise. Our second contribution is to extend the method to parametric space curves and, generally, to codimension > 1, by computing the equations of (hyper)surfaces intersecting precisely at the given object. By means of Chow forms, we propose a new, practical, randomized algorithm that always produces correct output but possibly with a non-minimal number of surfaces. For space curves, we obtain 3 surfaces whose polynomials are of near-optimal degree; in this case, computation reduces to a Sylvester resultant. We illustrate our algorithm through a series of examples and compare our Maple prototype with other methods implemented in Maple i.e., Groebner basis and implicit matrix representations. Our Maple prototype is not faster but yields fewer equations and seems more robust than Maple's implicitize; it is also comparable with the other methods for degrees up to 6.

Introduction

In manipulating geometric objects, it is essential to possess robust algorithms for changing representation. This paper considers three fundamental representations, namely implicit, parametric (possibly with base points), and point clouds. We offer two results on algorithms operating on the latter two representations. Our algorithms construct implicit equations or implicit matrix representations of parametric curves and (hyper)surfaces by means of some type of interpolation through points. In general, approaches to implicitization include resultants, Gröbner bases, moving lines and surfaces, and interpolation techniques.

Implicit matrix representations are quite robust, since they do not require developing the implicit equation; instead, they reduce geometric operations on the object to matrix algebra, e.g. [START_REF] Emiris | Geometric operations using sparse interpolation matrices[END_REF][START_REF] Busé | The surface/surface intersection problem by means of matrix based representations[END_REF]. Our matrix representation uses interpolation matrices. The method has been developed for plane curves, and (hyper)surfaces. The columns of the matrix are indexed by a superset of the implicit monomials, i.e., those in the implicit polynomial with non-zero coefficient [START_REF] Emiris | Sparse implicitization by interpolation: Characterizing non-exactness, and an application to computing discriminants[END_REF]. This monomial set is determined quite tightly for parametric models, by means of the sparse resultant of the parametric polynomials, thus exploiting the sparseness of the parametric and implicit polynomials. When the object is given as a point cloud, this superset is taken to be all monomials with total degree ≤ δ , where δ is an educated guess of the total degree of the implicit polynomial.

Our contribution is twofold. First, we reduce certain geometric operations, most notably ray shooting, to univariate solving, after preprocessing by means of matrix operations. To check whether the computed intersections of the ray with the implicit surface lie on the given parametric patch, we employ inversion by means of bivariate solving. Our method is not sensitive to intersections of high multiplicity, nor to surfaces with base points, and it can be applied to objects given either parametrically or as a point cloud. But it exhibits numerical issues when the interpolation matrix dimension is large, say of a few hundreds. On the other hand, the straightforward method of solving the system by equating the parametric expressions to the ray coordinates typically exhibits robustness issues, especially at the presence of base points. We present experiments based on an open-source, Maple implementation.

Second, we extend the method to varieties of codimension higher than 1. Given a variety in parametric form, and starting with the powerful and classic theory of Chow forms, which generalizes resultant theory, we design an original, practical algorithm that uses randomization in order to compute implicit hypersurfaces containing the variety. The question is reduced to computing a (sparse) resultant. For space curves in 3D-space, our method yields a small number of implicit polynomials defining the curve set-theoretically. Their degree is near to the optimal degree and the algorithm is significantly simpler and is expected to be much faster than computing the optimal-degree polynomials. Self-intersections do not affect the degree of the output.

Previous work

A direct method to reduce implicitization to linear algebra is to construct a square matrix M, indexed by all possible monomials in the implicit equation (columns) and different values (rows) at which all monomials get evaluated. The vector of coefficients of the implicit equation is in the kernel of M, even in the presence of base points. This idea has been extensively used [START_REF] Corless | Numerical implicitization of parametric hypersurfaces with linear algebra[END_REF][START_REF] Dokken | Approximate implicitization[END_REF][START_REF] Marco | Implicitization of rational surfaces by means of polynomial interpolation[END_REF][START_REF] Bajaj | Hermite interpolation of rational space curves using real algebraic surfaces[END_REF]. In [START_REF] Emiris | Sparse implicitization by interpolation: Characterizing non-exactness, and an application to computing discriminants[END_REF][START_REF] Emiris | Geometric operations using sparse interpolation matrices[END_REF], sparse elimination theory is employed to predict the implicit monomials and build the interpolation matrix.

Resultants, and their matrix formulae, have been used to express the implicit surface equation, e.g., in [START_REF] Manocha | The implicit representation of rational parametric surfaces[END_REF], under the assumption of no base points. This technique reduces ray shooting (and surface-surface intersection) to substituting the parameterization of one object into the resultant matrix, then solving the resulting univariate polynomial.

Our method, unlike previous methods, works also for objects defined by point clouds. Moreover, the interpolation matrices are mostly numeric, except for one symbolic row that captures the second object, namely the ray or the second surface. This allows for several preprocessing operations to be executed numerically. In [START_REF] Emiris | Geometric operations using sparse interpolation matrices[END_REF], Membership and Sidedness predicates are expressed in terms of operations on the interpolation matrix. This approach is extended to ray shooting here.

A modern method for representing implicit equations by matrices is based on the theory of syzygies, e.g., [START_REF] Sederberg | Implicitization using moving curves and surfaces[END_REF][START_REF] Busé | Implicit matrix representations of rational Bézier curves and surfaces[END_REF], and handles base points. It yields matrices whose entries are linear polynomials in the implicit variables and they indirectly represent implicit objects: their rank drops exactly on the curve or surface. They allow for geometric operations, such as surface-surface intersection [START_REF] Busé | The surface/surface intersection problem by means of matrix based representations[END_REF] and, more recently, ray shooting [START_REF] Shen | A line/trimmed NURBS surface intersection algorithm using matrix representations[END_REF], to be executed by linear algebra operations. Their advantage is that the matrices are much smaller than interpolation matrices, and allow for inversion by solving an eigenproblem on these matrices.

Chow forms have been studied in computational algebra, in particular for varieties of codimension ≥ 2, since they provide a method to describe the variety by a single polynomial, e.g. [START_REF] Gelfand | Discriminants, Resultants and Multidimensional Determinants[END_REF][START_REF] Dalbec | Introduction to chow forms[END_REF]. Given a d-dimensional irreducible variety V in n-dimensional space, and d + 1 linear forms H i on V , the Chow form R V of V is a polynomial whose variables are the coefficients of the linear forms and that vanishes whenever the linear forms vanish on V .

For example, the Chow form of a space curve in projective 3D-space is a polynomial in the indeterminates u 0 , u 1 , u 2 , u 3 , v 0 , v 1 , v 2 , v 3 that vanishes whenever the planes

H 0 = u 0 x 0 + u 1 x 1 + u 2 x 2 + u 3 x 3 = 0, (1)
H 1 = v 0 x 0 + v 1 x 1 + v 2 x 2 + v 3 x 3 = 0, intersect on the curve.
There are many ways to describe R V . The most common is by using the maximal minors of the (d + 1) × (n + 1) matrix whose rows are the normals to the hyperplanes H i , see Sec. 4.1. These are known as Plücker coordinates or brackets and are the variables of R V . Equivalently, the dual Plücker coordinates or dual brackets can be used; these are the maximal minors of a (nd) × (n + 1) matrix whose rows are nd points that span the intersection L of the hyperplanes H i . There are algorithms to recover the implicit or affine equations of V from R V [START_REF] Gelfand | Discriminants, Resultants and Multidimensional Determinants[END_REF][START_REF] Sturmfels | Algorithms in Invariant Theory[END_REF]. Here we adopt a practical method that avoids conversion but may yield a multiple of R V . Our algorithm requires a parametric representation of V .

Interpolation matrices

This section describes implicitization by an interpolation matrix M for a plane curve or a (hyper)surface given either in parametric form, or, as a point cloud. If the input object is affinely parameterized by rational functions

x i = f i (t), i = 1, . . . , n, t = (t 1 , . . . ,t n-1 ) (2)
then it is possible to predict the implicit monomials. This set is included in the predicted (implicit) polytope computed by software ResPol1 [START_REF] Emiris | An oracle-based, outputsensitive algorithm for projections of resultant polytopes[END_REF]. If the input is a point cloud, we consider a coarse estimation of the monomial set by guessing the total degree of the variety and taking all the monomials of that degree or lower. The set S of monomials is used to construct a numerical matrix M, expressing a linear system whose unknowns are the coefficients c i of the monomials with exponents in S in the implicit polynomial, as discussed above. To obtain the linear system in the c i we substitute x i by its rational parametric expression

x i = f i (t) into equation ∑ |S| i=1 c i x a i = 0, x a i := x a i1 1 • • • x a in n
, and evaluate the parameters t at generic points (randomized in practice) τ k ∈ C n-1 , k = 1, . . . , µ, µ ≥ |S|, avoiding values that make the denominators of the parametric expressions close to 0.

Letting m i = m i (t) denote the monomial x a i after substituting x i by its parametric expression in (2), and m i | t=τ k its evaluation at t = τ k , we end up with a matrix M of the form:

M =     m 1 | t=τ 1 • • • m |S| | t=τ 1 . . . • • • . . . m 1 | t=τ µ • • • m |S| | t=τ µ     . (3)
Typically µ = |S| for performing exact kernel computation, and µ = 2|S| for approximate numeric computation. When constructing matrix M it is assumed that the parametric hypersurface is sampled sufficiently generically by evaluating the parametric expressions at random points τ k ∈ C n-1 . Let us recall some further properties. Let M be the (|S| -1) × |S| numeric matrix obtained by evaluating the monomials S at |S| -1 points τ k , k = 1 . . . , |S| -1. We obtain the |S| × |S| matrix M(x), which is numeric except for its last row, by appending the row of monomials S to matrix M :

M(x) = M S(x) . (4)
Notice that matrices M , M and M(q), for q a point lying on the hypersurface have the same kernel of corank r. Matrix M(x) has an important property: Lemma 1. [EKK15, Lem.7] Assuming M is of full rank, then det M(x) equals the implicit polynomial p(x) up to a constant.

Ray shooting

This section examines the geometric operation of ray shooting implemented by means of linear algebra so as to avoid computation of the implicit equation. Our main motivation is to complete the properties of interpolation matrices used in expressing the implicit form of objects initially given either by a parameterization or a point cloud. Our experimental results are preliminary but show that our approach can be competitive to well-established methods possessing optimized implementations.

In particular, the representation by interpolation matrices reduces ray shooting at a surface patch to operations on the interpolation matrix and univariate equation solving. An analogous approach applies to computing surface-surface intersections by determining the implicit curve of the intersection in the plane of the parameters of one surface, discussed in Sec. 5.

Given a parameterization as in (2), the (hyper)surface is represented by an interpolation matrix of the form (4). The standard case of surfaces is n = 3. We substitute the parametric ray into the monomials S(x) of the surface's implicit support. This defines the last row of the interpolation matrix. Let the ray be given by

x i = r i (ρ), i = 1, . . . , n, ρ > 0,
where the r i ∈ R[ρ] are linear polynomials in ρ. The ray may lie in arbitrary dimension but again, the standard 3D case is n = 3. Starting with the square matrix M in expression (4), we substitute x i by r i (ρ) in the monomials of the last row S(x):

M(ρ) = M S(ρ) , (5)
where S(ρ) is a row vector of monomials in ρ. The intersection points correspond to the real roots of polynomial equation p(ρ) = det M(ρ) = 0, which expresses the implicit equation by Lem. 1, evaluated at x i = r i (ρ). Of course, if the ray is totally included in the (hyper)surface, then p is identically 0. In the following, we assume that it is not the case. Lastly, our algorithm uses inversion of the parameterization to check whether the intersection points lie on the patch.

The structure of M(ρ) allows us to perform numeric operations so as to obtain another matrix with same determinant, that is easier to compute. Our goal is to use some preprocessing time, without knowledge of the ray, so as to accelerate the actual ray shooting computation. This preprocessing can be considered as a better representation of the surface, amenable to geometric operations. We apply a PLU-decomposition or QR-decomposition on M .

Computation

It shall be easier to discuss matrix operations on the transpose matrices, namely M T and (M ) T , although mathematically equivalent with operations on the original matrices. Below, we employ the PLU-decomposition, but we have also experimented with the QR-decomposition with analogous results.

First, consider the standard PLU-decomposition on |S| × |S| -1 matrix (M ) T :

(M ) T = PLU,
where P, L are square |S| × |S| permutation and lower-triangular matrices, respectively, where L has a unit diagonal, and U is a numeric rectangular |S| × |S| -1 upper triangular matrix, with the last row equal to zero. We use the above decomposition to derive an analogous decomposition of M(ρ) T . Clearly, P -1 = P T , hence we can write square |S| × |S| matrix L -1 P T M(ρ) T as follows:

L -1 P T [ (M ) T | S(ρ) T ] = [U | L -1 P T S(ρ) T ], (6)
where [ • | • ] indicates concatenation of matrices, and S(ρ) is a row vector, hence the right-hand side is clearly a square |S| × |S| matrix. Now, det M(ρ) = det(L -1 P T M(ρ) T ) because det L = 1. To compute the determinant, we use the right-hand side of equation ( 6), which is a square upper triangular matrix. Its diagonal is the diagonal of U with one last element equal to the last element of column L -1 P T S(ρ) T . Since U is numeric, the determinant equals, up to a constant nonzero factor, the last element on the column L -1 P T S(ρ) T . Hence, it suffices to store the last row of L -1 P T , which depends on the surface, and compute its inner product with vector S(ρ) T , for any given ray. This discussion establishes the following:

Lemma 2. Given M(ρ) T = [(M ) T | S(ρ)]
as defined above, consider the PLU-decomposition (M ) T = PLU, which can be computed at preprocessing, and does not require knowledge of the particular ray. It is possible to compute the determinant det M(ρ), up to a nonzero scalar factor, as the inner product of the last row of L -1 P T multiplied by column vector S(ρ) T . An analogous result follows if one computes the QR-decomposition (M ) T = QR.

Given p(ρ) = det M(ρ), we employ any state-of-the-art real solver to approximate all positive real roots ρ 0 . For every root ρ j > 0, for some j, we obtain the intersection point x i = r i (ρ j ) in implicit space. An alternative is to use solvers that only require values of the polynomial p(ρ) and, eventually, its derivative p (ρ), such as the Newton-Raphson method. We have experimented with this method and the results are comparable in terms of numerical accuracy and speed.

The last step consists in computing the corresponding vector (t 1 ,t 2 , . . . ,t n-1 ) ∈ R n-1 that maps to this implicit point, by inverting the parameterization. Hence we check whether the point lies on the patch of interest. Assuming a normal and proper parameterization, there is a unique solution (t 1 ,t 2 , . . . ,t n-1 ) ∈ R n-1 such that r i (ρ j ) = f i (t 1 ,t 2 , . . . ,t n-1 ). Since we are dealing with patches here, we may not want to consider the missing points of the parameterization if we do not have a normal (surjective) parameterization. In that situation, the normal parameterization hypothesis can be forgotten.

There are several ways to solve the inversion problem, including methods for solving the overconstrained system r i (ρ j ) = f i (t), i = 1, . . . , n for every ρ j . Alternatively, we have used real solving methods applied to the following well-constrained system of n polynomials in t 1 , . . . ,t n-1 : r i (ρ j )f i (t) = 0, for i = 1, . . . , n. Then, we substitute all real roots into the last polynomial r n (ρ j )f n (t), and keep the root that satisfies it, or best satisfies it if approximate numerical computation is involved. This approach works satisfactorily for n = 3 in our experiments.

Asymptotic and practical complexity

Let us first consider the asymptotic bit complexity, if we omit preprocessing. Let D be the maximum total degree of any monomial in S, hence polynomial p(ρ) has degree D. Real solving of this polynomial costs O * (D 3 |S|), e.g. [START_REF] Tsigaridas | On the complexity of real root isolation using continued fractions[END_REF]. Inversion is then required to check which solutions belong to the relevant patch; it is reduced to bivariate solving and evaluation of the polynomials r i (ρ j )f i (t) defined by the parametric equations, for every real root ρ j of p(ρ) = 0. If these polynomials have total degree δ in s,t and bitsize τ, then each system is solved in bit complexity O * (δ 12 + δ 10 τ2 ), see e.g. [START_REF] Diochnos | On the complexity of real solving bivariate systems[END_REF], and this is repeated for each of the O(D) roots ρ j . It turns out that this is the bottleneck in practice, which can be explained by the analysis since D is only a few times larger than d. For example, the bicubic surface has D = 18 and δ = 6.

We now experiment with code developed on top of the basic method 2 , running Maple 14 on a laptop with a 3 GHz Intel Core i5 processor. Runtimes are averaged over at least 3 runs.

Example 1. The standard benchmark of a bicubic surface has an implicit equation of total degree 18, with a support of 715 monomials. Here, the predicted support is optimal. However, matrix size makes computations such as PLU-or QR-decompositions over Q infeasible. We use floating-point numbers, but then numerical error becomes significant. A concrete ray we experimented with was x 1 = ρ -13, x 2 = 12 + ρ, x 3 = 3 -5ρ. If we omit preprocessing, the total runtime is up to 0.12 sec, and the bottleneck is inversion.

Our method may be juxtaposed to the direct approach of solving the 3×3 system f i (t 1 ,t 2 ) = r i (ρ), i = 1, 2, 3, which takes about 0.02 sec on Maple and is substantially more accurate numerically, on this example. But numerical issues arise for the direct approach when the parameterization has base points. Both direct and our method handle intersection points with high multiplicity easily.

Let us also compare with the matrix representation based on syzygies and employed for ray shooting by [START_REF] Shen | A line/trimmed NURBS surface intersection algorithm using matrix representations[END_REF]. This method is very robust and can handle bicubic surfaces with preprocessing and ray-shooting runtimes both below 1 msec, using a C++ implementation and state-of-the-art libraries such as LAPACK and Eigen. The second runtime may increase in the case of multiple intersection points, but base points pose no limitation. We have focused on methods potentially offering exact solutions if given enough precision; hence, we have not compared to methods relying on a triangular mesh since they are completely different and may need to refine candidate solutions against the parametric patch.

To avoid numerical issues, our approach may employ exact rational computation for the PLUor QR-decompositions. This is feasible with medium-size inputs: e.g., for the Moebious surface, a predicted support of size 21 allows for QR decomposition over the rationals and solving following Lem.2 in 2 sec. Note that this surface has complex base points which make certain real solvers to fail for simple rays. On simpler examples, e.g., the crossed surface

x 1 = t 1 , x 2 = t 2 , x 3 = t 2 1 t 2 2 , with implicit equation x 3 -x 2 1 x 2 2
, the direct method is up to 2 times faster and equally accurate as ours.

Beyond hypersurfaces

This section examines the case of varieties of codimension greater than 1. A fundamental tool in algebraic elimination is the Chow form of a variety.

Definition 1. Let V ⊂ P n be a d-dimensional irreducible variety and H 0 , . . . , H d be linear forms where

H i = u i0 x 0 + • • • + u in x n , i = 0, . . . , d (7)
and u i j are new variables, 0

≤ i ≤ d, 0 ≤ j ≤ n. The Chow form R V of V is a polynomial in the variables u i j such that R V (u i j ) = 0 ⇔ V ∩ {H 0 = 0, . . . , H d = 0} = / 0.
The intersection of the d + 1 hyperplanes H i defined in equation ( 7) is a (nd -1)-dimensional linear subspace L of P n . The set of k-dimensional linear subspaces L of P n is the Grassmannian G(k, n). Those subspaces L that intersect V form a hypersurface Z(V ) on the Grassmannian G(nd -1, n) by [GKZ94, Prop.2.2,p.99]. The Chow form is the unique polynomial, up to a scalar factor, that defines Z(V ). To recover V from Z(V ) we have the following: Proposition 3. [GKZ94, Prop.2.5,p.102] A d-dimensional irreducible subvariety V ⊂ P n is uniquely determined by its Chow form. More precisely, a point ξ ∈ P n lies in V if and only if any (nd -1)dimensional plane containing ξ belongs to the Chow form.

Consider a variety V as in Prop. 3 parameterized as

x i = f i (t), i = 0, . . . , n, t = (t 0 : • • • : t d ), and d + 1 hyperplanes H 0 = • • • = H d = 0
, where H i is defined as in equation ( 7). Substituting x i = f i (t) in every equation H i = 0, we obtain an overdetermined system of equations in the parameters t which we can eliminate by using resultants. This reduces the computation of any Chow form to the computation of a multivariate resultant: Corollary 4. Consider any V ⊂ P n of dimension d, with parameterization x i = f i (t), i = 0, . . . , n. Then, the Chow form R V is the resultant of the hyperplane equations H i , where one eliminates t.

Recall that R V is a polynomial in the Plücker coordinates, i.e., the maximal minors of the (d + 1)×(n+1) matrix whose rows are the normals to the hyperplanes H i . To obtain a representation for V as intersection of implicit hypersurfaces from its Chow form one can apply [GKZ94, Cor.2.6,p.102], or [DS95, Prop.3.1]. Both methods reduce to applying a rewriting algorithm such as the one in [START_REF] Sturmfels | Algorithms in Invariant Theory[END_REF] and typically yield more implicit polynomials than necessary, all of degree equal to the degree of V .

To illustrate the approach in [START_REF] Dalbec | Introduction to chow forms[END_REF], let us focus on varieties of codimension 2 in P 3 , i.e., space curves. Consider the planes H 0 , H 1 in equation (1). The Chow form R V is a polynomial in the brackets [i, j], where [i, j] denotes the maximal minor indexed by the columns 0 ≤ i, j ≤ 3, of the matrix

U := u 0 u 1 u 2 u 3 v 0 v 1 v 2 v 3 .
R V is then rewritten as a polynomial in the dual brackets [[i, j]]. The dual brackets are substituted by the determinant u i v ju j v i of the corresponding minor of U. Finally, the result is expanded as a polynomial whose variables are polynomials in the v 0 , v 1 , v 2 , v 3 and its coefficients are polynomials in the u 0 , u 1 , u 2 , u 3 . The latter polynomials are all of degree equal to the degree of V and form a system of implicit equations of V .

Implicit representations for space curves

We now derive implicit representations by avoiding complex computations, such as the rewriting algorithm. Suppose that we have a space curve parameterized as

x i = f i (t), i = 0, . . . , 3, t = (t 0 : t 1 ).

Let line L be defined by a symbolic point ξ = (ξ 0 : • • • : ξ 3 ) and a sufficiently generic point G ∈ V . Define two planes H 0 , H 1 that intersect along L, by choosing one more random point P 0 and P 1 , respectively. Let their equations be as in (1), where the coefficients now are linear polynomials in ξ .

The Sylvester resultant of this system, where we set x i = f i (t), eliminates t and returns a polynomial in ξ which vanishes on V (but not only on V ), thus offering a necessary but not sufficient condition, see Algorithm 1.

Lemma 5. Let δ = max{deg f i (t), i = 0, . . . , 3} and R G V be the Sylvester resultant of H 0 ( f 0 (t) : . . . : f 3 (t)) and H 1 ( f 0 (t) : . . . : f 3 (t)).

Then R G V factors into a degree δ polynomial defining a surface S G i V containing the space curve V , and a polynomial E δ L , where E L is a linear polynomial defining the plane passing through points G, P 0 , P 1 .

Proof. The degree of the Sylvester resultant in the coefficients of each of the H 0 , H 1 , is δ , and its total degree is 2δ ; ξ is involved (linearly) in the coefficients of both H 0 and H 1 , since it is taken to lie in the intersection of the two planes. Hence the degree of the sought polynomial in ξ is 2δ . It vanishes only in two cases: if ξ belongs to the plane defined by G, P 0 and P 1 , or if ξ belongs to a line passing by G and intersecting V . Hence we can divide (possibly several times) the sought polynomial in ξ by the equation E L of the plane defined by G, P 0 and P 1 , thus obtaining an equation of the conical surface of vertex G and directrix V . Since such a conical surface is of degree δ , its equation is Res(H 0 , H 1 )/E δ L .

Theorem 6. Let p : P 1 (R) → P 3 (R) be a rational parameterization of a space curve V and S G i V , i = 1, 2, 3 be three conical surfaces obtained by the method above with 3 different random points G i / ∈ V . We distinguish two cases.

1. If V is not planar and the points G i are not collinear, then

V = S G 1 V ∩ S G 2 V ∩ S G 3 V .
2. If V is contained in a plane P, at least two of the points G i are not in P and if the G i are not collinear, then V = ∩ i=1,2,3 S G i V .

Proof. Case (1). Since S G i V are three different cones (or cylinders when G i are at infinity), they have no 2-dimensional component in common. We first reduce the problem to the case G 1 = P x := (1 : 0 : 0 : 0), G 2 = P y := (0 : 1 : 0 : 0), G 3 = P z := (0 : 0 : 1 : 0).

We shall prove that an algebraic space curve is the intersection of the 3 cylinders spanned by the curve itself and of directions x 1 , x 2 and x 3 respectively.

Let C be an irreducible component of (S

G 1 V ∩ S G 2 V ∩ S G 3 V ).
Since G 1 , G 2 , G 3 are not collinear, there exists a map φ ∈ PGL(4, R) that sends G 1 to P x 1 , G 2 to P x 2 and G 3 to P x 3 . By linearity of φ (in particular, φ preserves alignment), we have φ (S G 1 V ) = S P x 1 φ (V ) and similar equalities for φ (S G 2 V ) and φ (S G 3 V ). Also, φ (V ) is a rational variety parameterized by φ • p.

Thus, if C ⊂ S G 3 V , then φ (C ) ⊂ S P z φ (V )
. By this argument, we only have to prove that there is no rational space curve φ (V ) for which

S P x 1 φ (V ) ∩ S P x 2 φ (V ) ∩ S P x 3 φ (V ) contains another curve different than φ (V ). For convenience, φ (V ) is denoted by V and φ (C ) is denoted by C in what follows.
Remark 1. Since V does not lie in the plane at infinity, we can switch to the affine setting. We now have "nicer" expressions for our surfaces:

• S P x 1 V = {(x 1 , p 2 (t 1 ), p 3 (t 1 )) | x 1 ,t 1 ∈ R}, • S P x 2 V = {(p 1 (t 1 ), x 2 , p 3 (t 1 )) | x 2 ,t 1 ∈ R}, • S P x 3 V = {(p 1 (t 1 ), p 2 (t 1 ), x 3 ) | x 3 ,t 1 ∈ R}. Since C ⊂ S P x 1
V , there is a (not necessarily rational) parameterization of C given by q : t 1 ∈ R → (q 1 (t 1 ), p 2 (ϕ(t 1 )), p 3 (ϕ(t 1 ))), where q 1 and ϕ are continuous piecewise smooth maps.

Remark 2. We see that ϕ (resp. q 1 ) is not locally constant: otherwise, a part of C would be included in a straight line (resp. a plane), which contradicts the fact that V is not planar.

We can thus pick a small interval I := (a, b) ⊂ R such that q |I is injective and so, without loss of generality, we assume that ϕ |I (t 1 ) = t 1 .

Also, since V is not planar and p is rational, the sets of non regular points of the three maps π x 1 : t 1 → (p 2 (t 1 ), p 3 (t 1 )), π x 2 : t 1 → (p 1 (t 1 ), p 3 (t 1 )) and π x 3 : t 1 → (p 1 (t 1 ), p 2 (t 1 )) are finite. Shrinking I if necessary, we assume there is no such non regular point in I. Now, we have a local curve q(I) = {(q 1 (t 1 ), p 2 (t 1 ), p 3 (t 1 )) | t 1 ∈ I} included in (S

P x 1 V ∩ S P x 2 V ∩ S P x 3 V ).
Using the fact it lies on S P x 2 V , we have another parameterization of q(I) given by r(t 1 ) = (r 1 (t 1 ), p 2 (t 1 ), r 3 (t 1 )) with r i injective on I and r i (I) ⊂ p i (R). Similarly, q(I) ⊂ S P x 3 V gives a third parameterization s(t) = (s 1 (t 1 ), s 2 (t 1 ), p 3 (t 1 )) with s i injective on I and s i (I) ⊂ p i (R).

Comparing s with q and r, we have s(t 1 ) = (s 1 (t 1 ), p 2 (t 1 ), p 3 (t 1 )). Lastly, since π x 1 is regular on I, there is only one branch in π x 1 (s(I)) and so s -1 1 (s(I)) = p -1 1 (s(I)). The curve C is thus locally contained in V ; it follows that C = V .

Case (2). Let us look more closely at Remarks 1 and 2 where we used the non-planarity hypothesis:

1. At Remark 1 we used the fact that V is not contained in the plane containing G 1 , G 2 and G 3 .

2. At Remark 2 we used the fact that the points G i do not belong to a plane containing V , i = 1, . . . , 3.

These facts are still valid under the assumption of G i ∈ P, hence the proof of Case (1) holds.

The general case

Let us generalize the construction above. Let V ⊂ P n be a d-dimensional variety parameterized as

x i = f i (t), i = 0, . . . , n, t = (t 0 : • • • : t d ),
and G be a set of n-d -1 sufficiently generic points G i not in V . Let L be the (n-d -1)-dimensional linear subspace defined by a symbolic point ξ = (ξ 0 : • • • : ξ n ) and the points in G. Chose d + 1 sets of d random points P i = {P i1 , . . . , P id }, P i j / ∈ V . Let H i , i = 0, . . . , d be the hyperplane defined as the span of the points ξ , G, P i . Substitute

x i = f i (t) in each H i to obtain H 0 ( f 0 (t) : . . . : f n (t)) = 0, . . . , H d ( f 0 (t) : . . . : f n (t)) = 0. ( 8
)
The resultant of the polynomials in (8) eliminates t and returns a polynomial R G V in ξ which vanishes on V (but not only on V ), thus offering a necessary but not sufficient condition, see Algorithm 1. To achieve the hypothesis of Prop. 3, one may iterate for a few distinct pointsets G thus obtaining implicit hypersurfaces R G V = 0 whose intersection is V .

To compute the resultant one option is to use interpolation in conjunction with information on its support. This might be obtained from degree bounds on R G V , as explained below, or by the computation of the monomials of the polynomials in (8) using software ResPol from [START_REF] Emiris | An oracle-based, outputsensitive algorithm for projections of resultant polytopes[END_REF]. The latter is analogous to the basic approach for defining the interpolation matrix by the support obtained from the resultant polytope, see Sec. 2.1.

In the general setting we do not have a result similar to Lem. 5. Our method gives rise to extraneous factors, whose expressions and powers are less clear than in the case of space curves. From experiments, it appears that there are only extraneous factors of degree codim(V ). Of course, there are other extraneous factors depending on the chosen method for computing the resultant, e.g., interpolation, or as the determinant of the Macaulay matrix.

To find the degree of R G V in ξ , first note that the degree in t of every H i ( f 0 (t) : . . . : f n (t)) is δ = max{deg t f j (t), j = 0, . . . , n}, and that the coefficients of the H i 's are linear polynomials in ξ . The resultant of these polynomials has degree in the coefficients of each H i equal to δ d , therefore total degree (d + 1)δ d , see e.g. [START_REF] Cox | Using Algebraic Geometry[END_REF]Thm.3.1].

A tighter bound on the degree of R G V can be obtained by considering sparse resultants and mixed volumes [CLO05, Ch.7]. Denote

MV -i = MV (H 0 , . . . , H i-1 , H i+1 , . . . , H d ), 0 ≤ i ≤ d,
the mixed volume of all polynomials excluding H i . The degree of the sparse resultant in the coefficients of H i is known to equal MV -i , therefore its total degree equals ∑ d i=0 MV -i . Both degree bounds above are typically much larger than the degree of V .

An issue arises of course at sampling: all generated points ξ lie on V , whereas we are trying to compute a hypersurface containing V . We expect that the kernel of M will have large dimension and among the kernel vectors we may choose one or more "small" vectors to define the implicit equation. Here "small" may refer to the number of non-zero vector entries, or to the total degree of the monomials corresponding to its non-zero entries.

We can avoid interpolation if there exists a rational or, even better, determinantal formula, for the resultant, such as Sylvester's for curves in P n , or D'Andrea's formula for sparse resultants [D'A02], which are optimal for generic coefficients. For arbitrary coefficients, an infinitesimal perturbation may be applied.

Algorithm 1: Implicit equation Input : V ⊂ P d , parameterized by x i = f i (t), i = 0, . . . , n. Output :

A polynomial vanishing on V .

(1) Define (nd -1)-dimensional linear subspace L by affinely independent random points G = {G 1 , . . . , G n-d-1 } : G i ∈ V , and consider symbolic point ξ = (ξ 0 , . . . , ξ n ).

(2) Define d + 1 hyperplanes H i through G, ξ , P i , for random points P i affinely independent from points in G.

(3) Set x i = f i (t) in the H i : their resultant, where we eliminate t is the sought polynomial in ξ . Compute it by formula or by interpolation as in Sec. 2.1.

Examples and complexity

In this section we illustrate our method through examples of curves in 3D and 4D space and a surface in 4D space.

We switch from the projective to the affine setting and set ξ = (ξ 1 , . . . , ξ n ) ≡ (x 1 , . . . , x n ), for emphasizing these are the implicit variables. For curves in any ambient dimension we compute the resultant directly using the Sylvester matrix, and also by interpolation, employing degree bounds and the expected extraneous linear factor raised to the power 4. This yields a surface containing V but not of minimal degree. As in the previous example, repeating the procedure 3 times, the ideal of the resulting polynomials equals the ideal defined from the polynomials in (10). Note that all iterations produce polynomials of degree 8 which all factor into a degree 4 polynomial vanishing on V and a degree 1 polynomial raised to the power 4, as predicted by Lem. 5.

Alternatively, we interpolate the Sylvester resultant above using as support the lattice points in a 3-simplex of size 8. The 165 × 165 matrix constructed has a kernel of dimension 133. The degrees of the corresponding 133 kernel polynomials vary from 2 to 8. Among them, there are two polynomials of degree 2. They coincide with those in (10).

The equations above were obtained by choosing random points with integral coordinates and relatively close to the curve. When we increase the range or allow for non-integral points, in order to have better chances to avoid the non-generic points, the size of the coefficients increase. For example, using random points with integral coordinates in a box of size 100 around the curve yields equations as below, where we omit most terms:

3983438998535755975578507593x 4 + • • • -6421697880560981054975490304
Example 4. We illustrate Lem. 5 in two more examples. For the space curve of parametric degree 4 in [START_REF] Rueda | An algorithm to parametrize approximately space curves[END_REF], we compute the 3 implicit equations of degree 4 in 0.171 sec. For the Viviani curve which has parametric degree 4 and implicit equations:

x 2 1 + x 2 2 + x 2 3 = 4 a 2 , (x 1 -a) 2 + x 2 2 =
a 2 , we compute 3 implicit equations of total implicit degree 4, and degree 74 in the parameter a, in 7.72 sec.

Example 5. Consider the curve in R 4 with parameterization:

(x 1 , x 2 , x 3 , x 4 ) = (t 2 -t -1, t 3 + 2t 2 -t, t 2 + t -1, t 3 -2t + 3)
The curve is defined by 5 implicit equations of degrees 1,2,2,2 and 3. Our algorithm computes 5 equations of degree 6 in 0.06 sec. These equations contain linear extraneous factors raised to the power 3. When these are divided out we obtain 5 degree 3 equations which define the curve set theoretically.

Example 6. We tested our method on a surface in R 4 with parametric degree 2. It has an implicit representation defined by 7 equations all of degree 3. The resultant of the hyperplane equations is of degree 12, while the polynomial computed as the determinant of the Macaulay matrix is of degree 15 and factors into 4 polynomials: an irreducible polynomial of degree 4 that contains the surface, two polynomials due to the Macaulay resultant, and a power of a linear polynomial due to our method. A total of 5 polynomials define the surface set-theoretically.

In all examples in this section, we used our Maple implementation of Algorithm 1. We compared to Maple's native command algcurves[implicitize] on a laptop with a 2GHz Intel Celeron processor running Maple 18. Our implementation is slower per output equation for curves of degree larger than 6, see Fig. 2. For curves of moderate and high degrees, the Sylvester resultant computation (or interpolation) represents more than 99% of our algorithm's computation time.

Though the runtime per equation is lower, the number of implicit equations that implicitize outputs cannot be predicted and it seems to be only dependent on the input degree: for inputs of the same degree it always returns the same number of implicit equations. For degrees larger than 8, their number is greater than 100. In contrast, our algorithm can return the implicit equations one by one, and Thm. 6 guarantees that 3 of them define the curve set-theoretically. For parameterizations of degree up to 4, implicitize computes a small number of implicit equations whose total degree is lower than the degree of those computed by our method. On the other hand, our algorithm always outputs equations of the same total degree equal to the degree of the parameterization, by Lem. 5. For instance, when the degree of the parameterization is 4, implicitize typically returns 7 equations of degree 3 and 11 equations of degree 4, whereas our algorithm returns 3 equations of degree 4. For Examp. 3, implicitize computes the two equations in (10). This cannot be done by our algorithm unless we know which apex G for the conical surface should be used: for this example these are points at infinity.

The largest example that we computed is of a parametric space curve of degree 18 in less than 3 min.

We also compared our method against implicitization using Groebner basis and the method of matrix representations based on syzygies [BLB12] using the authors' Maple code. For parameterizations of degree up to 6 all methods are comparable in running time. For larger degrees Groebner basis outperform our method by a factor that depends on the degree and the number of monomials in the parametric polynomials and ranges from 10 to 10 3 . The method based on syzygies yields a matrix M which has the drop of rank property and, moreover, its minors of size Rank(M) give a system implicit equations of the variety. However, computing these minors, unlike the computation of M, is very inefficient and it is not in the spirit of this method. On the upside, our algorithm always outputs 3 implicit equations, while in both other methods their number depends on the input and is typically much larger than 3.

Our algorithm was also implemented in Sage with similar runtimes as in Maple.

Future work

Interesting mathematical questions arise in studying the complexity of our method, namely how to bound the degree and the number of computed conical implicit surfaces in the general case by results similar to Thm. 6. This requires to study the extraneous factors of the computed equations and to distinguish them into those coming from the chosen method to compute the resultant, and those inherited in our method. Interpolation matrices extend to other geometric operations such as surface-surface intersection, assuming one surface is either parameterized or given by a point cloud, and the second is parameterized. The first surface is captured by matrix M through a point sample. The parameterization of the second surface is substituted into the last row S(x). Now det M is an implicit equation of the intersection curve, defined in the parametric plane of the second surface. We may avoid developing the determinant but instead we may operate on the matrix. For instance, we can sample the intersection curve by fixing values of one parameter and solving the resulting univariate polynomial det[M |S(t)] in the other parameter t, much like we did with ray shooting.

The interpolation matrices seem to extend beyond the codimension 1 case: given a variety of codimension > 1 (i.e., in (2) we have t = (t 1 , . . . ,t d ) with d ≤ n -2), we can construct a matrix M(x) as in Sect.2.1 that still "represents" the variety, under some assumptions on the support S used: point q ∈ C n lies in the variety iff Rank(M(q)) = Rank(M ). However, in that situation, matrix M is never of full rank. The implicit equations of the variety can be found among the non-zero minors of M(x) of size Rank(M ), but this computation is very expensive. Similar properties hold for the matrix representations in [START_REF] Busé | Implicit matrix representations of rational Bézier curves and surfaces[END_REF] but these matrices are smaller in size.

Figure 1 :

 1 Figure 1: Left: The space curve of Exam. 3. Right: The two implicit surfaces defining the curve.

Figure 2 :

 2 Figure 2: Runtime (sec) per output equation vs parametric degree. Our method's times are in (light) blue on both graphs, Maple's native command's are in (dark) red.

http://sourceforge.net/projects/respol

Maple implementation at http://ergawiki.di.uoa.gr/experiments/simpl.mpl

Acknowledgments

IE, CK, CL belong to team AROMATH, joint between INRIA Sophia-Antipolis and NKUA; their work is partially supported as part of Project ARCADES that has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 675789. We thank L.Busé and J.R.Sendra for discussions.

relying on mixed volume. For d > 1 we use the Macaulay (or the toric resultant) matrix. Interpolating the resultant leads to matrices with very large kernels: on the upside, the polynomials we obtain from the kernel vectors are of degree no greater than those of the former method. Moreover, among them we can find a number of polynomials of small degree and, often, smaller than the degree predicted by degree bounds: these polynomials define the variety set-theoretically, in other words, as the intersection of a few surfaces.

Example 2. Consider the twisted cubic curve affinely parameterized as:

The implicit equations of V are

Let L be the line passing through symbolic point ξ = (x 1 , x 2 , x 3 ), and generic point G ∈ V (here the set G consists of only one point also denoted G). We define two random planes

), intersecting at L by considering additional random points P 1 , P 2 / ∈ L, respectively. Then, the Sylvester resultant of H 1 (t,t 2 ,t 3 ) = H 2 (t,t 2 ,t 3 ) = 0 is a polynomial of degree 6 in ξ which factors into the degree 3 polynomial

and the expected extraneous linear factor raised to the power 3. This yields a surface containing V but not of minimal degree. Repeating the procedure 3 times, the ideal of the resulting polynomials equals the ideal defined from the polynomials in (9).

Alternatively, we may interpolate the Sylvester resultant above. We take as predicted support the lattice points in a 3-simplex of size 6. The 84 × 84 matrix constructed has a kernel of dimension 65. The corresponding 65 kernel polynomials are of degrees from 2 to 6. Among them, the three polynomials of minimal degree 2: x 1 x 3x 2 2 , x 1 x 2x 3 ,x 2 + x 2 1 , define the twisted cubic exactly. In contrast, [DS95, Sec. 3.3] gives 16 (homogeneous) implicit equations, all of degree 3.

Example 3. Consider the space curve in Fig. 1 affinely parameterized as:

It is the intersection of two cylinders:

Let line L be defined from the symbolic point ξ = (x 1 , x 2 , x 3 ), and "generic" point G ∈ V . Define two random planes H 1 and H 2 that intersect at L, by choosing random points P 1 , P 2 / ∈ L, respectively. Then, the Sylvester resultant of

is a polynomial of degree 8 in ξ which factors into the following degree 4 polynomial:

3 + 56x 3 2 + 49x 4 2 + 25x 4 1 -48x 2 3 x 2 -72x 2 3 x 1 + +12x 2 3 x 2 2 + 28x 3 x 3 2 + 12x 2 3 x 2 1 -168x 3 2 x 1 + 222x 2 2 x 2 1 --30x 3 x 3 1 -120x 2 x 3 1 -48x 3 x 2 x 1 -102x 3 x 2 2 x 1 + 76x 3 x 2 1 x 2


Characteristic Constraint Modes for Component Mode Synthesis

Introduction

W HEN modeling the dynamics of a complex structure, it is often impractical to perform a nite element analysis of the entire structure. In some cases, the nite element model (FEM) of the full structure has so many degrees of freedom (DOF) that a global nite element analysisis impossibledue to computermemory constraints or node limits in the software. Also, it is common for component structures to be designed or redesigned independently, which makes it more convenientto perform a separate nite element analysis for each component.

Component mode synthesis (CMS) was developed as a practical and ef cient approach to modeling and analyzing the dynamics of the global structure in such circumstances. [START_REF] Hurty | Dynamic Analysis of Structural Systems Using Component Modes[END_REF][START_REF] Craig | Coupling of Substructures for Dynamic Analyses[END_REF] In CMS, the dynamics of a structure are described by selected sets of normal modes of the individual component structures, plus a set of static vectors that account for the coupling at each interface where component structures are connected. Thus, each componentstructure is dynamically reduced by a separate modal analysis before being coupled at the system level, yielding savings in nite element costs as well as providing signi cant order reduction for the CMS model relative to the FEM of the complex structure. For a review of the CMS literature, the reader is referred to the work by Craig 3 and Seshu. [START_REF] Seshu | Substructuring and Component Mode Synthesis[END_REF] One of the most accurate and widely used CMS methods is the Craig-Bampton method. [START_REF] Craig | Coupling of Substructures for Dynamic Analyses[END_REF][START_REF] Craig | Substructure Methods in Vibration[END_REF] In this method, the component normal modes are calculated with the interface between the component structures held xed. Attachment at the interface is achieved by a set of constraint modes. A constraint mode shape is the static deection induced in the structure by applying a unit displacement to one interface DOF while all of the other interface DOF are held xed. Thus, the motion at the interface is completely described by the constraint modes. However, there is a computational cost associated with these modes because the CMS model must have one DOF for each DOF in the interface of the FEM. This means that, whereas dramatic order reduction may be achieved for the individual components (normal-mode DOF), there is no order reduction for the interface(constraint-modeDOF). In fact, if the nite element mesh is suf ciently ne, the size of the CMS model may be dominated by the constraint-modeDOF. This can make the CMS model cumbersome to use, for example, in analysis of large-scale structures with many components or in probabilistic analysis.

In this paper, a new technique is presented for reducing the size of a CMS model by performing an eigenanalysison the constraintmode partitions of the mass and stiffness matrices. The resultant eigenvectors are called the characteristic constraint ( ) modes. These modes may then be truncated to yield a highly reducedorder model. Therefore, the method presented in this paper offers increased computational ef ciency for performing analyses with a CMS model. Furthermore, the modes evidence the characteristic motion of the interface, so that this method provides signi cant physical insight into the transmission of vibration energy between substructures. In particular, a -mode-based model is well suited to the ef cient calculation of power ow in complex structures, due to the compact representation of interface motion.

To the authors' knowledge, the only similar method is that of Bourquin and d'Hennezel. [START_REF] Bourquin | Intrinsic Component Mode Synthesis and Plate Vibrations[END_REF][START_REF] Bourquin | Numerical Study of an Intrinsic Component Mode Synthesis Method[END_REF] They introduced a xed-interface CMS method "based on a non-conventional choice of constraint modes tied to the normal modes of the Poincaré-Steklov operator associated with the interface between the substructures." [START_REF] Bourquin | Intrinsic Component Mode Synthesis and Plate Vibrations[END_REF] They called these alternative constraint modes the "coupling modes," and they proposed algorithms for handling the various numerical computations needed to nd the coupling modes. [START_REF] Bourquin | Numerical Study of an Intrinsic Component Mode Synthesis Method[END_REF] In contrast, the method proposed here is based on the classic Craig-Bampton method [START_REF] Craig | Coupling of Substructures for Dynamic Analyses[END_REF][START_REF] Craig | Substructure Methods in Vibration[END_REF] and employs a straightforward eigenanalysis to nd the characteristic constraint modes. Therefore, this new method is quite simple, both conceptually and computationally.

This paper is organized as follows. In the second section, the Craig-Bampton method of CMS is brie y reviewed, and the constraint modes are discussed. In the third section, the modes are introduced as a means of reducing the size of the CMS model, and the advantages of a -mode-based representationare covered. For illustration, traditional constraint modes and modes are compared for two different nite element models of a two-span beam on simple supports. In the fourth section, the new CMS method is applied to an example system of a cantilever plate. The mode shapes are shown, and the accuracy in calculating the natural frequencies with various sets of modes is discussed. This work is summarized in the nal section.

Component Mode Synthesis and Constraint Modes

Consider a FEM of a complex structurethat is partitionedinto two substructures, that is, component structures. The boundary where the component structures are connected will be referred to as the interface. The group of DOF in the interface will be denoted by C , and the number of interface DOF by N C . The DOF of substructure 1 (2) that are not in C will be denoted by a (b ), and the number of these DOF will be denoted by N a (N b ). Then the mass and stiffness matrices in nite element coordinatesmay be partitionedas follows:

M m C 0 0 0 m a 0 0 0 m b , K k C k C a k C b k a C k a 0 k b C 0 k b (1)
where subscripts denote the associated DOF and double subscripts indicate coupling terms. The partitioning scheme of Eq. ( 1) is chosen for simplicity of presentation. In general, one may choose to have other active DOF, such as points where forcing is applied or where a damper is attached. These DOF would then be grouped with the interface DOF. Also, CMS formulationsare usually cast in terms of the substructure matrices. However, for the purposesof this study, it is more straightforward to partition the matrices of the full complex structure into interface and noninterface DOF. Furthermore, only two component structures are considered, but this is not restrictive.

The Craig-Bampton method 2 of component mode synthesis employs xed-interface component normal modes, as well as a set of vectors called constraint modes. As mentioned earlier, a constraint mode is the static de ection induced in the structure by applying a unit displacement to one interface DOF while all other interface DOF are held xed. Therefore, the interface partition of the constraint modes is an identity matrix of dimension N C , and the set of constraint modes W C is of the form

W C I W C a W C b (2)
where W C a and W C b are the static shapes induced in substructures 1 and 2, respectively, by the unit displacements at the interface. These shapes are determined by posing the following statics problem:

k C k C a k C b k a C k a 0 k b C 0 k b I W C a W C b f 0 0 ( 3
)
where f is a vector of forces that would be needed to impose the successive unit displacements at the interface DOF. The two sets of homogeneous equations in Eq. ( 3) yield

W C a (k a ) 1 k a C , W C b (k b ) 1 k b C (4)
Using these solutions, the constraint modes of Eq. ( 2) are known.

Next, the component normal modes are obtained by performing a modal for each component structure. Because the Craig-Bampton method uses xed-interface modes, they are determined from the following eigenvalue problems:

k a Á 1 k m a Á 1 , k b Á 2 k m b Á 2 (5)
In general, these modes are truncatedto reduce the size of the model.

For example, if U i is the full eigenvector matrix for substructurei and N i modes are selected,

U 1 Na Na select modes U 1 Na N1 , U 2
Nb Nb select modes U 2 Nb N2 (6)
Then, by the use of the constraint modes and selected sets of component normal modes, the transformation from nite element coordinates x to CMS generalized coordinates y is

x I 0 0 W C a U 1 0 W C b 0 U 2 y ( 7
)
Applying the transformationof Eq. ( 7), the global mass and stiffness matrices for the CMS model have the form 3

MCMS m

C m C 1 m C 2 m 1C m1 0 m 2C 0 m2 , KCMS k C 0 0 0 k1 0 0 0 k2 (8)
Note that this transformation yields only inertial coupling between the normal modes and the constraint modes.

The number of normal-mode DOF depends, of course, on how many normal modes are selected for each component structure. However, the number of constraint-mode DOF is equal to N C , the number of FEM DOF that are in the interface. Thus, the size of the constraint-mode partitions is determined by the nite element mesh. If there are many nite element nodes in the interface region, then the constraint-mode partitions of the CMS matrices may be relatively large. Furthermore, as the number of constraint modes increases, the individual constraint modes may have limited physical meaning.

To illustrate the last point, consider the simply supported twospan beam presented in Fig. 1. This system is treated as a complex structure consisting of two component structures with the interface at the center support. Each span is treated as a substructure,referred to as beam 1 and beam 2. (The modeling of power ow in this system was examined in detail in an earlier work by the authors. [START_REF] Tan | Approximations of Power Flow Between Two Coupled Beams Using Statistical Energy Methods[END_REF] ) Two different nite element models were generated for this system: one consisting of Timoshenko beam elements and a second, more generalFEM consistingof solid elements.All nite element analysis in this study was performed using the commercial code ABAQUS.

For the beam-element FEM, the displacement at the center support is constrained by the pin, so that there is only one (rotational) DOF at the interface. Thus, there is only one constraint mode in the corresponding CMS model. This single constraint mode, which is obtained by imposing a static angular displacement about the center support, is shown in Fig. 2. Only beam 1 is shown. The interface is at the right end.

In contrast, the solid-element FEM has 25 nodes at the interface. Note that one horizontalline of nodes in the center of the interface is xed to represent a pinned condition. This leaves 20 unconstrained nite element nodes (60 DOF) in the interface of the solid-element FEM, leading to 60 constraint modes in the CMS model. One of these constraintmodes is shown in Fig. 3. Only beam 1 is shown. The interface node at which a unit displacement is imposed is marked by a dot. Signi cant local deformation is evident in this constraint mode. It can also be seen that the offset of the imposed displacement from the pinned line of nodes leads to bending of the beam in the XY plane. It is easy to imagine that a linear combination of similar constraint modes could yield a more natural exural motion, such as that of the lone constraint mode for the beam-element FEM. Furthermore, a complete set of these linear combinations could be de ned as a new set of constraint modes that describe the principal ways in which the interfacedeforms.This idea is consideredin detail in the next section.

Modes

It is now suggestedthat the number of necessaryconstraintmodes may be reduced by seeking a new set of modes that correspond to more natural physical motion at the interface. This is posed as an eigenvalue problem for the constraint-mode partitions of the CMS matrices:

k C Á k m C Á (9)
Suppose that these eigenvectorsare then truncated,as in a traditional modal analysis:

U NC NC select modes

U

NC N (10) This selected set of eigenvectorsmay be used to transform the mass and stiffness matrices to yield a reduced-order model (ROM). The The mass and stiffness matrices in ROM generalized coordinates are, thus,

MROM U T m C U U T m C 1 U T m C 2 m 1C U m1 0 m 2C U 0 m2 m m 1 m 2 m1 m1 0 m2 0 m2 (12) K ROM U T k C U 0 0 0 k1 0 0 0 k2 k 0 0 0 k1 0 0 0 k2 (13)
Now, compared to the mass and stiffness matrices of Eq. ( 1), the size of every matrix partition has been reduced. Finally, an eigenvector from Eq. ( 9) may be transformed from constraint-mode coordinates into nite element coordinates:

à W C Á (14)
Note that à is an eigenvector-based linear combination of constraint modes. Therefore, it captures some characteristic physical motion in the interface region, and accordingly, it is called a characteristic constraint ( ) mode. The modes provide, at least in an approximate sense, the principal modes of deformation for the interface. In fact, the calculation and selection of the modes is essentially a secondary modal analysis, and the bene ts are similar to those of a traditional modal analysis. For example, a re nement of the nite element mesh would lead to an increase in the accuracy of the ROM without requiring an increase in the size of the ROM because the same number of modes could be retained. Note that, in certain cases, it might be more bene cial to perform an eigenanalysison the full matrices of the CMS model. This would yield a set of system modes for the complex structure. This idea has been investigatedby the authorsin the contextof modelingthe vibration of turbomachinery rotors, 9,10 but it will not be considered here.

To illustrate the characteristic motion that is captured by the modes, the rst mode for the solid-element FEM of the beam system is shown in Fig. 4b. Only beam 1 is shown. The interface is at the right end. Note that it appears to be identical to the single constraint mode for the beam-element FEM, which is shown again in Fig. 4a for comparison. It is clear that the rst mode captures the XY-plane exural motion at the interface. The second mode for the solid-elementFEM is shown in Fig. 5. It is apparent that the second mode corresponds to exural motion in the XZ plane. Only beam 1 is shown. Figure 5 is a rotated view with the interface shown in the foreground.

Thus, depending on the type of motion and frequency range of interest, the modes may be truncated to yield a new CMS model with a greatly reduced number of DOF. The ef ciency and accuracy of this method will be investigated in the next section.

Example

The cantilever plate shown in Fig. 6 is now considered as an example structure. The dimensions and material properties are chosen to be the same as those used by Craig and Bampton 2 : 2024-T3 aluminum,Young's modulus E 10.5 10 6 psi (72 GPa), Poisson's ratio m 0.33, density q 0.101 lb/in. [START_REF] Craig | Structural Dynamics: An Introduction to Computer Methods[END_REF] (2800 kg/m 3 ), and thickness t 0.125 in. (0.3175 cm). The plate is considered to consist of two substructures as shown. This partitioning is also the same as that of Ref. 2. A FEM of this plate was constructed using a 48 24 1 mesh of ABAQUS shell elements (element S4R, four nodes/element, six DOF/node). The analyticalresultsfor the rst 10 naturalfrequencies of the cantilever plate are given by Gorman, 11 and these were used as a benchmark. The nite element mesh is suf ciently ne that the rst 10 FEM natural frequenciesare all within 0.7% of the analytical results.

A Craig-Bampton CMS model 2 of the plate was then generated. The rst 10 normal modes were retained for each component structure. In addition, there were 144 constraint modes, yielding a CMS model with 164 DOF. For reference, the rst 15 natural of this CMS model are all within 0.9% of the corresponding FEM natural frequencies.

Next, modes were calculated, and 10 different ROMs were constructed by keeping only the lowest one to 10 modes. Because a total of 20 component normal modes were selected for the CMS model, the smallest ROM considered contains 21 DOF and the largest contains 30 DOF.

The error for the rst 15 natural frequenciesof the cantileverplate is shown vs the number of selected modes in Fig. 7. This error Fig. 8 First six characteristic constraint modes for the cantilever plate.

is calculated relative to the 164-DOF CMS model because a mode-based model will be no more accurate than its parent CMS model. In fact, if all of the modes are included,no approximation is made. Note that only 10 modes are needed to drive the error below 0.01% for the rst 15 natural frequencies of the cantilever plate. Therefore, using the present method, a 30-DOF CMS model is highly accurate relative to the 164-DOF CMS model for the frequency range considered. Note that if the FEM mesh were re ned, the full CMS model size would increase. Yet a ROM calculated using the selected set of modes would become more accurate without a penalty of any additional DOF.

Insight into the results presented in Fig. 7 may be gained by comparing the modes to the global modes of the cantileverplate. The rst six modes for plate 2 are shown in Fig. 8. (Only plate 2 is shown. Views are rotated such that the interface is seen at the lower left edge, along the y axis.) Consider, for example, the second mode. It is clear that the second mode captures much of the interface-inducedmotion seen in the second and third global modes, which are shown in Fig. 9. This explains why the frequency error decreases by more than two orders of magnitude for these global modes (Fig. 7) when mode 2 is included in the ROM. In Fig. 9, for reference, the interface between plates 1 and 2 is marked by a dark line. Returning to Fig. 8, note that the wavelength of the interface motion decreases with higher modes, as would be expected. Furthermore, there is an attendant decrease in the in uence of the interface motion at the far end of the component structure. Thus, it appearsthat the lower modes capture most of the global in uence of the interface motion. Higher modes, then, would be needed only to model the details of the local interface motion.

Conclusions

This paper introduceda techniquefor reducingthe size of a Craig-Bampton CMS 2,5 model by improving the representation of the interface between component structures. This new method is based on computing modes. The modes are found by performing an eigenanalysis of the partitions of the CMS mass and stiffness matrices that correspond to the Craig-Bampton constraint modes.

It was seen that the modes capture the characteristic motion of the interface, and thus, they may be truncated as if they were traditional modes of vibration. This truncation leads to a highly reduced-orderCMS model. In addition, the modes allow signicant insight into the physical mechanisms of vibration transmission between the component structures. This information could be used, for example, to determine the design parameters that have a critical impact on power ow.

The calculation and selection of the modes is essentially a secondary modal analysis. Therefore, the bene ts are the same as those of a traditional modal analysis. For instance, re ning thenite element mesh would increase the accuracy of a ROM without introducing any additional DOF.

In general, this investigation showed that the use of modes in CMS allows highly ef cient modeling of the dynamics of complex structures.Also, this techniqueis capable of providingphysical insight into how the vibration energy ows through the structure. Whereas it is believed that this method has wide applicability, it seems especially suited to predicting power ow in complex structures.

  Matthew P. Castanier, Yung-Chang Tan, and Christophe Pierre University of Michigan, Ann Arbor, Michigan 48109-2125 A technique is presented for reducing the size of a model generated by the Craig-Bampton method (Craig, R. R., and Bampton, M. C. C., "Coupling of Substructures for Dynamic Analyses," AIAA Journal, Vol. 6, No. 7, 1968, pp. 1313-1319) of component mode synthesis (CMS). An eigenanalysis is performed on the partitions of the CMS mass and stiffness matrices that correspond to the so-called constraint modes. The resultant eigenvectors are referred to as characteristic constraint modes because they represent the characteristic motion of the interface between the component structures. When the characteristic constraint modes are truncated, a CMS model with a highly reduced number of degrees of freedom may be obtained. An example of a cantilever plate is considered. It is shown that relatively few characteristic constraint modes are needed to yield accurate approximationsof the lower natural frequencies. Furthermore, this method yields physical insightinto the mechanisms of vibration transmission in complex structures, and it provides an excellent framework for the ef cient calculation of power ow. Nomenclature K = stiffness matrix k = stiffness matrix partition k = stiffness matrix partition that has been reduced in both dimensions (rows and columns) by coordinate transformations using selected sets of modes M = mass matrix m = mass matrix partition m = mass matrix partition that has been reduced in both dimensions (rows and columns) by coordinate transformations using selected sets of modes N = number of physical or modal degrees of freedom (DOF) x = displacements in physical coordinates, for example, nite element DOF y = generalized coordinates of the component mode synthesis model z = generalized coordinates of the reduced-order model k = eigenvalue U = full set of characteristic constraint eigenvectors U = selected set of characteristic constraint eigenvectors U i = full set of normal modes for substructure i U i = selected set of normal modes for substructure i Á = characteristic constraint eigenvector (in constraint-mode DOF) Á i = normal mode for substructure i W C = full set of Craig-Bampton constraint modes 2 Ã = characteristic constraint mode (in physical coordinates) Subscripts C = constraint-mode DOF = characteristic-constraint-mode DOF CMS = component mode synthesis model ROM = reduced-ordermodel a = physical DOF of substructure 1, excluding the interface DOF a C , b C , = coupling terms between physical DOF C a , C b of a substructure and the interface b = physical DOF of substructure 2, excluding the interface DOF C = physical DOF of the interface between substructures 1 = normal-mode DOF of substructure 1 2 = normal-mode DOF of substructure 2 1C, 2C, = coupling terms between normal-mode C1, C 2 and constraint-mode DOF

Fig. 1 Fig. 2

 12 Fig. 1 Two-span beam on simple supports.

Fig. 3

 3 Fig. 3 One of the 60 constraint modes for the solid-element FEM: a) side view of beam 1, as in Fig. 1, with the interface at the right end and b) rotated view of beam 1 with the interface shown in the foreground.

Fig. 4 Fig. 5

 45 Fig. 4 Comparison of a) single constraint mode for the beam-element FEM and b) rst mode for the solid-element FEM.

Fig. 6 Fig. 7

 67 Fig. 6 Cantilever plate partitioned into two component structures, referred to as plate 1 and plate 2.

Fig. 9

 9 Fig. 9 Second and third global modes for the cantilever plate.

,

, = coupling terms between normal-mode 1, 2 and characteristic-constraint-mode DOF