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/TrapezoidalRule.html">
24 * Trapezoid Rule</a> for integration of real univariate functions.
25 *
26 * See <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, chapter 3.
27 *
28 * <p>
29 * The function should be integrable.
30 *
31 * <p>
32 * <em>Caveat:</em> At each iteration, the algorithm refines the estimation by
33 * evaluating the function twice as many times as in the previous iteration;
34 * When specifying a {@link #integrate(int,UnivariateFunction,double,double)
35 * maximum number of function evaluations}, the caller must ensure that it
36 * is compatible with the {@link #TrapezoidIntegrator(int,int) requested
37 * minimal number of iterations}.
38 *
39 * @since 1.2
40 */
41 public class TrapezoidIntegrator extends BaseAbstractUnivariateIntegrator {
42 /** Maximum number of iterations for trapezoid. */
43 private static final int TRAPEZOID_MAX_ITERATIONS_COUNT = 30;
44
45 /** Intermediate result. */
46 private double s;
47
48 /**
49 * Build a trapezoid integrator with given accuracies and iterations counts.
50 * @param relativeAccuracy relative accuracy of the result
51 * @param absoluteAccuracy absolute accuracy of the result
52 * @param minimalIterationCount minimum number of iterations
53 * @param maximalIterationCount maximum number of iterations
54 * @throws org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException
55 * if {@code minimalIterationCount <= 0}.
56 * @throws org.apache.commons.math4.legacy.exception.NumberIsTooSmallException
57 * if {@code maximalIterationCount < minimalIterationCount}.
58 * is lesser than or equal to the minimal number of iterations
59 * @throws NumberIsTooLargeException if {@code maximalIterationCount > 30}.
60 */
61 public TrapezoidIntegrator(final double relativeAccuracy,
62 final double absoluteAccuracy,
63 final int minimalIterationCount,
64 final int maximalIterationCount) {
65 super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
66 if (maximalIterationCount > TRAPEZOID_MAX_ITERATIONS_COUNT) {
67 throw new NumberIsTooLargeException(maximalIterationCount,
68 TRAPEZOID_MAX_ITERATIONS_COUNT, false);
69 }
70 }
71
72 /**
73 * Build a trapezoid integrator with given iteration counts.
74 * @param minimalIterationCount minimum number of iterations
75 * @param maximalIterationCount maximum number of iterations
76 * @throws org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException
77 * if {@code minimalIterationCount <= 0}.
78 * @throws org.apache.commons.math4.legacy.exception.NumberIsTooSmallException
79 * if {@code maximalIterationCount < minimalIterationCount}.
80 * is lesser than or equal to the minimal number of iterations
81 * @throws NumberIsTooLargeException if {@code maximalIterationCount > 30}.
82 */
83 public TrapezoidIntegrator(final int minimalIterationCount,
84 final int maximalIterationCount) {
85 super(minimalIterationCount, maximalIterationCount);
86 if (maximalIterationCount > TRAPEZOID_MAX_ITERATIONS_COUNT) {
87 throw new NumberIsTooLargeException(maximalIterationCount,
88 TRAPEZOID_MAX_ITERATIONS_COUNT, false);
89 }
90 }
91
92 /**
93 * Construct a trapezoid integrator with default settings.
94 */
95 public TrapezoidIntegrator() {
96 super(DEFAULT_MIN_ITERATIONS_COUNT, TRAPEZOID_MAX_ITERATIONS_COUNT);
97 }
98
99 /**
100 * Compute the n-th stage integral of trapezoid rule. This function
101 * should only be called by API <code>integrate()</code> in the package.
102 * To save time it does not verify arguments - caller does.
103 * <p>
104 * The interval is divided equally into 2^n sections rather than an
105 * arbitrary m sections because this configuration can best utilize the
106 * already computed values.</p>
107 *
108 * @param baseIntegrator integrator holding integration parameters
109 * @param n the stage of 1/2 refinement, n = 0 is no refinement
110 * @return the value of n-th stage integral
111 * @throws org.apache.commons.math4.legacy.exception.TooManyEvaluationsException if the maximal number of evaluations
112 * is exceeded.
113 */
114 double stage(final BaseAbstractUnivariateIntegrator baseIntegrator, final int n) {
115 if (n == 0) {
116 final double max = baseIntegrator.getMax();
117 final double min = baseIntegrator.getMin();
118 s = 0.5 * (max - min) *
119 (baseIntegrator.computeObjectiveValue(min) +
120 baseIntegrator.computeObjectiveValue(max));
121 return s;
122 } else {
123 final long np = 1L << (n-1); // number of new points in this stage
124 double sum = 0;
125 final double max = baseIntegrator.getMax();
126 final double min = baseIntegrator.getMin();
127 // spacing between adjacent new points
128 final double spacing = (max - min) / np;
129 double x = min + 0.5 * spacing; // the first new point
130 for (long i = 0; i < np; i++) {
131 sum += baseIntegrator.computeObjectiveValue(x);
132 x += spacing;
133 }
134 // add the new sum to previously calculated result
135 s = 0.5 * (s + sum * spacing);
136 return s;
137 }
138 }
139
140 /** {@inheritDoc} */
141 @Override
142 protected double doIntegrate() {
143 double oldt = stage(this, 0);
144 iterations.increment();
145 while (true) {
146 final int i = iterations.getCount();
147 final double t = stage(this, i);
148 if (i >= getMinimalIterationCount()) {
149 final double delta = JdkMath.abs(t - oldt);
150 final double rLimit =
151 getRelativeAccuracy() * (JdkMath.abs(oldt) + JdkMath.abs(t)) * 0.5;
152 if (delta <= rLimit || delta <= getAbsoluteAccuracy()) {
153 return t;
154 }
155 }
156 oldt = t;
157 iterations.increment();
158 }
159 }
160 }