Неверный учет вторых производных (интерполяция многочленами Ньютона с кратными узлами)

Необходимо реализовать интерполяцию по заданным значениям функции и узлов (кратных) многочленами Ньютона любого порядка до l = 10. Из-за кратности узлов соответствующие разностные отношения нужно заменять производными. Мой код решает задачу в следующем порядке:

  • считает матрицу разностных отношений (по схеме Эйткена)
  • выбирает узлы интерполяции путем минимизации остаточного члена
  • и, собственно, интерполирует

Код в целом все делает верно:

  1. значения многочлена равны значениям функции в узлах
  2. первые производные совпадают с данными таблицы

НО (и в этом проблема) почему-то неверно учитываются вторые производные. Проверку осуществил в origin, попросив его дважды продифференцировать график, посмотрев значения вторых производных и сравнив их с табличными.

В чем может быть проблема?

Прикрепляю значения функции, узлов и производных, а также сам код.

Таблица

P.S. Программа способна интерполировать функцию многочленами Ньютона любого порядка до l = 10, исходя из этого матрица разностных отношений имеет соответствующую размерность. В первом столбце матрицы стоят значения функции для удобства дальнейших расчетов

static void Main(string[] args)
    {

        double[] fxx = { 0.0341, 0, 0, 0, -0.0103, 0, 0, 0, -0.0085, 0, 0 }; //первые производные из таблицы
        double[] fxxx = { -0.00164, 0, 0, 0, -0.001, 0, 0, 0, 0, 0, 0 }; //вторые производные из таблицы
        double[] X = { 1, 1, 1, 1.3, 1.7, 1.7, 1.7, 2.0, 2.4, 2.4, 3.0}; //узлы
        double[] fxArray = { 2.598, 2.598, 2.598, 2.643, 2.619, 2.619, 2.619, 2.586, 2.546, 2.546, 2.509 }; //значения функции

        //создаем и готовим зубчатый массив для разностных отношений
        int k = fxArray.Length - 1;
        double[][] fid = new double[k + 1][];
        fid[0] = fxArray; 
        int i = 1;
        do
        {
            fid[i] = new double[k];
            k = k - 1;
            i++;
        } while (k != 0);

        //вычисляем разностные отношения
        for (i = 0; i < (fid.Length - 1); i++)
        {
            for (int j = 0; j < (fid[i + 1].Length); j++)
            {
                if ((X[j] == X[j + 1]) && (i == 0))
                {
                    if ((X[j] == X[0]))
                    {
                        fid[i + 1][j] = fxx[0];
                    }
                    else if ((X[j] == X[4]))
                    {
                        fid[i + 1][j] = fxx[4];
                    }
                    else
                    {
                        fid[i + 1][j] = fxx[8];
                    }
                }
                else if ((j != X.Length - 2) && (X[j] == X[j + 2]) && (i == 1))
                {
                    if ((X[j] == X[0]))
                    {
                        fid[i + 1][j] = fxxx[0] / 2; 
                    }
                    else
                    {
                        fid[i + 1][j] = fxxx[4] / 2;
                    }
                }
                else
                {
                    fid[i + 1][j] = (fid[i][j + 1] - fid[i][j]) / (X[i + j + 1] - X[j]);
                }
            }
        }

        //проверяем результат (вывод матрицы разностных отношений)
        for (int m = 0; m < fid.Length; m++)
        {
            for (int n = 0; n < fid[m].Length; n++)
            {
                Console.Write($"{fid[n][m]} \t");
            }
            Console.WriteLine();
        }

        
        Console.WriteLine("Укажите порядок полинома: ");
        int l = Convert.ToInt32(Console.ReadLine());

        //количество сдвижек для поиска минимального остаточного члена
        int v = (X.Length - (l + 1));
        int w = 0;
        double[] arr = new double[v];

        double L = 0;
        double C = 1;
        double x = X[0];
        double h = 0.01;
        int g = (int)(2 / h) + 1;
        int p = 0;
        double[] Ln = new double[g];
        double[] Xi = new double[g];
        for (i = 0; i < Xi.Length; i++)  
        {
            //перебор и выделение узлов, обеспечивающих минимизацию остаточного члена
            for (w = 0; w <= (v - 1); w++)
            {
                double R = 1;
                for (int j = 0; j <= l; j++)
                {
                    double r = (x - X[j + w]);
                    R = R * r;
                }
                arr[w] = R;
            }
            for (int n = 1; n < arr.Length; n++)
            {
                if ((Math.Abs(arr[n - 1]) < Math.Abs(arr[n])))
                {
                    w = n - 1;
                    break;
                }
            }

            Console.WriteLine($"Кол-во сдвижек узлов на {i+1} шаге: {w}");

            //вычисление многочлена
            for (int n = 0; n <= l; n++)
            {
                double ln = fid[p][w];
                for (int j = 0; j < n; j++)
                {
                    C = C * (x - X[j + w]);
                }
                L = L + C * ln;
                C = 1;
                p++;
            }
            Ln[i] = L;
            Xi[i] = x;
            x = x + h;
            L = 0;
            p = 0;
        }

        StreamWriter sw = new StreamWriter(@"C:\Users\79872\Desktop\Interpol");
        sw.WriteLine("x\tLn(x)");
        for (int n = 0; n < Xi.Length; n++)
        {
            sw.WriteLine(string.Format("{0:0.0000E000}\t{1:0.0000E000}", Xi[n], Ln[n]));
        }
        sw.Close();
    }

Ответы (0 шт):