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/RombergIntegration.html"> 28 * Romberg Algorithm</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 * Romberg integration employs k successive refinements of the trapezoid 33 * rule to remove error terms less than order O(N^(-2k)). Simpson's rule 34 * is a special case of k = 2.</p> 35 * 36 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $ 37 * @since 1.2 38 */ 39 public class RombergIntegrator extends UnivariateRealIntegratorImpl { 40 41 /** 42 * Construct an integrator for the given function. 43 * 44 * @param f function to integrate 45 * @deprecated as of 2.0 the integrand function is passed as an argument 46 * to the {@link #integrate(UnivariateRealFunction, double, double)}method. 47 */ 48 @Deprecated RombergIntegrator(UnivariateRealFunction f)49 public RombergIntegrator(UnivariateRealFunction f) { 50 super(f, 32); 51 } 52 53 /** 54 * Construct an integrator. 55 */ RombergIntegrator()56 public RombergIntegrator() { 57 super(32); 58 } 59 60 /** {@inheritDoc} */ 61 @Deprecated integrate(final double min, final double max)62 public double integrate(final double min, final double max) 63 throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException { 64 return integrate(f, min, max); 65 } 66 67 /** {@inheritDoc} */ integrate(final UnivariateRealFunction f, final double min, final double max)68 public double integrate(final UnivariateRealFunction f, final double min, final double max) 69 throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException { 70 71 final int m = maximalIterationCount + 1; 72 double previousRow[] = new double[m]; 73 double currentRow[] = new double[m]; 74 75 clearResult(); 76 verifyInterval(min, max); 77 verifyIterationCount(); 78 79 TrapezoidIntegrator qtrap = new TrapezoidIntegrator(); 80 currentRow[0] = qtrap.stage(f, min, max, 0); 81 double olds = currentRow[0]; 82 for (int i = 1; i <= maximalIterationCount; ++i) { 83 84 // switch rows 85 final double[] tmpRow = previousRow; 86 previousRow = currentRow; 87 currentRow = tmpRow; 88 89 currentRow[0] = qtrap.stage(f, min, max, i); 90 for (int j = 1; j <= i; j++) { 91 // Richardson extrapolation coefficient 92 final double r = (1L << (2 * j)) - 1; 93 final double tIJm1 = currentRow[j - 1]; 94 currentRow[j] = tIJm1 + (tIJm1 - previousRow[j - 1]) / r; 95 } 96 final double s = currentRow[i]; 97 if (i >= minimalIterationCount) { 98 final double delta = FastMath.abs(s - olds); 99 final double rLimit = relativeAccuracy * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5; 100 if ((delta <= rLimit) || (delta <= absoluteAccuracy)) { 101 setResult(s, i); 102 return result; 103 } 104 } 105 olds = s; 106 } 107 throw new MaxIterationsExceededException(maximalIterationCount); 108 } 109 110 /** {@inheritDoc} */ 111 @Override verifyIterationCount()112 protected void verifyIterationCount() throws IllegalArgumentException { 113 super.verifyIterationCount(); 114 // at most 32 bisection refinements due to higher order divider 115 if (maximalIterationCount > 32) { 116 throw MathRuntimeException.createIllegalArgumentException( 117 LocalizedFormats.INVALID_ITERATIONS_LIMITS, 118 0, 32); 119 } 120 } 121 } 122