/***********************************************/ /* T.Kouya's GSL sample program collection */ /* Interpolation of Discrete Points */ /* Written by Tomonori Kouya */ /* */ /* Version 0.01: 2007-10-04 */ /***********************************************/ #include #include #include #include /* 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; }