How can I calculate the area between two curves?

Hello everyone. I am trying to calculate the area between two curves. How can I do this having only the x and y coordinates for each curve? Thanks in advance.

[Muhammad Usman]

Please Let me know if this is of any help...

area = 0;
for i=1:size(x1,2)-1
wid = ( sqrt((x1(i+1) - x1(i) )^2 + (y1(i+1) - y1(i) )^2)  + sqrt((x2(i+1) - x2(i))^2 + (y2(i+1) - y2(i))^2) )/2;
height = ( sqrt((y2(i) - y1(i))^2 + (x2(i) - x1(i))^2)  +  sqrt( ( y2(i+1) - y1(i+1) )^2 + (x2(i+1) - x1(i+1))^2 ))/2
area(i) = wid*height;

end
sum(area)

Answer 1

[Star Strider]

Use the trapz (or cumtrapz if their x-coordinates are the same and you want to subtact them element-wise) function to integrate each one, and then subtract the integral of one (calculated by trapz) from the other.

[Star Strider]

It depends on why your data are not strictily monotonic. If you have duplicate consecutive values in your independent variable, I add increasing values of eps (for example eps, 2*eps ...) to the subsequent duplicates.

If they are random but plot as continuous points, just sort them (using sortrows) by the independent variable.

[Paschalis Garouniatis]

I tried this code and got me the same error. If i wrote it right then interpolation should be done between x(1,1) and x(2,1). As for the yq, x should be my x-data and y should be my y-data right? I am still trying to figure out what is wrong. Any ideas?

 [x0,y0]=intersections(d1{1,10}(:,1),d1{1,10}(:,2),d2{1,10}(:,1),d2{1,10}(:,2));
sortrows(d1{1,10})
x=d1{1,10}(:,1)
y=d1{1,10}(:,2)
xq=x(1,1):0.00005:x(2,1)
yq=interp1(x,y,xq)

[Star Strider]

I have no idea how to interpret that code snippet.

[Star Strider]

Yours are not easy data to analyse. In the end, I solved the interpolation problem by first getting unique rows, then deleting rows with duplicate independent variables (this probably could have been done in one step rather than two), then doing the interpolation. The plotting and integral calculations were straightforward after that.

The code:

C1 = load('Paschalis Garouniatis coordinates1.mat');
C2 = load('Paschalis Garouniatis coordinates2.mat');
A1 = C1.coordinates;
A2 = C2.alpha;

Au1 = unique(A1, 'rows');                                           % Unique Rows (Remove Duplicates)
Au2 = unique(A2, 'rows');                                           % Unique Rows (Remove Duplicates)

D1 = diff([0; Au1(:,1)]);                                           % Eliminate Rows With Duplicate Independent Variable Values
Au1(D1(2:end)==0,:) = [];                                           % Eliminate Rows With Duplicate Independent Variable Values
D2 = diff([0; Au2(:,1)]);                                           % Eliminate Rows With Duplicate Independent Variable Values
Au2(D2(2:end)==0,:) = [];                                           % Eliminate Rows With Duplicate Independent Variable Values

xrange = [min([Au1(:,1);Au2(:,1)])  max([Au1(:,1);Au2(:,1)])];      % Independent Variable Extremes
xi = linspace(xrange(1), xrange(2), 200);                           % Interpolation Vector
Au1i = interp1(Au1(:,1), Au1(:,2), xi, 'linear', 'extrap');         % Interpolate ‘Au1’, Allow For Linear Extrapolation
Au2i = interp1(Au2(:,1), Au2(:,2), xi, 'linear', 'extrap');         % Interpolate ‘Au2’, Allow For Linear Extrapolation

figure(1)
plot(xi, Au1i)
hold on
plot(xi, Au2i)
hold off
grid

IntAu1 = cumtrapz(xi,Au1i);
IntAu2 = cumtrapz(xi,Au2i);

figure(2)
plot(xi, IntAu1-IntAu2, '-b', 'LineWidth',1.5)
hold on
plot(xi, IntAu1+1, '-m', 'LineWidth',1)
plot(xi, IntAu2+1, '--c', 'LineWidth',1)
hold off
grid
set(gca, 'YScale','log')
xlabel('X')
ylabel('Integral')
legend('Area Difference', 'Area 1 (‘coordinates’)', 'Area 2 (‘alpha’)', 'Location','SE')

The plot:

[Paschalis Garouniatis]

Thanks a lot for your answer and your time Star Strider. Didn't comment sooner because I wanted to understand the code and not just fit it in my script. So this question is answered.

[Star Strider]

My pleasure.

[Vanika Sharma]

I know the value of area of intersection and data points for one of the curve. Is there any way to get the second curve or points for that?

Answer 2

[John D'Errico]

If the two curves are defined over the same support, so identical x values for each curve, then the simple answer is to use trapz on the difference if the y-values. So...

trapz(x,y1-y2)

This is equivalent to computing trapz for each curve, then subtracting the two results, since integration is a linear operator. So, IF they are defined on the same interval in x...

trapz(x1,y1) - trapz(x2-y2)

If the two curves cross at one or more points, then, assuming you want to compute the area of each section as being a positive area, summing them up, then you COULD use trapz on the absolute difference.

trapz(x,abs(y1-y2))

the problem is this fails to be fully accurate because at the intersection point, if that falls between two of your sampled data points, the abs introduces an error.

So a better scheme would be to identify the crossing points using linear interpolation. Insert new values into each curve at those locations, then trapz will have higher accuracy. Doug Schwarz has a nice tool to locate those intersections on the file exchange, here:

intersections

If the two curves are not defined over the same support, then you will best use some interpolation method to put the two curves onto the same support, then apply one of the schemes I describe above.

[Paschalis Garouniatis]

Thanks a lot for your response John. I found the intersection points via the file you suggested. I am still having problem though with the interpolation. I tried interp1 but I am getting an error that my data is not strictly monotonic. I understand why this errors pops up but don't know how to overcome this. Any suggestions?

[Image Analyst]

You'll get faster and better answers if you just simply post your data so people can try things and show you.

[Paschalis Garouniatis]

Thank you Image Analyst for your suggestion. I included 3 files, coordinates1.mat is the original data file which contains pairs of x and y coordinates for the first curve, coordinates2.mat for the second curve and intersection.mat contains the intersection points between them. Any help is most welcome.

[John D'Errico]

Oh well, I knew that by not being complete there in show how to do the interpolation, I'd be forced to add that part anyway. :)

It is easy enough. intersections does the interpolation for you.

I'll keep the points with the same support foreach function, so x will be used for both curves.

x = linspace(0,2*pi,10);
y1 = sin(x);
y2 = cos(x);

plot(x1,y1,'r-o',x2,y2,'b-o')

[xnew,ynew] = intersections(x,y1,x,y2)
Warning: NARGCHK will be removed in a future release. Use NARGINCHK or NARGOUTCHK instead. 
> In intersections (line 91) 
xnew =
         0.790220733000934
          3.92363040454986
ynew =
         0.687902741912839
        -0.667001165195401

(It looks like Doug needs to fix intersections to remove that warning message.) Those are the new points we need to add to each curve.

x = [x,xnew'];
y1 = [y1,ynew'];
y2 = [y2,ynew'];

xold = x; % save for a later plot
[x,tags] = sort(x);
y1 = y1(tags);
y2 = y2(tags);

plot(x,y1,'r-o',x,y2,'b-o')

Now, I happen to know that the integral of abs(sin(x)-cos(x)), over the limits [0,2*pi] is 4*sqrt(2).

First, compare the two absolute differences.

plot(xold,abs(sin(xold) - cos(xold)),'r--',x,abs(y1-y2),'b:')

As you can see, just taking the abs of the difference WITHOUT inserting the points of intersection of the curves, I will get a poor result.

trapz(x,abs(y1-y2))
ans =
          5.43086292018145

Test that result:

syms u
int(abs(sin(u) - cos(u)),u,[0,2*pi])
ans =
4*2^(1/2)

vpa(ans)
ans =
5.6568542494923801952067548968388

IF the two functions did not live on the same support, then use interp1 to put them on the same support.

Answer 3

[Zeeshan Salam]

Can anyon solve this question?

How i find area between two curves?

clear;

clc;

%define x for both curves

x_ref = 0:0.01:1;

x = 0:0.01:1;

%define the reference curve and the other curve

c_ref = sin(x_ref*2*pi);

y = 0.5+0.5*sin(x*2*pi);

%plot both curves

figure(1)

plot (x_ref,c_ref,'r');

hold on

plot(x,y,'b');

axis([min(x_ref)-1,max(x_ref)+1,min(c_ref)-1,max((c_ref))+1]);

legend('c_{ref}','y')

%the relative x&y positions

x_rel = abs(x-min(x)/(max(x)-min(x)));

y_rel = abs((y-min(y))/(max(y)-min(y)));

%% Parameter before shifting

%corner points (in p1 matrix)

p1(:,1) = [min(x); max(x); min(x); max(x)]; % x colomn

p1(:,2) = [min(y); min(y); max(y); max(y)]; % y Spalte

for i = 1:size(p1,1)

text(p1(i,1),p1(i,2),[' p',num2str(i)])

end

%plot corner Points as circles

plot(p1(:,1),p1(:,2),'ko')

% hold off

%finding Xup, Xlow, Yright, Yleft

x_up_temp = (p1(4,1)-p1(3,1))*x_rel;

x_up1 = x_up_temp;

x_low_temp = (p1(2,1)-p1(1,1))*x_rel;

x_low1 = x_low_temp;

y_right_temp = (p1(4,2)-p1(2,2))*y_rel;

y_right1 = y_right_temp;

y_left_temp = (p1(3,2)-p1(1,2))*y_rel;

y_left1 = y_left_temp;

%% Parameter after shifting the first corner

%p1 = before shifting

%p2 = first point shifting

for i=1:3

p = [i-2 0; 0 0; 0 0; 0 0];

for j=1:3

p2 = p1 + p + [0 j-2; 0 0; 0 0; 0 0];

%finding new Xup, Xlow, Yright, Yleft

f1 = p2(4,1)-p2(3,1);

x_up_temp = f1*x_rel;

x_up2 = x_up_temp;

f2 = p2(2,1)-p2(1,1);

x_low_temp = f2*x_rel;

x_low2 = x_low_temp;

f3 = p2(4,2)-p2(2,2);

y_right_temp = f3*y_rel;

y_right2 = y_right_temp;

f4 = p2(3,2)-p2(1,2);

y_left_temp = f4*y_rel;

y_left2 = y_left_temp;

%calculating delta_x & delta_y

x_up_d = x_up2 - x_up1;

x_low_d = x_low2 - x_low1;

delta_x = x_up_d - x_low_d;

delta_y = (y_right2 - y_right1) - (y_left2 - y_left1);

%% the new x&y values are

xn = x + delta_x .* y_rel;

yn = y + delta_y .* x_rel;

cn = [xn; yn];% should drow this as a new curve

% calculate the Area

%P2(i;j) = flaeche();

end

end

%find the minimum area and apply the shifting

%a = find(flaeche == min(min(flaeche)));

%P2{a}

%% p3, p4, p5 would be the same

[Steven Lord]

Rather than posting this as an answer on a more than three year old question, you should ask this as a separate question of its own.