Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
area under curve

Subject: area under curve

From: huijing

Date: 25 May, 2009 10:27:01

Message: 1 of 8

Hello,
I am trying to get area value under curve using integration functions, like trapz or cumtrapz,
e.g.
x=[1:15];
y=[10 9 4 5 6 -1 -5 3 7 8 -3 -4 -1 5 0];
area(y)

the problem is how can I just get the area above zero in a simpler way.

Thank you!

Subject: area under curve

From: John D'Errico

Date: 25 May, 2009 11:01:02

Message: 2 of 8

"huijing " <denghuijing820@hotmail.com> wrote in message <gvdrpl$bli$1@fred.mathworks.com>...
> Hello,
> I am trying to get area value under curve using integration functions, like trapz or cumtrapz,
> e.g.
> x=[1:15];
> y=[10 9 4 5 6 -1 -5 3 7 8 -3 -4 -1 5 0];
> area(y)
>
> the problem is how can I just get the area above zero in a simpler way.
>
> Thank you!

So what would this do?

trapz(max(y,0))

What is un-simple about it?

John

Subject: area under curve

From: tpl@eng.cam.ac.uk (Tim Love)

Date: 25 May, 2009 11:12:05

Message: 3 of 8

"John D'Errico" <woodchips@rochester.rr.com> writes:


>So what would this do?

>trapz(max(y,0))

For
x=[0 1 2]
y=[1 -1 1]

I think the original poster would want the answer 0.5, not 1.

Subject: area under curve

From: John D'Errico

Date: 25 May, 2009 12:46:02

Message: 4 of 8

tpl@eng.cam.ac.uk (Tim Love) wrote in message <gvdue5$ee4$1@gemini.csx.cam.ac.uk>...
> "John D'Errico" <woodchips@rochester.rr.com> writes:
>
>
> >So what would this do?
>
> >trapz(max(y,0))
>
> For
> x=[0 1 2]
> y=[1 -1 1]
>
> I think the original poster would want the answer 0.5, not 1.

Not necessarily. The use of trapz implies an approximation,
as does ANY numerical integration. Your implied answer is
just as much of an approximation, just a different one. In
fact, there are many approximations for any such problem
in mathematics. All are equally valid in a sense, since we
cannot know the true underlying function that the data
points have been sampled from, therefore we cannot know
the true integral to be computed.

So what you "think" is irrelevant. I will note that while you
told us all what you think, you did not bother to actually
tell us what the "correct" solution is. Since there is NO
correct solution here, only various choices of approximation,
that explains why you did not do so, but only try to shoot
down the option that I posed.

My point is, while you might have offered a constructive
alternative, you did not bother to do so.

John

Subject: area under curve

From: huijing

Date: 25 May, 2009 13:58:01

Message: 5 of 8

> So what would this do?
>
> trapz(max(y,0))
>
> What is un-simple about it?
>
> John

Thank you John. I was trying to do in this way:
m=length(y);
dt =0.5;
ss=dt .* (y(1:m-1,:) + y(2:m,:));
ss=ss(ss>0);
area=sum(ss);

I am not sure about it. and the result is 48. while I have tried with your suggestion, and the result is 52.
Please tell me where the problem of my try is from the mathematic perspective?

Subject: area under curve

From: Steven Lord

Date: 25 May, 2009 14:28:43

Message: 6 of 8


"huijing " <denghuijing820@hotmail.com> wrote in message
news:gve859$hde$1@fred.mathworks.com...
>> So what would this do?
>>
>> trapz(max(y,0))
>>
>> What is un-simple about it?
>>
>> John
>
> Thank you John. I was trying to do in this way:
> m=length(y);
> dt =0.5;
> ss=dt .* (y(1:m-1,:) + y(2:m,:));
> ss=ss(ss>0);
> area=sum(ss);
>
> I am not sure about it. and the result is 48. while I have tried with your
> suggestion, and the result is 52.
> Please tell me where the problem of my try is from the mathematic
> perspective?

There are several different ways to modify the data so it stays positive.


x = [1 2 3];
y = [1 -1 1];

subplot(3, 1, 1)
x1 = x;
y1 = y;
plot(x1, y1, '-', x, y, 'r--')
title('Unmodified data')

subplot(3, 1, 2)
x2 = x;
y2 = max(y, 0);
plot(x2, y2, '-', x, y, 'r--')
title(Set any points in the original data with y < 0 to 0')

subplot(3, 1, 3)
x3 = [1 1.5 2 2.5 3];
y3 = [1 0 0 0 1];
plot(x3, y3, '-', x, y, 'r--')
title('Linear interpolate the function and cut it off at y = 0')


Note that both the latter two subplots satisfy the requirement that the
function not go below y = 0, but they do so using two different methods.

--
Steve Lord
slord@mathworks.com

Subject: area under curve

From: John D'Errico

Date: 25 May, 2009 15:28:01

Message: 7 of 8

"huijing " <denghuijing820@hotmail.com> wrote in message <gve859$hde$1@fred.mathworks.com>...
> > So what would this do?
> >
> > trapz(max(y,0))
> >
> > What is un-simple about it?
> >
> > John
>
> Thank you John. I was trying to do in this way:
> m=length(y);
> dt =0.5;
> ss=dt .* (y(1:m-1,:) + y(2:m,:));
> ss=ss(ss>0);
> area=sum(ss);
>
> I am not sure about it. and the result is 48. while I have tried with your suggestion, and the result is 52.
> Please tell me where the problem of my try is from the mathematic perspective?

You have supplied only a list of points sampled
from a function. What function will you presume
to have generated the points? Is it truly piecewise
linear? Probably not, but that need not be a poor
approximation.

x=[1:15];
y=[10 9 4 5 6 -1 -5 3 7 8 -3 -4 -1 5 0];
plot(x,y,'o-')
grid on

Most likely, the original function was not truly
piecewise linear, but I'd not be surprised to learn
that there is noise in your data too. This complicates
things a bit more.

You wish to compute the area above y == 0,
but below the curve y(x). Look back at the
approximation you have tried. What does it do?

m=length(y);
dt =0.5;
ss=dt .* (y(1:m-1) + y(2:m));

Take just these first few lines. What do they do?
This just trapezoidal rule. You have taken each
pair of points, and formed a trapezoid.

ss =
  Columns 1 through 7
          9.5 6.5 4.5 5.5 2.5 -3 -1
  Columns 8 through 14
            5 7.5 2.5 -3.5 -2.5 2 2.5

See that some of those trapezoids have negative
"areas". A negative area here just means that
this part of the curve fell below the line y == 0,
at least in part.

Then, at the end, you delete those trapezoidal
segments that gave negative areas. In fact, this
may not be the best choice, although as I suggested,
any choice is an approximation. What matters is
that you understand what approximation you
have made, and choose intelligently from among
the choices possible.

ss=ss(ss>0);
area=sum(ss);

When you delete a trapezoid that yields a negative
area from the sum, there may have been a significant
component above the y = 0 line. The alternative
that I offered was to use

trapz(max(y,0))

This too is an approximation, a quite simple one.
But you did ask for a simple solution in your post. It
adjusts those trapezoids that had one negative edge,
turning them effectively into triangles, now with a
fully non-negative area. You can appreciate the
difference from this plot:

plot(x,y,'bo-',x,max(y,0),'r-*')
grid on

See that the red curve will have a somewhat higher
area above the line y == 0 than the blue curve.
It will also be different from the trapezoidal
variant that you tried, where you simply deleted
some trapezoids. In fact, if we presumed a
piecewise linear interpolant from this data, then
the integral over that domain will be different from
either solution. A better solution might have us
using interpolation to find that point where the
curve actually crosses the line y == 0, breaking
the curve there. This would be exact, at least if
we presume an originally piecewise linear function.

The downside is that better solution will be more
complicated. You originally asked for "simple". So
you need to decide whether your need for accuracy
is more important than your desire for simplicity.

And, of course, if we were to presume a smooth
interpolant for this data, perhaps a spline, then
the integral above y == 0 will be a bit messier yet.

Can we find a reasonable solution here? Suppose
I posed a simple problem. Compute the area above
the line y == 0 for the linear function on the
interval [0,1], where y(0) = h, and y(1) = k.

When k > 0, the area is simple, and is given by
the trapezoidal area. For a unit width trapezoid,
we get

   A = (h + k)/2

Suppose k is negative? What is the area? We only
want that portion above the line y == 0. Where
does the line intersect? The solution comes from
writing the equation of this line. It is

  y = h + (k-h)*x

So we cross the x axis at

  x_k = h/(h-k)

and the area above the line y == 0 is

  A = h/(2*(h-k))

We could use this to get the integral of the
piecewise linear function, at least that portion
above y == 0. Or, we could insert new data
points in this curve at the intersections as
derived.

Either one will work. A reasonable solution might
have us start out not unlike your approach.

  ind = 1:(length(x)-1)
  dx = diff(x);
  ipos = find((y(ind) >= 0) & (y(ind+1) >= 0));
  A = sum(dx(ipos).*(y(ipos) + y(ipos+1))/2);

This takes only the segments where both edges
are positive. Then add in those pieces where we
crossed zero.

 ineg = find(((y(ind) > 0) & (y(ind+1) < 0)) | ...
      ((y(ind) < 0) & (y(ind+1) > 0)));

 A = A + sum(dx(ineg).*max(y(ineg)+abs(y(ineg+1)))./ ...
       (abs(y(ineg)+abs(y(ineg+1)))))/2

A =
       46.411

If I've done my back of the envelope algebra correctly,
this would be exact, but only for the piecewise linear
case. It is also not as simple as trapz(max(y,0)). Simple
is not always terribly good. But good is also not
always terribly simple. The case where the function
is not presumed to be linear will be messier yet.

You get what you pay for.

John

Subject: area under curve

From: huijing

Date: 26 May, 2009 12:35:01

Message: 8 of 8

John, I really appreciate for your reply. It's quite clear and logical. Also, I have tried your suggestion on my data, and it worked though not precise, but enough.
In fact, the data I am dealing with is from biological experiments. There is no function, just recording data, and in order to respect the original experimental result, I am not allowed to fit the curve. The curve of the data is kind of like the curve generated by the vectors I made up in the poster.
I'd like to show you the real data, but I don't know how to present it here.
Anyway, yesterday after reading your suggestion carefully, I tried to revise my code like this:

trapz(x(a:b),max(y(a:b),0));

a,b are indices of selected two points on graph( b>a), I wanna get area between these two pints, of course, omitting the negative area.

And it seems that it works out, because the result makes sense. it's bigger than the previous value which I got from
trapz(x(a:b),y(a:b)); which has calculated the negative area in addition to positive area.

By the way, I am reading the code you wrote years ago about how to select data points graphically, though the recent MATLAB has provided the Brush and Datalink, your code still helped to get better understand. Thank you again, you are big help to me.

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us