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.math.analysis.integration;
018
019 import org.apache.commons.math.FunctionEvaluationException;
020 import org.apache.commons.math.MathRuntimeException;
021 import org.apache.commons.math.MaxIterationsExceededException;
022 import org.apache.commons.math.analysis.UnivariateRealFunction;
023 import org.apache.commons.math.exception.util.LocalizedFormats;
024 import org.apache.commons.math.util.FastMath;
025
026 /**
027 * Implements the <a href="http://mathworld.wolfram.com/SimpsonsRule.html">
028 * Simpson's 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 * This implementation employs basic trapezoid rule as building blocks to
033 * calculate the Simpson's rule of alternating 2/3 and 4/3.</p>
034 *
035 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
036 * @since 1.2
037 */
038 public class SimpsonIntegrator extends UnivariateRealIntegratorImpl {
039
040 /**
041 * Construct an integrator for the given function.
042 *
043 * @param f function to integrate
044 * @deprecated as of 2.0 the integrand function is passed as an argument
045 * to the {@link #integrate(UnivariateRealFunction, double, double)}method.
046 */
047 @Deprecated
048 public SimpsonIntegrator(UnivariateRealFunction f) {
049 super(f, 64);
050 }
051
052 /**
053 * Construct an integrator.
054 */
055 public SimpsonIntegrator() {
056 super(64);
057 }
058
059 /** {@inheritDoc} */
060 @Deprecated
061 public double integrate(final double min, final double max)
062 throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
063 return integrate(f, min, max);
064 }
065
066 /** {@inheritDoc} */
067 public double integrate(final UnivariateRealFunction f, final double min, final double max)
068 throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
069
070 clearResult();
071 verifyInterval(min, max);
072 verifyIterationCount();
073
074 TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
075 if (minimalIterationCount == 1) {
076 final double s = (4 * qtrap.stage(f, min, max, 1) - qtrap.stage(f, min, max, 0)) / 3.0;
077 setResult(s, 1);
078 return result;
079 }
080 // Simpson's rule requires at least two trapezoid stages.
081 double olds = 0;
082 double oldt = qtrap.stage(f, min, max, 0);
083 for (int i = 1; i <= maximalIterationCount; ++i) {
084 final double t = qtrap.stage(f, min, max, i);
085 final double s = (4 * t - oldt) / 3.0;
086 if (i >= minimalIterationCount) {
087 final double delta = FastMath.abs(s - olds);
088 final double rLimit =
089 relativeAccuracy * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5;
090 if ((delta <= rLimit) || (delta <= absoluteAccuracy)) {
091 setResult(s, i);
092 return result;
093 }
094 }
095 olds = s;
096 oldt = t;
097 }
098 throw new MaxIterationsExceededException(maximalIterationCount);
099 }
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 }