Символьное вычисление полинома из алгоритма Невилла

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

// fill C (allocated if null) with params for interpolating polynomial
// use params with interp_poly_eval
// !! these are NOT polynomial coefficients. 

double* interp_poly( Int deg, const double* x, const double* y, double* restrict C)
{
double* c = C ?: calloc( deg+1, sizeof *c);
Int i, j;
    memcpy( c, y, (deg+1)*sizeof *y);
    for (i=1; i<=deg; i++) 
    {   for (j=deg; j>=i; j--)
        {   c[j] = (c[j]-c[j-1]) / (x[j]-x[j-i]);
        }
    }
    return c;
}


double  interp_poly_eval( Int deg, const double* c, const double* x, double X)
{
double  p = c[deg];
Int i = deg;
    while( --i >= 0)
    {   p = c[i] + (X-x[i])*p;
    }
    return p;
}

// as above but also returns derivative of the polynomial through pdp
double  interp_poly_eval_d( Int deg, const double* c, const double* x, double X, double* pdp)
{
double  p = c[deg];
double  dp = 0.0;
Int i = deg;
    while( --i >= 0)
    {   dp = (X-x[i])*dp + p;
        p = c[i] + (X-x[i])*p;
    }   
    *pdp = dp;
    return p;
}

Вам нужно иметь класс/структуру данных, которая может представлять одномерный (по одной переменной) полином. Тогда легко вычислить фактический многочлен (т.е. его коэффициенты), а не только отдельные точки, используя сам алгоритм Невилла, т.е.

  P[0,0] = y[0]   ; constant
  P[1,1] = y[1]   ; constant

а потом рекурсивно

  P[a,b] = P[a,b-1] * x[b]/C + P[a,b-1] * -X/C + P[a+1,b] * -x[a]/C + P[a+1,b] * X/C

где X — полином «x» (т. е. моном первой степени), а C — константа X[b]-X[a].

Это дает вам фактический полином, когда вы, как обычно, следуете рекурсии в алгоритме.

Обратите внимание, что выше все арифметические действия относятся к полиному, т.е. P[a,b-1]*x[b]/C означает многочлен P[a,b-1], умноженный на константу x[b]/C (т.е. x-координата b-й точки, деленная на C).

Если вы хотите получить точный результат, используйте пакет рациональной арифметики произвольной точности (например, GMP для C/C++). В качестве альтернативы, вычисляйте с использованием чисел с плавающей запятой, но вы можете столкнуться с ошибками округления, которые повлияют на вас.