Last Update: Oct. 4, 2007
/***********************************************/
/* T.Kouya's GSL sample program collection */
/* Interpolation of Discrete Points */
/* Written by Tomonori Kouya */
/* */
/* Version 0.01: 2007-10-04 */
/***********************************************/
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
/* Number of Discrete Points */
#define DIM 10
/* exact func(x) */
/* sin(x) */
double func(double x)
{
return sin(x);
}
int main(void)
{
int i, times, status;
gsl_interp *workspace;
gsl_interp_accel *accel;
double h;
double x[DIM], y[DIM], xtmp, ytmp;
/* Polynomial Interpolation */
/* Set discrete points */
x[0] = 0.0;
x[DIM - 1] = 2 * M_PI;
h = (x[DIM - 1] - x[0]) / DIM;
for(i = 0; i < DIM; i++)
{
x[i] = x[0] + (double)i * h;
y[i] = func(x[i]);
}
printf("Discrete Points for Interpolation (sin(x))\n");
printf(" i x y\n");
for(i = 0; i < DIM; i++)
printf("%3d: %25.17e %25.17e\n", i, x[i], y[i]);
printf("\n");
/* Define type of Interpolation */
workspace = gsl_interp_alloc(gsl_interp_polynomial, DIM);
accel = gsl_interp_accel_alloc();
/* Initialize */
gsl_interp_init(workspace, x, y, DIM);
/* Evaluation */
h /= 3.0;
printf(" x interpolated y abs_error \n");
for(i = 0; i < 3 * DIM; i++)
{
xtmp = x[0] + h * (double)i;
ytmp = gsl_interp_eval(workspace, x, y, xtmp, NULL);
printf("%25.17e %25.17e %10.3e\n", xtmp, ytmp, fabs(ytmp - func(xtmp)));
}
/* free */
gsl_interp_free(workspace);
return 0;
}
$ ./interpolation
Discrete Points for Interpolation (sin(x))
i x y
0: 0.00000000000000000e+00 0.00000000000000000e+00
1: 6.28318530717958623e-01 5.87785252292473137e-01
2: 1.25663706143591725e+00 9.51056516295153531e-01
3: 1.88495559215387587e+00 9.51056516295153642e-01
4: 2.51327412287183449e+00 5.87785252292473248e-01
5: 3.14159265358979312e+00 1.22464679914735321e-16
6: 3.76991118430775174e+00 -5.87785252292473026e-01
7: 4.39822971502571036e+00 -9.51056516295153531e-01
8: 5.02654824574366899e+00 -9.51056516295153642e-01
9: 5.65486677646162761e+00 -5.87785252292473359e-01
x interpolated y abs_error
0.00000000000000000e+00 0.00000000000000000e+00 0.000e+00
2.09439510239319532e-01 2.07860795644523644e-01 5.090e-05
4.18879020478639064e-01 4.06712580221638642e-01 2.406e-05
6.28318530717958623e-01 5.87785252292473137e-01 0.000e+00
8.37758040957278127e-01 7.43151997300360856e-01 7.172e-06
1.04719755119659763e+00 8.66029786417928271e-01 4.383e-06
1.25663706143591725e+00 9.51056516295153420e-01 1.110e-16
1.46607657167523664e+00 9.94519932551033659e-01 1.963e-06
1.67551608191455625e+00 9.94520472133828548e-01 1.423e-06
1.88495559215387587e+00 9.51056516295153642e-01 0.000e+00
2.09439510239319526e+00 8.66026266881472884e-01 8.631e-07
2.30383461263251466e+00 7.43145543942779097e-01 7.185e-07
2.51327412287183449e+00 5.87785252292473248e-01 0.000e+00
2.72271363311115389e+00 4.06736078624623798e-01 5.645e-07
2.93215314335047328e+00 2.07911159104951504e-01 5.317e-07
3.14159265358979312e+00 3.48786849800863176e-16 2.263e-16
3.35103216382911251e+00 -2.07911159104950893e-01 5.317e-07
3.56047167406843190e+00 -4.06736078624623076e-01 5.645e-07
3.76991118430775174e+00 -5.87785252292473248e-01 2.220e-16
3.97935069454707113e+00 -7.43145543942779541e-01 7.185e-07
4.18879020478639053e+00 -8.66026266881472884e-01 8.631e-07
4.39822971502571036e+00 -9.51056516295152754e-01 7.772e-16
4.60766922526502931e+00 -9.94520472133827882e-01 1.423e-06
4.81710873550434915e+00 -9.94519932551033659e-01 1.963e-06
5.02654824574366899e+00 -9.51056516295155085e-01 1.443e-15
5.23598775598298793e+00 -8.66029786417928715e-01 4.383e-06
5.44542726622230777e+00 -7.43151997300363854e-01 7.172e-06
5.65486677646162761e+00 -5.87785252292475802e-01 2.442e-15
5.86430628670094656e+00 -4.06712580221640141e-01 2.406e-05
6.07374579694026639e+00 -2.07860795644525476e-01 5.090e-05