/* nag_1d_cheb_fit (e02adc) Example Program. * * Copyright 1998 Numerical Algorithms Group. * * Mark 5, 1998. * Mark 8 revised, 2004. * */ #include #include #include #include int main(void) { #define A(I,J) a[(I)*tda + J] Integer exit_status=0, i, iwght, j, k, m, r, tda; NagError fail; double *a=0, *ak=0, d1, fit, *s=0, *w=0, *x=0, x1, xarg, xcapr, xm, *y=0; INIT_FAIL(fail); Vprintf("nag_1d_cheb_fit (e02adc) Example Program Results \n"); /* Skip heading in data file */ Vscanf("%*[^\n]"); while ((scanf("%ld",&m)) != EOF) { if (m >=2) { if ( !( x = NAG_ALLOC(m, double)) || !( y = NAG_ALLOC(m, double)) || !( w = NAG_ALLOC(m, double)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } } else { Vprintf("Invalid m.\n"); exit_status = 1; return exit_status; } Vscanf("%ld",&k); if (k>=1) { if ( !( a = NAG_ALLOC((k+1)*(k+1), double)) || !( s = NAG_ALLOC(k+1, double)) || !( ak = NAG_ALLOC(k+1, double)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } tda = k+1; } else { Vprintf("Invalid k.\n"); exit_status = 1; return exit_status; } Vscanf("%ld",&iwght); for (r = 0; r < m; ++r) { if (iwght != 1) { Vscanf("%lf",&x[r]); Vscanf("%lf",&y[r]); Vscanf("%lf",&w[r]); } else { Vscanf("%lf",&x[r]); Vscanf("%lf",&y[r]); w[r] = 1.0; } } /* nag_1d_cheb_fit (e02adc). * Computes the coefficients of a Chebyshev series * polynomial for arbitrary data */ nag_1d_cheb_fit(m, k+1, tda, x, y, w, a, s, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_1d_cheb_fit (e02adc).\n%s\n", fail.message); exit_status = 1; goto END; } for (i = 0; i <= k; ++i) { Vprintf("\n"); Vprintf(" %s%4ld%s%12.2e\n", "Degree", i, " R.M.S. residual =", s[i]); Vprintf("\n J Chebyshev coeff A(J) \n"); for (j = 0; j < i+1; ++j) Vprintf(" %3ld%15.4f\n", j+1, A(i,j)); } for (j = 0; j < k+1; ++j) ak[j] = A(k,j); x1 = x[0]; xm = x[m-1]; Vprintf("\n %s%4ld\n\n", "Polynomial approximation and residuals for degree",k); Vprintf(" R Abscissa Weight Ordinate Polynomial Residual \n"); for (r = 1; r <= m; ++r) { xcapr = (x[r-1] - x1 - (xm - x[r-1])) / (xm - x1); /* nag_1d_cheb_eval (e02aec). * Evaluates the coefficients of a Chebyshev series * polynomial */ nag_1d_cheb_eval(k+1, ak, xcapr, &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; } d1 = fit - y[r-1]; Vprintf(" %3ld%11.4f%11.4f%11.4f%11.4f%11.2e\n", r, x[r-1], w[r-1], y[r-1], fit, d1); if (r < m) { xarg = (x[r-1] + x[r]) * 0.5; xcapr = (xarg - x1 - (xm - xarg)) / (xm - x1); /* nag_1d_cheb_eval (e02aec), see above. */ nag_1d_cheb_eval(k+1, ak, xcapr, &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(" %11.4f %11.4f\n",xarg,fit); } } END: if (a) NAG_FREE(a); if (x) NAG_FREE(x); if (y) NAG_FREE(y); if (w) NAG_FREE(w); if (s) NAG_FREE(s); if (ak) NAG_FREE(ak); } return exit_status; }