/* nag_numdiff_1d_real_eval (d04bac) Example Program. * * Copyright 2011, Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #include #include #include int main(void) { /* DER_COMP is used to store results in column major format only. */ #define DER_COMP(I, J, K) der_comp[I + J*n_hbase + K*n_hbase*n_der_comp] Integer exit_status = 0; double x_0; Integer n_der_comp, n_display, n_hbase; double hbase; Integer i, j, k; char title[51]; double der[14], erest[14], fval[21], xval[21]; double *actder = 0, *der_comp = 0; char *clabs = 0, **clabsc = 0; NagError fail; INIT_FAIL(fail); printf("nag_numdiff_1d_real_eval (d04bac) Example Program Results\n"); printf("\n" "Find the derivatives of the polygamma (psi) function\n" "using function values generated by nag_real_polygamma (s14aec).\n\n" "Demonstrate the effect of successively reducing hbase.\n\n"); fflush(stdout); n_der_comp = 3; n_display = 3; n_hbase = 4; if (!(clabs = NAG_ALLOC(n_der_comp*10, char)) || !(clabsc = NAG_ALLOC(n_der_comp, char *)) || !(actder = NAG_ALLOC(n_display, double)) || !(der_comp = NAG_ALLOC(n_hbase*n_der_comp*14, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } x_0 = 0.05; /* Select an initial separation distance hbase. */ hbase = 0.0025; /* Compute the actual derivatives at target location x_0 * using nag_real_polygamma (s14aec), for comparison. */ for (j = 0; j < n_display; j++) { /* nag_real_polygamma (s14aec). * Derivative of the psi function psi(x). */ actder[j] = nag_real_polygamma(x_0, j+1, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_real_polygamma (s14aec).\n%s\n", fail.message); exit_status = 1; goto END; } } /* Attempt n_hbase approximations, reducing hbase by factor 0.1 each time. */ for (j = 0; j < n_hbase; j++) { /* nag_numdiff_1d_real_absci (d04bbc). * Generates sample points for function evaluations * by nag_numdiff_1d_real_eval (d04bac). */ nag_numdiff_1d_real_absci(x_0, hbase, xval); /* Calculate the corresponding objective function values. */ for (k = 0; k < 21; k++) { fval[k] = nag_real_polygamma(xval[k], 0, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_real_polygamma (s14aec).\n%s\n", fail.message); exit_status = 1; goto END; } } /* Calculate the derivative estimates. */ /* nag_numdiff_1d_real_eval (d04bac). * Numerical differentiation, user-supplied function values, * derivatives up to order 14, * derivatives with respect to one real variable. */ nag_numdiff_1d_real_eval(xval, fval, der, erest, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_numdiff_1d_real_eval (d04bac).\n%s\n", fail.message); exit_status = 1; goto END; } /* Store results in DER_COMP in column major format. */ for (i = 0; i < 14; i++) { DER_COMP(j, 0, i) = hbase; DER_COMP(j, 1, i) = der[i]; DER_COMP(j, 2, i) = erest[i]; } /* Decrease hbase for next loop. */ hbase = hbase * 0.1; } /* Display results for first n_display derivatives. */ for (j = 0; j < n_display; j++) { printf(" Derivative (%ld) calculated using " "nag_real_polygamma (s14aec): %13.4e\n", j+1, actder[j]); fflush(stdout); sprintf(title, "Derivative and error estimates for derivative (%ld)", j+1); sprintf(&clabs[0], "hbase"); sprintf(&clabs[10], "der[%ld]", j); sprintf(&clabs[20], "erest[%ld]", j); for (k = 0; k < n_der_comp; k++) clabsc[k] = &clabs[k*10]; /* nag_gen_real_mat_print_comp (x04cbc). * Print real general matrix (comprehensive). */ nag_gen_real_mat_print_comp(Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, n_hbase, n_der_comp, &DER_COMP(0, 0, j), n_hbase, NULL, title, Nag_NoLabels, NULL, Nag_CharacterLabels, (const char **) clabsc, 0, 0, NULL, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\n"); } END: NAG_FREE(actder); NAG_FREE(der_comp); NAG_FREE(clabs); NAG_FREE(clabsc); return exit_status; }