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.