/* nag_1d_cheb_interp_fit (e02afc) Example Program. * * Copyright 1998 Numerical Algorithms Group. * * Mark 5, 1998. * Mark 8 revised, 2004. * */ #include #include #include #include #include #include int main(void) { Integer exit_status = 0; double *an = 0, d1, *f = 0, fit, pi, piby2n, *xcap = 0; Integer i, j, n; Integer r; NagError fail; INIT_FAIL(fail); printf("nag_1d_cheb_interp_fit (e02afc) Example Program Results \n"); /* Skip heading in data file */ scanf("%*[^\n]"); /* nag_pi (x01aac). * pi */ pi = nag_pi; while ((scanf("%ld", &n)) != EOF) { if (n > 0) { if (!(an = NAG_ALLOC(n+1, double)) || !(f = NAG_ALLOC(n+1, double)) || !(xcap = NAG_ALLOC(n+1, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } } else { printf("Invalid n.\n"); exit_status = 1; return exit_status; } piby2n = pi * 0.5 / (double) n; for (r = 0; r < n+1; ++r) scanf("%lf", &f[r]); for (r = 0; r < n+1; ++r) { i = r; /* The following method of evaluating xcap = cos(pi*i/n) * ensures that the computed value has a small relative error * and, moreover, is bounded in modulus by unity for all * i = 0, 1, ..., n. (It is assumed that the sine routine * produces a result with a small relative error for values * of the argument between -PI/4 and PI/4). */ if (2*i <= n) { d1 = sin(piby2n * i); xcap[i] = 1.0 - d1 * d1 * 2.0; } else if (2*i > n * 3) { d1 = sin(piby2n * (n - i)); xcap[i] = d1 * d1 * 2.0 - 1.0; } else { xcap[i] = sin(piby2n * (n - 2*i)); } } /* nag_1d_cheb_interp_fit (e02afc). * Computes the coefficients of a Chebyshev series * polynomial for interpolated data */ nag_1d_cheb_interp_fit(n+1, f, an, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_1d_cheb_interp_fit (e02afc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\n Chebyshev \n"); printf(" J coefficient A(J) \n"); for (j = 0; j < n+1; ++j) printf(" %3ld%14.7f\n", j+1, an[j]); printf("\n R Abscissa Ordinate Fit \n"); for (r = 0; r < n+1; ++r) { /* nag_1d_cheb_eval (e02aec). * Evaluates the coefficients of a Chebyshev series * polynomial */ nag_1d_cheb_eval(n+1, an, xcap[r], &fit, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_1d_cheb_eval (e02aec).\n%s\n", fail.message); exit_status = 1; goto END; } printf(" %3ld%11.4f%11.4f%11.4f\n", r+1, xcap[r], f[r], fit); } END: NAG_FREE(an); NAG_FREE(f); NAG_FREE(xcap); } return exit_status; }