/* nag_quad_1d_gauss_wrec (d01tdc) Example Program.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.1, 2017.
*/
#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd01.h>
#include <nags.h>
#include <nagx01.h>
#include <nagx04.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_1d_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_gamma (s14aac).
* Gamma function Gamma(x)
*/
fa1 = nag_gamma(a1, &fail);
fb1 = nag_gamma(b1, &fail1);
fab2 = nag_gamma(ab2, &fail2);
if (fail.code != NE_NOERROR || fail1.code != NE_NOERROR ||
fail2.code != NE_NOERROR) {
printf("Error from nag_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_gamma (s14aac).
* gamma(x)
*/
muzero = nag_gamma(a1, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_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_1d_gauss_wrec (d01tdc).
* Compute weights and abscissae for a Gaussian quadrature rule
* governed by a three-term recurrence relation.
*/
nag_quad_1d_gauss_wrec(n, a, b, c, muzero, weight, abscis, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_quad_1d_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;
}