Raymond | dee0849 | 2015-04-02 10:43:13 -0700 | [diff] [blame] | 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.math.analysis.integration; |
| 18 | |
| 19 | import org.apache.commons.math.FunctionEvaluationException; |
| 20 | import org.apache.commons.math.MathRuntimeException; |
| 21 | import org.apache.commons.math.MaxIterationsExceededException; |
| 22 | import org.apache.commons.math.analysis.UnivariateRealFunction; |
| 23 | import org.apache.commons.math.exception.util.LocalizedFormats; |
| 24 | import org.apache.commons.math.util.FastMath; |
| 25 | |
| 26 | /** |
| 27 | * Implements the <a href="http://mathworld.wolfram.com/TrapezoidalRule.html"> |
| 28 | * Trapezoidal Rule</a> for integration of real univariate functions. For |
| 29 | * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, |
| 30 | * chapter 3. |
| 31 | * <p> |
| 32 | * The function should be integrable.</p> |
| 33 | * |
| 34 | * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $ |
| 35 | * @since 1.2 |
| 36 | */ |
| 37 | public class TrapezoidIntegrator extends UnivariateRealIntegratorImpl { |
| 38 | |
| 39 | /** Intermediate result. */ |
| 40 | private double s; |
| 41 | |
| 42 | /** |
| 43 | * Construct an integrator for the given function. |
| 44 | * |
| 45 | * @param f function to integrate |
| 46 | * @deprecated as of 2.0 the integrand function is passed as an argument |
| 47 | * to the {@link #integrate(UnivariateRealFunction, double, double)}method. |
| 48 | */ |
| 49 | @Deprecated |
| 50 | public TrapezoidIntegrator(UnivariateRealFunction f) { |
| 51 | super(f, 64); |
| 52 | } |
| 53 | |
| 54 | /** |
| 55 | * Construct an integrator. |
| 56 | */ |
| 57 | public TrapezoidIntegrator() { |
| 58 | super(64); |
| 59 | } |
| 60 | |
| 61 | /** |
| 62 | * Compute the n-th stage integral of trapezoid rule. This function |
| 63 | * should only be called by API <code>integrate()</code> in the package. |
| 64 | * To save time it does not verify arguments - caller does. |
| 65 | * <p> |
| 66 | * The interval is divided equally into 2^n sections rather than an |
| 67 | * arbitrary m sections because this configuration can best utilize the |
| 68 | * alrealy computed values.</p> |
| 69 | * |
| 70 | * @param f the integrand function |
| 71 | * @param min the lower bound for the interval |
| 72 | * @param max the upper bound for the interval |
| 73 | * @param n the stage of 1/2 refinement, n = 0 is no refinement |
| 74 | * @return the value of n-th stage integral |
| 75 | * @throws FunctionEvaluationException if an error occurs evaluating the function |
| 76 | */ |
| 77 | double stage(final UnivariateRealFunction f, |
| 78 | final double min, final double max, final int n) |
| 79 | throws FunctionEvaluationException { |
| 80 | |
| 81 | if (n == 0) { |
| 82 | s = 0.5 * (max - min) * (f.value(min) + f.value(max)); |
| 83 | return s; |
| 84 | } else { |
| 85 | final long np = 1L << (n-1); // number of new points in this stage |
| 86 | double sum = 0; |
| 87 | final double spacing = (max - min) / np; // spacing between adjacent new points |
| 88 | double x = min + 0.5 * spacing; // the first new point |
| 89 | for (long i = 0; i < np; i++) { |
| 90 | sum += f.value(x); |
| 91 | x += spacing; |
| 92 | } |
| 93 | // add the new sum to previously calculated result |
| 94 | s = 0.5 * (s + sum * spacing); |
| 95 | return s; |
| 96 | } |
| 97 | } |
| 98 | |
| 99 | /** {@inheritDoc} */ |
| 100 | @Deprecated |
| 101 | public double integrate(final double min, final double max) |
| 102 | throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException { |
| 103 | return integrate(f, min, max); |
| 104 | } |
| 105 | |
| 106 | /** {@inheritDoc} */ |
| 107 | public double integrate(final UnivariateRealFunction f, final double min, final double max) |
| 108 | throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException { |
| 109 | |
| 110 | clearResult(); |
| 111 | verifyInterval(min, max); |
| 112 | verifyIterationCount(); |
| 113 | |
| 114 | double oldt = stage(f, min, max, 0); |
| 115 | for (int i = 1; i <= maximalIterationCount; ++i) { |
| 116 | final double t = stage(f, min, max, i); |
| 117 | if (i >= minimalIterationCount) { |
| 118 | final double delta = FastMath.abs(t - oldt); |
| 119 | final double rLimit = |
| 120 | relativeAccuracy * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5; |
| 121 | if ((delta <= rLimit) || (delta <= absoluteAccuracy)) { |
| 122 | setResult(t, i); |
| 123 | return result; |
| 124 | } |
| 125 | } |
| 126 | oldt = t; |
| 127 | } |
| 128 | throw new MaxIterationsExceededException(maximalIterationCount); |
| 129 | } |
| 130 | |
| 131 | /** {@inheritDoc} */ |
| 132 | @Override |
| 133 | protected void verifyIterationCount() throws IllegalArgumentException { |
| 134 | super.verifyIterationCount(); |
| 135 | // at most 64 bisection refinements |
| 136 | if (maximalIterationCount > 64) { |
| 137 | throw MathRuntimeException.createIllegalArgumentException( |
| 138 | LocalizedFormats.INVALID_ITERATIONS_LIMITS, |
| 139 | 0, 64); |
| 140 | } |
| 141 | } |
| 142 | } |