001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math4.legacy.analysis.integration;
018
019import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
020import org.apache.commons.math4.legacy.analysis.integration.gauss.GaussIntegrator;
021import org.apache.commons.math4.legacy.analysis.integration.gauss.GaussIntegratorFactory;
022import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException;
023import org.apache.commons.math4.legacy.exception.TooManyEvaluationsException;
024import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
025import org.apache.commons.math4.core.jdkmath.JdkMath;
026
027/**
028 * This algorithm divides the integration interval into equally-sized
029 * sub-interval and on each of them performs a
030 * <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
031 * Legendre-Gauss</a> quadrature.
032 * Because of its <em>non-adaptive</em> nature, this algorithm can
033 * converge to a wrong value for the integral (for example, if the
034 * function is significantly different from zero toward the ends of the
035 * integration interval).
036 * In particular, a change of variables aimed at estimating integrals
037 * over infinite intervals as proposed
038 * <a href="http://en.wikipedia.org/w/index.php?title=Numerical_integration#Integrals_over_infinite_intervals">
039 *  here</a> should be avoided when using this class.
040 *
041 * @since 3.1
042 */
043
044public class IterativeLegendreGaussIntegrator
045    extends BaseAbstractUnivariateIntegrator {
046    /** Factory that computes the points and weights. */
047    private static final GaussIntegratorFactory FACTORY
048        = new GaussIntegratorFactory();
049    /** Number of integration points (per interval). */
050    private final int numberOfPoints;
051
052    /**
053     * Builds an integrator with given accuracies and iterations counts.
054     *
055     * @param n Number of integration points.
056     * @param relativeAccuracy Relative accuracy of the result.
057     * @param absoluteAccuracy Absolute accuracy of the result.
058     * @param minimalIterationCount Minimum number of iterations.
059     * @param maximalIterationCount Maximum number of iterations.
060     * @throws NotStrictlyPositiveException if minimal number of iterations
061     * or number of points are not strictly positive.
062     * @throws org.apache.commons.math4.legacy.exception.NumberIsTooSmallException
063     * if the maximal number of iterations is smaller than or equal to the
064     * minimal number of iterations.
065     */
066    public IterativeLegendreGaussIntegrator(final int n,
067                                            final double relativeAccuracy,
068                                            final double absoluteAccuracy,
069                                            final int minimalIterationCount,
070                                            final int maximalIterationCount) {
071        super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
072        if (n <= 0) {
073            throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_POINTS, n);
074        }
075       numberOfPoints = n;
076    }
077
078    /**
079     * Builds an integrator with given accuracies.
080     *
081     * @param n Number of integration points.
082     * @param relativeAccuracy Relative accuracy of the result.
083     * @param absoluteAccuracy Absolute accuracy of the result.
084     * @throws NotStrictlyPositiveException if {@code n < 1}.
085     */
086    public IterativeLegendreGaussIntegrator(final int n,
087                                            final double relativeAccuracy,
088                                            final double absoluteAccuracy) {
089        this(n, relativeAccuracy, absoluteAccuracy,
090             DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT);
091    }
092
093    /**
094     * Builds an integrator with given iteration counts.
095     *
096     * @param n Number of integration points.
097     * @param minimalIterationCount Minimum number of iterations.
098     * @param maximalIterationCount Maximum number of iterations.
099     * @throws NotStrictlyPositiveException if minimal number of iterations
100     * is not strictly positive.
101     * @throws org.apache.commons.math4.legacy.exception.NumberIsTooSmallException
102     * if the maximal number of iterations is smaller than or equal to the
103     * minimal number of iterations.
104     * @throws NotStrictlyPositiveException if {@code n < 1}.
105     */
106    public IterativeLegendreGaussIntegrator(final int n,
107                                            final int minimalIterationCount,
108                                            final int maximalIterationCount) {
109        this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY,
110             minimalIterationCount, maximalIterationCount);
111    }
112
113    /** {@inheritDoc} */
114    @Override
115    protected double doIntegrate() {
116        // Compute first estimate with a single step.
117        double oldt = stage(1);
118
119        int n = 2;
120        while (true) {
121            // Improve integral with a larger number of steps.
122            final double t = stage(n);
123
124            // Estimate the error.
125            final double delta = JdkMath.abs(t - oldt);
126            final double limit =
127                JdkMath.max(getAbsoluteAccuracy(),
128                             getRelativeAccuracy() * (JdkMath.abs(oldt) + JdkMath.abs(t)) * 0.5);
129
130            // check convergence
131            if (iterations.getCount() + 1 >= getMinimalIterationCount() &&
132                delta <= limit) {
133                return t;
134            }
135
136            // Prepare next iteration.
137            final double ratio = JdkMath.min(4, JdkMath.pow(delta / limit, 0.5 / numberOfPoints));
138            n = JdkMath.max((int) (ratio * n), n + 1);
139            oldt = t;
140            iterations.increment();
141        }
142    }
143
144    /**
145     * Compute the n-th stage integral.
146     *
147     * @param n Number of steps.
148     * @return the value of n-th stage integral.
149     * @throws TooManyEvaluationsException if the maximum number of evaluations
150     * is exceeded.
151     */
152    private double stage(final int n)
153        throws TooManyEvaluationsException {
154        // Function to be integrated is stored in the base class.
155        final UnivariateFunction f = new UnivariateFunction() {
156                /** {@inheritDoc} */
157                @Override
158                public double value(double x) {
159                    return computeObjectiveValue(x);
160                }
161            };
162
163        final double min = getMin();
164        final double max = getMax();
165        final double step = (max - min) / n;
166
167        double sum = 0;
168        for (int i = 0; i < n; i++) {
169            // Integrate over each sub-interval [a, b].
170            final double a = min + i * step;
171            final double b = a + step;
172            final GaussIntegrator g = FACTORY.legendreHighPrecision(numberOfPoints, a, b);
173            sum += g.integrate(f);
174        }
175
176        return sum;
177    }
178}