View Javadoc

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  
24  /**
25   * Implements the <a href="http://mathworld.wolfram.com/SimpsonsRule.html">
26   * Simpson's Rule</a> for integration of real univariate functions. For
27   * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
28   * chapter 3.
29   * <p>
30   * This implementation employs basic trapezoid rule as building blocks to
31   * calculate the Simpson's rule of alternating 2/3 and 4/3.</p>
32   *  
33   * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
34   * @since 1.2
35   */
36  public class SimpsonIntegrator extends UnivariateRealIntegratorImpl {
37  
38      /**
39       * Construct an integrator for the given function.
40       * 
41       * @param f function to integrate
42       * @deprecated as of 2.0 the integrand function is passed as an argument
43       * to the {@link #integrate(UnivariateRealFunction, double, double)}method.
44       */
45      @Deprecated
46      public SimpsonIntegrator(UnivariateRealFunction f) {
47          super(f, 64);
48      }
49  
50      /**
51       * Construct an integrator.
52       */
53      public SimpsonIntegrator() {
54          super(64);
55      }
56  
57      /** {@inheritDoc} */
58      @Deprecated
59      public double integrate(final double min, final double max)
60          throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
61          return integrate(f, min, max);
62      }
63  
64      /** {@inheritDoc} */
65      public double integrate(final UnivariateRealFunction f,
66                              final double min, final double max)
67          throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
68          
69          int i = 1;
70          double s, olds, t, oldt;
71          
72          clearResult();
73          verifyInterval(min, max);
74          verifyIterationCount();
75  
76          TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
77          if (minimalIterationCount == 1) {
78              s = (4 * qtrap.stage(f, min, max, 1) - qtrap.stage(f, min, max, 0)) / 3.0;
79              setResult(s, 1);
80              return result;
81          }
82          // Simpson's rule requires at least two trapezoid stages.
83          olds = 0;
84          oldt = qtrap.stage(f, min, max, 0);
85          while (i <= maximalIterationCount) {
86              t = qtrap.stage(f, min, max, i);
87              s = (4 * t - oldt) / 3.0;
88              if (i >= minimalIterationCount) {
89                  final double delta = Math.abs(s - olds);
90                  final double rLimit =
91                      relativeAccuracy * (Math.abs(olds) + Math.abs(s)) * 0.5; 
92                  if ((delta <= rLimit) || (delta <= absoluteAccuracy)) {
93                      setResult(s, i);
94                      return result;
95                  }
96              }
97              olds = s;
98              oldt = t;
99              i++;
100         }
101         throw new MaxIterationsExceededException(maximalIterationCount);
102     }
103 
104     /** {@inheritDoc} */
105     @Override
106     protected void verifyIterationCount() throws IllegalArgumentException {
107         super.verifyIterationCount();
108         // at most 64 bisection refinements
109         if (maximalIterationCount > 64) {
110             throw MathRuntimeException.createIllegalArgumentException(
111                     "invalid iteration limits: min={0}, max={1}",
112                     0, 64);
113         }
114     }
115 }