Package gov.nih.mipav.model.algorithms

Class Integration2

java.lang.Object
gov.nih.mipav.model.algorithms.Integration2
Direct Known Subclasses:
AlgorithmCircleGeneration.erfcModel2, AlgorithmCircleGeneration.IntTorquato95ModelMean2, AlgorithmFRAP.FitFullIntModel2i, AlgorithmFRAP.FitFullIntModel2p, AlgorithmFRAP.FitFullIntModel2s, AlgorithmFRAP.IntModel2, AlgorithmFRAP.IntModelBessel, AlgorithmFRAP.IntModelI0NuclearArea, AlgorithmSM2.Integration2All, AlgorithmSphereGeneration.IntModelMean2, AlgorithmSphereGeneration.IntModelMeanSquared2, AlgorithmSphereGeneration.IntTorquato95ModelMean2, AlgorithmSphereGeneration.IntTorquatoModelMean2, SymmsIntegralMapping.qaws
public abstract class Integration2 extends Object
This is a port of FORTRAN numerical integration routines in QUADPACK found at http://www.netlib.org/quadpack Reference: R. Piessens, E. deDoncker-Kapenga, C. Uberhuber, D. Kahaner Quadpack: a Subroutine Package for Automatic Integration Springer Verlag, 1983. Series in Computational Mathematics v. 1 The original dqage, dqagie, dqagpe, dqagse, dqawce, dqawfe, dqawoe, adn dqawse routines were written by Robert Piessens and Elise de Doncker. The original dqng routine was written by Robert Piessens and Elise de Doncker and modified by David Kahaner. The original dqc25c, dqc25f, dqc25s, dqcheb, dqelg, dqk15, dqk15i, dqk15w, dqk21, dqk31, dqk41, dqk51, dqk61, dqmomo, dqpsrt, dqwgtc, dqwgtf, and dqwgts routines were written by Robert Piessens and Elise de Doncker. The simple contents of dqwgtc, dqwgtf, and dqwgts are incorporated into dqk15w. The linpack routine dgtsl was written by Jack Dongarra of the Argonne National Laboratory. Self tests 1 to 15 come from quadpack_prb.f90 by John Burkardt. Porting was performed by William Gandler.
The QUADPACK code is made freely available, with the following restrictions: 
 Re-distributions of the source code should credit Quadpack and retain the notice with the names of the original authors;
 the names of the copyright holders or contributors may not be used to endorse or promote any derived software without prior permission;
 the Quadpack software is provided "as is", without warranties;
 Quadpack and the authors deny any liability for situations resulting from the use of this software.

  • Field Summary

    Fields
    Modifier and Type
    Field
    Description
    private final double[]
    abserr
    Estimate of the absolute value of the error, which should equal or exceed abs(actual integral - result).
    private double
    actualAnswer
     
    private double
    alfa
    parameter in the weight function, alfa.gt.(-1) if alfa.le.(-1), the routine will end with ier[0] = 6.
    private double[]
    alist
    The left end points of the subintervals in the partition of the given integration range (lower, upper).
    private double
    beta
    parameter in the weight function, beta.gt.(-1) if beta.le.(-1), the routine will end with ier[0] = 6.
    private double[]
    blist
    The right end points of the subintervals in the partition of the given integration range (lower, upper).
    private double
    bound
    finite bound of integration range used in dqagie (has no meaning if interval is doubly-infinite).
    private double[]
    breakPoints
    Used in dqagpe This array must always be >= 2 in length.
    private double
    c
    Parameter in the weight function, c !
    private double[][]
    chebmo
    array of dimension (maxp1,25) containing the chebyshev moments
    protected static final int
    DQAGE
    Integral without singularities, discontinuities, or infinite bounds.
    protected static final int
    DQAGIE
    Integral over (bound, +infinity), (-infinity, bound), or (-infinity, +infinity).
    protected static final int
    DQAGPE
    Integral over (lower, upper) with user specified break points.
    protected static final int
    DQAGSE
    Integral over (lower, upper) which can handle end point singularities.
    protected static final int
    DQAWCE
    Cauchy principal value
    protected static final int
    DQAWFE
    Automatic integrator, special-purpose, fourier integrals
    protected static final int
    DQAWOE
    Automatic integrator, integrand with oscillatory cosine or sine factor
    protected static final int
    DQAWSE
    Adaptive integrator for integrands with algebraic or logarithmic singularities at endpoints.
    protected static final int
    DQK15
    Gauss-Kronod quadrature rule
    protected static final int
    DQK21
    Gauss-Kronod quadrature rule
    protected static final int
    DQK31
    Gauss-Kronod quadrature rule
    protected static final int
    DQK41
    Gauss-Kronod quadrature rule
    protected static final int
    DQK51
    Gauss-Kronod quadrature rule
    protected static final int
    DQK61
    Gauss-Kronod quadrature rule
    protected static final int
    DQNG
    Non-adaptive integration.
    private double[]
    elist
    Estimates of the absolute errors on the subintervals.
    private double
    epmach
    epmach = D1MACH(4) Machine epmach is the smallest positive epmach such that (1.0 + epmach) !
    private double
    epsabs
    If epsabs invalid input: '<'= 0.0 and epsrel invalid input: '<' Math.max(50*relative machine accuracy, 5.0E-29), the routine will exit with an error message.
    private double
    epsrel
    DOCUMENT ME!
    private double[]
    erlst
    In dqawfe vector of dimension at least limlst erlst[k] contains the error estimate corresponding wiht rstlst[k].
    private double
    gamma
     
    private int
    icall
    if dqawoe is to be used only once, icall must be set to 1. assume that during this call, the c chebyshev moments (for clenshaw-curtis integration c of degree 24) have been computed for intervals of c lenghts (abs(b-a))*2**(-l), l=0,1,2,...momcom-1. c if icall.gt.1 this means that dqawoe has been c called twice or more on intervals of the same c length abs(b-a). the chebyshev moments already c computed are then re-used in subsequent calls. c if icall.lt.1, the routine will end with ier[0] = 6.
    private final int[]
    ier
    ier[0] = 0 for normal and reliable termination of the routine.
    private int[]
    ierlst
    In dqawfe vector of dimension at least limlst ierlst[k] contains the error flag corresponding with rstlst[k].
    private int
    inf
    In dqagie indicates the kind of integration range involved inf = 1 corresponds to (bound, +infinity) inf = -1 corresponds to (-infinity, bound) inf = 2 corresponds to (-infinity, +infinity).
    private int
    integr
    With dqawoe: indicates which of the weight functions is to be used integr = 1 w(x) = cos(omega*x) integr = 2 w(x) = sin(omega*x) if integr.ne.1 and integr.ne.2, the routine will end with ier[0] = 6.
    private int[]
    iord
    The first k elements of iord are pointers to the error estimates over the subintervals, such that elist[iord[0]], ..., elist[iord[k-1]] form a decreasing sequence, with k = last if last invalid input: '<'= (limit/2 + 2) , and k = limit + 1 - last otheriwse.
    private int
    key
    In dqage key selects the local integration rule A gauss-kronrod pair is used with 7 - 15 points if key invalid input: '<' 2 10 - 21 points if key = 2 15 - 31 points if key = 3 20 - 41 points if key = 4 25 - 51 points if key = 5 30 - 61 points if key > 5.
    private int
    last
    The number of subintervals actually called in the subdivision process.
    private int[]
    level
    level contains the subdivision levels of the subinterval, i.e., if (aa, bb) is a subinterval of (p1, p2) where p1 as well as p2 is a user-provided break point or integration limit, then (aa, bb) has level 1 if abs(bb - aa) = abs(p2 - p1) * 2**(-1).
    private int
    limit
    Gives an upper bound on the number of subintervals in the partition of lower, upper.
    private int
    limlst
    In dqawfe limlst gives an upper bound on the number of cycles, limlst >= 1.
    private double
    lower
    DOCUMENT ME!
    private int
    lst
    In dqawfe number of subintervals needed for the integration.
    private int
    maxp1
    gives an upper bound on the number of chebyshev c moments which can be stored, i.e. for the c intervals of lenghts abs(b-a)*2**(-l), c l=0,1, ..., maxp1-2, maxp1.ge.1. c if maxp1.lt.1, the routine will end with ier[0] = 6.
    private int
    momcom
    indicating that the chebyshev moments c have been computed for intervals of lengths c (abs(b-a))*2**(-l), l=0,1,2, ..., momcom-1, c momcom.lt.maxp1
    private int[]
    ndin
    After the first integration over the intervals (pts[i], pts[i+1]), i = 0, 1, ..., npts2 - 2, the error estimates over some of the intervals may have been increased artificially, in order to put their subdivision forward.
    private final int[]
    neval
    Number of integrand evaluations.
    private int[]
    nnlog
    vector of dimension at least limit, containing the c subdivision levels of the subintervals, i.e. c iwork(i) = l means that the subinterval c numbered i is of length abs(b-a)*2**(1-l)
    private int
    npts
    npts = npts2 - 2.
    private int
    npts2
    npts2 = breakPoints.length.
    private final double
    oflow
    D1MACH(2) = 2**(emax)*(1 - 2**(-doubleDigits)) = 2**1024*(1 - 2**-53) D1MACH(2) = Double.MAX_VALUE.
    private double
    omega
    parameter in the integrand weight function
    private double[]
    pts
    Integration limits and the break points of the interval in ascending sequence.
    private final double[]
    resabs
     
    private final double[]
    resasc
     
    private final double[]
    result
    Approximation to the integral.
    private double[]
    rlist
    The integral approximations of the subintervals.
    private int
    routine
    DOCUMENT ME!
    private double[]
    rslst
    In dqawfe vector of dimension at least limlst. rslst[k-1] contains the integral contribution over the interval (a+(k-1)c, a+kc), where c = (2*int(abs(omega))+1)*pi/abs(omega), k = 1, 2, ..., lst.
    private boolean
    selfTest
     
    private int
    testCase
     
    private final double
    uflow
    DOCUMENT ME!
    private double
    upper
    DOCUMENT ME!

  • Constructor Summary

    Constructors
    Constructor
    Description
    Integration2()
    Constructor for running self tests.
    Integration2(double lower, double upper, double c, int routine, double epsabs, double epsrel, int limit)
    Constructor for Integration2 Used with routine = DQAWCE.
    Integration2(double lower, double upper, int routine, double[] breakPoints, double epsabs, double epsrel, int limit)
    Constructor for Integration Used with routine = DQAGPE.
    Integration2(double lower, double upper, int routine, double epsabs, double epsrel)
    Constructor for Integration Used with routine = DQNG.
    Integration2(double lower, double upper, int routine, double epsabs, double epsrel, int limit)
    Constructor for Integration Used with routine = DQAGSE.
    Integration2(double lower, double upper, int routine, double alfa, double beta, int integr, double epsabs, double epsrel, int limit)
    This constructor must be used with routine DQAWSE.
    Integration2(double lower, double upper, int routine, double omega, int integr, double epsabs, double epsrel, int limit, int icall, int maxp1)
    Constructor must be used with routine dqawoe
    Integration2(double lower, double upper, int routine, int key, double epsabs, double epsrel, int limit)
    Constructor for Integration Used with routine = DQAGE.
    Integration2(double lower, int routine, double omega, int integr, double epsabs, int limlst, int limit, int maxp1)
    Constructor must be used with routine dqawfe
    Integration2(double bound, int routine, int inf, double epsabs, double epsrel, int limit)
    Constructor used with routine = DQAGIE.

  • Method Summary

    Modifier and Type
    Method
    Description
    private void
    dgtsl(int n, double[] c, double[] d, double[] e, double[] b, int[] info)
    This is a port of the FORTRAN routine whose header is given below: c c dgtsl given a general tridiagonal matrix and a right hand c side will find the solution. c c on entry c c n integer c is the order of the tridiagonal matrix. c c c double precision(n) c is the subdiagonal of the tridiagonal matrix. c c(2) through c(n) should contain the subdiagonal. c on output c is destroyed. c c d double precision(n) c is the diagonal of the tridiagonal matrix. c on output d is destroyed. c c e double precision(n) c is the superdiagonal of the tridiagonal matrix. c e(1) through e(n-1) should contain the superdiagonal. c on output e is destroyed. c c b double precision(n) c is the right hand side vector. c c on return c c b is the solution vector. c c info integer c = 0 normal value. c = k if the k-th element of the diagonal becomes c exactly zero. the subroutine returns when c this is detected. c c linpack. this version dated 08/14/78 . c jack dongarra, argonne national laboratory.
    private void
    dqage()
    This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqage c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a1 c***keywords automatic integrator, general-purpose, c integrand examinator, globally adaptive, c gauss-kronrod c***author piessens,robert,appl. math.
    void
    dqagie(double bound, int inf, double epsrel)
    This is a port of the orginal FORTRAN code whose header is given below: c***begin prologue dqagie c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a3a1,h2a4a1 c***keywords automatic integrator, infinite intervals, c general-purpose, transformation, extrapolation, c globally adaptive c***author piessens,robert,appl. math invalid input: '&' progr. div - k.u.leuven c de doncker,elise,appl. math invalid input: '&' progr. div - k.u.leuven c***purpose the routine calculates an approximation result to a given c integral i = integral of f over (bound,+infinity) c or i = integral of f over (-infinity,bound) c or i = integral of f over (-infinity,+infinity), c hopefully satisfying following claim for accuracy c abs(i-result).le.max(epsabs,epsrel*abs(i)) c***description c c integration over infinite intervals c standard fortran subroutine
    void
    dqagpe()
    dqagpe.
    private void
    dqagse()
    This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqagse c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a1 c***keywords automatic integrator, general-purpose, c (end point) singularities, extrapolation, c globally adaptive c***author piessens,robert,appl. math.
    private void
    dqawce()
    This is the port of the original FORTRAN routine whose header is given below: c***begin prologue dqawce c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a1,j4 c***keywords automatic integrator, special-purpose, c cauchy principal value, clenshaw-curtis method c***author piessens,robert,appl.
    private void
    dqawfe()
    This is the port of the original FORTRAN routine whose header is given below: c***begin prologue dqawfe c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a3a1 c***keywords automatic integrator, special-purpose, c fourier integrals, c integration between zeros with dqawoe, c convergence acceleration with dqelg c***author piessens,robert,appl. math.
    private void
    dqawoe(double a, double b, double epsabs, double epsrel, int icall, double[] result, double[] abserr, int[] neval, int[] ier)
    This is the port of the original FORTRAN routine whose header is given below: c***begin prologue dqawoe c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a1 c***keywords automatic integrator, special-purpose, c integrand with oscillatory cos or sin factor, c clenshaw-curtis method, (end point) singularities, c extrapolation, globally adaptive c***author piessens,robert,appl. math.
    private void
    dqawse()
    This is the port of the original FORTRAN routine whose header is given below: c***begin prologue dqawse c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a1 c***keywords automatic integrator, special-purpose, c algebraico-logarithmic end point singularities, c clenshaw-curtis method c***author piessens,robert,appl. math.
    private void
    dqc25c(double a, double b, double c, double[] result, double[] abserr, int[] krul, int[] neval)
    This is the port of the original FORTRAN routine whose header is given below: c***begin prologue dqc25c c***date written 810101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a2,j4 c***keywords 25-point clenshaw-curtis integration c***author piessens,robert,appl. math.
    private void
    dqc25f(double a, double b, double omega, int nrmom, int ksave, double[] result, double[] abserr, int[] neval, double[] resabs, double[] resasc)
    This is a port of the original FORTRAN whose header is given below: c***begin prologue dqc25f c***date written 810101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a2 c***keywords integration rules for functions with cos or sin c factor, clenshaw-curtis, gauss-kronrod c***author piessens,robert,appl. math.
    private void
    dqc25s(double bl, double br, double[] ri, double[] rj, double[] rg, double[] rh, double[] result, double[] abserr, double[] resasc, int[] nev)
    This is a port of the original FORTRAN whose header is given below: c***begin prologue dqc25s c***date written 810101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a2 c***keywords 25-point clenshaw-curtis integration c***author piessens,robert,appl. math.
    private void
    dqcheb(double[] x, double[] fval, double[] cheb12, double[] cheb24)
    This is a port of the original FORTRAN whose header is given below: c***begin prologue dqcheb c***refer to dqc25c,dqc25f,dqc25s c***routines called (none) c***revision date 830518 (yymmdd) c***keywords chebyshev series expansion, fast fourier transform c***author piessens,robert,appl. math.
    private void
    dqelg(int[] n, double[] epstab, double[] result, double[] abserr, double[] res31a, int[] nres)
    This is a port of the original FORTRAN whose header is given below: c***begin prologue dqelg c***refer to dqagie,dqagoe,dqagpe,dqagse c***routines called d1mach c***revision date 830518 (yymmdd) c***keywords epsilon algorithm, convergence acceleration, c extrapolation c***author piessens,robert,appl. math.
    private void
    dqk15(double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
    The is a port of the original FORTRAN whose header is given below: c***begin prologue dqk15 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 15-point gauss-kronrod rules c***author piessens,robert,appl. math.
    private void
    dqk15i(double boun, int inf, double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
    This is a port of the original FORTRAN whose header is given below: c***begin prologue dqk15i c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a3a2,h2a4a2 c***keywords 15-point transformed gauss-kronrod rules c***author piessens,robert,appl. math.
    private void
    dqk15w(double p1, double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
    This a a port of the original FORTRAN whose header is given below: c***begin prologue dqk15w c***date written 810101 (yymmdd) c***revision date 830518 (mmddyy) c***category no. h2a2a2 c***keywords 15-point gauss-kronrod rules c***author piessens,robert,appl. math.
    private void
    dqk21(double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
    This is a port of the original FORTRAN whose header is given below: c***begin prologue dqk21 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 21-point gauss-kronrod rules c***author piessens,robert,appl. math.
    private void
    dqk31(double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
    This is a port of the original FORTRAN whose header is given below: c***begin prologue dqk31 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 31-point gauss-kronrod rules c***author piessens,robert,appl. math.
    private void
    dqk41(double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
    This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqk41 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 41-point gauss-kronrod rules c***author piessens,robert,appl. math.
    private void
    dqk51(double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
    This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqk51 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 51-point gauss-kronrod rules c***author piessens,robert,appl. math.
    private void
    dqk61(double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
    This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqk61 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 61-point gauss-kronrod rules c***author piessens,robert,appl. math.
    private void
    dqmomo(double[] ri, double[] rj, double[] rg, double[] rh)
    This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqmomo c***date written 820101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a1,c3a2 c***keywords modified chebyshev moments c***author piessens,robert,appl. math.
    private void
    dqng()
    This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqng c***date written 800101 (yymmdd) c***revision date 810101 (yymmdd) c***category no. h2a1a1 c***keywords automatic integrator, smooth integrand, c non-adaptive, gauss-kronrod(patterson) c***author piessens,robert,appl. math.
    private void
    dqpsrt(int limit, int last, int[] maxErr, double[] ermax, double[] elist, int[] iord, int[] nrmax)
    This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqpsrt c***refer to dqage,dqagie,dqagpe,dqawse c***routines called (none) c***revision date 810101 (yymmdd) c***keywords sequential sorting c***author piessens,robert,appl. math.
    void
    driver()
    DOCUMENT ME!
    double
    getAbserr()
    DOCUMENT ME!
    int
    getErrorStatus()
    DOCUMENT ME!
    double
    getIntegral()
    DOCUMENT ME!
    int
    getNeval()
    DOCUMENT ME!
    abstract double
    intFunc(double x)
    DOCUMENT ME!
    private double
    intFuncTest(double x)
     

    Methods inherited from class java.lang.Object

    clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

  • Field Details

    • DQAGPE

      protected static final int DQAGPE
      Integral over (lower, upper) with user specified break points.
      See Also:

    • DQAGIE

      protected static final int DQAGIE
      Integral over (bound, +infinity), (-infinity, bound), or (-infinity, +infinity).
      See Also:

    • DQAGSE

      protected static final int DQAGSE
      Integral over (lower, upper) which can handle end point singularities.
      See Also:

    • DQAGE

      protected static final int DQAGE
      Integral without singularities, discontinuities, or infinite bounds.
      See Also:

    • DQNG

      protected static final int DQNG
      Non-adaptive integration.
      See Also:

    • DQK15

      protected static final int DQK15
      Gauss-Kronod quadrature rule
      See Also:

    • DQK21

      protected static final int DQK21
      Gauss-Kronod quadrature rule
      See Also:

    • DQK31

      protected static final int DQK31
      Gauss-Kronod quadrature rule
      See Also:

    • DQK41

      protected static final int DQK41
      Gauss-Kronod quadrature rule
      See Also:

    • DQK51

      protected static final int DQK51
      Gauss-Kronod quadrature rule
      See Also:

    • DQK61

      protected static final int DQK61
      Gauss-Kronod quadrature rule
      See Also:

    • DQAWCE

      protected static final int DQAWCE
      Cauchy principal value
      See Also:

    • DQAWOE

      protected static final int DQAWOE
      Automatic integrator, integrand with oscillatory cosine or sine factor
      See Also:

    • DQAWSE

      protected static final int DQAWSE
      Adaptive integrator for integrands with algebraic or logarithmic singularities at endpoints.
      See Also:

    • DQAWFE

      protected static final int DQAWFE
      Automatic integrator, special-purpose, fourier integrals
      See Also:

    • abserr

      private final double[] abserr
      Estimate of the absolute value of the error, which should equal or exceed abs(actual integral - result).

    • alfa

      private double alfa
      parameter in the weight function, alfa.gt.(-1) if alfa.le.(-1), the routine will end with ier[0] = 6.

    • alist

      private double[] alist
      The left end points of the subintervals in the partition of the given integration range (lower, upper).

    • beta

      private double beta
      parameter in the weight function, beta.gt.(-1) if beta.le.(-1), the routine will end with ier[0] = 6.

    • blist

      private double[] blist
      The right end points of the subintervals in the partition of the given integration range (lower, upper).

    • bound

      private double bound
      finite bound of integration range used in dqagie (has no meaning if interval is doubly-infinite).

    • breakPoints

      private double[] breakPoints
      Used in dqagpe This array must always be >= 2 in length. The first breakPoints.length - 2 points are user provided break points. If these points are not in an ascending sequence, they will automatically be sorted.

    • c

      private double c
      Parameter in the weight function, c != lower, c != upper. If c == a or c == b, the routine will end with errorStatus = 6

    • chebmo

      private double[][] chebmo
      array of dimension (maxp1,25) containing the chebyshev moments

    • elist

      private double[] elist
      Estimates of the absolute errors on the subintervals.

    • epmach

      private double epmach
      epmach = D1MACH(4) Machine epmach is the smallest positive epmach such that (1.0 + epmach) != 1.0. epmach = 2**(1 - doubleDigits) = 2**(1 - 53) = 2**(-52) epmach = 2.224460e-16 epmach is called the largest relative spacing

    • epsabs

      private double epsabs
      If epsabs invalid input: '<'= 0.0 and epsrel invalid input: '<' Math.max(50*relative machine accuracy, 5.0E-29), the routine will exit with an error message.

    • epsrel

      private double epsrel
      DOCUMENT ME!

    • erlst

      private double[] erlst
      In dqawfe vector of dimension at least limlst erlst[k] contains the error estimate corresponding wiht rstlst[k].

    • gamma

      private double gamma

    • icall

      private int icall
      if dqawoe is to be used only once, icall must be set to 1. assume that during this call, the c chebyshev moments (for clenshaw-curtis integration c of degree 24) have been computed for intervals of c lenghts (abs(b-a))*2**(-l), l=0,1,2,...momcom-1. c if icall.gt.1 this means that dqawoe has been c called twice or more on intervals of the same c length abs(b-a). the chebyshev moments already c computed are then re-used in subsequent calls. c if icall.lt.1, the routine will end with ier[0] = 6.

    • ier

      private final int[] ier
      ier[0] = 0 for normal and reliable termination of the routine. It is assumed that the requested accuracy has been achieved. ier[0] > 0 for abnormal termination of the routine. The estimates for the integral and error are less reliable. It is assumed that the requested accuracy has not been achieved. ier[0] = 1 maximum number of subdivisions allowed has been achieved. One can allow more subdivisions by increasing the value of limit (and taking the according dimension adjustments into account). However, if this yields no improvement, it is advised to analyze the integrand in order to determine the integration difficulties. If the position of a local difficulty can be determined( i.e. singularity, discontinuity within the interval), it should be supplied to the routine as an element of the breakPoints array. If necessary, an appropriate special purpose integrator must be used, which is designed for handling the type of difficulty involved. ier[0] = 2 The occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be under-estimated. ier[0] = 3 Extremely bad integrand behavior occurs at some points of the integration interval. ier[0] = 4 The algorithm does not converge. Roundoff error is detected in the extrapolation table. It is presumed that the requested tolerance cannot be achieved, and that the returned result is the best that can be obtained. ier[0] = 5 The integral is probably divergent, or slowly convergent. It must be noted that divergence can occur with any other value of ier[0] > 0. ier[0] = 6 The input is invalid because (epsabs invalid input: '<'= 0 and epsrel invalid input: '<' max(50.0 * relative machine accuracy, 5.0E-29) or in the case of dqagpe can also happen because npts2 invalid input: '<' 2, the breakPoints are specified outside the integration range, or , or limit invalid input: '<'= npts.

    • ierlst

      private int[] ierlst
      In dqawfe vector of dimension at least limlst ierlst[k] contains the error flag corresponding with rstlst[k]. For the meaning of the local error flags see description of the output parameter ier.

    • inf

      private int inf
      In dqagie indicates the kind of integration range involved inf = 1 corresponds to (bound, +infinity) inf = -1 corresponds to (-infinity, bound) inf = 2 corresponds to (-infinity, +infinity).

    • integr

      private int integr
      With dqawoe: indicates which of the weight functions is to be used integr = 1 w(x) = cos(omega*x) integr = 2 w(x) = sin(omega*x) if integr.ne.1 and integr.ne.2, the routine will end with ier[0] = 6. With dqawse: indicates which weight function is to be used = 1 (x-lower)**alfa*(upper-x)**beta = 2 (x-lower)**alfa*(upper-x)**beta*log(x-lower) = 3 (x-lower)**alfa*(upper-x)**beta*log(upper-x) = 4 (x-lower)**alfa*(upper-x)**beta*log(x-lower)*log(upper-x) if integr.lt.1 or integr.gt.4, the routine will end with ier[0] = 6.

    • iord

      private int[] iord
      The first k elements of iord are pointers to the error estimates over the subintervals, such that elist[iord[0]], ..., elist[iord[k-1]] form a decreasing sequence, with k = last if last invalid input: '<'= (limit/2 + 2) , and k = limit + 1 - last otheriwse.

    • key

      private int key
      In dqage key selects the local integration rule A gauss-kronrod pair is used with 7 - 15 points if key invalid input: '<' 2 10 - 21 points if key = 2 15 - 31 points if key = 3 20 - 41 points if key = 4 25 - 51 points if key = 5 30 - 61 points if key > 5.

    • last

      private int last
      The number of subintervals actually called in the subdivision process.

    • level

      private int[] level
      level contains the subdivision levels of the subinterval, i.e., if (aa, bb) is a subinterval of (p1, p2) where p1 as well as p2 is a user-provided break point or integration limit, then (aa, bb) has level 1 if abs(bb - aa) = abs(p2 - p1) * 2**(-1).

    • limit

      private int limit
      Gives an upper bound on the number of subintervals in the partition of lower, upper. In dqagpe limit must be >= npts2. If limit invalid input: '<' npts2 in dqagpe, the routine exits with an error message. In dqawse limit must be >= 2. If limit invalid input: '<' 2 in dqawse, the routine exits with ier[0] = 6.

    • limlst

      private int limlst
      In dqawfe limlst gives an upper bound on the number of cycles, limlst >= 1. If limlst invalid input: '<' 3, the routine will end with ier[0] = 6.

    • lower

      private double lower
      DOCUMENT ME!

    • lst

      private int lst
      In dqawfe number of subintervals needed for the integration. If omega = 0 then lst is set to 1.

    • maxp1

      private int maxp1
      gives an upper bound on the number of chebyshev c moments which can be stored, i.e. for the c intervals of lenghts abs(b-a)*2**(-l), c l=0,1, ..., maxp1-2, maxp1.ge.1. c if maxp1.lt.1, the routine will end with ier[0] = 6.

    • momcom

      private int momcom
      indicating that the chebyshev moments c have been computed for intervals of lengths c (abs(b-a))*2**(-l), l=0,1,2, ..., momcom-1, c momcom.lt.maxp1

    • ndin

      private int[] ndin
      After the first integration over the intervals (pts[i], pts[i+1]), i = 0, 1, ..., npts2 - 2, the error estimates over some of the intervals may have been increased artificially, in order to put their subdivision forward. If this happens for the subinterval numbered k, ndin[k] is put to 1, otherwise ndin[k] = 0.

    • neval

      private final int[] neval
      Number of integrand evaluations.

    • nnlog

      private int[] nnlog
      vector of dimension at least limit, containing the c subdivision levels of the subintervals, i.e. c iwork(i) = l means that the subinterval c numbered i is of length abs(b-a)*2**(1-l)

    • npts

      private int npts
      npts = npts2 - 2.

    • npts2

      private int npts2
      npts2 = breakPoints.length.

    • oflow

      private final double oflow
      D1MACH(2) = 2**(emax)*(1 - 2**(-doubleDigits)) = 2**1024*(1 - 2**-53) D1MACH(2) = Double.MAX_VALUE.
      See Also:

    • omega

      private double omega
      parameter in the integrand weight function

    • pts

      private double[] pts
      Integration limits and the break points of the interval in ascending sequence.

    • resabs

      private final double[] resabs

    • resasc

      private final double[] resasc

    • result

      private final double[] result
      Approximation to the integral.

    • rlist

      private double[] rlist
      The integral approximations of the subintervals.

    • routine

      private int routine
      DOCUMENT ME!

    • rslst

      private double[] rslst
      In dqawfe vector of dimension at least limlst. rslst[k-1] contains the integral contribution over the interval (a+(k-1)c, a+kc), where c = (2*int(abs(omega))+1)*pi/abs(omega), k = 1, 2, ..., lst. Note that if omega = 0, rslst[0] contains the value of the integral over (a, infinity).

    • uflow

      private final double uflow
      DOCUMENT ME!

    • upper

      private double upper
      DOCUMENT ME!

    • selfTest

      private boolean selfTest

    • testCase

      private int testCase

    • actualAnswer

      private double actualAnswer

  • Constructor Details

    • Integration2

      public Integration2()
      Constructor for running self tests.

    • Integration2

      public Integration2(double lower, double upper, int routine, double epsabs, double epsrel)
      Constructor for Integration Used with routine = DQNG.
      Parameters:
      lower - DOCUMENT ME!
      upper - DOCUMENT ME!
      routine - DOCUMENT ME!
      epsabs - DOCUMENT ME!
      epsrel - DOCUMENT ME!

    • Integration2

      public Integration2(double bound, int routine, int inf, double epsabs, double epsrel, int limit)
      Constructor used with routine = DQAGIE.
      Parameters:
      bound - double
      routine - int
      inf - int
      epsabs - double
      epsrel - double
      limit - int

    • Integration2

      public Integration2(double lower, double upper, int routine, double epsabs, double epsrel, int limit)
      Constructor for Integration Used with routine = DQAGSE.
      Parameters:
      lower - DOCUMENT ME!
      upper - DOCUMENT ME!
      routine - DOCUMENT ME!
      epsabs - DOCUMENT ME!
      epsrel - DOCUMENT ME!
      limit - DOCUMENT ME!

    • Integration2

      public Integration2(double lower, double upper, double c, int routine, double epsabs, double epsrel, int limit)
      Constructor for Integration2 Used with routine = DQAWCE.
      Parameters:
      lower - DOCUMENT ME!
      upper - DOCUMENT ME!
      c -
      routine - DOCUMENT ME!
      epsabs - DOCUMENT ME!
      epsrel - DOCUMENT ME!
      limit - DOCUMENT ME!

    • Integration2

      public Integration2(double lower, int routine, double omega, int integr, double epsabs, int limlst, int limit, int maxp1)
      Constructor must be used with routine dqawfe
      Parameters:
      lower -
      routine -
      omega -
      integr -
      epsabs -
      limlst -
      limit -
      maxp1 -

    • Integration2

      public Integration2(double lower, double upper, int routine, double omega, int integr, double epsabs, double epsrel, int limit, int icall, int maxp1)
      Constructor must be used with routine dqawoe
      Parameters:
      lower -
      upper -
      routine -
      omega -
      integr -
      epsabs -
      epsrel -
      limit -
      icall -
      maxp1 -

    • Integration2

      public Integration2(double lower, double upper, int routine, double alfa, double beta, int integr, double epsabs, double epsrel, int limit)
      This constructor must be used with routine DQAWSE.
      Parameters:
      lower -
      upper -
      routine -
      alfa -
      beta -
      integr -
      epsabs -
      epsrel -
      limit -

    • Integration2

      public Integration2(double lower, double upper, int routine, double[] breakPoints, double epsabs, double epsrel, int limit)
      Constructor for Integration Used with routine = DQAGPE.
      Parameters:
      lower - DOCUMENT ME!
      upper - DOCUMENT ME!
      routine - DOCUMENT ME!
      breakPoints - DOCUMENT ME!
      epsabs - DOCUMENT ME!
      epsrel - DOCUMENT ME!
      limit - DOCUMENT ME!

    • Integration2

      public Integration2(double lower, double upper, int routine, int key, double epsabs, double epsrel, int limit)
      Constructor for Integration Used with routine = DQAGE.
      Parameters:
      lower - DOCUMENT ME!
      upper - DOCUMENT ME!
      routine - DOCUMENT ME!
      key - DOCUMENT ME!
      epsabs - DOCUMENT ME!
      epsrel - DOCUMENT ME!
      limit - DOCUMENT ME!

  • Method Details

    • intFunc

      public abstract double intFunc(double x)
      DOCUMENT ME!
      Parameters:
      x - DOCUMENT ME!
      Returns:
      DOCUMENT ME!

    • intFuncTest

      private double intFuncTest(double x)

    • dgtsl

      private void dgtsl(int n, double[] c, double[] d, double[] e, double[] b, int[] info)
      This is a port of the FORTRAN routine whose header is given below: c c dgtsl given a general tridiagonal matrix and a right hand c side will find the solution. c c on entry c c n integer c is the order of the tridiagonal matrix. c c c double precision(n) c is the subdiagonal of the tridiagonal matrix. c c(2) through c(n) should contain the subdiagonal. c on output c is destroyed. c c d double precision(n) c is the diagonal of the tridiagonal matrix. c on output d is destroyed. c c e double precision(n) c is the superdiagonal of the tridiagonal matrix. c e(1) through e(n-1) should contain the superdiagonal. c on output e is destroyed. c c b double precision(n) c is the right hand side vector. c c on return c c b is the solution vector. c c info integer c = 0 normal value. c = k if the k-th element of the diagonal becomes c exactly zero. the subroutine returns when c this is detected. c c linpack. this version dated 08/14/78 . c jack dongarra, argonne national laboratory.
      Parameters:
      n -
      c -
      d -
      e -
      b -

    • dqagie

      public void dqagie(double bound, int inf, double epsrel)
      This is a port of the orginal FORTRAN code whose header is given below: c***begin prologue dqagie c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a3a1,h2a4a1 c***keywords automatic integrator, infinite intervals, c general-purpose, transformation, extrapolation, c globally adaptive c***author piessens,robert,appl. math invalid input: '&' progr. div - k.u.leuven c de doncker,elise,appl. math invalid input: '&' progr. div - k.u.leuven c***purpose the routine calculates an approximation result to a given c integral i = integral of f over (bound,+infinity) c or i = integral of f over (-infinity,bound) c or i = integral of f over (-infinity,+infinity), c hopefully satisfying following claim for accuracy c abs(i-result).le.max(epsabs,epsrel*abs(i)) c***description c c integration over infinite intervals c standard fortran subroutine

    • dqagpe

      public void dqagpe()
      dqagpe. This is a port of the original FORTRAN code whose header is given below: c***begin prologue dqagie c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a3a1,h2a4a1 c***keywords automatic integrator, infinite intervals, c general-purpose, transformation, extrapolation, c globally adaptive c***author piessens,robert,appl. math invalid input: '&' progr. div - k.u.leuven c de doncker,elise,appl. math invalid input: '&' progr. div - k.u.leuven c***purpose the routine calculates an approximation result to a given c integral i = integral of f over (bound,+infinity) c or i = integral of f over (-infinity,bound) c or i = integral of f over (-infinity,+infinity), c hopefully satisfying following claim for accuracy c abs(i-result).le.max(epsabs,epsrel*abs(i)) c***description c c integration over infinite intervals c standard fortran subroutine

    • driver

      public void driver()
      DOCUMENT ME!

    • getAbserr

      public double getAbserr()
      DOCUMENT ME!
      Returns:
      abserr

    • getErrorStatus

      public int getErrorStatus()
      DOCUMENT ME!
      Returns:
      ier[0]

    • getIntegral

      public double getIntegral()
      DOCUMENT ME!
      Returns:
      result

    • getNeval

      public int getNeval()
      DOCUMENT ME!
      Returns:
      neval

    • dqage

      private void dqage()
      This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqage c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a1 c***keywords automatic integrator, general-purpose, c integrand examinator, globally adaptive, c gauss-kronrod c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose the routine calculates an approximation result to a given c definite integral i = integral of f over (a,b), c hopefully satisfying following claim for accuracy c abs(i-reslt).le.max(epsabs,epsrel*abs(i)). c***description c c computation of a definite integral c standard fortran subroutine c double precision version Calculates an integral over (lower, upper), hopefully satisfying the abs(actual integral - result) invalid input: '<'= max(epsabs, epsresl * abs(actual integral)).

    • dqagse

      private void dqagse()
      This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqagse c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a1 c***keywords automatic integrator, general-purpose, c (end point) singularities, extrapolation, c globally adaptive c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose the routine calculates an approximation result to a given c definite integral i = integral of f over (a,b), c hopefully satisfying following claim for accuracy c abs(i-result).le.max(epsabs,epsrel*abs(i)). c***description c c computation of a definite integral c standard fortran subroutine c double precision version This routine calculates the integral of intFunc over (lower, upper) and can handle end-point singularities at lower and upper. It tries to satisfy the claim for accuracy abs(actual integral - result) invalid input: '<'= max(epsabs, esprel * abs(actual integral))

    • dqawce

      private void dqawce()
      This is the port of the original FORTRAN routine whose header is given below: c***begin prologue dqawce c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a1,j4 c***keywords automatic integrator, special-purpose, c cauchy principal value, clenshaw-curtis method c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c*** purpose the routine calculates an approximation result to a c cauchy principal value i = integral of f*w over (a,b) c (w(x) = 1/(x-c), (c.ne.a, c.ne.b), hopefully satisfying c following claim for accuracy c abs(i-result).le.max(epsabs,epsrel*abs(i)) c***description c c computation of a cauchy principal value c standard fortran subroutine c double precision version

    • dqawfe

      private void dqawfe()
      This is the port of the original FORTRAN routine whose header is given below: c***begin prologue dqawfe c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a3a1 c***keywords automatic integrator, special-purpose, c fourier integrals, c integration between zeros with dqawoe, c convergence acceleration with dqelg c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c dedoncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose the routine calculates an approximation result to a c given fourier integal c i = integral of f(x)*w(x) over (a,infinity) c where w(x)=cos(omega*x) or w(x)=sin(omega*x), c hopefully satisfying following claim for accuracy c abs(i-result).le.epsabs. c***description c c computation of fourier integrals c standard fortran subroutine c double precision version c c parameters c on entry c f - double precision c function subprogram defining the integrand c function f(x). the actual name for f needs to c be declared e x t e r n a l in the driver program. c c a - double precision c lower limit of integration c c omega - double precision c parameter in the weight function c c integr - integer c indicates which weight function is used c integr = 1 w(x) = cos(omega*x) c integr = 2 w(x) = sin(omega*x) c if integr.ne.1.and.integr.ne.2, the routine will c end with ier = 6. c c epsabs - double precision c absolute accuracy requested, epsabs.gt.0 c if epsabs.le.0, the routine will end with ier = 6. c c limlst - integer c limlst gives an upper bound on the number of c cycles, limlst.ge.1. c if limlst.lt.3, the routine will end with ier = 6. c c limit - integer c gives an upper bound on the number of subintervals c allowed in the partition of each cycle, limit.ge.1 c each cycle, limit.ge.1. c c maxp1 - integer c gives an upper bound on the number of c chebyshev moments which can be stored, i.e. c for the intervals of lengths abs(b-a)*2**(-l), c l=0,1, ..., maxp1-2, maxp1.ge.1 c c on return c result - double precision c approximation to the integral x c c abserr - double precision c estimate of the modulus of the absolute error, c which should equal or exceed abs(i-result) c c neval - integer c number of integrand evaluations c c ier - ier = 0 normal and reliable termination of c the routine. it is assumed that the c requested accuracy has been achieved. c ier.gt.0 abnormal termination of the routine. the c estimates for integral and error are less c reliable. it is assumed that the requested c accuracy has not been achieved. c error messages c if omega.ne.0 c ier = 1 maximum number of cycles allowed c has been achieved., i.e. of subintervals c (a+(k-1)c,a+kc) where c c = (2*int(abs(omega))+1)*pi/abs(omega), c for k = 1, 2, ..., lst. c one can allow more cycles by increasing c the value of limlst (and taking the c according dimension adjustments into c account). c examine the array iwork which contains c the error flags on the cycles, in order to c look for eventual local integration c difficulties. if the position of a local c difficulty can be determined (e.g. c singularity, discontinuity within the c interval) one will probably gain from c splitting up the interval at this point c and calling appropriate integrators on c the subranges. c = 4 the extrapolation table constructed for c convergence acceleration of the series c formed by the integral contributions over c the cycles, does not converge to within c the requested accuracy. as in the case of c ier = 1, it is advised to examine the c array iwork which contains the error c flags on the cycles. c = 6 the input is invalid because c (integr.ne.1 and integr.ne.2) or c epsabs.le.0 or limlst.lt.3. c result, abserr, neval, lst are set c to zero. c = 7 bad integrand behaviour occurs within one c or more of the cycles. location and type c of the difficulty involved can be c determined from the vector ierlst. here c lst is the number of cycles actually c needed (see below). c ierlst(k) = 1 the maximum number of c subdivisions (= limit) has c been achieved on the k th c cycle. c = 2 occurrence of roundoff error c is detected and prevents the c tolerance imposed on the c k th cycle, from being c achieved. c = 3 extremely bad integrand c behaviour occurs at some c points of the k th cycle. c = 4 the integration procedure c over the k th cycle does c not converge (to within the c required accuracy) due to c roundoff in the c extrapolation procedure c invoked on this cycle. it c is assumed that the result c on this interval is the c best which can be obtained. c = 5 the integral over the k th c cycle is probably divergent c or slowly convergent. it c must be noted that c divergence can occur with c any other value of c ierlst(k). c if omega = 0 and integr = 1, c the integral is calculated by means of dqagie c and ier = ierlst(1) (with meaning as described c for ierlst(k), k = 1). c c rslst - double precision c vector of dimension at least limlst c rslst(k) contains the integral contribution c over the interval (a+(k-1)c,a+kc) where c c = (2*int(abs(omega))+1)*pi/abs(omega), c k = 1, 2, ..., lst. c note that, if omega = 0, rslst(1) contains c the value of the integral over (a,infinity). c c erlst - double precision c vector of dimension at least limlst c erlst(k) contains the error estimate corresponding c with rslst(k). c c ierlst - integer c vector of dimension at least limlst c ierlst(k) contains the error flag corresponding c with rslst(k). for the meaning of the local error c flags see description of output parameter ier. c c lst - integer c number of subintervals needed for the integration c if omega = 0 then lst is set to 1. c c alist, blist, rlist, elist - double precision c vector of dimension at least limit, c c iord, nnlog - integer c vector of dimension at least limit, providing c space for the quantities needed in the subdivision c process of each cycle c c chebmo - double precision c array of dimension at least (maxp1,25), providing c space for the chebyshev moments needed within the c cycles c c***references (none) c***routines called d1mach,dqagie,dqawoe,dqelg c***end prologue dqawfe c

    • dqawoe

      private void dqawoe(double a, double b, double epsabs, double epsrel, int icall, double[] result, double[] abserr, int[] neval, int[] ier)
      This is the port of the original FORTRAN routine whose header is given below: c***begin prologue dqawoe c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a1 c***keywords automatic integrator, special-purpose, c integrand with oscillatory cos or sin factor, c clenshaw-curtis method, (end point) singularities, c extrapolation, globally adaptive c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose the routine calculates an approximation result to a given c definite integral c i = integral of f(x)*w(x) over (a,b) c where w(x) = cos(omega*x) or w(x)=sin(omega*x), c hopefully satisfying following claim for accuracy c abs(i-result).le.max(epsabs,epsrel*abs(i)). c***description c c computation of oscillatory integrals c standard fortran subroutine c double precision version c c parameters c on entry c f - double precision c function subprogram defining the integrand c function f(x). the actual name for f needs to be c declared e x t e r n a l in the driver program. c c a - double precision c lower limit of integration c c b - double precision c upper limit of integration c c omega - double precision c parameter in the integrand weight function c c integr - integer c indicates which of the weight functions is to be c used c integr = 1 w(x) = cos(omega*x) c integr = 2 w(x) = sin(omega*x) c if integr.ne.1 and integr.ne.2, the routine c will end with ier = 6. c c epsabs - double precision c absolute accuracy requested c epsrel - double precision c relative accuracy requested c if epsabs.le.0 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), c the routine will end with ier = 6. c c limit - integer c gives an upper bound on the number of subdivisions c in the partition of (a,b), limit.ge.1. c c icall - integer c if dqawoe is to be used only once, icall must c be set to 1. assume that during this call, the c chebyshev moments (for clenshaw-curtis integration c of degree 24) have been computed for intervals of c lenghts (abs(b-a))*2**(-l), l=0,1,2,...momcom-1. c if icall.gt.1 this means that dqawoe has been c called twice or more on intervals of the same c length abs(b-a). the chebyshev moments already c computed are then re-used in subsequent calls. c if icall.lt.1, the routine will end with ier = 6. c c maxp1 - integer c gives an upper bound on the number of chebyshev c moments which can be stored, i.e. for the c intervals of lenghts abs(b-a)*2**(-l), c l=0,1, ..., maxp1-2, maxp1.ge.1. c if maxp1.lt.1, the routine will end with ier = 6. c c on return c result - double precision c approximation to the integral c c abserr - double precision c estimate of the modulus of the absolute error, c which should equal or exceed abs(i-result) c c neval - integer c number of integrand evaluations c c ier - integer c ier = 0 normal and reliable termination of the c routine. it is assumed that the c requested accuracy has been achieved. c - ier.gt.0 abnormal termination of the routine. c the estimates for integral and error are c less reliable. it is assumed that the c requested accuracy has not been achieved. c error messages c ier = 1 maximum number of subdivisions allowed c has been achieved. one can allow more c subdivisions by increasing the value of c limit (and taking according dimension c adjustments into account). however, if c this yields no improvement it is advised c to analyze the integrand, in order to c determine the integration difficulties. c if the position of a local difficulty can c be determined (e.g. singularity, c discontinuity within the interval) one c will probably gain from splitting up the c interval at this point and calling the c integrator on the subranges. if possible, c an appropriate special-purpose integrator c should be used which is designed for c handling the type of difficulty involved. c = 2 the occurrence of roundoff error is c detected, which prevents the requested c tolerance from being achieved. c the error may be under-estimated. c = 3 extremely bad integrand behaviour occurs c at some points of the integration c interval. c = 4 the algorithm does not converge. c roundoff error is detected in the c extrapolation table. c it is presumed that the requested c tolerance cannot be achieved due to c roundoff in the extrapolation table, c and that the returned result is the c best which can be obtained. c = 5 the integral is probably divergent, or c slowly convergent. it must be noted that c divergence can occur with any other value c of ier.gt.0. c = 6 the input is invalid, because c (epsabs.le.0 and c epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) c or (integr.ne.1 and integr.ne.2) or c icall.lt.1 or maxp1.lt.1. c result, abserr, neval, last, rlist(1), c elist(1), iord(1) and nnlog(1) are set c to zero. alist(1) and blist(1) are set c to a and b respectively. c c last - integer c on return, last equals the number of c subintervals produces in the subdivision c process, which determines the number of c significant elements actually in the c work arrays. c alist - double precision c vector of dimension at least limit, the first c last elements of which are the left c end points of the subintervals in the partition c of the given integration range (a,b) c c blist - double precision c vector of dimension at least limit, the first c last elements of which are the right c end points of the subintervals in the partition c of the given integration range (a,b) c c rlist - double precision c vector of dimension at least limit, the first c last elements of which are the integral c approximations on the subintervals c c elist - double precision c vector of dimension at least limit, the first c last elements of which are the moduli of the c absolute error estimates on the subintervals c c iord - integer c vector of dimension at least limit, the first k c elements of which are pointers to the error c estimates over the subintervals, c such that elist(iord(1)), ..., c elist(iord(k)) form a decreasing sequence, with c k = last if last.le.(limit/2+2), and c k = limit+1-last otherwise. c c nnlog - integer c vector of dimension at least limit, containing the c subdivision levels of the subintervals, i.e. c iwork(i) = l means that the subinterval c numbered i is of length abs(b-a)*2**(1-l) c c on entry and return c momcom - integer c indicating that the chebyshev moments c have been computed for intervals of lengths c (abs(b-a))*2**(-l), l=0,1,2, ..., momcom-1, c momcom.lt.maxp1 c c chebmo - double precision c array of dimension (maxp1,25) containing the c chebyshev moments c c***references (none) c***routines called d1mach,dqc25f,dqelg,dqpsrt c***end prologue dqawoe c

    • dqawse

      private void dqawse()
      This is the port of the original FORTRAN routine whose header is given below: c***begin prologue dqawse c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a1 c***keywords automatic integrator, special-purpose, c algebraico-logarithmic end point singularities, c clenshaw-curtis method c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose the routine calculates an approximation result to a given c definite integral i = integral of f*w over (a,b), c (where w shows a singular behaviour at the end points, c see parameter integr). c hopefully satisfying following claim for accuracy c abs(i-result).le.max(epsabs,epsrel*abs(i)). c***description c c integration of functions having algebraico-logarithmic c end point singularities c standard fortran subroutine c double precision version c c parameters c on entry c f - double precision c function subprogram defining the integrand c function f(x). the actual name for f needs to be c declared e x t e r n a l in the driver program. c c a - double precision c lower limit of integration c c b - double precision c upper limit of integration, b.gt.a c if b.le.a, the routine will end with ier = 6. c c alfa - double precision c parameter in the weight function, alfa.gt.(-1) c if alfa.le.(-1), the routine will end with c ier = 6. c c beta - double precision c parameter in the weight function, beta.gt.(-1) c if beta.le.(-1), the routine will end with c ier = 6. c c integr - integer c indicates which weight function is to be used c = 1 (x-a)**alfa*(b-x)**beta c = 2 (x-a)**alfa*(b-x)**beta*log(x-a) c = 3 (x-a)**alfa*(b-x)**beta*log(b-x) c = 4 (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x) c if integr.lt.1 or integr.gt.4, the routine c will end with ier = 6. c c epsabs - double precision c absolute accuracy requested c epsrel - double precision c relative accuracy requested c if epsabs.le.0 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), c the routine will end with ier = 6. c c limit - integer c gives an upper bound on the number of subintervals c in the partition of (a,b), limit.ge.2 c if limit.lt.2, the routine will end with ier = 6. c c on return c result - double precision c approximation to the integral c c abserr - double precision c estimate of the modulus of the absolute error, c which should equal or exceed abs(i-result) c c neval - integer c number of integrand evaluations c c ier - integer c ier = 0 normal and reliable termination of the c routine. it is assumed that the requested c accuracy has been achieved. c ier.gt.0 abnormal termination of the routine c the estimates for the integral and error c are less reliable. it is assumed that the c requested accuracy has not been achieved. c error messages c = 1 maximum number of subdivisions allowed c has been achieved. one can allow more c subdivisions by increasing the value of c limit. however, if this yields no c improvement, it is advised to analyze the c integrand in order to determine the c integration difficulties which prevent the c requested tolerance from being achieved. c in case of a jump discontinuity or a local c singularity of algebraico-logarithmic type c at one or more interior points of the c integration range, one should proceed by c splitting up the interval at these c points and calling the integrator on the c subranges. c = 2 the occurrence of roundoff error is c detected, which prevents the requested c tolerance from being achieved. c = 3 extremely bad integrand behaviour occurs c at some points of the integration c interval. c = 6 the input is invalid, because c b.le.a or alfa.le.(-1) or beta.le.(-1), or c integr.lt.1 or integr.gt.4, or c (epsabs.le.0 and c epsrel.lt.max(50*rel.mach.acc.,0.5d-28), c or limit.lt.2. c result, abserr, neval, rlist(1), elist(1), c iord(1) and last are set to zero. alist(1) c and blist(1) are set to a and b c respectively. c c alist - double precision c vector of dimension at least limit, the first c last elements of which are the left c end points of the subintervals in the partition c of the given integration range (a,b) c c blist - double precision c vector of dimension at least limit, the first c last elements of which are the right c end points of the subintervals in the partition c of the given integration range (a,b) c c rlist - double precision c vector of dimension at least limit,the first c last elements of which are the integral c approximations on the subintervals c c elist - double precision c vector of dimension at least limit, the first c last elements of which are the moduli of the c absolute error estimates on the subintervals c c iord - integer c vector of dimension at least limit, the first k c of which are pointers to the error c estimates over the subintervals, so that c elist(iord(1)), ..., elist(iord(k)) with k = last c if last.le.(limit/2+2), and k = limit+1-last c otherwise form a decreasing sequence c c last - integer c number of subintervals actually produced in c the subdivision process c c***references (none) c***routines called d1mach,dqc25s,dqmomo,dqpsrt c***end prologue dqawse c

    • dqc25c

      private void dqc25c(double a, double b, double c, double[] result, double[] abserr, int[] krul, int[] neval)
      This is the port of the original FORTRAN routine whose header is given below: c***begin prologue dqc25c c***date written 810101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a2,j4 c***keywords 25-point clenshaw-curtis integration c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose to compute i = integral of f*w over (a,b) with c error estimate, where w(x) = 1/(x-c) c***description c c integration rules for the computation of cauchy c principal value integrals c standard fortran subroutine c double precision version c c c a - double precision c left end point of the integration interval c c b - double precision c right end point of the integration interval, b.gt.a c c c - double precision c parameter in the weight function c c result - double precision c approximation to the integral c result is computed by using a generalized c clenshaw-curtis method if c lies within ten percent c of the integration interval. in the other case the c 15-point kronrod rule obtained by optimal addition c of abscissae to the 7-point gauss rule, is applied. c c abserr - double precision c estimate of the modulus of the absolute error, c which should equal or exceed abs(i-result) c c krul - integer c key which is decreased by 1 if the 15-point c gauss-kronrod scheme has been used c c neval - integer c number of integrand evaluations c c....................................................................... c***references (none) c***routines called dqcheb,dqk15w,dqwgtc c***end prologue dqc25c c
      Parameters:
      a -
      b -
      c -
      result -
      abserr -
      krul -
      neval -

    • dqc25f

      private void dqc25f(double a, double b, double omega, int nrmom, int ksave, double[] result, double[] abserr, int[] neval, double[] resabs, double[] resasc)
      This is a port of the original FORTRAN whose header is given below: c***begin prologue dqc25f c***date written 810101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a2 c***keywords integration rules for functions with cos or sin c factor, clenshaw-curtis, gauss-kronrod c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose to compute the integral i=integral of f(x) over (a,b) c where w(x) = cos(omega*x) or w(x)=sin(omega*x) and to c compute j = integral of abs(f) over (a,b). for small value c of omega or small intervals (a,b) the 15-point gauss-kronro c rule is used. otherwise a generalized clenshaw-curtis c method is used. c***description c c integration rules for functions with cos or sin factor c standard fortran subroutine c double precision version c c parameters c on entry c f - double precision c function subprogram defining the integrand c function f(x). the actual name for f needs to c be declared e x t e r n a l in the calling program. c c a - double precision c lower limit of integration c c b - double precision c upper limit of integration c c omega - double precision c parameter in the weight function c c integr - integer c indicates which weight function is to be used c integr = 1 w(x) = cos(omega*x) c integr = 2 w(x) = sin(omega*x) c c nrmom - integer c the length of interval (a,b) is equal to the length c of the original integration interval divided by c 2**nrmom (we suppose that the routine is used in an c adaptive integration process, otherwise set c nrmom = 0). nrmom must be zero at the first call. c c maxp1 - integer c gives an upper bound on the number of chebyshev c moments which can be stored, i.e. for the c intervals of lengths abs(bb-aa)*2**(-l), c l = 0,1,2, ..., maxp1-2. c c ksave - integer c key which is one when the moments for the c current interval have been computed c c on return c result - double precision c approximation to the integral i c c abserr - double precision c estimate of the modulus of the absolute c error, which should equal or exceed abs(i-result) c c neval - integer c number of integrand evaluations c c resabs - double precision c approximation to the integral j c c resasc - double precision c approximation to the integral of abs(f-i/(b-a)) c c on entry and return c momcom - integer c for each interval length we need to compute the c chebyshev moments. momcom counts the number of c intervals for which these moments have already been c computed. if nrmom.lt.momcom or ksave = 1, the c chebyshev moments for the interval (a,b) have c already been computed and stored, otherwise we c compute them and we increase momcom. c c chebmo - double precision c array of dimension at least (maxp1,25) containing c the modified chebyshev moments for the first momcom c momcom interval lengths c c ...................................................................... c***references (none) c***routines called d1mach,dgtsl,dqcheb,dqk15w,dqwgtf c***end prologue dqc25f c
      Parameters:
      a -
      b -
      omega -
      nrmom -
      ksave -
      result -
      abserr -
      neval -
      resabs -
      resasc -

    • dqc25s

      private void dqc25s(double bl, double br, double[] ri, double[] rj, double[] rg, double[] rh, double[] result, double[] abserr, double[] resasc, int[] nev)
      This is a port of the original FORTRAN whose header is given below: c***begin prologue dqc25s c***date written 810101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a2 c***keywords 25-point clenshaw-curtis integration c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose to compute i = integral of f*w over (bl,br), with error c estimate, where the weight function w has a singular c behaviour of algebraico-logarithmic type at the points c a and/or b. (bl,br) is a part of (a,b). c***description c c integration rules for integrands having algebraico-logarithmic c end point singularities c standard fortran subroutine c double precision version c c parameters c f - double precision c function subprogram defining the integrand c f(x). the actual name for f needs to be declared c e x t e r n a l in the driver program. c c a - double precision c left end point of the original interval c c b - double precision c right end point of the original interval, b.gt.a c c bl - double precision c lower limit of integration, bl.ge.a c c br - double precision c upper limit of integration, br.le.b c c alfa - double precision c parameter in the weight function c c beta - double precision c parameter in the weight function c c ri,rj,rg,rh - double precision c modified chebyshev moments for the application c of the generalized clenshaw-curtis c method (computed in subroutine dqmomo) c c result - double precision c approximation to the integral c result is computed by using a generalized c clenshaw-curtis method if b1 = a or br = b. c in all other cases the 15-point kronrod c rule is applied, obtained by optimal addition of c abscissae to the 7-point gauss rule. c c abserr - double precision c estimate of the modulus of the absolute error, c which should equal or exceed abs(i-result) c c resasc - double precision c approximation to the integral of abs(f*w-i/(b-a)) c c integr - integer c which determines the weight function c = 1 w(x) = (x-a)**alfa*(b-x)**beta c = 2 w(x) = (x-a)**alfa*(b-x)**beta*log(x-a) c = 3 w(x) = (x-a)**alfa*(b-x)**beta*log(b-x) c = 4 w(x) = (x-a)**alfa*(b-x)**beta*log(x-a)* c log(b-x) c c nev - integer c number of integrand evaluations c***references (none) c***routines called dqcheb,dqk15w c***end prologue dqc25s c
      Parameters:
      bl -
      br -
      ri -
      rj -
      rg -
      rh -
      result -
      abserr -
      resasc -
      nev -

    • dqcheb

      private void dqcheb(double[] x, double[] fval, double[] cheb12, double[] cheb24)
      This is a port of the original FORTRAN whose header is given below: c***begin prologue dqcheb c***refer to dqc25c,dqc25f,dqc25s c***routines called (none) c***revision date 830518 (yymmdd) c***keywords chebyshev series expansion, fast fourier transform c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose this routine computes the chebyshev series expansion c of degrees 12 and 24 of a function using a c fast fourier transform method c f(x) = sum(k=1,..,13) (cheb12(k)*t(k-1,x)), c f(x) = sum(k=1,..,25) (cheb24(k)*t(k-1,x)), c where t(k,x) is the chebyshev polynomial of degree k. c***description c c chebyshev series expansion c standard fortran subroutine c double precision version c c parameters c on entry c x - double precision c vector of dimension 11 containing the c values cos(k*pi/24), k = 1, ..., 11 c c fval - double precision c vector of dimension 25 containing the c function values at the points c (b+a+(b-a)*cos(k*pi/24))/2, k = 0, ...,24, c where (a,b) is the approximation interval. c fval(1) and fval(25) are divided by two c (these values are destroyed at output). c c on return c cheb12 - double precision c vector of dimension 13 containing the c chebyshev coefficients for degree 12 c c cheb24 - double precision c vector of dimension 25 containing the c chebyshev coefficients for degree 24 c c***end prologue dqcheb c
      Parameters:
      x -
      fval -
      cheb12 -
      cheb24 -

    • dqelg

      private void dqelg(int[] n, double[] epstab, double[] result, double[] abserr, double[] res31a, int[] nres)
      This is a port of the original FORTRAN whose header is given below: c***begin prologue dqelg c***refer to dqagie,dqagoe,dqagpe,dqagse c***routines called d1mach c***revision date 830518 (yymmdd) c***keywords epsilon algorithm, convergence acceleration, c extrapolation c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math invalid input: '&' progr. div. - k.u.leuven c***purpose the routine determines the limit of a given sequence of c approximations, by means of the epsilon algorithm of c p.wynn. an estimate of the absolute error is also given. c the condensed epsilon table is computed. only those c elements needed for the computation of the next diagonal c are preserved. c***description c c epsilon algorithm c standard fortran subroutine c double precision version This routine determines the limit of a given sequence of approximations, by means of the epsilon algorithm of P. Wynn. An estimate of the absolute error is also given. The condensed epsilon table is computed. Only those elements needed for the computation of the next diagonal are preserved.
      Parameters:
      n - epstab[n[0]-1] contains the new element in the first column of the epsilon table
      epstab - 52 element array containing the elements of the two lower diagonals of the triangular epsilon table. The elements are numbered starting at the right-hand corner of the triangle
      result - resulting approximation to the integral The element in the new diagonal with the least value of error.
      abserr - estimate of the absolute error computed from result and the 3 previous results
      res31a - array containing the last 3 results
      nres - number of calls to the routine (should be zero at first call)

    • dqk15

      private void dqk15(double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
      The is a port of the original FORTRAN whose header is given below: c***begin prologue dqk15 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 15-point gauss-kronrod rules c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div - k.u.leuven c***purpose to compute i = integral of f over (a,b), with error c estimate c j = integral of abs(f) over (a,b) c***description c c integration rules c standard fortran subroutine c double precision version This routine computes theintegral of the intFunc over (a, b) and the integral of the abs(intFunc) over (a, b).
      Parameters:
      a - lower limit for integration
      b - upper limit for integration
      result - Approximation to the integral. The result is computed by applying the 15-point kronrod rule (resk) obtained by optimal addition of abscissae to the 7-point gauss rule (resg).
      abserr - Estimate of the modulus of the absolute error, which should equal or exceed abs(integral of intFunc over (a,b) - result)
      resabs - Approximation to integral of abs(intFunc)
      resasc - Approximation to the integral of abs(intFunc - integral of intFunc/(b - a)) over (a, b)

    • dqk15i

      private void dqk15i(double boun, int inf, double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
      This is a port of the original FORTRAN whose header is given below: c***begin prologue dqk15i c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a3a2,h2a4a2 c***keywords 15-point transformed gauss-kronrod rules c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose the original (infinite integration range is mapped c onto the interval (0,1) and (a,b) is a part of (0,1). c it is the purpose to compute c i = integral of transformed integrand over (a,b), c j = integral of abs(transformed integrand) over (a,b). c***description c c integration rule c standard fortran subroutine c double precision version The original infinite integration range is mapped onto the interval (0, 1) and (a, b) is a part of (0, 1). This routine computes the integral of the transformed integrand over (a, b) and the integral of the abs(transformed integrand) over (a, b).
      Parameters:
      boun - finite bound of original integration range (set to zero if inf = +2)
      inf - If inf = -1, the original interval is (-infinity, bound) If inf = +1, the original interval is (bound, +infinity) If inf = +2, the original interval is (-infinity, +infinity) and the integral is computed as the sum of two integrals, one over (-infinity, 0) and one over (0, +infinity)
      a - lower limit for integration over subrange of (0, 1)
      b - upper limit for integration over subrange of (0, 1)
      result - Approximation to the integral. The result is computed by applying the 15-point kronrod rule (resk) obtained by optimal addition of abscissae to the 7-point gauss rule (resg).
      abserr - Estimate of the modulus of the absolute error, which should equal or exceed abs(integral of transformed integrand - result)
      resabs - Approximation to integral of abs(transformed integrand)
      resasc - Approximation to the integral of abs(transformed integrand - integral of transformed integrand/(b - a)) over (a, b)

    • dqk15w

      private void dqk15w(double p1, double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
      This a a port of the original FORTRAN whose header is given below: c***begin prologue dqk15w c***date written 810101 (yymmdd) c***revision date 830518 (mmddyy) c***category no. h2a2a2 c***keywords 15-point gauss-kronrod rules c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose to compute i = integral of f*w over (a,b), with error c estimate c j = integral of abs(f*w) over (a,b) c***description c c integration rules c standard fortran subroutine c double precision version c c parameters c on entry c c p1, p2, p3, p4 - double precision c parameters in the weight function c c kp - integer c key for indicating the type of weight function c c a - double precision c lower limit of integration c c b - double precision c upper limit of integration c c on return c result - double precision c approximation to the integral i c result is computed by applying the 15-point c kronrod rule (resk) obtained by optimal addition c of abscissae to the 7-point gauss rule (resg). c c abserr - double precision c estimate of the modulus of the absolute error, c which should equal or exceed abs(i-result) c c resabs - double precision c approximation to the integral of abs(f) c c resasc - double precision c approximation to the integral of abs(f-i/(b-a)) c
      Parameters:
      p1 -
      a -
      b -
      result -
      abserr -
      resabs -
      resasc -
      kp -

    • dqk21

      private void dqk21(double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
      This is a port of the original FORTRAN whose header is given below: c***begin prologue dqk21 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 21-point gauss-kronrod rules c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose to compute i = integral of f over (a,b), with error c estimate c j = integral of abs(f) over (a,b) c***description c c integration rules c standard fortran subroutine c double precision version Computes integral of function over (a,b), with error estimate and computes integral of abs(function) over (a,b).
      Parameters:
      a - lower limit of integration
      b - upper limit of integration
      result - Approximation to the integral. Result is computed by applying the 21-point kronrod rule (resk) obtained by optimal addition of abscissae to the 10-point gauss rule (resg).
      abserr - Estimate of the modulus of the absolute error, which should not exceed abs(actual integral - result)
      resabs - Approximation to the integral of abs(function) over (a,b).
      resasc - Approximation to the integral of abs(function - actual integral/(b-a)) over (a,b)

    • dqk31

      private void dqk31(double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
      This is a port of the original FORTRAN whose header is given below: c***begin prologue dqk31 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 31-point gauss-kronrod rules c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose to compute i = integral of f over (a,b) with error c estimate c j = integral of abs(f) over (a,b) c***description c c integration rules c standard fortran subroutine c double precision version This routine computes theintegral of the intFunc over (a, b) and the integral of the abs(intFunc) over (a, b).
      Parameters:
      a - lower limit for integration
      b - upper limit for integration
      result - Approximation to the integral. The result is computed by applying the 31-point kronrod rule (resk) obtained by optimal addition of abscissae to the 15-point gauss rule (resg).
      abserr - Estimate of the modulus of the absolute error, which should equal or exceed abs(integral of intFunc over (a,b) - result)
      resabs - Approximation to integral of abs(intFunc)
      resasc - Approximation to the integral of abs(intFunc - integral of intFunc/(b - a)) over (a, b)

    • dqk41

      private void dqk41(double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
      This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqk41 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 41-point gauss-kronrod rules c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose to compute i = integral of f over (a,b), with error c estimate c j = integral of abs(f) over (a,b) c***description c c integration rules c standard fortran subroutine c double precision version This routine computes theintegral of the intFunc over (a, b) and the integral of the abs(intFunc) over (a, b).
      Parameters:
      a - lower limit for integration
      b - upper limit for integration
      result - Approximation to the integral. The result is computed by applying the 41-point kronrod rule (resk) obtained by optimal addition of abscissae to the 20-point gauss rule (resg).
      abserr - Estimate of the modulus of the absolute error, which should equal or exceed abs(integral of intFunc over (a,b) - result)
      resabs - Approximation to integral of abs(intFunc)
      resasc - Approximation to the integral of abs(intFunc - integral of intFunc/(b - a)) over (a, b)

    • dqk51

      private void dqk51(double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
      This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqk51 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 51-point gauss-kronrod rules c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math invalid input: '&' progr. div. - k.u.leuven c***purpose to compute i = integral of f over (a,b) with error c estimate c j = integral of abs(f) over (a,b) c***description c c integration rules c standard fortran subroutine c double precision version This routine computes theintegral of the intFunc over (a, b) and the integral of the abs(intFunc) over (a, b).
      Parameters:
      a - lower limit for integration
      b - upper limit for integration
      result - Approximation to the integral. The result is computed by applying the 51-point kronrod rule (resk) obtained by optimal addition of abscissae to the 25-point gauss rule (resg).
      abserr - Estimate of the modulus of the absolute error, which should equal or exceed abs(integral of intFunc over (a,b) - result)
      resabs - Approximation to integral of abs(intFunc)
      resasc - Approximation to the integral of abs(intFunc - integral of intFunc/(b - a)) over (a, b)

    • dqk61

      private void dqk61(double a, double b, double[] result, double[] abserr, double[] resabs, double[] resasc)
      This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqk61 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 61-point gauss-kronrod rules c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose to compute i = integral of f over (a,b) with error c estimate c j = integral of dabs(f) over (a,b) c***description c c integration rule c standard fortran subroutine c double precision version This routine computes the integral of the intFunc over (a, b) and the integral of the abs(intFunc) over (a, b).
      Parameters:
      a - lower limit for integration
      b - upper limit for integration
      result - Approximation to the integral. The result is computed by applying the 61-point kronrod rule (resk) obtained by optimal addition of abscissae to the 30-point gauss rule (resg).
      abserr - Estimate of the modulus of the absolute error, which should equal or exceed abs(integral of intFunc over (a,b) - result)
      resabs - Approximation to integral of abs(intFunc)
      resasc - Approximation to the integral of abs(intFunc - integral of intFunc/(b - a)) over (a, b)

    • dqmomo

      private void dqmomo(double[] ri, double[] rj, double[] rg, double[] rh)
      This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqmomo c***date written 820101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a1,c3a2 c***keywords modified chebyshev moments c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose this routine computes modified chebsyshev moments. the k-th c modified chebyshev moment is defined as the integral over c (-1,1) of w(x)*t(k,x), where t(k,x) is the chebyshev c polynomial of degree k. c***description c c modified chebyshev moments c standard fortran subroutine c double precision version c c parameters c alfa - double precision c parameter in the weight function w(x), alfa.gt.(-1) c c beta - double precision c parameter in the weight function w(x), beta.gt.(-1) c c ri - double precision c vector of dimension 25 c ri(k) is the integral over (-1,1) of c (1+x)**alfa*t(k-1,x), k = 1, ..., 25. c c rj - double precision c vector of dimension 25 c rj(k) is the integral over (-1,1) of c (1-x)**beta*t(k-1,x), k = 1, ..., 25. c c rg - double precision c vector of dimension 25 c rg(k) is the integral over (-1,1) of c (1+x)**alfa*log((1+x)/2)*t(k-1,x), k = 1, ..., 25. c c rh - double precision c vector of dimension 25 c rh(k) is the integral over (-1,1) of c (1-x)**beta*log((1-x)/2)*t(k-1,x), k = 1, ..., 25. c c integr - integer c input parameter indicating the modified c moments to be computed c integr = 1 compute ri, rj c = 2 compute ri, rj, rg c = 3 compute ri, rj, rh c = 4 compute ri, rj, rg, rh c c***references (none) c***routines called (none) c***end prologue dqmomo c
      Parameters:
      ri -
      rj -
      rg -
      rh -

    • dqng

      private void dqng()
      This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqng c***date written 800101 (yymmdd) c***revision date 810101 (yymmdd) c***category no. h2a1a1 c***keywords automatic integrator, smooth integrand, c non-adaptive, gauss-kronrod(patterson) c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl math invalid input: '&' progr. div. - k.u.leuven c kahaner,david,nbs - modified (2/82) c***purpose the routine calculates an approximation result to a c given definite integral i = integral of f over (a,b), c hopefully satisfying following claim for accuracy c abs(i-result).le.max(epsabs,epsrel*abs(i)). c***description c c non-adaptive integration c standard fortran subroutine c double precision version Calculates an integral over (lower, upper), hopefully satisfying the abs(actual integral - result) invalid input: '<'= max(epsabs, epsresl * abs(actual integral)) Result is obtained by applying the 21-point gauss-kronrod rule (res21) obtained by optimal addition of abscissae to the 10-point gauss rule (res10), or by applying the 43-point rule (res43) obtained by optimal addition of abscissae to the 21-point gauss-kronrod rule, or by applying the 87-point rule (res87) obtained by optimal addition of abscissae to the 43-point rule.

    • dqpsrt

      private void dqpsrt(int limit, int last, int[] maxErr, double[] ermax, double[] elist, int[] iord, int[] nrmax)
      This is a port of the original FORTRAN routine whose header is given below: c***begin prologue dqpsrt c***refer to dqage,dqagie,dqagpe,dqawse c***routines called (none) c***revision date 810101 (yymmdd) c***keywords sequential sorting c***author piessens,robert,appl. math. invalid input: '&' progr. div. - k.u.leuven c de doncker,elise,appl. math. invalid input: '&' progr. div. - k.u.leuven c***purpose this routine maintains the descending ordering in the c list of the local error estimated resulting from the c interval subdivision process. at each call two error c estimates are inserted using the sequential search c method, top-down for the largest error estimate and c bottom-up for the smallest error estimate. c***description c c ordering routine c standard fortran subroutine c double precision version This routine maintains the descending ordering in the list of the local error estimated resulting from the interval subdivision process. At each call two error estimates are inserted using the sequential search method, top-down for the largest error estimate and bottom-up for the smallest error estimate.
      Parameters:
      limit - maximum number of error estimates the list can contain
      last - number of error estimates currently in the list
      maxErr - maxErr points to the nrmax-th largest estimate currently in the list
      ermax - nrmax-th largest error estimate ermax[0] = elist[maxErr[0]]
      elist - the error estimates
      iord - Array of size last , the first k elements of which contain pointers to the error estimates, such that elist[iord[0]],...,elist[iord[k-1]] form a decreasing sequence, with k = last if last invalid input: '<'= (limit/2 + 2), and k = limit + 1 - last otherwise
      nrmax - maxErr[0] = iord[nrmax[0] - 1]