MATLAB 中具有任意数量交叉点的曲线之间的区域

Area between curves with arbitary number of intersections in MATLAB

我想计算 matlab 中具有不同 x 范围的两条曲线之间的面积。有两种情况。当一条曲线(图中的红色)高于另一条曲线时,我希望该区域为正,当该曲线低于蓝色曲线时,该区域为负。我意识到这通常是一个简单的问题,通常会找到两条曲线之间的所有交点,但在我的例子中,交点的数量可能非常大。我的困难是,根据所讨论的数据集,两条曲线将相交任意次数或根本不相交。

示例:

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);

有什么建议

您计算交叉点之间面积的方法是合理的。您可以使用 this library to calculate intersection points. Direct Link 压缩文件。

假设您有交叉点,您可以使用 trapz 函数遍历该交叉点以减去两个间隔之间的面积。如果曲线 1 高于曲线 2,唯一棘手的部分就是面积符号。为此,您可以设置一个阈值 epsilon 并在 x2 - epsilon 处计算 y1y2其中 (x1, x2)X 交点。 y1 > y2 意味着你必须做 Area 1 - Area2.

综上所述,代码如下所示:

%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);