Using trapz just for positive areas

조회 수: 26(최근 30일)
Simone Speltoni
Simone Speltoni 2021년 6월 9일
댓글: Star Strider 2021년 6월 9일
Hi, I already tried to solve my problem by looking at this answer https://it.mathworks.com/matlabcentral/answers/254940-how-to-find-the-area-between-the-x-axis-and-positive-section-of-a-graph, but without success, I guess because I'm missing how to link y-data to x-data. My goal consists of getting just the portion of area under a curve (which is not a function) and above the x-axis, then a positive area. For example:
clear all
close all
xA=0:4;
A=[-2 -1 2 3 4]; %if one considers this as a polygonal chain,
%positive area is equal to 6.66, negative area -1.665
plot(xA,A,'o');
grid on
hold on
Valpos = find(A>=0);
posArea = trapz(xA(Valpos),A(Valpos)); %(positive area) =6
integral=trapz(A); %(positive area)-(negative area)=5
%polyfit to build a curve with "A" data
[pA,~,muA] = polyfit(xA, A, 4);
fA = polyval(pA,xA,[],muA);
plot(xA,fA);
trapzinterpolaz=trapz(fA); %(positive area)-(negative area)=5
I understand how trapz works, and it's right to get these results (6 as positive area and 5 if one considers negative trapezoids too). Anyway, what I would like to have is considering data in A as a function, because in this case I would get 6.66 as positive area, and -1.665 as negative area. Basically, I'm missing the positive area between x=1 and x=2 IF the A data were a function. Any advice? thanks!

채택된 답변

Star Strider
Star Strider 2021년 6월 9일
I would just use interp1 to interpolate the data, then use trapz to integrate it. The resolution of ‘xi’ determines the precision of the integrated result. (I use the default length of 100 here.)
xA=0:4;
A=[-2 -1 2 3 4]; %if one considers this as a polygonal chain,
%positive area is equal to 6.66, negative area -1.665
plot(xA,A,'o');
grid on
hold on
Valpos = find(A>=0);
posArea = trapz(xA(Valpos),A(Valpos)); %(positive area) =6
integral=trapz(A); %(positive area)-(negative area)=5
%polyfit to build a curve with "A" data
[pA,~,muA] = polyfit(xA, A, 4);
fA = polyval(pA,xA,[],muA);
plot(xA,fA);
trapzinterpolaz=trapz(fA); %(positive area)-(negative area)=5
xi = linspace(min(xA), max(xA)); % Interpolation Vector
Ai = interp1(xA,A,xi,'linear'); % Interpolated Data
posIdx = Ai>0; % Logical Index Of Positive 'A' Values
areaPos = trapz(xi(posIdx), Ai(posIdx)) % Positive Area
areaPos = 6.6638
areaNeg = trapz(xi(~posIdx), Ai(~posIdx)) % Negative ARea
areaNeg = -1.6664
.
  댓글 수: 2
Star Strider
Star Strider 2021년 6월 9일
As always, my pleasure!

댓글을 달려면 로그인하십시오.

추가 답변(1개)

Joel Lynch
Joel Lynch 2021년 6월 9일
You need to use more fine-grid X data in polyval to get test points at y=0.
xA_fine = linspace(xA(1),xA(end),100);
fA_fine = polyval(pA,xA_fine,[],muA);
% Full Integral
trapz( xA_fine, fA_fine )
% Just positive Values (logical indicies)
idx = fA_fine>=0;
trapz( xA_fine(idx), fA_fine(idx) )
% Just positive Values (using a starting index)
i = find(fA_fine>=0,1)
trapz( xA_fine(i:end), fA_fine(i:end) )
Also note: multiple ways to get values greater than zero. The first is to use logical indices, the second is to find the first place and integrate from there. The second only works when the function is always increasing. The first may run into problems if there are multiple y=0 cross-overs, depending on how xA_fine is defined.
  댓글 수: 1
Simone Speltoni
Simone Speltoni 2021년 6월 9일
Thank you Joel, this approach is very useful too!

댓글을 달려면 로그인하십시오.

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by