Intersection d'une surface reglée avec une sphère

Version imprimable

10/04/2007, 12h18
IXaMaXI

Intersection d'une surface reglée avec une sphère
Salut à tous,

J'utilise Matlab 7.0.0.19920 (R14) et j'essaie de calculer et représenter l'intersection entre une surface reglée et une sphère.

Ma surface a comme équations paramétriques :
x = u;
y = v;
z = (u² * v)/16;

Comme sphère j'ai choisi :
x = 2 * pi * sin(q) * sin(s);
y = 2 * pi * sin(q) * sin(s);
z = 2 * pi * cos(q);

J'ai donc dit que l'intersection devait être :
u = 2 * pi * sin(q) * sin(s);
v = 2 * pi * sin(q) * sin(s);
(u² * v)/16 = 2 * pi * cos(q);

Mais là déjà mathématiquement ça me paraît bizarre :/ M'enfin bon, je change ma dernière équation et j'obtiens :
(u² * v) / (32 * pi * cos(q)) = 1;

C'est bof...

Dans matlab j'ai mis ça :
Code:

1 2 3 4 5 6 7 8 [q,s] = meshgrid(0:pi./32:pi,0:pi./32:2.*pi); u=2.*pi.*sin(q).*cos(s); v=2.*pi.*sin(q).*sin(s); hold on; mesh(u, v, (u.^2).*v./(32.*cos(q))); hold off;
Donc on a notre sphère et notre surface reglée :
www.ixamaxi.be/images/matlab.png

Puis j'avais écrit ça :
Code:

1 2 3 4 5 6 7 8 [q,s] = meshgrid(0:pi./32:pi,0:pi./32:2.*pi); u=2.*pi.*sin(q).*cos(s); v=2.*pi.*sin(q).*sin(s); hold on; mesh(u, v, (u.^2).*v./16); hold off;
Ce qui me donnait :
www.ixamaxi.be/images/matlab1.png

Et ensuite avec le code donné plus haut :
www.ixamaxi.be/images/matlab2.png

Ce qui me perturbe c'est qu'on a une surface et pas une courbe. On s'attendait à voir une sorte de sinusoïde enroulée autour d'un cercle. Et c'est ce qu'on a dans le dernier screenshot mais on aimerait ne pas voir le cercle :D

Donc il doit y avoir soit une erreur dans notre routine matlab soit une erreur dans nos équations paramétriques.

J'imagine que les gens qui manipulent matlab depuis plus longtemps que 3 heures, doivent aussi savoir répondre au problème des équations paramétriques donc j'espère avoir une réponse :)

D'avance merci,
@micalement,
`ixM

Je ne suis pas un spécialiste des ces surfaces mais... les points doivent se trouver sur la sphère, donc :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
[q,s] = meshgrid(0:pi./32:pi,0:pi./32:2.*pi);
u=2.*pi.*sin(q).*cos(s);
v=2.*pi.*sin(q).*sin(s);
 
% Transformations des coordonnées cartésiennes en polaires
r=sqrt(u(:).^2+v(:).^2);
 
% Recherche des indices des valeurs de r égale à 2*pi (points sur la sphère)
% Comme 2*pi est un nombre réél, on ne fait pas r==2*pi mais abs(r-2*pi)<epsilon
idx=find(abs(r-2*pi)<eps);
 
% Affichage des indices trouvés uniquement
plot3(u(idx),v(idx),(u(idx).^2).*v(idx)./(32.*cos(q(idx))))

Je ne suis pas sur :?

10/04/2007, 14h07
IXaMaXI

Ok super merci !

Je venais de capter l'histoire du rayon (enfin rien de compliqué) et le truc de prendre que les points. Mais j'y serais jamais arrivé seul comme je connaissais pas toutes ces fonctions ;)

Merci beaucoup!

@micalement,
`ixM
10/04/2007, 14h58
IXaMaXI

Encore une dernière question :)

C'est un peu plus technique mais je me demandais comment matlab trouvait les points.

Selon moi, ce que find retourne (j'ai regardé) c'est une liste d'id qui correspondent à des points de la sphère qui correspondent eux-mêmes à une valeur de q et s. Ensuite on injecte cette référence dans u et v :)

C'est ça ? :bug:

Merci beaucoup pour ton aide ;)

@micalement,
`ixM
10/04/2007, 15h06
Jerome Briot
Citation:

Envoyé par IXaMaXI

C'est ça ?

On peut voir cela comme ça ;)

Voici un exemple détaillé de ce qui se passe :
Code:

1 2 3 4 a=rand(1,5) b=a>0.5 idx=find(b) a(idx)
ce qui peut directement s'écrire :
Code:

1 2 a=rand(1,5) a(a>0.5)

Re-salut :p

Bon, en fait on s'est rendu compte qu'il y avait une erreur dans la routine que tu nous avais fournie. Le problème venait d'abord du /(32.*pi.*cos(q)) qui faisait renvoyer des valeurs énormes à z (/0).

Ensuite on s'est rendu compte que la condition n'était pas bonne, l'intersection était celle d'un cylindre avec notre surface reglée.

Après on a réfléchi dans l'autre sens, on a essayé de trouver tous les points sur la surface reglée qui se trouvaient à approximativement 2pi du centre. En jouant sur la sensibilité et le nombre de pas, on est arrivé à avoir une belle courbe. Seulement c'est pas bon puisqu'on a pas vraiment récupéré les points logiquement sans raisonnement mathématique.

Donc on a recommencé :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 [x,y] = meshgrid(-7:0.25:7, -7:0.25:7); [q,s] = meshgrid(0:pi./32:pi, 0:pi./32:2.*pi); u = 2.*pi.*sin(q).*cos(s); v = 2.*pi.*sin(q).*sin(s); z = 2.*pi.*cos(q); cond = abs((u.^2).*v./16)-z; new_param = find(cond < eps); hold on; %surf(u, v, z); %surf(x, y, (x.^2).*y./16); plot3(u(new_param),v(new_param),(u(new_param).^2).*v(new_param)./16, '+k', 'Color', 'black', 'Linewidth', 1); hold off;
Là, on obtient la bonne forme mais on a de nouveau tous ces foutus points au milieu.

u = 2 * pi * sin(q) * sin(s);
v = 2 * pi * sin(q) * sin(s);
(u² * v)/16 = 2 * pi * cos(q);

C'est donc notre système pour trouver l'intersection. Les deux premiers ne posent pas de problème mais la dernière oui. On a donc repris l'idée de la condition et ça nous a donné le résultat décrit plus haut.

www.ixamaxi.be/images/matlab4.png

Et donc nous voilà de nouveau bloqués :/

Merci d'avance :aie:
@micalement,
`ixM

10/04/2007, 19h30
millie

Citation:

Selon moi, ce que find retourne (j'ai regardé) c'est une liste d'id qui correspondent à des points de la sphère qui correspondent eux-mêmes à une valeur de q et s.

Peut être est-ce pour ça que vous n'avez que des points. find retourne peut être un espace discret... (bon, je dis ça, mais j'y connais pas grand chose :aie:)

Salut :)

En fait c'est ce que nous chercons à avoir :p, rien que des points qui correspondraient à l'intersection de notre surface avec la sphère. Là pour l'instant, on a une sorte de surface qui vient de je sais pas trop où :/

A mon avis ça provient encore d'une erreur de paramétrisation mais cette fois-ci je m'imagine vraiment pas où :(

Voici toutefois la "mauvaise" méthode dont j'avais parlé tout-à-l'heure. On peut voir la surface reglée en décommentant une ligne.
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 [x,y] = meshgrid(-7:0.25:7, -7:0.25:7); [x1,y1] = meshgrid(-7:0.01:7, -7:0.01:7); [q,s] = meshgrid(0:pi./32:pi, 0:pi./32:2.*pi); u = x1; v = y1; z = (u.^2).*v./16; new_param = find(abs(sqrt(u.^2+v.^2+z.^2)-2.*pi) < 0.001); hold on; surf(2.*pi.*sin(q).*cos(s), 2.*pi.*sin(q).*sin(s), 2.*pi.*cos(q)); %surf(x, y, (x.^2).*y./16); plot3(u(new_param),v(new_param),z(new_param), '+k', 'linewidth', 1); hold off;
Alors, aussi entretemps, je viens de découvrir que tous les points que vous voyez représentés sur le screenshot de mon précédent post sont situés à 2pi du centre de la sphère au moment de la condition. Ce qui me pousse à confirmer l'hypothèse de la paramétrisation foireuse. Normalement, il ne devrait prendre que l'intersection et là il récupère des points qui ne peuvent pas appartenir à l'intersection. Ceux du centre du graphe par exemple, ne peuvent pas appartenir puisque (j'ai vérifié ça aussi) la sphère est vide (on le voit si on la trace en plot3 avec des croix) donc, il n'y a pas d'intersection entre la courbe et la boule.

J'espère que je ne déborde pas du cadre de ce forum en posant ce genre de question mais comme je suis assez peu expérimenté (ahah) avec matlab, j'arrive pas à dire si c'est une erreur de matlab ou une erreur de paramétrisation.

Si quelqu'un peut apporter une petite pierre à l'édifice de ma réflexion, je lui en serait très reconnaissant.

D'avance merci,
@micalement,
`ixM

P.S. Je ne demande pas de solution complète mais juste un élément de réponse :) Le reste je me débrouillerai.

11/04/2007, 08h08
LordPeterPan2
Citation:
Envoyé par IXaMaXI

Salut à tous,

J'utilise Matlab 7.0.0.19920 (R14) et j'essaie de calculer et représenter l'intersection entre une surface reglée et une sphère.

Ma surface a comme équations paramétriques :
x = u;
y = v;
z = (u² * v)/16;

Comme sphère j'ai choisi :
x = 2 * pi * sin(q) * sin(s);
y = 2 * pi * sin(q) * sin(s);
z = 2 * pi * cos(q);

J'ai donc dit que l'intersection devait être :
u = 2 * pi * sin(q) * sin(s);
v = 2 * pi * sin(q) * sin(s);
(u² * v)/16 = 2 * pi * cos(q);

Mais là déjà mathématiquement ça me paraît bizarre :/ M'enfin bon, je change ma dernière équation et j'obtiens :
(u² * v) / (32 * pi * cos(q)) = 1;

C'est bof...

Dans matlab j'ai mis ça :

Code:

1 2 3 4 5 6 7 8 [q,s] = meshgrid(0:pi./32:pi,0:pi./32:2.*pi); u=2.*pi.*sin(q).*cos(s); v=2.*pi.*sin(q).*sin(s); hold on; mesh(u, v, (u.^2).*v./(32.*cos(q))); hold off;

Donc on a notre sphère et notre surface reglée :
www.ixamaxi.be/images/matlab.png

Puis j'avais écrit ça :

Code:

1 2 3 4 5 6 7 8 [q,s] = meshgrid(0:pi./32:pi,0:pi./32:2.*pi); u=2.*pi.*sin(q).*cos(s); v=2.*pi.*sin(q).*sin(s); hold on; mesh(u, v, (u.^2).*v./16); hold off;

Ce qui me donnait :
www.ixamaxi.be/images/matlab1.png

Et ensuite avec le code donné plus haut :
www.ixamaxi.be/images/matlab2.png

Ce qui me perturbe c'est qu'on a une surface et pas une courbe. On s'attendait à voir une sorte de sinusoïde enroulée autour d'un cercle. Et c'est ce qu'on a dans le dernier screenshot mais on aimerait ne pas voir le cercle :D

Donc il doit y avoir soit une erreur dans notre routine matlab soit une erreur dans nos équations paramétriques.

J'imagine que les gens qui manipulent matlab depuis plus longtemps que 3 heures, doivent aussi savoir répondre au problème des équations paramétriques donc j'espère avoir une réponse :)

D'avance merci,
@micalement,
`ixM
Il suffit de résoudre à la main :
Code:

1 2 3 4 x = 2*pi*sin(q)*cos(s) y = 2*pi*sin(q)*sin(s) z = (x^2*y)/16 = 2 * pi * cos(q)
Non ???
11/04/2007, 10h46
IXaMaXI

Salut,

En fait on nous demande de trouver les équations paramétriques de l'intersection. Donc il n'y a pas vraiment grand chose à résoudre. On a déjà tenté de mettre dans l'équation z = x² * y / 16, x et y mais ça ne donne rien de bien.

C'est en fait après avoir écrit x² * y / 16 = 2 * pi * cos(q) que nous avons décidé d'écrire la dernière condition comme une restriction de points de la sphère grâce à find en écrivant

Code:

new_param = find((abs((u.^2).*v./16)-abs(z)) < eps);

Ce qui me semblait être la manière d'écrire dans matlab x² * y / 16 = 2 * pi * cos(q).

Est-ce c'est de là que pourrait venir l'erreur ? J'ai beau chercher, je ne trouve pas. J'ai passé douze (vraies) heures dessus hier avec mon groupe et on arrive pas à s'en sortir.

Notre problème se situe au niveau du fait que matlab trouve, avec ce qu'on a écrit, des intersections entre la surface et la sphère à des endroits où on voit bien qu'il n'y en a pas quand on n'affiche que la surface et la sphère.
Ca doit donc être un problème dans notre routine. :bug:

Merci :)
@micalement,
`ixM

Tiens essaye çà :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
 
Q = [] ;
S = [] ;
pas = 1000 ;
for q = -pi:pi/pas:pi
    X=bombelli(-1,4*cos(q)/(pi^2*(sin(q))^3));
    X=X(imag(X)==0) ;
    if (abs(X)<=1)
        for xx=X
            s = asin(xx)    ;
            Q = [q Q    q] ;
            S = [s S pi-s] ;
        end;
    end;
end;
 
[S,I] = sort(S);
Q=Q(I);
 
x=2*pi.*sin(Q).*cos(S);
y=2*pi.*sin(Q).*sin(S);
z=x.^2.*y/16;
 
[xx,yy,zz]=sphere(10);
xx=2*pi*xx;
yy=2*pi*yy;
zz=2*pi*zz;
 
figure;
surf(xx,yy,zz);
hold on;
plot3(x,y,z,'r+');

Avec :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
 
function [x] = bombelli(p,q)
% [x] = bombelli(p,q) solve the equation x^3+px+q=0
 
if (p==0)
    x=(-q)^(1/3) ;
else
    if (q==0)
        x = [-sqrt(-p) 0 sqrt(-p)];
    else
        DELTA = (27*q^2+4*p^3)/27 ;
        switch sign(DELTA)
            case  1
                temp = sqrt(DELTA)/2 ;
                x1 = power_lto(-q/2-temp,3)+power_lto(-q/2+temp,3);
                delta = x1^2+4*q/x1 ;
                x = [ (-x1-sqrt(delta))/2 x1 (-x1+sqrt(delta))/2 ] ;
            case  0
                x = sort([power_lto(-4*q,3) power_lto(q/2,3)]) ;
            case -1
                X1 = (-q-sqrt(DELTA))/2 ;
                u1 = power_lto(X1,3) ;
                u2 = exp(2*pi*i/3)*u1 ;
                u3 = exp(4*pi*i/3)*u1 ;
                v1 = -p/(3*u1) ;
                v2 = -p/(3*u2) ;
                v3 = -p/(3*u3) ;
                x = sort([(u1+v1) (u2+v2) (u3+v3)]);
        end;
    end;
end;
 
function z=power_lto(x,n)
 
if isreal(x)
    z = sign(x)*abs(x)^(1/n);
else
    z = x^(1/n);
end;

Ce code viens de la page web http://homeomath.imingo.net/equa31.htm
Et du raisonnement suivant :
L'intersection de ta sphère et de ta surface vérifie :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
 
x = 2*pi*sin(q)*cos(s)
y = 2*pi*sin(q)*sin(s)
z = (x^2*y)/16 = 2*pi*cos(q)
 
Donc :
 
8*pi^3*sin(q)^3*sin(s)*cos(s)^2/16 = 2*pi*cos(q)
 
Ce qui donne en simplifiant et en posant X = sin(s) :
X^3-X-4*cos(q)/(pi^2*sin(q)^3)=0
Equation du troisième ordre résolvable pour q donné d'après la méthode de bombelie.
Ensuite Si |X|<=1 s= acos(x) et s=pi-acos(x) conviennent

12/04/2007, 13h04
IXaMaXI

Salut,

Merci pour ta réponse. Ca nous a ouvert les yeux. On a laissé l'équation de la sphère en cartésien et la surface en paramétrique. On injecte et paf, voilà l'intersection :D

J'ai quand même une petite question. Comment est-ce qu'on fait pour utiliser ton code ? On peut coller la fonction avant l'autre code dans le même fichier ?

Merci beaucoup !
@micalement,
`ixM
12/04/2007, 17h24
LordPeterPan2

Moi j'ai fait deux fichiers :
un fichier script et un fichier pour la fonction bombelli qui contient à la fin la fonction power_lto
Les deux code que j'ai mis, correspondent aux deux fichiers que j'ai utilisé.