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