Pascal Ochem
email: ochem@lirmm.fr

Matthieu Rosenfeld
email: matthieu.rosenfeld@ens-lyon.fr

Avoidability of Formulas with two Variables

Keywords: Word, Pattern avoidance, -AAB, ABA, ABB, BBA, BAB

In combinatorics on words, a word w over an alphabet Σ is said to avoid a pattern p over an alphabet ∆ of variables if there is no factor f of w such that f = h(p) where h : ∆ * → Σ * is a non-erasing morphism. A pattern p is said to be k-avoidable if there exists an infinite word over a k-letter alphabet that avoids p. We consider the patterns such that at most two variables appear at least twice, or equivalently, the formulas with at most two variables. For each such formula, we determine whether it is 2-avoidable.

Introduction

A pattern p is a non-empty finite word over an alphabet ∆ = {A, B, C, . . .} of capital letters called variables. An occurrence of p in a word w is a nonerasing morphism h : ∆ * → Σ * such that h(p) is a factor of w. The avoidability index λ(p) of a pattern p is the size of the smallest alphabet Σ such that there exists an infinite word over Σ containing no occurrence of p. Bean, Ehrenfeucht, and McNulty [START_REF] Bean | Avoidable patterns in strings of symbols[END_REF] and Zimin [START_REF] Zimin | Blocking sets of terms[END_REF] characterized unavoidable patterns, i.e., such that λ(p) = ∞. We say that a pattern p is t-avoidable if λ(p) t. For more informations on pattern avoidability, we refer to Chapter 3 of Lothaire's book [START_REF] Lothaire | Algebraic Combinatorics on Words[END_REF].

A variable that appears only once in a pattern is said to be isolated. Following Cassaigne [START_REF] Cassaigne | Motifs évitables et régularité dans les mots[END_REF], we associate to a pattern p the formula f obtained by replacing every isolated variable in p by a dot. The factors between the dots are called fragments.

An occurrence of f in a word w is a non-erasing morphism h : ∆ * → Σ * such that the h-image of every fragment of f is a factor of w. As for patterns, the avoidability index λ(f ) of a formula f is the size of the smallest alphabet allowing an infinite word containing no occurrence of p. Clearly, every word avoiding f also avoids p, so λ(p) λ(f ). Recall that an infinite word is recurrent if every finite factor appears infinitely many times. If there exists an infinite word over Σ avoiding p, then there there exists an infinite recurrent word over Σ avoiding p. This recurrent word also avoids f , so that λ(p) = λ(f ). Without loss of generality, a formula is such that no variable is isolated and no fragment is a factor of another fragment.

Cassaigne [START_REF] Cassaigne | Motifs évitables et régularité dans les mots[END_REF] began and Ochem [START_REF] Ochem | A generator of morphisms for infinite words[END_REF] finished the determination of the avoidability index of every pattern with at most 3 variables. A doubled pattern contains every variable at least twice. Thus, a doubled pattern is a formula with exactly one fragment. Every doubled pattern is 3-avoidable [START_REF] Ochem | Doubled patterns are 3-avoidable[END_REF]. A formula is said to be binary if it has at most 2 variables. In this paper, we determine the avoidability index of every binary formula.

We say that a formula f is divisible by a formula f if f does not avoid f , that is, there is a non-erasing morphism such that the image of any fragment of f by h is a factor of a fragment of f . If f is divisible by f , then every word avoiding f also avoids f and thus λ(f ) λ(f ). For example, the fact that ABA.AABB is 2-avoidable implies that ABAABB and ABAB.BBAA are 2avoidable. Moreover, the reverse f R of a formula f satisfies λ(f R ) = λ(f ). See Cassaigne [START_REF] Cassaigne | Motifs évitables et régularité dans les mots[END_REF] and Clark [START_REF] Clark | Avoidable formulas in combinatorics on words[END_REF] for more information on formulas and divisibility.

First, we check that every avoidable binary formula is 3-avoidable. Since λ(AA) = 3, every formula containing a square is 3-avoidable. Then, the only square free avoidable binary formula is ABA.BAB with avoidability index 3 [START_REF] Cassaigne | Motifs évitables et régularité dans les mots[END_REF]. Thus, we have to distinguish between avoidable binary formulas with avoidability index 2 and 3. A binary formula is minimally 2-avoidable if it is 2-avoidable and is not divisible by any other 2-avoidable binary formula. A binary formula f is maximally 2-unavoidable if it is 2-unavoidable and every other binary formula that is divisible by f is 2-avoidable.

Theorem 1. Up to symmetry, the maximally 2-unavoidable binary formulas are:

-AABA.BABA -AAA -ABA.BAAB.BAB -AABA.ABAA.BAB -AABA.ABAA.BAAB -ABAAB
To obtain the 2-unavoidability of the formulas in the first part of Theorem 1, we use a standard backtracking algorithm. In the rest of the paper, we consider the 2-avoidable formulas in the second part of Theorem 1. Fig. 1 gives the maximal length and number of binary words avoiding each maximally 2-unavoidable formula.

We show in Section 3 that the first three of these formulas are avoided by polynomially many binary words only. The proof uses a technical lemma given in Section 2. Then we show in Section 4 that the other formulas are avoided by exponentially many binary words.

The Useful Lemma

Let us define the following words:

-b 2 is the fixed point of 0 → 01, 1 → 10. -b 3 is the fixed point of 0 → 012, 1 → 02, 2 → 1. -b 4 is the fixed point of 0 → 01, 1 → 03, 2 → 21, 3 → 23. -b 5 is the fixed point of 0 → 01, 1 → 23, 2 → 4, 3 → 21, 4 → 0.
Let w and w be infinite (right infinite or bi-infinite) words. We say that w and w are equivalent if they have the same set of finite factors. We write w ∼ w if w and w are equivalent. A famous result of Thue [START_REF] Thue | Über unendliche Zeichenreihen[END_REF] can be stated as follows: Theorem 2. [START_REF] Thue | Über unendliche Zeichenreihen[END_REF] Every bi-infinite ternary word avoiding 010, 212, and squares is equivalent to b 3 .

Given an alphabet Σ and forbidden structures S, we say that a finite set W of infinite words over Σ essentially avoids S if every word in W avoids S and every bi-infinite words over Σ avoiding S is equivalent to one of the words in S.

If W contains only one word w, we denote the set W by w instead of {w}. Then we can restate Theorem 2: b 3 essentially avoids 010, 212, and squares

The results in the next section involve b 3 . We have tried without success to prove them by using Theorem 2. We need the following stronger property of b 3 : Then there is no possible suffix letter for P : 0 gives 010, 1 gives 11, and 2 gives 212. This rules out the case N = P 1. So N = P 0 and w contains 2P 02P 0. This factor cannot extend to 02P 02P 0, since this contains the large square 02P 02P and this case is covered by the first item. Thus w contains 12P 02P 0. If P = Q1, then w contains 12Q102Q10. Since {22, 1021} ⊂ F , the factor 102Q implies that Q = 0R, so that w contains 120R1020R10. Then there is no possible prefix letter for R: 0 gives 00, 1 gives 1201, and 2 gives 0202. This rules out the case P = Q1. So P = Q2 and w contains 12Q202Q20. The factor Q202 implies that Q = R1 and w contains 12R1202R120. Since {00, 1201} ⊂ F , w contains 12R1202R1202, which contains 2Y Y with Y = R1202.

Formulas Avoided by Few Binary Words

The first three 2-avoidable formulas in Theorem 1 are not avoided by exponentially many binary words:

-{g x (b 3 ), g y (b 3 ), g z (b 3 ), g z (b 3 )} essentially avoids AA.ABA.ABBA.

-{g x (b 3 ), g t (b 3 )} essentially avoids ABA.AABB.

g x (b 3 ) essentially avoids AABA.ABB.BBA.

The words avoiding these formulas are morphic images of b 3 by the morphisms given below. Let w denote the word obtained from the (finite or bi-infinite) binary word w by exchanging 0 and 1. Obviously, if w avoids a given formula, then so does w. A (bi-infinite) binary word w is self-complementary if w ∼ w.

The words g x (b 3 ), g y (b 3 ), and g t (b 3 ) are self-complementary. Since the frequency of 0 in g z (b 3 ) is 5 9 , g z (b 3 ) is not self-complementary. Then g z is obtained from g z by exchanging 0 and 1, so that g z (b 3 ) = g z (b 3 ).

g x (0) = 01110, g x (1) = 0110, g x (2) = 0. g y (0) = 0111, g y (1) = 01, g y (2) = 00. g z (0) = 0001, g z (1) = 001, g z (2) = 11. g t (0) = 01011011010, g t (1) = 01011010, g t (2) = 010.
To prove the avoidability, we have implemented Cassaigne's algorithm that decides, under mild assumptions, whether a morphic word avoids a formula [START_REF] Cassaigne | Motifs évitables et régularité dans les mots[END_REF]. For the first two formulas, we have to explain how the long enough binary words split into 4 or 2 distinct incompatible types. A similar phenomenon has been described for AABB.ABBA [START_REF] Ochem | Binary words avoiding the pattern AABBCABBA[END_REF].

First, consider any infinite binary word w avoiding AA.ABA.ABBA. A computer check shows by backtracking that w must contain the factor 01110001110. In particular, w contains 00. Thus, w cannot contain both 010 and 0110, since it would produce an occurrence of AA.ABA.ABBA. Moreover, a computer check shows by backtracking that w cannot avoid both 010 and 0110. So, w must contain either 010 or 0110 (this is an exclusive or). Similarly, w must contain either 101 or 1001. There are thus at most 4 possibilities for w, depending on which subset of {010, 0110, 101, 1001} appears among the factors of w, see Figure 2a. Now, consider any infinite binary word w avoiding ABA.AABB. Notice that w cannot contain both 010 and 0011. Also, a computer check shows by backtracking that w cannot avoid both 010 and 1100. By symmetry, there are thus at most 2 possibilities for w, depending on which subset of {010, 0011, 101, 1100} appears among the factors of w, see Figure 2b. Let us first prove that g y (b 3 ) essentially avoids AA.ABA.ABBA, 0110, and 1001. We check that the set of prolongable binary words of length 100 avoiding AA.ABA.ABBA, 0110, and 1001 is exactly the set of factors of length 100 of g y (b 3 ). Using Cassaigne's notion of circular morphism [START_REF] Cassaigne | Motifs évitables et régularité dans les mots[END_REF], this is sufficient to prove that every bi-infinite binary word of this type is the g y -image of some biinfinite ternary word w 3 . It also ensures that w 3 and b 3 have the same set of small factors. Suppose for contradiction that w 3 = b 3 . By Lemma 3, w 3 contains 2Y Y . Then w 3 contains 2Y Y a with a ∈ Σ 3 . Notice that 0 is a prefix of the g y -image of every letter. So g y (w 3 ) contains g y (2Y Y a) = 000U 0U 0V with U, V ∈ Σ + 3 , which contains an occurrence of AA.ABA.ABBA with A = 0 and B = 0U . This shows that w 3 ∼ b 3 , and thus g y (w 3 ) ∼ g y (b 3 ). Thus g y (b 3 ) essentially avoids AA.ABA.ABBA, 0110, and 1001. The argument is similar for the other types and we only detail the final contradiction:

-Since 1 is a suffix of the g z -image of every letter, g z (2Y Y ) = 11U 1U 1 contains an occurrence of AA.ABA.ABBA with A = 1 and B = 1U . -Since 010 is a prefix and a suffix of the g t -image of every letter, g t (u2Y Y ) = V 010010010U 010010U 010 contains an occurrence of ABA.AABB with A = 010 and B = 010U 010. -Since 0 is a prefix and a suffix of the g x -image of every letter, g x (u2Y Y a) = V 000U 00U 00W contains an occurrence of AABA.AABBA with A = 0 and B = 0U 0. Therefore, g x (u2Y Y a) contains an occurrence of AA.ABA.ABBA, ABA.AABB, and AABA.ABB.BBA.

Formulas Avoided by Exponentially Many Binary Words

The other 2-avoidable formulas in Theorem 1 are avoided by exponentially many binary words. For every such formula f , we give below a uniform morphism g that maps every ternary square free word to a binary word avoiding f . If possible, we simultaneously avoid the reverse formula f R of f . We also avoid large squares. Let SQ t denote the pattern corresponding to squares of period at least t, that is, SQ 1 = AA, SQ 2 = ABAB, SQ 3 = ABCABC, and so on. The morphism g produces words avoiding SQ t with t as small as possible.

-AA.ABA.BABB is avoided with its reverse by the following 22-uniform morphism which also avoids SQ 6 : Previous papers [START_REF] Ochem | A generator of morphisms for infinite words[END_REF][START_REF] Ochem | Binary words avoiding the pattern AABBCABBA[END_REF] have considered a 102-uniform morphism to avoid AABB.ABBA and SQ 27 . No infinite binary word avoids AABB.ABBA and SQ 15 . -ABAB.BABA is self-reverse, 50-uniform morphism, avoids SQ 3 , see [START_REF] Ochem | A generator of morphisms for infinite words[END_REF]:

0 → 0001101101110011100011 1 →
0 → 00011001011000111001011001110001011100101100010111 1 → 00011001011000101110010110011100010110001110010111 2 → 00011001011000101110010110001110010111000101100111
Notice that a binary word avoiding ABAB.BABA and SQ 3 contains only the squares 00, 11, and 0101 (or 00, 11, and 1010).

-AABA.BABA: A case analysis of the small factors shows that a recurrent binary word avoids AABA.BABA, ABAA.ABAB, and SQ 3 if and only if it contains only the squares 00, 11, and 0101 (or 00, 11, and 1010). We thus obtain the same morphism as for ABAB.BABA. -AAA is self-reverse, 32-uniform morphism, avoids SQ 4 :

0 → 00101001101101001011001001101011 1 → 00101001101100101101001001101011 2 → 00100101101001001101101001011011
-ABA.BAAB.BAB is self-reverse, 10-uniform morphism, avoids SQ 3 :

0 → 0001110101 1 → 0001011101 2 → 0001010111
-AABA.ABAA.BAB is self-reverse, 57-uniform morphism, avoids SQ 6 :

0 → 000101011100010110010101100010111001011000101011100101011 1 → 000101011100010110010101100010101110010110001011100101011 2 → 000101011100010110010101100010101110010101100010111001011
-AABA.ABAA.BAAB is self-reverse, 30-uniform morphism, avoids SQ 3 :

0 → 000101110001110101000101011101 1 → 000101110001110100010101110101 2 → 000101110001010111010100011101
-ABAAB is avoided with its reverse, 10-uniform morphism, avoids SQ 3 , see [START_REF] Ochem | A generator of morphisms for infinite words[END_REF]:

0 → 0001110101 1 → 0000111101 2 → 0000101111
For every q-uniform morphism g above, we say that a binary word is an sqf-g-image if it is the g-image of a ternary square free word. Let us show that for every minimally 2-avoidable formula f and corresponding morphism g, every sqf-g-image avoids f . We start by checking that every morphism is synchronizing, that is, for every letters a, b, c ∈ Σ 3 , the factor g(a) only appears as a prefix or a suffix in g(bc).

For every morphism g, the sqf-g-images are claimed to avoid SQ t with 2t < q. Let us prove that SQ t is avoided. We first check exhaustively that the sqf-gimages contain no square uu such that t |u| < 2q -1. Now suppose for contradiction that an sqf-g-image contains a square uu with |u| 2q -1. The condition |u| 2q -1 implies that u contains a factor g(a) with a ∈ Σ 3 . This factor g(a) only appears as the g-image of the letter a because g is synchronizing. Thus the distance between any two factors u in an sqf-g-image is a multiple of q. Since uu is a factor of an sqf-g-image, we have q | |u|. Also, the center of the square uu cannot lie between the g-images of two consecutive letters, since otherwise there would be a square in the pre-image. The only remaining possibility is that the ternary square free word contains a factor aXbXc with a, b, c ∈ Σ 3 and X ∈ Σ + 3 such that g(aXbXc) = bsY psY pe contains the square uu = sY psY p, where g(X) = Y , g(a) = bs, g(b) = ps, g(c) = pe. Then, we also have a = b and b = c since aXbXc is square free. Then abc is square free and g(abc) = bspspe contains a square with period |s| + |p| = |g(a)| = q. This is a contradiction since the sqf-g-images contain no square with period q.

Notice that f is not square free, since the only avoidable square free binary formula is ABA.BAB, which is not 2-avoidable. Now, we distinguish two kinds of formula. A formula is easy if every appearing variable is contained in at least one square. Every potential occurrence of an easy formula then satisfies |A| < t and |B| < t since SQ t is avoided. The longest fragment of every easy formula has length 4. So, to check that the sqf-g-images avoids an easy formula, it is sufficient to consider the set of factors of the sqf-g-images with length at most 4(t -1).

A tough formula is such that one of the variables is not contained in any square. The tough formulas have been named so that this variable is B. The tough formulas are ABA.BAAB.BAB, ABAAB, AABA.ABAA.BAAB, and AABA.ABAA.BAB. As before, every potential occurrence of a tough formula satisfies |A| < t since SQ t is avoided. Suppose for contradiction that |B| 2q-1. By previous discussion, the distance between any two occurrences of B in an sqf-g-image is a multiple of q. The case of ABA.BAAB.BAB can be settled as follows. The factor BAAB implies that q | |BAA| and the factor BAB implies that q | |BA|. This implies that q | |A|, which contradicts |A| < t. For the other formulas, only one fragment contains B twice. This fragment is said to be important. Since |A| < t, the important fragment is a repetition which is "almost" a square. The important fragment is BAB for AABA.ABAA.BAB, BAAB for AABA.ABAA.BAAB, and ABAAB for ABAAB. Informally, this almost square implies a factor aXbXc in the ternary pre-image, such that |a| = |c| = 1 and 1 |b| 2. If |X| is small, then |B| is small and we check exhaustively that there exists no small occurrence of f . If |X| is large, there would exist a ternary square free factor aY bY c with |Y | small, such that g(aY bY c) contains the important fragment of an occurrence of f if and only if g(aXbXc) contains the important fragment of a smaller occurrence of f .

Concluding Remarks

From our results, every minimally 2-avoidable binary formula, and thus every 2-avoidable binary formula, is avoided by some morphic image of b 3 .

What can we forbid so that there exists only few infinite avoiding words ? The known examples from the literature [START_REF] Badkobeh | Characterization of some binary words with few squares[END_REF][START_REF] Baker | Growth problems for avoidable words[END_REF][START_REF] Thue | Über unendliche Zeichenreihen[END_REF]

Now we can extend this list:

one formula over two variables:

• g x (b 3 ) essentially avoids AAB.BAA.BBAB.

• {g x (b 3 ), g t (b 3 )} essentially avoids ABA.AABB.

• {g x (b 3 ), g y (b 3 ), g z (b 3 ), g z (b 3 )} essentially avoids AA.ABA.ABBA. one pattern over three variables: ABACAABB (same as ABA.AABB).

Fig. 1 :

 1 Fig.1:The number and maximal length of binary words avoiding the maximally 2-unavoidable formulas.

Lemma 3 .

 3 b 3 essentially avoids 010, 212, XX with 1 |X| 3, and 2Y Y with |Y | 4. Proof. We start by checking by computer that b 3 has the same set of factors of length 100 as every bi-infinite ternary word avoiding 010, 212, XX with 1 |X| 3, and 2Y Y with |Y | 4. The set of the forbidden factors of b 3 of length at most 4 is F = {00, 11, 22, 010, 212, 0202, 2020, 1021, 1201}. To finish the proof, we use Theorem 2 and we suppose for contradiction that w is a biinfinite ternary word that contains a large square M M and avoids both F and large factors of the form 2Y Y . -Case M = 0N . Then w contains M M = 0N 0N . Since 00 ∈ F and 2Y Y is forbidden, w contains 10N 0N . Since {11, 010} ⊂ F , w contains 210N 0N . If N = P 1, then w contains 210P 10P 1, which contains 2Y Y with Y = 10P . So N = P 2 and w contains 210P 20P 2. If P = Q1, then w contains 210Q120Q12. Since {11, 212} ⊂ F , the factor Q12 implies that Q = R0 and w contains 210R0120R012. Moreover, since {00, 1201} ⊂ F , the factor 120R implies that R = 2S and w contains 2102S01202S012. Then there is no possible prefix letter for S: 0 gives 2020, 1 gives 1021, and 2 gives 22. This rules out the case P = Q1. So P = Q0 and w contains 210Q020Q02. The factor Q020Q implies that Q = 1R1, so that w contains 2101R10201R102. Since {11, 010} ⊂ F , the factor 01R implies that R = 2S, so that w contains 21012S102012S102. The only possible right extension with respect to F of 102 is 102012. So w contains 21012S102012S102012, which contains 2Y Y with Y = S102012. -Case M = 1N . Then w contains M M = 1N 1N . In order to avoid 11 and 2Y Y , w must contain 01N 1N . If N = P 0, then w contains 01P 01P 0. So w contains the large square 01P 01P and this case is covered by the previous item. So N = P 2 and w contains 01P 21P 2. Then there is no possible prefix letter for P : 0 gives 010, 1 gives 11, and 2 gives 212. -Case M = 2N . Then w contains M M = 2N 2N . If N = P 1, then w contains 2P 12P 1. This factor cannot extend to 2P 12P 12, since this is 2Y Y with Y = P 12. So w contains 2P 12P 10.

  The four bi-infinite binary words avoiding AA.ABA.ABBA. The two bi-infinite binary words avoiding ABA.AABB.

Fig. 2

 2 Fig.2

  are: one pattern and two factors: • b 3 essentially avoids AA, 010, and 212. • A morphic image of b 5 essentially avoids AA, 010, and 020. • A morphic image of b 5 essentially avoids AA, 121, and 212. • b 2 essentially avoids ABABA, 000, and 111. two patterns: b 2 essentially avoids ABABA and AAA. one formula over three variables: b 4 and two words from b 4 obtained by letter permutation essentially avoid AB.AC.BA.BC.CA.

-

  0001101101110001100011 2 → 0001101101100011100111 Notice that {AA.ABA.BABB, AA.ABA.BBAB, SQ 5 } is 2-unavoidable. However, {AA.ABA.BABB, SQ 4 } is 2-avoidable: AA.ABB.BBAB is avoided with its reverse, 60-uniform morphism, avoids SQ 11 : Notice that {AA.ABB.BBAB, SQ 10 } is 2-unavoidable. -AA.ABAB.BB is self-reverse, 11-uniform morphism, avoids SQ 4 :

	0 → 00010010011000111001001100010011100100100111
	1 → 00010010011000100111001001100011100100100111
	2 → 00010010011000100111001001001100011100100111
	0 → 000110011100011001110011000111000110011100011100110001110011
	1 → 000110011100011001110001110011000111000110011100110001110011
	2 → 000110011100011001110001100111000111001100011100110001110011
	0 → 00100110111
	1 → 00100110001
	2 → 00100011011

-

  AA.ABBA.BAB is self-reverse, 30-uniform morphism, avoids SQ 6 :

	0 → 000110001110011000110011100111
	1 → 000110001100111001100011100111
	2 → 000110001100011001110011100111

-

  AAB.ABB.BBAA is self-reverse, 30-uniform morphism, avoids SQ 5 :

	0 → 000100101110100010110111011101
	1 → 000100101101110100010111011101
	2 → 000100010001011101110111010001

-

  AAB.ABBA.BAA is self-reverse, 38-uniform morphism, avoids SQ 5 : AABB.ABBA is unavoidable with its reverse, 193-uniform morphism, avoids SQ 16 :

	-0 → 00010001011011101100010110111000101101110111000101100010001011
	011101100010110111011100010110111011000101101110001011011101110001
	01100010001011011100010110111011100010110111011000101101110001011
	1 → 00010001011011101100010110111000101101110111000101100010001011
	011100010110111011100010110111011000101101110001011011101110001011
	00010001011011101100010110111011100010110111011000101101110001011
	2 → 00010001011011100010110111011100010110001000101101110110001011
	011101110001011011101100010110111000101101110111000101100010001011
	01110110001011011100010110111011100010110111011000101101110001011
	0 → 00010001000101110111010001011100011101
	1 → 00010001000101110100011100010111011101
	2 → 00010001000101110001110100010111011101


Robert E Kalaba

Harriet H Natsuyama

Firdaus E Udwadia

Y Y Fan

Reflections on the Gauss Principle of Least Constraint

Keywords: Gauss principle of least constraint, nonideal nonholonomic constraints, generalized inverses of matrices

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

Introduction

In his epochal paper of 1829, Gauss (Ref. 1) began by remarking that the D'Alembert principle reduced all of dynamics to statics and that the principle of virtual works reduced all of statics to a mathematical problem. Thus, there could be no new principle of mechanics that is not included already in those two. Yet, he observed that every new principle is not without merit, especially if it can shed new light on mechanical processes and perhaps render the solution of certain problems simpler to obtain. He then went on to state his own new principle, the principle of least constraint reducing all of mechanics, dynamics, and statics, to a single principle. Since its enunciation, it has been a cornerstone of analytical dynamics. His own derivation, which relies on the aforementioned two principles and the law of cosines, has often been presented but never surpassed for brevity and clarity.

This paper, besides being a tribute to Gauss in the 175th anniversary of the principle bearing his name, will present a rederivation of the Gauss principle and then an extension of his principle to cases in which the standard principle of virtual works is not applicable. These cases include those in which sliding friction (for example) is significant and cannot be neglected.

Gauss Principle of Least Constraint

Consider a system of n particles. The ith particle has mass m i and its position vector is x i in an inertial Cartesian frame of reference. Its velocity is ẋi and its acceleration is ẍi . It is subjected to an impressed force f i ; if no constraints are present, its free motion acceleration is a i = f i /m i .

The system mass matrix M is of dimension 3n by 3n and is a nonsingular diagonal matrix with the masses down the main diagonal in sets of three and zeroes elsewhere. The system displacement vector x i , of dimension 3n by 1, has the vectors x i stacked in the usual fashion. Similarly, we construct the system velocity and acceleration vectors ẋ and ẍ. Likewise for the system free acceleration vector a.

It is assumed that the system is subjected to m consistent equality constraint conditions of the form

f i (x, ẋ, t) = 0, i= 1, 2, . . . , m.
Note that this constraint form is general and includes both holonomic or nonholonomic constraints. Through use of the chain rule of differentiation, these constraints assume the form

A ẍ = b,
where A is a matrix of dimension m by 3n and b is a vector of dimension m by 1. Both may depend upon x, ẋ, t. Thus, the constraint on the acceleration is linear. Together with the initial conditions x(t 0 ) and ẋ(t 0 ), this constraint is equivalent to the original form

f i = 0, i = 1, 2, . . . , m.
For some years, it has been known (Ref. 2) that the actual system acceleration vector is given by the explicit formula

ẍ = a + M -1/2 (AM -1/2 ) + (b -Aa), (1)
where (AM -1/2 ) + denotes the usual pseudoinverse of the matrix AM -1/2 . At the present moment, we need not enter into the details of this explicit formula for the system acceleration, because later, in Eq. ( 17), we shall produce a more comprehensive formula from a more general point of view. Our aim now is to derive the Gauss principle of least constraint from formula [START_REF] Gauss | Ueber ein neues algemeines Grundgesetz der Mechanik[END_REF] It is convenient to make the substitutions

D = AM -1/2 , ( 2
)
e = b -Aa, ( 3
) ẍs = M 1/2 ẍ, (4)
a s = M 1/2 a. ( 5
.t. (AM -1/2 )( ẍs -a s ) = b -Aa. ( (9) s
) 10
Reverting now to the original variables, we see that the variational problem becomes

min ẍ ( ẍ -a) T M( ẍ -a), (11) s
.t. A ẍ = b. ( 12
)
Equations ( 11) and ( 12) constitute a form of the Gauss principle of least constraint.

General Equation of Motion

Next, we wish to view matters from a more general point of view. In fact, we wish to determine all the equations of motion that are compatible with the constraint A ẍ = b with no physical assumptions being made at all. Experience has shown the importance of the matrix AM -1/2 , so we shall rewrite the above equation in the form

AM -1/2 (M -1/2 ẍ) = b. ( 13
)
The theory of generalized inverses shows us that the general solution of this constraint set of linear algebraic equations for ẍ is

M 1/2 ẍ = (AM -1/2 ) + b + [I -(AM -1/2 ) + (AM -1/2 )]z, ( 14
)
where z is an arbitrary vector of dimension 3n by 1. Furthermore, it is most revealing to write the arbitrary vector z in the form

z = M 1/2 a + M -1/2 C, ( 15
)
where a is the free motion acceleration vector and C is an arbitrary vector, both being of dimension 3n by 1. It follows that

ẍ = a + M -1/2 (AM -1/2 ) + (b -Aa) +M -1/2 [I -(AM -1/2 ) + (AM -1/2 )]M -1/2 C, ( 16
)
or

M ẍ = Ma + M 1/2 (AM -1/2 ) + (b -Aa) +M 1/2 [I -(AM -1/2 ) + (AM -1/2 )]M -1/2 C. ( 17
)
Thus, Equation (17) is the most general possible equation of motion that is compatible with the constraint relation A ẍ = b. No physical assumptions were employed in the derivation of this equation. We notice that Eq. ( 1) is a special case of Eq. ( 17). It is convenient to rewrite Eq. (17) in the form

M ẍ = F N + F L + F C , ( 18
)
where

F N = Ma, ( 19
)
F L = M 1/2 (AM -1/2 ) + (b -Aa), ( 20
)
F C = M 1/2 I -(AM -1/2 ) + AM -1/2 M -1/2 C. ( 21
)
The notation bears the names of Newton, Lagrange, and Coulomb. Later, we shall elaborate on this. Note that only two essential mathematical ideas have entered the analysis: the chain rule of differentiation and generalized inverses of matrices. Modern computing environments, such as Matlab, have built-in commands for calculating the generalized inverse of a matrix, so it makes the approach highly suitable for numerical studies.

Extended Gauss Principle

Once again, we emphasize that Eq. ( 17) is the most general possible equation of motion that is compatible with the constraint condition A ẍ = b, of course assuming that the matrix M is nonsingular. Now, we wish to obtain an extended Gauss principle that leads to this equation of motion. Using the same notation as earlier, plus the additional notation

c s = M -1/2 C,
we find that Eq. ( 17) becomes

ẍ = a s + D + e + [I -D + D]c s . ( 22
)
We rearrange this to yield

ẍs -a s -c s = D + (e -Dc s ). ( 23
)
Recalling Eqs. ( 7) and (8), we find that ẍs solves the problem

min ẍs [ ẍs -(a s + c s )] T [ ẍs -(a s + c s )], (24) s
.t AM -1/2 [ ẍs -(a s + c s )] = b -Aa -AM -1 C. ( 25
)
In terms of the original variables, this problem is

min ẍ [ ẍ -a -M -1 C] T M[ ẍ -a -M -1 C], (26)
s.t. A ẍ = b. ( 27
)
Equations ( 26) and ( 27) constitute an extended Gauss principle of least constraint. This principle covers the cases of nonideal and of course nonholonomic constraints.

Physical Significance of the Terms in Eq. (17)

Let us now comment further on the terms in the general equation of motion, Eq. (17) or Eq. (18). The special vector a is chosen to be the acceleration that the system would have if there were no constraints. This is why in Eq. ( 19) we have denoted Ma as being F N , the Newtonian impressed force vector. Notice that Eq. (17) involves only the concepts of mass, distance, and time. It is in the physical interpretation of F N as Ma that we introduce a further physical concept, that of force. Now, we consider Eq. (20). We introduce the notion of a virtual displacement vector ν to be any vector ν that satisfies the homogeneous equation Av = 0. Thus, if ẍ is the actual acceleration at time t of the system, then ẍ + v is a possible acceleration, since

A( ẍ + v) = A ẍ + Av = A ẍ = b. ( 28
)
The equation Av = 0 is homogeneous, so that the vector v can be multiplied by any dimensional constant and it remains a solution. Thus, v is not necessarily a displacement. Finally, we introduce the notion of the work done by the force F L in a virtual displacement ν to be ν T F L . Since ν satisfies the equation

(AM -1/2 )M 1/2 v = 0, (29)
we have

M 1/2 v = [I -(AM -1/2 ) + AM -1/2 ]q, ( 30
)
where q is an arbitrary vector of dimension 3n. Thus,

v T = q T [I -(AM -1/2 ) + AM -1/2 ]M -1/2 . ( 31
)
It follows that

v T F L = q T [I -(AM -1/2 ) + AM -1/2 ]M -1/2 M 1/2 (AM -1/2 ) + (b -Aa) = 0. ( 32
)
Thus, F L is a constraint force that does no work on the system in a virtual displacement ν. We see that the other constraint force F C would be the null vector if C = 0, so F L must be the constraint force that maintains the constraints while doing no work on the system in any virtual displacement. It is the force that Lagrange denoted A T λ, where λ is an m-dimensional vector of Lagrange multipliers. It should be remembered that the columns of A + and A T span the same space.

There remains the third term on the right in Eq. ( 17),

F C = M 1/2 [I -(AM -1/2 ) + AM -1/2 ]M -1/2 C. ( 33
)
How much work does it do in a virtual displacement? We have

v T F C = q T [I -(AM -1/2 ) + AM -1/2 ]M -1/2 ×M 1/2 [I -(AM -1/2 ) + (AM -1/2 )]M -1/2 C = q T [I -(AM -1/2 ) + AM -1/2 ]M -1/2 C = v T C. ( 34
)
From the last equation, we see that a specification of the vector C is a specification of the work done by the constraint force F C in a virtual displacement v. Since the total force of constraint is F L + F C and since the force F L does no work in a virtual displacement, we have

v T (F L + F C ) = v T C. ( 35
)
This equation is a generalization of the classical equation for the principle of virtual work. That principle, a cornerstone of analytical mechanics since it was enunciated by Lagrange, states that the constraint forces do no work in any virtual displacement. This is referred to as an ideal constraint.

A constraint force such as sliding friction is ruled out of consideration, as Pars and Goldstein have stated in their textbooks. The new principle in Eq. ( 35) allows the constraint forces to do work in virtual displacements.

Toward a Theory of Underdetermined Systems

As general as the previous discussion has been, it is possible to go even beyond the bounds of mechanics. In fact, we can even begin to see the outlines of a theory of underdetermined systems. Now, let x = x(t) be the state vector of a system S and let the dimension of x be n. Next, suppose that the theory suggests that certain relations must be satisfied. We may also desire certain relations to be fulfilled. Thus, at least in the formative states of a theory, we are led to m relations of the form

f i (x, ẋ, t) = 0, i= 1, 2, . . . , m. ( 36
)
Notice that the number of relations (constraints) may be less than, equal to, or greater than the dimension of x. Through use of the chain rule of differentiation, we arrive at the equation

A ẍ = b, ( 37
)
where the matrix A is of dimension m by n and the vector b is of dimension m by 1.

Again, the theory of generalized inverses of matrices lets us write the general solution of this consistent system of linear algebraic equations in the form

ẍ = A + b + (I -A + A)z, ( 38
)
where A + is the generalized inverse of A and z is an arbitrary n-dimensional vector. We may write the arbitrary vector z in the form z = s + c, where s is a special vector that is up to the modeler choice and c is an arbitrary vector. Then, Eq. (38) becomes

ẍ = A + b + (I -A + A)(s + c) (39) or ẍ = s + A + (b -As) + (I -A + A)c. ( 40
)
This equation may be written as

ẍ = s + v 1 + v 2 , ( 41
)
where

v 1 = A + (b -As), ( 42
)
v 2 = (I -A + A)c. ( 43
)
The second-order ordinary differential equation above becomes the equation of motion of system S, based on the partial knowledge that we possess, as expressed in the m equations of constraint. The terms s, v 1 , v 2 may now be investigated somewhat more. If there were no constraints, then A = b = 0, so that ẍ = z. This indicates that the arbitrary vector z in Equation ( 38) is what ẍ would be if there were no constraints present, as is intuitively obvious. When there are constraints, we could follow the hints provided by mechanics. We let v be a vector in the null space of A, for which Av = 0. The vector v must have the form

v = (I -A + A)w,
where w is an arbitrary vector. Since I -A + A is symmetric and

A + AA + = A + , it follows that v T v 1 = 0. ( 44
)
Lastly, with regard to the vector ν 2 , we see that

v T ν 2 = w T (I -A + A)(I -A + A)c = w T (I -A + A)c, (45)
since the matrix I -A + A is also idempotent. But this implies that

ν T ν 2 = ν T c. ( 46
)
In addtion, from their definitions, we see that the vectors ν 1 and ν 2 are orthogonal,

ν T 2 ν 1 = 0. (47)
In summary, we may say that the underdetermined system S satisfies the second-order ordinary differential equation (33). The arbitrary vector z is what ẍ would be if there were no constraints. The matrix A and the vector b come from the constraints. This shows the nature of the vector ν l . Finally, a specification of the vector c is a prescription of

ν T (ν l + ν 2 ) = ν T ν 2 (48) or ν T ( ẍ -s) = ν T c, (49)
where

z = s + c.
If c is in the null space of A, so that Ac = 0, then ν 2 = c. If c is in the range space of A T , then ν 2 = 0.

Discussion

The entire analysis presented in this paper was started directly with the constraints on the system, which emphasizes again the importance of the constraints in a mechanical system. In fact, it is the constraints that make a set of point masses and rigid bodies into a system. Our previous analysis shows that the Lagrange principle of virtual works is a sagacious hypothesis that F c = 0. This assumption is what makes analytical dynamics into pure mathematics and makes it possible to eliminate all thermodynamic considerations. This assumption works well because most practical mechanical systems, through design and the use of lubricants, do minimize the effects of constraint forces that do work on a system. However, to be more precise, nonideal constraint forces must also be considered.

This paper considers mechanical systems including constraints that are holonomic or nonholonomic and also ideal or nonideal. It also looks into more general underdetermined systems. This widens the applicability of the magnificent contribution of the Gauss principle to human thought.


Charlotte Danino

Redundancy and Analogy: a Cognitive Discourse Analysis look at as in spoken data

Keywords:

Context: online description of a live event. Necessity to describe and explain events as perceived and understood on the go.

Research Question: why did speakers clearly prefer AS, and not while, when, because, since? As= polysemous and rather vague. Importance of cotext Objective: criticism of approach in terms of identification.

• costly cognitive operation (plausibility?)

• not consider the value of as but the value of the links between informational content 1 and informational content 2; • information cannot be identical to be compared or assimilated

Limitations

For as much as analyses are clearly essential to our research, this work remains primarily concerned by theoretical problems. This implies that the corpus is considered as a way to explore the epistemological and methodological problems singled out by our reflections. In other words, our approach is both corpus-based and corpus-driven, both inductive and deductive.

Before coming back to the particular status of the corpus in this study, it is necessary to present a detailed model of it and of the idea of constraint.

Situation, discourse and linguistic productions

The "funnel" model presents the type of constraints weighing on/in the speakers' verbalizations, which the dissertation targets. The principle is that a superordinate level

  (Langacker 1987, Talmy 2000): link between linguistic and conceptual structures -Cognitive Discourse Analysis (Tenbrink 2015): language activity reflects current state of thought. Speakers verbalize from a set of options they choose from Key concepts Meaning= instructional dynamic construction (Col et al 2012) Analogy= cognitive process (Hofstadter/Sanders 2013) Redundancy= same information coded several times Research framework -First 4 hours of CNN's live broadcast on 9/11 (70 speakers) -Highly constrained speech situation = ''natural'' control parameters Corpus Results Conclusion Redundancy and Analogy: a Cognitive Discourse Analysis look at as in spoken data Charlotte Danino (Université de Versailles Saint-Quentin, France) charlotte.danino@hotmail.fr -Corpus study: annotation of 146 conjunctive uses -Analec software: Semantic annotation + geometric representation of grammatical features and semantic repertories-Case study: must be tested on other data -Based on conjunctive uses: must be applied in greater detail to comparative and prepositional uses (as a journalist, as tall as he is, etc.)

  Figure 10 -The Funnel Heuristic: modeling the corpus


L Cortés

J.-M Mencik

S Méo

Interpolatory model reduction for component mode synthesis analysis of structures involving substructures with frequency-dependent parameters

An interpolatory model order reduction (MOR) strategy is proposed to compute the harmonic forced response of structures built up of substructures with frequency-dependent parameters. In this framework, the Craig-Bampton (CB) method is used for modeling each substructure by means of static modes and a reduced number of fixed-interface modes which are interpolated between several master frequencies. Emphasis is on the analysis of several substructures which can vibrate at different scales and, as such, do not need to be modeled with the same sets of interpolation points, depending on whether their modal density is low or high. For this purpose, an error indicator is developed to determine, through greedy algorithm procedure, the optimal number of interpolation points needed for each substructure. Additional investigations concern the selection of the fixed-interface modes which need to be retained for each substructure. Numerical experiments are carried out to highlight the relevance of the proposed approach, in terms of computational saving and accuracy.

Introduction

The problem of predicting the dynamic response of structures which are composed of substructures with frequency-dependent parameters is addressed in the present work. Those substructures may refer, for instance, to structural components and viscoelastic junctions connected to each other, where it is understood that the material characteristics of the junctions vary over the frequency domain. Also, the substructures are likely to vibrate at different scales, i.e., their modal densities can be largely disparate. Assessing the frequency response of this kind of multi-physics and multi-scale vibrating systems is the motivation behind the present work. Emphasis is placed on the development of efficient finite element (FE) based MOR approaches which are accurate and low-time consuming.

Recent works have demonstrated the potentiality of interpolatory MOR methods for analyzing dynamic models with parameter-dependent matrices [START_REF] Panzer | Parametric Model Order Reduction by Matrix Interpolation[END_REF][START_REF] Amsallem | An Online Method for Interpolating Linear Parametric Reduced-Order Models[END_REF]. These techniques consist in computing a finite set of reduced-order bases at some interpolation points -i.e., at some discrete values of the varying parametersand interpolating reduced matrices -i.e., once projected in the subspaces spanned by the reduced basesbetween the interpolation points. One of the key idea in these methods is the consideration of a Modal Assurance Criterion (MAC) -based congruence transformation to ensure the compatibility among the coordinates systems of the reduced bases, i.e., between the interpolation points [START_REF] Panzer | Parametric Model Order Reduction by Matrix Interpolation[END_REF][START_REF] Amsallem | An Online Method for Interpolating Linear Parametric Reduced-Order Models[END_REF]. For frequency-dependent matrices, another approach has been proposed in [START_REF] Ohayon | Clarification about Component Mode Synthesis Methods for Substructures with Physical Flexible Interfaces[END_REF] which interpolates reduced bases of vibration modes, rather than projected matrices, between master frequencies. In this framework, the Craig-Bampton (CB) method [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF] is used for modeling substructures by means of static modes and reduced sets of fixed-interface modes which are thus interpolated, e.g., using Lagrange polynomials. The static modes are being kept constant and are computed from the stiffness matrices of the substructures expressed at a certain frequency f , e.g., f = 0. The interesting feature of this approach is that it involves interpolating a matrix of fixed-interface modes only, instead of interpolating several matrices and vectors in different ways which may penalize the determination of the subspaces between the interpolation points. This technique will be considered and improved in the present paper for modeling structures involving several substructures with frequency-dependent parameters. The originality of the present work lies in the consideration of several substructures which can vibrate at different scales, i.e., they do not need to be modeled with the same sets of interpolation points. Here, an error indicator is proposed to determine, through greedy algorithm procedure, the optimal number of interpolation points -i.e., the number of master frequencies -for each substructure. The development of this error indicator follows from the analysis of the force balance equation on each interface between two connected substructures. It can be shown that this procedure is much more relevant than the usual one when a residual force vector is analyzed at a few excitation points, regardless of the coupling conditions between the substructures. Additional investigations are also proposed to determine the number of fixed-interface modes for modeling each substructure.

The rest of the paper is organized as follows. In Section 2, the basics of the CB method, for modeling a set of coupled substructures by means of static modes and fixed-interface modes, are recalled. In Section 3, the interpolatory MOR strategy is presented. The matrices of fixed-interface modes of the substructures are assessed and post-treated according to the MAC criterion proposed in [START_REF] Panzer | Parametric Model Order Reduction by Matrix Interpolation[END_REF][START_REF] Amsallem | An Online Method for Interpolating Linear Parametric Reduced-Order Models[END_REF]. Also, the error indicator used to select the number of interpolation points for each substructure is derived. In Section 4, numerical experiments are carried out concerning an assembly of plates with two different frequency-dependent substructures and a frequency-independent substructure. The efficiency of the proposed approach is highlighted in comparison with the conventional CB method. Additional discussions are brought about the selection of the number of fixed-interfaces modes per substructure.

Craig Bampton method

The Craig Bampton (CB) method is one of the most popular Component Mode Synthesis (CMS) techniques which combine substructuring and model order reduction to solve large-scale dynamic systems in an efficient way. The method has been initiated in [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF] and is recalled hereafter of the sake of clarity.

Consider the FE model of a structure which is composed of several substructures, see Figure 1. Each substructure is assumed to be linear, elastic and dissipative. Also, for the sake of simplicity, assumption is made that the FE meshes of the substructures are compatible across their coupling interfaces.

Consider a given substructure s. In the frequency domain, its dynamic equilibrium equation is given by

-ω 2 M s + iωC s + K s q s = F s , (1)
where M s , C s and K s are the mass, damping and stiffness matrices of the substructure, respectively, which are supposed to be real symmetric; q s and F s are the vectors of nodal displacements and nodal forces, respectively; also, ω is the angular frequency. Eq. ( 1) can be partitioned into boundary (B) and internal (I) degrees of freedom (DOFs) -i.e., those which belong and do not belong to the coupling interface, respectively -, as follows:

-

ω 2 M s BB M s BI M s IB M s II + iω C s BB C s BI C s IB C s II + K s BB K s BI K s IB K s II q s B q s I = F s B F s I . (2)
Denote as n s B and n s I the number of boundary and internal DOFs, respectively. Within the CB framework, the vector of nodal displacements q s I is expressed in terms of static modes {χ s j }

q s I ≈ X s qs B + m s I j=1 φs j αs j where qs B ≈ q s B , (3)
where qs B is the approximate vector of boundary DOFs -i.e., which results from the resolution of the CB model -, while αs j refers to generalized coordinates; also, X s is the matrix of static modes, which is defined by:

X s = [χ s 1 • • • χ s n s B ] = -(K s II ) -1 K s IB . (4)
On the other hand, the fixed-interface modes are the eigenvectors of the matrix pencil (K s II , M s II ). In matrix form, Eq. ( 3) is written by

q s I ≈ X s qs B + Φs αs , (5)
where Φs = [ φs

1 • • • φs m s I ] and αs = [α s 1 • • • αs m s I ]
T are the n s I × m s I matrix of fixed-interface modes and the m s I × 1 vector of generalized coordinates, respectively. As a result, the vector of nodal displacements q s is expressed by

q s ≈ qs B qs I = Ts qs B αs , (6)
where Ts = I 0 X s Φs .

(7)

As it turns out, the dynamic equilibrium equation (1) can be approximated through Galerkin procedure, as follows:

-ω 2 Ms + iω Cs + Ks qs

B αs = Fs , (8)
where Ms , Cs and Ks are reduced matrices of size (n s B + m s I ) × (n s B + m s I ). The matrices Ms and Cs are written as

Ms = ( Ts ) T M s Ts = M s BB + M s BI X s + (X s ) T M s IB + (X s ) T M s II X s M s BI Φs + (X s ) T M s II Φs ( Φs ) T M s IB + ( Φs ) T M s II X s ( Φs ) T M s II Φs , (9)
and

Cs = ( Ts ) T C s Ts = C s BB + C s BI X s + (X s ) T C s IB + (X s ) T C s II X s C s BI Φs + (X s ) T C s II Φs ( Φs ) T C s IB + ( Φs ) T C s II X s ( Φs ) T C s II Φs .
(10) Also, from Eq. ( 4), the matrix Ks is written as

Ks = ( Ts ) T K s Ts = K s BB + K s BI X s 0 0 ( Φs ) T K s II Φs . (11)
Finally, in Eq. ( 8), Fs denotes the (n s B + m s I ) × 1 vector of excitations:

Fs = ( Ts ) T F s = F s B + (X s ) T F s I ( Φs ) T F s I . ( 12
)
The CB modeling of an assembly of N substructures, connected to each other through several coupling interfaces, can be readily addressed through classic FE assembly procedure. Define the global vector of generalized coordinates α:

α =    α1 . . . αN    . (13)
Also, denote as qB ≈ q B the approximate global vector of boundary DOFs of the whole structure. The local and global vectors are linked as follows:

qs B αs = L s qB α s = 1, . . . , N, (14)
where L s denotes a Boolean localization matrix for the substructure s. Hence, the Galerkin approximation of the dynamic equilibrium equation of the whole structure is expressed by

-ω 2 M + iω C + K qB α = F, (15)
where

M = N s=1 (L s ) T Ms L s , C = N s=1 (L s ) T Cs L s , K = N s=1 (L s ) T Ks L s . ( 16
)
Also, F is the vector of external forces which are applied to the structure:

F = N s=1 (L s ) T Fs . ( 17
)
3 Interpolatory MOR

Introduction

In dynamic analysis, the mass, damping and stiffness matrices are often assumed to be frequency-independent and, consequently, a single CB reduced modeling is enough to analyze the dynamic response of a structure over a certain frequency range. Things are not so straightforward when materials having frequencydependent parameters are dealt with. In this case, the mass, damping and stiffness matrices are no longer constant, i.e., M s = M s (ω), C s = C s (ω) and/or K s = K s (ω) for some substructure s composing a whole structure. Also, the matrices of fixed-interface modes and static modes depend on the frequency -i.e., Φ s = Φ s (ω) and X s = X s (ω) -, meaning that they have to be recomputed at each discrete frequency considered in the frequency band of interest. This leads to computational costs which may largely exceed those involved when solving the standard unreduced FE problem -e.g., using a sparse solver -, which makes the CB method inefficient. To solve this issue, an interpolatory MOR strategy is proposed which aims at computing the matrices of fixed-interface modes and static modes of the substructures at a small number of master frequencies, only. The key idea is to linearly interpolate those matrices between the master frequencies so as to obtain approximate expressions -i.e., at any intermediary point between two consecutive master frequencies -, which can be done in a very fast way.

However, due to the frequency-dependency of the matrices, reduced bases of fixed-interface modes computed at different master frequencies span different subspaces with "coordinates systems" which are not necessarily compatible to each other. The issue may be viewed as to track each mode along the "master" frequency range, i.e., among all the reduced bases. As a result, directly interpolating the reduced bases of fixedinterface modes, without any mode tracking procedure, can be awkward. To solve this issue, a MAC-like procedure has been proposed in [START_REF] Panzer | Parametric Model Order Reduction by Matrix Interpolation[END_REF][START_REF] Amsallem | An Online Method for Interpolating Linear Parametric Reduced-Order Models[END_REF]. In this framework, the reduced bases are first transformed into a set of generalized coordinates with respect to a reference subspace of the same order spanned by a reference basis.

The proposed procedure enforces maximum correlation between the transformed bases and the reference basis. The choice of the reference basis is also addressed in [START_REF] Panzer | Parametric Model Order Reduction by Matrix Interpolation[END_REF] and it greatly determines the accuracy of the interpolation approximation. One possibility is to arbitrarily choose one of the master bases -i.e., those computed at the master frequencies -as being the reference subspace and transform the remaining master bases accordingly. However, there is no certainty that the chosen basis is the one that approximates the most important dynamics of the system best. A more convenient approach consists in computing the singular value decomposition (SVD) of the set of master bases in order to compute a meaningful reference basis, as proposed in [START_REF] Panzer | Parametric Model Order Reduction by Matrix Interpolation[END_REF]. Such an approach will be applied here to interpolate frequency-dependent reduced mode bases (see next subsection).

Mode basis interpolation

Consider a given substructure s and a set of M master angular frequencies {Ω s p } M p=1 for which M matrices of fixed-interface modes { Φs p } M p=1 -i.e., Φs p = Φs (Ω p ) for p ∈ {1, . . . , M } -are computed. Also, assume that the matrices Φs p have the same size n s I × m s I . Consider a so-called reference matrix of modes Ψs . As proposed in [START_REF] Panzer | Parametric Model Order Reduction by Matrix Interpolation[END_REF], the choice of a suitable reference basis is made by considering a singular value decomposition (SVD), as follows:

Φs 1 • • • Φs M = U s Σ s (V s ) T . ( 18
)
The reference matrix of modes Ψs is then determined by considering the first m s I left singular vectors, i.e., those associated with the m s I highest singular values:

Ψs = U s ( : , 1 : m s I ). (19)
As explained in [START_REF] Panzer | Parametric Model Order Reduction by Matrix Interpolation[END_REF], each matrix of modes needs to be tightly connected with the reference matrix Ψs , which is done by considering the following modal coordinates transformation:

Φs p = Φs p Rs p p = 1, . . . , M, (20)
Rs p = Ûs p ( Vs p ) T p = 1, . . . , M. (22)
In the present work, a linear interpolation scheme of the matrix of modes is dealt with. Hence, for a substructure s having frequency-dependent parameters, a matrix Ts p (see Eq. ( 7)) can be defined at each master frequency Ω p as

Ts p = I 0 X s p Φs p , (23)
where X s p = X s (Ω p ) is the matrix of static modes at Ω p :

X s p = -K s II (Ω p ) -1 K s IB (Ω p ). ( 24
)
If one assumes that the master frequencies {Ω s p } M p=1 are sorted in ascending order, the interpolated transformed matrix of modes, at a given frequency ω between two consecutive master frequencies Ω s p and Ω s p+1 , is given by

Ts (ω) = Ω s p+1 -ω Ω s p+1 -Ω s p Ts p + ω -Ω s p Ω s p+1 -Ω s p Ts p+1 for Ω s p ≤ ω < Ω s p+1 p = 1, . . . , M -1. (25)
As it turns out, the following alternative mode expansion can be proposed to assess the vector of nodal displacements of the substructure:

q s (ω) ≈ qs B qs I = Ts (ω) qs B αs for Ω s p ≤ ω < Ω s p+1 p = 1, . . . , M -1, (26)
where qs B ≈ q s B , while αs = αs (ω) is the m s I × 1 vector of generalized coordinates associated to the interpolated reduced matrix of modes Φs (ω). As a result, a modified Galerkin approximation of the dynamic equilibrium equation of the whole structure is to be considered, as follows:

-ω 2 M(ω) + iω Ĉ(ω) + K(ω) qB α = F, (27)
where matrices M(ω), K(ω) and Ĉ(ω) are given by:

Ms (ω) = Ts (ω) T M s (ω) Ts (ω), (28) Ĉs (ω) = Ts (ω) T C s (ω) Ts (ω),
(29) Ks (ω) = Ts (ω) T K s (ω) Ts (ω).

(30)

Selection of master frequencies

The efficiency of the interpolatory MOR strategy is mainly determined by the number of master frequencies considered, since the computation of the matrices of modes at many master frequencies is likely to be cumbersome. On the other hand, the accuracy of the approximate solutions is proportional to the number of master frequencies. In other words, the smaller the number of master frequencies, the larger the frequency gap between two master bases and, consequently, the greater the approximation error in interpolating those master bases. Thus, there is a tradeoff between computational cost and accuracy with respect to the number of master frequencies considered.

The analysis of a structure through "substructured" interpolatory MOR provides further benefits. In this framework, the interpolatory reduction is performed at the substructure level, i.e., the number of master frequencies may vary from one substructure to another. This offers the possibility to adapt the computational burden among the substructures. For instance, a stiff (frequency-dependent) substructure may be reduced using a very small number of master frequencies, whereas a more flexible (frequency-dependent) substructure would require a larger number of master frequencies for obtaining the same interpolation accuracy in a given frequency range. Therefore, frequency-dependent components with different numbers of master frequencies may coexist in a same CB model, together with frequency-independent components.

Such a strategy is proposed here. Notice that the number of master frequencies is defined at the substructure level and it determines the accuracy of the approximate solution. Then, the selection of the number of master frequencies per substructure appears to be a crucial task, which requires the use of an efficient error indicator (see next subsections).

A greedy-based selection algorithm which determines the required number of master frequencies per substructure is proposed. The procedure starts by defining, for each substructure s, a coarse set of equally-spaced master frequencies {Ω s p } M s p=1 . The following steps are then considered:

1. Compute the matrix of transformed modes Φs p at each master frequency Ω s p , see Eq. (20).

2. For each intermediary frequency ω considered within the frequency band [Ω s p , Ω s p+1 ), compute the approximate solution, see Eq. ( 27).

3. Set as a new master frequency any intermediary frequency ω for which the error indicator ¯ s (ω) exceeds the expected tolerance threshold, see Eq. (36).

The greedy algorithm runs until the error indicator ¯ s (ω), for each substructure s, is under the expected tolerance threshold at each intermediary frequency ω. Notice that the matrices of fixed-interface modes { Φs p } p and static modes {X s p } p computed (and stored) at previous iterations do not need to be recomputed at subsequent iterations.

The success of the greedy algorithm greatly relies on the interpolatory error estimate. Computing the true error of the approximate solution is unfeasible since the exact solution would be required for that purpose. Instead, a computationally cheap error indicator must be considered (see next subsection). The error indicator must be tight enough so that the greedy algorithm does not overestimate the number of necessary master frequencies.

Error indicator

Assessing the efficiency of a MOR approach in structural dynamics usually involves ensuring that the vectors of input forces are well described by the approximate displacement solutions [START_REF] Hetmaniuk | An adaptive scheme for a class of interpolatory model reduction methods for frequency response problems[END_REF]. The issue hence consists in analyzing a set of residual force vectors at a few excitation points. One major drawback of this approach is that it is rather local and discards the analysis of the coupling conditions between the substructures, which is a crucial point when coupled systems are dealt with.

An alternative error analysis is proposed here to solve this issue. In this framework, the force balance equation on each interface between two (or more) connected substructures is invoked, which appears to be much more precise than the usual local error analysis based on residual forces. The proposed strategy is derived as follows.

Eq. ( 1) can be rewritten as

D s q s B q s I = F s B F s I , (31)
where D s = -ω 2 M s (ω) + iωC s (ω) + K s (ω) is the dynamic stiffness matrix of a substructure s which, in matrix form, is expressed by

D s = D s BB D s BI D s IB D s II . (32)
In Eq. (31), F s B is the vector of boundary forces, which is given by F s B = D s BB q s B + D s BI q s I . On the other hand, by considering the approximate vectors of nodal displacements qs B and qs I (see Eq. ( 26)), the vector of boundary forces is rather expressed as

D s BB qs B + D s BI qs I .
As it turns out, the error made at the level of the substructure s is expressed by:

Rs B = F s B -(D s BB qs B + D s BI qs I ) . (33)
Also, the error made at the level of a whole structure, made up of N substructures, is expressed by

RB = N s=1 (L s B ) T Rs B = F B - N s=1 (L s B ) T (D s BB qs B + D s BI qs I ) , (34)
where L s B is a Boolean localization matrix concerning the boundary DOFs of a substructure s, and

F B = N s=1 (L s B ) T F s B
is the vector of boundary forces -i.e., the vector of input forces which are applied to the coupling interfaces between the substructures -, which is supposed to be known. Notice that L s B (L s B ) T = 0 for two substructures s and s coupled through a common interface. In this sense, a relative error made when expressing the force balance equation on the coupling interface between a substructure s and one (or more) connected substructure can be derived, as follows:

s (ω) = L s B RB 2 F 2 . ( 35
)
where • 2 denotes the 2-norm, while F is the vector of input forces applied to the structure, which is known.

Scaled error indicator

The error indicator s is expressed in terms of boundary forces, and, as such, does not reflect the true displacement error. To solve this issue, a scaled error indicator may be proposed as follows:

¯ s (ω) = ˜ HR q R s (ω), (36)
where ˜ HR q is an estimate of the displacement error issued from the CB modeling:

˜ HR q = q -qHR 2 qHR 2 , ( 37
)
where qHR is the displacement vector issued from a "high resolution" CB modeling, i.e., with a large number of fixed-interface modes for each substructure (e.g., twice as many as those usually required). Also, R is the error indicator for a given reference substructure R (R ∈ {1, . . . , N }) which may be arbitrarily defined.

For practical purpose, both error terms ˜ HR q and R can be set as constant, i.e., they need to be computed at one single master frequency only, e.g., the maximum frequency within the frequency band of interest. By considering the scaled error indicator ¯ s , the difference between the boundary force error and the true displacement error is greatly reduced. In this way, a meaningful error tolerance threshold -e.g., 10% as this is usually done when estimating displacement errors -may be proposed, which can be directly compared to the scaled error indicator ¯ s at a low computational cost. Here, 90, 200 and 90 plate elements are respectively used to model substructures 1, 2 and 3, which yields 1287 DOFs for the whole structure. The structure has simply supported boundary conditions on one edge of each substructure, as shown in Figure 2. A 10N vertical harmonic point force is applied at location (x, y, z) = (0.71 m, 1 m, 0 m) in substructure 1, while the harmonic response of the structure is assessed at a measurement point at the interface between substructures 1 and 2 (Figure 2). The material properties of substructures 1, 2 and 3 are listed in Table 1. Here, the Young's moduli of substructures 1 and 3 are assumed to be frequency-dependent, i.e.:

E s (ω, x) = E s 0 + ωE s 1 (x), ( 38
)
where E s 1 (x) is the linear variation rate, which is here assumed to vary randomly in space (with mean value Ēs 1 ), and also, to be different between substructures 1 and 3. Also, the damping matrix of each substructure is assumed to be proportional to its stiffness matrix, as displayed in Table 1. The frequency response function (FRF) of the structure is assessed on a frequency band of [0 Hz, 50 Hz] by considering 251 equally-spaced sample frequencies. For each sample frequency, the harmonic response is computed through CB modeling. The selection of the number of fixed-interface modes for each substructure is addressed in the following subsection.

Mass

Selection of fixed-interface modes

The accuracy of the CB solution is linked to the number of fixed-interface modes which are retained for modeling each substructure. As a rule of thumb, the selected modes are those for which the related eigenfrequencies are below a certain cut-off frequency, say twice the maximum frequency in the frequency band of interest, i.e., 2f max = 2 × 50 = 100 Hz in the present case. Notice that for substructures having frequencydependent parameters, the eigenfrequencies are frequency-dependent as well. Hence, the aforementioned selection criterion is to be applied by considering the eigenfrequencies computed at a certain frequency only, e.g., f = f max = 50 Hz. This yields 3, 25 and 7 fixed-interface modes for substructures 1, 2 and 3, respectively, that is, those for which the eigenfrequencies computed at f max are below 2f max .

As a first test, the FRF issued from the CB method (without interpolation) is compared to that issued from the FE method as shown in Figure 3. Also, the relative errors induced when computing the displacement and force vectors are displayed. Those displacement and force errors are defined as follows:

˜ q (ω) = q -q 2 q 2 , ˜ F (ω) = Dq -F 2 F 2 , ( 39
)
where D is the dynamic stiffness matrix of the whole structure. Also, the errors ˜ q and ˜ F are compared to the error indicator s (see Eq. ( 35)) for s = 1, 3, i.e., for the frequency-dependent substructures.

˜ q (ω); - ˜ F (ω); -1 (ω); -3 (ω).
As it can be seen in Figure 3, the CB-based FRF appears to be close to the FE solution. Significant differences occur at high frequencies, which are due to the mode truncation process of the CB modeling. Regarding the errors, the true displacement error ˜ q is shown to be mostly greater than 10% over the frequency band. Notice that the force error ˜ F does not appear to be a good indicator of the true error ˜ q because it shows values exceeding 100% over the whole frequency band. On the contrary, the proposed error indicators 1 and 2 provide almost the same trends, in terms of curve shapes, than the displacement error.

As it turns out, the displacement error ˜ q is large especially at high frequencies. Hence, a higher cut-off frequency is required so as to include a larger number of fixed-interface modes for modeling each substructure. In Figure [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF], the displacement error at f max is plotted as a function of the cut-off frequency used to select the fixed-interface modes. As it can be seen, the cut-off frequency has to be much above 2f max so as to lower the errors at a satisfactory level. Thus, a cut-off frequency of 10f max = 500 Hz is chosen which provides a displacement error under 1%. This corresponds to 23, 137 and 45 modes for substructures 1, 2 and 3, respectively.

It may be further observed in Figure 4 that the gap between the error indicators 1 (ω) and 3 (ω), and the true error ˜ q , broadens as the cut-off frequency increases. This is even exacerbated for the force error ˜ F (ω). This makes it necessary to use a scaled error indicator, see Eq. (36). -˜ q (ω); -˜ F (ω); -1 (ω); -3 (ω).

Interpolatory MOR

The FRF of the structure is assessed using the interpolatory MOR strategy. Recall that a cut-off frequency of 500 Hz is considered to select the number of fixed-interface modes required for each substructure (see previous subsection). Substructure 2 is frequency-independent and, as such, does not requires interpolation of its modes. Regarding substructures 1 and 3, three equidistant master frequencies are initially considered, i.e., 0 Hz, 25 Hz and 50 Hz. Notice that the variation rate of the Young's modulus (Eq. ( 38)) of substructure 1 is ten times larger than that of substructure 3, so that different numbers of master frequencies may be required between these substructures. In fact, the greedy algorithm working with the scaled error indicator ¯ s and an error threshold of 10% (see Section 3.3) is used to determine the adequate number of master frequencies per substructure. (bottom) relative errors: -ˆ q (ω); -¯ 1 (ω); -¯ 3 (ω); --error threshold.

In Figure [START_REF] Hetmaniuk | An adaptive scheme for a class of interpolatory model reduction methods for frequency response problems[END_REF], the FRF is shown at the first and last iteration of the greedy algorithm, along with the true displacement error ˆ q : ˆ q (ω) = qq 2 q 2 . (40) Also, the scaled error indicators ¯ 1 (ω) and ¯ 3 (ω) are displayed.

At the first iteration where 3 master frequencies per substructure are considered, the computed FRF is clearly incorrect between the two first master frequencies. Both displacement error ˆ q and scaled error indicator ¯ 1 are greater than the threshold, while ¯ 3 lays under the threshold. Consequently, two new master frequencies have been added for modeling substructure 1 (only) at the second iteration. The greedy algorithm ends once ¯ 1 and ¯ 3 are under the threshold at the intermediary points between the master frequencies. This occurs at the fourth iteration where 9 and 3 master frequencies are respectively considered for substructures 1 and 3.

As expected, the true displacement error ˆ q appears to be smaller than the threshold at the fourth iteration as shown in Figure 5.

Table 2 shows the execution times required at each iteration of the greedy algorithm to select the master frequencies and compute the FRF. Hence, the total execution time is 49.2 s. In comparison, it would have taken 724.6 s to compute the FRF with the CB method, i.e., by computing the fixed-interface modes and static modes of substructures 1 and 3 at every sampling frequency. This means more than 93% time saving in benefit of the proposed interpolatory MOR strategy. Also, the relevance of the proposed scaled error indicator ¯ s is analyzed when compared to the one based on residual forces. For this purpose, the following scaled force error indicator ¯ F is considered:

¯ F (ω) = ˜ HR q R Dq -F 2 F 2 , (41)
The FRFs and errors at the first and last iteration of the greedy algorithm are shown in Figure 6. As it can be seen, the error indicator ¯ F differs much more from the displacement error ˆ q than the proposed error indicators ¯ 1 and ¯ 3 , see Figure 5. In this case, the greedy algorithm stops at the second iteration where 4 master frequencies are considered for both substructures 1 and 3. Notice that the error indicator ¯ F is global, i.e., the same number of master frequencies are added to each frequency-dependent substructure at each iteration of the greedy algorithm. In contrast, the proposed strategy enables the number of master frequencies to be optimized between the substructures. As a second drawback of the error indicator ¯ F , the displacement error cannot be lowered under the tolerance threshold at all sampling frequencies.

As it can be seen in Figure 6, the force error indicator ¯ F exhibits a rather flat shape over the frequency band. This appears to be a consequence of the CB modeling in the sense that a modal reduction is performed at the DOF where the loading force is applied. The modal reduction ensures a good representation of the displacement field in substructures; however, the spatial distribution of the loading is not necessarily well captured by the same modal basis. Hence, significant errors may be introduced when approximating the loading term in the dynamic equilibrium equation. This may be observed in Figures 3 and4, where the non-scaled force error ˜ F is around or even above 100% over the whole frequency band.

To clarify this point, consider the reciprocal problem in which the loading point and the measurement point are interchanged, as shown in Figure 7. In this case the point force is located on the coupling interface master frequencies for substructure 1; × master frequencies for substructure 3. (bottom) relative errors: -ˆ q (ω); -¯ F (ω); --error threshold.

between substructures 1 and 2, so that no mode truncation approximation is introduced when expressing the loading term. Figure 8 shows that the force error ˜ F decreases as the number of fixed-interface modes increases in a similar manner as the displacement error ˜ q and the error indicators 1 and 3 . In this sense, the force error indicator ¯ F (Eq. ( 41)) is likely to be much more accurate when compared to Figure 6.

Again, the FRFs and errors are computed as shown in Figure 9. By using either the force error indicator ¯ F , or the proposed indicators ¯ 1 and ¯ 3 , the displacement error ˆ q is reduced under the threshold of 10% after 5 iterations. However, the numbers of master frequencies are significantly different, i.e., 10 and 4 master frequencies for substructures 1 and 3 (respectively) with the proposed approach, against 11 master frequencies for both substructures 1 and 3 with the residual force approach. This clearly demonstrates the fact that the proposed MOR strategy works in a more efficient way to select the master frequencies.

Conclusions

An interpolatory MOR strategy has been proposed to compute the forced response of structures composed of substructures with frequency-dependent parameters. The CB method has been considered for modeling the substructures, and a linear interpolation scheme has been proposed to approximate the matrices of substructure modes between a set of master frequencies. The main advantage of the proposed MOR strategy lies in the consideration of an error indicator which enables the number of master frequencies to be modulated and optimized among the substructures. The selection of those master frequencies is automatically achieved through greedy algorithm. Several test cases have been carried out which clearly demonstrate the efficiency of the proposed approach.

Figure 1 :

 1 Figure 1: FE model of an abstract structure composed of four substructures.

  The interpolatory MOR strategy is used to analyze the frequency forced response of an assembly of three plates -say, substructures 1, 2 and 3 -, as shown in Figure 2. Substructures 1 and 3 are made of steel material and have similar dimensions 1.6 m × 4 m × 0.01 m; on the other hand, substructure 2 consists of an aluminum plate with dimensions 2 m × 4 m × 0.01 m. The FE meshes of the substructures are composed of 4-node isotropic thin plate elements with 3 DOFs per node, and are compatible across the coupling interfaces.

Figure 2 :

 2 Figure 2: FE mesh of a three coupled plates structure with a 10N vertical point force in substructure 1 and a response measurement point at the interface between substructures 1 and 2.

Figure 3 :

 3 Figure 3: (top) FRF of the structure: -FE solution; • • • CB solution. (bottom) relative errors: -˜ q (ω); -˜ F (ω); -1 (ω); -3 (ω).

Figure 4 :

 4 Figure 4: Errors at f max as functions of the cut-off frequency.-˜ q (ω); -˜ F (ω); -1 (ω); -3 (ω).

Figure 5 :

 5 Figure 5: (top) FRF at the first (left) and last (right) iterations of the greedy algorithm: -FE solution; • • • interpolated solution; master frequencies for substructure 1; × master frequencies for substructure 3.(bottom) relative errors: -ˆ q (ω); -¯ 1 (ω); -¯ 3 (ω); --error threshold.

Figure 6 :

 6 Figure 6: (top) FRF at first (left) and last (right) iterations of the greedy algorithm: -FE solution; • • • interpolated solution;master frequencies for substructure 1; × master frequencies for substructure 3. (bottom) relative errors: -ˆ q (ω); -¯ F (ω); --error threshold.

Figure 7 :

 7 Figure 7: FE mesh of the three coupled plates structure with a 10N vertical point force acting on the interface between substructures 1 and 2, and a response measurement point in substructure 1.

Figure 8 :

 8 Figure 8: Errors at f max (Figure 7) as functions of the cutoff frequency. -˜ q (ω); -˜ F (ω); -1 (ω); -3 (ω).

Figure 9 :

 9 Figure 9: (top) FRF at the last iteration of the greedy algorithm, using the error indicator ˜ F (left) and the indicators 1 and 3 (right): -FE solution; • • • interpolated solution; master frequencies for substructure 1;× master frequencies for substructure 3. (bottom) relative errors: -ˆ q (ω); -¯ F (ω); -¯ 1 (ω); -¯ 3 (ω); --error threshold.

Table 1 :

 1 Material properties of substructures 1 (S1), 2 (S2) and 3 (S3).

			Elastic modulus	Poisson's	Damping
		density	E s 0	Ēs 1	ratio	proportionality
	S1 (steel)	7800 kg/m 3 210 GPa 2.1 GPa.s/rad	0.31	4 × 10 -5
	S2 (aluminum) 2700 kg/m 3 70 GPa	-	0.34	2 × 10 -5
	S3 (steel)	7800 kg/m 3 210 GPa 0.21 GPa.s/rad	0.31	4 × 10 -5

Table 2 :

 2 Execution times involved in the interpolatory MOR strategy.

Acknowledgements

The authors express their thanks for the financial support provided by "Le Commissariat à l'énergie atomique et aux énergies alternatives (CEA)" and "La Région Centre".