该函数无法检测两条曲线的交点

Question

下面的源代码应该找到两条曲线之间的交点。

但是，该功能无法检测交点。

我该如何修复它？

using MathNet.Numerics.Interpolation;
using System.Collections.Generic;
using System.Linq;

public static Vec2 FindIntersectionOfTwoCurves(List<double> xList1, List<double> yList1,
                                               List<double> xList2, List<double> yList2)
{
    if (xList1 == null || yList1 == null || xList2 == null || yList2 == null ||
        xList1.Count != yList1.Count || xList2.Count != yList2.Count)
        return null;

    IInterpolation interpolation1 = Interpolate.Linear(xList1, yList1);
    IInterpolation interpolation2 = Interpolate.Linear(xList2, yList2);

    double lowerBound = Math.Max(xList1.Min(), xList2.Min());
    double upperBound = Math.Min(xList1.Max(), xList2.Max());

    double step = (upperBound - lowerBound) / 1000; // adjust the step size as needed
    for (double x = lowerBound; x <= upperBound; x += step)
    {
        double y1 = interpolation1.Interpolate(x);
        double y2 = interpolation2.Interpolate(x);

        if (Math.Abs(y1 - y2) < 1e-7)
        {
            return new Vec2(x, y1);
        }
    }

    return null;
}

public class Vec2
{
    public double X { get; set; }
    public double Y { get; set; }

    public Vec2(double x, double y)
    {
        X = x;
        Y = y;
    }
}

Answer 1

由于您试图找到两条曲线与非连续一阶导数的交点，因此您应该使用 二分法。

一个快速的C#实现如下，二分查找的实现来自这里

namespace ConsoleApp2
{
    public class Interpolator
    {
        private static int BinarySearch<T>(IList<T> list, T value)
        {
            if (list == null)
                throw new ArgumentNullException("list");
            var comp = Comparer<T>.Default;
            int lo = 0, hi = list.Count - 1;
            while (lo < hi)
            {
                int m = (hi + lo) / 2;  // this might overflow; be careful.
                if (comp.Compare(list[m], value) < 0) lo = m + 1;
                else hi = m - 1;
            }
            if (comp.Compare(list[lo], value) < 0) lo++;
            return lo;
        }

        private static int FindFirstIndexGreaterThanOrEqualTo<T>
                                  (List<T> sortedList, T key)
        {
            return BinarySearch(sortedList, key);
        }

        List<double> x_values;
        List<double> y_values;
        public Interpolator(List<double> x, List<double> y)
        {
            List<int> indicies = x.AsEnumerable().Select((v,i) => new {obj = v, index = i}).OrderBy(c=> c.obj).Select(c => c.index).ToList();
            
            x_values = indicies.Select(i => x[i]).ToList();
            y_values = indicies.Select(i => y[i]).ToList();
        }

        public double Interpolate(double x)
        {
            int index = FindFirstIndexGreaterThanOrEqualTo(x_values, x);
            if (index == x_values.Count-1)
            {
                return y_values[^1];
            }
            double y1 = y_values[index];
            double y2 = y_values[index + 1];
            double x1 = x_values[index];
            double x2 = x_values[index+1];
            return (x - x1) / (x2 - x1) * (y2 - y1) + y1;
        }

    }

    class IntersectionFinder
    {
        public static Nullable<(double,double)> FindIntersectionOfTwoCurves(List<double> xList1, List<double> yList1,
                                                   List<double> xList2, List<double> yList2)
        {
            if (xList1 == null || yList1 == null || xList2 == null || yList2 == null ||
                xList1.Count != yList1.Count || xList2.Count != yList2.Count)
                return null;

            Interpolator interpolator1 = new Interpolator(xList1, yList1);
            Interpolator interpolator2 = new Interpolator(xList2, yList2);

            double lowerBound = Math.Max(xList1.Min(), xList2.Min());
            double upperBound = Math.Min(xList1.Max(), xList2.Max());

            double diff_start = interpolator1.Interpolate(lowerBound) - interpolator2.Interpolate(lowerBound);
            double diff_end = interpolator1.Interpolate(upperBound) - interpolator2.Interpolate(upperBound);
            if ((diff_start > 0 && diff_end > 0) || (diff_start < 0 && diff_end < 0)) // make sure they actually intersect
            {
                return null;
            }

            double mid = (lowerBound + upperBound) / 2;
            double diff_mid = interpolator1.Interpolate(mid) - interpolator2.Interpolate(mid);

            while (diff_mid > 1e-7)
            {
                if (diff_start > diff_end) // list1 is higher
                {
                    if (diff_mid > 0) // mid is also higher, intersection in right side
                    {
                        lowerBound = mid;
                        diff_start = diff_mid;
                    }
                    else // mid is lower, intersection in left side
                    {
                        upperBound = mid;
                        diff_end = diff_mid;
                    }
                }
                else // list 2 is higher
                {
                    if (diff_mid < 0) 
                    {
                        lowerBound = mid;
                        diff_start = diff_mid;
                    }
                    else 
                    {
                        upperBound = mid;
                        diff_end = diff_mid;
                    }
                }
                mid = (lowerBound + upperBound) / 2;
                diff_mid = interpolator1.Interpolate(mid) - interpolator2.Interpolate(mid);
            }

            return new (mid, interpolator1.Interpolate(mid));
        }
    }

    internal class Program
    {
        static void Main(string[] args)
        {
            List<double> xList1 = [ 1, 2, 3, 4, 5 ];
            List<double> xList2 = [ 1, 2, 3, 4, 5 ];
            List<double> yList1 = [ 1, 2, 3, 4, 5 ];
            List<double> yList2 = [ 5, 4, 3, 2, 1 ];
            Console.WriteLine(IntersectionFinder.FindIntersectionOfTwoCurves(xList1, yList1, xList2, yList2).ToString());
        }
    }
}