f77 [ flags ] file(s) ... -L/usr/local/lib -lgjlC (K&R, 89, 99), C++ (98):SUBROUTINE gjqfd(x, w, deltaw, deltax, alpha, beta, nquad, ierr) DOUBLE PRECISION alpha, beta, deltax(*) DOUBLE PRECISION deltaw(*), w(*), x(*) INTEGER ierr, nquad
cc [ flags ] -I/usr/local/include file(s) ... -L/usr/local/lib -lgjl
Use#include <gjl.h>to get this prototype:void gjqfd(fortran_double_precision x_[], fortran_double_precision w_[], fortran_double_precision deltaw_[], fortran_double_precision deltax_[], const fortran_double_precision * alpha_, const fortran_double_precision * beta_, const fortran_integer * nquad_, fortran_integer * ierr_);
NB: The definition of C/C++ data types fortran_ xxx, and the mapping of Fortran external names to C/C++ external names, is handled by the C/C++ header file. That way, the same function or subroutine name can be used in C, C++, and Fortran code, independent of compiler conventions for mangling of external names in these programming languages.
Compute the nodes and weights for the evaluation of the integral
\int_{\(mi1}^{1} (1\(mix)^\alpha (1+x)^\beta \ln(1+x) f(x) dx (\alpha > \(mi1, \beta > \(mi1)
as the quadrature sum
\sum_{i=1}^{N}[\deltaW_i(\alpha,\beta) f(x_i(\alpha,\beta)) + \deltax_i(\alpha,\beta) f'(x_i(\alpha,\beta))]
The nonlogarithmic integral
\int_{\(mi1}^{1} (1\(mix)^\alpha (1+x)^\beta f(x) dx (\alpha > \(mi1, \beta > \(mi1)
can be computed from the quadrature sum
\sum_{i=1}^{N}W_i(\alpha, \beta) f(x_i(\alpha, \beta)).
The quadrature is exact to machine precision for f(x) of polynomial order less than or equal to 2*nquad \(mi 1.
This form of the quadrature requires values of the function and its derivative at N (== nquad) points. For a derivative-free quadrature at 2N points, see the companion routine, gjqf().
On entry:
- alpha
- Power of (1\(mix) in the integrand (alpha > \(mi1).
- beta
- Power of (1+x) in the integrand (beta > \(mi1).
- nquad
- Number of quadrature points to compute. It must be less than the limit MAXPTS defined in the header file, maxpts.inc. The default value chosen there should be large enough for any realistic application.
On return:
- x(1..nquad)
- Nodes of the Jacobi quadrature, denoted x_i(\alpha,\beta) above.
- w(1..nquad)
- Weights of the Jacobi quadrature, denoted W_i(\alpha,\beta) above.
- deltaw(1..nquad)
- Weights of the quadrature, denoted \deltaW_i(\alpha,\beta) above.
- deltax(1..nquad)
- Weights of the quadrature, denoted \deltax_i(\alpha,\beta) above.
- ierr
- Error indicator:
= 0 (success), 1 (eigensolution could not be obtained), 2 (destructive overflow), 3 (nquad out of range), 4 (alpha out of range), 5 (beta out of range).
The logarithmic integral can then be computed by code like this:
sum = 0.0d+00 do 10 i = 1,nquad sum = sum + deltaw(i)*f(x(i)) + deltax(i)*fprime(x(i)) 10 continue
where fprime(x(i)) is the derivative of the function f(x) with respect to x, evaluated at x = x(i).
The nonlogarithmic integral can be computed by:
sum = 0.0d+00 do 20 i = 1,nquad sum = sum + w(i)*f(x(i)) 20 continue
Fast Gaussian Quadrature for Two Classes of Logarithmic Weight Functionsin ACM Transactions on Mathematical Software, Volume ??, Number ??, Pages ????--???? and ????--????, 20xx, by
andNelson H. F. Beebe Center for Scientific Computing University of Utah Department of Mathematics, 110 LCB 155 S 1400 E RM 233 Salt Lake City, UT 84112-0090 Tel: +1 801 581 5254 FAX: +1 801 581 4148 Email: beebe@math.utah.edu, beebe@acm.org, beebe@computer.org WWW URL: http://www.math.utah.edu/~beebe
James S. Ball University of Utah Department of Physics Salt Lake City, UT 84112-0830 USA Tel: +1 801 581 8397 FAX: +1 801 581 6256 Email: ball@physics.utah.edu WWW URL: http://www.physics.utah.edu/people/faculty/ball.html