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