Les nouvelles technologies pour l’enseignement des mathématiques
Intégration des TICE dans l’enseignement des mathématiques

MathémaTICE, première revue en ligne destinée à promouvoir les TICE à travers l’enseignement des mathématiques.

Géométrie dynamique et visualisation des systèmes différentiels
Article mis en ligne le 30 avril 2017
dernière modification le 16 avril 2021

par Alain Busser, Batoul Saoud, Patrice Debrabant

A l’origine, les logiciels de géométrie dynamique ne sont pas conçus (pensés) pour le calcul différentiel. Mais ces logiciels ont évolué et sont maintenant dotés de fonctionnalités séquentielles qui leur permettent d’investir ce champ d’investigation. C’est d’autant plus vrai depuis l’intégration (récente) de la 3D.
On va s’intéresser en particulier à trois logiciels :

  • CaRMetal
  • GeoGebra
  • DGPad

On verra en particulier ce que la géométrie dynamique peut apporter à la « visualisation du chaos » (approche des attracteurs étranges).
On traitera différentes équations différentielles ordinaires et une équation aux dérivées partielles, celle de la chaleur.

Illustration : visualisation de l’attracteur étrange de Lorenz

gif réalisé à partir d’une animation CaRMetal

Note : L’article comprend plusieurs figures DGPad manipulables en ligne. Pour afficher certaines d’entre elles, il peut être nécessaire de réafficher la page à l’aide du navigateur (firefox conseillé). Les fichiers audio présentés ne sont pas intégrés correctement par le navigateur Safari (firefox ou Chrome conseillé).

Cet article peut être librement diffusé et son contenu réutilisé pour une utilisation non commerciale (contacter les auteurs pour une utilisation commerciale) suivant la licence CC-by-nc-sa

Plan de l’article :

  • Introduction
  • Equations différentielles ordinaires
    • Loi de refroidissement de Newton
    • Les oscillateurs
      • Le pendule simple
      • L’oscillateur de Van der Pol
      • L’oscillateur de Duffing
      • L’oscillateur de Lorenz
      • L’oscillateur de Chua
    • Equations de Lotka-Volterra
  • Equations aux dérivées partielles
    • L’équation de la chaleur

Introduction

La résolution approchée des systèmes d’équations du calcul différentiel passe par une discrétisation d’intervalles continus et la création automatisée de plusieurs points. En géométrie dynamique, ceux-ci doivent rester dynamiques par rapport aux conditions initiales et aux fonctions définissant le système (en pratique à tous leurs paramètres potentiels si l’on veut garder de la manipulation directe).

Dans CaRMetal, ces points sont créés par une boucle Javascript, et leur dynamisme est assuré par la composition dynamique des expressions de CaRMetal.

Dans GeoGebra et DGPad, on peut utiliser l’objet liste de points (appelé séquence de points dans GGB) qui est un objet propre du logiciel et s’actualise en conséquence.
Dans le cas de GeoGebra, on a une syntaxe particulière, de type récursif.
Dans le cas de DGPad, on a une syntaxe limpide de type Javascript (ou Blockly, voire tortue Blockly).

Avec GGB, on peut utiliser aussi l’ingénieux tableur, mais cette utilisation reste un peu anecdotique (dans ce cadre) car le tableur fait perdre un peu de dynamisme. On n’en parlera pas dans cet article.
On n’évoquera pas non plus la notation de Leibniz qui, en tant qu’outil implémenté, n’existe que dans CaRMetal [1]

Remarque : Dans cet article, pour résoudre les équations différentielles ordinaires on a toujours utilisé la méthode d’Euler. Les scripts peuvent être adaptés pour mettre en œuvre les méthodes de Runge-Kutta. Dans certains cas, on y gagnera une meilleure stabilité.



Petit échauffement : la loi de refroidissement de Newton y’=-Ky

La loi de refroidissement de Newton alimente régulièrement les sujets du BAC STI2D et des BTS Biologie.

Une recette de grand-mère permet de résoudre l’équation différentielle associée.
Voici par exemple un exercice du sujet du BAC STI2D Polynésie de juin 2013.

Le mioche ferait bien d’écouter sa grand-mère car elle a raison.

L’équation différentielle est $y’=0,04(20-y)$

La solution du problème est $80e^{-0,04t}+20$

On va pouvoir la comparer à la solution approchée dynamique obtenue par la méthode d’Euler.

En utilisant les infinitésimaux de Newton, toujours lui, on obtient $dy=0,04(20-y)dt$.

Avec CaRMetal

On crée un curseur dt entre 0 et 1. En réduisant dt, on remontera le temps tout en augmentant la précision.

Le script est le suivant :

On crée un point (nommé automatiquement par CaRMetal) et on retient son nom dans la variable Javascript P .
Puis on boucle :

  • Créer le point suivant et retenir son nom dans Q ;
  • Créer le segment (=relier P et Q) ;
  • P=Q (se décaler en Q).

On notera l’utilisation de l’underscore (conception d’Eric Hackenholz) qui remplace une variable par son contenu et permet d’augmenter considérablement la lisibilité du code.

On peut ensuite ajouter quelques éléments cosmétiques et s’intéresser à l’erreur relative par rapport à la valeur exacte.
La courbe exacte est tracée en vert.
La courbe de la valeur approchée est tracée en rouge.
On ajoute la courbe de l’erreur relative en bleu (en prenant une échelle non sensible au zoom de 100 pixel pour une erreur relative de 1 %). Cette erreur est très faible, l’équation différentielle a tendance à corriger les erreurs d’approximation !

Le dynamisme est obtenu par l’utilisation des expressions : les coordonnées étant écrites entre guillemets, elles sont automatiquement transformées en expressions (qui sont des objets dynamiques) par CaRMetal.

Avec GeoGebra

On ne va pas utiliser le calcul formel. Le but, c’est de préparer les outils que l’on utilisera dans une situation plus compliquée.

On crée le curseur dt et on utilise le Javascript global du curseur dt.

  1. function calculerErreur() {
  2.     ggbApplet.evalCommand("P0=(0,100)");
  3.     ggbApplet.setVisible("P0",false);
  4.         nbPoints=1000;
  5.         ggbApplet.setValue("er",0);
  6.         for (var j=1;j<nbPoints;j++) {
  7.                 ggbApplet.evalCommand("P"+j+"=("+j+"*dt,y(P"+(j-1)+")+0.04*(20-y(P"+(j-1)+"))*dt)");
  8.                 nerr=Math.abs((80*Math.exp(-0.04*ggbApplet.getXcoord("P"+j))+20-ggbApplet.getYcoord("P"+j)))/(80*Math.exp(-0.04*ggbApplet.getXcoord("P"+j))+20)*100;
  9.                 if (ggbApplet.getValue("er")<nerr) {
  10.                         ggbApplet.setValue("er",nerr);
  11.                 }
  12.                 ggbApplet.setVisible("P"+j,false);
  13.                 ggbApplet.evalCommand("s"+j+"=Segment[P"+(j-1)+",P"+j+"]");
  14.                 ggbApplet.setLabelVisible("s"+j,false);
  15.                 ggbApplet.setColor("s"+j,210,10,10);
  16.                 }
  17. }
  18. ggbApplet.registerObjectClickListener("Bouton1", "calculerErreur");

Télécharger

La fonction calculerErreur calcule l’erreur relative et trace la courbe par la méthode d’Euler. Cette fonction est lancée au clic sur le bouton Go !

Dans cet exemple, l’erreur n’est pas dynamique, et on a moins de puissance et de fluidité qu’avec CaRMetal. Au delà de 1000 points, le logiciel commence à souffrir.

Avec DGPad

Avec DGPad, on peut utiliser les listes de points et en pratique on obtient plus de dynamisme et un peu plus de puissance qu’avec CaRMetal (en outre, la figure dynamique peut être publiée sur une page web).
Ce dynamisme supplémentaire donne beaucoup de liberté mais doit être employé de façon raisonnée sous peine de demander des taches surhumaines à l’ordinateur. Par exemple, on peut créer un nombre de points proportionnel à 1/dt, mais ce nombre peut devenir trop grand. On a fait le choix ici de garder un nombre fixe de points, ce qui garantit une certaine stabilité. Mais on pourrait aussi faire un curseur dynamique pour le nombre de points.

Ouvrir la figure dynamique dans un autre onglet

On veut construire trois courbes de couleurs différentes, donc on construit trois listes de points.
Le code « classique » (sans utiliser les propriété Blockly) de la liste de points pour tracer la courbe exacte est le suivant :

  1. qab=[];
  2. qab.push([0,100]);
  3. nbPoints=5000;
  4. for (var i=1;i<nbPoints;i++) {
  5. qab.push([i*dt,20+80*exp(-0.04*i*dt)]);
  6. };
  7. qab

Télécharger

La liste de points est un tableau de tableaux de 2 coordonnées construit en Javascript et qui prend vie dans DGPad (sous forme d’une suite de points) après le dernier point-virgule.
Il s’agit d’un objet dynamique de DGPad (c’est une expression particulière).

On peut remarquer que les expressions (ou curseurs) présents sur la figure peuvent être lus sans passer par GetExpressionValue() (qui serait plus logique), au risque d’une collision avec le Javascript. C’est pratique, mais prudence...

La mise en forme se fait ultérieurement pour toute la suite de points (couleur, épaisseur, on relie les poins ou pas, etc).

Le code pour la courbe approchée est :

  1. ab=[];
  2. P=[0,100];
  3. tab.push(P);
  4. nbPoints=4000;
  5. for (var j=1;j<nbPoints;j++) {
  6. Q=[j*dt,P[1]+0.04*(20-P[1])*dt];
  7. tab.push(Q);
  8. P=Q;
  9. }
  10. tab

Télécharger

Le code pour la courbe approchée et la courbe de l’erreur relative est :

  1. tab=[];
  2. tabE=[];
  3. P=[0,100];
  4. PE=[0,0];
  5. tab.push(P);
  6. tabE.push(PE);
  7. nbPoints=4000;
  8. SetExpressionValue("er",0);
  9. for (var j=1;j<nbPoints;j++) {
  10. Q=[j*dt,P[1]+0.04*(20-P[1])*dt];
  11. QE=[j*dt,10000/pixel()*Math.abs(20+80*exp(-0.04*Q[0])-Q[1])/(20+80*exp(-0.04*Q[0]))];
  12. if (GetExpressionValue("er")<100*Math.abs(20+80*exp(-0.04*Q[0])-Q[1])/(20+80*exp(-0.04*Q[0]))) {
  13.         SetExpressionValue("er",100*Math.abs(20+80*exp(-0.04*Q[0])-Q[1])/(20+80*exp(-0.04*Q[0])));
  14.         }
  15. tab.push(Q);
  16. tabE.push(QE);
  17. P=Q;
  18. PE=QE;
  19. }
  20. tabE

Télécharger


La version avec les listes de points construites par les propriétés Blockly est :

Ouvrir la figure dynamique dans un autre onglet

La programmation visuelle par blocs demande un peu de patience (et de pratique). On a choisi de faire 3 listes.

Le script visuel de celle pour la courbe exacte :

Le script visuel pour la courbe construite par la méthode d’Euler :

Le script visuel pour la courbe de l’erreur relative :



La grand-mère a aussi l’habitude de couper le chauffage la nuit. Voici un autre exercice extrait du sujet de BAC STI2D (Métropole juin 2017) qui s’intéresse à la question.


Le gradient de température est plus faible, on peut comparer l’erreur relative.



Conclusion partielle : on voit que les trois logiciels fonctionnent et donnent le même résultat.
Pour comparer leur efficacité, dans l’exemple suivant on pendra un nombre de points proportionnel à 1/dt et on observera à partir de quelle valeur ça plante.


La grand-mère est aussi préoccupée par le taux d’alcool de son petit fils qui siphonne sa Bénédictine.
Voici un exercice du sujet de BAC STI2D (Polynésie septembre 2014) :

L’équation différentielle est $y’=2,5e^{-t}-y$

La solution du problème est $2,5te^{-t}$

Cette fois, l’erreur relative est croissante (mais elle reste faible).

On va prendre un nombre de points égal à 5/dt (de quoi tenir 5 h).

Avec CaRMetal

On obtient progressivement la courbe, ce qui est agréable et garantit contre le plantage. Pour obtenir tous les points (jusqu’à t=5) en pratique il ne faut pas descendre en dessous de dt = 0.002

Avec GeoGebra

On obtient la courbe entière au bout d’un certain temps. Et là encore, en pratique il ne faut pas descendre en dessous de dt = 0.002, sinon le logiciel se met à planter.

Avec DGPad

Ouvrir la figure dynamique dans un autre onglet

On peut aller jusqu’à dt = 0.0002

Avec la tortue de DGPad

Pour dessiner une courbe solution de l’équation différentielle $\frac{dy}{dx}=2,5e^{-x}-y$, la tortue va dessiner un polygone approchant la courbe. En zoomant sur un pas de tortue, on voit que la trigonométrie permet de calculer l’angle a (orientation de la tortue) et la longueur à parcourir (en rouge) :

En effet le triangle est rectangle et la trigonométrie permet de calculer l’angle (on « connaît » le côté adjacent dx et le côté opposé dy) comme arctangente de y’ (lui-même calculé par l’équation différentielle), puis l’hypoténuse rouge puisqu’on connaît l’angle et le côté adjacent dx.

Le script tortue donne ceci :

Voici le fichier obtenu. On peut non seulement modifier dx avec le curseur, mais aussi choisir par mouvement du point bleu, la quantité initiale d’éthanol, normalement nulle (selon l’énoncé) :

Bien entendu, cette technique peut aussi être utilisée pour les autres équations différentielles du premier ordre traitées dans cet article.



Verdict : du point de vue de la puissance, DGPad s’impose largement car il permet largement d’aller jusqu’à dt = 0.0002 (soit 10 fois plus bas) et permet d’obtenir une erreur relative dynamique et un nombre de points dynamique.
CaRMetal est moins puissant mais donne une fluidité qui le rend agréable à utiliser.
GeoGebra fait le job, mais semble moins commode à utiliser dans ce cas précis (il est possible que cela soit dû, au moins en partie, au manque d’aisance des auteurs avec le scripting dans ce logiciel).



A partir de maintenant, l’article sera plus libre et on traitera chaque nouvel exemple avec un seul logiciel selon l’inspiration du moment.

Oscillateurs

Le pendule simple

Lors des soirées pluvieuses, après la Bénédictine, la grand-mère et son petit-fils s’adonnent à la radiésthésie : Le pendule simple est un oscillateur mécanique très classique.

On va s’intéresser ici au modèle simplifié pour les petites oscillations pour lequel on connaît une solution exacte, ce qui permet de comparer la valeur exacte à la valeur approchée obtenue par la méthode d’Euler.

$\dfrac{d^2 \theta(t)}{dt^2}+\omega \theta^2 (t)=0$

On met en oeuvre la méthode d’Euler sur le vecteur $(\theta,\theta’)$ :

$(\theta,\theta’)’=(\theta’,-\omega^2 \theta)$

On trouvera des exemples intéressants de pendule 2D réalisés par Christian Mercat avec CaRMetal sur cette page de i2geo.

On va présenter ici une modélisation 3D du pendule qui exploite certains des nouveaux outils 3D de CaRMetal.

Avec CaRMetal

Le pendule bleu correspond à la valeur exacte, le pendule rouge à la valeur approchée par la méthode d’Euler.
Pour le fun, on a modélisé des pendules en forme de cube.
D’un point de vue technique, on code les sommets du cube en base 2 pour pouvoir factoriser le code de leur déplacement.

On pourra noter que l’on a recréé une fonction Move en 3D pour accélérer le mouvement. Cette fonction, nommée MoveEspace(), supprime le point avant de le recréer, ce qui oblige à recréer aussi les objets dépendant.

Le fichier n’est pas très réactif. On a préféré ne pas rendre dynamiques les coefficients oméga et dt (ce serait très simple de le faire, il suffirait d’entrer les GetExpressionValue dans la boucle principale).

On voit que les pendules évoluent de façon très proche.

gif réalisé à partir d’une figure CaRMetal

Réaliser un pendule 2D serait plus simple, et on bénéficierait de plus de dynamisme.


La notion de portrait de phase

Non, Phase n’est pas le petit nom de la grand-mère, qui n’a d’ailleurs jamais posé pour un portrait : Le portrait de phase est la (des) courbe(s) de variation de la vitesse θ’(t) en fonction de θ(t).
La géométrie dynamique permet d’en donner une nouvelle définition :
(pour nous) le portrait de phase est la courbe dynamique de θ’ en fonction de θ avec manipulation directe du point donnant les conditions initiales.

Quelques remarques :

1) Si une trajectoire de phase (= un état de la courbe dynamique de phase) est fermée, cela traduit un mouvement périodique, durant indéfiniment dans le temps.

2) Dans le cas d’un pendule simple, les trajectoires de phase sont des ellipses (des cercles pour $\omega=1$).

3) On peut interpréter les différentes trajectoires de phase comme des courbes de niveau constant (niveau d’énergie en général).
Dans le cas de l’oscillateur simple, les ellipses sont concentriques, ce qui est due à la linéarité de l’équation.

Influence de l’amortissement sur le portrait de phase :

Si on ajoute un terme d’amortissement $\gamma v$ à un pendule simple, le pendule n’oscille plus indéfiniment, mais il s’arrête au bout d’un temps dépendant de $\gamma$.
Sa trajectoire de phase n’est plus fermée. Et le temps n’est plus réversible (le système n’est plus invariant par changement de signe du temps).
L’amplitude de l’oscillation décroît avec le temps, car l’énergie mécanique du pendule diminue et tend vers 0, à cause de la dissipation sous forme d’énergie thermique due aux frottements.
La trajectoire de phase n’est plus fermée mais tend vers un point qui correspond à l’état stable du pendule. On appelle ce point un attracteur. [2]

Théorie de Rayleigh : Si le frottement est de type visqueux (force opposée au mouvement, en général proportionnelle à la vitesse), les oscillations sont amorties (voir le pendule). Pour qu’il y ait oscillateur, il faut donc que, pour les petits mouvements, on ait

  • ou bien un frottement négatif quand la position est petite (Van der Pol, 1920 ; typiquement v(1-x2))
  • ou bien une force de rappel positive quand la position est petite (Duffing, 1918 ; typiquement en x-x3)

On va examiner de plus près ces deux oscillateurs.


L’oscillateur de Van der Pol

Le grand-père est un pionnier de l’électronique : Il faisait des oscillateurs à triode selon la recette de van der Pol : L’oscillateur de Van der Pol est un oscillateur de nature électronique développé par le le physicien néerlandais Balthasar van der Pol. Son équation est :

$\dfrac{d^2 x(t)}{dt^2}+\gamma \dfrac{dx(t)}{dt}-\epsilon\omega(1-x^2(t))\dfrac{dx(t)}{dt}+\omega^2 x(t)=0$

On peut représenter son évolution avec DGPad :

code de la liste de points

  1. var tab=[];
  2. var x0=xinit;
  3. var v0=vinit;
  4. tab.push([0,x0]);
  5. nbPoints=1000;
  6. for (var j=1;j<nbPoints;j++) {
  7. x0+=v0*dt;
  8. v0+=(-x0*omega*omega+epsilon*omega*(1-x0*x0)*v0-gamma*v0)*dt;
  9. tab.push([j*dt,x0]);
  10. };
  11. tab

Télécharger

On voit bien l’influence des différents paramètres.

Portrait de phase (avec manipulation directe du point rouge) :

Pour $\epsilon \neq 0$ on observe une dynamique régulière vers un cycle limite (= attracteur).

Lassé de voir osciller le pendule de sa grand-mère, le petit-fils décide de se mettre à la clarinette : Le fait que les formes d’onde obtenues ressemblent à celle de la pression à l’intérieur du bec d’une clarinette, laisse en effet penser que l’équation de van der Pol pourrait être un modèle simple de l’anche de la clarinette. Pour le vérifier, on a fait varier le paramètre (ε ci-dessus, m ci-dessous) progressivement depuis 0 (silence) jusqu’à 5 (anche battante) puis à nouveau 0. Est-ce que le son obtenu ressemble à celui d’une clarinette ? Pour en juger, cliquer sur le bouton « play » ci-dessous :

En fait, pour le timbre et la longueur de l’attaque c’est plutôt réaliste, mais une vraie clarinette ne sonne pas aussi grave si on souffle plus fort : Pour la modéliser correctement il faudrait introduire un retard qui synchroniserait l’oscillateur sur le retard : C’est le rôle du tuyau auquel est attachée l’anche.

van der Pol et Python

Pour fabriquer le son ci-dessus, on a simplement reprogrammé la résolution de l’équation différentielle x’’=m(1-x2)x’-x en faisant varier progressivement le paramètre m. Il est nécessaire d’ajouter un faible bruit de fond (random.random()/32000.0) pour échapper à la stabilité de la solution (x=0,x’=0) donnée comme condition initiale : La présence de fluctuations dans le souffle du clarinettiste est peut-être nécessaire au fonctionnement de la clarinette. Voici le script :

  1. import wave, struct, math, random
  2.  
  3. fe = 44100.0 # hertz
  4. duree = 10.0       # secondes
  5. dt = 0.04
  6. x, y = 0.0,0.0
  7. m = 5.0
  8.  
  9.  
  10. wavef = wave.open('vanderpol.wav','w')
  11. wavef.setnchannels(1) # mono
  12. wavef.setsampwidth(2) #2 octets par piste
  13. wavef.setframerate(fe)
  14.  
  15. for i in range(int(duree * fe)):
  16.         m = 0.5-abs(i-duree*fe/2.0)/duree/fe
  17.         m *= 10.0
  18.         x += random.random()/32000.0
  19.         y += random.random()/32000.0
  20.         x,y = x+dt*y, y+dt*(-x+m*(1-x**2)*y)
  21.         value = int(2000*x)
  22.         data = struct.pack('<h', value)
  23.         wavef.writeframesraw( data )
  24.  
  25. wavef.writeframes('')
  26. wavef.close()

Télécharger

Et voici le script de la variante « van der Pol forcé » ci-dessous :

  1. import wave, struct, math, random
  2.  
  3. fe = 44100.0 # hertz
  4. duree = 10.0       # secondes
  5. freq = 110.0    # hertz
  6. dt = 0.04
  7. x, y = 0.0,0.0
  8. m = 5.0
  9.  
  10.  
  11. wavef = wave.open('vanderpol2.wav','w')
  12. wavef.setnchannels(1) # mono
  13. wavef.setsampwidth(2) #2 octets par piste
  14. wavef.setframerate(fe)
  15.  
  16. for i in range(int(duree * fe)):
  17.         m = 0.5-abs(i-duree*fe/2.0)/duree/fe
  18.         m *= 10.0
  19.         x += random.random()/32000.0
  20.         y += 0.02*math.cos(freq*math.pi*float(i)/float(fe))
  21.         x,y = x+dt*y, y+dt*(-x+m*(1-x**2)*y)
  22.         value = int(2000*x)
  23.         data = struct.pack('<h', value)
  24.         wavef.writeframesraw( data )
  25.  
  26. wavef.writeframes('')
  27. wavef.close()

Télécharger

Il s’agit ici de chanter la note « la grave » (à 110 Hz) dans la « clarinette ».

Par contre, l’équation de van der Pol forcée permet de simuler le growl (musique) à la clarinette : Le clarinettiste chante dans la clarinette au lieu de simplement souffler dedans ’avec la Bénédictine c’est plus facile). Avec le son glissant, on constate à quel point l’effet est rauque vers le milieu de l’enregistrement (anche franchement battante) :

La présence de chaos dans les équations différentielles forcées (voir ci-dessous) permet aussi de modéliser les sons rauques produits par les violonistes débutants (pression de l’archet trop grande). En effet, une variante de l’équation de van der Pol a été proposée par Rayleigh à la fin du XIXe siècle pour modéliser les vibrations dues aux frottements.


L’oscillateur de Duffing

L’oscillateur de Duffing est un modèle d’oscillateur non linéaire. Il peut être simulé électroniquement lui aussi.

Version de base de l’oscillateur :

L’équation de cet oscillateur sans frottement est :

$\dfrac{d^2 x(t)}{dt^2}- x(t)+x^3(t)=0$

On obtient cette courbe de x :

Et le portrait de phase est le suivant, caractéristique de l’oscillation (l’attracteur est une courbe fermée) :

On va considérer maintenant la version de cet oscillateur amorti par un frottement fluide, puis celle où l’on compense cet amortissement par un forçage sinusoïdal.

Version de l’oscillateur amorti par un frottement fluide :

$\dfrac{d^2 x(t)}{dt^2}+0,1 \dfrac{dx(t)}{dt}- x(t)+x^3(t)=0$

Le portrait de phase est particulièrement intéressant :

Il s’agit d’un oscillateur à deux puits : on a des points fixes pour $– x + x^3 = 0$, soit $x=$ ±1 qui sont des centres pour l’oscillateur non amorti et des puits si celui-ci est amorti, et $x=$ 0, qui est un point-selle.

Du fait des frottements, l’énergie mécanique diminue, le point tourne dans le sens des aiguilles d’une montre en se rapprochant ; à un instant donné, il entre dans l’un des deux puits et tourne autour de celui-ci jusqu’à sa position finale au fond du puits. Il n’y a pas de chaos possible (pas d’attracteur étrange), sauf si on « force » l’oscillateur, ce qui revient à lui donner un second membre non nul.

Version de l’oscillateur amorti par un frottement fluide compensé par un forçage sinusoïdal :

On introduit un forçage sinusoïdal r.sin(t).
L’équation devient :

$\dfrac{d^2 x(t)}{dt^2}+0,1 \dfrac{dx(t)}{dt}- x(t)+x^3(t)=r.sin(t)$

On obtient cette courbe de variation de x :

Le portrait de phase est :

La trajectoire est alors chaotique. L’oscillateur tourne tantôt autour d’un puits, tantôt autour de l’autre, ou autour des deux, sans parvenir à se fixer. On voit émerger l’attracteur étrange de Duffing.


Les attracteurs étranges sont des ensembles difficiles à définir. Et on ne va pas le faire ici, c’est étrange...

Selon Wikipédia :
La forme d’un attracteur étrange « n’est pas une courbe ni une surface et n’est même pas continue mais reconstituée point par point de manière discontinue par la dynamique qui, bien qu’apparemment désordonnée, reconstitue ce type spécial d’ordre ». C’est un objet mathématique abstrait (c’est-à-dire qu’il ne peut être observé dans la nature) qui modélise un ou des paramètres du système à l’étude. Même si la forme est dite « étrange », elle permet d’étudier des phénomènes apparemment désordonnés qui sont influencés par des contraintes déterministes. La stabilité de cet attracteur est la conséquence de la structure sous-jacente du système. L’attracteur étrange sert à « élucider les mécanismes fondamentaux de la turbulence, les réactions chimiques, les prévisions météorologiques, la génétique des populations bactériennes ».


L’oscillateur de Lorenz

La grand-mère ressent une douleur dans ses articulations. Pour savoir si c’est le signe d’une pluie à venir, elle consulte le bulletin météo. L’oscillateur de Lorenz (ou attracteur de Lorenz) est un oscillateur d’un type particulier, tout comme l’oscillateur de Chua (voir plus loin).
Il modélise le chaos. Dans une certaine mesure seulement...

Historique :
L’attracteur de Lorenz et les équations associées ont été rendues publiques en 1963 par Edward Lorenz.
Depuis quelques années, Lorenz étudiait des modèles simplifiés du système de Navier Stokes de la mécanique des fluides (ce système d’équations était beaucoup trop compliqué à résoudre numériquement pour les ordinateurs existant au temps de Lorenz) décrivant le mouvement de l’atmosphère à travers des équations différentielles ordinaires dépendant d’un petit nombre de variables.

Correctif :
En fait, le système de Lorenz décrit les mouvements de convection (dûs à un gradient de température) alors que le système de Navier Stokes décrit les mouvements provoqués par un gradient de pression. Mais les équations étant analogues, on peut considérer que, formellement, c’est un cas particulier de Navier Stokes.

Une équation différentielle typique présentant à la fois viscosité et forçage est de la forme suivante :
$\dfrac{dx_i}{dt}=\sum_{\substack{i,k}}a_{i,j,k}x_jx_k-\sum_{}b_{i,j}x_j+c_i$

Lorenz observe que le champ de vecteurs est transverse aux sphères de grands rayons si bien que les trajectoires qui pénètrent dans une grande boule y restent confinées. Il peut alors discuter de diverses notions de stabilité, familières aux mathématiciens d’aujourd’hui ; périodicité, quasi-périodicité, stabilité ...
Puis, Lorenz prend l’exemple du phénomène de convection. Une couche d’un fluide visqueux est placée entre deux plans horizontaux, qui sont à deux températures différentes et il s’agit de décrire le mouvement qui en résulte. Les parties hautes du fluide sont plus froides et donc plus denses ; elles ont donc tendance à descendre par gravité et elles sont réchauffées lorsqu’elles parviennent dans les zones basses. La circulation du fluide qui en résulte est complexe.
En simplifiant, Lorenz parvient à obtenir une équation différentielle ordinaire décrivant l’évolution :

$\dfrac{dx}{dt}=-\sigma c + \sigma y$

$\dfrac{dy}{dt}=-xz+rx-y$

$\dfrac{dz}{dt}=xy-bz$

Ici $x$ représente l’intensité du mouvement de convection, $y$ représente la différence de température entre les courants ascendants et descendants et $z$ est proportionnel à la « distortion of the vertical temperature profile from linearity, a positive value indicating that the strongest gradients occur near the boundaries ».

En réalité, il ne faut pas chercher dans cette équation différentielle ordinaire une représentation fidèle du phénomène physique. La constante $σ$ est le nombre de Prandtl. Guidé par des considérations physiques, Lorenz est amené à choisir les valeurs numériques $r = 28$, $σ = 10$, $b = \dfrac{8}{3}$ ; c’était un bon choix et ces valeurs sont restées traditionnelles.
On peut alors résoudre numériquement ces équations et observer par exemple quelques orbites.

La nouveauté ici sera de construire des orbites dynamiques par rapport aux coefficients et à la position initiale.

Avec CaRMetal

Programmer la méthode d’Euler en 3D dans CaRMetal ne pose pas de problème particulier.

Tous les paramètres sont dynamiques :
 on peut modifier la couleur de la courbe (codes R,V,B) ;
 on peut modifier les coefficients (c, r, et e) ;
 on peut zoomer et modifier le repère tournant par clic droit ;
 En diminuant dt, on remonte dans le passé en augmentant la précision ;
 Le point initial est un point 3D libre que l’on peut déplacer dans le plan frontal (cette manipulation directe d’un point 3D n’existe que dans CaRMetal).

Quand on sent que le programme ralentit, on peut « tuer toutes les taches » pour interrompre le script. On obtient alors une figure plus réactive, le nombre de points créés n’évoluant plus.

Gif

réalisé avec des copies d’écran de CaRMetal

Comme on le voit sur le gif, pour que le papillon de Lorenz batte des ailes, on tourne le repère (c’est un effet d’optique, le papillon est figé).

C’est en 1972, lors d’une conférence intitulée « Predictability : does the flap of a butterfly’s wings in Brazil set off a tornado in Texas ? » que Lorenz popularisa « l’effet papillon ».
Au delà du principe général, qui est devenu un lieu commun pour bobos, on peut s’arrêter sur cette phrase : « More generally, I am proposing that over the years minuscule disturbances nei-ther increase nor decrease the frequency of occurrence of various weather events such as tornados ; the most they may do is to modify the sequence in which these events occur. »
On peut considérer que cette phrase a un véritable contenu scientifique car elle suggère (et invite à démontrer) que les statistiques décrivant le futur d’un système dynamique pourraient être insensibles aux conditions initiales.
C’est ce que l’on propose d’illustrer maintenant.

Pour commencer, on crée un cube avec l’outil correspondant de la palette.
Le centre du cube est un point 3D libre.
Ce cube est livré avec une expression k qui pilote la longueur des arêtes. On transforme cette expression en curseur.

L’idée est ensuite de construire la trajectoire issue de chaque sommet du cube et d’observer le comportement du faisceau.

Gif

gif réalisé à partir d’une figure CaRMetal

Tout d’abord, on voit bien la divergence très rapide du point jaune et du point blanc, cette divergence se faisant en restant dans le papillon.

Ensuite on recentre le papillon et on observe l’évolution en fonction de k (à t= 4,6)

Gif

gif réalisé à partir d’une figure CaRMetal

L’oscillateur de Chua

Le petit-fils a hérité de l’intérêt de son grand-père pour l’électronique ; mais lui n’utilise plus de triode, il préfère des « amplis-ops » [3] : L’oscillateur de Chua est un circuit électronique qui présente un attracteur étrange.

Pour l’explorer, on peut utiliser le dynamisme et la manipulation directe en 3D de CaRMetal. La méthode est la même que celle utilisée pour l’attracteur de Lorentz :

Sur la vidéo suivante, on a animé par script la rotation du repère et l’augmentation de dt, et superposé les portraits de phases (en 2D) :

  • en x (rouge)
  • en y (cyan)
  • en z (bleu)


Les équations de Lotka-Volterra

Les équations de Lotka-Volterra sont un classique des équations différentielles qui décrivent la dynamique de systèmes biologiques.

On se contentera ici de reporter le lecteur aux travaux d’Yves Martin sur la question (avec DGPad).

Les heureux (ou malheureux) possesseurs d’un iPad ou d’un ordinateur Apple pourront consulter la version iBook :

(Faire une recherche iBooks sur Yves Martin Volterra).



L’équation de la chaleur

L’équation de la chaleur est une équation aux dérivées partielles qui permet de décrire la conduction thermique et les évolutions de température qui en découlent.

En dimension 1

En dimension 1, cette équation s’écrit $\dfrac{\partial T}{\partial t}=K\dfrac{\partial^2 T}{{\partial x}^2}$

où T est la température (en °C), t le temps (en s), et K est un coefficient positif de conduction thermique.

Exemple avec des conditions de Dirichlet : on plonge un oeuf rectiligne de 1 cm de rayon à la température de 8°C dans de l’eau à 95°C et on veut calculer quand on a un oeuf dur.

On voudrait pouvoir considérer que l’on peut prendre des conditions de Dirichlet, à savoir que la température aux bords ne varie pas.
Mais mémé, par continuité cela implique qu’à t=0 l’intérieur de l’oeuf n’est pas exactement à 8°C. Il y a un état transitoire de durée 0, la condition initiale ne peut pas être $T_i(x)=8$, c’est impossible.
La question est donc de savoir quelle est cette condition initiale $T_i(x)$.

Les solutions générales de l’équation aux dérivées partielles sont de la forme

$d.e^{-a.t}cos(b.x) \quad $ et $\quad d.e^{-a.t}sin(b.x) \quad$ avec $\quad \dfrac{a}{b^2}=K$

Mais comme la fonction doit être paire en x (on centre l’oeuf dans le repère), on ne garde que les

$d.e^{-a.t}cos(b.x)$

Pour $x=R$ la température doit être constante, donc $b=\dfrac{(2n+1)\pi}{2R}=\dfrac{\pi}{2}+n\pi$

On va considérer que la conduction de la chaleur à l’intérieur du corps est plus lente qu’à sa surface. Et comme condition initiale, on va prendre la solution pour n=0.

On admet que $K=\dfrac{1}{677}$ (donné par les conditions physiques).

Par conséquent, la solution du problème est $T=95+(8-95)e^{-\frac{t}{274}}cos(\dfrac{\pi}{2}x)$

Quand le centre de l’oeuf est à 65°C, l’oeuf est cuit dur.
Un calcul simple montre que cela arrive en 4 min 52s (c’est un tantinet fantaisiste, mais depuis le départ ;)

On peut alors calculer une solution approchée en utilisant la méthode des différences finies, et la comparer à la solution exacte $T=95-87e^{-\frac{t}{274}}cos(\dfrac{\pi}{2}x)$

L’intervalle [0,1] est discrétisé en N + 1 noeuds de coordonnées $x_i$ (i variant de 0 à N ) régulièrement espacés.
On note ∆x le pas d’espace. Le temps est discrétisé en intervalles de pas constant ∆t.
On note $T_i^n$ la température au noeud $x_i = i\Delta x$ et à l’instant $t = n\Delta t$.

La résolution se fait en $t$ à partir de la condition initiale.
Par conséquent on utilise un schéma avant d’ordre 1 pour la dérivée en $t$ et un schéma centré d’ordre 2 pour la dérivée seconde en x :

$\left(\dfrac{\partial T}{\partial t}\right)_i^n=\dfrac{T_i^{n+1}-T_i^n}{\Delta t}$

$\left(\dfrac{\partial^2 T}{{\partial x}^2}\right)_i^n=\dfrac{T_{i+1}^n-2T_i^n+T_{i-1}^n}{\Delta x^2}$

On représente en vert la solution exacte et en rouge la solution approchée.

Avec CaRMetal

On doit renoncer à un peu de dynamisme, sinon la charge de travail est trop importante.
L’échelle étant radicalement différente pour x et pour T, on ne peut pas tracer les courbes en repère orthonormé. On peut

  1. soit tracer les courbes de T/100 ;
  2. soit utiliser un repère flottant pour tracer les courbes de T.

1) Dans le premier cas, le code est le suivant :

sol est la fonction solution à 2 variables, qui est entrée directement dans l’interface.

2) Dans le deuxième cas, c’est un peu plus technique, et l’utilisation d’un repère flottant dynamique ralentit le script.
On utilise l script fourni sur le site de CaRMetal pour tracer le repère flottant. Ensuite, on entre les formules.

Pour $\Delta t > 4$ (s), on observe un phénomène qui peut sembler étonnant.


Mais c’est un phénomène bien connu : il s’agit de la condition de Courant-Friedrichs-Lewy (CFL) : $\Delta t \leq \dfrac{(\Delta x)^2}{2K} $
En calculant, on obtient $\Delta t \leq 3,38$ ce qui est très probant ici.


Exemple avec des conditions de Neumann : Soit une barre de longueur $l=\pi$ parfaitement isolée thermiquement (= aucun échange de chaleur aux extrémités de la tige).
Mathématiquement, cela se traduit par des conditions de Neumann aux extrémités : $\dfrac{\partial T}{\partial x}(0,t)=0 \quad ; \quad \dfrac{\partial T}{\partial x}(l,t)=0 \quad\forall t \geq 0$

$t$ est le temps écoulé (en minutes), $x$ la longueur de la barre (en dm).

La condition initiale est $\quad T(x,0)=0,3 \pi x^2-0,2x^3 \quad $ pour $x \in [0, \pi ]$, condition qui vérifie les conditions de Neumann aux extrémités en t=0.

$\dfrac{\partial T}{\partial t}=K\dfrac{\partial^2 T}{{\partial x}^2}$

On admet que $K=1$ avec les unités choisies.

En utilisant la méthode de séparation de variables, on trouve que la solution est
$T(x,t)=\dfrac{a_0}{2}+\sum_{n=1}^{\infty}a_n.cos(nx)e^{-n^2t} \quad $ où

$a_0=\dfrac{2}{\pi} \int_0^{\pi}(0,3\pi x^2-0,2x^3) dx = \dfrac{\pi^3}{10} $

et $\quad a_n=\dfrac{2}{\pi} \int_0^{\pi}(0,3\pi x^2-0,2x^3).cos(nx)dx = \dfrac{0,2}{\pi}\left(\dfrac{(\pi^3.n^3+6\pi n)sin(n\pi)+12cos(n\pi)}{n^4}-\dfrac{12}{n^4}\right)$

On peut alors comparer la solution approchée obtenue avec les différences finies à la solution quasi-exacte constituée des 100 premiers termes de la série de Fourier.

La condition de Courant-Friedrichs-Lewy est : $\Delta t \leq \dfrac{(\Delta x)^2}{2} $


En dimension 2

En dimension 2, l’équation devient $\dfrac{\partial T}{\partial t}(x,y,t)=K(\dfrac{\partial^2 T}{{\partial x}^2}(x,y,t)+\dfrac{\partial^2 T}{{\partial y}^2}(x,y,t))$

On va s’intéresser à un problème avec conditions de Dirichlet :

$(x,y) \in \quad ]-0,5 ;0,5[^2$

La température est nulle et reste nulle à la frontière du domaine.

La condition initiale est la suivante :

En guise d’aperçu, voici un gif réalisé à partir de copies d’écran de la figure obtenue (avec CaRMetal) :

Gif

gif réalisé à partir d’une figure CaRMetal

Avec CaRMetal

On crée une nouvelle figure 3D.
Puis on peut construire sous forme d’une fonction CaRMetal à deux variables la condition initiale. Ici, il faut utiliser des if(,,) qui permettent de construire cette fonction (des mécanismes du même genre existent dans GeoGebra). Cette fonction est

if(abs(x)<0,4;if(abs(y)<0,4;1;(0,5-abs(y))/0,1);if(abs(y)<0,4;(0,5-abs(x))/0,1;(0,5-abs(x))*(0,5-abs(y))/0,01))

(Remarque : si cette condition devenait trop compliquée, on pourrait aussi la programer plutôt en Javascript dans le script.)

Ensuite, on adapte le script. Après initialisation, la boucle principale est :

  1. pour j allant de 1 à 10000 { // boucle temps
  2.         pour i allant de 1 à nbX-1 {
  3.                 pour p allant de 1 à nbX-1 {
  4.                         T[i][p]=exT[i][p]+K*dt/(dx*dx)*(exT[i-1][p]-2*exT[i][p]+exT[i+1][p]+exT[i][p-1]-2*exT[i][p]+exT[i][p+1]);
  5.                         Déplacer(P[i][p],a+i*dx,a+p*dx,T[i][p]);
  6.                
  7.                 }
  8.         }
  9.         pour i allant de 1 à nbX-1 {
  10.                 pour p allant de 1 à nbX-1 {
  11.                         exT[i][p]=T[i][p];
  12.                 }
  13.         }
  14.         MettreValeurExpression("tps",j*dt);
  15. }

Télécharger

Et voici le fichier CaRMetal :