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    *      https://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  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.analysis.integration;
23  
24  import org.hipparchus.exception.LocalizedCoreFormats;
25  import org.hipparchus.exception.MathIllegalArgumentException;
26  import org.hipparchus.exception.MathIllegalStateException;
27  import org.hipparchus.util.FastMath;
28  
29  /**
30   * Implements <a href="http://mathworld.wolfram.com/SimpsonsRule.html">
31   * Simpson's Rule</a> for integration of real univariate functions. For
32   * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
33   * chapter 3.
34   * <p>
35   * This implementation employs the basic trapezoid rule to calculate Simpson's
36   * rule.</p>
37   *
38   */
39  public class SimpsonIntegrator extends BaseAbstractUnivariateIntegrator {
40  
41      /** Maximal number of iterations for Simpson. */
42      public static final int SIMPSON_MAX_ITERATIONS_COUNT = 64;
43  
44      /**
45       * Build a Simpson integrator with given accuracies and iterations counts.
46       * @param relativeAccuracy relative accuracy of the result
47       * @param absoluteAccuracy absolute accuracy of the result
48       * @param minimalIterationCount minimum number of iterations
49       * @param maximalIterationCount maximum number of iterations
50       * (must be less than or equal to {@link #SIMPSON_MAX_ITERATIONS_COUNT})
51       * @exception MathIllegalArgumentException if minimal number of iterations
52       * is not strictly positive
53       * @exception MathIllegalArgumentException if maximal number of iterations
54       * is lesser than or equal to the minimal number of iterations
55       * @exception MathIllegalArgumentException if maximal number of iterations
56       * is greater than {@link #SIMPSON_MAX_ITERATIONS_COUNT}
57       */
58      public SimpsonIntegrator(final double relativeAccuracy,
59                               final double absoluteAccuracy,
60                               final int minimalIterationCount,
61                               final int maximalIterationCount)
62          throws MathIllegalArgumentException {
63          super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
64          if (maximalIterationCount > SIMPSON_MAX_ITERATIONS_COUNT) {
65              throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED,
66                                                     maximalIterationCount, SIMPSON_MAX_ITERATIONS_COUNT);
67          }
68      }
69  
70      /**
71       * Build a Simpson integrator with given iteration counts.
72       * @param minimalIterationCount minimum number of iterations
73       * @param maximalIterationCount maximum number of iterations
74       * (must be less than or equal to {@link #SIMPSON_MAX_ITERATIONS_COUNT})
75       * @exception MathIllegalArgumentException if minimal number of iterations
76       * is not strictly positive
77       * @exception MathIllegalArgumentException if maximal number of iterations
78       * is lesser than or equal to the minimal number of iterations
79       * @exception MathIllegalArgumentException if maximal number of iterations
80       * is greater than {@link #SIMPSON_MAX_ITERATIONS_COUNT}
81       */
82      public SimpsonIntegrator(final int minimalIterationCount,
83                               final int maximalIterationCount)
84          throws MathIllegalArgumentException {
85          super(minimalIterationCount, maximalIterationCount);
86          if (maximalIterationCount > SIMPSON_MAX_ITERATIONS_COUNT) {
87              throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED,
88                                                     maximalIterationCount, SIMPSON_MAX_ITERATIONS_COUNT);
89          }
90      }
91  
92      /**
93       * Construct an integrator with default settings.
94       * (max iteration count set to {@link #SIMPSON_MAX_ITERATIONS_COUNT})
95       */
96      public SimpsonIntegrator() {
97          super(DEFAULT_MIN_ITERATIONS_COUNT, SIMPSON_MAX_ITERATIONS_COUNT);
98      }
99  
100     /** {@inheritDoc} */
101     @Override
102     protected double doIntegrate()
103         throws MathIllegalStateException {
104 
105         TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
106         if (getMinimalIterationCount() == 1) {
107             return (4 * qtrap.stage(this, 1) - qtrap.stage(this, 0)) / 3.0;
108         }
109 
110         // Simpson's rule requires at least two trapezoid stages.
111         double olds = 0;
112         double oldt = qtrap.stage(this, 0);
113         while (true) {
114             final double t = qtrap.stage(this, iterations.getCount());
115             iterations.increment();
116             final double s = (4 * t - oldt) / 3.0;
117             if (iterations.getCount() >= getMinimalIterationCount()) {
118                 final double delta = FastMath.abs(s - olds);
119                 final double rLimit =
120                     getRelativeAccuracy() * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5;
121                 if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) {
122                     return s;
123                 }
124             }
125             olds = s;
126             oldt = t;
127         }
128 
129     }
130 
131 }