Discussion :

[jokfun]

Bonjour,

Après des recherches infructueuses sur le web j'ai décidé de poster un message ici. Je possède une matrice correspondant à mes coordonnées z. Ainsi lorsque je trace celle ci a l'aide de la fonction surf(z) j'obtiens bien ma surface. Maintenant je souhaiterais calculer le volume sous cette surface. La fonction dblquad me permettrait de le faire facilement mais je ne connais pas l'équation de ma fonction. Existe-t-il une manipulation spécifique pour utiliser dblquad directement avec ma matrice? Ou une autre fonction qui me permettrait d'effectuer ce calcul?

voici l'aperçu d'une surface obtenue :



Merci d'avance

jokfun

[phryte]

Salut.

Ou une autre fonction qui me permettrait d'effectuer ce calcul?

Peux-être une idée :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
clear
T=0.1;
N=30;
X=1:T:N;
Y=1:T:N;
Z=ones(length(X));;
surfc(X,Y,Z)
colormap hsv
axis([-30 30 -30 30 0 5])
c=0;
for n = X(1):T:X(30)
    c=c+cumtrapz(Y,Z(1,:));
end
c(end)

[mr_samurai]

Salut,

Je te propose cette méthode :

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
 
d1 = 0.2;
d2 = 0.1;
xx = -3:d1:3;
yy = -2:d2:2;
[X, Y] = meshgrid(xx, yy);
 
%MyFun = @(u,v)(exp(-(u.^2 + v.^2)));
MyFun = @(u,v)(u.^2 + v.^2);
 
Z = MyFun(X,Y) ;
 
%--- Approximative value
s = size(Z);
Ztmp = (Z(2:s(1),2:s(2)) + Z(1:s(1)-1,1:s(2)-1))/2;  % Moyenne des pavés
Iapprox = sum(Ztmp(:))*d1*d2; % Approximation en escaliers 
 
%--- Exact Value
 
Iexact = dblquad(MyFun,-3,3,-2,2);
 
%---- Compare
fprintf(' ## Approx = %f  \n ## Exact = %f  \n ## Error = %f \n',...
    Iapprox,Iexact, abs(Iapprox - Iexact)/Iexact)

Si ton maillage n'est pas trop grossier, les résultats devront collés .

++ bonne chance.

[casier]

Salut.

Peux-être une idée :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
clear
T=0.1;
N=30;
X=1:T:N;
Y=1:T:N;
Z=ones(length(X));;
surfc(X,Y,Z)
colormap hsv
axis([-30 30 -30 30 0 5])
c=0;
for n = X(1):T:X(30)
    c=c+cumtrapz(Y,Z(1,:));
end
c(end)

Tu utilises une intégration monodimensionnelle que tu répète pour chaque couche, ça ne peut pas faire une intégration bidimensionnelle, il va manquer un dy dans ton calcul si je me trompe pas non ?

@jokfun
Pourquoi ne pas bêtement intégrer ça en découpant le volume en cubes ? Tu découpes les axes x et y en N points chacun, et le volume que tu cherches peut être approché par la somme des delta_x*delta_y*Z(xi,yj) pour i=1..N et j=1..N tout simplement C'est codé en une ligne et ça marche très bien

Code :

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

1
2
3
4
5
6
7
8
9
dx=0.01;
dy=0.01;
X=0:dx:1;
Y=0:dy:1;
Z=ones(length(X),length(Y));
mesh(X,Y,Z) % voila une surface simple de volume 1
 
% calcul du volume
volume=dx*dy*sum(sum(Z(1:end-1,1:end-1)))

Par contre, sur ton image le maillage sur un des axes est assez grossier, le pic risque d'induire pas mal d'erreur car il n'est pas du tout discrétisé (il est constitué d'un seul point). Si tu peux, recalcule ta surface avec plus de points...

edit : je vois que mr samurai a posté une solution qui ressemble à celle-là mais avec des fonctions savantes

[phryte]

Salut.

il va manquer un dy dans ton calcul si je me trompe pas non ?

Pourtant mon calcul est bon !

[mr_samurai]

Salut,

Il manque des choses à ton code :

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
clear
T=0.1;
N=30;
X=1:T:N;
Y=1:T:N;
Z=ones(length(X));;
surfc(X,Y,Z)
colormap hsv
axis([-30 30 -30 30 0 5])
c=0;
for n = 1:length(X)
    c=c+cumtrapz(Y,Z(n,:));
end

c = c *T ; % sinon le pas suivant Y est égale à 1.
c(end)

[phryte]

Salut.

ll manque des choses à ton code

OK mr_samurai, le calcul est bon avec ta modif.
Merci

[mr_samurai]

edit : je vois que mr samurai a posté une solution qui ressemble à celle-là mais avec des fonctions savantes

Absolument pas . Je prend juste la moyenne des deux extrémités des pavés. Le reste du code est pour comparer pour des fonctions particulières...

Sinon, un vecteur de N points, contient N-1 segments :

Code :

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

1
2
3
4
5
6
7
dx = (X(2)-X(1));
dy = (Y(2)-Y(1));
 
% Ceci :
volume =  dx * y *sum(sum(Z(1:end-1,1:end-1)))
% Est une meilleur Approximation que : 
volume=  dx * dy * sum(sum(Z))

++

[casier]

mince, j'avais commencé à faire avec des length(X)-1 et length(Y)-1 dans Z, mais c'était un peu lourd alors j'ai décidé de faire avec Nx et Ny comme paramètres de départ et j'ai oublié de mettre end-1, pis j'avais aussi oublié de définir dx et dy du coup, ça m'apprendra Puisque c'est ça je remet l'ancienne version

[jokfun]

Merci beaucoup pour ces nombreuses réponses. Cependant voila je débute sur Matlab donc j'ai quelques soucis pour appliquer vos codes à ma situation...Je m'explique : ma matrice Z comporte 20 colonnes et 51 lignes ce qui ne donne pas un graph très précis mais je ne peux pas y toucher...voici le code de casier appliqué à mon cas :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
dx=1;
dy=1;
X=0:dx:19;
Y=0:dy:50;
mesh(X,Y,Z)
 
% calcul du volume
volume=dx*dy*sum(sum(Z(1:end-1,1:end-1)))
 
volume =
 
  7.4608e+003

Ainsi je suis obligé d'avoir un pas de 1 ce qui ne me semble pas très précis...est ce que je me trompe ou non?

J'ai aussi essayé la méthode de mr_samurai mais je n'ai pas réussit à la faire marcher avec mon exemple...je ne comprends pas vraiment toutes les lignes comme la fonction Myfun par exemple ou bien où rentre en compte ma matrice Z.

Edit : Ok donc je comprends un peu mieux ton code...cependant je ne vois vraiment pas comment je pourrais remplacer ta fonction MyFun u^2 + v^2 par ma matrice Z...

Enfin je me demande s'il n'est pas possible d'utiliser dblquad en posant une fonction Myfun qui traduirait ce que l'on voit sur le graphique?

[mr_samurai]

rebonjour,

tu peux poster un bout de ton code ?

++

[jokfun]

Justement c'est ce que j'explique dans mon post précédent je n'arrive pas à voir comment je pourrais réécrire ton code en l'appliquant à mon cas...je ne comprends pas l'utilisation de la fonction MyFun avec les puissances 2...elle ne correspond pas du tout à mon graphique...

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
d1 = 0.2;
d2 = 0.1;
xx = -3:d1:3;
yy = -2:d2:2;
[X, Y] = meshgrid(xx, yy);

%MyFun = @(u,v)(exp(-(u.^2 + v.^2)));
MyFun = @(u,v)(u.^2 + v.^2);
% => d'où vient cette fonction, pourquoi celle-ci?

Z = MyFun(X,Y) ;

%--- Approximative value
s = size(Z);
Ztmp = (Z(2:s(1),2:s(2)) + Z(1:s(1)-1,1:s(2)-1))/2;  % Moyenne des pavés
Iapprox = sum(Ztmp(:))*d1*d2; % Approximation en escaliers 

%--- Exact Value

Iexact = dblquad(MyFun,-3,3,-2,2);

%---- Compare
fprintf(' ## Approx = %f  \n ## Exact = %f  \n ## Error = %f \n',...
    Iapprox,Iexact, abs(Iapprox - Iexact)/Iexact)

[mr_samurai]

salut,

il n'y a que 3 lignes qui t'intéresse dans le code que j'ai posté, le reste est la pour comparer le calcul exact (DBQUAD) et l'approximation:

Code :

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

1
2
3
4
5
6
7
8
% - Variables à définir avant 
% d1 : le pas suivant l'axe des x
% d2 : le pas suivant l'axe des y 
% Z : la matrice de la surface.
 
s = size(Z);
Ztmp = (Z(2:s(1),2:s(2)) + Z(1:s(1)-1,1:s(2)-1))/2;  % Moyenne des pavés
I = sum(Ztmp(:))*d1*d2;   % I est le volume en dessous de la surface Z

++

[jokfun]

Bonjour,

et bien avec les dernières informations de mr_samurai cela fonctionne a merveille. Je tiens a préciser que le code de casier marche aussi mais il est nettement moins précis (je ne comprends pas pourquoi d'ailleurs...j'indique le même pas dans les deux cas...enfin bon ce n'est pas très important!)

Merci à tous pour votre aide vous m'enlevez une belle épine du pied!

jokfun

[mr_samurai]

Salut,

Les deux méthodes reposent sur le même principe : Approximer la fonction sur chaque rectangle de cotés dx et dy par une constante C. Le volume est donc C*dx*dy.

Pour les indices (i,j)
- Méthode 1 : (mon code)
C = [ Z(i,j) + Z(i+1,j+1) ] / 2

- Méthode 2 : (code de casier)
C = Z(i,j)

J'espère avoir était clair.

Résolu ?

++ bonne continuation

[jokfun]

Ok c'est noté je comprends mieux pourquoi ton calcul est plus précis maintenant (parce que mes points sont assez espacés en fait).

Oui problème résolu merci encore.