Неверный учет вторых производных (интерполяция многочленами Ньютона с кратными узлами)
Необходимо реализовать интерполяцию по заданным значениям функции и узлов (кратных) многочленами Ньютона любого порядка до l = 10. Из-за кратности узлов соответствующие разностные отношения нужно заменять производными. Мой код решает задачу в следующем порядке:
- считает матрицу разностных отношений (по схеме Эйткена)
- выбирает узлы интерполяции путем минимизации остаточного члена
- и, собственно, интерполирует
Код в целом все делает верно:
- значения многочлена равны значениям функции в узлах
- первые производные совпадают с данными таблицы
НО (и в этом проблема) почему-то неверно учитываются вторые производные. Проверку осуществил в 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();
}
