/* nag_quad_dim1_gauss_wrec (d01tdc) Example Program.
*
* Copyright 2020 Numerical Algorithms Group.
*
* Mark 27.1, 2020.
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
int main(void) {
/* Scalars */
Integer exit_status = 0;
Integer n, i;
double muzero, alpha, beta, a1, b1, ab1, ab2, d, ri, fa1, fb1, fab2;
/* Arrays */
char nag_enum_arg[40];
double *a = 0, *b = 0, *c = 0, *abscis = 0, *weight = 0;
const char *str_type;
/* Nag Types */
Nag_QuadType quadtype;
NagError fail, fail1, fail2;
INIT_FAIL(fail);
printf("nag_quad_dim1_gauss_wrec (d01tdc) Example Program Results\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Input number of abscissae required, n */
scanf("%" NAG_IFMT "%*[^\n] ", &n);
/* Input quadrature rule to simulate, Nag_QuadType */
scanf("%39s%*[^\n] ", nag_enum_arg);
quadtype = (Nag_QuadType)nag_enum_name_to_value(nag_enum_arg);
/* Allocate coefficient, weight and abscissae arrays */
if (!(a = NAG_ALLOC(n, double)) || !(b = NAG_ALLOC(n, double)) ||
!(c = NAG_ALLOC(n, double)) || !(weight = NAG_ALLOC(n, double)) ||
!(abscis = NAG_ALLOC(n, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Convert QuadType to string using
* nag_enum_value_to_name (x04nbc).
* Converts NAG enum member value to its name
*/
str_type = nag_enum_value_to_name(quadtype);
printf("\nQuadrature type = %s\n\n", str_type);
/* Generate recursion coefficients a, b, c from quadtype */
switch (quadtype) {
case Nag_Quad_Gauss_Legendre: {
muzero = 2.0;
for (i = 0; i < n; ++i) {
ri = (double)(i);
a[i] = (2.0 * ri + 1.0) / (ri + 1.0);
b[i] = 0.0;
c[i] = ri / (ri + 1.0);
}
} break;
case Nag_Quad_Gauss_Chebyshev_first: {
muzero = X01AAC;
for (i = 0; i < n; ++i) {
a[i] = 2.0;
b[i] = 0.0;
c[i] = 1.0;
}
a[0] = 1.0;
} break;
case Nag_Quad_Gauss_Chebyshev_second: {
muzero = 0.5 * X01AAC;
for (i = 0; i < n; ++i) {
a[i] = 2.0;
b[i] = 0.0;
c[i] = 1.0;
}
} break;
case Nag_Quad_Gauss_Jacobi: {
/* Input quadrature paramaters alpha and beta */
scanf("%lf%lf\n", &alpha, &beta);
printf("Using parameters alpha = %10.5f, beta = %10.5f\n\n", alpha, beta);
a1 = alpha + 1.0;
b1 = beta + 1.0;
ab1 = alpha + beta + 1.0;
ab2 = a1 + b1;
INIT_FAIL(fail1);
INIT_FAIL(fail2);
/* nag_specfun_gamma (s14aac).
* Gamma function Gamma(x)
*/
fa1 = nag_specfun_gamma(a1, &fail);
fb1 = nag_specfun_gamma(b1, &fail1);
fab2 = nag_specfun_gamma(ab2, &fail2);
if (fail.code != NE_NOERROR || fail1.code != NE_NOERROR ||
fail2.code != NE_NOERROR) {
printf("Error from nag_specfun_gamma (s14aac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
muzero = pow(2.0, ab1) * fa1 * fb1 / fab2;
a[0] = 0.5 * ab2;
b[0] = 0.5 * (alpha - beta);
c[0] = 0.0;
for (i = 1; i < n; ++i) {
ri = (double)i;
d = (2.0 * ri + 2.0) * (ri + ab1);
a[i] = (2.0 * ri + ab1) * (2.0 * ri + ab2) / d;
d = (2.0 * ri + alpha + beta) * d;
b[i] = (2.0 * ri + ab1) * (alpha * alpha - beta * beta) / d;
c[i] = 2.0 * (ri + alpha) * (ri + beta) * (2.0 * ri + ab2) / d;
}
} break;
case Nag_Quad_Gauss_Laguerre: {
/* Input quadrature paramaters alpha */
scanf("%lf\n", &alpha);
printf("Using parameter alpha = %10.5f\n\n", alpha);
a1 = alpha + 1.0;
/* nag_specfun_gamma (s14aac).
* gamma(x)
*/
muzero = nag_specfun_gamma(a1, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_specfun_gamma (s14aac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
for (i = 0; i < n; ++i) {
ri = (double)i;
a[i] = -1.0 / (ri + 1.0);
b[i] = (2.0 * ri + a1) / (ri + 1.0);
c[i] = (ri + alpha) / (ri + 1.0);
}
} break;
case Nag_Quad_Gauss_Hermite: {
muzero = sqrt(X01AAC);
for (i = 0; i < n; ++i) {
a[i] = 2.0;
b[i] = 0.0;
c[i] = (double)2 * i;
}
} break;
default: {
exit_status = 1;
printf("The quadrature type %s is not handled here\n", str_type);
goto END;
}
}
/* nag_quad_dim1_gauss_wrec (d01tdc).
* Compute weights and abscissae for a Gaussian quadrature rule
* governed by a three-term recurrence relation.
*/
nag_quad_dim1_gauss_wrec(n, a, b, c, muzero, weight, abscis, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_quad_dim1_gauss_wrec (d01tdc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
printf(" %3" NAG_IFMT " points\n\n Abscissae Weights\n\n", n);
for (i = 0; i < n; i++) {
printf("%15.6e", abscis[i]);
printf("%15.6e\n", weight[i]);
}
printf("\n");
END:
NAG_FREE(a);
NAG_FREE(b);
NAG_FREE(c);
NAG_FREE(abscis);
NAG_FREE(weight);
return exit_status;
}