Matlab trapz but splitting areas above and below the x-axis

Question

I'm trying to find the area under a curve in the region both above the x-axis (from y = 0 to +inf) and below the x-axis, (y = 0 to -inf), separately.

As far as I can figure out, I can use

int = trapz(x,y)

to get the area under the curve, but not as 2 separate values.

For example, if we have x = [0,1,1,2] & y = [7,7,-4,-4]

Trapz will work out "3" as the area ((7.*1) + (-4.*1)), but I want the values of the area above the x-axis and below the x-axis as 2 separate numbers, 7 & -4.

Trapz must fit the trapezoids, calculate the area of the trapezoids and then sum over all. I want a way to fit & find the areas of the trapezoids, but then sort them in to +ve & -ve and finds the sum's of the +ve and -ve areas.

The documentation for trapz: https://au.mathworks.com/help/matlab/ref/trapz.html

Answer 1

You just need to be careful about how you handle the segments where a zero-crossing happens. Picchiolu's answer is much simpler to implement, but ignores this case and you will introduce errors with that approach.

The clearest way is probably to just write your own implementation of the trapezium rule and handle each segment accordingly, with particular focus on handling segments which cross y=0.

See the comments for details

% Demo data
x = [1 2 3 6 7 9 12 13 14 16];
y = sin(3*pi*(x/max(x)));

% Trapz with sign change
aPos = 0; % positive area
aNeg = 0; % negative area
% Loop over point pairs
for ii = 2:numel(x)
    idx = (ii-1:ii);    % two points we care about
    dx = diff(x(idx));  % change in x
    my = sum(y(idx))/2; % mean y
    sgn = sign(y(idx));
    if prod(sgn) >= 0 % Both y values are same sign or one/both equals 0
        % two y points have same sign, simple trapezium area added to +/-
        if max(sgn) > 0 % at least one y is positive
            aPos = aPos + (dx * my);
        else % 0 or -1 sign
            aNeg = aNeg + (dx * my);
        end
    else % change in sign
        pct = 1 - y(ii) / diff(y(idx)); % percentage into interval  
        if sgn(1) < 0
            pct = 1 - pct;      
        end
        aPos = aPos + sum( y(idx)/2 .* pct .* dx .* (y(idx)>0) );
        aNeg = aNeg + sum( y(idx)/2 .* (1-pct) .* dx .* (y(idx)<0) );
    end
end

% Total area
aTotal = aPos + aNeg;

Plot to visualise the integral. Note that the data points do not include the zero-crossing points, but the integral creates distinct segments in the above loop

As a sense check, using the above gives the same output for aTotal as trapz(x,y).