两组离散数据求交点

在做课题研究时，经常需要求解两组离散数据的交点。常见的做法是，先采用多项式模型或者指数模型等函数模型对两组离散数据数据进行拟合，得到拟合函数再使用solve求解两条曲线的精确解。但这样做有一个很明显的局限，当两组离散数据不满足上述任何一种模型时，那么使用上述模型拟合出来的曲线就会与实际的离散数据坐标点严重偏离，造成求解结果的误差非常大。因此，本文对两组离散数据（可以看做是两个分段线性函数）的交点进行求解总结出另一种求解方法，可直接求得两组离散数据折线（线性）的精确交点。

本文所述的求解两组离散数据折线的交点的求解精度由用户给定的两组离散数据的密集程度决定。用户所给的两组离散数据越密集，分段线性函数（离散数据折线）也就看起来也就越光滑， 所求得的解的精度自然就越高。

（1）两组离散数据（相同序列X）求交点

假定有这样两组离散数据（相同序列X）

第1组：

X = [1, 2, 3, 4, 5, 6, 7]
Y1 = [6, 5, 3, 10, 5, 6, 8]

第2组：

X = [1, 2, 3, 4, 5, 6, 7]
Y2 = [2, 5, 5, 4, 3, 5, 8]

两组两组离散数据（相同序列X）的坐标点如图所示：

由上图可以看出，两组离散数据的总共应该有3个交点，于是我们对两组离散数据求交点，代码如下：

numpy_get_roots_with_two_piecewise_linear_by_same_X.py：

import numpy as np
import matplotlib.pyplot as plt
# 支持中文
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号



#### 两组离散数据（同X）求交点 ####
def numpy_get_roots_with_two_piecewise_linear_by_same_X(
    Y1 = [6, 5, 3, 10, 5, 6, 8], 
    Y2 = [2, 5, 5, 4, 3, 5, 8], 
    X = None, 
    plot = 0
    ):
    # 注1：Y1与Y2长度相等，每条直线段最多只能有一个交点（不可能有两个）
    # 注1：X的值必须由小到大（X须为单调递增序列）
    if X == None:
        X = [i+1 for i in range(len(Y1))]
    else:
        pass
    
    roots = []
    for i in range(len(Y1)-1):
        Y_i_1 = [Y1[i], Y1[i+1]]
        Y_i_2 = [Y2[i], Y2[i+1]]
        X_i = [X[i], X[i+1]]

        poly_i_1 = np.poly1d(np.polyfit(X_i, Y_i_1, 1))
        poly_i_2 = np.poly1d(np.polyfit(X_i, Y_i_2, 1))
        poly_delta = poly_i_1 - poly_i_2
        
        # 非水平直线函数poly_delta必定有且只有一个解
        root_i = poly_delta.roots[0]
        y = poly_i_1(root_i)
        
        # 减少精度，找回因数据精度损失造成节点处丢失的根
        root_i = round(root_i, 6)
        y = round(y, 6)
        #print((root_i, y), X_i)
        if root_i >=  X_i[0] and root_i <= X_i[1]:
            #print((root_i, y), X_i)
            roots.append((root_i, y))

    if plot == 1:
        # 显示曲线
        plt.figure(figsize=(12, 6))
        plt.title('两组离散数据(相同序列X)求交点')
        plt.xlabel('X')
        plt.ylabel('Y')

        plt.plot(
            X, 
            Y1, 
            marker='*',
            label='X,Y1离散数据')

        plt.plot(
            X, 
            Y2, 
            marker='*',
            label='X,Y2离散数据')

        plt.legend()
        plt.show()

    roots = list(set(roots))
    roots.sort()
    return roots


if __name__ == '__main__':

    roots = numpy_get_roots_with_two_piecewise_linear_by_same_X(Y1=[6, 5, 3, 10, 5, 6, 8], Y2=[2, 5, 5, 4, 3, 5, 8], X=None, plot=1)
    print(roots)

输出两组离散数据（相同序列X）的交点坐标：

[(2.0, 5.0), (3.25, 4.75), (7.0, 8.0)]

（2）两组离散数据（不同序列X）求交点

假定有这样两组离散数据（不同序列X）

第1组：

X1 = [1, 1.8, 3, 4, 5]
Y1 = [5, 5, 9, 10, 5]

第2组：

X2 = [2, 3.2, 4.5, 5.1, 6, 8, 11]
Y2 = [8, 5, 5, 8, 3, 5, 4]

两组两组离散数据（不同序列X）的坐标点如图所示：

由上图可以看出，两组离散数据的总共应有2个交点。

先将不同序列X的两组离散数据转化为相同序列X的两组离散数据如下图

再调用numpy_get_roots_with_two_piecewise_linear_by_same_X.py中的函数对两组不同序列X离散数据求交点，代码如下：

import numpy as np
import matplotlib.pyplot as plt
# 支持中文
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号



from numpy_get_roots_with_two_piecewise_linear_by_same_X import *


# 求交集范围
def X1X2ab(X1ab=(1, 6), X2ab=(5,9)):
    '''
    xa, xb = X1X2ab(X1ab=(X1[0], X1[-1]), X2ab=(X2[0], X2[-1]))
    '''
    X1a = X1ab[0]
    X1b = X1ab[1]
    if X2ab == None:
        X2a = X1a
        X2b = X1b
    else:
        X2a = X2ab[0]
        X2b = X2ab[1]
        if X2a == None:
            X2a = X1a
        if X2b == None:
            X2b = X1b


    # 先确定X2a的位置 , r_left
    if X2a <= X1a:
        r_left = X1a
    elif X1a < X2a and X2a <X1b:
        r_left = X2a
    elif X2a >= X1b:
        return None
    # 再确定X2b的位置, r_right
    if X2b <= X1a:
       return None
    elif X1a < X2b and X2b <X1b:
        r_right = X2b
    elif X2b >= X1b:
        r_right = X1b
    xab = (r_left, r_right)
    return xab





def numpy_get_roots_with_two_piecewise_linear_any(

    X1 = [1, 1.8, 3, 4, 5],
    Y1 = [5, 5, 9, 10, 5],

    X2 = [2, 3.2, 4.5, 5.1, 6, 8, 11],
    Y2 = [8, 5, 5, 8, 3, 5, 4],

    plot = 0

    ):

    # 构造新的X_new
    X_ = X1 + X2
    X_.sort()
    X_new = []
    # 共同定义域
    a, b = X1X2ab(X1ab=(X1[0], X1[-1]), X2ab=(X2[0], X2[-1]))
    for i in X_:
        if i >= a and i <= b:
            X_new.append(i)
    X_new = list(set(X_new)) # 去重
    X_new.sort()
    #print(X_new)

    # 构造新的Y1_new、Y2_new
    Y1_new = np.interp(X_new, X1, Y1)
    Y2_new = np.interp(X_new, X2, Y2)

    roots = numpy_get_roots_with_two_piecewise_linear_by_same_X(Y1=Y1_new, Y2=Y2_new, X=X_new, plot=plot)
    return roots




if __name__ == '__main__':

    X1 = [1, 1.8, 3, 4, 5]
    Y1 = [5, 5, 9, 10, 5]

    X2 = [2, 3.2, 4.5, 5.1, 6, 8, 11]
    Y2 = [8, 5, 5, 8, 3, 5, 4]


    roots = numpy_get_roots_with_two_piecewise_linear_any(X1=X1, Y1=Y1, X2=X2, Y2=Y2, plot=1)
    print(roots)


    plt.figure(figsize=(12, 6))
    plt.title('两组离散数据(不同序列X)求交点')
    plt.xlabel('X')
    plt.ylabel('Y')

    plt.plot(
        X1, 
        Y1, 
        marker='*',
        label='X1,X1离散数据')

    plt.plot(
        X2, 
        Y2, 
        marker='*',
        label='X2,X2离散数据')

    plt.legend()
    plt.show()

输出两组离散数据（不同序列X）的交点坐标：

[(2.4, 7.0), (4.75, 6.25)]