How to calculate area between curve and horizontal line?

조회 수: 17 (최근 30일)
Hamzah Mahmood
Hamzah Mahmood 2020년 8월 21일
댓글: esat gulhan 2020년 8월 23일
Apologies if this has already been answered in another form elsewhere on the forums, however I couldn't find any answers that helped me directly with my code. I currently have code that maps the profile and area of a channel. I wish to have a horizontal y value across my channel width to simulate the fluctuating water level, the example line is shown in the image as the green line. So for example I would want to calculate and shade the area underneath the greenline only.
With the function would I be able to loop the same function for a range of values for the horizontal values, to provide an output I could then plot? A vector form would be ideal if possible.
I've already tried to use the area and fill commands unsuccesfully, or maybe I was not using the correct format for them.
I've also attached the relevant code below, any help would be much appreciated.
x = [ 0 1 2 3 4 5 6 7 8 9 10 11 12];
y = -[0 -0.5 -0.8 -0.8 -1 -1.1 -1.2 -1.2 -1.4 -1.2 -1.1 -1 0]; % y values in image were produced by rand function
subplot(2, 1, 1);
area(x,-y)
grid on;
hold on;
for k = 1 : length(x)
xline(x(k), 'Color', 'k');
end
yt = yticks;
for k = 1 : length(yt)
yline(yt(k), 'Color', 'k');
end
plot (x, -y, 'r.-', 'LineWidth', 3, 'MarkerSize', 30)
hold on
yline(-0.5,'g',11);
xlabel('Chainage (m)', 'FontSize', 11);
ylabel('Depth (m)', 'FontSize', 11);
An = (trapz(x,y)); % Overall area.
% Compute areas within each pair of x locations.
midpoints = (y(1:end-1) + y(2:end))/2;
deltax = diff(x);
Ani = deltax .* midpoints;
subplot(2, 1, 2);
bar(x(1:end-1)-x(1)+ 0.5, Ani);
title('Areas of Slices', 'FontSize', 14);
xlabel('x', 'FontSize', 11);
ylabel('Area (m^2)', 'FontSize', 11);
grid on;

채택된 답변

esat gulhan
esat gulhan 2020년 8월 21일
clc;clear;
x = [ 0 1 2 3 4 5 6 7 8 9 10 11 12];
y = [0 -0.5 -0.8 -0.8 -1 -1.1 -1.2 -1.2 -1.4 -1.2 -1.1 -1 0];
h0= -0.5; %the axis in the graph
h=ones(size(x))*h0;
Int=pchip(x,y);
xx=linspace(x(1),x(end),500);
yy=y-h0;
yyy=pchip(x,yy,xx);
gt=yyy<0;
theArea =-trapz(xx(gt), (yyy(gt)))
hold off
area(xx, min(ppval(Int,xx), h(1)), h(1), 'EdgeColor', 'none', 'FaceColor', 'b'),grid,
hold on
plot(xx,ppval(Int,xx),'k',xx,h0,'Linewidth',2)
set(gca, 'XLim', [x(1) x(end)], 'YLim', [min(y) 0]);
title(strcat('BlueArea=', num2str(theArea)))
xlabel('Chainage (m)', 'FontSize', 11);
ylabel('Depth (m)', 'FontSize', 11);
You can change h0 which is horizontal line. when h0=-1
when h0=-0.5
when h0=0
I RESHAPED YOUR X AND Y DATAS 500 INTERVALS WITH PCHIP INTERPOLATION, IF NOT THE DATA IS NOT ENOUGH FOR GOOD SOLUTION.
  댓글 수: 4
Hamzah Mahmood
Hamzah Mahmood 2020년 8월 22일
I'll look into that, cheers.
Another thing has come up, how would I be able to then keep these values as an array, im trying this however it isn't working.
clc;
clear;
x = [ 0 1 2 3 4 5 6 7 8 9 10 11 12];
y = [0 -0.5 -0.8 -0.8 -1 -1.1 -1.2 -1.2 -1.4 -1.2 -1.1 -1 0];
for stage= [-2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0] %the axis in the graph
stage_dim= length(h)
h=ones(size(x))*stage;
Int=pchip(x,y);
xx=linspace(x(1),x(end),500);
yy=y-stage;
yyy=pchip(x,yy,xx);
gt=yyy<0;
Area =zeros(stage_dim);
for i = 1:stage_dim
Area(i) =-trapz(xx(gt), (yyy(gt)))
end
hold off
area(xx, min(ppval(Int,xx), h(1)), h(1), 'EdgeColor', 'none', 'FaceColor', 'b'),grid,
hold on
plot(xx,ppval(Int,xx),'k',xx,stage,'Linewidth',2)
set(gca, 'XLim', [x(1) x(end)], 'YLim', [min(y) 0]);
title(strcat('Area=', num2str(Area)))
xlabel('Chainage (m)', 'FontSize', 11);
ylabel('Depth (m)', 'FontSize', 11);
end
esat gulhan
esat gulhan 2020년 8월 23일
Sorry, i dont understand your question. If you ask this question in new stage, other users can see and answer...

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

추가 답변 (2개)

esat gulhan
esat gulhan 2020년 8월 21일
편집: esat gulhan 2020년 8월 21일
x = [ 0 1 2 3 4 5 6 7 8 9 10 11 12]
y = -[0 -0.5 -0.8 -0.8 -1 -1.1 -1.2 -1.2 -1.4 -1.2 -1.1 -1 0]
Int=pchip(x,y)
plot(x,ppval(Int,x))
MIItemp = (fnint(Int))
Area = fnval(MIItemp,12)%Area between 0 12
If you want to find area between 0-5 , Area = fnval(MIItemp,5)
There are other ways but this is the easiest way .
Note:If your matlab is latetest version you can write makima(x,y),instead of phcip(x,y). If your lines are strict use makima
  댓글 수: 1
Hamzah Mahmood
Hamzah Mahmood 2020년 8월 21일
Hi Esat,
I think you may have misinterpered which axes I'm looking to cacluate the area from. The -0.5 value is a line that shifts up and down the y-axis rather than along the x-axis like so (the shade blue area).
Would your code allow me to calculate the area n the verticle directon if moodified slightly( I tried for a box with y = ones(1,13) however this didnt work). I'm relatively new to Matlab, so your code has helped introduce me to some new commands I hadn't come across before too.
Thanks.

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


esat gulhan
esat gulhan 2020년 8월 22일
Did you see my second answer. I think it helps.
If you add pause to the code it will seem like an animation.
clc;clear;
x = [ 0 1 2 3 4 5 6 7 8 9 10 11 12];
y = [0 -0.5 -0.8 -0.8 -1 -1.1 -1.2 -1.2 -1.4 -1.2 -1.1 -1 0];
for h0= 0:-0.04:-2; %the axis in the graph
h=ones(size(x))*h0;Int=pchip(x,y);xx=linspace(x(1),x(end),500);
yy=y-h0;yyy=pchip(x,yy,xx);gt=yyy<0;theArea =-trapz(xx(gt), (yyy(gt)))
hold off
area(xx, min(ppval(Int,xx), h(1)), h(1), 'EdgeColor', 'none', 'FaceColor', 'b'),grid,
hold on
plot(xx,ppval(Int,xx),'k',xx,h0,'Linewidth',2)
set(gca, 'XLim', [x(1) x(end)], 'YLim', [min(y) 0]);
title(strcat('BlueArea=', num2str(theArea)))
xlabel('Chainage (m)', 'FontSize', 11);ylabel('Depth (m)', 'FontSize', 11);pause(0.00001)
end

카테고리

Help CenterFile Exchange에서 Visual Exploration에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by