Home | History | Annotate | Download | only in integration 
      1 /*
      2  * Licensed to the Apache Software Foundation (ASF) under one or more
      3  * contributor license agreements.  See the NOTICE file distributed with
      4  * this work for additional information regarding copyright ownership.
      5  * The ASF licenses this file to You under the Apache License, Version 2.0
      6  * (the "License"); you may not use this file except in compliance with
      7  * the License.  You may obtain a copy of the License at
      8  *
      9  *      http://www.apache.org/licenses/LICENSE-2.0
     10  *
     11  * Unless required by applicable law or agreed to in writing, software
     12  * distributed under the License is distributed on an "AS IS" BASIS,
     13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
     14  * See the License for the specific language governing permissions and
     15  * limitations under the License.
     16  */
     17 package org.apache.commons.math.analysis.integration;
     18 
     19 import org.apache.commons.math.FunctionEvaluationException;
     20 import org.apache.commons.math.MathRuntimeException;
     21 import org.apache.commons.math.MaxIterationsExceededException;
     22 import org.apache.commons.math.analysis.UnivariateRealFunction;
     23 import org.apache.commons.math.exception.util.LocalizedFormats;
     24 import org.apache.commons.math.util.FastMath;
     25 
     26 /**
     27  * Implements the <a href="http://mathworld.wolfram.com/TrapezoidalRule.html">
     28  * Trapezoidal Rule</a> for integration of real univariate functions. For
     29  * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
     30  * chapter 3.
     31  * <p>
     32  * The function should be integrable.</p>
     33  *
     34  * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 fvr. 2011) $
     35  * @since 1.2
     36  */
     37 public class TrapezoidIntegrator extends UnivariateRealIntegratorImpl {
     38 
     39     /** Intermediate result. */
     40     private double s;
     41 
     42     /**
     43      * Construct an integrator for the given function.
     44      *
     45      * @param f function to integrate
     46      * @deprecated as of 2.0 the integrand function is passed as an argument
     47      * to the {@link #integrate(UnivariateRealFunction, double, double)}method.
     48      */
     49     @Deprecated
     50     public TrapezoidIntegrator(UnivariateRealFunction f) {
     51         super(f, 64);
     52     }
     53 
     54     /**
     55      * Construct an integrator.
     56      */
     57     public TrapezoidIntegrator() {
     58         super(64);
     59     }
     60 
     61     /**
     62      * Compute the n-th stage integral of trapezoid rule. This function
     63      * should only be called by API <code>integrate()</code> in the package.
     64      * To save time it does not verify arguments - caller does.
     65      * <p>
     66      * The interval is divided equally into 2^n sections rather than an
     67      * arbitrary m sections because this configuration can best utilize the
     68      * alrealy computed values.</p>
     69      *
     70      * @param f the integrand function
     71      * @param min the lower bound for the interval
     72      * @param max the upper bound for the interval
     73      * @param n the stage of 1/2 refinement, n = 0 is no refinement
     74      * @return the value of n-th stage integral
     75      * @throws FunctionEvaluationException if an error occurs evaluating the function
     76      */
     77     double stage(final UnivariateRealFunction f,
     78                  final double min, final double max, final int n)
     79         throws FunctionEvaluationException {
     80 
     81         if (n == 0) {
     82             s = 0.5 * (max - min) * (f.value(min) + f.value(max));
     83             return s;
     84         } else {
     85             final long np = 1L << (n-1);           // number of new points in this stage
     86             double sum = 0;
     87             final double spacing = (max - min) / np; // spacing between adjacent new points
     88             double x = min + 0.5 * spacing;    // the first new point
     89             for (long i = 0; i < np; i++) {
     90                 sum += f.value(x);
     91                 x += spacing;
     92             }
     93             // add the new sum to previously calculated result
     94             s = 0.5 * (s + sum * spacing);
     95             return s;
     96         }
     97     }
     98 
     99     /** {@inheritDoc} */
    100     @Deprecated
    101     public double integrate(final double min, final double max)
    102         throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
    103         return integrate(f, min, max);
    104     }
    105 
    106     /** {@inheritDoc} */
    107     public double integrate(final UnivariateRealFunction f, final double min, final double max)
    108         throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
    109 
    110         clearResult();
    111         verifyInterval(min, max);
    112         verifyIterationCount();
    113 
    114         double oldt = stage(f, min, max, 0);
    115         for (int i = 1; i <= maximalIterationCount; ++i) {
    116             final double t = stage(f, min, max, i);
    117             if (i >= minimalIterationCount) {
    118                 final double delta = FastMath.abs(t - oldt);
    119                 final double rLimit =
    120                     relativeAccuracy * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5;
    121                 if ((delta <= rLimit) || (delta <= absoluteAccuracy)) {
    122                     setResult(t, i);
    123                     return result;
    124                 }
    125             }
    126             oldt = t;
    127         }
    128         throw new MaxIterationsExceededException(maximalIterationCount);
    129     }
    130 
    131     /** {@inheritDoc} */
    132     @Override
    133     protected void verifyIterationCount() throws IllegalArgumentException {
    134         super.verifyIterationCount();
    135         // at most 64 bisection refinements
    136         if (maximalIterationCount > 64) {
    137             throw MathRuntimeException.createIllegalArgumentException(
    138                     LocalizedFormats.INVALID_ITERATIONS_LIMITS,
    139                     0, 64);
    140         }
    141     }
    142 }
    143