[trapz] résultats abérrants

**hittie** · 24/12/2011, 18h59

Bonsoir à tous,

je dois résoudre une intégrale :

mais je n'arrive pas à obtenir de résultat cohérent. Mon code est le suivant

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
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';
THETA=THETA';
RADIUS=RADIUS';
%%rechnen von der geschwindigkeit
[l,m]=size(C);
i=1;
for a=1:l
    density=C(a,40)/(287*(C(a,42)+273.15));
    Density(a,1)=density;
    B(a,1)=density;
    for b=39:-1:18
        geschwindigkeit(i,1)=sqrt((2*(C(a,b)-C(a,40)))/density)
        %massenstrom(i,1)=density*geschwindigkeit(i,1);
        i=i+1;
    end
end
density=mean(Density(:,1));
% 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;
[x,y] = pol2cart(THETA,RADIUS); 
% figure(1)  % Affichage valeurs "connues" dans repere cartesien
% Grille interpolation dans repere polaire
[R,T] = meshgrid(linspace(rint,rext,300),linspace(-pi,pi,100));
 
% Interpolation dans repere polaire (domaine convexe) 
zi = griddata(RADIUS,THETA,z,R,T,'cubic');
A=density.*zi(:,:);
for i=1:300
    B(i,1)=trapz(A(:,i));
end
for i=1:100
    massenstrom=trapz(B(:,1)*2*pi);
end
 
[X,Y] = pol2cart(T,R); 
 
figure(1)
contourf(X,Y,zi,40)
% surf(X,Y,zeros(size(X)),zi,'edgecolor','none')
colorbar 
view(2)
axis equal tight

La variable A contient le rho*c (matrice 100*300) pour chacuns des points (colonne=position radiale croissante et ligne=position angulaire). Comment je peux calculer la double intégrale sachant que j'ai essayé trapz mais j'obtient un résultat complètement abérrant?

merci beaucoup pour votre aide

**Jean Dumoncel** · 25/12/2011, 22h35

Bonjour,

la boucle for :

Code :

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

1
2
3
for i=1:100
    massenstrom=trapz(B(:,1)*2*pi);
end

est inutile : tu fais 100 fois le même calcul...

Et je pense que tu devrais utiliser la forme Z = trapz(X,Y) sinon trapz ne sait pas comment est définie ta fonction. Ce qui donnerait :

Code :

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

1
2
B = trapz(THETA,A);
massenstrom = trapz(RADIUS,B*2*pi); % pourquoi multiplies-tu par 2pi?

**hittie** · 28/12/2011, 11h32

merci magelan mais dès que j'essaie cette méthode j'obtient :

Code :

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

1
2
3
4
5
6
??? Error using ==> trapz at 59
LENGTH(X) must equal the length of the first non-singleton dimension of
Y.
 
Error in ==> calculus at 56
B=trapz(THETA,A);

ma matrice A est une 100x300 et THETA (comme RADIUS) une 198x1. Ce qui est bizarre c'est que ma matrice B est une 9x1.

**Jean Dumoncel** · 28/12/2011, 12h25

Oui, je me suis trompé dans les noms des variables, mais tu aurais pu rectifier tout seul :

Code :

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

1
2
B = trapz(T(:,1),A);
massenstrom = trapz(R(1,:),B*2*pi); % pourquoi multiplies-tu par 2pi?

**hittie** · 28/12/2011, 17h55

j'ai déjà testé mais j'obtient toujours la même erreur...Je sais que je dois intégrer la matrice A d'abord sur les colonnes pour avoir une valeur pour un rayon donné sur un tour complet (une colonne )puis après intégrer sur les colonnes pour avoir l'intégrale sur le rayon mais je ne ois pas comment faire. A chaque fois j'ai comme erreur que j'ai la dimension qui coince.

Invité · 28/12/2011, 18h18

Bonjour,

Aurais-tu une nouvelle version du code?
Car avec ces seules lignes:

Code :

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

1
2
3
4
5
6
7
A = rand(100,300);  % Génération d'une matrice 100x300 pour le test
% Les bornes du domaine
rint =0; rext = 0.220; 
% Grille interpolation dans repère polaire
[R,T] = meshgrid(linspace(rint,rext,300),linspace(-pi,pi,100));
B = trapz(T(:,1),A);
massenstrom = trapz(R(1,:),B*2*pi);

Je n'obtiens aucune erreur

**hittie** · 28/12/2011, 20h41

en fait j'ai un résultat qui sort qui est assez aberrant car j'ai pas l'impression qu'il fasse bien la double intégrale .

Invité · 28/12/2011, 20h56

Tu n'as donc plus de message d'erreur, mais un résultat non attendu

Donc encore une fois, nous montrer le nouveau code du calcul de ton intégrale double pourrait aider.

**hittie** · 28/12/2011, 21h03

voici le code

la matrice zi contient les infos sur la vitesse (ligne =position angulaire et colonne =position radiale). J'ai donc 100 angles et 300 positions radiales au total

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
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';
THETA=THETA';
RADIUS=RADIUS';
%%rechnen von der geschwindigkeit
[l,m]=size(C);
i=1;
for a=1:l
    density=C(a,40)/(287*(C(a,42)+273.15));
    Density(a,1)=density;
    B(a,1)=density;
    for b=39:-1:18
        geschwindigkeit(i,1)=sqrt((2*(C(a,b)-C(a,40)))/density)
        %massenstrom(i,1)=density*geschwindigkeit(i,1);
        i=i+1;
    end
end
density=mean(Density(:,1));
% 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;
[x,y] = pol2cart(THETA,RADIUS); 
% figure(1)  % Affichage valeurs "connues" dans repere cartesien
% Grille interpolation dans repere polaire
[R,T] = meshgrid(linspace(rint,rext,300),linspace(-pi,pi,100));
 
% Interpolation dans repere polaire (domaine convexe) 
zi = griddata(RADIUS,THETA,z,R,T,'cubic');
A=density.*zi(:,:);
B=trapz(T(:,1),A);
massenstrom=trapz(R(1,:),B);
[X,Y] = pol2cart(T,R); 
figure(1)
contourf(X,Y,zi,40)
% surf(X,Y,zeros(size(X)),zi,'edgecolor','none')
colorbar 
view(2)
axis equal tight

merci beaucoup pour l'aide

**Jean Dumoncel** · 29/12/2011, 12h54

Qu'appelles-tu un résultat aberrant?

As-tu vérifié les valeurs contenu dans A? Elles te paraissent cohérentes?

Le code avec trapz permet de faire la double intégrale sur les valeurs numériques contenues dans A. Tu obtiens donc une approximation, mais l'erreur ne devrait pas être très grande. Donc si le résultat est très éloigné de ce que tu attendais, le problème se situe peut-être dans le calcul de cette matrice A.

**hittie** · 29/12/2011, 13h23

comme résultat aberrant j'obtient un débit massique gigantesque.
ma matrice A est correcte car j'ai des valeurs correctes pour zi. En fait le débit que je dois calculer correspond à rho*c*Aire mais comme j'ai la vitesse qui varie en chaque point(la variable c) je dois intégrer sur l'ensemble de la surface connaissant rho (valeur moyenne). J'ai vérifié les valeurs de zi (la matrice qui contient les vitesses en chaque point) et elle est correcte.

[trapz] résultats abérrants

MATLAB

Discussions similaires

Partager

Partager