/* 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); Vprintf("nag_1d_cheb_interp_fit (e02afc) Example Program Results \n"); /* Skip heading in data file */ Vscanf("%*[^\n]"); /* nag_pi (x01aac). * pi */ pi = X01AAC; 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)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } } else { Vprintf("Invalid n.\n"); exit_status = 1; return exit_status; } piby2n = pi * 0.5 / (double) n; for (r = 0; r < n+1; ++r) Vscanf("%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) { Vprintf("Error from nag_1d_cheb_interp_fit (e02afc).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("\n Chebyshev \n"); Vprintf(" J coefficient A(J) \n"); for (j = 0; j < n+1; ++j) Vprintf(" %3ld%14.7f\n",j+1,an[j]); Vprintf("\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) { Vprintf("Error from nag_1d_cheb_eval (e02aec).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf(" %3ld%11.4f%11.4f%11.4f\n",r+1,xcap[r],f[r],fit); } END: if (an) NAG_FREE(an); if (f) NAG_FREE(f); if (xcap) NAG_FREE(xcap); } return exit_status; }