Question de proba/simulation

Plop,

J'ai une envie de simulation mais je suis face à un problème théorique dont j'arrive pas à trouver une solution satisfaisante …

Alors voilà le cadre :

J'ai une fonction $f$ holomorphe et je cherche des cycles attractifs. C'est-à-dire des ensembles finis de valeurs qui forment un cycle (chacun s'envoie sur le suivant) et dont le produit de leurs dérivés est plus petit que 1 en valeur absolue.

Je suis gentil, je n'en demande qu'un et j'accepte une marge d'erreur contrôlable si ça peut se faire en un temps raisonnable.

Jusque là j'ai pensé à faire tourner des points au hasard pendant 30-40 itérations et essayer d'examiner les valeurs pour essayer de croiser un cycle … mais bof. Qu'est-ce que vous en pensez ? Des idées ?

Ty

mais bof. Qu'est-ce que vous en pensez ? Des idées ?

Bof pourquoi ? Faire un simulation te prendra du temps si tu fais plus de 10⁹ itérations. Donc 30-40, pas de problème =). Tu peut tester pas mal de points à l'avance comme ça. Si en plus tu as une restriction sur les points initiaux (le coup de la dérivée), c'est encore plus simple.

Après, est-ce que tes cycles sont aussi attractifs si tu passes "proche" d'un d'entre eux ? Et connais-tu quelque chose à l'avance sur leur taille ?

Le sens "attractif" n'est pas trompeur, localement les orbites s'accumulent près d'un tel cycle. D'où la possibilité d'avoir une certaine marge d'erreur.

Pour la taille, pas grand chose non …

Le problème est probablement là, j'ai dis 30-40 itérations en pensant que ce serait suffisant pour converger vers un cycle, mais par contre rien n'est moins sûr quant à ce que ce soit suffisant pour le détecter …

Effectivement, il n'est à ma connaissance pas du tout impossible que la convergence vers les cycles ne se fasse qu'à partir d'un certain rang. Avec le théorème de Montel, on devrait réussir à montrer qu'il existe des suites $(f_n)_{\in\mathbb N} de fonctions holomorphes qui convergent uniformément au voisinage des points attracteurs. Mais là, à froid je ne vois pas comment construire explicitement de telles suites.

Je suis content que tu parles du théorème de Montel, il est exactement question ici d'ensembles de Julia et Fatou.

Après au niveau des probabilités ont a de la chance : l'ensemble de Julia est de mesure nulle. Ainsi en prenant un point au hasard, il est statistiquement impossible d'avoir une orbite positive qui ne finit pas par s'accumuler près d'un cycle attractif !

Du coup j'aurai envie de prendre un premier point, construire un graphe de l'orbite où les noeuds sont à une marge d'erreur près et d'y chercher les cycles. À priori ça devrait marcher … peut-être que pour vérifier on pourrait prendre quelques points en plus qui seraient proches d'un cycle et dont on vérifierait la convergence vers ce cycle.

Réactions ?

Bonjour, Je pense avoir compris la question, chouette ! Mais dans ce cas, le nombre d'itérations retenues est insuffisant.

Imaginons la fonction :

Ro = (Ro+1) / 2

Theta= Theta + 2* PI / 100

Elle est holomorphe , et elle a une infinité de cycles attractifs, mais il faut 100 pts pour faire un cycle. A priori indétectable en 50 itérations.

Pour combiner performance et exhaustivité, je ferais ainsi : Fixer un nombre d'itérations à 40. Tester 5 points au hasard. Si pas de solution trouvée, augmenter le nombre d'itérations à 45, refaire une dizaine de tests avec d'autres points de départ (ou avec les 5 premiers, plus un 6ème pris a hasard, pour limiter les calculs), Et ainsi un certain temps , en augmentant régulièrement la profondeur d'analyse, et le nombre de points de départs analysés.

$f$ n'est pas seulement holomorphe, c'est aussi une fraction rationnelle. J'arrive pas à exprimer ta fonction comme une fraction rationnelle et du coup je tends à penser que ce n'est pas un bon exemple …

Je n'ai lu qu'en diagonale mais je tiens tout de même à dire que tu vas devoir te pencher sur la question de l'arrêt optimal. En effet, je ne sais pas quelle forme a ta fonction et la manière dont tu la calcules numériquement mais une fois tombé dans ton cycle, il va être possible d'en sortir à cause de la propagation des erreurs en arithmétique flottante, ce qui peut, si tu ne le détectes pas, t'empêcher de voir le cycle ou te faire penser que ce n'en était pas un.

C'est un problème récurrent sur les attracteurs.

La solution overkill serait CESTAC par exemple. Maintenant tu peux simplement te contenter de vérifications visuelles et d'une marge d'erreur variable pour tester le comportement face à ce problème (tu commences avec une marge d'erreur large, et tu la réduis jusqu'à voir si ta la trajectoire sort du cycle).

Mais si c'est un cycle attractif, je n'aurai pas ce problème, si ?

Si je comprends bien tu me dis que par les erreurs de calcul je pourrais m'éloigner d'un possible cycle, mais d'après l'attractivité, ça devrait être compensé non ?

Je connais pas CESTAC, des détails ?

Si ton cycle est attractif, effectivement tu y reviendras toujours en temps fini. Mais si je ne me trompe pas, la remarque de Höd reste valable : il n'est pas impossible que tu rentres dans un cycle, puis que tu en sortes au bout d'un certain temps à cause des erreurs de calcul. Tu y retourneras ensuite en temps fini, mais tu ne sais pas au bout de combien d'itérations…

En fait j'arrive pas à savoir ce que vous entendez par "sortir" du cycle ..

On y rentre jamais vraiment, ou du moins avec une proba nulle, et c'est bien sûr la composante connexe du bassin d'attraction qui m'intéresse ici. Une fois qu'on est dans cette composante, toutes les orbites s'accumulent près du cycle, et même si on s'en écarte un peu à un moment donné, on y revient puisque la composante est en principe de mesure suffisamment grande.

Après si je comprends bien, vous dîtes que je pourrais commettre une erreur tout juste suffisante pour que je crois que l'orbite n'induit pas de cycle. Alors je sais pas … quels sont les ordres de grandeurs de telles erreurs ?

quels sont les ordres de grandeurs de telles erreurs ?

Pas si important que ça, mais elles s'accumulent. Donc il faut voir comment le propagateur agit sur les erreurs avant de conclure.

Une idée d'implémentation du truc : générer des trajectoires pour un grand nombre de points de départ, puis analyser ces trajectoires pour y trouver les cycles à partir de propriétés topologiques de ces cycles. Chaque trajectoire peut faire facilement 10⁴ points, ce qui est un bon départ pour trouver les cycles. Ensuite pour l'analyse, il doit y avoir un moyen de faire ça plus intelligemment que juste en cherchant si la trajectoire repasse par un point particulier.

C'est peut-être de la déformation professionnelle, mais je vois une grosse analogie entre ces trajectoires et les trajectoires dans l'espace des phases pour un système à 1 degré de liberté. Peut-être une voie à creuser.

Sinon, elles ressemblent à quoi tes fonctions $f$ ? Ça aiderai peut-être à orienter la réflexion.

Il faut bien comprendre que lorsque tu travailles sur ordinateur, tu ne travailles pas sur le corps des réels mais sur le corps des flottants. Par la force des choses, l'ensemble des réels ne sont pas représentables informatiquement car il faudra une mémoire infinie. De fait, il y a des erreurs d'affections dû à la représentation même des nombres réels sur machine. A savoir par exemple que des nombres aussi simples que 0.1 ne sont PAS représentables en machine.

Je te donne un exemple très simplifié pour comprendre le problème. Imagine que tu es une fonction $f$ et $x_0 = 0.1$ et $x_1 = 0$ telle que $f(x_0) = x_1$ et $f(x_1) = x_0$. Tu as donc un cycle extrêmement basique. Sauf qu'au niveau informatique, $x_0$ n'est pas représentable et sera représentée par $X_0 = x_0 + \delta$, de fait que $f(X_0)$ ne vaudra pas exactement $x_1$. Selon la sensibilité de la fonction (au sens numérique notamment), tu pourras t'éloigner très rapidement du cycle pour ne jamais y revenir (puisque les erreurs se cumulent et se propagent). Tu pourras également obtenir des orbites totalement différentes de celles espérées (des lemniscates au lieu d’ellipses par exemple).

Voici un début d'introduction plus formel sur l'arithmétique stochastique qui est à la base de la méthode CESTAC (c'est plus pour la culture parce qu'à moins d'avoir de réels besoins, c'est lourd à mettre en place - utilisé dans l'industrie critique comme l'aéronautique par exemple -):

Le calcul numérique en nombre flottant est connu pour poser un certain nombre de problèmes de qualité des solutions obtenues, par plusieurs phénomènes :

• Entâchement d'erreur par la représentation du corps des réels (non dénombrable) par le corps des flottants (dénombrable), impliquant une erreur d'arrondi ou de troncature.
• La propagation de ces erreurs tout au long de la manipulation des flottants.

En résulte une précision détériorée et parfois des résultats qui n'ont plus aucune signification. S'il existe bon nombre de moyens à la portée de n'importe quel numéricien pour éviter une propagation trop rapide, il subsiste nécessairement une perte de précision, qui, selon le domaine, peut s'avérer dramatique, notamment pour les systèmes critiques. Ainsi, nous aimerions détecter et quantifier cette erreur (ou tout du moins la borner) et c'est ce que fait la méthode de Contrôle et Estimation STochastique des Arrondis de Calculs (CESTAC).

Soit $b$ la base de l'arithmétique dans laquelle nous travaillerons. Généralement, $b=2$ ou $b=16$ pour les unités de calculs. Alors, tout nombre $x\in \mathbb{R}$ est représenté par la machine de manière normalisée par un flottant $X$ de sorte que : $$X = \pm Mb^E$$ Avec $\frac{1}{b} \leq M < 1$ et $M$ la mantisse codée sur un nombre $n$ fixe de digits et $E$ l'exposant codé également sur un nombre fixe de digits.
On peut donc écrire $M$ en base $b$ : $\sum_{i=1}^n M_ib^{-i}, 0 \leq M_i < b$

A l'instar, on peut également décomposer $x$ en base $b$ : $$X = \pm mb^E$$ Avec cette fois une mantisse $m$ possédant un nombre (possiblement) illimité de chiffres.

On considère l'opération d'affectation ($\rightarrow$) : $X \rightarrow x, \forall x \in \mathbb{R}$. Alors, l'erreur relative d'affectation sera $\alpha = \frac{X - x}{X}$.

Si l'on pose $r = m - M$ le résidu, c'est à dire la partie perdue de la mantisse, on a $\alpha = \frac{-r}{M}$ avec un encadrement du résidu selon que l'on se trouve en arithmétique d'arrondi ou de troncature :

• $0 \leq r < \frac{1}{b^n}$ en arithmétique de troncature.
• $\frac{-1}{2b^n} \leq r < \frac{1}{2b^n}$ en arithmétique d'arrondi.

Partant du principe que tout résultat d'une opération arithmétique qui n'est pas un flottant est borné par deux flottants, $R^-$ et $R^+$, l'arithmétique stochastique consiste à retenir aléatoirement une de ces deux bornes avec une probabilité égale.
L'objectif de l'arithmétique stochastique n'est pas d'améliorer la précision mais d'estimer l'erreur commise et le nombre de chiffres significatifs en effectuant $N$ fois le même calcul.

De fait, avec CESTAC le principe est d'effectuer plusieurs fois (en pratique 3), chaque opération élémentaire (addition, multiplication, soustraction, division et affectation) avec un mode d'arrondi particulier, pour estimer avec un test de Student le nombre de chiffres significatifs.

J'ai un très bon article dessus mais malheureusement il est payant . Ceci dit, on peut s'arranger par e-mail.

Pas si important que ça, mais elles s'accumulent. Donc il faut voir comment le propagateur agit sur les erreurs avant de conclure.

Rien ne te permet de l'affirmer à priori. Il y a des systèmes où la propagation est énorme en quelques itérations et d'autres qui sont bien plus stable. Dans le cas de système linéaire on peut regarder le conditionnement du système par exemple. Je cherche un exemple parlant et j'édite ou poste un nouveau message avec.

Rien ne te permet de l'affirmer à priori. Il y a des systèmes où la propagation est énorme en quelques itérations et d'autres qui sont bien plus stable. Dans le cas de système linéaire on peut regarder le conditionnement du système par exemple. Je cherche un exemple parlant et j'édite ou poste un nouveau message avec.

Ce que je voulais dire, c'est que c'est plus la propagation des erreurs qui est gênante que l’existence d'erreurs. Et cette propagation dépends à la fois de :

• La fonction (mathématique) de propagation ;
• L'implémentation du calcul de cette propagation.

Il faudrait donc avoir la fonction (ou la classe de fonctions) étudiées avant de pouvoir être plus précis.

Effectivement, mais je pense qu'on est un peu hors sujet pour son problème. Je voulais juste attirer l'attention sur ce problème qui est souvent inconnu de ceux qui n'ont pas fait beaucoup d'implémentation. J'en ai vu s'arracher les cheveux pour comprendre pourquoi leur code ne connait pas les résultats escomptés alors qu'au final leur code était juste mais la stabilité numérique du schéma était à revoir.

Bref, si jamais à l'implémentation il trouve des choses bizarres il faudra qu'il se penche sur ces problèmes. Jusque là je pense qu'il peut se consacrer au problème initial.

Merci Höd pour ces résultat, je vais essayer de digérer ça demain mais dans les grandes lignes on s'est compris. Pour l'article … je verrais si j'en ressens le besoin et si je trouve pas d'autres articles sur le sujet (si t'as des liens mathscinet ce sera plus facile).

Pour Luthaf : toutes mes fonctions sont des fractions rationnelles à coefficients complexes (puisque ce sont les seules holomorphes sur $\bar{\mathbf{C}}$…).

Après j'ai envie de dire que les orbites sont en général très stables, une petite perturbation ne change pas le comportement asymptotique : convergence vers un cycle attractif (sauf au voisinage de l'ensemble de Julia, ce qui n'a pas de raison d'arriver ici).

Pour Luthaf : toutes mes fonctions sont des fractions rationnelles à coefficients complexes (puisque ce sont les seules holomorphes sur C¯…).

Ça fait longtemps que je n'ai pas fait de maths … Un truc de ce genre ?

$$ f(z) = \frac{az + b}{cz + d} $$

Yup, un rapport de deux polynômes (ce que tu donnes là c'est une fraction rationnelle de degré 1). Le degré d'une fraction rationnelle c'est le max des degrés du numérateur et dénominateur. On suppose en plus que les polynômes sont premiers entre eux (normal).

Dans le contexte, quand on parle de stabilité, on ne parlait pas de stabilité au sens mathématique mais au sens numérique. Or ici, ce genre de fonction va être plus ou moins stable numériquement en fonction des coefficients et de la manière dont tu la calcules.

C'est une propriété indépendante de la stabilité mathématique de ton attracteur.

Bon j'ai avancé dans tout ça. Au final :

• je tire des points au hasard (jusqu'à 100) ;
• pour chacun je cherche un cycle en faisant tourner le point (pendant 200 itérations au plus) et en cherchant une correspondance ;
• si je trouve rien j'indique que je regarde l'infini (mais bon, c'est bif bof je suis d'accord).

On a quelques premiers résultats qui sont sympas. Le tout marche assez lentement puisque maintenant je suis passé à des fractions rationnelles (les polynômes c'est moins de calculs). D'ailleurs si vous avez des méthodes pour rapidement calculer des valeurs selon une fraction rationnelle je dis pas non

Exemple de résultat avec $1 - 1/z^2$ :