r5.01.02 Résolution du problème modal

r5.01.02 quadratique (QEP)

Résumé

L’étude de la stabilité dynamique d’une structure amortie et / ou tournante conduit à la résolution d’un problème modal plus élevé que les traditionnels problèmes modaux standards (SEP) ou généralisé (GEP).

Pour les appréhender, on propose des méthodes via l’opérateur [CALC_MODES] : puissances inverses et méthode de Müller, Lanczos, IRA et QZ. Elles ont chacune leur périmètre d’utilisation, leur avantages / inconvénients et leur historique de développement.

Dans la première partie du document nous résumons la problématique de résolution d’un problème quadratique et sa déclinaison dans l’architecture générale d’un calcul modal. Puis nous détaillons les aspects numériques, informatiques et fonctionnels de chacune des approches disponibles dans le code. Les différents résultats, algorithmes ou paramètres abordés dans ce document s’appuient souvent sur les méthodes modales d’ordres inférieures (SEP et GEP) décrites dans le document [Solveurs modaux et résolution du problème généralisé]. La lecture de ce dernier est donc un pré-requis conseillé.

Annexe_1__Interpr_tation_des_valeurs_propres_complexes

Contexte

Problématique

Nous considérons le problème modal quadratique (QEP)

\[\begin{split}\begin{array}{c}\text{Trouver }(\lambda, u) \text{ tel que :}\hfill \\ ({\lambda}^{2}\tensTwo{B}+\lambda \tensTwo{C}+\tensTwo{A})\vector{u} = \vectorZero \text{ avec } \vector{u} \ne \vectorZero \end{array}\end{split}\]

où \(\tensTwo{A}\), \(\tensTwo{B}\) et \(\tensTwo{C}\) ont des matrices carrées de taille \(n\), à coefficients réels ou complexes, symétriques ou non. Ce type de problème correspond, en mécanique, notamment à l’étude des vibrations libres d’une structure amortie et/ou tournante. Pour cette structure, on recherche les valeurs propres \({\lambda}_{i}\) (et leurs vecteurs propres associés \({\vector{u}}_{i}\) ) les plus proches, dans le plan complexe, d’une valeur de référence donnée (le «shift» \(\sigma\) ) pour savoir si une force excitatrice peut créer une résonance. Dans ce cas standard:

• La matrice \(\tensTwo{A}\) est la matrice de rigidité, notée \(\tensTwo{K}\), symétrique réelle (éventuellement augmentée de la matrice symétrique complexe, notée \(\tensTwo{E}_{\text{hyst}}\), si la structure présente un amortissement hystérétique) soit \(\tensTwo{A}=\tensTwo{K}+\tensTwo{E}_{\text{hyst}}\). Donc \(\tensTwo{A}\) est symétrique réelle ou complexe.

• La matrice \(\tensTwo{B}\) est la matrice de masse ou d’inertie notée \(\tensTwo{M}\) (symétrique réelle).

• La matrice \(\tensTwo{C}\) regroupe, elle, les éventuels effets gyroscopiques et ceux d’amortissement visqueux via la combinaison:

\[\tensTwo{C} := \tensTwo{E}_{\text{visq}} + \xi \tensTwo{G}\]

avec \(\tensTwo{E}_{\text{visq}}\), matrice (symétrique réelle) d’amortissement induit par des forces dissipatives, \(\tensTwo{G}\) matrice de gyroscopie (anti-symétrique réelle) et \(\xi\) un paramètre réel représentatif de la vitesse de rotation. Donc \(\tensTwo{C}\) est potentiellement non symétrique réelle.

Ce type de problématique est activé par le mot-clé TYPE_RESU='DYNAMIQUE'. Dans le cas général du QEP et en n’exploitant que leur propriété de symétrie et d’arithmétique des matrices, on peut décomposer le périmètre des opérateurs modaux comme suit.

Tableau 48 Périmètre d’utilisation de [CALC_MODES] selon les méthodes d’analyse, en fonction de propriétés des matrices \(\tensTwo{A}\) et \(\tensTwo{C}\) du QEP. \(\tensTwo{B}\) est réelle quelconque.

\(\tensTwo{A}\) et \(\tensTwo{C}\)	Réelle symétrique	Réelle non symétrique	Complexe
Réelle symétrique	Cas le plus courant SIMULT et INV sans restriction sur les méthodes; Modes réels ou complexes \(\left( \lambda, \tilde{\lambda} \right)\) .	SIMULT avec SORENSEN/QZ; Modes réels ou complexes \((\lambda, \tilde{\lambda}\) .	Cas non traité
Réelle non symétrique	SIMULT avec SORENSEN/QZ; Modes réels ou complexes \((\lambda, \tilde{\lambda}\)	SIMULT avec SORENSEN/QZ; Modes réels ou complexes \(\left( \lambda, \tilde{\lambda} \right)\) .	Cas non traité
Complexe symétrique	SIMULT avec SORENSEN/QZ; Modes réels, complexes quelconques ou \(\left( \lambda, \tilde{\lambda} \right)\)	Cas non traité	Cas non traité
Autres complexes (hermitien, non symétrique,…)	Cas non traité	Cas non traité	Cas non traité

En absence d’amortissement hystérétique, le QEP standard à résoudre est donc constitué de matrices réelles symétriques. Ces valeurs propres sont, soit réelles, soit complexes conjuguées par paire \(\left( \lambda, \tilde{\lambda} \right)\). Les vecteurs propres sont alors potentiellement à composantes complexes.

En présence d’amortissement hystérétique, le QEP devient complexe et perd la propriété précédente (Pour la recouvrer, il faut que les matrices soient, par exemple, complexes hermitiennes). Ces modes alors sont réelles, complexes par paire ou dépareillés.

Remarques

• Contrairement au GEP [Solveurs modaux et résolution du problème généralisé], le mot-clé TYPE_RESU='FLAMBEMENT' n’est pas licite dans le cadre des QEP.

• Toutes les options de calculs fonctionnent quel que soit le type de prise en compte des conditions limites: dualisation ou élimination (voir les commandes [AFFE_CHAR_MECA] et [AFFE_CHAR_CINE]). Toutefois, celles-ci ne doivent pas être redondantes: par exemple, utiliser la dualisation et l’élimination sur une même inconnue.

• Informatiquement, le périmètre de des options CENTRE, PLUS_PETITE et TOUT a été étendu aux matrices réelles non symétriques \(\tensTwo{A}\), \(\tensTwo{B}\) et \(\tensTwo{C}\) du QEP. Ce n’est pas encore le cas pour les options PROCHE et AJUSTE qui se limitent, pour l’instant, à des matrices réelles symétriques. Ces dernières ne peuvent donc pas prendre en compte d’effets gyroscopiques.

• La modélisation de l’amortissement [Modélisation de l’amortissement en dynamique linéaire] peut se décomposer en deux classes: amortissement visqueux proportionnel (dit de Rayleigh \(\tensTwo{E}_{\text{visq}}:=\alpha \tensTwo{K}+\beta \tensTwo{M}\)) ou amortissement hystérétique (\(\tensTwo{E}_{\text{hyst}}:=(1+i \, \eta )\tensTwo{K}\)). Chacun se décline au niveau global de la structure ou au moyen d’amortissements localisés (groupes de mailles, éléments discrets ad hoc). Ces derniers permettant de mieux représenter l’hétérogénéité de la structure par rapport à l’amortissement. Reste cependant la question délicate de l’identification des coefficients \((\alpha ,\beta, \eta )\) et de leur influence sur le résultat.

Propriétés des modes propres

Le tableau précédent ne tient pas compte des cas particuliers (gyroscopie pure, structure sur-amortie…). Pour mémoire, on va en énoncer quelques-uns que l’on retrouve souvent dans la littérature. Mais il faut bien garder en mémoire que le contexte des modélisations avec l’usage quasi-systématique de double Lagrange et les tris de modes propres programmés dans les opérateurs modaux ne tiennent pas compte de ces spécificités. Actuellement en QEP, on ne retient que les modes couplés \(\left( \lambda, \tilde{\lambda} \right)\) et on conserve celui à partie imaginaire positive. Suivant les cas de figure, la présence de valeurs propres d’autres natures (réelle ou complexe non conjuguée) est signalée par une alarme ou à titre informatif.

On a donc (en reprenant les éléments et les notations de [Solveurs modaux et résolution du problème généralisé] et du papier [TM01]), la distribution des propriétés suivantes. Celles-ci sont parfois cumulatives.

	Propriétés des matrices	Propriétés des valeurs propres	Propriétés des vecteurs propres	Exemple de liens avec les modélisations
1	\(\tensTwo{B}\) non singulier	\(2n\) valeurs propres \(\lambda\) finies		Structure non bloquée ou bloquée par élimination
2	\(\tensTwo{B}\) singulier	\(2n\) valeurs propres \(\lambda\) finies ou infinies		Structure bloquée par dualisation
3	\(\tensTwo{A}\), \(\tensTwo{B}\) et \(\tensTwo{C}\) réelles	Réelles ou par paire \(\left( \lambda, \tilde{\lambda} \right)\)	Si \(\vector{u}\) vecteur propre à droite de \(\lambda\) , alors \(\tilde{\vector{u}}\) est celui de \(\tilde{\lambda}\)	QEP classique» (pas d’amortissement hystérétique)
4	\(\tensTwo{A}\), \(\tensTwo{B}\) et \(\tensTwo{C}\) hermitiennes	Réelles ou par paire \(\left( \lambda, \tilde{\lambda} \right)\)	Si \(\vector{u}\) vecteur propre à droite de \(\lambda\), alors il est aussi vecteur propre à gauche de \(\tilde{\lambda}\)	QEP sans gyroscopie et avec amortissement visqueux
5	\(\tensTwo{A}\), \(\tensTwo{B}\) et \(\tensTwo{C}\) hermitiennes avec \(\tensTwo{B}>0\), \(\tensTwo{A} \ge 0\) et \(\tensTwo{C}\ge 0\)	\(\Re(\lambda)\le 0\)		Structure libre ou bloquée par élimination, avec un amortissement (dissipatif) de Rayleigh et sans gyroscopie
6	\(\tensTwo{A}\), \(\tensTwo{B}\) et \(\tensTwo{C}\) symétriques avec \(\tensTwo{B}>0\), \(\tensTwo{C} \ge 0\) et \(\tensTwo{A} \ge 0\) + Condition de sur-amortissement [2]_	Spectre réel négatif décomposé en deux paquets de taille \(n\)	\(n\) vecteurs propres linéairement indépendants sont associés aux \(n\) plus grandes (ou plus petites) valeurs propres	Structure libre ou bloquée par élimination, avec un amortissement ad hoc et sans gyroscopie
7	\(\tensTwo{A}\) et \(\tensTwo{B}\) hermitiennes et \(\tensTwo{C}=-{\tensTwo{C}}^{\star}\) (anti-hermitienne) avec \(\tensTwo{B}>0\)	Imaginaire pure ou par paire \(\left( \lambda ,-\tilde{\lambda} \right)\)	Si \(\vector{u}\) vecteur propre à droite de \(\lambda\), alors il est aussi vecteur propre à gauche de \(-\tilde{\lambda}\)	Structure libre ou bloquée par élimination, avec gyroscopie
8	\(\tensTwo{A}\) et \(\tensTwo{B}\) symétriques réelles strictement positives et \(\tensTwo{C}=-\transpose{{\tensTwo{C}}}\) (anti-symétrique)	Imaginaire pure

Les figures suivantes illustrent, sur des cas canoniques, certains cas de figure.

Fig. 167 Exemples du spectre d’un QEP \((n=8)\) d’un modèle simplifié d’enceinte nucléaire [TM01]. Figure de gauche, on est dans le cas n°3 du tableau. Celle de droite, on rajoute à \(\tensTwo{K}\) de l’amortissement hystérétique donc apparaissent des modes dépareillés (cas général n°1)

Fig. 168 Exemples du spectre d’un QEP \((n=50)\) d’un système masses-ressorts amorti[TM01]. Figure de gauche, on est dans le cas n°5 du tableau. Celle de droite, on rajoute la contrainte de sur-amortissement et le spectre devient réel négatif en deux paquets (cas n°6).

Spécificités des QEP et éléments de théorie

Une spécificité algébrique majeure entre les QEP et les autres classes de problèmes modaux de niveau inférieure (GEP et SEP) est l’existence de deux fois plus de modes propres (\(2n\)) que la taille du problème discrétisé. Les vecteurs propres ne peuvent plus être linéairement indépendants, de surcroît, les valeurs propres peuvent finies ou infinies et tout «ce petit monde» vit généralement dans le plan complexe [3]_. Il faut alors définir des heuristiques robustes pour filtrer les valeurs propres souhaitées par l’utilisateur et distinguer les réelles, des complexes, les complexes dépareillés de ceux couplés…

Une autre complication inhérente au QEP est d’ordre algorithmique. Il n’existe pas de décomposition de Schur (resp. Schur généralisé) comme pour les SEP (resp. GEP) sur laquelle va pouvoir s’appuyer l’algorithme de résolution. Par exemple, pour le SEP \(\tensTwo{A}\vector{u}=\lambda \vector{u}\), cette décomposition nous assure l’existence d’une matrice unitaire (donc bien conditionnée et facilement inversible) \(\tensTwo{U}\), permettant la réécriture de la matrice de travail \(\tensTwo{A}\) sous une forme plus facile à manipuler [4]: la matrice triangulaire supérieure \(\tensTwo{T}\) .

\[\tensTwo{U}\tensTwo{A}{\tensTwo{U}}^{\star} = \tensTwo{T}\]

Les QEP constituent donc une classe très particulière et importante [5] des problèmes non linéaires. Leur résolution est moins routinière que pour les autres classes de problèmes modaux. En particulier, peu de méthodes permettent de résoudre le problème directement, il faut souvent transiter par un GEP ( via une linéarisation ad hoc ) puis un SEP (transformation spectrale) ce qui induit une perte d’information spectrale et des instabilités numériques (propagation d’erreurs d’arrondis).

Dans la littérature théorique, plutôt qu’au QEP (2.1-1), on s’intéresse au polynôme matriciel (quadratique en \(\lambda\) )

\[\tensTwo{Q}(\lambda) := \left( {\lambda}^{2}\tensTwo{B} + \lambda \tensTwo{C} + \tensTwo{A} \right )\]

On l’appelle aussi \(\lambda\)-matrice ou matrice dynamique et on cherche son spectre noté \(\spectre(\tensTwo{Q}):=\lbrace \lambda \in C/det\tensTwo{Q}(\lambda)=0\rbrace\). Lorsque ce déterminant est nul quel que soit \(\lambda\), la \(\lambda\)-matrice est dite singulière. Dans le cas contraire, elle est dite régulière (on se place en général dans ce cadre).

Le polynôme caractéristique du problème s’écrit donc sous la forme \(\det\tensTwo{Q}(\lambda)=\det(\tensTwo{B}){\lambda}^{2n}\) plus des termes d’ordre inférieur en \(\lambda\). Ainsi, dès que \(\tensTwo{B}=\tensTwo{M}\) est singulière, ce polynôme admet \(r<2n\) racines finies auxquelles il faut rajouter \(2n-r\) infinies. D’autre part, ces valeurs propres distinctes peuvent partager les mêmes vecteurs propres. C’est la «situation duale» des valeurs propres multiples du SEP où une même valeur propre est rattachée à plusieurs vecteurs propres. Ainsi, pour le QEP canonique

\[\begin{split}\tensTwo{A} = {\tensTwoUnit}_{3} \tensTwo{B} = \left[\begin{array}{ccc} 0& 6& 0\\ 0& 6& 0\\ 0& 0& 1 \end{array}\right] \text{ et } \tensTwo{C} = \left[\begin{array}{ccc} 1& -6& 0\\ 2& -7& 0\\ 0& 0& 0 \end{array}\right]\end{split}\]

on a le spectre suivant

\[\begin{split}\begin{array}{ccccccc} i& 1& 2& 3& 4& 5& 6\\ {\lambda}_{i}& 1/3& 1/2& 1& i& -i& \infty \\ {\vector{u}}_{i} & \left[\begin{array}{c}1\\ 1\\ 0\end{array}\right]& \left[\begin{array}{c}1\\ 1\\ 0\end{array}\right]& \left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]& \left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]& \left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]& \left[\begin{array}{c}1\\ 0\\ 0\end{array}\right] \end{array}\end{split}\]

De la singularité de \(\tensTwo{B}\) provient la valeur propre infinie. On note aussi le partage de vecteurs propres par les deux premières valeurs propres et les quatrième et cinquième.

Pour appréhender plus facilement les valeurs propres infinies de \(\tensTwo{Q}(\lambda)\), on les associe aux valeurs propres nulles du polynôme inverse

\[\rev \tensTwo{Q}(\lambda) := {\lambda}^{2} \tensTwo{Q}(\frac{1}{\lambda}) = \left( \tensTwo{B}+\lambda \tensTwo{C}+{\lambda}^{2}\tensTwo{A} \right)\]

Ces valeurs propres finies de \(\rev\tensTwo{Q}(\lambda)\) sont plus faciles à calculer et partagent les mêmes caractéristiques spectrales (multiplicités géométrique et algébrique, espace propre,…) que les valeurs infinies initiales.

Pour de plus amples informations sur ces aspects théoriques on pourra consulter, par exemple, les ouvrages de P.Lancaster [GLR82][LT85] et la thèse de D.S.Mackey[Mac06] (et ses nombreuses références bibliographiques).

Stratégie de linéarisation

Introduction

L’approche standard pour résoudre un problème modal est de le transformer pour faire apparaître des matrices canoniques révélant les modes recherchés. Par exemple, la décomposition de Schur pour les SEP ou celle généralisée pour les GEP. Malheureusement, cette pratique ne peut pas se généraliser aux problèmes modaux d’ordre supérieur. Deux solutions s’offrent alors:

• Résoudre le problème non linéaire directement (racine de polynôme, factorisation de \(\lambda\)-matrice, méthodes de type Newton,…).

• Le linéariser en un GEP, puis traiter ce dernier via les méthodes adaptées dans [Solveurs modaux et résolution du problème généralisé].

La première stratégie est déclinée dans les options PROCHE et AJUSTE (méthode de Müller-Traub) pour capturer une première estimation des valeurs propres. On la panache ensuite avec la stratégie suivante, car pour affiner les estimations et accélérer la convergence, on fait appel à une méthode de type puissances inverses (variante de Jennings) sur le QEP linéarisé. Certains auteurs proposent des approches directes complètes. Elles s’appuient souvent sur des variantes de Newton. Mais leur convergence est lente (un mode à la fois) et pas toujours assurée (cf. Kublanovskaya 1970, G.Peters & J.H.Wilkinson 1979 et A.Ruhe 1973).

La stratégie de linéarisation est donc souvent privilégiée. Elle double la taille du problème modal à traiter, change sa nature (au risque de perdre de l’information spectrale [6]) et introduit des instabilités numériques (sensibilités aux erreurs d’arrondis [7] et aux paramètres numériques [8]). Mais elle présente l’immense avantage de réutiliser toute «l’artillerie numérique» déjà déployée pour les GEP [9].

On peut ainsi capturer simultanément tout le spectre (cf. méthode QZ utilisable avec les options CENTRE, PLUS_PETITE et TOUT) ou seulement le partie centrée sur une zone d’intérêt (cf. IRAM ou Lanczos) en recourant à la transformation spectrale idoine.

Principe

Maintenant que le contexte a été brossé, pour se conforter à la terminologie employée dans la littérature, on va désormais simplifier les notations en considérant le QEP associé à la \(\lambda\)-matrice \(\tensTwo{Q}(\lambda):=({\lambda}^{2}\tensTwo{M}+\lambda \tensTwo{C}+\tensTwo{K})\) . Cette matrice est aussi appelée «matrice dynamique» dans le jargon du dynamicien des structures . C’est le pendant de celle décrite dans la note dédiée au GEP: \(\tensTwo{Q}(\lambda):=(\tensTwo{K}-\lambda \tensTwo{M})\) .

Définition 1

D’un point de vue théorique, linéariser \(\tensTwo{Q}(\lambda)\) revient à trouver une autre \(\lambda\)-matrice (complexe carrée de taille \(2n\) ) du type \(\tensTwo{A}-\lambda \tensTwo{B}\) telle qu’il existe deux autres \(\lambda\)-matrices (\(\tensTwo{E}(\lambda)\) et \(\tensTwo{F}(\lambda)\) ) à déterminant constant non nul, vérifiant la relation

\[\begin{split}\left[\begin{array}{cc} \tensTwo{Q}(\lambda)& {0}_{n}\\ {0}_{n}& {\tensTwoUnit}_{n} \end{array}\right] = \tensTwo{E}(\lambda)(\tensTwo{A}-\lambda \tensTwo{B})\tensTwo{F}(\lambda)\end{split}\]

Les matrices \(\tensTwo{A}\) et \(\tensTwo{B}\) sont appelées les matrices «compagnon» de la linéarisation du QEP.

Clairement, les spectres du QEP initial \((\tensTwo{Q}(\lambda))\) et celui de sa linéarisation \((\tensTwo{A}-\lambda \tensTwo{B})\) coïncident. Cependant cette décomposition n’est pas unique et il faut choisir, si possible, celle qui:

• Préserve les propriétés des matrices initiales (arithmétique, symétrie,…),

• Présente le moins de sensibilité aux erreurs d’arrondis («backward stable»: matrices équilibrées, bien conditionnées, inversibles,…),

• Manipule des matrices canoniques(identité, triangulaire,…),

• S’adapte à un large périmètre de QEP (amortissement visqueux, gyroscopie, sur-amortissement,…).

Une manière plus intuitive d’appréhender les techniques de linéarisation, consiste à introduire un changement de variable du type \(\vector{v}=\lambda \vector{u}\) et à réarranger sous forme matricielle:

\[\begin{split}\left( \underbrace{{\left[\begin{array}{cc}{0}_{n}& {\tensTwoUnit}_{n}\\ -\tensTwo{K}& -\tensTwo{C}\end{array}\right]}}_{\tensTwo{A}} - \lambda \underbrace{{\left[\begin{array}{cc}{\tensTwoUnit}_{n}& {0}_{n}\\ {0}_{n}& \tensTwo{M}\end{array}\right]}}_{\tensTwo{B}} \right) \left[\begin{array}{c} \vector{u}\\ \vector{v} \end{array}\right] = {0}_{2n}\end{split}\]

Cette linéarisation «très classique» rentre dans le cadre de la définition 1 si on utilise les \(\lambda\)-matrices :

\[\begin{split}\tensTwo{E}(\lambda) := \left[\begin{array}{cc} -(\tensTwo{C}+\lambda \tensTwo{M})& -{\tensTwoUnit}_{n}\\ {\tensTwoUnit}_{n}& {0}_{n} \end{array}\right] \text{ et } \tensTwo{F}(\lambda) := \left[\begin{array}{cc} {\tensTwoUnit}_{n}& {0}_{n}\\ \lambda {\tensTwoUnit}_{n}& {\tensTwoUnit}_{n} \end{array}\right]\end{split}\]

En fait, la linéarisation d’un QEP épouse le cadre général de celle des polynômes matriciels. En posant \(\tensTwo{Q}(\lambda):=\sum_{i=1}^{k}{\lambda}_{i}{\tensTwo{X}}_{i}\) (ici \(k=2,{\tensTwo{X}}_{i}=\tensTwo{M}\) etc), on est toujours assuré de trouver deux linéarisations, dites «compagnons», de la forme \({\tensTwo{A}}_{1}-\lambda {\tensTwo{B}}_{1}\) et \({\tensTwo{A}}_{2}-\lambda {\tensTwo{B}}_{2}\) avec

\[\begin{split}{\tensTwo{B}}_{1} = {\tensTwo{B}}_{2} := \left[\begin{array}{cc} -{\tensTwo{X}}_{k} & {0}_{(k-1)n}\\ {0}_{(k-1)n}& -{\tensTwoUnit}_{(k-1)n} \end{array}\right]\end{split}\]

\[\begin{split}{\tensTwo{A}}_{1} := \left[\begin{array}{cccc} {\tensTwo{X}}_{k-1}& {\tensTwo{X}}_{k-2}& \cdots & {\tensTwo{X}}_{0}\\ {\tensTwoUnit}_{n}& {0}_{n}& \cdots & {0}_{n}\\ \cdots & \cdots & \cdots & \cdots\\ {0}_{n}& \cdots & -{\tensTwoUnit}_{n}& {0}_{n} \end{array}\right]\end{split}\]

\[\begin{split}{\tensTwo{A}}_{2} := \left[\begin{array}{cccc} {\tensTwo{X}}_{k-1}& -{\tensTwoUnit}_{n}& \cdots & {0}_{n}\\ {\tensTwo{X}}_{k-2}& {0}_{n}& \cdots & \cdots\\ \cdots & \cdots & \cdots & -{\tensTwoUnit}_{n}\\ {\tensTwo{X}}_{0}& {0}_{n}& \cdots & {0}_{n} \end{array}\right]\end{split}\]

En pratique, on leur préfère des linéarisations plus adaptées au cas traité (cf. paragraphe suivant).

Quelques exemples

Dans la littérature, il est souvent question des linéarisations suivantes (avec \(\tensTwo{N}\) une matrice régulière carrée complexe de taille \(n\) ):

\[\begin{split}{L}_{1} : \left( \underbrace{ \left[\begin{array}{cc} {0}_{n} & \tensTwo{N}\\ -\tensTwo{K} & -\tensTwo{C} \end{array}\right] }_{\tensTwo{A}} - \lambda \underbrace{ \left[\begin{array}{cc} \tensTwo{N}& {0}_{n}\\ {0}_{n}& \tensTwo{M} \end{array}\right] }_{\tensTwo{B}} \right) \left[\begin{array}{c} \vector{u}\\ \lambda \vector{u} \end{array}\right] = {0}_{2n}\end{split}\]

\[\begin{split}{L}_{2} : \left( \underbrace{ \left[\begin{array}{cc} -\tensTwo{K}& {0}_{n}\\ {0}_{n}& \tensTwo{N} \end{array}\right] }_{\tensTwo{A}} - \lambda \underbrace{ \left[\begin{array}{cc} \tensTwo{C}& \tensTwo{M}\\ \tensTwo{N}& {0}_{n} \end{array}\right] }_{\tensTwo{B}} \right) \left[\begin{array}{c} \vector{u}\\ \lambda \vector{u} \end{array}\right] = {0}_{2n}\end{split}\]

Le choix dépend de la singularité potentielle des matrices \(\tensTwo{K}\) et \(\tensTwo{M}\) . On souhaite pouvoir inverser les termes diagonaux, donc si \(\tensTwo{K}\) est singulière (resp. \(\tensTwo{M}\) ), on utilise \({L}_{1}\) (resp. \({L}_{2}\)). Pour équilibrer les termes des matrices et faciliter les manipulations, on prend souvent \(\tensTwo{N}=\alpha {\tensTwoUnit}_{n}\) avec \(\alpha\) facteur réel d’équilibrage (cf. [Solveurs modaux et résolution du problème généralisé], annexe 1) construit à partir des normes matricielles des autres termes \(\alpha =\alpha \left( \Vert \tensTwo{K}\Vert ,\Vert \tensTwo{M}\Vert ,\Vert \tensTwo{C}\Vert \right)\) . Pour travailler avec un GEP symétrique (voire défini ou semi-défini positif), on peut aussi initialiser la matrice auxiliaire \(\tensTwo{N}\) à \(\alpha \tensTwo{K}\) ou \(\beta \tensTwo{M}\) avec \((\alpha , \beta) \in \setR\).

Pour résoudre plus efficacement les QEP avec effets gyroscopiques, on conseille dans la littérature les pendants «Hamiltoniens [11]» des linéarisations précédentes. Soient les matrices «compagnons»

\[\begin{split}{L}_{3} : \left( \underbrace{ \left[\begin{array}{cc} \tensTwo{K}& {0}_{n}\\ \tensTwo{C}& \tensTwo{N} \end{array}\right] }_{\tensTwo{A}} - \lambda \underbrace{ \left[\begin{array}{cc} {0}_{n}& \tensTwo{K}\\ -\tensTwo{M}& {0}_{n} \end{array}\right] }_{\tensTwo{B}} \right) \left[\begin{array}{c} \vector{u}\\ \lambda \vector{u} \end{array}\right] = {0}_{2n}\end{split}\]

\[\begin{split}{L}_{4} : \left( \underbrace{ \left[\begin{array}{cc} {0}_{n}& -\tensTwo{K}\\ \tensTwo{M}& \tensTwo{N} \end{array}\right] }_{\tensTwo{A}} - \lambda \underbrace{ \left[\begin{array}{cc} \tensTwo{M}& \tensTwo{C}\\ {0}_{n}& \tensTwo{M} \end{array}\right] }_{\tensTwo{B}} \right) \left[\begin{array}{c} \vector{u}\\ \lambda \vector{u} \end{array}\right] = {0}_{2n}\end{split}\]

Ici l’objectif n’est pas de respecter une propriété de symétrie mais plutôt celles «Hamiltoniennes [12] » du problème. Si \(\tensTwo{K}\) est singulière (resp. \(\tensTwo{M}\)), on utilise \({L}_{4}\) (resp. \({L}_{3}\)).

Dans le cas d’un QEP dit «T-palindromique» [Mac06] \((\tensTwo{K}=\transpose{\tensTwo{M}}\) et \(\tensTwo{C}=\transpose{\tensTwo{C}}\), on préfère la linéarisation

\[\begin{split}{L}_{3} : \left( \underbrace{ \left[\begin{array}{cc} \tensTwo{K} & \tensTwo{K}\\ \tensTwo{C}-\tensTwo{M} & \tensTwo{K} \end{array}\right] }_{\tensTwo{A}} - \lambda \underbrace{ \left[\begin{array}{cc} -\tensTwo{M} & \tensTwo{K}-\tensTwo{C} \\ -\tensTwo{M} & -\tensTwo{M} \end{array}\right] }_{\tensTwo{B}} \right) \left[\begin{array}{c} \vector{u}\\ \lambda \vector{u} \end{array}\right] = {0}_{2n}\end{split}\]

(\({L}_{5}\) ) \((\underset{\tensTwo{A}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\tensTwo{K}& \tensTwo{K}\\ \tensTwo{C}-\tensTwo{M}& \tensTwo{K}\end{array}\right]}}-\lambda \underset{\tensTwo{B}}{\underset{\underbrace{}}{\left[\begin{array}{cc}-\tensTwo{M}& \tensTwo{K}-\tensTwo{C}\\ -\tensTwo{M}& -\tensTwo{M}\end{array}\right]}})\left[\begin{array}{c}\vector{u}\\ \lambda \vector{u}\end{array}\right]={0}_{2n}\) (2.4-9)

Stratégies retenues

On a donc une très grande variété de choix pour adapter la linéarisation du QEP à chaque situation. On ne propose pas à l’utilisateur de paramètre permettant de moduler uniquement la linéarisation. Celle-ci est figée pour un solveur modal donné:

• Avec OPTION parmi PROCHE et AJUSTE (seconde phase) le passage QEP/GEP s’effectue via les matrices compagnons de \({L}_{1}\), en posant \(\tensTwo{N}={\tensTwoUnit}_{n}.\) Le système linéarisé n’est pas symétrique mais ce n’est pas un pré-requis pour appliquer une méthode de type puissances inverses.

\[\begin{split}{L}_{1}^{\prime} : \left( \underbrace{ \left[\begin{array}{cc} {0}_{n}& {\tensTwoUnit}_{n}\\ -\tensTwo{K}& -\tensTwo{C} \end{array}\right] }_{\tensTwo{A}} - \lambda \underbrace{ \left[\begin{array}{cc} {\tensTwoUnit}_{n}& {0}_{n}\\ {0}_{n}& \tensTwo{M} \end{array}\right] }_{\tensTwo{B}} \right) \left[\begin{array}{c} \vector{u}\\ \lambda \vector{u} \end{array}\right] = {0}_{2n}\end{split}\]

• Avec OPTION parmi CENTRE et PLUS_PETITE et METHODE=’TRI_DIAG' ou METHODE=’SORENSEN', le passage QEP/GEP s’effectue via les matrices compagnons de \({L}_{2}\) en posant \(\tensTwo{N}=\tensTwo{M}\). Si les matrices de rigidité, de masse et d’amortissement sont symétriques réelles, on retombe dans le périmètre des GEP standards (symétriques réels). Les méthodes de Lanczos et de Sorensen sont alors accessibles. Si math:tensTwo{K} devient complexe ou si l’une d’entre elles est non symétrique, seul IRAM reste en course (avec QZ cf. item suivant).

\[\begin{split}{L}_{2}^{\prime} : \left( \underbrace{ \left[\begin{array}{cc} -\tensTwo{K}& {0}_{n}\\ {0}_{n}& \tensTwo{M} \end{array}\right] }_{\tensTwo{A}} - \lambda \underbrace{ \left[\begin{array}{cc} \tensTwo{C}& \tensTwo{M} \\ \tensTwo{M}& {0}_{n} \end{array}\right] }_{\tensTwo{B}} \right) \left[\begin{array}{c} \vector{u}\\ \lambda \vector{u} \end{array}\right] = {0}_{2n}\end{split}\]

Toutefois, lorsque \(\tensTwo{M}\) n’est pas inversible (du fait par exemple de Lagrange de blocage), l’inversion du terme diagonal de \(\tensTwo{A}\) n’est plus assurée. Pour pallier ce problème, on introduit une «régularisation» de la matrice \(\tensTwo{M}\), notée \(\tensTwo{M}_{R}\) , dans le terme diagonal. Celle-ci a la propriété d’annuler, par multiplication, toutes les composantes du noyau de \(\tensTwo{M}\): \({\tensTwo{M}}_{R}\vector{u}=0\text{ si }\vector{u}\in \kernel(\tensTwo{M})\) . Ainsi, on peut définir \(\inverse{{\tensTwo{M}}_{R}}\) qui vérifie formellement la propriété \(\inverse{{\tensTwo{M}}_{R}} \tensTwo{M}\approx \tensTwo{M}\inverse{{\tensTwo{M}}_{R}} \approx {\tensTwoUnit}_{n}\) sur les composantes «hors noyau de \(\tensTwo{M}\) » des vecteurs manipulés. Du fait de la structure particulière de la matrice dualisée \(\tensTwo{M}\) (cf. [Solveurs modaux et résolution du problème généralisé]), cette manipulation revient à proscrire les composantes liées aux Lagrange et aux degrés de liberté bloqués (les \({\vector{v}}_{i}^{2,3}\) de la démonstration de la propriété n°4 [Solveurs modaux et résolution du problème généralisé]). Ce qui ne gêne en rien le processus puisqu’on ne s’intéresse qu’aux modes physiques et non à ces artefacts de modélisations [13] .

\[\begin{split}{L}_{2}^{\prime \prime} : \left( \underbrace{ \left[\begin{array}{cc} -\tensTwo{K}& {0}_{n}\\ {0}_{n}& {\tensTwo{M}}_{R} \end{array}\right] }_{\tensTwo{A}} - \lambda \underbrace{ \left[\begin{array}{cc} \tensTwo{C}& \tensTwo{M}\\ \tensTwo{M}& {0}_{n} \end{array}\right] }_{\tensTwo{B}} \right) \left[\begin{array}{c} \vector{u}\\ \lambda \vector{u} \end{array}\right] = {0}_{2n}\end{split}\]

• Avec OPTION parmi CENTRE, PLUS_PETITE et TOUT et avec METHODE=’QZ’, on a préféré éviter ces complications. L’objectif étant de privilégier la robustesse numérique pour traiter de petits cas (moins de mille degrés de liberté). Le scénario \({L}_{2}\) a été empiriquement retenu en posant \(\tensTwo{N}=\alpha {\tensTwoUnit}_{n}\) avec \(\alpha :=\frac{1}{3n}(\Vert \tensTwo{K}\Vert +\Vert \tensTwo{M}\Vert +\Vert \tensTwo{C}\Vert )\)

\[\begin{split}{L}_{2}^{\star} : \left( \underbrace{ \left[\begin{array}{cc} -\tensTwo{K}& {0}_{n}\\ {0}_{n}& \alpha {\tensTwoUnit}_{n} \end{array}\right] }_{\tensTwo{A}} - \lambda \underbrace{ \left[\begin{array}{cc} \tensTwo{C}& \tensTwo{M}\\ \alpha {\tensTwoUnit}_{n}& {0}_{n} \end{array}\right] }_{\tensTwo{B}} \right) \left[\begin{array}{c} \vector{u}\\ \lambda \vector{u} \end{array}\right] = {0}_{2n}\end{split}\]

Implantation

Le déroulement d’un calcul modal quadratique est, par construction, très proche de celui d’un problème généralisé. On va donc se rattacher au schéma fonctionnel de résolution des GEP décrit dans [Solveurs modaux et résolution du problème généralisé]. Ses étapes se décomposent comme suit:

Pré-traitements (linéarisation, calcul du shift)

Détermination de la forme linéarisée et de la transformation spectrale associée suivant la méthode retenue (puissance, Lanczos, IRAM et QZ) et l’option de recherche AJUSTE``PROCHE`` ou PLUS_PETITE/CENTRE.

A noter que les modes propres étant dans le plan complexe, le test de Sturm actuellement implanté est inopérant. Donc l’heuristique associée à l’option PROCHE et AJUSTE n’est pas basée sur les positions modales, et l’option BANDE est proscrite. De même pour la post-vérification basée sur le test de Sturm [Boi12].

Pour les méthodes de sous-espace (OPTION parmi CENTRE et PLUS_PETITE et METHODE='TRI_DIAG' ou METHODE='SORENSEN'), une fois le problème linéarisé \({L}_{2}^{\prime \prime}\) mis à jour, on applique une transformation spectrale de type «shift-and-invert» (cf. [Solveurs modaux et résolution du problème généralisé]). Cela va permettre d’achever la transformation en un SEP [14] \({S}_{1}\) , d’orienter la recherche du spectre dans une zone d’intérêt pour l’utilisateur et de mieux séparer les valeurs propres.

\[\begin{split}{S}_{1} : \underbrace{\inverse{\left( \tensTwo{A}-\sigma \tensTwo{B} \right)}}_{\tensTwo{A}_{\sigma}} \tensTwo{B} \underbrace{ {\left[\begin{array}{c} \vector{u}\\ \lambda \vector{u} \end{array}\right]} }_{\tensTwo{W}} = \underbrace{ \frac{1}{\lambda -\sigma} }_{\mu} \underbrace{ {\left[\begin{array}{c} \vector{u}\\ \lambda \vector{u} \end{array}\right]} }_{\tensTwo{W}}\end{split}\]

La méthode des puissances inverses applique elle une transformation spectrale particulière au problème linéarisé \({L}_{1}^{\prime}\) .

Quant à la méthode QZ (OPTION parmi CENTRE, PLUS_PETITE et TOUT et METHODE='QZ'), cette transformation spectrale est superflue. L’algorithme traite directement le problème linéarisé \({L}_{2}^{\star}\) et calcule tous ses modes. Suivant les souhaits de l’utilisateur (OPTION , CALC_FREQ=_F(NMAX_FREQ,...)), il filtre ensuite les résultats.

Remarques

• Comme pour les techniques de linéarisation, aucun paramètre utilisateur ne permet encore de changer de transformation spectrale. En quadratique, il serait intéressant d’utiliser des transformations plus adaptées au plan complexe que le traditionnel «shift-and-invert». Par exemple, une transformation de type Cayley pour faire converger simultanément les couples de valeurs propres conjuguées \(\left( \lambda, \tilde{\lambda} \right)\)

\[{S}_{1} : \underbrace{ \left( \tensTwo{A} - \sigma \tensTwo{B} \right) \inverse{\tensTwo{B}} \left( \tensTwo{A}-\tilde{\sigma}\tensTwo{B} \right) }_{\tensTwo{A}_{\sigma}} \inverse{\tensTwo{B}} \tensTwo{W} = \underbrace{ \frac{1}{(\lambda-\sigma)(\tilde{\lambda}-\sigma)} }_{\mu} {\tensTwo{W}}\]

• Cette stratégie de Cayley se généralise pour traiter la gyroscopie (cf. V.Mehrman et D.Watkins 2001) et faire converger simultanément les quadruplets \((\lambda, \tilde{\lambda}, -\lambda , -\tilde{\lambda})\)

• Certains auteurs proposent aussi d’intervertir la phase de linéarisation et celle de transformation spectrale [Mee07].

Factorisation de matrices dynamiques

Comme on l’a déjà souligné, ces factorisations de matrices dynamiques du type \(\tensTwo{Q}(\lambda)\) ne servent pas à mettre en œuvre un test de Sturm. Ce dernier étant illicite pour les QEP. Cependant elles restent un ingrédient très présent et très coûteux [15] pour:

• Estimer la fonction coût de la partie heuristique (méthode de Müller-Traub) avec OPTION=’AJUSTE'.

• Lorsqu’on doit manipuler une matrice de travail comportant des inverses. Pour être plus efficace, on cherche alors une formulation ne faisant intervenir que \(\inverse{\tensTwo{Q}(\lambda)}\) (OPTION parmi PROCHE et AJUSTE et SOLVEUR_MODAL=_F(OPTION_INV+’DIRECT'), ou bien OPTION parmi CENTRE et PLUS_PETITE et SOLVEUR_MODAL=_F(METHODE=’TRI_DIAG'/’SORENSEN')).

L’approche globale (OPTION parmi CENTRE, PLUS_PETITE et TOUT et SOLVEUR_MODAL=_F(METHODE=’QZ')) n’est pas concernée par cette factorisation préliminaire de la matrice dynamique.

Calcul modal

Le calcul modal proprement dit est effectué: on résout le problème standard SEP, on revient au GEP intermédiaire, puis au QEP initial. A chacune de ces conversions, on filtre et transcrit les résultats du calcul précédent. Concernant les valeurs propres, on ne retient que les \(\lambda\) à partie imaginaire positive des modes couplés \(\left( \lambda, \tilde{\lambda} \right).\)

Exemple 1. Impressions dans le fichier message d’alarmes récapitulant le nombre de modes non retenus (ici deux réels). Extrait du cas-test sdll123a.:

votre problème est fortement amorti.

valeur(s) propre(s) réelle(s) : 2

valeur(s) propre(s) complexe(s) avec conjuguée : 107

valeur(s) propre(s) complexe(s) sans conjuguée : 0

Pour les vecteurs propres, deux possibilités restent ouvertes:

• Prendre directement les \(n\) premières composantes de \(w\) (partie dite «haute»),

• Choisir les \(n\) dernières (partie dite «basse») éventuellement divisées par le même scalaire.

C’est la deuxième option qui a été retenue. Ce choix n’est pas toujours anodin [16] notamment vis-à-vis de la qualité des modes et de leurs sensibilités aux erreurs d’arrondis.

Avec OPTION parmi PROCHE et AJUSTE, une seule méthode est disponible, une variante de la méthode des itérations inverses due à Jennings (OPTION_INV=‘DIRECT'). Elle affine les valeurs propres préalablement détectées par la méthode de Müller-Traub (AJUSTE) ou les estimations fournies par l’utilisateur (PROCHE). Pour ce qui est de OPTION parmi CENTRE, PLUS_PETITE et TOUT, on peut utiliser trois méthodes distinctes: la méthode de Lanczos (TRI_DIAG ), celle de Sorensen (SORENSEN) et QZ (QZ). Seules les deux dernières sont disponibles avec des matrices non symétriques réelles ou complexes symétriques (pour \(\tensTwo{K}\) seulement).

Remarques

• Chacune de ces méthodes possède des tests d’arrêts internes. Sans compter que les méthodes de projection emploient des méthodes modales auxiliaires (QR/QL) pour Lanczos et IRAM. Elles nécessitent aussi des tests d’arrêts. L’utilisateur a souvent accès à ces paramètres, bien qu’il soit chaudement recommandé, au moins dans un premier temps, de conserver leurs valeurs par défaut.

• Contrairement au GEP, on ne force pas la \(\tensTwo{M}\)-orthonormalisation des vecteurs propres. Elle n’a aucun sens en QEP. De plus, en appliquant le résultat de la proposition 2 de [Solveurs modaux et résolution du problème généralisé] au problème linéarisé, il est clair que les vecteurs propres sont \(\tensTwo{A}\) et \(B\)-orthogonaux. Mais du fait de la complexité des réductions linéaires utilisées, cela ne conduit pas (sauf cas particulier) à des vecteurs propres \(\tensTwo{K}\), \(\tensTwo{M}\) ou \(\tensTwo{C}\)-orthogonaux.

Post-traitements de vérification

Cette dernière partie regroupe le seul post-traitement qui vérifie le bon déroulement du calcul (en l’absence de test de Sturm adapté au cas complexe). On calcule l’erreur relative sur la norme [17] du résidudu QEP initial.

1. \(\vector{u} \Leftarrow \frac{\vector{u}}{{\Vert\vector{u}\Vert}_{\infty}}\)

2. Si \(\mid \lambda \mid >\text{SEUIL\_FREQ}\), alors \(\frac{{\Vert{\lambda}^{2}\tensTwo{M} \vector{u}+\lambda \tensTwo{C} \vector{u}+\tensTwo{K} \vector{u}\Vert}_{2}}{{\Vert\tensTwo{K} \vector{u}\Vert}_{2}} < \text{SEUIL}\)

3. Sinon \({\Vert{\lambda}^{2}\tensTwo{M} \vector{u}+\lambda \tensTwo{C} \vector{u}+\tensTwo{K} \vector{u}\Vert}_{2} <\text{SEUIL}\)

Cette séquence est paramétrée par les mot-clés SEUIL et SEUIL_FREQ, appartenant respectivement aux mot-clés facteurs VERI_MODE et CALC_FREQ. D’autre part, ce post-traitement est activé par l’initialisation à 'OUI' (valeur par défaut) de STOP_ERREUR dans le mot-clé facteur VERI_MODE. Lorsque cette règle est activée et non-respectée, le calcul s’arrête, sinon l’erreur est juste signalée par une alarme. On resommande de ne pas désactiver ce paramètre.

Affichage dans le fichier message

Dans le fichier message sont mentionnées les informations relatives aux modes retenus. Par exemple, dans le cas d’un QEP, on précise le solveur modal utilisé et la liste des valeurs propres \({\lambda}_{i}\) retenues par ordre croissant de partie imaginaire \(\frac{\Im({\lambda}_{i})}{2\pi }\) (colonne FREQUENCE). Puis on précise son amortissement \({\xi}_{i}=\frac{-\Re({\lambda}_{i})}{{\lambda}_{i}}\) et sa norme d’erreur.

Exemple 2. Impressions dans le fichier message des valeurs propre s retenues lors de QEP. Extrait du cas-test sdll123a:

CALCUL MODAL: METHODE GLOBALE DE TYPE QR

ALGORITHME **QZ_EQUI**

NUMERO **FREQUENCE (HZ)** **AMORTISSEMENT** **NORME D'ERREUR**

1 1.23915E+02 1.55604E-08 9.93686E-02

2 1.24546E+02 -3.76772E-10 1.81854E-01

3 4.97033E+02 -5.88927E-10 1.41043E-03

4 4.99575E+02 1.41315E-12 2.14401E-03

5 1.11531E+03 4.75075E-12 1.17669E-04

Maintenant que le contexte des QEP a été brossé, nous allons nous intéresser plus particulièrement aux méthodes de type puissance et à leur implantation dans l’OPTION PROCHE et AJUSTE.

[2]

Condition de sur-amortissement: \(\gamma (\tensTwo{A},\tensTwo{B},\tensTwo{C}):=\underset{{\Vert \vector{x}\Vert }_{2}=1}{min}\left[{({\vector{x}}^{\star}\tensTwo{C}\vector{x})}^{2}-4({\vector{x}}^{\star}\tensTwo{M}\vector{x})({\vector{x}}^{\star}\tensTwo{K}\vector{x})\right]>0\).

[3]

C’est déjà le cas pour certaines classes de GEP: non symétrique ou à matrice complexe.

[4]

Les valeurs propres du problème se retrouvent sur la diagonale principale de \(\tensTwo{T}\) et on peut déduire facilement des vecteurs de Schur (les vecteurs colonnes de \(\tensTwo{Q}\)), leurs vecteurs propres associés.

[5]

Du fait notamment de ces applications croissantes dans l’industrie: analyse vibratoire, vibro-accoustique, mécanique des fluides avec turbulence, électronique (VLSI),…

[6]

On trouve la solution exacte d’un «GEP approché» mais celle-ci ne correspond pas forcément à la solution exacte d’un «QEP approché» car les solveurs GEP ne respectent pas la structure particulière du problème linéarisé (cf. [TM01] §5).

[7]

L’application de solveur modaux de GEP à un QEP linéarisé n’est pas «backward stable». C’est-à-dire que cette technique a une propension à transmettre/amplifier les erreurs d’arrondis.

[8]

Les résultats peuvent être sensibles aux différents paramètres des solveurs modaux utilisés in fine, mais aussi et surtout, à ceux des deux pré-traitements: linéarisation et transformation spectrale.

[9]

On a fait de même pour la transition GEP/SEP ([Boi09] §3.7).

[11]

Une matrice carrée d’ordre \(2n\) \(\tensTwo{A}\), est dite Hamiltonienne (resp. anti-Hamiltonienne) si elle vérifie la relation \({(\tensTwo{A} \tensTwo{J})}^{\tensTwo{T}}=\tensTwo{A} \tensTwo{J}\) (resp. \({(\tensTwo{A} \tensTwo{J})}^{\tensTwo{T}}=-\tensTwo{A} \tensTwo{J}\) avec \(\tensTwo{J}\) la matrice unité anti-symétrique \(\tensTwo{J}:=\left[\begin{array}{cc}{0}_{n}& {\tensTwoUnit}_{n}\\ -{\tensTwoUnit}_{n}& {0}_{n}\end{array}\right]\).

[12]

C’est-à-dire que les solutions sont symétriques par rapport aux axes réel et imaginaire:\((\lambda ,\tilde{\lambda},-\lambda ,-\tilde{\lambda}.\)

[13]

Un procédé identique est souvent utilisé pour résoudre les GEP. Il consiste à initialiser l’algorithme de résolution par un vecteur de l’ensemble image de \(\tensTwo{M}\) (filtrant ainsi les éléments du noyaux). Au lieu du vecteur aléatoire \({\vector{u}}_{i}\) on prend \({\tensTwo{M} \vector{u}}_{i}.\)

[14]

Sur lequel s’appliquent les solveurs modaux classiques.

[15]

Complexité calcul en \(O(bn^2)\) avec \(b\) la largeur de bande et \(n\) la taille du problème.

[16]

La variante de la méthode QZ proposée par F.Tisseur [TM01] se base sur les valeurs de «l’erreur inverse» de chaque mode, pour choisir la partie haute ou basse des vecteurs propres.

[17]

Norme hermitienne et uniquement sur les composantes des degrés de liberté physiques.

Méthode des puissances inverses (OPTION parmi [“PROCHE”,”AJUSTE”])

Introduction

Comme pour le traitement des GEP, en QEP cet opérateur fonctionne en deux parties :

• La localisation des valeurs propres . On détermine une approximation de chaque valeur propre contenue dans un intervalle donné par l’algorithme de Müller-Traub [Mul56][DB08] ou on prend une estimation fournie par l’utilisateur.

• L’amélioration de ces estimations et le calcul de leurs vecteurs propres associés par une méthode d’itérations inverses (variante de Jennings[Jen77]).

La recherche d’une valeur approchée pour chaque valeur propre considérée est sélectionnée avec le mot-clé OPTION:

• Si OPTION=”AJUSTE”, dans chaque intervalle de fréquences définies par le mot-clé CALC_FREQ=_F(FREQ), une valeur approchée de chaque valeur propre contenue dans cet intervalle est calculée en utilisant la méthode de Müller-Traub (cf. §3.2).

• Soit on ne retient que les NMAX_FREQ (sous le mot-clé facteur CALC_FREQ) fréquences les plus basses contenues dans l’intervalle maximal spécifié par l’utilisateur, soit on calcule toutes les valeurs de cet intervalle (si NMAX_FREQ=0).

• Si OPTION=”PROCHE”, les fréquences données par le mot-clé CALC_FREQ=_F(FREQ), sont considérées comme les valeurs approchées des valeurs propres cherchées.

Remarques:

• Bien sûr, comme on l’a déjà précisé, ces options sont à utiliser plutôt pour déterminer ou affiner quelques valeurs propres. Pour une recherche plus étendue il faut utiliser les options [“CENTRE”,”PLUS_PETITE”,”TOUT”] .

• Avec l’option “ PROCHE “ on ne peut pas calculer de modes multiples.

• C’est un algorithme coûteux car il fait beaucoup appel aux factorisations d’une matrice dynamique de travail \(Q(\lambda ).\)

Méthode de Müller-Traub

Comme en QEP, l’astucieux test de Sturm ne peut plus s’appliquer, l’heuristique de localisation des valeurs propres doit être basée sur une autre stratégie. On se «rabat» alors sur la classique recherche des zéros du polynôme caractéristique associé au QEP \(p(\lambda )=\text{dét}\mathrm{Q}(\lambda )\) dans le plan complexe.

Comme on ne dispose alors plus de relation d’ordre, il n’est pas possible d’appliquer la dichotomie mise en œuvre pour les GEP. Une généralisation à trois points de la méthode de la sécante permet toutefois de s’en sortir. C’est la méthode proposée par D.E.Muller[Mul56][DB08] en 1956. C’est une méthode itérative utilisant comme courbe d’interpolation une parabole à axe horizontal. Alors que la méthode de la sécante s’appuie sur une interpolation linéaire entre les deux derniers itérés \(\left\lbrace ({\lambda}_{k-1},p({\lambda}_{k-1})),({\lambda}_{k},p({\lambda}_{k}))\right\rbrace ,\) la méthode dite de «Müller-Traub» cherche à construire l’interpolée quadratique passant par les trois derniers

\(\left\lbrace ({\lambda}_{k-2},p({\lambda}_{k-2})),({\lambda}_{k-1},p({\lambda}_{k-1})),({\lambda}_{k},p({\lambda}_{k}))\right\rbrace .\)

En écrivant que l’interpolée quadratique de Lagrange du polynôme caractéristique

\(p(\lambda )\approx {a}_{0}{\lambda}^{2}+{a}_{1}\lambda +{a}_{2}\) (3.2-1)

doit passer par les trois derniers points cités, on construit un système de trois équations à trois inconnues. De sa résolution analytique, on déduit la racine la plus proche du dernier itéré

\(\begin{array}{ccc}{\rho}_{k+1}& =& \frac{-\mathrm{2p}({\lambda}_{k}){\delta}_{k}}{{g}_{k}\pm \sqrt{{g}_{k}^{2}-\mathrm{4p}({\lambda}_{k}){\delta}_{k}{\rho}_{k}\left[p({\lambda}_{k-2}){\rho}_{k}-p({\lambda}_{k-1}){\delta}_{k}+p({\lambda}_{k})\right]}}\\ \text{avec}& & {g}_{k}:=p({\lambda}_{k-2}){\rho}_{k}^{2}-p({\lambda}_{k-1}){\delta}_{i}^{2}+p({\lambda}_{k})({\rho}_{k}+{\delta}_{k})\\ & & {\rho}_{k}:=\frac{{\lambda}_{k}-{\lambda}_{k-1}}{{\lambda}_{k-1}-{\lambda}_{k-2}}\end{array}\) (3.2-2)

De \({\rho}_{k}\) , on déduit facilement alors le nouvel itéré \({\lambda}_{k+1}\) . Le signe du dénominateur est choisi de manière à minimiser \(\mid {\rho}_{k+1}\mid\) ainsi on privilégie le choix «conservatif» du zéro le plus proche de l’itéré précédent \({\lambda}_{k}\) . On considère que la méthode a convergé, lorsque le critère relatif suivant a été atteint (avant \(k+1<\) NMAX_ITER_AJUSTE itérations)

\(\frac{∣{\lambda}_{k+1}-{\lambda}_{k}∣}{∣{\lambda}_{k}∣}<{\varepsilon}_{\text{PREC\_AJUSTE}}\) (3.2-3)

Cette méthode est relativement aisée à mettre en oeuvre mais elle se prête mal aux recherches de zéros de fonctions réelles à racines réelles, car elle plonge l’interpolation dans le plan complexe et ceci même en partant de valeurs réelles. Son intérêt est lié à la classe de cette méthode dite «par courbes d’interpolation», à savoir :

• La sûreté de la méthode de bissection, puisque la recherche s’effectue dans une boule d’approximés \({({\lambda}_{k})}_{k}\) se réduisant progressivement,

• Seul le calcul de la fonction est requis (pas de calcul de dérivée comme par exemple dans la méthode de Newton-Raphson).

• La convergence s’apparente à une convergence quadratique ( q ~1.89).

En pratique, l’évaluation se fait non pas directement sur le déterminant du polynôme \(p(\lambda ):=\text{dét}\mathrm{Q}(\lambda )\) caractéristique mais sur son «déflaté» pour tenir compte de l’information spectrale déjà exhumée (sinon on reconverge toujours vers le même mode). D’autre part, après avoir calculé \(k\) valeurs propres \({({\lambda}_{i})}_{i=1}^{k}\) , on déflate aussi en tenant compte de leur conjugué [18]_ pour ne pas les recalculer

\(p(\lambda ):=\frac{\text{dét}\mathrm{Q}(\lambda )}{\prod_{i=1}^{k}(\lambda -{\lambda}_{i})(\lambda -\stackrel{ˉ}{{\lambda}_{i}})}\) (3.2-4)

L’évaluation du dénominateur est triviale. Celle du numérateur, \(\text{dét}\mathrm{Q}(\lambda )\) , s’effectue en factorisant [19]_ sous forme \({\mathit{LDL}}^{T}\) la matrice dynamique \(\mathrm{Q}(\lambda )\) et en effectuant le produit des termes diagonaux de la matrice diagonale \(\mathrm{D}\) .

Cette factorisation de matrice dynamique est très coûteuse avec une complexité calcul en \(O(\mathrm{bn²})\) où \(b\) est la largeur de bande et \(n\) la taille du problème. C’est la partie la gourmande en capacité de traitement et en place mémoire de l’algorithme. Elle est à reproduire \(p(i+2)\) fois avec \(p\) le nombre de modes souhaités et \(i\) le nombre d’itérations moyen à convergence. L’algorithme est donc plutôt à réserver aux QEP de petite taille (<104degrés de liberté) ou aux benchmarks de méthodes.

Remarques:

• Müller[Mul56] propose une procédure particulière et efficace pour initialiser le processus dans le cas où la fonction à interpoler \(p(\lambda )\) est un polynôme d’ordre \(n\) dont les coefficients sont connus. Le cas de figure qui nous intéresse ne rentre malheureusement pas dans cette classe de problème. L’initialisation de l’algorithme requiert ainsi 2 évaluations préliminaires de \(\mathrm{p.}\)

• Cette méthode ne calcule pas de position modale, elle gère donc mal les multiplicités. Contrairement à la dichotomie mise en œuvre pour les GEP. La notion de seuil en deçà duquel deux valeurs numériques sont supposées multiples est lié au paramètre \({\varepsilon}_{\text{PREC\_AJUSTE}}\)

• On ne retient dans la formule (3.2-2) que le complexe à partie imaginaire positive.

Méthode des puissances inverses (variante de Jennings)

La seconde partie consiste à affiner les estimations des valeurs propres calculées par la méthode de Müller-Traub et à calculer leurs vecteurs propres associés. Pour ce faire, on linéarise le QEP suivant la stratégie

(L1)’(cf. §2.4) et on applique une transformation spectrale particulière due à Jennings [Jen77](1977) . Son objectif est de ne faire apparaître, comme inverse de matrice, que celui de la matrice dynamique \(\mathrm{Q}{(\lambda )}^{-1}\) . L’algorithme mis en place dans le code se découpe comme suit:

• Initialisation de la valeur propre à partir de l’estimation de la première phase: \({\lambda}_{0}\) .

• Construction de la matrice dynamique \(\mathrm{Q}({\lambda}_{0})\) et de sa factorisation.

• Initialisation du processus itératif par les vecteurs: \({\mathrm{x}}_{0}\text{aléatoire,}{\mathrm{y}}_{0}={\lambda}_{0}{\mathrm{x}}_{0}.\)

• Pour \(k=0,\text{NMAX\_ITER-1}\) faire

• Normaliser \({\mathrm{x}}_{k}\leftarrow \frac{{\mathrm{x}}_{k}}{{({\stackrel{ˉ}{\mathrm{x}}}_{k})}^{t}{\mathrm{Mx}}_{k}},{\mathrm{y}}_{k}\leftarrow \frac{{\mathrm{y}}_{k}}{{({\stackrel{ˉ}{\mathrm{y}}}_{k})}^{t}{\mathrm{My}}_{k}}\)

• Résoudre \(\mathrm{Q}({\lambda}_{0}){\mathrm{x}}_{k+1}={\mathrm{Cx}}_{k}+{\lambda}_{0}{\mathrm{Mx}}_{k-1}+{\mathrm{My}}_{k-1}\)

• Calculer \({\mathrm{y}}_{k+1}=-{\lambda}_{0}{\mathrm{x}}_{k+1}+{\mathrm{x}}_{k}\)

• Evaluer \({\lambda}_{k}=\frac{{(\stackrel{ˉ}{{\mathrm{x}}_{k}})}^{T}{\mathrm{Cx}}_{k}+\sqrt{{({(\stackrel{ˉ}{{\mathrm{x}}_{k}})}^{T}{\mathrm{Cx}}_{k})}^{2}-4({(\stackrel{ˉ}{{\mathrm{x}}_{k}})}^{T}{\mathrm{Mx}}_{k})({(\stackrel{ˉ}{{\mathrm{x}}_{k}})}^{T}{\mathrm{Kx}}_{k})}}{2({(\stackrel{ˉ}{{\mathrm{x}}_{k}})}^{T}{\mathrm{Mx}}_{k})}\) (3.3-1)

• Si \(\frac{∣{\lambda}_{k+1}-{\lambda}_{k}∣}{∣{\lambda}_{k}∣}<\text{PREC}\) et \(\max(\mid \text{Re}({\lambda}_{k}-{\lambda}_{0})\mid ,\mid \text{Im}({\lambda}_{k}-{\lambda}_{0})\mid )<\text{PREC}\) alors (3.3-2)

• Le mode propre solution est \(({\lambda}_{k},{\mathrm{x}}_{k})\) , exit.

Algorithme 2. Méthode des puissances inverses (variante de Jennings).

Le critère d’arrêt PRECet le nombre maximal d’itérations autorisées NMAX_ITER sont des arguments du mot-clé facteur CALC_MODE.

Remarques:

• Le critère d’arrêt est basé sur l’extension du quotient de Rayleigh (cf. [Solveurs modaux et résolution du problème généralisé] § 4.3) au QEP. Connaissant un vecteur propre \({\mathrm{x}}_{k}\) une bonne estimation de sa valeur propre \({\lambda}_{k}\) est fournie par la formule (3.3-1). On stoppe l’algorithme dès que les variations relatives (en complexe) et absolues (sur les parties réelle et imaginaire séparément), par rapport à l’estimée initiale, sont inférieures à une certaine valeur. Ainsi conjointement, on affine l’estimée initiale de la valeur propre et on détermine son vecteur propre associé.

• Malgré la linéarisation du problème initial, la formulation «astucieuse» de Jennings ne manipule que des entités de taille \(n\) . En réalité, l’algorithme 2 peut se schématiser sous la forme vectorielle

• \(\left[\begin{array}{cc}\mathrm{Q}({\lambda}_{0})& {0}_{n}\\ -{\lambda}_{0}{\mathrm{I}}_{n}& {\mathrm{I}}_{n}\end{array}\right]\left\lbrace \begin{array}{c}{\mathrm{x}}_{k+1}\\ {\mathrm{y}}_{k+1}\end{array}\right\rbrace =\left[\begin{array}{cc}-C-{\lambda}_{0}\mathrm{M}& -\mathrm{M}\\ -{\mathrm{I}}_{n}& {0}_{n}\end{array}\right]\left\lbrace \begin{array}{c}{\mathrm{x}}_{k}\\ {\mathrm{y}}_{k}\end{array}\right\rbrace\) (3.3-3)

Périmètre d’utilisation

QEP à matrices réelles symétriques.

L’utilisateur ne peut spécifier que TYPE_RESU=’DYNAMIQUE’ comme classe d’appartenance de son calcul (pas de flambement). On renseigne alors éventuellement le vecteur FREQ (et pas CHAR_CRIT).

Affichage dans le fichier message

Dans le fichier message, les résultats sont affichés sous forme de tableau

---

LE NOMBRE DE DDL

TOTAL EST: 56

DE LAGRANGE EST: 32

LE NOMBRE DE DDL ACTIFS EST: 8

---

CALCUL MODAL: METHODE D’ITERATION INVERSE

MULLER INVERSE

NUMERO FREQUENCE(HZ) AMORTISSEMENT NB_ITER PRECISION NB_ITER PRECISION NORME D’ERREUR

1 5.52915E+00 1.52090E-02 5 9.49522E-07 1 4.47353E-17 2.17748E-15

2 1.08959E+01 2.87575E-02 3 1.17993E-05 1 3.24203E-17 4.83597E-15

3 1.59270E+01 3.95645E-02 5 2.53170E-05 1 8.86844E-18 7.12414E-16

…

---

VERIFICATION A POSTERIORI DES MODES

---

Exemple 3. Trace de CALC_MODES(QEP) dans le fichier message.Extrait du cas-test sdld27b.

Pour chaque valeur propre (représentée sous la forme FREQUENCE= \(\frac{\text{Im}({\lambda}_{i})}{2\pi }\) et AMORTISSEMENT= \(\frac{-\text{Re}({\lambda}_{i})}{\mid {\lambda}_{i}\mid }\) cf §2.5), on trace le nombre d’itérations et la précision obtenue [20]_ des deux phases de l’opérateur modal. Avec l’option “PROCHE”, les colonnes concernant la méthode de Müller n’apparaissent évidemment pas. La dernière colonne, NORME D’ERREUR, reprend la norme d’erreur du résidu déterminé suivant l’algorithme n°1 (§2.5).

Récapitulatif du paramétrage

Récapitulons maintenant le paramétrage de l’opérateur CALC_MODES avec OPTION parmi [“PROCHE”,”AJUSTE”].

Opérande	Mot-clé	Valeur par défaut	Références
	TYPE_RESU=”DYNAMIQUE”	“DYNAMIQUE”	§2.1
	OPTION =‘AJUSTE’	“AJUSTE”	§3.1
	‘PROCHE’		§3.1
CALC_FREQ	FREQ (siPROCHE)		§3.1
	NMAX_FREQ	0	§3.1
SOLVEUR_MODAL	NMAX_ITER_AJUSTE	15	§3.2
	PREC_AJUSTE	1.E-04	§3.2
	OPTION_INV=‘DIRECT’	“DIRECT”	§3.3
	‘RAYLEIGH’(même traitement)		§3.3
	PREC_INV	1.E-05	§3.3
	NMAX_ITER_INV	30	§3.3
VERI_MODE	STOP_ERREUR=”OUI”	“OUI”	[Procédure de dénombrement de valeurs propres] §3.2
	“NON”
	SEUIL	1.E-02	[Procédure de dénombrement de valeurs propres] §3.2

Tableau 3.6-1. Récapitulatif du paramétrage de CALC_MODES avec OPTIONparmi [“PROCHE”,”AJUSTE”] (QEP).

Remarques:

• On retrouve toute la « tripaille » de paramètres liée aux post-traitements de vérification ( SEUIL_FREQ, VERI_MODE ).

• Lors des premiers passages, il est fortement conseillé de ne modifier que les paramètres principaux notés en gras. Les autres concernent plus les arcanes de l’algorithme et ils ont été initialisés empiriquement à des valeurs standards.

[1]

C’est aussi pour cette raison que le périmètre d’utilisation de CALC_MODES avec OPTION parmi [“PROCHE”,”AJUSTE”], est limité au QEP standard (symétrique réel). On s’attend à trouver uniquement des valeurs propres conjuguées par paire.

[2]

On s’attend à pouvoir factoriser sans encombre (sans avoir recours au pivotage) et à manipuler une factorisée symétrique. C’est une autre raison qui limite le périmètre de CALC_MODES avec OPTION parmi [“PROCHE”,”AJUSTE”], au QEP symétrique réel.

[3]

Précision de la méthode de Müller au sens de la formule (3.2-3), celle de la méthode de Jennings reprend la formule (3.3-2).

Méthode de sous-espace (METHODE=’TRI_DIAG’/’SORENSEN’)

Du QEP au SEP

Passage QEP/GEP: linéarisation

Comme on l’a expliqué au §2.4, le passage du QEP au GEP s’effectue via la linéarisation suivante

(L2)’’ \((\underset{\mathrm{A}}{\underset{\underbrace{}}{\left[\begin{array}{cc}-\mathrm{K}& {0}_{n}\\ {0}_{n}& {\mathrm{M}}_{R}\end{array}\right]}}-\lambda \underset{\mathrm{B}}{\underset{\underbrace{}}{\left[\begin{array}{cc}-\mathrm{C}& \mathrm{M}\\ \mathrm{M}& {0}_{n}\end{array}\right]}})\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]={0}_{\mathrm{2n}}\) (4.1-1)

Le GEP obtenu est structurellement symétrique réel et hérite des propriétés des matrices \(\mathrm{K}\) , \(\mathrm{M}\) et \(\mathrm{C}\) . Dès que l’une d’elles est non symétrique ou complexe (possible uniquement avec \(\mathrm{K}\) ), le GEP devient non symétrique ou complexe. L’aspect symétrique ou non du GEP n’est en fait pas bloquant car chaque méthode va s’employer à trouver le couple «opérateur de travail/produit scalaire» qui lui permet de fonctionner . D’ailleurs, seul Lanczos a besoin d’un couple symétrique. La méthode IRAM peut-elle très bien fonctionner en non symétrique et avec souvent plus de robustesse [21]_ .

Passage GEP/SEP: transformation spectrale

Puis une transformation spectrale (cf. [Solveurs modaux et résolution du problème généralisé] §3.7/§5) permet d’achever la phase de pré-traitement du QEP en transformant le GEP \(({L}_{2})’’\) précédente en un SEP. Lorsqu’on manipule des modes complexes, deux cas de figures sont envisageables:

• Travailler complètement en arithmétique complexe ,

• Rester le plus souvent en arithmétique réelle en isolant les contributions réelle et imaginaire de la transformation spectrale classique.

La première stratégie utilise classiquement la transformation spectrale dite de «shift-and-invert» et métamorphose \(({L}_{2})’’\) en

\(({S}_{1})\) \(\begin{array}{c}\underset{{\mathrm{A}}_{\sigma}}{\underset{\underbrace{}}{{(\mathrm{A}-\sigma \mathrm{B})}^{-1}\mathrm{B}}}\mathrm{w}=\mu \mathrm{w}\\ \text{avec}\mu :=\frac{1}{\lambda -\sigma },\mathrm{w}:=\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\end{array}\) (4.1-2)

Elle est mise en œuvre en activant l’approche complexe d’IRAM (‘SORENSEN’+APPROCHE=’COMPLEXE’). C’est la seule alternative lorsqu’on manipule une matrice \(\mathrm{K}\) complexe [22]_ . Elle est utilisable avec \(\mathrm{K},\mathrm{M}\) et \(\mathrm{C}\) symétriques ou non.

La seconde stratégie repose sur la transformation spectrale à «double shifts somme» introduite en 1987 par B.N.Parlett & Y.Saad[PS87], qui conduit à manipuler deux types de SEP

\(({S}_{2})\) \(\begin{array}{c}\text{Re}\underset{{\mathrm{A}}_{\sigma}^{+}}{\underset{\underbrace{}}{({(\mathrm{A}-\sigma \mathrm{B})}^{-1}\mathrm{B})}}\mathrm{w}={\mu}^{+}\mathrm{w}\\ \text{avec}{\mu}^{+}:=\frac{1}{2}(\frac{1}{\lambda -\sigma }+\frac{1}{\lambda -\stackrel{ˉ}{\sigma}}),\mathrm{w}:=\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\end{array}\) (4.1-3)

\(({S}_{3})\) \(\begin{array}{c}\text{Im}\underset{{\mathrm{A}}_{\sigma}^{-}}{\underset{\underbrace{}}{({(\mathrm{A}-\sigma \mathrm{B})}^{-1}\mathrm{B})}}\mathrm{w}={\mu}^{-}\mathrm{w}\\ \text{avec}{\mu}^{-}:=\frac{1}{\mathrm{2i}}(\frac{1}{\lambda -\sigma }+\frac{1}{\lambda -\stackrel{ˉ}{\sigma}}),\mathrm{w}:=\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\end{array}\) (4.1-4)

Ces deux variantes sont mises en oeuvre avec Lanczos et IRAM via les valeurs, respectivement, ‘REEL’ et ‘IMAG’ du mot-clé APPROCHE. Avec Lanczos, le périmètre d’utilisation est limité aux QEP classiques (symétriques réels), tandis qu’IRAM est moins restreint et permet aussi le non symétrique.

Remarques:

• L’origine de la dénomination «double shifts somme» de la transformation spectrale ci-dessus apparaît plus clairement si on reformule les opérateurs comme suit

\(\begin{array}{c}{\mathrm{A}}_{\sigma}^{+}:=\frac{1}{2}\left[{(\mathrm{A}-\sigma \mathrm{B})}^{-1}+{(\mathrm{A}-\stackrel{ˉ}{\sigma}\mathrm{B})}^{-1}\right]\mathrm{B}\\ {\mathrm{A}}_{\sigma}^{-}:=\frac{1}{\mathrm{2i}}\left[{(\mathrm{A}-\sigma \mathrm{B})}^{-1}+{(\mathrm{A}-\stackrel{ˉ}{\sigma}\mathrm{B})}^{-1}\right]\mathrm{B}\end{array}\) (4.1-5)

• La terminologie «somme» est par opposition au «double shifts produit» proposé initialement par J.C.Francis et utilisé d’ailleurs, en sous-main, par la méthode QZ et les redémarrages «implicites» d’IRAM (cf. [Solveurs modaux et résolution du problème généralisé] Annexe 1):

\({\mathrm{A}}_{\sigma}:={\left[(\mathrm{A}-\sigma \mathrm{B})+(\mathrm{A}-\stackrel{ˉ}{\sigma}\mathrm{B})\right]}^{-1}\mathrm{B}\) (4.1-6)

Celle-ci a l’inconvénient majeur de produire des matrices de travail très remplies (donc plus coûteuses à manipuler et à stocker) par rapport à la stratégie de Parlett & Saad.

Une fois déterminé l’opérateur de travail, il reste à choisir son produit scalaire associé . Pour IRAM, on peut se contenter d’un couple «opérateur/produit scalaire» non symétrique (cf. [Solveurs modaux et résolution du problème généralisé] §6), mais pour Lanczos, celui-ci doit obligatoirement être symétrique (cf. [Solveurs modaux et résolution du problème généralisé] §7). Plusieurs choix de (pseudo)produit-scalaires matriciels permettent cette symétrie:

\(\begin{array}{c}{(\mathrm{w},\mathrm{z})}_{\mathrm{B}}:={\mathrm{z}}^{T}\mathrm{Bw}\\ {(\mathrm{w},\mathrm{z})}_{+}:={\mathrm{z}}^{T}\text{Re}{\left[{(\mathrm{A}-\sigma \mathrm{B})}^{-1}\right]}^{-1}\mathrm{w}\\ {(\mathrm{w},\mathrm{z})}_{-}:={\mathrm{z}}^{T}\text{Im}{\left[{(\mathrm{A}-\sigma \mathrm{B})}^{-1}\right]}^{-1}\mathrm{w}\end{array}\) (4.1-7)

Le premier fonctionne avec les deux formulations (S2) et (S3), tandis que le deuxième (resp. le troisième) ne rend symétrique que \({\mathrm{A}}_{\sigma}^{+}\) (resp. \({\mathrm{A}}_{\sigma}^{-}\) ). Dans l’implantation effective dans Code _Aste r, ce sont ces deux derniers qui sont retenus pour les approches réelle et imaginaire de Lanczos.

Remarque:

• Ces produits scalaires sont les extensions «naturelles» du (pseudo) produit scalaire matriciel utilisé par la variante de Newmann-Pipano en GEP:

\((\mathrm{w},\mathrm{z}):={\mathrm{z}}^{T}(\mathrm{A}-\sigma \mathrm{B})\mathrm{w}\) (4.1-8)

Implantation dans Code_Aster

Choix du décalage spectral

En GEP il y’a quatre façons de choisir ce shift (cf. [Solveurs modaux et résolution du problème généralisé]§5.4). En QEP, l’option ‘BANDE’ étant proscrite du fait de l’absence de test de Sturm adapté, ce chiffre est réduit à trois:

• \(\sigma =0,\) on cherche les plus petites valeurs propres du problème de départ. Ceci correspond à OPTION=”PLUS_PETITE’sous le mot-clé facteur CALC_FREQ.

• \(\sigma ={\sigma}_{0}\text{avec}{\sigma}_{0}={(2\pi {f}_{0})}^{2},\) on cherche les fréquences proches de la fréquence FREQ= \({f}_{0}\) (OPTION=”CENTRE”).

• \(\sigma ={\sigma}_{0}\text{avec}{\sigma}_{0}={(2\pi {f}_{0})}^{2}\text{et}{\sigma}_{1}=\eta {(2\pi {f}_{0})}^{2}/2\) on cherche les fréquences proches de la fréquence FREQ= \({f}_{0}\) (OPTION=”CENTRE”) et de l’amortissement réduit AMOR_REDUIT= \(\eta\) (OPTION=”CENTRE”).

Le nombre de fréquences à calculer est donné en général par l’utilisateur à l’aide de NMAX_FREQ sous le mot-clé facteur CALC_FREQ.

Calcul de l’opérateur et du produit scalaire de travail

On récapitule dans le tableau suivant l’ensemble des couples (opérateur de travail, produit scalaire) possibles, en fonction des options choisies par l’utilisateur.

APPROCHE/METHODE	‘TRI_DIAG’	‘SORENSEN’ (par défaut)
‘COMPLEXE’	Indisponible	\(({A}_{\sigma},{L}^{2})\)
‘REEL’ ( par défaut)	\(({A}_{\sigma}^{+},{()}_{+})\)	\(({A}_{\sigma}^{+},{L}^{2})\)
‘IMAG’	\(({A}_{\sigma}^{-},{()}_{-})\)	\(({A}_{\sigma}^{-},{L}^{2})\)

Tableau 4.2-1. Ensemble des couples (opérateur de travail, produit scalaire) possibles en fonction des options choisies par l’utilisateur.

Lecalcul des opérateurs de travail \({\mathrm{A}}_{\sigma}^{\pm}\) peut se ramener à de simples produits matrice-vecteur et à une inversion de la matrice dynamique \(\mathrm{Q}(\lambda )\) (préalablement factorisée) via la formulation (astucieuse) suivante

\({(\mathrm{A}-\sigma \mathrm{B})}^{-1}\mathrm{B}=\left[\begin{array}{cc}{0}_{n}& {0}_{n}\\ {\mathrm{M}}_{R}^{-1}\mathrm{M}& {0}_{n}\end{array}\right]-\left[\begin{array}{cc}\mathrm{Q}{(\sigma )}^{-1}& {0}_{n}\\ {0}_{n}& \sigma {\mathrm{M}}_{R}^{-1}\mathrm{M}\mathrm{Q}{(\sigma )}^{-1}\end{array}\right]\left[\begin{array}{cc}\mathrm{C}+\sigma \mathrm{M}& \mathrm{M}\\ \mathrm{C}+\sigma \mathrm{M}& \mathrm{M}\end{array}\right]\) (4.2-1)

D’autre part, cette formulation a le «bon goût» de permettre des traitements en pure arithmétique réelle et de conserver la structure creuse des matrices. On peut ainsi construire le mode opératoire suivant:

• Préparation dans \(C\) [23]_

(une seule fois, en phase d’initialisation de l’algorithme, cf.

§2.5):

Former \(\mathrm{Q}(\sigma ):=({\sigma}^{2}\mathrm{M}+\sigma \mathrm{C}+\mathrm{K})\)

Factoriser \(\mathrm{Q}(\sigma )\) sous forme LDLT.

• Calcul de \(\begin{array}{c}\\ {\mathrm{A}}_{{\sigma}^{\pm}}\underset{\mathrm{w}}{\underset{\underbrace{}}{\left[\begin{array}{c}\mathrm{u}\\ \mathrm{v}\end{array}\right]}}\end{array}\) (plusieurs fois par itération de Lanczos ou d’IRAM):

Former dans \(R{\mathrm{u}}_{1}:=\mathrm{Cu},{\mathrm{u}}_{2}:=\mathrm{Mu},{\mathrm{u}}_{3}:=\mathrm{Mv},\)

Calculer dans \(C{\mathrm{u}}_{4}:=\mathrm{Q}{(\sigma )}^{-1}+\sigma {\mathrm{u}}_{2}+{\mathrm{u}}_{3}\)

• Suivant le choix de l’approche \(\begin{array}{c}{\mathrm{A}}_{\sigma}^{+}\mathrm{w}=\left[\begin{array}{c}-\text{Re}({\mathrm{u}}_{4})\\ {\mathrm{M}}_{R}^{-1}\mathrm{M}(\mathrm{u}-\text{Re}(\sigma {\mathrm{u}}_{4}))\end{array}\right]\\ {\mathrm{A}}_{\sigma}^{-}\mathrm{w}=\left[\begin{array}{c}-\text{Im}({\mathrm{u}}_{4})\\ -{\mathrm{M}}_{R}^{-1}\mathrm{M}(\mathrm{u}-\text{Im}(\sigma {\mathrm{u}}_{4}))\end{array}\right]\end{array}\)

Algorithme 3. Mode opératoire pour calculer les opérateurs de travail \({\mathrm{A}}_{\sigma}^{\pm}\) .

En approche complexe, le calcul se fait directement via la formule (4.1-2) en arithmétique complexe.

Concernant les produits scalaires matriciels, leurs manipulations dans Code_Aster s’appuient sur les formulations suivantes:

\(\begin{array}{c}{(\mathrm{w},\mathrm{z})}_{+}:={\mathrm{z}}^{T}\left[(\mathrm{A}-\text{Re}(\sigma )\mathrm{B})+{\text{Im}}^{2}(\sigma )\mathrm{B}{(\mathrm{A}-\text{Re}(\sigma )\mathrm{B})}^{-1}\mathrm{B}\right]\mathrm{w}\\ {(\mathrm{w},\mathrm{z})}_{-}:={\mathrm{z}}^{T}\left[\frac{1}{\text{Im}(\sigma )}(\mathrm{A}-\text{Re}(\sigma )\mathrm{B}){\mathrm{B}}^{-1}(\mathrm{A}-\text{Re}(\sigma )\mathrm{B})+{\text{Im}}^{2}(\sigma )\mathrm{B}\right]\mathrm{w}\end{array}\) (4.2-2)

Remarque:

• Le calcul des produits scalaires matriciels peut s’effectuer uniquement en arithmétique réelle. Pour le premier produit scalaire on a besoin d’un matrice dynamique supplémentaire \(\mathrm{Q}(\text{Re}(\sigma ))\) et de sa factorisée (d’où un surcoût calcul et mémoire).

Périmètre d’utilisation

QEP à matrices symétriques réelles pour Lanczos, éventuellement non symétriques (et \(\mathrm{K}\) peut être symétrique complexe) pour IRAM.

L’utilisateur ne peut spécifier que TYPE_RESU=’DYNAMIQUE’ comme classe d’appartenance de son calcul (pas de flambement). On renseigne alors éventuellement le vecteur FREQ (et pas CHAR_CRIT).

Affichage dans le fichier message

Dans le fichier message, les résultats sont affichés sous forme de tableau

---

LE NOMBRE DE DDL

TOTAL EST: 130

DE LAGRANGE EST: 16

LE NOMBRE DE DDL ACTIFS EST: 106

---

CALCUL MODAL: METHODE D’ITERATION SIMULTANEE

METHODE DE SORENSEN

NUMERO FREQUENCE (HZ) AMORTISSEMENT NORME D’ERREUR

1 1.24163E+02 -1.89229E-02 1.16550E-09

2 1.24164E+02 1.89229E-02 1.04882E-09

…

6 1.07321E+03 1.63495E-01 1.36092E-13

NORME D’ERREUR MOYENNE: 0.36947E-09

Exemple 4. Trace de CALC_MODES avec OPTIONparmi [“CENTRE”,”PLUS_PETITE”] et METHODE=‘SORENSEN’(QEP)dans le fichier message. Extrait du cas-test sdll123b.

Pour chaque valeur propre (représentée sous la forme FREQUENCE= \(\frac{\text{Im}({\lambda}_{i})}{2\pi }\) et AMORTISSEMENT= \(\frac{-\text{Re}({\lambda}_{i})}{\mid {\lambda}_{i}\mid }\) cf §2.5), on trace la norme d’erreur du résidu déterminé suivant l’algorithme n°1 (§2.5).

Récapitulatif du paramétrage

Récapitulons maintenant le paramétrage de l’opérateur CALC_MODES avec OPTION parmi [“CENTRE”,”PLUS_PETITE”] pour les méthodes de Lanczos et de Sorensen.

Opérande	Mot-clé	Valeur par défaut	Références
	TYPE_RESU =”DYNAMIQUE”	“DYNAMIQUE”	§2.1
	OPTION =‘PLUS_PETITE’	‘PLUS_PETITE’	§4.2
	‘CENTRE’		§4.2
CALC_FREQ	FREQ (si’CENTRE”)		§4.2
	AMOR_REDUIT (si’CENTRE”)		§4.2
	NMAX_FREQ	10	§4.2
SOLVEUR_MODAL	METHODE =‘TRI_DIAG’	‘SORENSEN’	[Solveurs modaux et résolution du problème généralisé] §6
	‘SORENSEN’		[Solveurs modaux et résolution du problème généralisé] §7
	APPROCHE = ’REEL’	’REEL’	§4.1
	‘IMAG’		§4.1
	‘ COMPLEXE’ (si “SORENSEN”)		§4.1
	DIM_SOUS_ESPACE COEF_DIM_ESPACE	Calculé	[Solveurs modaux et résolution du problème généralisé] §6/7
	Paramétrage usuel de ‘TRI_DIAG’/’SORENSEN’		[Solveurs modaux et résolution du problème généralisé] §6/7

Tableau 4.5-1. Récapitulatif du paramétrage de CALC_MODES avec OPTIONparmi [“CENTRE”,”PLUS_PETITE”] et METHODE=‘TRI_DIAG’ou ’SORENSEN’ (QEP) .

Remarques:

• Lors des premiers passages, il est fortement conseillé de ne modifier que les paramètres principaux notés en gras. Les autres concernent plus les arcanes de l’algorithme et ils ont été initialisés empiriquement à des valeurs standards.

• En particulier, pour améliorer la qualité d’un mode, le paramètre fondamental est la dimension du sous-espace de projection du SEP, DIM_SOUS_ESPACE .

[1]

Au prix toutefois d’un surplus de calculs et de stockages.

[2]

Dès que :math:`mathrm{K}`est complexe symétrique et que l’on a choisi IRAM, quelque soit l’approche choisie dans le fichier de commande, c’est l’approche ‘COMPLEXE’ qui est prise par défaut.

[3]

Car le shift (ainsi que:math:`mathrm{K}`avec IRAM) peut être complexe.

Méthode globale QZ (METHODE=’QZ’)

Du QEP au GEP

Contrairement aux méthodes vues précédemment, la méthode QZ est capable de gérer directement un GEP. On a donc qu’une seule transformation à appliquer au QEP initiale, c’est la linéarisation . Celle empiriquement retenue et détaillée au §2.4 est

(L2)* \((\underset{\mathrm{A}}{\underset{\underbrace{}}{\left[\begin{array}{cc}-\mathrm{K}& {0}_{n}\\ {0}_{n}& \alpha {\mathrm{I}}_{n}\end{array}\right]}}-\lambda \underset{\mathrm{B}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{C}& \mathrm{M}\\ \alpha {\mathrm{I}}_{n}& {0}_{n}\end{array}\right]}})\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]={0}_{\mathrm{2n}}\) (5.1-1)

Le GEP obtenu n’est pas structurellement symétrique (ou hermitien) mais ce n’est pas préjudiciable au calcul car les drivers de LAPACK adaptés ne travaillent qu’en non symétrique. Pour plus de robustesse et pour simplifier la structuration du code , on formule et résout le GEP (L2)* dans le plan complexe . Comme pour les GEP à modes complexes on utilise les routines LAPACK standard (ZGGEV) ou expert (ZGGEVX) suivant le paramètre TYPE_QZ=’QZ_SIMPLE’/’QZ_EQUI’.

A l’issu du calcul modal, QZ renvoie à Code_Aster les modes calculés et ces derniers vont être vérifiés, filtrés et ordonnés suivant les mêmes préconisations que pour les modes complexes des GEP (cf. [Solveurs modaux et résolution du problème généralisé] §9.4).

Périmètre d’utilisation

QEP à matrices symétriques ou non (et \(\mathrm{K}\) peut être symétrique complexe).

L’utilisateur ne peut spécifier que TYPE_RESU=’DYNAMIQUE’ comme classe d’appartenance de son calcul (pas de flambement). On renseigne alors éventuellement le vecteur FREQ (et pas CHAR_CRIT).

Impressions dans le fichier message

Un extrait commenté du cas-test sdll113a est repris au §2.5.

Récapitulatif du paramétrage

Récapitulons maintenant le paramétrage de l’opérateur CALC_MODES avec OPTION parmi [“CENTRE”,”PLUS_PETITE”,”TOUT”] pour la méthode QZ.

Opérande	Mot-clé	Valeur par défaut	Références
	TYPE_RESU=”DYNAMIQUE”	“DYNAMIQUE”	§2.1
	OPTION =‘PLUS_PETITE’	‘PLUS_PETITE’	[Solveurs modaux et résolution du problème généralisé] §9.4
	“CENTRE’ “TOUT”		[Solveurs modaux et résolution du problème généralisé] §9.4
CALC_FREQ	FREQ (siCENTRE)		[Solveurs modaux et résolution du problème généralisé] §9.4
	AMOR_REDUIT (siCENTRE)		[Solveurs modaux et résolution du problème généralisé] §9.4
	NMAX_FREQ	10	[Solveurs modaux et résolution du problème généralisé] §9.4
	APPROCHE=…	Sans objet
SOLVEUR_MODAL	METHODE =‘QZ’	‘SORENSEN’	[Solveurs modaux et résolution du problème généralisé] §9
	Paramétrage usuel de ‘QZ’		[Solveurs modaux et résolution du problème généralisé] §9

Tableau 5.4-1. Récapitulatif du paramétrage de CALC_MODES avec OPTION parmi [“”CENTRE”,”PLUS_PETITE”,”TOUT”] avec METHODE=‘QZ’ (QEP).

Remarque:

• Lors des premiers passages, il est fortement conseillé de ne modifier que les paramètres principaux notés en gras. Les autres concernent plus les détails de l’algorithme et ils ont été initialisés empiriquement à des valeurs standards.

Bibliographie

[DB08] G.Dahlquist & A.Bjorck. Numerical methods in scientific computing . SIAM (2008).

[GLR82] I.Gohberg, P.Lancaster et L.Rodman. Matrix polynomials . Ed. Academic Press (1982).

[Jen77] A.Jennings. Matrix computation for engineers and scientists . Ed. Wisley (1977).

[LT85] P.Lancaster et M.Tsimenetsky. The theory of matrices . Ed. Academic Press (1985).

[Mac06] D.S.Mackey. Structured linearizations for matrix polynomials . Thèse de l’Université de Manchester (2006).

[Mee07] K.Meerbergen. The quadratic Arnoldi method for the solution of the QEP . Rapport interne de l’Université de Leuven (2007).

[Mul56] D.E.Müller. A method for solving algebraic equations using an automatic computer . Math. Tab. Wash., 10, pp208-215 (1956).

[PS87] B.N.Parlett & Y.Saad. Complex shift and invert strategies for real matrices . Linear algebra and its applications, 88-89, pp575-595 (1987).

[TM01] F.Tisseur & K.Meerbergen. The quadratics eigenvalue problem . SIAM review, vol. 43, n°2, pp235-286 (2001).