'trapz' to find area under curve not working

.fig file attached.

For all positive y values, I want to find the area under the graph.

if y>0
    for x = length(y)
    int = trapz(x, y)
    end
end

However that gives me int = 0 and that doesn't seem right.

Any idea what I am doing wrong?

Thanks!

Answer 1

[Walter Roberson]

x=length(y) is a scalar value, just the length. You then pass that one length as the X coordinates and pass all of y as the Y coordinates. When you pass a scalar as the X coordinates then trapz treats it as the difference between successive X values. It seems unlikely that the difference between X values just happens to be the same as the number of Y values.

Your for loop is processing only that one calculation because the =length(y) part is a scalar.

If you were to change that to be a range like 1:length(y) then you would be overwriting the int variable each iteration.

Your y is a vector. When you test a vector in an if statement the test is true only if all of vector values are nonzero simultaneously. Your body of the if will not fire once for each y that is greater than 0: it will fire once total if ALL of the y are greater than 0 and not at all if even one y is not positive.

Suppose your y is 1 2 3 -1 5 6 then how do you want to process only the positive values? Do you want to trapz 1 2 3 and add trapz 5 6 (each block of positive values independently) or do you want to process 1 2 3 5 6 as if the 5 were immediately beside the 3, or do you want to process like 1 2 3 0 5 6 (negative values become 0) or do you want to "hold" the last positive value through the gap, like 1 2 3 3 5 6?

[JD]

Walter,

Thank you for responding. I would like to calculate trapz 1 2 3 and then add that to traps 5 6 because I believe that is what will give me what I'm looking for (area under the curve above y=0)

Is the code below doing what I want it to?

int = trapz(x, y>0)

[JD]

Ok that is definitely not doing the right thing. However, I don't know how to go about writing the code to do this. Any suggestions?

[JD]

I cannot construct a for loop, because for y > 0 is not a valid input, which is why I used if. But if I understand what you said above, that will not compute if it is not true for all y.

Following another method you suggested above would be to create a new array with only positive y. To do that, I can execute the following code

y_pos = y(y>0);

But now how do I get the x only corresponding to the above y's?

[Walter Roberson]

mask = y>0;
y_pos = y(mask);
x_pos = x(mask);

However, this would be the equivalent of holding the last positive value, like 1 2 3 3 5 6

I think you need to refine your needs a bit. Consider that taking trapz([1 2 3]) + trapz([5 6]) is like assuming that the input "suddenly starts" at 1, goes to 2 and then 3, and then "suddenly stops", and then "suddenly starts" at 5 and goes to 6 and "suddenly stops". Whereas if you were to think about 1 2 3 -1 5 6, it does not suddenly stop at 3, it drops to -1 passing through 0 as it does, and from -1 it goes to 5 passing through 0 as it does, so the area above 0 should include that drop from 3 to 0 and the rise from 0 to 5. You can calculate those intersection points without a lot of trouble and inject them into the (x,y) pairs.

Meanwhile, if you do want trapz([1 2 3]) + trapz([5 6]) then:

starts = strfind( [false, y>0], [0 1] );
stops  = strfind( [y>0, false], [1 0] );
tot_area = 0;
for K = 1 : length(starts)
    from = starts(K); to = stops(K);
    tot_area = tot_area + trapz(x(from:to), y(from:to));
end

[JD]

That make sense! Thought it would be easy to calculate the area under curve, but I guess it's not as straight forward.

Thanks so much for the help!