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
31 /**
32 * Implements the <a href="http://en.wikipedia.org/wiki/Midpoint_method">
33 * Midpoint Rule</a> for integration of real univariate functions. For
34 * reference, see <b>Numerical Mathematics</b>, ISBN 0387989595,
35 * chapter 9.2.
36 * <p>
37 * The function should be integrable.</p>
38 * @param <T> Type of the field elements.
39 * @since 2.0
40 */
41 public class FieldMidPointIntegrator<T extends CalculusFieldElement<T>> extends BaseAbstractFieldUnivariateIntegrator<T> {
42
43 /** Maximum number of iterations for midpoint. */
44 public static final int MIDPOINT_MAX_ITERATIONS_COUNT = 64;
45
46 /**
47 * Build a midpoint integrator with given accuracies and iterations counts.
48 * @param field field to which function argument and value belong
49 * @param relativeAccuracy relative accuracy of the result
50 * @param absoluteAccuracy absolute accuracy of the result
51 * @param minimalIterationCount minimum number of iterations
52 * @param maximalIterationCount maximum number of iterations
53 * (must be less than or equal to {@link #MIDPOINT_MAX_ITERATIONS_COUNT}
54 * @exception MathIllegalArgumentException if minimal number of iterations
55 * is not strictly positive
56 * @exception MathIllegalArgumentException if maximal number of iterations
57 * is lesser than or equal to the minimal number of iterations
58 * @exception MathIllegalArgumentException if maximal number of iterations
59 * is greater than {@link #MIDPOINT_MAX_ITERATIONS_COUNT}
60 */
61 public FieldMidPointIntegrator(final Field<T> field,
62 final double relativeAccuracy,
63 final double absoluteAccuracy,
64 final int minimalIterationCount,
65 final int maximalIterationCount)
66 throws MathIllegalArgumentException {
67 super(field, relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
68 if (maximalIterationCount > MIDPOINT_MAX_ITERATIONS_COUNT) {
69 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED,
70 maximalIterationCount, MIDPOINT_MAX_ITERATIONS_COUNT);
71 }
72 }
73
74 /**
75 * Build a midpoint integrator with given iteration counts.
76 * @param field field to which function argument and value belong
77 * @param minimalIterationCount minimum number of iterations
78 * @param maximalIterationCount maximum number of iterations
79 * (must be less than or equal to {@link #MIDPOINT_MAX_ITERATIONS_COUNT}
80 * @exception MathIllegalArgumentException if minimal number of iterations
81 * is not strictly positive
82 * @exception MathIllegalArgumentException if maximal number of iterations
83 * is lesser than or equal to the minimal number of iterations
84 * @exception MathIllegalArgumentException if maximal number of iterations
85 * is greater than {@link #MIDPOINT_MAX_ITERATIONS_COUNT}
86 */
87 public FieldMidPointIntegrator(final Field<T> field,
88 final int minimalIterationCount,
89 final int maximalIterationCount)
90 throws MathIllegalArgumentException {
91 super(field, minimalIterationCount, maximalIterationCount);
92 if (maximalIterationCount > MIDPOINT_MAX_ITERATIONS_COUNT) {
93 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED,
94 maximalIterationCount, MIDPOINT_MAX_ITERATIONS_COUNT);
95 }
96 }
97
98 /**
99 * Construct a midpoint integrator with default settings.
100 * @param field field to which function argument and value belong
101 * (max iteration count set to {@link #MIDPOINT_MAX_ITERATIONS_COUNT})
102 */
103 public FieldMidPointIntegrator(final Field<T> field) {
104 super(field, DEFAULT_MIN_ITERATIONS_COUNT, MIDPOINT_MAX_ITERATIONS_COUNT);
105 }
106
107 /**
108 * Compute the n-th stage integral of midpoint rule.
109 * This function should only be called by API <code>integrate()</code> in the package.
110 * To save time it does not verify arguments - caller does.
111 * <p>
112 * The interval is divided equally into 2^n sections rather than an
113 * arbitrary m sections because this configuration can best utilize the
114 * already computed values.</p>
115 *
116 * @param n the stage of 1/2 refinement. Must be larger than 0.
117 * @param previousStageResult Result from the previous call to the
118 * {@code stage} method.
119 * @param min Lower bound of the integration interval.
120 * @param diffMaxMin Difference between the lower bound and upper bound
121 * of the integration interval.
122 * @return the value of n-th stage integral
123 * @throws MathIllegalStateException if the maximal number of evaluations
124 * is exceeded.
125 */
126 private T stage(final int n,
127 T previousStageResult,
128 T min,
129 T diffMaxMin)
130 throws MathIllegalStateException {
131
132 // number of new points in this stage
133 final long np = 1L << (n - 1);
134 T sum = getField().getZero();
135
136 // spacing between adjacent new points
137 final T spacing = diffMaxMin.divide(np);
138
139 // the first new point
140 T x = min.add(spacing.multiply(0.5));
141 for (long i = 0; i < np; i++) {
142 sum = sum.add(computeObjectiveValue(x));
143 x = x.add(spacing);
144 }
145 // add the new sum to previously calculated result
146 return previousStageResult.add(sum.multiply(spacing)).multiply(0.5);
147 }
148
149
150 /** {@inheritDoc} */
151 @Override
152 protected T doIntegrate()
153 throws MathIllegalArgumentException, MathIllegalStateException {
154
155 final T min = getMin();
156 final T diff = getMax().subtract(min);
157 final T midPoint = min.add(diff.multiply(0.5));
158
159 T oldt = diff.multiply(computeObjectiveValue(midPoint));
160
161 while (true) {
162 iterations.increment();
163 final int i = iterations.getCount();
164 final T t = stage(i, oldt, min, diff);
165 if (i >= getMinimalIterationCount()) {
166 final double delta = FastMath.abs(t.subtract(oldt)).getReal();
167 final double rLimit = FastMath.abs(oldt).add(FastMath.abs(t)).multiply(0.5 * getRelativeAccuracy()).getReal();
168 if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) {
169 return t;
170 }
171 }
172 oldt = t;
173 }
174
175 }
176
177 }