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