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.numbers.core.ArithmeticUtils;
20 import org.apache.commons.math4.legacy.exception.NumberIsTooLargeException;
21 import org.apache.commons.math4.core.jdkmath.JdkMath;
22
23 /**
24 * Implements the <a href="https://en.wikipedia.org/wiki/Riemann_sum#Midpoint_rule">
25 * Midpoint Rule</a> for integration of real univariate functions. For
26 * reference, see <b>Numerical Mathematics</b>, ISBN 0387989595,
27 * chapter 9.2.
28 * <p>
29 * The function should be integrable.</p>
30 *
31 * @since 3.3
32 */
33 public class MidPointIntegrator extends BaseAbstractUnivariateIntegrator {
34
35 /** Maximum number of iterations for midpoint. 39 = floor(log_3(2^63)), the
36 * maximum number of triplings allowed before exceeding 64-bit bounds.
37 */
38 private static final int MIDPOINT_MAX_ITERATIONS_COUNT = 39;
39
40 /**
41 * Build a midpoint integrator with given accuracies and iterations counts.
42 * @param relativeAccuracy relative accuracy of the result
43 * @param absoluteAccuracy absolute accuracy of the result
44 * @param minimalIterationCount minimum number of iterations
45 * @param maximalIterationCount maximum number of iterations
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 39.
52 */
53 public MidPointIntegrator(final double relativeAccuracy,
54 final double absoluteAccuracy,
55 final int minimalIterationCount,
56 final int maximalIterationCount) {
57 super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
58 if (maximalIterationCount > MIDPOINT_MAX_ITERATIONS_COUNT) {
59 throw new NumberIsTooLargeException(maximalIterationCount,
60 MIDPOINT_MAX_ITERATIONS_COUNT, false);
61 }
62 }
63
64 /**
65 * Build a midpoint integrator with given iteration counts.
66 * @param minimalIterationCount minimum number of iterations
67 * @param maximalIterationCount maximum number of iterations
68 * @exception org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException if minimal number of iterations
69 * is not strictly positive
70 * @exception org.apache.commons.math4.legacy.exception.NumberIsTooSmallException if maximal number of iterations
71 * is lesser than or equal to the minimal number of iterations
72 * @exception NumberIsTooLargeException if maximal number of iterations
73 * is greater than 39.
74 */
75 public MidPointIntegrator(final int minimalIterationCount,
76 final int maximalIterationCount) {
77 super(minimalIterationCount, maximalIterationCount);
78 if (maximalIterationCount > MIDPOINT_MAX_ITERATIONS_COUNT) {
79 throw new NumberIsTooLargeException(maximalIterationCount,
80 MIDPOINT_MAX_ITERATIONS_COUNT, false);
81 }
82 }
83
84 /**
85 * Construct a midpoint integrator with default settings.
86 * (max iteration count set to {@link #MIDPOINT_MAX_ITERATIONS_COUNT})
87 */
88 public MidPointIntegrator() {
89 super(DEFAULT_MIN_ITERATIONS_COUNT, MIDPOINT_MAX_ITERATIONS_COUNT);
90 }
91
92 /**
93 * Compute the n-th stage integral of midpoint rule.
94 * This function should only be called by API <code>integrate()</code> in the package.
95 * To save time it does not verify arguments - caller does.
96 * <p>
97 * The interval is divided equally into 3^n sections rather than an
98 * arbitrary m sections because this configuration can best utilize the
99 * already computed values.</p>
100 *
101 * @param n the stage of 1/3 refinement. Must be larger than 0.
102 * @param previousStageResult Result from the previous call to the
103 * {@code stage} method.
104 * @param min Lower bound of the integration interval.
105 * @param diffMaxMin Difference between the lower bound and upper bound
106 * of the integration interval.
107 * @return the value of n-th stage integral
108 * @throws org.apache.commons.math4.legacy.exception.TooManyEvaluationsException if the maximal number of evaluations
109 * is exceeded.
110 */
111 private double stage(final int n,
112 double previousStageResult,
113 double min,
114 double diffMaxMin) {
115 // number of points in the previous stage. This stage will contribute
116 // 2*3^{n-1} more points.
117 final long np = ArithmeticUtils.pow(3L, n - 1);
118 double sum = 0;
119
120 // spacing between adjacent new points
121 final double spacing = diffMaxMin / np;
122 final double leftOffset = spacing / 6;
123 final double rightOffset = 5 * leftOffset;
124
125 double x = min;
126 for (long i = 0; i < np; i++) {
127 // The first and second new points are located at the new midpoints
128 // generated when each previous integration slice is split into 3.
129 //
130 // |--------x--------|
131 // |--x--|--x--|--x--|
132 sum += computeObjectiveValue(x + leftOffset);
133 sum += computeObjectiveValue(x + rightOffset);
134 x += spacing;
135 }
136 // add the new sum to previously calculated result
137 return (previousStageResult + sum * spacing) / 3.0;
138 }
139
140
141 /** {@inheritDoc} */
142 @Override
143 protected double doIntegrate() {
144 final double min = getMin();
145 final double diff = getMax() - min;
146 final double midPoint = min + 0.5 * diff;
147
148 double oldt = diff * computeObjectiveValue(midPoint);
149
150 while (true) {
151 iterations.increment();
152 final int i = iterations.getCount();
153 final double t = stage(i, oldt, min, diff);
154 if (i >= getMinimalIterationCount()) {
155 final double delta = JdkMath.abs(t - oldt);
156 final double rLimit =
157 getRelativeAccuracy() * (JdkMath.abs(oldt) + JdkMath.abs(t)) * 0.5;
158 if (delta <= rLimit || delta <= getAbsoluteAccuracy()) {
159 return t;
160 }
161 }
162 oldt = t;
163 }
164 }
165 }