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 <a href="http://mathworld.wolfram.com/SimpsonsRule.html">
024 * Simpson's 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 * This implementation employs the basic trapezoid rule to calculate Simpson's
030 * rule.
031 *
032 * <p>
033 * <em>Caveat:</em> At each iteration, the algorithm refines the estimation by
034 * evaluating the function twice as many times as in the previous iteration;
035 * When specifying a {@link #integrate(int,UnivariateFunction,double,double)
036 * maximum number of function evaluations}, the caller must ensure that it
037 * is compatible with the {@link #SimpsonIntegrator(int,int) requested minimal
038 * number of iterations}.
039 *
040 * @since 1.2
041 */
042public class SimpsonIntegrator extends BaseAbstractUnivariateIntegrator {
043    /** Maximal number of iterations for Simpson. */
044    private static final int SIMPSON_MAX_ITERATIONS_COUNT = 30;
045
046    /**
047     * Build a Simpson integrator with given accuracies and iterations counts.
048     * @param relativeAccuracy relative accuracy of the result
049     * @param absoluteAccuracy absolute accuracy of the result
050     * @param minimalIterationCount Minimum number of iterations.
051     * @param maximalIterationCount Maximum number of iterations.
052     * It must be less than or equal to 30.
053     * @throws org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException
054     * if {@code minimalIterationCount <= 0}.
055     * @throws org.apache.commons.math4.legacy.exception.NumberIsTooSmallException
056     * if {@code maximalIterationCount < minimalIterationCount}.
057     * is lesser than or equal to the minimal number of iterations
058     * @throws NumberIsTooLargeException if {@code maximalIterationCount > 30}.
059     */
060    public SimpsonIntegrator(final double relativeAccuracy,
061                             final double absoluteAccuracy,
062                             final int minimalIterationCount,
063                             final int maximalIterationCount) {
064        super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
065        if (maximalIterationCount > SIMPSON_MAX_ITERATIONS_COUNT) {
066            throw new NumberIsTooLargeException(maximalIterationCount,
067                                                SIMPSON_MAX_ITERATIONS_COUNT, false);
068        }
069    }
070
071    /**
072     * Build a Simpson integrator with given iteration counts.
073     * @param minimalIterationCount Minimum number of iterations.
074     * @param maximalIterationCount Maximum number of iterations.
075     * It must be less than or equal to 30.
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 SimpsonIntegrator(final int minimalIterationCount,
084                             final int maximalIterationCount) {
085        super(minimalIterationCount, maximalIterationCount);
086        if (maximalIterationCount > SIMPSON_MAX_ITERATIONS_COUNT) {
087            throw new NumberIsTooLargeException(maximalIterationCount,
088                                                SIMPSON_MAX_ITERATIONS_COUNT, false);
089        }
090    }
091
092    /**
093     * Construct an integrator with default settings.
094     */
095    public SimpsonIntegrator() {
096        super(DEFAULT_MIN_ITERATIONS_COUNT, SIMPSON_MAX_ITERATIONS_COUNT);
097    }
098
099    /** {@inheritDoc} */
100    @Override
101    protected double doIntegrate() {
102        // Simpson's rule requires at least two trapezoid stages.
103        // So we set the first sum using two trapezoid stages.
104        final TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
105
106        final double s0 = qtrap.stage(this, 0);
107        double oldt = qtrap.stage(this, 1);
108        double olds = (4 * oldt - s0) / 3.0;
109        while (true) {
110            // The first iteration is the first refinement of the sum.
111            iterations.increment();
112            final int i = getIterations();
113            final double t = qtrap.stage(this, i + 1); // 1-stage ahead of the iteration
114            final double s = (4 * t - oldt) / 3.0;
115            if (i >= getMinimalIterationCount()) {
116                final double delta = JdkMath.abs(s - olds);
117                final double rLimit = getRelativeAccuracy() * (JdkMath.abs(olds) + JdkMath.abs(s)) * 0.5;
118                if (delta <= rLimit ||
119                    delta <= getAbsoluteAccuracy()) {
120                    return s;
121                }
122            }
123            olds = s;
124            oldt = t;
125        }
126    }
127}