    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     */
    017    package org.apache.commons.math3.analysis.integration;
    018
    019    import org.apache.commons.math3.exception.MaxCountExceededException;
    020    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    021    import org.apache.commons.math3.exception.NumberIsTooLargeException;
    022    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    023    import org.apache.commons.math3.exception.TooManyEvaluationsException;
    024    import org.apache.commons.math3.util.FastMath;
    025
    026    /**
    027     * Implements the <a href="http://mathworld.wolfram.com/TrapezoidalRule.html">
    028     * Trapezoid Rule</a> for integration of real univariate functions. For
    029     * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
    030     * chapter 3.
    031     * <p>
    032     * The function should be integrable.</p>
    033     *
    034     * @version $Id: TrapezoidIntegrator.java 1364387 2012-07-22 18:14:11Z tn $
    035     * @since 1.2
    036     */
    037    public class TrapezoidIntegrator extends BaseAbstractUnivariateIntegrator {
    038
    039        /** Maximum number of iterations for trapezoid. */
    040        public static final int TRAPEZOID_MAX_ITERATIONS_COUNT = 64;
    041
    042        /** Intermediate result. */
    043        private double s;
    044
    045        /**
    046         * Build a trapezoid integrator with given accuracies and iterations counts.
    047         * @param relativeAccuracy relative accuracy of the result
    048         * @param absoluteAccuracy absolute accuracy of the result
    049         * @param minimalIterationCount minimum number of iterations
    050         * @param maximalIterationCount maximum number of iterations
    051         * (must be less than or equal to {@link #TRAPEZOID_MAX_ITERATIONS_COUNT}
    052         * @exception NotStrictlyPositiveException if minimal number of iterations
    053         * is not strictly positive
    054         * @exception NumberIsTooSmallException if maximal number of iterations
    055         * is lesser than or equal to the minimal number of iterations
    056         * @exception NumberIsTooLargeException if maximal number of iterations
    057         * is greater than {@link #TRAPEZOID_MAX_ITERATIONS_COUNT}
    058         */
    059        public TrapezoidIntegrator(final double relativeAccuracy,
    060                                   final double absoluteAccuracy,
    061                                   final int minimalIterationCount,
    062                                   final int maximalIterationCount)
    063            throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {
    064            super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
    065            if (maximalIterationCount > TRAPEZOID_MAX_ITERATIONS_COUNT) {
    066                throw new NumberIsTooLargeException(maximalIterationCount,
    067                                                    TRAPEZOID_MAX_ITERATIONS_COUNT, false);
    068            }
    069        }
    070
    071        /**
    072         * Build a trapezoid integrator with given iteration counts.
    073         * @param minimalIterationCount minimum number of iterations
    074         * @param maximalIterationCount maximum number of iterations
    075         * (must be less than or equal to {@link #TRAPEZOID_MAX_ITERATIONS_COUNT}
    076         * @exception NotStrictlyPositiveException if minimal number of iterations
    077         * is not strictly positive
    078         * @exception NumberIsTooSmallException if maximal number of iterations
    079         * is lesser than or equal to the minimal number of iterations
    080         * @exception NumberIsTooLargeException if maximal number of iterations
    081         * is greater than {@link #TRAPEZOID_MAX_ITERATIONS_COUNT}
    082         */
    083        public TrapezoidIntegrator(final int minimalIterationCount,
    084                                   final int maximalIterationCount)
    085            throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {
    086            super(minimalIterationCount, maximalIterationCount);
    087            if (maximalIterationCount > TRAPEZOID_MAX_ITERATIONS_COUNT) {
    088                throw new NumberIsTooLargeException(maximalIterationCount,
    089                                                    TRAPEZOID_MAX_ITERATIONS_COUNT, false);
    090            }
    091        }
    092
    093        /**
    094         * Construct a trapezoid integrator with default settings.
    095         * (max iteration count set to {@link #TRAPEZOID_MAX_ITERATIONS_COUNT})
    096         */
    097        public TrapezoidIntegrator() {
    098            super(DEFAULT_MIN_ITERATIONS_COUNT, TRAPEZOID_MAX_ITERATIONS_COUNT);
    099        }
    100
    101        /**
    102         * Compute the n-th stage integral of trapezoid rule. This function
    103         * should only be called by API <code>integrate()</code> in the package.
    104         * To save time it does not verify arguments - caller does.
    105         * <p>
    106         * The interval is divided equally into 2^n sections rather than an
    107         * arbitrary m sections because this configuration can best utilize the
    108         * already computed values.</p>
    109         *
    110         * @param baseIntegrator integrator holding integration parameters
    111         * @param n the stage of 1/2 refinement, n = 0 is no refinement
    112         * @return the value of n-th stage integral
    113         * @throws TooManyEvaluationsException if the maximal number of evaluations
    114         * is exceeded.
    115         */
    116        double stage(final BaseAbstractUnivariateIntegrator baseIntegrator, final int n)
    117            throws TooManyEvaluationsException {
    118
    119            if (n == 0) {
    120                final double max = baseIntegrator.getMax();
    121                final double min = baseIntegrator.getMin();
    122                s = 0.5 * (max - min) *
    123                          (baseIntegrator.computeObjectiveValue(min) +
    124                           baseIntegrator.computeObjectiveValue(max));
    125                return s;
    126            } else {
    127                final long np = 1L << (n-1);           // number of new points in this stage
    128                double sum = 0;
    129                final double max = baseIntegrator.getMax();
    130                final double min = baseIntegrator.getMin();
    131                // spacing between adjacent new points
    132                final double spacing = (max - min) / np;
    133                double x = min + 0.5 * spacing;    // the first new point
    134                for (long i = 0; i < np; i++) {
    135                    sum += baseIntegrator.computeObjectiveValue(x);
    136                    x += spacing;
    137                }
    138                // add the new sum to previously calculated result
    139                s = 0.5 * (s + sum * spacing);
    140                return s;
    141            }
    142        }
    143
    144        /** {@inheritDoc} */
    145        @Override
    146        protected double doIntegrate()
    147            throws TooManyEvaluationsException, MaxCountExceededException {
    148
    149            double oldt = stage(this, 0);
    150            iterations.incrementCount();
    151            while (true) {
    152                final int i = iterations.getCount();
    153                final double t = stage(this, i);
    154                if (i >= getMinimalIterationCount()) {
    155                    final double delta = FastMath.abs(t - oldt);
    156                    final double rLimit =
    157                        getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5;
    158                    if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) {
    159                        return t;
    160                    }
    161                }
    162                oldt = t;
    163                iterations.incrementCount();
    164            }
    165
    166        }
    167
    168    }