001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 package org.apache.commons.math.analysis.integration; 018 019 import org.apache.commons.math.FunctionEvaluationException; 020 import org.apache.commons.math.MathRuntimeException; 021 import org.apache.commons.math.MaxIterationsExceededException; 022 import org.apache.commons.math.analysis.UnivariateRealFunction; 023 import org.apache.commons.math.exception.util.LocalizedFormats; 024 import org.apache.commons.math.util.FastMath; 025 026 /** 027 * Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html"> 028 * Romberg Algorithm</a> for integration of real univariate functions. For 029 * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, 030 * chapter 3. 031 * <p> 032 * Romberg integration employs k successive refinements of the trapezoid 033 * rule to remove error terms less than order O(N^(-2k)). Simpson's rule 034 * is a special case of k = 2.</p> 035 * 036 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $ 037 * @since 1.2 038 */ 039 public class RombergIntegrator extends UnivariateRealIntegratorImpl { 040 041 /** 042 * Construct an integrator for the given function. 043 * 044 * @param f function to integrate 045 * @deprecated as of 2.0 the integrand function is passed as an argument 046 * to the {@link #integrate(UnivariateRealFunction, double, double)}method. 047 */ 048 @Deprecated 049 public RombergIntegrator(UnivariateRealFunction f) { 050 super(f, 32); 051 } 052 053 /** 054 * Construct an integrator. 055 */ 056 public RombergIntegrator() { 057 super(32); 058 } 059 060 /** {@inheritDoc} */ 061 @Deprecated 062 public double integrate(final double min, final double max) 063 throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException { 064 return integrate(f, min, max); 065 } 066 067 /** {@inheritDoc} */ 068 public double integrate(final UnivariateRealFunction f, final double min, final double max) 069 throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException { 070 071 final int m = maximalIterationCount + 1; 072 double previousRow[] = new double[m]; 073 double currentRow[] = new double[m]; 074 075 clearResult(); 076 verifyInterval(min, max); 077 verifyIterationCount(); 078 079 TrapezoidIntegrator qtrap = new TrapezoidIntegrator(); 080 currentRow[0] = qtrap.stage(f, min, max, 0); 081 double olds = currentRow[0]; 082 for (int i = 1; i <= maximalIterationCount; ++i) { 083 084 // switch rows 085 final double[] tmpRow = previousRow; 086 previousRow = currentRow; 087 currentRow = tmpRow; 088 089 currentRow[0] = qtrap.stage(f, min, max, i); 090 for (int j = 1; j <= i; j++) { 091 // Richardson extrapolation coefficient 092 final double r = (1L << (2 * j)) - 1; 093 final double tIJm1 = currentRow[j - 1]; 094 currentRow[j] = tIJm1 + (tIJm1 - previousRow[j - 1]) / r; 095 } 096 final double s = currentRow[i]; 097 if (i >= minimalIterationCount) { 098 final double delta = FastMath.abs(s - olds); 099 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 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 }