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.math4.legacy.analysis.integration;
18
19 import org.apache.commons.math4.legacy.exception.NumberIsTooLargeException;
20 import org.apache.commons.math4.core.jdkmath.JdkMath;
21
22 /**
23 * Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html">
24 * Romberg Algorithm</a> for integration of real univariate functions. For
25 * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
26 * chapter 3.
27 * <p>
28 * Romberg integration employs k successive refinements of the trapezoid
29 * rule to remove error terms less than order O(N^(-2k)). Simpson's rule
30 * is a special case of k = 2.</p>
31 *
32 * @since 1.2
33 */
34 public class RombergIntegrator extends BaseAbstractUnivariateIntegrator {
35
36 /** Maximal number of iterations for Romberg. */
37 public static final int ROMBERG_MAX_ITERATIONS_COUNT = 32;
38
39 /**
40 * Build a Romberg integrator with given accuracies and iterations counts.
41 * @param relativeAccuracy relative accuracy of the result
42 * @param absoluteAccuracy absolute accuracy of the result
43 * @param minimalIterationCount minimum number of iterations
44 * @param maximalIterationCount maximum number of iterations
45 * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
46 * @exception org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException if minimal number of iterations
47 * is not strictly positive
48 * @exception org.apache.commons.math4.legacy.exception.NumberIsTooSmallException if maximal number of iterations
49 * is lesser than or equal to the minimal number of iterations
50 * @exception NumberIsTooLargeException if maximal number of iterations
51 * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT}
52 */
53 public RombergIntegrator(final double relativeAccuracy,
54 final double absoluteAccuracy,
55 final int minimalIterationCount,
56 final int maximalIterationCount) {
57 super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
58 if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) {
59 throw new NumberIsTooLargeException(maximalIterationCount,
60 ROMBERG_MAX_ITERATIONS_COUNT, false);
61 }
62 }
63
64 /**
65 * Build a Romberg integrator with given iteration counts.
66 * @param minimalIterationCount minimum number of iterations
67 * @param maximalIterationCount maximum number of iterations
68 * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
69 * @exception org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException if minimal number of iterations
70 * is not strictly positive
71 * @exception org.apache.commons.math4.legacy.exception.NumberIsTooSmallException if maximal number of iterations
72 * is lesser than or equal to the minimal number of iterations
73 * @exception NumberIsTooLargeException if maximal number of iterations
74 * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT}
75 */
76 public RombergIntegrator(final int minimalIterationCount,
77 final int maximalIterationCount) {
78 super(minimalIterationCount, maximalIterationCount);
79 if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) {
80 throw new NumberIsTooLargeException(maximalIterationCount,
81 ROMBERG_MAX_ITERATIONS_COUNT, false);
82 }
83 }
84
85 /**
86 * Construct a Romberg integrator with default settings
87 * (max iteration count set to {@link #ROMBERG_MAX_ITERATIONS_COUNT}).
88 */
89 public RombergIntegrator() {
90 super(DEFAULT_MIN_ITERATIONS_COUNT, ROMBERG_MAX_ITERATIONS_COUNT);
91 }
92
93 /** {@inheritDoc} */
94 @Override
95 protected double doIntegrate() {
96 final int m = iterations.getMaximalCount() + 1;
97 double[] previousRow = new double[m];
98 double[] currentRow = new double[m];
99
100 TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
101 currentRow[0] = qtrap.stage(this, 0);
102 iterations.increment();
103 double olds = currentRow[0];
104 while (true) {
105
106 final int i = iterations.getCount();
107
108 // switch rows
109 final double[] tmpRow = previousRow;
110 previousRow = currentRow;
111 currentRow = tmpRow;
112
113 currentRow[0] = qtrap.stage(this, i);
114 iterations.increment();
115 for (int j = 1; j <= i; j++) {
116 // Richardson extrapolation coefficient
117 final double r = (1L << (2 * j)) - 1;
118 final double tIJm1 = currentRow[j - 1];
119 currentRow[j] = tIJm1 + (tIJm1 - previousRow[j - 1]) / r;
120 }
121 final double s = currentRow[i];
122 if (i >= getMinimalIterationCount()) {
123 final double delta = JdkMath.abs(s - olds);
124 final double rLimit = getRelativeAccuracy() * (JdkMath.abs(olds) + JdkMath.abs(s)) * 0.5;
125 if (delta <= rLimit || delta <= getAbsoluteAccuracy()) {
126 return s;
127 }
128 }
129 olds = s;
130 }
131 }
132 }