Calculer un volume à partir d'un graphique

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 :

http://nsa02.casimages.com/img/2008/...3729279176.jpg

Merci d'avance

jokfun

Salut.

Citation:

---

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

---

Peux-être une idée :

Code:

---

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)

---

Salut,

Je te propose cette méthode :

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

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.

Citation:

---

Envoyé par phryte

Salut.

Peux-être une idée :

Code:

---

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:

---

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

Salut.

Citation:

---

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

---

Pourtant mon calcul est bon !

Salut,

Il manque des choses à ton code :

Code:

---

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)

---

Salut.

Citation:

---

ll manque des choses à ton code

---

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

Citation:

---

Envoyé par casier

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

---

Absolument pas :aie: . 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:

---

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

---

++

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 :aie: Puisque c'est ça je remet l'ancienne version :yaisse2:

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:

---

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?

rebonjour,

tu peux poster un bout de ton code ?

++

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:

---

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)

---

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:

---

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

---

++

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

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

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.