Широко признано - и было указано в комментариях - что матрица Вандермонда часто плохо обусловлена. Например, стандартный учебник по численному анализу, D. Kincaid и W. Cheney, "Numerical Analysis, 2nd ed.", Brooks / Cole Publishing 1996, утверждает на странице 336:
Матрица коэффициентов в Уравнение (12) называется матрицей Вандермонда . Это не единственное число, потому что система имеет уникальное решение для любого выбора y 0 , y 1 , ..., y n [...]. Следовательно, определитель матрицы Вандермонда отличен от нуля для различных узлов x 0 , x 1 , ..., x п . [...] Однако матрица Вандермонда часто плохо обусловлена , и поэтому коэффициенты a i могут быть определены неточно путем решения Системы (12). (См. Gautschi [1984]). [...] Следовательно, этот подход не рекомендуется .
Проверка матрицы, возвращаемой GETRF, уже сообщает нам, что у нас проблемы с полиномами довольно низкого градусов. При степени 4 величина наименьшего релевантного элемента после разложения LU составляет порядка 10 -12 , а при степени 5 наименьший такой элемент составляет 10 -14 . Исходя из того факта, что все элементы исходной матрицы были близки к единице, и зная основные c шагов процесса разложения LU, ясно, что крошечные элементы в результате GETRF получаются путем вычитания.
Величина элементов, приближающихся к эпсилону двойной точности (≈ 10 -16 ), говорит нам, что большая часть точности была потеряна. Для полиномов степени 6 и выше код в основном работает на чистом шуме.
Мы можем дополнительно подтвердить это, сравнив вычисленные коэффициенты интерполяционного полинома с более надежным эталоном. Для простоты я использовал Wolfram Alpha для этого сравнения. Для полинома степени 4 код на Фортране вычисляет коэффициенты с точностью примерно до трех десятичных цифр, а для полинома степени 5 это сжимается до единственного правильного десятичного числа di git.
В терминах простой и более надежный алгоритм для генерации коэффициентов полинома интерполяции, я обнаружил следующее:
JN Lyness и C .B. Молер, "Системы Ван-дер-Монд и численное дифференцирование". Numerische Mathematik 8, 458-464 (1966)
Я перевел код Algol в статье на ISO-C99, и он, кажется, дает разумные результаты до степени 8 по сравнению с Wolfram Alpha . Wolfram Alpha отказывается от оценок выше 8, и у меня нет под рукой других справочных материалов. Даже с этим более надежным алгоритмом, похоже, наблюдается значительная потеря точности для более высоких степеней, с совпадением только 6 десятичных цифр между Wolfram Alpha и алгоритмом Luness / Moler.
#include <stdio.h>
#include <stdlib.h>
/* J. N. Lyness and C.B. Moler, "Van Der Monde Systems and Numerical
Differentiation." Numerische Mathematik 8, 458-464 (1966)
*/
void update (int k, int p, double *x, double fxk, double *C)
{
int d, s, m, n;
double xk, xkd;
xk = x[k];
m = k*(k+3)/2;
C[m] = fxk;
for (d = 1; d <= k; d++) {
xkd = xk - x[k-d];
for (s = 0; s <= ((d > p) ? p : d); s++) {
m = m - 1;
n = m + d;
if (s == 0) {
C[m] = C[n] + xk * (C[n-k-1] - C[n]) / xkd;
} else if (s == d) {
C[m] = (C[n+1] - C[n-k]) / xkd;
} else {
C[m] = C[n] + (xk * (C[n-k-1] - C[n]) + (C[n+1] - C[n-k]))/ xkd;
}
if (d > p) {
m = m - d + p;
}
}
}
}
/* Solve (k+1) x (k+1) system Vc = f, where V[i,j] = x[i]**j */
void vandal (int k, double *x, double *f, double *c)
{
double C [k*(k+3)/2+1];
for (int i = 0; i <= k; i++) {
update (i, k, x, f[i], C);
}
for (int i = 0; i <= k; i++) {
c[i] = C[k-i];
}
}
#define k 10
int main (void)
{
double x[k+1] = {1.000, 1.001, 1.002, 1.003, 1.004, 1.005,
1.006, 1.007, 1.008, 1.009, 1.010};
double f[k+1] = {1.00000000000, 1.01004512021, 1.02018096336,
1.03040825707, 1.04072773401, 1.05114013204,
1.06164619412, 1.07224666847, 1.08294230847,
1.09373387280, 1.10462212541};
double c[k+1] = {0};
vandal (k, x, f, c);
for (int i = 0; i <= k; i++) {
printf ("c[%2d] = % 23.16e\n", i, c[i]);
}
return EXIT_SUCCESS;
}