Как получить полиномиальные коэффициенты интерполяции с помощью gsl_interp?

Question

Итак, у меня есть код ниже. Он отлично вычисляет все y-точки многочлена (и печатает их для построения с помощью gnuplot), но как получить полученный полином (1-x² в этом случае)?

void twoDegreePoly() { 
    int n = 3; 
    double x[n],y[n]; 
    printf ("#m=0,S=16\n"); 
    for (int i=0; i<n ;i++) { 
     x[i] = ((double)2*i)/2 -1; 
     y[i] = f(x[i]); 
     printf ("%g %g\n", x[i], y[i]); 
    } 
    printf ("#m=1,S=0\n"); 

    gsl_interp_accel *acc = gsl_interp_accel_alloc(); 
    const gsl_interp_type *t = gsl_interp_polynomial; 
    gsl_interp* poly = gsl_interp_alloc(t,n); 
    gsl_interp_init (poly, x, y,n); 
    for (double xi=x[0]; xi<x[n-1]; xi+= 0.01) { 
     double yi = gsl_interp_eval (poly, x, y, xi, acc); 
     printf ("%g %g\n", xi, yi); 
    } 
} 

Answer 1

После быстрого сканирования над documentation, это не кажется, что такая функция доступна в GSL. Это может быть вызвано двумя причинами: во-первых, получение полиномиальных коэффициентов является особым для этого метода интерполяции, который не вписывается в общий дизайн (который может обрабатывать произвольные функции). Во-вторых, ссылаясь на числовые рецепты:

Пожалуйста, будьте уверены, что коэффициенты - это то, что вам нужно. Как правило, коэффициенты интерполяционного многочлена могут быть определены гораздо менее точно, чем его значение при желаемой абсциссе. При этом нецелесообразно определять коэффициенты только для использования при вычислении интерполяционных значений. Значения вычисленных таким образом не будет проходить точно через табличные точки, например, ...

Причина этого заключается в том, что в принципе вычисление коэффициентов предполагает решение линейной системы с матрицей Вандермонда, которая очень плохо -conditioned.

Тем не менее, в числовых рецептах приводится процедура polcoe, с помощью которой вы можете получить интерполяционный многочлен. Вы можете найти его в главе 3.5. в free second edition.

Answer 2

Я сделал что-то подобное с интерполяцией Акимы. Во-первых, определить состояние как GSL сделать:

typedef struct 
{ 
    double *b; 
    double *c; 
    double *d; 
    double *_m; 
}akima_state_t; 

Затем создайте интерполянт

spline = gsl_spline_alloc (gsl_interp_akima, M_size); 
gsl_spline_init (spline, x, y, M_size); 

и после этого, вы можете сделать:

const akima_state_t *state = (const akima_state_t *) (spline -> interp -> state); 
    double _b,_c,_d; 
    for (int i = 0; i < M_size; i++) 
    { 
    _b = state->b[i]; 
    _c = state->c[i]; 
    _d = state->d[i]; 
    std::cout << "(x>"<<x[i]<<")*(x<"<<x[i+1]<<")*("<<y[i]<< "+ (x-"<< x[i]<<")*("<<_b<<"+(x-"<< x[i]<<")*("<<_c<<"+"<<_d<<"*(x-"<<x[i]<<")))) + "; 
    } 

Я не попробовал с полиномиальная интерполяция, но здесь структура состояния для многочлена, она должна быть хорошей отправной точкой.

typedef struct 
{ 
    double *d; 
    double *coeff; 
    double *work; 
} 
polynomial_state_t;