"huijing " <denghuijing820@hotmail.com> wrote in message <gve859$hde$1@fred.mathworks.com>...
> > So what would this do?
> >
> > trapz(max(y,0))
> >
> > What is unsimple about it?
> >
> > John
>
> Thank you John. I was trying to do in this way:
> m=length(y);
> dt =0.5;
> ss=dt .* (y(1:m1,:) + 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:m1) + 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 nonnegative 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 + (kh)*x
So we cross the x axis at
x_k = h/(hk)
and the area above the line y == 0 is
A = h/(2*(hk))
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
