jeudi 28 février 2019

polexpr et xint : applications (6)

Comment utiliser la macro \listplot de pst-plot avec xint (1)

Le texte ci-après est de Jean-François Burnol, l'auteur des extensions polexpr et xint :

Récemment, Jürgen Gilg et Thomas Söll m'ont sollicité pour voir comment utiliser la macro \listplot de pst-plot avec xint. Il se trouve que cette macro reçoit des options et comme argument obligatoire simplement des coordonnées (horizontale puis verticale) séparées par un caractère espace. De plus, ce matériel peut être obtenu par expansion à partir des macros, du moment que cela marcherait à l'intérieur d'un \edef.
Comme xintexpr procède de manière purement développable, ses macros sont bien adaptées à ce contexte, d'autant plus qu'il a une syntaxe pour générer des valeurs en progression arithmétique (par exemple t = 0..[3]..300) avec une variable t (la lettre est arbitraire) qui peut alors être utilisée comme paramètre dans des fonctions, disons par exemple seq((t, t^2),t=-10..[0.1]..10).
Un problème se pose de transformer les virgules en les caractères espaces attendus par \listplot, j'ai donc fait une macro \xintthelistplotpoints, c'est une cousine de la macro \xintthecoords de la documentation de xint (elle est plus simple car \listplot contrairement à TikZ ne limite pas le nombre maximal d'expansions).
Ensuite, avec Thomas et Jürgen on a commencé à regarder l'utilisation d'un pas variable puisqu'on peut parfaitement faire seq((t, t^2), t=-10..-5, -5..[0.1]..5, 5..10), et on peut aussi dans la même expression pour xint se faire succéder plusieurs tels seq((X(t),Y(t)), t=), simplement en les séparant par des virgules.
C'est alors que j'ai exploré une idée qui a en fait une riche mathématique sous-jacente, que je n'évoquerai qu'à peine ici : supposons qu'on en est à la valeur T du paramètre, comme choisir le suivant T + ΔT, d'une manière gouvernée par la courbure de la courbe ? Clairement plus la courbure est élevée (rayon de courbure petit) plus il faudrait un ΔT petit, et plus la courbure est plate plus on peut se contenter d'un ΔT (assez) grand.
L'idée suivie a été d'imposer que le nouveau point P + ΔP vérifie la contrainte que |ΔP|≤ η R avec un paramètre η, a priori sans dimensions, mais c'est plus subtil comme nous verrons sous peu. Bien sûr on ne va pas s'amuser à faire des tests jusqu'à trouver un ΔT qui marche donc on fait l'approximation |ΔP| = V ΔT et nous avons notre formule :
ΔT = η R/V.
Comme xint n'est pas un outil de calcul formel (on pourrait utiliser polexpr pour des paramétrisations polynomiales, voire rationnelles avec un peu de travail), l'utilisateur fournit sous formes de fonctions explicites (dans la syntaxe de xint qui est la notation "infix" des opérateurs) la paramétrisation de la courbe, le vecteur vitesse et le vecteur accélération ; ensuite xint calculera tout seul R et le module de la vitesse V. Bien sûr si on a fait tous les calculs à la main et simplifie les formules alors on peut les donner directement à xint.
Les images ci-jointes illustrent les résultats. L'idée de départ était simplement de choisir plus de points là où il faut pour que \listplot (avec plotstyle = curve) fasse une jolie interpolation par des quadriques de Bézier, et moins là où la courbe est plutôt plate.




Dans la pratique, les points d'inflexion (ou plus généralement les points où la courbure s'annule) posent un problème car le ΔT risque d'être grand et du coup la paramètre ira d'un seul coup peut-être au-delà d'un intervalle où la courbe tourne à nouveau. Donc un deuxième paramètre α est utilisé. Il a les dimensions d'une distance, et grosso modo si les zones plates de la courbe parcourt une distance D, il faudrait que α soit au plus D/2 ou D/4. Ça reste un peu artisanal. Il est utilisé via la contrainte ΔT≤ α/V, qui s'ajoute à la formule principale  ΔT ≤ η R/V.
On pourrait aussi simplement prendre un ΔT constant mais il serait dommage de le choisir inutilement grand...
La mathématique sous-jacent est très intéressante et on est obligé d'y réfléchir lorsque l'on fait des animations. L'animation aura N images par seconde. Chaque image passe d'un T au suivant T + ΔT. Numérotons les images par k = 0, 1, 2, .... À l'instant k/N nous sommes au paramètre T_k correspondant au point P_k, et à l'instant k/N + 1/N nous serons au point P_{k+1} et le vecteur de source P_k et de but P_{k+1} est à peu près le vecteur vitesse V_k multiplié par ΔT_k = η R/|V_k|. Cela veut dire que dans la nouvelle paramétrisation par les instants k/N, le nouveau vecteur vitesse est NηR fois le vecteur tangent unitaire. Autrement dit, la nouvelle vitesse est dans un rapport constant (égal à Nη) avec le rayon de courbure !
Cela confirme nos problèmes avec les points où la courbure s'annule, car la vitesse devient infinie.
Un petit calcul mathématique montre qu'avec cette nouvelle paramétrisation, le vecteur tangent unitaire tourne à une vitesse angulaire constante, qui est égale à Nη. Par exemple, avec η=0.1 et N = 10fps, la direction du vecteur vitesse tourne à 1 rad/s. Il peut aussi tourner à -1 rad/s bien sûr, si la concavité de la courbe change, ce changement brusque de la direction de rotation confirme les problèmes de notre paramétrisation aux points de courbure nulle.
Et au passage on voit que η, dans cette interprétation a des dimensions : il est en radians par image. Un η de 0,1 signifie qu'à chaque image, le vecteur vitesse aura tourné d'environ 0,1 radian. Je donnerai plus de détails dans un texte séparé, car la façon dont une variable sans dimension η se retrouve soudainement affecté d'une dimension est un phénomène de renormalisation, qui est tacite à toute notre discussion.
Dans la pratique, l'action du second paramètre α est de ralentir le point dans les zone plates, donc au total le temps de l'animation pour des courbes qui bouclent sur elles-mêmes peut être légèrement supérieur à la rotation absolue complète du vecteur vitesse divisé par la vitesse angulaire théorique de 1 rad/s.
Mais si vous contemplez les animations ci-jointes et jouez au chronométreur, vous obtiendrez des durées qui divisées par 2π auront des significations intéressantes...
Dans ces animations les pointillés verts indiquent la direction du vecteur accélération. En ce qui concerne le code utilisé il sera fourni prochainement car il nécessite un peu de nettoyage.
Pour les fonctions sind() et cosd() on a utilisé des macros ad-hoc, que l'on a regroupées dans un  fichier poormantrig.tex. Car actuellement xint n'a pas encore d'implémentation officielle.
Jean-François Burnol, Mercredi 27 février 2019




Comme le titre indique (1), il y aura une suite : Jean-François Burnol rédigera prochainement la présentation détaillée des problèmes (et des solutions) évoqués ci-dessus avec d'autres exemples et les fichiers sources. 

01/03/2019 : voici la présentation détaillée promise et écrite par Jean-François Burnol :

Le choix des points que nous envoyons à PSTricks via sa macro \listplot a été décrit plus haut :  on procède par des incréments de la paramétrisation de T en T + ∆T avec un ∆T donné par la formule η R/v où R est le rayon de courbure au point de la courbe de paramètre T et v la vitesse de parcours en ce même point. Il se trouve que ceci est une approximation à un procédé « idéal » qui utiliserait des pas constants unités dans une reparamétrisation de la courbe pour laquelle le (nouveau) vecteur vitesse tourne exactement à la vitesse angulaire de η rad/s ou - η rad/s (on change de signe aux points d'inflexion). Le texte mathématique ci-joint fait la théorie sous-jacente.

Paramétrisation par la (variation absolue de la) direction de déplacement

02/03/2019 : Le document  précédent étudiait la paramétrisation d'une courbe par la (variation absolue de la) direction de déplacement, dans le cas général. Voici une application dédiée aux ellipses, écrite par Jean-François Burnol, dans le document  suivant :


 04/03/2019 : les fichiers sources :
Tous les fichiers concernant les applications des extensions polexpr et xint de Jean-François Burnol présentées dans ce blog et des exemples d'utilisation de celles-ci avec PSTricks sont disponibles aussi sur :

mercredi 27 février 2019

Nouvelle version de pst-marble 1.4

La nouvelle version de pst-marble 1.4 est disponible sur le serveur du CTAN :
Je rappelle que l'auteur en est Aubrey Jaffer et que cette extension PSTricks a été adaptée par Juergen Gilg et Manuel Luque sous la supervision de l'auteur. Elle permet la réalisation de papiers marbrés en simulant différentes techniques utilisées par les artisans.
La page de l'auteur : http://people.csail.mit.edu/jaffer/Marbling/
La documentation vous renseignera sur les modifications et ajouts qui ont été faits par rapport à la version précédente. Voici quelques images extraites de la nouvelle documentation et des exemples qui lui sont joints.
 Sur les pages de ce blog vous trouverez de nombreux autres exemples, qu'il faudra peut-être adapter à cette nouvelle version, en faisant une recherche sur "pst-marble".

jeudi 21 février 2019

polexpr et xint : applications (5)

Le polynôme de Rump avec xint

Le polynôme de Rump est bien connu  de tous ceux qui se sont intéressés à la précision des calculs avec des nombres à virgule flottante et il est maintes fois cité sur tous les documents qu'on peut lire sur le sujet.
Voici le polynôme de Rump revu par Jean-François Burnol avec une présentation inédite, ce qui fait qu'on comprend, enfin très simplement, pourquoi ce polynôme est typique des problèmes rencontrés lors des calculs en virgule flottante. Jean-François Burnol s'est aussi intéressé  à l'ordre des opérations lorsqu'on fait le calcul, il a un rôle.
Il détermine aussi la période de la valeur exacte de f (77617, 33096)=-54767/66192.

Il a aussi créé une macro purement développable avec la syntaxe de xintexpr permettant de déterminer la longueur de la période décimale d'une fraction. La procédure est expliquée en détail dans le document. En plus de la valeur du polynôme de Rump, un autre exemple est calculé.
Vous pourrez  calculer avec 1000 000 de chiffres après la virgule 54767/66192, et enregistrer le résultat sur le disque, 3 lignes suffisent :

\immediate\openout\out=Rump-digits.txt
\immediate\write\out{\xintXTrunc{1000000}{\xinteval{Rump(77617,33096)}}

\immediate\closeout\out

Pour terminer cette présentation, je pense qu'après l'avoir lu, vous serez aussi admiratif que moi sur la qualité typographique du document qui avec sa pédagogie et sa rigueur mathématique en font  un document exceptionnel.

En voici le début  :

Son évaluation en x = 77617, y = 33096 est un exemple qui illustre l’un des problèmes typiques des calculs sur les nombres flottants, les « soustractions catastrophiques » : la soustraction de deux nombres très proches l’un de l’autre donnera zéro si l’arrondi de ces nombres à des mantisses de P chiffres est le même, et même si non nul aura en tout cas strictement moins que P chiffres. Or en virgule flottante, tout nombre flottant x doit être vu comme la classe d’équivalence des nombres réels x′ dont l’arrondi à P chiffres significatifs est x.
Le problème est que les x′−y′ exacts peuvent alors correspondre à beaucoup de classes d’équivalence mais que l’ordinateur ayant déjà arrondi x′ et y′ n’en obtiendra qu’une seule. Pour illustrer ce problème écrivons l’expression ci-dessus en mettant 4y en dénominateur commun, et en indiquant les valeurs numériques exactes :
On peut constater que les deux valeurs entières ont un nombre très important de chiffres en commun : 37 sur 43 ! Si l’on transforme ces entiers en des nombres flottants avec P < 37 chiffres, ils deviennent opposés et l’addition donne zéro. Nos calculs en virgule flottante n’utiliseront pas cette réduction à un dénominateur commun (comme le ferait éventuellement un logiciel de calcul formel) et le résultat final ne sera pas nécessairement zéro, mais en tout cas il sera très probablement « aberrant » si P n’est pas au moins 37.
 
 Fichier source et pdf sont accessibles ici :
 Les calculs ont été effectués avec l'extension xint, on appréciera dans ce document leur efficacité, elle est téléchargeable sur le CTAN :    https://ctan.org/pkg/xint

Tous les fichiers concernant les applications des extensions polexpr et xint de Jean-François Burnol présentées dans ce blog et des exemples d'utilisation de celles-ci avec PSTricks sont disponibles aussi sur :



mercredi 20 février 2019

polexpr et xint : applications (4)

Étude de suites avec XINT

Juergen Gilg propose une étude de suites en utilisant l'extension xint de Jean-François Burnol, on peut télécharger ce package sur le serveur du CTAN :
il contient différents modules qui sont utilisés dans ce document.
Il faut aussi récupérer l'extension poormanlog qui permet le calcul des logarithmes en base 10 et des puissances fractionnaires de 10 avec environ 9 chiffres de précision.
La documentation de XINT décrit toutes les commandes nécessaires permettant l'étude d'une grande variété de suites avec des paramètres permettant d'affiner cette étude ainsi de nombreux exemples.
Le document de Juergen Gilg illustre quelques cas  :
  • suites de valeurs de fonction dont les suites imbriquées avec nombres entiers ou décimaux avec une application du calcul du nombre d'Euler ;
  • suites définies par récurrence avec diverses applications comme les suites constantes, la méthode de Héron pour calculer la racine carrée. Voici le tableau obtenu pour sqrt(2) :
 et l'animation pour sqrt(10) réalisée avec le package animate d'Alexander Grahn :
 l'animation a été réalisée dans l'environnement PSTricks et les calculs effectués avec XINT.
D'autres applications complètent ce document comme l'algorithme de Brent-Salamin pour calculer π et des suites récursives.
Le document suites-xint.zip est dans le répertoire :
xint
Tous les fichiers concernant les applications des extensions polexpr et xint de Jean-François Burnol présentées dans ce blog et des exemples d'utilisation de celles-ci avec PSTricks sont disponibles aussi sur :

lundi 18 février 2019

polexpr et xint : applications (3)

La date de Pâques change d'une année à l'autre. De quelle manière est-elle déterminée ? Le site de la Libre Belgique :
en rappelle la définition suivante établie par le Concile de Nicée en 325 :

<< ce serait le premier dimanche après la pleine lune qui suit l'équinoxe de printemps, le 21 mars. Dès lors, Pâques tombe au plus tôt le 22 mars et au plus tard le 25 avril. À titre d'exemple, en 2014, Pâques a eu lieu le 20 avril; en 2015, elle a lieu le 5 avril; et en 2016, le 27 mars. >>
Il existe plusieurs méthodes pour calculer cette date, les auteurs de cet article(Jean-François Burnol et Juergen Gilg)  ont choisi la méthode du grand mathématicien Gauss et pour faire les calculs l'extension xint :
Une occasion pour ceux qui ne connaissent pas xint  de se familiariser avec sa syntaxe.
Les fichiers JouerAvecXint.tex et JouerAvecXint.pdf sont accessibles ici :
le tableau suivant, qui conclut le document, établit  la date de Pâques des années 2016 à 2019 ainsi que des fêtes qui en dépendent.
Tous les fichiers concernant les applications des extensions polexpr et xint de Jean-François Burnol présentées dans ce blog et des exemples d'utilisation de celles-ci avec PSTricks sont disponibles aussi sur :

dimanche 17 février 2019

L'encensoir de Botamufeiro (Juergen Gilg)

« Suspendu par une corde attachée au haut de la croisée de la cathédrale, l’encensoir gigantesque est dévié, d’une poussée, de sa position de repos verticale. Pendant qu’il se balance, huit hommes tirent sur une corde qui soulève l’encensoir lorsque celui-ci passe par la verticale et relâchent la corde  quand l’encensoir est au plus haut. Les tireurs, sous les ordres d’un conducteur, amplifient ainsi les oscillations de l’encensoir, jusqu’à ce que celui-ci monte à une hauteur de 21 mètres, décrit un arc de 65 mètres de longueur, et passe en vrombissant à la vitesse de 68 kilomètres à l’heure en un point situé à ras du sol. »
« . . . le dispositif actuel a, en son centre, deux tambours de châtaignier, de 58 centimètres et 29 centimètres de diamètre, dont l’axe commun repose sur la carcasse métallique. . . . L’encensoir actuel est en laiton argenté. »
« La corde dont une extrémité est dans les mains des tireurs s’enroule d’abord autour du tambour de petit diamètre, puis autour du tambour de grand diamètre, et redescend vers le sol ; la seconde extrémité est nouée à l’encensoir. Quand l’encensoir et la corde passent par la position verticale, les hommes tirent sur l’autre extrémité de la corde qui fait tourner le petit tambour de plus d’un tour et demi, lequel entraîne le grand tambour et la même corde soulève ainsi l’encensoir d’environ trois mètres. Ce dispositif à double tambour amplifie les déplacements. Les tireurs relâchent la même longueur de corde lorsque l’encensoir atteint le point d’amplitude maximale. Après la poussée initiale, l’amplitude angulaire est d’environ 13 degrés. En 80 secondes et 17 demi-périodes d’oscillation de l’encensoir, l’amplitude maximale de 82 degrés est atteinte, l’encensoir atteignant un point situé à un demi-mètre sous la voûte.»
Cité d’un article “Pour la science” numéro 155, de Juan Sanmartin Losada.
La vidéo d'une cérémonie du Botafumeiro à la cathédrale de Saint-Jacques-de-Compostelle :
Une animation et quelques figures de l'encensoir de Botafumeiro par Juergen Gilg.
Les fichiers sont dans le répertoire :


jeudi 14 février 2019

polexpr et xint : applications (2)

Le polynôme de Wilkinson 

Le polynôme de Wilkinson est d'apparence très simple, sous sa forme factorisée, il est défini par :
Il y a une très belle étude sur Wikipedia :
Le mathématicien  James H. Wilkinson qui a donné son nom au polynôme a titré son article : "The perfidious polynomial". Pour ceux qui souhaitent lire l'article original, voici le lien où on peut le télécharger.
Les extensions polexpr et xint de Jean-François Burnol ont déjà été présentées dans un précédent article :
elles sont disponibles sur le CTAN :

C'est polexpr qui a permis de réaliser la magnifique étude illustrée de ce polynôme, où on verra pourquoi J. H. Wilkinson le qualifia de ``perfide''. Pour les calculs avec les logarithmes une extension ``poormanlog.sty'' est nécessaire, elle est présente dans le répertoire et dans le fichier zippé.  Tous les fichiers sont dans ce répertoire :
Ce document est l’œuvre de Jean-François Burnol, Juergen Gilg et Thomas Soell. J'ai très modestement aidé à mettre au point les illustrations avec PSTricks, mais les calculs à partir de la forme développée des polynômes (voir le document) nécessitent une précision qui dépasse les possibilités de Postscript, la précision des calculs est discutée en détail dans le document. L'extension polexpr a permis de lever les limitations de Postscript et de dessiner les courbes.
Si un lecteur de cette page est aussi un utilisateur de Tickz et s'il souhaite en donner une version et la faire partager, il sera le bienvenu.
Je ne donne ici que quelques images, il faut lire le document et regarder le listing pour apprécier la puissance et la simplicité de l'extension polexpr. Par exemple, la syntaxe utilisée pour définir le polynôme s'écrit :
\poldef W(x):= mul(x-i, i=1..20);
L’appel $\PolTypeset{W}$ permet de l’afficher sous forme développée :
W(x) =
x20− 210x19 + 20615x18
− 1256850x17 + 53327946x16 − 1672280820x15
 + 40171771630x14 − 756111184500x13  +11310276995381x12
-135585182899530x11+1307535010540395x10
−10142299865511450x9+63030812099294896x8
−311333643161390640x7+1206647803780373360x6
-3599979517947607200x5+8037811822645051776x4
−12870931245150988800x3+13803759753640704000x2
−8752948036761600000x+2432902008176640000

 Wilkinson's polynomial
 sgn(W(x)) log(1 + |W(x)|)
Calculs effectués avec postscript

 Calculs effectués avec polexpr et poormanlog.
La particularité du polynôme de Wilkinson est d’être extrêmement sensible aux plus infimes variations des coefficients. On le voit sur le graphique ci-dessus, où la nouvelle courbe est celle des (logarithmes des) valeurs du polynôme obtenu par la modification du coefficient du terme en x19 en lui ajoutant −2−23. Cette minime variation sur le coefficient sous-dominant a eu un effet dévastateur à partir de x = 9 : les racines du polynôme, sauf celle proche de x = 20 sont devenues complexes !
L’extension polexpr ne sait pas faire les calculs avec des complexes mais est capable de déterminer avec une précision arbitraire tous les zéros réels et leurs multiplicités :
(mult. 1) Z1 = 0.99999999999999999999 . . .
(mult. 1) Z2 = 2.00000000000000000976 . . .
(mult. 1) Z3 = 2.99999999999980523297 . . .
(mult. 1) Z4 = 4.00000000026102318914 . . .
(mult. 1) Z5 = 4.99999992755153790956 . . .
(mult. 1) Z6 = 6.00000694395229570720 . . .
(mult. 1) Z7 = 6.99969723393601394867 . . .
(mult. 1) Z8 = 8.00726760345037685489 . . .
(mult. 1) Z9 = 8.91725024851707049429 . . .
(mult. 1) Z10 = 20.84690810148225691492 . . .

PS : le répertoire contient aussi un petit fichier de calculs de logarithmes avec l'extension poormanlog écrit par Juergen Gilg qui pourra éventuellement servir de mode d'emploi (tests-poormanlog.tex et pdf).

Tous les fichiers concernant les applications des extensions polexpr et xint de Jean-François Burnol présentées dans ce blog et des exemples d'utilisation de celles-ci avec PSTricks sont disponibles aussi sur :



jeudi 7 février 2019

polexpr et xint : applications (1)

polexpr et xint sont deux packages créés par Jean-François Burnol, disponibles sur le serveur du CTAN :

 xint qui permet les calculs sur les très grands nombres est constitué de plusieurs modules.
polexpr permet de faire du calcul formel et numérique avec des polynômes ayant des coefficients rationnels (ou en notation scientifique) :  calcul de dérivées, de PGCD, calcul du « contenu de Gauss » (générateur de l'idéal fractionnaire engendré par les coefficients), possibilité de trouver toutes les racines réelles avec une précision arbitrairement grande (et leurs multiplicités) et d'identifier au passage toutes les racines rationnelles, et de produire le polynôme obtenu en supprimant toutes les racines rationnelles, ou en supprimant toutes les multiplicités.
Pour le moment, polexpr (qui en est à la version 0.7.4) a les limitations suivantes :
  1. il ne connaît pas les nombres complexes,
  2. il ne connaît pas les fractions rationnelles (mais on peut ajouter des macros, comme dans le travail de Jürgen Gilg et Thomas Söll),
  3. il ne sait manipuler qu'une seule variable formelle à la fois,
  4. il représente intérieurement les polynômes par la suite de leurs coefficients, ce qui rend peu efficace la manipulation de polynômes de grands degrés comme x^1000-1.
 Les outils fournis ces packages sont très puissants et leur exécution très rapide. Jürgen Gilg et Thomas Söll séduits et enthousiasmés par les possibilités offertes pas ces packages ont écrit des applications à partir du travail de Jean-François Burnol et ont fait des liens, quand il  était utile d'avoir des représentations graphiques, avec PSTricks.
 Les commandes écrites par Jürgen Gilg et Thomas Söll sont dans le préambule des documents, ils envisagent d'en faire un package avec la documentation.
 Les documents, source .tex et pdf sont dans le répertoire :
Le fichier zippé les contient tous.
 Il m'est difficile de donner un aperçu de cette première série d'applications(j'espère que Jürgen Gilg et Thomas Söll en enverront d'autres) dans cette page, car elle n'est pas configurée pour afficher des équations mathématiques.
Dans polynom-example.tex, la courbe d'équation f'(x)=x^3+7x^2+14x+8 est étudiée, ses dérivées première et deuxième calculées ainsi que l'équation de la tangente en un point, PSTricks en fait les représentations graphiques.
 Dans le fichier Calculations-matrices.tex, vous y verrez une démonstration de la puissance du calcul matriciel exploitant les extraordinaires possibilités des packages polexpr et xint mises en œuvre par Jürgen Gilg et Thomas Söll.
 Dans le fichier RationalFunction.tex, ils étudient une fonction rationnelle, en font calculer les dérivées jusqu'au troisième ordre :

et dessiner les représentations graphiques par PSTricks :

Tous les fichiers concernant les applications des extensions polexpr et xint de Jean-François Burnol présentées dans ce blog et des exemples d'utilisation de celles-ci avec PSTricks sont disponibles aussi sur :