numpy 测量两条曲线之间的积分(线性函数和任意曲线)

eufgjt7s  于 2022-11-29  发布在  其他
关注(0)|答案(2)|浏览(136)

在下面的图中,我的目标是在1 / 2 / 3区域内求积分,这样我就知道直线下方有多少区域(1 / 3区域),直线上方有多少区域(2区域)。
我不是在寻找精确的积分,只是一个近似值来测量。一个近似值,它将以同样的方式工作,为其他版本的曲线,我已经表示。
y1:蓝线为线性函数y= -0.148x + 1301.35
y2:黄线是任意曲线
两条曲线共享相同的x轴。
image of curves linear & arbitrary curve
我试过几种方法,在这里找到了关于stackoverflow的方法,主要是这两种方法引起了我的注意:
https://stackoverflow.com/a/57827807
&
https://stackoverflow.com/a/25447819
他们给予我完全相同的输出为整个地区,我的问题是分开它上面/下面。
我的最佳尝试示例:(https://stackoverflow.com/a/25447819/20441461的修改版本)
y1 / y2 / x -是用于上图中曲线的数据

y1 = [1298.54771845, 1298.40019417, 1298.2526699, 1298.10514563, 
     1297.95762136,1297.81009709, 1297.66257282, 1297.51504854]

y2 = [1298.59, 1297.31, 1296.04, 1297.31, 1296.95, 1299.18, 1297.05, 1297.45]

x = np.arange(len(y1))

z = y1-y2
dx = x[1:] - x[:-1]
cross_test = np.sign(z[:-1] * z[1:])

x_intersect = x[:-1] - dx / (z[1:] - z[:-1]) * z[:-1]
dx_intersect = - dx / (z[1:] - z[:-1]) * z[:-1]

areas_pos = abs(z[:-1] + z[1:]) * 0.5 * dx # signs of both z are same   
areas_neg = 0.5 * dx_intersect * abs(z[:-1]) + 0.5 * (dx - dx_intersect) * abs(z[1:])

negatives = np.where(cross_test < 0)
negative_sum = np.sum(x_intersect[negatives])

positives = np.where(cross_test >= 0)
positive_sum = np.sum(x_intersect[positives])`

就是给予我这个结果
负积分= 10.15
正积分= 9.97
仅仅从图片上看,我就可以看出这不可能是正确的值。(直线下面的面积比上面的大得多。)
我已经花了大量时间在这上面了,而且很困--任何建议或建议都是受欢迎的。

efzxgjgh

efzxgjgh1#

下面是一段代码,它 * 精确 * 计算所有区域,并以 * 矢量化 * 的方式(快速)完成计算:

def areas(x, y1, y2, details=None):
    dy = y1 - y2
    b0 = dy[:-1]
    b1 = dy[1:]
    b = np.c_[b0, b1]
    r = np.abs(b0) / (np.abs(b0) + np.abs(b1))
    rr = np.c_[r, 1-r]
    dx = np.diff(x)
    h = rr * dx[:, None]
    br = (b * rr[:, ::-1]).sum(1)
    a = (b + br[:, None]) * h / 2
    result = np.sum(a[a > 0]), np.sum(a[a < 0])
    if details is not None:
        details.update(locals())  # for debugging
    return result

示例:

x = np.array([0,1,2])
y1 = np.array([1,0,3])
y2 = np.array([0,3,4])

>>> areas(x, y1, y2)
(0.125, -3.125)

您的原始示例:

y1 = np.array([
    1298.54771845, 1298.40019417, 1298.2526699, 1298.10514563, 
    1297.95762136,1297.81009709, 1297.66257282, 1297.51504854])

y2 = np.array([1298.59, 1297.31, 1296.04, 1297.31, 1296.95, 1299.18, 1297.05, 1297.45])

x = np.arange(len(y1))

>>> areas(x, y1, y2)
(5.228440802728334, -0.8687563377282332)
说明

为了理解它是如何工作的,让我们考虑四点的四边形:[a, b, c, d],其中ab位于相同的x位置,cd也位于相同的x位置。如果没有边相交,则可以是"直的",否则可以是"扭曲的"。在这两种情况下,我们都考虑扭曲版本相交的交点的x位置,以及在该位置的四边形的实际垂直截面(如果是扭曲的,则为0,或者如果是直的,则为垂直边的加权平均)。
假设我们把垂直距离称为b0b1,它们是有方向的(如果y1 > y2,则为正),交点在x坐标x + r * dx处,其中r = |b0| / (|b0| + |b1|)和是0和1之间的一个因子。
对于绞合四边形,左侧(三角形)区域为b0*r*dx/2,右侧区域为b1*(1-r)*dx/2
对于一个直四边形,左边区域(梯形)是(b0 + br)/2 * r * dx,右边是(b1 + br) / 2 * (1 - r) * dx,其中brr水平比例的高度,br = b0 * (1 - r) + b1 * r
一般来说,我们在计算中总是使用br。对于扭曲四边形,它是0,我们可以使用与直四边形相同的表达式。这是消除任何测试并产生纯矢量化函数的关键。
剩下的是一点numpy表达式来高效地计算所有这些值。

示例详细信息
def plot_details(details, ax=None):
    x, y1, y2, dx, r, a = [details[k] for k in 'x y1 y2 dx r a'.split()]
    ax = ax if ax else plt.gca()
    ax.plot(x, y1, 'b.-')
    ax.plot(x, y2, 'r.-')
    xmid = x[:-1] + dx * r
    [ax.axvline(xi) for xi in xmid]
    xy1 = np.c_[x, y1]
    xy2 = np.c_[x, y2]
    for A,B,C,D,r,(a0,a1) in zip(xy1, xy2, xy1[1:], xy2[1:], r, a):
        ACmid = A*(1-r) + C*r
        BDmid = B*(1-r) + D*r
        q0 = np.c_[A,ACmid,BDmid,B]
        q1 = np.c_[ACmid,C,D,BDmid]
        ax.fill(*q0, alpha=.2)
        ax.fill(*q1, alpha=.2)
        ax.text(*q0.mean(1), f'{a0:.3f}', ha='center')
        ax.text(*q1.mean(1), f'{a1:.3f}', ha='center')

以前面的例子为例:

x = np.array([0,1,2])
y1 = np.array([1,0,3])
y2 = np.array([0,3,4])

details = {}
>>> areas(x, y1, y2, details)
(0.125, -3.125)

>>> details
{'x': array([0, 1, 2]),
 'y1': array([1, 0, 3]),
 'y2': array([0, 3, 4]),
 'details': {...},
 'dy': array([ 1, -3, -1]),
 'b0': array([ 1, -3]),
 'b1': array([-3, -1]),
 'b': array([[ 1, -3],
        [-3, -1]]),
 'r': array([0.25, 0.75]),
 'rr': array([[0.25, 0.75],
        [0.75, 0.25]]),
 'dx': array([1, 1]),
 'h': array([[0.25, 0.75],
        [0.75, 0.25]]),
 'br': array([ 0. , -1.5]),
 'a': array([[ 0.125 , -1.125 ],
        [-1.6875, -0.3125]]),
 'result': (0.125, -3.125)}

还有:

plot_details(details)


中的一个

vddsk6oq

vddsk6oq2#

也许您可以对两个数组的绝对差进行积分:

>>> np.trapz(np.abs(y2 - y1))
7.1417718350001

相关问题