Area between curves with arbitary number of intersections in MATLAB

Question

I am looking to calculate the area between two curves in matlab that have different x ranges. There are two cases. I want the area to be positive when one curve (red in graph) is above the other and negative when that curve is below the blue curve. I realise this is typically an easy problem and one would typically find all points of intersection between the two curves but in my case the number of intersections could be quite large. My difficulty is that depending on the dataset in question, the two curves will intersect any number of times or not at all.

Example:

n=4
x = n*pi*rand(100,1);
x=[0;x;pi];
x=sort(x);
y=.5+exp(-.5*x).*sin(x);
x2 = n*pi*rand(1000,1);
x2=[0;x2;pi];
x2=sort(x2);
y2=.5+exp(-.5*x2).*0.6.*sin(x2/1.2);

figure,hold on 
plot(x,y,'r-')
plot(x2,y2,'b-')

Area1  = trapz(x, y);
Area2  = trapz(x2, y2);

Any suggestions ?

Answer 1

Your approach to calculate area between intersections is reasonable. You can use this library to calculate intersection points. Direct Link to zip file.

Assuming you have the intersections you can loop through that to subtract area between two intervals using trapz function. The only tricky part remains is of the sign of area if curve 1 is above curve 2. For that you can have a threshold epsilon and calculate y1 and y2 at x2 - epsilon where (x1, x2) were the X intersection points. y1 > y2 will imply you have to do Area 1 - Area2.

With all the above, the code would look like below:

%Find all intersections
[xIntersect,yIntersect] = intersections(x,y,x2,y2,1);
%Number of intersections
numIntersections = size(xIntersect, 1);
%Threshold to calculate sign of area
epsilon = 0.1;

curveArea = 0; %Initial value
lastX1Index = 1; % End X1 Index processes in last iteration
lastX2Index = 1; % End X2 Index processes in last iteration

%Loop over intersections processing pair of points hence until
%numIntersections - 1
for i = 1:numIntersections-1
%     startX = xIntersect(i); %Start of interval 
%     startY = yIntersect(i); 
    endX = xIntersect(i+1); % Next Intersection point
%     endY = yIntersect(i+1);
    xMinus = endX - epsilon; % get a X value before intersection

    %Sign of area is positive if y1Minus > y2Minus
    y1Minus = .5+exp(-.5*xMinus).*sin(xMinus); %Convert this to function1
    y2Minus = .5+exp(-.5*xMinus).*0.6.*sin(xMinus/1.2); %Convert this to function 2

    x1Index = find(x < xIntersect(i+1), 1, 'last' ); % Find last X1 element before intersection
    x2Index = find(x2 < xIntersect(i+1), 1, 'last' ); % Find last X2 element before intersection

    %Calculate area for interval
    Area1 = trapz(x(lastX1Index:x1Index), y(lastX1Index:x1Index));
    Area2 = trapz(x2(lastX2Index:x2Index), y2(lastX2Index:x2Index));

    %Store last X1, X2 Index for next run
    lastX1Index = x1Index+1;
    lastX2Index = x2Index+1;

    %Save curveArea
    if y1Minus > y2Minus
        curveArea = curveArea + (Area1 - Area2);
    else
        curveArea = curveArea + (Area2 - Area1);
    end
end

fprintf('curveArea= %f \n',curveArea);