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/SimpsonsRule.html">
     28  * Simpson's 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  * This implementation employs basic trapezoid rule as building blocks to
     33  * calculate the Simpson's rule of alternating 2/3 and 4/3.</p>
     34  *
     35  * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 fvr. 2011) $
     36  * @since 1.2
     37  */
     38 public class SimpsonIntegrator extends UnivariateRealIntegratorImpl {
     39 
     40     /**
     41      * Construct an integrator for the given function.
     42      *
     43      * @param f function to integrate
     44      * @deprecated as of 2.0 the integrand function is passed as an argument
     45      * to the {@link #integrate(UnivariateRealFunction, double, double)}method.
     46      */
     47     @Deprecated
     48     public SimpsonIntegrator(UnivariateRealFunction f) {
     49         super(f, 64);
     50     }
     51 
     52     /**
     53      * Construct an integrator.
     54      */
     55     public SimpsonIntegrator() {
     56         super(64);
     57     }
     58 
     59     /** {@inheritDoc} */
     60     @Deprecated
     61     public double integrate(final double min, final double max)
     62         throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
     63         return integrate(f, min, max);
     64     }
     65 
     66     /** {@inheritDoc} */
     67     public double integrate(final UnivariateRealFunction f, final double min, final double max)
     68         throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
     69 
     70         clearResult();
     71         verifyInterval(min, max);
     72         verifyIterationCount();
     73 
     74         TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
     75         if (minimalIterationCount == 1) {
     76             final double s = (4 * qtrap.stage(f, min, max, 1) - qtrap.stage(f, min, max, 0)) / 3.0;
     77             setResult(s, 1);
     78             return result;
     79         }
     80         // Simpson's rule requires at least two trapezoid stages.
     81         double olds = 0;
     82         double oldt = qtrap.stage(f, min, max, 0);
     83         for (int i = 1; i <= maximalIterationCount; ++i) {
     84             final double t = qtrap.stage(f, min, max, i);
     85             final double s = (4 * t - oldt) / 3.0;
     86             if (i >= minimalIterationCount) {
     87                 final double delta = FastMath.abs(s - olds);
     88                 final double rLimit =
     89                     relativeAccuracy * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5;
     90                 if ((delta <= rLimit) || (delta <= absoluteAccuracy)) {
     91                     setResult(s, i);
     92                     return result;
     93                 }
     94             }
     95             olds = s;
     96             oldt = t;
     97         }
     98         throw new MaxIterationsExceededException(maximalIterationCount);
     99     }
    100 
    101     /** {@inheritDoc} */
    102     @Override
    103     protected void verifyIterationCount() throws IllegalArgumentException {
    104         super.verifyIterationCount();
    105         // at most 64 bisection refinements
    106         if (maximalIterationCount > 64) {
    107             throw MathRuntimeException.createIllegalArgumentException(
    108                     LocalizedFormats.INVALID_ITERATIONS_LIMITS,
    109                     0, 64);
    110         }
    111     }
    112 }
    113