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 }