samedi 27 octobre 2012

simulation de trajectoires d'un billard polygonal

J'ai lu cette semaine un article consacré aux billards polygonaux sur le blog images des maths. L'article traite surtout de l'existence de trajectoires périodiques pour un billard donné et propose un petite animation d'une telle trajectoire. Ça m'a donné envie d'écrire un petit programme, avec Scilab, pour simuler de telles trajectoires :





Le problème est assez simple à modéliser (plus simple que le problème des trajectoires planétaires par exemple) il dépend de 4 données :

  • M=la position du mobile au départ
  • V=vitesse du mobile au départ
  • T= durée de la simulation
  • Bord=la liste des points délimitant le billard

Le mobile avance en ligne droite à vitesse constante jusqu'à entrer en collision avec le bord du billard, ce qui va modifier sa vitesse jusqu'à la prochaine collision ... l'algorithme s'écrit simplement :

 t=0, dt=pas de temps pour la simulation
tant_que t<T
          si collision(Bord,M,V,dt)  alors 
                          M=point de collision avec le bord
                          V=vitesse après la collision
                    sinon M=M+V*dt 
          fin_si
         afficher la position M
          t=t+dt
fin_tant_que


 le Bord du billard va être décrit dans scilab par une liste de points Bord=$[A_1;A_2; \dots;A_n]$ chacun des $A_i$ étant décrit  par 2 coordonnées x,y, voici deux exemples de billard : 

Carre=[-1 -1;1,-1;1,1;-1,1];
Triangle=[0,1;-1,-0.5;1 -0.5];

une fois identifié le morceau $[A_i,A_{i+1}]$ du bord du billard où va avoir lieu une collision, la modification de la vitesse résulte d'une simple réflexion spéculaire par rapport à ce bord, donc si on note U un vecteur normé orthogonal au segment $[A_i,A_{i+1}]$ alors la vitesse V devient  V-2<V,U>U après la collision

function V=modifie_vitesse(A, B, V)
    A=A(:),B=B(:),V=V(:)
    //U= vecteur  normé orthogonal à [AB]
    U=A-B
    d=sqrt(sum(U.^2))
    U=[U(2); -U(1)]/d
    //vitesse après réflexion V=V-2<V,U>U
    ps=sum(V.*U)
    V=V-2*ps*U
endfunction
 

La partie la plus difficile est de détecter la collision, cela revient à trouver un point d'intersection entre un morceau du bord $[A_i,A_{i+1}]$  et la droite définie par le point M et la vitesse V, ce qui s'obtient en résolvant un système linéaire de deux équations pour chacun des segments du bord. Là c'est un peu plus compliqué à programmer :

function [X]=intersection(A, B, M, V)
    //X= intersection entre (AB) et {M+Vt|t dans R}
    //   [ ] sinon
    xA=A(1);yA=A(2)
    xB=B(1);yB=B(2)
    //matrice du système 
    N=[yA-yB, -(xA-xB); V(2), -V(1)]
    Y=[xA*(yA-yB)-yA*(xA-xB); M(1)*V(2)-M(2)*V(1)]
    if abs(det(N))>10^(-10) then 
        X=linsolve(N,-Y)
        if (X(1)<min(xA,xB)-0.1)|(X(1)>max(xA,xB)+0.1)|(X(2)<min(yA,yB)-0.1)|(X(2)>max(yA,yB)+0.1) then
            X=[]// le point de collision n'est pas sur le segment!
        end
    else X=[]//si droites // pas d'intersection
    end
endfunction

function [bool, X]=collision(A, B, M, V, dt)
    M=M(:);V=V(:);// M et V vecteurs de même taille
    [X]=intersection(A,B,M,V)
    if X==[] then bool=%F
    else // si intersection vérifier que
         d=sqrt((X(1)-M(1))^2+(X(2)-M(2))^2)// distance à la collision
         s=sum((X-M).*V)// s >0 la collision à lieu dans le futur
        if (d<=1.1*dt*sqrt(sum(V.^2)))&(s>1D-8) then //collision dans moins de 1 pas de temps
            bool=%T,X=X(:)'
            else bool=%F,X=[]
        end
    end
endfunction

function [X, A, B]=detecte_collision(L, M, V, dt)
    // X = point de collision
    // A,B= bord du segment où à lieu la collision
    n=size(L,1)
    L($+1,:)=L(1,:)
    X=[],A=[],B=[]
    k=0
    while k<n //on teste tous les bords
        k=k+1
        [bool,X]=collision(L(k,:),L(k+1,:),M,V,dt) 
        A=L(k,:),B=L(k+1,:)
        if bool|(k==n) then break
        end
    end
endfunction

Enfin le programme principal s'écrit :

function billard(T, L, M, V)
    rect=[min(L(:,1)),min(L(:,2));max(L(:,1)),max(L(:,2))]
    t=0
    dt=0.1
    seuil=0.1*dt*sqrt(sum(V.^2))
    X=M(1),Y=M(2)
    while t<T
            [Z,A,B]=detecte_collision(L,M,V,dt)
            if Z<>[] then V0=V
                V=modifie_vitesse(A,B,V)
                M=Z
                if (abs(sqrt(sum((A-Z).^2)))<seuil)|(abs(sqrt(sum((B-Z).^2)))<seuil) then
                    V=-V0.*(0.95+0.1*grand(V0,'def'))// si collision dans un coin!!!
                end
            else //pas de collision on avance  
                M=M(:)',V=V(:)',M=M+V*dt 
            end
        X($+1)=M(1)
        Y($+1)=M(2)
        t=t+dt
        drawlater()
        clf
        xpoly(L(:,1),L(:,2),"lines",1);
        E=gce();E.thickness=2;
        plot(X,Y,'-c');
        E=gce();E.children.thickness=2;
        A=gca();A.data_bounds=rect;A.axes_visible=["off" "off" "off"];A.isoview="on";
        xinfo('t='+string(t))
        sleep(10)
        drawnow()
    end
endfunction

on obtient la simulation de départ enchargent les différentes fonctions dans Scilab puis en lançant  la commande

-->billard(37.2,Triangle,[0,0],[1,0.1])

On remarquera qu'il y a des complications liées à plusieurs cas particuliers :
  • quand la trajectoire est trop rasant par rapport au bord les calculs peuvent être erronés, on peut le détecter facilement en regardant si le déterminant du système est petit
  • il faut vérifier que la collision a bien lieu dans le futur  (en regardant si le produit scalaire entre V et le vecteur de M vers le point de collision est positif ) et dans au plus un pas de temps (facile en comparant  la distance entre M et le point de collision  et la vitesse )
  • quand la collision est trop proche  d'un coin du billard le calcul de la vitesse après la collision est en général faux ... dans ce cas là, pour plus de réalisme j'ai pris pour vitesse après la collision  l'opposé de V+ une perturbation aléatoire
le script devrait pouvoir marcher avec des billards plus complexes :
  • des billards non-convexes comme  celui ci-dessous Mais le programme actuel ne traite pas bien les réflexion sur les coins "non-convexes" ...
    Banane=[2 1; -2 1; -2 -2; -1 -2; -1 -1; 1 -1; 1 -2; 2 -2];
  • des billards "courbes" comme un disque, mais dans ce cas là il y a un très grand nombre de petits segments dans le bord ce qui ralentit beaucoup la fonction (et peux créer des problèmes de calcul des points d'intersection) ... il faudrait optimiser la recherche des collisions aux seuls segments proches du point courant.

11 commentaires:

  1. Bonjour, étant un débutant sur scilab, et voulant réaliser de tel simulation, serait-il possible de mettre ce programme en fichier .sci téléchargeable ? Ou du moins l'entièreté du code que l'on ai juste a copier coller, car j'ai pris les fonctions "bêtement" je les ai mis bout à bout mais le programme ne lance rien (une erreur de ma part j'en suis sur)

    Merci !

    RépondreSupprimer
    Réponses
    1. Pardon ! Suffit d'énoncer le problème pour qu'il se sublime ^^

      Ma démarche, faites un copié collé des codes !! dans l'ordre !!,

      enregistrez,

      puis dans la console de scilab faites un Get('/chemin/de/votre/fichier'), (chemin que vous trouverez facilement en retourntant dans SciNote et en cliquant droit sur l'onglet qui contient le programme),

      puis faites dans la console billard(37.2,[0,1;-1,-0.5;1 -0.5],[0,0],[1,0.1])

      Enjoy ! Et merci à toi car ça va m'aider ton blog ;)

      Supprimer
    2. tant mieux si ça t'a aidé, il faudrait que je fasse un billet sur le billard de Sinaï que j'ai programmé en modifiant un peu le programme de ce billet .

      Supprimer
  2. Bonjour, nous travaillons sur votre programme. Néanmoins, impossible pour nous de comprendre la fonction modifie vitesse.Serait il possible que vous nous l'expliquiez ?
    Merci d'avance

    RépondreSupprimer
  3. A chaque fois que la particule rencontre le bord du billard il faut modifier la direction de sa vitesse V, pour un choc élastique il suffit de connaître un vecteur normal unitaire U au point d'impact car alors la vitesse V est modifiée en V'=V-2U où :

    =V(1)U(1)+V(2)U(2)

    est le produit scalaire des vecteurs U et V. Dans la modélisation que j'ai utilisé le bord du billard est représenté par une liste de points, par exemple :

    Triangle=[0,1;-1,-0.5;1 -0.5];

    chaque point est représentés par 2 coordonnées sur 1 ligne de cette matrice, on a donc les points (0,1); (-1,-0.5) et (1,-0.5) ... Deux points consécutifs de cette liste (A et B) forme un segment du bord du billard . En soustrayant leurs coordonnées (U=A-B) on obtient un vecteur parallèle au bord. Un vecteur normal est alors donné par [U(2),-U(1)] qu'il faut diviser par sa longueur sqrt(sum(U.^2)) pour le rendre unitaire.

    RépondreSupprimer
    Réponses
    1. je me rends compte qu'une partie du texte à été mal interprété (les inférieur/supérieur sont pris pour des balises html !!!), j'aurais du écrire en LaTeX :

      la vitesse modifiée est $$V'=V-2\langle V,U\rangle U$$

      où on a écrit le produit scalaire usuel $$\langle V,U\rangle =V(1)U(1)+V(2)U(2) $$

      Supprimer
  4. Merci énormément, cela nous a beaucoup aidée !!!
    Bonne journée :-)

    RépondreSupprimer
  5. Bonjour, je travaille sur votre programme mais je n'arrive pas à comprendre votre function intersection. Pourriez-vous me l'expliquer ?

    Cordialement

    RépondreSupprimer
    Réponses
    1. 1)la fonction X=intersection(A, B, M, V) prend en entré :

      - un segment [A,B] délimitant bord du billard
      - l'état du point mouvant M=position et V=vitesse

      ces 4 A,B,M,V données sont donc des listes de 2 coordonnées.

      2)On veut savoir si la trajectoire (la droite M+V*t avec t>0) coupe le bord [A,B]. Pour ça on commence par regarder si la trajectoire n'est pas parallèle au bord en calculant le déterminant de V et du vecteur AB (qui forment la matrice N).

      3)Si le déterminant de N est non nul (à 10^-10 près pour tenir compte des erreurs d'arrondi) alors il y a un point d'intersection entre la droite M+V*t et la droite (AB). Pour trouver les coordonnées de l'intersection on résout le système d'équations linéaires correspondant au point d'intersection ( qui est N*X=Y=> X=linsolve(N,-Y) en scilab ).

      4) Cependant le point d'intersection n'appartient peut être pas au segment [A,B] ce qui veut dire qu'il n'y aurait pas de collision avec ce bord! On vérifie donc que les coordonnées de X sont bien dans le segment [A,B] (avec une tolérance de +- 0.1 pour éviter de perdre un point de collision) et on renvoie X, sinon on renvoie X=[]

      Supprimer
    2. Merci beaucoup ! Ça m'a énormément aidé

      Supprimer
  6. si vous cherchez d'autres exemples de Billards intéressants à simuler et d'autres exemple de code scilab voir les 2 autres billets :

    le Billard de Sinaï

    le stade de Bunimovich

    RépondreSupprimer

Pour écrire des formules mathématiques vous pouvez utiliser la syntaxe latex en mettant vos formules entre des "dollars" $ \$....\$ $ par exemple :
- $\sum_{n=1}^\infty {1\over n^2}={\pi^2\over 6}$ s'obtient avec \sum_{n=1}^\infty {1\over n^2}={\pi^2\over 6}
- $\mathbb R$ s'obtient avec {\mathbb R} et $\mathcal D$ s'obtient avec {\mathcal D}
- pour les crochets $\langle .,. \rangle$ dans les commentaires utilisez \langle .,. \rangle
vous pouvez écrire du html dans les commentaires :
- italique <i> ... </i> gras <b> ... </b>
- lien <a href="http://adresse "> .... </a>