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://mathworld.wolfram.com/TrapezoidalRule.html">
33 * Trapezoid Rule</a> for integration of real univariate functions. For
34 * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
35 * chapter 3.
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 FieldTrapezoidIntegrator<T extends CalculusFieldElement<T>> extends BaseAbstractFieldUnivariateIntegrator<T> {
42
43 /** Maximum number of iterations for trapezoid. */
44 public static final int TRAPEZOID_MAX_ITERATIONS_COUNT = 64;
45
46 /** Intermediate result. */
47 private T s;
48
49 /**
50 * Build a trapezoid 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 #TRAPEZOID_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 #TRAPEZOID_MAX_ITERATIONS_COUNT}
63 */
64 public FieldTrapezoidIntegrator(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 > TRAPEZOID_MAX_ITERATIONS_COUNT) {
72 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED,
73 maximalIterationCount, TRAPEZOID_MAX_ITERATIONS_COUNT);
74 }
75 }
76
77 /**
78 * Build a trapezoid 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 #TRAPEZOID_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 #TRAPEZOID_MAX_ITERATIONS_COUNT}
89 */
90 public FieldTrapezoidIntegrator(final Field<T> field,
91 final int minimalIterationCount,
92 final int maximalIterationCount)
93 throws MathIllegalArgumentException {
94 super(field, minimalIterationCount, maximalIterationCount);
95 if (maximalIterationCount > TRAPEZOID_MAX_ITERATIONS_COUNT) {
96 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED,
97 maximalIterationCount, TRAPEZOID_MAX_ITERATIONS_COUNT);
98 }
99 }
100
101 /**
102 * Construct a trapezoid integrator with default settings.
103 * @param field field to which function argument and value belong
104 * (max iteration count set to {@link #TRAPEZOID_MAX_ITERATIONS_COUNT})
105 */
106 public FieldTrapezoidIntegrator(final Field<T> field) {
107 super(field, DEFAULT_MIN_ITERATIONS_COUNT, TRAPEZOID_MAX_ITERATIONS_COUNT);
108 }
109
110 /**
111 * Compute the n-th stage integral of trapezoid rule. This function
112 * should only be called by API <code>integrate()</code> in the package.
113 * To save time it does not verify arguments - caller does.
114 * <p>
115 * The interval is divided equally into 2^n sections rather than an
116 * arbitrary m sections because this configuration can best utilize the
117 * already computed values.</p>
118 *
119 * @param baseIntegrator integrator holding integration parameters
120 * @param n the stage of 1/2 refinement, n = 0 is no refinement
121 * @return the value of n-th stage integral
122 * @throws MathIllegalStateException if the maximal number of evaluations
123 * is exceeded.
124 */
125 T stage(final BaseAbstractFieldUnivariateIntegrator<T> baseIntegrator, final int n)
126 throws MathIllegalStateException {
127
128 if (n == 0) {
129 final T max = baseIntegrator.getMax();
130 final T min = baseIntegrator.getMin();
131 s = max.subtract(min).multiply(0.5).
132 multiply(baseIntegrator.computeObjectiveValue(min).
133 add(baseIntegrator.computeObjectiveValue(max)));
134 return s;
135 } else {
136 final long np = 1L << (n-1); // number of new points in this stage
137 T sum = getField().getZero();
138 final T max = baseIntegrator.getMax();
139 final T min = baseIntegrator.getMin();
140 // spacing between adjacent new points
141 final T spacing = max.subtract(min).divide(np);
142 T x = min.add(spacing.multiply(0.5)); // the first new point
143 for (long i = 0; i < np; i++) {
144 sum = sum.add(baseIntegrator.computeObjectiveValue(x));
145 x = x.add(spacing);
146 }
147 // add the new sum to previously calculated result
148 s = s.add(sum.multiply(spacing)).multiply(0.5);
149 return s;
150 }
151 }
152
153 /** {@inheritDoc} */
154 @Override
155 protected T doIntegrate()
156 throws MathIllegalArgumentException, MathIllegalStateException {
157
158 T oldt = stage(this, 0);
159 iterations.increment();
160 while (true) {
161 final int i = iterations.getCount();
162 final T t = stage(this, i);
163 if (i >= getMinimalIterationCount()) {
164 final double delta = FastMath.abs(t.subtract(oldt)).getReal();
165 final double rlimit = FastMath.abs(oldt).add(FastMath.abs(t)).multiply(0.5 * getRelativeAccuracy()).getReal();
166 if ((delta <= rlimit) || (delta <= getAbsoluteAccuracy())) {
167 return t;
168 }
169 }
170 oldt = t;
171 iterations.increment();
172 }
173
174 }
175
176 }