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 the <a href="http://mathworld.wolfram.com/RombergIntegration.html">
31   * Romberg Algorithm</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   * Romberg integration employs k successive refinements of the trapezoid
36   * rule to remove error terms less than order O(N^(-2k)). Simpson's rule
37   * is a special case of k = 2.</p>
38   *
39   */
40  public class RombergIntegrator extends BaseAbstractUnivariateIntegrator {
41  
42      /** Maximal number of iterations for Romberg. */
43      public static final int ROMBERG_MAX_ITERATIONS_COUNT = 32;
44  
45      /**
46       * Build a Romberg integrator with given accuracies and iterations counts.
47       * @param relativeAccuracy relative accuracy of the result
48       * @param absoluteAccuracy absolute accuracy of the result
49       * @param minimalIterationCount minimum number of iterations
50       * @param maximalIterationCount maximum number of iterations
51       * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
52       * @exception MathIllegalArgumentException if minimal number of iterations
53       * is not strictly positive
54       * @exception MathIllegalArgumentException if maximal number of iterations
55       * is lesser than or equal to the minimal number of iterations
56       * @exception MathIllegalArgumentException if maximal number of iterations
57       * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT}
58       */
59      public RombergIntegrator(final double relativeAccuracy,
60                               final double absoluteAccuracy,
61                               final int minimalIterationCount,
62                               final int maximalIterationCount)
63          throws MathIllegalArgumentException {
64          super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
65          if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) {
66              throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED,
67                                                     maximalIterationCount, ROMBERG_MAX_ITERATIONS_COUNT);
68          }
69      }
70  
71      /**
72       * Build a Romberg integrator with given iteration counts.
73       * @param minimalIterationCount minimum number of iterations
74       * @param maximalIterationCount maximum number of iterations
75       * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
76       * @exception MathIllegalArgumentException if minimal number of iterations
77       * is not strictly positive
78       * @exception MathIllegalArgumentException if maximal number of iterations
79       * is lesser than or equal to the minimal number of iterations
80       * @exception MathIllegalArgumentException if maximal number of iterations
81       * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT}
82       */
83      public RombergIntegrator(final int minimalIterationCount,
84                               final int maximalIterationCount)
85          throws MathIllegalArgumentException {
86          super(minimalIterationCount, maximalIterationCount);
87          if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) {
88              throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED,
89                                                     maximalIterationCount, ROMBERG_MAX_ITERATIONS_COUNT);
90          }
91      }
92  
93      /**
94       * Construct a Romberg integrator with default settings
95       * (max iteration count set to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
96       */
97      public RombergIntegrator() {
98          super(DEFAULT_MIN_ITERATIONS_COUNT, ROMBERG_MAX_ITERATIONS_COUNT);
99      }
100 
101     /** {@inheritDoc} */
102     @Override
103     protected double doIntegrate()
104         throws MathIllegalStateException {
105 
106         final int m = iterations.getMaximalCount() + 1;
107         double previousRow[] = new double[m];
108         double currentRow[]  = new double[m];
109 
110         TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
111         currentRow[0] = qtrap.stage(this, 0);
112         iterations.increment();
113         double olds = currentRow[0];
114         while (true) {
115 
116             final int i = iterations.getCount();
117 
118             // switch rows
119             final double[] tmpRow = previousRow;
120             previousRow = currentRow;
121             currentRow = tmpRow;
122 
123             currentRow[0] = qtrap.stage(this, i);
124             iterations.increment();
125             for (int j = 1; j <= i; j++) {
126                 // Richardson extrapolation coefficient
127                 final double r = (1L << (2 * j)) - 1;
128                 final double tIJm1 = currentRow[j - 1];
129                 currentRow[j] = tIJm1 + (tIJm1 - previousRow[j - 1]) / r;
130             }
131             final double s = currentRow[i];
132             if (i >= getMinimalIterationCount()) {
133                 final double delta  = FastMath.abs(s - olds);
134                 final double rLimit = getRelativeAccuracy() * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5;
135                 if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) {
136                     return s;
137                 }
138             }
139             olds = s;
140         }
141 
142     }
143 
144 }