le stade de Bunimovich et autres billards

mon dernier billet   sur le  billard de Sinaï m'a donné l'envie de continuer mon exploration de ce domaine, par exemple avec le stade de Bunimovich .  

le stade de Bunimovich

Comment simuler un billard?

Toutes ces animation sont générées avec scilab en modifiant un peu les fonctions utilisées pour simuler des billards polygonaux ou le billard de Sinaï  . Le problème dépend de 4 données :

• M= [M(1) M(2)] la position de départ
• V = [V(1) V(2)] la vitesse de départ
• T la durée de la simulation
• et bien sûr la forme du billard

 Le principe de l'algorithme de simulation est toujours le même :

fonction   Billard(M,V,T)
t=0, dt=pas de temps pour la simulation (0.1  à 0.02 suivant les cas)  
tant_que t<T
          si  detecte_collision(M,V,dt) alors 
                          calculer
                           V=vitesse modifiée après la collision
                          M= position au temps t+dt (donc après la collision)
                    sinon On avance en ligne droite
                              M=M+V*dt 
          fin_si
         afficher la position M ou le segment ...
          t=t+dt
fin_tant_que

la fonction qui détecte les collisions  va dépendre de la forme du billard  qui va être décrit par une (ou plusieurs) équation(s)  disons $f(M)\geq 0$ , dans ce cas on aura une structure de la forme :

fonction  detecte_collision(M,V,dt)
si f(M+V*dt)>0  alors pas de collision
sinon chercher le point de collision
X=M+s*V tel que f(X)=0  ( avec 0<s<dt)
fin_si

On peut se débrouiller pour que la fonction qui détecte les collisions  renvoie  les informations nécessaire pour calculer la nouvelle position et la nouvelle vitesse après le changement de trajectoire. Pour le changement  de vitesse il s'agit toujours d'une réflexion parfaite sur le bord  qui ne dépend  que de la normale  au point de collision. Si on veut que seule la fonction detecte_collision dépende du bord  (pour simplifier l'étude des différents billards) on peut utiliser  une fonction modifie_vitesse  qui calcule cette réflexion  à partir d'un vecteur  normal  (ou tangent) au bord :

fonction   modifie_vitesse(V,N)
N= vecteur normal unitaire (sinon diviser par ||N||)
ps=  V.N= V(1)* N(1) + V(2)*N(2) // produit scalaire de V et N
V=V-2 *ps*N

le stade de Bunimovich

le bord de ce billard est défini  par 2 demi-cercles dont les extrémités sont reliés par deux segments parallèles. Comme le billard de Sinaï ce billard  est dispersif (la plupart des trajectoire finisse par recouvrir tout le domaine spatial comme ci-dessus)et possède des propriétés de nature chaotiques. Pour le voir  on peut par exemple   lancer plusieurs particules de positions assez proches et voir qu'elles finissent assez rapidement par se répartir sur tout le domaine.

Pour programmer ce billard j'ai choisi de prendre pour le bord

• les arcs de cercle de rayon 1 et de centres (-1,0) et (1,0)
• les deux segments qui  relient les points (-1,1) à (1,1) et (-1,-1) à (1,-1)

la détection des collisions se fait alors suivant le schéma :

fonction   detecte_collision(M,V,dt)
X=M+dt*V
si |X(2)|>1  alors  collision sur un des 2 segments
   sinon si (X(1)-1)^2+X(2)^2>1 alors collision sur l'arc de cercle droit
   sinon si (X(1)+1)^2+X(2)^2>1 alors collision sur l'arc de cercle gauche
   sinon  pas de collision
fin_si

le billard circulaire

Ce billard peut être vu comme un billard  de Bunimovich pour lequel les deux segments parallèles sont de longueur nulle. On pourrait croire qu'il possède des propriétés analogues, pourtant il n'en est rien! Ici la plupart   des trajectoires recouvrent  une couronne  délimité par le bord extérieur et le cercle de même centre et tangent au premier segment de la trajectoire :

Il existe bien sûr tout une classe de trajectoires périodiques (par exemples beaucoup  forment un polygone régulier inscrit dans le cercle ) mais comme pour les autres billards une très faible erreur sur la position ou la vitesse de départ redonne une trajectoire couvrant une couronne . Dans l'exemple ci-dessous  le cercle est centré en (0,0) et de rayon R=1, de sorte que le carré de centre (0,0) et d’arête $\sqrt{2}$  est une trajectoire périodique  du billard  qui part du point $(\sqrt{2}, \sqrt{2})$ avec la vitesse $(-1,0)$. Si pour simuler une erreur d'arrondi on change  la coordonnée du point de départ en $(\sqrt{2}, 1.41)$ , la trajectoire ressemble à un carré mais qui semble tourner petit à petit  de sorte qu'il recouvre une couronne interne au cercle de départ.

le fait que ce comportement disparaisse totalement dans le stade de Bunimovich (quelque soit la longueur $>0$ des 2 segments  qui relie les arcs de cercle) montre l'importance de la régularité du bord. Dans le cas du cercle le bord est partout  $C^\infty$  alors que pour le stade il sont seulement $C^1$ et même pas $C^2$ et ceci uniquement à cause  points de raccordement demi-cercles/segments. La fonction de détection des collision est très simple à écrire si on prend comme billard le disque centré en (0,0) et de rayon R=1 :

fonction   detecte_collision(M,V,dt)
X=M+dt*V
si  X(1)^2+X(2)^2>1  alors  collision
 s= solution de  || M+s*V||^2=1 // trinôme du second degré
 X=M+s*V 
 N= X  est un vecteur normal  au point de collision  
 V'=modifie_vitesse(V,N)  // vitesse après collision
 X'= X+(dt-s)* V'  //  déplacement après collision
sinon pas de collision
fin_si

les billards elliptiques

Ce sont des généralisations du billard circulaire. On y retrouve  les trajectoires recouvrant une "couronne" (mais entre 2 ellipses cette fois)  lorsque le segment de départ ne coupe pas le segment joignant les deux foyers de l'ellipse.

Mais dans le cas contraire la trajectoire va couvrir une zone contenue  entre les branches d'une ellipse  de mêmes foyers que le bord elliptique :

Dans le cas du cercle cela correspondrait à une trajectoire passant par le centre du cercle ,  qui se réduirait au final à un rayon du cercle parcouru indéfiniment ... Pour une ellipse de demi-grand axes a et b l'équation du domaine est donnée par :
$$ {x^2\over a^2} + {y^2\over b^2}\leq 1
\Leftrightarrow
f(x,y)=1-\left({x^2\over a^2} + {y^2\over b^2}\right)\geq 0$$
 ce qui permet d'exprimer  la normale en un point du bord par :
$${\bf  N}=\nabla f (x,y)=\left(2{x\over a^2} ; 2 {y\over b^2}\right)$$
pour écrire la fonction qui détecte les collisions.