Discussion :

[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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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

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

[Jerome Briot]

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

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[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

[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 ?

Merci beaucoup pour ton aide

@micalement,
`ixM

[Jerome Briot]

C'est ça ?

On peut voir cela comme ça

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

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
a=rand(1,5)
b=a>0.5
idx=find(b)
a(idx)

ce qui peut directement s'écrire :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[IXaMaXI]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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
@micalement,
`ixM

[millie]

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 )

[IXaMaXI]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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.

[LordPeterPan2]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 ???

[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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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.

Merci
@micalement,
`ixM

[LordPeterPan2]

Tiens essaye çà :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[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

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

[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é.