# Continuation interactive de diagrammes de bifurcation complexes

Bruno Cochelin

<figure><img src="image_1.png" /></figure>

## To cite this version:

Bruno Cochelin. Continuation interactive de diagrammes de bifurcation complexes. 17e Congrès Français de Mécanique, Sep 2005, Troyes, France. pp.1-9. hal-04089459

HAL Id: hal-04089459

https://hal.science/hal-04089459

Submitted on 4 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Continuation interactive de diagrammes de bifurcation complexes

**Bruno Cochelin**

Laboratoire de Mécanique et d'Acoustique (LMA)
Ecole Généraliste d'Ingénieur de Marseille (EGIM)
Technopôle de Château-Gombert 13381 Marseille cedex 13
bruno.cochelin@egim-mrs.fr

**Résumé :**

L'objet de cette communication est de présenter une méthode et un logiciel de continuation interactif destinés à tracer des diagrammes de bifurcation complexes. La méthode de continuation repose sur des développements en séries des branches de solutions : la méthode asymptotique numerique. Cette technique est très bien adaptée à la détection des bifurcations et au branchement par perturbation. Une interface graphique permet de gérer interactivement le pilotage de la continuation vers telle ou telle branche, tout en offrant des possibilités de visualisation et d'analyse définies par l'utilisateur. Lors de la présentation, on s'appuiera sur des exemples concrets, généralement issus de la mécanique des structures, pour dégager l'intérêt et les limitations de ces techniques par rapport aux méthodes et aux logiciels de prédictions-corrections.

**Abstract :**

In this communication, we describe a method and an interactif software for the path-following of complex bifurcation diagrams. The continuation method is the so-called Méthode Asymptotique Numerique which relies on high order power series expansions of the branches of solutions. A graphical interface has been developed in order to analyse interactively the solutions and to take the suitable decisions for the path following. During the oral presentation, many examples from structural mechanics will be considered to show the advantages and drawbacks as compared to classical predictor corrector algorithms.

**Mots-clefs :**

Continuation ; bifurcation ; Branchement ; Perturbations ; Méthode Asymptotique Numérique ;

## 1 Méthodes de continuation

Les méthodes de continuation par rapport à un paramètre sont très couramment utilisées en science de l'ingénieur pour comprendre et simuler le comportement des systèmes physiques non linéaires qui sont modélisés par des équations algébriques, des équations différentielles ou aux dérivées partielles. Elles ont pour objet de trouver les "branches" de solutions de l'équation $R(U) = 0$, où $R : \mathbb{R}^{N+1} \rightarrow \mathbb{R}^N$, est une fonction régulière de $U$. Pour les équations différentielles ou aux dérivées partielles, on se ramène à cette forme en appliquant au préalable une méthode de discrétisation de type éléments finis ou différences finies. En mécanique des structures, l'équation $R(U) = 0$ représente en général des équations d'équilibre aux noeuds d'un maillage et des relations de comportement aux points de Gauss des éléments. Le vecteur inconnu $U$ contient typiquement des champs de déplacements, de déformations ou de contraintes. Il comporte également un paramètre de continuation tel que : l'intensité du chargement, la fréquence d'une oscillation, une élévation de température, un coefficient matériau, etc.

Parmi les diverses méthodes de continuation qui ont été développées pour résoudre l'équation $R(U) = 0$, la plus populaire est certainement la méthode de prédiction-corrections [1] [2] [3],
[4] que l'on dénomme également par Newton-Raphson ou Incrémentale-itératives en mécanique des structures [5] [6] [7]. Cette technique consiste à générer une suite de points solutions pour décrire les branches de solutions. De nombreuses variantes existent, selon le choix du prédicteur, du correcteur, et de la stratégie d'adaptation des pas, c'est à dire de la gestion de la distance entre les points solutions. Des logiciels de continuation généraux tels que AUTO2000 (E. Doedel) [8] où BIFPACK (R. Seydel) [9] sont disponibles sur internet. Les mécaniciens des structures ont aménagé ces algorithmes généraux pour leurs modèles éléments finis qui sont généralement de grande taille, et qui peuvent comporter des non linéarités matérielles (hyper-élasticité, plasticité, endommagement), des conditions de contact-frottement, et des non linéairités géométriques (post-flambement), etc. Ce sont aujourd'hui ces techniques de prédiction-corrections qui équipent les logiciels de calcul de structures standards du commerce.

Quelques méthodes alternatives aux prédiction-corrections ont également été développées. Il s'agit par exemple des méthodes "piecewise-linear" qui reposent sur la recherche de chemins polygonaux voisins des branches de solutions de $R(U)$. L'algorithme s'appuie sur une "triangulation" de l'espace $R^{N+1}$, des interpolations linéaires de $R(U)$ entre les sommets de la triangulation et des pas continuation qui s'apparentent à une méthode de Simplex[10]. Ces techniques ne réclament pas de régularité sur les équations, elles ont été appliquées avec succès pour le flambement de poutres confinées entre deux murs par Holmes et al. [11]. Pour la mécanique des structures, on citera sans chercher à être exhaustif deux autres méthodes alternatives, la méthode à grand imcrément de temps (LATIN) [12] due à Ladevèze et al. et la méthode asymptotique numérique (MAN) due à Damil, Potier-Ferry et Al., dont on trouvera les détails dans [13]. La méthode LATIN repose sur des corrections globales des solutions pour un intervalle donné du paramètre de continuation, et non pour une valeur unique du paramètre. Elle s'appuie sur une décomposition astucieuse en équations globales linéaires et locales non linéaires, avec une représentation adéquate des inconnues sur l'intervalle. Cette méthode est particulièrement bien adaptée aux relations de comportement complexes, aux analyses multi-échelle et à la sous-structuration. La méthode MAN repose quand à elle sur la représentation des branches sous forme de série entière à des ordres élevés. Ses points forts sont la facilité du pilotage et la robustesse de la continuation. Nous allons montrer, dans cette communication, qu'elle est bien adaptée pour les diagrammes de bifurcation qui présentent des bifurcations.

## 2 Continuation de diagrammes de bifurcation complexes

Dans certaines applications, le diagramme de bifurcation de $R(U) = 0$, peut se révéler particulièrement riche. Il est formé de nombreuses branches qui se connectent en divers points de bifurcation. C'est par exemple le cas typique d'une structure que l'on a optimisée vis à vis du flambement linéaire. Après optimisation, les premières charges critiques sont pratiquement confondues et les branches bifurquées primaires peuvent présenter des bifurcations secondaires avec des branches sous critiques. La structure devient alors sensible aux imperfections et elle flambe en pratique avant sa première charge critique. C'est le phénomène bien connu de "l'érosion" de l'optimisation de la charge critique [14]. Un autre exemple concerne la recherche des solutions périodiques forcées dans les problèmes de vibrations non linéaires. La figure 1 présente un diagramme de réponses forcées pour une poutre bi-encastrée [15]. On peut y voir de nombreuses branches correpondant aux résonnances principales et super-harmoniques autour des modes non linéaires du système.

Pour tracer des diagrammes de bifurcation compliqués, il faut être capable de detecter les bifurcations et de tourner sur les branches bifurquées. Ces opérations, qui reposent sur des
<figure><img src="image_3.png" /><figcaption>FIG. 1 – Exemple de diagramme avec de nombreuses bifurcations. A gauche : vue générale du diagramme. A droite : vue de détail</figcaption></figure>

concepts théoriques bien établis [16], sont disponibles dans les logiciels généraux de continuation les plus aboutis. On les trouve beaucoup plus rarement dans les programmes que les chercheurs et les ingénieurs écrivent pour leurs applications, ou dans les codes de calcul de structures du commerce. Il faut dire que la facilité de mise en oeuvre et la fiabilité de ces techniques doivent encore être améliorées pour un usage courant. On notera que ces opérations de branchement gagnent à être faites de façon interactive, avec un utilisateur qui prend les décisions au fur et à mesure qu'il découvre les branches de solutions. Un logiciel de continuation qui aurait pour ambition de tracer automatiquement l'ensemble des solutions est d'une part, très difficile à réaliser, d'autre part, d'un intérêt limité. Le logiciel risque en effet de passer son temps sur des réseaux de branches qui n'intéressent pas forcément l'utilisateur.

L'objet de cette communication est de dégager l'intérêt des méthodes de continuation MAN pour les diagrammes complexes. Nous allons montrer que la technique des perturbations aléatoires permet de faire la détection et le branchement d'une façon simple et fiable. La robustesse de la procédure tient au fait que, grâce aux séries entières, on dispose d'une moyen simple de régler l'intensité de la perturbation pour assurer le branchement ou le non branchement. Dans la suite, nous allons commencer par rappeler le principe de continuation MAN, puis nous aborderons le problème du branchement, et enfin nous présenterons le logiciel de continuation interactif Manlab.

### 3 Le principe de continuation d'une branche par la MAN

Soit $U_0$ un point solution régulier vérifiant $||R(U_0)|| < \varepsilon_R$ où $\varepsilon_R$ est un critère de tolérance donné. Soit $U_t$ le vecteur tangent (normé et orienté) en $U_0$, et $a$ le paramètre de pseudo-longueur d'arc défini par $a = U_t^T \cdot (U - U_0)$. On représente la branche qui passe en $U_0$ sous la forme d'une série entière en $a$ tronquée à l'ordre $N : U(a) = U_0 + aU_1 + a^2U_2 + \dots + a^N U_N$. Les vecteurs $U_i, i = 1, N$, sont solutions de systèmes linéaires récurrents à matrice invariante. Pour des ordres de troncature élevés, disons de 15 à 30, la série tronquée fournit une description locale continue de la branche qui est de très bonne qualité pour les valeurs du paramètre $a$ qui sont à l'intérieur du rayon de convergence de la série. On peut facilement contrôler la qualité de cette série tronquée eu utilisant l'approximation asymptotique suivante pour le résidu :

$$R(U_0 + aU_1 + \dots + a^N U_N) \approx R(U_0) + a^{N+1} R^{N+1}$$
où $\mathbf{R}(\mathbf{U}_0)$ est le résidu au point de départ et $\mathbf{R}^{N+1}$ le vecteur second membre à l'ordre $N + 1$. Si on tolère un accroissement de l'erreur de $\varepsilon_R$ par rapport à l'erreur initiale $||\mathbf{R}(\mathbf{U}_0)||$, c'est à dire si $||a^{N+1}\mathbf{R}^{N+1}|| < \varepsilon_R$, on obtient la valeur maximale admissible du paramètre $a$ qui est $a_{max} = \left( \frac{\varepsilon_R}{||\mathbf{R}^{N+1}||} \right)^{\frac{1}{N+1}}$. En résumé, sur l'intervalle $[0, a_{max}]$, la série tronquée fournit une représentation continue de la branche avec une erreur maximale de l'ordre<fnref n="1" /> de $||\mathbf{R}(\mathbf{U}_0)|| + \varepsilon_R$. Le point final $\mathbf{U}(a_{max})$ sert de point de départ pour une nouvelle représentation en série. On progresse ainsi le long de la branche qui est décrit par une succession de réprésentations en série valables localement. Cet algorithme de continuation s'apparente à une méthode de prédiction-correction, où le prédicteur est d'un ordre très élévé. On notera qu'il fournit une description locale continue de la branche, au lieu d'un seul point de prédiction. L'utilisation d'un correcteur est facultatif. Il est néanmoins utile pour corriger les erreurs qui s'accumulent au cours des pas et pour contrôler le résidu du point de départ en vue du branchement, comme nous allons le voir plus loin. Cet algorithme de continuation n'a que deux paramètres : l'ordre $N$ et le critère de tolérance $\varepsilon_R$. L'ordre de troncature $N$ n'ayant pas d'influence sur les résultats et on le fixe à 20 par défaut en général. Le choix d'un ordre adapté ne se justifie que pour les très gros systèmes où l'on cherche à optimiser le temps de calcul. En pratique, l'utilisateur ne fournit que de la tolérance $\varepsilon_R$. La longueur des pas de continuation est gérée automatiquement par l'algorithme, via la valeur de $a_{max}$. Des exemples de continuation seront présentés lors de la présentation orale.

## 4 Detection des bifurcations et branchement

La détection des bifurcations à été l'objet d'une littérature abondante. L'idée principale est de définir une fonction scalaire le long de la branche qui change de signe en chaque point de bifurcation. Il suffit alors de surveiller les valeurs de cette fonction (indicateur de bifurcation) lors de la continuation. On trouvera les indicateurs les plus classiques dans [4] et des versions adaptées à la MAN dans [17].

Le branchement de la continuation vers une branche bifurquée<fnref n="2" /> est également un sujet classique, mais cette possibilité est plus rarement disponible dans les logiciels de calcul de structures. On distingue deux grandes classes de méthodes. La première consiste à introduire une perturbation $\mathbf{P}$ dans le système d'équations qui devient $\mathbf{R}_{pert}(\mathbf{U}) = \mathbf{R}(\mathbf{U}) + \mathbf{P}$. Si le vecteur de perturbation $\mathbf{P}$ est aléatoire, alors avec une probabilité de un point de bifurcation disparait comme cela est montré sur la figure 2. Dans la pratique, on peut choisir un vecteur de perturbation adapté au lieu d'un vecteur aléatoire. C'est souvent le cas dans les problèmes de flambement où le vecteur de perturbation est choisi proportionnellement à un mode de flambement linéaire. Cette technique est très simple à mettre en oeuvre, puisque l'on est ramené à de la continuation d'une branche sans bifurcation pour effectuer le branchement.

Cependant, avec une méthode de prédiction-correction, le choix de l'amplitude de la perturbation réclame un certain "doigté". La perturbation doit être suffisamment petite pour ne pas modifier trop fortement le système. Elle ne doit pas être trop petite non plus, sous peine de voir l'algorithme sauter sur la branche d'en face, et rater ainsi le branchement. La seconde façon de faire le branchement s'appuie sur la théorie de la bifurcation [16]. Elle repose sur un calcul précis du point de bifurcation puis des deux tangentes à chaque branche en ce point. On peut dès lors générer des points de prédiction adéquats pour embarquer la continuation sur l'une ou

<footnote n="1">¹On rappelle que le calcul de $a_{max}$ repose sur une approximation asymptotique du résidu.</footnote>

<footnote n="2">²On se restreint pour simplifier aux bifurcations simples.</footnote>
<figure><img src="image_4.png" /></figure>

FIG. 2 – Perturbation d’une bifurcation simple : les deux branches (en pointillés) qui se croisent au point de bifurcation sont remplacées par des branches sans bifurcation (en trait plein).

l'autre branche. Ces techniques sont efficaces, mais elles sont évidemment plus lourdes à mettre en oeuvre.

<figure><img src="image_5.png" /></figure>

FIG. 3 – Branchement de la continuation sur la branche bifurquée au moyen d’une petite perturbation $p = 10^{-8}$. L’algorithme MAN adapte automatiquement le pas pour tourner sur la bifurquée après une série de petit pas.

La MAN permet de redécouvrir la méthode des perturbations aléatoires comme un moyen fiable et robuste de faire à la fois la détection et le branchement. Même avec des perturbations très petites, l'algorithme de continuation MAN arrive à "négocier le virage" sans sauter sur l'autre branche. Nous allons d'abord illustrer cette propriété remarquable sur un exemple. On considère pour celà le système $R(u, \lambda) = (u - \lambda^2)(\lambda - (u - 2)^2) = 0$ dont les solutions sont formées de deux paraboles qui s'intersectent dans le plan $u,\lambda$. On ajoute une très petite perturbation $p = 10^{-8}$ dans ce système qui devient $R_{pert}(u, \lambda) = R(u, \lambda) + p = 0$. On choisit
une tolérance $\varepsilon_R = 10^{-10}$ sur le résidu qui est inférieure à cette perturbation $p$. Le résultat de continuation à partir du point de départ $(u, \lambda) = (4, -2)$ est montré sur la figure 3. On constate que l’algorithme de continuation raccourcit drastiquement la longueur de ces pas à l’approche de la zone de bifurcation, là où la branche change brutalement de direction. Avec des perturbations plus fortes ($p = 10^{-4}$, $\varepsilon_R = 10^{-6}$) ou encore plus petite ($p = 10^{-12}$, $\varepsilon_R = 10^{-14}$) on obtient le même type de résultat. Simplement, plus la perturbation est petite, plus le virage est serré et plus le nombre de petits pas augmente. En définitive, quelque soit l’intensité de la perturbation, tant que la tolérance demandée est inférieure à la perturbation, l’algorithme de continuation effectue le virage sans sauter sur la branche d’en face. Si maintenant, on supprime la perturbation ($p = 0$, $\varepsilon_R = 10^{-6}$), l’algorithme de continuation traverse la bifurcation sans difficulté<fnref n="3" /> comme cela est montré sur la figure 4.

<figure><img src="image_6.png" /><figcaption>FIG. 4 – Sans perturbation, et avec des points de départ dont le résidu est inférieur à $\varepsilon_R$, l’algorithme de continuation MAN passe les bifurcations</figcaption></figure>

Cet exemple est révélateur du fait qu'une perturbation, même très petite, introduit une situation dans les séries et peut modifier de façon importante le rayon de convergence, ainsi que le $a_{max}$ [18]. Nous allons préciser maintenant comment choisir l'intensité de la perturbation par rapport à $\varepsilon_R$ pour assurer le branchement sur la branche bifurquée comme à la figure 3, ou pour rester sur la branche fondamentale au delà du point de bifurcation comme à la figure 4. Nous utilisons pour cela l'exemple élémentaire d'une barre articulée soumise à une force verticale (figure 5). On note $u$ l'angle de rotation et $\lambda$ la force. L'équation d'équilibre de la barre est $R(u, \lambda) = u(1 - \lambda) = 0$. On ajoute une perturbation $P$ dans le système, on considère le point de départ (exact) $(u_0, \lambda_0) = (p, 0)$, et on prend le paramètre $\lambda$ comme paramètre de chemin.

Avec ces données, on obtient la série suivante $u(\lambda) = P + P\lambda + P\lambda^2 + ...$ dont le rayon de convergence est égal à un<fnref n="4" />. Le résidu de la série tronquée à l'ordre $N$ vaut exactement $R(P + ... + P\lambda^N) = P\lambda^{N+1}$ et la valeur maximale de $\lambda$ qui permet d'assurer la tolérance $\varepsilon_R$

<footnote n="3">³nous allons voir plus loin qu'il faut également maintenir le résidu des points de départ inférieur à $\varepsilon_R$</footnote>

<footnote n="4">⁴On notera que c'est exactement la distance du point de départ au point de bifurcation</footnote>
<figure><img src="image_7.png" /><figcaption>FIG. 5 – Exemple élémentaire d’une barre articulée en compression : diagramme de bifurcation avec et sans perturbation</figcaption></figure>

est $\lambda_{max} = (\frac{\epsilon_R}{P})^{\frac{1}{N+1}}$. On voit que si $p$ est plus grand que $\epsilon_R$, alors $\lambda_{max}$ est plus petit que 1 et la fin du premier pas est avant la bifurcation. Dans ce cas l’algorithme de continuation va décrire le virage. En revanche, si $p$ est plus petit que $\epsilon_R$, alors la fin du premier pas de continuation est déja au delà du point de bifurcation, et l’algorithme continue sur la branche verticale. on comprend sur cet exemple élémentaire que l’algorithme reste sur la branche fondamentale ou bien tourne sur la branche bifurquée, selon que l’intensité de la perturbation est plus petite ou plus grande que la précision demandée $\epsilon_R$.

On notera pour terminer qu’au cours de la continuation, les points de départs successifs ne sont pas des solutions exactes mais seulement des solutions approchées. Le résidu $||R(U_0)||$ au point de départ joue exactement le même rôle que le vecteur de perturbation **P**. Aussi, lorsque l’on fait de la continuation MAN de base, c’est à dire sans correction et sans perturbation, le résidu des points de départ devient de l’ordre de quelques fois la valeur de $\epsilon_R$ et l’algorithme tourne spontanément au point de bifurcation [18]. Si l’utilisateur ne veut pas tourner, il suffit de faire des corrections pour maintenir le résidu $||R(U_0)||$ à une valeur inférieure à $\epsilon_R$.

## 5 Le logiciel interactif MANLAB

Ce logiciel, développé par R. Arquier et B. Cochelin, comprend un moteur de calcul MAN et une interface graphique qui permet de réaliser de façon interactive les opérations nécessaires à la continuation et au branchement par perturbation. Le logiciel et sa notice sont disponibles sur le site web [?](https://manlab). Le logiciel est écrit dans l'environnement Matlab avec une programmation totalement objet. Il permet de traiter des systèmes algébriques réguliers quelconques de petite taille, disons typiquement de 10 à 100 équations. Il est également interfacé avec un code éléments finis "maison" pour aborder les problèmes de vibrations et de flambement des structures élastiques en nonlinéaire géométrique.

## 6 Conclusions

Dans ce papier, nous présentons une méthode (MAN) et un logiciel interactif de continuation dédié aux diagrammes de bifurcation complexes. Au travers d'exemples qui seront traités avec le logiciel lors de la présentation orale, nous soulignerons l'intérêt de ces techniques par rapport aux méthodes de prédiction-corrections classiques.
<figure><img src="image_8.png" /><figcaption>FIG. 6 – Le logiciel comporte une boîte à boutons, le diagramme de bifurcation, et des fenêtres de visualisation des solutions</figcaption></figure>

Références

[1] H. B. Keller *Lectures on numerical methods in bifurcation problems* Spinger Verlag, 1987.

[2] E.L. Allgower, K. Georg. *Numerical continuation methods, an introduction* Volumes 13 of Springer series in Computational Mathematics, Springer-Verlag, 1990.

[3] E. Doedel, H.B. Keller, J.P. Kernevez Bifurcation in finite dimension. *Int. J. of Bifurcation and chaos*, 1:493–520, 1991.

[4] R. Seydel *From equilibrium to chaos. Practical bifurcation and stability, analysis* Second Edition, Springer, New-York, 1994.

[5] Riks E. The application of Newton's method to the problem of elastic stability. *Journal of Applied Mechanics, Transaction of ASME* 1972; **39**:1060–1066.
[6] Ramm E. Strategies for tracing non-linear responses near limit points. In *Non-linear finite element analysis in structural mechanics*, Stein E., Wunderlich W., Bathe K.J. (eds). Springer-Verlag : New-York, 1981 ; 68-89.

[7] M.A. Crisfield *Non-linear finite element analysis of solids and structures*; Volume 1, John Wiley and Sons, New York, 1991.

[8] https://sourceforge.net/projects/auto2000/

[9] http://www.bifurcation.de/software.html

[10] E.L. Allgower, K. Georg. Continuation and path following. *Acta Numerica* 2:1-64, 1993.

[11] P. Holmes, G. Domokos, J. Schmitt, I. Szeberenyi. Constrained euler buckling : an interplay of computation and analysis *Comp. Methods. Applied Mech. Engrg* 170:175-207, 1999.

[12] P. Ladeveze *Nonlinear computational Structural Mechanics- New approaches and non-incremental methods of calculation*; Springer-Verlag, 1999.

[13] *Méthodes Asymptotiques Numériques* ; Numéro spécial de la Revue Européenne des Elements Finis, Vol 13, 2004.

[14] J.M.T. Thompson, G. Hunt *A general theory of elastic stability*; John Wiley and Sons, London, 1973.

[15] F. Perignon Vibrations forcées de structures minces élastiques non linéaires. *Thèse de l'Université de la méditerranée* ; Marseille, 2004.

[16] H. B. Keller Numerical solution of bifurcation and nonlinear eigenvalue problems. In *Application of bifurcation theory* P.H. Rabinowitz editor, Academic Press, New-York, London, 359-384, 1977.

[17] E. H. Boutyour, H. Zahrouni, M. Potier-Ferry, M. Boudi Asymptotic-numericam method for buckling analysis of shell structures with larges rotations *J. of Comput. and App. Math.* 168:77-85, 2004.

[18] On the behaviour of the ANM continuation in the presence of bifurcations *Comm. in Numer. Methods in Engng* 19:459-471, 2004.

[19] R. Arquier Manlab : Logiciel de continuation interactif (Notice) http://www.lma.cnrs-mrs.fr/manlab/, 2005.