Area between two curves without intersection

Hi everyone.

I would like to calculate the area between two curves (see attachment).

I have used the trapz and polyarea function but these do not work because the curves are superimposed on each other without intersection.

Could someone please help me?

x=[0,-1.66128688049154,-3.71843384492024,-6.03903044153544,-8.52179344691878,-11.0684783490837,-13.5546470968919,-15.8324222826000,-17.7183932651871,-19.0241469744085,-19.5928500159198,-19.2715861063891,-18.0007426019886,-15.8645181604802,-13.1050363442789,-10.1023488165208,-7.36091774112053,-5.20617947547990,-3.73171179000825,-2.80200973068434,-2.17303027950468,-1.64687923319577,-1.14657359693582,-0.683525804975150,-0.261378364022391,0.167766706401400,0.646305657504070,1.13828938838504,1.49690490570653,1.50928713174259,0.991312332299208];

y=[0,3.80978510632932,8.62533289690098,13.7146367945814,18.4791657883180,22.5818350036559,26.0354921546852,29.1006123032565,31.9144606689929,33.8206797037834,33.4764658683987,29.5025727265991,21.5010392613040,10.9350767496541,13.3409505060801,18.4784450622125,22.0669438223010,22.8065275890501,20.8225751873008,17.4709163499800,14.5397350670517,13.2179835701845,13.4862976982229,14.3068759442281,14.3461752980756,12.9833980247855,10.3519556168192,7.26333666757206,4.38583660003191,2.05034193641872,0.472418853310666];

Answer 1

[darova]

Make sure curves have the same start and end

xx = linspace(x1(1),x1(end),100);   % new mesh
y11 = interp1(x1,y1,xx);            % interpolate curve1
y22 = interp1(x2,y2,xx);            % interpolate curve2
A = trapz(xx,abs(y22-y11));         % calculate positive area

[Torsten]

x and y are not two different curves, but one curve given in a (x,y) representation (like e.g. (x,y) = (cos(t),sin(t)) for a circle)

[Crocola Cool]

#Darova

Thanks for your feedback but it doesn't work with your proposal. You should not separate the x and y data. this said, (x,y).

#code

x=[0,-1.66128688049154,-3.71843384492024,-6.03903044153544,-8.52179344691878,-11.0684783490837,-13.5546470968919,-15.8324222826000,-17.7183932651871,-19.0241469744085,-19.5928500159198,-19.2715861063891,-18.0007426019886,-15.8645181604802,-13.1050363442789,-10.1023488165208,-7.36091774112053,-5.20617947547990,-3.73171179000825,-2.80200973068434,-2.17303027950468,-1.64687923319577,-1.14657359693582,-0.683525804975150,-0.261378364022391,0.167766706401400,0.646305657504070,1.13828938838504,1.49690490570653,1.50928713174259,0.991312332299208];

y=[0,3.80978510632932,8.62533289690098,13.7146367945814,18.4791657883180,22.5818350036559,26.0354921546852,29.1006123032565,31.9144606689929,33.8206797037834,33.4764658683987,29.5025727265991,21.5010392613040,10.9350767496541,13.3409505060801,18.4784450622125,22.0669438223010,22.8065275890501,20.8225751873008,17.4709163499800,14.5397350670517,13.2179835701845,13.4862976982229,14.3068759442281,14.3461752980756,12.9833980247855,10.3519556168192,7.26333666757206,4.38583660003191,2.05034193641872,0.472418853310666];

time=[10,39,69,99,129,158,188,218,248,277,307,337,367,397,426,456,486,516,545,575,605,635,665,694,724,754,784,813,843,873,903];

figure(1)

plot(x,y,'-x');

tq=min(time):1:max(time);

interp_x = interp1(time,x,tq);

interp_y= interp1(time,y,tq);

figure(2)

plot(interp_x,interp_y,'-O')

A=trapz(tq,abs(interp_y-interp_x))

Answer 2

[Paul]

I think this is what you're looking for:

x=[0,-1.66128688049154,-3.71843384492024,-6.03903044153544,-8.52179344691878,-11.0684783490837,-13.5546470968919,-15.8324222826000,-17.7183932651871,-19.0241469744085,-19.5928500159198,-19.2715861063891,-18.0007426019886,-15.8645181604802,-13.1050363442789,-10.1023488165208,-7.36091774112053,-5.20617947547990,-3.73171179000825,-2.80200973068434,-2.17303027950468,-1.64687923319577,-1.14657359693582,-0.683525804975150,-0.261378364022391,0.167766706401400,0.646305657504070,1.13828938838504,1.49690490570653,1.50928713174259,0.991312332299208];
y=[0,3.80978510632932,8.62533289690098,13.7146367945814,18.4791657883180,22.5818350036559,26.0354921546852,29.1006123032565,31.9144606689929,33.8206797037834,33.4764658683987,29.5025727265991,21.5010392613040,10.9350767496541,13.3409505060801,18.4784450622125,22.0669438223010,22.8065275890501,20.8225751873008,17.4709163499800,14.5397350670517,13.2179835701845,13.4862976982229,14.3068759442281,14.3461752980756,12.9833980247855,10.3519556168192,7.26333666757206,4.38583660003191,2.05034193641872,0.472418853310666];
plot(x,y,'-+')

p=polyshape(x,y);

Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.

plot(p)

p.area

ans = 202.0491

[Crocola Cool]

No I'm not looking for that. I've done this before but I'm really looking for the area between the two curves . No I'm not looking for that. I have done this before but I am really looking for the area between the two curves. Because Polyarea and trapz don't work. It gives me a very low value.

[Torsten]

The answer

area = 202.0491

comes close to what I expect by approximating the upper part by a triangle and the lower part by a rectangle.

[Crocola Cool]

#Torsten

Thank you for your feedback.

Your suggestion is correct but it doesn't work for the other cases. Here is an example.

Area= 44.3634 with your method and I find Area= 42.3287 with polyarea.

clear all

x=[0,-1.58143962981297,-3.03392588365747,-4.31293028513336,-5.43316939640765,-6.43266612032772,-7.34751230187266,-8.19610716571275,-8.97053149511526,-9.63779510279898,-10.1519335567333,-10.4585616901778,-10.5342395749024,-10.3758946306706,-9.98217498098106,-9.34944176705073,-8.47251867761138,-7.35956551245431,-6.05095668830785,-4.62502866942295,-3.18355609435549,-1.82611854945020,-0.611557284260889,0.477054737792964,1.47779860348098,2.40953293828952,3.22796956291444,3.81329910953911,4.00379047888411,3.66637669160807,2.77294873115999,1.43501593689118];

y=[0,-0.0914589891686976,-0.169258183169858,-0.137641932843007,0.115356308463203,0.650321396134900,1.44445426860180,2.40398635583359,3.39379566406802,4.25865376154095,4.84684849930363,5.03502448735570,4.83627564548227,4.37126677838120,3.80309691225722,3.29016826646019,2.94696153078223,2.83659552454877,2.98508831160315,3.38247203682421,3.96326780273134,4.58673760893439,5.02975831885881,5.09870485313676,4.70518969374756,3.90891301840727,2.89444462683932,1.89004771352925,1.07492003950258,0.526156770959942,0.222079396840722,0.0776605885521645];

figure(1)

plot(x,y)

p=polyshape(x,y);

figure (2)

plot(p)

Aire=p.area

poly=polyarea(x,y)

[Torsten]

"My" method is just an estimate for the magnitude of the area ( a means to decide whether the area returned by another program is realistic or not ). What polyshape gives under p.area will be far more exact. But the two numbers don't differ that much for your last example - so I'd trust in what polyshape returns.

[Paul]

I think polyarea is keeping track of the direction of the points. So the big blob in the lower right is a positive area, and the little blob in the upper left is a negative area. Whereas polyshape() by default fixes things so that all the vertices traverse in a consistent direction around the boundary. Consider a rectangle made up of two unit squares:

xy = [0 1;1 1;1 0;0 0;0 1;-1 1;-1 0;0 0];
plot(xy(:,1),xy(:,2),'-x');
axis([-2 2 -0.5 1.5]);

polyarea(xy(:,1),xy(:,2))

ans = 0

p1 = polyshape(xy(:,1),xy(:,2)); % default, note the message about adjusting the points

Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.

plot(p1)

p1.area

ans = 2

p2 = polyshape(xy(:,1),xy(:,2),'Simplify',false);
plot(p2)

p2.area

ans = 0

So with actual example data:

x=[0,-1.58143962981297,-3.03392588365747,-4.31293028513336,-5.43316939640765,-6.43266612032772,-7.34751230187266,-8.19610716571275,-8.97053149511526,-9.63779510279898,-10.1519335567333,-10.4585616901778,-10.5342395749024,-10.3758946306706,-9.98217498098106,-9.34944176705073,-8.47251867761138,-7.35956551245431,-6.05095668830785,-4.62502866942295,-3.18355609435549,-1.82611854945020,-0.611557284260889,0.477054737792964,1.47779860348098,2.40953293828952,3.22796956291444,3.81329910953911,4.00379047888411,3.66637669160807,2.77294873115999,1.43501593689118];
y=[0,-0.0914589891686976,-0.169258183169858,-0.137641932843007,0.115356308463203,0.650321396134900,1.44445426860180,2.40398635583359,3.39379566406802,4.25865376154095,4.84684849930363,5.03502448735570,4.83627564548227,4.37126677838120,3.80309691225722,3.29016826646019,2.94696153078223,2.83659552454877,2.98508831160315,3.38247203682421,3.96326780273134,4.58673760893439,5.02975831885881,5.09870485313676,4.70518969374756,3.90891301840727,2.89444462683932,1.89004771352925,1.07492003950258,0.526156770959942,0.222079396840722,0.0776605885521645];
plot(x,y,'-x')

p1=polyshape(x,y);

Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.

plot(p1)

p2 = polyshape(x,y,'Simplify',false)

p2 = 
  polyshape with properties:

      Vertices: [32×2 double]
    NumRegions: 1
      NumHoles: 0

plot(p2)

Aire1=p1.area

Aire1 = 44.3634

Aire2=p2.area

Aire2 = 42.3287

poly=polyarea(x,y)

poly = 42.3287

[Crocola Cool]

thanks Paul.

In example 1: the area is 202.0491

In example 2: the area is 44.3634.

According to the visual analysis, the area of this example must be greater than the area of the first one.

Do you agree with me?

[Torsten]

I'd estimate the length of the big region as 10 and its height as 4, and 10x4 = 40. So no, the area of this example will be much smaller than the area for the first one.

[Paul]

Plotting both shows that area in example 2 is much smaller than in example 1.

x1=[0,-1.66128688049154,-3.71843384492024,-6.03903044153544,-8.52179344691878,-11.0684783490837,-13.5546470968919,-15.8324222826000,-17.7183932651871,-19.0241469744085,-19.5928500159198,-19.2715861063891,-18.0007426019886,-15.8645181604802,-13.1050363442789,-10.1023488165208,-7.36091774112053,-5.20617947547990,-3.73171179000825,-2.80200973068434,-2.17303027950468,-1.64687923319577,-1.14657359693582,-0.683525804975150,-0.261378364022391,0.167766706401400,0.646305657504070,1.13828938838504,1.49690490570653,1.50928713174259,0.991312332299208];
y1=[0,3.80978510632932,8.62533289690098,13.7146367945814,18.4791657883180,22.5818350036559,26.0354921546852,29.1006123032565,31.9144606689929,33.8206797037834,33.4764658683987,29.5025727265991,21.5010392613040,10.9350767496541,13.3409505060801,18.4784450622125,22.0669438223010,22.8065275890501,20.8225751873008,17.4709163499800,14.5397350670517,13.2179835701845,13.4862976982229,14.3068759442281,14.3461752980756,12.9833980247855,10.3519556168192,7.26333666757206,4.38583660003191,2.05034193641872,0.472418853310666];
x2=[0,-1.58143962981297,-3.03392588365747,-4.31293028513336,-5.43316939640765,-6.43266612032772,-7.34751230187266,-8.19610716571275,-8.97053149511526,-9.63779510279898,-10.1519335567333,-10.4585616901778,-10.5342395749024,-10.3758946306706,-9.98217498098106,-9.34944176705073,-8.47251867761138,-7.35956551245431,-6.05095668830785,-4.62502866942295,-3.18355609435549,-1.82611854945020,-0.611557284260889,0.477054737792964,1.47779860348098,2.40953293828952,3.22796956291444,3.81329910953911,4.00379047888411,3.66637669160807,2.77294873115999,1.43501593689118];
y2=[0,-0.0914589891686976,-0.169258183169858,-0.137641932843007,0.115356308463203,0.650321396134900,1.44445426860180,2.40398635583359,3.39379566406802,4.25865376154095,4.84684849930363,5.03502448735570,4.83627564548227,4.37126677838120,3.80309691225722,3.29016826646019,2.94696153078223,2.83659552454877,2.98508831160315,3.38247203682421,3.96326780273134,4.58673760893439,5.02975831885881,5.09870485313676,4.70518969374756,3.90891301840727,2.89444462683932,1.89004771352925,1.07492003950258,0.526156770959942,0.222079396840722,0.0776605885521645];
plot(x1,y1,'-x',x2,y2,'-o'),grid