/* nag_numdiff_1d_real_eval (d04bac) Example Program.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.1, 2017.
*/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd04.h>
#include <nags.h>
int main(void)
{
/* DER_COMP is used to store results for printing */
#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;
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. */
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 (%" NAG_IFMT ") calculated using "
"nag_real_polygamma (s14aec): %13.4e\n", j + 1, actder[j]);
printf(" Derivative and error estimates for derivative (%" NAG_IFMT ")\n",
j + 1);
printf(" hbase der[%" NAG_IFMT "] erest[%"
NAG_IFMT "]\n", j, j);
for (k = 0; k < n_hbase; k++)
printf("%15.4e %13.4e %13.1e\n", DER_COMP(k, 0, j), DER_COMP(k, 1, j),
DER_COMP(k, 2, j));
printf("\n");
}
END:
NAG_FREE(actder);
NAG_FREE(der_comp);
NAG_FREE(clabs);
NAG_FREE(clabsc);
return exit_status;
}