Кусочно-квадратичная интерполяция, полином Лагранжа

class Program
{
    delegate double Foo(double x);
 
    static void Main(string[] args)
    {
        double[,] grid = { { 1, 4.1 }, { 2, 3.4 }, { 3, 2.9 }, { 4, 2.6 }, { 5, 2.5 } };
        PiecewiseInterpolant l = new Lagrange(grid), q = new Quadratic(grid);
        Foo L = t => l.Calc(t); // Интерполянта Лагранжа
        Foo Q = t => q.Calc(t); // Квадратичная интерполянта
 
        // Строим табличку для сравнения
        int n = grid.GetLength(0);
        double h = (grid[n - 1, 0] - grid[0, 0]) / ((n-1) * 4);
        string pattern = "{0, 6} | {1, 6} | {2, 6} | {3, 6}";
        Console.WriteLine(pattern, "x", "L(x)", "Q(x)", "Abs(L-Q)");
        for (double x = grid[0, 0]; x <= grid[n - 1, 0]; x += h )
        {
            double ly = L(x), qy = Q(x);
            Console.WriteLine(pattern,
                x,
                ly.ToString("F4"), 
                qy.ToString("F4"), 
                Math.Abs(ly - qy).ToString("F4"));
        }
 
        Console.ReadLine();
    }
}
 
/// <summary>
/// Базовый класс всех интерполянт
/// </summary>
abstract class PiecewiseInterpolant
{
    protected double[] x, y;
 
    public PiecewiseInterpolant(double[,] Points)
    {
        if (Points == null) throw new ArgumentNullException();
 
        int size = Points.GetLength(0);
        if (size < 3 || Points.GetLength(1) != 2)
            throw new ArgumentException("Wrong format");
 
        var points = Enumerable.Range(0, size)
            .Select(i => new { x = Points[i, 0], y = Points[i, 1] })
            .OrderBy(p => p.x)          // сортируем
            .GroupBy(p => p.x)          // Проверяем, 
            .Select(g => g.First())     // чтобы x' не совпадали
            .ToArray();
 
        if (points.Length != size)
            throw new ArgumentException("values X must be unique");
 
        x = points.Select(p => p.x).ToArray();
        y = points.Select(p => p.y).ToArray();
    }
 
    public double Calc(double X)
    {
        // Смотрим, вдруг, X совпадет с одним из значений на сетке
        int idx = Array.BinarySearch(x, X); // Тут нужен компаратор 
        if (idx >= 0) return y[idx];
 
        // Если не совпал, проверим, чтобы лежал внутри наименьшего 
        // и наибольшего значения
        idx = -(idx + 1);
        if (idx == 0 || idx == x.Length)
            throw new ArgumentException(String.Format(
                "Argument must be in intereval: {0}..{1}",
                x[0], x[x.Length -1]));
 
        // Если все хорошо, начнем вычислять интерполянту 
        // на данном промежутке [x{i}, x{i+1}]
        return InnerCalc(idx - 1, X);
    }
 
    protected abstract double InnerCalc(int idx, double x);
}
 
class Quadratic : PiecewiseInterpolant
{
    double[,] a;
    public Quadratic(double[,] Points) : base(Points)
    {
        // Вычисляем коэффициенты парабол для всех интервалов,
        // здесь не будем пропускать среднюю точку, 
        // для последнего интервала сдвинем параболу влево
        a = new double[x.Length, 3];
        for (int i = 0; i < x.Length - 1; i++)
        {
            int d = i == x.Length - 2 ? 1 : 0;
            Fill(i, i - d, i + 1 - d, i + 2 - d);
        }
    }
 
    /// <summary>
    /// Вычисляет коэффициенты параболы по трем точкам
    /// </summary>
    private void Fill(int idx, int i1, int i2, int i3)
    {
        double x1 = x[i1], x2 = x[i2], x3 = x[i3];
        double y1 = x[i1], y2 = y[i2], y3 = y[i3];
 
        // просто формулы из учебника
        a[idx, 2] = (y3 - y1) / ((x3 - x1) * (x3 - x2));
        a[idx, 2] -= (y2 - y1) / ((x2 - x1) * (x3 - x2));
        a[idx, 1] = (y2 - y1) / (x2 - x1) - a[idx, 2] * (x2 + x1);
        a[idx, 0] = y1 - a[idx, 1] * x1 - a[idx, 2] * x1 * x1;
    }
 
    protected override double InnerCalc(int idx, double X)
    {
        return a[idx, 0] + (a[idx, 1] + a[idx, 2] * X) * X;
    }
}
 
class Lagrange: PiecewiseInterpolant
{
    double[] l;
 
    public Lagrange(double[,] Points) : base(Points)
    {
        l = new double[x.Length];
        for (int i = 0; i < x.Length; i++)
        {
            l[i] = y[i];
            for (int j = 0; j < x.Length; j++)
                if (i != j) l[i] /= x[i] - x[j];
        }
    }
 
    protected override double InnerCalc(int idx, double X)
    {
        double y = 0;
        for (int i = 0; i < x.Length; i++)
        {
            int j = 0;
            double t = l[i];
            while (j < i) t *= X - x[j++];
            while (++j < x.Length) t *= X - x[j];
            y += t;
        }
 
        return y;
    }
}