C#借助开源数值分析库求解两条离散曲线交点前最后一点

原逻辑问题说明 #

你之前的判断逻辑存在三个明显缺陷：

• 最后一个判断条件y1 == Y2没有物理意义，不符合曲线相交的判定规则
• 未处理浮点计算的精度误差，直接用等号判断会漏过大量临界场景
• 仅判断了采样点上相交的情况，没有覆盖两个采样点间隔之间发生相交的场景

基于MathNet.Numerics的实现方案 #

MathNet.Numerics是.NET生态应用最广泛的开源数值计算库，可以直接通过NuGet安装使用，不需要额外配置。

核心思路 #

由于两条曲线是相同t坐标采样的，我们可以通过差值的符号变化判断相交：

1. 构造差值序列delta[i] = y[i] - Y[i]
2. 遍历差值序列，当相邻两个差值的乘积小于等于0时，说明该区间内两条曲线发生相交，区间左端点就是相交前的最后一个采样点
3. 如果需要获取精确交点坐标，可以通过线性插值/样条插值计算

代码实现 #

基础版（线性插值，性能最优，适合超大点集）

using System;
using MathNet.Numerics;

// 输入数据：等长的t坐标序列、第一条曲线y值序列、第二条曲线Y值序列
double[] t = /* 你的t坐标数组 */;
double[] y = /* 第一条曲线y值数组 */;
double[] Y = /* 第二条曲线y值数组 */;

int? lastPointIndex = null;
double intersectionT = double.NaN;
double intersectionY = double.NaN;
const double precision = 1e-9; // 浮点精度阈值，可根据业务需求调整

for (int i = 0; i < y.Length - 1; i++)
{
    double deltaCurr = y[i] - Y[i];
    double deltaNext = y[i + 1] - Y[i + 1];

    // 差值符号变化，说明区间内存在交点
    if (deltaCurr * deltaNext <= precision)
    {
        lastPointIndex = i;
        // 计算精确交点
        double t1 = t[i], t2 = t[i + 1];
        double y1 = y[i], y2 = y[i + 1];
        double Y1 = Y[i], Y2 = Y[i + 1];

        double denominator = (y2 - y1) - (Y2 - Y1);
        if (Math.Abs(denominator) < precision)
        {
            // 区间内两条线平行，直接取当前点为交点
            intersectionT = t1;
            intersectionY = y1;
        }
        else
        {
            double u = (-deltaCurr) / denominator;
            intersectionT = t1 + u * (t2 - t1);
            intersectionY = y1 + u * (y2 - y1);
        }
        break;
    }
}

// 结果输出
if (lastPointIndex.HasValue)
{
    Console.WriteLine($"相交前最后一个采样点索引：{lastPointIndex.Value}");
    Console.WriteLine($"采样点坐标：(y={y[lastPointIndex.Value]}, t={t[lastPointIndex.Value]})");
    Console.WriteLine($"精确交点坐标：(y={intersectionY}, t={intersectionT})");
}
else
{
    Console.WriteLine("给定区间内两条曲线无交点");
}

进阶版（三次样条插值，适合平滑曲线场景）

如果你的曲线是平滑曲线，需要按样条插值计算交点，可以使用以下代码：

using System;
using MathNet.Numerics.Interpolation;
using MathNet.Numerics.RootFinding;

double[] t = /* 你的t坐标数组 */;
double[] y = /* 第一条曲线y值数组 */;
double[] Y = /* 第二条曲线y值数组 */;
const double precision = 1e-9;

// 构造两条曲线的三次自然样条插值函数
var splineY = CubicSpline.InterpolateNatural(t, y);
var splineY2 = CubicSpline.InterpolateNatural(t, Y);

// 构造差值函数 f(t) = y(t) - Y(t)
Func<double, double> diffFunc = (x) => splineY.Interpolate(x) - splineY2.Interpolate(x);

// 查找差值为0的根，即交点的t坐标
if (Brent.TryFindRoot(diffFunc, t[0], t[^1], precision, out double rootT))
{
    // 查找小于交点t的最后一个采样点索引
    int searchResult = Array.BinarySearch(t, rootT);
    int lastPointIndex = searchResult < 0 ? (~searchResult - 1) : (searchResult - 1);
    double rootY = splineY.Interpolate(rootT);
    
    Console.WriteLine($"相交前最后一个采样点索引：{lastPointIndex}");
    Console.WriteLine($"采样点坐标：(y={y[lastPointIndex]}, t={t[lastPointIndex]})");
    Console.WriteLine($"精确交点坐标：(y={rootY}, t={rootT})");
}
else
{
    Console.WriteLine("给定区间内两条曲线无交点");
}

内容的提问来源于stack exchange，提问作者Leon