/* nag_ode_bvp_ps_lin_quad_weights (d02uyc) Example Program.
*
* Copyright 2024 Numerical Algorithms Group.
*
* Mark 30.2, 2024.
*/
#include <math.h>
#include <nag.h>
#ifdef __cplusplus
extern "C" {
#endif
static double NAG_CALL exact(double x);
#ifdef __cplusplus
}
#endif
int main(void) {
/* Scalars */
Integer exit_status = 0;
Integer i, n;
double a = -1.0, b = 3.0;
double integ, scale, uerr;
double teneps = 10.0 * nag_machine_precision;
/* Arrays */
double *f = 0, *w = 0, *x = 0;
/* NAG types */
Nag_Boolean reqerr = Nag_FALSE, reqwgt = Nag_FALSE;
NagError fail;
INIT_FAIL(fail);
printf("nag_ode_bvp_ps_lin_quad_weights (d02uyc) "
"Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%" NAG_IFMT "%*[^\n] ", &n);
if (!(f = NAG_ALLOC((n + 1), double)) || !(w = NAG_ALLOC((n + 1), double)) ||
!(x = NAG_ALLOC((n + 1), double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Set up solution grid:
* nag_ode_bvp_ps_lin_cgl_grid (d02ucc).
* Chebyshev Gauss-Lobatto grid generation.
*/
nag_ode_bvp_ps_lin_cgl_grid(n, a, b, x, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_ode_bvp_ps_lin_cgl_grid (d02ucc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Set up problem right hand sides for grid. */
for (i = 0; i < n + 1; i++)
f[i] = exact(x[i]);
scale = 0.5 * (b - a);
/* Solve on equally spaced grid:
* nag_ode_bvp_ps_lin_quad_weights (d02uyc).
* Clenshaw-Curtis quadrature weights for integration using computed
* Chebyshev coefficients.
*/
nag_ode_bvp_ps_lin_quad_weights(n, w, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_ode_bvp_ps_lin_quad_weights (d02uyc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Apply the weights, w, to the function values, f, and scale. */
integ = 0.0;
for (i = 0; i < n + 1; i++)
integ = integ + w[i] * f[i];
integ = scale * integ;
/* Print function values and weights if required. */
if (reqwgt) {
printf("f(x) and integral weights\n\n");
printf(" x f(x) w\n");
for (i = 0; i < n + 1; i++) {
printf("%10.4f %10.4f %10.4f\n", x[i], f[i], w[i]);
}
printf("\n");
}
/* Print approximation to integral. */
printf("Integral of f(x) from %6.1f to %6.1f = %13.5f\n", a, b, integ);
if (reqerr) {
uerr = fabs(integ - 28.0);
printf("Integral is within a multiple ");
printf("%8" NAG_IFMT " ", 10 * ((Integer)(uerr / teneps) + 1));
printf(" of machine precision.\n");
}
END:
NAG_FREE(f);
NAG_FREE(w);
NAG_FREE(x);
return exit_status;
}
static double NAG_CALL exact(double x) { return 3.0 * pow(x, 2); }