求两条水平重叠曲线之间的积分

Finding the integral between two horizontally overlapping curves

我有一个代码可以绘制文本文件中的数据。我已经能够得到每条数据曲线下的积分(虽然我不知道它们中的任何一个的功能,所以我只是使用 integral = s = simps(y, x))。其中两条曲线水平重叠,我想找到它们重叠之间的区域。我在这方面遇到了一些麻烦,因为它们水平重叠而不是垂直重叠。我想找到的区域是这个 graph 上的红线和蓝线之间。数据在 numpy 数组中,例如 bxby 代表蓝线 x 和 y,rxry 代表红线 x 和 y。我没有代表它们的功能,但它们在那些数组中。很难知道如何进行,因为它们水平重叠而不是垂直重叠。

我对此深有体会,如有任何帮助,我们将不胜感激!

更新:JohanC has worked in one case but not another. The non-working plot initially looks like this (I want the area only of that middle overlapping area), but when that did not work either, I had tried limiting the x range to hopefully help the code narrow in on the area where I would like the overlapping integral to be calculated (that is where that plot with the smaller range came from). The code I am currently using is the one JohanC提供的代码如下所示:

def get_overlap_integral(bx, by, rx, ry):
    fig, ax = plt.subplots()
    ax.plot(bx, by, color='b')
    ax.plot(rx, ry, color='r')

    # create a common x range
    gx = np.linspace(min(bx.min(), rx.min()), max(bx.max(), rx.max()), 1000)
    gby = np.interp(gx, bx, by)  # set bx,by onto the common x range
    gry = np.interp(gx, rx, ry)  # set rx,ry onto the common x range
    gy = np.minimum(gby, gry)  # let gx,gy be the intersection of the two curves
    ax.fill_between(gx, gy, color='g', alpha=0.3)

    area = np.trapz(gy, gx)  # calculate the green area
    ax.text(0.05, 0.95, f'Overlap: {area:.3f}', color='g', ha='left', va='top', transform=ax.transAxes)
    plt.show()
    return area

文本 bx、by、rx 和 ry 的值如下:

BX: [999.5 999.  998.5 ... 201.  200.5 200. ]
BY: [-3.786867e-05 -4.366451e-05 -4.745308e-05 ...  1.068685e-05  1.555391e-05
 -8.949840e-06]
RX: [999.5 999.  998.5 ... 201.  200.5 200. ]
RY: [-5.865443e-05 -7.808241e-05 -5.887286e-05 ... -1.556630e-06 -3.473830e-06
 -6.367470e-06]

我不确定为什么这个功能在一种情况下不起作用,但在另一种情况下却可以完美地工作。任何帮助将不胜感激!

np.interp(new_x, old_x, old_y) 可以计算新 x 范围的曲线。当设置在公共 x 范围内时,np.minimum 找到两条正值曲线的公共区域。 np.trapz 计算面积。

下面是一些示例代码,从创建测试数据开始:

from matplotlib import pyplot as plt
import numpy as np

# create some test data
Nb = 300
bx = np.linspace(320, 750, Nb)
by = np.random.randn(Nb).cumsum()
by -= by[0]  # start at zero
by -= np.linspace(0, 1, Nb) * by[-1]  # let end in zero
by = (by ** 2) ** .3  # positive curve starting and ending with zero
Nr = 200
rx = np.linspace(610, 1080, Nr)
ry = np.random.randn(Nr).cumsum()
ry -= ry[0]  # start at zero
ry -= np.linspace(0, 1, Nr) * ry[-1]  # let end in zero
ry = (ry ** 2) ** .3  # positive curve starting and ending with zero

fig, ax = plt.subplots()
ax.plot(bx, by, color='b')
ax.plot(rx, ry, color='r')

# create a common x range
gx = np.linspace(min(bx.min(), rx.min()), max(bx.max(), rx.max()), 1000)
gby = np.interp(gx, bx, by)  # set bx,by onto the common x range
gry = np.interp(gx, rx, ry)  # set rx,ry onto the common x range
gy = np.minimum(gby, gry)  # let gx,gy be the intersection of the two curves
ax.fill_between(gx, gy, color='g', alpha=0.3)

area = np.trapz(gy, gx)  # calculate the green area
ax.text(0.05, 0.95, f'Overlap: {area:.3f}', color='g', ha='left', va='top', transform=ax.transAxes)
plt.show()