kronrodintegral.cpp

有很多的函数库

                                                     Real relativeAccuracy)
    : Integrator(absoluteAccuracy, maxEvaluations),
      relativeAccuracy_(relativeAccuracy) {}

    Real
    GaussKronrodNonAdaptive::integrate(const boost::function<Real (Real)>& f,
                                       Real a,
                                       Real b) const {
        Real result;
        //Size neval;
        Real fv1[5], fv2[5], fv3[5], fv4[5];
        Real savfun[21];  /* array of function values which have been computed */
        Real res10, res21, res43, res87;    /* 10, 21, 43 and 87 point results */
        Real err;
        Real resAbs; /* approximation to the integral of abs(f) */
        Real resasc; /* approximation to the integral of abs(f-i/(b-a)) */
        int k ;

        QL_REQUIRE(a<b, "b must be greater than a)");

        const Real halfLength = 0.5 * (b - a);
        const Real center = 0.5 * (b + a);
        const Real fCenter = f(center);

        // Compute the integral using the 10- and 21-point formula.

        res10 = 0;
        res21 = w21b[5] * fCenter;
        resAbs = w21b[5] * std::fabs(fCenter);

        for (k = 0; k < 5; k++) {
            Real abscissa = halfLength * x1[k];
            Real fval1 = f(center + abscissa);
            Real fval2 = f(center - abscissa);
            Real fval = fval1 + fval2;
            res10 += w10[k] * fval;
            res21 += w21a[k] * fval;
            resAbs += w21a[k] * (std::fabs(fval1) + std::fabs(fval2));
            savfun[k] = fval;
            fv1[k] = fval1;
            fv2[k] = fval2;
        }

        for (k = 0; k < 5; k++) {
            Real abscissa = halfLength * x2[k];
            Real fval1 = f(center + abscissa);
            Real fval2 = f(center - abscissa);
            Real fval = fval1 + fval2;
            res21 += w21b[k] * fval;
            resAbs += w21b[k] * (std::fabs(fval1) + std::fabs(fval2));
            savfun[k + 5] = fval;
            fv3[k] = fval1;
            fv4[k] = fval2;
        }

        result = res21 * halfLength;
        resAbs *= halfLength ;
        Real mean = 0.5 * res21;
        resasc = w21b[5] * std::fabs(fCenter - mean);

        for (k = 0; k < 5; k++)
            resasc += (w21a[k] * (std::fabs(fv1[k] - mean)
                        + std::fabs(fv2[k] - mean))
                        + w21b[k] * (std::fabs(fv3[k] - mean)
                        + std::fabs(fv4[k] - mean)));

        err = rescaleError ((res21 - res10) * halfLength, resAbs, resasc) ;
        resasc *= halfLength ;

        // test for convergence.
        if (err < absoluteAccuracy() || err < relativeAccuracy() * std::fabs(result)){
            setAbsoluteError(err);
            setNumberOfEvaluations(21);
            return result;
        }

        /* compute the integral using the 43-point formula. */

        res43 = w43b[11] * fCenter;

        for (k = 0; k < 10; k++)
            res43 += savfun[k] * w43a[k];

        for (k = 0; k < 11; k++){
            Real abscissa = halfLength * x3[k];
            Real fval = (f(center + abscissa)
                + f(center - abscissa));
            res43 += fval * w43b[k];
            savfun[k + 10] = fval;
            }

        // test for convergence.

        result = res43 * halfLength;
        err = rescaleError ((res43 - res21) * halfLength, resAbs, resasc);

       if (err < absoluteAccuracy() || err < relativeAccuracy() * std::fabs(result)){
            setAbsoluteError(err);
            setNumberOfEvaluations(43);
            return result;
        }

        /* compute the integral using the 87-point formula. */

        res87 = w87b[22] * fCenter;

        for (k = 0; k < 21; k++)
            res87 += savfun[k] * w87a[k];

        for (k = 0; k < 22; k++){
            Real abscissa = halfLength * x4[k];
            res87 += w87b[k] * (f(center + abscissa)
                + f(center - abscissa));
        }

        // test for convergence.
        result = res87 * halfLength ;
        err = rescaleError ((res87 - res43) * halfLength, resAbs, resasc);

        setAbsoluteError(err);
        setNumberOfEvaluations(87);
        return result;
    }

    Real
    GaussKronrodAdaptive::integrate(const boost::function<Real (Real)>& f,
                                    Real a,
                                    Real b) const {
        return integrateRecursively(f, a, b, absoluteAccuracy());
    }

    // weights for 7-point Gauss-Legendre integration
    // (only 4 values out of 7 are given as they are symmetric)
    static const Real g7w[] = { 0.417959183673469,
                                0.381830050505119,
                                0.279705391489277,
                                0.129484966168870 };
    // weights for 15-point Gauss-Kronrod integration
    static const Real k15w[] = { 0.209482141084728,
                                 0.204432940075298,
                                 0.190350578064785,
                                 0.169004726639267,
                                 0.140653259715525,
                                 0.104790010322250,
                                 0.063092092629979,
                                 0.022935322010529 };
    // abscissae (evaluation points)
    // for 15-point Gauss-Kronrod integration
    static const Real k15t[] = { 0.000000000000000,
                                 0.207784955007898,
                                 0.405845151377397,
                                 0.586087235467691,
                                 0.741531185599394,
                                 0.864864423359769,
                                 0.949107912342758,
                                 0.991455371120813 };

    Real GaussKronrodAdaptive::integrateRecursively(
                                    const boost::function<Real (Real)>& f,
                                    Real a,
                                    Real b,
                                    Real tolerance) const {

            Real halflength = (b - a) / 2;
            Real center = (a + b) / 2;

            Real g7; // will be result of G7 integral
            Real k15; // will be result of K15 integral

            Real t, fsum; // t (abscissa) and f(t)
            Real fc = f(center);
            g7 = fc * g7w[0];
            k15 = fc * k15w[0];

            // calculate g7 and half of k15
            Integer j, j2;
            for (j = 1, j2 = 2; j < 4; j++, j2 += 2) {
                t = halflength * k15t[j2];
                fsum = f(center - t) + f(center + t);
                g7  += fsum * g7w[j];
                k15 += fsum * k15w[j2];
            }

            // calculate other half of k15
            for (j2 = 1; j2 < 8; j2 += 2) {
                t = halflength * k15t[j2];
                fsum = f(center - t) + f(center + t);
                k15 += fsum * k15w[j2];
            }

            // multiply by (a - b) / 2
            g7 = halflength * g7;
            k15 = halflength * k15;

            // 15 more function evaluations have been used
            increaseNumberOfEvaluations(15);

            // error is <= k15 - g7
            // if error is larger than tolerance then split the interval
            // in two and integrate recursively
            if (std::fabs(k15 - g7) < tolerance) {
                return k15;
            } else {
                QL_REQUIRE(numberOfEvaluations()+30 <=
                           maxEvaluations(),
                           "maximum number of function evaluations "
                           "exceeded");
                return integrateRecursively(f, a, center, tolerance/2)
                    + integrateRecursively(f, center, b, tolerance/2);
            }
        }


    GaussKronrodAdaptive::GaussKronrodAdaptive(Real absoluteAccuracy,
                                               Size maxEvaluations)
    : Integrator(absoluteAccuracy, maxEvaluations) {
        QL_REQUIRE(maxEvaluations >= 15,
                   "required maxEvaluations (" << maxEvaluations <<
                   ") not allowed. It must be >= 15");
    }
}