area under and over the curve

19 ビュー (過去 30 日間)
jackie
jackie 2018 年 8 月 12 日
コメント済み: Angelavtc 2020 年 4 月 28 日
hello I hope you can help me!
as per diagram below, I have a curve (blue) and a straight line (pink) say y=1.
I will need to find the unit area under the curve (yellow area) and "over" the curve (green area).
also the curve dont have many data points, I hope this wouldn't be a huge problem!
Please advise and many thanks in advance!

採用された回答

Image Analyst
Image Analyst 2018 年 8 月 12 日
Try this:
xLast = find(y > 1, 1, 'last');
areaUnderCurve = sum(abs(y(1:xLast) - 1))
This is basically assuming the digitized y is a bar chart and you're summing the heights of the rectangular bars. Of use trapz() if you want trapezoids instead of rectangular bars.
  2 件のコメント
jackie
jackie 2018 年 8 月 12 日
編集済み: jackie 2018 年 8 月 12 日
hey Image Analyst, thank you for the prompt reply!
if possible , Can you quickly explain the 2 line of code please if possible?
and also how to use trapz()?
thanks!
Image Analyst
Image Analyst 2018 年 8 月 12 日
y starts below 1, then rises above 1, then finally drops below 1 again. It's this x coordinate, where it drops below 1 that we need to find because that defines the rightmost place where you want to integrate.
If you do y-1, that will give you the area when y > 1, but that won't work on the left because it will give you negative areas. You'd want 1-y there. But 1-y is simply abs(y-1) so that now allows us to use the same formula everywhere because that will also work for the right area because abs(a positive number) is equal to the positive number.
trapz() is explained in the documentation. Do you really want me to just repeat what's already in the documentation? Basically instead of assuming the areas of each point are rectangular bars, it draws a line from one data point to the adjacent one and this gives you bars with sloping tops rather than flat tops. It's just a different definition and way of doing it. Use whichever you want.

サインインしてコメントする。

その他の回答 (3 件)

jackie
jackie 2018 年 8 月 13 日
編集済み: jackie 2018 年 8 月 13 日
ok I figured it out, thanks for you help!
x=[0,5,10,20,50,100,150,250,350,450];
b=[4.01E-01,8.30E-01,9.49E-01,1.08E+00,1.21E+00,1.23E+00
,1.20E+00,1.12E+00,1.11E+00,9.69E-01];
line=[1,1,1,1,1,1,1,1,1,1];
x_fine=0:10:450
b_fine = interp1(x,b,x_fine,'spline');
figure
plot(x,b,'o',x_fine,b_fine,'r', x,line,'b');
xFirst = find(b_fine > 1, 1, 'first');
xLast = find(b_fine > 1, 1, 'last');
area_total = trapz(x_fine,b_fine-1)
area_Loss = trapz(x_fine(1:xFirst),b_fine(1:xFirst)-1)
area_Gain = trapz(x_fine(xFirst:xLast),b_fine(xFirst:xLast)-1)
  1 件のコメント
Angelavtc
Angelavtc 2020 年 4 月 27 日
But using this formula, you are not caculating the yellow area but also all the area that is below y=1

サインインしてコメントする。


Image Analyst
Image Analyst 2020 年 4 月 27 日
This is how I'd do it:
%--------------------------------------------------------------------------------------------------------------------------------------------------------
% CLEAN UP - INITIALIZATION STEPS
clc; % Clear the command window.
fprintf('Beginning to run %s.m.\n', mfilename);
close all; % Close all figures (except those of imtool.)
clearvars; % Erase all existing variables.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
%--------------------------------------------------------------------------------------------------------------------------------------------------------
% TRAINING DATA PREPARATION
x = [0,5,10,20,50,100,150,250,350,450];
b = [4.01E-01,8.30E-01,9.49E-01,1.08E+00,1.21E+00,1.23E+00,...
1.20E+00,1.12E+00,1.11E+00,9.69E-01];
line = ones(1, length(x));
x_fine = linspace(min(x), max(x), 1000); % a thouisand points.
b_fine = interp1(x,b,x_fine,'spline');
line_fine = ones(1, length(x_fine));
hFig = figure
plot(x, b, 'bo', 'LineWidth', 2);
hold on;
plot(x, line, 'b', 'LineWidth', 2);
plot(x_fine, b_fine, 'r-', 'LineWidth', 2);
grid on;
xlabel('X', 'FontSize', fontSize);
ylabel('b', 'FontSize', fontSize);
title('trapz() demo', 'FontSize', fontSize);
indexLeft = find(b_fine > 1, 1, 'first');
indexRight = find(b_fine > 1, 1, 'last');
% Put up vertical lines there.
xline(x_fine(indexLeft), 'Color', 'm', 'LineWidth', 2);
xline(x_fine(indexRight), 'Color', 'm', 'LineWidth', 2);
%--------------------------------------------------------------------------------------------------------------------------------------------------------
% COMPUTE AREAS BETWEEN CURVE AND FLAT BLUE LINE
% Find that area below the line at 1.
y = line_fine(1 : indexLeft - 1) - b_fine(1 : indexLeft - 1);
area_Loss = trapz(x_fine(1 : indexLeft - 1), y)
% To double check, compute area by summing also.
area_Loss2 = sum(y) * abs(x_fine(2) - x_fine(1)) % Should be fairly close.
str = sprintf('Area below blue line = %.3f', area_Loss);
text(20, 0.86, str, 'Color', [0, 0.5, 0], 'FontSize', fontSize);
% Find that area above the line at 1.
y = b_fine(indexLeft : indexRight) - line_fine(indexLeft : indexRight);
area_Gain = trapz(x_fine(indexLeft : indexRight), b_fine(indexLeft : indexRight)-1)
% To double check, compute area by summing also.
area_Gain2 = sum(y) * abs(x_fine(2) - x_fine(1)) % Should be fairly close.
str = sprintf('Area above blue line = %.3f', area_Gain);
text(140, 1.05, str, 'Color', [0, 0.5, 0], 'FontSize', fontSize);
area_total = area_Loss + area_Gain
str = sprintf('Sum of both areas = %.3f', area_total);
text(220, 0.75, str, 'Color', [0, 0.5, 0], 'FontSize', fontSize);
hFig.WindowState = 'maximized';
fprintf('Done running %s.m.\n', mfilename);

Angelavtc
Angelavtc 2020 年 4 月 28 日
This is wonderful @Image Analyst, with this, you helped me to answer a question that I posted here: https://www.mathworks.com/matlabcentral/answers/521056-how-to-find-the-are-between-two-graphs?s_tid=mlc_ans_email_view#comment_836228
if you want you can post this answer and I will accept it.
Thank you!
  1 件のコメント
Angelavtc
Angelavtc 2020 年 4 月 28 日
So basically, the point was to calculate the difference between both areas (blue line vs red spline curve)

サインインしてコメントする。

Community Treasure Hunt

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

Start Hunting!

Translated by