double intégration d'une matrice

Version imprimable

bonjour à tous,

je dois effecteur une double intégrale sur des valeurs stockées dans une matrice. la situation est la suivante :
ma matrice zi2 j'ai des valeurs de débit (une ligne correspond à une position angulaire et une colonne à une position radiale ). j'aimerai effectuer la double intégration selon theta (entre 0 et 2pi) puis entre 0 et rayon extérieur.

J'ai essayé la fonction quad mais ne disposant pas de la fonction je ne peux pas l'utiliser. Comment faire?

Voici mon code

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
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
clear all
 
 
%% file wird gelesen und bearbeitet
fichier='massenstrom_tutorium_gruppe2.txt';
fid=fopen(fichier,'r'); % oeffnet den file
X=fread(fid); %liesst den file
fclose(fid); %schliess den file
X=strrep(X,',','.'); %matlab erkennt nur die "." deswegen muss man alle kommas veraendern 
fichier2='massenstrom_tutorium_gruppe2_b.txt'; %neuen file korrigiert 
fid=fopen(fichier2,'w');
fwrite(fid,X); 
fclose(fid); 
fichier='massenstrom_tutorium_gruppe2_b.txt';
fid=fopen(fichier,'r');
C = textscan(fid,'%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s','headerLines',1);%den file wird gelesen und alle Schrifte werden gesprungen
fclose(fid); 
C=cell2mat( cellfun(@str2double,C,'UniformOutput',false) ); %die fruehere matrix besteht aus Zellen und wir mussen sie in eine numerische matrix verwandeln
%% von polar koordinaten--->kartesianische koordinaten
THETA=[-pi -pi -pi -pi -pi -pi -pi -pi -pi -pi -pi -pi -pi -pi -pi -pi -pi -pi -pi -pi -pi -pi -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -(3*pi)/4 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 -pi/4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/4 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 pi/2 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 (3*pi)/4 pi pi pi pi pi pi pi pi pi pi pi pi pi pi pi pi pi pi pi pi pi pi];
RADIUS=[0.220 0.2 0.18 0.16 0.140 0.130 0.120 0.110 0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.025 0.02 0.015 0.01 0.005 0 0.220 0.2 0.18 0.16 0.140 0.130 0.120 0.110 0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.025 0.02 0.015 0.01 0.005 0 0.220 0.2 0.18 0.16 0.140 0.130 0.120 0.110 0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.025 0.02 0.015 0.01 0.005 0 0.220 0.2 0.18 0.16 0.140 0.130 0.120 0.110 0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.025 0.02 0.015 0.01 0.005 0 0.220 0.2 0.18 0.16 0.140 0.130 0.120 0.110 0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.025 0.02 0.015 0.01 0.005 0 0.220 0.2 0.18 0.16 0.140 0.130 0.120 0.110 0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.025 0.02 0.015 0.01 0.005 0 0.220 0.2 0.18 0.16 0.140 0.130 0.120 0.110 0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.025 0.02 0.015 0.01 0.005 0 0.220 0.2 0.18 0.16 0.140 0.130 0.120 0.110 0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.025 0.02 0.015 0.01 0.005 0 0.220 0.2 0.18 0.16 0.140 0.130 0.120 0.110 0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.025 0.02 0.015 0.01 0.005 0];
%[x,y]=pol2cart(THETA,RADIUS);
%%von einen linien vektor zu einen Spalten vektor
%x=x';
%y=y';
%%rechnen von der geschwindigkeit
[l,m]=size(C);
i=1;
for a=1:l
    density=C(a,40)/(287*(C(a,42)+273.15));
    B(a,1)=density;
    for b=18:39
        geschwindigkeit(i,1)=sqrt((2*(C(a,b)-C(a,40)))/density)
        massenstrom(i,1)=density*geschwindigkeit(i,1);
        i=i+1;
    end
end
% Les bornes du domaine
rint = 0; rext = 0.220; 
[xext,yext] =pol2cart(linspace(0,2*pi,100),repmat(rext,1,100)); 
[xint,yint] = pol2cart(linspace(0,2*pi,100),repmat(rext-rint,1,100)); 
% Quelques valeurs "connues"
z = geschwindigkeit;
z2=massenstrom;
[x,y] = pol2cart(THETA,RADIUS); 
% figure(1)  % Affichage valeurs "connues" dans repere cartesien 
% subplot(221)
% plot(xext,yext,'b-',xint,yint,'b-')
% hold on 
% scatter(x,y,z,z)
% grid on 
% title('1 - Valeurs connues - Repere cartesien')
% colorbar
% axis equal tight  
% % Affichage valeurs "connues" dans repere polaire
% subplot(222)
% scatter(RADIUS,THETA,z,z)
% grid on 
% axis tight
% title('2 - Valeurs connues - Repere polaire') 
% Grille interpolation dans repere polaire
[R,T] = meshgrid(linspace(rint,rext,100),linspace(-pi,pi,100),'linerar');
% Interpolation dans repere polaire (domaine convexe) 
zi2 = griddata(RADIUS,THETA,z2,R,T);
zi = griddata(RADIUS,THETA,z,R,T);
% Affichage valeurs interpolees dans repere polaire
%subplot(223)
%scatter(R(:),T(:),zi(:),zi(:)) 
%grid on 
%axis tight
%title('3 - Valeurs interpolees - Repere polaire') 
% Affichage valeurs interpolees dans repere cartesien
[X,Y] = pol2cart(T,R); 
%subplot(224)
pcolor(X,Y,zi)
%surfc(X,Y,zeros(size(X)),zi,'edgecolor','none')
shading interp
colorbar view(2)
axis equal tight
title('4 - Valeurs interpolees - Repere cartesien')
[a,b]=size(zi2);
Massenstrom=quad(zi2,0,rext,0,2*pi)

merci beaucoup de l'aide

19/12/2011, 13h10
Jean Dumoncel

Bonjour,

si tu veux faire une intégration numérique basée sur un vecteur de ta fonction, tu peux utiliser trapz.
19/12/2011, 13h20
hittie
merci magelan mais quand j'essaie (pour effectuer la première intégrale (celle entre 0 et 2 pi)
Code:

1 2 3 4 5 6 [a,b]=size(zi2); j=1; for i=1:a Massenstrom(i,j)=trapz(zi2(i,j)) j=j+1 end
j'obtient une matrice nulle. Je ne vois pas comment je peux d'abord effectuer une intégration sur l'angle (une ligne) pour chaque point du rayon (une colonne) puis après intégrer l'ensemble sur le rayon total (l'ensemble des colonnes)

merci encore pour l'aide
19/12/2011, 13h37
Jean Dumoncel

C'est normal que tu obtiennes 0 : tu intègres unes seule valeur donc une constante. Tu peux intégrer les données d'une ligne en faisant trapz(zi2(i,:)) pour la ième ligne par exemple
19/12/2011, 14h05
FR119492
Salut!

Citation:

une double intégrale sur des valeurs stockées dans une matrice

Deux remarques:
1. "matrice" est un terme d'algèbre linéaire; dans ton cas, il s'agit d'un tableau et non d'une matrice.
2. On ne peut pas calculer l'intégrale d'un tableau, mais seulement d'une fonction définie dans un certain domaine, probablement un cercle. Tes données sont les valeurs échantillonnées de cette fonction.
La première étape que tu dois franchir est donc de choisir l'algorithme qui te permettra de calculer ce que tu cherches (la valeur de l'intégrale) à partir de ce que tu as (les valeurs stockées dans ton tableau). Ensuite, tu pourras commencer à programmer cet algorithme.
Jean-Marc Blanc
19/12/2011, 14h12
hittie
effectivement j'avais pas fait gaffe :oops: sinon j'ai une question concernant trapz: quand je fais
Code:

1 2 3 4 5 6 7 8 for k=1:a for i=1:a Massenstrom(c,d)=trapz(zi2(:,i)); c=c+1; end c=1; d=d+1; end
j'obtient une matrice (normal) mais j'aimerai après avoir une valeur globale (un nombre et non une matrice) .

Merci beaucoup à vous deux :)
19/12/2011, 14h19
Invité
Bonjour,

Ce que tu sembles ne pas avoir compris, c'est que
- Au début tu as un tableau 2D.
- Après intégration sur la première dimension, il te reste un vecteur 1D et non 2D comme tu sembles vouloir le calculer avec Massenstrom(c,d)
- Après intégration sur la deuxième, il te reste un scalaire
19/12/2011, 15h02
hittie

j'ai toujours un drole de résultat à savoir que ma matrice zi2 contient des valeurs comprises entre 30.3 et 50.3 mais quand je calcule ma variable MAssenstrom avec trapz j'obtient des valeurs de l'ordre de 10^5 et je ne comprend pas d'où cela peut venir...

encore merci

Un exemple:
Code:

1 2 3 4 x = 0:0.01:1; y = x.^2; trapz(x,y) trapz(y)

20/12/2011, 00h59
hittie

mais la question que je me pose c'est pourquoi en intégrant sur des valeurs qui sont inférieures à 50 j'obtient du 10^3

sinon merci beaucoup pour les explications :)
20/12/2011, 16h01
FR119492

Salut!

Citation:

pourquoi en intégrant sur des valeurs qui sont inférieures à 50 j'obtient du 10^3

C'est tout à fait normal: par exemple, si ta fonction est une constante f(x)=50, son intégrale de 0 à 20 vaut 1000.
Jean-Marc Blanc
20/12/2011, 16h30
hittie
désolé de reposer la question mais j'arrive pas à voir comme intégrer sur la colonne :
Code:

1 2 3 4 [a b]=size (zi2) for i=1:a I=trapz(zi2(i,:) end
j'obtient bien un vecteur colonne car il intègre ligne par ligne mais après quand je fais

Code:

Massenstrom=trapz(I)

j'obtient tjs 0 alors que normalement selon l'aide de Matlab

Code:

If Y is a vector, trapz(Y) is the integral of Y.

donc pour un vecteur non nul je devrais avoir une valeur différente de 0

J'ai testé avec le vecteur V=[1 2 3; 4 5 6;6 7 8]

et j'obtient 0...
comment je peux intégrer sur la colonne?

merci beaucoup
20/12/2011, 16h56
Jean Dumoncel
Citation:

Envoyé par hittie

j'obtient bien un vecteur colonne car il intègre ligne par ligne

Si tu avais vérifié tu aurais vu que I est un scalaire car tu écrases la valeur à chaque itération, tu obtiens uniquement l'intégrale de la dernière ligne.
Il aurait fallu écrire :
Code:

1 2 3 4 [a b]=size (zi2) for i=1:a I(i)=trapz(zi2(i,:) end
ce qui peut s'écrire plus simplement :

Code:

trapz(zi2,2)
20/12/2011, 17h05
hittie

merci magelan, Winjerome et FR119492 pour votre aide précieuse :)